Merge branch 'convergenceTestsMonodomain' into 'feature'

Convergence tests monodomain

See merge request !163

.gitlab-ci.yml

@@ -2,76 +2,63 @@ stages:
   - build
   - test
   - debug
+  - benchmark
   - deploy
 
-variables:
-  OS: "ubuntu"
-  OS_VERSION: "20.04" 
-  IMAGE_NAME_CARDMECH: "cardmech-${CI_COMMIT_SHORT_SHA}-${OS}${OS_VERSION}-${CI_PIPELINE_ID}"
-  DOCKER_TLS_CERTDIR: "/certs"
+include:
+  - project: 'mpp/mpp'
+    ref: feature
+    file: '.mpp-ci-templates.yml'
 
-services:
-   - docker:dind
+variables:
+  extends: .global-variables
 
 
-build cardmech:
-  stage: build
+build:
+  extends: .build_template
   variables:
-    GIT_SUBMODULE_STRATEGY: recursive
-  image: docker:latest
-  before_script:
-    - docker info 
-    - docker login -u $MPP_REGISTRY_USER -p $MPP_REGISTRY_PASS ${MPP_REGISTRY}
-  script:
-    - sed s/REGISTRY/${MPP_REGISTRY}\\/${MPP_REGISTRY_REPO_RELEASE}/g docker/cardmech.baseimage >
-      docker/cardmech_registry.baseimage
-    - sed s/UBUNTUVERSION/${OS_VERSION}/g docker/cardmech_registry.baseimage >
-      docker/${IMAGE_NAME_CARDMECH}.baseimage
-    - docker build
-      --build-arg UPSTREAM_COMMIT=${UPSTREAM_COMMIT}
-      --no-cache -t ${MPP_REGISTRY}/${MPP_REGISTRY_REPO_DEV}/${IMAGE_NAME_CARDMECH}
-      -f docker/${IMAGE_NAME_CARDMECH}.baseimage .
-    - docker push ${MPP_REGISTRY}/${MPP_REGISTRY_REPO_DEV}/${IMAGE_NAME_CARDMECH}
-  dependencies: [ ]
-  tags: [ docker ]
+    PROJECT_NAME: 'cardmech'
 
 
 test cardmech 1/2:
-  stage: test
+  extends: .test_template
   variables:
-    GIT_STRATEGY: none
-  image: ${MPP_REGISTRY}/${MPP_REGISTRY_REPO_DEV}/${IMAGE_NAME_CARDMECH}
+    PROJECT_NAME: 'cardmech'
   script:
-    - cd /mpp/build
     - ./cardmech/test/cellmodels/TestElphyModels
     - ./cardmech/test/elasticity/TestMaterial
-  dependencies: [ "build cardmech" ]
-  tags: [ docker ]
 
 
 test cardmech 2/2:
-  stage: test
+  extends: .test_template
   variables:
-    GIT_STRATEGY: none
-  image: ${MPP_REGISTRY}/${MPP_REGISTRY_REPO_DEV}/${IMAGE_NAME_CARDMECH}
+    PROJECT_NAME: 'cardmech'
   script:
-    - cd /mpp/build
     - python3 mppyrun.py --mpi_tests=1 --mute=0
-  dependencies: [ "build cardmech" ]
-  tags: [ docker ]
 
 
-debug cardmech:
-  stage: debug
-  when: on_failure
+#debug cardmech:
+#  stage: debug
+#  when: on_failure
+#  variables:
+#    GIT_STRATEGY: none
+#  image: ${MPP_REGISTRY}/${MPP_REGISTRY_REPO_DEV}/${IMAGE_NAME_CARDMECH}
+#  script:
+#    - cd /mpp/build
+#    - ctest
+#  dependencies: [ "build cardmech" ]
+#  tags: [ docker ]
+
+
+ventricle-test:
+  extends: .benchmark_on_horeka
   variables:
