Merge pull request #34 from finitewave/develop

Develop

Author: TiNezlobinsky

CHANGELOG.md

@@ -1,5 +1,32 @@
 # Changelog
 
+## [0.9.3] May 2026
+
+### Added
+
+- Initial conditions support. Model variables can now be initialized either as scalars (uniform across the tissue)
+  or as arrays (node-wise values). All initial condition fields use the `init_*` prefix.
+  See the **change_initial_conditions.py** example.
+
+- Courant–Friedrichs–Lewy (CFL) condition check to warn about potential instability of the user-defined
+  numerical parameters.
+
+- Propagation tests for validating model conduction velocity.
+
+- **Reentry and Spiral waves** tutorial.
+
+### Fixes
+
+- Adjusted diffusion coefficients (`D_model`) in Fenton–Karma and Mitchell–Schaeffer models
+  to ensure physiological conduction velocities.
+
+- Removed the unnecessary variable `irel` in the Courtemanche model (now treated as a parameter).
+
+- Standardized variable and parameter names in the ten-Tusscher–Panfilov 2006 model.
+
+- Fixed an API issue and a potential deadlock in AnimationBuilder2D/3D.
+
+
 ## [0.9.0] March 2026
 
 ### Added

README.md

@@ -147,6 +147,8 @@ plt.tight_layout()
 plt.show()
 ```
 
+![Alt text](images/quick_start_output.png)
+
 Now, let's move on to a detailed explanation. 
 
 ## Table of Contents

Tutorials/Anisotropy2D.ipynb

@@ -0,0 +1,106 @@
+"""
+Change initial conditions in the Luo-Rudy model 
+======================================
+
+Overview:
+---------
+This example demonstrates how to set spatially varying initial conditions 
+in the Luo-Rudy 1991 cardiac model. Here we create a non-conductive obstacle in the center of the tissue 
+and initialize an anatomical reentry by keeping the right half of the tissue refractory.
+
+Simulation Setup:
+-----------------
+- Tissue Grid: A 256×256 cardiac tissue domain.
+- Non-conductive Obstacle: A circular non-conductive region is created in the center of 
+    the tissue to simulate an anatomical obstacle.
+- Stimulation: A localized stimulus is applied to the left half of the upper border.
+- Initial Conditions: The right half of the tissue is initialized in a refractory state by setting the gating 
+    variables 'h' and 'j' to 0, while the left half is initialized to 0.9 (non-refractory).
+- Time and Space Resolution:
+    - Temporal step (dt): 0.01 
+    - Spatial resolution (dr): 0.4 
+    - Total simulation time (t_max): 200
+
+Execution:
+----------
+1. Create a 2D cardiac tissue grid and define a non-conductive obstacle.
+2. Apply a stimulus to the left half of the upper border.
+3. Set up and initialize the Luo-Rudy 1991 model with spatially varying initial conditions.
+4. Run the simulation to observe wave propagation and interaction with the obstacle.
+5. Visualize the final membrane potential distribution, highlighting the non-conductive region.
+
+How to use initial conditions:
+------------------------------
+Initial conditions can be set by defining attributes that start with 'init_' 
+followed by the name of the model variable. In this example, we set 'init_h' and 'init_j' 
+to create a refractory region in the right half of the tissue. This allows us to control 
+the initial state of the tissue and observe how it affects wave propagation.
+"""
+
+
+import numpy as np
+import matplotlib.pyplot as plt
+import finitewave as fw
+
+import numpy as np
+import matplotlib.pyplot as plt
+import finitewave as fw
+
+
+n = 300
+tissue = fw.CardiacTissue([n, n])
+
+# Create a circular non-conductive obstacle in the center
+x, y = np.meshgrid(np.arange(n), np.arange(n), indexing="ij")
+
+cx, cy = n // 2, n // 2
+radius = 70
+
+obstacle = (x - cx) ** 2 + (y - cy) ** 2 < radius ** 2
+
+tissue.mesh[obstacle] = 0 # non-conductive
+
+
+stim_sequence = fw.StimSequence()
+
+# Stimulate only the left half of the upper border
+stim_sequence.add_stim(
+    fw.StimVoltageCoord(
+        time=0,
+        volt_value=20,
+        x1=0, x2=n // 2,
+        y1=0, y2=n // 2,
+    )
+)
+
+model = fw.LuoRudy91()
+model.dt = 0.01
+model.dr = 0.5
+model.t_max = 200
+
+model.cardiac_tissue = tissue
+model.stim_sequence = stim_sequence
+
+# Here we use initial conditions to keep the right half of the tissue refractory, 
+# which will prevent wave propagation in that region until it is released.
+#
+# Variables starts with init_ are used as initial conditions for the corresponding model variables.
+# You can set them to scalar values (all elements will have the same initial value) or to spatially varying arrays. 
+# In this case, we set 'h' and 'j' to arrays that are 0 in the left half and 1 in the right half, 
+# which corresponds to the refractory state in the Luo-Rudy model.
+model.init_h = np.full_like(tissue.mesh, 0.9, dtype=np.float64)
+model.init_h[n // 2:, :] = 0.0
+model.init_j = np.full_like(tissue.mesh, 0.9, dtype=np.float64)
+model.init_j[n // 2:, :] = 0.0
+
+model.run()
+
+# Show the potential map at the end of calculations.
+# We also overlay the obstacle in gray to visualize its location relative to the wave propagation.
+fig, axs = plt.subplots(ncols=1)
+plt.rcParams["axes.titley"] = -0.25
+obstacle_overlay = np.ma.masked_where(~obstacle, obstacle)
+axs.imshow(model.u)
+axs.imshow(obstacle_overlay, cmap='gray', alpha=1.0, vmin=0, vmax=1)
+plt.tight_layout()
+plt.show()

Tutorials/Reentry and Spiral waves.ipynb

@@ -141,6 +141,7 @@ def run(self, initialize=True, num_of_threads=None):
             Whether to (re)initialize the model before running the simulation.
             Default is True.
         """
+        self._check_cfl_condition()
         if initialize:
             self.initialize()
 
@@ -241,6 +242,52 @@ def clone(self):
         """
         return copy.deepcopy(self)
 
+    def reset_variables_to_defaults(self):
+        """
+        Resets the state variables to their default initial conditions.
+        """
+        for name, value in self.default_variables.items():
+            setattr(self, f"init_{name}", value)
+
+    def reset_parameters_to_defaults(self):
+        """
+        Resets the parameters to their default values.
+        """
+        for name, value in self.default_parameters.items():
+            setattr(self, name, value)
+
+    def _check_cfl_condition(self, safety_factor=0.75):
+        """
+        Checks the CFL stability condition for the explicit diffusion scheme and issues a warning 
+        if it is likely to be violated.
+        """
+        if self.D_model is None:
+            return
+
+        if self.D_model <= 0:
+            return
+
+        dt_limit = self.dr**2 / (2 * self.cardiac_tissue.dimensions * self.D_model)
+        dt_recommended = safety_factor * dt_limit
+
+        if self.dt > dt_limit:
+            warnings.warn(
+                (
+                    "The selected time step may violate the CFL stability "
+                    "condition for the explicit diffusion scheme.\n\n"
+                    f"Current values:\n"
+                    f"  dt = {self.dt}\n"
+                    f"  dr = {self.dr}\n"
+                    f"  D_max = {self.D_model}\n"
+                    f"  dimensions = {self.cardiac_tissue.dimensions}\n\n"
+                    f"Estimated stable dt <= {dt_limit:.6g}\n"
+                    f"Recommended dt <= {dt_recommended:.6g}\n\n"
+                    "Consider decreasing dt, increasing dr, or decreasing D."
+                ),
+                RuntimeWarning,
+                stacklevel=2
+            )
+        
     def _initialize_variables_and_parameters(self, ops):
         self.default_parameters = ops.get_parameters()
         self.default_variables = ops.get_variables()
@@ -249,22 +296,19 @@ def _initialize_variables_and_parameters(self, ops):
         self.state_pars = list(self.default_parameters.keys())
 
         # expose parameters as direct attributes (scalar or array)
-        for name, value in self.default_parameters.items():
-            setattr(self, name, value)
+        self.reset_parameters_to_defaults()
 
         # expose initial conditions as init_*
-        for name, value in self.default_variables.items():
-            setattr(self, f"init_{name}", value)
-
-        # declare arrays (optional, for readability/debug)
-        for name in self.default_variables.keys():
-            setattr(self, name, np.ndarray)      
+        self.reset_variables_to_defaults() 
 
     def _allocate_state_arrays(self):
         # allocate state arrays
         for name in self.default_variables.keys():
             init_val = getattr(self, f"init_{name}")
-            setattr(self, name, init_val * np.ones_like(self.u, dtype=self.npfloat))
+            if isinstance(init_val, np.ndarray):
+                setattr(self, name, init_val)
+            else:
+                setattr(self, name, init_val * np.ones_like(self.u, dtype=self.npfloat))
             if name == 'u':
                 self.u_new = self.u.copy()
 

Tutorials/RestitutionCurve.ipynb Tutorials/Restitution curve.ipynbTutorials/RestitutionCurve.ipynb renamed to Tutorials/Restitution curve.ipynb

@@ -31,7 +31,7 @@ def __init__(self):
             "u", "nai", "ki", "cai",
             "m", "h", "j",
             "oa", "oi", "ua", "ui", "xr", "xs", "d", "f", "fca",
-            "urel", "vrel", "wrel", "irel",
+            "urel", "vrel", "wrel",
             "caup", "carel",
             "nao", "ko", "cao", "R", "T", "F", "Cm",
             "gna", "gk1", "gto", "gcal", "gnab", "gcab",
@@ -67,7 +67,6 @@ def generate_body(self) -> str:
 
         urel_old = {model['urel']}
         vrel_old = {model['vrel']}
-        irel_old = {model['irel']}
         wrel_old = {model['wrel']}
 
         caup_old = {model['caup']}
@@ -200,7 +199,13 @@ def generate_body(self) -> str:
             {model['ipcamax']}, cai_old, {model['Cm']}
         )
 
-        Fn = calc_Fn(irel_old, ical, inaca, {model['F']}, {model['Vrel']})
+        irel = calc_irel(
+            urel_old, vrel_old, wrel_old,
+            {model['krel']},
+            carel_old, cai_old
+        )
+
+        Fn = calc_Fn(irel, ical, inaca, {model['F']}, {model['Vrel']})
 
         tau_urel = calc_tau_urel()
         urel_inf = calc_urel_inf(Fn)
@@ -214,12 +219,6 @@ def generate_body(self) -> str:
         wrel_inf = calc_wrel_inf(u_old)
         wrel_new = calc_gating_variable_rush_larsen(wrel_old, wrel_inf, tau_wrel, dt)
 
-        irel_new = calc_irel(
-            urel_old, vrel_old, irel_old, wrel_old,
-            {model['krel']},
-            carel_old, cai_old
-        )
-
         itr = calc_itr(caup_old, carel_old)
         iup = calc_iup({model['iupmax']}, cai_old, {model['kup']})
         iupleak = calc_iupleak(caup_old, {model['caupmax']}, {model['iupmax']})
@@ -241,15 +240,15 @@ def generate_body(self) -> str:
 
         dcai = calc_dcai(
             cai_old, inaca, ipca, ical, ibca,
-            iup, iupleak, irel_old,
+            iup, iupleak, irel,
             {model['Vrel']}, {model['Vup']},
             {model['trpnmax']}, {model['kmtrpn']},
             {model['cmdnmax']}, {model['kmcmdn']},
             {model['F']}, {model['Vj']}
         )
 
         dcarel = calc_dcarel(
-            carel_old, itr, irel_old,
+            carel_old, itr, irel,
             {model['csqnmax']}, {model['kmcsqn']}
         )
 
@@ -270,7 +269,6 @@ def generate_body(self) -> str:
 
         {model['urel']} = urel_new
         {model['vrel']} = vrel_new
-        {model['irel']} = irel_new
         {model['wrel']} = wrel_new
 
         {model['caup']} = caup_old + dt * dcaup

Tutorials/SpiralWaves2D.ipynb

@@ -130,7 +130,7 @@ def __init__(self):
         Initializes the Fenton-Karma instance with default parameters.
         """
         super().__init__()
-        self.D_model = 1.
+        self.D_model = 0.1
         self.npfloat    = 'float64'
 
         self._initialize_variables_and_parameters(ops)