GCC Code Coverage Report

Directory: ./
File: Solver/schafSolver.cpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 0 68 0.0%
Branches: 0 106 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: A. Collin
13 //
14
15 // System includes
16
17 // External includes
18
19 // Project includes
20 #include "Solver/schafSolver.hpp"
21
22 // I_{ion} = reac_amp* (w/tau_in) (u-v_min) * (u-v_min) *(v_max - u)/(v_min - v_max) - (u- v_min) /(tau_out * (v_min- v_max))
23 //
24 // ((v_max-v_min)^{-2} - w )/tau_open if u < vcrit
25 // dw/dt ={
26 // -w/tau_close if u > vcrit
27 //
28 //w_0 = (v_max-v_min)^{-2}
29
30 namespace felisce {
31 ✗ SchafSolver::SchafSolver(FelisceTransient::Pointer fstransient):
32 ✗ IonicSolver(fstransient)
33 {}
34
35
36 ✗ SchafSolver::~SchafSolver() {
37 ✗ m_tauClose.clear();
38 ✗ m_tauOut.clear();
39 ✗ m_tauIn.clear();
40 }
41
42
43 ✗ void SchafSolver::computeRHS() {
44 ✗ double& vGate = FelisceParam::instance().vGate;
45 ✗ double& vMin = FelisceParam::instance().vMin;
46 ✗ double& vMax = FelisceParam::instance().vMax;
47 ✗ double& dt = m_fstransient->timeStep;
48 ✗ double& tauOpen = FelisceParam::instance().tauOpen;
49
50 ✗ m_bdf.computeRHSTime(dt, m_RHS);
51
52 felInt pos;
53 double value_uExtrap;
54 double value_RHS;
55
56 ✗ for (felInt i = 0; i < m_size; i++) {
57 ✗ ISLocalToGlobalMappingApply(m_localDofToGlobalDof,1,&i,&pos);
58 ✗ m_uExtrap.getValues(1,&pos,&value_uExtrap);//value_uExtrap = m_uExtrap(i)
59 ✗ if( value_uExtrap < vGate) {
60 ✗ value_RHS = (1./((vMax - vMin)*(vMax - vMin) * tauOpen));
61 } else
62 ✗ value_RHS = 0.;
63 ✗ m_RHS.setValue(pos,value_RHS, ADD_VALUES);
64 }
65 ✗ m_RHS.assembly();
66
67 }
68
69 ✗ void SchafSolver::solveEDO() {
70 ✗ double& tauOpen = FelisceParam::instance().tauOpen;
71 ✗ double& vGate = FelisceParam::instance().vGate;
72 ✗ double& tauClose = FelisceParam::instance().tauClose;
73 ✗ double& coeffDeriv = m_bdf.coeffDeriv0();
74 ✗ double& dt = m_fstransient->timeStep;
75
76 felInt pos;
77 double value_uExtrap;
78 double value_RHS;
79 double valuem_solEDO;
80
81 ✗ for (felInt i = 0; i < m_size; i++) {
82 ✗ ISLocalToGlobalMappingApply(m_localDofToGlobalDof,1,&i,&pos);
83 ✗ m_uExtrap.getValues(1,&pos,&value_uExtrap);//value_uExtrap = m_uExtrap(i)
84 ✗ m_RHS.getValues(1,&pos,&value_RHS);//value_RHS = m_RHS(i)
85 ✗ if( value_uExtrap < vGate) {
86 ✗ valuem_solEDO = value_RHS * 1./(coeffDeriv/dt + 1./tauOpen);
87 } else {
88 ✗ if (FelisceParam::instance().hasHeteroTauClose) {
89 ✗ valuem_solEDO = value_RHS * 1./(coeffDeriv/dt + 1./m_tauClose[pos]);
90 } else {
91 ✗ valuem_solEDO = value_RHS * 1./(coeffDeriv/dt + 1./tauClose);
92 }
93 }
94 ✗ m_solEDO.setValue(pos,valuem_solEDO, INSERT_VALUES);
95 }
96 ✗ m_solEDO.assembly();
97 }
98
99 ✗ void SchafSolver::computeIon() {
100 ✗ if (FelisceParam::instance().stateFilter) {
101 ✗ if(!m_stabInit) {
102 ✗ m_stabTerm.duplicateFrom(m_ion);
103 ✗ m_stabInit=true;
104 }
105 ✗ m_stabTerm.zeroEntries();
106 }
107
108 ✗ double& tauOut = FelisceParam::instance().tauOut;
109 ✗ double& tauIn = FelisceParam::instance().tauIn;
110 ✗ double& vMin = FelisceParam::instance().vMin;
111 ✗ double& vMax = FelisceParam::instance().vMax;
112
113 felInt pos;
114 double value_uExtrap;
115 double valuem_solEDO;
116 double value_ion;
117 ✗ felInt sizeVent = 45580;
118 ✗ if (FelisceParam::instance().typeOfAppliedCurrent == "ellipseheart")
119 ✗ sizeVent = 32320;
120
121 ✗ for (felInt i = 0; i < m_size; i++) {
122 ✗ ISLocalToGlobalMappingApply(m_localDofToGlobalDof,1,&i,&pos);
123 ✗ felInt meshId = pos;
124 ✗ AOPetscToApplication(m_ao,1,&meshId);
125
126 ✗ m_solEDO.getValues(1,&pos,&valuem_solEDO);//valuem_solEDO = m_solEDO(i)
127 ✗ m_uExtrap.getValues(1,&pos,&value_uExtrap);//value_uExtrap = m_uExtrap(i)
128 ✗ if (FelisceParam::instance().hasInfarct) {
129 ✗ value_ion = valuem_solEDO/tauIn*(value_uExtrap-vMin)*(value_uExtrap-vMin)*(vMax-value_uExtrap)/(vMax-vMin)-(value_uExtrap-vMin)/(m_tauOut[pos] *(vMax-vMin));
130 ✗ if (FelisceParam::instance().typeOfIonicModel == "courtAtriaSchafVent") {
131 ✗ if (meshId > (sizeVent-1) )
132 ✗ value_ion = 0.0;
133 }
134 }
135 else {
136 ✗ value_ion = valuem_solEDO/tauIn*(value_uExtrap-vMin)*(value_uExtrap-vMin)*(vMax-value_uExtrap)/(vMax-vMin)-(value_uExtrap-vMin)/(tauOut *(vMax-vMin));
137 ✗ if (FelisceParam::instance().typeOfIonicModel == "courtAtriaSchafVent") {
138 ✗ if (meshId > (sizeVent-1) )
139 ✗ value_ion = 0.0;
140 }
141 }
142 ✗ m_ion.setValue(pos,value_ion, INSERT_VALUES);
143 }
144
145 ✗ m_ion.assembly();
146 }
147
148 }
149