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: A. Collin |
13 |
|
|
// |
14 |
|
|
|
15 |
|
|
// System includes |
16 |
|
|
|
17 |
|
|
// External includes |
18 |
|
|
|
19 |
|
|
// Project includes |
20 |
|
|
#include "Model/bidomainModel.hpp" |
21 |
|
|
#include "Solver/CVGraphInterface.hpp" |
22 |
|
|
|
23 |
|
|
|
24 |
|
|
namespace felisce { |
25 |
|
✗ |
BidomainModel::BidomainModel(): |
26 |
|
|
Model(), |
27 |
|
✗ |
m_rom(nullptr), |
28 |
|
✗ |
m_verbose(0), |
29 |
|
✗ |
m_buildIonic(false), |
30 |
|
✗ |
m_buildCourt(false), |
31 |
|
✗ |
m_ionic(nullptr), |
32 |
|
✗ |
m_schaf(nullptr), |
33 |
|
✗ |
m_fhn(nullptr), |
34 |
|
✗ |
m_court(nullptr), |
35 |
|
✗ |
m_mv(nullptr), |
36 |
|
✗ |
m_paci(nullptr), |
37 |
|
✗ |
m_bcl(nullptr), |
38 |
|
✗ |
m_observations(nullptr), |
39 |
|
✗ |
m_extrapolateIsCreated(false) { |
40 |
|
✗ |
for(std::size_t i=0; i<m_linearProblemBidomain.size(); i++) { |
41 |
|
✗ |
m_linearProblemBidomain[i] = nullptr; |
42 |
|
|
} |
43 |
|
✗ |
m_name = "Bidomain"; |
44 |
|
|
} |
45 |
|
|
|
46 |
|
✗ |
BidomainModel::~BidomainModel() { |
47 |
|
✗ |
for(std::size_t i=0; i<m_linearProblemBidomain.size(); i++) { |
48 |
|
✗ |
m_linearProblemBidomain[i] = nullptr; |
49 |
|
|
} |
50 |
|
✗ |
if(FelisceParam::instance().useROM) { |
51 |
|
✗ |
delete m_rom; |
52 |
|
✗ |
m_rom=nullptr; |
53 |
|
|
} |
54 |
|
✗ |
if (m_buildIonic) { |
55 |
|
✗ |
delete m_ionic; |
56 |
|
✗ |
m_ionic=nullptr; |
57 |
|
|
} |
58 |
|
✗ |
if (m_buildCourt) { |
59 |
|
✗ |
delete m_court; |
60 |
|
✗ |
m_court = nullptr; |
61 |
|
|
} |
62 |
|
✗ |
if(m_observations) { |
63 |
|
✗ |
delete m_observations; |
64 |
|
✗ |
m_observations=nullptr; |
65 |
|
|
} |
66 |
|
|
|
67 |
|
✗ |
for (std::size_t i=0; i<m_W_0.size(); i++) { |
68 |
|
✗ |
m_W_0[i].destroy(); |
69 |
|
|
} |
70 |
|
✗ |
if ( (FelisceParam::instance().typeOfIonicModel == "courtAtriaSchafVent") || (FelisceParam::instance().typeOfIonicModel == "courtAtriaMvVent") ) { |
71 |
|
✗ |
m_m.destroy(); |
72 |
|
✗ |
m_h.destroy(); |
73 |
|
✗ |
m_j.destroy(); |
74 |
|
✗ |
m_ao.destroy(); |
75 |
|
✗ |
m_io.destroy(); |
76 |
|
✗ |
m_ua.destroy(); |
77 |
|
✗ |
m_ui.destroy(); |
78 |
|
✗ |
m_xr.destroy(); |
79 |
|
✗ |
m_xs.destroy(); |
80 |
|
✗ |
m_d.destroy(); |
81 |
|
✗ |
m_f.destroy(); |
82 |
|
✗ |
m_fca.destroy(); |
83 |
|
✗ |
m_urel.destroy(); |
84 |
|
✗ |
m_vrel.destroy(); |
85 |
|
✗ |
m_wrel.destroy(); |
86 |
|
✗ |
m_nai.destroy(); |
87 |
|
✗ |
m_nao.destroy(); |
88 |
|
✗ |
m_cao.destroy(); |
89 |
|
✗ |
m_ki.destroy(); |
90 |
|
✗ |
m_ko.destroy(); |
91 |
|
✗ |
m_cai.destroy(); |
92 |
|
✗ |
m_naiont.destroy(); |
93 |
|
✗ |
m_kiont.destroy(); |
94 |
|
✗ |
m_caiont.destroy(); |
95 |
|
✗ |
m_ileak.destroy(); |
96 |
|
✗ |
m_iup.destroy(); |
97 |
|
✗ |
m_itr.destroy(); |
98 |
|
✗ |
m_irel.destroy(); |
99 |
|
✗ |
m_cmdn.destroy(); |
100 |
|
✗ |
m_trpn.destroy(); |
101 |
|
✗ |
m_nsr.destroy(); |
102 |
|
✗ |
m_jsr.destroy(); |
103 |
|
✗ |
m_csqn.destroy(); |
104 |
|
|
} |
105 |
|
|
|
106 |
|
✗ |
if (m_extrapolateIsCreated) |
107 |
|
✗ |
m_extrapolate.destroy(); |
108 |
|
✗ |
delete m_iApp; |
109 |
|
|
|
110 |
|
✗ |
if (FelisceParam::instance().hasAppliedExteriorCurrent) |
111 |
|
✗ |
delete m_iAppExt; |
112 |
|
|
|
113 |
|
✗ |
if(FelisceParam::instance().stateFilter) { |
114 |
|
✗ |
m_luenbFlter.destroy(); |
115 |
|
✗ |
m_Observ.destroy(); |
116 |
|
|
} |
117 |
|
|
|
118 |
|
✗ |
if (FelisceParam::instance().restartSolution) { |
119 |
|
✗ |
if (FelisceParam::instance().orderBdfEdp > 0) |
120 |
|
✗ |
m_bdf_sol_n.destroy(); |
121 |
|
✗ |
if (FelisceParam::instance().orderBdfEdp > 1) |
122 |
|
✗ |
m_bdf_sol_n_1.destroy(); |
123 |
|
✗ |
if (FelisceParam::instance().orderBdfEdp > 2) |
124 |
|
✗ |
m_bdf_sol_n_2.destroy(); |
125 |
|
✗ |
if (FelisceParam::instance().orderBdfIonic > 0) { |
126 |
|
✗ |
for (std::size_t i=0; i<m_ionic_sol_n.size(); i++) { |
127 |
|
✗ |
m_ionic_sol_n[i].destroy(); |
128 |
|
|
} |
129 |
|
|
} |
130 |
|
✗ |
if (FelisceParam::instance().orderBdfIonic > 1) { |
131 |
|
✗ |
for (std::size_t i=0; i<m_ionic_sol_n_1.size(); i++) { |
132 |
|
✗ |
m_ionic_sol_n_1[i].destroy(); |
133 |
|
|
} |
134 |
|
|
} |
135 |
|
✗ |
if (FelisceParam::instance().orderBdfIonic > 2) { |
136 |
|
✗ |
for (std::size_t i=0; i<m_ionic_sol_n_2.size(); i++) { |
137 |
|
✗ |
m_ionic_sol_n_2[i].destroy(); |
138 |
|
|
} |
139 |
|
|
} |
140 |
|
✗ |
if ( (FelisceParam::instance().typeOfIonicModel == "courtAtriaSchafVent") || (FelisceParam::instance().typeOfIonicModel == "courtAtriaMvVent") ) { |
141 |
|
✗ |
for (int iCourt=0 ; iCourt<29; iCourt++) { |
142 |
|
✗ |
m_courtVarInit[iCourt].destroy(); |
143 |
|
|
} |
144 |
|
|
} |
145 |
|
|
} |
146 |
|
|
} |
147 |
|
|
|
148 |
|
✗ |
void BidomainModel::initializeDerivedLinearProblem() { |
149 |
|
✗ |
for (std::size_t i=0; i<m_linearProblem.size(); i++) { |
150 |
|
✗ |
m_linearProblemBidomain.push_back(static_cast<LinearProblemBidomain*>(m_linearProblem[i])); |
151 |
|
|
} |
152 |
|
|
} |
153 |
|
|
|
154 |
|
|
|
155 |
|
|
|
156 |
|
|
//! Define m_ionic and m_bdfEdp (used in linearProblem). |
157 |
|
✗ |
void BidomainModel::initializeDerivedModel() { |
158 |
|
✗ |
m_verbose = FelisceParam::verbose(); |
159 |
|
|
|
160 |
|
✗ |
if ( FelisceParam::instance().withCVG ){ |
161 |
|
✗ |
m_linearProblemBidomain[0]->initPetscVectors(); |
162 |
|
|
} |
163 |
|
|
|
164 |
|
|
//Allocate m_ionic. |
165 |
|
|
// Choose type of ionic model to be used. |
166 |
|
✗ |
if ((FelisceParam::instance().typeOfIonicModel == "schaf") || (FelisceParam::instance().typeOfIonicModel == "courtAtriaSchafVent")) { |
167 |
|
|
// Initialize Mitchell and Schaeffer model. |
168 |
|
✗ |
m_schaf = new SchafSolver(m_fstransient); |
169 |
|
✗ |
m_ionic = m_schaf; |
170 |
|
|
|
171 |
|
✗ |
if (FelisceParam::instance().typeOfIonicModel == "courtAtriaSchafVent") { |
172 |
|
|
// Initialize Courtemanche model |
173 |
|
✗ |
m_court = new CourtemancheSolver(m_fstransient); |
174 |
|
|
} |
175 |
|
✗ |
} else if (FelisceParam::instance().typeOfIonicModel == "schafRev") { |
176 |
|
|
// Initialize Mitchell and Schaeffer Revised model. |
177 |
|
✗ |
m_schaf = new SchafRevisedSolver(m_fstransient); |
178 |
|
✗ |
m_ionic = m_schaf; |
179 |
|
✗ |
} else if (FelisceParam::instance().typeOfIonicModel == "fhn") { |
180 |
|
|
// Initialize Fitzhugh-Nagumo model |
181 |
|
✗ |
m_fhn = new FhNSolver(m_fstransient); |
182 |
|
✗ |
m_ionic =m_fhn; |
183 |
|
✗ |
} else if ( (FelisceParam::instance().typeOfIonicModel == "MV") || (FelisceParam::instance().typeOfIonicModel == "courtAtriaMvVent") ) { |
184 |
|
|
// Initialize MV model |
185 |
|
✗ |
m_mv = new MVSolver(m_fstransient); |
186 |
|
✗ |
m_ionic = m_mv; |
187 |
|
|
|
188 |
|
✗ |
if (FelisceParam::instance().typeOfIonicModel == "courtAtriaMvVent") { |
189 |
|
|
// Initialize Courtemanche model |
190 |
|
✗ |
m_court = new CourtemancheSolver(m_fstransient); |
191 |
|
|
} |
192 |
|
✗ |
} else if (FelisceParam::instance().typeOfIonicModel == "Paci") { |
193 |
|
|
// Initialize Paci model |
194 |
|
✗ |
m_paci = new PaciSolver(m_fstransient); |
195 |
|
✗ |
m_ionic = m_paci; |
196 |
|
✗ |
} else if (FelisceParam::instance().typeOfIonicModel == "BCL") { |
197 |
|
|
// Initialize BCL model |
198 |
|
✗ |
m_bcl = new BCLSolver(m_fstransient); |
199 |
|
✗ |
m_ionic = m_bcl; |
200 |
|
|
} |
201 |
|
|
|
202 |
|
✗ |
userDefinedIonicModel(); |
203 |
|
|
|
204 |
|
|
//Fix order of bdf for ionic solver. |
205 |
|
✗ |
if ( (FelisceParam::instance().typeOfIonicModel == "MV") || (FelisceParam::instance().typeOfIonicModel == "courtAtriaMvVent") ) { |
206 |
|
✗ |
m_ionic->defineOrderBdf(FelisceParam::instance().orderBdfIonic,3); |
207 |
|
✗ |
m_W_0.resize(3); |
208 |
|
✗ |
} else if(FelisceParam::instance().typeOfIonicModel == "Paci") { |
209 |
|
✗ |
m_ionic->defineOrderBdf(FelisceParam::instance().orderBdfIonic,42); |
210 |
|
✗ |
m_W_0.resize(42); |
211 |
|
✗ |
} else if(FelisceParam::instance().typeOfIonicModel == "BCL") { |
212 |
|
✗ |
m_ionic->defineOrderBdf(FelisceParam::instance().orderBdfIonic,4); |
213 |
|
✗ |
m_W_0.resize(4); |
214 |
|
|
} else |
215 |
|
✗ |
m_ionic->defineOrderBdf(FelisceParam::instance().orderBdfIonic,1); |
216 |
|
|
|
217 |
|
|
//Define order for bdf used by EDP (linearProblem). |
218 |
|
✗ |
m_bdfEdp.defineOrder(FelisceParam::instance().orderBdfEdp); |
219 |
|
✗ |
initializeAppCurrent(); |
220 |
|
✗ |
if (FelisceParam::instance().hasAppliedExteriorCurrent) |
221 |
|
✗ |
initializeAppCurrentExt(); |
222 |
|
✗ |
if(FelisceParam::instance().stateFilter) { |
223 |
|
✗ |
m_observations = new EnsightCase(); |
224 |
|
✗ |
m_observations->read(FelisceParam::instance().observDir,FelisceParam::instance().observFileName); |
225 |
|
|
} |
226 |
|
|
|
227 |
|
|
} |
228 |
|
|
|
229 |
|
✗ |
void BidomainModel::initializeAppCurrent() { |
230 |
|
✗ |
std::unordered_map<std::string, int> mapOfType; |
231 |
|
✗ |
mapOfType["ellibi"] = 1; |
232 |
|
✗ |
mapOfType["zygote"] = 2; |
233 |
|
|
|
234 |
|
✗ |
bool warning_Iapp = true; |
235 |
|
|
|
236 |
|
✗ |
switch(mapOfType[FelisceParam::instance().typeOfAppliedCurrent]) { |
237 |
|
✗ |
case 1 : |
238 |
|
✗ |
m_iApp = new AppCurrent(); |
239 |
|
✗ |
warning_Iapp = false; |
240 |
|
✗ |
break; |
241 |
|
✗ |
case 2 : |
242 |
|
✗ |
m_iApp = new zygoteIapp(); |
243 |
|
✗ |
warning_Iapp = false; |
244 |
|
✗ |
break; |
245 |
|
|
} |
246 |
|
✗ |
if (warning_Iapp) { |
247 |
|
✗ |
std::cout << "WARNING !\nApplied current not defined : used non pathological case for ellibi geometry.\n"; |
248 |
|
✗ |
m_iApp = new AppCurrent(); |
249 |
|
|
} |
250 |
|
✗ |
m_iApp->initialize(m_fstransient); |
251 |
|
|
} |
252 |
|
|
|
253 |
|
✗ |
void BidomainModel::initializeAppCurrentExt() { |
254 |
|
✗ |
m_iAppExt = new AppCurrent(); |
255 |
|
✗ |
m_iAppExt->initialize(m_fstransient); |
256 |
|
|
} |
257 |
|
|
|
258 |
|
|
|
259 |
|
✗ |
void BidomainModel::evalIapp() { |
260 |
|
✗ |
m_iAppValue.clear(); |
261 |
|
✗ |
if ( (FelisceParam::instance().typeOfAppliedCurrent == "zygote") || (FelisceParam::instance().typeOfAppliedCurrent == "heart") || (FelisceParam::instance().typeOfAppliedCurrent == "ellipseheart") ) { |
262 |
|
|
// Warning: extra-data m_linearProblemBidomain[0]->EndocardiumDistance() size has to be numDof (not number of mesh points). |
263 |
|
✗ |
m_linearProblemBidomain[0]->evalFunctionOnDof(*m_iApp,m_fstransient->time,m_iAppValue,m_linearProblemBidomain[0]->EndocardiumDistance()); |
264 |
|
|
} |
265 |
|
|
else { |
266 |
|
✗ |
m_linearProblemBidomain[0]->evalFunctionOnDof(*m_iApp,m_fstransient->time,m_iAppValue); |
267 |
|
|
} |
268 |
|
|
} |
269 |
|
|
|
270 |
|
✗ |
void BidomainModel::evalIappExt() { |
271 |
|
✗ |
m_iAppValueExt.clear(); |
272 |
|
✗ |
m_linearProblemBidomain[0]->evalFunctionOnDof(*m_iAppExt,m_fstransient->time,m_iAppValueExt); |
273 |
|
|
} |
274 |
|
|
|
275 |
|
|
|
276 |
|
|
//Calculate RHS std::vector to be multiplied by mass matrix. |
277 |
|
✗ |
void BidomainModel::finalizeRHSDP(int iProblem) |
278 |
|
|
{ |
279 |
|
✗ |
int size = m_linearProblemBidomain[iProblem]->numDofLocalPerUnknown(potTransMemb); |
280 |
|
✗ |
double* valueForPetsc = new double[size]; |
281 |
|
✗ |
double* obsForPetsc = new double[size]; |
282 |
|
✗ |
felInt* idLocalValue = new felInt[size]; |
283 |
|
✗ |
felInt* idGlobalValue = new felInt[size]; |
284 |
|
|
|
285 |
|
✗ |
double coefAm = FelisceParam::instance().Am; |
286 |
|
✗ |
double coefCm = FelisceParam::instance().Cm; |
287 |
|
|
|
288 |
|
✗ |
double* observ=nullptr; |
289 |
|
✗ |
if (FelisceParam::instance().stateFilter ) { |
290 |
|
|
//RHS of Luenb. filter |
291 |
|
✗ |
m_Observ.set(0.); |
292 |
|
✗ |
felInt sizeVar = m_linearProblemBidomain[iProblem]->numDofPerUnknown(potTransMemb); |
293 |
|
✗ |
int indexTime = m_observations->timeIndex(m_fstransient->time); |
294 |
|
✗ |
int idVariable = m_observations->variableIndex("potTransMemb"); |
295 |
|
✗ |
observ=new double[sizeVar]; |
296 |
|
✗ |
m_observations->readVariable(idVariable, indexTime,observ,sizeVar); |
297 |
|
|
} |
298 |
|
|
|
299 |
|
|
// RHS = Am*Cm*BDF.RHS (Warning : bdf.RHS containes both problem variables). |
300 |
|
✗ |
m_linearProblemBidomain[iProblem]->vector().axpy(1.,m_bdfEdp.vector()); |
301 |
|
✗ |
m_linearProblemBidomain[iProblem]->vector().scale(coefAm*coefCm); |
302 |
|
|
|
303 |
|
|
// RHS = RHS + Am*Iapp |
304 |
|
✗ |
for (felInt i = 0; i < size; i++) { |
305 |
|
✗ |
idLocalValue[i] = i; |
306 |
|
|
} |
307 |
|
✗ |
ISLocalToGlobalMappingApply(m_linearProblemBidomain[iProblem]->mappingDofLocalToDofGlobal(potTransMemb),size,&idLocalValue[0],&idGlobalValue[0]); |
308 |
|
✗ |
for (felInt i = 0; i < size; i++) { |
309 |
|
✗ |
felInt indexGlobal=idGlobalValue[i]; |
310 |
|
|
#ifdef NDEBUG |
311 |
|
|
valueForPetsc[i] = coefAm*m_iAppValue[indexGlobal]; |
312 |
|
|
#else |
313 |
|
✗ |
valueForPetsc[i] = coefAm*m_iAppValue.at(indexGlobal); |
314 |
|
|
#endif |
315 |
|
✗ |
if (FelisceParam::instance().stateFilter ) { |
316 |
|
✗ |
AOPetscToApplication(m_linearProblemBidomain[iProblem]->ao(),1,&indexGlobal); |
317 |
|
✗ |
obsForPetsc[i] = observ[indexGlobal]; |
318 |
|
|
} |
319 |
|
|
} |
320 |
|
✗ |
m_linearProblemBidomain[iProblem]->vector().setValues(size,&idGlobalValue[0],&valueForPetsc[0],ADD_VALUES); |
321 |
|
|
|
322 |
|
✗ |
m_linearProblemBidomain[iProblem]->vector().assembly(); |
323 |
|
|
|
324 |
|
|
// For the second equation : RHS = 0 ( or applied exterior current ) |
325 |
|
✗ |
int sizePotExtraCell = m_linearProblemBidomain[iProblem]->numDofLocalPerUnknown(potExtraCell); |
326 |
|
✗ |
double* valueForPetscPotExtraCell = new double[sizePotExtraCell]; |
327 |
|
✗ |
felInt* idLocalValuePotExtraCell = new felInt[sizePotExtraCell]; |
328 |
|
✗ |
felInt* idGlobalValuePotExtraCell = new felInt[sizePotExtraCell]; |
329 |
|
|
|
330 |
|
✗ |
for (felInt i = 0; i < sizePotExtraCell; i++) { |
331 |
|
✗ |
idLocalValuePotExtraCell[i] = i; |
332 |
|
|
} |
333 |
|
✗ |
ISLocalToGlobalMappingApply(m_linearProblemBidomain[iProblem]->mappingDofLocalToDofGlobal(potExtraCell),sizePotExtraCell,&idLocalValuePotExtraCell[0],&idGlobalValuePotExtraCell[0]); |
334 |
|
✗ |
for (felInt i = 0; i < sizePotExtraCell; i++) { |
335 |
|
|
// For an applied exterior current |
336 |
|
✗ |
if (FelisceParam::instance().hasAppliedExteriorCurrent) { |
337 |
|
✗ |
valueForPetscPotExtraCell[i] = m_iAppValueExt[idGlobalValuePotExtraCell[i]]; |
338 |
|
|
} |
339 |
|
|
else { |
340 |
|
✗ |
valueForPetscPotExtraCell[i] = 0.0; |
341 |
|
|
} |
342 |
|
|
} |
343 |
|
✗ |
m_linearProblemBidomain[iProblem]->vector().setValues(sizePotExtraCell ,&idGlobalValuePotExtraCell[0],&valueForPetscPotExtraCell[0],INSERT_VALUES); |
344 |
|
|
|
345 |
|
✗ |
m_linearProblemBidomain[iProblem]->vector().assembly(); |
346 |
|
|
|
347 |
|
|
// RHS = RHS + Am*Iion |
348 |
|
✗ |
if ( (FelisceParam::instance().typeOfIonicModel == "courtAtriaSchafVent") || (FelisceParam::instance().typeOfIonicModel == "courtAtriaMvVent") ) { |
349 |
|
✗ |
m_linearProblemBidomain[iProblem]->vector().axpy(coefAm,m_ionic->ion()); |
350 |
|
✗ |
m_linearProblemBidomain[iProblem]->vector().axpy(-coefAm*0.001,m_court->ion()); |
351 |
|
|
} else { |
352 |
|
✗ |
m_linearProblemBidomain[iProblem]->vector().axpy(coefAm,m_ionic->ion()); |
353 |
|
|
} |
354 |
|
|
|
355 |
|
|
// Elisa, Dec. 2014 - Luenberger filter |
356 |
|
✗ |
if (FelisceParam::instance().electrodeControl) { |
357 |
|
✗ |
m_linearProblemBidomain[iProblem]->addElectrodeCondtrol(); |
358 |
|
|
} |
359 |
|
|
|
360 |
|
|
// RHS = RHS + Am*A2*A3*H*V_m |
361 |
|
✗ |
if (FelisceParam::instance().stateFilter ) { |
362 |
|
✗ |
m_Observ.setValues(size,&idGlobalValue[0],&obsForPetsc[0],INSERT_VALUES); |
363 |
|
✗ |
m_Observ.assembly(); |
364 |
|
✗ |
pointwiseMult(m_Observ,m_Observ,m_ionic->stabTerm()); |
365 |
|
✗ |
m_linearProblemBidomain[iProblem]->vector().axpy(coefAm,m_Observ); |
366 |
|
|
} |
367 |
|
|
|
368 |
|
|
// m_linearProblemBidomain[iProblem]->applyUserDefinedStimulation(m_fstransient->time); |
369 |
|
|
|
370 |
|
✗ |
delete [] valueForPetsc; |
371 |
|
✗ |
delete [] obsForPetsc; |
372 |
|
✗ |
delete [] idLocalValue; |
373 |
|
✗ |
delete [] idGlobalValue; |
374 |
|
✗ |
delete [] valueForPetscPotExtraCell; |
375 |
|
✗ |
delete [] idLocalValuePotExtraCell; |
376 |
|
✗ |
delete [] idGlobalValuePotExtraCell; |
377 |
|
|
|
378 |
|
|
} |
379 |
|
|
|
380 |
|
✗ |
void BidomainModel::preAssemblingMatrixRHS(std::size_t iProblem) { |
381 |
|
|
//Initialize vectors for ionic models and m_bdfEdp (with m_U_0) with a std::vector with same structure as m_RHS in linearProblem. |
382 |
|
|
//(in this case give PotTransMemb with value Vmin and 0 for PotExtraCell at time t == 0). |
383 |
|
|
|
384 |
|
|
// 1) initialize Vectors for ionic solver |
385 |
|
✗ |
if ( (FelisceParam::instance().typeOfIonicModel == "MV") || (FelisceParam::instance().typeOfIonicModel == "courtAtriaMvVent") ) { |
386 |
|
|
double valueByDofW0[3]; |
387 |
|
✗ |
valueByDofW0[0] = 1.; |
388 |
|
✗ |
valueByDofW0[1] = 1.; |
389 |
|
✗ |
valueByDofW0[2] = 0.; |
390 |
|
✗ |
for (felInt i=0; i<3; i++) { |
391 |
|
✗ |
m_W_0[i].duplicateFrom(m_linearProblemBidomain[iProblem]->vector()); |
392 |
|
✗ |
m_W_0[i].set(valueByDofW0[i]); |
393 |
|
|
} |
394 |
|
|
} |
395 |
|
✗ |
else if(FelisceParam::instance().typeOfIonicModel == "Paci") { |
396 |
|
|
double valueByDofW0[42]; |
397 |
|
|
|
398 |
|
✗ |
valueByDofW0[0] = 0.75; |
399 |
|
✗ |
valueByDofW0[1] = 0.75; |
400 |
|
✗ |
valueByDofW0[2] = 0.; |
401 |
|
✗ |
valueByDofW0[3] = 0.; |
402 |
|
✗ |
valueByDofW0[4] = 1.; |
403 |
|
✗ |
valueByDofW0[5] = 1.; |
404 |
|
✗ |
valueByDofW0[6] = 1.; |
405 |
|
✗ |
valueByDofW0[7] = 0.; |
406 |
|
✗ |
valueByDofW0[8] = 1.; |
407 |
|
✗ |
valueByDofW0[9] = 0.; |
408 |
|
✗ |
valueByDofW0[10] = 1.; |
409 |
|
✗ |
valueByDofW0[11] = 0.; |
410 |
|
✗ |
valueByDofW0[12] = 0.1; |
411 |
|
✗ |
valueByDofW0[13] = 1.; |
412 |
|
✗ |
valueByDofW0[14] = 0.0002; |
413 |
|
✗ |
valueByDofW0[15] = 0.3; |
414 |
|
✗ |
valueByDofW0[16] = 10.; |
415 |
|
✗ |
valueByDofW0[17] = 0.; |
416 |
|
✗ |
valueByDofW0[18] = 0.; |
417 |
|
✗ |
valueByDofW0[19] = 0.; |
418 |
|
✗ |
valueByDofW0[20] = 0.; |
419 |
|
✗ |
valueByDofW0[21] = 0.; |
420 |
|
✗ |
valueByDofW0[22] = 0.; |
421 |
|
✗ |
valueByDofW0[23] = 0.; |
422 |
|
✗ |
valueByDofW0[24] = 0.; |
423 |
|
✗ |
valueByDofW0[25] = 0.; |
424 |
|
✗ |
valueByDofW0[26] = 0.; |
425 |
|
✗ |
valueByDofW0[27] = 0.; |
426 |
|
✗ |
valueByDofW0[28] = 0.; |
427 |
|
✗ |
valueByDofW0[29] = 0.; |
428 |
|
✗ |
valueByDofW0[30] = 0.; |
429 |
|
✗ |
valueByDofW0[31] = 0.; |
430 |
|
✗ |
valueByDofW0[32] = 0.; |
431 |
|
✗ |
valueByDofW0[33] = 0.; |
432 |
|
✗ |
valueByDofW0[34] = 0.; |
433 |
|
✗ |
valueByDofW0[35] = 0.; |
434 |
|
✗ |
valueByDofW0[36] = 0.0074621; |
435 |
|
✗ |
valueByDofW0[37] = 0.692477; |
436 |
|
✗ |
valueByDofW0[38] = 0.692574; |
437 |
|
✗ |
valueByDofW0[39] = 0.692591; |
438 |
|
✗ |
valueByDofW0[40] = 0.496116; |
439 |
|
✗ |
valueByDofW0[41] = 0.000194015; |
440 |
|
|
|
441 |
|
|
|
442 |
|
✗ |
for (felInt i=0; i<42; i++) { |
443 |
|
✗ |
m_W_0[i].duplicateFrom(m_linearProblemBidomain[iProblem]->vector()); |
444 |
|
✗ |
m_W_0[i].set(valueByDofW0[i]); |
445 |
|
|
} |
446 |
|
|
} |
447 |
|
✗ |
else if (FelisceParam::instance().typeOfIonicModel == "BCL") { |
448 |
|
|
double valueByDofW0[4]; |
449 |
|
✗ |
valueByDofW0[0] = 1.; |
450 |
|
✗ |
valueByDofW0[1] = 1.; |
451 |
|
✗ |
valueByDofW0[2] = 0.; |
452 |
|
✗ |
valueByDofW0[3] = 1.; |
453 |
|
|
|
454 |
|
✗ |
for (felInt i=0; i<4; i++) { |
455 |
|
✗ |
m_W_0[i].duplicateFrom(m_linearProblemBidomain[iProblem]->vector()); |
456 |
|
✗ |
m_W_0[i].set(valueByDofW0[i]); |
457 |
|
|
} |
458 |
|
|
} |
459 |
|
|
// "schaf", "schafRev", "courtAtriaSchafVent", "fhn" |
460 |
|
|
else { |
461 |
|
✗ |
m_W_0.resize(1); |
462 |
|
✗ |
m_W_0[0].duplicateFrom(m_linearProblemBidomain[iProblem]->vector()); |
463 |
|
✗ |
if (FelisceParam::instance().typeOfIonicModel == "fhn") { |
464 |
|
✗ |
m_W_0[0].set(0.); |
465 |
|
|
} else { |
466 |
|
✗ |
double valueByDofW_0 = 1./((FelisceParam::instance().vMax - FelisceParam::instance().vMin)*(FelisceParam::instance().vMax - FelisceParam::instance().vMin)); |
467 |
|
✗ |
m_W_0[0].set(valueByDofW_0); |
468 |
|
|
} |
469 |
|
|
} |
470 |
|
|
|
471 |
|
✗ |
if ( (FelisceParam::instance().typeOfIonicModel == "courtAtriaSchafVent") || (FelisceParam::instance().typeOfIonicModel == "courtAtriaMvVent") ) { |
472 |
|
✗ |
double valueByDofm = 0.00291; |
473 |
|
✗ |
double valueByDofh = 0.965; |
474 |
|
✗ |
double valueByDofj = 0.978; |
475 |
|
✗ |
double valueByDofao = 0.0304; |
476 |
|
✗ |
double valueByDofio = 0.999; |
477 |
|
✗ |
double valueByDofua = 0.00496; |
478 |
|
✗ |
double valueByDofui = 0.999; |
479 |
|
✗ |
double valueByDofxr = 0.0000329; |
480 |
|
✗ |
double valueByDofxs = 0.0187; |
481 |
|
✗ |
double valueByDofd = 0.000137; |
482 |
|
✗ |
double valueByDoff = 0.999837; |
483 |
|
✗ |
double valueByDoffca = 0.775; |
484 |
|
✗ |
double valueByDofvrel = 1.00; |
485 |
|
✗ |
double valueByDofwrel = 0.999; |
486 |
|
✗ |
double valueByDofnai = 11.2; |
487 |
|
✗ |
double valueByDofnao = 140.0; |
488 |
|
✗ |
double valueByDofcao = 1.8; |
489 |
|
✗ |
double valueByDofki = 139.0; |
490 |
|
✗ |
double valueByDofko = 4.5; |
491 |
|
✗ |
double valueByDofcai = 0.000102; |
492 |
|
✗ |
double valueByDofcmdn = 0.00205; |
493 |
|
✗ |
double valueByDoftrpn = 0.0118; |
494 |
|
✗ |
double valueByDofnsr = 1.49; |
495 |
|
✗ |
double valueByDofjsr = 1.49; |
496 |
|
✗ |
double valueByDofcsqn = 6.51; |
497 |
|
|
|
498 |
|
✗ |
m_m.duplicateFrom(m_linearProblemBidomain[iProblem]->vector()); |
499 |
|
✗ |
m_m.set(valueByDofm); |
500 |
|
✗ |
m_h.duplicateFrom(m_linearProblemBidomain[iProblem]->vector()); |
501 |
|
✗ |
m_h.set(valueByDofh); |
502 |
|
✗ |
m_j.duplicateFrom(m_linearProblemBidomain[iProblem]->vector()); |
503 |
|
✗ |
m_j.set(valueByDofj); |
504 |
|
✗ |
m_ao.duplicateFrom(m_linearProblemBidomain[iProblem]->vector()); |
505 |
|
✗ |
m_ao.set(valueByDofao); |
506 |
|
✗ |
m_io.duplicateFrom(m_linearProblemBidomain[iProblem]->vector()); |
507 |
|
✗ |
m_io.set(valueByDofio); |
508 |
|
✗ |
m_ua.duplicateFrom(m_linearProblemBidomain[iProblem]->vector()); |
509 |
|
✗ |
m_ua.set(valueByDofua); |
510 |
|
✗ |
m_ui.duplicateFrom(m_linearProblemBidomain[iProblem]->vector()); |
511 |
|
✗ |
m_ui.set(valueByDofui); |
512 |
|
✗ |
m_xr.duplicateFrom(m_linearProblemBidomain[iProblem]->vector()); |
513 |
|
✗ |
m_xr.set(valueByDofxr); |
514 |
|
✗ |
m_xs.duplicateFrom(m_linearProblemBidomain[iProblem]->vector()); |
515 |
|
✗ |
m_xs.set(valueByDofxs); |
516 |
|
✗ |
m_d.duplicateFrom(m_linearProblemBidomain[iProblem]->vector()); |
517 |
|
✗ |
m_d.set(valueByDofd); |
518 |
|
✗ |
m_f.duplicateFrom(m_linearProblemBidomain[iProblem]->vector()); |
519 |
|
✗ |
m_f.set(valueByDoff); |
520 |
|
✗ |
m_fca.duplicateFrom(m_linearProblemBidomain[iProblem]->vector()); |
521 |
|
✗ |
m_fca.set(valueByDoffca); |
522 |
|
✗ |
m_urel.duplicateFrom(m_linearProblemBidomain[iProblem]->vector()); |
523 |
|
✗ |
m_urel.set(0.); |
524 |
|
✗ |
m_vrel.duplicateFrom(m_linearProblemBidomain[iProblem]->vector()); |
525 |
|
✗ |
m_vrel.set(valueByDofvrel); |
526 |
|
✗ |
m_wrel.duplicateFrom(m_linearProblemBidomain[iProblem]->vector()); |
527 |
|
✗ |
m_wrel.set(valueByDofwrel); |
528 |
|
✗ |
m_nai.duplicateFrom(m_linearProblemBidomain[iProblem]->vector()); |
529 |
|
✗ |
m_nai.set(valueByDofnai); |
530 |
|
✗ |
m_nao.duplicateFrom(m_linearProblemBidomain[iProblem]->vector()); |
531 |
|
✗ |
m_nao.set(valueByDofnao); |
532 |
|
✗ |
m_cao.duplicateFrom(m_linearProblemBidomain[iProblem]->vector()); |
533 |
|
✗ |
m_cao.set(valueByDofcao); |
534 |
|
✗ |
m_ki.duplicateFrom(m_linearProblemBidomain[iProblem]->vector()); |
535 |
|
✗ |
m_ki.set(valueByDofki); |
536 |
|
✗ |
m_ko.duplicateFrom(m_linearProblemBidomain[iProblem]->vector()); |
537 |
|
✗ |
m_ko.set(valueByDofko); |
538 |
|
✗ |
m_cai.duplicateFrom(m_linearProblemBidomain[iProblem]->vector()); |
539 |
|
✗ |
m_cai.set(valueByDofcai); |
540 |
|
✗ |
m_naiont.duplicateFrom(m_linearProblemBidomain[iProblem]->vector()); |
541 |
|
✗ |
m_naiont.set(0.); |
542 |
|
✗ |
m_kiont.duplicateFrom(m_linearProblemBidomain[iProblem]->vector()); |
543 |
|
✗ |
m_kiont.set(0.); |
544 |
|
✗ |
m_caiont.duplicateFrom(m_linearProblemBidomain[iProblem]->vector()); |
545 |
|
✗ |
m_caiont.set(0.); |
546 |
|
✗ |
m_ileak.duplicateFrom(m_linearProblemBidomain[iProblem]->vector()); |
547 |
|
✗ |
m_ileak.set(0.); |
548 |
|
✗ |
m_iup.duplicateFrom(m_linearProblemBidomain[iProblem]->vector()); |
549 |
|
✗ |
m_iup.set(0.); |
550 |
|
✗ |
m_itr.duplicateFrom(m_linearProblemBidomain[iProblem]->vector()); |
551 |
|
✗ |
m_itr.set(0.); |
552 |
|
✗ |
m_irel.duplicateFrom(m_linearProblemBidomain[iProblem]->vector()); |
553 |
|
✗ |
m_irel.set(0.); |
554 |
|
✗ |
m_cmdn.duplicateFrom(m_linearProblemBidomain[iProblem]->vector()); |
555 |
|
✗ |
m_cmdn.set(valueByDofcmdn); |
556 |
|
✗ |
m_trpn.duplicateFrom(m_linearProblemBidomain[iProblem]->vector()); |
557 |
|
✗ |
m_trpn.set(valueByDoftrpn); |
558 |
|
✗ |
m_nsr.duplicateFrom(m_linearProblemBidomain[iProblem]->vector()); |
559 |
|
✗ |
m_nsr.set(valueByDofnsr); |
560 |
|
✗ |
m_jsr.duplicateFrom(m_linearProblemBidomain[iProblem]->vector()); |
561 |
|
✗ |
m_jsr.set(valueByDofjsr); |
562 |
|
✗ |
m_csqn.duplicateFrom(m_linearProblemBidomain[iProblem]->vector()); |
563 |
|
✗ |
m_csqn.set(valueByDofcsqn); |
564 |
|
|
} |
565 |
|
|
|
566 |
|
|
|
567 |
|
|
|
568 |
|
|
// 2) copy structure in m_extrapolate (contains extrapolate values of linearProblem solution). |
569 |
|
✗ |
m_extrapolate.duplicateFrom(m_linearProblemBidomain[iProblem]->vector()); |
570 |
|
✗ |
m_extrapolate.set(0.); |
571 |
|
✗ |
m_extrapolateIsCreated = true; |
572 |
|
|
|
573 |
|
|
// Lun. filter |
574 |
|
✗ |
if(FelisceParam::instance().stateFilter) { |
575 |
|
✗ |
m_Observ.duplicateFrom(m_linearProblemBidomain[iProblem]->vector()); |
576 |
|
✗ |
m_Observ.set(0.); |
577 |
|
|
} |
578 |
|
|
|
579 |
|
|
// Read fibers and distance mapp (for Zygote geometry). |
580 |
|
✗ |
if(FelisceParam::instance().testCase == 1) { // Any testCase with fibers. |
581 |
|
✗ |
m_linearProblemBidomain[iProblem]->readData(*io()); |
582 |
|
|
} |
583 |
|
|
|
584 |
|
|
//Initialize heterogeneous parameters for schaf solver. |
585 |
|
✗ |
if ( (FelisceParam::instance().typeOfIonicModel == "schaf") || (FelisceParam::instance().typeOfIonicModel == "schafRev") || (FelisceParam::instance().typeOfIonicModel == "courtAtriaSchafVent") ) { |
586 |
|
✗ |
initializeSchafParameter(); |
587 |
|
✗ |
if (FelisceParam::instance().typeOfIonicModel == "courtAtriaSchafVent") |
588 |
|
✗ |
initializeCourtParameter(); |
589 |
|
|
} |
590 |
|
|
|
591 |
|
|
//Initialize heterogeneous parameters for MV solver. |
592 |
|
✗ |
if ( (FelisceParam::instance().typeOfIonicModel == "MV") || (FelisceParam::instance().typeOfIonicModel == "courtAtriaMvVent") ) { |
593 |
|
✗ |
initializeMVParameter(); |
594 |
|
✗ |
if (FelisceParam::instance().typeOfIonicModel == "courtAtriaMvVent") |
595 |
|
✗ |
initializeCourtParameter(); |
596 |
|
|
} |
597 |
|
|
|
598 |
|
|
//Initialize heterogeneous parameters for BCL solver. |
599 |
|
✗ |
if (FelisceParam::instance().typeOfIonicModel == "BCL") { |
600 |
|
✗ |
initializeBCLParameter(); |
601 |
|
|
} |
602 |
|
|
|
603 |
|
|
//Initialize heterogeneous parameters for FhN solver. |
604 |
|
✗ |
if (FelisceParam::instance().typeOfIonicModel == "fhn") { |
605 |
|
✗ |
initializeFhNParameter(); |
606 |
|
|
} |
607 |
|
|
|
608 |
|
✗ |
if(FelisceParam::instance().useROM) { |
609 |
|
✗ |
std::cout << "m_rom->initializeRom\n"; |
610 |
|
✗ |
m_rom = new RomPETSC(); |
611 |
|
✗ |
m_rom->initializeRom(m_fstransient, MpiInfo::petscComm(), m_linearProblemBidomain[iProblem]->ao()); |
612 |
|
✗ |
m_linearProblemBidomain[iProblem]->initializeROM(m_rom); |
613 |
|
✗ |
m_rom->initializeRomBasis(); |
614 |
|
|
} |
615 |
|
|
|
616 |
|
|
} |
617 |
|
|
|
618 |
|
✗ |
void BidomainModel::postAssemblingMatrixRHS(std::size_t iProblem) { |
619 |
|
|
|
620 |
|
|
// Evaluate applied current std::vector. |
621 |
|
✗ |
if (iProblem == 0) { |
622 |
|
✗ |
evalIapp(); |
623 |
|
✗ |
if (FelisceParam::instance().hasAppliedExteriorCurrent) |
624 |
|
✗ |
evalIappExt(); |
625 |
|
|
} |
626 |
|
|
|
627 |
|
|
|
628 |
|
|
// Calculate RHS std::vector (to be multiplied by mass matrix). |
629 |
|
|
// RHS = Am*Ion + Am*Iapp + ... |
630 |
|
✗ |
finalizeRHSDP(iProblem); |
631 |
|
|
|
632 |
|
|
// Calculate m_RHS : |
633 |
|
|
// if model == 'Bidomain' => m_RHS = mass * RHS. |
634 |
|
|
// if model == 'BidomainDecoupled' : |
635 |
|
|
// iProblem==0 => m_RHS=mass*RHS+Ksigmai*u_e |
636 |
|
|
// iProblem==1 => m_RHS=Ksigmai*V_m |
637 |
|
✗ |
m_linearProblemBidomain[iProblem]->addMatrixRHS(); |
638 |
|
|
|
639 |
|
|
//RHS of Luenb. filter (Cesare) |
640 |
|
✗ |
if (FelisceParam::instance().stateFilter) { |
641 |
|
✗ |
if (iProblem == 0) { |
642 |
|
✗ |
double coefAm = FelisceParam::instance().Am; |
643 |
|
✗ |
m_luenbFlter.zeroEntries(); |
644 |
|
✗ |
m_luenbFlter.copyFrom(m_linearProblemBidomain[iProblem]->matrix(1),SAME_NONZERO_PATTERN); |
645 |
|
✗ |
m_luenbFlter.assembly(MAT_FINAL_ASSEMBLY); |
646 |
|
|
// copy diagonal stabiliz. matrix in luenberger matrix |
647 |
|
✗ |
m_luenbFlter.diagonalScale(PetscVector::null(),m_ionic->stabTerm()); |
648 |
|
✗ |
m_luenbFlter.scale(coefAm); |
649 |
|
✗ |
m_linearProblemBidomain[iProblem]->matrix(0).axpy(1,m_luenbFlter,SAME_NONZERO_PATTERN); |
650 |
|
|
} |
651 |
|
|
} |
652 |
|
|
} |
653 |
|
|
|
654 |
|
|
// Re-define class to change clearMatrixRHS call |
655 |
|
✗ |
void BidomainModel::updateTime(const FlagMatrixRHS flagMatrixRHS) { |
656 |
|
✗ |
for(std::size_t i=0; i<m_linearProblem.size(); i++) { |
657 |
|
✗ |
m_linearProblemBidomain[i]->clearMatrixRHS(flagMatrixRHS); |
658 |
|
|
} |
659 |
|
✗ |
m_fstransient->iteration++; |
660 |
|
✗ |
m_fstransient->time +=m_fstransient->timeStep; |
661 |
|
|
} |
662 |
|
|
|
663 |
|
✗ |
void BidomainModel::initializeSchafParameter() { |
664 |
|
|
// Initialize TauClose. |
665 |
|
✗ |
if (FelisceParam::instance().hasHeteroTauClose) { |
666 |
|
✗ |
m_schaf->fctTauClose().initialize(m_fstransient); |
667 |
|
✗ |
if ( (FelisceParam::instance().typeOfAppliedCurrent == "zygote") || (FelisceParam::instance().typeOfAppliedCurrent == "heart") || (FelisceParam::instance().typeOfAppliedCurrent == "ellipseheart") ) { |
668 |
|
|
// Warning: extra-data m_linearProblemBidomain[0]->EndocardiumDistance() size has to be numDof (not number of mesh points). |
669 |
|
✗ |
m_linearProblemBidomain[0]->evalFunctionOnDof(m_schaf->fctTauClose(), m_schaf->tauClose(),m_linearProblemBidomain[0]->EndocardiumDistance()); |
670 |
|
|
} else { |
671 |
|
✗ |
m_linearProblemBidomain[0]->evalFunctionOnDof(m_schaf->fctTauClose(), m_schaf->tauClose()); |
672 |
|
|
} |
673 |
|
|
} |
674 |
|
|
// Initialize TauOut : heterogeneous in case of infarct. |
675 |
|
✗ |
if (FelisceParam::instance().hasInfarct) { |
676 |
|
✗ |
m_schaf->fctTauOut().initialize(m_fstransient); |
677 |
|
✗ |
m_linearProblemBidomain[0]->evalFunctionOnDof(m_schaf->fctTauOut(), m_schaf->tauOut()); |
678 |
|
|
} |
679 |
|
|
} |
680 |
|
|
|
681 |
|
✗ |
void BidomainModel::initializeMVParameter() { |
682 |
|
✗ |
if (FelisceParam::instance().CellsType == "hetero") { |
683 |
|
✗ |
m_mv->MVCoeff().initialize(m_fstransient); |
684 |
|
|
// Warning: extra-data m_linearProblemBidomain[0]->EndocardiumDistance() size has to be numDof (not number of mesh points). |
685 |
|
✗ |
m_linearProblemBidomain[0]->evalFunctionOnDof(m_mv->MVCoeff(), m_mv->cellType(),m_linearProblemBidomain[0]->EndocardiumDistance()); |
686 |
|
|
} |
687 |
|
|
// Initialize gfi |
688 |
|
✗ |
if (FelisceParam::instance().hasInfarct) { |
689 |
|
✗ |
m_mv->fctTauOut().initialize(m_fstransient); |
690 |
|
✗ |
m_linearProblemBidomain[0]->evalFunctionOnDof(m_mv->fctTauOut(), m_mv->tauOut()); |
691 |
|
|
} |
692 |
|
|
} |
693 |
|
|
|
694 |
|
✗ |
void BidomainModel::initializeBCLParameter() { |
695 |
|
✗ |
if (FelisceParam::instance().CellsType == "hetero") { |
696 |
|
✗ |
m_bcl->MVCoeff().initialize(m_fstransient); |
697 |
|
|
// Warning: extra-data m_linearProblemBidomain[0]->EndocardiumDistance() size has to be numDof (not number of mesh points). |
698 |
|
✗ |
m_linearProblemBidomain[0]->evalFunctionOnDof(m_bcl->MVCoeff(), m_bcl->cellType(),m_linearProblemBidomain[0]->EndocardiumDistance()); |
699 |
|
|
} |
700 |
|
|
} |
701 |
|
|
|
702 |
|
✗ |
void BidomainModel::initializeFhNParameter() { |
703 |
|
✗ |
if (FelisceParam::instance().hasInfarct) { |
704 |
|
✗ |
m_fhn->fctSpar().initialize(m_fstransient); |
705 |
|
✗ |
m_linearProblemBidomain[0]->evalFunctionOnDof(m_fhn->fctSpar(), m_fhn->f0Par()); |
706 |
|
|
} |
707 |
|
|
} |
708 |
|
|
|
709 |
|
✗ |
void BidomainModel::initializeCourtParameter() { |
710 |
|
|
// Initialize gToMax |
711 |
|
✗ |
m_court->fctCourtCondIto().initialize(m_fstransient); |
712 |
|
|
// Initialize gCaLMax |
713 |
|
✗ |
m_court->fctCourtCondICaL().initialize(m_fstransient); |
714 |
|
✗ |
if (FelisceParam::instance().hasHeteroCourtPar) { |
715 |
|
✗ |
m_linearProblemBidomain[0]->evalFunctionOnDof(m_court->fctCourtCondIto(), m_court->courtCondIto()); |
716 |
|
✗ |
m_linearProblemBidomain[0]->evalFunctionOnDof(m_court->fctCourtCondICaL(), m_court->courtCondICaL()); |
717 |
|
|
} |
718 |
|
|
// Initialize multiplicative coefficient |
719 |
|
✗ |
m_court->fctCourtCondMultCoeff().initialize(m_fstransient); |
720 |
|
|
// Warning: extra-data m_linearProblemBidomain[0]->Reference() size has to be numDof (not number of mesh points). |
721 |
|
✗ |
if (FelisceParam::instance().testCase == 1) |
722 |
|
✗ |
m_linearProblemBidomain[0]->evalFunctionOnDof(m_court->fctCourtCondMultCoeff(), m_court->courtCondMultCoeff(), m_linearProblemBidomain[0]->Reference()); |
723 |
|
|
else |
724 |
|
✗ |
m_linearProblemBidomain[0]->evalFunctionOnDof(m_court->fctCourtCondMultCoeff(), m_court->courtCondMultCoeff()); |
725 |
|
|
} |
726 |
|
|
|
727 |
|
|
// Pay attention on the call of this : |
728 |
|
|
// call BidomainModel::writeSolution() function instead of (non-virtual) Model::writeSolution() function |
729 |
|
✗ |
void BidomainModel::writeSolution() { |
730 |
|
✗ |
if (MpiInfo::rankProc() == 0) { |
731 |
|
✗ |
if (m_meshIsWritten == false) writeMesh(); |
732 |
|
|
} |
733 |
|
|
|
734 |
|
✗ |
if( (m_fstransient->iteration % FelisceParam::instance().frequencyWriteSolution == 0) or m_hasFinished) { |
735 |
|
✗ |
if(FelisceParam::verbose() > 1) PetscPrintf(MpiInfo::petscComm(),"Write solutions\n"); |
736 |
|
✗ |
for (std::size_t ipb = 0; ipb < m_linearProblem.size(); ipb++) { |
737 |
|
|
|
738 |
|
✗ |
m_linearProblemBidomain[ipb]->writeSolution(MpiInfo::rankProc(), m_ios, m_fstransient->time, m_fstransient->iteration); |
739 |
|
|
|
740 |
|
✗ |
if(FelisceParam::instance().printIonicVar) { |
741 |
|
|
int rankProc; |
742 |
|
✗ |
MPI_Comm_rank(MpiInfo::petscComm(),&rankProc); |
743 |
|
✗ |
if ( m_fstransient->iteration == 0 ) { |
744 |
|
✗ |
double * ionicSolToSave = new double[m_linearProblemBidomain[ipb]->numDofPerUnknown(0)]; |
745 |
|
✗ |
m_linearProblemBidomain[ipb]->fromVecToDoubleStar(ionicSolToSave, m_W_0[0], rankProc, 1, m_linearProblemBidomain[ipb]->numDofPerUnknown(0)); |
746 |
|
✗ |
if (rankProc == 0) { |
747 |
|
✗ |
m_linearProblemBidomain[ipb]->writeEnsightScalar(ionicSolToSave, static_cast<int>(m_fstransient->iteration/FelisceParam::instance().frequencyWriteSolution), "ionicVar"); |
748 |
|
✗ |
m_linearProblemBidomain[ipb]->writeEnsightCase(static_cast<int>(m_fstransient->iteration/FelisceParam::instance().frequencyWriteSolution)+1,"ionicVar"); |
749 |
|
|
} |
750 |
|
✗ |
delete [] ionicSolToSave; |
751 |
|
|
|
752 |
|
✗ |
if(FelisceParam::instance().typeOfIonicModel == "MV") { |
753 |
|
|
//MV |
754 |
|
✗ |
double * ionicSolToSaveV = new double[m_linearProblemBidomain[ipb]->numDofPerUnknown(0)]; |
755 |
|
✗ |
double * ionicSolToSaveW = new double[m_linearProblemBidomain[ipb]->numDofPerUnknown(0)]; |
756 |
|
✗ |
double * ionicSolToSaveS = new double[m_linearProblemBidomain[ipb]->numDofPerUnknown(0)]; |
757 |
|
|
|
758 |
|
✗ |
std::vector<std::string> varNames; |
759 |
|
✗ |
varNames.emplace_back("v"); |
760 |
|
✗ |
varNames.emplace_back("w"); |
761 |
|
✗ |
varNames.emplace_back("s"); |
762 |
|
|
|
763 |
|
✗ |
m_linearProblemBidomain[ipb]->fromVecToDoubleStar(ionicSolToSaveV, m_W_0[0], rankProc, 1, m_linearProblemBidomain[ipb]->numDofPerUnknown(0)); |
764 |
|
✗ |
m_linearProblemBidomain[ipb]->fromVecToDoubleStar(ionicSolToSaveW, m_W_0[1], rankProc, 1, m_linearProblemBidomain[ipb]->numDofPerUnknown(0)); |
765 |
|
✗ |
m_linearProblemBidomain[ipb]->fromVecToDoubleStar(ionicSolToSaveS, m_W_0[2], rankProc, 1, m_linearProblemBidomain[ipb]->numDofPerUnknown(0)); |
766 |
|
|
|
767 |
|
✗ |
if (MpiInfo::rankProc() == 0) { |
768 |
|
✗ |
m_linearProblemBidomain[ipb]->writeEnsightScalar(ionicSolToSaveV, static_cast<int>(m_fstransient->iteration/FelisceParam::instance().frequencyWriteSolution), varNames[0]); |
769 |
|
✗ |
m_linearProblemBidomain[ipb]->writeEnsightScalar(ionicSolToSaveW, static_cast<int>(m_fstransient->iteration/FelisceParam::instance().frequencyWriteSolution), varNames[1]); |
770 |
|
✗ |
m_linearProblemBidomain[ipb]->writeEnsightScalar(ionicSolToSaveS, static_cast<int>(m_fstransient->iteration/FelisceParam::instance().frequencyWriteSolution), varNames[2]); |
771 |
|
|
|
772 |
|
✗ |
m_linearProblemBidomain[ipb]->writeEnsightCaseCurrents(static_cast<int>(m_fstransient->iteration/FelisceParam::instance().frequencyWriteSolution)+1,varNames); |
773 |
|
|
|
774 |
|
|
} |
775 |
|
|
|
776 |
|
✗ |
delete [] ionicSolToSaveV; |
777 |
|
✗ |
delete [] ionicSolToSaveW; |
778 |
|
✗ |
delete [] ionicSolToSaveS; |
779 |
|
|
} |
780 |
|
|
|
781 |
|
✗ |
else if(FelisceParam::instance().typeOfIonicModel == "BCL") { |
782 |
|
|
//BCL |
783 |
|
✗ |
double * ionicSolToSaveH = new double[m_linearProblemBidomain[ipb]->numDofPerUnknown(0)]; |
784 |
|
✗ |
double * ionicSolToSaveF = new double[m_linearProblemBidomain[ipb]->numDofPerUnknown(0)]; |
785 |
|
✗ |
double * ionicSolToSaveR = new double[m_linearProblemBidomain[ipb]->numDofPerUnknown(0)]; |
786 |
|
✗ |
double * ionicSolToSaveS = new double[m_linearProblemBidomain[ipb]->numDofPerUnknown(0)]; |
787 |
|
|
|
788 |
|
✗ |
std::vector<std::string> varNames; |
789 |
|
✗ |
varNames.emplace_back("H"); |
790 |
|
✗ |
varNames.emplace_back("F"); |
791 |
|
✗ |
varNames.emplace_back("R"); |
792 |
|
✗ |
varNames.emplace_back("S"); |
793 |
|
|
|
794 |
|
✗ |
m_linearProblemBidomain[ipb]->fromVecToDoubleStar(ionicSolToSaveH, m_W_0[0], rankProc, 1, m_linearProblemBidomain[ipb]->numDofPerUnknown(0)); |
795 |
|
✗ |
m_linearProblemBidomain[ipb]->fromVecToDoubleStar(ionicSolToSaveF, m_W_0[1], rankProc, 1, m_linearProblemBidomain[ipb]->numDofPerUnknown(0)); |
796 |
|
✗ |
m_linearProblemBidomain[ipb]->fromVecToDoubleStar(ionicSolToSaveR, m_W_0[2], rankProc, 1, m_linearProblemBidomain[ipb]->numDofPerUnknown(0)); |
797 |
|
✗ |
m_linearProblemBidomain[ipb]->fromVecToDoubleStar(ionicSolToSaveS, m_W_0[3], rankProc, 1, m_linearProblemBidomain[ipb]->numDofPerUnknown(0)); |
798 |
|
|
|
799 |
|
✗ |
if (MpiInfo::rankProc() == 0) { |
800 |
|
✗ |
m_linearProblemBidomain[ipb]->writeEnsightScalar(ionicSolToSaveH, static_cast<int>(m_fstransient->iteration/FelisceParam::instance().frequencyWriteSolution), varNames[0]); |
801 |
|
✗ |
m_linearProblemBidomain[ipb]->writeEnsightScalar(ionicSolToSaveF, static_cast<int>(m_fstransient->iteration/FelisceParam::instance().frequencyWriteSolution), varNames[1]); |
802 |
|
✗ |
m_linearProblemBidomain[ipb]->writeEnsightScalar(ionicSolToSaveR, static_cast<int>(m_fstransient->iteration/FelisceParam::instance().frequencyWriteSolution), varNames[2]); |
803 |
|
✗ |
m_linearProblemBidomain[ipb]->writeEnsightScalar(ionicSolToSaveS, static_cast<int>(m_fstransient->iteration/FelisceParam::instance().frequencyWriteSolution), varNames[3]); |
804 |
|
|
|
805 |
|
✗ |
m_linearProblemBidomain[ipb]->writeEnsightCaseCurrents(static_cast<int>(m_fstransient->iteration/FelisceParam::instance().frequencyWriteSolution)+1,varNames); |
806 |
|
|
|
807 |
|
|
} |
808 |
|
|
|
809 |
|
✗ |
delete [] ionicSolToSaveH; |
810 |
|
✗ |
delete [] ionicSolToSaveF; |
811 |
|
✗ |
delete [] ionicSolToSaveR; |
812 |
|
✗ |
delete [] ionicSolToSaveS; |
813 |
|
|
|
814 |
|
|
} |
815 |
|
|
|
816 |
|
|
} else { |
817 |
|
✗ |
if(FelisceParam::instance().typeOfIonicModel == "MV") { |
818 |
|
✗ |
std::vector<std::string> varNames; |
819 |
|
✗ |
varNames.emplace_back("v"); |
820 |
|
✗ |
varNames.emplace_back("w"); |
821 |
|
✗ |
varNames.emplace_back("s"); |
822 |
|
|
|
823 |
|
✗ |
double * ionicSolToSave = new double[m_linearProblemBidomain[ipb]->numDofPerUnknown(0)]; |
824 |
|
|
|
825 |
|
✗ |
for (int iMvVar=0; iMvVar<3; iMvVar++) { |
826 |
|
✗ |
m_linearProblemBidomain[ipb]->fromVecToDoubleStar(ionicSolToSave, m_ionic->vec_sol()[iMvVar], rankProc, 1, m_linearProblemBidomain[ipb]->numDofPerUnknown(0)); |
827 |
|
✗ |
if (MpiInfo::rankProc() == 0) |
828 |
|
✗ |
m_linearProblemBidomain[ipb]->writeEnsightScalar(ionicSolToSave, static_cast<int>(m_fstransient->iteration/FelisceParam::instance().frequencyWriteSolution), varNames[iMvVar]); |
829 |
|
|
} |
830 |
|
|
|
831 |
|
✗ |
if (MpiInfo::rankProc() == 0) |
832 |
|
✗ |
m_linearProblemBidomain[ipb]->writeEnsightCaseCurrents(static_cast<int>(m_fstransient->iteration/FelisceParam::instance().frequencyWriteSolution)+1,varNames); |
833 |
|
|
|
834 |
|
✗ |
delete [] ionicSolToSave; |
835 |
|
✗ |
} else if(FelisceParam::instance().typeOfIonicModel == "BCL") { |
836 |
|
✗ |
std::vector<std::string> varNames; |
837 |
|
✗ |
varNames.emplace_back("H"); |
838 |
|
✗ |
varNames.emplace_back("F"); |
839 |
|
✗ |
varNames.emplace_back("R"); |
840 |
|
✗ |
varNames.emplace_back("S"); |
841 |
|
|
|
842 |
|
✗ |
double * ionicSolToSave = new double[m_linearProblemBidomain[ipb]->numDofPerUnknown(0)]; |
843 |
|
|
|
844 |
|
✗ |
for (int iMvVar=0; iMvVar<4; iMvVar++) { |
845 |
|
✗ |
m_linearProblemBidomain[ipb]->fromVecToDoubleStar(ionicSolToSave, m_ionic->vec_sol()[iMvVar], rankProc, 1, m_linearProblemBidomain[ipb]->numDofPerUnknown(0)); |
846 |
|
✗ |
if (MpiInfo::rankProc() == 0) |
847 |
|
✗ |
m_linearProblemBidomain[ipb]->writeEnsightScalar(ionicSolToSave, static_cast<int>(m_fstransient->iteration/FelisceParam::instance().frequencyWriteSolution), varNames[iMvVar]); |
848 |
|
|
} |
849 |
|
|
|
850 |
|
✗ |
if (MpiInfo::rankProc() == 0) |
851 |
|
✗ |
m_linearProblemBidomain[ipb]->writeEnsightCaseCurrents(static_cast<int>(m_fstransient->iteration/FelisceParam::instance().frequencyWriteSolution)+1,varNames); |
852 |
|
|
|
853 |
|
✗ |
delete [] ionicSolToSave; |
854 |
|
✗ |
} else if(FelisceParam::instance().typeOfIonicModel == "Paci") { |
855 |
|
✗ |
std::vector<std::string> varNames; |
856 |
|
✗ |
varNames.emplace_back("h"); |
857 |
|
✗ |
varNames.emplace_back("j"); |
858 |
|
✗ |
varNames.emplace_back("m"); |
859 |
|
✗ |
varNames.emplace_back("d"); |
860 |
|
✗ |
varNames.emplace_back("fCa"); |
861 |
|
✗ |
varNames.emplace_back("f1"); |
862 |
|
✗ |
varNames.emplace_back("f2"); |
863 |
|
✗ |
varNames.emplace_back("r"); |
864 |
|
✗ |
varNames.emplace_back("q"); |
865 |
|
✗ |
varNames.emplace_back("xr1"); |
866 |
|
✗ |
varNames.emplace_back("xr2"); |
867 |
|
✗ |
varNames.emplace_back("xs"); |
868 |
|
✗ |
varNames.emplace_back("xf"); |
869 |
|
✗ |
varNames.emplace_back("g"); |
870 |
|
✗ |
varNames.emplace_back("Cai"); |
871 |
|
✗ |
varNames.emplace_back("Casr"); |
872 |
|
✗ |
varNames.emplace_back("Nai"); |
873 |
|
✗ |
varNames.emplace_back("df1"); |
874 |
|
✗ |
varNames.emplace_back("INaCa"); |
875 |
|
✗ |
varNames.emplace_back("IpCa"); |
876 |
|
✗ |
varNames.emplace_back("IbCa"); |
877 |
|
✗ |
varNames.emplace_back("ENa"); |
878 |
|
✗ |
varNames.emplace_back("EKs"); |
879 |
|
✗ |
varNames.emplace_back("ECa"); |
880 |
|
✗ |
varNames.emplace_back("INa"); |
881 |
|
✗ |
varNames.emplace_back("ICaL"); |
882 |
|
✗ |
varNames.emplace_back("Ito"); |
883 |
|
✗ |
varNames.emplace_back("IKr"); |
884 |
|
✗ |
varNames.emplace_back("IKs"); |
885 |
|
✗ |
varNames.emplace_back("IK1"); |
886 |
|
✗ |
varNames.emplace_back("If"); |
887 |
|
✗ |
varNames.emplace_back("INaK"); |
888 |
|
✗ |
varNames.emplace_back("Irel"); |
889 |
|
✗ |
varNames.emplace_back("Iup"); |
890 |
|
✗ |
varNames.emplace_back("Ileak"); |
891 |
|
✗ |
varNames.emplace_back("IbNa"); |
892 |
|
|
|
893 |
|
✗ |
double * ionicSolToSave = new double[m_linearProblemBidomain[ipb]->numDofPerUnknown(0)]; |
894 |
|
|
|
895 |
|
✗ |
for (int iPaciVar=0; iPaciVar<42; iPaciVar++) { |
896 |
|
✗ |
m_linearProblemBidomain[ipb]->fromVecToDoubleStar(ionicSolToSave, m_ionic->vec_sol()[iPaciVar], rankProc, 1, m_linearProblemBidomain[ipb]->numDofPerUnknown(0)); |
897 |
|
✗ |
if (MpiInfo::rankProc() == 0) |
898 |
|
✗ |
m_linearProblemBidomain[ipb]->writeEnsightScalar(ionicSolToSave,static_cast<int>(m_fstransient->iteration/FelisceParam::instance().frequencyWriteSolution),varNames[iPaciVar]); |
899 |
|
|
} |
900 |
|
✗ |
if (MpiInfo::rankProc() == 0) |
901 |
|
✗ |
m_linearProblemBidomain[ipb]->writeEnsightCaseCurrents(static_cast<int>(m_fstransient->iteration/FelisceParam::instance().frequencyWriteSolution)+1,varNames); |
902 |
|
|
|
903 |
|
✗ |
delete [] ionicSolToSave; |
904 |
|
|
|
905 |
|
✗ |
} else { |
906 |
|
✗ |
double * ionicSolToSave = new double[m_linearProblemBidomain[ipb]->numDofPerUnknown(0)]; |
907 |
|
✗ |
m_linearProblemBidomain[ipb]->fromVecToDoubleStar(ionicSolToSave, m_ionic->sol(), rankProc, 1, m_linearProblemBidomain[ipb]->numDofPerUnknown(0)); |
908 |
|
✗ |
if (rankProc == 0) { |
909 |
|
✗ |
m_linearProblemBidomain[ipb]->writeEnsightScalar(ionicSolToSave, static_cast<int>(m_fstransient->iteration/FelisceParam::instance().frequencyWriteSolution), "ionicVar"); |
910 |
|
✗ |
m_linearProblemBidomain[ipb]->writeEnsightCase(static_cast<int>(m_fstransient->iteration/FelisceParam::instance().frequencyWriteSolution)+1,"ionicVar"); |
911 |
|
|
} |
912 |
|
✗ |
delete [] ionicSolToSave; |
913 |
|
|
} |
914 |
|
|
} |
915 |
|
|
} |
916 |
|
✗ |
if (MpiInfo::rankProc() == 0) { |
917 |
|
✗ |
for(std::size_t iio=0; iio<m_ios.size(); ++iio) { |
918 |
|
✗ |
if ( m_ios[iio]->typeOutput == 1 ) // 1 : ensight{ |
919 |
|
✗ |
m_ios[iio]->postProcess(m_fstransient->time/*, !m_sameGeometricMeshForVariable*/); |
920 |
|
|
} |
921 |
|
|
} |
922 |
|
|
} |
923 |
|
|
} |
924 |
|
|
|
925 |
|
✗ |
for (std::size_t ipb = 0; ipb < m_linearProblem.size(); ipb++) { |
926 |
|
✗ |
const int iMesh = m_linearProblem[ipb]->currentMesh(); |
927 |
|
|
|
928 |
|
✗ |
if( m_solutionBackup[ipb].hasToBackup(m_fstransient->iteration /*+ std::max(FelisceParam::instance().orderBdfEdp,FelisceParam::instance().orderBdfIonic)*/, FelisceParam::instance().frequencyBackup) && ( m_fstransient->iteration > ( std::max(FelisceParam::instance().orderBdfEdp,FelisceParam::instance().orderBdfIonic) - 1 ) ) ) { |
929 |
|
✗ |
if(FelisceParam::verbose() > 1) |
930 |
|
✗ |
PetscPrintf(MpiInfo::petscComm(),"Backup the solutions\n"); |
931 |
|
|
|
932 |
|
✗ |
double * supportVariable1 = nullptr; |
933 |
|
✗ |
double * supportVariable2 = nullptr; |
934 |
|
✗ |
double * supportVariable3 = nullptr; |
935 |
|
✗ |
double * supportVariable4 = nullptr; |
936 |
|
✗ |
double * supportVariable5 = nullptr; |
937 |
|
✗ |
double * supportVariable6 = nullptr; |
938 |
|
✗ |
std::vector<double*> supportVariableCourt; |
939 |
|
|
|
940 |
|
|
// Backup of linear problem variables |
941 |
|
✗ |
m_linearProblemBidomain[ipb]->writeSolutionBackup(m_initialCondition.listVariable(), MpiInfo::rankProc(), m_solutionBackup[ipb].listIO(), m_fstransient->time, m_fstransient->iteration); |
942 |
|
✗ |
felInt sizeSupportVariable = 0; |
943 |
|
|
|
944 |
|
|
// Backup of linear problem support variables (bdfEdp, bdfIonicModel, IonicModel variables, ... ) |
945 |
|
✗ |
if ( FelisceParam::instance().orderBdfEdp > 0 ) { |
946 |
|
✗ |
m_bdfEdp.sol_n().getSize(&sizeSupportVariable); |
947 |
|
✗ |
supportVariable1 = new double[sizeSupportVariable]; |
948 |
|
✗ |
m_linearProblemBidomain[ipb]->fromVecToDoubleStar(supportVariable1, m_bdfEdp.sol_n(), MpiInfo::rankProc(), 1); |
949 |
|
✗ |
m_linearProblem[ipb]->writeSupportVariableBackup("bdf_sol_n", supportVariable1, sizeSupportVariable, 1, MpiInfo::rankProc(), m_solutionBackup[ipb].io(iMesh), m_fstransient->time, m_fstransient->iteration); |
950 |
|
|
} |
951 |
|
|
|
952 |
|
✗ |
if ( FelisceParam::instance().orderBdfEdp > 1 ) { |
953 |
|
✗ |
m_bdfEdp.sol_n_1().getSize(&sizeSupportVariable); |
954 |
|
✗ |
supportVariable2 = new double[sizeSupportVariable]; |
955 |
|
✗ |
m_linearProblemBidomain[ipb]->fromVecToDoubleStar(supportVariable2, m_bdfEdp.sol_n_1(), MpiInfo::rankProc(), 1); |
956 |
|
✗ |
m_linearProblem[ipb]->writeSupportVariableBackup("bdf_sol_n_1", supportVariable2, sizeSupportVariable, 1, MpiInfo::rankProc(), m_solutionBackup[ipb].io(iMesh), m_fstransient->time, m_fstransient->iteration); |
957 |
|
|
} |
958 |
|
|
|
959 |
|
✗ |
if ( FelisceParam::instance().orderBdfEdp > 2 ) { |
960 |
|
✗ |
m_bdfEdp.sol_n_2().getSize(&sizeSupportVariable); |
961 |
|
✗ |
supportVariable3 = new double[sizeSupportVariable]; |
962 |
|
✗ |
m_linearProblemBidomain[ipb]->fromVecToDoubleStar(supportVariable3, m_bdfEdp.sol_n_2(), MpiInfo::rankProc(), 1); |
963 |
|
✗ |
m_linearProblem[ipb]->writeSupportVariableBackup("bdf_sol_n_2", supportVariable3, sizeSupportVariable, 1, MpiInfo::rankProc(), m_solutionBackup[ipb].io(iMesh), m_fstransient->time, m_fstransient->iteration); |
964 |
|
|
} |
965 |
|
|
|
966 |
|
✗ |
if ( (FelisceParam::instance().typeOfIonicModel == "MV") || (FelisceParam::instance().typeOfIonicModel == "courtAtriaMvVent") ) { |
967 |
|
|
//MV |
968 |
|
✗ |
double * ionicSolToSaveV = nullptr; |
969 |
|
✗ |
double * ionicSolToSaveW = nullptr; |
970 |
|
✗ |
double * ionicSolToSaveS = nullptr; |
971 |
|
✗ |
if ( FelisceParam::instance().orderBdfIonic > 0 ) { |
972 |
|
✗ |
m_ionic->vec_sol(0).getSize(&sizeSupportVariable); |
973 |
|
✗ |
ionicSolToSaveV = new double[sizeSupportVariable]; |
974 |
|
✗ |
m_ionic->vec_sol(1).getSize(&sizeSupportVariable); |
975 |
|
✗ |
ionicSolToSaveW = new double[sizeSupportVariable]; |
976 |
|
✗ |
m_ionic->vec_sol(2).getSize(&sizeSupportVariable); |
977 |
|
✗ |
ionicSolToSaveS = new double[sizeSupportVariable]; |
978 |
|
|
|
979 |
|
✗ |
m_linearProblemBidomain[ipb]->fromVecToDoubleStar(ionicSolToSaveV, m_ionic->vec_sol(0), MpiInfo::rankProc(), 1); |
980 |
|
✗ |
m_linearProblem[ipb]->writeSupportVariableBackup("ion_sol_v_n", ionicSolToSaveV, sizeSupportVariable, 1, MpiInfo::rankProc(), m_solutionBackup[ipb].io(iMesh), m_fstransient->time, m_fstransient->iteration); |
981 |
|
|
|
982 |
|
✗ |
m_linearProblemBidomain[ipb]->fromVecToDoubleStar(ionicSolToSaveW, m_ionic->vec_sol(1), MpiInfo::rankProc(), 1); |
983 |
|
✗ |
m_linearProblem[ipb]->writeSupportVariableBackup("ion_sol_w_n", ionicSolToSaveW, sizeSupportVariable, 1, MpiInfo::rankProc(), m_solutionBackup[ipb].io(iMesh), m_fstransient->time, m_fstransient->iteration); |
984 |
|
|
|
985 |
|
✗ |
m_linearProblemBidomain[ipb]->fromVecToDoubleStar(ionicSolToSaveS, m_ionic->vec_sol(2), MpiInfo::rankProc(), 1); |
986 |
|
✗ |
m_linearProblem[ipb]->writeSupportVariableBackup("ion_sol_s_n", ionicSolToSaveS, sizeSupportVariable, 1, MpiInfo::rankProc(), m_solutionBackup[ipb].io(iMesh), m_fstransient->time, m_fstransient->iteration); |
987 |
|
|
} |
988 |
|
|
|
989 |
|
✗ |
if ( FelisceParam::instance().orderBdfIonic > 1 ) { |
990 |
|
✗ |
m_linearProblemBidomain[ipb]->fromVecToDoubleStar(ionicSolToSaveV, m_ionic->vec_sol_n_1()[0], MpiInfo::rankProc(), 1); |
991 |
|
✗ |
m_linearProblem[ipb]->writeSupportVariableBackup("ion_sol_v_n_1", ionicSolToSaveV, sizeSupportVariable, 1, MpiInfo::rankProc(), m_solutionBackup[ipb].io(iMesh), m_fstransient->time, m_fstransient->iteration); |
992 |
|
|
|
993 |
|
✗ |
m_linearProblemBidomain[ipb]->fromVecToDoubleStar(ionicSolToSaveW, m_ionic->vec_sol_n_1()[1], MpiInfo::rankProc(), 1); |
994 |
|
✗ |
m_linearProblem[ipb]->writeSupportVariableBackup("ion_sol_w_n_1", ionicSolToSaveW, sizeSupportVariable, 1, MpiInfo::rankProc(), m_solutionBackup[ipb].io(iMesh), m_fstransient->time, m_fstransient->iteration); |
995 |
|
|
|
996 |
|
✗ |
m_linearProblemBidomain[ipb]->fromVecToDoubleStar(ionicSolToSaveS, m_ionic->vec_sol_n_1()[2], MpiInfo::rankProc(), 1); |
997 |
|
✗ |
m_linearProblem[ipb]->writeSupportVariableBackup("ion_sol_s_n_1", ionicSolToSaveS, sizeSupportVariable, 1, MpiInfo::rankProc(), m_solutionBackup[ipb].io(iMesh), m_fstransient->time, m_fstransient->iteration); |
998 |
|
|
} |
999 |
|
|
|
1000 |
|
✗ |
if ( FelisceParam::instance().orderBdfIonic > 2 ) { |
1001 |
|
✗ |
m_linearProblemBidomain[ipb]->fromVecToDoubleStar(ionicSolToSaveV, m_ionic->vec_sol_n_2()[0], MpiInfo::rankProc(), 1); |
1002 |
|
✗ |
m_linearProblem[ipb]->writeSupportVariableBackup("ion_sol_v_n_2", ionicSolToSaveV, sizeSupportVariable, 1, MpiInfo::rankProc(), m_solutionBackup[ipb].io(iMesh), m_fstransient->time, m_fstransient->iteration); |
1003 |
|
|
|
1004 |
|
✗ |
m_linearProblemBidomain[ipb]->fromVecToDoubleStar(ionicSolToSaveW, m_ionic->vec_sol_n_2()[1], MpiInfo::rankProc(), 1); |
1005 |
|
✗ |
m_linearProblem[ipb]->writeSupportVariableBackup("ion_sol_w_n_2", ionicSolToSaveW, sizeSupportVariable, 1, MpiInfo::rankProc(), m_solutionBackup[ipb].io(iMesh), m_fstransient->time, m_fstransient->iteration); |
1006 |
|
|
|
1007 |
|
✗ |
m_linearProblemBidomain[ipb]->fromVecToDoubleStar(ionicSolToSaveS, m_ionic->vec_sol_n_2()[2], MpiInfo::rankProc(), 1); |
1008 |
|
✗ |
m_linearProblem[ipb]->writeSupportVariableBackup("ion_sol_s_n_2", ionicSolToSaveS, sizeSupportVariable, 1, MpiInfo::rankProc(), m_solutionBackup[ipb].io(iMesh), m_fstransient->time, m_fstransient->iteration); |
1009 |
|
|
} |
1010 |
|
✗ |
if (ionicSolToSaveV) delete[] ionicSolToSaveV; |
1011 |
|
✗ |
if (ionicSolToSaveW) delete[] ionicSolToSaveW; |
1012 |
|
✗ |
if (ionicSolToSaveS) delete[] ionicSolToSaveS; |
1013 |
|
|
} |
1014 |
|
|
else { |
1015 |
|
✗ |
if ( FelisceParam::instance().orderBdfIonic > 0 ) { |
1016 |
|
✗ |
m_ionic->sol().getSize(&sizeSupportVariable); |
1017 |
|
✗ |
supportVariable4 = new double[sizeSupportVariable]; |
1018 |
|
✗ |
m_linearProblemBidomain[ipb]->fromVecToDoubleStar(supportVariable4, m_ionic->sol(), MpiInfo::rankProc(), 1); |
1019 |
|
✗ |
m_linearProblem[ipb]->writeSupportVariableBackup("ion_sol_n", supportVariable4, sizeSupportVariable, 1, MpiInfo::rankProc(), m_solutionBackup[ipb].io(iMesh), m_fstransient->time, m_fstransient->iteration); |
1020 |
|
|
} |
1021 |
|
|
|
1022 |
|
✗ |
if ( FelisceParam::instance().orderBdfIonic > 1 ) { |
1023 |
|
✗ |
m_ionic->sol_n_1().getSize(&sizeSupportVariable); |
1024 |
|
✗ |
supportVariable5 = new double[sizeSupportVariable]; |
1025 |
|
✗ |
m_linearProblemBidomain[ipb]->fromVecToDoubleStar(supportVariable5, m_ionic->sol_n_1(), MpiInfo::rankProc(), 1); |
1026 |
|
✗ |
m_linearProblem[ipb]->writeSupportVariableBackup("ion_sol_n_1", supportVariable5, sizeSupportVariable, 1, MpiInfo::rankProc(), m_solutionBackup[ipb].io(iMesh), m_fstransient->time, m_fstransient->iteration); |
1027 |
|
|
} |
1028 |
|
|
|
1029 |
|
✗ |
if ( FelisceParam::instance().orderBdfIonic > 2 ) { |
1030 |
|
✗ |
m_ionic->sol_n_2().getSize(&sizeSupportVariable); |
1031 |
|
✗ |
supportVariable6 = new double[sizeSupportVariable]; |
1032 |
|
✗ |
m_linearProblemBidomain[ipb]->fromVecToDoubleStar(supportVariable6, m_ionic->sol_n_2(), MpiInfo::rankProc(), 1); |
1033 |
|
✗ |
m_linearProblem[ipb]->writeSupportVariableBackup("ion_sol_n_2", supportVariable6, sizeSupportVariable, 1, MpiInfo::rankProc(), m_solutionBackup[ipb].io(iMesh), m_fstransient->time, m_fstransient->iteration); |
1034 |
|
|
} |
1035 |
|
|
} |
1036 |
|
|
|
1037 |
|
✗ |
if ( (FelisceParam::instance().typeOfIonicModel == "courtAtriaSchafVent") || (FelisceParam::instance().typeOfIonicModel == "courtAtriaMvVent") ) { |
1038 |
|
✗ |
std::string ionicVarFileName; |
1039 |
|
✗ |
for (int iCourt=0 ; iCourt<29; iCourt++) { |
1040 |
|
✗ |
supportVariableCourt.push_back(new double[sizeSupportVariable]); |
1041 |
|
✗ |
m_linearProblemBidomain[ipb]->fromVecToDoubleStar(supportVariableCourt[iCourt], m_court->ionicVariable(iCourt), MpiInfo::rankProc(), 1); |
1042 |
|
✗ |
ionicVarFileName = "court_ionic_var" + std::to_string(iCourt); |
1043 |
|
✗ |
m_linearProblem[ipb]->writeSupportVariableBackup(ionicVarFileName, supportVariableCourt[iCourt], sizeSupportVariable, 1, MpiInfo::rankProc(), m_solutionBackup[ipb].io(iMesh), m_fstransient->time, m_fstransient->iteration); |
1044 |
|
|
} |
1045 |
|
|
} |
1046 |
|
|
|
1047 |
|
✗ |
m_solutionBackup[ipb].backup(m_fstransient->time, m_mesh, MpiInfo::rankProc()); |
1048 |
|
|
|
1049 |
|
✗ |
if (supportVariable1 != nullptr) |
1050 |
|
✗ |
delete [] supportVariable1; |
1051 |
|
✗ |
if (supportVariable2 != nullptr) |
1052 |
|
✗ |
delete [] supportVariable2; |
1053 |
|
✗ |
if (supportVariable3 != nullptr) |
1054 |
|
✗ |
delete [] supportVariable3; |
1055 |
|
✗ |
if (supportVariable4 != nullptr) |
1056 |
|
✗ |
delete [] supportVariable4; |
1057 |
|
✗ |
if (supportVariable5 != nullptr) |
1058 |
|
✗ |
delete [] supportVariable5; |
1059 |
|
✗ |
if (supportVariable6 != nullptr) |
1060 |
|
✗ |
delete [] supportVariable6; |
1061 |
|
✗ |
if ( (FelisceParam::instance().typeOfIonicModel == "courtAtriaSchafVent") || (FelisceParam::instance().typeOfIonicModel == "courtAtriaMvVent") ) { |
1062 |
|
✗ |
for (int iCourt=0 ; iCourt<29; iCourt++) { |
1063 |
|
✗ |
if (supportVariableCourt[iCourt] != nullptr) |
1064 |
|
✗ |
delete [] supportVariableCourt[iCourt]; |
1065 |
|
|
} |
1066 |
|
|
} |
1067 |
|
|
} |
1068 |
|
|
} |
1069 |
|
|
} |
1070 |
|
|
|
1071 |
|
✗ |
void BidomainModel::userDefinedInitialConditions(std::size_t ipb) { |
1072 |
|
✗ |
if(FelisceParam::verbose() > 0) PetscPrintf(MpiInfo::petscComm(),"Restore solutions with previous simulation.\n"); |
1073 |
|
✗ |
double * supportVariable = nullptr; |
1074 |
|
|
// for (std::size_t ipb = 0; ipb < m_linearProblem.size(); ipb++) { |
1075 |
|
✗ |
int iterationToRead = FelisceParam::instance().restartSolutionIter + 1 - std::max(FelisceParam::instance().orderBdfEdp, FelisceParam::instance().orderBdfIonic); |
1076 |
|
|
|
1077 |
|
✗ |
int count = 0; |
1078 |
|
✗ |
supportVariable = new double[m_linearProblemBidomain[ipb]->numDof()]; |
1079 |
|
✗ |
m_initialCondition.readVariable(count, iterationToRead, supportVariable, m_linearProblemBidomain[ipb]->numDofPerUnknown(0) ); |
1080 |
|
✗ |
count++; |
1081 |
|
✗ |
for (felInt idof=0 ; idof<m_linearProblemBidomain[ipb]->numDofPerUnknown(0); idof++) { |
1082 |
|
✗ |
supportVariable[m_linearProblemBidomain[ipb]->numDofPerUnknown(0)+idof] = 0.0; |
1083 |
|
|
} |
1084 |
|
✗ |
m_linearProblemBidomain[ipb]->fromDoubleStarToVec(supportVariable,&m_U_0, m_linearProblemBidomain[ipb]->numDof()); |
1085 |
|
✗ |
if(FelisceParam::verbose() >= 40) |
1086 |
|
✗ |
m_U_0.view(); |
1087 |
|
✗ |
delete [] supportVariable; |
1088 |
|
✗ |
supportVariable = nullptr; |
1089 |
|
|
|
1090 |
|
✗ |
if ( FelisceParam::instance().orderBdfEdp > 0 ) { |
1091 |
|
✗ |
supportVariable = new double[m_linearProblemBidomain[ipb]->numDof()]; |
1092 |
|
✗ |
m_initialCondition.readVariable(count, iterationToRead, supportVariable, m_linearProblemBidomain[ipb]->numDofPerUnknown(0) ); |
1093 |
|
✗ |
for (felInt idof=0 ; idof<m_linearProblemBidomain[ipb]->numDofPerUnknown(0); idof++) { |
1094 |
|
✗ |
supportVariable[m_linearProblemBidomain[ipb]->numDofPerUnknown(0)+idof] = 0.0; |
1095 |
|
|
} |
1096 |
|
✗ |
count++; |
1097 |
|
✗ |
m_linearProblemBidomain[ipb]->fromDoubleStarToVec(supportVariable,&m_bdf_sol_n, m_linearProblemBidomain[ipb]->numDof()); |
1098 |
|
✗ |
if(FelisceParam::verbose() >= 40) |
1099 |
|
✗ |
m_bdf_sol_n.view(); |
1100 |
|
✗ |
delete [] supportVariable; |
1101 |
|
✗ |
supportVariable = nullptr; |
1102 |
|
|
} |
1103 |
|
|
|
1104 |
|
✗ |
if ( FelisceParam::instance().orderBdfEdp > 1 ) { |
1105 |
|
✗ |
supportVariable = new double[m_linearProblemBidomain[ipb]->numDof()]; |
1106 |
|
✗ |
m_initialCondition.readVariable(count, iterationToRead, supportVariable, m_linearProblemBidomain[ipb]->numDofPerUnknown(0) ); |
1107 |
|
✗ |
count++; |
1108 |
|
✗ |
for (felInt idof=0 ; idof<m_linearProblemBidomain[ipb]->numDofPerUnknown(0); idof++) { |
1109 |
|
✗ |
supportVariable[m_linearProblemBidomain[ipb]->numDofPerUnknown(0)+idof] = 0.0; |
1110 |
|
|
} |
1111 |
|
✗ |
m_linearProblemBidomain[ipb]->fromDoubleStarToVec(supportVariable,&m_bdf_sol_n_1, m_linearProblemBidomain[ipb]->numDof()); |
1112 |
|
✗ |
if(FelisceParam::verbose() >= 40) |
1113 |
|
✗ |
m_bdf_sol_n_1.view(); |
1114 |
|
✗ |
delete [] supportVariable; |
1115 |
|
✗ |
supportVariable = nullptr; |
1116 |
|
|
} |
1117 |
|
|
|
1118 |
|
✗ |
if ( FelisceParam::instance().orderBdfEdp > 2 ) { |
1119 |
|
✗ |
supportVariable = new double[m_linearProblemBidomain[ipb]->numDof()]; |
1120 |
|
✗ |
m_initialCondition.readVariable(count, iterationToRead, supportVariable, m_linearProblemBidomain[ipb]->numDofPerUnknown(0) ); |
1121 |
|
✗ |
count++; |
1122 |
|
✗ |
for (felInt idof=0 ; idof<m_linearProblemBidomain[ipb]->numDofPerUnknown(0); idof++) { |
1123 |
|
✗ |
supportVariable[m_linearProblemBidomain[ipb]->numDofPerUnknown(0)+idof] = 0.0; |
1124 |
|
|
} |
1125 |
|
✗ |
m_linearProblemBidomain[ipb]->fromDoubleStarToVec(supportVariable,&m_bdf_sol_n_2, m_linearProblemBidomain[ipb]->numDof()); |
1126 |
|
✗ |
if(FelisceParam::verbose() >= 40) |
1127 |
|
✗ |
m_bdf_sol_n_2.view(); |
1128 |
|
✗ |
delete [] supportVariable; |
1129 |
|
✗ |
supportVariable = nullptr; |
1130 |
|
|
} |
1131 |
|
|
|
1132 |
|
✗ |
if ( (FelisceParam::instance().typeOfIonicModel == "MV") || (FelisceParam::instance().typeOfIonicModel == "courtAtriaMvVent") ) { |
1133 |
|
✗ |
if ( FelisceParam::instance().orderBdfIonic > 0 ) { |
1134 |
|
✗ |
m_ionic_sol_n.resize(3); |
1135 |
|
✗ |
for (int iVar=0; iVar<3; iVar++) { |
1136 |
|
✗ |
supportVariable = new double[m_linearProblemBidomain[ipb]->numDof()]; |
1137 |
|
✗ |
m_initialCondition.readVariable(count, iterationToRead, supportVariable, m_linearProblemBidomain[ipb]->numDofPerUnknown(0) ); |
1138 |
|
✗ |
count++; |
1139 |
|
✗ |
for (felInt idof=0 ; idof<m_linearProblemBidomain[ipb]->numDofPerUnknown(0); idof++) { |
1140 |
|
✗ |
supportVariable[m_linearProblemBidomain[ipb]->numDofPerUnknown(0)+idof] = 0.0; |
1141 |
|
|
} |
1142 |
|
✗ |
m_linearProblemBidomain[ipb]->fromDoubleStarToVec(supportVariable,&m_ionic_sol_n[iVar], m_linearProblemBidomain[ipb]->numDof()); |
1143 |
|
✗ |
if(FelisceParam::verbose() >= 40) |
1144 |
|
✗ |
m_ionic_sol_n[iVar].view(); |
1145 |
|
✗ |
delete [] supportVariable; |
1146 |
|
✗ |
supportVariable = nullptr; |
1147 |
|
|
} |
1148 |
|
|
} |
1149 |
|
✗ |
if ( FelisceParam::instance().orderBdfIonic > 1 ) { |
1150 |
|
✗ |
m_ionic_sol_n_1.resize(3); |
1151 |
|
✗ |
for (int iVar=0; iVar<3; iVar++) { |
1152 |
|
✗ |
supportVariable = new double[m_linearProblemBidomain[ipb]->numDof()]; |
1153 |
|
✗ |
m_initialCondition.readVariable(count, iterationToRead, supportVariable, m_linearProblemBidomain[ipb]->numDofPerUnknown(0) ); |
1154 |
|
✗ |
count++; |
1155 |
|
✗ |
for (felInt idof=0 ; idof<m_linearProblemBidomain[ipb]->numDofPerUnknown(0); idof++) { |
1156 |
|
✗ |
supportVariable[m_linearProblemBidomain[ipb]->numDofPerUnknown(0)+idof] = 0.0; |
1157 |
|
|
} |
1158 |
|
✗ |
m_linearProblemBidomain[ipb]->fromDoubleStarToVec(supportVariable,&m_ionic_sol_n_1[iVar], m_linearProblemBidomain[ipb]->numDof()); |
1159 |
|
✗ |
if(FelisceParam::verbose() >= 40) |
1160 |
|
✗ |
m_ionic_sol_n_1[iVar].view(); |
1161 |
|
✗ |
delete [] supportVariable; |
1162 |
|
✗ |
supportVariable = nullptr; |
1163 |
|
|
} |
1164 |
|
|
} |
1165 |
|
✗ |
if ( FelisceParam::instance().orderBdfIonic > 2 ) { |
1166 |
|
✗ |
m_ionic_sol_n_2.resize(3); |
1167 |
|
✗ |
for (int iVar=0; iVar<3; iVar++) { |
1168 |
|
✗ |
supportVariable = new double[m_linearProblemBidomain[ipb]->numDof()]; |
1169 |
|
✗ |
m_initialCondition.readVariable(count, iterationToRead, supportVariable, m_linearProblemBidomain[ipb]->numDofPerUnknown(0) ); |
1170 |
|
✗ |
count++; |
1171 |
|
✗ |
for (felInt idof=0 ; idof<m_linearProblemBidomain[ipb]->numDofPerUnknown(0); idof++) { |
1172 |
|
✗ |
supportVariable[m_linearProblemBidomain[ipb]->numDofPerUnknown(0)+idof] = 0.0; |
1173 |
|
|
} |
1174 |
|
✗ |
m_linearProblemBidomain[ipb]->fromDoubleStarToVec(supportVariable,&m_ionic_sol_n_2[iVar], m_linearProblemBidomain[ipb]->numDof()); |
1175 |
|
✗ |
if(FelisceParam::verbose() >= 40) |
1176 |
|
✗ |
m_ionic_sol_n_2[iVar].view(); |
1177 |
|
✗ |
delete [] supportVariable; |
1178 |
|
✗ |
supportVariable = nullptr; |
1179 |
|
|
} |
1180 |
|
|
} |
1181 |
|
|
} |
1182 |
|
|
else { |
1183 |
|
✗ |
if ( FelisceParam::instance().orderBdfIonic > 0 ) { |
1184 |
|
✗ |
m_ionic_sol_n.resize(1); |
1185 |
|
✗ |
supportVariable = new double[m_linearProblemBidomain[ipb]->numDof()]; |
1186 |
|
✗ |
m_initialCondition.readVariable(count, iterationToRead, supportVariable, m_linearProblemBidomain[ipb]->numDofPerUnknown(0) ); |
1187 |
|
✗ |
count++; |
1188 |
|
✗ |
for (felInt idof=0 ; idof<m_linearProblemBidomain[ipb]->numDofPerUnknown(0); idof++) { |
1189 |
|
✗ |
supportVariable[m_linearProblemBidomain[ipb]->numDofPerUnknown(0)+idof] = 0.0; |
1190 |
|
|
} |
1191 |
|
✗ |
m_linearProblemBidomain[ipb]->fromDoubleStarToVec(supportVariable,&m_ionic_sol_n[0], m_linearProblemBidomain[ipb]->numDof()); |
1192 |
|
✗ |
if(FelisceParam::verbose() >= 40) |
1193 |
|
✗ |
m_ionic_sol_n[0].view(); |
1194 |
|
✗ |
delete [] supportVariable; |
1195 |
|
✗ |
supportVariable = nullptr; |
1196 |
|
|
} |
1197 |
|
|
|
1198 |
|
✗ |
if ( FelisceParam::instance().orderBdfIonic > 1 ) { |
1199 |
|
✗ |
m_ionic_sol_n_1.resize(1); |
1200 |
|
✗ |
supportVariable = new double[m_linearProblemBidomain[ipb]->numDof()]; |
1201 |
|
✗ |
m_initialCondition.readVariable(count, iterationToRead, supportVariable, m_linearProblemBidomain[ipb]->numDofPerUnknown(0) ); |
1202 |
|
✗ |
count++; |
1203 |
|
✗ |
for (felInt idof=0 ; idof<m_linearProblemBidomain[ipb]->numDofPerUnknown(0); idof++) { |
1204 |
|
✗ |
supportVariable[m_linearProblemBidomain[ipb]->numDofPerUnknown(0)+idof] = 0.0; |
1205 |
|
|
} |
1206 |
|
✗ |
m_linearProblemBidomain[ipb]->fromDoubleStarToVec(supportVariable,&m_ionic_sol_n_1[0], m_linearProblemBidomain[ipb]->numDof()); |
1207 |
|
✗ |
if(FelisceParam::verbose() >= 40) |
1208 |
|
✗ |
m_ionic_sol_n_1[0].view(); |
1209 |
|
✗ |
delete [] supportVariable; |
1210 |
|
✗ |
supportVariable = nullptr; |
1211 |
|
|
} |
1212 |
|
|
|
1213 |
|
✗ |
if ( FelisceParam::instance().orderBdfIonic > 2 ) { |
1214 |
|
✗ |
m_ionic_sol_n_2.resize(1); |
1215 |
|
✗ |
supportVariable = new double[m_linearProblemBidomain[ipb]->numDof()]; |
1216 |
|
✗ |
m_initialCondition.readVariable(count, iterationToRead, supportVariable, m_linearProblemBidomain[ipb]->numDofPerUnknown(0) ); |
1217 |
|
✗ |
count++; |
1218 |
|
✗ |
for (felInt idof=0 ; idof<m_linearProblemBidomain[ipb]->numDofPerUnknown(0); idof++) { |
1219 |
|
✗ |
supportVariable[m_linearProblemBidomain[ipb]->numDofPerUnknown(0)+idof] = 0.0; |
1220 |
|
|
} |
1221 |
|
✗ |
m_linearProblemBidomain[ipb]->fromDoubleStarToVec(supportVariable,&m_ionic_sol_n_2[0], m_linearProblemBidomain[ipb]->numDof()); |
1222 |
|
✗ |
if(FelisceParam::verbose() >= 40) |
1223 |
|
✗ |
m_ionic_sol_n_2[0].view(); |
1224 |
|
✗ |
delete [] supportVariable; |
1225 |
|
✗ |
supportVariable = nullptr; |
1226 |
|
|
} |
1227 |
|
|
} |
1228 |
|
|
|
1229 |
|
✗ |
if ( (FelisceParam::instance().typeOfIonicModel == "courtAtriaSchafVent") || (FelisceParam::instance().typeOfIonicModel == "courtAtriaMvVent") ) { |
1230 |
|
✗ |
std::string ionicVarFileName; |
1231 |
|
✗ |
m_courtVarInit.resize(29); |
1232 |
|
✗ |
for (int iCourt=0 ; iCourt<29; iCourt++) { |
1233 |
|
✗ |
supportVariable = new double[m_linearProblemBidomain[ipb]->numDof()]; |
1234 |
|
✗ |
m_initialCondition.readVariable(count, iterationToRead, supportVariable, m_linearProblemBidomain[ipb]->numDof()); |
1235 |
|
✗ |
count++; |
1236 |
|
✗ |
m_linearProblemBidomain[ipb]->fromDoubleStarToVec(supportVariable,&m_courtVarInit[iCourt], m_linearProblemBidomain[ipb]->numDof()); |
1237 |
|
✗ |
if(FelisceParam::verbose() >= 40) |
1238 |
|
✗ |
m_courtVarInit[iCourt].view(); |
1239 |
|
✗ |
delete [] supportVariable; |
1240 |
|
✗ |
supportVariable = nullptr; |
1241 |
|
|
} |
1242 |
|
|
} |
1243 |
|
|
|
1244 |
|
|
// } |
1245 |
|
|
|
1246 |
|
|
} |
1247 |
|
|
|
1248 |
|
✗ |
void BidomainModel::forward() |
1249 |
|
|
{ |
1250 |
|
|
|
1251 |
|
|
|
1252 |
|
|
//Write solution with ensight. |
1253 |
|
✗ |
BidomainModel::writeSolution(); |
1254 |
|
|
|
1255 |
|
|
//Advance time step. |
1256 |
|
✗ |
if (FelisceParam::instance().stateFilter) { |
1257 |
|
|
// clear RHS, and m_Matrix[0]. |
1258 |
|
✗ |
updateTime(FlagMatrixRHS::matrix_and_rhs); |
1259 |
|
|
} |
1260 |
|
|
else { |
1261 |
|
|
// clear only RHS, NOT m_Matrix[0]. |
1262 |
|
✗ |
updateTime(FlagMatrixRHS::only_rhs); |
1263 |
|
|
} |
1264 |
|
✗ |
printNewTimeIterationBanner(); |
1265 |
|
|
|
1266 |
|
|
|
1267 |
|
✗ |
if ( m_fstransient->iteration == 1 && !FelisceParam::instance().restartSolution) { |
1268 |
|
|
|
1269 |
|
|
//Initialization of bdf solver. |
1270 |
|
✗ |
m_bdfEdp.initialize(m_linearProblemBidomain[0]->solution()); |
1271 |
|
✗ |
m_linearProblemBidomain[0]->initializeTimeScheme(&m_bdfEdp); |
1272 |
|
|
|
1273 |
|
|
//m_extrapolate = initial solution. |
1274 |
|
✗ |
m_bdfEdp.extrapolate(m_extrapolate); |
1275 |
|
|
|
1276 |
|
|
//Initialization of Schaf solver. |
1277 |
|
✗ |
m_ionic->defineSizeAndMappingOfIonicProblem(m_linearProblemBidomain[0]->numDofLocalPerUnknown(potTransMemb), m_linearProblemBidomain[0]->mappingDofLocalToDofGlobal(potTransMemb), m_linearProblemBidomain[0]->ao()); |
1278 |
|
|
|
1279 |
|
✗ |
m_ionic->initializeExtrap(m_extrapolate); |
1280 |
|
✗ |
m_ionic->initialize(m_W_0); |
1281 |
|
✗ |
m_buildIonic = true; |
1282 |
|
✗ |
if ( (FelisceParam::instance().typeOfIonicModel == "courtAtriaSchafVent") || (FelisceParam::instance().typeOfIonicModel == "courtAtriaMvVent") ) { |
1283 |
|
✗ |
m_court->defineSizeAndMappingOfIonicProblem(m_linearProblemBidomain[0]->numDofLocalPerUnknown(potTransMemb), m_linearProblemBidomain[0]->mappingDofLocalToDofGlobal(potTransMemb), m_linearProblemBidomain[0]->ao()); |
1284 |
|
✗ |
m_court->initializeExtrap(m_extrapolate); |
1285 |
|
✗ |
m_court->initializeConcentrationsAndGateConditions(m_m, m_h, m_j, m_ao, m_io, m_ua, m_ui, m_xr, m_xs, m_d, m_f, m_fca, m_urel, m_vrel, m_wrel, m_nai, m_nao, m_cao, m_ki, m_ko, m_cai, m_naiont, m_kiont, m_caiont, m_ileak, m_iup, m_itr, m_irel, m_cmdn, m_trpn, m_nsr, m_jsr, m_csqn); |
1286 |
|
✗ |
m_buildCourt = true; |
1287 |
|
|
} |
1288 |
|
|
|
1289 |
|
|
} |
1290 |
|
✗ |
else if ( m_fstransient->iteration == 1 && FelisceParam::instance().restartSolution) { |
1291 |
|
|
// Read vectors: |
1292 |
|
|
//Initialization of bdf solver. |
1293 |
|
✗ |
if (FelisceParam::instance().orderBdfEdp == 1) { |
1294 |
|
✗ |
m_bdfEdp.initialize(m_bdf_sol_n); |
1295 |
|
|
} |
1296 |
|
✗ |
else if (FelisceParam::instance().orderBdfEdp == 2) { |
1297 |
|
✗ |
m_bdfEdp.initialize(m_bdf_sol_n_1,m_bdf_sol_n); |
1298 |
|
|
} |
1299 |
|
✗ |
else if (FelisceParam::instance().orderBdfEdp == 2) { |
1300 |
|
✗ |
m_bdfEdp.initialize(m_bdf_sol_n_2,m_bdf_sol_n_1,m_bdf_sol_n); |
1301 |
|
|
} |
1302 |
|
✗ |
m_linearProblemBidomain[0]->initializeTimeScheme(&m_bdfEdp); |
1303 |
|
|
//m_extrapolate = initial solution. |
1304 |
|
✗ |
m_bdfEdp.extrapolate(m_extrapolate); |
1305 |
|
|
|
1306 |
|
|
|
1307 |
|
|
//Initialization of ionic solver. |
1308 |
|
✗ |
m_ionic->defineSizeAndMappingOfIonicProblem(m_linearProblemBidomain[0]->numDofLocalPerUnknown(potTransMemb), m_linearProblemBidomain[0]->mappingDofLocalToDofGlobal(potTransMemb), m_linearProblemBidomain[0]->ao()); |
1309 |
|
|
|
1310 |
|
✗ |
if (FelisceParam::instance().orderBdfIonic == 1) { |
1311 |
|
✗ |
m_ionic->initialize(m_ionic_sol_n); |
1312 |
|
|
} |
1313 |
|
✗ |
else if (FelisceParam::instance().orderBdfIonic == 2) { |
1314 |
|
✗ |
m_ionic->initialize(m_ionic_sol_n, m_ionic_sol_n_1); |
1315 |
|
|
} |
1316 |
|
✗ |
else if (FelisceParam::instance().orderBdfIonic == 3) { |
1317 |
|
✗ |
m_ionic->initialize(m_ionic_sol_n, m_ionic_sol_n_1, m_ionic_sol_n_2); |
1318 |
|
|
} |
1319 |
|
✗ |
m_ionic->initializeExtrap(m_extrapolate); |
1320 |
|
✗ |
m_buildIonic = true; |
1321 |
|
|
|
1322 |
|
✗ |
if ( (FelisceParam::instance().typeOfIonicModel == "courtAtriaSchafVent") || (FelisceParam::instance().typeOfIonicModel == "courtAtriaMvVent") ) { |
1323 |
|
✗ |
m_court->defineSizeAndMappingOfIonicProblem(m_linearProblemBidomain[0]->numDofLocalPerUnknown(potTransMemb), m_linearProblemBidomain[0]->mappingDofLocalToDofGlobal(potTransMemb), m_linearProblemBidomain[0]->ao()); |
1324 |
|
✗ |
m_court->initializeExtrap(m_extrapolate); |
1325 |
|
✗ |
m_court->initializeConcentrationsAndGateConditions(m_courtVarInit[0], m_courtVarInit[1], m_courtVarInit[2], m_courtVarInit[3], m_courtVarInit[4], m_courtVarInit[5], m_courtVarInit[6], m_courtVarInit[7], m_courtVarInit[8], m_courtVarInit[9], m_courtVarInit[10], m_courtVarInit[11], m_courtVarInit[12], m_courtVarInit[13], m_courtVarInit[14], m_courtVarInit[15], m_courtVarInit[16], m_courtVarInit[17], m_courtVarInit[18], m_courtVarInit[19], m_courtVarInit[20], m_courtVarInit[21], m_courtVarInit[22], m_courtVarInit[23], m_courtVarInit[23], m_courtVarInit[23], m_courtVarInit[23], m_courtVarInit[23], m_courtVarInit[24], m_courtVarInit[25], m_courtVarInit[26], m_courtVarInit[27], m_courtVarInit[28]); |
1326 |
|
✗ |
m_buildCourt = true; |
1327 |
|
|
} |
1328 |
|
|
|
1329 |
|
|
} |
1330 |
|
|
else { |
1331 |
|
✗ |
m_bdfEdp.update(m_linearProblemBidomain[0]->solution()); |
1332 |
|
✗ |
m_bdfEdp.extrapolate(m_extrapolate); |
1333 |
|
✗ |
m_ionic->update(m_extrapolate); |
1334 |
|
|
|
1335 |
|
✗ |
if ( (FelisceParam::instance().typeOfIonicModel == "courtAtriaSchafVent") || (FelisceParam::instance().typeOfIonicModel == "courtAtriaMvVent") ) { |
1336 |
|
✗ |
m_court->updateExtrap(m_extrapolate); |
1337 |
|
|
} |
1338 |
|
|
} |
1339 |
|
|
|
1340 |
|
|
// Long QT segment during period [torsadeTimeBegin,torsadeTimeEnd] |
1341 |
|
✗ |
if (FelisceParam::instance().torsade) { |
1342 |
|
✗ |
if ( std::abs( m_fstransient->time - FelisceParam::instance().torsadeTimeBegin ) < 2*FelisceParam::instance().timeStep ) { |
1343 |
|
|
// Warning: extra-data m_linearProblemBidomain[0]->EndocardiumDistance() size has to be numDof (not number of mesh points). |
1344 |
|
✗ |
m_linearProblemBidomain[0]->evalFunctionOnDof(m_schaf->fctTauClose(), m_schaf->tauClose(),m_linearProblemBidomain[0]->EndocardiumDistance(), m_fstransient->time); |
1345 |
|
|
} |
1346 |
|
✗ |
if ( std::abs( m_fstransient->time - FelisceParam::instance().torsadeTimeEnd ) < 2*FelisceParam::instance().timeStep ) { |
1347 |
|
|
// Warning: extra-data m_linearProblemBidomain[0]->EndocardiumDistance() size has to be numDof (not number of mesh points). |
1348 |
|
✗ |
m_linearProblemBidomain[0]->evalFunctionOnDof(m_schaf->fctTauClose(), m_schaf->tauClose(),m_linearProblemBidomain[0]->EndocardiumDistance(), m_fstransient->time); |
1349 |
|
|
} |
1350 |
|
|
} |
1351 |
|
|
|
1352 |
|
|
//Solve ionic model(s) and calculate Iion. |
1353 |
|
|
|
1354 |
|
✗ |
m_ionic->computeRHS(); |
1355 |
|
✗ |
m_ionic->solveEDO(); |
1356 |
|
✗ |
m_ionic->computeIon(); |
1357 |
|
✗ |
if ( (FelisceParam::instance().typeOfIonicModel != "BCL") && (FelisceParam::instance().typeOfIonicModel != "Paci")) { |
1358 |
|
✗ |
m_ionic->updateBdf(); |
1359 |
|
|
} |
1360 |
|
|
|
1361 |
|
|
|
1362 |
|
|
|
1363 |
|
✗ |
if ( (FelisceParam::instance().typeOfIonicModel == "courtAtriaSchafVent") || (FelisceParam::instance().typeOfIonicModel == "courtAtriaMvVent") ) { |
1364 |
|
✗ |
m_court->computeIon(); |
1365 |
|
|
} |
1366 |
|
|
|
1367 |
|
✗ |
m_bdfEdp.computeRHSTime(m_fstransient->timeStep); |
1368 |
|
|
|
1369 |
|
|
// Assembling of matrices at first iteration (because of needs some coefficients not defined at iteration 0). |
1370 |
|
|
|
1371 |
|
✗ |
if ( m_fstransient->iteration == 1 ) { |
1372 |
|
✗ |
if (FelisceParam::instance().hasMoveMesh) { |
1373 |
|
✗ |
m_linearProblemBidomain[0]->readDataDisplacement(m_ios,m_fstransient->iteration); |
1374 |
|
✗ |
m_linearProblemBidomain[0]->meshLocal()->moveMesh(m_linearProblemBidomain[0]->Displacement(),1.0); |
1375 |
|
✗ |
m_linearProblemBidomain[0]->assembleMatrixRHS(MpiInfo::rankProc()); |
1376 |
|
|
} |
1377 |
|
|
else { |
1378 |
|
✗ |
if ( FelisceParam::instance().hasCoupledAtriaVent ) { |
1379 |
|
✗ |
m_linearProblemBidomain[0]->assembleMatrixRHSCurrentAndCurvilinearElement(MpiInfo::rankProc()); |
1380 |
|
|
//m_linearProblemBidomain[0]->initAndStoreMassMatrix(); |
1381 |
|
|
} |
1382 |
|
|
else { |
1383 |
|
✗ |
m_linearProblemBidomain[0]->assembleMatrixRHS(MpiInfo::rankProc()); |
1384 |
|
|
} |
1385 |
|
|
|
1386 |
|
✗ |
if(FelisceParam::instance().stateFilter) { |
1387 |
|
✗ |
m_luenbFlter.duplicateFrom(m_linearProblemBidomain[0]->matrix(1),MAT_COPY_VALUES); |
1388 |
|
✗ |
m_luenbFlter.assembly(MAT_FINAL_ASSEMBLY); |
1389 |
|
✗ |
m_luenbFlter.zeroEntries(); |
1390 |
|
|
} |
1391 |
|
|
} |
1392 |
|
|
} |
1393 |
|
|
else { |
1394 |
|
✗ |
if (FelisceParam::instance().hasMoveMesh) { |
1395 |
|
✗ |
m_linearProblemBidomain[0]->clearMatrixRHS(FlagMatrixRHS::matrix_and_rhs); |
1396 |
|
✗ |
m_linearProblemBidomain[0]->clearMatrix(1); |
1397 |
|
✗ |
m_linearProblemBidomain[0]->readDataDisplacement(m_ios,m_fstransient->iteration); |
1398 |
|
✗ |
if (FelisceParam::instance().hasMoveFibers) |
1399 |
|
✗ |
m_linearProblemBidomain[0]->updateFibers(*io(),m_fstransient->iteration*m_fstransient->timeStep); |
1400 |
|
✗ |
m_linearProblemBidomain[0]->meshLocal()->moveMesh(m_linearProblemBidomain[0]->Displacement(),1.0); |
1401 |
|
✗ |
m_linearProblemBidomain[0]->assembleMatrixRHS(MpiInfo::rankProc()); |
1402 |
|
|
|
1403 |
|
|
} |
1404 |
|
|
} |
1405 |
|
|
|
1406 |
|
|
|
1407 |
|
|
|
1408 |
|
|
//Specific operations before solve the system : calculate m_RHS std::vector used in solve(). |
1409 |
|
✗ |
postAssemblingMatrixRHS(); |
1410 |
|
|
|
1411 |
|
|
|
1412 |
|
|
//Simulation of MEA experimentation |
1413 |
|
✗ |
if(FelisceParam::instance().bc_potExtraCell == "MEA_grounds") { |
1414 |
|
✗ |
m_linearProblem[0]->finalizeEssBCTransient(); |
1415 |
|
✗ |
m_linearProblem[0]->applyBC(FelisceParam::instance().essentialBoundaryConditionsMethod, MpiInfo::rankProc()); |
1416 |
|
|
} |
1417 |
|
|
|
1418 |
|
|
|
1419 |
|
|
//Solve the system m_Matrix[0]*m_sol=m_RHS. |
1420 |
|
|
//std::cout << "MPI INFO : " << MpiInfo::rankProc() << " " << MpiInfo::numProc() << std::endl; |
1421 |
|
✗ |
m_linearProblemBidomain[0]->solve(MpiInfo::rankProc(), MpiInfo::numProc()); |
1422 |
|
|
|
1423 |
|
|
// Choose zero mean value solution. |
1424 |
|
✗ |
if(FelisceParam::instance().bc_potExtraCell != "MEA_grounds") { |
1425 |
|
✗ |
m_linearProblemBidomain[0]->removeAverageFromSolution(potExtraCell, MpiInfo::rankProc()); |
1426 |
|
|
} |
1427 |
|
|
//m_linearProblemBidomain[0]->removeAverageForPotExtraCell(); |
1428 |
|
|
|
1429 |
|
|
|
1430 |
|
|
} |
1431 |
|
|
|
1432 |
|
✗ |
void BidomainModel::forwardROM() { |
1433 |
|
|
//Write solution with ensight. |
1434 |
|
✗ |
BidomainModel::writeSolution(); |
1435 |
|
|
//Advance time step. |
1436 |
|
✗ |
updateTime(FlagMatrixRHS::only_rhs); |
1437 |
|
✗ |
printNewTimeIterationBanner(); |
1438 |
|
|
|
1439 |
|
✗ |
if ( m_fstransient->iteration == 1) { |
1440 |
|
|
//Initialization of bdf solver. |
1441 |
|
✗ |
m_bdfEdp.initialize(m_linearProblemBidomain[0]->solution()); |
1442 |
|
✗ |
m_linearProblemBidomain[0]->initializeTimeScheme(&m_bdfEdp); |
1443 |
|
|
//m_extrapolate = initial solution. |
1444 |
|
✗ |
m_bdfEdp.extrapolate(m_extrapolate); |
1445 |
|
|
//Initialization of Schaf solver. |
1446 |
|
✗ |
m_ionic->defineSizeAndMappingOfIonicProblem(m_linearProblemBidomain[0]->numDofLocalPerUnknown(potTransMemb), m_linearProblemBidomain[0]->mappingDofLocalToDofGlobal(potTransMemb), m_linearProblemBidomain[0]->ao()); |
1447 |
|
✗ |
m_ionic->initializeExtrap(m_extrapolate); |
1448 |
|
✗ |
m_ionic->initialize(m_W_0[0]); |
1449 |
|
✗ |
m_buildIonic = true; |
1450 |
|
|
} else { |
1451 |
|
✗ |
m_bdfEdp.update(m_linearProblemBidomain[0]->solution()); |
1452 |
|
✗ |
m_bdfEdp.extrapolate(m_extrapolate); |
1453 |
|
✗ |
m_ionic->update(m_extrapolate); |
1454 |
|
|
} |
1455 |
|
|
|
1456 |
|
|
//Solve Mithcell and Schaeffer model and calculate Iion. |
1457 |
|
✗ |
m_ionic->computeRHS(); |
1458 |
|
✗ |
m_ionic->solveEDO(); |
1459 |
|
✗ |
m_ionic->computeIon(); |
1460 |
|
✗ |
m_ionic->updateBdf(); |
1461 |
|
|
|
1462 |
|
✗ |
m_bdfEdp.computeRHSTime(m_fstransient->timeStep); |
1463 |
|
|
|
1464 |
|
|
// Assembling of matrices at first iteration (because of needs some coefficients not defined at iteration 0). |
1465 |
|
✗ |
if ( m_fstransient->iteration == 1 ) { |
1466 |
|
✗ |
m_linearProblemBidomain[0]->assembleMatrixRHS(MpiInfo::rankProc()); |
1467 |
|
✗ |
if(FelisceParam::instance().stateFilter) { |
1468 |
|
✗ |
m_luenbFlter.duplicateFrom(m_linearProblemBidomain[0]->matrix(1),MAT_COPY_VALUES); |
1469 |
|
✗ |
m_luenbFlter.assembly(MAT_FINAL_ASSEMBLY); |
1470 |
|
✗ |
m_luenbFlter.zeroEntries(); |
1471 |
|
|
} |
1472 |
|
|
} |
1473 |
|
|
//Specific operations before solve the system : calculate m_RHS std::vector used in solve(). |
1474 |
|
✗ |
postAssemblingMatrixRHS(); |
1475 |
|
|
//Solve the system m_Matrix[0]*m_sol=m_RHS. |
1476 |
|
✗ |
m_rom->solveROM(m_linearProblemBidomain[0]->solution()); |
1477 |
|
|
|
1478 |
|
✗ |
m_linearProblemBidomain[0]->removeAverageFromSolution(potExtraCell, MpiInfo::rankProc()); |
1479 |
|
|
|
1480 |
|
|
} |
1481 |
|
|
|
1482 |
|
✗ |
MVSolver* BidomainModel::getMVSolver() { |
1483 |
|
✗ |
return m_mv; |
1484 |
|
|
} |
1485 |
|
|
#ifdef FELISCE_WITH_CVGRAPH |
1486 |
|
✗ |
void BidomainModel::forwardLight() { |
1487 |
|
|
|
1488 |
|
|
//Advance time step. |
1489 |
|
✗ |
std::stringstream msg; |
1490 |
|
✗ |
msg<<"Starting a new time step, preparing forward"<<std::endl; |
1491 |
|
✗ |
PetscPrintf(MpiInfo::petscComm(),"%s",msg.str().c_str()); |
1492 |
|
✗ |
BidomainModel::writeSolution(); |
1493 |
|
✗ |
updateTime(FlagMatrixRHS::only_rhs); |
1494 |
|
✗ |
printNewTimeIterationBanner(); |
1495 |
|
|
|
1496 |
|
|
// Initialize a bunch of stuff. Should be moved elsewhere. |
1497 |
|
✗ |
if (m_fstransient->iteration == 1) { |
1498 |
|
|
//Initialization of bdf solver. |
1499 |
|
✗ |
m_bdfEdp.initialize(m_linearProblemBidomain[0]->solution()); |
1500 |
|
✗ |
m_linearProblemBidomain[0]->initializeTimeScheme(&m_bdfEdp); |
1501 |
|
|
//m_extrapolate = initial solution. |
1502 |
|
✗ |
m_bdfEdp.extrapolate(m_extrapolate); |
1503 |
|
|
//Initialization of ionic solver. |
1504 |
|
✗ |
m_ionic->defineSizeAndMappingOfIonicProblem(m_linearProblemBidomain[0]->numDofLocalPerUnknown(potTransMemb), m_linearProblemBidomain[0]->mappingDofLocalToDofGlobal(potTransMemb), m_linearProblemBidomain[0]->ao()); |
1505 |
|
✗ |
m_ionic->initializeExtrap(m_extrapolate); |
1506 |
|
✗ |
m_ionic->initialize(m_W_0); |
1507 |
|
✗ |
m_buildIonic = true; |
1508 |
|
|
} else { |
1509 |
|
✗ |
m_bdfEdp.update(m_linearProblemBidomain[0]->solution()); |
1510 |
|
✗ |
m_bdfEdp.extrapolate(m_extrapolate); |
1511 |
|
✗ |
m_ionic->update(m_extrapolate); |
1512 |
|
|
} |
1513 |
|
|
|
1514 |
|
✗ |
m_ionic->computeRHS(); |
1515 |
|
✗ |
m_ionic->solveEDO(); |
1516 |
|
✗ |
m_ionic->computeIon(); |
1517 |
|
✗ |
m_ionic->updateBdf(); |
1518 |
|
✗ |
m_bdfEdp.computeRHSTime(m_fstransient->timeStep); |
1519 |
|
|
|
1520 |
|
✗ |
if (m_fstransient->iteration == 1){ |
1521 |
|
✗ |
m_linearProblemBidomain[0]->assembleMatrixRHSCurrentAndCurvilinearElement(MpiInfo::rankProc()); |
1522 |
|
|
//m_linearProblem[0]->createAndCopyMatrixRHSWithoutBC(FlagMatrixRHS::only_matrix); |
1523 |
|
✗ |
m_linearProblemBidomain[0]->initAndStoreMassMatrix(); |
1524 |
|
|
} |
1525 |
|
|
|
1526 |
|
|
//Specific operations before solve the system : calculate m_RHS std::vector used in solve(). |
1527 |
|
✗ |
postAssemblingMatrixRHS(); |
1528 |
|
|
|
1529 |
|
|
// Apply the BC |
1530 |
|
✗ |
m_linearProblem[0]->finalizeEssBCConstantInT(); |
1531 |
|
✗ |
m_linearProblem[0]->finalizeEssBCTransient(); |
1532 |
|
✗ |
if (m_fstransient->iteration == 1){ |
1533 |
|
✗ |
m_linearProblem[0]->applyBC(FelisceParam::instance().essentialBoundaryConditionsMethod, MpiInfo::rankProc(),FlagMatrixRHS::matrix_and_rhs, FlagMatrixRHS::matrix_and_rhs, 0, true, ApplyNaturalBoundaryConditions::no); |
1534 |
|
|
} else { |
1535 |
|
✗ |
m_linearProblem[0]->applyBC(FelisceParam::instance().essentialBoundaryConditionsMethod, MpiInfo::rankProc(),FlagMatrixRHS::only_rhs, FlagMatrixRHS::only_rhs, 0, true, ApplyNaturalBoundaryConditions::no); |
1536 |
|
|
} |
1537 |
|
|
|
1538 |
|
|
|
1539 |
|
|
//Solve the system m_Matrix[0]*m_sol=m_RHS. |
1540 |
|
✗ |
m_linearProblemBidomain[0]->solve(MpiInfo::rankProc(), MpiInfo::numProc()); |
1541 |
|
|
// Remove average |
1542 |
|
✗ |
m_linearProblemBidomain[0]->removeAverageForPotExtraCell(); |
1543 |
|
|
//m_linearProblemBidomain[0]->removeAverageFromSolution(potExtraCell, MpiInfo::rankProc()); |
1544 |
|
|
|
1545 |
|
|
} |
1546 |
|
|
|
1547 |
|
|
|
1548 |
|
✗ |
void BidomainModel::startIterationCVG() { |
1549 |
|
✗ |
cvgMainSlave* slave=m_linearProblem[0]->slave(); |
1550 |
|
✗ |
if (slave==nullptr) { |
1551 |
|
✗ |
FEL_ERROR("slave not initialized"); |
1552 |
|
|
} |
1553 |
|
|
//Advance time step. |
1554 |
|
✗ |
if( slave->newTimeStep() ) { |
1555 |
|
✗ |
m_kount = 0; |
1556 |
|
✗ |
std::stringstream msg; |
1557 |
|
✗ |
msg<<"Starting a new time step, preparing forward"<<std::endl; |
1558 |
|
✗ |
PetscPrintf(MpiInfo::petscComm(),"%s",msg.str().c_str()); |
1559 |
|
✗ |
BidomainModel::writeSolution(); |
1560 |
|
✗ |
updateTime(FlagMatrixRHS::only_rhs); |
1561 |
|
✗ |
printNewTimeIterationBanner(); |
1562 |
|
✗ |
} else if( slave->redoTimeStep() ){ |
1563 |
|
✗ |
m_kount++; |
1564 |
|
✗ |
std::stringstream msg; |
1565 |
|
✗ |
msg << "......................................" << std::endl; |
1566 |
|
✗ |
msg << "Redo time step at time t= " << m_fstransient->time << std::endl; |
1567 |
|
✗ |
PetscPrintf(MpiInfo::petscComm(),"%s",msg.str().c_str()); |
1568 |
|
✗ |
m_linearProblemBidomain[0]->clearMatrixRHS(FlagMatrixRHS::only_rhs); |
1569 |
|
✗ |
} else { |
1570 |
|
✗ |
slave->msgInfo().print(1); |
1571 |
|
✗ |
FEL_ERROR("Error: strange timeStatus"); |
1572 |
|
|
} |
1573 |
|
|
|
1574 |
|
|
|
1575 |
|
|
// Initialize a bunch of stuff. Should be moved elsewhere. |
1576 |
|
✗ |
if ( ( m_fstransient->iteration == 1) && slave->newTimeStep()) { |
1577 |
|
|
//Initialization of bdf solver. |
1578 |
|
✗ |
m_bdfEdp.initialize(m_linearProblemBidomain[0]->solution()); |
1579 |
|
✗ |
m_linearProblemBidomain[0]->initializeTimeScheme(&m_bdfEdp); |
1580 |
|
|
//m_extrapolate = initial solution. |
1581 |
|
✗ |
m_bdfEdp.extrapolate(m_extrapolate); |
1582 |
|
|
//Initialization of ionic solver. |
1583 |
|
✗ |
m_ionic->defineSizeAndMappingOfIonicProblem(m_linearProblemBidomain[0]->numDofLocalPerUnknown(potTransMemb), m_linearProblemBidomain[0]->mappingDofLocalToDofGlobal(potTransMemb), m_linearProblemBidomain[0]->ao()); |
1584 |
|
✗ |
m_ionic->initializeExtrap(m_extrapolate); |
1585 |
|
✗ |
m_ionic->initialize(m_W_0); |
1586 |
|
✗ |
m_buildIonic = true; |
1587 |
|
✗ |
} else if (slave->newTimeStep()) { |
1588 |
|
✗ |
m_bdfEdp.update(m_linearProblemBidomain[0]->solution()); |
1589 |
|
✗ |
m_bdfEdp.extrapolate(m_extrapolate); |
1590 |
|
✗ |
m_ionic->update(m_extrapolate); |
1591 |
|
|
} |
1592 |
|
|
|
1593 |
|
✗ |
if (slave->newTimeStep()){ |
1594 |
|
✗ |
m_ionic->computeRHS(); |
1595 |
|
✗ |
m_ionic->solveEDO(); |
1596 |
|
✗ |
m_ionic->computeIon(); |
1597 |
|
✗ |
m_ionic->updateBdf(); |
1598 |
|
✗ |
m_bdfEdp.computeRHSTime(m_fstransient->timeStep); |
1599 |
|
|
} |
1600 |
|
|
|
1601 |
|
|
|
1602 |
|
✗ |
if ( (m_fstransient->iteration == 1) && slave->newTimeStep() ){ |
1603 |
|
✗ |
m_linearProblemBidomain[0]->assembleMatrixRHSCurrentAndCurvilinearElement(MpiInfo::rankProc()); |
1604 |
|
✗ |
m_linearProblem[0]->createAndCopyMatrixRHSWithoutBC(FlagMatrixRHS::only_matrix); |
1605 |
|
✗ |
m_linearProblemBidomain[0]->initAndStoreMassMatrix(); |
1606 |
|
|
} |
1607 |
|
|
|
1608 |
|
|
//Specific operations before solve the system : calculate m_RHS std::vector used in solve(). |
1609 |
|
✗ |
postAssemblingMatrixRHS(); |
1610 |
|
✗ |
m_linearProblemBidomain[0]->applyUserDefinedStimulation(m_fstransient->time); |
1611 |
|
✗ |
m_linearProblem[0]->createAndCopyMatrixRHSWithoutBC(FlagMatrixRHS::only_rhs); |
1612 |
|
|
|
1613 |
|
|
// Apply the BC |
1614 |
|
✗ |
m_linearProblem[0]->finalizeEssBCConstantInT(); |
1615 |
|
✗ |
m_linearProblem[0]->finalizeEssBCTransient(); |
1616 |
|
✗ |
if (( m_fstransient->iteration == 1) && slave->newTimeStep() ){ |
1617 |
|
|
//std::cout << "BidomainModel::startIterationCVG(): FIRST TIME STEP" <<std::endl; |
1618 |
|
✗ |
if (slave->thereIsAtLeastOneRobinCondition()){ |
1619 |
|
|
//std::cout << " CVG Robin: we apply the robin condition in the MATRIX ONLY. The RHS parts are handled in cvgAddExternalResidualToRHS() and cvgAddRobinToRHS()." << std::endl; |
1620 |
|
✗ |
m_linearProblem[0]->applyBC(FelisceParam::instance().essentialBoundaryConditionsMethod,MpiInfo::rankProc(),FlagMatrixRHS::matrix_and_rhs, FlagMatrixRHS::only_matrix, 0, true, ApplyNaturalBoundaryConditions::yes); |
1621 |
|
|
} else { |
1622 |
|
|
//std::cout << " CVG Neumann: Neumann from other compartment handled in cvgAddExternalResidualToRHS()" << std::endl; |
1623 |
|
✗ |
m_linearProblem[0]->applyBC(FelisceParam::instance().essentialBoundaryConditionsMethod, MpiInfo::rankProc(),FlagMatrixRHS::matrix_and_rhs, FlagMatrixRHS::matrix_and_rhs, 0, true, ApplyNaturalBoundaryConditions::no); |
1624 |
|
|
} |
1625 |
|
|
} else { |
1626 |
|
✗ |
m_linearProblem[0]->applyBC(FelisceParam::instance().essentialBoundaryConditionsMethod, MpiInfo::rankProc(),FlagMatrixRHS::only_rhs, FlagMatrixRHS::only_rhs, 0, true, ApplyNaturalBoundaryConditions::no); |
1627 |
|
|
} |
1628 |
|
|
|
1629 |
|
✗ |
if (slave->thereIsAtLeastOneRobinCondition()){ |
1630 |
|
|
//std::cout << " CVG Robin: Add remaining term in RHS: beta_torso*massBD*u_torso" << std::endl; |
1631 |
|
✗ |
m_linearProblemBidomain[0]->cvgAddRobinToRHS(); |
1632 |
|
|
} |
1633 |
|
|
// Add current from other compartment |
1634 |
|
|
//m_linearProblemBidomain[0]->addCurrentToRHS(); |
1635 |
|
|
//Solve the system m_Matrix[0]*m_sol=m_RHS. |
1636 |
|
|
double tic, toc, tsolve; |
1637 |
|
✗ |
tic = clock(); |
1638 |
|
✗ |
m_linearProblemBidomain[0]->solve(MpiInfo::rankProc(), MpiInfo::numProc()); |
1639 |
|
✗ |
toc = clock(); |
1640 |
|
✗ |
tsolve = (toc-tic)/CLOCKS_PER_SEC; |
1641 |
|
|
// Remove average |
1642 |
|
|
//double trmvav, trmvav2; |
1643 |
|
|
//m_linearProblemBidomain[0]->removeAverageFromSolution(potExtraCell, MpiInfo::rankProc()); |
1644 |
|
|
//toc = clock(); |
1645 |
|
|
//trmvav = (toc-tic)/CLOCKS_PER_SEC; |
1646 |
|
|
//tic = clock(); |
1647 |
|
|
|
1648 |
|
✗ |
m_linearProblemBidomain[0]->removeAverageForPotExtraCell(); |
1649 |
|
|
|
1650 |
|
|
//toc = clock(); |
1651 |
|
|
//trmvav = (toc-tic)/CLOCKS_PER_SEC; |
1652 |
|
✗ |
std::cout << "Tsolve = " << tsolve << std::endl;// ",Trmvav2=" << trmvav2<< std::endl; |
1653 |
|
|
} |
1654 |
|
|
#endif |
1655 |
|
|
|
1656 |
|
|
} |
1657 |
|
|
|