Monodomain project

.github/workflows/ci_pipeline.yml

@@ -144,7 +144,53 @@ jobs:
           path: |
             data_libpressio
             coverage_libpressio_3.10.dat
-
+
+  user_monodomain_tests_linux:
+    runs-on: ubuntu-latest
+
+    defaults:
+      run:
+        shell: bash -l {0}
+
+    steps:
+      - name: Checkout
+        uses: actions/checkout@v3
+
+      - name: Install Conda environment with Micromamba
+        uses: mamba-org/setup-micromamba@v1
+        with:
+          environment-file: "pySDC/projects/Monodomain/etc/environment-monodomain.yml"
+          create-args: >-
+              python=3.10
+
+      - name: Compile C++ ionic models
+        env:
+          IONIC_MODELS_PATH: "pySDC/projects/Monodomain/problem_classes/ionicmodels/cpp"
+        run: |
+          c++ -O3 -Wall -shared -std=c++11 -fPIC -fvisibility=hidden $(python3 -m pybind11 --includes) ${IONIC_MODELS_PATH}/bindings_definitions.cpp -o ${IONIC_MODELS_PATH}/ionicmodels$(python3-config --extension-suffix)
+
+      - name: Run pytest for CPU stuff
+        run: |
+          echo "print('Loading sitecustomize.py...')
+          import coverage
+          coverage.process_startup() " > sitecustomize.py
+          coverage run -m pytest --continue-on-collection-errors -v --durations=0 pySDC/tests -m monodomain
+
+      - name: Make coverage report
+        run: |
+          mv data data_monodomain
+          coverage combine
+          mv .coverage coverage_monodomain_3.10.dat
+
+      - name: Uploading artifacts
+        uses: actions/upload-artifact@v3
+        with:
+          name: cpu-test-artifacts
+          path: |
+            data_monodomain
+            coverage_monodomain_3.10.dat
+
+
 #  user_cpu_tests_macos:
 #    runs-on: macos-12
 #
@@ -206,6 +252,7 @@ jobs:
       - lint
       - user_cpu_tests_linux
       - user_libpressio_tests
+      - user_monodomain_tests_linux
 #       - wait_for_gitlab
 
     defaults:

docs/source/index.rst

@@ -52,6 +52,7 @@ Projects
    projects/DAE.rst
    projects/compression.rst
    projects/second_order.rst
+   projects/monodomain.rst
 
 
 API documentation

docs/source/projects/monodomain.rst

@@ -0,0 +1 @@
+.. include:: /../../pySDC/projects/Monodomain/README.rst

pySDC/projects/Monodomain/README.rst

