GCC Code Coverage Report

Directory: ./
File: Model/ALPModel.cpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 49 175 28.0%
Branches: 29 135 21.5%
Line Branch Exec Source
1 // ______ ______ _ _ _____ ______
2 // | ____| ____| | (_)/ ____| | ____|
3 // | |__ | |__ | | _| (___ ___| |__
4 // | __| | __| | | | |\___ \ / __| __|
5 // | | | |____| |____| |____) | (__| |____
6 // |_| |______|______|_|_____/ \___|______|
7 // Finite Elements for Life Sciences and Engineering
8 //
9 // License: LGL2.1 License
10 // FELiScE default license: LICENSE in root folder
11 //
12 // Main authors: E. Schenone
13 //
14
15 // System includes
16
17 // External includes
18
19 // Project includes
20 #include "Model/ALPModel.hpp"
21
22 namespace felisce {
23 1 ALPModel::ALPModel():
24 1 Model() {
25
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
 1 for(std::size_t i=0; i<m_eigenProblem.size(); i++) {
26 ✗ m_eigenProblem[i] = nullptr;
27 }
28
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 m_name = "ALP";
29 1 }
30
31 2 ALPModel::~ALPModel() {
32
2/2
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 1 times.
 4 for(std::size_t i=0; i<m_eigenProblem.size(); i++) {
33
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 2 if (m_eigenProblem[i])
34
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 2 delete m_eigenProblem[i];
35 }
36 }
37
38 1 void ALPModel::initializeEigenProblem(std::vector<EigenProblemALP*> eigenPb) {
39
40
2/2
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 1 times.
 2 for (std::size_t i=0; i<eigenPb.size(); i++) {
41 1 m_eigenProblem.push_back(eigenPb[i]);
42 }
43
44
2/2
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 1 times.
 2 for (std::size_t ipb = 0; ipb < eigenPb.size(); ipb++) {
45 //Define linear problem
46 //=======================
47
1/2
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
 1 m_eigenProblem[ipb]->initialize(mesh(), m_fstransient, MpiInfo::petscComm());
48
1/2
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
 1 if (FelisceParam::instance().solver.size() < eigenPb.size())
49 ✗ m_eigenProblem[ipb]->fixIdOfTheProblemSolver(0);
50 else
51 1 m_eigenProblem[ipb]->fixIdOfTheProblemSolver(ipb);
52
53 //Compute Degrees of freedom (DOF) the problem
54 //=========================================
55 1 m_eigenProblem[ipb]->computeDof(MpiInfo::numProc(), MpiInfo::rankProc());
56 }
57 // Specific initalization of the user model
58 1 initializeDerivedModel();
59
60 //Define initial conditions
61 //=========================================
62
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
 1 if (FelisceParam::instance().hasInitialCondition) {
63 ✗ for (std::size_t ipb = 0; ipb < eigenPb.size(); ipb++) {
64 ✗ for (std::size_t ii = 0; ii < FelisceParam::instance().valueInitCond.size(); ii++) {
65 ✗ m_initialCondition.addVariable(*m_eigenProblem[ipb]->listVariable().getVariable(FelisceParam::instance().nameVariableInitCond[ii]));
66 }
67 }
68 ✗ m_initialCondition.print(FelisceParam::verbose(),std::cout);
69 }
70
71
2/2
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 1 times.
 2 for (std::size_t ipb = 0; ipb < eigenPb.size(); ipb++) {
72 //Degrees of freedom partitionning with ParMetis
73 //===============================================
74 1 m_eigenProblem[ipb]->cutMesh();
75 //Allocate memory for the matrix _Matrix and std::vector _RHS in the linearProblem class
76 //============================================================
77
78 1 m_eigenProblem[ipb]->allocateMatrix();
79
80 }
81 1 setExternalVec(); // make the connection between the different linear pb (to be defined in the derived class)
82 1 std::cout << std::endl;
83
2/2
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 1 times.
 2 for (std::size_t ipb = 0; ipb < eigenPb.size(); ipb++) {
84 //Apply specific operations before assembling
85 //===========================================
86 1 preAssemblingMatrixRHS(ipb);
87
88 // Gather to write initial solution.
89 1 m_eigenProblem[ipb]->gatherSolution();
90
91 }
92
93
2/2
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 1 times.
 2 for(std::size_t iio=0; iio<m_ios.size(); ++iio)
94 1 m_ios[iio]->initializeOutput();
95
96
2/2
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 1 times.
 2 for(std::size_t ipb = 0; ipb < m_eigenProblem.size(); ipb++) {
97 1 m_eigenProblem[ipb]->clearMatrix();
98 }
99 1 setInitialCondition();
100
101 1 m_fstransient->iteration=0;
102
103 1 m_eigenProblemIsInitialized = true;
104 1 }
105
106
107 ✗ void ALPModel::preAssembleMatrix(const int iProblem) {
108 ✗ std::unordered_map<std::string, int> mapOfType;
109 ✗ mapOfType["EXPLICIT_EULER"] = 0;
110 ✗ mapOfType["EXPLICIT_EULER_MONOLITHIC"] = 1;
111 ✗ mapOfType["EXPLICIT_EULER_MONOLITHIC2"] = 2;
112 ✗ mapOfType["IMPLICIT_EULER"] = 3;
113 ✗ mapOfType["BACKWARD_DF_2"] = 4;
114
115 ✗ m_method = mapOfType[FelisceParam::instance().integrationTimeMethod];
116
117 ✗ m_eigenProblem[iProblem]->setIntegrationMethod(m_method);
118
119 ✗ if (FelisceParam::instance().hasInfarct) {
120 ✗ HeteroSparFHN heterof0;
121 ✗ std::vector<double> valuef0;
122 ✗ heterof0.initialize(m_fstransient);
123 ✗ m_eigenProblem[iProblem]->evalFunctionOnDof(heterof0, valuef0);
124 ✗ m_eigenProblem[iProblem]->setFhNf0(valuef0);
125 }
126 }
127
128 // Re-define
129 ✗ void ALPModel::updateTime(const FlagMatrixRHS flagMatrixRHS) {
130 IGNORE_UNUSED_FLAG_MATRIX_RHS;
131 ✗ m_fstransient->iteration++;
132 ✗ m_fstransient->time +=m_fstransient->timeStep;
133 }
134
135 // Pay attention on the call of this :
136 // call ALPModel::writeSolution() function instead of (non-virtual) Model::writeSolution() function
137 ✗ void ALPModel::writeSolution() {
138 ✗ if (MpiInfo::rankProc() == 0) {
139 ✗ if (m_meshIsWritten == false) writeMesh();
140 }
141
142 ✗ if( (m_fstransient->iteration % FelisceParam::instance().frequencyWriteSolution == 0)
143 ✗ or m_hasFinished) {
144 ✗ if(FelisceParam::verbose() > 1) PetscPrintf(MpiInfo::petscComm(),"Write solutions\n");
145 ✗ for (std::size_t ipb = 0; ipb < m_eigenProblem.size(); ipb++) {
146 // Calculate FE solution (projection of alpha coeff)
147 ✗ m_eigenProblem[ipb]->writeEnsightSolution(static_cast<int>(m_fstransient->iteration/FelisceParam::instance().frequencyWriteSolution));
148 }
149 }
150 }
151
152 ✗ void ALPModel::forward() {
153 ✗ for (std::size_t ipb = 0; ipb < m_eigenProblem.size(); ipb++) {
154 ✗ if ( m_fstransient->iteration == 0 ) {
155 // Initialize Rom object and calculate reduced basis
156
157 // Read initial data
158 ✗ m_eigenProblem[ipb]->readData(*io());
159 ✗ preAssembleMatrix(ipb);
160 // Create _Matrix[0] and _Matrix[1] of m_eigenProblem
161 ✗ m_eigenProblem[ipb]->assembleMatrix();
162
163 // Read basis from ensight files
164 ✗ if (FelisceParam::instance().readBasisFromFile) {
165 // Build basis reading vectors from ensight files
166 ✗ m_eigenProblem[ipb]->initializeROM();
167 }
168 // Calculate basis functions solving Schrodinger equation
169 else {
170 ✗ if (MpiInfo::rankProc() == 0) {
171 ✗ if (m_meshIsWritten == false) writeMesh();
172 }
173 // Initialize Slepc solver
174 ✗ m_eigenProblem[ipb]->buildSolver();
175 // Solve with slepc the generilized eigen problem: _Matrix[0] v = m_matrix[1] lambda v
176 ✗ m_eigenProblem[ipb]->solve();
177 // Writes modes (eigenvectors) in ensight format
178 ✗ m_eigenProblem[ipb]->writeMode();
179 ✗ m_eigenProblem[ipb]->initializeSolution();
180 }
181
182 ✗ if (FelisceParam::instance().hasSource)
183 ✗ postAssembleMatrix(ipb);
184
185 ✗ m_eigenProblem[ipb]->computeMatrixZeta();
186
187 ✗ m_eigenProblem[ipb]->computeTensor();
188
189 ✗ m_eigenProblem[ipb]->computeGamma();
190
191 ✗ m_eigenProblem[ipb]->computeMatrixM();
192
193 ✗ switch (m_method) {
194 ✗ case 0: // EXPLICIT_EULER
195 ✗ break;
196 ✗ case 1: // EXPLICIT_EULER_MONOLITHIC
197 case 2: // EXPLICIT_EULER_MONOLITHIC2
198 case 3: // IMPLICIT_EULER
199 case 4: // BACKWARD_DF_2
200 ✗ m_eigenProblem[ipb]->initializeSystemSolver();
201 ✗ break;
202 ✗ default:
203 ✗ FEL_ERROR("This integration method is not implemented.");
204 ✗ break;
205 }
206
207 ✗ if (FelisceParam::instance().writeECG) {
208 ✗ m_eigenProblem[ipb]->initializeECG();
209 ✗ m_eigenProblem[ipb]->writeECG(m_fstransient->iteration);
210 }
211
212 }
213 }
214
215 //Write solution with ensight.
216 ✗ ALPModel::writeSolution();
217 //Advance time step.
218 ✗ updateTime();
219 ✗ printNewTimeIterationBanner();
220
221
222
223 ✗ for (std::size_t ipb = 0; ipb < m_eigenProblem.size(); ipb++) {
224 ✗ switch (m_method) {
225 ✗ case 0: // EXPLICIT_EULER
226 ✗ m_eigenProblem[ipb]->updateBasis();
227
228 ✗ m_eigenProblem[ipb]->updateBeta();
229
230 ✗ m_eigenProblem[ipb]->updateEigenvalue();
231
232 ✗ m_eigenProblem[ipb]->updateTensor();
233
234 ✗ m_eigenProblem[ipb]->computeGamma();
235
236 // Compute Matrix M
237 ✗ m_eigenProblem[ipb]->computeMatrixM();
238
239 ✗ break;
240 ✗ case 1: // EXPLICIT_EULER_MONOLITHIC
241 ✗ m_eigenProblem[ipb]->updateBasis();
242
243 ✗ m_eigenProblem[ipb]->updateBetaMonolithic();
244
245 ✗ m_eigenProblem[ipb]->updateEigenvalue();
246
247 ✗ m_eigenProblem[ipb]->updateTensor();
248
249 ✗ m_eigenProblem[ipb]->computeGamma();
250
251 // Compute Matrix M
252 ✗ m_eigenProblem[ipb]->computeMatrixM();
253
254 ✗ break;
255 ✗ case 2: // EXPLICIT_EULER_MONOLITHIC2
256 ✗ m_eigenProblem[ipb]->updateBasis();
257
258 ✗ m_eigenProblem[ipb]->updateBetaMonolithic();
259
260 ✗ m_eigenProblem[ipb]->updateTensor();
261
262 ✗ m_eigenProblem[ipb]->updateEigenvalue();
263
264 ✗ m_eigenProblem[ipb]->computeGamma();
265
266 // Compute Matrix M
267 ✗ m_eigenProblem[ipb]->computeMatrixM();
268
269 ✗ break;
270 ✗ case 3: // IMPLICIT_EULER
271 ✗ m_eigenProblem[ipb]->updateBasis();
272
273 ✗ m_eigenProblem[ipb]->updateBetaEI();
274
275 ✗ m_eigenProblem[ipb]->updateTensor();
276
277 ✗ m_eigenProblem[ipb]->updateEigenvalue();
278
279 ✗ m_eigenProblem[ipb]->computeGamma();
280
281 // Compute Matrix M
282 ✗ m_eigenProblem[ipb]->computeMatrixM();
283
284 ✗ break;
285 ✗ case 4: // BACKWARD_DF_2
286 ✗ m_eigenProblem[ipb]->updateBasis();
287
288 ✗ m_eigenProblem[ipb]->computeGammaExtrap();
289
290 ✗ m_eigenProblem[ipb]->updateBetaBdf2();
291
292 ✗ m_eigenProblem[ipb]->updateTensor();
293
294 ✗ m_eigenProblem[ipb]->updateEigenvalue();
295
296 ✗ m_eigenProblem[ipb]->computeGamma();
297
298 // Compute Matrix M
299 ✗ m_eigenProblem[ipb]->computeMatrixM();
300
301 ✗ break;
302 ✗ default:
303 ✗ FEL_ERROR("This integration method is not implemented.");
304 ✗ break;
305 }
306
307 ✗ if (FelisceParam::instance().writeECG) {
308 ✗ m_eigenProblem[ipb]->updateEcgOperator();
309 ✗ m_eigenProblem[ipb]->writeECG(m_fstransient->iteration);
310 }
311
312 }
313
314 }
315
316 1 void ALPModel::solveEigenProblem() {
317 //Write solution with ensight.
318
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 if (MpiInfo::rankProc() == 0) {
319
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 if (m_meshIsWritten == false) writeMesh();
320 }
321
2/2
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 1 times.
 2 for (std::size_t ipb = 0; ipb < m_eigenProblem.size(); ipb++) {
322 // Read initial data
323 1 m_eigenProblem[ipb]->readData(*io());
324 // Create _Matrix[0] and _Matrix[1] of m_eigenProblem
325 1 m_eigenProblem[ipb]->assembleMatrix();
326 // Initialize Slepc solver
327 1 m_eigenProblem[ipb]->buildSolver();
328 // Solve with slepc the generilized eigen problem: _Matrix[0] v = m_matrix[1] lambda v
329 1 m_eigenProblem[ipb]->solve();
330
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 if (! FelisceParam::instance().optimizePotential) {
331 // Writes modes (eigenvectors) in ensight format
332 1 m_eigenProblem[ipb]->writeMode();
333 }
334 }
335
336
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 if (! FelisceParam::instance().optimizePotential) {
337
2/2
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 1 times.
 2 for (std::size_t ipb = 0; ipb < m_eigenProblem.size(); ipb++)
338 1 m_eigenProblem[ipb]->checkBasis();
339 }
340
341 1 }
342
343 ✗ void ALPModel::readInitialData()
344 {
345 ✗ for (std::size_t ipb = 0; ipb < m_eigenProblem.size(); ipb++)
346 {
347 ✗ m_eigenProblem[ipb]->readData(*io());
348
349 // Create _Matrix[0] and _Matrix[1] of m_eigenProblem
350 ✗ m_eigenProblem[ipb]->assembleMatrix();
351
352 // Read basis from ensight files
353 ✗ if (FelisceParam::instance().readBasisFromFile && FelisceParam::instance().optimizePotential) {
354 // Build basis reading vectors from ensight files
355 ✗ m_eigenProblem[ipb]->EigenProblemALP::readBasis();
356 }
357 }
358 }
359
360 ✗ void ALPModel::optimizePotential()
361 {
362 ✗ for (std::size_t ipb = 0; ipb < m_eigenProblem.size(); ipb++) {
363 ✗ m_eigenProblem[ipb]->setPotential();
364 // Writes modes (eigenvectors) in ensight format
365 ✗ m_eigenProblem[ipb]->writeMode();
366 }
367
368 }
369
370
371 }
372