GCC Code Coverage Report

Directory: ./
File: Solver/BCLSolver.cpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 0 316 0.0%
Branches: 0 186 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: F. Raphel
13 //
14
15 // System includes
16
17 // External includes
18
19 // Project includes
20 #include "Solver/BCLSolver.hpp"
21
22 #ifdef FELISCE_WITH_SUNDIALS
23
24 //For CVODE
25 #define Ith(v,i) NV_Ith_S(v,i) /* Ith numbers components 1..NEQ */
26 #define IJth(A,i,j) DENSE_ELEM(A,i-1,j-1) /* IJth numbers rows,cols 1..NEQ */
27
28 static constexpr double RTOL = RCONST(1.0e-8);
29 static constexpr double ATOL = RCONST(1.0e-4); /* std::vector absolute tolerance components */
30
31 #endif
32
33 namespace felisce {
34 ✗ BCLSolver::BCLSolver(FelisceTransient::Pointer fstransient):
35 ✗ IonicSolver(fstransient) {
36 ✗ m_nComp = 4;
37
38 ✗ if (FelisceParam::instance().CellsType != "hetero") {
39 ✗ tauhp.resize(1);
40 ✗ tauhm.resize(1);
41 ✗ taufp.resize(1);
42 ✗ taufm.resize(1);
43 ✗ taurp.resize(1);
44 ✗ taurm.resize(1);
45 ✗ tausp.resize(1);
46 ✗ tausm.resize(1);
47 ✗ gfi.resize(1);
48 ✗ Vfi.resize(1);
49 ✗ gso.resize(1);
50 ✗ gsi.resize(1);
51 ✗ beta1.resize(1);
52 ✗ beta2.resize(1);
53 ✗ V1.resize(1);
54 ✗ V2.resize(1);
55 ✗ gto.resize(1);
56 ✗ Vc.resize(1);
57 ✗ Vs.resize(1);
58
59 ✗ if (FelisceParam::instance().CellsType == "endo") {
60 ✗ tauhp[0] = 10.8;
61 ✗ tauhm[0] = 10.8;
62 ✗ taufp[0] = 355.;
63 ✗ taufm[0] = 52.3;
64 ✗ taurp[0] = 7.54;
65 ✗ taurm[0] = 6.07;
66 ✗ tausp[0] = 29.1;
67 ✗ tausm[0] = 10.4;
68 ✗ gfi[0] = 1.72;
69 ✗ Vfi[0] = 1.24;
70 ✗ gso[0] = 0.00891;
71 ✗ gsi[0] = 0.414;
72 ✗ beta1[0] = 22.8;
73 ✗ beta2[0] = 2.95;
74 ✗ V1[0] = 0.522;
75 ✗ V2[0] = 0.596;
76 ✗ gto[0] = 0.300;
77 ✗ Vc[0] = 0.130;
78 ✗ Vs[0] = 0.6;
79 ✗ } else if (FelisceParam::instance().CellsType == "M") {
80 ✗ tauhp[0] = 11.3;
81 ✗ tauhm[0] = 1.88;
82 ✗ taufp[0] = 101.;
83 ✗ taufm[0] = 228.;
84 ✗ taurp[0] = 2.15;
85 ✗ taurm[0] = 0.371;
86 ✗ tausp[0] = 4.67;
87 ✗ tausm[0] = 1.75;
88 ✗ gfi[0] = 2.62;
89 ✗ Vfi[0] = 1.60;
90 ✗ gso[0] = 0.0278;
91 ✗ gsi[0] = 0.103;
92 ✗ beta1[0] = 27.1;
93 ✗ beta2[0] = 6.12;
94 ✗ V1[0] = 0.668;
95 ✗ V2[0] = 1.08;
96 ✗ gto[0] = 1.36;
97 ✗ Vc[0] = 0.130;
98 ✗ Vs[0] = 0.3;
99 ✗ } else if (FelisceParam::instance().CellsType == "epi") {
100 ✗ tauhp[0] = 17.9;
101 ✗ tauhm[0] = 11.4;
102 ✗ taufp[0] = 123.;
103 ✗ taufm[0] = 183.;
104 ✗ taurp[0] = 2.51;
105 ✗ taurm[0] = 2.00;
106 ✗ tausp[0] = 57.0;
107 ✗ tausm[0] = 10.6;
108 ✗ gfi[0] = 4.00;
109 ✗ Vfi[0] = 1.46;
110 ✗ gso[0] = 0.0161;
111 ✗ gsi[0] = 0.176;
112 ✗ beta1[0] = 3.99;
113 ✗ beta2[0] = 1.56;
114 ✗ V1[0] = 0.529;
115 ✗ V2[0] = 0.386;
116 ✗ gto[0] = 2.10;
117 ✗ Vc[0] = 0.130;
118 ✗ Vs[0] = 0.3;
119 ✗ } else if (FelisceParam::instance().CellsType == "LRd") {
120 ✗ tauhp[0] = 25.8;
121 ✗ tauhm[0] = 0.950;
122 ✗ taufp[0] = 488.;
123 ✗ taufm[0] = 25.4;
124 ✗ taurp[0] = 5.10;
125 ✗ taurm[0] = 13.1;
126 ✗ tausp[0] = 300.;
127 ✗ tausm[0] = 7.11;
128 ✗ gfi[0] = 10.0;
129 ✗ Vfi[0] = 1.20;
130 ✗ gso[0] = 0.0316;
131 ✗ gsi[0] = 2.32;
132 ✗ beta1[0] = 6.90;
133 ✗ beta2[0] = 6.36;
134 ✗ V1[0] = 0.453;
135 ✗ V2[0] = 0.828;
136 ✗ gto[0] = 2.50;
137 ✗ Vc[0] = 0.380;
138 ✗ Vs[0] = 0.6;
139 ✗ } else if (FelisceParam::instance().CellsType == "TNNP") {
140 ✗ tauhp[0] = 9035;
141 ✗ tauhm[0] = 6.61;
142 ✗ taufp[0] = 43.1;
143 ✗ taufm[0] = 181.;
144 ✗ taurp[0] = 13.5;
145 ✗ taurm[0] = 2.20;
146 ✗ tausp[0] = 99.9;
147 ✗ tausm[0] = 4.34;
148 ✗ gfi[0] = 14.0;
149 ✗ Vfi[0] = 1.18;
150 ✗ gso[0] = 0.0498;
151 ✗ gsi[0] = 0.138;
152 ✗ beta1[0] = 11.8;
153 ✗ beta2[0] = 8.05;
154 ✗ V1[0] = 0.200;
155 ✗ V2[0] = 1.02;
156 ✗ gto[0] = 9.82;
157 ✗ Vc[0] = 0.350;
158 ✗ Vs[0] = 0.6;
159 }
160
161 } else {
162
163 ✗ tauhp.resize(4);
164 ✗ tauhm.resize(4);
165 ✗ taufp.resize(4);
166 ✗ taufm.resize(4);
167 ✗ taurp.resize(4);
168 ✗ taurm.resize(4);
169 ✗ tausp.resize(4);
170 ✗ tausm.resize(4);
171 ✗ gfi.resize(4);
172 ✗ Vfi.resize(4);
173 ✗ gso.resize(4);
174 ✗ gsi.resize(4);
175 ✗ beta1.resize(4);
176 ✗ beta2.resize(4);
177 ✗ V1.resize(4);
178 ✗ V2.resize(4);
179 ✗ gto.resize(4);
180 ✗ Vc.resize(4);
181 ✗ Vs.resize(4);
182
183 // RV
184 ✗ tauhp[0] = 17.9;
185 ✗ tauhm[0] = 11.4;
186 ✗ taufp[0] = 123.;
187 ✗ taufm[0] = 183.;
188 ✗ taurp[0] = 2.51;
189 ✗ taurm[0] = 2.00;
190 ✗ tausp[0] = 57.0;
191 ✗ tausm[0] = 10.6;
192 ✗ gfi[0] = 4.00;
193 ✗ Vfi[0] = 1.46;
194 ✗ gso[0] = 0.0161;
195 ✗ gsi[0] = 0.176;
196 ✗ beta1[0] = 3.99;
197 ✗ beta2[0] = 1.56;
198 ✗ V1[0] = 0.529;
199 ✗ V2[0] = 0.386;
200 ✗ gto[0] = 2.10;
201 ✗ Vc[0] = 0.130;
202 ✗ Vs[0] = 0.3;
203
204 // Endo
205 ✗ tauhp[1] = 10.8;
206 ✗ tauhm[1] = 10.8;
207 ✗ taufp[1] = 355.;
208 ✗ taufm[1] = 52.3;
209 ✗ taurp[1] = 7.54;
210 ✗ taurm[1] = 6.07;
211 ✗ tausp[1] = 29.1;
212 ✗ tausm[1] = 10.4;
213 ✗ gfi[1] = 1.72;
214 ✗ Vfi[1] = 1.24;
215 ✗ gso[1] = 0.00891;
216 ✗ gsi[1] = 0.414;
217 ✗ beta1[1] = 22.8;
218 ✗ beta2[1] = 2.95;
219 ✗ V1[1] = 0.522;
220 ✗ V2[1] = 0.596;
221 ✗ gto[1] = 0.300;
222 ✗ Vc[1] = 0.130;
223 ✗ Vs[1] = 0.6;
224
225 // M-Cells
226 ✗ tauhp[2] = 11.3;
227 ✗ tauhm[2] = 1.88;
228 ✗ taufp[2] = 101.;
229 ✗ taufm[2] = 228.;
230 ✗ taurp[2] = 2.15;
231 ✗ taurm[2] = 0.371;
232 ✗ tausp[2] = 4.67;
233 ✗ tausm[2] = 1.75;
234 ✗ gfi[2] = 2.62;
235 ✗ Vfi[2] = 1.60;
236 ✗ gso[2] = 0.0278;
237 ✗ gsi[2] = 0.103;
238 ✗ beta1[2] = 27.1;
239 ✗ beta2[2] = 6.12;
240 ✗ V1[2] = 0.668;
241 ✗ V2[2] = 1.08;
242 ✗ gto[2] = 1.36;
243 ✗ Vc[2] = 0.130;
244 ✗ Vs[2] = 0.3;
245
246 // Epi
247 ✗ tauhp[3] = 17.9;
248 ✗ tauhm[3] = 11.4;
249 ✗ taufp[3] = 123.;
250 ✗ taufm[3] = 183.;
251 ✗ taurp[3] = 2.51;
252 ✗ taurm[3] = 2.00;
253 ✗ tausp[3] = 57.0;
254 ✗ tausm[3] = 10.6;
255 ✗ gfi[3] = 4.00;
256 ✗ Vfi[3] = 1.46;
257 ✗ gso[3] = 0.0161;
258 ✗ gsi[3] = 0.176;
259 ✗ beta1[3] = 3.99;
260 ✗ beta2[3] = 1.56;
261 ✗ V1[3] = 0.529;
262 ✗ V2[3] = 0.386;
263 ✗ gto[3] = 2.10;
264 ✗ Vc[3] = 0.130;
265 ✗ Vs[3] = 0.3;
266 }
267 }
268
269 ✗ BCLSolver::~BCLSolver()
270 = default;
271
272 ✗ void BCLSolver::initialize(std::vector<PetscVector>& sol_0) {
273 #ifdef FELISCE_WITH_SUNDIALS
274 int flag;//Informations for CVODE (if flag==0: OK)
275 #else
276 IGNORE_UNUSED_ARGUMENT(sol_0);
277 #endif
278 ✗ m_userData.resize(m_size);
279 ✗ m_cvode_mem.resize(m_size);
280 #ifdef FELISCE_WITH_SUNDIALS
281 m_initvalue.resize(m_size);
282 m_abstol.resize(m_size);
283 #endif
284
285 double value_uExtrap;
286 ✗ double Vmin = FelisceParam::instance().vMin;
287 ✗ double Vmax = FelisceParam::instance().vMax;
288
289 ✗ int iType = 0;
290
291 felInt pos;
292 ✗ for (felInt i = 0; i < m_size; i++) {
293 // Set m_userData.u[i] with transmembrane potential extrapolate (rescaled on [0,1])
294 ✗ ISLocalToGlobalMappingApply(m_localDofToGlobalDof,1,&i,&pos);
295 ✗ m_uExtrap.getValues(1,&pos,&value_uExtrap);
296 ✗ value_uExtrap = (value_uExtrap-Vmin)/(Vmax-Vmin);
297 ✗ m_userData[i].u = value_uExtrap;
298
299 ✗ if (FelisceParam::instance().CellsType == "hetero")
300 ✗ iType = m_cellType[pos];
301
302 ✗ m_userData[i].iType = iType;
303 ✗ m_userData[i].pos = i;
304
305 ✗ m_cvode_mem[i] = nullptr;
306 #ifdef FELISCE_WITH_SUNDIALS
307 m_initvalue[i] = N_VNew_Serial(m_nComp);
308 m_abstol[i] = N_VNew_Serial(m_nComp);
309
310 double valuem_solEDO;//0=h, 1=f, 2=r and 3=s
311 double w0value;
312 for (std::size_t j=0; j<m_nComp; j++) {
313 sol_0[j].getValues(1,&pos,&w0value);
314 Ith(m_initvalue[i],j) = w0value;//h
315 Ith(m_abstol[i],j) = ATOL;
316 }
317
318 m_cvode_mem[i] = CVodeCreate(CV_BDF,CV_NEWTON);//CV_BDF or CV_ADAMS and CV_NEWTON or CV_FUNCTIONAL || (CV_BDF,CV_NEWTON) for stiff pb
319 m_data = &m_userData[i];
320
321 flag = CVodeSetUserData(m_cvode_mem[i],m_data);
322 flag = CVodeInit(m_cvode_mem[i],M_f,0.0,m_initvalue[i]);
323 flag = CVodeSVtolerances(m_cvode_mem[i],RTOL,m_abstol[i]);
324 flag = CVDense(m_cvode_mem[i],m_nComp);
325
326 //CVODE Options:
327 //
328 flag = CVodeSetMaxOrd(m_cvode_mem[i],5);//default: 5, order max of BDF
329 flag = CVodeSetMaxNumSteps(m_cvode_mem[i],500);//default: 500
330 flag = CVodeSetMaxStep(m_cvode_mem[i],FelisceParam::instance().timeStep);//0.1 is the step interpolation
331 #endif
332 }
333 }
334
335
336 ✗ void BCLSolver::solveEDO() {
337
338 #ifdef FELISCE_WITH_SUNDIALS
339 int flag;//Informations for CVODE (if flag==0: OK)
340 #endif
341 double value_uExtrap;
342 ✗ double Vmin = FelisceParam::instance().vMin;
343 ✗ double Vmax = FelisceParam::instance().vMax;
344
345
346 felInt pos;
347 ✗ for (felInt i = 0; i < m_size; i++) {
348 // Set m_userData.u[i] with transmembrane potential extrapolate (rescaled on [0,1])
349 ✗ ISLocalToGlobalMappingApply(m_localDofToGlobalDof,1,&i,&pos);
350 ✗ m_uExtrap.getValues(1,&pos,&value_uExtrap);
351 ✗ value_uExtrap = (value_uExtrap-Vmin)/(Vmax-Vmin);
352 ✗ m_userData[i].u = value_uExtrap;
353
354 #ifdef FELISCE_WITH_SUNDIALS
355 int iType = 0;
356 if (FelisceParam::instance().CellsType == "hetero")
357 iType = m_cellType[pos];
358
359 double time = m_fstransient->time;
360 double valuem_solEDO;//0=h, 1=f, 2=r and 3=s
361 realtype tout;
362 if (iType < 4) {
363 flag = CVode(m_cvode_mem[i], time, m_initvalue[i], &tout, CV_NORMAL);
364 if ( flag != CV_SUCCESS ) {
365 FEL_ERROR("CVODE solver failed ");
366 }
367
368 for(std::size_t j=0; j<m_nComp; j++) {
369 valuem_solEDO = Ith(m_initvalue[i],j);
370 m_vecSolEDO[j].setValue(pos,valuem_solEDO, INSERT_VALUES);
371 }
372 } else {
373 valuem_solEDO = 0.;
374 for(std::size_t j=0; j<m_nComp; j++) {
375 m_vecSolEDO[j].setValue(pos,valuem_solEDO, INSERT_VALUES);
376 }
377 }
378 #endif
379 }
380 ✗ for(std::size_t j=0; j<m_nComp; j++) {
381 ✗ m_vecSolEDO[j].assembly();
382 }
383
384 }
385
386 ✗ void BCLSolver::computeIon() {
387 ✗ std::cout << " void BCLSolver::computeIon()" << std::endl;
388 double value_uExtrap;
389 ✗ std::vector<double> valuem_solEDO;
390 ✗ valuem_solEDO.resize(4);
391 double value_ion;
392
393 double Jfi;
394 double Jsi;
395 double Jto;
396 double Jso;
397
398 ✗ double Vmin = FelisceParam::instance().vMin;
399 ✗ double Vmax = FelisceParam::instance().vMax;
400
401 ✗ felInt iType = 0;
402
403 felInt pos;
404 ✗ for (felInt i = 0; i < m_size; i++) {
405 ✗ ISLocalToGlobalMappingApply(m_localDofToGlobalDof,1,&i,&pos);
406 ✗ m_uExtrap.getValues(1,&pos,&value_uExtrap);//value_uExtrap = m_uExtrap(i)
407 ✗ value_uExtrap = (value_uExtrap-Vmin)/(Vmax-Vmin);
408
409 ✗ for(int j=0; j<4; j++) {
410 ✗ m_vecSolEDO[j].getValues(1,&pos,&valuem_solEDO[j]);
411 }
412
413 //Calculation of minf,dinf,fpinf and kinf-----fpinf=f'm_inf
414 ✗ if (FelisceParam::instance().CellsType == "hetero")
415 ✗ iType = m_cellType[pos];
416
417 ✗ if (iType < 4) {
418 double hinf;
419 double finf;
420 double rinf;
421 double tauh;
422 double tauf;
423 double taur;
424
425 ✗ hgate(hinf, tauh, value_uExtrap, iType);
426 ✗ sgate(finf, tauf, value_uExtrap, iType);
427 ✗ rgate(rinf, taur, value_uExtrap, iType);
428
429 double H;
430 double minf;
431 double dinf;
432 double kinf;
433 double fpinf;
434
435 ✗ H = Heaviside(value_uExtrap,Vc[iType]);
436
437 ✗ minf = (value_uExtrap-Vc[iType])*H;
438 ✗ dinf = /*0.5**/H*(1.+std::tanh(beta1[iType]*(value_uExtrap-V1[iType])));
439 ✗ kinf = (value_uExtrap/Vc[iType])+H*(1.-(value_uExtrap/Vc[iType]));
440 ✗ fpinf = /*0.5**/(1.-std::tanh(beta2[iType]*(value_uExtrap-V2[iType])));
441
442 ✗ Jfi = -gfi[iType]*valuem_solEDO[0]*minf*(Vfi[iType]-value_uExtrap);
443 ✗ Jsi = -gsi[iType]*dinf*valuem_solEDO[1]*fpinf;
444 ✗ Jto = gto[iType]*valuem_solEDO[2]*valuem_solEDO[3]*(value_uExtrap);//U-Vto, but Vto=0
445 ✗ Jso = gso[iType]*kinf;
446 } else {
447 ✗ Jfi = 0.;
448 ✗ Jsi = 0.;
449 ✗ Jto = 0.;
450 ✗ Jso = 0.;
451 }
452
453 //Currents:
454 ✗ const double Cm = FelisceParam::instance().Cm; //about 1 uF cm-2------1 avant
455 ✗ const double dV = Vmax-Vmin;//mV----100mV
456
457 ✗ value_ion = - ( Jfi + Jsi + Jto + Jso )*Cm*dV;
458
459 ✗ m_ion.setValue(pos,value_ion, INSERT_VALUES);
460
461 }
462 ✗ m_ion.assembly();
463
464 }
465
466
467 //******************************************************************************************************************
468 // Gates + Heaviside
469 //******************************************************************************************************************
470
471 ✗ void BCLSolver::hgate(double& hinf,double& tauh, double u, int iType) {
472
473 ✗ if( iType < 4 ) {
474 ✗ const double H = BCLSolver::Heaviside(u,Vc[iType]);
475 ✗ hinf = 1.-H;
476 ✗ tauh = tauhp[iType]-H*(tauhp[iType]-tauhm[iType]);
477 } else {
478 ✗ hinf = 0.;
479 ✗ tauh = 1.e+06;
480 }
481
482 }
483
484 ✗ void BCLSolver::fgate(double& finf,double& tauf, double u, int iType) {
485
486 ✗ if( iType < 4 ) {
487 ✗ const double H = BCLSolver::Heaviside(u,Vc[iType]);
488 ✗ finf = 1.-H;
489 ✗ tauf = taufp[iType]-H*(taufp[iType]-taufm[iType]);
490 } else {
491 ✗ finf = 0.;
492 ✗ tauf = 1.e+06;
493 }
494
495 }
496
497 ✗ void BCLSolver::rgate(double& rinf,double& taur, double u, int iType) {
498 // const double H = BCLSolver::Heaviside(u,0.6);//Vr=0.6
499 ✗ const double H = BCLSolver::Heaviside(u,Vs[iType]); // V_r = V_s (typo in the paper)
500 ✗ rinf = H;
501
502
503 ✗ if( iType < 4 ) {
504 ✗ taur = taurp[iType]-H*(taurp[iType]-taurm[iType]);
505 } else {
506 ✗ taur = 0.;
507 ✗ rinf = 1.e+06;
508 }
509 }
510
511 ✗ void BCLSolver::sgate(double& sinf,double& taus, double u, int iType) {
512
513
514 ✗ if( iType < 4 ) {
515 ✗ const double H = BCLSolver::Heaviside(u,Vc[iType]);
516 ✗ sinf = 1.-H;
517 ✗ taus = tausp[iType]-H*(tausp[iType]-tausm[iType]);
518 } else {
519 ✗ sinf = 0.;
520 ✗ taus = 1.e+06;
521 }
522
523 }
524
525 ✗ double BCLSolver::Heaviside(double V,double Vsubstr) {
526 ✗ double H = 0.;
527 ✗ const double substr = V-Vsubstr;
528
529 ✗ if(substr>=0.) {
530 ✗ H = 1.;
531 }
532 ✗ return H;
533 }
534
535
536 //******************************************************************************************************************
537 //******************************************************************************************************************
538
539
540 //******************************************************************************************************************
541 // For CVODE
542 //******************************************************************************************************************
543
544 #ifdef FELISCE_WITH_SUNDIALS
545
546 int BCLSolver::M_f(realtype t, N_Vector y, N_Vector ydot, void *f_data) {
547 (void) t;
548
549 BCLUserData* data = static_cast< BCLUserData* >(f_data);
550 const int pos = data->pos;
551 const double u = data->u;
552 const int iType = data->iType;
553
554 double hinf;
555 double tauh;
556 BCLSolver::hgate(hinf,tauh,u,iType);
557
558 double finf;
559 double tauf;
560 BCLSolver::fgate(finf,tauf,u,iType);
561
562 double rinf;
563 double taur;
564 BCLSolver::rgate(rinf,taur,u,iType);
565
566 double sinf;
567 double taus;
568 BCLSolver::sgate(sinf,taus,u,iType);
569
570
571 Ith(ydot,0) = (hinf-Ith(y,0))/tauh;//hdot
572 Ith(ydot,1) = (finf-Ith(y,1))/tauf;//fdot
573 Ith(ydot,2) = (rinf-Ith(y,2))/taur;//rdot
574 Ith(ydot,3) = (sinf-Ith(y,3))/taus;//sdot
575
576 return 0;
577 }
578 #endif
579
580 //******************************************************************************************************************
581 //******************************************************************************************************************
582 }
583
584
585