-    GIT_STRATEGY: none
-  image: ${MPP_REGISTRY}/${MPP_REGISTRY_REPO_DEV}/${IMAGE_NAME_CARDMECH}
+    CMAKE_ARGS: '-DMPP_BUILD_TYPE=MppRelease'
   script:
-    - cd /mpp/build
-    - ctest
-  dependencies: [ "build cardmech" ]
-  tags: [ docker ]
+    - salloc -p cpuonly -t 00:10:00 -n 4 mpirun Elphy-M++ BiVentricleTest/start
+
+
+
 
 deploy cardmech:
   stage: deploy

cardmech/conf/BiVentricleTest/cellmodels.conf

+# === Electrophysioligcal Models: =================================
+# = TenTusscher                                                   =
+# =                                                               =
+# === Tension Models: =============================================
+# = Land                                                          =
+# =                                                               =
+# =================================================================
+ElphyModelClass = MElphyModel
+ElphyModelName = BeelerReuter
+#CellModelName = FitzHughNagumo
+#CellModelScheme=LinearImplicit
+CellModelScheme=ExponentialIntegrator
+ThresholdVForActivation=-40.0
+

cardmech/conf/BiVentricleTest/solvers.conf

+# === Iterative Solvers: ==========================================
+# = LS, CG, CGX, CGNE, BiCGStab, BiCGStab2, GMRES, MINRES, TestDG =
+# =                                                               =
+# === Direct Solvers: =============================================
+# = SuperLU, SuperLU_local, LIB_PS, PS                            =
+# =                                                               =
+# === Preconditioners: ============================================
+# = Jacobi, PointBlockJacobi, PointBlockJacobi_dG, DampedJacobi,  =
+# = GaussSeidel, PointBlockGaussSeidel, SSOR, SGS, UMF            =
+# =================================================================
+
+
+EulerReAssemble = 0;
+
+SetSolution = 0;
+
+# === Linear Solver === #
+LinearSolver = gmres;
+LinearEpsilon = 1e-8;
+LinearReduction = 1e-8;
+LinearReduction = 1e-11;
+LinearSteps = 10000;
+BasePreconditioner = GaussSeidel;
+Preconditioner = GaussSeidel;
+Preconditioner = LIB_PS;
+
+# =================================================================
+# === Gating Solver ===============================================
+GatingSolver = gmres;
+GatingPreconditioner = GaussSeidel;
+GatingEpsilon = 1e-8;
+GatingReduction = 1e-12;
+
+# === Elphy Solver ================================================
+ElphySolver = gmres;
+ElphyPreconditioner =Jacobi# WICHTIG: Mit SuperLU funktioniert die Parallelisierung nicht richtig!;
+#ElphyPreconditioner LIB_PS;
+ElphyEpsilon = 1e-8;
+ElphyReduction = 1e-12;
+
+# === Stretch Solver ==============================================
+StretchSolver = gmres;
+StretchPreconditioner = GaussSeidel;
+StretchEpsilon = 1e-8;
+StretchReduction = 1e-12;
+
+# === Mechanics Solver ============================================
+MechSolver = gmres
+MechPreconditioner = GaussSeidel;
+MechEpsilon = 1e-8;
+MechReduction = 1e-12;
+MechSteps = 10000;
+
+# === Newton Method === #
+NewtonEpsilon = 1e-6;
+NewtonReduction = 1e-12;
+NewtonSteps = 30;
+NewtonLineSearchSteps = 5;
+NewtonLinearizationReduction = 1e-12;
+
+
+
+Smoother = GaussSeidel;
+#Smoother = SuperLU;
+#SmootherDamp = 0.8;
+
+presmoothing = 3;
+postsmoothing = 3;
+
+Overlap = NoOverlap
+Overlap_Distribution = 0
+
+DebugLevel = 0;
+TimeLevel = 0;
+#gnuplot  = 1;
+precision = 10;
\ No newline at end of file

cardmech/conf/BiVentricleTest/start.conf