@@ -0,0 +1,94 @@
+Exponential SDC for the Monodomain Equation in Cardiac Electrophysiology
+==============================================================================
+This project implements the exponential spectral deferred correction (ESDC) method for the monodomain equation in cardiac electrophysiology.
+The method proposed here is an adaptation of the `ESDC method proposed by T. Buvoli  <https://doi.org/10.1137/19M1256166>`_ to the monodomain equation. 
+In particular, the implicit-explicit Rush-Larsen method is used as correction scheme. Moreover, not all model components have exponential terms, therefore the resulting method is an hybrid between ESDC and the standard SDC method.
+
+Monodomain equation
+-------------------
+The monodomain equation models the electrical activity in the heart. It is a reaction-diffusion equation coupled with an ordinary differential equation for the ionic model and is given by
+
+.. math::
+    \begin{align}
+    \chi (C_m\frac{\partial V}{\partial t}+I_{ion}(V,z_E,z_e, t)) &= \nabla \cdot (D \nabla V) & \quad \text{in } &\Omega \times (0,T), \\
+    \frac{\partial z_E}{\partial t} &= g_E(V,z_E,z_e) & \quad \text{in } &\Omega \times (0,T), \\
+    \frac{\partial z_e}{\partial t} &= \Lambda_e(V)(z_e-z_{e,\infty}(V)) & \quad \text{in } &\Omega \times (0,T), \\
+    \end{align}
+
+plus the boundary conditions, where :math:`V(t,x)\in\mathbb{R}` is the transmembrane potential and :math:`z_E(t,x)\in\mathbb{R}^n`, :math:`z_e(t,x)\in\mathbb{R}^m` are the ionic model state variables. 
+The ionic model right-hand side :math:`g_E` is a general nonlinear term, while :math:`\Lambda_e` is a diagonal matrix. The typical range for the number of unknowns :math:`N=1+n+m` is :math:`N\in [4,50]` and depends on the ionic model of choice. 
+
+Spatial discretization yields a system of ODEs which can be written in compact form as
+
+.. math::
+    \mathbf y'=f_I(\mathbf y)+f_E(\mathbf y)+f_e(\mathbf y),
+
+where :math:`\mathbf y(t)\in\mathbb{R}^{M N}` is the vector of unknowns and :math:`M` the number of mesh nodes. 
+Concerning the right-hand sides, :math:`f_I` is a linear term for the discrete diffusion, :math:`f_E` is a nonlinear but non-stiff term for :math:`I_{ion},g_E`, and :math:`f_e` is a severely stiff term for :math:`\Lambda_e(V)(z_e-z_{e,\infty}(V))`.
+
+The standard (serial) way of integrating the monodomain equation is by using a splitting method, where :math:`f_I` is integrated implicitly, :math:`f_E` explicitly, and :math:`f_e` using the exponential Euler method (which is inexpensive due to the diagonal structure of :math:`\Lambda_e`). We denote this method as IMEXEXP.
+
+The ESDC method for the monodomain equation
+-------------------------------------------
+A possible way to parallelize the integration of the monodomain equation is by employing the SDC method in combination with the IMEXEXP approach for the correction scheme (preconditioner).
+However, this approach is unstable due to the severe stiffness of :math:`f_e`. 
+Therefore we propose a hybrid method, where we employ SDC for the :math:`f_I,f_E` terms and ESDC for the :math:`f_e` term. For the correcttion scheme we still use the IMEXEXP method. 
+The resulting method can be seen as a particular case of ESDC and will be denoted by ESDC in the next figures, for simplicity.
+
+Running the code
+----------------
+Due to their complexity, ionic models are coded in C++ and wrapped to Python. Therefore, before running any example you need to compile the ionic models by running the following command in the root folder:
+
+.. code-block::
+
+   export IONIC_MODELS_PATH=pySDC/projects/Monodomain/problem_classes/ionicmodels/cpp
+   c++ -O3 -Wall -shared -std=c++11 -fPIC -fvisibility=hidden $(python3 -m pybind11 --includes) ${IONIC_MODELS_PATH}/bindings_definitions.cpp -o ${IONIC_MODELS_PATH}/ionicmodels$(python3-config --extension-suffix)
+
+Then an example can be run:
+
+.. code-block::
+
+   cd pySDC/projects/Monodomain/run_scripts
+   mpirun -n 4 python run_MonodomainODE_cli.py --dt 0.05 --end_time 0.2 --num_nodes 6,3 --domain_name cube_1D --refinements 0 --ionic_model_name TTP --truly_time_parallel --n_time_ranks 4
+
+Stability
+---------
+We display here the stability domain of the ESDC and SDC methods, both with IMEXEXP as correction scheme, applied to the test problem 
+
+.. math:: 
+    y'=\lambda_I y+\lambda_E y+\lambda_e y, 
+
+with :math:`\lambda_I,\lambda_E,\lambda_e` representing :math:`f_I,f_E,f_e`, respectively.
+We fix :math:`\lambda_E=-1` and vary the stiff terms :math:`\lambda_I,\lambda_e` only. We see that the ESDC method is stable for all tested values of :math:`\lambda_I,\lambda_e`, while SDC is not.
+
+.. image:: ../../../data/stability_domain_IMEXEXP_EXPRK.png
+   :scale: 60 %
+.. image:: ../../../data/stability_domain_IMEXEXP.png
+   :scale: 60 %
+
+Convergence
+-----------
+Here we verify convergence of the ESDC method for the monodomain equation. 
+We fix the number of collocation nodes to :math:`m=6` and perform a convergence experiment fixing the number of sweeps to either :math:`k=3` or :math:`k=6`.
+We use the ten Tusscher-Panfilov ionic model, which is employed in practical applications.
+We see that we gain one order of accuracy per sweep, as expected. 
+
+.. image:: ../../../data/convergence_ESDC_fixed_iter.png
+   :scale: 100 %
+
+
+Iterations
+----------
+Here we consider three methods:
+
+* ESDC: with :math:`m=6` collocation nodes.
+* MLESDC: This is a multilevel version of ESDC with :math:`m=6` collocation nodes on the fine level and :math:`m=3` nodes on the coarse level.
+* PFASST: Combination of the PFASST parallelization method with MLESDC, using 24 processors.
+
+We display the number of iterations required by each method to reach a given tolerance and the residual at convergence. As ionic model we use again the ten Tusscher-Panfilov model.
+We see that PFASST requires a reasonalbly small number of iterations, comparable to the serial counterparts ESDC and MLESDC.
+
+.. image:: ../../../data/niter_VS_time.png
+   :scale: 100 %
+.. image:: ../../../data/res_VS_time.png
+   :scale: 100 %

