GCC Code Coverage Report

Directory: ./
File: Model/PODModel.cpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 0 115 0.0%
Branches: 0 150 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 D.Lombardi
13 //
14
15 // System includes
16
17 // External includes
18
19 // Project includes
20 #include "Model/PODModel.hpp"
21
22 namespace felisce {
23 ✗ PODModel::PODModel():
24 ✗ Model() {
25 ✗ for(std::size_t i=0; i<m_eigenProblem.size(); i++) {
26 ✗ m_eigenProblem[i] = nullptr;
27 }
28 ✗ m_name = "POD";
29 }
30
31 ✗ PODModel::~PODModel() {
32 ✗ for(std::size_t i=0; i<m_eigenProblem.size(); i++) {
33 ✗ if (m_eigenProblem[i])
34 ✗ delete m_eigenProblem[i];
35 }
36 }
37
38 ✗ void PODModel::initializeEigenProblem(std::vector<EigenProblemPOD*> eigenPb) {
39 ✗ std::unordered_map<std::string, int> mapOfType;
40 ✗ mapOfType["EXPLICIT_EULER"] = 0;
41 ✗ mapOfType["IMPLICIT_EULER"] = 1;
42 ✗ mapOfType["SEMI_IMPLICIT_EULER"] = 2;
43
44 ✗ m_method = mapOfType[FelisceParam::instance().integrationTimeMethod];
45
46 ✗ for (std::size_t i=0; i<eigenPb.size(); i++) {
47 ✗ m_eigenProblem.push_back(eigenPb[i]);
48 }
49
50 ✗ for (std::size_t ipb = 0; ipb < eigenPb.size(); ipb++) {
51 //Define linear problem
52 //=======================
53 ✗ m_eigenProblem[ipb]->initialize(mesh(), m_fstransient, MpiInfo::petscComm());
54 ✗ if (FelisceParam::instance().solver.size() < eigenPb.size()) {
55 ✗ m_eigenProblem[ipb]->fixIdOfTheProblemSolver(0);
56 } else
57 ✗ m_eigenProblem[ipb]->fixIdOfTheProblemSolver(ipb);
58
59 //Compute Degrees of freedom (DOF) the problem
60 //=========================================
61 ✗ m_eigenProblem[ipb]->computeDof(MpiInfo::numProc(), MpiInfo::rankProc());
62 }
63 // Specific initalization of the user model
64 ✗ initializeDerivedModel();
65
66 //Define initial conditions
67 //=========================================
68 ✗ if (FelisceParam::instance().hasInitialCondition) {
69 ✗ for (std::size_t ipb = 0; ipb < eigenPb.size(); ipb++) {
70 ✗ for (std::size_t ii = 0; ii < FelisceParam::instance().valueInitCond.size(); ii++) {
71 ✗ m_initialCondition.addVariable(*m_eigenProblem[ipb]->listVariable().getVariable(FelisceParam::instance().nameVariableInitCond[ii]));
72 }
73 }
74 ✗ m_initialCondition.print(FelisceParam::verbose(),std::cout);
75 }
76
77 ✗ for (std::size_t ipb = 0; ipb < eigenPb.size(); ipb++) {
78 //Degrees of freedom partitionning with ParMetis
79 //===============================================
80 ✗ m_eigenProblem[ipb]->cutMesh();
81 //Allocate memory for the matrix _Matrix and std::vector _RHS in the linearProblem class
82 //============================================================
83
84 ✗ m_eigenProblem[ipb]->allocateMatrix();
85
86 }
87 ✗ setExternalVec(); // make the connection between the different linear pb (to be defined in the derived class)
88 ✗ std::cout << std::endl;
89 ✗ for (std::size_t ipb = 0; ipb < eigenPb.size(); ipb++) {
90 //Apply specific operations before assembling
91 //===========================================
92 ✗ preAssemblingMatrixRHS(ipb);
93
94 // Gather to write initial solution.
95 ✗ m_eigenProblem[ipb]->gatherSolution();
96
97 }
98 ✗ for(std::size_t iio=0; iio<m_ios.size(); ++iio)
99 ✗ m_ios[iio]->initializeOutput();
100 ✗ for(std::size_t ipb = 0; ipb < m_eigenProblem.size(); ipb++) {
101 ✗ m_eigenProblem[ipb]->clearMatrix();
102 }
103 ✗ setInitialCondition();
104
105 ✗ m_fstransient->iteration=0;
106
107 ✗ m_eigenProblemIsInitialized = true;
108 }
109
110
111 ✗ void PODModel::preAssembleMatrix(const int iProblem) {
112 ✗ if (FelisceParam::instance().hasInfarct) {
113 ✗ HeteroSparFHN heterof0;
114 ✗ std::vector<double> valuef0;
115 ✗ heterof0.initialize(m_fstransient);
116 ✗ m_eigenProblem[iProblem]->evalFunctionOnDof(heterof0, valuef0);
117 ✗ m_eigenProblem[iProblem]->setFhNf0(valuef0);
118 }
119 }
120
121 // Re-define
122 ✗ void PODModel::updateTime(const FlagMatrixRHS flagMatrixRHS) {
123 IGNORE_UNUSED_FLAG_MATRIX_RHS;
124 ✗ m_fstransient->iteration++;
125 ✗ m_fstransient->time +=m_fstransient->timeStep;
126 }
127
128 // Pay attention on the call of this :
129 // call PODModel::writeSolution() function instead of (non-virtual) Model::writeSolution() function
130 ✗ void PODModel::writeSolution() {
131 ✗ if (MpiInfo::rankProc() == 0) {
132 ✗ if (m_meshIsWritten == false) writeMesh();
133 }
134
135 ✗ if( (m_fstransient->iteration % FelisceParam::instance().frequencyWriteSolution == 0)
136 ✗ or m_hasFinished) {
137 ✗ if(FelisceParam::verbose() > 1) PetscPrintf(MpiInfo::petscComm(),"Write solutions\n");
138 ✗ for (std::size_t ipb = 0; ipb < m_eigenProblem.size(); ipb++) {
139 // Calculate FE solution (projection of alpha coeff)
140 ✗ m_eigenProblem[ipb]->writeEnsightSolution(static_cast<int>(m_fstransient->iteration/FelisceParam::instance().frequencyWriteSolution));
141 }
142 }
143 }
144
145 ✗ void PODModel::forward() {
146 ✗ for (std::size_t ipb = 0; ipb < m_eigenProblem.size(); ipb++) {
147 ✗ if ( m_fstransient->iteration == 0 ) {
148 // Initialize Rom object and calculate reduced basis
149
150 // Read initial data
151 ✗ m_eigenProblem[ipb]->readData(*io());
152 ✗ preAssembleMatrix(ipb);
153 // Create Matrix
154 ✗ m_eigenProblem[ipb]->assembleMatrix();
155
156 // Read basis from ensight files
157 ✗ if (FelisceParam::instance().readBasisFromFile) {
158 // Build basis reading vectors from ensight files
159 ✗ m_eigenProblem[ipb]->initializeROM();
160 }
161
162 else {
163 ✗ if (MpiInfo::rankProc() == 0) {
164 ✗ if (m_meshIsWritten == false) writeMesh();
165 }
166 // Initialize Slepc solver
167 ✗ m_eigenProblem[ipb]->buildSolver();
168 // Solve with slepc the generilized eigen problem: _Matrix[0] v = m_matrix[1] lambda v
169 ✗ m_eigenProblem[ipb]->solve();
170 // Writes modes (eigenvectors) in ensight format
171 ✗ m_eigenProblem[ipb]->writeMode();
172 ✗ m_eigenProblem[ipb]->initializeSolution();
173 }
174
175 ✗ if (FelisceParam::instance().hasSource)
176 ✗ postAssembleMatrix(ipb);
177
178 ✗ m_eigenProblem[ipb]->computeMatrixZeta();
179
180 ✗ m_eigenProblem[ipb]->computeTensor();
181
182 ✗ m_eigenProblem[ipb]->computeGamma();
183
184
185 ✗ if (FelisceParam::instance().writeECG) {
186 ✗ m_eigenProblem[ipb]->initializeECG();
187 ✗ m_eigenProblem[ipb]->writeECG(m_fstransient->iteration);
188 }
189
190 }
191 }
192
193 //Write solution with ensight.
194 ✗ PODModel::writeSolution();
195 //Advance time step.
196 ✗ updateTime();
197 ✗ printNewTimeIterationBanner();
198
199
200
201 ✗ for (std::size_t ipb = 0; ipb < m_eigenProblem.size(); ipb++) {
202 ✗ switch (m_method) {
203 ✗ case 0: // EXPLICIT_EULER
204
205 ✗ m_eigenProblem[ipb]->updateBeta();
206 ✗ m_eigenProblem[ipb]->computeGamma();
207
208
209 ✗ break;
210 ✗ case 1: // SEMI_IMPLICIT_EULER using fix point method
211 // updateTensor();
212 // semiImplicitFixPt();
213 ✗ break;
214 ✗ case 2: // IMPLICIT_EULER
215 // implicitFixPt();
216 ✗ break;
217 ✗ default:
218 ✗ FEL_ERROR("This integration method is not implemented.");
219 ✗ break;
220 }
221
222 ✗ if (FelisceParam::instance().writeECG) {
223 ✗ m_eigenProblem[ipb]->updateEcgOperator();
224 ✗ m_eigenProblem[ipb]->writeECG(m_fstransient->iteration);
225 }
226
227 }
228
229 }
230
231 ✗ void PODModel::solveEigenProblem() {
232 //Write solution with ensight.
233 ✗ if (MpiInfo::rankProc() == 0) {
234 ✗ if (m_meshIsWritten == false) writeMesh();
235 }
236 ✗ for (std::size_t ipb = 0; ipb < m_eigenProblem.size(); ipb++) {
237 // Read initial data
238 ✗ m_eigenProblem[ipb]->assembleMatrix();
239 // Read snapshots and compute autocorrelation:
240 ✗ m_eigenProblem[ipb]->readSnapshot(*io());
241 ✗ m_eigenProblem[ipb]->computeAutocorr();
242
243 // Initialize Slepc solver
244 ✗ m_eigenProblem[ipb]->buildSolver();
245
246
247 // Solve with slepc the generilized eigen problem: _Matrix[0] v = m_matrix[1] lambda v
248 ✗ m_eigenProblem[ipb]->solve();
249
250 // Writes modes (eigenvectors) in ensight format
251 ✗ m_eigenProblem[ipb]->writeMode();
252 }
253
254 ✗ for (std::size_t ipb = 0; ipb < m_eigenProblem.size(); ipb++)
255 ✗ m_eigenProblem[ipb]->checkOrthonormality();
256 // m_eigenProblem[ipb]->checkBasis();
257
258 }
259
260
261 }
262