Niederer benchmark

Niederer benchmark#

In this example we will use the same setup as in the Niederer benchmark [LGA+15].

from pathlib import Path
import json

from mpi4py import MPI

# import beat.conductivities
import dolfinx
import scifem
import numpy as np
import numpy.typing as npt

import pyvista
import gotranx
import matplotlib.pyplot as plt

import beat

def setup_initial_conditions() -> npt.NDArray:
    ic = {
        "V": -85.23,  # mV
        "Xr1": 0.00621,
        "Xr2": 0.4712,
        "Xs": 0.0095,
        "m": 0.00172,
        "h": 0.7444,
        "j": 0.7045,
        "d": 3.373e-05,
        "f": 0.7888,
        "f2": 0.9755,
        "fCass": 0.9953,
        "s": 0.999998,
        "r": 2.42e-08,
        "Ca_i": 0.000126,  # millimolar
        "R_prime": 0.9073,
        "Ca_SR": 3.64,  # millimolar
        "Ca_ss": 0.00036,  # millimolar
        "Na_i": 8.604,  # millimolar
        "K_i": 136.89,  # millimolar
    }
    values = model.init_state_values(**ic)
    return values

comm = MPI.COMM_WORLD

dx = 0.5
dt = 0.05
# Increase T to 100 to reproduce Niederer benchmark
T = 20.0

here = Path.cwd()

model_path = Path("tentusscher_panfilov_2006_epi_cell.py")
if not model_path.is_file():
    here = Path.cwd()
    ode = gotranx.load_ode(
        here
        / ".."
        / "odes"
        / "tentusscher_panfilov_2006"
        / "tentusscher_panfilov_2006_epi_cell.ode",
    )
    code = gotranx.cli.gotran2py.get_code(
        ode,
        scheme=[gotranx.schemes.Scheme.generalized_rush_larsen],
    )
    model_path.write_text(code)

import tentusscher_panfilov_2006_epi_cell as model

fun = model.generalized_rush_larsen
init_states = setup_initial_conditions()
parameters = model.init_parameter_values(stim_amplitude=0.0)

mesh_unit = "mm"

Lx = 20.0 * beat.units.ureg("mm").to(mesh_unit).magnitude
Ly = 7 * beat.units.ureg("mm").to(mesh_unit).magnitude
Lz = 3 * beat.units.ureg("mm").to(mesh_unit).magnitude
dx_mm = dx * beat.units.ureg("mm").to(mesh_unit).magnitude

geo = beat.geometry.get_3D_slab_geometry(
    comm=comm,
    Lx=Lx,
    Ly=Ly,
    Lz=Lz,
    dx=dx_mm,
)
ode_space = dolfinx.fem.functionspace(geo.mesh, ("Lagrange", 1))

pyvista.start_xvfb()
plotter = pyvista.Plotter()
grid = pyvista.UnstructuredGrid(*dolfinx.plot.vtk_mesh(geo.mesh))
plotter.add_mesh(grid, show_edges=True)
plotter.show_grid()
plotter.add_axes(line_width=5)
plotter.show_axes()
plotter.view_xy()

if not pyvista.OFF_SCREEN:
    plotter.show()
else:
    figure = plotter.screenshot("niederer_mesh.png")

2025-03-24 07:53:29 [info     ] Load ode /__w/fenicsx-beat/fenicsx-beat/demos/../odes/tentusscher_panfilov_2006/tentusscher_panfilov_2006_epi_cell.ode

2025-03-24 07:53:29 [info     ] Num states 19                 

2025-03-24 07:53:29 [info     ] Num parameters 53             

error: XDG_RUNTIME_DIR is invalid or not set in the environment.

# Surface to volume ratio
conductivities = beat.conductivities.default_conductivities("Niederer")
# # Membrane capacitance
C_m = 1.0 * beat.units.ureg("uF/cm**2")

time = dolfinx.fem.Constant(geo.mesh, 0.0)
L = 1.5 * beat.units.ureg("mm").to(mesh_unit).magnitude
S1_marker = 1
L = 1.5


tol = 1.0e-10


def S1_subdomain(x):
    return np.logical_and(
        np.logical_and(x[0] <= L + tol, x[1] <= L + tol),
        x[2] <= L + tol,
    )


cells = dolfinx.mesh.locate_entities(geo.mesh, geo.mesh.topology.dim, S1_subdomain)
S1_markers = dolfinx.mesh.meshtags(
    geo.mesh,
    geo.mesh.topology.dim,
    cells,
    np.full(len(cells), S1_marker, dtype=np.int32),
)