+# = Removes all files in the data folder upon start ======
+ClearData = True
+# ========================================================
+# = Location of .geo files ===============================
+GeoPath = ../geo/
+# ========================================================
+# = Parallelization for cardiac meshes need (plevel = 0) =
+plevel = 0
+# ========================================================
+
+# ======================================
+# = Choose Configuration for Problems  =
+# ======================================
+# === Electrophysiology ================
+loadconf = BiVentricleTest/test.conf
+# ======================================
+# = Default Solver Configuration Files =
+loadconf = BiVentricleTest/solvers.conf;
+loadconf = BiVentricleTest/cellmodels.conf
+loadconf = BiVentricleTest/verbose.conf
+# ======================================

cardmech/conf/BiVentricleTest/test.conf

+ClearData = 1 # Removes existing plots from data folder.
+Distribution=RCB
+Model = Monodomain;
+ElphyModel = SemiImplicit
+ElphyProblem =BiVentricleProblem
+# Spielt für diesen Test keine Rolle, nur beim Splitting (Iext splitted in PDE =1, Iext in ODE =0)
+IextInPDE=1
+ReassembleRHSonce=1
+OrderStaggeredScheme=wcV
+
+ExternalCurrentSmoothingInTime = Arctan;
+#Glättung Anregung im Raum: discrete, linear
+ExternalCurrentSmoothingInSpace =linear;
+
+
+######################
+#       Meshes       #
+######################
+
+ElphyPolynomialDegree = 1
+
+StartTime = 0.0
+EndTime=0.1
+
+DeltaTime = 0.0001
+DeltaTimeMin = 0.00001
+DeltaTimeMax = 0.001
+
+######################
+#  Elphy Parameters  #
+#####################
+ElphyLevel = 0
+ElphySubsteps = 1
+
+ElphyMeshPart = ventriclesonly;
+
+# Monodomain Parameters
+SurfaceToVolumeRatio= 140;.
+MembraneCapacitance=0.01;
+
+IntraLongitudinal = 0.170;
+ExtraLongitudinal = 0.620;
+IntraTransversal = 0.019;
+ExtraTransversal = 0.240;
+
+CrankNicolsonTheta = 0.5
+
+# Gerade nicht implementiert:
+#can be set to ExplicitEuler, AdamsBashforth (second order)
+#NumericalMethodForElphy=AdamsBashforth;
+ElphyIntegrator=ExplicitEuler;
+
+######################
+#      Plotting      #
+######################
+PlotStress = 0
+PlotVertexPotential = 1
+#PlotCellPotential = 1
+
+PrintVSteps =1
+# vtk Dateien werden erstellt
+PlotVTK  = 0
+PlottingSteps = 10000
+
+TimeSeries = uniformabaqus;
+
+
+EulerReAssemble = 0;
+
+
+
+
+

cardmech/conf/BiVentricleTest/verbose.conf

+ConfigVerbose = 1
+MeshVerbose = 2
+AssembleVerbose = 1
+
+# = Solver Verbose =
+MultigridVerbose = 0
+BaseSolverVerbose = -1
+LinearVerbose = -5
+NewtonVerbose = 1
+EulerVerbose = 0
+LinearImplicitVerbose = 5
+ElphyLinearSolverVerbose =-1
+MonodomainVerbose =1
+NonLinearVerbose = -5
+ElphySolverVerbose = 1
+CoupledSolverVerbose = 2
+CellModelVerbose = 10
+
+# = Main Verbose =
+ElphyVerbose = 1
+MechVerbose = 1
+
+# = Plot Verbose =
+ElphySolverVTK = 1
+CoupledSolverVTK = 1
+
+PrestressPlotVerbose = 0
+PressurePlotVerbose = 1
+DynamicPlotVerbose = 1
+PlotVTK = 1

cardmech/conf/TetraTest/solvers.conf

@@ -52,8 +52,8 @@ MechReduction = 1e-12;
 MechSteps = 10000;
 
 # === Newton Method === #
-NewtonEpsilon = 1e-3;
-NewtonReduction = 1e-10;
+NewtonEpsilon = 1e-6;
+NewtonReduction = 1e-12;
 NewtonSteps = 30;
 NewtonLineSearchSteps = 5;
 NewtonLinearizationReduction = 1e-12;

