GCC Code Coverage Report

Directory: ./
File: Solver/schafRevisedSolver.cpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 0 75 0.0%
Branches: 0 94 0.0%
Line Branch Exec Source
1 // ______ ______ _ _ _____ ______
2 // | ____| ____| | (_)/ ____| | ____|
3 // | |__ | |__ | | _| (___ ___| |__
4 // | __| | __| | | | |\___ \ / __| __|
5 // | | | |____| |____| |____) | (__| |____
6 // |_| |______|______|_|_____/ \___|______|
7 // Finite Elements for Life Sciences and Engineering
8 //
9 // License: LGL2.1 License
10 // FELiScE default license: LICENSE in root folder
11 //
12 // Main authors: E. Schenone
13 //
14
15 // System includes
16
17 // External includes
18
19 // Project includes
20 #include "Solver/schafRevisedSolver.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 // dw/dt = a(v) * (1-w) - b(v) * w
25 //
26 // a(v) = (1-f) / (tau_open + (tau_open-tau_close)*f)
27 //
28 // b(v) = f / (tau_open + (tau_open-tau_close)*f)
29 //
30 // f(v) = 1/2 * (1 + tanh(k(u-v_g))) ~ 1/2 * (1 + k(u-v_g) - 1/3 * k^3*(u-v_g)^3 )
31 //
32 //w_0 = (v_max-v_min)^{-2}
33
34 namespace felisce {
35 ✗ SchafRevisedSolver::SchafRevisedSolver(FelisceTransient::Pointer fstransient):
36 SchafSolver(fstransient),
37 ✗ m_fTanh(nullptr),
38 ✗ m_fA(nullptr),
39 ✗ m_fB(nullptr)
40 {}
41
42
43 ✗ SchafRevisedSolver::~SchafRevisedSolver() {
44 ✗ delete [] m_fTanh;
45 ✗ delete [] m_fA;
46 ✗ delete [] m_fB;
47 }
48
49 ✗ void SchafRevisedSolver::computeTanhFunction(double* u) {
50 ✗ double& vGate = FelisceParam::instance().vGate;
51 ✗ double& k = FelisceParam::instance().kTanhMSR;
52
53 ✗ if (m_fTanh == nullptr) {
54 ✗ m_fTanh = new double[m_size];
55 }
56 ✗ for (felInt i=0; i<m_size; i++) {
57 ✗ m_fTanh[i] = 0.5 * (1 + k*(u[i]-vGate) - 1/3 * k*k*k * (u[i]-vGate)*(u[i]-vGate)*(u[i]-vGate));
58 }
59 }
60
61 ✗ void SchafRevisedSolver::computeAFunction() {
62 ✗ double& tauOpen = FelisceParam::instance().tauOpen;
63 ✗ double& tauClose = FelisceParam::instance().tauClose;
64
65 ✗ if (m_fA == nullptr) {
66 ✗ m_fA = new double[m_size];
67 }
68
69 ✗ for (felInt i=0; i<m_size; i++) {
70 ✗ if (FelisceParam::instance().hasHeteroTauClose) {
71 ✗ m_fA[i] = (1-m_fTanh[i])/(tauOpen+(m_tauClose[i]-tauOpen)*m_fTanh[i]);
72
73 } else {
74 ✗ m_fA[i] = (1-m_fTanh[i])/(tauOpen+(tauClose-tauOpen)*m_fTanh[i]);
75 }
76 }
77 }
78
79 ✗ void SchafRevisedSolver::computeBFunction() {
80 ✗ double& tauOpen = FelisceParam::instance().tauOpen;
81 ✗ double& tauClose = FelisceParam::instance().tauClose;
82
83 ✗ if (m_fB == nullptr) {
84 ✗ m_fB = new double[m_size];
85 }
86
87 ✗ for (felInt i=0; i<m_size; i++) {
88 ✗ if (FelisceParam::instance().hasHeteroTauClose) {
89 ✗ m_fB[i] = m_fTanh[i]/(tauOpen+(m_tauClose[i]-tauOpen)*m_fTanh[i]);
90 } else {
91 ✗ m_fB[i] = m_fTanh[i]/(tauOpen+(tauClose-tauOpen)*m_fTanh[i]);
92 }
93 }
94 }
95
96 ✗ void SchafRevisedSolver::computeRHS() {
97 ✗ double& dt = m_fstransient->timeStep;
98
99 ✗ m_bdf.computeRHSTime(dt, m_RHS);
100
101 felInt pos;
102 ✗ double* value_uExtrap = new double[m_size];
103 ✗ for (felInt i = 0; i < m_size; i++) {
104 ✗ ISLocalToGlobalMappingApply(m_localDofToGlobalDof,1,&i,&pos);
105 ✗ m_uExtrap.getValues(1,&pos,&value_uExtrap[i]);//value_uExtrap = m_uExtrap(i)
106 }
107 ✗ computeTanhFunction(value_uExtrap);
108 ✗ computeAFunction();
109 ✗ computeBFunction();
110
111 ✗ for (felInt i = 0; i < m_size; i++) {
112 ✗ ISLocalToGlobalMappingApply(m_localDofToGlobalDof,1,&i,&pos);
113 ✗ m_RHS.setValue(pos,m_fA[i], ADD_VALUES);
114 }
115 ✗ m_RHS.assembly();
116 }
117
118 ✗ void SchafRevisedSolver::solveEDO() {
119 ✗ double& coeffDeriv = m_bdf.coeffDeriv0();
120 ✗ double& dt = m_fstransient->timeStep;
121
122 felInt pos;
123 double value_RHS;
124 double valuem_solEDO;
125
126 ✗ for (felInt i = 0; i < m_size; i++) {
127 ✗ ISLocalToGlobalMappingApply(m_localDofToGlobalDof,1,&i,&pos);
128 ✗ m_RHS.getValues(1,&pos,&value_RHS);//value_RHS = m_RHS(i)
129 ✗ valuem_solEDO = value_RHS * 1./(coeffDeriv/dt + m_fA[i] + m_fB[i]);
130 ✗ m_solEDO.setValue(pos,valuem_solEDO, INSERT_VALUES);
131 }
132 ✗ m_solEDO.assembly();
133 }
134
135 ✗ void SchafRevisedSolver::computeIon() {
136 ✗ double& tauOut = FelisceParam::instance().tauOut;
137 ✗ double& tauIn = FelisceParam::instance().tauIn;
138 ✗ double& vMin = FelisceParam::instance().vMin;
139 ✗ double& vMax = FelisceParam::instance().vMax;
140
141 felInt pos;
142 double value_uExtrap;
143 double value_m;
144 double valuem_solEDO;
145 double value_ion;
146
147 ✗ for (felInt i = 0; i < m_size; i++) {
148 ✗ ISLocalToGlobalMappingApply(m_localDofToGlobalDof,1,&i,&pos);
149 ✗ m_solEDO.getValues(1,&pos,&valuem_solEDO);//valuem_solEDO = m_solEDO(i)
150 ✗ m_uExtrap.getValues(1,&pos,&value_uExtrap);//value_uExtrap = m_uExtrap(i)
151 ✗ value_m = value_uExtrap;
152 ✗ if (value_uExtrap < vMin) {
153 ✗ value_m = vMin;
154 ✗ } else if (value_uExtrap > vMax) {
155 ✗ value_m = vMax;
156 }
157 ✗ if (FelisceParam::instance().hasInfarct) {
158 ✗ value_ion = valuem_solEDO/tauIn*(value_m-vMin)*(value_m-vMin)*(vMax-value_uExtrap)/(vMax-vMin)-(value_uExtrap-vMin)/(m_tauOut[pos] *(vMax-vMin));
159 } else {
160 ✗ value_ion = valuem_solEDO/tauIn*(value_m-vMin)*(value_m-vMin)*(vMax-value_uExtrap)/(vMax-vMin)-(value_uExtrap-vMin)/(tauOut *(vMax-vMin));
161 }
162 ✗ m_ion.setValue(pos,value_ion, INSERT_VALUES);
163 }
164 ✗ m_ion.assembly();
165 }
166
167 }
168