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-05-04 12:06:48 [info     ] Load ode /__w/fenicsx-beat/fenicsx-beat/demos/../odes/tentusscher_panfilov_2006/tentusscher_panfilov_2006_epi_cell.ode

2025-05-04 12:06:49 [info     ] Num states 19                 

2025-05-04 12:06:49 [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.36600032177187), v.min() =np.float64(-88.12451569538197)
{'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.82232230263223), v.min() =np.float64(-88.8787068710618)
{'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.351494101563965), v.min() =np.float64(-89.7372833110324)
{'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.1662476429743), v.min() =np.float64(-89.80801491982093)
{'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.41090392613747), v.min() =np.float64(-89.6123960030701)
{'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.190371141128367), v.min() =np.float64(-88.84572609961091)
{'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.430908674172642), v.min() =np.float64(-89.1471275162617)
{'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.99167142971675), v.min() =np.float64(-89.02304335390977)
{'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.774564716058055), v.min() =np.float64(-89.57966329961306)
{'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.62406249770202), v.min() =np.float64(-89.63379823891792)
{'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.541739226283976), v.min() =np.float64(-89.71751110950103)
{'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.885913233917588), v.min() =np.float64(-90.14831791112933)
{'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.21798498977755), v.min() =np.float64(-90.2330556455518)
{'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.73356122942338), v.min() =np.float64(-91.026273042075)
{'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.35006888285377), v.min() =np.float64(-91.24173909092258)
{'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.55111224293659), v.min() =np.float64(-91.81472987278796)
{'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.011458715856676), v.min() =np.float64(-91.84013094253714)
{'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.6520022275682), v.min() =np.float64(-92.23203435105722)
{'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.33648005334444), v.min() =np.float64(-92.73718674873417)
{'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.56291356894891), v.min() =np.float64(-92.84698850952601)
{'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.