cardmech/conf/cellmodels.conf

@@ -10,9 +10,11 @@ ElphyModelClass = MElphyModel
 ElphyModelName = BeelerReuter
 #CellModelName = FitzHughNagumo
 #CellModelScheme=LinearImplicit
-CellModelScheme=RungeKutta2
+#CellModelScheme=ExplicitEuler
+#CellModelScheme=RungeKutta2
 CellModelScheme = ExponentialIntegrator
 
+ThresholdVForActivation=-40.0
 TensionModel = IBT;
 
 ElphyModel = TenTusscher2;

cardmech/conf/elphy.conf

@@ -2,26 +2,29 @@ ClearData = 1 # Removes existing plots from data folder.
 #Distribution =Stripes
 Distribution = RCB
 
+PrintMeshOnly =0;
 #Model: OneCell
 Model = Monodomain;
 
 #ElphyModel = Diffusion
 #ElphyModel = Splitting
 #ElphyModel = ImplictEuler
-#ElphyModel = SemiImplicit
+ElphyModel = SemiImplicit
 #ElphyModel = SemiImplicitNodes
 ElphyModel = LinearImplicit
 ElphyModel = LinearImplicitNodes
 #ElphyModel = LinearImplicitQuadrature
-OrderStaggeredScheme=wVc
+OrderStaggeredScheme=wcV
 ReassembleRHSonce=1
 
-ElphySplittingMethod = Strang
+ElphySplittingMethod =Godunov  #Strang, Stein,Godunov
 
 #ElphyProblem: VariableMeshProblem,TetraTestProblem,HexaTestProblem, HexaFineTestProblem, HexaLinExTestProblem, OneCellProblem,Hexa2dTestProblem,WithoutDiffusion,EigenvalueProblem, EllipsoidProblem,FullHeartProblem
-#ElphyProblem=FullHeartProblem
+#ElphyProblem =FullHeartProblem
+
 #ElphyProblem=TetraTestProblem
-#ElphyProblem=EllipsoidProblem
+ElphyProblem=BiVentricleProblem
+ElphyProblem=EllipsoidProblem
 #ElphyProblem =HexaTestProblem
 
 #ElphyProblem=Hexa2dTestProblem
@@ -40,6 +43,10 @@ ExternalCurrentSmoothingInSpace =linear;
 ######################
 #       Meshes       #
 ######################
+Mesh=TestEllipsoid3dExcitation
+Mesh=EllipsoidExcitationApex
+Mesh=ChamberDebug2DQuad
+Mesh=ChamberCube0252DQuad
 Mesh = BenchmarkBeam3DTet
 Mesh = RossiBlockSmall3DTet
 #Mesh=KovachevaHeart
@@ -48,24 +55,27 @@ Mesh=SmallTestCube053DQuad
 Mesh=TestChamberCube05WithActivation2DQuad
 #Mesh=ChamberCube052DQuad
 
-ElphyLevel = 0
 ElphyPolynomialDegree = 1
+#HexaQuadExactupTo =6
 
 StartTime = 0.0
+EndTime=0.0001
 EndTime = 0.4
 
 #EndTime = 0.02
 #DeltaTime=0.000003125
 DeltaTime = 0.0001
+
 DeltaTimeMin = 0.00001
 DeltaTimeMax = 0.001
 
 ######################
 #  Elphy Parameters  #
 #####################
+ElphyLevel = 0
 ElphySubsteps = 1
 
-#ElphyMeshPart = ventriclesonly;
+ElphyMeshPart = ventriclesonly;
 #ElphyMeshPart = atriaonly;
 ElphyMeshPart = fourchamber;
 
@@ -98,7 +108,7 @@ PlotVertexPotential = 1
 PrintVSteps =1
 # vtk Dateien werden erstellt
 PlotVTK  = 1
-PlottingSteps = 50
+PlottingSteps = 10
 
 TimeSeries = uniformabaqus;
 

