Hands-on example

In this section, we define and run a simple 2D simulation of a spiral wave in a two-variable model that breaks up to spiral chaos in a circular domain with no-flux boundary conditions. With this example, it is shown in which steps a Pigreads simulation is usually defined.

First, define the coordinates of the geometry to be used. In this example, we use a 2D plane with 200 points in both \(x\) and \(y\):

import pigreads as pig
import numpy as np

R = 10
z, y, x = np.mgrid[0:1, -R:R:200j, -R:R:200j]
Nz, Ny, Nx = x.shape
dz, dy, dx = pig.deltas(z, y, x)

Pigreads is optimised for three-dimensional space. For lower-dimensional simulations, set the number of points in additional dimensions to one, as done above for the \(z\)-dimension. Note that np.mgrid is used to define a dense multi-dimensional grid where the coordinates, in this case, are defined as the range of integers from 0 up to but excluding 1 in z, i.e., \(z \in \{0\}\), and -R to R divided into 200 equal steps including both end points in y and x. See the documentation of NumPy for details. The function pigreads.helper.deltas() returns the grid spacing dz, dy, dx that is used in the arrays z, y, x, which each have a shape of (Nz, Ny, Nx). Space in Pigreads is periodic by default, such that the points at index ix = 0 are neighbours with the points at ix = 1 and ix = Nx - 1, and likewise for the other dimensions.

The grid spacing dz, dy, dx must be chosen small enough to resolve the dynamics of the system, but large enough to keep the computational cost manageable. For instance, at insufficient spatial resolution, the conduction is slower along grid axes than diagonal to them. One may use this effect to determine an as large as possible grid spacing dx such that the difference between the fastest and slowest stimulation at some distance is less than a given threshold, for instance 5%.

The integer field inhom defines which points are outside (inhom == 0) and inside the domain \(\heartsuit\) (inhom > 0). Pigreads implements no-flux boundary conditions on \(\partial\heartsuit\). In our example, we define the domain \(\heartsuit\) as a disk of radius \(R\):

inhom = np.ones(x.shape, dtype=int)
r = np.linalg.norm((x, y, z), axis=0)
inhom[r >= R] = 0

Values of inhom larger than zero can be used to select one or multiple models, i.e., reaction terms \(\underline{r}\). For an inhom value of 1, models[0] is used; and models[1] for a value of 2, etc. One or more of the available models can be selected using an instance of the pigreads.models.Models class:

models = pig.Models()
models.add("marcotte2017dynamical", beta=1.389)
Nv = models.Nv

The key of a model – an identifying string – is here used to choose a model and keyword arguments are used to set model parameters to different values than the defaults; the function add may be called multiple times to add more models. When using multiple models, the largest number of variables Nv is used.

Various models are pre-defined in Pigreads, and defining another model is straight-forward, see Models. A list of the available models can also be obtained programmatically using:

for key in pig.Models.available.keys():
    print(key)
# aliev1996simple
# barkley1991model
# beeler1977reconstruction
# ...

For some models, there are multiple sets of parameters available, they may be accessed via the additional meta data of a model definition:

key = "tentusscher2006alternans"
model_def = pig.Models.available[key]
endo = model_def.meta["parameter sets"]["endo"]
for parameter, value in endo.items():
    print(f"{parameter} = {value}")
# g_Ks = 0.392
# g_to = 0.073
# s_offset = 28.0
# s_variant = 1.0
models_ = pig.Models()
models_.add(key, **endo)

In the next step, memory needs to be allocated to store the state variables of the model. For this, we typically use a 5D array of shape (Nfr, Nz, Ny, Nx, Nv) called states in Pigreads. It consists of the values of all Nv variables at each point in space on Nfr so-called frames in time.

One way to initialise a states array is using the function pigreads.models.Models.resting_states(): It creates an array of the correct shape and fills its first frame with the appropriate resting values for each model according to the value of inhom:

Nfr = 100
states = models.resting_states(inhom, Nframes=Nfr)
assert states.shape == (Nfr, Nz, Ny, Nx, Nv)

Note that this allocates a potentially large amount of memory, as states is a 5D array of shape (Nfr, Nz, Ny, Nx, Nv). For large simulations, it may be better to only store a few frames and overwrite older frames, or to save frames to disk and not keep them in memory, for instance via memory mapping:

states = np.lib.format.open_memmap(
    "states.npy",
    mode="w+",
    dtype=np.float32,
    shape=(Nfr, Nz, Ny, Nx, Nv),
)
states[:1] = models.resting_states(inhom, Nframes=1)

If the full history over time does not need to be kept, it is also possible to create a states array of only a few frames, continually overwriting them in an infinite loop. We use this scheme in some of the interactive examples that can be found in the code repository.

Initial conditions can then be set in the frame with ifr = 0:

states[0, x < -8, 0] = 1
states[0, y < 0, 1] = 2

Or equivalently expressed in mathematical notation:

\[\begin{split}u_0(t = 0, \symbfit x) &= \begin{cases} 1 & \text{if}\, x < -8 \\ u_{0,\,\text{rest}} & \text{else} \end{cases} \\ u_1(t = 0, \symbfit x) &= \begin{cases} 2 & \text{if}\, y < 0 \\ u_{1,\,\text{rest}} & \text{else} \end{cases}\end{split}\]

The calculation of the diffusion term \(\underline{P} \nabla \cdot \symbfit D \nabla \underline{u}\) is implemented as a weighted sum of neighbouring points. The weights can be calculated using the function pigreads.models.Models.weights(), which also requires the diffusivity \(\symbfit D\) as input:

diffusivity = pig.diffusivity_matrix(Df=0.03)
weights = models.weights(dz, dy, dx, inhom, diffusivity)

See also

Documentation of the diffusivity matrix pigreads.diffusivity.diffusivity_matrix().

Calculating the weights is an expensive operation. Once they are calculated for a given geometry, it is possible to store and reuse them. Together with inhom, they fully encode all the necessary information about geometry needed for the numerical scheme below.

Finally, the simulation can be started using pigreads.models.Models.run() to advance the simulation from one frame to the next, which is done in a loop over the number of frames Nfr. The function run does Nt forward Euler steps and only returns the final states after those steps:

Nt = 200
dt = 0.025
for ifr in range(Nfr - 1):
    states[ifr + 1] = models.run(
        inhom,
        weights,
        states[ifr],
        Nt=Nt,
        dt=dt,
    )

Note that for numerical stability, the time step dt needs to be chosen small enough to at least meet the Courant-Friedrichs-Lewy (CFL) condition, which informally can be summarised as: “The wave must travel less than one grid length per time step.” A smaller time step will lead to higher accuracy at higher computational cost. We usually choose dt as large as possible while still obtaining a solution that does not change significantly for smaller dt.

Additional arguments to pigreads.models.Models.run() may be used to add stimulus currents at specific times and locations.

Now that the simulation is done, the 5D array states containing the result can be analysed and visualised, for instance with Matplotlib:

import matplotlib.pyplot as plt

plt.imshow(states[-1, 0, :, :, 0])
plt.show()

Or as a movie using FFmpeg:

from pigreads.plot import movie

movie("example.mp4", states[:, 0, :, :, 0])

Full examples with more sophisticated plotting, a progress bar, and stability checks can be found in the `examples folder in the Git repository of this project <https://gitlab.com/pigreads/pigreads/-/tree/main/examples>`__.