I_s = beat.stimulation.define_stimulus(
    mesh=geo.mesh,
    chi=conductivities["chi"],
    time=time,
    subdomain_data=S1_markers,
    marker=S1_marker,
    mesh_unit=mesh_unit,
    amplitude=50_000.0,
)


M = beat.conductivities.define_conductivity_tensor(
    f0=geo.f0,
    **conductivities,
)

params = {"preconditioner": "sor", "use_custom_preconditioner": False}

pde = beat.MonodomainModel(
    time=time,
    mesh=geo.mesh,
    M=M,
    I_s=I_s,
    params=params,
    C_m=C_m.to(f"uF/{mesh_unit}**2").magnitude,
    dx=I_s.dZ,
)
ode = beat.odesolver.DolfinODESolver(
    v_ode=dolfinx.fem.Function(ode_space),
    v_pde=pde.state,
    fun=fun,
    init_states=init_states,
    parameters=parameters,
    num_states=len(init_states),
    v_index=model.state_index("V"),
)

solver = beat.MonodomainSplittingSolver(pde=pde, ode=ode)
output_dir = Path("results-niederer-benchmark")
output_dir.mkdir(exist_ok=True)
filename = output_dir / f"results-{dt}-{dx}.xdmf"
filename.unlink(missing_ok=True)
filename.with_suffix(".h5").unlink(missing_ok=True)

vtx = dolfinx.io.VTXWriter(
    comm,
    "niederer_benchmark.bp",
    [solver.pde.state],
    engine="BP4",
)

points = {
    "P1": (0, 0, 0),
    "P2": (0.0, Ly, 0.0),
    "P3": (Lx, 0.0, 0.0),
    "P4": (Lx, Ly, 0.0),
    "P5": (0.0, 0.0, Lz),
    "P6": (0.0, Ly, Lz),
    "P7": (Lx, 0.0, Lz),
    "P8": (Lx, Ly, Lz),
    "P9": (Lx / 2, Ly / 2, Lz / 2),
}

activation_times = {p: -1.0 for p in points}
save_freq = int(1.0 / dt)
i = 0
plotter_voltage = pyvista.Plotter()
viridis = plt.get_cmap("viridis")
grid.point_data["V"] = solver.pde.state.x.array
grid.set_active_scalars("V")
renderer = plotter_voltage.add_mesh(
    grid,
    show_edges=True,
    lighting=False,
    cmap=viridis,
    clim=[-90.0, 40.0],
)
gif_file = Path("niederer_benchmark.gif")
gif_file.unlink(missing_ok=True)
plotter_voltage.open_gif(gif_file.as_posix())

T = 20
# T = 100  # Change to 100 to reproduce Niederer benchmark
t = 0.0
while t < T + 1e-12 and any(at < 0.0 for at in activation_times.values()):
    v = solver.pde.state.x.array
    if i % save_freq == 0:
        print(f"Solve for {t=:.2f}, {v.max() =}, {v.min() =}")
        print(activation_times)
        vtx.write(t)
        grid.point_data["V"] = solver.pde.state.x.array
        plotter_voltage.write_frame()
    solver.step((t, t + dt))

    for p in points:
        value = scifem.evaluate_function(solver.pde.state, [points[p]]).squeeze()
        if value > 0.0 and activation_times[p] < 0.0:
            activation_times[p] = t
    i += 1
    t += dt

plotter_voltage.close()

Solve for t=0.00, v.max() =np.float64(-85.23), v.min() =np.float64(-85.23)
{'P1': -1.0, 'P2': -1.0, 'P3': -1.0, 'P4': -1.0, 'P5': -1.0, 'P6': -1.0, 'P7': -1.0, 'P8': -1.0, 'P9': -1.0}

error: XDG_RUNTIME_DIR is invalid or not set in the environment.

Solve for t=1.00, v.max() =np.float64(-27.365497242149793), v.min() =np.float64(-88.12835714910064)
{'P1': -1.0, 'P2': -1.0, 'P3': -1.0, 'P4': -1.0, 'P5': -1.0, 'P6': -1.0, 'P7': -1.0, 'P8': -1.0, 'P9': -1.0}

Solve for t=2.00, v.max() =np.float64(61.81849085607121), v.min() =np.float64(-88.87815768175166)
{'P1': 1.2500000000000004, 'P2': -1.0, 'P3': -1.0, 'P4': -1.0, 'P5': -1.0, 'P6': -1.0, 'P7': -1.0, 'P8': -1.0, 'P9': -1.0}