cardmech/conf/m++.conf

@@ -11,13 +11,13 @@ GeoPath = ../geo/
 # = Choose Configuration for Problems  =
 # ======================================
 # === Electrophysiology ================
-#loadconf = elphy.conf
+loadconf = elphy.conf
 #loadconf = monodomain.conf;
 # === Cardiac Mechanics ================
 #loadconf = coupled.conf
 #loadconf = klotz.conf
 #loadconf = elasticity_benchmarks/test.conf
-loadconf = elasticity_benchmarks/linear.conf
+#loadconf = elasticity_benchmarks/linear.conf
 #loadconf = elasticity_benchmarks/beam.conf
 #loadconf = elasticity_benchmarks/ellipsoid.conf
 #loadconf = elasticity.conf

cardmech/conf/solvers.conf

@@ -34,7 +34,7 @@ GatingReduction = 1e-12;
 # === Elphy Solver ================================================
 ElphySolver = gmres;
 #ElphyPreconditioner =GaussSeidel# WICHTIG: Mit SuperLU funktioniert die Parallelisierung nicht richtig!;
-ElphyPreconditioner LIB_PS;
+ElphyPreconditioner=Jacobi;
 ElphyEpsilon = 1e-8;
 ElphyReduction = 1e-12;
 
@@ -52,7 +52,7 @@ MechReduction = 1e-12
 MechSteps = 20000
 
 # === Newton Method === #
-NewtonEpsilon = 1e-3
+NewtonEpsilon = 1e-6
 NewtonReduction = 1e-12;
 NewtonSteps = 300
 NewtonLineSearchSteps = 10;

cardmech/conf/verbose.conf

 ConfigVerbose = 1
 MeshVerbose = 2
-AssembleVerbose = 2
+AssembleVerbose = 1
 
 # = Solver Verbose =
 MultigridVerbose = 0
@@ -8,28 +8,23 @@ BaseSolverVerbose = -1
 LinearVerbose = -5
 NewtonVerbose = 1
 EulerVerbose = 0
-LinearImplicitVerbose = -5
-MonodomainVerbose =-1
-NonLinearVerbose = -5
-ElphySolverVerbose = -1
-MechVerbose = -1
+LinearImplicitVerbose = 5
 ElphyLinearSolverVerbose =-1
-PressureSolverVerbose = -3
-DynamicSolverVerbose = -3
-CoupledSolverVerbose = 5
-CellModelVerbose = 5 #10
+MonodomainVerbose =1
+NonLinearVerbose = -5
+ElphySolverVerbose = 1
+CoupledSolverVerbose = 2
+CellModelVerbose = 10
 
 # = Main Verbose =
 ElphyVerbose = 1
-MainMechVerbose = 2
+MechVerbose = 1
 
 # = Plot Verbose =
-ElphySolverVTK = -1
-PressureSolverVTK = 1
-DynamicSolverVTK = 1
+ElphySolverVTK = 1
 CoupledSolverVTK = 1
 
-PrestressPlotVerbose = -1
+PrestressPlotVerbose = 0
 PressurePlotVerbose = 1
 DynamicPlotVerbose = 1
 PlotVTK = 1

cardmech/src/CardiacProblem.hpp

