added possibility to use different time integration schemes for ion...

added possibility to use different time integration schemes for ion concentration for exponential integrator

cardmech/conf/BiVentricleTest/cellmodels.conf

@@ -10,5 +10,6 @@ ElphyModelName = BeelerReuter
 #CellModelName = FitzHughNagumo
 #CellModelScheme=LinearImplicit
 CellModelScheme=ExponentialIntegrator
+IonScheme = ExplicitEuler
 ThresholdVForActivation=-40.0
 

cardmech/conf/TetraTest/cellmodels.conf

@@ -10,5 +10,6 @@ ElphyModelName = BeelerReuter
 #CellModelName = FitzHughNagumo
 #CellModelScheme=LinearImplicit
 CellModelScheme=ExponentialIntegrator
-#CellModelScheme = ExplicitEuler
+IonScheme = ExplicitEuler
+ThresholdVForActivation=-40.0
 

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;

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]);
+  }
+}
+

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{