Solve for t=3.00, v.max() =np.float64(41.34981696963944), v.min() =np.float64(-89.7364447528398)
{'P1': 1.2500000000000004, 'P2': -1.0, 'P3': -1.0, 'P4': -1.0, 'P5': -1.0, 'P6': -1.0, 'P7': -1.0, 'P8': -1.0, 'P9': -1.0}

Solve for t=4.00, v.max() =np.float64(30.161805433955415), v.min() =np.float64(-89.80917134206945)
{'P1': 1.2500000000000004, 'P2': -1.0, 'P3': -1.0, 'P4': -1.0, 'P5': -1.0, 'P6': -1.0, 'P7': -1.0, 'P8': -1.0, 'P9': -1.0}

Solve for t=5.00, v.max() =np.float64(32.40973734877997), v.min() =np.float64(-89.61056456259517)
{'P1': 1.2500000000000004, 'P2': -1.0, 'P3': -1.0, 'P4': -1.0, 'P5': -1.0, 'P6': -1.0, 'P7': -1.0, 'P8': -1.0, 'P9': -1.0}

Solve for t=6.00, v.max() =np.float64(31.184536665183742), v.min() =np.float64(-88.84192460298502)
{'P1': 1.2500000000000004, 'P2': -1.0, 'P3': -1.0, 'P4': -1.0, 'P5': -1.0, 'P6': -1.0, 'P7': -1.0, 'P8': -1.0, 'P9': -1.0}

Solve for t=7.00, v.max() =np.float64(29.419247394071423), v.min() =np.float64(-89.14510444472783)
{'P1': 1.2500000000000004, 'P2': -1.0, 'P3': -1.0, 'P4': -1.0, 'P5': -1.0, 'P6': -1.0, 'P7': -1.0, 'P8': -1.0, 'P9': -1.0}

Solve for t=8.00, v.max() =np.float64(33.97757153482185), v.min() =np.float64(-89.01983433203414)
{'P1': 1.2500000000000004, 'P2': -1.0, 'P3': -1.0, 'P4': -1.0, 'P5': -1.0, 'P6': -1.0, 'P7': -1.0, 'P8': -1.0, 'P9': -1.0}

Solve for t=9.00, v.max() =np.float64(28.772046673645136), v.min() =np.float64(-89.58268722004469)
{'P1': 1.2500000000000004, 'P2': -1.0, 'P3': -1.0, 'P4': -1.0, 'P5': -1.0, 'P6': -1.0, 'P7': -1.0, 'P8': -1.0, 'P9': -1.0}

Solve for t=10.00, v.max() =np.float64(31.61169662091426), v.min() =np.float64(-89.63493241557963)
{'P1': 1.2500000000000004, 'P2': -1.0, 'P3': -1.0, 'P4': -1.0, 'P5': -1.0, 'P6': -1.0, 'P7': -1.0, 'P8': -1.0, 'P9': -1.0}

Solve for t=11.00, v.max() =np.float64(30.53473548348199), v.min() =np.float64(-89.72261684217871)
{'P1': 1.2500000000000004, 'P2': -1.0, 'P3': -1.0, 'P4': -1.0, 'P5': -1.0, 'P6': -1.0, 'P7': -1.0, 'P8': -1.0, 'P9': -1.0}

Solve for t=12.00, v.max() =np.float64(27.869778477038494), v.min() =np.float64(-90.14736054688737)
{'P1': 1.2500000000000004, 'P2': -1.0, 'P3': -1.0, 'P4': -1.0, 'P5': -1.0, 'P6': -1.0, 'P7': -1.0, 'P8': -1.0, 'P9': -1.0}

Solve for t=13.00, v.max() =np.float64(36.231030844372924), v.min() =np.float64(-90.23053287575775)
{'P1': 1.2500000000000004, 'P2': -1.0, 'P3': -1.0, 'P4': -1.0, 'P5': -1.0, 'P6': -1.0, 'P7': -1.0, 'P8': -1.0, 'P9': -1.0}

Solve for t=14.00, v.max() =np.float64(38.72168598131495), v.min() =np.float64(-91.00630332024203)
{'P1': 1.2500000000000004, 'P2': -1.0, 'P3': -1.0, 'P4': -1.0, 'P5': -1.0, 'P6': -1.0, 'P7': -1.0, 'P8': -1.0, 'P9': -1.0}