@@ -86,6 +86,7 @@ enum ProblemType {
   ExactDiffusionOne, ExactDiffusionTwo,
   Ellipsoid,
   ChamberCube2dProblem,
+  BiVentricle,
   FullHeart,
   RobinEllipsoid, DynTractionTestPol, DynDirichletTestPol, DynDirichletTestExp, DynTractionTestExp,
   DynDirichletTestZero, DynNeumannTestPol, DynNeumannTestExp,
@@ -140,6 +141,7 @@ static std::map<std::string, ProblemType> problemIndex{
     {"OneCellProblem",           OneCell},
     {"EigenvalueProblem",        EigenvaluesCellmodel},
     {"EllipsoidProblem",         Ellipsoid},
+    {"BiVentricleProblem",         BiVentricle},
     {"ChamberCube2d",            ChamberCube2dProblem},
     {"FullHeartProblem",         FullHeart},
     {"DynamicDirichletTestPol",  DynDirichletTestPol},

cardmech/src/cellmodels/solvers/MCellSolver.cpp

@@ -25,6 +25,23 @@ void RungeKutte2Step(MCellModel &cellModel, double t, double dt, vector<double>
     VW[i] += dt * g[1][i];
   }
 }
+void HeunStep(MCellModel &cellModel, double t, double dt, vector<double> &VW,
+              vector<std::vector<double>> &g){
+  std::vector<double> k(VW.size());
+  std::uninitialized_copy(VW.begin(), VW.end(), k.begin());
+
+  double gSize = g[0].size();
+
+  cellModel.EvaluateRHS(t, VW, g[0]);
+  for (int i = 0; i < gSize; ++i) {
+    k[i] = VW[i] +  dt * g[0][i];
+  }
+
+  cellModel.EvaluateRHS(t +  dt, k, g[1]);
+  for (int i = 0; i < gSize; ++i) {
+    VW[i] += 0.5*dt * (g[1][i]+g[0][i]);
+  }
+}
 
 void RungeKutta4Step(MCellModel &cellModel, double t, double dt, vector<double> &VW,
                      vector<std::vector<double>> &g) {

cardmech/src/cellmodels/solvers/MCellSolver.hpp

@@ -11,6 +11,8 @@ void ExplicitEulerStep(MCellModel &cellModel, double t, double dt, std::vector<d
 
 void RungeKutte2Step(MCellModel &cellModel, double t, double dt, std::vector<double> &VW,
                      std::vector<std::vector<double>> &g);
+void HeunStep(MCellModel &cellModel, double t, double dt, std::vector<double> &VW,
+                     std::vector<std::vector<double>> &g);
 
 void RungeKutta4Step(MCellModel &cellModel, double t, double dt, std::vector<double> &VW,
                      std::vector<std::vector<double>> &g);

cardmech/src/cellmodels/solvers/MElphySolver.cpp

 #include "MElphySolver.hpp"
 #include "Logging.hpp"
 #include "Spectrum.hpp"
+#include <solvers/SplittingSolver.hpp>
 
 
 void MElphySolver::initialize(Vector &u) {
@@ -164,26 +165,34 @@ void MElphySolver::PrintMaxMinGating() const {
 }
 
 void MElphySolver::selectScheme(const string &solvingScheme) {
-  if (solvingScheme == "ExplicitEuler") {
-    steps = 1;
-    SolveStep = ExplicitEulerStep;
-  } else if (solvingScheme == "RungeKutta2") {
+  if(getSplitting()==STEIN){
     steps = 2;
-    SolveStep = RungeKutte2Step;
-  } else if (solvingScheme == "RungeKutta4") {
-    steps = 4;
-    SolveStep = RungeKutta4Step;
-  } else if (solvingScheme == "ExponentialIntegrator") {
-    steps = 1;
-    SolveStep = ExponentialIntegratorStep;
-  }else if(solvingScheme=="EIwithNewtonStep"){
-    steps=1;
-    SolveStep=EIWithNewtonStep;
-  }else if (solvingScheme == "Eigenvalue") {
-    steps = 1;
-    SolveStep = EigenvalueStep;
-  } else {
-    Exit("CellModelScheme \"" + solvingScheme + "\" unknown.")
+    SolveStep = HeunStep;
+  }else {
+    if (solvingScheme == "ExplicitEuler") {
+      steps = 1;
+      SolveStep = ExplicitEulerStep;
+    } else if (solvingScheme == "RungeKutta2") {
+      steps = 2;
+      SolveStep = RungeKutte2Step;
+    }else if (solvingScheme == "Heun") {
+      steps = 2;
+      SolveStep = HeunStep;
+    } else if (solvingScheme == "RungeKutta4") {
+      steps = 4;
+      SolveStep = RungeKutta4Step;