GCC Code Coverage Report

Directory: ./
File: Model/ALPDeimModel.cpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 0 80 0.0%
Branches: 0 65 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/ALPDeimModel.hpp"
21
22 namespace felisce {
23 ✗ ALPDeimModel::ALPDeimModel():
24 ✗ ALPModel() {
25 ✗ for(std::size_t i=0; i<m_eigenProblemDeim.size(); i++) {
26 ✗ m_eigenProblemDeim[i] = nullptr;
27 }
28 }
29
30 ✗ void ALPDeimModel::preAssembleMatrix(const int iProblem) {
31
32 ✗ m_eigenProblemDeim.push_back(static_cast<EigenProblemALPDeim*>(m_eigenProblem[iProblem]));
33
34 ✗ std::unordered_map<std::string, int> mapOfType;
35 ✗ mapOfType["EXPLICIT_EULER"] = 0;
36 ✗ mapOfType["RUNGE_KUTTA_2"] = 1;
37
38 ✗ m_method = mapOfType[FelisceParam::instance().integrationTimeMethod];
39
40 ✗ m_eigenProblemDeim[iProblem]->setIntegrationMethod(m_method);
41
42 // Read collocation points from file (mesh nodes)
43 ✗ m_eigenProblemDeim[iProblem]->readCollocationPts();
44
45 ✗ if (FelisceParam::instance().hasInfarct) {
46 ✗ HeteroSparFHN heterof0;
47 ✗ std::vector<double> valuef0;
48 ✗ heterof0.initialize(m_fstransient);
49 ✗ m_eigenProblemDeim[iProblem]->evalFunctionOnDof(heterof0, valuef0);
50 ✗ m_eigenProblemDeim[iProblem]->setFhNf0(valuef0);
51 }
52 }
53
54 // Pay attention on the call of this :
55 // call ALPModel::writeSolution() function instead of (non-virtual) Model::writeSolution() function
56 ✗ void ALPDeimModel::writeSolution() {
57 ✗ if (MpiInfo::rankProc() == 0) {
58 ✗ if (m_meshIsWritten == false) writeMesh();
59 }
60
61 ✗ if( (m_fstransient->iteration % FelisceParam::instance().frequencyWriteSolution == 0)
62 ✗ or m_hasFinished) {
63 ✗ if(FelisceParam::verbose() > 1) PetscPrintf(MpiInfo::petscComm(),"Write solutions\n");
64 ✗ for (std::size_t ipb = 0; ipb < m_eigenProblem.size(); ipb++) {
65 // Calculate FE solution (projection of alpha coeff)
66 ✗ m_eigenProblemDeim[ipb]->writeEnsightSolution(static_cast<int>(m_fstransient->iteration/FelisceParam::instance().frequencyWriteSolution));
67 }
68 }
69 }
70
71 ✗ void ALPDeimModel::forward() {
72 ✗ for (std::size_t ipb = 0; ipb < m_eigenProblem.size(); ipb++) {
73 ✗ if ( m_fstransient->iteration == 0 ) {
74 // Initialize Rom object and calculate reduced basis
75
76 // Read initial data
77 ✗ m_eigenProblem[ipb]->readData(*io());
78 ✗ preAssembleMatrix(ipb);
79 // Create _Matrix[0] and _Matrix[1] of m_eigenProblem
80 ✗ m_eigenProblem[ipb]->assembleMatrix();
81
82 // Read basis from ensight files
83 ✗ if (FelisceParam::instance().readBasisFromFile) {
84 // Build basis reading vectors from ensight files
85 ✗ m_eigenProblemDeim[ipb]->initializeROM();
86 }
87 // Calculate basis functions solving Schrodinger equation
88 else {
89 ✗ if (MpiInfo::rankProc() == 0) {
90 ✗ if (m_meshIsWritten == false) writeMesh();
91 }
92 // Initialize Slepc solver
93 ✗ m_eigenProblem[ipb]->buildSolver();
94 // Solve with slepc the generilized eigen problem: _Matrix[0] v = m_matrix[1] lambda v
95 ✗ m_eigenProblem[ipb]->solve();
96 // Writes modes (eigenvectors) in ensight format
97 ✗ m_eigenProblem[ipb]->writeMode();
98 // m_massDeim = mass matrix on V
99 ✗ m_eigenProblemDeim[ipb]->computeMassDeim();
100 ✗ m_eigenProblemDeim[ipb]->initializeSolution();
101 }
102
103 ✗ if (FelisceParam::instance().hasSource)
104 ✗ postAssembleMatrix(ipb);
105
106 // m_hatF and M_hatG
107 ✗ m_eigenProblemDeim[ipb]->computeRHSDeim();
108
109 // m_theta = \Phi^T G F * \Phi
110 ✗ m_eigenProblemDeim[ipb]->computeTheta();
111
112 // M_ij = \chi / (\lambda_j - \lambda_i) \theta_ij
113 ✗ m_eigenProblemDeim[ipb]->computeMatrixM();
114
115 // matrices E and Q for improved recontruction and/or bidomain problem
116 ✗ m_eigenProblemDeim[ipb]->computeTensor();
117
118 // \Pi_V = V W G
119 ✗ m_eigenProblemDeim[ipb]->computeProjectionPiV();
120
121 ✗ if (FelisceParam::instance().writeECG) {
122 ✗ m_eigenProblemDeim[ipb]->initializeECG();
123 ✗ m_eigenProblemDeim[ipb]->writeECG(m_fstransient->iteration);
124 }
125
126 }
127 }
128
129 //Write solution with ensight.
130 ✗ ALPDeimModel::writeSolution();
131 //Advance time step.
132 ✗ updateTime();
133 ✗ printNewTimeIterationBanner();
134
135 ✗ for (std::size_t ipb = 0; ipb < m_eigenProblem.size(); ipb++) {
136 ✗ switch (m_method) {
137 ✗ case 0: // EXPLICIT_EULER
138 ✗ m_eigenProblemDeim[ipb]->updateBasis();
139
140 ✗ m_eigenProblemDeim[ipb]->updateBeta();
141
142 ✗ m_eigenProblemDeim[ipb]->updateEigenvalue();
143
144 ✗ m_eigenProblemDeim[ipb]->computeHatU();
145
146 ✗ m_eigenProblemDeim[ipb]->computeRHSDeim();
147 ✗ m_eigenProblemDeim[ipb]->computeTheta();
148 ✗ m_eigenProblemDeim[ipb]->computeMatrixM();
149
150 ✗ break;
151 ✗ case 1: // RUNGE_KUTTA_2
152 ✗ m_eigenProblemDeim[ipb]->updateBasis();
153
154 ✗ m_eigenProblemDeim[ipb]->updateBeta();
155
156 ✗ m_eigenProblemDeim[ipb]->updateEigenvalue();
157
158 ✗ m_eigenProblemDeim[ipb]->computeHatU();
159
160 ✗ m_eigenProblemDeim[ipb]->computeRHSDeim();
161 ✗ m_eigenProblemDeim[ipb]->computeTheta();
162 ✗ m_eigenProblemDeim[ipb]->computeMatrixM();
163 ✗ break;
164 ✗ default:
165 ✗ FEL_ERROR("This integration method is not implemented.");
166 ✗ break;
167 }
168
169 ✗ if (FelisceParam::instance().writeECG) {
170 ✗ m_eigenProblemDeim[ipb]->updateEcgOperator();
171 ✗ m_eigenProblemDeim[ipb]->writeECG(m_fstransient->iteration);
172 }
173
174 }
175
176 }
177
178
179 }
180