GCC Code Coverage Report

Directory: ./
File: Solver/eigenProblemALPDeim.cpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 0 1597 0.0%
Branches: 0 2468 0.0%
Line Branch Exec Source
1 // ______ ______ _ _ _____ ______
2 // | ____| ____| | (_)/ ____| | ____|
3 // | |__ | |__ | | _| (___ ___| |__
4 // | __| | __| | | | |\___ \ / __| __|
5 // | | | |____| |____| |____) | (__| |____
6 // |_| |______|______|_|_____/ \___|______|
7 // Finite Elements for Life Sciences and Engineering
8 //
9 // License: LGL2.1 License
10 // FELiScE default license: LICENSE in root folder
11 //
12 // Main authors: E. Schenone
13 //
14
15 // System includes
16
17 // External includes
18 // Project includes
19 #include "Solver/eigenProblemALPDeim.hpp"
20 #include "Core/felisceTools.hpp"
21
22 namespace felisce {
23 /*! Construtor
24 */
25 ✗ EigenProblemALPDeim::EigenProblemALPDeim():
26 EigenProblemALP(),
27 ✗ m_resPoint(false),
28 ✗ m_fePotExtraCell(nullptr),
29 ✗ m_dimCollocationPts(0),
30 ✗ m_collocationNode(nullptr),
31 ✗ m_basisDeimColloc(nullptr),
32 ✗ m_basisDeim(nullptr),
33 ✗ m_orthCompDeim(nullptr),
34 ✗ m_improvedRecType(0),
35 ✗ m_hatU(nullptr),
36 ✗ m_hatW(nullptr),
37 ✗ m_hatUe(nullptr),
38 ✗ m_potential(nullptr),
39 ✗ m_hatBasis(nullptr),
40 ✗ m_theta(nullptr),
41 ✗ m_thetaOrth(nullptr),
42 ✗ m_massDeim(nullptr),
43 ✗ m_stiffDeim(nullptr),
44 ✗ m_invGK(nullptr),
45 ✗ m_resU(nullptr),
46 ✗ m_hatF(nullptr),
47 ✗ m_hatG(nullptr),
48 ✗ m_matrixA(nullptr),
49 ✗ m_matrixB(nullptr),
50 ✗ m_sumG(nullptr),
51 ✗ m_coefFhNs(nullptr),
52 ✗ m_coefTauClose(nullptr),
53 ✗ m_coefTauOut(nullptr)
54 {}
55
56 ✗ EigenProblemALPDeim::~EigenProblemALPDeim() {
57 ✗ m_fePotExtraCell = nullptr;
58 ✗ if (m_collocationNode) {
59 ✗ delete [] m_collocationNode;
60 }
61 ✗ if (m_hatU) {
62 ✗ delete[] m_hatU;
63 }
64 ✗ if (m_hatW) {
65 ✗ delete[] m_hatW;
66 }
67 ✗ if (m_hatUe) {
68 ✗ delete[] m_hatUe;
69 }
70 ✗ if (m_potential) {
71 ✗ delete [] m_potential;
72 }
73 ✗ if (m_hatF) {
74 ✗ delete[] m_hatF;
75 }
76 ✗ if (m_hatG) {
77 ✗ delete[] m_hatG;
78 }
79 ✗ if (m_matrixA) {
80 ✗ for (int i=0; i< m_dimRomBasis+1; i++) {
81 ✗ delete [] m_matrixA[i];
82 }
83 ✗ delete [] m_matrixA;
84 }
85 ✗ if (m_matrixB) {
86 ✗ for (int i=0; i< m_dimRomBasis; i++) {
87 ✗ delete [] m_matrixB[i];
88 }
89 ✗ delete [] m_matrixB;
90 }
91 ✗ if (m_sumG) {
92 ✗ delete [] m_sumG;
93 }
94 ✗ if (m_basisDeim) {
95 ✗ for (int i=0; i< m_dimCollocationPts; i++) {
96 ✗ delete [] m_basisDeim[i];
97 }
98 ✗ delete [] m_basisDeim;
99 }
100 ✗ if (m_basisDeimColloc) {
101 ✗ for (int i=0; i< m_dimCollocationPts; i++) {
102 ✗ delete [] m_basisDeimColloc[i];
103 }
104 ✗ delete [] m_basisDeimColloc;
105 }
106 ✗ if (m_orthCompDeim) {
107 ✗ for (int i=0; i< m_dimOrthComp; i++) {
108 ✗ delete [] m_orthCompDeim[i];
109 }
110 ✗ delete [] m_orthCompDeim;
111 }
112 ✗ if (m_hatBasis) {
113 ✗ for (int i=0; i< m_dimCollocationPts; i++) {
114 ✗ delete [] m_hatBasis[i];
115 }
116 ✗ delete[] m_hatBasis;
117 }
118 ✗ if (m_theta) {
119 ✗ for (int i=0; i< m_dimRomBasis; i++) {
120 ✗ delete [] m_theta[i];
121 }
122 ✗ delete[] m_theta;
123 }
124 ✗ if (m_thetaOrth) {
125 ✗ for (int i=0; i< m_dimOrthComp; i++) {
126 ✗ delete [] m_thetaOrth[i];
127 }
128 ✗ delete[] m_thetaOrth;
129 }
130 ✗ if (m_massDeim) {
131 ✗ for (int i=0; i< m_dimCollocationPts; i++) {
132 ✗ delete [] m_massDeim[i];
133 }
134 ✗ delete[] m_massDeim;
135 }
136 ✗ if (m_stiffDeim) {
137 ✗ for (int i=0; i< m_dimCollocationPts; i++) {
138 ✗ delete [] m_stiffDeim[i];
139 }
140 ✗ delete[] m_stiffDeim;
141 }
142
143 ✗ if (m_invGK) {
144 ✗ for (int i=0; i< m_dimCollocationPts; i++) {
145 ✗ delete [] m_invGK[i];
146 }
147 ✗ delete [] m_invGK;
148 }
149
150 ✗ if (m_resU) {
151 ✗ delete [] m_resU;
152 }
153
154 ✗ if (FelisceParam::instance().ROMmethod == "ALP-DEIM") {
155 ✗ for (int i=0; i<m_dimCollocationPts; i++) {
156 ✗ m_basis[i].destroy();
157 }
158 }
159
160 ✗ if (m_projectorPiV.size() > 0) {
161 ✗ for (int i=0; i<m_dimCollocationPts; i++) {
162 ✗ m_projectorPiV[i].destroy();
163 }
164 }
165
166 ✗ if (m_coefFhNs) {
167 ✗ delete [] m_coefFhNs;
168 }
169 ✗ if (m_coefTauClose) {
170 ✗ delete [] m_coefTauClose;
171 }
172 ✗ if (m_coefTauOut) {
173 ✗ delete [] m_coefTauOut;
174 }
175 }
176
177 ✗ void EigenProblemALPDeim::initialize(const GeometricMeshRegion::Pointer& mesh, FelisceTransient::Pointer fstransient, MPI_Comm& comm) {
178 ✗ if (FelisceParam::verbose() > 10 ) PetscPrintf(PETSC_COMM_WORLD,"\nEigenProblemALPDeim::initialize\n");
179
180 ✗ EigenProblem::initialize(mesh, fstransient, comm);
181
182 ✗ std::unordered_map<std::string, int> mapOfPbm;
183 ✗ mapOfPbm["FKPP"] = 0;
184 ✗ mapOfPbm["Monodomain"] = 1;
185 ✗ mapOfPbm["Bidomain"] = 2;
186 ✗ mapOfPbm["BidomainCurv"] = 2;
187 ✗ mapOfPbm["BidomainCoupled"] = 3;
188 ✗ m_problem = mapOfPbm[FelisceParam::instance().model];
189
190 ✗ m_dimRomBasis = FelisceParam::instance().dimRomBasis;
191 ✗ m_dimCollocationPts = FelisceParam::instance().dimCollocationPts;
192 ✗ m_useImprovedRec = FelisceParam::instance().useImprovedRec;
193 ✗ m_dimOrthComp = FelisceParam::instance().dimOrthComp;
194
195 ✗ m_numSource = FelisceParam::instance().numberOfSource;
196
197 ✗ int id=0;
198 ✗ m_idL = id;
199 ✗ id++;
200 ✗ m_idG = id;
201 ✗ id++;
202 ✗ if (!FelisceParam::instance().solveEigenProblem) {
203 ✗ if (m_useImprovedRec) {
204 ✗ m_idK = id;
205 ✗ id++;
206 }
207 }
208 ✗ m_numberOfMatrix = id;
209
210 ✗ m_coefChi = FelisceParam::instance().chiSchrodinger;
211 ✗ std::vector<PhysicalVariable> listVariable;
212 ✗ std::vector<std::size_t> listNumComp;
213 ✗ switch (m_problem) {
214 ✗ case 0: {
215 ✗ nameOfTheProblem = "Problem Fisher-Kolmogorov-Petrovski-Piskunov equation (FKPP)";
216 ✗ listVariable.push_back(potTransMemb);
217 ✗ listNumComp.push_back(1);
218 //define unknown of the linear system.
219 ✗ m_listUnknown.push_back(potTransMemb);
220 ✗ break;
221 }
222 ✗ case 1: {
223 ✗ nameOfTheProblem = "Problem cardiac Monodomain";
224 ✗ listVariable.push_back(potTransMemb);
225 ✗ listNumComp.push_back(1);
226 //define unknown of the linear system.
227 ✗ m_listUnknown.push_back(potTransMemb);
228 ✗ break;
229 }
230 ✗ case 2: {
231 ✗ nameOfTheProblem = "Problem cardiac Bidomain";
232 ✗ listVariable.push_back(potTransMemb);
233 ✗ listNumComp.push_back(1);
234 //define unknown of the linear system.
235 ✗ m_listUnknown.push_back(potTransMemb);
236 ✗ break;
237 }
238 ✗ case 3: {
239 ✗ nameOfTheProblem = "Problem cardiac Bidomain";
240 ✗ listVariable.push_back(potTransMemb);
241 ✗ listNumComp.push_back(1);
242 ✗ listVariable.push_back(potExtraCell);
243 ✗ listNumComp.push_back(1);
244 //define unknown of the linear system.
245 ✗ m_listUnknown.push_back(potTransMemb);
246 ✗ m_listUnknown.push_back(potExtraCell);
247 ✗ break;
248 }
249 ✗ default:
250 ✗ FEL_ERROR("Model not defined for ALP solver.");
251 ✗ break;
252 }
253 ✗ definePhysicalVariable(listVariable,listNumComp);
254
255 }
256
257 // Define Physical Variable associate to the problem
258 ✗ void EigenProblemALPDeim::initPerElementType() {
259 ✗ switch(m_problem) {
260 ✗ case 0: {
261 //Init pointer on Finite Element, Variable or idVariable
262 ✗ m_ipotTransMemb = m_listVariable.getVariableIdList(potTransMemb);
263 ✗ m_fePotTransMemb = m_listCurrentFiniteElement[m_ipotTransMemb];
264 ✗ m_elemFieldU0.initialize(DOF_FIELD,*m_fePotTransMemb);
265 ✗ break;
266 }
267 ✗ case 1: {
268 //Init pointer on Finite Element, Variable or idVariable
269 ✗ m_ipotTransMemb = m_listVariable.getVariableIdList(potTransMemb);
270 ✗ m_fePotTransMemb = m_listCurrentFiniteElement[m_ipotTransMemb];
271 ✗ m_elemFieldU0.initialize(DOF_FIELD,*m_fePotTransMemb);
272 ✗ break;
273 }
274 ✗ case 2: {
275 //Init pointer on Finite Element, Variable or idVariable
276 ✗ m_ipotTransMemb = m_listVariable.getVariableIdList(potTransMemb);
277 ✗ m_fePotTransMemb = m_listCurrentFiniteElement[m_ipotTransMemb];
278 ✗ m_elemFieldU0.initialize(DOF_FIELD,*m_fePotTransMemb);
279 ✗ break;
280 }
281 ✗ case 3: {
282 //Init pointer on Finite Element, Variable or idVariable
283 ✗ m_ipotTransMemb = m_listVariable.getVariableIdList(potTransMemb);
284 ✗ m_ipotExtraCell = m_listVariable.getVariableIdList(potExtraCell);
285 ✗ m_fePotTransMemb = m_listCurrentFiniteElement[m_ipotTransMemb];
286 ✗ m_fePotExtraCell = m_fePotTransMemb;
287 ✗ m_elemFieldU0.initialize(DOF_FIELD,*m_fePotTransMemb);
288 ✗ m_elemFieldUe0.initialize(DOF_FIELD,*m_fePotExtraCell);
289 ✗ break;
290 }
291 ✗ default:
292 ✗ FEL_ERROR("Model not defined for ALP solver.");
293 ✗ break;
294 }
295
296 ✗ if (FelisceParam::instance().hasInfarct && FelisceParam::instance().typeOfIonicModel == "fhn") {
297 ✗ m_elemFieldFhNf0.initialize(DOF_FIELD,*m_fePotTransMemb);
298 }
299 }
300
301 // Assemble Matrix
302 ✗ void EigenProblemALPDeim::computeElementArray(const std::vector<Point*>& elemPoint, const std::vector<felInt>& elemIdPoint, felInt& iel, FlagMatrixRHS flagMatrixRHS) {
303 IGNORE_UNUSED_ELEM_ID_POINT;
304 IGNORE_UNUSED_FLAG_MATRIX_RHS;
305 ✗ if (FelisceParam::verbose() > 50 ) PetscPrintf(PETSC_COMM_WORLD,"\nEigenProblemALPDeim::computeElementArray\n");
306 ✗ switch(m_problem) {
307 ✗ case 0: { // FKPP
308 ✗ m_fePotTransMemb->updateMeasQuadPt(0, elemPoint);
309 ✗ m_fePotTransMemb->updateFirstDeriv(0, elemPoint);
310
311 // m_Matrix[0] = grad(phi_i) * grad(phi_j)
312 ✗ m_elementMat[m_idL]->grad_phi_i_grad_phi_j(1.,*m_fePotTransMemb,0,0,1);
313
314 ✗ if (FelisceParam::instance().hasInitialCondition) {
315 // m_Matrix[0] += - m_coefChi * V * phi_i * phi_j
316 ✗ m_elemFieldU0.setValue(m_U_0_seq, *m_fePotTransMemb, iel, m_ipotTransMemb, m_ao, m_dof);
317 ✗ m_elementMat[m_idL]->a_phi_i_phi_j(- m_coefChi,m_elemFieldU0,*m_fePotTransMemb,0,0,1);
318 }
319
320 // Matrix[1] = phi_i * phi_j
321 ✗ m_elementMat[m_idG]->phi_i_phi_j(1.,*m_fePotTransMemb,0,0,1);
322
323 ✗ if (!FelisceParam::instance().solveEigenProblem) {
324 ✗ if (m_useImprovedRec) {
325 // m_Matrix[m_numberOfMatrix-1] = grad(phi_i) * grad(phi_j)
326 ✗ m_elementMat[m_idK]->grad_phi_i_grad_phi_j(1.,*m_fePotTransMemb,0,0,1);
327 }
328 }
329 ✗ break;
330 }
331 ✗ case 1:
332 case 2: { // Monodomain
333 ✗ m_fePotTransMemb->updateMeasQuadPt(0, elemPoint);
334 ✗ m_fePotTransMemb->updateFirstDeriv(0, elemPoint);
335
336 // m_Matrix[0] = grad(phi_i) * grad(phi_j)
337 ✗ m_elementMat[m_idL]->grad_phi_i_grad_phi_j(FelisceParam::instance().intraTransvTensor,*m_fePotTransMemb,0,0,1);
338
339 ✗ if (FelisceParam::instance().testCase == 1) {
340 ✗ std::vector<double> elemFiber;
341 ✗ getFiberDirection(iel,m_ipotTransMemb, elemFiber);
342 // m_Matrix[0] += (\sigma_i^l-\sigma_i^t) a vec a grad(phi_i) * grad(phi_j)
343 ✗ m_elementMat[m_idL]->tau_grad_phi_i_tau_grad_phi_j(FelisceParam::instance().intraFiberTensor-FelisceParam::instance().intraTransvTensor,elemFiber,*m_fePotTransMemb,0,0,1);
344 }
345
346 ✗ if (FelisceParam::instance().hasInitialCondition) {
347 // m_Matrix[0] += - m_coefChi * V * phi_i * phi_j
348 ✗ m_elemFieldU0.setValue(m_U_0_seq, *m_fePotTransMemb, iel, m_ipotTransMemb, m_ao, m_dof);
349 ✗ m_elementMat[m_idL]->a_phi_i_phi_j(- m_coefChi,m_elemFieldU0,*m_fePotTransMemb,0,0,1);
350 }
351
352 // m_Matrix[1] = phi_i * phi_j
353 ✗ m_elementMat[m_idG]->phi_i_phi_j(1.,*m_fePotTransMemb,0,0,1);
354
355 ✗ if (!FelisceParam::instance().solveEigenProblem) {
356 ✗ if (m_useImprovedRec) {
357 // m_Matrix[2] = E_ij = grad(phi_i) * grad(phi_j)
358 ✗ m_elementMat[m_idK]->grad_phi_i_grad_phi_j(FelisceParam::instance().intraTransvTensor,*m_fePotTransMemb,0,0,1);
359 ✗ if (FelisceParam::instance().testCase == 1) {
360 ✗ std::vector<double> elemFiber;
361 ✗ getFiberDirection(iel,m_ipotTransMemb, elemFiber);
362 // m_Matrix[2] += (\sigma_i^l-\sigma_i^t) a vec a grad(phi_i) * grad(phi_j)
363 ✗ m_elementMat[m_idK]->tau_grad_phi_i_tau_grad_phi_j(FelisceParam::instance().intraFiberTensor-FelisceParam::instance().intraTransvTensor,elemFiber,*m_fePotTransMemb,0,0,1);
364 }
365 }
366 }
367 ✗ break;
368 }
369 ✗ case 3: { // Bidomain (coupled)
370 ✗ m_fePotTransMemb->updateMeasQuadPt(0, elemPoint);
371 ✗ m_fePotTransMemb->updateFirstDeriv(0, elemPoint);
372
373 // m_Matrix[0] = [sigma_i grad(phi_i) * grad(phi_j) sigma_i grad(phi_i) * grad(phi_j) ]
374 // [sigma_i grad(phi_i) * grad(phi_j) (sigma_i+sigma_e) grad(phi_i) * grad(phi_j)]
375 ✗ m_elementMat[m_idL]->grad_phi_i_grad_phi_j(FelisceParam::instance().intraTransvTensor,*m_fePotTransMemb,0,0,1);
376 ✗ m_elementMat[m_idL]->grad_phi_i_grad_phi_j(FelisceParam::instance().intraTransvTensor,*m_fePotExtraCell,0,1,1);
377 ✗ m_elementMat[m_idL]->grad_phi_i_grad_phi_j(FelisceParam::instance().intraTransvTensor,*m_fePotTransMemb,1,0,1);
378 ✗ m_elementMat[m_idL]->grad_phi_i_grad_phi_j(FelisceParam::instance().intraTransvTensor+FelisceParam::instance().extraTransvTensor,*m_fePotExtraCell,1,1,1);
379
380 ✗ if (FelisceParam::instance().testCase == 1) {
381 ✗ std::vector<double> elemFiber;
382 ✗ getFiberDirection(iel,m_ipotTransMemb, elemFiber);
383 // m_Matrix[0] += (\sigma_i^l-\sigma_i^t) a vec a grad(phi_i) * grad(phi_j)
384 ✗ m_elementMat[m_idL]->tau_grad_phi_i_tau_grad_phi_j(FelisceParam::instance().intraFiberTensor-FelisceParam::instance().intraTransvTensor,elemFiber,*m_fePotTransMemb,0,0,1);
385 ✗ m_elementMat[m_idL]->tau_grad_phi_i_tau_grad_phi_j(FelisceParam::instance().intraFiberTensor-FelisceParam::instance().intraTransvTensor,elemFiber,*m_fePotExtraCell,0,1,1);
386 ✗ m_elementMat[m_idL]->tau_grad_phi_i_tau_grad_phi_j(FelisceParam::instance().intraFiberTensor-FelisceParam::instance().intraTransvTensor,elemFiber,*m_fePotTransMemb,1,0,1);
387 ✗ m_elementMat[m_idL]->tau_grad_phi_i_tau_grad_phi_j(FelisceParam::instance().intraFiberTensor-FelisceParam::instance().intraTransvTensor+FelisceParam::instance().extraFiberTensor-FelisceParam::instance().extraTransvTensor,elemFiber,*m_fePotExtraCell,1,1,1);
388 }
389
390 ✗ if (FelisceParam::instance().hasInitialCondition) {
391 // m_Matrix[0] += - m_coefChi * V * phi_i * phi_j
392 ✗ m_elemFieldU0.setValue(m_U_0_seq, *m_fePotTransMemb, iel, m_ipotTransMemb, m_ao, m_dof);
393 ✗ m_elementMat[m_idL]->a_phi_i_phi_j(- m_coefChi,m_elemFieldU0,*m_fePotTransMemb,0,0,1);
394 ✗ m_elemFieldUe0.setValue(m_Ue_0_seq, *m_fePotExtraCell, iel, m_ipotExtraCell, m_ao, m_dof);
395 ✗ m_elementMat[m_idL]->a_phi_i_phi_j(- m_coefChi,m_elemFieldUe0,*m_fePotExtraCell,1,1,1);
396 }
397
398 // m_Matrix[1] = phi_i * phi_j
399 ✗ m_elementMat[m_idG]->phi_i_phi_j(1.,*m_fePotTransMemb,0,0,1);
400 ✗ m_elementMat[m_idG]->phi_i_phi_j(1.,*m_fePotExtraCell,1,1,1);
401
402 ✗ break;
403 }
404 ✗ default:
405 ✗ FEL_ERROR("Model not defined for ALP solver.");
406 ✗ break;
407 }
408 }
409
410 // Projection functions from FE space to RO one
411 ✗ void EigenProblemALPDeim::projectOnBasis(double* vectorDof, double* vectorBasis, bool flagMass) {
412 ✗ if (FelisceParam::verbose() > 10 ) PetscPrintf(PETSC_COMM_WORLD,"\nEigenProblemALPDeim::projectOnBasis\n");
413
414 ✗ if (flagMass) {
415 ✗ for (int iBasis=0; iBasis<m_dimRomBasis; iBasis++) {
416 ✗ vectorBasis[iBasis] =0.;
417 ✗ for (int jPts=0; jPts<m_dimCollocationPts; jPts++) {
418 ✗ double tmpVal=0;
419 ✗ for (int kPts=0; kPts<m_dimCollocationPts; kPts++) {
420 ✗ tmpVal += m_massDeim[jPts][kPts]*m_basisDeim[kPts][iBasis];
421 }
422 ✗ vectorBasis[iBasis] += tmpVal*vectorDof[jPts];
423 }
424 }
425 } else {
426 ✗ for (int iBasis=0; iBasis<m_dimRomBasis; iBasis++) {
427 ✗ vectorBasis[iBasis] =0.;
428 ✗ for (int jPts=0; jPts<m_dimCollocationPts; jPts++) {
429 ✗ vectorBasis[iBasis] += m_basisDeim[jPts][iBasis]*vectorDof[jPts];
430 }
431 }
432 }
433
434 }
435
436 // Projection functions from RO space to FE one
437 ✗ void EigenProblemALPDeim::projectOnDof(double* vectorBasis, double* vectorDof, int size) {
438 ✗ if (FelisceParam::verbose() > 10 ) PetscPrintf(PETSC_COMM_WORLD,"\nEigenProblemALPDeim::projectOnDof\n");
439 ✗ for (int jPts=0; jPts<m_dimCollocationPts; jPts++) {
440 ✗ vectorDof[jPts] = 0.;
441 ✗ for (int iBasis=0; iBasis<size; iBasis++) {
442 ✗ vectorDof[jPts] += m_basisDeim[jPts][iBasis]*vectorBasis[iBasis];
443 }
444 }
445 }
446
447 ✗ void EigenProblemALPDeim::readData(IO& io) {
448 ✗ if (FelisceParam::verbose() > 0) PetscPrintf(PETSC_COMM_WORLD,"\nEigenProblemALPDeim::readData.\n");
449 ✗ m_Matrix[m_idG].getVecs(m_U_0,nullPetscVector);
450 ✗ m_U_0.setFromOptions();
451
452 ✗ m_U_0_seq.createSeq(PETSC_COMM_SELF,m_numDof);
453 ✗ m_U_0_seq.setFromOptions();
454
455 ✗ m_W_0.duplicateFrom(m_U_0);
456 ✗ m_W_0.copyFrom(m_U_0);
457
458 ✗ m_W_0_seq.createSeq(PETSC_COMM_SELF,m_numDof);
459 ✗ m_W_0_seq.setFromOptions();
460
461 ✗ m_Ue_0.duplicateFrom(m_U_0);
462 ✗ m_Ue_0.copyFrom(m_U_0);
463
464 ✗ m_Ue_0_seq.createSeq(PETSC_COMM_SELF,m_numDof);
465 ✗ m_Ue_0_seq.setFromOptions();
466
467 ✗ if (FelisceParam::instance().hasInitialCondition) {
468 // Read initial solution file (*.00000.scl)
469 ✗ int idInVar = 0;
470 // potTransMemb
471 ✗ double* initialSolution = nullptr;
472 ✗ initialSolution = new double[m_mesh->numPoints()];
473 ✗ io.readVariable(idInVar, 0.,initialSolution, m_mesh->numPoints());
474 ✗ idInVar ++;
475 // ionicVar
476 ✗ double* ionicSolution = nullptr;
477 ✗ if ( (m_problem == 1) || (m_problem == 2) ) {
478 ✗ ionicSolution = new double[m_mesh->numPoints()];
479 ✗ io.readVariable(idInVar, 0.,ionicSolution, m_mesh->numPoints());
480 ✗ idInVar ++;
481 }
482 // potExtraCell
483 ✗ double* extraCellSolution = nullptr;
484 ✗ if (m_problem == 2) {
485 ✗ extraCellSolution = new double[m_mesh->numPoints()];
486 ✗ io.readVariable(idInVar, 0.,extraCellSolution, m_mesh->numPoints());
487 ✗ idInVar ++;
488 }
489 // initial potential
490 ✗ double* initPot = nullptr;
491 ✗ if (FelisceParam::instance().optimizePotential) {
492 ✗ m_initPot.createSeq(PETSC_COMM_SELF,m_numDof);
493 ✗ m_initPot.setFromOptions();
494 ✗ initPot = new double[m_mesh->numPoints()];
495 ✗ io.readVariable(idInVar, 0.,initPot, m_mesh->numPoints());
496 ✗ idInVar ++;
497 }
498 // initialize (parallel) initial solution vectors
499 ✗ fromDoubleStarToVec(initialSolution, m_U_0, m_mesh->numPoints());
500 ✗ if ( (m_problem == 1) || (m_problem == 2) )
501 ✗ fromDoubleStarToVec(ionicSolution, m_W_0, m_mesh->numPoints());
502 ✗ if (m_problem == 2 )
503 ✗ fromDoubleStarToVec(extraCellSolution, m_Ue_0, m_mesh->numPoints());
504
505 // initialize SEQUENTIAL initial solution vectors (in order to use elemField and define collocation pts solution)
506 felInt ia;
507 ✗ for (felInt i=0; i< m_numDof; i++) {
508 ✗ ia = i;
509 ✗ AOApplicationToPetsc(m_ao,1,&ia);
510 ✗ m_U_0_seq.setValues(1,&ia,&initialSolution[i],INSERT_VALUES);
511 ✗ if ( (m_problem == 1) || (m_problem == 2) )
512 ✗ m_W_0_seq.setValues(1,&ia,&ionicSolution[i],INSERT_VALUES);
513 ✗ if (m_problem == 2)
514 ✗ m_Ue_0_seq.setValues(1,&ia,&extraCellSolution[i],INSERT_VALUES);
515 ✗ if (FelisceParam::instance().optimizePotential)
516 ✗ m_initPot.setValues(1,&ia,&initPot[i],INSERT_VALUES);
517 }
518 ✗ m_U_0_seq.assembly();
519 ✗ if ( (m_problem == 1) || (m_problem == 2) ) {
520 ✗ m_W_0_seq.assembly();
521 }
522 ✗ if (m_problem == 2) {
523 ✗ m_Ue_0_seq.assembly();
524 }
525 ✗ if (FelisceParam::instance().optimizePotential) {
526 ✗ m_initPot.assembly();
527 }
528
529 ✗ if (initialSolution) delete [] initialSolution;
530 ✗ if (ionicSolution) delete [] ionicSolution;
531 ✗ if (extraCellSolution) delete [] extraCellSolution;
532 ✗ if (initPot) delete [] initPot;
533
534 // Read fibers file (*.00000.vct)
535 ✗ if (FelisceParam::instance().testCase == 1) {
536 ✗ if (m_fiber == nullptr)
537 ✗ m_fiber = new double[m_mesh->numPoints()*3];
538 ✗ io.readVariable(idInVar, 0.,m_fiber, m_mesh->numPoints()*3);
539 ✗ idInVar ++;
540 }
541 }
542 else {
543 ✗ PetscPrintf(PETSC_COMM_WORLD, "U_0 not read from file (zero vector).");
544 ✗ m_U_0.set( 0.);
545 ✗ m_U_0.assembly();
546 ✗ if ( (m_problem == 1) || (m_problem == 2) ) {
547 ✗ if (FelisceParam::instance().typeOfIonicModel == "schaf") {
548 ✗ m_W_0.set( 1.);
549 ✗ m_W_0_seq.set( 1.);
550 }
551 else {
552 ✗ m_W_0.set( 0.);
553 ✗ m_W_0_seq.set( 0.);
554 }
555 ✗ m_W_0.assembly();
556 ✗ m_W_0_seq.assembly();
557 }
558 ✗ m_U_0_seq.set( 0.);
559 ✗ m_U_0_seq.assembly();
560 ✗ if (m_problem == 2) {
561 ✗ m_Ue_0.set( 0.);
562 ✗ m_Ue_0.assembly();
563 ✗ m_Ue_0_seq.set( 0.);
564 ✗ m_Ue_0_seq.assembly();
565 }
566 }
567
568 }
569
570 ✗ void EigenProblemALPDeim::readCollocationPts() {
571 ✗ if (FelisceParam::verbose() > 10 ) PetscPrintf(PETSC_COMM_WORLD,"\nEigenProblemALPDeim::readCollocationPts\n");
572 ✗ std::string collocationPtsFileName;
573 ✗ collocationPtsFileName = FelisceParam::instance().inputDirectory + "/" + FelisceParam::instance().collocationPtsFile;// load the collocation points file
574
575 felInt numCollocationNode;
576 ✗ std::ifstream collocationPtsFile(collocationPtsFileName.c_str());
577 ✗ if(!collocationPtsFile) {
578 ✗ FEL_ERROR("Fatal error: cannot open match file " + collocationPtsFileName );
579 }
580 ✗ collocationPtsFile >> numCollocationNode;
581
582 ✗ if (numCollocationNode != m_dimCollocationPts) {
583 ✗ FEL_WARNING("Wrong number of collocation points in ECGnode file.");
584 }
585
586 ✗ if (m_collocationNode == nullptr) {
587 ✗ m_collocationNode = new felInt[m_dimCollocationPts];
588
589 ✗ for(felInt j=0; j<m_dimCollocationPts; j++) {
590 ✗ collocationPtsFile >> m_collocationNode[j];
591 }
592 }
593 ✗ collocationPtsFile.close();
594 }
595
596
597 ✗ void EigenProblemALPDeim::initializeROM() {
598 ✗ if (FelisceParam::verbose() > 10 ) PetscPrintf(PETSC_COMM_WORLD,"\nEigenProblemALPDeim::initializeROM()\n");
599
600 // Basis in collocation points
601 ✗ if (m_basisDeimColloc == nullptr) {
602 ✗ m_basisDeimColloc = new double* [m_dimCollocationPts];
603 ✗ for (int i=0 ; i< m_dimCollocationPts; i++) {
604 ✗ m_basisDeimColloc[i] = new double[m_dimCollocationPts];
605 }
606 }
607 ✗ if (m_basisDeim == nullptr) {
608 ✗ m_basisDeim = new double* [m_dimCollocationPts];
609 ✗ for (int i=0 ; i< m_dimCollocationPts; i++) {
610 ✗ m_basisDeim[i] = new double[m_dimRomBasis];
611 }
612 }
613
614
615 ✗ if (m_useImprovedRec) {
616 ✗ if (m_orthCompDeim == nullptr){
617 ✗ if (m_dimRomBasis+m_dimOrthComp > m_dimCollocationPts) {
618 ✗ FEL_WARNING("Dimension of orthogonal component exceeds");
619 ✗ m_dimOrthComp = m_dimCollocationPts - m_dimRomBasis;
620 }
621 ✗ m_orthCompDeim = new double* [m_dimCollocationPts];
622 ✗ for (int i=0 ; i< m_dimCollocationPts; i++) {
623 ✗ m_orthCompDeim[i] = new double[m_dimOrthComp];
624 }
625 }
626 }
627
628 ✗ if (FelisceParam::instance().solveEigenProblem)
629 {
630 int rankProc;
631 ✗ MPI_Comm_rank(m_petscComm,&rankProc);
632 ✗ for (int i=0; i<m_dimCollocationPts; i++) {
633 ✗ double* vec = new double[m_numDof];
634 ✗ fromVecToDoubleStar(vec, m_basis[i], rankProc, 1, m_numDof);
635 // Collocation of basis : m_basisDeim = P(m_basis)
636 ✗ for (int j=0; j<m_dimCollocationPts; j++) {
637 ✗ int k=m_collocationNode[j]-1;
638 ✗ m_basisDeimColloc[j][i] = vec[k];
639 ✗ if (i<m_dimRomBasis) {
640 ✗ m_basisDeim[j][i] = m_basisDeimColloc[j][i];
641 }
642 }
643 ✗ delete [] vec;
644 ✗ if (FelisceParam::verbose() > 40)
645 ✗ m_basis[i].view();
646 }
647 }
648 else
649 {
650 // Basis in mesh points
651 ✗ m_basis.resize(m_dimCollocationPts);
652 ✗ for (int i=0; i<m_dimCollocationPts; i++) {
653 ✗ m_Matrix[m_idG].getVecs(m_basis[i],nullPetscVector);
654 ✗ m_basis[i].setFromOptions();
655 ✗ double* vec = new double[m_numDof];
656 ✗ readEnsightFile(vec,"basis",i,m_numDof);
657 ✗ fromDoubleStarToVec(vec, m_basis[i], m_numDof);
658 // Collocation of basis : m_basisDeim = P(m_basis)
659 ✗ for (int j=0; j<m_dimCollocationPts; j++) {
660 ✗ int k=m_collocationNode[j]-1;
661 ✗ m_basisDeimColloc[j][i] = vec[k];
662 ✗ if (i<m_dimRomBasis) {
663 ✗ m_basisDeim[j][i] = m_basisDeimColloc[j][i];
664 }
665 }
666 ✗ delete [] vec;
667 ✗ if (FelisceParam::verbose() > 40)
668 ✗ m_basis[i].view();
669 }
670
671 ✗ readEigenValueFromFile();
672 ✗ computeMassDeim();
673 ✗ initializeSolution();
674
675 }
676
677 ✗ if (m_useImprovedRec && m_improvedRecType==0) {
678 ✗ for (int i=0; i<m_dimCollocationPts; i++) {
679 ✗ for (int j=0; j<m_dimOrthComp; j++) {
680 ✗ m_orthCompDeim[i][j] = m_basisDeimColloc[i][j+m_dimRomBasis];
681 }
682 }
683 }
684
685 }
686
687 ✗ void EigenProblemALPDeim::initializeSolution() {
688 ✗ if (FelisceParam::verbose() > 10 ) PetscPrintf(PETSC_COMM_WORLD,"\nEigenProblemALPDeim::initializeSolution()\n");
689
690 ✗ if ( !m_solutionInitialized) {
691
692 ✗ if (m_hatU == nullptr)
693 ✗ m_hatU = new double[m_dimCollocationPts];
694 // m_hatU = P(m_U_0)
695 ✗ for (int j=0; j<m_dimCollocationPts; j++) {
696 ✗ felInt iout = m_collocationNode[j]-1;
697 ✗ AOApplicationToPetsc(m_ao, 1, &iout);
698 double valueVec;
699 ✗ m_U_0_seq.getValues(1, &iout, &valueVec);
700 ✗ m_hatU[j] = valueVec;
701 }
702
703 // m_beta = \Phi*m_hatU
704 ✗ if (m_beta == nullptr)
705 ✗ m_beta = new double[m_dimRomBasis];
706 //EigenProblemALP::
707 ✗ projectOnBasis(m_U_0,m_beta,true);
708
709 // Ionic variable W
710 ✗ if ( (m_problem == 1) || (m_problem == 2) ) {
711 ✗ if (m_hatW == nullptr)
712 ✗ m_hatW = new double[m_dimCollocationPts];
713 // m_hatW = P(m_W_0)
714 ✗ for (int j=0; j<m_dimCollocationPts; j++) {
715 ✗ felInt iout = m_collocationNode[j]-1;
716 ✗ AOApplicationToPetsc(m_ao, 1, &iout);
717 double valueVec;
718 ✗ m_W_0_seq.getValues(1, &iout, &valueVec);
719 ✗ m_hatW[j] = valueVec;
720 }
721
722 }
723
724 // Extra-cell potential
725 ✗ if (m_problem == 2) {
726 ✗ if (m_hatUe == nullptr)
727 ✗ m_hatUe = new double[m_dimCollocationPts];
728 // m_hatUe = P(m_Ue_0)
729 ✗ for (int j=0; j<m_dimCollocationPts; j++) {
730 ✗ felInt iout = m_collocationNode[j]-1;
731 ✗ AOApplicationToPetsc(m_ao, 1, &iout);
732 double valueVec;
733 ✗ m_Ue_0_seq.getValues(1, &iout, &valueVec);
734 ✗ m_hatUe[j] = valueVec;
735 }
736 // m_xi = \Phi*m_hatUe
737 ✗ if (m_xi == nullptr)
738 ✗ m_xi = new double[m_dimRomBasis];
739 //EigenProblemALP::
740 ✗ projectOnBasis(m_Ue_0,m_xi,true);
741 }
742
743 ✗ if (FelisceParam::instance().optimizePotential) {
744 ✗ if (m_potential == nullptr)
745 ✗ m_potential = new double[m_dimCollocationPts];
746 // m_potential = P(m_initPot)
747 ✗ for (int j=0; j<m_dimCollocationPts; j++) {
748 ✗ felInt iout = m_collocationNode[j]-1;
749 ✗ AOApplicationToPetsc(m_ao, 1, &iout);
750 double valueVec;
751 ✗ m_initPot.getValues(1, &iout, &valueVec);
752 ✗ m_potential[j] = valueVec;
753 }
754 }
755
756 ✗ if (m_resPoint) {
757 ✗ if (m_resU == nullptr) {
758 ✗ m_resU = new double[m_dimCollocationPts];
759 }
760
761 ✗ for (int i=0; i<m_dimCollocationPts; i++) {
762 ✗ m_resU[i] = m_hatU[i];
763 ✗ for (int j=0; j<m_dimRomBasis; j++) {
764 ✗ m_resU[i] -= m_basisDeim[i][j]*m_beta[j];
765 }
766 }
767
768 ✗ double maxRes = 0.;
769 ✗ for (int i=0; i<m_dimCollocationPts; i++) {
770 ✗ if (m_resU[i] > maxRes) {
771 ✗ maxRes = m_resU[i];
772 }
773 }
774 ✗ std::cout << "maxRes = " << maxRes << std::endl;
775
776 }
777
778 ✗ m_solutionInitialized = true;
779 }
780
781 }
782
783 // Modified Graam-Schmidt: orthonormalization (in the sense of m_massDeim) of A
784 ✗ void EigenProblemALPDeim::MGS(double** A, double** R, int size) {
785 ✗ if (FelisceParam::verbose() > 10 ) PetscPrintf(PETSC_COMM_WORLD,"\nEigenProblemALPDeim::MGS\n");
786
787 ✗ for(int i=0; i<size; i++) {
788
789 // tmp = G phi_i
790 ✗ double tmp[m_dimCollocationPts];
791 ✗ for (int j=0; j<m_dimCollocationPts; j++) {
792 ✗ tmp[j] = 0.;
793 ✗ for (int k=0; k<m_dimCollocationPts; k++)
794 ✗ tmp[j] += m_massDeim[j][k]*A[k][i];
795 }
796
797 // Diagonal of R
798 ✗ R[i][i] =0.;
799 ✗ for (int k=0; k<m_dimCollocationPts; k++)
800 ✗ R[i][i] += A[k][i] * tmp[k];
801
802 ✗ R[i][i] = std::sqrt(R[i][i]);
803 ✗ for (int k=0; k<m_dimCollocationPts; k++)
804 ✗ A[k][i] = A[k][i]/R[i][i];
805
806 // Extradiag R
807 ✗ for (int j=i+1; j<size; j++) {
808 ✗ R[i][j] =0.;
809 ✗ for (int k=0; k<m_dimCollocationPts; k++)
810 ✗ R[i][j] += A[k][j]*tmp[k];
811 ✗ R[j][i] = 0.;
812 // A
813 ✗ for (int k=0; k<m_dimCollocationPts; k++)
814 ✗ A[k][j] -= R[i][j] * A[k][i];
815 }
816 }
817
818 }
819
820 ✗ void EigenProblemALPDeim::factorizationQR(double** A, double** Q, double** R, int sizeRow, int sizeCol) {
821 ✗ if (FelisceParam::verbose() > 10 ) PetscPrintf(PETSC_COMM_WORLD,"\nEigenProblemALPDeim::factorizationQR\n");
822 // Copy A into W
823 ✗ double** W = new double* [sizeRow];
824 ✗ for (int i=0; i<sizeRow; i++) {
825 ✗ W[i] = new double [sizeCol];
826 ✗ for (int j=0; j<sizeCol; j++)
827 ✗ W[i][j] = A[i][j];
828 }
829
830 ✗ for(int i=0; i<sizeCol; i++) {
831 // Diagonal of R
832 ✗ R[i][i] =0.;
833 ✗ for (int j=0; j<sizeRow; j++)
834 ✗ R[i][i] += W[j][i] * W[j][i];
835 ✗ R[i][i] = std::sqrt(R[i][i]);
836
837 // Q
838 ✗ for (int j=0; j<sizeRow; j++)
839 ✗ Q[j][i] = W[j][i] / R[i][i];
840
841 // Extradiag R
842 ✗ for (int k=i+1; k<sizeCol; k++) {
843 ✗ R[i][k] =0.;
844 ✗ for (int j=0; j<sizeRow; j++)
845 ✗ R[i][k] += Q[j][i] * W[j][k];
846 ✗ R[k][i] = 0.;
847 // W
848 ✗ for (int j=0; j<sizeRow; j++)
849 ✗ W[j][k] -= R[i][k] * Q[j][i];
850 }
851 }
852
853 ✗ for (int i=0; i<sizeRow; i++)
854 ✗ delete [] W[i];
855 ✗ delete [] W;
856
857 }
858
859 ✗ void EigenProblemALPDeim::inverseTriangular(double** A, double** Ainv, int sizeA) {
860 ✗ if (FelisceParam::verbose() > 10 ) PetscPrintf(PETSC_COMM_WORLD,"\nEigenProblemALPDeim::inverseTriangular\n");
861 ✗ double det=1.;
862 ✗ for (int l=0; l<sizeA; l++) {
863 ✗ det = det*A[l][l];
864 }
865 ✗ if (det < 1.e-7) {
866 ✗ std::cout << "det(R) = " << det << std::endl;
867 ✗ FEL_ERROR("W is singular.");
868 }
869
870 ✗ for (int l=0; l<sizeA; l++) {
871 double* tmpX;
872 ✗ tmpX = new double[l+1];
873 ✗ tmpX[l] = 1.0/(A[l][l]);
874 ✗ for (int j=l-1; j>=0; j--) {
875 ✗ tmpX[j] = 0.0;
876 ✗ for (int k=j+1; k<=l; k++) {
877 ✗ tmpX[j] -=A[j][k]*tmpX[k];
878 }
879 ✗ tmpX[j] = tmpX[j]/A[j][j];
880 }
881
882 ✗ for (int j=0; j<l+1; j++) {
883 ✗ Ainv[j][l] = tmpX[j];
884 ✗ if (j != l)
885 ✗ Ainv[l][j] = 0.;
886 }
887
888 ✗ delete [] tmpX;
889 }
890
891 }
892
893 // Projection on collocation points of mass matrix
894 ✗ void EigenProblemALPDeim::computeMassDeim() {
895 ✗ if (FelisceParam::verbose() > 10 ) PetscPrintf(PETSC_COMM_WORLD,"\nEigenProblemALPDeim::computeMassDeim\n");
896
897 ✗ if (m_massDeim == nullptr) {
898 ✗ m_massDeim = new double* [m_dimCollocationPts];
899 ✗ for (int i=0; i< m_dimCollocationPts; i++) {
900 ✗ m_massDeim[i] = new double [m_dimCollocationPts];
901 }
902 }
903 ✗ if (m_useImprovedRec || m_resPoint) {
904 ✗ if (m_stiffDeim == nullptr) {
905 ✗ m_stiffDeim = new double* [m_dimCollocationPts];
906 ✗ for (int i=0; i< m_dimCollocationPts; i++) {
907 ✗ m_stiffDeim[i] = new double [m_dimCollocationPts];
908 }
909 }
910 }
911
912 double** matrixQ;
913 ✗ matrixQ = new double* [m_dimCollocationPts];
914 ✗ for (int i=0; i< m_dimCollocationPts; i++) {
915 ✗ matrixQ[i] = new double [m_dimCollocationPts];
916 }
917 double** matrixR;
918 ✗ matrixR = new double* [m_dimCollocationPts];
919 ✗ for (int i=0; i< m_dimCollocationPts; i++) {
920 ✗ matrixR[i] = new double [m_dimCollocationPts];
921 }
922 double** matrixRinv;
923 ✗ matrixRinv = new double* [m_dimCollocationPts];
924 ✗ for (int i=0; i< m_dimCollocationPts; i++) {
925 ✗ matrixRinv[i] = new double [m_dimCollocationPts];
926 }
927 double** matrixH;
928 ✗ matrixH = new double* [m_dimCollocationPts];
929 ✗ for (int i=0; i< m_dimCollocationPts; i++) {
930 ✗ matrixH[i] = new double [m_dimCollocationPts];
931 }
932
933 ✗ factorizationQR(m_basisDeimColloc, matrixQ, matrixR, m_dimCollocationPts, m_dimCollocationPts);
934 #ifndef NDEBUG
935 ✗ double checkQR = 0.;
936 ✗ for (int i=0; i<m_dimCollocationPts; i++) {
937 ✗ for (int j=0; j<m_dimCollocationPts; j++) {
938 ✗ checkQR = m_basisDeimColloc[i][j];
939 ✗ for (int k=0; k<m_dimCollocationPts; k++) {
940 ✗ checkQR -= matrixQ[i][k]*matrixR[k][j];
941 }
942 ✗ if (std::abs(checkQR) > 1.e-3) {
943 ✗ FEL_WARNING("Error in QR factorization of matrix W.");
944 }
945 }
946 }
947 #endif
948 ✗ inverseTriangular(matrixR, matrixRinv, m_dimCollocationPts);
949 #ifndef NDEBUG
950 ✗ double checkRinv = 0.;
951 ✗ for (int i=0; i<m_dimCollocationPts; i++) {
952 ✗ for (int j=0; j<m_dimCollocationPts; j++) {
953 ✗ checkRinv = 0.;
954 ✗ if ( i==j ) checkRinv = 1.;
955 ✗ for (int k=0; k<m_dimCollocationPts; k++) {
956 ✗ checkRinv -= matrixR[i][k]*matrixRinv[k][j];
957 }
958 ✗ if (std::abs(checkRinv) > 1.e-3) {
959 ✗ FEL_WARNING("Error in inverse of triangular matrix R.");
960 }
961 }
962 }
963 #endif
964 // H = R^-1 Q^T
965 ✗ for (int iG=0; iG< m_dimCollocationPts; iG++) {
966 ✗ for (int jG=0; jG< m_dimCollocationPts; jG++) {
967 ✗ matrixH[iG][jG] =0.;
968 ✗ for (int kSum=0; kSum< m_dimCollocationPts; kSum++)
969 ✗ matrixH[iG][jG] += matrixRinv[iG][kSum] * matrixQ[jG][kSum];
970 }
971 }
972
973 // G = H^T H
974 ✗ for (int iG=0; iG< m_dimCollocationPts; iG++) {
975 ✗ for (int jG=0; jG< m_dimCollocationPts; jG++) {
976 ✗ m_massDeim[iG][jG] = 0.;
977 ✗ for (int kSum=0; kSum< m_dimCollocationPts; kSum++)
978 ✗ m_massDeim[iG][jG] += matrixH[kSum][iG] * matrixH[kSum][jG] ;
979 }
980 }
981
982 ✗ if (m_useImprovedRec || m_resPoint) {
983 // K : W^T K W = V^T K_fem V
984 ✗ std::vector<PetscVector> matrixX;
985 ✗ std::vector<PetscVector> matrixY;
986 ✗ matrixX.resize(m_dimCollocationPts);
987 ✗ matrixY.resize(m_dimCollocationPts);
988 ✗ for (int i=0; i< m_dimCollocationPts; i++) {
989 ✗ matrixX[i].duplicateFrom(m_basis[0]);
990 ✗ matrixY[i].duplicateFrom(m_basis[0]);
991 }
992 // X = V H
993 // V = m_basis
994 ✗ for (int i=0; i<m_dimCollocationPts; i++) {
995 ✗ for (int j=0; j<m_dimCollocationPts; j++) {
996 ✗ matrixX[i].axpy(matrixH[j][i],m_basis[j]);
997 }
998 }
999 // Y = K_fem X
1000 // K_fem = m_Matrix[m_idK];
1001 ✗ for (int i=0; i<m_dimCollocationPts; i++) {
1002 ✗ mult(m_Matrix[m_idK],matrixX[i],matrixY[i]);
1003 }
1004 // K = X^T K_fem X = X^T Y
1005 ✗ for (int i=0; i<m_dimCollocationPts; i++) {
1006 ✗ for (int j=0; j<m_dimCollocationPts; j++) {
1007 double valueK;
1008 ✗ dot(matrixX[i],matrixY[j], &valueK);
1009 ✗ m_stiffDeim[i][j] = valueK;
1010 }
1011 }
1012 ✗ for (int i=0; i< m_dimCollocationPts; i++) {
1013 ✗ matrixX[i].destroy();
1014 ✗ matrixY[i].destroy();
1015 }
1016 }
1017
1018 ✗ if (m_resPoint) {
1019 ✗ m_invGK = new double* [m_dimCollocationPts];
1020 ✗ for (int i=0; i<m_dimCollocationPts; i++) {
1021 ✗ m_invGK[i] = new double[m_dimCollocationPts];
1022 }
1023
1024 ✗ double** invG = new double* [m_dimCollocationPts];
1025 ✗ for (int i=0; i<m_dimCollocationPts; i++) {
1026 ✗ invG[i] = new double[m_dimCollocationPts];
1027 }
1028
1029 ✗ for (int i=0; i<m_dimCollocationPts; i++) {
1030 ✗ for (int j=i; j<m_dimCollocationPts; j++) {
1031 ✗ invG[i][j] = 0.;
1032 ✗ for (int k=0; k<m_dimCollocationPts; k++) {
1033 ✗ invG[i][j] += m_basisDeimColloc[i][k]*m_basisDeimColloc[j][k];
1034 }
1035 ✗ invG[j][i] = invG[i][j];
1036 }
1037 }
1038
1039 ✗ for (int i=0; i<m_dimCollocationPts; i++) {
1040 ✗ for (int j=i; j<m_dimCollocationPts; j++) {
1041 ✗ m_invGK[i][j] = 0.;
1042 ✗ for (int k=0; k<m_dimCollocationPts; k++) {
1043 ✗ m_invGK[i][j] += invG[i][k]*m_stiffDeim[k][j];
1044 }
1045 ✗ m_invGK[j][i] = m_invGK[i][j];
1046 }
1047 }
1048
1049 ✗ for (int i=0; i< m_dimCollocationPts; i++) {
1050 ✗ delete [] invG[i];
1051 }
1052 ✗ delete [] invG;
1053
1054 }
1055
1056
1057 ✗ for (int i=0; i< m_dimCollocationPts; i++) {
1058 ✗ delete [] matrixQ[i];
1059 ✗ delete [] matrixR[i];
1060 ✗ delete [] matrixRinv[i];
1061 ✗ delete [] matrixH[i];
1062 }
1063 ✗ delete [] matrixQ;
1064 ✗ delete [] matrixR;
1065 ✗ delete [] matrixRinv;
1066 ✗ delete [] matrixH;
1067 }
1068
1069 ✗ void EigenProblemALPDeim::computeProjectionPiV() {
1070 ✗ if (FelisceParam::verbose() > 10 ) PetscPrintf(PETSC_COMM_WORLD,"\nEigenProblemALPDeim::computeProjectionPiV\n");
1071
1072 ✗ if (m_projectorPiV.size() < 1) {
1073 ✗ std::vector<PetscVector> projectorPiW;
1074 ✗ projectorPiW.resize(m_dimCollocationPts);
1075 ✗ for (int i=0; i<m_dimCollocationPts; i++) {
1076 ✗ m_Matrix[m_idG].getVecs(projectorPiW[i],nullPetscVector);
1077 ✗ projectorPiW[i].setFromOptions();
1078 ✗ for (int k=0; k<m_dimCollocationPts; k++)
1079 ✗ projectorPiW[i].axpy( m_basisDeimColloc[i][k], m_basis[k]);
1080 }
1081
1082 ✗ m_projectorPiV.resize(m_dimCollocationPts);
1083 ✗ for (int i=0; i<m_dimCollocationPts; i++) {
1084 ✗ m_Matrix[m_idG].getVecs(m_projectorPiV[i],nullPetscVector);
1085 ✗ m_projectorPiV[i].setFromOptions();
1086 ✗ for (int k=0; k<m_dimCollocationPts; k++)
1087 ✗ m_projectorPiV[i].axpy( m_massDeim[k][i], projectorPiW[k]);
1088 }
1089 ✗ for (int i=0; i<m_dimCollocationPts; i++) {
1090 ✗ projectorPiW[i].destroy();
1091 }
1092 }
1093
1094 }
1095
1096 ✗ void EigenProblemALPDeim::writeEnsightSolution(const int iter) {
1097 ✗ if (FelisceParam::verbose() > 10 ) PetscPrintf(PETSC_COMM_WORLD,"\nEigenProblemALPDeim::writeEnsightSolution\n");
1098 ✗ PetscVector solU;
1099 ✗ solU.duplicateFrom(m_U_0);
1100 ✗ PetscVector solUe;
1101 ✗ solUe.duplicateFrom(m_U_0);
1102 ✗ PetscVector solW;
1103 ✗ solW.duplicateFrom(m_W_0);
1104
1105 ✗ if (m_projectorPiV.size() < 1)
1106 ✗ computeProjectionPiV();
1107
1108 ✗ if (iter > 0) {
1109 ✗ for (int k=0; k<m_dimCollocationPts; k++) {
1110 ✗ solU.axpy( m_hatU[k], m_projectorPiV[k]);
1111 }
1112 ✗ if ( (m_problem == 1) || (m_problem == 2)) {
1113 ✗ if (FelisceParam::instance().printIonicVar) {
1114 ✗ for (int k=0; k<m_dimCollocationPts; k++) {
1115 ✗ solW.axpy( m_hatW[k], m_projectorPiV[k]);
1116 }
1117 }
1118 }
1119 ✗ if (m_problem == 2) {
1120 ✗ for (int k=0; k<m_dimCollocationPts; k++) {
1121 ✗ solUe.axpy( m_hatUe[k], m_projectorPiV[k]);
1122 }
1123 }
1124 } else {
1125 ✗ solU.copyFrom(m_U_0);
1126 ✗ if ( (m_problem == 1) || (m_problem == 2)) {
1127 ✗ if (FelisceParam::instance().printIonicVar) {
1128 ✗ solW.copyFrom(m_W_0);
1129 }
1130 }
1131 ✗ if (m_problem == 2) {
1132 ✗ solUe.copyFrom(m_Ue_0);
1133 }
1134 }
1135
1136 int rankProc;
1137 ✗ MPI_Comm_rank(m_petscComm,&rankProc);
1138
1139 ✗ double* tmpSolToSave = nullptr;
1140 ✗ tmpSolToSave = new double[m_numDof];
1141
1142 ✗ if (m_problem == 0) {
1143 ✗ fromVecToDoubleStar(tmpSolToSave, solU, rankProc, 1);
1144 ✗ writeEnsightVector(tmpSolToSave, iter, "sol");
1145 }
1146
1147 ✗ if ( (m_problem == 1) || (m_problem == 2)) {
1148 ✗ fromVecToDoubleStar(tmpSolToSave, solU, rankProc, 1);
1149 ✗ writeEnsightVector(tmpSolToSave, iter, "potTransMemb");
1150 ✗ if (FelisceParam::instance().printIonicVar) {
1151 ✗ fromVecToDoubleStar(tmpSolToSave, solW, rankProc, 1);
1152 ✗ writeEnsightVector(tmpSolToSave, iter, "ionicVar");
1153 }
1154 }
1155
1156 ✗ if (m_problem == 2) {
1157 ✗ fromVecToDoubleStar(tmpSolToSave, solUe, rankProc, 1);
1158 ✗ writeEnsightVector(tmpSolToSave, iter, "potExtraCell");
1159 }
1160
1161 ✗ delete [] tmpSolToSave;
1162 ✗ solU.destroy();
1163 ✗ solUe.destroy();
1164 ✗ solW.destroy();
1165
1166 ✗ std::vector<std::string> varNames;
1167 ✗ if (m_problem == 0) {
1168 ✗ varNames.resize(1);
1169 ✗ varNames[0] = "sol";
1170 ✗ } else if (m_problem == 1) {
1171 ✗ if (FelisceParam::instance().printIonicVar) {
1172 ✗ varNames.resize(2);
1173 ✗ varNames[0] = "potTransMemb";
1174 ✗ varNames[1] = "ionicVar";
1175 } else {
1176 ✗ varNames.resize(1);
1177 ✗ varNames[0] = "potTransMemb";
1178 }
1179 ✗ } else if (m_problem == 2) {
1180 ✗ if (FelisceParam::instance().printIonicVar) {
1181 ✗ varNames.resize(3);
1182 ✗ varNames[0] = "potTransMemb";
1183 ✗ varNames[1] = "potExtraCell";
1184 ✗ varNames[2] = "ionicVar";
1185 } else {
1186 ✗ varNames.resize(2);
1187 ✗ varNames[0] = "potTransMemb";
1188 ✗ varNames[1] = "potExtraCell";
1189 }
1190 }
1191 ✗ if (rankProc == 0)
1192 ✗ writeEnsightCase(iter+1, FelisceParam::instance().timeStep * FelisceParam::instance().frequencyWriteSolution,varNames, FelisceParam::instance().time);
1193
1194 }
1195
1196 ✗ void EigenProblemALPDeim::computeRHSDeim() {
1197 ✗ if (FelisceParam::verbose() > 10 ) PetscPrintf(PETSC_COMM_WORLD,"\nEigenProblemALPDeim::computeRHSDeim\n");
1198
1199 ✗ if (m_hatF == nullptr) {
1200 ✗ m_hatF = new double[m_dimCollocationPts];
1201 }
1202
1203 ✗ switch (m_problem) {
1204 ✗ case 0: {
1205 // FKPP:
1206 // TO DO
1207 ✗ FEL_ERROR("FKPP RHS NOT DEFINED");
1208 ✗ break;
1209 }
1210 ✗ case 1:
1211 case 2: {
1212 ✗ if (m_hatG == nullptr) {
1213 ✗ m_hatG = new double[m_dimCollocationPts];
1214 }
1215
1216 // MONODOMAIN or BIDOMAIN:
1217 ✗ double Am = FelisceParam::instance().Am;
1218 ✗ double Cm = FelisceParam::instance().Cm;
1219 // FhN
1220 ✗ double a = FelisceParam::instance().alpha;
1221 ✗ double s = FelisceParam::instance().f0;
1222 ✗ double eps = FelisceParam::instance().epsilon;
1223 ✗ double gam = FelisceParam::instance().gammaEl;
1224 // Schaf
1225 ✗ double& tauIn = FelisceParam::instance().tauIn;
1226 ✗ double& tauOut = FelisceParam::instance().tauOut;
1227 ✗ double& tauOpen = FelisceParam::instance().tauOpen;
1228 ✗ double& tauClose = FelisceParam::instance().tauClose;
1229
1230
1231 ✗ for (int i=0; i<m_dimCollocationPts; i++) {
1232 ✗ m_hatF[i] = 0.;
1233
1234 ✗ if(m_resPoint) {
1235 ✗ for (int k=0; k<m_dimCollocationPts; k++) {
1236 ✗ m_hatF[i] += m_invGK[i][k]*m_resU[k];
1237 }
1238 }
1239
1240 ✗ for (int j=0; j<m_dimRomBasis; j++) {
1241 ✗ double coef = m_beta[j];
1242 ✗ if (m_problem == 2)
1243 ✗ coef += m_xi[j];
1244 ✗ if (FelisceParam::instance().optimizePotential) {
1245 ✗ m_hatF[i] += coef * (- m_eigenValue[j] - m_coefChi*m_potential[i]) * m_basisDeim[i][j];
1246 }
1247 else {
1248 ✗ m_hatF[i] += coef * (- m_eigenValue[j] - m_coefChi*m_hatU[i]) * m_basisDeim[i][j];
1249 }
1250 }
1251 ✗ if (FelisceParam::instance().typeOfIonicModel == "fhn") {
1252 ✗ if (FelisceParam::instance().hasInfarct)
1253 ✗ s = m_coefFhNs[i];
1254 ✗ m_hatF[i] += s * Am * ( - m_hatU[i]*m_hatU[i]*m_hatU[i] + (1+a) * m_hatU[i]*m_hatU[i] - a * m_hatU[i] ) - Am * m_hatW[i];
1255 ✗ m_hatG[i] = eps * ( gam * m_hatU[i] - m_hatW[i] );
1256 }
1257 ✗ else if (FelisceParam::instance().typeOfIonicModel == "schaf") {
1258 ✗ if (FelisceParam::instance().hasInfarct)
1259 ✗ tauOut = m_coefTauOut[i];
1260 ✗ m_hatF[i] += Am * ( m_hatW[i] * m_hatU[i] * m_hatU[i] * (1 - m_hatU[i]) / tauIn - m_hatU[i] / tauOut);
1261
1262 ✗ if (m_hatU[i] < FelisceParam::instance().vGate) {
1263 ✗ m_hatG[i] = (1 - m_hatW[i]) / tauOpen;
1264 }
1265 else {
1266 ✗ if (FelisceParam::instance().hasHeteroTauClose) {
1267 ✗ tauClose = m_coefTauClose[i];
1268 }
1269 ✗ m_hatG[i] = - m_hatW[i] / tauClose;
1270 }
1271 }
1272
1273 ✗ if (FelisceParam::instance().hasSource) {
1274 ✗ double delay[m_numSource];
1275 ✗ double tPeriod=fmod(m_fstransient->time,FelisceParam::instance().timePeriod);
1276 ✗ for (int iS=0; iS<m_numSource; iS++) {
1277 ✗ delay[iS] = FelisceParam::instance().delayStim[iS];
1278 ✗ if ((tPeriod>=delay[iS])&&(tPeriod<=(delay[iS]+FelisceParam::instance().stimTime)))
1279 ✗ m_hatF[i] += Am*m_source[iS][i];
1280 }
1281 }
1282
1283 ✗ m_hatF[i] = m_hatF[i] / (Am*Cm);
1284 }
1285 }
1286 }
1287
1288 }
1289
1290 ✗ void EigenProblemALPDeim::computeTheta() {
1291 ✗ if (FelisceParam::verbose() > 10 ) PetscPrintf(PETSC_COMM_WORLD,"\nEigenProblemALPDeim::computeTheta\n");
1292 ✗ if (m_theta == nullptr) {
1293 ✗ m_theta = new double* [m_dimRomBasis];
1294 ✗ for (int i=0; i< m_dimRomBasis; i++) {
1295 ✗ m_theta[i] = new double [m_dimRomBasis];
1296 }
1297 ✗ if(m_useImprovedRec && m_improvedRecType==0) {
1298 ✗ m_thetaOrth = new double* [m_dimOrthComp];
1299 ✗ for (int i=0; i< m_dimOrthComp; i++) {
1300 ✗ m_thetaOrth[i] = new double [m_dimRomBasis];
1301 }
1302 }
1303 }
1304
1305 ✗ double* tmp = new double [m_dimCollocationPts];
1306
1307 // \theta = ( \Phi^T G ) (F * \Phi)
1308
1309 ✗ double* extraF = nullptr;
1310 ✗ if (m_problem == 0) {
1311 ✗ double Cm = FelisceParam::instance().Cm;
1312 ✗ extraF = new double[m_dimCollocationPts];
1313 ✗ for (int k=0; k< m_dimCollocationPts; k++) {
1314 ✗ extraF[k] = Cm * m_hatU[k];
1315 ✗ for (int id=0; id< m_dimRomBasis; id++) {
1316 ✗ extraF[k] += - m_eigenValue[id] * m_beta[id] * m_basisDeim[k][id];
1317 }
1318 }
1319 }
1320
1321 ✗ for (int id=0; id< m_dimRomBasis; id++) {
1322 ✗ for (int jd=0; jd< m_dimCollocationPts; jd++) {
1323 ✗ tmp[jd] = 0.;
1324 ✗ for (int k=0; k< m_dimCollocationPts; k++)
1325 ✗ tmp[jd] += m_basisDeim[k][id] * m_massDeim[k][jd];
1326 }
1327
1328 ✗ for (int jd=0; jd< id+1; jd++) {
1329 ✗ m_theta[id][jd] = 0.;
1330 ✗ for (int k=0; k< m_dimCollocationPts; k++) {
1331 ✗ m_theta[id][jd] += tmp[k] * m_hatF[k] * m_basisDeim[k][jd];
1332 ✗ if (m_problem == 0)
1333 ✗ m_theta[id][jd] += tmp[k] * extraF[k] * m_basisDeim[k][jd];
1334 }
1335 ✗ if (jd != id)
1336 ✗ m_theta[jd][id] = m_theta[id][jd];
1337 }
1338 }
1339 ✗ if(m_useImprovedRec && m_improvedRecType==0) {
1340 // \theta^orth = ( \Psi^T G ) (F * \Phi)
1341 ✗ for (int id=0; id< m_dimOrthComp; id++) {
1342 ✗ for (int jd=0; jd< m_dimCollocationPts; jd++) {
1343 ✗ tmp[jd] = 0.;
1344 ✗ for (int k=0; k< m_dimCollocationPts; k++)
1345 ✗ tmp[jd] += m_orthCompDeim[k][id] * m_massDeim[k][jd];
1346 }
1347
1348 ✗ for (int jd=0; jd< m_dimRomBasis; jd++) {
1349 ✗ m_thetaOrth[id][jd] = 0.;
1350 ✗ for (int k=0; k< m_dimCollocationPts; k++)
1351 ✗ m_thetaOrth[id][jd] += tmp[k] * m_hatF[k] * m_basisDeim[k][jd];
1352 }
1353 }
1354 }
1355
1356 ✗ if (extraF) delete [] extraF;
1357 ✗ delete [] tmp;
1358
1359 }
1360
1361 ✗ void EigenProblemALPDeim::computeMatrixM() {
1362 ✗ if (FelisceParam::verbose() > 10 ) PetscPrintf(PETSC_COMM_WORLD,"\nEigenProblemALPDeim::computeMatrixM()\n");
1363 ✗ if ( m_matrixM == nullptr ) {
1364 ✗ m_matrixM = new double*[m_dimRomBasis];
1365 ✗ for (int i=0; i<m_dimRomBasis; i++) {
1366 ✗ m_matrixM[i] = new double[m_dimRomBasis];
1367 }
1368 }
1369
1370 ✗ double epsilonLambda = 5.e-01; //1.e-01;
1371 ✗ for (int i=0; i<m_dimRomBasis; i++) {
1372 ✗ m_matrixM[i][i] = 0.;
1373 ✗ for (int j=0; j<i; j++) {
1374 ✗ m_matrixM[i][j] = 0.;
1375 ✗ if ( std::abs(m_eigenValue[i] - m_eigenValue[j] ) > epsilonLambda ) {
1376 ✗ m_matrixM[i][j] = m_coefChi / ( m_eigenValue[j] - m_eigenValue[i] ) * m_theta[i][j];
1377 }
1378 ✗ m_matrixM[j][i] = (-1.0) * m_matrixM[i][j];
1379 }
1380 }
1381
1382 }
1383
1384 ✗ void EigenProblemALPDeim::computeHatU() {
1385 ✗ if (FelisceParam::verbose() > 10 ) PetscPrintf(PETSC_COMM_WORLD,"\nEigenProblemALPDeim::computeHatU\n");
1386
1387 ✗ projectOnDof(m_beta,m_hatU,m_dimRomBasis);
1388
1389 ✗ if (m_problem == 0) { // FKPP (fix max and min u)
1390 ✗ for (int i=0; i<m_dimCollocationPts; i++) {
1391 ✗ if (m_hatU[i] < 0.)
1392 ✗ m_hatU[i] = 0.;
1393 ✗ else if (m_hatU[i] > 1.)
1394 ✗ m_hatU[i] = 1.;
1395 }
1396 }
1397
1398 ✗ if (m_resPoint) {
1399 ✗ double maxRes = 0.;
1400 ✗ for (int i=0; i<m_dimCollocationPts; i++) {
1401 ✗ if (m_resU[i] > maxRes) {
1402 ✗ maxRes = m_resU[i];
1403 }
1404 ✗ m_hatU[i] += m_resU[i];
1405 }
1406 ✗ std::cout << "maxRes = " << maxRes << std::endl;
1407 }
1408
1409
1410 ✗ if ( (m_problem == 1) || (m_problem == 2)) {
1411 ✗ double dt = FelisceParam::instance().timeStep;
1412 ✗ for (int i=0; i<m_dimCollocationPts; i++) {
1413 ✗ m_hatW[i] += dt*m_hatG[i];
1414 }
1415 }
1416
1417 ✗ if (m_problem == 2)
1418 ✗ projectOnDof(m_xi,m_hatUe,m_dimRomBasis);
1419
1420 // Potential evolution
1421 // Hyp: dt(U) = dt(u) = hat(F(u))
1422 ✗ if (FelisceParam::instance().optimizePotential) {
1423 ✗ double dt = FelisceParam::instance().timeStep;
1424 ✗ for (int i=0; i<m_dimCollocationPts; i++) {
1425 ✗ m_potential[i] += dt*m_hatF[i]; // Warning! Not work like that!
1426 }
1427 }
1428
1429 }
1430
1431 ✗ void EigenProblemALPDeim::computeTensor() {
1432 ✗ if (FelisceParam::verbose() > 10 ) PetscPrintf(PETSC_COMM_WORLD,"\nEigenProblemALPDeim::computeTensor()\n");
1433
1434 ✗ if (m_problem ==2) {
1435 // If \sigma_I = \alpha \sigma_E
1436 // then A \xi + B \beta = 0
1437 // A_ij = \sum_{h,l}^{N_p} (\alpha + 1)(\chi m_hatU[l] + lambda_i) \phi_j(x_h) G_hl \phi_i(x_h)
1438 // "+" \int u_e = 0
1439 ✗ if (m_matrixA == nullptr) {
1440 ✗ m_matrixA = new double* [m_dimRomBasis+1];
1441 ✗ for (int i=0; i<m_dimRomBasis+1; i++) {
1442 ✗ m_matrixA[i] = new double[m_dimRomBasis+1];
1443 }
1444 }
1445 // B_ij = \sum_{h,l}^{N_p} (\chi m_hatU[l] + lambda_i) \phi_j(x_h) G_hl \phi_i(x_h)
1446 ✗ if (m_matrixB == nullptr) {
1447 ✗ m_matrixB = new double* [m_dimRomBasis];
1448 ✗ for (int i=0; i<m_dimRomBasis; i++) {
1449 ✗ m_matrixB[i] = new double[m_dimRomBasis];
1450 }
1451 }
1452 // In order to compute \int u_e = 0 :
1453 // m_sumG = [1 ... 1] * G
1454 ✗ if (m_sumG == nullptr) {
1455 ✗ m_sumG = new double[m_dimCollocationPts];
1456 ✗ for (int i=0; i<m_dimCollocationPts; i++) {
1457 ✗ m_sumG[i] = 0.;
1458 ✗ for (int j=0; j<m_dimCollocationPts; j++) {
1459 ✗ m_sumG[i] += m_massDeim[j][i];
1460 }
1461 }
1462 }
1463 }
1464 }
1465
1466 ✗ void EigenProblemALPDeim::setHeteroTauClose(std::vector<double>& valueTauClose) {
1467
1468 ✗ if (m_coefTauClose == nullptr)
1469 ✗ m_coefTauClose = new double[m_dimCollocationPts];
1470 // m_coefFhNs = P(valuef0)
1471 ✗ for (int j=0; j<m_dimCollocationPts; j++) {
1472 ✗ m_coefTauClose[j] = valueTauClose[m_collocationNode[j]-1];
1473 }
1474
1475 }
1476
1477 ✗ void EigenProblemALPDeim::setHeteroTauOut(std::vector<double>& valueTauOut) {
1478
1479 ✗ if (m_coefTauOut == nullptr)
1480 ✗ m_coefTauOut = new double[m_dimCollocationPts];
1481 // m_coefFhNs = P(valuef0)
1482 ✗ for (int j=0; j<m_dimCollocationPts; j++) {
1483 ✗ m_coefTauOut[j] = valueTauOut[m_collocationNode[j]-1];
1484 }
1485
1486 }
1487
1488 ✗ void EigenProblemALPDeim::setFhNf0(std::vector<double>& valuef0) {
1489
1490 ✗ if (m_coefFhNs == nullptr)
1491 ✗ m_coefFhNs = new double[m_dimCollocationPts];
1492 // m_coefFhNs = P(valuef0)
1493 ✗ for (int j=0; j<m_dimCollocationPts; j++) {
1494 ✗ m_coefFhNs[j] = valuef0[m_collocationNode[j]-1];
1495 }
1496
1497 }
1498
1499 ✗ void EigenProblemALPDeim::setIapp(std::vector<double>& iApp, int& idS) {
1500 ✗ if (m_source == nullptr) {
1501 ✗ m_source = new double*[m_numSource];
1502 ✗ for (int i=0; i<m_numSource; i++) {
1503 ✗ m_source[i] = new double[m_dimCollocationPts];
1504 }
1505 }
1506
1507 // m_source = P(iApp)
1508 ✗ for (int j=0; j<m_dimCollocationPts; j++) {
1509 ✗ m_source[idS][j] = iApp[m_collocationNode[j]-1];
1510 }
1511 }
1512
1513 //Update functions
1514 ✗ void EigenProblemALPDeim::updateBasis() {
1515 ✗ if (FelisceParam::verbose() > 10 ) PetscPrintf(PETSC_COMM_WORLD,"\nEigenProblemALPDeim::updateBasis()\n");
1516
1517 ✗ double deltaPhi[m_dimCollocationPts][m_dimRomBasis];
1518 ✗ double dt = FelisceParam::instance().timeStep;
1519
1520 ✗ double** psiC = nullptr;
1521 ✗ double** res = nullptr;
1522 // if (m_useImprovedRec==true) => d\Phi/dt = \Phi*M + R
1523 // if improvedRecType==0 => R = \Psi*C
1524 // if improvedRecType==1 => R = (-L-lambda)^{-1} t
1525 ✗ if (m_useImprovedRec)
1526 {
1527 ✗ switch (m_improvedRecType) {
1528 ✗ case 0: {
1529 // 1) Compute \Psi = (I-Phi^T Phi G)W_\orth
1530 ✗ double** phiphiT = new double* [m_dimCollocationPts];
1531 ✗ for (int i=0; i<m_dimCollocationPts; i++)
1532 ✗ phiphiT[i] = new double[m_dimCollocationPts];
1533 ✗ for (int i=0; i<m_dimCollocationPts; i++) {
1534 ✗ for (int j=0; j<m_dimCollocationPts; j++) {
1535 ✗ phiphiT[i][j] = 0.;
1536 ✗ for (int k=0; k<m_dimRomBasis; k++) {
1537 ✗ phiphiT[i][j] += m_basisDeim[i][k] * m_basisDeim[j][k];
1538 }
1539 }
1540 }
1541 ✗ for (int i=0; i<m_dimCollocationPts; i++) {
1542 ✗ double* vec = new double[m_dimCollocationPts];
1543 ✗ for (int k=0; k<m_dimCollocationPts; k++) {
1544 ✗ if (k == i)
1545 ✗ vec[k] = 1.;
1546 else
1547 ✗ vec[k] = 0.;
1548 ✗ for (int l=0; l<m_dimCollocationPts; l++) {
1549 ✗ vec[k] -= phiphiT[i][l] * m_massDeim[l][k];
1550 }
1551 }
1552 ✗ for (int j=0; j<m_dimOrthComp; j++) {
1553 ✗ m_orthCompDeim[i][j] = 0.;
1554 ✗ for (int k=0; k<m_dimCollocationPts; k++) {
1555 ✗ m_orthCompDeim[i][j] += vec[k] * m_basisDeimColloc[k][j+m_dimRomBasis];
1556 }
1557 }
1558 ✗ delete [] vec;
1559 }
1560 ✗ for (int i=0; i<m_dimCollocationPts; i++)
1561 ✗ delete [] phiphiT[i];
1562 ✗ delete [] phiphiT;
1563
1564 // 2) Compute "C" and \Psi*C
1565 ✗ psiC = new double* [m_dimCollocationPts];
1566 ✗ for (int i=0; i<m_dimCollocationPts; i++)
1567 ✗ psiC[i] = new double[m_dimRomBasis];
1568 ✗ double* vectorC = new double [m_dimOrthComp];
1569 ✗ double** matrixHtmp = new double* [m_dimCollocationPts];
1570 ✗ for (int i=0; i<m_dimCollocationPts; i++)
1571 ✗ matrixHtmp[i] = new double[m_dimOrthComp];
1572
1573 // H_tmp = K \Psi - \chi G (hat{u}*\Psi)
1574 ✗ for (int j=0; j<m_dimCollocationPts; j++) {
1575 ✗ for (int k=0; k<m_dimOrthComp; k++) {
1576 ✗ matrixHtmp[j][k] = 0.;
1577 ✗ for (int l=0; l<m_dimCollocationPts; l++) {
1578 ✗ matrixHtmp[j][k] += (m_stiffDeim[j][l] - m_coefChi * m_massDeim[j][l] * m_hatU[l] ) * m_orthCompDeim[l][k];
1579 }
1580 }
1581 }
1582
1583 // solve (H_const+h_i) c_i = \chi \Psi^T G f*\Phi \forall i
1584 ✗ for (int i=0; i<m_dimRomBasis; i++) {
1585 ✗ double checkOrth = 0.;
1586 ✗ double* rhsC = new double[m_dimOrthComp];
1587 ✗ for (int j=0; j<m_dimOrthComp; j++) {
1588 ✗ rhsC[j] = m_coefChi * m_thetaOrth[j][i];
1589 ✗ vectorC[j] = 1.;
1590 ✗ checkOrth += m_thetaOrth[j][i]*m_thetaOrth[j][i];
1591 }
1592
1593 ✗ if (FelisceParam::verbose()) {
1594 ✗ std::cout << "checkOrth = " << checkOrth << std::endl;
1595 }
1596 ✗ if (checkOrth > 1.e-4)
1597 {
1598 ✗ double** matrixH = new double* [m_dimOrthComp];
1599 ✗ for (int j=0; j<m_dimOrthComp; j++)
1600 ✗ matrixH[j] = new double[m_dimOrthComp];
1601
1602 // H = \Psi^T ( H_tmp + h_i )
1603 // h_i = - \lambda_i G \Psi
1604
1605 ✗ for (int j=0; j<m_dimOrthComp; j++) {
1606 ✗ double* tmpVec = new double[m_dimCollocationPts];
1607 ✗ for (int k=0; k<m_dimCollocationPts; k++) {
1608 ✗ tmpVec[k] = 0.;
1609 ✗ for (int l=0; l<m_dimCollocationPts; l++) {
1610 ✗ tmpVec[k] += m_eigenValue[i] * m_massDeim[k][l] * m_orthCompDeim[l][j];
1611 }
1612 }
1613 ✗ for (int k=0; k<m_dimOrthComp; k++) {
1614 ✗ matrixH[k][j] = 0.;
1615 ✗ for (int l=0; l<m_dimCollocationPts; l++) {
1616 ✗ matrixH[k][j] += m_orthCompDeim[l][k] * ( matrixHtmp[l][j] - tmpVec[l]);
1617 }
1618 }
1619 ✗ delete [] tmpVec;
1620 }
1621
1622 ✗ ConjGrad(matrixH, rhsC, vectorC, 1.e-7, 1000, m_dimOrthComp, true);
1623 ✗ for (int j=0; j<m_dimCollocationPts; j++) {
1624 ✗ psiC[j][i] = 0.;
1625 ✗ for (int k=0; k<m_dimOrthComp; k++) {
1626 ✗ psiC[j][i] += m_orthCompDeim[j][k]*vectorC[k];
1627 }
1628 }
1629 ✗ for (int j=0; j<m_dimOrthComp; j++)
1630 ✗ delete [] matrixH[j];
1631 ✗ delete [] matrixH;
1632 }
1633 else {
1634 ✗ for (int j=0; j<m_dimCollocationPts; j++) {
1635 ✗ psiC[j][i] = 0.;
1636 }
1637 }
1638 ✗ delete [] rhsC;
1639 }
1640
1641 ✗ delete [] vectorC;
1642 ✗ for (int i=0; i<m_dimOrthComp; i++)
1643 ✗ delete [] matrixHtmp[i];
1644 ✗ delete [] matrixHtmp;
1645
1646 ✗ break;
1647 }
1648 ✗ case 1: {
1649 ✗ res = new double* [m_dimRomBasis];
1650 ✗ for (int i=0; i<m_dimRomBasis; i++) {
1651 ✗ res[i] = new double[m_dimCollocationPts];
1652 }
1653 // 1) Compute H = K - \chi diag(\hat{u}^1/2) G diag(\hat{u}^1/2) - G \lambda_i = Htmp - G \lambda_i
1654 ✗ double** matrixH = new double* [m_dimCollocationPts];
1655 ✗ double** matrixHtmp = new double* [m_dimCollocationPts];
1656 ✗ for (int j=0; j<m_dimCollocationPts; j++) {
1657 ✗ matrixH[j] = new double[m_dimCollocationPts];
1658 ✗ matrixHtmp[j] = new double[m_dimCollocationPts];
1659 }
1660 // Htmp = K - 1/2 (diag(\hat{U}^{1/2}) G diag(\hat{U}^{1/2})
1661 ✗ for (int i=0; i<m_dimCollocationPts; i++) {
1662 ✗ for (int j=0; j<m_dimCollocationPts; j++) {
1663 ✗ matrixHtmp[i][j] = m_stiffDeim[i][j];
1664 ✗ if (FelisceParam::instance().optimizePotential) {
1665 ✗ matrixHtmp[i][j] -= m_coefChi * (m_potential[i]/std::fabs(m_potential[i])) * std::sqrt(std::fabs(m_potential[i])) * m_massDeim[i][j] * std::sqrt(std::fabs(m_potential[j]));
1666 }
1667 else {
1668 ✗ matrixHtmp[i][j] -= m_coefChi * (m_hatU[i]/std::fabs(m_hatU[i])) * std::sqrt(std::fabs(m_hatU[i])) * m_massDeim[i][j] * std::sqrt(std::fabs(m_hatU[j]));
1669 }
1670 }
1671 }
1672 // 2) Compute t_i = \chi G (\hat{f}*\phi_i - \sum_j \theta_ji \phi_j)
1673 ✗ double* vectorT = new double[m_dimCollocationPts];
1674 // H*r_i = t_i
1675 ✗ for (int i=0; i<m_dimRomBasis; i++) {
1676 // H = Htmp - G \lambda_i
1677 ✗ for (int j=0; j<m_dimCollocationPts; j++) {
1678 ✗ for (int k=0; k<m_dimCollocationPts; k++) {
1679 ✗ matrixH[j][k] = matrixHtmp[j][k] - m_eigenValue[i]*m_massDeim[j][k];
1680 }
1681 }
1682 // tmpVec = \sum_j \theta_ji \phi_j
1683 ✗ double tmpVec[m_dimCollocationPts];
1684 ✗ for (int k=0; k<m_dimCollocationPts; k++) {
1685 ✗ tmpVec[k] = 0.;
1686 ✗ for (int j=0; j<m_dimRomBasis; j++) {
1687 ✗ tmpVec[k] += m_theta[j][i]*m_basisDeim[k][j];
1688 }
1689 }
1690 // tmpVec2 = \hat{f} * \phi_i - tmpVec
1691 ✗ double tmpVec2[m_dimCollocationPts];
1692 ✗ for (int k=0; k<m_dimCollocationPts; k++) {
1693 ✗ tmpVec2[k] = m_hatF[k]*m_basisDeim[k][i] - tmpVec[k];
1694 }
1695 // t_i = \chi G (\hat{f} * \phi_i - \sum_j \theta_ji \phi_j) = \chi G (\hat{f} * \phi_i - tmpVec) = \chi G tmpVec2
1696 ✗ for (int j=0; j<m_dimCollocationPts; j++) {
1697 ✗ vectorT[j] = 0.;
1698 ✗ for (int k=0; k<m_dimCollocationPts; k++) {
1699 ✗ vectorT[j] += m_coefChi * m_massDeim[j][k] * tmpVec2[k];
1700 }
1701 }
1702
1703
1704 ✗ double Ht[m_dimCollocationPts];
1705 ✗ double tTHt = 0.;
1706 ✗ double tTt = 0.;
1707 ✗ for (int j=0; j<m_dimCollocationPts; j++) {
1708 ✗ Ht[j] = 0.;
1709 ✗ for (int k=0; k<m_dimCollocationPts; k++) {
1710 ✗ Ht[j] += matrixH[j][k]*vectorT[k];
1711 }
1712 }
1713 ✗ for (int k=0; k<m_dimCollocationPts; k++) {
1714 ✗ tTHt += vectorT[k]*Ht[k];
1715 ✗ tTt += vectorT[k]*vectorT[k];
1716 }
1717 double alpha;
1718 ✗ alpha = tTt/tTHt;
1719 // std::cout << "alpha " << alpha << std::endl;
1720 // double resT = 0.;
1721 ✗ double resNorm = 0.;
1722 ✗ for (int j=0; j<m_dimCollocationPts; j++) {
1723 // resT += vectorT[j]*vectorT[j];
1724 ✗ resNorm += (vectorT[j]-alpha*Ht[j])*(vectorT[j]-alpha*Ht[j]);
1725 }
1726 // std::cout << "res = " << resT << " - " << resNorm << std::endl;
1727 ✗ if (resNorm > 1.e-3) {
1728 ✗ for (int j=0; j<m_dimCollocationPts; j++) {
1729 ✗ res[i][j] = alpha*vectorT[j];
1730 }
1731 ✗ ConjGrad(matrixH, vectorT, res[i], 1.e-2, 10, m_dimCollocationPts, true);
1732 }
1733 else {
1734 ✗ for (int j=0; j<m_dimCollocationPts; j++) {
1735 ✗ res[i][j] = 0.;
1736 }
1737 }
1738
1739 }
1740
1741 ✗ delete [] vectorT;
1742 ✗ for (int i=0; i<m_dimCollocationPts; i++) {
1743 ✗ delete [] matrixH[i];
1744 ✗ delete [] matrixHtmp[i];
1745 }
1746 ✗ delete [] matrixH;
1747 ✗ delete [] matrixHtmp;
1748
1749 ✗ break;
1750 }
1751 ✗ default:
1752 ✗ break;
1753 }
1754 }
1755
1756 ✗ for (int i=0; i<m_dimCollocationPts; i++) {
1757 ✗ for (int j=0; j<m_dimRomBasis; j++) {
1758 ✗ deltaPhi[i][j] = 0.;
1759 ✗ for (int k=0; k<m_dimRomBasis; k++) {
1760 ✗ deltaPhi[i][j] += m_basisDeim[i][k]*m_matrixM[j][k];
1761 }
1762 // 3) d\Phi += \Psi*C
1763 ✗ if(m_useImprovedRec) {
1764 ✗ if (m_improvedRecType==0) {
1765 ✗ deltaPhi[i][j] += psiC[i][j];
1766 }
1767 ✗ else if (m_improvedRecType==1) {
1768 ✗ deltaPhi[i][j] += res[j][i];
1769 }
1770 }
1771 }
1772 }
1773 ✗ if(m_useImprovedRec) {
1774 ✗ if (m_improvedRecType==0) {
1775 ✗ for (int i=0; i<m_dimCollocationPts; i++)
1776 ✗ delete [] psiC[i];
1777 ✗ delete [] psiC;
1778 }
1779 ✗ else if (m_improvedRecType==1) {
1780 ✗ for (int i=0; i<m_dimRomBasis; i++)
1781 ✗ delete [] res[i];
1782 ✗ delete [] res;
1783 }
1784 }
1785
1786 ✗ for (int i=0; i<m_dimCollocationPts; i++)
1787 ✗ for (int j=0; j<m_dimRomBasis; j++)
1788 ✗ m_basisDeim[i][j] += dt*deltaPhi[i][j];
1789
1790
1791 ✗ if (m_resPoint) {
1792 ✗ for (int i=0; i<m_dimCollocationPts; i++) {
1793 ✗ m_resU[i] = dt*m_hatF[i];
1794 ✗ for (int j=0; j<m_dimRomBasis; j++) {
1795 ✗ m_resU[i] -= dt*deltaPhi[i][j]*m_beta[j];
1796 }
1797 }
1798 }
1799
1800
1801 double** matrixR;
1802 ✗ matrixR = new double* [m_dimRomBasis];
1803 ✗ for (int i=0; i< m_dimRomBasis; i++) {
1804 ✗ matrixR[i] = new double [m_dimRomBasis];
1805 }
1806
1807 ✗ MGS(m_basisDeim, matrixR, m_dimRomBasis);
1808
1809 ✗ for (int i=0; i< m_dimRomBasis; i++) {
1810 ✗ delete [] matrixR[i];
1811 }
1812 ✗ delete [] matrixR;
1813
1814 }
1815
1816 ✗ void EigenProblemALPDeim::setPotential()
1817 {
1818 ✗ if (FelisceParam::verbose() > 10 ) PetscPrintf(PETSC_COMM_WORLD,"\nEigenProblemALPDeim::setPotential()\n");
1819 ✗ int Nm = FelisceParam::instance().numApproxMode;
1820 ✗ int nCutOff = FelisceParam::instance().nCutOff;
1821
1822 ✗ m_initPot.duplicateFrom(m_basis[0]);
1823 ✗ m_initPot.copyFrom(m_U_0);
1824
1825 ✗ int sizePb = -1;
1826 ✗ std::vector<PetscVector> lagMult;
1827 ✗ std::vector<PetscVector> errorRep;
1828 ✗ if (m_problem < 2)
1829 { // FKPP or Monodomain
1830 ✗ sizePb = 1;
1831 }
1832 ✗ else if (m_problem == 2)
1833 { // Bidomain
1834 ✗ sizePb = 2;
1835 }
1836 else
1837 ✗ FEL_ERROR("Problem not found.");
1838 ✗ lagMult.resize(sizePb);
1839 ✗ errorRep.resize(sizePb);
1840
1841 ✗ for (int i=0; i<sizePb; i++) {
1842 ✗ lagMult[i].duplicateFrom(m_basis[0]);
1843 ✗ errorRep[i].duplicateFrom(m_basis[0]);
1844 }
1845
1846 ✗ PetscVector tmpVec;
1847 ✗ tmpVec.duplicateFrom(m_basis[0]);
1848
1849 ✗ double errorU = 1.;
1850 ✗ double errorUe = 1.;
1851
1852 ✗ std::vector<double*> coeffBeta;
1853 ✗ coeffBeta.resize(sizePb);
1854
1855 //EigenProblemALP::
1856 ✗ if (m_beta == nullptr) {
1857 ✗ m_beta = new double[Nm];
1858 ✗ projectOnBasis(m_U_0,m_beta,true,Nm);
1859 }
1860 ✗ m_solutionInitialized = true;
1861 ✗ coeffBeta[0] = m_beta;
1862
1863 ✗ errorRep[0].copyFrom(m_U_0);
1864 ✗ for (int i=0; i<Nm; i++) {
1865 ✗ double coef = -m_beta[i];
1866 ✗ errorRep[0].axpy(coef,m_basis[i]);
1867 }
1868 ✗ mult(m_Matrix[m_idG],errorRep[0],tmpVec);
1869 ✗ dot(errorRep[0],tmpVec,&errorU);
1870 // errorU = std::sqrt(errorU)/Nm;
1871
1872 ✗ if (m_problem == 2) {
1873 //EigenProblemALP::
1874 ✗ if (m_xi == nullptr) {
1875 ✗ m_xi = new double[Nm];
1876 ✗ projectOnBasis(m_Ue_0,m_xi,true,Nm);
1877 }
1878 ✗ coeffBeta[1] = m_xi;
1879 ✗ errorRep[1].copyFrom(m_Ue_0);
1880 ✗ for (int i=0; i<Nm; i++) {
1881 ✗ double coef = -m_xi[i];
1882 ✗ errorRep[1].axpy(coef,m_basis[i]);
1883 }
1884 ✗ mult(m_Matrix[m_idG],errorRep[1],tmpVec);
1885 ✗ dot(errorRep[1],tmpVec,&errorUe);
1886 // errorUe = std::sqrt(errorUe)/Nm;
1887 }
1888
1889 ✗ for (int i=0; i<sizePb; i++) {
1890 ✗ lagMult[i].copyFrom(errorRep[i]);
1891 ✗ lagMult[i].scale(-1.);
1892 }
1893
1894 ✗ PetscVector dLdU;
1895 ✗ dLdU.duplicateFrom(m_basis[0]);
1896
1897 ✗ PetscVector phi_ij;
1898 ✗ phi_ij.duplicateFrom(m_basis[0]);
1899
1900 ✗ PetscVector dUphi;
1901 ✗ dUphi.duplicateFrom(m_basis[0]);
1902
1903 ✗ PetscVector dU;
1904 ✗ dU.duplicateFrom(m_basis[0]);
1905
1906 ✗ PetscVector dUphiTG;
1907 ✗ dUphiTG.duplicateFrom(m_basis[0]);
1908 ✗ std::vector<PetscVector> lKeps;
1909 ✗ lKeps.resize(sizePb);
1910 ✗ for (int i=0; i<sizePb; i++) {
1911 ✗ lKeps[i].duplicateFrom(m_basis[0]);
1912 }
1913
1914 ✗ int numIter = 0;
1915 ✗ double step = 0.25; // 1.e-1;
1916 double k[2];
1917 ✗ k[0] = 50.;
1918 ✗ k[1] = 10.;
1919
1920 ✗ double** q = new double* [sizePb];
1921 ✗ for (int iPb=0; iPb<sizePb; iPb++) {
1922 ✗ q[iPb] = new double[nCutOff];
1923 }
1924 ✗ double** theta = new double* [nCutOff];
1925 ✗ for (int i=0; i<nCutOff; i++) {
1926 ✗ theta[i] = new double[nCutOff];
1927 }
1928
1929 ✗ PetscPrintf(PETSC_COMM_WORLD,"numIter = %i \n",numIter);
1930 ✗ PetscPrintf(PETSC_COMM_WORLD,"error = %e, %e \n",errorU,errorUe);
1931
1932 ✗ while ( (errorU > 1.e-08) && (errorUe > 1.e-08) && (numIter < 5000) )
1933 {
1934 // q[iPb]_j = <lagMult[iPb]+k*eps[iPb], phi_j>
1935 ✗ for (int iPb=0; iPb<sizePb; iPb++) {
1936 ✗ lKeps[iPb].copyFrom(lagMult[iPb]);
1937 ✗ lKeps[iPb].axpy(k[iPb],errorRep[iPb]);
1938 //EigenProblemALP::
1939 ✗ projectOnBasis(lKeps[iPb], q[iPb], true, nCutOff);
1940 }
1941
1942 // dL/dU = (U - u0)
1943 ✗ dLdU.copyFrom(m_initPot);
1944 ✗ dLdU.axpy(-1.,m_U_0);
1945 ✗ double reg = 1.;
1946 ✗ dLdU.scale(reg);
1947 // dL/dU -= \chi * \sum \sum ( sum_pb(beta_i q_j) / (lambda_j - lambda_i) phi_i phi_j
1948 ✗ for (int i=0; i<Nm; i++) {
1949 ✗ for (int j=0; j<nCutOff; j++) {
1950 ✗ if ( std::fabs(m_eigenValue[j] - m_eigenValue[i]) > 1.e-1 ) {
1951 ✗ pointwiseMult(phi_ij, m_basis[i], m_basis[j]);
1952 ✗ double coef = 0.;
1953 ✗ for (int iPb=0; iPb<sizePb; iPb++) {
1954 ✗ coef += coeffBeta[iPb][i]*q[iPb][j];
1955 }
1956 ✗ coef = - m_coefChi * coef/(m_eigenValue[j] - m_eigenValue[i]);
1957 ✗ dLdU.axpy(coef,phi_ij);
1958 }
1959 }
1960 }
1961
1962 // dU = -s dL/dU
1963 ✗ dU.set(0.);
1964 ✗ dU.axpy(-step,dLdU);
1965
1966 // theta_ij = <dU phi_i, phi_j>
1967 ✗ for (int i=0; i<nCutOff; i++) {
1968 ✗ pointwiseMult(dUphi, dU, m_basis[i]);
1969 ✗ mult(m_Matrix[m_idG], dUphi, dUphiTG);
1970 ✗ for (int j=0; j<=i; j++) {
1971 ✗ dot(dUphiTG, m_basis[j], &theta[i][j]);
1972 ✗ if (i != j)
1973 ✗ theta[j][i] = theta[i][j];
1974 }
1975 }
1976
1977 // Update phi
1978 ✗ for (int i=0; i<nCutOff; i++) {
1979 ✗ for (int j=0; j<nCutOff; j++) {
1980 ✗ if ( std::fabs(m_eigenValue[j] - m_eigenValue[i]) > 1.e-1 ) {
1981 ✗ double coef = m_coefChi * theta[i][j]/(m_eigenValue[j] - m_eigenValue[i]);
1982 ✗ m_basis[i].axpy(coef, m_basis[j]);
1983 }
1984 }
1985 }
1986 //EigenProblemALP::
1987 ✗ MGS(nCutOff);
1988
1989 // Update lambda
1990 ✗ for (int i=0; i<nCutOff; i++) {
1991 ✗ m_eigenValue[i] -= m_coefChi * theta[i][i];
1992 }
1993
1994 // Update beta
1995 //EigenProblemALP::
1996 ✗ projectOnBasis(m_U_0,m_beta,true,Nm);
1997 // Update xi
1998 ✗ if (m_problem == 2) {
1999 //EigenProblemALP::
2000 ✗ projectOnBasis(m_Ue_0,m_xi,true,Nm);
2001 }
2002
2003 // Compute error
2004 ✗ errorRep[0].copyFrom(m_U_0);
2005 ✗ for (int i=0; i<Nm; i++) {
2006 ✗ errorRep[0].axpy(-1.*m_beta[i],m_basis[i]);
2007 }
2008 ✗ mult(m_Matrix[m_idG],errorRep[0],tmpVec);
2009 ✗ dot(errorRep[0],tmpVec,&errorU);
2010 // errorU = std::sqrt(errorU)/Nm;
2011 ✗ if (m_problem == 2) {
2012 ✗ errorRep[1].copyFrom(m_Ue_0);
2013 ✗ for (int i=0; i<Nm; i++) {
2014 ✗ errorRep[1].axpy(-1.*m_xi[i],m_basis[i]);
2015 }
2016 ✗ mult(m_Matrix[m_idG],errorRep[1],tmpVec);
2017 ✗ dot(errorRep[1],tmpVec,&errorUe);
2018 // errorUe = std::sqrt(errorUe)/Nm;
2019 }
2020
2021 // Update Lagrangian Multiplicator
2022 ✗ for (int iPb=0; iPb<sizePb; iPb++) {
2023 ✗ lagMult[iPb].axpy(step, errorRep[iPb]);
2024 }
2025
2026 // Update initial potential
2027 ✗ m_initPot.axpy(1.,dU);
2028
2029 ✗ numIter++;
2030
2031 ✗ PetscPrintf(PETSC_COMM_WORLD,"numIter = %i \n",numIter);
2032 ✗ PetscPrintf(PETSC_COMM_WORLD,"error = %e, %e \n",errorU,errorUe);
2033 }
2034
2035 ✗ for (int iPb=0; iPb<sizePb; iPb++) {
2036 ✗ lagMult[iPb].destroy();
2037 ✗ errorRep[iPb].destroy();
2038 ✗ lKeps[iPb].destroy();
2039 }
2040 ✗ dLdU.destroy();
2041 ✗ phi_ij.destroy();
2042 ✗ dUphi.destroy();
2043 ✗ dU.destroy();
2044 ✗ dUphiTG.destroy();
2045 ✗ tmpVec.destroy();
2046
2047 ✗ for (int i=0; i<sizePb; i++) {
2048 ✗ delete [] q[i];
2049 }
2050 ✗ delete [] q;
2051 ✗ for (int i=0; i<nCutOff; i++) {
2052 ✗ delete [] theta[i];
2053 }
2054 ✗ delete [] theta;
2055
2056 // Write initial potential in Ensight file
2057 ✗ double* tmpSolToSave = nullptr;
2058 ✗ tmpSolToSave = new double[m_numDof];
2059 int rankProc;
2060 ✗ MPI_Comm_rank(m_petscComm,&rankProc);
2061 ✗ fromVecToDoubleStar(tmpSolToSave, m_initPot, rankProc, 1);
2062 ✗ writeEnsightVector(tmpSolToSave, 0, "initialPotential");
2063 ✗ delete [] tmpSolToSave;
2064
2065 ✗ std::string fileName = FelisceParam::instance().resultDir + "/eigenValue";
2066 ✗ viewALP(m_eigenValue,m_dimRomBasis,fileName);
2067
2068 ✗ checkBasis(Nm);
2069
2070 }
2071
2072 ✗ void EigenProblemALPDeim::updateBeta() {
2073 ✗ if (FelisceParam::verbose() > 10 ) PetscPrintf(PETSC_COMM_WORLD,"\nEigenProblemALPDeim::updateBeta()\n");
2074
2075 ✗ double dt = FelisceParam::instance().timeStep;
2076 ✗ double deltaBeta[m_dimRomBasis];
2077
2078 ✗ double Cm = FelisceParam::instance().Cm;
2079
2080 // beta^{n+1} = beta^n - dt * M^n * beta^n + dt * (\Phi^T G F^)^n
2081
2082 ✗ double PhiTG[m_dimRomBasis][m_dimCollocationPts];
2083 ✗ for (int i=0; i<m_dimRomBasis; i++) {
2084 ✗ for (int j=0; j<m_dimCollocationPts; j++) {
2085 ✗ PhiTG[i][j] = 0.;
2086 ✗ for (int k=0; k<m_dimCollocationPts; k++) {
2087 ✗ PhiTG[i][j] += m_basisDeim[k][i] * m_massDeim[k][j];
2088 }
2089 }
2090 }
2091
2092 // delta_beta^n
2093 ✗ for (int i=0; i<m_dimRomBasis; i++) {
2094 ✗ deltaBeta[i] = 0.;
2095
2096 // M^n * beta^n
2097 ✗ for (int j=0; j<m_dimRomBasis; j++) {
2098 ✗ deltaBeta[i] -= m_beta[j] * m_matrixM[j][i];
2099 }
2100 // (\Phi^T G F^)^n
2101 ✗ for (int k=0; k<m_dimCollocationPts; k++) {
2102 ✗ deltaBeta[i] += PhiTG[i][k] * m_hatF[k];
2103 }
2104 }
2105 // beta^{n+1} = beta^n - dt delta_beta^n
2106 ✗ for (int i=0; i<m_dimRomBasis; i++) {
2107 ✗ if (m_problem == 0) { // FKPP -> semi-implicit scheme
2108 // beta^{n+1} = beta^n - dt * M^n * beta^n + dt * (\Phi^T G F^)^n + dt (Cm - lambda) beta^n, F^ = - (chi + Cm) u^.2
2109 ✗ m_beta[i] = (1. - m_eigenValue[i]/2.*dt + Cm/2.*dt)/(1. + m_eigenValue[i]/2.*dt- Cm/2.*dt) * m_beta[i] + dt/(1. +m_eigenValue[i]/2.*dt - Cm/2.*dt) * deltaBeta[i];
2110 }
2111 else { // Explicit Euler
2112 ✗ m_beta[i] += dt * deltaBeta[i];
2113 }
2114 }
2115
2116 ✗ if (m_resPoint) {
2117 ✗ for (int i=0; i<m_dimCollocationPts; i++) {
2118 ✗ for (int j=0; j<m_dimRomBasis; j++) {
2119 ✗ m_resU[i] -= dt*m_basisDeim[i][j]*deltaBeta[j];
2120 }
2121 }
2122 }
2123
2124 ✗ if (m_problem == 2) {
2125 // If \sigma_E = \alpha \sigma_I
2126 // A \xi + B \beta = 0
2127 // A_ij = \sum_{h,l}^{N_p} (\alpha + 1)(\chi m_hatU[l] + lambda_i) \phi_j(x_h) G_hl \phi_i(x_l)
2128 // B_ij = \sum_{h,l}^{N_p} (\chi m_hatU[l] + lambda_i) \phi_j(x_h) G_hl \phi_i(x_h)
2129
2130 ✗ double alpha = FelisceParam::instance().extraTransvTensor / FelisceParam::instance().intraTransvTensor;
2131 ✗ if (FelisceParam::instance().testCase == 1) {
2132 ✗ double beta = FelisceParam::instance().extraFiberTensor / FelisceParam::instance().intraFiberTensor;
2133 ✗ if ( !felisce::Tools::almostEqual(alpha,beta,1.e-3) ) {
2134 ✗ FEL_ERROR("Intra and extra-cellular conductivity must have the same anysotropy");
2135 }
2136 }
2137
2138 ✗ for (int i=0; i<m_dimRomBasis; i++) {
2139 ✗ for (int j=0; j<m_dimRomBasis; j++) {
2140 ✗ m_matrixA[i][j] = .0;
2141 ✗ m_matrixB[i][j] = .0;
2142 ✗ for (int l=0; l<m_dimCollocationPts; l++) {
2143 ✗ if (FelisceParam::instance().optimizePotential) {
2144 ✗ m_matrixA[i][j] += (alpha+1)*(m_coefChi * m_potential[l] + m_eigenValue[i])*PhiTG[j][l]*m_basisDeim[l][i];
2145 ✗ m_matrixB[i][j] += (m_coefChi * m_potential[l] + m_eigenValue[i])*PhiTG[j][l]*m_basisDeim[l][i];
2146 }
2147 else {
2148 ✗ m_matrixA[i][j] += (alpha+1)*(m_coefChi * m_hatU[l] + m_eigenValue[i])*PhiTG[j][l]*m_basisDeim[l][i];
2149 ✗ m_matrixB[i][j] += (m_coefChi * m_hatU[l] + m_eigenValue[i])*PhiTG[j][l]*m_basisDeim[l][i];
2150 }
2151 }
2152 }
2153 ✗ m_matrixA[i][m_dimRomBasis] = 0.;
2154 ✗ for (int l=0; l<m_dimCollocationPts; l++) {
2155 ✗ m_matrixA[i][m_dimRomBasis] += m_sumG[l] * m_basisDeim[l][i];
2156 }
2157 ✗ m_matrixA[m_dimRomBasis][i] = m_matrixA[i][m_dimRomBasis];
2158 }
2159 ✗ m_matrixA[m_dimRomBasis][m_dimRomBasis] = 0.;
2160 ✗ double* newXi = new double[m_dimRomBasis+1];
2161 ✗ double* rhsXi = new double[m_dimRomBasis+1];
2162 ✗ for (int i=0; i<m_dimRomBasis; i++) {
2163 ✗ newXi[i] = m_xi[i];
2164 ✗ rhsXi[i] = 0.;
2165 ✗ for (int j=0; j<m_dimRomBasis; j++) {
2166 ✗ rhsXi[i] += - m_matrixB[i][j] * m_beta[j];
2167 }
2168 }
2169 ✗ newXi[m_dimRomBasis] = 0.;
2170 ✗ rhsXi[m_dimRomBasis] = 0.;
2171
2172 ✗ ConjGrad(m_matrixA, rhsXi, newXi, 1.e-09, 1000, m_dimRomBasis+1);
2173
2174 ✗ for (int i=0; i<m_dimRomBasis; i++) {
2175 ✗ m_xi[i] = newXi[i];
2176 }
2177 ✗ delete [] newXi;
2178 ✗ delete [] rhsXi;
2179 }
2180
2181 }
2182
2183 ✗ void EigenProblemALPDeim::ConjGrad(double** A, double* b, double* x, double tol, int maxIter, int n, bool sym) {
2184
2185 ✗ if (sym) {
2186 ✗ double* res = new double[n];
2187 ✗ double* p = new double[n];
2188 ✗ double* Ax = new double[n];
2189
2190 ✗ for (int i=0; i<n; i++) {
2191 ✗ Ax[i] = 0.;
2192 ✗ for (int j=0; j<n; j++) {
2193 ✗ Ax[i] += A[i][j]*x[j];
2194 }
2195 }
2196 ✗ for (int i=0; i<n; i++) {
2197 ✗ res[i] = b[i] - Ax[i];
2198 ✗ p[i] = res[i];
2199 }
2200 ✗ double res2 = 0.;
2201 ✗ double newRes2 = 0.;
2202 ✗ for (int i=0; i<n; i++) {
2203 ✗ res2 += res[i]*res[i];
2204 }
2205
2206 ✗ double* Ap = new double[n];
2207 double pAp;
2208 double alpha;
2209 double beta;
2210 ✗ bool isConv = false;
2211 ✗ int iter = 1;
2212
2213 ✗ if (res2 < tol) {
2214 ✗ isConv = true;
2215 ✗ if (FelisceParam::verbose()) {
2216 ✗ PetscPrintf(PETSC_COMM_WORLD,"Conjugate Gradient converged in %i iterations, residual: %1.10f \n",0,res2);
2217 }
2218 }
2219
2220 ✗ while (!isConv) {
2221 ✗ for (int i=0; i<n; i++) {
2222 ✗ Ap[i] = 0.;
2223 ✗ for (int j=0; j<n; j++) {
2224 ✗ Ap[i] += A[i][j]*p[j];
2225 }
2226 }
2227 ✗ pAp = 0.;
2228 ✗ for (int i=0; i<n; i++)
2229 ✗ pAp += p[i]*Ap[i];
2230 ✗ alpha = res2/pAp;
2231
2232 ✗ for (int i=0; i<n; i++) {
2233 ✗ x[i] += alpha*p[i];
2234 ✗ res[i] -= alpha*Ap[i];
2235 ✗ newRes2 += res[i]*res[i];
2236 }
2237
2238 ✗ if ( (newRes2 < tol) || (iter >= maxIter) ) {
2239 ✗ if (newRes2 < tol) {
2240 ✗ isConv = true;
2241 ✗ if (FelisceParam::verbose()) {
2242 ✗ PetscPrintf(PETSC_COMM_WORLD,"Conjugate Gradient converged in %i iterations, residual: %1.10f \n",iter,newRes2);
2243 }
2244 } else {
2245 ✗ PetscPrintf(PETSC_COMM_WORLD,"Conjugate Gradient: %i iterations, residual: %1.10f \n",iter,newRes2);
2246 ✗ FEL_WARNING("Conjugate Gradient did not converge.");
2247 }
2248 ✗ break;
2249 }
2250 else {
2251 ✗ beta = newRes2/res2;
2252 ✗ res2 = newRes2;
2253 ✗ newRes2 = 0.;
2254 ✗ for (int i=0; i<n; i++) {
2255 ✗ p[i] = res[i] + beta*p[i];
2256 }
2257 ✗ iter++;
2258 }
2259 }
2260
2261 ✗ delete [] res;
2262 ✗ delete [] p;
2263 ✗ delete [] Ap;
2264 ✗ delete [] Ax;
2265
2266 }
2267 else
2268 {
2269 ✗ double** AtA = new double*[n];
2270 ✗ for (int i=0; i<n; i++)
2271 ✗ AtA[i] = new double[n];
2272
2273 ✗ double* Atb = new double[n];
2274
2275 ✗ for (int i=0; i<n; i++) {
2276 ✗ Atb[i] = 0.;
2277 ✗ for (int j=0; j<n; j++) {
2278 ✗ Atb[i] += A[j][i]*b[j];
2279 ✗ AtA[i][j] = 0.;
2280 ✗ for (int k=0; k<n; k++) {
2281 ✗ AtA[i][j] += A[k][i]*A[k][j];
2282 }
2283 }
2284 }
2285
2286 ✗ double* res = new double[n];
2287 ✗ double* p = new double[n];
2288 ✗ double* Ax = new double[n];
2289
2290 ✗ for (int i=0; i<n; i++) {
2291 ✗ Ax[i] = 0.;
2292 ✗ for (int j=0; j<n; j++) {
2293 ✗ Ax[i] += AtA[i][j]*x[j];
2294 }
2295 }
2296 ✗ for (int i=0; i<n; i++) {
2297 ✗ res[i] = Atb[i] - Ax[i];
2298 ✗ p[i] = res[i];
2299 }
2300 ✗ double res2 = 0.;
2301 ✗ double newRes2 = 0.;
2302 ✗ for (int i=0; i<n; i++) {
2303 ✗ res2 += res[i]*res[i];
2304 }
2305
2306 ✗ double* Ap = new double[n];
2307 double pAp;
2308 double alpha;
2309 double beta;
2310 ✗ bool isConv = false;
2311 ✗ int iter = 1;
2312
2313 ✗ if (res2 < tol) {
2314 ✗ isConv = true;
2315 ✗ if (FelisceParam::verbose()) {
2316 ✗ PetscPrintf(PETSC_COMM_WORLD,"Conjugate Gradient converged in %i iterations, residual: %1.10f \n",0,res2);
2317 }
2318 }
2319
2320 ✗ while (!isConv) {
2321 ✗ for (int i=0; i<n; i++) {
2322 ✗ Ap[i] = 0.;
2323 ✗ for (int j=0; j<n; j++) {
2324 ✗ Ap[i] += AtA[i][j]*p[j];
2325 }
2326 }
2327 ✗ pAp = 0.;
2328 ✗ for (int i=0; i<n; i++)
2329 ✗ pAp += p[i]*Ap[i];
2330 ✗ alpha = res2/pAp;
2331
2332 ✗ for (int i=0; i<n; i++) {
2333 ✗ x[i] += alpha*p[i];
2334 ✗ res[i] -= alpha*Ap[i];
2335 ✗ newRes2 += res[i]*res[i];
2336 }
2337
2338 ✗ if ( (newRes2 < tol) || (iter >= maxIter) ) {
2339 ✗ if (newRes2 < tol) {
2340 ✗ isConv = true;
2341 ✗ if (FelisceParam::verbose()) {
2342 ✗ PetscPrintf(PETSC_COMM_WORLD,"Conjugate Gradient converged in %i iterations, residual: %1.10f \n",iter,newRes2);
2343 }
2344 } else {
2345 ✗ PetscPrintf(PETSC_COMM_WORLD,"Conjugate Gradient: %i iterations, residual: %1.10f \n",iter,newRes2);
2346 ✗ FEL_ERROR("Conjugate Gradient did not converge.");
2347 }
2348 ✗ break;
2349 } else {
2350 ✗ beta = newRes2/res2;
2351 ✗ res2 = newRes2;
2352 ✗ newRes2 = 0.;
2353 ✗ for (int i=0; i<n; i++) {
2354 ✗ p[i] = res[i] + beta*p[i];
2355 }
2356 ✗ iter++;
2357 }
2358 }
2359
2360 ✗ delete [] res;
2361 ✗ delete [] p;
2362 ✗ delete [] Ap;
2363 ✗ delete [] Ax;
2364
2365 ✗ for (int i=0; i<n; i++)
2366 ✗ delete AtA[i];
2367 ✗ delete [] AtA;
2368 ✗ delete [] Atb;
2369 }
2370
2371 }
2372
2373 ✗ void EigenProblemALPDeim::updateEigenvalue() {
2374 ✗ if (FelisceParam::verbose() > 10 ) PetscPrintf(PETSC_COMM_WORLD,"\nEigenProblemALPDeim::updateEigenvalue()\n");
2375
2376 ✗ double dt = FelisceParam::instance().timeStep;
2377 // lambda_i^{n+1} = lambda_i^n - dt * m_coefChi * m_theta_ii
2378 ✗ for (int i=0; i<m_dimRomBasis; i++) {
2379 ✗ m_eigenValue[i] -= dt * m_coefChi * m_theta[i][i];
2380 }
2381
2382 }
2383
2384 ✗ void EigenProblemALPDeim::checkBasis(int size) {
2385 ✗ if (FelisceParam::verbose() > 10 ) PetscPrintf(PETSC_COMM_WORLD,"\nEigenProblemALPDeim::checkBasis()\n");
2386
2387 ✗ if (size == 0)
2388 ✗ size = m_dimRomBasis;
2389
2390 int rankProc;
2391 ✗ MPI_Comm_rank(m_petscComm,&rankProc);
2392 ✗ PetscVector tmpSol;
2393 ✗ tmpSol.duplicateFrom(m_U_0);
2394
2395 ✗ double* tmpSolToSave = nullptr;
2396 ✗ tmpSolToSave = new double[m_numDof];
2397
2398 // Check convergence on m_U_0
2399 ✗ double convergence[size];
2400 ✗ for (int iBasis=0; iBasis < size; iBasis++) {
2401 ✗ double tmpVec[iBasis+1];
2402 ✗ for (int jBasis=0; jBasis<iBasis+1; jBasis++) {
2403 ✗ tmpVec[jBasis] = m_beta[jBasis];
2404 }
2405
2406 //EigenProblemALP::
2407 ✗ projectOnDof(tmpVec,tmpSol,iBasis+1);
2408
2409 ✗ checkSolution(tmpSol,m_U_0,convergence[iBasis]);
2410
2411 ✗ fromVecToDoubleStar(tmpSolToSave, tmpSol, rankProc, 1);
2412 ✗ if (m_problem == 0) {
2413 ✗ writeEnsightVector(tmpSolToSave, iBasis, "uConv");
2414 }
2415 ✗ if ( (m_problem == 1) || (m_problem == 2) ) {
2416 ✗ writeEnsightVector(tmpSolToSave, iBasis, "potTransMembConv");
2417 }
2418
2419 ✗ tmpSol.axpy(-1.,m_U_0);
2420 ✗ tmpSol.abs();
2421 double maxU0;
2422 ✗ m_U_0.max(&maxU0);
2423 double minU0;
2424 ✗ m_U_0.min(&minU0);
2425 ✗ double deltaU0 = maxU0-minU0;
2426 ✗ tmpSol.scale(1./deltaU0);
2427 ✗ fromVecToDoubleStar(tmpSolToSave, tmpSol, rankProc, 1);
2428 ✗ if (m_problem == 0) {
2429 ✗ writeEnsightVector(tmpSolToSave, iBasis, "uError");
2430 }
2431 ✗ if ( (m_problem == 1) || (m_problem == 2) ) {
2432 ✗ writeEnsightVector(tmpSolToSave, iBasis, "potTransMembError");
2433 }
2434
2435 ✗ PetscPrintf(PETSC_COMM_WORLD, "convergence[%d] = %e \n", iBasis, convergence[iBasis]);
2436
2437 }
2438
2439
2440 ✗ std::string fileName = FelisceParam::instance().resultDir + "/convergence";
2441 ✗ viewALP(convergence,size,fileName);
2442
2443 // Check convergence on m_Ue_0
2444 ✗ if (m_problem == 2) {
2445 ✗ double convergenceUe[size];
2446 ✗ for (int iBasis=0; iBasis < size; iBasis++) {
2447 ✗ double tmpVec[iBasis+1];
2448 ✗ for (int jBasis=0; jBasis<iBasis+1; jBasis++) {
2449 ✗ tmpVec[jBasis] = m_xi[jBasis];
2450 }
2451
2452 //EigenProblemALP::
2453 ✗ projectOnDof(tmpVec,tmpSol,iBasis+1);
2454
2455 ✗ checkSolution(tmpSol,m_Ue_0,convergenceUe[iBasis]);
2456
2457 ✗ fromVecToDoubleStar(tmpSolToSave, tmpSol, rankProc, 1);
2458 ✗ writeEnsightVector(tmpSolToSave, iBasis, "potExtraCellConv");
2459
2460 ✗ tmpSol.axpy(-1.,m_Ue_0);
2461 ✗ tmpSol.abs();
2462 double maxU0;
2463 ✗ m_Ue_0.max(&maxU0);
2464 double minU0;
2465 ✗ m_Ue_0.min(&minU0);
2466 ✗ double deltaU0 = maxU0-minU0;
2467 ✗ tmpSol.scale(1./deltaU0);
2468 ✗ fromVecToDoubleStar(tmpSolToSave, tmpSol, rankProc, 1);
2469 ✗ writeEnsightVector(tmpSolToSave, iBasis, "potExtraCellError");
2470
2471 ✗ PetscPrintf(PETSC_COMM_WORLD, "convergenceUe[%d] = %e \n", iBasis, convergenceUe[iBasis]);
2472
2473 }
2474
2475 ✗ fileName = FelisceParam::instance().resultDir + "/convergenceUe";
2476 ✗ viewALP(convergenceUe,size,fileName);
2477 }
2478
2479 ✗ std::vector<std::string> varNames;
2480 ✗ if (m_problem == 0) {
2481 ✗ varNames.resize(2);
2482 ✗ varNames[0] = "uConv";
2483 ✗ varNames[1] = "uError";
2484 }
2485 ✗ else if (m_problem == 1) {
2486 ✗ varNames.resize(2);
2487 ✗ varNames[0] = "potTransMembConv";
2488 ✗ varNames[1] = "potTransMembError";
2489 }
2490 ✗ else if (m_problem == 2) {
2491 ✗ varNames.resize(4);
2492 ✗ varNames[0] = "potTransMembConv";
2493 ✗ varNames[1] = "potExtraCellConv";
2494 ✗ varNames[2] = "potTransMembError";
2495 ✗ varNames[3] = "potExtraCellError";
2496 }
2497 ✗ if (rankProc == 0)
2498 ✗ writeEnsightCase(size, 1.,varNames);
2499
2500 ✗ delete [] tmpSolToSave;
2501 ✗ tmpSol.destroy();
2502 }
2503
2504 ✗ void EigenProblemALPDeim::checkBasisCollocation() {
2505 ✗ if (FelisceParam::verbose() > 10 ) PetscPrintf(PETSC_COMM_WORLD,"\nEigenProblemALPDeim::checkBasisCollocation()\n");
2506
2507 ✗ if (!m_solutionInitialized) {
2508 ✗ readCollocationPts();
2509 ✗ initializeROM();
2510 ✗ computeMassDeim();
2511 ✗ initializeSolution();
2512 ✗ computeProjectionPiV();
2513 }
2514
2515 int rankProc;
2516 ✗ MPI_Comm_rank(m_petscComm,&rankProc);
2517 ✗ PetscVector tmpSol;
2518 ✗ tmpSol.duplicateFrom(m_U_0);
2519
2520 ✗ double* tmpSolToSave = nullptr;
2521 ✗ tmpSolToSave = new double[m_numDof];
2522
2523 // Check convergence on m_U_0
2524 ✗ double convergence[m_dimRomBasis];
2525 ✗ for (int iBasis=1; iBasis < m_dimRomBasis; iBasis++) {
2526 ✗ tmpSol.set(0.);
2527 ✗ double* tmpVec = new double[iBasis];
2528 ✗ for (int jBasis=0; jBasis<iBasis; jBasis++) {
2529 ✗ tmpVec[jBasis] = m_beta[jBasis];
2530 }
2531
2532 ✗ projectOnDof(tmpVec,m_hatU,iBasis);
2533
2534 ✗ for (int k=0; k<m_dimCollocationPts; k++) {
2535 ✗ tmpSol.axpy( m_hatU[k], m_projectorPiV[k]);
2536 }
2537
2538 ✗ checkSolution(tmpSol,m_U_0,convergence[iBasis]);
2539
2540 ✗ fromVecToDoubleStar(tmpSolToSave, tmpSol, rankProc, 1);
2541 ✗ if (m_problem == 0) {
2542 ✗ writeEnsightVector(tmpSolToSave, iBasis, "uConv");
2543 }
2544 ✗ if ( (m_problem == 1) || (m_problem == 2) ) {
2545 ✗ writeEnsightVector(tmpSolToSave, iBasis, "potTransMembConv");
2546 }
2547
2548 ✗ tmpSol.axpy(-1.,m_U_0);
2549 ✗ tmpSol.abs();
2550 double maxU0;
2551 ✗ m_U_0.max(&maxU0);
2552 double minU0;
2553 ✗ m_U_0.min(&minU0);
2554 ✗ double deltaU0 = maxU0-minU0;
2555 ✗ tmpSol.scale(1./deltaU0);
2556 ✗ fromVecToDoubleStar(tmpSolToSave, tmpSol, rankProc, 1);
2557 ✗ if (m_problem == 0) {
2558 ✗ writeEnsightVector(tmpSolToSave, iBasis, "uError");
2559 }
2560 ✗ if ( (m_problem == 1) || (m_problem == 2) ) {
2561 ✗ writeEnsightVector(tmpSolToSave, iBasis, "potTransMembError");
2562 }
2563
2564 ✗ PetscPrintf(PETSC_COMM_WORLD, "convergence[%d] = %e \n", iBasis, convergence[iBasis]);
2565
2566 ✗ delete [] tmpVec;
2567 }
2568
2569
2570 ✗ std::string fileName = FelisceParam::instance().resultDir + "/convergence";
2571 ✗ viewALP(convergence,m_dimRomBasis,fileName);
2572
2573 // Check convergence on m_W_0
2574 // if ( (m_problem == 1) || (m_problem == 2) ) {
2575 // double convergenceW[m_dimRomBasis];
2576 // for (int iBasis=1; iBasis < m_dimRomBasis; iBasis++) {
2577 // tmpSol.set(0.);
2578 // double* tmpVec = new double[iBasis];
2579 // for (int jBasis=0; jBasis<iBasis; jBasis++) {
2580 // tmpVec[jBasis] = m_mu[jBasis];
2581 // }
2582 //
2583 // projectOnDof(tmpVec,m_hatW,iBasis);
2584 //
2585 // for (int k=0; k<m_dimCollocationPts; k++) {
2586 // tmpSol.axpy( m_hatW[k], m_projectorPiV[k]);
2587 // }
2588 //
2589 // checkSolution(tmpSol,m_W_0,convergenceW[iBasis]);
2590 //
2591 // fromVecToDoubleStar(tmpSolToSave, tmpSol, rankProc, 1);
2592 // writeEnsightVector(tmpSolToSave, iBasis, "ionicVarConv");
2593 //
2594 // PetscPrintf(PETSC_COMM_WORLD, "convergenceW[%d] = %e \n", iBasis, convergenceW[iBasis]);
2595 //
2596 // delete [] tmpVec;
2597 // }
2598 // if (rankProc == 0)
2599 // writeEnsightCase(m_dimRomBasis, 1., "ionicVarConv");
2600 // fileName = FelisceParam::instance().resultDir + "/convergenceW";
2601 // viewALP(convergenceW,m_dimRomBasis,fileName);
2602 // }
2603
2604 // Check convergence on m_Ue_0
2605 ✗ if (m_problem == 2) {
2606 ✗ double convergenceUe[m_dimRomBasis];
2607 ✗ for (int iBasis=1; iBasis < m_dimRomBasis; iBasis++) {
2608 ✗ tmpSol.set(0.);
2609 ✗ double* tmpVec = new double[iBasis];
2610 ✗ for (int jBasis=0; jBasis<iBasis; jBasis++) {
2611 ✗ tmpVec[jBasis] = m_xi[jBasis];
2612 }
2613
2614 ✗ projectOnDof(tmpVec,m_hatUe,iBasis);
2615
2616 ✗ for (int k=0; k<m_dimCollocationPts; k++) {
2617 ✗ tmpSol.axpy( m_hatUe[k], m_projectorPiV[k]);
2618 }
2619
2620 ✗ checkSolution(tmpSol,m_Ue_0,convergenceUe[iBasis]);
2621
2622 ✗ fromVecToDoubleStar(tmpSolToSave, tmpSol, rankProc, 1);
2623 ✗ writeEnsightVector(tmpSolToSave, iBasis, "potExtraCellConv");
2624
2625 ✗ tmpSol.axpy(-1.,m_Ue_0);
2626 ✗ tmpSol.abs();
2627 double maxU0;
2628 ✗ m_Ue_0.max(&maxU0);
2629 double minU0;
2630 ✗ m_Ue_0.min(&minU0);
2631 ✗ double deltaU0 = maxU0-minU0;
2632 ✗ tmpSol.scale(1./deltaU0);
2633 ✗ fromVecToDoubleStar(tmpSolToSave, tmpSol, rankProc, 1);
2634 ✗ writeEnsightVector(tmpSolToSave, iBasis, "potExtraCellError");
2635
2636 ✗ PetscPrintf(PETSC_COMM_WORLD, "convergenceUe[%d] = %e \n", iBasis, convergenceUe[iBasis]);
2637
2638 ✗ delete [] tmpVec;
2639 }
2640
2641 ✗ fileName = FelisceParam::instance().resultDir + "/convergenceUe";
2642 ✗ viewALP(convergenceUe,m_dimRomBasis,fileName);
2643 }
2644
2645 ✗ std::vector<std::string> varNames;
2646 ✗ if (m_problem == 0) {
2647 ✗ varNames.resize(2);
2648 ✗ varNames[0] = "uConv";
2649 ✗ varNames[1] = "uError";
2650 }
2651 ✗ else if (m_problem == 1) {
2652 ✗ varNames.resize(2);
2653 ✗ varNames[0] = "potTransMembConv";
2654 ✗ varNames[1] = "potTransMembError";
2655 }
2656 ✗ else if (m_problem == 2) {
2657 ✗ varNames.resize(4);
2658 ✗ varNames[0] = "potTransMembConv";
2659 ✗ varNames[1] = "potExtraCellConv";
2660 ✗ varNames[2] = "potTransMembError";
2661 ✗ varNames[3] = "potExtraCellError";
2662 }
2663 ✗ if (rankProc == 0)
2664 ✗ writeEnsightCase(m_dimRomBasis, 1.,varNames);
2665
2666 ✗ delete [] tmpSolToSave;
2667 ✗ tmpSol.destroy();
2668 }
2669
2670
2671 }
2672