Solve for t=15.00, v.max() =np.float64(36.337892961807235), v.min() =np.float64(-91.26670861686922)
{'P1': 1.2500000000000004, 'P2': -1.0, 'P3': -1.0, 'P4': -1.0, 'P5': 14.000000000000064, 'P6': -1.0, 'P7': -1.0, 'P8': -1.0, 'P9': -1.0}

Solve for t=16.00, v.max() =np.float64(39.54544934028266), v.min() =np.float64(-91.82584090545595)
{'P1': 1.2500000000000004, 'P2': -1.0, 'P3': -1.0, 'P4': -1.0, 'P5': 14.000000000000064, 'P6': -1.0, 'P7': -1.0, 'P8': -1.0, 'P9': -1.0}

Solve for t=17.00, v.max() =np.float64(35.00655295067564), v.min() =np.float64(-91.82813057493661)
{'P1': 1.2500000000000004, 'P2': -1.0, 'P3': -1.0, 'P4': -1.0, 'P5': 14.000000000000064, 'P6': -1.0, 'P7': -1.0, 'P8': -1.0, 'P9': -1.0}

Solve for t=18.00, v.max() =np.float64(34.6482815390785), v.min() =np.float64(-92.22360466706088)
{'P1': 1.2500000000000004, 'P2': -1.0, 'P3': -1.0, 'P4': -1.0, 'P5': 14.000000000000064, 'P6': -1.0, 'P7': -1.0, 'P8': -1.0, 'P9': -1.0}

Solve for t=19.00, v.max() =np.float64(32.12200578162016), v.min() =np.float64(-92.73689682071492)
{'P1': 1.2500000000000004, 'P2': -1.0, 'P3': -1.0, 'P4': -1.0, 'P5': 14.000000000000064, 'P6': -1.0, 'P7': -1.0, 'P8': -1.0, 'P9': -1.0}

Solve for t=20.00, v.max() =np.float64(32.560529686553984), v.min() =np.float64(-92.84581253319969)
{'P1': 1.2500000000000004, 'P2': -1.0, 'P3': -1.0, 'P4': -1.0, 'P5': 14.000000000000064, 'P6': -1.0, 'P7': -1.0, 'P8': -1.0, 'P9': -1.0}

_

# Save activation times
activation_times["dx"] = dx
activation_times["dt"] = dt
at_file_name = output_dir / "activation_times.json"
if at_file_name.is_file():
    all_at = json.loads(at_file_name.read_text())
else:
    all_at = []
all_at.append(activation_times)
at_file_name.write_text(json.dumps(all_at, indent=2))

213

The activation times are saved in the file output-niederer-benchmark/activation_times.json. The file contains a list of dictionaries, each dictionary contains the activation times for a specific dx and dt. Here are the activation times for the different dx and dt:

dx

dt

P1

P2

P3

P4

P5

P6

P7

P8

P9

0

0.5

0.05

1.25

51.1

34.9

58.9

14.1

49.5

34

56.65

26.05

1

0.5

0.01

1.22

50.85

33.96

58.05

13.98

49.36

33.07

55.91

25.64

2

0.5

0.005

1.215

50.775

33.825

57.96

13.97

49.345

32.945

55.825

25.595

3

0.2

0.05

1.25

29.7

32.9

40.2

9.55

30

32.95

39.9

18.9

4

0.2

0.01

1.24

29.09

31.25

38.66

9.34

29.4

31.29

38.42

18.14

5

0.2

0.005

1.235

29.015

31.05

38.475

9.315

29.32

31.08

38.235

18.045

6

0.1

0.05

1.25

26.85

33.3

40.35

8.4

27.5

33.85

40.55

18.95

7

0.1

0.01

1.23

25.64

31.46

38.08

8.03

26.24

31.94

38.21

17.95

8

0.1

0.005

1.225

25.5

31.26

37.81

7.99

26.09

31.72

37.93

17.835
[LGA+15]

Sander Land, Viatcheslav Gurev, Sander Arens, Christoph M Augustin, Lukas Baron, Robert Blake, Chris Bradley, Sebastian Castro, Andrew Crozier, Marco Favino, and others. Verification of cardiac mechanics software: benchmark problems and solutions for testing active and passive material behavior. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 471(2184):20150641, 2015.

[ZLH+21]

Zhaoyang Zhang, Michael B Liu, Xiaodong Huang, Zhen Song, and Zhilin Qu. Mechanisms of premature ventricular complexes caused by qt prolongation. Biophysical journal, 120(2):352–369, 2021.