diff --git a/cardmech/conf/BiVentricleTest/cellmodels.conf b/cardmech/conf/BiVentricleTest/cellmodels.conf
index 2fe237f67949e300ce19d955a846008b81d90cf7..c2cbad731c2c7541e95661baea207dcae6e3a22c 100644
--- a/cardmech/conf/BiVentricleTest/cellmodels.conf
+++ b/cardmech/conf/BiVentricleTest/cellmodels.conf
@@ -10,5 +10,6 @@ ElphyModelName = BeelerReuter
#CellModelName = FitzHughNagumo
#CellModelScheme=LinearImplicit
CellModelScheme=ExponentialIntegrator
+IonScheme = ExplicitEuler
ThresholdVForActivation=-40.0
diff --git a/cardmech/conf/TetraTest/cellmodels.conf b/cardmech/conf/TetraTest/cellmodels.conf
index 352619b454cf5c0eff76064657da736e94832b8c..c2cbad731c2c7541e95661baea207dcae6e3a22c 100644
--- a/cardmech/conf/TetraTest/cellmodels.conf
+++ b/cardmech/conf/TetraTest/cellmodels.conf
@@ -10,5 +10,6 @@ ElphyModelName = BeelerReuter
#CellModelName = FitzHughNagumo
#CellModelScheme=LinearImplicit
CellModelScheme=ExponentialIntegrator
-#CellModelScheme = ExplicitEuler
+IonScheme = ExplicitEuler
+ThresholdVForActivation=-40.0
diff --git a/cardmech/conf/cellmodels.conf b/cardmech/conf/cellmodels.conf
index 5c5f8b028bb278d563ac8d948a70222aea0f856e..2a6b46b0180150efc5d32607fa8ea355050a92c3 100644
--- a/cardmech/conf/cellmodels.conf
+++ b/cardmech/conf/cellmodels.conf
@@ -13,6 +13,7 @@ ElphyModelName = BeelerReuter
#CellModelScheme=ExplicitEuler
#CellModelScheme=RungeKutta2
CellModelScheme = ExponentialIntegrator
+IonScheme = ExplicitEuler # ExplicitEuler,RungeKutta2,RungeKutta4
ThresholdVForActivation=-40.0
TensionModel = IBT;
diff --git a/cardmech/src/electrophysiology/solvers/LinearImplicitSolver.cpp b/cardmech/src/electrophysiology/solvers/LinearImplicitSolver.cpp
index 3ccda8e5d5d8aa1b667ddf0c7b6a2d009da99371..adf6c7b51c6a4d9e8408130fe5b2ae42e3c24daa 100644
--- a/cardmech/src/electrophysiology/solvers/LinearImplicitSolver.cpp
+++ b/cardmech/src/electrophysiology/solvers/LinearImplicitSolver.cpp
@@ -109,26 +109,7 @@ void LinearImplicitSolver::SolveGating( double t, double dt) {
}
}
-void LinearImplicitSolver::SolveConcentration( double t, double dt) {
- std::vector<double> vcw(M.Size());
- std::vector<double> G(M.Size());
- int position=M.getGatingIndex();
- for (row r = (*gating)[0].rows(); r != (*gating)[0].rows_end(); ++r) {
- for (int j = 0; j < r.n(); ++j) {
- for (int i = 0; i < M.Size(); ++i) {
- vcw[i] = (*gating)[i](r, j);
- }
- M.EvaluateConcentration( t,vcw, G);
- for(int i=1;i<position;i++) {
- vcw[i] += dt * M.TimeScale() * G[i];
- (*gating)[i](r, j) = vcw[i];
- }
- if (vcw[1] < 0.0) {
- mout << "concentration of Ca is less than zero " << vcw[1] << endl;
- }
- }
- }
-}
+
void LinearImplicitSolver::SolveTension(Vectors &values,double t, double dt) {
std::vector<double> vcw(M.Size()+1);
@@ -407,10 +388,87 @@ void LinearImplicitSolver::selectOrder(std::string o) {
};
}
}
-/*void LiniearImplicitSolverQuad::SolveStep( double t, double dt){
- A.SystemMatrix((*gating)[0], *K);
- (*S)(*K);
- A.RHS((*gating), *b);
- (*gating)[0] = (*S) * (*b);
- //updateMaxMin(VCW);
-}*/
+void LinearImplicitSolver::selectIonScheme(const string &solvingIonScheme) {
+ if (solvingIonScheme == "ExplicitEuler") {
+ ionSteps = 1;
+ SolveConcentrationStep = ExplicitEulerIonStep;
+ } else if (solvingIonScheme == "RungeKutta2") {
+ ionSteps = 2;
+ SolveConcentrationStep = RungeKutta2IonStep;
+ }else if (solvingIonScheme == "RungeKutta4") {
+ ionSteps = 4;
+ SolveConcentrationStep = RungeKutta4IonStep;
+ } else {
+ Exit("CellModelScheme \"" + solvingIonScheme + "\" unknown.")
+ }
+}
+
+void LinearImplicitSolver::SolveConcentration(double t, double dt){
+ std::vector<double> vcw(M.Size());
+ std::vector<vector<double>> G(ionSteps);
+ for (int i = 0; i < ionSteps; ++i) {
+ G[i] = std::vector<double>(M.Size());
+ }
+ for (row r = (*gating)[0].rows(); r != (*gating)[0].rows_end(); ++r) {
+ for (int j = 0; j < r.n(); ++j) {
+ for (int i = 0; i < M.Size(); ++i) {
+ vcw[i] = (*gating)[i](r, j);
+ }
+ SolveConcentrationStep(M,t,M.TimeScale()*dt,vcw,G);
+ for(int i=1;i<M.getGatingIndex();i++) {
+ (*gating)[i](r, j) = vcw[i];
+ }
+
+ if (vcw[1] < 0.0) {
+ mout << "concentration of Ca is less than zero " << vcw[1] << endl;
+ }
+ }
+ }
+}
+void ExplicitEulerIonStep(MCellModel &cellModel,double t, double dt, std::vector<double> &vcw,
+ std::vector<std::vector<double>> &g){
+ cellModel.EvaluateConcentration( t,vcw, g[0]);
+ for(int i=1;i<cellModel.getGatingIndex();i++) {
+ vcw[i] += dt * g[0][i];
+ }
+}
+
+void RungeKutta2IonStep( MCellModel &cellModel,double t, double dt, std::vector<double> &vcw,
+ std::vector<std::vector<double>> &g) {
+ std::vector<double> k(vcw.size());
+ std::uninitialized_copy(vcw.begin(), vcw.end(), k.begin());
+ int position=cellModel.getGatingIndex();
+ cellModel.EvaluateConcentration(t, vcw, g[0]);
+ for (int i = 1; i < position; ++i) {
+ k[i] = vcw[i] + 0.5 * dt * g[0][i];
+ }
+ cellModel.EvaluateConcentration(t + 0.5 * dt, k, g[1]);
+ for (int i = 1; i < position; ++i) {
+ vcw[i] += dt * g[1][i];
+ }
+}
+
+void RungeKutta4IonStep(MCellModel &cellModel,double t, double dt, std::vector<double> &vcw,
+ std::vector<std::vector<double>> &g) {
+ std::vector<double> k(vcw.size());
+ std::uninitialized_copy(vcw.begin(), vcw.end(), k.begin());
+ int position=cellModel.getGatingIndex();
+
+ cellModel.EvaluateConcentration(t, k, g[0]);
+ for (int i = 1; i < position; ++i) {
+ k[i] = vcw[i] + 0.5 * dt * g[0][i];
+ }
+ cellModel.EvaluateConcentration(t + 0.5 * dt, k, g[1]);
+ for (int i = 1; i < position; ++i) {
+ k[i] = vcw[i] + 0.5 * dt * g[1][i];
+ }
+ cellModel.EvaluateConcentration(t + 0.5 * dt, k, g[2]);
+ for (int i = 1; i < position; ++i) {
+ k[i] = vcw[i] + dt * g[2][i];
+ }
+ cellModel.EvaluateConcentration(t + dt, k, g[3]);
+ for (int i = 1; i < position; ++i) {
+ vcw[i] += (dt / 6.0) * (g[0][i] + 2.0 * g[1][i] + 2.0 * g[2][i] + g[3][i]);
+ }
+}
+
diff --git a/cardmech/src/electrophysiology/solvers/LinearImplicitSolver.hpp b/cardmech/src/electrophysiology/solvers/LinearImplicitSolver.hpp
index a23ebb5360c414fb83eef1652498777a168de947..35f9558464d745a050530452fc14561bc19a0d08 100644
--- a/cardmech/src/electrophysiology/solvers/LinearImplicitSolver.hpp
+++ b/cardmech/src/electrophysiology/solvers/LinearImplicitSolver.hpp
@@ -8,6 +8,15 @@
#include <MCellModel.hpp>
#include "Newton.hpp"
+
+void ExplicitEulerIonStep( MCellModel &cellModel,double t, double dt, std::vector<double> &VW,
+ std::vector<std::vector<double>> &g);
+
+void RungeKutta2IonStep( MCellModel &cellModel,double t, double dt, std::vector<double> &VW,
+ std::vector<std::vector<double>> &g);
+
+void RungeKutta4IonStep( MCellModel &cellModel,double t, double dt, std::vector<double> &VW,
+ std::vector<std::vector<double>> &g);
class LinearImplicitSolver : public ElphySolver {
int ReferenceCalculation;
int max_estep;
@@ -33,6 +42,10 @@ protected:
bool iextInPDE{false};
bool plotVTK{false};
int sizeElphyModel{0};
+ int ionSteps{1};
+ void selectIonScheme(const string &solvingIonScheme);
+ std::function<void( MCellModel &cellModel,double t, double dt, std::vector<double> &VW,
+ std::vector<std::vector<double>> &g)> SolveConcentrationStep;
public:
explicit LinearImplicitSolver(IElphyAssemble &assemble, MCellModel &cellModel) :
ElphySolver(assemble),
@@ -47,6 +60,9 @@ public:
if (cellModelType==TENSION){
sizeElphyModel=M.ElphySize();
}
+ string solvingIonScheme{"ExplicitEuler"};
+ config.get("IonScheme", solvingIonScheme);
+ selectIonScheme(solvingIonScheme);
config.get("PlotVTK",plotVTK);
config.get("LinearImplicitVerbose", verbose);
@@ -69,9 +85,10 @@ public:
virtual void Solve(Vector &V, TimeSeries &timeSeries);
+ void SolveConcentration(double t, double dt);
+
void SolveGating( double t, double dt);
- void SolveConcentration( double t, double dt);
void SolveTension(Vectors &values,double t,double dt);
virtual void SolvePDE();
@@ -89,6 +106,8 @@ public:
void updateMaxMin(Vectors &VCW);
void PrintMaxMinGating();
+
+
};
/*class LiniearImplicitSolverQuad: public LinearImplicitSolver{