pySDC/projects/Monodomain/datatype_classes/my_mesh.py

@@ -0,0 +1,5 @@
+from pySDC.implementations.datatype_classes.mesh import MultiComponentMesh
+
+
+class imexexp_mesh(MultiComponentMesh):
+    components = ['impl', 'expl', 'exp']

pySDC/projects/Monodomain/etc/environment-monodomain.yml

@@ -0,0 +1,24 @@
+name: pySDC_monodomain
+channels:
+  - conda-forge
+  - defaults
+dependencies:
+  - python
+  - numpy
+  - scipy>=0.17.1
+  - matplotlib>=3.0
+  - sympy>=1.0
+  - numba>=0.35
+  - dill>=0.2.6
+  - pytest
+  - pytest-benchmark
+  - pytest-timeout
+  - pytest-order
+  - coverage[toml]
+  - sphinx
+  - numdifftools
+  - pybind11
+  - mpi4py
+  - mpich
+  - tqdm
+  - pymp-pypi

pySDC/projects/Monodomain/hooks/HookClass_pde.py

@@ -0,0 +1,34 @@
+from pySDC.core.Hooks import hooks
+
+
+class pde_hook(hooks):
+    """
+    Hook class to write the solution to file.
+    """
+
+    def __init__(self):
+        super(pde_hook, self).__init__()
+
+    def pre_run(self, step, level_number):
+        """
+        Overwrite default routine called before time-loop starts
+        It calls the default routine and then writes the initial value to file.
+        """
+        super(pde_hook, self).pre_run(step, level_number)
+
+        L = step.levels[level_number]
+        P = L.prob
+        if level_number == 0 and L.time == P.t0:
+            P.write_solution(L.u[0], P.t0)
+
+    def post_step(self, step, level_number):
+        """
+        Overwrite default routine called after each step.
+        It calls the default routine and then writes the solution to file.
+        """
+        super(pde_hook, self).post_step(step, level_number)
+
+        if level_number == 0:
+            L = step.levels[level_number]
+            P = L.prob
+            P.write_solution(L.uend, L.time + L.dt)

pySDC/projects/Monodomain/hooks/HookClass_post_iter_info.py

@@ -0,0 +1,34 @@
+import time
+from pySDC.core.Hooks import hooks
+
+
+class post_iter_info_hook(hooks):
+    """
+    Hook class to write additional iteration information to the command line.
+    It is used to print the final residual, after u[0] has been updated with the new value from the previous step.
+    This residual is the one used to check the convergence of the iteration and when running in parallel is different from
+    the one printed at IT_FINE.
+    """
+
+    def __init__(self):
+        super(post_iter_info_hook, self).__init__()
+
+    def post_iteration(self, step, level_number):
+        """
+        Overwrite default routine called after each iteration.
+        It calls the default routine and then writes the residual to the command line.
+        We call this the residual at IT_END.
+        """
+        super().post_iteration(step, level_number)
+        self.__t1_iteration = time.perf_counter()
+
+        L = step.levels[level_number]
+
+        self.logger.info(
+            "Process %2i on time %8.6f at stage %15s: ----------- Iteration: %2i --------------- " "residual: %12.8e",
+            step.status.slot,
+            L.time,
+            "IT_END",
+            step.status.iter,
+            L.status.residual,
+        )