GCC Code Coverage Report

Directory: ./
File: Model/ALPCurvModel.cpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 0 137 0.0%
Branches: 0 131 0.0%
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/ALPCurvModel.hpp"
21
22 namespace felisce {
23 ✗ ALPCurvModel::ALPCurvModel():
24 ✗ ALPModel()
25 {}
26
27 ✗ ALPCurvModel::~ALPCurvModel()
28 = default;
29
30 ✗ void ALPCurvModel::initializeEigenProblem(std::vector<EigenProblemALPCurv*> eigenPb) {
31 ✗ std::unordered_map<std::string, int> mapOfType;
32 ✗ mapOfType["EXPLICIT_EULER"] = 0;
33 ✗ mapOfType["EXPLICIT_EULER_MONOLITHIC"] = 1;
34 ✗ mapOfType["EXPLICIT_EULER_MONOLITHIC2"] = 2;
35 ✗ mapOfType["IMPLICIT_EULER"] = 3;
36 ✗ mapOfType["BACKWARD_DF_2"] = 4;
37
38 ✗ m_method = mapOfType[FelisceParam::instance().integrationTimeMethod];
39
40 ✗ for (std::size_t i=0; i<eigenPb.size(); i++) {
41 ✗ m_eigenProblem.push_back(eigenPb[i]);
42 }
43
44 ✗ for (std::size_t ipb = 0; ipb < eigenPb.size(); ipb++) {
45 //Define linear problem
46 //=======================
47 ✗ m_eigenProblem[ipb]->initialize(mesh(), m_fstransient, MpiInfo::petscComm());
48 ✗ if (FelisceParam::instance().solver.size() < eigenPb.size()) {
49 ✗ m_eigenProblem[ipb]->fixIdOfTheProblemSolver(0);
50 } else
51 ✗ m_eigenProblem[ipb]->fixIdOfTheProblemSolver(ipb);
52
53 //Compute Degrees of freedom (DOF) the problem
54 //=========================================
55 ✗ m_eigenProblem[ipb]->computeDof(MpiInfo::numProc(), MpiInfo::rankProc());
56 }
57 // Specific initalization of the user model
58 ✗ initializeDerivedModel();
59
60 //Define initial conditions
61 //=========================================
62 ✗ 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 ✗ for (std::size_t ipb = 0; ipb < eigenPb.size(); ipb++) {
72 //Degrees of freedom partitionning with ParMetis
73 //===============================================
74 ✗ m_eigenProblem[ipb]->cutMesh();
75 //Allocate memory for the matrix m_Matrix and std::vector m_RHS in the linearProblem class
76 //============================================================
77
78 ✗ m_eigenProblem[ipb]->allocateMatrix();
79
80 }
81 ✗ setExternalVec(); // make the connection between the different linear pb (to be defined in the derived class)
82 ✗ std::cout << std::endl;
83 ✗ for (std::size_t ipb = 0; ipb < eigenPb.size(); ipb++) {
84 //Apply specific operations before assembling
85 //===========================================
86 ✗ preAssemblingMatrixRHS(ipb);
87
88 // Gather to write initial solution.
89 ✗ m_eigenProblem[ipb]->gatherSolution();
90
91 }
92
93 ✗ for(std::size_t iio=0; iio<m_ios.size(); ++iio)
94 ✗ m_ios[iio]->initializeOutput();
95
96 ✗ for(std::size_t ipb = 0; ipb < m_eigenProblem.size(); ipb++) {
97 ✗ m_eigenProblem[ipb]->clearMatrix();
98 }
99 ✗ setInitialCondition();
100
101 ✗ for(std::size_t ipb = 0; ipb < m_eigenProblem.size(); ipb++) {
102 ✗ m_eigenProblem[ipb]->setIntegrationMethod(m_method);
103 }
104
105 ✗ m_fstransient->iteration=0;
106
107 ✗ m_eigenProblemIsInitialized = true;
108 }
109
110 ✗ void ALPCurvModel::forward() {
111 ✗ for (std::size_t ipb = 0; ipb < m_eigenProblem.size(); ipb++) {
112 ✗ if ( m_fstransient->iteration == 0 ) {
113 // Initialize Rom object and calculate reduced basis
114
115 // Read initial data
116 ✗ m_eigenProblem[ipb]->readData(*io());
117 ✗ preAssembleMatrix(ipb);
118 // Create _Matrix[0] and _Matrix[1] of m_eigenProblem
119 ✗ m_eigenProblem[ipb]->assembleMatrixBD();
120
121 // Read basis from ensight files
122 ✗ if (FelisceParam::instance().readBasisFromFile) {
123 // Build basis reading vectors from ensight files
124 ✗ m_eigenProblem[ipb]->initializeROM();
125 }
126 // Calculate basis functions solving Schrodinger equation
127 else {
128 ✗ if (MpiInfo::rankProc() == 0) {
129 ✗ if (m_meshIsWritten == false) writeMesh();
130 }
131 // Initialize Slepc solver
132 ✗ m_eigenProblem[ipb]->buildSolver();
133 // Solve with slepc the generilized eigen problem: _Matrix[0] v = m_matrix[1] lambda v
134 ✗ m_eigenProblem[ipb]->solve();
135 // Writes modes (eigenvectors) in ensight format
136 ✗ m_eigenProblem[ipb]->writeMode();
137 ✗ m_eigenProblem[ipb]->initializeSolution();
138 }
139
140 ✗ if (FelisceParam::instance().hasSource)
141 ✗ postAssembleMatrix(ipb);
142
143 ✗ m_eigenProblem[ipb]->computeMatrixZeta();
144
145 ✗ m_eigenProblem[ipb]->computeTensor();
146
147 ✗ m_eigenProblem[ipb]->computeGamma();
148
149 ✗ m_eigenProblem[ipb]->computeMatrixM();
150
151 ✗ switch (m_method) {
152 ✗ case 0: // EXPLICIT_EULER
153 ✗ break;
154 ✗ case 1: // EXPLICIT_EULER_MONOLITHIC
155 case 2: // EXPLICIT_EULER_MONOLITHIC2
156 case 3: // IMPLICIT_EULER
157 case 4: // BACKWARD_DF_2
158 ✗ m_eigenProblem[ipb]->initializeSystemSolver();
159 ✗ break;
160 ✗ default:
161 ✗ FEL_ERROR("This integration method is not implemented.");
162 ✗ break;
163 }
164
165 ✗ if (FelisceParam::instance().writeECG) {
166 ✗ m_eigenProblem[ipb]->initializeECG();
167 ✗ m_eigenProblem[ipb]->writeECG(m_fstransient->iteration);
168 }
169
170 }
171 }
172
173 //Write solution with ensight.
174 ✗ ALPModel::writeSolution();
175 //Advance time step.
176 ✗ updateTime();
177 ✗ printNewTimeIterationBanner();
178
179
180 ✗ for (std::size_t ipb = 0; ipb < m_eigenProblem.size(); ipb++) {
181 ✗ switch (m_method) {
182 ✗ case 0: // EXPLICIT_EULER
183 ✗ m_eigenProblem[ipb]->updateBasis();
184
185 //m_eigenProblem[ipb]->computeMatrixZeta();
186
187 ✗ m_eigenProblem[ipb]->updateBeta();
188
189 ✗ m_eigenProblem[ipb]->updateEigenvalue();
190
191 ✗ m_eigenProblem[ipb]->updateTensor();
192
193 ✗ m_eigenProblem[ipb]->computeGamma();
194
195 // Compute Matrix M
196 ✗ m_eigenProblem[ipb]->computeMatrixM();
197
198 ✗ break;
199 ✗ case 1: // EXPLICIT_EULER_MONOLITHIC
200 ✗ m_eigenProblem[ipb]->updateBasis();
201
202 //m_eigenProblem[ipb]->computeMatrixZeta();
203
204 ✗ m_eigenProblem[ipb]->updateBetaMonolithic();
205
206 ✗ m_eigenProblem[ipb]->updateEigenvalue();
207
208 ✗ m_eigenProblem[ipb]->updateTensor();
209
210 ✗ m_eigenProblem[ipb]->computeGamma();
211
212 // Compute Matrix M
213 ✗ m_eigenProblem[ipb]->computeMatrixM();
214
215 ✗ break;
216 ✗ case 2: // EXPLICIT_EULER_MONOLITHIC2
217 ✗ m_eigenProblem[ipb]->updateBasis();
218
219 //m_eigenProblem[ipb]->computeMatrixZeta();
220
221 ✗ m_eigenProblem[ipb]->updateBetaMonolithic();
222
223 ✗ m_eigenProblem[ipb]->updateTensor();
224
225 ✗ m_eigenProblem[ipb]->updateEigenvalue();
226
227 ✗ m_eigenProblem[ipb]->computeGamma();
228
229 // Compute Matrix M
230 ✗ m_eigenProblem[ipb]->computeMatrixM();
231
232 ✗ break;
233 ✗ case 3: // IMPLICIT_EULER
234 ✗ m_eigenProblem[ipb]->updateBasis();
235
236 //m_eigenProblem[ipb]->computeMatrixZeta();
237
238 ✗ m_eigenProblem[ipb]->updateBetaEI();
239
240 ✗ m_eigenProblem[ipb]->updateTensor();
241
242 ✗ m_eigenProblem[ipb]->updateEigenvalue();
243
244 ✗ m_eigenProblem[ipb]->computeGamma();
245
246 // Compute Matrix M
247 ✗ m_eigenProblem[ipb]->computeMatrixM();
248
249 ✗ break;
250 ✗ case 4: // BACKWARD_DF_2
251 ✗ m_eigenProblem[ipb]->updateBasis();
252
253 //m_eigenProblem[ipb]->computeMatrixZeta();
254
255 ✗ m_eigenProblem[ipb]->computeGammaExtrap();
256
257 ✗ m_eigenProblem[ipb]->updateBetaBdf2();
258
259 ✗ m_eigenProblem[ipb]->updateTensor();
260
261 ✗ m_eigenProblem[ipb]->updateEigenvalue();
262
263 ✗ m_eigenProblem[ipb]->computeGamma();
264
265 // Compute Matrix M
266 ✗ m_eigenProblem[ipb]->computeMatrixM();
267
268 ✗ break;
269 ✗ default:
270 ✗ FEL_ERROR("This integration method is not implemented.");
271 ✗ break;
272 }
273
274 ✗ if (FelisceParam::instance().writeECG) {
275 ✗ m_eigenProblem[ipb]->updateEcgOperator();
276 ✗ m_eigenProblem[ipb]->writeECG(m_fstransient->iteration);
277 }
278
279 }
280
281 }
282
283 ✗ void ALPCurvModel::solveEigenProblem() {
284 //Write solution with ensight.
285 ✗ if (MpiInfo::rankProc() == 0) {
286 ✗ if (m_meshIsWritten == false) writeMesh();
287 }
288 ✗ for (std::size_t ipb = 0; ipb < m_eigenProblem.size(); ipb++) {
289 // Read initial data
290 ✗ m_eigenProblem[ipb]->readData(*io());
291 // Create m_Matrix[0] and m_Matrix[1] of m_eigenProblem
292 ✗ m_eigenProblem[ipb]->assembleMatrixBD();
293 // Initialize Slepc solver
294 ✗ m_eigenProblem[ipb]->buildSolver();
295 // Solve with slepc the generilized eigen problem: m_Matrix[0] v = m_matrix[1] lambda v
296 ✗ m_eigenProblem[ipb]->solve();
297 // Writes modes (eigenvectors) in ensight format
298 ✗ m_eigenProblem[ipb]->writeMode();
299 }
300
301 ✗ for (std::size_t ipb = 0; ipb < m_eigenProblem.size(); ipb++)
302 ✗ m_eigenProblem[ipb]->checkBasis();
303
304 }
305
306
307 }
308