GCC Code Coverage Report

Directory: ./
File: Solver/eigenProblemALP.cpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 334 1961 17.0%
Branches: 292 3241 9.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 "Solver/eigenProblemALP.hpp"
21 #include "FiniteElement/elementMatrix.hpp"
22 #include "FiniteElement/elementField.hpp"
23 #include "Core/felisceParam.hpp"
24 #include "Core/felisceTools.hpp"
25
26 namespace felisce
27 {
28 /*! Construtor
29 */
30 1 EigenProblemALP::EigenProblemALP():
31
14/28
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 1 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 1 times.
✗ Branch 12 not taken.
✓ Branch 14 taken 1 times.
✗ Branch 15 not taken.
✓ Branch 17 taken 1 times.
✗ Branch 18 not taken.
✓ Branch 20 taken 1 times.
✗ Branch 21 not taken.
✓ Branch 23 taken 1 times.
✗ Branch 24 not taken.
✓ Branch 26 taken 1 times.
✗ Branch 27 not taken.
✓ Branch 29 taken 1 times.
✗ Branch 30 not taken.
✓ Branch 35 taken 1 times.
✗ Branch 36 not taken.
✓ Branch 38 taken 1 times.
✗ Branch 39 not taken.
✓ Branch 41 taken 1 times.
✗ Branch 42 not taken.
✓ Branch 44 taken 1 times.
✗ Branch 45 not taken.
 1 EigenProblem()
32 1 {}
33
34 2 EigenProblemALP::~EigenProblemALP() {
35 2 m_fePotTransMemb = nullptr;
36
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 2 delete [] m_eigenValue;
37
3/6
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
 2 if ( (FelisceParam::instance().ROMmethod == "ALP") || (FelisceParam::instance().solveEigenProblem) ) {
38
2/2
✓ Branch 0 taken 25 times.
✓ Branch 1 taken 1 times.
 52 for (int i=0; i<m_dimRomBasis; i++) {
39 50 m_basis[i].destroy();
40 }
41 }
42
2/6
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 1 times.
 2 if ( (m_useImprovedRec) && (FelisceParam::instance().ROMmethod == "ALP") ) {
43 ✗ for (int i=0; i<m_dimOrthComp; i++) {
44 ✗ m_orthComp[i].destroy();
45 }
46 }
47
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 2 if (m_matrixM) {
48 ✗ for (int i=0; i<m_dimRomBasis; i++) {
49 ✗ delete [] m_matrixM[i];
50 }
51 ✗ delete [] m_matrixM;
52 }
53
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
 2 if (m_zeta.size()) {
54 ✗ for (int i=0; i<m_size2; i++) {
55 ✗ m_zeta[i].destroy();
56 }
57 }
58
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 2 if (m_tensorB) delete [] m_tensorB;
59
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 2 if (m_tensorE) delete [] m_tensorE;
60
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 2 if (m_tensorEs) delete [] m_tensorEs;
61
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 2 if (m_tensorQ) delete [] m_tensorQ;
62
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 2 if (m_tensorT) delete [] m_tensorT;
63
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 2 if (m_tensorTs) delete [] m_tensorTs;
64
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 2 if (m_tensorU) {
65 ✗ for (int i=0; i<m_dimRomBasis; i++) {
66 ✗ delete [] m_tensorU[i];
67 }
68 ✗ delete [] m_tensorU;
69 }
70
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 2 if (m_tensorY) delete [] m_tensorY;
71 2 m_U_0.destroy();
72 2 m_U_0_seq.destroy();
73 2 m_W_0.destroy();
74 2 m_W_0_seq.destroy();
75 2 m_Ue_0.destroy();
76 2 m_Ue_0_seq.destroy();
77
2/6
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 1 times.
 2 if ( (FelisceParam::instance().hasInfarct) && (FelisceParam::instance().ROMmethod == "ALP") ) {
78 ✗ m_FhNf0.destroy();
79 }
80
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 2 if (m_solutionInitialized) {
81
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 2 delete [] m_beta;
82
2/4
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 2 if (m_mu) delete [] m_mu;
83
2/4
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 2 if (m_xi) delete [] m_xi;
84
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 2 if (m_gamma) delete [] m_gamma;
85
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 2 if (m_eta) delete [] m_eta;
86
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 2 if (m_method == 4) {
87 ✗ delete [] m_beta_n_1;
88 ✗ delete [] m_mu_n_1;
89 ✗ delete [] m_xi_n_1;
90 ✗ delete [] m_gamma_extrap;
91 ✗ delete [] m_eta_extrap;
92 }
93
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
 2 if (FelisceParam::instance().hasSource) {
94 ✗ for (int i=0; i<m_numSource; i++) {
95 ✗ delete [] m_source[i];
96 }
97 ✗ delete [] m_source;
98 }
99
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 2 delete [] m_modeMean;
100 }
101
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 2 if (m_fiber) {
102 ✗ delete [] m_fiber;
103 }
104
105
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
 2 if (FelisceParam::instance().optimizePotential) {
106 ✗ m_initPot.destroy();
107 }
108
109
3/10
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✓ Branch 12 taken 1 times.
 2 if ( (m_method != 0 ) && (m_method != -1 ) && (!FelisceParam::instance().solveEigenProblem) && (FelisceParam::instance().ROMmethod == "ALP")) {
110 ✗ m_solutionRom.destroy();
111 ✗ m_rhsRom.destroy();
112 }
113
114 }
115
116 1 void EigenProblemALP::initialize(const GeometricMeshRegion::Pointer& mesh, FelisceTransient::Pointer fstransient, MPI_Comm& comm) {
117
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 EigenProblem::initialize(mesh, fstransient, comm);
118
119 1 std::unordered_map<std::string, int> mapOfPbm;
120
2/4
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
 1 mapOfPbm["FKPP"] = 0;
121
2/4
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
 1 mapOfPbm["Monodomain"] = 1;
122
2/4
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
 1 mapOfPbm["Bidomain"] = 2;
123
2/4
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
 1 mapOfPbm["BidomainCurv"] = 2;
124
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 1 m_problem = mapOfPbm[FelisceParam::instance().model];
125
126
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 m_dimRomBasis = FelisceParam::instance().dimRomBasis;
127
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 m_useImprovedRec = FelisceParam::instance().useImprovedRec;
128
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 m_dimOrthComp = FelisceParam::instance().dimOrthComp;
129
130
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 m_tensorOrder = FelisceParam::instance().orderALP;
131
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 m_numSource = FelisceParam::instance().numberOfSource;
132
133 1 int id=0;
134 1 m_idL = id;
135 1 id++;
136 1 m_idG = id;
137 1 id++;
138
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
 1 if (!FelisceParam::instance().solveEigenProblem) {
139 ✗ m_idGs = id;
140 ✗ id++;
141 ✗ if (FelisceParam::instance().model == "Monodomain") {
142 ✗ if (m_tensorOrder == 3 ) {
143 ✗ m_idK = id;
144 ✗ id++;
145 ✗ m_idKs = id;
146 ✗ id++;
147 }
148 ✗ } else if ( (FelisceParam::instance().model == "Bidomain") || (FelisceParam::instance().model == "BidomainCurv") ) {
149 ✗ m_idK = id;
150 ✗ id++;
151 ✗ m_idKie = id;
152 ✗ id++;
153 ✗ if (m_tensorOrder == 3 ) {
154 ✗ m_idKs = id;
155 ✗ id++;
156 }
157 }
158
159 ✗ if ( (m_useImprovedRec) & (m_idK == 0) ) {
160 ✗ m_idK = id;
161 ✗ id++;
162 }
163 }
164 1 m_numberOfMatrix = id;
165
166
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 m_coefChi = FelisceParam::instance().chiSchrodinger;
167 1 std::vector<PhysicalVariable> listVariable;
168 1 std::vector<std::size_t> listNumComp;
169
1/4
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 switch (m_problem) {
170 ✗ case 0: {
171 ✗ nameOfTheProblem = "Problem Fisher-Kolmogorov-Petrovski-Piskunov equation (FKPP)";
172 ✗ listVariable.push_back(potTransMemb);
173 ✗ listNumComp.push_back(1);
174 //define unknown of the linear system.
175 ✗ m_listUnknown.push_back(potTransMemb);
176 ✗ break;
177 }
178 ✗ case 1: {
179 ✗ nameOfTheProblem = "Problem cardiac Monodomain";
180 ✗ listVariable.push_back(potTransMemb);
181 ✗ listNumComp.push_back(1);
182 //define unknown of the linear system.
183 ✗ m_listUnknown.push_back(potTransMemb);
184 ✗ break;
185 }
186 1 case 2: {
187
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 nameOfTheProblem = "Problem cardiac Bidomain";
188
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 listVariable.push_back(potTransMemb);
189
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 listNumComp.push_back(1);
190 //define unknown of the linear system.
191
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 m_listUnknown.push_back(potTransMemb);
192 1 break;
193 }
194 ✗ default:
195 ✗ FEL_ERROR("Model not defined for ALP solver.");
196 ✗ break;
197 }
198
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 definePhysicalVariable(listVariable,listNumComp);
199
200 1 }
201
202 ✗ void EigenProblemALP::initializeSystemSolver() {
203 ✗ if (m_problem != 2) {
204 ✗ FEL_ERROR("Method not implemented for this problem.");
205 }
206 ✗ felInt* nnz = new felInt[3*m_dimRomBasis+1];
207 ✗ for ( felInt i = 0; i<m_dimRomBasis; i++) {
208 ✗ if (m_method == 1 || m_method == 2) {
209 ✗ nnz[i] = 1;
210 ✗ nnz[i+m_dimRomBasis] = 1;
211 ✗ } else if (m_method == 3 || m_method == 4) {
212 ✗ nnz[i] = m_dimRomBasis;
213 ✗ nnz[i+m_dimRomBasis] = m_dimRomBasis;
214 }
215 ✗ nnz[i+2*m_dimRomBasis] = 2*m_dimRomBasis+1;
216 }
217 ✗ nnz[3*m_dimRomBasis] = m_dimRomBasis+1;
218 ✗ m_matrixRom.createSeqAIJ(PETSC_COMM_SELF, 3*m_dimRomBasis+1, 3*m_dimRomBasis+1, 0, nnz);
219 ✗ m_matrixRom.setFromOptions();
220 ✗ m_matrixRom.setOption(MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE);
221 ✗ delete [] nnz;
222
223 ✗ m_matrixRom.getVecs(m_solutionRom,m_rhsRom);
224
225 int ii;
226 ✗ for (int i=0; i<m_dimRomBasis; i++) {
227 ✗ ii=i;
228 ✗ m_solutionRom.setValues(1,&ii,&m_beta[i],INSERT_VALUES);
229 ✗ ii=i+m_dimRomBasis;
230 ✗ m_solutionRom.setValues(1,&ii,&m_mu[i],INSERT_VALUES);
231 ✗ ii=i+2*m_dimRomBasis;
232 ✗ m_solutionRom.setValues(1,&ii,&m_xi[i],INSERT_VALUES);
233 }
234 ✗ m_solutionRom.assembly();
235
236 ✗ if (FelisceParam::verbose() > 40)
237 ✗ m_solutionRom.view();
238
239 ✗ m_kspRom->init();
240 ✗ m_kspRom->setKSPandPCType(KSPGMRES,PCASM);
241 ✗ m_kspRom->setTolerances(1.e-10,1.e-10,1000,200);
242 ✗ m_kspRom->setKSPOptions(KSPGMRES,true);
243 }
244
245 // Define Physical Variable associate to the problem
246 1 void EigenProblemALP::initPerElementType() {
247
1/4
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 switch(m_problem) {
248 ✗ case 0: {
249 //Init pointer on Finite Element, Variable or idVariable
250 ✗ m_ipotTransMemb = m_listVariable.getVariableIdList(potTransMemb);
251 ✗ m_fePotTransMemb = m_listCurrentFiniteElement[m_ipotTransMemb];
252 ✗ m_elemFieldU0.initialize(DOF_FIELD,*m_fePotTransMemb);
253 ✗ break;
254 }
255 ✗ case 1: {
256 //Init pointer on Finite Element, Variable or idVariable
257 ✗ m_ipotTransMemb = m_listVariable.getVariableIdList(potTransMemb);
258 ✗ m_fePotTransMemb = m_listCurrentFiniteElement[m_ipotTransMemb];
259 ✗ m_elemFieldU0.initialize(DOF_FIELD,*m_fePotTransMemb);
260 ✗ break;
261 }
262 1 case 2: {
263 //Init pointer on Finite Element, Variable or idVariable
264 1 m_ipotTransMemb = m_listVariable.getVariableIdList(potTransMemb);
265 1 m_fePotTransMemb = m_listCurrentFiniteElement[m_ipotTransMemb];
266 1 m_elemFieldU0.initialize(DOF_FIELD,*m_fePotTransMemb);
267 1 break;
268 }
269 ✗ default:
270 ✗ FEL_ERROR("Model not defined for ALP solver.");
271 ✗ break;
272 }
273
274
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
 1 if (FelisceParam::instance().hasInfarct) {
275 ✗ m_elemFieldFhNf0.initialize(DOF_FIELD,*m_fePotTransMemb);
276 }
277 1 }
278
279 // Assemble Matrix
280 8192 void EigenProblemALP::computeElementArray(const std::vector<Point*>& elemPoint, const std::vector<felInt>& elemIdPoint, felInt& iel, FlagMatrixRHS flagMatrixRHS) {
281 IGNORE_UNUSED_ELEM_ID_POINT;
282 IGNORE_UNUSED_FLAG_MATRIX_RHS;
283
1/4
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✓ Branch 2 taken 8192 times.
✗ Branch 3 not taken.
 8192 switch(m_problem) {
284 ✗ case 0: { // FKPP
285 ✗ m_fePotTransMemb->updateMeasQuadPt(0, elemPoint);
286 ✗ m_fePotTransMemb->updateFirstDeriv(0, elemPoint);
287
288
289 // m_Matrix[0] = Am * grad(phi_i) * grad(phi_j)
290 ✗ m_elementMat[m_idL]->grad_phi_i_grad_phi_j(FelisceParam::instance().Am,*m_fePotTransMemb,0,0,1);
291
292 ✗ if (FelisceParam::instance().hasInitialCondition || FelisceParam::instance().schroFilter) {
293 // m_Matrix[0] += - m_coefChi * V * phi_i * phi_j
294 ✗ m_elemFieldU0.setValue(m_U_0_seq, *m_fePotTransMemb, iel, m_ipotTransMemb, m_ao, m_dof);
295 ✗ m_elementMat[m_idL]->a_phi_i_phi_j(- m_coefChi,m_elemFieldU0,*m_fePotTransMemb,0,0,1);
296 }
297
298 // Matrix[1] = phi_i * phi_j
299 ✗ m_elementMat[m_idG]->phi_i_phi_j(1.,*m_fePotTransMemb,0,0,1);
300
301 ✗ if (!FelisceParam::instance().solveEigenProblem) {
302 ✗ if (m_useImprovedRec) {
303 // m_Matrix[m_numberOfMatrix-1] = grad(phi_i) * grad(phi_j)
304 ✗ m_elementMat[m_idK]->grad_phi_i_grad_phi_j(1.,*m_fePotTransMemb,0,0,1);
305 }
306 }
307 ✗ break;
308 }
309 ✗ case 1: { // Monodomain
310 ✗ m_fePotTransMemb->updateMeasQuadPt(0, elemPoint);
311 ✗ m_fePotTransMemb->updateFirstDeriv(0, elemPoint);
312
313
314 // m_Matrix[0] = grad(phi_i) * grad(phi_j)
315 ✗ m_elementMat[m_idL]->grad_phi_i_grad_phi_j(FelisceParam::instance().intraTransvTensor,*m_fePotTransMemb,0,0,1);
316
317 ✗ if (FelisceParam::instance().testCase == 1) {
318 ✗ std::vector<double> elemFiber;
319 ✗ getFiberDirection(iel,m_ipotTransMemb, elemFiber);
320 // m_Matrix[0] += (\sigma_i^l-\sigma_i^t) a vec a grad(phi_i) * grad(phi_j)
321 ✗ m_elementMat[m_idL]->tau_grad_phi_i_tau_grad_phi_j(FelisceParam::instance().intraFiberTensor-FelisceParam::instance().intraTransvTensor,elemFiber,*m_fePotTransMemb,0,0,1);
322 }
323
324 ✗ if (FelisceParam::instance().hasInitialCondition || FelisceParam::instance().schroFilter) {
325 // m_Matrix[0] += - m_coefChi * V * phi_i * phi_j
326 ✗ m_elemFieldU0.setValue(m_U_0_seq, *m_fePotTransMemb, iel, m_ipotTransMemb, m_ao, m_dof);
327 ✗ m_elementMat[m_idL]->a_phi_i_phi_j(- m_coefChi,m_elemFieldU0,*m_fePotTransMemb,0,0,1);
328 }
329
330 // m_Matrix[1] = phi_i * phi_j
331 ✗ m_elementMat[m_idG]->phi_i_phi_j(1.,*m_fePotTransMemb,0,0,1);
332
333 ✗ if (!FelisceParam::instance().solveEigenProblem) {
334
335 ✗ if (FelisceParam::instance().hasInfarct) {
336 // m_Matrix[1] = B_ij = s(x) * phi_i * phi_j
337 ✗ m_elemFieldFhNf0.setValue(m_FhNf0, *m_fePotTransMemb, iel, m_ipotTransMemb, m_ao, m_dof);
338 ✗ m_elementMat[m_idGs]->a_phi_i_phi_j(1.,m_elemFieldFhNf0,*m_fePotTransMemb,0,0,1);
339 } else {
340 ✗ double s = FelisceParam::instance().f0;
341 // m_Matrix[1] = B_ij = s * phi_i * phi_j
342 ✗ m_elementMat[m_idGs]->phi_i_phi_j(s,*m_fePotTransMemb,0,0,1);
343 }
344
345 ✗ if ( ( (m_tensorOrder == 4 ) & (m_useImprovedRec) ) || (m_tensorOrder == 3 ) ) {
346 // m_Matrix[2] = E_ij = grad(phi_i) * grad(phi_j)
347 ✗ m_elementMat[m_idK]->grad_phi_i_grad_phi_j(FelisceParam::instance().intraTransvTensor,*m_fePotTransMemb,0,0,1);
348 ✗ if (FelisceParam::instance().testCase == 1) {
349 ✗ std::vector<double> elemFiber;
350 ✗ getFiberDirection(iel,m_ipotTransMemb, elemFiber);
351 // m_Matrix[2] += (\sigma_i^l-\sigma_i^t) a vec a grad(phi_i) * grad(phi_j)
352 ✗ m_elementMat[m_idK]->tau_grad_phi_i_tau_grad_phi_j(FelisceParam::instance().intraFiberTensor-FelisceParam::instance().intraTransvTensor,elemFiber,*m_fePotTransMemb,0,0,1);
353 }
354 }
355
356 ✗ if (m_tensorOrder == 3 ) {
357 ✗ if (FelisceParam::instance().hasInfarct) {
358 ✗ if (FelisceParam::instance().testCase == 2) {
359 // m_Matrix[3] = E^s_ij = s(x) * \sigma_i^t grad(phi_i), grad(phi_j)
360 ✗ m_elemFieldFhNf0.setValue(m_FhNf0, *m_fePotTransMemb, iel, m_ipotTransMemb, m_ao, m_dof);
361 ✗ m_elementMat[m_idKs]->a_grad_phi_i_grad_phi_j(FelisceParam::instance().intraTransvTensor,m_elemFieldFhNf0,*m_fePotTransMemb,0,0,1);
362 ✗ } else if (FelisceParam::instance().testCase == 1) {
363 ✗ FEL_ERROR("Error: infarct in fiber directed conductivity case not yet implemented.");
364 }
365 } else {
366 ✗ double s = FelisceParam::instance().f0;
367 ✗ if (FelisceParam::instance().testCase == 2) {
368 // m_Matrix[3] = E^s_ij = s*\sigma_i^t grad(phi_i), grad(phi_j)
369 ✗ m_elementMat[m_idKs]->grad_phi_i_grad_phi_j(s*FelisceParam::instance().intraTransvTensor,*m_fePotTransMemb,0,0,1);
370 ✗ } else if (FelisceParam::instance().testCase == 1) {
371 ✗ std::vector<double> elemFiber;
372 ✗ getFiberDirection(iel,m_ipotTransMemb, elemFiber);
373 // m_Matrix[3] = E^s_ij += s*(\sigma_i^l-\sigma_i^t) a vec a grad(phi_i) * grad(phi_j)
374 ✗ m_elementMat[m_idKs]->tau_grad_phi_i_tau_grad_phi_j(s*(FelisceParam::instance().intraFiberTensor-FelisceParam::instance().intraTransvTensor),elemFiber,*m_fePotTransMemb,0,0,1);
375 }
376 }
377 }
378 }
379 ✗ break;
380 }
381 8192 case 2: {
382 8192 m_fePotTransMemb->updateMeasQuadPt(0, elemPoint);
383 8192 m_fePotTransMemb->updateFirstDeriv(0, elemPoint);
384
385
386 // m_Matrix[0] = grad(phi_i) * grad(phi_j)
387 8192 m_elementMat[m_idL]->grad_phi_i_grad_phi_j(FelisceParam::instance().intraTransvTensor,*m_fePotTransMemb,0,0,1);
388
389
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 8192 times.
 8192 if (FelisceParam::instance().testCase == 1) {
390 ✗ std::vector<double> elemFiber;
391 ✗ getFiberDirection(iel,m_ipotTransMemb, elemFiber);
392 // m_Matrix[0] = \sigma_i^t \grad V_m
393 ✗ m_elementMat[m_idL]->grad_phi_i_grad_phi_j(FelisceParam::instance().intraTransvTensor,*m_fePotTransMemb,0,0,1);
394 // m_Matrix[0] += (\sigma_i^l-\sigma_i^t) a vec a grad(phi_i) * grad(phi_j)
395 ✗ m_elementMat[m_idL]->tau_grad_phi_i_tau_grad_phi_j(FelisceParam::instance().intraFiberTensor-FelisceParam::instance().intraTransvTensor,elemFiber,*m_fePotTransMemb,0,0,1);
396 }
397
398
3/6
✓ Branch 1 taken 8192 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 8192 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 8192 times.
 8192 if (FelisceParam::instance().hasInitialCondition || FelisceParam::instance().schroFilter) {
399 // m_Matrix[0] += - m_coefChi * V * phi_i * phi_j
400 ✗ m_elemFieldU0.setValue(m_U_0_seq, *m_fePotTransMemb, iel, m_ipotTransMemb, m_ao, m_dof);
401 ✗ m_elementMat[m_idL]->a_phi_i_phi_j(- m_coefChi,m_elemFieldU0,*m_fePotTransMemb,0,0,1);
402 }
403
404 // m_Matrix[1] = phi_i * phi_j
405 8192 m_elementMat[m_idG]->phi_i_phi_j(1.,*m_fePotTransMemb,0,0,1);
406
407
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 8192 times.
 8192 if (!FelisceParam::instance().solveEigenProblem) {
408 ✗ if (FelisceParam::instance().hasInfarct) {
409 // m_Matrix[1] = B_ij = s(x) * phi_i * phi_j
410 ✗ m_elemFieldFhNf0.setValue(m_FhNf0, *m_fePotTransMemb, iel, m_ipotTransMemb, m_ao, m_dof);
411 ✗ m_elementMat[m_idGs]->a_phi_i_phi_j(1.,m_elemFieldFhNf0,*m_fePotTransMemb,0,0,1);
412 } else {
413 ✗ double s = FelisceParam::instance().f0;
414 // m_Matrix[1] = B_ij = s * phi_i * phi_j
415 ✗ m_elementMat[m_idGs]->phi_i_phi_j(s,*m_fePotTransMemb,0,0,1);
416 }
417
418 // Matrix[2] = E_ij = \sigma_i^t grad(phi_i), grad(phi_j)
419 ✗ m_elementMat[m_idK]->grad_phi_i_grad_phi_j(FelisceParam::instance().intraTransvTensor,*m_fePotTransMemb,0,0,1);
420 // Matrix[3] = Q_ij = (\sigma_i^t+\sigma_e^t) grad(phi_i), grad(phi_j)
421 ✗ m_elementMat[m_idKie]->grad_phi_i_grad_phi_j(FelisceParam::instance().intraTransvTensor+FelisceParam::instance().extraTransvTensor,*m_fePotTransMemb,0,0,1);
422 ✗ if (FelisceParam::instance().testCase == 1) {
423 ✗ std::vector<double> elemFiber;
424 ✗ getFiberDirection(iel,m_ipotTransMemb, elemFiber);
425 // Matrix[2] = E_ij = (\sigma_i^l-\sigma_i^t) a vec a grad(phi_i), grad(phi_j)
426 ✗ m_elementMat[m_idK]->tau_grad_phi_i_tau_grad_phi_j(FelisceParam::instance().intraFiberTensor-FelisceParam::instance().intraTransvTensor,elemFiber,*m_fePotTransMemb,0,0,1);
427 // Matrix[3] = Q_ij += (\sigma_i^l-\sigma_i^t+\sigma_e^l-\sigma_e^t) a vec a grad(phi_i) * grad(phi_j)
428 ✗ m_elementMat[m_idKie]->tau_grad_phi_i_tau_grad_phi_j(FelisceParam::instance().intraFiberTensor-FelisceParam::instance().intraTransvTensor+FelisceParam::instance().extraFiberTensor-FelisceParam::instance().extraTransvTensor,elemFiber,*m_fePotTransMemb,0,0,1);
429 }
430
431 ✗ if (m_tensorOrder == 3 ) {
432 ✗ if (FelisceParam::instance().hasInfarct) {
433 ✗ if (FelisceParam::instance().testCase == 2) {
434 // m_Matrix[3] = E^s_ij = s(x) * \sigma_i^t grad(phi_i), grad(phi_j)
435 ✗ m_elemFieldFhNf0.setValue(m_FhNf0, *m_fePotTransMemb, iel, m_ipotTransMemb, m_ao, m_dof);
436 ✗ m_elementMat[m_idKs]->a_grad_phi_i_grad_phi_j(FelisceParam::instance().intraTransvTensor,m_elemFieldFhNf0,*m_fePotTransMemb,0,0,1);
437 ✗ } else if (FelisceParam::instance().testCase == 1) {
438 ✗ FEL_ERROR("Error: infarct in fiber directed conductivity case not yet implemented.");
439 }
440 } else {
441 ✗ double s = FelisceParam::instance().f0;
442 ✗ if (FelisceParam::instance().testCase == 2) {
443 // m_Matrix[3] = E^s_ij = s*\sigma_i^t grad(phi_i), grad(phi_j)
444 ✗ m_elementMat[m_idKs]->grad_phi_i_grad_phi_j(s*FelisceParam::instance().intraTransvTensor,*m_fePotTransMemb,0,0,1);
445 ✗ } else if (FelisceParam::instance().testCase == 1) {
446 ✗ std::vector<double> elemFiber;
447 ✗ getFiberDirection(iel,m_ipotTransMemb, elemFiber);
448 // m_Matrix[3] = E^s_ij += s*(\sigma_i^l-\sigma_i^t) a vec a grad(phi_i) * grad(phi_j)
449 ✗ m_elementMat[m_idKs]->tau_grad_phi_i_tau_grad_phi_j(s*(FelisceParam::instance().intraFiberTensor-FelisceParam::instance().intraTransvTensor),elemFiber,*m_fePotTransMemb,0,0,1);
450 }
451 }
452 }
453 }
454 8192 break;
455 }
456 ✗ default:
457 ✗ FEL_ERROR("Model not defined for ALP solver.");
458 ✗ break;
459 }
460 8192 }
461
462 // Visualize modes as a variable
463 1 void EigenProblemALP::writeMode(const int iter) {
464 int rankProc;
465
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 MPI_Comm_rank(m_petscComm,&rankProc);
466
467 1 std::string fileName;
468
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 1 if (iter > 0)
469 ✗ fileName = "basis" + std::to_string(iter);
470 else
471
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 fileName = "basis";
472
473 1 double* tmpSolToSave = nullptr;
474
2/4
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 1 tmpSolToSave = new double[m_numDof];
475
476
2/2
✓ Branch 0 taken 25 times.
✓ Branch 1 taken 1 times.
 26 for (int iBasis=0; iBasis<m_dimRomBasis; iBasis++) {
477
1/2
✓ Branch 2 taken 25 times.
✗ Branch 3 not taken.
 25 fromVecToDoubleStar(tmpSolToSave, m_basis[iBasis], rankProc, 1);
478
2/4
✓ Branch 1 taken 25 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 25 times.
✗ Branch 5 not taken.
 25 writeEnsightVector(tmpSolToSave, iBasis, fileName);
479 }
480
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 delete [] tmpSolToSave;
481
482
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 if (rankProc == 0)
483
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 1 writeEnsightCase(m_dimRomBasis, 1.,fileName);
484
485 1 }
486
487 ✗ void EigenProblemALP::writeEnsightSolution(const int iter) {
488 ✗ PetscVector solUe;
489 ✗ solUe.duplicateFrom(m_U_0);
490 ✗ PetscVector solW;
491 ✗ solW.duplicateFrom(m_W_0);
492
493 ✗ if (iter > 0) {
494 ✗ projectOnDof(m_beta,m_sol,m_dimRomBasis);
495 ✗ if ( (m_problem == 1) || (m_problem == 2)) {
496 ✗ if (FelisceParam::instance().printIonicVar) {
497 ✗ projectOnDof(m_mu,solW,m_dimRomBasis);
498 }
499 }
500 ✗ if (m_problem == 2) {
501 ✗ projectOnDof(m_xi,solUe,m_dimRomBasis);
502 }
503 } else {
504 ✗ m_sol.copyFrom(m_U_0);
505 ✗ if ( (m_problem == 1) || (m_problem == 2)) {
506 ✗ if (FelisceParam::instance().printIonicVar) {
507 ✗ solW.copyFrom(m_W_0);
508 }
509 }
510 ✗ if (m_problem == 2) {
511 ✗ solUe.copyFrom(m_Ue_0);
512 }
513 }
514
515 int rankProc;
516 ✗ MPI_Comm_rank(m_petscComm,&rankProc);
517
518 ✗ double* tmpSolToSave = nullptr;
519 ✗ tmpSolToSave = new double[m_numDof];
520
521 ✗ if (m_problem == 0) {
522 ✗ fromVecToDoubleStar(tmpSolToSave, m_sol, rankProc, 1);
523 ✗ writeEnsightVector(tmpSolToSave, iter, "sol");
524 }
525
526 ✗ if ( (m_problem == 1) || (m_problem == 2)) {
527 ✗ fromVecToDoubleStar(tmpSolToSave, m_sol, rankProc, 1);
528 ✗ writeEnsightVector(tmpSolToSave, iter, "potTransMemb");
529 ✗ if (FelisceParam::instance().printIonicVar) {
530 ✗ fromVecToDoubleStar(tmpSolToSave, solW, rankProc, 1);
531 ✗ writeEnsightVector(tmpSolToSave, iter, "ionicVar");
532 }
533 }
534
535 ✗ if (m_problem == 2) {
536 ✗ fromVecToDoubleStar(tmpSolToSave, solUe, rankProc, 1);
537 ✗ writeEnsightVector(tmpSolToSave, iter, "potExtraCell");
538 }
539
540 ✗ delete [] tmpSolToSave;
541 ✗ solUe.destroy();
542 ✗ solW.destroy();
543
544 ✗ std::vector<std::string> varNames;
545 ✗ if (m_problem == 0) {
546 ✗ varNames.resize(1);
547 ✗ varNames[0] = "sol";
548 ✗ } else if (m_problem == 1) {
549 ✗ if (FelisceParam::instance().printIonicVar) {
550 ✗ varNames.resize(2);
551 ✗ varNames[0] = "potTransMemb";
552 ✗ varNames[1] = "ionicVar";
553 } else {
554 ✗ varNames.resize(1);
555 ✗ varNames[0] = "potTransMemb";
556 }
557 ✗ } else if (m_problem == 2) {
558 ✗ if (FelisceParam::instance().printIonicVar) {
559 ✗ varNames.resize(3);
560 ✗ varNames[0] = "potTransMemb";
561 ✗ varNames[1] = "potExtraCell";
562 ✗ varNames[2] = "ionicVar";
563 } else {
564 ✗ varNames.resize(2);
565 ✗ varNames[0] = "potTransMemb";
566 ✗ varNames[1] = "potExtraCell";
567 }
568 }
569 ✗ if (rankProc == 0)
570 ✗ writeEnsightCase(iter+1, FelisceParam::instance().timeStep * FelisceParam::instance().frequencyWriteSolution,varNames, FelisceParam::instance().time);
571
572 }
573
574 // Solve eigenproblem to build rom basis
575 1 void EigenProblemALP::buildSolver() {
576
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 if (FelisceParam::verbose()) PetscPrintf(PETSC_COMM_WORLD,"\nBuilding Slepc Solver.\n");
577
578 #ifdef FELISCE_WITH_SLEPC
579 // m_Matrix[0] * Phi = lambda * m_Matrix[1] * Phi
580 1 EPSSetOperators(m_eps,m_Matrix[m_idL].toPetsc(),m_Matrix[m_idG].toPetsc());
581 1 EPSSetDimensions(m_eps,2*m_dimRomBasis,PETSC_DECIDE,PETSC_DECIDE);
582 1 EPSSetProblemType(m_eps,EPS_GHEP);
583 1 EPSSetType(m_eps, EPSKRYLOVSCHUR);
584 1 EPSSetTolerances(m_eps,1.0e-8,1000); // 1.e-05, 200
585 1 EPSKrylovSchurSetRestart(m_eps,0.7);
586 1 EPSSetFromOptions(m_eps);
587 #endif
588 1 }
589
590 1 void EigenProblemALP::solve() {
591
3/6
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
 1 if (FelisceParam::verbose()) PetscPrintf(PETSC_COMM_WORLD,"\nSolve eigenvalue problem.\n");
592 #ifdef FELISCE_WITH_SLEPC
593 PetscInt nconv;
594 //First request a singular value from one end of the spectrum
595
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 EPSSetWhichEigenpairs(m_eps,EPS_SMALLEST_REAL);
596 //EPSSetInterval(m_eps,-1.e-1,200.);
597
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 1 if (m_verbose >= 10) {
598 ✗ m_Matrix[m_idL].view();
599 ✗ m_Matrix[m_idG].view();
600 ✗ EPSView(m_eps,PETSC_VIEWER_STDOUT_WORLD);
601 }
602
603
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 EPSSolve(m_eps);
604
605 //Get number of converged singular values
606
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 EPSGetConverged(m_eps,&nconv);
607
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 PetscPrintf(PETSC_COMM_WORLD," Number of converged values = %d.\n", nconv);
608
609
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 1 if ( m_dimRomBasis > m_numDof ) {
610 ✗ PetscPrintf(PETSC_COMM_WORLD," Warning! Number of function basis bigger than number of dof.\n");
611 ✗ m_dimRomBasis = m_numDof;
612 }
613
614
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 1 if ( m_dimRomBasis > nconv ) {
615 ✗ PetscPrintf(PETSC_COMM_WORLD," Warning! Number of function basis bigger than number of converged values.\n");
616 ✗ m_dimRomBasis = nconv;
617 }
618
619 double sigma;
620 //Get converged singular values: largest singular value is stored in sigma.
621
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 if (nconv > 0) {
622
2/2
✓ Branch 0 taken 25 times.
✓ Branch 1 taken 1 times.
 26 for (int i=0; i < m_dimRomBasis; i++) {
623
1/2
✓ Branch 1 taken 25 times.
✗ Branch 2 not taken.
 25 EPSGetEigenpair(m_eps,i,&sigma,FELISCE_PETSC_NULLPTR,FELISCE_PETSC_NULLPTR,FELISCE_PETSC_NULLPTR);
624 }
625 } else {
626 ✗ m_dimRomBasis = 0;
627 ✗ PetscPrintf(PETSC_COMM_WORLD," Unable to compute smallest singular value!\n\n");
628 ✗ exit(1);
629 }
630
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 PetscPrintf(PETSC_COMM_WORLD, "m_dimRomBasis = %d \n", m_dimRomBasis);
631
632
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 if (m_eigenValue == nullptr)
633
2/4
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 1 m_eigenValue = new double[m_dimRomBasis];
634
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 m_basis.resize(m_dimRomBasis);
635 double norm2;
636 double normM;
637
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 PetscVector tmpVec;
638
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 m_Matrix[m_idG].getVecs(nullPetscVector,tmpVec);
639 1 double imgpart = 0.;
640
2/2
✓ Branch 0 taken 25 times.
✓ Branch 1 taken 1 times.
 26 for (int i=0; i < m_dimRomBasis; i++) {
641
1/2
✓ Branch 3 taken 25 times.
✗ Branch 4 not taken.
 25 m_Matrix[m_idG].getVecs(nullPetscVector,m_basis[i]);
642
2/4
✓ Branch 2 taken 25 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 25 times.
✗ Branch 6 not taken.
 25 EPSGetEigenpair(m_eps,i,&m_eigenValue[i],&imgpart,m_basis[i].toPetsc(),FELISCE_PETSC_NULLPTR);
643
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 25 times.
 25 if (Tools::notEqual(imgpart,0.0)) {
644 ✗ PetscPrintf(PETSC_COMM_WORLD, "Warning! Im(m_eigenValue[%d]) = %e \n", i, imgpart);
645 }
646
1/2
✓ Branch 2 taken 25 times.
✗ Branch 3 not taken.
 25 m_basis[i].norm(NORM_2,&norm2);
647
648
1/2
✓ Branch 3 taken 25 times.
✗ Branch 4 not taken.
 25 mult(m_Matrix[m_idG], m_basis[i], tmpVec);
649
1/2
✓ Branch 2 taken 25 times.
✗ Branch 3 not taken.
 25 dot(m_basis[i],tmpVec,&normM);
650
1/2
✓ Branch 0 taken 25 times.
✗ Branch 1 not taken.
 25 if (m_verbose >= 1) {
651
1/2
✓ Branch 1 taken 25 times.
✗ Branch 2 not taken.
 25 PetscPrintf(PETSC_COMM_WORLD, "m_eigenValue[%d] = %e \n", i, m_eigenValue[i]);
652
1/2
✓ Branch 1 taken 25 times.
✗ Branch 2 not taken.
 25 PetscPrintf(PETSC_COMM_WORLD, "norm2[%d] = %e \n", i, norm2);
653
1/2
✓ Branch 1 taken 25 times.
✗ Branch 2 not taken.
 25 PetscPrintf(PETSC_COMM_WORLD, "normM[%d] = %e \n", i, normM);
654 }
655 }
656
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 tmpVec.destroy();
657
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 1 std::string fileName = FelisceParam::instance().resultDir + "/eigenValue";
658
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 1 viewALP(m_eigenValue,m_dimRomBasis,fileName);
659 #endif
660 1 }
661
662 // Reduced model functions
663 //========================
664 // Initialization
665 1 void EigenProblemALP::readData(IO& io) {
666
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 if (FelisceParam::verbose() > 0) PetscPrintf(PETSC_COMM_WORLD,"\nEigenProblemALP::readData.\n");
667 1 m_Matrix[m_idG].getVecs(m_U_0,nullPetscVector);
668 1 m_U_0.setFromOptions();
669
670 1 m_U_0_seq.createSeq(PETSC_COMM_SELF,m_numDof);
671 1 m_U_0_seq.setFromOptions();
672
673 1 m_W_0.duplicateFrom(m_U_0);
674 1 m_W_0.copyFrom(m_U_0);
675
676 1 m_W_0_seq.createSeq(PETSC_COMM_SELF,m_numDof);
677 1 m_W_0_seq.setFromOptions();
678
679 1 m_Ue_0.duplicateFrom(m_U_0);
680 1 m_Ue_0.copyFrom(m_U_0);
681
682 1 m_Ue_0_seq.createSeq(PETSC_COMM_SELF,m_numDof);
683 1 m_Ue_0_seq.setFromOptions();
684
685
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
 1 if (FelisceParam::instance().hasInitialCondition) {
686 // Read initial solution file (*.00000.scl)
687 ✗ int idInVar = 0;
688 // potTransMemb
689 ✗ double* initialSolution = nullptr;
690 ✗ initialSolution = new double[m_mesh->numPoints()];
691 ✗ io.readVariable(idInVar, 0.,initialSolution, m_mesh->numPoints());
692 ✗ idInVar ++;
693 // ionicVar
694 ✗ double* ionicSolution = nullptr;
695 ✗ if ( (m_problem == 1) || (m_problem == 2) ) {
696 ✗ ionicSolution = new double[m_mesh->numPoints()];
697 ✗ io.readVariable(idInVar, 0.,ionicSolution, m_mesh->numPoints());
698 ✗ idInVar ++;
699 }
700 // potExtraCell
701 ✗ double* extraCellSolution = nullptr;
702 ✗ if (m_problem == 2) {
703 ✗ extraCellSolution = new double[m_mesh->numPoints()];
704 ✗ io.readVariable(idInVar, 0.,extraCellSolution, m_mesh->numPoints());
705 ✗ idInVar ++;
706 }
707 // initialize (parallel) initial solution vectors
708 ✗ fromDoubleStarToVec(initialSolution, m_U_0, m_mesh->numPoints());
709 ✗ if ( (m_problem == 1) || (m_problem == 2) )
710 ✗ fromDoubleStarToVec(ionicSolution, m_W_0, m_mesh->numPoints());
711 ✗ if (m_problem == 2 )
712 ✗ fromDoubleStarToVec(extraCellSolution, m_Ue_0, m_mesh->numPoints());
713
714 // initialize SEQUENTIAL initial solution vectors (in order to use elemField and define collocation pts solution)
715 felInt ia;
716 ✗ for (felInt i=0; i< m_numDof; i++) {
717 ✗ ia = i;
718 ✗ AOApplicationToPetsc(m_ao,1,&ia);
719 ✗ m_U_0_seq.setValues(1,&ia,&initialSolution[i],INSERT_VALUES);
720 ✗ if ( (m_problem == 1) || (m_problem == 2) )
721 ✗ m_W_0_seq.setValues(1,&ia,&ionicSolution[i],INSERT_VALUES);
722 ✗ if (m_problem == 2)
723 ✗ m_Ue_0_seq.setValues(1,&ia,&extraCellSolution[i],INSERT_VALUES);
724 }
725 ✗ m_U_0_seq.assembly();
726 ✗ if ( (m_problem == 1) || (m_problem == 2) ) {
727 ✗ m_W_0_seq.assembly();
728 }
729 ✗ if (m_problem == 2) {
730 ✗ m_Ue_0_seq.assembly();
731 }
732
733 ✗ if (initialSolution) delete [] initialSolution;
734 ✗ if (ionicSolution) delete [] ionicSolution;
735 ✗ if (extraCellSolution) delete [] extraCellSolution;
736
737 // Read fibers file (*.00000.vct)
738 ✗ if (FelisceParam::instance().testCase == 1) {
739 ✗ if (m_fiber == nullptr)
740 ✗ m_fiber = new double[m_mesh->numPoints()*3];
741 ✗ io.readVariable(idInVar, 0.,m_fiber, m_mesh->numPoints()*3);
742 ✗ idInVar ++;
743 }
744 }
745 else {
746 1 PetscPrintf(PETSC_COMM_WORLD, "U_0 not read from file (zero vector).");
747
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 m_U_0.set( 0.);
748 1 m_U_0.assembly();
749
2/4
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 if ( (m_problem == 1) || (m_problem == 2) ) {
750
1/2
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
 1 if (FelisceParam::instance().typeOfIonicModel == "schaf") {
751 ✗ m_W_0.set( 1.);
752 ✗ m_W_0_seq.set( 1.);
753 }
754 else {
755
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 m_W_0.set( 0.);
756
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 m_W_0_seq.set( 0.);
757 }
758 1 m_W_0.assembly();
759 1 m_W_0_seq.assembly();
760 }
761
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 m_U_0_seq.set( 0.);
762 1 m_U_0_seq.assembly();
763
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 if (m_problem == 2) {
764
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 m_Ue_0.set( 0.);
765 1 m_Ue_0.assembly();
766
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 m_Ue_0_seq.set( 0.);
767 1 m_Ue_0_seq.assembly();
768 }
769 }
770
771 1 }
772
773 ✗ void EigenProblemALP::getFiberDirection(felInt iel, int iUnknown, std::vector<double>& elemFiber) {
774 ✗ int numSupport = static_cast<int>(m_supportDofUnknown[iUnknown].iEle()[iel+1] - m_supportDofUnknown[iUnknown].iEle()[iel]);
775 ✗ elemFiber.resize( numSupport*3,0.);
776 ✗ for (int i = 0; i < numSupport; i++) {
777 ✗ elemFiber[i*3] = m_fiber[3*(m_supportDofUnknown[iUnknown].iSupportDof()[m_supportDofUnknown[iUnknown].iEle()[iel]+i])];
778 ✗ elemFiber[i*3+1] = m_fiber[3*(m_supportDofUnknown[iUnknown].iSupportDof()[m_supportDofUnknown[iUnknown].iEle()[iel]+i])+1];
779 ✗ elemFiber[i*3+2] = m_fiber[3*(m_supportDofUnknown[iUnknown].iSupportDof()[m_supportDofUnknown[iUnknown].iEle()[iel]+i])+2];
780 }
781 }
782
783 ✗ void EigenProblemALP::readBasis() {
784 ✗ if (FelisceParam::verbose() > 0 ) PetscPrintf(PETSC_COMM_WORLD,"\nEigenProblemALP::readBasis()\n");
785
786 ✗ m_basis.resize(m_dimRomBasis);
787 ✗ for (int i=0; i<m_dimRomBasis; i++) {
788 ✗ m_Matrix[m_idG].getVecs(m_basis[i],nullPetscVector);
789 ✗ m_basis[i].setFromOptions();
790 ✗ double* vec = new double[m_numDof];
791 ✗ readEnsightFile(vec,"basis",i,m_numDof);
792 ✗ fromDoubleStarToVec(vec, m_basis[i], m_numDof);
793 ✗ delete [] vec;
794 ✗ if (FelisceParam::verbose() > 40) m_basis[i].view();
795 }
796
797 ✗ readEigenValueFromFile();
798
799 }
800
801 ✗ void EigenProblemALP::initializeROM() {
802 ✗ if (FelisceParam::verbose() > 10 ) PetscPrintf(PETSC_COMM_WORLD,"\nEigenProblemALP::initializeROM()\n");
803
804 ✗ m_basis.resize(m_dimRomBasis);
805 ✗ for (int i=0; i<m_dimRomBasis; i++) {
806 ✗ m_Matrix[m_idG].getVecs(m_basis[i],nullPetscVector);
807 ✗ m_basis[i].setFromOptions();
808 ✗ double* vec = new double[m_numDof];
809 ✗ readEnsightFile(vec,"basis",i,m_numDof);
810 ✗ fromDoubleStarToVec(vec, m_basis[i], m_numDof);
811 ✗ delete [] vec;
812 ✗ if (FelisceParam::verbose() > 40) m_basis[i].view();
813 }
814
815 ✗ if(m_useImprovedRec==true) {
816 ✗ m_orthComp.resize(m_dimOrthComp);
817 ✗ for(int i=0; i<m_dimOrthComp; i++) {
818 ✗ int fileInd = i + m_dimRomBasis;
819 ✗ m_Matrix[m_idG].getVecs(m_orthComp[i],nullPetscVector);
820 ✗ m_orthComp[i].setFromOptions();
821 ✗ double* vec = new double[m_numDof];
822 ✗ readEnsightFile(vec,"basis",fileInd,m_numDof);
823 ✗ fromDoubleStarToVec(vec, m_orthComp[i], m_numDof);
824 ✗ delete [] vec;
825 ✗ if (FelisceParam::verbose() > 40) m_orthComp[i].view();
826 }
827 }
828
829 ✗ readEigenValueFromFile();
830 ✗ initializeSolution();
831 }
832
833 ✗ void EigenProblemALP::readEnsightFile(double* vec, std::string varName, int idIter, felInt size) {
834 ✗ std::string iteration;
835 ✗ if (idIter < 10)
836 ✗ iteration = "0000" + std::to_string(idIter);
837 ✗ else if (idIter < 100)
838 ✗ iteration = "000" + std::to_string(idIter);
839 ✗ else if (idIter < 1000)
840 ✗ iteration = "00" + std::to_string(idIter);
841 ✗ else if (idIter < 10000)
842 ✗ iteration = "0" + std::to_string(idIter);
843 ✗ else if (idIter < 100000)
844 ✗ iteration = std::to_string(idIter);
845
846 ✗ std::string fileName = FelisceParam::instance().inputDirectory + "/" + varName + "." + iteration + ".scl";
847 ✗ if (FelisceParam::verbose() > 1) PetscPrintf(PETSC_COMM_WORLD, "Read file %s\n", fileName.c_str());
848 ✗ FILE* pFile = fopen (fileName.c_str(),"r");
849 ✗ if (pFile==nullptr) {
850 ✗ FEL_ERROR("Impossible to read "+fileName +".");
851 }
852 char str[80];
853 ✗ std::string display;
854 ✗ if ( fscanf(pFile,"%s", str) !=1 ) {
855 ✗ FEL_WARNING("Failed ot read str");
856 }
857 ✗ display = str;
858 ✗ if ( fscanf(pFile,"%s", str) !=1 ) {
859 ✗ FEL_WARNING("Failed ot read str");
860 }
861 ✗ display = str;
862 ✗ if ( fscanf(pFile,"%s", str) !=1) {
863 ✗ FEL_WARNING("Failed ot read str");
864 }
865 ✗ display = str;
866 //Read values
867 ✗ int count = 0;
868 ✗ for ( felInt i = 0; i < size; i++) {
869 ✗ if ( fscanf(pFile,"%lf", &vec[i]) !=1 ) {
870 ✗ FEL_WARNING("Failed ot read vec");
871 }
872 ✗ if ( count == 6 ) {
873 ✗ if ( fscanf( pFile, "\n" ) !=0 ) {
874 ✗ FEL_WARNING("Failed ot read");
875 }
876 ✗ count=0;
877 }
878 ✗ count++;
879 }
880 ✗ fclose(pFile);
881 }
882
883 ✗ void EigenProblemALP::readEigenValueFromFile() {
884 ✗ std::string fileName = FelisceParam::instance().inputDirectory + "/eigenValue";
885 ✗ if (FelisceParam::verbose() > 1) PetscPrintf(PETSC_COMM_WORLD, "Read file %s\n", fileName.c_str());
886 ✗ FILE* pFile = fopen (fileName.c_str(),"r");
887 ✗ if (pFile==nullptr) {
888 ✗ FEL_ERROR("Impossible to read "+fileName +".");
889 }
890 //Read values
891 ✗ m_eigenValue = new double[m_dimRomBasis];
892 ✗ for ( felInt i = 0; i < m_dimRomBasis; i++) {
893 ✗ if ( fscanf(pFile,"%lf", &m_eigenValue[i]) !=1 ) {
894 ✗ FEL_WARNING("Failed ot read vec");
895 }
896 ✗ if ( fscanf( pFile, "\n" ) !=0 ) {
897 ✗ FEL_WARNING("Failed ot read");
898 }
899 }
900 ✗ fclose(pFile);
901
902 ✗ if (FelisceParam::verbose() > 10) {
903 ✗ PetscPrintf(PETSC_COMM_WORLD, "EigenValues:\n");
904 ✗ for ( felInt i = 0; i < m_dimRomBasis; i++) {
905 ✗ PetscPrintf(PETSC_COMM_WORLD, "[%d] %e \n", i, m_eigenValue[i]);
906 }
907 }
908
909 }
910
911 ✗ void EigenProblemALP::computeGamma() {
912 ✗ if (FelisceParam::verbose() > 10 ) PetscPrintf(PETSC_COMM_WORLD,"\nEigenProblemALP::computeGamma()\n");
913
914 ✗ if (m_gamma == nullptr)
915 ✗ m_gamma = new double[m_dimRomBasis];
916
917 ✗ switch (m_problem) {
918 ✗ case 0: {
919 // FKPP:
920 ✗ double muFKPP = 10.;
921 ✗ double coef = - (m_coefChi + muFKPP);
922 // FKPP:
923 // m_gamma_p = (muFKPP - lambda_p) * beta_p - (m_coefChi + muFKPP) * sum_{ij} beta_i * beta_j * T_{ijp}
924 ✗ for (int p=0; p<m_dimRomBasis; p++) {
925
926 ✗ m_gamma[p] = (muFKPP - m_eigenValue[p]) * m_beta[p];
927 ✗ for (int i=0; i<m_dimRomBasis; i++) {
928 ✗ for (int j=0; j<m_dimRomBasis; j++) {
929 ✗ m_gamma[p] += coef * m_beta[i] * m_beta[j] * m_tensorT[thirdOrderGlobalIndex(i,j,p)];
930 }
931 }
932 }
933 ✗ break;
934 }
935 ✗ case 1:
936 case 2: {
937 // MONODOMAIN or BIDOMAIN:
938 ✗ double Am = FelisceParam::instance().Am;
939 ✗ double Cm = FelisceParam::instance().Cm;
940 ✗ double a = FelisceParam::instance().alpha;
941 ✗ double eps = FelisceParam::instance().epsilon;
942 ✗ double gam = FelisceParam::instance().gammaEl;
943
944 ✗ if (m_eta == nullptr)
945 ✗ m_eta = new double[m_dimRomBasis];
946
947 ✗ switch (m_tensorOrder) {
948 ✗ case 3: {
949 // MONODOMAIN:
950 // gamma_p = - Am mu_p - A_m a sum_{i} (a + (a+1)lambda_i / chi + lambda_i^2/chi^2 )beta_i B_{ip} - sum_{i} beta_i E_{ip} + Am sum_{i} ((a+1)/chi + lambda_i/chi^2) beta_i E^s_ip - Am/chi sum_{ij} beta_i beta_j U_ijp
951 // eta_p = eps * ro * beta_p - eps * mu_p
952
953 // BIDOMAIN:
954 // gamma_p = - Am mu_p - A_m a sum_{i} (a + (a+1)lambda_i / chi + lambda_i^2/chi^2 )beta_i B_{ip} - sum_{i} beta_i E_{ip} - sum_{i} xi_i E_{ip} + Am sum_{i} ((a+1)/chi + lambda_i/chi^2) beta_i E^s_ip - Am/chi sum_{ij} beta_i beta_j U_ijp
955 // eta_p = eps * ro * beta_p - eps * mu_p
956
957 ✗ for (int p=0; p<m_dimRomBasis; p++) {
958 ✗ m_gamma[p] = - Am * m_mu[p];
959 ✗ for (int i=0; i<m_dimRomBasis; i++) {
960 ✗ m_gamma[p] += - Am*( a + (1+a)*m_eigenValue[i]/m_coefChi + (m_eigenValue[i]*m_eigenValue[i])/(m_coefChi*m_coefChi) ) * m_beta[i] * m_tensorB[secondOrderGlobalIndex(i,p)];
961 ✗ m_gamma[p] += - m_beta[i]*m_tensorE[secondOrderGlobalIndex(i,p)];
962 ✗ if (m_problem == 2) { // bidomain
963 ✗ m_gamma[p] += - m_xi[i]*m_tensorE[secondOrderGlobalIndex(i,p)];
964 }
965 ✗ m_gamma[p] += Am*( (1+a)/m_coefChi + m_eigenValue[i]/(m_coefChi*m_coefChi) ) * m_beta[i]*m_tensorEs[secondOrderGlobalIndex(i,p)];
966 ✗ for (int j=0; j<m_dimRomBasis; j++) {
967 ✗ m_gamma[p] += - Am/m_coefChi * m_beta[i] * m_beta[j] * m_tensorU[i][secondOrderGlobalIndex(j,p)];
968 }
969 }
970 ✗ if (FelisceParam::instance().hasSource) {
971 ✗ double delay[m_numSource];
972 ✗ double tPeriod=fmod(m_fstransient->time,FelisceParam::instance().timePeriod);
973 ✗ for (int iS=0; iS<m_numSource; iS++) {
974 ✗ delay[iS] = FelisceParam::instance().delayStim[iS];
975 ✗ if ((tPeriod>=delay[iS])&&(tPeriod<=(delay[iS]+FelisceParam::instance().stimTime)))
976 ✗ m_gamma[p] += Am*m_source[iS][p];
977 }
978 }
979
980 ✗ m_gamma[p] = m_gamma[p] / (Am*Cm);
981
982 ✗ m_eta[p] = eps * ( gam * m_beta[p] - m_mu[p] );
983
984 }
985
986 ✗ break;
987 }
988 ✗ case 4: {
989 // MONODOMAIN:
990 // gamma_p = - lambda_p beta_p - Am mu_p - A_m a sum_{i} beta_i B_{ip} - \chi sum_{ij} beta_i beta_j T_{ijp} + (a+1)Am sum_{ij} beta_i beta_j T^s_ijp - Am sum_{ijk} beta_i beta_j beta_k Y_ijkp
991 // eta_p = eps * ro * beta_p - eps * mu_p
992
993 // BIDOMAIN:
994 // gamma_p = - lambda_p beta_p - lambda_p xi_p- Am mu_p - A_m a sum_{i} beta_i B_{ip} - \chi sum_{ij} beta_i beta_j T_{ijp} - \chi sum_{ij} beta_i xi_j T_{ijp} + (a+1)Am sum_{ij} beta_i beta_j T^s_ijp - Am sum_{ijk} beta_i beta_j beta_k Y_ijkp
995 // eta_p = eps * ro * beta_p - eps * mu_p
996
997 ✗ for (int p=0; p<m_dimRomBasis; p++) {
998 ✗ m_gamma[p] = - m_eigenValue[p] * m_beta[p];
999 ✗ m_gamma[p] += - Am * m_mu[p];
1000 ✗ if (m_problem == 2) { // bidomain
1001 ✗ m_gamma[p] += - m_eigenValue[p] * m_xi[p];
1002 }
1003 ✗ for (int i=0; i<m_dimRomBasis; i++) {
1004 ✗ m_gamma[p] += - a*Am * m_beta[i] * m_tensorB[secondOrderGlobalIndex(i,p)];
1005 ✗ for (int j=0; j<m_dimRomBasis; j++) {
1006 ✗ m_gamma[p] += (a+1)*Am * m_beta[i] * m_beta[j] * m_tensorTs[thirdOrderGlobalIndex(i,j,p)];
1007 ✗ m_gamma[p] += - m_coefChi * m_beta[i] * m_beta[j] * m_tensorT[thirdOrderGlobalIndex(i,j,p)];
1008 ✗ if (m_problem == 2) { // bidomain
1009 ✗ m_gamma[p] += - m_coefChi * m_beta[i] * m_xi[j] * m_tensorT[thirdOrderGlobalIndex(i,j,p)];
1010 }
1011 ✗ for (int k=0; k<m_dimRomBasis; k++) {
1012 ✗ m_gamma[p] += - Am * m_beta[i] * m_beta[j] * m_beta[k] * m_tensorY[fourthOrderGlobalIndex(i,j,k,p)];
1013 }
1014 }
1015 }
1016
1017 ✗ if (FelisceParam::instance().hasSource) {
1018 ✗ double delay[m_numSource];
1019 ✗ double tPeriod=fmod(m_fstransient->time,FelisceParam::instance().timePeriod);
1020 ✗ for (int iS=0; iS<m_numSource; iS++) {
1021 ✗ delay[iS] = FelisceParam::instance().delayStim[iS];
1022 ✗ if ((tPeriod>=delay[iS])&&(tPeriod<=(delay[iS]+FelisceParam::instance().stimTime)))
1023 ✗ m_gamma[p] += Am*m_source[iS][p];
1024 }
1025 }
1026
1027 ✗ m_gamma[p] = m_gamma[p] / (Am*Cm);
1028
1029 ✗ m_eta[p] = eps * ( gam * m_beta[p] - m_mu[p] );
1030
1031 }
1032 ✗ break;
1033 }
1034 ✗ default:
1035 ✗ break;
1036 }
1037 ✗ break;
1038 }
1039 ✗ default:
1040 ✗ FEL_ERROR("Model not defined for ALP solver.");
1041 ✗ break;
1042 }
1043 }
1044
1045 ✗ void EigenProblemALP::computeGammaExtrap() {
1046 ✗ if (FelisceParam::verbose() > 10 ) PetscPrintf(PETSC_COMM_WORLD,"\nEigenProblemALP::computeGamma()\n");
1047
1048 ✗ if (m_problem != 2) {
1049 ✗ FEL_ERROR("Time derivative method not implemented for this problem.");
1050 }
1051
1052 ✗ if (m_gamma_extrap == nullptr)
1053 ✗ m_gamma_extrap = new double[m_dimRomBasis];
1054 ✗ if (m_eta_extrap == nullptr)
1055 ✗ m_eta_extrap = new double[m_dimRomBasis];
1056
1057 ✗ double* beta_extrap = new double[m_dimRomBasis];
1058 ✗ double* mu_extrap = new double[m_dimRomBasis];
1059 ✗ double* xi_extrap = new double[m_dimRomBasis];
1060
1061 double coefExtrap[2];
1062 ✗ coefExtrap[0] = 2.;
1063 ✗ coefExtrap[1] = -1.;
1064
1065 ✗ for (int i=0; i<m_dimRomBasis; i++) {
1066 ✗ beta_extrap[i] = coefExtrap[0]*m_beta[i] + coefExtrap[1]*m_beta_n_1[i];
1067 ✗ mu_extrap[i] = coefExtrap[0]*m_mu[i] + coefExtrap[1]*m_mu_n_1[i];
1068 ✗ xi_extrap[i] = coefExtrap[0]*m_xi[i] + coefExtrap[1]*m_xi_n_1[i];
1069 }
1070
1071 // MONODOMAIN or BIDOMAIN:
1072 ✗ double Am = FelisceParam::instance().Am;
1073 ✗ double Cm = FelisceParam::instance().Cm;
1074 ✗ double a = FelisceParam::instance().alpha;
1075 ✗ double eps = FelisceParam::instance().epsilon;
1076 ✗ double gam = FelisceParam::instance().gammaEl;
1077
1078
1079 ✗ switch (m_tensorOrder) {
1080 ✗ case 3: {
1081 // MONODOMAIN:
1082 // gamma_p = - Am mu_p - A_m a sum_{i} (a + (a+1)lambda_i / chi + lambda_i^2/chi^2 )beta_i B_{ip} - sum_{i} beta_i E_{ip} + Am sum_{i} ((a+1)/chi + lambda_i/chi^2) beta_i E^s_ip - Am/chi sum_{ij} beta_i beta_j U_ijp
1083 // eta_p = eps * ro * beta_p - eps * mu_p
1084
1085 // BIDOMAIN:
1086 // gamma_p = - Am mu_p - A_m a sum_{i} (a + (a+1)lambda_i / chi + lambda_i^2/chi^2 )beta_i B_{ip} - sum_{i} beta_i E_{ip} - sum_{i} xi_i E_{ip} + Am sum_{i} ((a+1)/chi + lambda_i/chi^2) beta_i E^s_ip - Am/chi sum_{ij} beta_i beta_j U_ijp
1087 // eta_p = eps * ro * beta_p - eps * mu_p
1088
1089 ✗ for (int p=0; p<m_dimRomBasis; p++) {
1090 ✗ m_gamma_extrap[p] = - Am * mu_extrap[p];
1091 ✗ for (int i=0; i<m_dimRomBasis; i++) {
1092 ✗ m_gamma_extrap[p] += - Am*( a + (1+a)*m_eigenValue[i]/m_coefChi + (m_eigenValue[i]*m_eigenValue[i])/(m_coefChi*m_coefChi) ) * beta_extrap[i] * m_tensorB[secondOrderGlobalIndex(i,p)];
1093 ✗ m_gamma_extrap[p] += - beta_extrap[i]*m_tensorE[secondOrderGlobalIndex(i,p)];
1094 ✗ m_gamma_extrap[p] += - xi_extrap[i]*m_tensorE[secondOrderGlobalIndex(i,p)];
1095 ✗ m_gamma_extrap[p] += Am*( (1+a)/m_coefChi + m_eigenValue[i]/(m_coefChi*m_coefChi) ) * beta_extrap[i]*m_tensorEs[secondOrderGlobalIndex(i,p)];
1096 ✗ for (int j=0; j<m_dimRomBasis; j++) {
1097 ✗ m_gamma_extrap[p] += - Am/m_coefChi * beta_extrap[i] * beta_extrap[j] * m_tensorU[i][secondOrderGlobalIndex(j,p)];
1098 }
1099 }
1100 ✗ if (FelisceParam::instance().hasSource) {
1101 ✗ double delay[m_numSource];
1102 ✗ double tPeriod=fmod(m_fstransient->time,FelisceParam::instance().timePeriod);
1103 ✗ for (int iS=0; iS<m_numSource; iS++) {
1104 ✗ delay[iS] = FelisceParam::instance().delayStim[iS];
1105 ✗ if ((tPeriod>=delay[iS])&&(tPeriod<=(delay[iS]+FelisceParam::instance().stimTime)))
1106 ✗ m_gamma_extrap[p] += Am*m_source[iS][p];
1107 }
1108 }
1109
1110 ✗ m_gamma_extrap[p] = m_gamma_extrap[p] / (Am*Cm);
1111
1112 ✗ m_eta_extrap[p] = eps * ( gam * beta_extrap[p] - mu_extrap[p] );
1113
1114 }
1115
1116 ✗ break;
1117 }
1118 ✗ case 4: {
1119 // MONODOMAIN:
1120 // gamma_p = - lambda_p beta_p - Am mu_p - A_m a sum_{i} beta_i B_{ip} - \chi sum_{ij} beta_i beta_j T_{ijp} + (a+1)Am sum_{ij} beta_i beta_j T^s_ijp - Am sum_{ijk} beta_i beta_j beta_k Y_ijkp
1121 // eta_p = eps * ro * beta_p - eps * mu_p
1122
1123 // BIDOMAIN:
1124 // gamma_p = - lambda_p beta_p - lambda_p xi_p- Am mu_p - A_m a sum_{i} beta_i B_{ip} - \chi sum_{ij} beta_i beta_j T_{ijp} - \chi sum_{ij} beta_i xi_j T_{ijp} + (a+1)Am sum_{ij} beta_i beta_j T^s_ijp - Am sum_{ijk} beta_i beta_j beta_k Y_ijkp
1125 // eta_p = eps * ro * beta_p - eps * mu_p
1126
1127 ✗ for (int p=0; p<m_dimRomBasis; p++) {
1128 ✗ m_gamma_extrap[p] = - m_eigenValue[p] * beta_extrap[p];
1129 ✗ m_gamma_extrap[p] += - Am * m_mu[p];
1130 ✗ m_gamma_extrap[p] += - m_eigenValue[p] * xi_extrap[p];
1131 ✗ for (int i=0; i<m_dimRomBasis; i++) {
1132 ✗ m_gamma_extrap[p] += - a*Am * beta_extrap[i] * m_tensorB[secondOrderGlobalIndex(i,p)];
1133 ✗ for (int j=0; j<m_dimRomBasis; j++) {
1134 ✗ m_gamma_extrap[p] += (a+1)*Am * beta_extrap[i] * beta_extrap[j] * m_tensorTs[thirdOrderGlobalIndex(i,j,p)];
1135 ✗ m_gamma_extrap[p] += - m_coefChi * beta_extrap[i] * beta_extrap[j] * m_tensorT[thirdOrderGlobalIndex(i,j,p)];
1136 ✗ m_gamma_extrap[p] += - m_coefChi * beta_extrap[i] * xi_extrap[j] * m_tensorT[thirdOrderGlobalIndex(i,j,p)];
1137 ✗ for (int k=0; k<m_dimRomBasis; k++) {
1138 ✗ m_gamma_extrap[p] += - Am * beta_extrap[i] * beta_extrap[j] * beta_extrap[k] * m_tensorY[fourthOrderGlobalIndex(i,j,k,p)];
1139 }
1140 }
1141 }
1142
1143 ✗ if (FelisceParam::instance().hasSource) {
1144 ✗ double delay[m_numSource];
1145 ✗ double tPeriod=fmod(m_fstransient->time,FelisceParam::instance().timePeriod);
1146 ✗ for (int iS=0; iS<m_numSource; iS++) {
1147 ✗ delay[iS] = FelisceParam::instance().delayStim[iS];
1148 ✗ if ((tPeriod>=delay[iS])&&(tPeriod<=(delay[iS]+FelisceParam::instance().stimTime)))
1149 ✗ m_gamma_extrap[p] += Am*m_source[iS][p];
1150 }
1151 }
1152
1153 ✗ m_gamma_extrap[p] = m_gamma_extrap[p] / (Am*Cm);
1154
1155 ✗ m_eta_extrap[p] = eps * ( gam * beta_extrap[p] - mu_extrap[p] );
1156
1157 }
1158 ✗ break;
1159 }
1160 ✗ default:
1161 ✗ break;
1162 }
1163
1164 ✗ delete[] beta_extrap;
1165 ✗ delete[] mu_extrap;
1166 ✗ delete[] xi_extrap;
1167
1168 }
1169
1170 ✗ void EigenProblemALP::computeMatrixZeta() {
1171 ✗ if (FelisceParam::verbose() > 10 ) PetscPrintf(PETSC_COMM_WORLD,"\nEigenProblemALP::computeMatrixZeta()\n");
1172 ✗ if ( m_zeta.size() == 0 ) {
1173 ✗ m_size2 = static_cast<int>(m_dimRomBasis*(m_dimRomBasis+1)/2);
1174 ✗ m_zeta.resize(m_size2);
1175 ✗ for (int i=0; i<m_size2; i++)
1176 ✗ m_zeta[i].duplicateFrom(m_basis[0]);
1177 }
1178
1179 ✗ for (int i=0; i<m_dimRomBasis; i++) { // row
1180 ✗ for (int j=0; j<i+1; j++) { // column
1181 // secondOrderGlobalIndex(i,j) = i*(i+1)/2 + j
1182 ✗ pointwiseMult(m_zeta[secondOrderGlobalIndex(i,j)],m_basis[i],m_basis[j]);
1183 }
1184 }
1185
1186 }
1187
1188 ✗ void EigenProblemALP::computeMatrixM() {
1189 ✗ if (FelisceParam::verbose() > 10 ) PetscPrintf(PETSC_COMM_WORLD,"\nEigenProblemALP::computeMatrixM()\n");
1190 ✗ if ( m_matrixM == nullptr ) {
1191 ✗ m_matrixM = new double*[m_dimRomBasis];
1192 ✗ for (int i=0; i<m_dimRomBasis; i++) {
1193 ✗ m_matrixM[i] = new double[m_dimRomBasis];
1194 }
1195 }
1196
1197 double coef;
1198 ✗ double epsilonLambda = 5.e-01; //1.e-01;
1199 ✗ for (int i=0; i<m_dimRomBasis; i++) {
1200 ✗ m_matrixM[i][i] = 0.;
1201 ✗ for (int j=0; j<i; j++) {
1202 ✗ m_matrixM[i][j] = 0.;
1203 ✗ if ( std::abs(m_eigenValue[i] - m_eigenValue[j] ) > epsilonLambda ) {
1204 ✗ coef = m_coefChi / ( m_eigenValue[j] - m_eigenValue[i] );
1205 ✗ for (int k=0; k < m_dimRomBasis; k++) {
1206 ✗ m_matrixM[i][j] += coef * m_gamma[k] * m_tensorT[thirdOrderGlobalIndex(k,i,j)];
1207 }
1208 }
1209 ✗ m_matrixM[j][i] = (-1.0) * m_matrixM[i][j];
1210 }
1211 }
1212 }
1213
1214 ✗ int EigenProblemALP::secondOrderGlobalIndex(int& i, int& j) {
1215 int id;
1216
1217 ✗ if (i >= j) {
1218 ✗ id = static_cast<int>(i*(i+1)/2+j);
1219 } else {
1220 ✗ id = static_cast<int>(j*(j+1)/2+i);
1221 }
1222
1223 ✗ return id;
1224 }
1225
1226 ✗ int EigenProblemALP::thirdOrderGlobalIndex(int& i, int& j, int& k) {
1227 int id;
1228 int s[3];
1229 ✗ s[0] = i;
1230 ✗ s[1] = j;
1231 ✗ s[2] = k;
1232 int temp;
1233 ✗ for (int ii=0; ii<3; ii++) {
1234 ✗ for (int jj=0; jj<2; jj++) {
1235 ✗ if (s[jj] < s[jj+1]) {
1236 ✗ temp = s[jj];
1237 ✗ s[jj] = s[jj+1];
1238 ✗ s[jj+1] = temp;
1239 }
1240 }
1241 }
1242 ✗ id = static_cast<int>(s[0]*(s[0]+1)*(s[0]+2)/6);
1243 ✗ id += static_cast<int>(s[1]*(s[1]+1)/2);
1244 ✗ id += s[2];
1245 ✗ return id;
1246 }
1247
1248 ✗ int EigenProblemALP::fourthOrderGlobalIndex(int& i, int& j, int& k, int& h) {
1249 int id;
1250 int s[4];
1251 ✗ s[0] = i;
1252 ✗ s[1] = j;
1253 ✗ s[2] = k;
1254 ✗ s[3] = h;
1255 int temp;
1256 ✗ for (int ii=0; ii<4; ii++) {
1257 ✗ for (int jj=0; jj<3; jj++) {
1258 ✗ if (s[jj] < s[jj+1]) {
1259 ✗ temp = s[jj];
1260 ✗ s[jj] = s[jj+1];
1261 ✗ s[jj+1] = temp;
1262 }
1263 }
1264 }
1265 ✗ id = static_cast<int>(s[0]*(s[0]+1)*(s[0]*s[0]+5*s[0]+6)/24);
1266 ✗ id += static_cast<int>(s[1]*(s[1]+1)*(s[1]+2)/6);
1267 ✗ id += static_cast<int>(s[2]*(s[2]+1)/2);
1268 ✗ id += s[3];
1269 ✗ return id;
1270 }
1271
1272 ✗ void EigenProblemALP::computeTensor() {
1273 ✗ if (FelisceParam::verbose() > 10 ) PetscPrintf(PETSC_COMM_WORLD,"\nEigenProblemALP::computeTensor()\n");
1274 ✗ int n = m_dimRomBasis;
1275
1276 ✗ switch (m_problem) {
1277 ✗ case 0: {
1278 // T_ijk = < phi_i, zeta_jk >
1279 ✗ m_size3 = static_cast<int>(n*(n+1)*(n+2)/6);
1280 ✗ if (m_tensorT == nullptr)
1281 ✗ m_tensorT = new double [m_size3];
1282
1283 double valueT;
1284 ✗ PetscVector phiG;
1285 ✗ phiG.duplicateFrom(m_basis[0]);
1286 ✗ for (int i=0; i<m_dimRomBasis; i++) {
1287 ✗ for (int j=0; j<i+1; j++) {
1288 // phiG = zeta_ij^T * G
1289 ✗ mult(m_Matrix[m_idG],m_zeta[secondOrderGlobalIndex(i,j)],phiG);
1290 ✗ for (int k=0; k<j+1; k++) {
1291 // valueT = zeta_ij^T * G * phi_k
1292 ✗ dot(m_basis[k],phiG,&valueT);
1293 ✗ m_tensorT[thirdOrderGlobalIndex(i,j,k)] = valueT;
1294 }
1295 }
1296 }
1297 ✗ phiG.destroy();
1298 ✗ break;
1299 }
1300 ✗ case 1:
1301 case 2: {
1302 ✗ switch (m_tensorOrder) {
1303 ✗ case 3: {
1304 // B_ij = <s(x) phi_i, phi_j > = Phi^T G^s Phi
1305 ✗ if (m_tensorB == nullptr)
1306 ✗ m_tensorB = new double [m_size2];
1307 // E_ij = < sigma grad(phi_i), grad(phi_j) > = Phi^T K Phi
1308 ✗ if (m_tensorE == nullptr)
1309 ✗ m_tensorE = new double [m_size2];
1310 ✗ if (m_problem == 2) { // bidomain
1311 // Q_ij = < (sigma_I+sigma_E) grad(phi_i), grad(phi_j) > = Phi^T K_IE Phi
1312 ✗ if (m_tensorQ == nullptr)
1313 ✗ m_tensorQ = new double [secondOrderGlobalIndex(m_dimRomBasis,m_dimRomBasis)+1];
1314 }
1315 // Es_ij = < s(x) sigma grad(phi_i), grad(phi_j) > = Phi^T K^s Phi
1316 ✗ if (m_tensorEs == nullptr)
1317 ✗ m_tensorEs = new double [m_size2];
1318 // T_ijk = < phi_i, zeta_jk > = Phi^T G Z
1319 ✗ m_size3 = static_cast<int>(n*(n+1)*(n+2)/6);
1320 ✗ if (m_tensorT == nullptr)
1321 ✗ m_tensorT = new double [m_size3];
1322 // U_ijk = < s(x) sigma grad(phi_i), grad(zeta_jk) > = Phi^T K^s Z
1323 ✗ if (m_tensorU == nullptr) {
1324 ✗ m_tensorU = new double*[m_dimRomBasis];
1325 ✗ for (int i=0; i<m_dimRomBasis; i++) {
1326 ✗ m_tensorU[i] = new double[m_size2];
1327 }
1328 }
1329
1330 double value;
1331 ✗ for (int i=0; i<m_dimRomBasis; i++) {
1332 // phiG = phi_i^T * G
1333 ✗ PetscVector phiG;
1334 ✗ phiG.duplicateFrom(m_basis[0]);
1335 ✗ mult(m_Matrix[m_idG],m_basis[i],phiG);
1336
1337 // phiGs = phi_i^T * G^s
1338 ✗ PetscVector phiGs;
1339 ✗ phiGs.duplicateFrom(m_basis[0]);
1340 ✗ mult(m_Matrix[m_idGs],m_basis[i],phiGs);
1341
1342 // phiK = phi_i^T * K
1343 ✗ PetscVector phiK;
1344 ✗ phiK.duplicateFrom(m_basis[0]);
1345 ✗ mult(m_Matrix[m_idK],m_basis[i],phiK);
1346
1347 // phiKs = phi_i^T * K^s
1348 ✗ PetscVector phiKs;
1349 ✗ phiKs.duplicateFrom(m_basis[0]);
1350 ✗ mult(m_Matrix[m_idKs],m_basis[i],phiKs);
1351
1352 // phiKie = phi_i^T * K_ie
1353 ✗ PetscVector phiKie;
1354 ✗ phiKie.duplicateFrom(m_basis[0]);
1355 ✗ mult(m_Matrix[m_idKie],m_basis[i],phiKie);
1356
1357 ✗ for (int j=0; j<i+1; j++) {
1358 // valueB = phi_i^T * G^s * phi_j
1359 ✗ dot(m_basis[j],phiGs,&value);
1360 ✗ m_tensorB[secondOrderGlobalIndex(i,j)] = value;
1361
1362 // valueE = phi_i^T * K * phi_j
1363 ✗ dot(m_basis[j],phiK,&value);
1364 ✗ m_tensorE[secondOrderGlobalIndex(i,j)] = value;
1365
1366 ✗ if (m_problem == 2) { //bidomain
1367 // valueQ = phi_i^T * K_ie * phi_j
1368 ✗ dot(m_basis[j],phiKie,&value);
1369 ✗ m_tensorQ[secondOrderGlobalIndex(i,j)] = value;
1370 }
1371
1372 // valueEs = phi_i^T * K^s * phi_j
1373 ✗ dot(m_basis[j],phiKs,&value);
1374 ✗ m_tensorEs[secondOrderGlobalIndex(i,j)] = value;
1375
1376 ✗ for (int k=0; k<j+1; k++) {
1377 // valueT = phi_i^T * G * zeta_jk
1378 ✗ dot(m_zeta[secondOrderGlobalIndex(j,k)],phiG,&value);
1379 ✗ m_tensorT[thirdOrderGlobalIndex(i,j,k)] = value;
1380 }
1381 }
1382 ✗ for (int j=0; j<m_dimRomBasis; j++) {
1383 ✗ for (int k=0; k<j+1; k++) {
1384 // valueU = phi_i^T * K^s * zeta_jk
1385 ✗ dot(m_zeta[secondOrderGlobalIndex(j,k)],phiKs,&value);
1386 ✗ m_tensorU[i][secondOrderGlobalIndex(j,k)] = value;
1387 }
1388 }
1389 ✗ value = 0.;
1390 ✗ if (m_problem == 2) { //bidomain
1391 // in order to impose \int(u_E)=0
1392 // valueQ = sum_m phi_i,l^T * G_lm
1393 ✗ phiG.sum(&value);
1394 ✗ m_tensorQ[secondOrderGlobalIndex(i,m_dimRomBasis)] = value;
1395 }
1396
1397 ✗ phiG.destroy();
1398 ✗ phiGs.destroy();
1399 ✗ phiK.destroy();
1400 ✗ phiKie.destroy();
1401 ✗ phiKs.destroy();
1402 }
1403 ✗ if (m_problem == 2) //bidomain
1404 ✗ m_tensorQ[secondOrderGlobalIndex(m_dimRomBasis,m_dimRomBasis)] =0.;
1405 ✗ break;
1406 }
1407 ✗ case 4: {
1408 // B_ij = < s phi_i, phi_j > = Phi^T G^s Phi
1409 ✗ if (m_tensorB == nullptr)
1410 ✗ m_tensorB = new double [m_size2];
1411 ✗ if (m_problem == 2) { // bidomain
1412 // E_ij = < sigma grad(phi_i), grad(phi_j) > = Phi^T K Phi
1413 ✗ if (m_tensorE == nullptr)
1414 ✗ m_tensorE = new double [m_size2];
1415 // Q_ij = < (sigma_i+sigma_e) phi_i, phi_j > = Phi^T K_ie Phi
1416 ✗ if (m_tensorQ == nullptr)
1417 ✗ m_tensorQ = new double [secondOrderGlobalIndex(m_dimRomBasis,m_dimRomBasis)+1];
1418 }
1419 // T_ijk = < phi_i, zeta_jk > = Phi^T G Z
1420 ✗ m_size3 = static_cast<int>(n*(n+1)*(n+2)/6);
1421 ✗ if (m_tensorT == nullptr)
1422 ✗ m_tensorT = new double [m_size3];
1423 // Ts_ijk = < s(x) phi_i, zeta_jk > = Phi^T G^s Z
1424 ✗ m_size3 = static_cast<int>(n*(n+1)*(n+2)/6);
1425 ✗ if (m_tensorTs == nullptr)
1426 ✗ m_tensorTs = new double [m_size3];
1427 // Y_ijkh = < s(x) zeta_ij, zeta_kh > = Z^T G^s Z
1428 ✗ m_size4 = static_cast<int>(n*(n+1)*(n*n+5*n+6)/24);
1429 ✗ if (m_tensorY == nullptr)
1430 ✗ m_tensorY = new double [m_size4];
1431
1432 double value;
1433 ✗ for (int i=0; i<m_dimRomBasis; i++) {
1434 // phiG = phi_i^T * G
1435 ✗ PetscVector phiG;
1436 ✗ phiG.duplicateFrom(m_basis[0]);
1437 ✗ mult(m_Matrix[m_idG],m_basis[i],phiG);
1438
1439 // phiGs = phi_i^T * Gs
1440 ✗ PetscVector phiGs;
1441 ✗ phiGs.duplicateFrom(m_basis[0]);
1442 ✗ mult(m_Matrix[m_idGs],m_basis[i],phiGs);
1443
1444 // phiK = phi_i^T * K
1445 ✗ PetscVector phiK;
1446 ✗ phiK.duplicateFrom(m_basis[0]);
1447 ✗ if (m_problem == 2)
1448 ✗ mult(m_Matrix[m_idK],m_basis[i],phiK);
1449
1450 // phiKie = phi_i^T * K_ie
1451 ✗ PetscVector phiKie;
1452 ✗ phiKie.duplicateFrom(m_basis[0]);
1453 ✗ if (m_problem ==2)
1454 ✗ mult(m_Matrix[m_idKie],m_basis[i],phiKie);
1455
1456 ✗ for (int j=0; j<i+1; j++) {
1457 // valueB = phi_i^T * Gs * phi_j
1458 ✗ dot(m_basis[j],phiGs,&value);
1459 ✗ m_tensorB[secondOrderGlobalIndex(i,j)] = value;
1460
1461 ✗ if (m_problem == 2) { // bidomain
1462 // valueE = phi_i^T * K * phi_j
1463 ✗ dot(m_basis[j],phiK,&value);
1464 ✗ m_tensorE[secondOrderGlobalIndex(i,j)] = value;
1465 // valueQ = phi_i^T * K_ie * phi_j
1466 ✗ dot(m_basis[j],phiKie,&value);
1467 ✗ m_tensorQ[secondOrderGlobalIndex(i,j)] = value;
1468 }
1469
1470 // zetaG = zeta_ij^T * G
1471 ✗ PetscVector zetaG;
1472 ✗ zetaG.duplicateFrom(m_basis[0]);
1473 ✗ mult(m_Matrix[m_idG],m_zeta[secondOrderGlobalIndex(i,j)],zetaG);
1474
1475 // zetaGs = zeta_ij^T * G^s
1476 ✗ PetscVector zetaGs;
1477 ✗ zetaGs.duplicateFrom(m_basis[0]);
1478 ✗ mult(m_Matrix[m_idGs],m_zeta[secondOrderGlobalIndex(i,j)],zetaGs);
1479
1480 ✗ for (int k=0; k<j+1; k++) {
1481 // valueT = zeta_ij^T * G * phi_k
1482 ✗ dot(m_basis[k],zetaG,&value);
1483 ✗ m_tensorT[thirdOrderGlobalIndex(i,j,k)] = value;
1484
1485 // valueTs = zeta_ij^T * Gs * phi_k
1486 ✗ dot(m_basis[k],zetaGs,&value);
1487 ✗ m_tensorTs[thirdOrderGlobalIndex(i,j,k)] = value;
1488
1489 ✗ for (int h=0; h<k+1; h++) {
1490 // valueY = zeta_ij^T * Gs * zeta_kh
1491 ✗ dot(m_zeta[secondOrderGlobalIndex(k,h)],zetaGs,&value);
1492 ✗ m_tensorY[fourthOrderGlobalIndex(i,j,k,h)] = value;
1493 }
1494 }
1495 ✗ zetaG.destroy();
1496 ✗ zetaGs.destroy();
1497 }
1498 ✗ value = 0.;
1499 // in order to impose \int(u_E)=0
1500 ✗ if (m_problem == 2) { //bidomain
1501 // valueQ = sum_m phi_i,l^T * G_lm
1502 ✗ phiG.sum(&value);
1503 ✗ m_tensorQ[secondOrderGlobalIndex(i,m_dimRomBasis)] = value;
1504 }
1505
1506 ✗ phiG.destroy();
1507 ✗ phiGs.destroy();
1508 ✗ phiK.destroy();
1509 ✗ phiKie.destroy();
1510 }
1511 ✗ if (m_problem == 2) //bidomain
1512 ✗ m_tensorQ[secondOrderGlobalIndex(m_dimRomBasis,m_dimRomBasis)] =0.;
1513
1514 ✗ break;
1515 }
1516 ✗ default:
1517 ✗ break;
1518 }
1519 ✗ break;
1520 }
1521 ✗ default:
1522 ✗ FEL_ERROR("Model not defined for ALP solver.");
1523 ✗ break;
1524 }
1525
1526 }
1527
1528 1 void EigenProblemALP::initializeSolution() {
1529
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
 1 if (FelisceParam::verbose() > 10 ) PetscPrintf(PETSC_COMM_WORLD,"\nEigenProblemALP::initializeSolution()\n");
1530
1531
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 if ( !m_solutionInitialized) {
1532
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 if (m_beta == nullptr)
1533
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 m_beta = new double[m_dimRomBasis];
1534 1 projectOnBasis(m_U_0,m_beta,true);
1535
1536
2/4
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 if ( (m_problem == 1) || (m_problem == 2) ) {
1537
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 if (m_mu == nullptr)
1538
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 m_mu = new double[m_dimRomBasis];
1539 1 projectOnBasis(m_W_0,m_mu,true);
1540 }
1541
1542
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 if (m_problem == 2) {
1543
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 if (m_xi == nullptr)
1544
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 m_xi = new double[m_dimRomBasis];
1545 1 projectOnBasis(m_Ue_0,m_xi,true);
1546
1547 }
1548
1549
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 1 if (m_method == 4) {
1550 ✗ if (m_beta_n_1 == nullptr)
1551 ✗ m_beta_n_1 = new double[m_dimRomBasis];
1552 ✗ if (m_mu_n_1 == nullptr)
1553 ✗ m_mu_n_1 = new double[m_dimRomBasis];
1554 ✗ if (m_xi_n_1 == nullptr)
1555 ✗ m_xi_n_1 = new double[m_dimRomBasis];
1556 ✗ for (int i=0; i<m_dimRomBasis; i++) {
1557 ✗ m_beta_n_1[i] = m_beta[i];
1558 ✗ m_mu_n_1[i] = m_mu[i];
1559 ✗ m_xi_n_1[i] = m_xi[i];
1560 }
1561 }
1562
1563
1564 1 m_solutionInitialized = true;
1565 }
1566
1567 1 }
1568
1569 ✗ void EigenProblemALP::setHeteroTauClose(std::vector<double>& valueTauClose) {
1570 IGNORE_UNUSED_ARGUMENT(valueTauClose);
1571 ✗ FEL_WARNING("Hetero TauClose for Classical ALP not implemented.");
1572 }
1573
1574 ✗ void EigenProblemALP::setFhNf0(std::vector<double>& valuef0) {
1575 ✗ m_FhNf0.createSeq(PETSC_COMM_SELF,m_numDof);
1576 ✗ m_FhNf0.setFromOptions();
1577
1578 // initialize SEQUENTIAL vectors (in order to use elemField)
1579 felInt ia;
1580 ✗ for (felInt i=0; i< m_numDof; i++) {
1581 ✗ ia = i;
1582 ✗ AOApplicationToPetsc(m_ao,1,&ia);
1583 ✗ m_FhNf0.setValues(1,&ia,&valuef0[i],INSERT_VALUES);
1584 }
1585 ✗ m_FhNf0.assembly();
1586
1587 }
1588
1589 ✗ void EigenProblemALP::setIapp(std::vector<double>& iApp, int& idS) {
1590 ✗ if (m_source == nullptr) {
1591 ✗ m_source = new double*[m_numSource];
1592 ✗ for (int i=0; i<m_numSource; i++) {
1593 ✗ m_source[i] = new double[m_dimRomBasis];
1594 }
1595 }
1596
1597 ✗ PetscVector RHS;
1598 ✗ RHS.createSeq(PETSC_COMM_SELF,m_numDof);
1599 ✗ RHS.setFromOptions();
1600 // initialize SEQUENTIAL vectors (in order to use elemField)
1601 felInt ia;
1602 ✗ for (felInt i=0; i< m_numDof; i++) {
1603 ✗ ia = i;
1604 ✗ AOApplicationToPetsc(m_ao,1,&ia);
1605 ✗ RHS.setValues(1,&ia,&iApp[i],INSERT_VALUES);
1606 }
1607 ✗ RHS.assembly();
1608
1609 ✗ projectOnBasis(RHS, m_source[idS], true);
1610
1611 ✗ RHS.destroy();
1612 }
1613
1614
1615 75 void EigenProblemALP::checkSolution(PetscVector& solution, PetscVector& refSol, double & error) {
1616
2/6
✓ Branch 1 taken 75 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 75 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 75 if (FelisceParam::verbose() > 10 ) PetscPrintf(PETSC_COMM_WORLD,"\nEigenProblemALP::checkSolution\n");
1617
1618
1/2
✓ Branch 1 taken 75 times.
✗ Branch 2 not taken.
 75 PetscVector tmpVec1;
1619
1/2
✓ Branch 1 taken 75 times.
✗ Branch 2 not taken.
 75 tmpVec1.duplicateFrom(refSol);
1620
1621
1/2
✓ Branch 1 taken 75 times.
✗ Branch 2 not taken.
 75 PetscVector tmpVec2;
1622
1/2
✓ Branch 1 taken 75 times.
✗ Branch 2 not taken.
 75 tmpVec2.duplicateFrom(refSol);
1623
1624
1/2
✓ Branch 1 taken 75 times.
✗ Branch 2 not taken.
 75 waxpy(tmpVec1,-1.0,refSol,solution);
1625
1626
1/2
✓ Branch 2 taken 75 times.
✗ Branch 3 not taken.
 75 mult(m_Matrix[m_idG],tmpVec1,tmpVec2);
1627
1628
1/2
✓ Branch 1 taken 75 times.
✗ Branch 2 not taken.
 75 dot(tmpVec1,tmpVec2,&error);
1629
1630
1/2
✓ Branch 1 taken 75 times.
✗ Branch 2 not taken.
 75 tmpVec1.destroy();
1631
1/2
✓ Branch 1 taken 75 times.
✗ Branch 2 not taken.
 75 tmpVec2.destroy();
1632 75 }
1633
1634 1 void EigenProblemALP::checkBasis(int size) {
1635
1636 1 if (size == 0)
1637 1 size = m_dimRomBasis;
1638
1639
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 if (!m_solutionInitialized) {
1640
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 initializeSolution();
1641 }
1642
1643 int rankProc;
1644
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 MPI_Comm_rank(m_petscComm,&rankProc);
1645
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 PetscVector tmpSol;
1646
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 tmpSol.duplicateFrom(m_U_0);
1647
1648 1 double* tmpSolToSave = nullptr;
1649
2/4
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 1 tmpSolToSave = new double[m_numDof];
1650
1651 // Check convergence on m_U_0
1652 1 double convergence[size];
1653
2/2
✓ Branch 0 taken 25 times.
✓ Branch 1 taken 1 times.
 26 for (int iBasis=0; iBasis < size; iBasis++) {
1654 25 double tmpVec[iBasis+1];
1655
2/2
✓ Branch 0 taken 325 times.
✓ Branch 1 taken 25 times.
 350 for (int jBasis=0; jBasis<iBasis+1; jBasis++) {
1656 325 tmpVec[jBasis] = m_beta[jBasis];
1657 }
1658
1659
1/2
✓ Branch 1 taken 25 times.
✗ Branch 2 not taken.
 25 projectOnDof(tmpVec,tmpSol,iBasis+1);
1660
1661
1/2
✓ Branch 1 taken 25 times.
✗ Branch 2 not taken.
 25 checkSolution(tmpSol,m_U_0,convergence[iBasis]);
1662
1663
1/2
✓ Branch 1 taken 25 times.
✗ Branch 2 not taken.
 25 fromVecToDoubleStar(tmpSolToSave, tmpSol, rankProc, 1);
1664
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 25 times.
 25 if (m_problem == 0) {
1665 ✗ writeEnsightVector(tmpSolToSave, iBasis, "uConv");
1666 }
1667
2/4
✓ Branch 0 taken 25 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 25 times.
✗ Branch 3 not taken.
 25 if ( (m_problem == 1) || (m_problem == 2) ) {
1668
2/4
✓ Branch 2 taken 25 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 25 times.
✗ Branch 6 not taken.
 25 writeEnsightVector(tmpSolToSave, iBasis, "potTransMembConv");
1669 }
1670
1671
1/2
✓ Branch 1 taken 25 times.
✗ Branch 2 not taken.
 25 tmpSol.axpy(-1.,m_U_0);
1672
1/2
✓ Branch 1 taken 25 times.
✗ Branch 2 not taken.
 25 tmpSol.abs();
1673 double maxU0;
1674
1/2
✓ Branch 1 taken 25 times.
✗ Branch 2 not taken.
 25 m_U_0.max(&maxU0);
1675 double minU0;
1676
1/2
✓ Branch 1 taken 25 times.
✗ Branch 2 not taken.
 25 m_U_0.min(&minU0);
1677 25 double deltaU0 = maxU0-minU0;
1678
1/2
✓ Branch 1 taken 25 times.
✗ Branch 2 not taken.
 25 tmpSol.scale(1./deltaU0);
1679
1/2
✓ Branch 1 taken 25 times.
✗ Branch 2 not taken.
 25 fromVecToDoubleStar(tmpSolToSave, tmpSol, rankProc, 1);
1680
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 25 times.
 25 if (m_problem == 0) {
1681 ✗ writeEnsightVector(tmpSolToSave, iBasis, "uError");
1682 }
1683
2/4
✓ Branch 0 taken 25 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 25 times.
✗ Branch 3 not taken.
 25 if ( (m_problem == 1) || (m_problem == 2) ) {
1684
2/4
✓ Branch 2 taken 25 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 25 times.
✗ Branch 6 not taken.
 25 writeEnsightVector(tmpSolToSave, iBasis, "potTransMembError");
1685 }
1686
1687
1/2
✓ Branch 1 taken 25 times.
✗ Branch 2 not taken.
 25 PetscPrintf(PETSC_COMM_WORLD, "convergence[%d] = %e \n", iBasis, convergence[iBasis]);
1688
1689 25 }
1690
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 1 std::string fileName = FelisceParam::instance().resultDir + "/convergence";
1691
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 1 viewALP(convergence,size,fileName);
1692
1693 // Check convergence on m_W_0
1694
2/4
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 if ( (m_problem == 1) || (m_problem == 2) ) {
1695 1 double convergenceW[size];
1696
2/2
✓ Branch 0 taken 25 times.
✓ Branch 1 taken 1 times.
 26 for (int iBasis=0; iBasis < size; iBasis++) {
1697 25 double tmpVec[iBasis+1];
1698
2/2
✓ Branch 0 taken 325 times.
✓ Branch 1 taken 25 times.
 350 for (int jBasis=0; jBasis<iBasis+1; jBasis++) {
1699 325 tmpVec[jBasis] = m_mu[jBasis];
1700 }
1701
1702
1/2
✓ Branch 1 taken 25 times.
✗ Branch 2 not taken.
 25 projectOnDof(tmpVec,tmpSol,iBasis+1);
1703
1704
1/2
✓ Branch 1 taken 25 times.
✗ Branch 2 not taken.
 25 checkSolution(tmpSol,m_W_0,convergenceW[iBasis]);
1705
1706
1/2
✓ Branch 1 taken 25 times.
✗ Branch 2 not taken.
 25 fromVecToDoubleStar(tmpSolToSave, tmpSol, rankProc, 1);
1707
2/4
✓ Branch 2 taken 25 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 25 times.
✗ Branch 6 not taken.
 25 writeEnsightVector(tmpSolToSave, iBasis, "ionicVarConv");
1708
1709
1/2
✓ Branch 1 taken 25 times.
✗ Branch 2 not taken.
 25 tmpSol.axpy(-1.,m_W_0);
1710
1/2
✓ Branch 1 taken 25 times.
✗ Branch 2 not taken.
 25 tmpSol.abs();
1711 double maxU0;
1712
1/2
✓ Branch 1 taken 25 times.
✗ Branch 2 not taken.
 25 m_W_0.max(&maxU0);
1713 double minU0;
1714
1/2
✓ Branch 1 taken 25 times.
✗ Branch 2 not taken.
 25 m_W_0.min(&minU0);
1715 25 double deltaU0 = maxU0-minU0;
1716
1/2
✓ Branch 1 taken 25 times.
✗ Branch 2 not taken.
 25 tmpSol.scale(1./deltaU0);
1717
1/2
✓ Branch 1 taken 25 times.
✗ Branch 2 not taken.
 25 fromVecToDoubleStar(tmpSolToSave, tmpSol, rankProc, 1);
1718
2/4
✓ Branch 2 taken 25 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 25 times.
✗ Branch 6 not taken.
 25 writeEnsightVector(tmpSolToSave, iBasis, "ionicVarError");
1719
1720
1/2
✓ Branch 1 taken 25 times.
✗ Branch 2 not taken.
 25 PetscPrintf(PETSC_COMM_WORLD, "convergenceW[%d] = %e \n", iBasis, convergenceW[iBasis]);
1721 25 }
1722
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 1 fileName = FelisceParam::instance().resultDir + "/convergenceW";
1723
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 1 viewALP(convergenceW,size,fileName);
1724 1 }
1725
1726 // Check convergence on m_Ue_0
1727
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 if (m_problem == 2) {
1728 1 double convergenceUe[size];
1729
2/2
✓ Branch 0 taken 25 times.
✓ Branch 1 taken 1 times.
 26 for (int iBasis=0; iBasis < size; iBasis++) {
1730 25 double tmpVec[iBasis+1];
1731
2/2
✓ Branch 0 taken 325 times.
✓ Branch 1 taken 25 times.
 350 for (int jBasis=0; jBasis<iBasis+1; jBasis++) {
1732 325 tmpVec[jBasis] = m_xi[jBasis];
1733 }
1734
1735
1/2
✓ Branch 1 taken 25 times.
✗ Branch 2 not taken.
 25 projectOnDof(tmpVec,tmpSol,iBasis+1);
1736
1737
1/2
✓ Branch 1 taken 25 times.
✗ Branch 2 not taken.
 25 checkSolution(tmpSol,m_Ue_0,convergenceUe[iBasis]);
1738
1739
1/2
✓ Branch 1 taken 25 times.
✗ Branch 2 not taken.
 25 fromVecToDoubleStar(tmpSolToSave, tmpSol, rankProc, 1);
1740
2/4
✓ Branch 2 taken 25 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 25 times.
✗ Branch 6 not taken.
 25 writeEnsightVector(tmpSolToSave, iBasis, "potExtraCellConv");
1741
1742
1/2
✓ Branch 1 taken 25 times.
✗ Branch 2 not taken.
 25 tmpSol.axpy(-1.,m_Ue_0);
1743
1/2
✓ Branch 1 taken 25 times.
✗ Branch 2 not taken.
 25 tmpSol.abs();
1744 double maxU0;
1745
1/2
✓ Branch 1 taken 25 times.
✗ Branch 2 not taken.
 25 m_Ue_0.max(&maxU0);
1746 double minU0;
1747
1/2
✓ Branch 1 taken 25 times.
✗ Branch 2 not taken.
 25 m_Ue_0.min(&minU0);
1748 25 double deltaU0 = maxU0-minU0;
1749
1/2
✓ Branch 1 taken 25 times.
✗ Branch 2 not taken.
 25 tmpSol.scale(1./deltaU0);
1750
1/2
✓ Branch 1 taken 25 times.
✗ Branch 2 not taken.
 25 fromVecToDoubleStar(tmpSolToSave, tmpSol, rankProc, 1);
1751
2/4
✓ Branch 2 taken 25 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 25 times.
✗ Branch 6 not taken.
 25 writeEnsightVector(tmpSolToSave, iBasis, "potExtraCellError");
1752
1753
1/2
✓ Branch 1 taken 25 times.
✗ Branch 2 not taken.
 25 PetscPrintf(PETSC_COMM_WORLD, "convergenceUe[%d] = %e \n", iBasis, convergenceUe[iBasis]);
1754 25 }
1755
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 1 fileName = FelisceParam::instance().resultDir + "/convergenceUe";
1756
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 1 viewALP(convergenceUe,size,fileName);
1757 1 }
1758
1759 1 std::vector<std::string> varNames;
1760
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 1 if (m_problem == 0) {
1761 ✗ varNames.resize(2);
1762 ✗ varNames[0] = "uConv";
1763 ✗ varNames[1] = "uError";
1764
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 1 } else if (m_problem == 1) {
1765 ✗ if (FelisceParam::instance().printIonicVar) {
1766 ✗ varNames.resize(4);
1767 ✗ varNames[0] = "potTransMembConv";
1768 ✗ varNames[1] = "ionicVarConv";
1769 ✗ varNames[2] = "potTransMembError";
1770 ✗ varNames[3] = "ionicVarError";
1771 } else {
1772 ✗ varNames.resize(2);
1773 ✗ varNames[0] = "potTransMembConv";
1774 ✗ varNames[1] = "potTransMembError";
1775 }
1776
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 } else if (m_problem == 2) {
1777
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
 1 if (FelisceParam::instance().printIonicVar) {
1778
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 varNames.resize(6);
1779
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 varNames[0] = "potTransMembConv";
1780
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 varNames[1] = "potExtraCellConv";
1781
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 varNames[2] = "ionicVarConv";
1782
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 varNames[3] = "potTransMembError";
1783
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 varNames[4] = "potExtraCellError";
1784
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 varNames[5] = "ionicVarError";
1785 } else {
1786 ✗ varNames.resize(4);
1787 ✗ varNames[0] = "potTransMembConv";
1788 ✗ varNames[1] = "potExtraCellConv";
1789 ✗ varNames[3] = "potTransMembError";
1790 ✗ varNames[4] = "potExtraCellError";
1791 }
1792 }
1793
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 if (rankProc == 0)
1794
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 1 writeEnsightCase(size, 1.,varNames);
1795
1796
1797
1798
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 delete [] tmpSolToSave;
1799
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 tmpSol.destroy();
1800
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 2 }
1801
1802 //Update functions
1803 ✗ void EigenProblemALP::updateBasis() {
1804 ✗ if (FelisceParam::verbose() > 10 ) PetscPrintf(PETSC_COMM_WORLD,"\nEigenProblemALP::updateBasis()\n");
1805
1806 ✗ std::vector<PetscVector> newPhi;
1807 ✗ newPhi.resize(m_dimRomBasis);
1808 ✗ for (int i=0; i<m_dimRomBasis; i++) {
1809 ✗ newPhi[i].duplicateFrom(m_basis[0]);
1810 ✗ newPhi[i].set(0.);
1811 ✗ for (int j=0; j<m_dimRomBasis; j++) {
1812 ✗ newPhi[i].axpy(- m_matrixM[j][i], m_basis[j]);
1813 }
1814 }
1815
1816 ✗ std::vector<PetscVector> rPhi;
1817
1818 ✗ if(m_useImprovedRec==true) {
1819 ✗ rPhi.resize(m_dimRomBasis);
1820
1821 ✗ int psiSize = m_dimOrthComp*m_dimOrthComp;
1822 double* Psi1;
1823 ✗ Psi1 = new double [psiSize];
1824 double* Psi2;
1825 ✗ Psi2 = new double [psiSize];
1826 double* Psi;
1827 ✗ Psi = new double [psiSize];
1828
1829 ✗ double valP1=0.0;
1830 ✗ double valP2=0.0;
1831 ✗ PetscVector phiK;
1832 ✗ phiK.duplicateFrom(m_basis[0]);
1833 ✗ for(int i=0; i<m_dimOrthComp; i++) {
1834 ✗ mult(m_Matrix[m_idK],m_orthComp[i],phiK);
1835 ✗ for(int j=0; j<=i; j++) {
1836 ✗ dot(m_orthComp[j],phiK,&valP1);
1837 ✗ Psi1[i*(m_dimOrthComp)+j] = valP1;
1838 ✗ Psi1[j*(m_dimOrthComp)+i] = valP1;
1839 }
1840 }
1841 ✗ phiK.destroy();
1842
1843 ✗ PetscVector tmpU;
1844 ✗ tmpU.duplicateFrom(m_basis[0]);
1845 ✗ tmpU.set(0.);
1846 ✗ for(int i=0; i<m_dimRomBasis; i++) {
1847 ✗ tmpU.axpy(m_beta[i],m_basis[i]);
1848 }
1849 ✗ PetscVector psiG;
1850 ✗ psiG.duplicateFrom(m_basis[0]);
1851 ✗ mult(m_Matrix[m_idG],tmpU,psiG);
1852
1853
1854
1855 ✗ for(int i=0; i<m_dimOrthComp; i++) {
1856 ✗ for(int j=0; j<=i; j++) {
1857 ✗ PetscVector zetaP;
1858 ✗ zetaP.duplicateFrom(m_basis[0]);
1859 ✗ pointwiseMult(zetaP,m_orthComp[i],m_orthComp[j]);
1860 ✗ dot(zetaP,psiG,&valP2);
1861 ✗ Psi2[i*(m_dimOrthComp)+j] = valP2;
1862 ✗ Psi2[j*(m_dimOrthComp)+i] = valP2;
1863
1864 ✗ zetaP.destroy();
1865 }
1866 }
1867 ✗ psiG.destroy();
1868 ✗ tmpU.destroy();
1869
1870 ✗ PetscVector tmpF;
1871 ✗ tmpF.duplicateFrom(m_basis[0]);
1872 ✗ tmpF.set(0.);
1873 ✗ for(int i=0; i<m_dimRomBasis; i++) {
1874 ✗ tmpF.axpy(m_gamma[i],m_basis[i]);
1875 }
1876
1877 ✗ double tol = 1.0e-9;
1878 ✗ int maxIt = 2*m_dimOrthComp;
1879 double* tmpRhs;
1880 ✗ tmpRhs = new double[m_dimOrthComp];
1881 double* alpha;
1882 ✗ alpha = new double[m_dimOrthComp];
1883
1884
1885 ✗ for(int i=0; i<m_dimRomBasis; i++) {
1886 ✗ rPhi[i].duplicateFrom(m_basis[0]);
1887 ✗ rPhi[i].set(0.);
1888
1889 ✗ PetscVector FPhiI;
1890 ✗ FPhiI.duplicateFrom(m_basis[0]);
1891 ✗ pointwiseMult(FPhiI,m_basis[i],tmpF);
1892
1893 ✗ PetscVector FPhiG;
1894 ✗ FPhiG.duplicateFrom(m_basis[0]);
1895 ✗ mult(m_Matrix[m_idG],FPhiI,FPhiG);
1896
1897 ✗ for (int l=0; l<m_dimOrthComp; l++) {
1898 ✗ double rhsVal=0.0;
1899 ✗ dot(m_orthComp[l],FPhiG,&rhsVal);
1900 ✗ tmpRhs[l] = m_coefChi * rhsVal;
1901 }
1902
1903 ✗ for (int l=0; l<m_dimOrthComp; l++) {
1904 ✗ for(int m=0; m<m_dimOrthComp; m++) {
1905 ✗ Psi[l*m_dimOrthComp+m] = Psi1[l*m_dimOrthComp+m] - m_coefChi*Psi2[l*m_dimOrthComp+m];
1906 ✗ if(m==l) {
1907 ✗ Psi[l*m_dimOrthComp + l] = Psi[l*m_dimOrthComp + l] - m_eigenValue[i];
1908 }
1909 }
1910 }
1911 ✗ solveOrthComp(Psi, tmpRhs, alpha, tol, maxIt);
1912
1913 ✗ for(int l=0; l<m_dimOrthComp; l++) {
1914 ✗ rPhi[i].axpy(alpha[l],m_orthComp[l]);
1915 }
1916
1917 ✗ FPhiI.destroy();
1918 ✗ FPhiG.destroy();
1919 }
1920 ✗ tmpF.destroy();
1921
1922 ✗ delete [] Psi1;
1923 ✗ delete [] Psi2;
1924 ✗ delete [] Psi;
1925 ✗ delete [] tmpRhs;
1926 ✗ delete [] alpha;
1927
1928
1929 }
1930
1931
1932 ✗ double dt = FelisceParam::instance().timeStep;
1933 ✗ for (int i=0; i<m_dimRomBasis; i++) {
1934 ✗ m_basis[i].axpy(dt,newPhi[i]);
1935 ✗ newPhi[i].destroy();
1936
1937 ✗ if(m_useImprovedRec==true) {
1938 ✗ m_basis[i].axpy(dt,rPhi[i]);
1939 ✗ rPhi[i].destroy();
1940 }
1941 }
1942
1943 ✗ MGS();
1944 ✗ if(m_useImprovedRec==true) {
1945 ✗ MGS(m_orthComp);
1946 }
1947
1948 ✗ if (FelisceParam::instance().writeBasisEvolution)
1949 ✗ writeMode(m_fstransient->iteration);
1950
1951 }
1952
1953
1954 ✗ void EigenProblemALP::setPotential()
1955 {
1956 ✗ if (FelisceParam::verbose() > 10 ) PetscPrintf(PETSC_COMM_WORLD,"\nEigenProblemALP::setPotential()\n");
1957
1958 ✗ int Nm = FelisceParam::instance().numApproxMode;
1959 ✗ int nCutOff = FelisceParam::instance().nCutOff;
1960
1961 ✗ m_initPot.duplicateFrom(m_basis[0]);
1962 ✗ m_initPot.copyFrom(m_U_0);
1963
1964 ✗ int sizePb = -1;
1965 ✗ std::vector<PetscVector> lagMult;
1966 ✗ std::vector<PetscVector> errorRep;
1967 ✗ if (m_problem == 0)
1968 { // FKPP
1969 ✗ sizePb = 1;
1970 }
1971 ✗ else if (m_problem == 1)
1972 { // Monodomain
1973 ✗ sizePb = 2;
1974 }
1975 ✗ else if (m_problem == 2)
1976 { // Bidomain
1977 ✗ sizePb = 3;
1978 }
1979 else
1980 ✗ FEL_ERROR("Problem not found.");
1981 ✗ lagMult.resize(sizePb);
1982 ✗ errorRep.resize(sizePb);
1983
1984 ✗ for (int i=0; i<sizePb; i++) {
1985 ✗ lagMult[i].duplicateFrom(m_basis[0]);
1986 ✗ errorRep[i].duplicateFrom(m_basis[0]);
1987 }
1988
1989 ✗ PetscVector tmpVec;
1990 ✗ tmpVec.duplicateFrom(m_basis[0]);
1991
1992 ✗ double errorU = 1.;
1993 ✗ double errorW = 1.;
1994 ✗ double errorUe = 1.;
1995
1996 ✗ std::vector<double*> coeffBeta;
1997 ✗ coeffBeta.resize(sizePb);
1998
1999 ✗ if (m_beta == nullptr)
2000 ✗ m_beta = new double[Nm];
2001 ✗ projectOnBasis(m_U_0,m_beta,true,Nm);
2002 ✗ m_solutionInitialized = true;
2003 ✗ coeffBeta[0] = m_beta;
2004
2005 ✗ errorRep[0].copyFrom(m_U_0);
2006 ✗ for (int i=0; i<Nm; i++) {
2007 ✗ double coef = -m_beta[i];
2008 ✗ errorRep[0].axpy(coef,m_basis[i]);
2009 }
2010 ✗ mult(m_Matrix[m_idG],errorRep[0],tmpVec);
2011 ✗ dot(errorRep[0],tmpVec,&errorU);
2012 // errorU = std::sqrt(errorU)/Nm;
2013
2014 ✗ if (m_problem > 0) {
2015 ✗ if (m_mu == nullptr)
2016 ✗ m_mu = new double[Nm];
2017 ✗ projectOnBasis(m_W_0,m_mu,true,Nm);
2018 ✗ coeffBeta[1] = m_mu;
2019
2020 ✗ errorRep[1].copyFrom(m_W_0);
2021 ✗ for (int i=0; i<Nm; i++) {
2022 ✗ double coef = -m_mu[i];
2023 ✗ errorRep[1].axpy(coef,m_basis[i]);
2024 }
2025 ✗ mult(m_Matrix[m_idG],errorRep[1],tmpVec);
2026 ✗ dot(errorRep[1],tmpVec,&errorW);
2027 // errorW = std::sqrt(errorW)/Nm;
2028 }
2029 ✗ if (m_problem == 2) {
2030 ✗ if (m_xi == nullptr)
2031 ✗ m_xi = new double[Nm];
2032 ✗ projectOnBasis(m_Ue_0,m_xi,true,Nm);
2033 ✗ coeffBeta[2] = m_xi;
2034
2035 ✗ errorRep[2].copyFrom(m_Ue_0);
2036 ✗ for (int i=0; i<Nm; i++) {
2037 ✗ double coef = -m_xi[i];
2038 ✗ errorRep[2].axpy(coef,m_basis[i]);
2039 }
2040 ✗ mult(m_Matrix[m_idG],errorRep[2],tmpVec);
2041 ✗ dot(errorRep[2],tmpVec,&errorUe);
2042 // errorUe = std::sqrt(errorUe)/Nm;
2043 }
2044
2045 ✗ for (int i=0; i<sizePb; i++) {
2046 ✗ lagMult[i].copyFrom(errorRep[i]);
2047 ✗ lagMult[i].scale(-1.);
2048 }
2049
2050 ✗ PetscVector dLdU;
2051 ✗ dLdU.duplicateFrom(m_basis[0]);
2052
2053 ✗ PetscVector phi_ij;
2054 ✗ phi_ij.duplicateFrom(m_basis[0]);
2055
2056 ✗ PetscVector dUphi;
2057 ✗ dUphi.duplicateFrom(m_basis[0]);
2058
2059 ✗ PetscVector dU;
2060 ✗ dU.duplicateFrom(m_basis[0]);
2061
2062 ✗ PetscVector dUphiTG;
2063 ✗ dUphiTG.duplicateFrom(m_basis[0]);
2064 ✗ std::vector<PetscVector> lKeps;
2065 ✗ lKeps.resize(sizePb);
2066 ✗ for (int i=0; i<sizePb; i++) {
2067 ✗ lKeps[i].duplicateFrom(m_basis[0]);
2068 }
2069
2070 ✗ int numIter = 0;
2071 ✗ double step = 0.25; // 1.e-1;
2072 std::array<double, 3> k;
2073 ✗ k[0] = 50.;
2074 ✗ k[1] = 10.;
2075 ✗ k[2] = 10.;
2076
2077 ✗ double** q = new double* [sizePb];
2078 ✗ for (int iPb=0; iPb<sizePb; iPb++) {
2079 ✗ q[iPb] = new double[nCutOff];
2080 }
2081 ✗ double** theta = new double* [nCutOff];
2082 ✗ for (int i=0; i<nCutOff; i++) {
2083 ✗ theta[i] = new double[nCutOff];
2084 }
2085
2086 ✗ PetscPrintf(PETSC_COMM_WORLD,"numIter = %i \n",numIter);
2087 ✗ PetscPrintf(PETSC_COMM_WORLD,"error = %e, %e, %e \n",errorU,errorW,errorUe);
2088
2089 ✗ while ( (errorU > 1.e-08) && (errorW > 1.e-08) && (errorUe > 1.e-08) && (numIter < 5000) )
2090 {
2091 // q[iPb]_j = <lagMult[iPb]+k*eps[iPb], phi_j>
2092 ✗ for (int iPb=0; iPb<sizePb; iPb++) {
2093 ✗ lKeps[iPb].copyFrom(lagMult[iPb]);
2094 ✗ lKeps[iPb].axpy(k[iPb],errorRep[iPb]);
2095 ✗ projectOnBasis(lKeps[iPb], q[iPb], true, nCutOff);
2096 }
2097
2098 // dL/dU = (U - u0)
2099 ✗ dLdU.copyFrom(m_initPot);
2100 ✗ dLdU.axpy(-1.,m_U_0);
2101 ✗ double reg = 1.;
2102 ✗ dLdU.scale(reg);
2103 // dL/dU -= \chi * \sum \sum ( sum_pb(beta_i q_j) / (lambda_j - lambda_i) phi_i phi_j
2104 ✗ for (int i=0; i<Nm; i++) {
2105 ✗ for (int j=0; j<nCutOff; j++) {
2106 ✗ if ( std::fabs(m_eigenValue[j] - m_eigenValue[i]) > 1.e-1 ) {
2107 ✗ pointwiseMult(phi_ij, m_basis[i], m_basis[j]);
2108 ✗ double coef = 0.;
2109 ✗ for (int iPb=0; iPb<sizePb; iPb++) {
2110 ✗ coef += coeffBeta[iPb][i]*q[iPb][j];
2111 }
2112 ✗ coef = - m_coefChi * coef/(m_eigenValue[j] - m_eigenValue[i]);
2113 ✗ dLdU.axpy(coef,phi_ij);
2114 }
2115 }
2116 }
2117
2118 // dU = -s dL/dU
2119 ✗ dU.set(0.);
2120 ✗ dU.axpy(-step,dLdU);
2121
2122 // theta_ij = <dU phi_i, phi_j>
2123 ✗ for (int i=0; i<nCutOff; i++) {
2124 ✗ pointwiseMult(dUphi, dU, m_basis[i]);
2125 ✗ mult(m_Matrix[m_idG], dUphi, dUphiTG);
2126 ✗ for (int j=0; j<=i; j++) {
2127 ✗ dot(dUphiTG, m_basis[j], &theta[i][j]);
2128 ✗ if (i != j)
2129 ✗ theta[j][i] = theta[i][j];
2130 }
2131 }
2132
2133 // Update phi
2134 ✗ for (int i=0; i<nCutOff; i++) {
2135 ✗ for (int j=0; j<nCutOff; j++) {
2136 ✗ if ( std::fabs(m_eigenValue[j] - m_eigenValue[i]) > 1.e-1 ) {
2137 ✗ double coef = m_coefChi * theta[i][j]/(m_eigenValue[j] - m_eigenValue[i]);
2138 ✗ m_basis[i].axpy(coef, m_basis[j]);
2139 }
2140 }
2141 }
2142 ✗ MGS(nCutOff);
2143
2144 // Update lambda
2145 ✗ for (int i=0; i<nCutOff; i++) {
2146 ✗ m_eigenValue[i] -= m_coefChi * theta[i][i];
2147 }
2148
2149 // Update beta
2150 ✗ projectOnBasis(m_U_0,m_beta,true,Nm);
2151 ✗ if (m_problem > 0) {
2152 // Update mu
2153 ✗ projectOnBasis(m_W_0,m_mu,true,Nm);
2154 }
2155 ✗ if (m_problem == 2) {
2156 // Update xi
2157 ✗ projectOnBasis(m_Ue_0,m_xi,true,Nm);
2158 }
2159
2160 // Compute error
2161 ✗ errorRep[0].copyFrom(m_U_0);
2162 ✗ for (int i=0; i<Nm; i++) {
2163 ✗ errorRep[0].axpy(-1.*m_beta[i],m_basis[i]);
2164 }
2165 ✗ mult(m_Matrix[m_idG],errorRep[0],tmpVec);
2166 ✗ dot(errorRep[0],tmpVec,&errorU);
2167 // errorU = std::sqrt(errorU)/Nm;
2168 ✗ if (m_problem > 0) {
2169 ✗ errorRep[1].copyFrom(m_W_0);
2170 ✗ for (int i=0; i<Nm; i++) {
2171 ✗ errorRep[1].axpy(-1.*m_mu[i],m_basis[i]);
2172 }
2173 ✗ mult(m_Matrix[m_idG],errorRep[1],tmpVec);
2174 ✗ dot(errorRep[1],tmpVec,&errorW);
2175 // errorW = std::sqrt(errorW)/Nm;
2176 }
2177 ✗ if (m_problem == 2) {
2178 ✗ errorRep[2].copyFrom(m_Ue_0);
2179 ✗ for (int i=0; i<Nm; i++) {
2180 ✗ errorRep[2].axpy(-1.*m_xi[i],m_basis[i]);
2181 }
2182 ✗ mult(m_Matrix[m_idG],errorRep[2],tmpVec);
2183 ✗ dot(errorRep[2],tmpVec,&errorUe);
2184 // errorUe = std::sqrt(errorUe)/Nm;
2185 }
2186
2187 // Update Lagrangian Multiplicator
2188 ✗ for (int iPb=0; iPb<sizePb; iPb++) {
2189 ✗ lagMult[iPb].axpy(step, errorRep[iPb]);
2190 }
2191
2192 // Update initial potential
2193 ✗ m_initPot.axpy(1.,dU);
2194
2195 ✗ numIter++;
2196
2197 ✗ PetscPrintf(PETSC_COMM_WORLD,"numIter = %i \n",numIter);
2198 ✗ PetscPrintf(PETSC_COMM_WORLD,"error = %e, %e, %e \n",errorU,errorW,errorUe);
2199 }
2200
2201 ✗ for (int iPb=0; iPb<sizePb; iPb++) {
2202 ✗ lagMult[iPb].destroy();
2203 ✗ errorRep[iPb].destroy();
2204 ✗ lKeps[iPb].destroy();
2205 }
2206 ✗ dLdU.destroy();
2207 ✗ phi_ij.destroy();
2208 ✗ dUphi.destroy();
2209 ✗ dU.destroy();
2210 ✗ dUphiTG.destroy();
2211 ✗ tmpVec.destroy();
2212
2213 ✗ for (int i=0; i<sizePb; i++) {
2214 ✗ delete [] q[i];
2215 }
2216 ✗ delete [] q;
2217 ✗ for (int i=0; i<nCutOff; i++) {
2218 ✗ delete [] theta[i];
2219 }
2220 ✗ delete [] theta;
2221
2222 // Write initial potential in Ensight file
2223 ✗ double* tmpSolToSave = nullptr;
2224 ✗ tmpSolToSave = new double[m_numDof];
2225 int rankProc;
2226 ✗ MPI_Comm_rank(m_petscComm,&rankProc);
2227 ✗ fromVecToDoubleStar(tmpSolToSave, m_initPot, rankProc, 1);
2228 ✗ writeEnsightVector(tmpSolToSave, 0, "initialPotential");
2229 ✗ delete [] tmpSolToSave;
2230
2231 ✗ std::string fileName = FelisceParam::instance().resultDir + "/eigenValue";
2232 ✗ viewALP(m_eigenValue,m_dimRomBasis,fileName);
2233
2234 ✗ checkBasis(Nm);
2235 }
2236
2237 // Modified Gram-Schmidt (see Golub & Van Loan)
2238 ✗ void EigenProblemALP::MGS(int size)
2239 {
2240 ✗ if (FelisceParam::verbose() > 10 ) PetscPrintf(PETSC_COMM_WORLD,"\nEigenProblemALP::MGS(int)\n");
2241
2242 ✗ if (size == 0)
2243 ✗ size = m_dimRomBasis;
2244
2245 ✗ for(int iBasis=0; iBasis<size; iBasis++)
2246 {
2247 ✗ PetscVector tmpVec;
2248 ✗ tmpVec.duplicateFrom(m_basis[0]);
2249 ✗ mult(m_Matrix[m_idG],m_basis[iBasis],tmpVec);
2250 ✗ double rii = 0.0;
2251 ✗ dot(tmpVec,m_basis[iBasis],&rii);
2252 ✗ rii = 1.0/std::sqrt(rii);
2253 ✗ m_basis[iBasis].scale(rii);
2254
2255 ✗ for(int j=iBasis+1; j<size; j++) {
2256 ✗ double rij = 0.0;
2257 ✗ dot(tmpVec,m_basis[j],&rij);
2258 ✗ rij = -1.0 * rij;
2259 ✗ m_basis[j].axpy(rij,m_basis[iBasis]);
2260 }
2261 ✗ tmpVec.destroy();
2262 }
2263 }
2264
2265
2266 // Modified Gram-Schmidt (for an extra-space, orthonormal to the basis too)
2267 // 1 - orthonormalize w.r.t. basis, 2 - orthonormalize extraspace;
2268 ✗ void EigenProblemALP::MGS(std::vector<PetscVector>& extraSpace) {
2269 ✗ if (FelisceParam::verbose() > 10 ) PetscPrintf(PETSC_COMM_WORLD,"\nEigenProblemALP::MGS(std::vector<Vec>)\n");
2270
2271 ✗ for(int j=0; j<m_dimOrthComp; j++) {
2272 ✗ PetscVector tmpVec;
2273 ✗ tmpVec.duplicateFrom(m_basis[0]);
2274 ✗ mult(m_Matrix[m_idG],extraSpace[j],tmpVec);
2275 double* Sij;
2276 ✗ Sij = new double[m_dimRomBasis];
2277 ✗ for(int i =0; i<m_dimRomBasis; i++) {
2278 ✗ double sij = 0.0;
2279 ✗ dot(tmpVec,m_basis[i],&sij);
2280 ✗ Sij[i] = -1.0*sij;
2281 }
2282 double* Rjh;
2283 ✗ Rjh = new double[j];
2284 ✗ for(int h=0; h<j; h++) {
2285 ✗ double rjh = 0.0;
2286 ✗ dot(tmpVec,extraSpace[h],&rjh);
2287 ✗ Rjh[h] = -1.0*rjh;
2288 }
2289
2290 ✗ for(int i=0; i<m_dimRomBasis; i++) {
2291 ✗ double wheight = Sij[i];
2292 ✗ extraSpace[j].axpy(wheight,m_basis[i]);
2293 }
2294 ✗ for(int h=0; h<j; h++) {
2295 ✗ double wheight = Rjh[h];
2296 ✗ extraSpace[j].axpy(wheight,extraSpace[h]);
2297 }
2298
2299 ✗ PetscVector tmpNorm;
2300 ✗ tmpNorm.duplicateFrom(m_basis[0]);
2301 ✗ mult(m_Matrix[m_idG],extraSpace[j],tmpNorm);
2302 ✗ double rjj = 0.0;
2303 ✗ dot(tmpNorm,extraSpace[j],&rjj);
2304 ✗ rjj = 1.0/std::sqrt(rjj);
2305 ✗ extraSpace[j].scale(rjj);
2306
2307 ✗ tmpVec.destroy();
2308 ✗ tmpNorm.destroy();
2309 ✗ delete [] Rjh;
2310 ✗ delete [] Sij;
2311 }
2312 }
2313
2314 ✗ void EigenProblemALP::updateTensor() {
2315 ✗ if (FelisceParam::verbose() > 10 ) PetscPrintf(PETSC_COMM_WORLD,"\nEigenProblemALP::updateTensor()\n");
2316
2317 ✗ double dt = FelisceParam::instance().timeStep;
2318
2319 ✗ switch (m_problem) {
2320 ✗ case 0: {
2321 ✗ double deltaT[m_size3];
2322 // T^{n+1} = T^n + dt * [T^n, M^n]^{(3)}
2323 ✗ for (int i=0; i<m_dimRomBasis; i++) {
2324 ✗ for (int j=0; j<i+1; j++) {
2325 ✗ for (int k=0; k<j+1; k++) {
2326 ✗ deltaT[thirdOrderGlobalIndex(i,j,k)] = 0.;
2327 ✗ for (int l=0; l<m_dimRomBasis; l++) {
2328 // delta T_ijk = sum_l (M_il * T_ljk + M_jl * T_ilk + M_kl * T_ijl )
2329 ✗ deltaT[thirdOrderGlobalIndex(i,j,k)] += m_matrixM[i][l] * m_tensorT[thirdOrderGlobalIndex(l,j,k)];
2330 ✗ deltaT[thirdOrderGlobalIndex(i,j,k)] += m_matrixM[j][l] * m_tensorT[thirdOrderGlobalIndex(i,l,k)];
2331 ✗ deltaT[thirdOrderGlobalIndex(i,j,k)] += m_matrixM[k][l] * m_tensorT[thirdOrderGlobalIndex(i,j,l)];
2332 }
2333 }
2334 }
2335 }
2336 ✗ for (int idT=0; idT<m_size3; idT++)
2337 ✗ m_tensorT[idT] += dt * deltaT[idT];
2338 ✗ break;
2339 }
2340 ✗ case 1:
2341 case 2: {
2342 ✗ switch (m_tensorOrder) {
2343 ✗ case 3: {
2344 ✗ double deltaB[m_size2];
2345 // B^{n+1} = B^n + dt * [B^n, M^n]^{(2)}
2346 ✗ double deltaE[m_size2];
2347 // E^{n+1} = E^n + dt * [E^n, M^n]^{(2)}
2348 ✗ double deltaEs[m_size2];
2349 // Es^{n+1} = Es^n + dt * [Es^n, M^n]^{(2)}
2350 ✗ double deltaQ[m_size2+m_dimRomBasis];
2351 // Q^{n+1} = Q^n + dt * [Q^n, M^n]^{(2)}
2352 ✗ double deltaT[m_size3];
2353 // T^{n+1} = T^n + dt * [T^n, M^n]^{(3)}
2354 ✗ double** deltaU = new double*[m_dimRomBasis];
2355 ✗ for (int i=0; i<m_dimRomBasis; i++)
2356 ✗ deltaU[i] = new double[m_size2];
2357 // U^{n+1} = U^n + dt * [U^n, M^n]^{(3)}
2358
2359 ✗ for (int i=0; i<m_dimRomBasis; i++) {
2360 ✗ for (int j=0; j<i+1; j++) {
2361 ✗ deltaB[secondOrderGlobalIndex(i,j)] = 0.;
2362 ✗ deltaE[secondOrderGlobalIndex(i,j)] = 0.;
2363 ✗ deltaEs[secondOrderGlobalIndex(i,j)] = 0.;
2364 ✗ if (m_problem == 2) // bidomain
2365 ✗ deltaQ[secondOrderGlobalIndex(i,j)] = 0.;
2366
2367 ✗ for (int l=0; l<m_dimRomBasis; l++) {
2368 // delta B_ij = sum_l ( M_il * B_lj + M_jl * B_il )
2369 ✗ deltaB[secondOrderGlobalIndex(i,j)] += m_matrixM[i][l] * m_tensorB[secondOrderGlobalIndex(l,j)];
2370 ✗ deltaB[secondOrderGlobalIndex(i,j)] += m_matrixM[j][l] * m_tensorB[secondOrderGlobalIndex(i,l)];
2371 // delta E_ij = sum_l ( M_il * E_lj + M_jl * E_il )
2372 ✗ deltaE[secondOrderGlobalIndex(i,j)] += m_matrixM[i][l] * m_tensorE[secondOrderGlobalIndex(l,j)];
2373 ✗ deltaE[secondOrderGlobalIndex(i,j)] += m_matrixM[j][l] * m_tensorE[secondOrderGlobalIndex(i,l)];
2374 // delta Es_ij = sum_l ( M_il * Es_lj + M_jl * Es_il )
2375 ✗ deltaEs[secondOrderGlobalIndex(i,j)] += m_matrixM[i][l] * m_tensorEs[secondOrderGlobalIndex(l,j)];
2376 ✗ deltaEs[secondOrderGlobalIndex(i,j)] += m_matrixM[j][l] * m_tensorEs[secondOrderGlobalIndex(i,l)];
2377 ✗ if (m_problem == 2) { // bidomain
2378 // delta Q_ij = sum_l ( M_il * Q_lj + M_jl * Q_il )
2379 ✗ deltaQ[secondOrderGlobalIndex(i,j)] += m_matrixM[i][l] * m_tensorQ[secondOrderGlobalIndex(l,j)];
2380 ✗ deltaQ[secondOrderGlobalIndex(i,j)] += m_matrixM[j][l] * m_tensorQ[secondOrderGlobalIndex(i,l)];
2381 }
2382
2383 }
2384 ✗ for (int k=0; k<j+1; k++) {
2385 ✗ deltaT[thirdOrderGlobalIndex(i,j,k)] = 0.;
2386 ✗ for (int l=0; l<m_dimRomBasis; l++) {
2387 // delta T_ijk = sum_l (M_il * T_ljk + M_jl * T_ilk + M_kl * T_ijl )
2388 ✗ deltaT[thirdOrderGlobalIndex(i,j,k)] += m_matrixM[i][l] * m_tensorT[thirdOrderGlobalIndex(l,j,k)];
2389 ✗ deltaT[thirdOrderGlobalIndex(i,j,k)] += m_matrixM[j][l] * m_tensorT[thirdOrderGlobalIndex(i,l,k)];
2390 ✗ deltaT[thirdOrderGlobalIndex(i,j,k)] += m_matrixM[k][l] * m_tensorT[thirdOrderGlobalIndex(i,j,l)];
2391 }
2392 }
2393 }
2394 ✗ for (int j=0; j<m_dimRomBasis; j++) {
2395 ✗ for (int k=0; k<j+1; k++) {
2396 ✗ deltaU[i][secondOrderGlobalIndex(j,k)] = 0.;
2397 ✗ for (int l=0; l<m_dimRomBasis; l++) {
2398 // delta U_ijk = sum_l (M_il * U_ljk + M_jl * U_ilk + M_kl * U_ijl )
2399 ✗ deltaU[i][secondOrderGlobalIndex(j,k)] += m_matrixM[i][l] * m_tensorU[l][secondOrderGlobalIndex(j,k)];
2400 ✗ deltaU[i][secondOrderGlobalIndex(j,k)] += m_matrixM[j][l] * m_tensorU[i][secondOrderGlobalIndex(l,k)];
2401 ✗ deltaU[i][secondOrderGlobalIndex(j,k)] += m_matrixM[k][l] * m_tensorU[i][secondOrderGlobalIndex(j,l)];
2402 }
2403 }
2404 }
2405 ✗ if (m_problem == 2) {
2406 ✗ deltaQ[secondOrderGlobalIndex(m_dimRomBasis,i)] = 0.;
2407 ✗ for (int l=0; l<m_dimRomBasis; l++) {
2408 // delta Q_Nj = sum_l ( M_jl * Q_lj )
2409 ✗ deltaQ[secondOrderGlobalIndex(m_dimRomBasis,i)] += m_matrixM[i][l] * m_tensorQ[secondOrderGlobalIndex(m_dimRomBasis,l)];
2410 }
2411 }
2412 }
2413 ✗ for (int i=0; i<m_size2; i++) {
2414 ✗ m_tensorB[i] += dt * deltaB[i];
2415 ✗ m_tensorE[i] += dt * deltaE[i];
2416 ✗ m_tensorEs[i] += dt * deltaEs[i];
2417 }
2418 ✗ if (m_problem == 2) {
2419 ✗ for (int i=0; i<m_size2+m_dimRomBasis; i++)
2420 ✗ m_tensorQ[i] += dt * deltaQ[i];
2421 }
2422 ✗ for (int i=0; i<m_size3; i++)
2423 ✗ m_tensorT[i] += dt * deltaT[i];
2424 ✗ for (int i=0; i<m_dimRomBasis; i++) {
2425 ✗ for (int j=0; j<m_size2; j++)
2426 ✗ m_tensorU[i][j] += dt * deltaU[i][j];
2427 }
2428 ✗ for (int i=0; i<m_dimRomBasis; i++)
2429 ✗ delete [] deltaU[i];
2430 ✗ delete [] deltaU;
2431 ✗ break;
2432 }
2433 ✗ case 4: {
2434 ✗ double deltaB[m_size2];
2435 // B^{n+1} = B^n + dt * [B^n, M^n]^{(2)}
2436 ✗ double deltaE[m_size2];
2437 // E^{n+1} = E^n + dt * [E^n, M^n]^{(2)}
2438 ✗ double deltaQ[m_size2+m_dimRomBasis];
2439 // Q^{n+1} = Q^n + dt * [Q^n, M^n]^{(2)}
2440 ✗ double deltaT[m_size3];
2441 // T^{n+1} = T^n + dt * [T^n, M^n]^{(3)}
2442 ✗ double deltaTs[m_size3];
2443 // Ts^{n+1} = Ts^n + dt * [Ts^n, M^n]^{(3)}
2444 ✗ double deltaY[m_size4];
2445 // Y^{n+1} = Y^n + dt * [Y^n, M^n]^{(4)}
2446 ✗ for (int i=0; i<m_dimRomBasis; i++) {
2447 ✗ for (int j=0; j<i+1; j++) {
2448 ✗ deltaB[secondOrderGlobalIndex(i,j)] = 0.;
2449 ✗ if (m_problem == 2) { // bidomain
2450 ✗ deltaE[secondOrderGlobalIndex(i,j)] = 0.;
2451 ✗ deltaQ[secondOrderGlobalIndex(i,j)] = 0.;
2452 }
2453 ✗ for (int l=0; l<m_dimRomBasis; l++) {
2454 // delta B_ijk = sum_l (M_il * B_lj + M_jl * B_il )
2455 ✗ deltaB[secondOrderGlobalIndex(i,j)] += m_matrixM[i][l] * m_tensorB[secondOrderGlobalIndex(l,j)];
2456 ✗ deltaB[secondOrderGlobalIndex(i,j)] += m_matrixM[j][l] * m_tensorB[secondOrderGlobalIndex(i,l)];
2457 ✗ if (m_problem == 2) { // bidomain
2458 // delta E_ijk = sum_l (M_il * E_lj + M_jl * E_il )
2459 ✗ deltaE[secondOrderGlobalIndex(i,j)] += m_matrixM[i][l] * m_tensorE[secondOrderGlobalIndex(l,j)];
2460 ✗ deltaE[secondOrderGlobalIndex(i,j)] += m_matrixM[j][l] * m_tensorE[secondOrderGlobalIndex(i,l)];
2461 // delta Q_ijk = sum_l (M_il * Q_lj + M_jl * Q_il )
2462 ✗ deltaQ[secondOrderGlobalIndex(i,j)] += m_matrixM[i][l] * m_tensorQ[secondOrderGlobalIndex(l,j)];
2463 ✗ deltaQ[secondOrderGlobalIndex(i,j)] += m_matrixM[j][l] * m_tensorQ[secondOrderGlobalIndex(i,l)];
2464 }
2465 }
2466 ✗ for (int k=0; k<j+1; k++) {
2467 ✗ deltaT[thirdOrderGlobalIndex(i,j,k)] = 0.;
2468 ✗ deltaTs[thirdOrderGlobalIndex(i,j,k)] = 0.;
2469 ✗ for (int l=0; l<m_dimRomBasis; l++) {
2470 // delta T_ijk = sum_l (M_il * T_ljk + M_jl * T_ilk + M_kl * T_ijl )
2471 ✗ deltaT[thirdOrderGlobalIndex(i,j,k)] += m_matrixM[i][l] * m_tensorT[thirdOrderGlobalIndex(l,j,k)];
2472 ✗ deltaT[thirdOrderGlobalIndex(i,j,k)] += m_matrixM[j][l] * m_tensorT[thirdOrderGlobalIndex(i,l,k)];
2473 ✗ deltaT[thirdOrderGlobalIndex(i,j,k)] += m_matrixM[k][l] * m_tensorT[thirdOrderGlobalIndex(i,j,l)];
2474
2475 ✗ deltaTs[thirdOrderGlobalIndex(i,j,k)] += m_matrixM[i][l] * m_tensorTs[thirdOrderGlobalIndex(l,j,k)];
2476 ✗ deltaTs[thirdOrderGlobalIndex(i,j,k)] += m_matrixM[j][l] * m_tensorTs[thirdOrderGlobalIndex(i,l,k)];
2477 ✗ deltaTs[thirdOrderGlobalIndex(i,j,k)] += m_matrixM[k][l] * m_tensorTs[thirdOrderGlobalIndex(i,j,l)];
2478 }
2479 ✗ for (int h=0; h<k+1; h++) {
2480 ✗ deltaY[fourthOrderGlobalIndex(i,j,k,h)] = 0.;
2481 ✗ for (int l=0; l<m_dimRomBasis; l++) {
2482 // delta Y_ijkh = sum_l ( M_il * Y_ljkh + M_jl * Y_ilkh + M_kl * Y_ijlh + M_hl * Y_ijkl )
2483 ✗ deltaY[fourthOrderGlobalIndex(i,j,k,h)] += m_matrixM[i][l] * m_tensorY[fourthOrderGlobalIndex(l,j,k,h)];
2484 ✗ deltaY[fourthOrderGlobalIndex(i,j,k,h)] += m_matrixM[j][l] * m_tensorY[fourthOrderGlobalIndex(i,l,k,h)];
2485 ✗ deltaY[fourthOrderGlobalIndex(i,j,k,h)] += m_matrixM[k][l] * m_tensorY[fourthOrderGlobalIndex(i,j,l,h)];
2486 ✗ deltaY[fourthOrderGlobalIndex(i,j,k,h)] += m_matrixM[h][l] * m_tensorY[fourthOrderGlobalIndex(i,j,k,l)];
2487 }
2488 }
2489 }
2490 }
2491 ✗ if (m_problem == 2) {
2492 ✗ deltaQ[secondOrderGlobalIndex(m_dimRomBasis,i)] = 0.;
2493 ✗ for (int l=0; l<m_dimRomBasis; l++) {
2494 // delta Q_Nj = sum_l ( M_jl * Q_lj )
2495 ✗ deltaQ[secondOrderGlobalIndex(m_dimRomBasis,i)] += m_matrixM[i][l] * m_tensorQ[secondOrderGlobalIndex(m_dimRomBasis,l)];
2496 }
2497 }
2498 }
2499
2500 ✗ for (int id=0; id<m_size2; id++)
2501 ✗ m_tensorB[id] += dt * deltaB[id];
2502 ✗ if (m_problem == 2) {
2503 ✗ for (int id=0; id<m_size2; id++)
2504 ✗ m_tensorE[id] += dt * deltaE[id];
2505 ✗ for (int id=0; id<m_size2+m_dimRomBasis; id++)
2506 ✗ m_tensorQ[id] += dt * deltaQ[id];
2507 }
2508 ✗ for (int id=0; id<m_size3; id++) {
2509 ✗ m_tensorT[id] += dt * deltaT[id];
2510 ✗ m_tensorTs[id] += dt * deltaTs[id];
2511 }
2512 ✗ for (int id=0; id<m_size4; id++)
2513 ✗ m_tensorY[id] += dt * deltaY[id];
2514 ✗ break;
2515 }
2516 ✗ default:
2517 ✗ break;
2518 }
2519 ✗ break;
2520 }
2521 ✗ default:
2522 ✗ FEL_ERROR("Model not defined for ALP solver.");
2523 ✗ break;
2524 }
2525 }
2526
2527 ✗ void EigenProblemALP::updateBeta() {
2528 ✗ if (FelisceParam::verbose() > 10 ) PetscPrintf(PETSC_COMM_WORLD,"\nEigenProblemALP::updateBeta()\n");
2529
2530 ✗ double dt = FelisceParam::instance().timeStep;
2531 ✗ double deltaBeta[m_dimRomBasis];
2532 // beta^{n+1} = beta^n - dt * M^n * beta^n + dt * gamma^n
2533 ✗ for (int i=0; i<m_dimRomBasis; i++) {
2534 ✗ deltaBeta[i] = m_gamma[i];
2535
2536 ✗ for (int h=0; h<m_dimRomBasis; h++) {
2537 ✗ deltaBeta[i] -= m_beta[h] * m_matrixM[h][i];
2538 }
2539 }
2540 ✗ for (int i=0; i<m_dimRomBasis; i++) {
2541 ✗ m_beta[i] += dt * deltaBeta[i];
2542 }
2543
2544 ✗ if ( (m_problem == 1) || (m_problem == 2) ) {
2545 ✗ double deltaMu[m_dimRomBasis];
2546 // mu^{n+1} = mu^n - dt * M^n * mu^n + dt * eta^n
2547 ✗ for (int i=0; i<m_dimRomBasis; i++) {
2548 ✗ deltaMu[i] = m_eta[i];
2549
2550 ✗ for (int h=0; h<m_dimRomBasis; h++) {
2551 ✗ deltaMu[i] -= m_mu[h] * m_matrixM[h][i];
2552 }
2553 }
2554 ✗ for (int i=0; i<m_dimRomBasis; i++) {
2555 ✗ m_mu[i] += dt * deltaMu[i];
2556 }
2557 }
2558
2559 ✗ if (m_problem == 2) {
2560 ✗ double* newXi = new double[m_dimRomBasis+1];
2561 ✗ double* rhsXi = new double[m_dimRomBasis+1];
2562 ✗ for (int i=0; i<m_dimRomBasis; i++) {
2563 ✗ newXi[i] = m_xi[i];
2564 ✗ rhsXi[i] = 0.;
2565 ✗ for (int j=0; j<m_dimRomBasis; j++) {
2566 ✗ rhsXi[i] += - m_beta[j]*m_tensorE[secondOrderGlobalIndex(j,i)];
2567 }
2568 }
2569 ✗ newXi[m_dimRomBasis] = 0.;
2570 ✗ rhsXi[m_dimRomBasis] = 0.;
2571
2572 ✗ ConjGrad(m_tensorQ, rhsXi, newXi, 1.e-09, 1000, m_dimRomBasis+1);
2573
2574 ✗ for (int i=0; i<m_dimRomBasis; i++) {
2575 ✗ m_xi[i] = newXi[i];
2576 }
2577 ✗ delete [] newXi;
2578 ✗ delete [] rhsXi;
2579
2580 }
2581
2582 ✗ if (FelisceParam::instance().hasSource) {
2583 ✗ for (int idS=0; idS<m_numSource; idS++) {
2584 ✗ double* newSource = new double[m_dimRomBasis];
2585 // source^{n+1} = source^n - dt * M^n * source^n
2586 ✗ for (int i=0; i<m_dimRomBasis; i++) {
2587 ✗ newSource[i] = 0.;
2588
2589 ✗ for (int h=0; h<m_dimRomBasis; h++) {
2590 ✗ newSource[i] += m_source[idS][h] * m_matrixM[i][h];
2591 }
2592 }
2593 ✗ for (int i=0; i<m_dimRomBasis; i++) {
2594 ✗ m_source[idS][i] += dt * newSource[i];
2595 }
2596 ✗ delete [] newSource;
2597 }
2598 }
2599
2600 }
2601
2602 ✗ void EigenProblemALP::updateBetaMonolithic() {
2603 ✗ if (FelisceParam::verbose() > 10 )
2604 ✗ PetscPrintf(PETSC_COMM_WORLD,"\nEigenProblemALP::updateBetaMonolithic()\n");
2605
2606 ✗ double dt = FelisceParam::instance().timeStep;
2607 double value;
2608 double rhs1;
2609 double rhs2;
2610 felInt ia;
2611 ✗ m_rhsRom.set(0.);
2612
2613 // beta^{n+1} = beta^n + dt* gamma^n - dt * M * beta^n
2614 // mu^{n+1} = mu^n + dt * eta^n - dt * M * mu^n
2615 // E beta^{n+1} + Q xi^{n+1} = 0
2616 ✗ for (int i=0; i<m_dimRomBasis; i++) {
2617 ✗ rhs1 = m_beta[i] + dt * m_gamma[i];
2618 ✗ rhs2 = m_mu[i] + dt * m_eta[i];
2619
2620 ✗ for (int h=0; h<m_dimRomBasis; h++) {
2621 ✗ rhs1 -= dt * m_beta[h] * m_matrixM[h][i];
2622 ✗ rhs2 -= dt * m_mu[h] * m_matrixM[h][i];
2623 }
2624
2625 ✗ ia = i;
2626 ✗ m_rhsRom.setValue(ia,rhs1,INSERT_VALUES);
2627
2628 ✗ ia = i+m_dimRomBasis;
2629 ✗ m_rhsRom.setValue(ia,rhs2,INSERT_VALUES);
2630 }
2631 ✗ m_rhsRom.assembly();
2632
2633 // MatrixRom = [ I 0 0 ]
2634 // [ 0 I 0 ]
2635 // [ E 0 Q ]
2636
2637 felInt ii;
2638 felInt jj;
2639 ✗ for (int i=0; i<m_dimRomBasis; i++) {
2640 ✗ value = 1.;
2641
2642 ✗ ii = i;
2643 ✗ m_matrixRom.setValues(1,&ii,1,&ii,&value,INSERT_VALUES);
2644
2645 ✗ ii = i+m_dimRomBasis;
2646 ✗ m_matrixRom.setValues(1,&ii,1,&ii,&value,INSERT_VALUES);
2647
2648 ✗ for (int j=0; j<m_dimRomBasis; j++) {
2649 ✗ double valueMatrix = 0.;
2650
2651 ✗ valueMatrix = m_tensorE[secondOrderGlobalIndex(i,j)];
2652 ✗ ii = i+2*m_dimRomBasis;
2653 ✗ jj = j;
2654 ✗ m_matrixRom.setValues(1,&ii,1,&jj,&valueMatrix,INSERT_VALUES);
2655
2656 ✗ valueMatrix = m_tensorQ[secondOrderGlobalIndex(i,j)];
2657 ✗ jj = j+2*m_dimRomBasis;
2658 ✗ m_matrixRom.setValues(1,&ii,1,&jj,&valueMatrix,INSERT_VALUES);
2659 }
2660 }
2661 ✗ ii = 3*m_dimRomBasis;
2662 ✗ for (int j=0; j<m_dimRomBasis; j++) {
2663 ✗ jj = j+2*m_dimRomBasis;
2664 ✗ double valueMatrix = m_tensorQ[secondOrderGlobalIndex(j,m_dimRomBasis)];
2665 ✗ m_matrixRom.setValues(1,&ii,1,&jj,&valueMatrix,INSERT_VALUES);
2666 ✗ m_matrixRom.setValues(1,&jj,1,&ii,&valueMatrix,INSERT_VALUES);
2667 }
2668 ✗ jj = 3*m_dimRomBasis;
2669 ✗ double valueMatrix = 0.;
2670 ✗ m_matrixRom.setValues(1,&ii,1,&jj,&valueMatrix,INSERT_VALUES);
2671
2672 ✗ m_matrixRom.assembly(MAT_FINAL_ASSEMBLY);
2673
2674 ✗ m_solutionRom.set(0.);
2675
2676 // Solve
2677 ✗ m_kspRom->setKSPOperator(m_matrixRom, "");
2678 ✗ m_kspRom->solve(m_rhsRom, m_solutionRom, 0);
2679
2680 ✗ for (int i=0; i<m_dimRomBasis; i++) {
2681 ✗ ia = i;
2682 ✗ m_solutionRom.getValues(1,&ia,&value);
2683 ✗ m_beta[i] = value;
2684 ✗ ia = i+m_dimRomBasis;
2685 ✗ m_solutionRom.getValues(1,&ia,&value);
2686 ✗ m_mu[i] = value;
2687 ✗ ia = i+2*m_dimRomBasis;
2688 ✗ m_solutionRom.getValues(1,&ia,&value);
2689 ✗ m_xi[i] = value;
2690 }
2691
2692 }
2693
2694 ✗ void EigenProblemALP::updateBetaEI() {
2695 ✗ if (FelisceParam::verbose() > 10 )
2696 ✗ PetscPrintf(PETSC_COMM_WORLD,"\nEigenProblemALP::updateBetaEI()\n");
2697
2698 ✗ double dt = FelisceParam::instance().timeStep;
2699 double value;
2700 double rhs1;
2701 double rhs2;
2702 felInt ia;
2703 ✗ m_rhsRom.set(0.);
2704
2705 // (I + dt M) beta^{n+1} = beta^n + dt* gamma^n
2706 // (I + dt M) mu^{n+1} = mu^n + dt * eta^n
2707 // E beta^{n+1} + Q xi^{n+1} = 0
2708 ✗ for (int i=0; i<m_dimRomBasis; i++) {
2709 ✗ rhs1 = m_beta[i] + dt * m_gamma[i];
2710 ✗ rhs2 = m_mu[i] + dt * m_eta[i];
2711
2712 ✗ ia = i;
2713 ✗ m_rhsRom.setValue(ia,rhs1,INSERT_VALUES);
2714
2715 ✗ ia = i+m_dimRomBasis;
2716 ✗ m_rhsRom.setValue(ia,rhs2,INSERT_VALUES);
2717 }
2718 ✗ m_rhsRom.assembly();
2719
2720 // MatrixRom = [ I+dtM 0 0 ]
2721 // [ 0 I+dtM 0 ]
2722 // [ E 0 Q ]
2723
2724 felInt ii;
2725 felInt jj;
2726 ✗ for (int i=0; i<m_dimRomBasis; i++) {
2727 ✗ value = 1.;
2728
2729 ✗ ii = i;
2730 ✗ m_matrixRom.setValues(1,&ii,1,&ii,&value,INSERT_VALUES);
2731
2732 ✗ ii = i+m_dimRomBasis;
2733 ✗ m_matrixRom.setValues(1,&ii,1,&ii,&value,INSERT_VALUES);
2734
2735 ✗ for (int j=0; j<m_dimRomBasis; j++) {
2736 ✗ double valueMatrix = 0.;
2737
2738 ✗ valueMatrix = dt * m_matrixM[i][j];
2739 ✗ ii = i;
2740 ✗ jj = j;
2741 ✗ if (ii != jj) {
2742 ✗ m_matrixRom.setValues(1,&ii,1,&jj,&valueMatrix,INSERT_VALUES);
2743 }
2744 ✗ ii = i+m_dimRomBasis;
2745 ✗ jj = j+m_dimRomBasis;
2746 ✗ if (ii != jj) {
2747 ✗ m_matrixRom.setValues(1,&ii,1,&jj,&valueMatrix,INSERT_VALUES);
2748 }
2749
2750 ✗ valueMatrix = dt * m_tensorE[secondOrderGlobalIndex(i,j)];
2751 ✗ ii = i+2*m_dimRomBasis;
2752 ✗ jj = j;
2753 ✗ m_matrixRom.setValues(1,&ii,1,&jj,&valueMatrix,INSERT_VALUES);
2754
2755 ✗ valueMatrix = dt * m_tensorQ[secondOrderGlobalIndex(i,j)];
2756 ✗ jj = j+2*m_dimRomBasis;
2757 ✗ m_matrixRom.setValues(1,&ii,1,&jj,&valueMatrix,INSERT_VALUES);
2758 }
2759 }
2760 ✗ ii = 3*m_dimRomBasis;
2761 ✗ for (int j=0; j<m_dimRomBasis; j++) {
2762 ✗ jj = j+2*m_dimRomBasis;
2763 ✗ double valueMatrix = dt * m_tensorQ[secondOrderGlobalIndex(j,m_dimRomBasis)];
2764 ✗ m_matrixRom.setValues(1,&ii,1,&jj,&valueMatrix,INSERT_VALUES);
2765 ✗ m_matrixRom.setValues(1,&jj,1,&ii,&valueMatrix,INSERT_VALUES);
2766 }
2767 ✗ jj = 3*m_dimRomBasis;
2768 ✗ double valueMatrix = 0.;
2769 ✗ m_matrixRom.setValues(1,&ii,1,&jj,&valueMatrix,INSERT_VALUES);
2770
2771 ✗ m_matrixRom.assembly(MAT_FINAL_ASSEMBLY);
2772
2773 ✗ m_solutionRom.set(0.);
2774
2775 // Solve
2776 ✗ m_kspRom->setKSPOperator(m_matrixRom, "");
2777 ✗ m_kspRom->solve(m_rhsRom, m_solutionRom, 0);
2778
2779 ✗ for (int i=0; i<m_dimRomBasis; i++) {
2780 ✗ ia = i;
2781 ✗ m_solutionRom.getValues(1,&ia,&value);
2782 ✗ m_beta[i] = value;
2783 ✗ ia = i+m_dimRomBasis;
2784 ✗ m_solutionRom.getValues(1,&ia,&value);
2785 ✗ m_mu[i] = value;
2786 ✗ ia = i+2*m_dimRomBasis;
2787 ✗ m_solutionRom.getValues(1,&ia,&value);
2788 ✗ m_xi[i] = value;
2789 }
2790
2791 }
2792
2793 ✗ void EigenProblemALP::updateBetaBdf2() {
2794
2795 ✗ if (FelisceParam::verbose() > 10 )
2796 ✗ PetscPrintf(PETSC_COMM_WORLD,"\nEigenProblemALP::updateBetaEI()\n");
2797
2798 ✗ double dt = FelisceParam::instance().timeStep;
2799 double value;
2800 double rhs1;
2801 double rhs2;
2802 double coefTimeDer[3];
2803 ✗ coefTimeDer[0] = 3./2.;
2804 ✗ coefTimeDer[1] = 2.;
2805 ✗ coefTimeDer[2] = -1./2.;
2806 felInt ia;
2807 ✗ m_rhsRom.set(0.);
2808
2809 // (3/2 I + dt M) beta^{n+1} = 2 beta^n - 1/2 beta^{n-1}+ dt* gamma(beta_extrap = 2 beta^n - beta^{n-1})
2810 // (3/2 I + dt M) mu^{n+1} = 2 mu^n - 1/2 mu^{n-1} + dt * eta(mu_extrap = 2 mu^n - mu^{n-1))
2811 // E beta^{n+1} + Q xi^{n+1} = 0
2812 ✗ for (int i=0; i<m_dimRomBasis; i++) {
2813 ✗ rhs1 = coefTimeDer[1]*m_beta[i] + coefTimeDer[2]*m_beta_n_1[i] + dt * m_gamma_extrap[i];
2814 ✗ rhs2 = coefTimeDer[1]*m_mu[i] + coefTimeDer[2]*m_mu_n_1[i] + dt * m_eta_extrap[i];
2815
2816 ✗ ia = i;
2817 ✗ m_rhsRom.setValue(ia,rhs1,INSERT_VALUES);
2818
2819 ✗ ia = i+m_dimRomBasis;
2820 ✗ m_rhsRom.setValue(ia,rhs2,INSERT_VALUES);
2821 }
2822 ✗ m_rhsRom.assembly();
2823
2824 // MatrixRom = [ 3/2I+dtM 0 0 ]
2825 // [ 0 3/2I+dtM 0 ]
2826 // [ E 0 Q ]
2827
2828 felInt ii;
2829 felInt jj;
2830 ✗ for (int i=0; i<m_dimRomBasis; i++) {
2831 ✗ value = coefTimeDer[0];
2832
2833 ✗ ii = i;
2834 ✗ m_matrixRom.setValues(1,&ii,1,&ii,&value,INSERT_VALUES);
2835
2836 ✗ ii = i+m_dimRomBasis;
2837 ✗ m_matrixRom.setValues(1,&ii,1,&ii,&value,INSERT_VALUES);
2838
2839 ✗ for (int j=0; j<m_dimRomBasis; j++) {
2840 ✗ double valueMatrix = 0.;
2841
2842 ✗ valueMatrix = dt * m_matrixM[i][j];
2843 ✗ ii = i;
2844 ✗ jj = j;
2845 ✗ if (ii != jj) {
2846 ✗ m_matrixRom.setValues(1,&ii,1,&jj,&valueMatrix,INSERT_VALUES);
2847 }
2848 ✗ ii = i+m_dimRomBasis;
2849 ✗ jj = j+m_dimRomBasis;
2850 ✗ if (ii != jj) {
2851 ✗ m_matrixRom.setValues(1,&ii,1,&jj,&valueMatrix,INSERT_VALUES);
2852 }
2853
2854 ✗ valueMatrix = dt * m_tensorE[secondOrderGlobalIndex(i,j)];
2855 ✗ ii = i+2*m_dimRomBasis;
2856 ✗ jj = j;
2857 ✗ m_matrixRom.setValues(1,&ii,1,&jj,&valueMatrix,INSERT_VALUES);
2858
2859 ✗ valueMatrix = dt * m_tensorQ[secondOrderGlobalIndex(i,j)];
2860 ✗ jj = j+2*m_dimRomBasis;
2861 ✗ m_matrixRom.setValues(1,&ii,1,&jj,&valueMatrix,INSERT_VALUES);
2862 }
2863 }
2864 ✗ ii = 3*m_dimRomBasis;
2865 ✗ for (int j=0; j<m_dimRomBasis; j++) {
2866 ✗ jj = j+2*m_dimRomBasis;
2867 ✗ double valueMatrix = dt * m_tensorQ[secondOrderGlobalIndex(j,m_dimRomBasis)];
2868 ✗ m_matrixRom.setValues(1,&ii,1,&jj,&valueMatrix,INSERT_VALUES);
2869 ✗ m_matrixRom.setValues(1,&jj,1,&ii,&valueMatrix,INSERT_VALUES);
2870 }
2871 ✗ jj = 3*m_dimRomBasis;
2872 ✗ double valueMatrix = 0.;
2873 ✗ m_matrixRom.setValues(1,&ii,1,&jj,&valueMatrix,INSERT_VALUES);
2874
2875 ✗ m_matrixRom.assembly(MAT_FINAL_ASSEMBLY);
2876
2877 ✗ m_solutionRom.set(0.);
2878
2879 // Solve
2880 ✗ m_kspRom->setKSPOperator(m_matrixRom, "");
2881 ✗ m_kspRom->solve(m_rhsRom, m_solutionRom, 0);
2882
2883 ✗ for (int i=0; i<m_dimRomBasis; i++) {
2884 ✗ m_beta_n_1[i] = m_beta[i];
2885 ✗ m_mu_n_1[i] = m_mu[i];
2886 ✗ m_xi_n_1[i] = m_xi[i];
2887 }
2888
2889 ✗ for (int i=0; i<m_dimRomBasis; i++) {
2890 ✗ ia = i;
2891 ✗ m_solutionRom.getValues(1,&ia,&value);
2892 ✗ m_beta[i] = value;
2893 ✗ ia = i+m_dimRomBasis;
2894 ✗ m_solutionRom.getValues(1,&ia,&value);
2895 ✗ m_mu[i] = value;
2896 ✗ ia = i+2*m_dimRomBasis;
2897 ✗ m_solutionRom.getValues(1,&ia,&value);
2898 ✗ m_xi[i] = value;
2899 }
2900
2901 }
2902
2903 ✗ void EigenProblemALP::ConjGrad(double* A, double* b, double* x, double tol, int maxIter, int n) {
2904 ✗ double* res = new double[n];
2905 ✗ double* p = new double[n];
2906 ✗ double* Ax = new double[n];
2907
2908 ✗ for (int i=0; i<n; i++) {
2909 ✗ Ax[i] = 0.;
2910 ✗ for (int j=0; j<n; j++) {
2911 ✗ Ax[i] += A[secondOrderGlobalIndex(i,j)]*x[j];
2912 }
2913 }
2914 ✗ for (int i=0; i<n; i++) {
2915 ✗ res[i] = b[i] - Ax[i];
2916 ✗ p[i] = res[i];
2917 }
2918 ✗ double res2 = 0.;
2919 ✗ double newRes2 = 0.;
2920 ✗ for (int i=0; i<n; i++) {
2921 ✗ res2 += res[i]*res[i];
2922 }
2923
2924 ✗ double* Ap = new double[n];
2925 double pAp;
2926 double alpha;
2927 double beta;
2928 ✗ bool isConv = false;
2929 ✗ int iter = 1;
2930
2931 ✗ if (res2 < tol) {
2932 ✗ isConv = true;
2933 ✗ if (FelisceParam::verbose()) {
2934 ✗ PetscPrintf(PETSC_COMM_WORLD,"Conjugate Gradient converged in %i iterations, residual: %1.10f \n",0,res2);
2935 }
2936 }
2937
2938 ✗ while (!isConv) {
2939 ✗ for (int i=0; i<n; i++) {
2940 ✗ Ap[i] = 0.;
2941 ✗ for (int j=0; j<n; j++) {
2942 ✗ Ap[i] += A[secondOrderGlobalIndex(i,j)]*p[j];
2943 }
2944 }
2945 ✗ pAp = 0.;
2946 ✗ for (int i=0; i<n; i++)
2947 ✗ pAp += p[i]*Ap[i];
2948 ✗ alpha = res2/pAp;
2949
2950 ✗ for (int i=0; i<n; i++) {
2951 ✗ x[i] += alpha*p[i];
2952 ✗ res[i] -= alpha*Ap[i];
2953 ✗ newRes2 += res[i]*res[i];
2954 }
2955
2956 ✗ if ( (newRes2 < tol) || (iter >= maxIter) ) {
2957 ✗ if (newRes2 < tol) {
2958 ✗ isConv = true;
2959 ✗ if (FelisceParam::verbose()) {
2960 ✗ PetscPrintf(PETSC_COMM_WORLD,"Conjugate Gradient converged in %i iterations, residual: %1.10f \n",iter,newRes2);
2961 }
2962 } else {
2963 ✗ PetscPrintf(PETSC_COMM_WORLD,"Conjugate Gradient: %i iterations, residual: %1.10f \n",iter,newRes2);
2964 ✗ FEL_ERROR("Conjugate Gradient did not converge.");
2965 }
2966 ✗ break;
2967 } else {
2968 ✗ beta = newRes2/res2;
2969 ✗ res2 = newRes2;
2970 ✗ newRes2 = 0.;
2971 ✗ for (int i=0; i<n; i++) {
2972 ✗ p[i] = res[i] + beta*p[i];
2973 }
2974 ✗ iter++;
2975 }
2976 }
2977
2978 ✗ delete [] res;
2979 ✗ delete [] p;
2980 ✗ delete [] Ap;
2981 ✗ delete [] Ax;
2982
2983 }
2984
2985
2986
2987 ✗ void EigenProblemALP::solveOrthComp(double* A, double* b, double* x, double tol, int maxIter) {
2988 ✗ for (int i=0; i<m_dimOrthComp; i++) {
2989 ✗ x[i] = 0.0;
2990 }
2991
2992 ✗ double* res = new double[m_dimOrthComp];
2993 ✗ double* p = new double[m_dimOrthComp];
2994 ✗ double* Ax = new double[m_dimOrthComp];
2995
2996 ✗ for (int i=0; i<m_dimOrthComp; i++) {
2997 ✗ Ax[i] = 0.;
2998 ✗ for (int j=0; j<m_dimOrthComp; j++) {
2999 ✗ Ax[i] += A[i*(m_dimOrthComp)+j]*x[j];
3000 }
3001 }
3002 ✗ for (int i=0; i<m_dimOrthComp; i++) {
3003 ✗ res[i] = b[i] - Ax[i];
3004 ✗ p[i] = res[i];
3005 }
3006 ✗ double res2 = 0.;
3007 ✗ double newRes2 = 0.;
3008 ✗ for (int i=0; i<m_dimOrthComp; i++) {
3009 ✗ res2 += res[i]*res[i];
3010 }
3011
3012 ✗ double* Ap = new double[m_dimOrthComp];
3013 double pAp;
3014 double alpha;
3015 double beta;
3016 ✗ bool isConv = false;
3017 ✗ int iter = 1;
3018
3019 ✗ while (!isConv) {
3020 ✗ for (int i=0; i<m_dimOrthComp; i++) {
3021 ✗ Ap[i] = 0.;
3022 ✗ for (int j=0; j<m_dimOrthComp; j++) {
3023 ✗ Ap[i] += A[i*(m_dimOrthComp)+j]*p[j];
3024 }
3025 }
3026 ✗ pAp = 0.;
3027 ✗ for (int i=0; i<m_dimOrthComp; i++)
3028 ✗ pAp += p[i]*Ap[i];
3029 ✗ alpha = res2/pAp;
3030
3031 ✗ for (int i=0; i<m_dimOrthComp; i++) {
3032 ✗ x[i] += alpha*p[i];
3033 ✗ res[i] -= alpha*Ap[i];
3034 ✗ newRes2 += res[i]*res[i];
3035 }
3036
3037 ✗ if ( (newRes2 < tol) || (iter >= maxIter) ) {
3038 ✗ isConv = true;
3039 ✗ if (FelisceParam::verbose() > 10) {
3040 ✗ PetscPrintf(PETSC_COMM_WORLD,"Conjugate Gradient converged in %i iterations, residual: %1.10f \n",iter,newRes2);
3041 }
3042
3043 ✗ break;
3044 } else {
3045 ✗ beta = newRes2/res2;
3046 ✗ res2 = newRes2;
3047 ✗ newRes2 = 0.;
3048 ✗ for (int i=0; i<m_dimOrthComp; i++) {
3049 ✗ p[i] = res[i] + beta*p[i];
3050 }
3051 ✗ iter++;
3052 }
3053 }
3054
3055 ✗ for (int i=0; i<m_dimOrthComp; i++) {
3056 ✗ Ax[i] = 0.;
3057 ✗ for (int j=0; j<m_dimOrthComp; j++) {
3058 ✗ Ax[i] += A[i*(m_dimOrthComp)+j]*x[j];
3059 }
3060 }
3061
3062 ✗ delete [] res;
3063 ✗ delete [] p;
3064 ✗ delete [] Ap;
3065 ✗ delete [] Ax;
3066
3067 }
3068
3069
3070
3071 ✗ void EigenProblemALP::updateEigenvalue() {
3072 ✗ if (FelisceParam::verbose() > 10 ) PetscPrintf(PETSC_COMM_WORLD,"\nEigenProblemALP::updateEigenvalue()\n");
3073
3074 ✗ double dt = FelisceParam::instance().timeStep;
3075 ✗ double newLambda[m_dimRomBasis];
3076 // lambda_i^{n+1} = lambda_i^n - dt * m_coefChi * sum_h ( gamma_h * T_hii )
3077 ✗ for (int i=0; i<m_dimRomBasis; i++) {
3078 ✗ newLambda[i] = 0.;
3079 ✗ for (int h=0; h<m_dimRomBasis; h++) {
3080 ✗ newLambda[i] += m_gamma[h] * m_tensorT[thirdOrderGlobalIndex(h,i,i)];
3081 }
3082 }
3083 ✗ for (int i=0; i<m_dimRomBasis; i++) {
3084 ✗ m_eigenValue[i] -= dt * m_coefChi * newLambda[i];
3085 }
3086 }
3087
3088 // Projection functions from FE space to RO one
3089 3 void EigenProblemALP::projectOnBasis(PetscVector& vectorDof, double* vectorBasis, bool flagMass, int size)
3090 {
3091
1/2
✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
 3 if (size == 0)
3092 3 size = m_dimRomBasis;
3093
3094
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 PetscVector tmpVec;
3095
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 tmpVec.duplicateFrom(vectorDof);
3096
3097
1/2
✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
 3 if (flagMass) {
3098
2/2
✓ Branch 0 taken 75 times.
✓ Branch 1 taken 3 times.
 78 for (int iBasis=0; iBasis<size; iBasis++) {
3099
1/2
✓ Branch 3 taken 75 times.
✗ Branch 4 not taken.
 75 mult(m_Matrix[m_idG], m_basis[iBasis],tmpVec);
3100
3101
1/2
✓ Branch 1 taken 75 times.
✗ Branch 2 not taken.
 75 dot(tmpVec,vectorDof,&vectorBasis[iBasis]);
3102
3103 }
3104 } else {
3105 ✗ for (int iBasis=0; iBasis<size; iBasis++) {
3106 ✗ dot( m_basis[iBasis],vectorDof,&vectorBasis[iBasis]);
3107 }
3108 }
3109
3110
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 tmpVec.destroy();
3111 3 }
3112
3113 // Projection functions from RO space to FE one
3114 75 void EigenProblemALP::projectOnDof(double* vectorBasis, PetscVector& vectorDof, int size) {
3115
1/2
✓ Branch 1 taken 75 times.
✗ Branch 2 not taken.
 75 vectorDof.set(0.);
3116
2/2
✓ Branch 0 taken 975 times.
✓ Branch 1 taken 75 times.
 1050 for (int iBasis=0; iBasis<size; iBasis++) {
3117 975 vectorDof.axpy(vectorBasis[iBasis],m_basis[iBasis]);
3118 }
3119 75 }
3120
3121 // Auxiliary function to debug purposes.
3122 ✗ void EigenProblemALP::viewALP(double* vToPrint, int vSize, std::ofstream outStream) {
3123 int rankProc;
3124 ✗ MPI_Comm_rank(m_petscComm,&rankProc);
3125 ✗ if (rankProc == 0) {
3126 ✗ for(int i=0; i<vSize; i++) {
3127 ✗ outStream << vToPrint[i] << std::endl;
3128 }
3129 ✗ outStream.close();
3130 }
3131 }
3132
3133
3134 4 void EigenProblemALP::viewALP(double* vToPrint, int vSize, std::string fName) {
3135 int rankProc;
3136
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 MPI_Comm_rank(m_petscComm,&rankProc);
3137
1/2
✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
 4 if (rankProc == 0) {
3138
1/2
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
 4 std::ofstream outStream(fName.c_str());
3139
2/2
✓ Branch 0 taken 100 times.
✓ Branch 1 taken 4 times.
 104 for(int i=0; i<vSize; i++) {
3140
2/4
✓ Branch 1 taken 100 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 100 times.
✗ Branch 5 not taken.
 100 outStream << vToPrint[i] << std::endl;
3141 }
3142
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 outStream.close();
3143 4 }
3144 4 }
3145
3146 }
3147