Niederer benchmark

Niederer benchmark#

In this example we will use the same setup as in the Niederer benchmark [LGA+15]. It is a standard convergence/cross-code benchmark for monodomain solvers: a small cuboid slab of tissue with a detailed ionic current model (the ten Tusscher–Panfilov 2006 epicardial model), a conductivity tensor \(M\) built from the “Niederer” literature conductivity values, and a stimulus \(I_{stim}\) applied in one corner. The quantity of interest is the local activation time — the first time \(v\) crosses \(0\) mV — at nine fixed points in the slab, for a range of spatial and temporal resolutions \(dx\), \(dt\). See the mathematical background page for the underlying monodomain model and notation used below.

from pathlib import Path
import json
import time
import logging

from mpi4py import MPI

import os
import dolfinx
import scifem
import numpy as np
import numpy.typing as npt

logging.basicConfig(level=logging.INFO)
logger = logging.getLogger(__name__)

if os.environ.get("NO_PYVISTA", "0") == "1":
    pyvista = None
    logger.warning("Turn off pyvista")
else:
    try:
        import pyvista  # type: ignore[no-redef, assignment]
    except ImportError:
        pyvista = None  # type: ignore[no-redef]
        logger.warning("pyvista not installed, skipping visualization")

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()
    cell_ode = gotranx.load_ode(
        here
        / ".."
        / "odes"
        / "tentusscher_panfilov_2006"
        / "tentusscher_panfilov_2006_epi_cell.ode",
    )
    code = gotranx.cli.gotran2py.get_code(
        cell_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))

if pyvista is not None:
    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")

INFO:trame_server.utils.namespace:Translator(prefix=None)

INFO:wslink.backends.aiohttp:awaiting runner setup

INFO:wslink.backends.aiohttp:awaiting site startup

INFO:wslink.backends.aiohttp:Print WSLINK_READY_MSG

INFO:wslink.backends.aiohttp:Schedule auto shutdown with timeout 0

INFO:wslink.backends.aiohttp:awaiting running future

# Literature values for the intracellular/extracellular conductivities and the surface to volume
# ratio $\chi$, used below to build the conductivity tensor $M$.
conductivities = beat.conductivities.default_conductivities("Niederer")
# Membrane capacitance $C_m$
C_m = beat.units.ureg.Quantity(1.0, "uF/cm**2")

time_constant = 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),
)

# The stimulus current $I_{stim}$, applied to the S1 corner subdomain defined above.

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

# The conductivity tensor $M$, built from the fibre direction `geo.f0` and the conductivities above.

assert geo.f0 is not None
M = beat.conductivities.define_conductivity_tensor(
    f0=geo.f0,
    **conductivities,
)

params = {
    "petsc_options": {
        "ksp_type": "cg",
        "pc_type": "hypre",
        "pc_hypre_type": "boomeramg",
    },
}

save_freq = int(1.0 / dt)
# Create a monitor to log the solver time and performance metrics every `save_freq` time steps
monitor = beat.PerformanceMonitor(log_frequency=save_freq)


pde = beat.MonodomainModel(
    time=time_constant,
    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,
    monitor=monitor,
)
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"),
    monitor=monitor,
)

# Combine the PDE and ODE solver, using the default $\theta = 1$ (Godunov splitting).

INFO:beat.conductivities:Get harmonic mean conductivity: g_il=Quantity(0.17, "siemens / meter") g_it=Quantity(0.019, "siemens / meter") g_el=Quantity(0.62, "siemens / meter") g_et=Quantity(0.24, "siemens / meter") chi=Quantity(1400.0, "1 / centimeter")

INFO:beat.conductivities:Harmonic mean conductivities sigma_l=Quantity(0.13341772151898734, "siemens / meter") sigma_t=Quantity(0.017606177606177605, "siemens / meter")

INFO:beat.conductivities:Scaled harmonic mean conductivities s_l=0.0009529837251356239 s_t=0.00012575841147269718

INFO:beat.conductivities:Define conductivity tensor s_l=0.0009529837251356239 s_t=0.00012575841147269718 dim=3

solver = beat.MonodomainSplittingSolver(pde=pde, ode=ode, monitor=monitor)
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.v],
    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
if pyvista is not None:
    plotter_voltage = pyvista.Plotter()
    viridis = plt.get_cmap("viridis")
    grid.point_data["V"] = solver.pde.v.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
times: list[float] = []
while t < T + 1e-12 and any(at < 0.0 for at in activation_times.values()):
    v = solver.pde.v.x.array
    if i % save_freq == 0:
        logger.info(f"Solve for {t=:.2f}, {v.max() =}, {v.min() =}")
        if len(times) > 0:
            logger.info(f"Average solver time per time step: {np.mean(times):.5f} s")
        logger.info(activation_times)
        vtx.write(t)
        if pyvista is not None:
            grid.point_data["V"] = solver.pde.v.x.array
            plotter_voltage.write_frame()
    t0 = time.perf_counter()
    solver.step((t, t + dt))
    times.append(time.perf_counter() - t0)
    for p in points:
        value = scifem.evaluate_function(solver.pde.v, [points[p]]).squeeze()
        if value > 0.0 and activation_times[p] < 0.0:
            activation_times[p] = t
    i += 1
    t += dt

if pyvista is not None:
    plotter_voltage.close()

monitor.display_summary()
monitor.save_summary(output_dir / "performance_summary.json")

INFO:__main__:Solve for t=0.00, v.max() =np.float64(-85.23), v.min() =np.float64(-85.23)

INFO:__main__:{'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}

INFO:beat.telemetry:PDE step timing step=20, t=(0.45000, 0.50000), ksp_iterations=5, ksp_residual_norm=6.221910e-03, ksp_converged_reason=2, ode_function_call=0.026692s, ode_state_update=0.000108s, ode_total_step=0.026860s, ode_step=0.029692s, ode_to_dolfin=0.000071s, ode_to_pde=0.000066s, pde_assign_previous_before=0.000029s, pde_set_time=0.000043s, pde_update_matrices=0.010126s, pde_update_rhs=0.001134s, pde_linear_solve=0.005864s, pde_scatter_forward=0.000027s, pde_total_step=0.017327s, pde_step=0.017353s, pde_to_ode=0.000064s, ode_from_dolfin=0.000020s, pde_assign_previous_after=0.000026s, total_step=0.047751s

INFO:beat.telemetry:PDE step timing step=40, t=(0.95000, 1.00000), ksp_iterations=5, ksp_residual_norm=6.202186e-03, ksp_converged_reason=2, ode_function_call=0.052675s, ode_state_update=0.000212s, ode_total_step=0.052997s, ode_step=0.058934s, ode_to_dolfin=0.000128s, ode_to_pde=0.000111s, pde_assign_previous_before=0.000051s, pde_set_time=0.000077s, pde_update_matrices=0.010126s, pde_update_rhs=0.002043s, pde_linear_solve=0.010631s, pde_scatter_forward=0.000040s, pde_total_step=0.023155s, pde_step=0.023202s, pde_to_ode=0.000123s, ode_from_dolfin=0.000039s, pde_assign_previous_after=0.000050s, total_step=0.083450s

INFO:__main__:Solve for t=1.00, v.max() =np.float64(-27.608542903240274), v.min() =np.float64(-88.1310920950999)

INFO:__main__:Average solver time per time step: 0.00424 s

INFO:__main__:{'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}

INFO:beat.telemetry:PDE step timing step=60, t=(1.45000, 1.50000), ksp_iterations=5, ksp_residual_norm=6.204528e-03, ksp_converged_reason=2, ode_function_call=0.078966s, ode_state_update=0.000326s, ode_total_step=0.079455s, ode_step=0.088577s, ode_to_dolfin=0.000192s, ode_to_pde=0.000157s, pde_assign_previous_before=0.000076s, pde_set_time=0.000115s, pde_update_matrices=0.010126s, pde_update_rhs=0.003043s, pde_linear_solve=0.015359s, pde_scatter_forward=0.000056s, pde_total_step=0.029047s, pde_step=0.029117s, pde_to_ode=0.000181s, ode_from_dolfin=0.000058s, pde_assign_previous_after=0.000075s, total_step=0.119639s

INFO:beat.telemetry:PDE step timing step=80, t=(1.95000, 2.00000), ksp_iterations=5, ksp_residual_norm=6.212893e-03, ksp_converged_reason=2, ode_function_call=0.104676s, ode_state_update=0.000429s, ode_total_step=0.105316s, ode_step=0.117618s, ode_to_dolfin=0.000244s, ode_to_pde=0.000200s, pde_assign_previous_before=0.000097s, pde_set_time=0.000146s, pde_update_matrices=0.010126s, pde_update_rhs=0.003870s, pde_linear_solve=0.020001s, pde_scatter_forward=0.000068s, pde_total_step=0.034656s, pde_step=0.034746s, pde_to_ode=0.000237s, ode_from_dolfin=0.000075s, pde_assign_previous_after=0.000099s, total_step=0.154885s

INFO:__main__:Solve for t=2.00, v.max() =np.float64(60.29388459312645), v.min() =np.float64(-88.87722531862953)

INFO:__main__:Average solver time per time step: 0.00393 s

INFO:__main__:{'P1': 1.3000000000000005, 'P2': -1.0, 'P3': -1.0, 'P4': -1.0, 'P5': -1.0, 'P6': -1.0, 'P7': -1.0, 'P8': -1.0, 'P9': -1.0}

INFO:beat.telemetry:PDE step timing step=100, t=(2.45000, 2.50000), ksp_iterations=5, ksp_residual_norm=6.198655e-03, ksp_converged_reason=2, ode_function_call=0.130672s, ode_state_update=0.000534s, ode_total_step=0.131469s, ode_step=0.147058s, ode_to_dolfin=0.000308s, ode_to_pde=0.000245s, pde_assign_previous_before=0.000121s, pde_set_time=0.000182s, pde_update_matrices=0.010126s, pde_update_rhs=0.005034s, pde_linear_solve=0.024731s, pde_scatter_forward=0.000083s, pde_total_step=0.040713s, pde_step=0.040825s, pde_to_ode=0.000296s, ode_from_dolfin=0.000095s, pde_assign_previous_after=0.000124s, total_step=0.191031s

INFO:beat.telemetry:PDE step timing step=120, t=(2.95000, 3.00000), ksp_iterations=5, ksp_residual_norm=6.199786e-03, ksp_converged_reason=2, ode_function_call=0.156286s, ode_state_update=0.000637s, ode_total_step=0.157234s, ode_step=0.175952s, ode_to_dolfin=0.000359s, ode_to_pde=0.000288s, pde_assign_previous_before=0.000143s, pde_set_time=0.000215s, pde_update_matrices=0.010126s, pde_update_rhs=0.005877s, pde_linear_solve=0.029404s, pde_scatter_forward=0.000095s, pde_total_step=0.046373s, pde_step=0.046506s, pde_to_ode=0.000351s, ode_from_dolfin=0.000113s, pde_assign_previous_after=0.000148s, total_step=0.226164s

INFO:__main__:Solve for t=3.00, v.max() =np.float64(40.813109935474685), v.min() =np.float64(-89.5941091366833)

INFO:__main__:Average solver time per time step: 0.00383 s

INFO:__main__:{'P1': 1.3000000000000005, 'P2': -1.0, 'P3': -1.0, 'P4': -1.0, 'P5': -1.0, 'P6': -1.0, 'P7': -1.0, 'P8': -1.0, 'P9': -1.0}

INFO:beat.telemetry:PDE step timing step=140, t=(3.45000, 3.50000), ksp_iterations=5, ksp_residual_norm=6.202101e-03, ksp_converged_reason=2, ode_function_call=0.182376s, ode_state_update=0.000742s, ode_total_step=0.183481s, ode_step=0.205342s, ode_to_dolfin=0.000428s, ode_to_pde=0.000336s, pde_assign_previous_before=0.000168s, pde_set_time=0.000254s, pde_update_matrices=0.010126s, pde_update_rhs=0.006914s, pde_linear_solve=0.034204s, pde_scatter_forward=0.000111s, pde_total_step=0.052378s, pde_step=0.052536s, pde_to_ode=0.000412s, ode_from_dolfin=0.000133s, pde_assign_previous_after=0.000174s, total_step=0.262240s

INFO:beat.telemetry:PDE step timing step=160, t=(3.95000, 4.00000), ksp_iterations=5, ksp_residual_norm=6.204700e-03, ksp_converged_reason=2, ode_function_call=0.208401s, ode_state_update=0.000846s, ode_total_step=0.209662s, ode_step=0.234646s, ode_to_dolfin=0.000487s, ode_to_pde=0.000382s, pde_assign_previous_before=0.000191s, pde_set_time=0.000294s, pde_update_matrices=0.010126s, pde_update_rhs=0.008018s, pde_linear_solve=0.038997s, pde_scatter_forward=0.000125s, pde_total_step=0.058442s, pde_step=0.058622s, pde_to_ode=0.000472s, ode_from_dolfin=0.000152s, pde_assign_previous_after=0.000199s, total_step=0.298265s

INFO:__main__:Solve for t=4.00, v.max() =np.float64(29.926891976603628), v.min() =np.float64(-89.78076447962178)

INFO:__main__:Average solver time per time step: 0.00379 s

INFO:__main__:{'P1': 1.3000000000000005, 'P2': -1.0, 'P3': -1.0, 'P4': -1.0, 'P5': -1.0, 'P6': -1.0, 'P7': -1.0, 'P8': -1.0, 'P9': -1.0}

INFO:beat.telemetry:PDE step timing step=180, t=(4.45000, 4.50000), ksp_iterations=5, ksp_residual_norm=6.199545e-03, ksp_converged_reason=2, ode_function_call=0.247253s, ode_state_update=0.000956s, ode_total_step=0.248680s, ode_step=0.276886s, ode_to_dolfin=0.000555s, ode_to_pde=0.000431s, pde_assign_previous_before=0.000216s, pde_set_time=0.000333s, pde_update_matrices=0.010126s, pde_update_rhs=0.009090s, pde_linear_solve=0.043839s, pde_scatter_forward=0.000141s, pde_total_step=0.064524s, pde_step=0.064726s, pde_to_ode=0.000533s, ode_from_dolfin=0.000171s, pde_assign_previous_after=0.000225s, total_step=0.347277s

INFO:beat.telemetry:PDE step timing step=200, t=(4.95000, 5.00000), ksp_iterations=5, ksp_residual_norm=6.198219e-03, ksp_converged_reason=2, ode_function_call=0.287523s, ode_state_update=0.001061s, ode_total_step=0.289106s, ode_step=0.320458s, ode_to_dolfin=0.000607s, ode_to_pde=0.000475s, pde_assign_previous_before=0.000238s, pde_set_time=0.000365s, pde_update_matrices=0.010126s, pde_update_rhs=0.009927s, pde_linear_solve=0.048634s, pde_scatter_forward=0.000153s, pde_total_step=0.070300s, pde_step=0.070523s, pde_to_ode=0.000590s, ode_from_dolfin=0.000189s, pde_assign_previous_after=0.000250s, total_step=0.397219s

INFO:__main__:Solve for t=5.00, v.max() =np.float64(31.362948246090834), v.min() =np.float64(-89.4201789074145)

INFO:__main__:Average solver time per time step: 0.00403 s

INFO:__main__:{'P1': 1.3000000000000005, 'P2': -1.0, 'P3': -1.0, 'P4': -1.0, 'P5': -1.0, 'P6': -1.0, 'P7': -1.0, 'P8': -1.0, 'P9': -1.0}

INFO:beat.telemetry:PDE step timing step=220, t=(5.45000, 5.50000), ksp_iterations=5, ksp_residual_norm=6.198736e-03, ksp_converged_reason=2, ode_function_call=0.313563s, ode_state_update=0.001167s, ode_total_step=0.315304s, ode_step=0.349836s, ode_to_dolfin=0.000671s, ode_to_pde=0.000521s, pde_assign_previous_before=0.000262s, pde_set_time=0.000404s, pde_update_matrices=0.010126s, pde_update_rhs=0.010958s, pde_linear_solve=0.053434s, pde_scatter_forward=0.000169s, pde_total_step=0.076297s, pde_step=0.076544s, pde_to_ode=0.000651s, ode_from_dolfin=0.000209s, pde_assign_previous_after=0.000276s, total_step=0.433282s

INFO:beat.telemetry:PDE step timing step=240, t=(5.95000, 6.00000), ksp_iterations=5, ksp_residual_norm=6.187648e-03, ksp_converged_reason=2, ode_function_call=0.340861s, ode_state_update=0.001284s, ode_total_step=0.342793s, ode_step=0.380974s, ode_to_dolfin=0.000759s, ode_to_pde=0.000580s, pde_assign_previous_before=0.000296s, pde_set_time=0.000458s, pde_update_matrices=0.010126s, pde_update_rhs=0.012402s, pde_linear_solve=0.058491s, pde_scatter_forward=0.000189s, pde_total_step=0.083017s, pde_step=0.083299s, pde_to_ode=0.000738s, ode_from_dolfin=0.000231s, pde_assign_previous_after=0.000304s, total_step=0.472041s

INFO:__main__:Solve for t=6.00, v.max() =np.float64(31.406417554027634), v.min() =np.float64(-89.17776874023409)

INFO:__main__:Average solver time per time step: 0.00400 s

INFO:__main__:{'P1': 1.3000000000000005, 'P2': -1.0, 'P3': -1.0, 'P4': -1.0, 'P5': -1.0, 'P6': -1.0, 'P7': -1.0, 'P8': -1.0, 'P9': -1.0}

INFO:beat.telemetry:PDE step timing step=260, t=(6.45000, 6.50000), ksp_iterations=5, ksp_residual_norm=6.195862e-03, ksp_converged_reason=2, ode_function_call=0.379415s, ode_state_update=0.001402s, ode_total_step=0.381517s, ode_step=0.422887s, ode_to_dolfin=0.000822s, ode_to_pde=0.000626s, pde_assign_previous_before=0.000319s, pde_set_time=0.000497s, pde_update_matrices=0.010126s, pde_update_rhs=0.013621s, pde_linear_solve=0.063282s, pde_scatter_forward=0.000209s, pde_total_step=0.089200s, pde_step=0.089507s, pde_to_ode=0.000798s, ode_from_dolfin=0.000250s, pde_assign_previous_after=0.000328s, total_step=0.520832s

INFO:beat.telemetry:PDE step timing step=280, t=(6.95000, 7.00000), ksp_iterations=5, ksp_residual_norm=6.181814e-03, ksp_converged_reason=2, ode_function_call=0.418949s, ode_state_update=0.001507s, ode_total_step=0.421207s, ode_step=0.465660s, ode_to_dolfin=0.000874s, ode_to_pde=0.000669s, pde_assign_previous_before=0.000342s, pde_set_time=0.000531s, pde_update_matrices=0.010126s, pde_update_rhs=0.014523s, pde_linear_solve=0.067984s, pde_scatter_forward=0.000222s, pde_total_step=0.094955s, pde_step=0.095285s, pde_to_ode=0.000853s, ode_from_dolfin=0.000268s, pde_assign_previous_after=0.000352s, total_step=0.569968s

INFO:__main__:Solve for t=7.00, v.max() =np.float64(29.585741063385573), v.min() =np.float64(-88.31355090839853)

INFO:__main__:Average solver time per time step: 0.00413 s

INFO:__main__:{'P1': 1.3000000000000005, 'P2': -1.0, 'P3': -1.0, 'P4': -1.0, 'P5': -1.0, 'P6': -1.0, 'P7': -1.0, 'P8': -1.0, 'P9': -1.0}

INFO:beat.telemetry:PDE step timing step=300, t=(7.45000, 7.50000), ksp_iterations=5, ksp_residual_norm=6.192413e-03, ksp_converged_reason=2, ode_function_call=0.445525s, ode_state_update=0.001617s, ode_total_step=0.447953s, ode_step=0.495806s, ode_to_dolfin=0.000953s, ode_to_pde=0.000721s, pde_assign_previous_before=0.000370s, pde_set_time=0.000578s, pde_update_matrices=0.010126s, pde_update_rhs=0.015798s, pde_linear_solve=0.072858s, pde_scatter_forward=0.000241s, pde_total_step=0.101302s, pde_step=0.101662s, pde_to_ode=0.000919s, ode_from_dolfin=0.000302s, pde_assign_previous_after=0.000381s, total_step=0.607278s

INFO:beat.telemetry:PDE step timing step=320, t=(7.95000, 8.00000), ksp_iterations=5, ksp_residual_norm=6.272298e-03, ksp_converged_reason=2, ode_function_call=0.471502s, ode_state_update=0.001722s, ode_total_step=0.474087s, ode_step=0.525214s, ode_to_dolfin=0.001011s, ode_to_pde=0.000767s, pde_assign_previous_before=0.000395s, pde_set_time=0.000623s, pde_update_matrices=0.010126s, pde_update_rhs=0.016962s, pde_linear_solve=0.077624s, pde_scatter_forward=0.000256s, pde_total_step=0.107416s, pde_step=0.107803s, pde_to_ode=0.000982s, ode_from_dolfin=0.000322s, pde_assign_previous_after=0.000407s, total_step=0.643538s

INFO:__main__:Solve for t=8.00, v.max() =np.float64(33.040237216333026), v.min() =np.float64(-89.1236498084329)

INFO:__main__:Average solver time per time step: 0.00408 s

INFO:__main__:{'P1': 1.3000000000000005, 'P2': -1.0, 'P3': -1.0, 'P4': -1.0, 'P5': -1.0, 'P6': -1.0, 'P7': -1.0, 'P8': -1.0, 'P9': -1.0}

INFO:beat.telemetry:PDE step timing step=340, t=(8.45000, 8.50000), ksp_iterations=5, ksp_residual_norm=6.364465e-03, ksp_converged_reason=2, ode_function_call=0.509901s, ode_state_update=0.001826s, ode_total_step=0.512651s, ode_step=0.566927s, ode_to_dolfin=0.001082s, ode_to_pde=0.000832s, pde_assign_previous_before=0.000422s, pde_set_time=0.000669s, pde_update_matrices=0.010126s, pde_update_rhs=0.018288s, pde_linear_solve=0.082436s, pde_scatter_forward=0.000275s, pde_total_step=0.113772s, pde_step=0.114191s, pde_to_ode=0.001050s, ode_from_dolfin=0.000343s, pde_assign_previous_after=0.000434s, total_step=0.692401s

INFO:beat.telemetry:PDE step timing step=360, t=(8.95000, 9.00000), ksp_iterations=5, ksp_residual_norm=6.382525e-03, ksp_converged_reason=2, ode_function_call=0.549315s, ode_state_update=0.001928s, ode_total_step=0.552220s, ode_step=0.609605s, ode_to_dolfin=0.001145s, ode_to_pde=0.000878s, pde_assign_previous_before=0.000446s, pde_set_time=0.000709s, pde_update_matrices=0.010126s, pde_update_rhs=0.019341s, pde_linear_solve=0.087116s, pde_scatter_forward=0.000291s, pde_total_step=0.119680s, pde_step=0.120122s, pde_to_ode=0.001105s, ode_from_dolfin=0.000362s, pde_assign_previous_after=0.000458s, total_step=0.741660s

INFO:__main__:Solve for t=9.00, v.max() =np.float64(28.104588035486408), v.min() =np.float64(-89.51789132775824)

INFO:__main__:Average solver time per time step: 0.00418 s

INFO:__main__:{'P1': 1.3000000000000005, 'P2': -1.0, 'P3': -1.0, 'P4': -1.0, 'P5': -1.0, 'P6': -1.0, 'P7': -1.0, 'P8': -1.0, 'P9': -1.0}

INFO:beat.telemetry:PDE step timing step=380, t=(9.45000, 9.50000), ksp_iterations=5, ksp_residual_norm=6.391491e-03, ksp_converged_reason=2, ode_function_call=0.588649s, ode_state_update=0.002035s, ode_total_step=0.591750s, ode_step=0.652354s, ode_to_dolfin=0.001235s, ode_to_pde=0.000954s, pde_assign_previous_before=0.000488s, pde_set_time=0.000763s, pde_update_matrices=0.010126s, pde_update_rhs=0.021186s, pde_linear_solve=0.092312s, pde_scatter_forward=0.000316s, pde_total_step=0.126969s, pde_step=0.127452s, pde_to_ode=0.001190s, ode_from_dolfin=0.000385s, pde_assign_previous_after=0.000501s, total_step=0.792709s

INFO:beat.telemetry:PDE step timing step=400, t=(9.95000, 10.00000), ksp_iterations=5, ksp_residual_norm=6.436577e-03, ksp_converged_reason=2, ode_function_call=0.629468s, ode_state_update=0.002147s, ode_total_step=0.632779s, ode_step=0.696566s, ode_to_dolfin=0.001320s, ode_to_pde=0.001019s, pde_assign_previous_before=0.000532s, pde_set_time=0.000820s, pde_update_matrices=0.010126s, pde_update_rhs=0.023020s, pde_linear_solve=0.097577s, pde_scatter_forward=0.000341s, pde_total_step=0.134321s, pde_step=0.134845s, pde_to_ode=0.001277s, ode_from_dolfin=0.000409s, pde_assign_previous_after=0.000531s, total_step=0.845281s

INFO:__main__:Solve for t=10.00, v.max() =np.float64(32.027472850692284), v.min() =np.float64(-89.54034271038506)

INFO:__main__:Average solver time per time step: 0.00429 s

INFO:__main__:{'P1': 1.3000000000000005, 'P2': -1.0, 'P3': -1.0, 'P4': -1.0, 'P5': -1.0, 'P6': -1.0, 'P7': -1.0, 'P8': -1.0, 'P9': -1.0}

INFO:beat.telemetry:PDE step timing step=420, t=(10.45000, 10.50000), ksp_iterations=5, ksp_residual_norm=6.442955e-03, ksp_converged_reason=2, ode_function_call=0.668587s, ode_state_update=0.002260s, ode_total_step=0.672112s, ode_step=0.739091s, ode_to_dolfin=0.001416s, ode_to_pde=0.001088s, pde_assign_previous_before=0.000578s, pde_set_time=0.000876s, pde_update_matrices=0.010126s, pde_update_rhs=0.024911s, pde_linear_solve=0.102794s, pde_scatter_forward=0.000367s, pde_total_step=0.141685s, pde_step=0.142250s, pde_to_ode=0.001366s, ode_from_dolfin=0.000433s, pde_assign_previous_after=0.000561s, total_step=0.896199s

INFO:beat.telemetry:PDE step timing step=440, t=(10.95000, 11.00000), ksp_iterations=5, ksp_residual_norm=6.445296e-03, ksp_converged_reason=2, ode_function_call=0.725716s, ode_state_update=0.002371s, ode_total_step=0.729586s, ode_step=0.800507s, ode_to_dolfin=0.001500s, ode_to_pde=0.001312s, pde_assign_previous_before=0.000700s, pde_set_time=0.000964s, pde_update_matrices=0.010126s, pde_update_rhs=0.027921s, pde_linear_solve=0.110314s, pde_scatter_forward=0.000394s, pde_total_step=0.152523s, pde_step=0.153154s, pde_to_ode=0.001457s, ode_from_dolfin=0.000458s, pde_assign_previous_after=0.000593s, total_step=0.969872s

INFO:__main__:Solve for t=11.00, v.max() =np.float64(27.530393148255257), v.min() =np.float64(-89.60004871887763)

INFO:__main__:Average solver time per time step: 0.00447 s

INFO:__main__:{'P1': 1.3000000000000005, 'P2': -1.0, 'P3': -1.0, 'P4': -1.0, 'P5': -1.0, 'P6': -1.0, 'P7': -1.0, 'P8': -1.0, 'P9': -1.0}

INFO:beat.telemetry:PDE step timing step=460, t=(11.45000, 11.50000), ksp_iterations=5, ksp_residual_norm=6.445145e-03, ksp_converged_reason=2, ode_function_call=0.752215s, ode_state_update=0.002495s, ode_total_step=0.756284s, ode_step=0.830913s, ode_to_dolfin=0.001595s, ode_to_pde=0.001375s, pde_assign_previous_before=0.000737s, pde_set_time=0.001022s, pde_update_matrices=0.010126s, pde_update_rhs=0.029675s, pde_linear_solve=0.115344s, pde_scatter_forward=0.000419s, pde_total_step=0.159555s, pde_step=0.160225s, pde_to_ode=0.001534s, ode_from_dolfin=0.000481s, pde_assign_previous_after=0.000624s, total_step=1.008280s

INFO:beat.telemetry:PDE step timing step=480, t=(11.95000, 12.00000), ksp_iterations=5, ksp_residual_norm=6.431985e-03, ksp_converged_reason=2, ode_function_call=0.779029s, ode_state_update=0.002612s, ode_total_step=0.783310s, ode_step=0.861776s, ode_to_dolfin=0.001702s, ode_to_pde=0.001445s, pde_assign_previous_before=0.000779s, pde_set_time=0.001088s, pde_update_matrices=0.010126s, pde_update_rhs=0.031556s, pde_linear_solve=0.120470s, pde_scatter_forward=0.000446s, pde_total_step=0.166828s, pde_step=0.167542s, pde_to_ode=0.001625s, ode_from_dolfin=0.000505s, pde_assign_previous_after=0.000656s, total_step=1.047495s

INFO:__main__:Solve for t=12.00, v.max() =np.float64(29.2850063219413), v.min() =np.float64(-89.88407948626902)

INFO:__main__:Average solver time per time step: 0.00443 s

INFO:__main__:{'P1': 1.3000000000000005, 'P2': -1.0, 'P3': -1.0, 'P4': -1.0, 'P5': -1.0, 'P6': -1.0, 'P7': -1.0, 'P8': -1.0, 'P9': -1.0}

INFO:beat.telemetry:PDE step timing step=500, t=(12.45000, 12.50000), ksp_iterations=5, ksp_residual_norm=6.424940e-03, ksp_converged_reason=2, ode_function_call=0.817873s, ode_state_update=0.002723s, ode_total_step=0.822345s, ode_step=0.904088s, ode_to_dolfin=0.001785s, ode_to_pde=0.001505s, pde_assign_previous_before=0.000816s, pde_set_time=0.001140s, pde_update_matrices=0.010126s, pde_update_rhs=0.033205s, pde_linear_solve=0.125480s, pde_scatter_forward=0.000469s, pde_total_step=0.173726s, pde_step=0.174478s, pde_to_ode=0.001701s, ode_from_dolfin=0.000528s, pde_assign_previous_after=0.000684s, total_step=1.097611s

INFO:beat.telemetry:PDE step timing step=520, t=(12.95000, 13.00000), ksp_iterations=5, ksp_residual_norm=6.428446e-03, ksp_converged_reason=2, ode_function_call=0.858549s, ode_state_update=0.002832s, ode_total_step=0.863236s, ode_step=0.948207s, ode_to_dolfin=0.001871s, ode_to_pde=0.001571s, pde_assign_previous_before=0.000859s, pde_set_time=0.001198s, pde_update_matrices=0.010126s, pde_update_rhs=0.035100s, pde_linear_solve=0.130632s, pde_scatter_forward=0.000496s, pde_total_step=0.181040s, pde_step=0.181847s, pde_to_ode=0.001786s, ode_from_dolfin=0.000552s, pde_assign_previous_after=0.000713s, total_step=1.150059s

INFO:__main__:Solve for t=13.00, v.max() =np.float64(31.858037866789587), v.min() =np.float64(-89.89345659238829)

INFO:__main__:Average solver time per time step: 0.00449 s

INFO:__main__:{'P1': 1.3000000000000005, 'P2': -1.0, 'P3': -1.0, 'P4': -1.0, 'P5': -1.0, 'P6': -1.0, 'P7': -1.0, 'P8': -1.0, 'P9': -1.0}

INFO:beat.telemetry:PDE step timing step=540, t=(13.45000, 13.50000), ksp_iterations=5, ksp_residual_norm=6.421607e-03, ksp_converged_reason=2, ode_function_call=0.896686s, ode_state_update=0.002937s, ode_total_step=0.901538s, ode_step=0.989734s, ode_to_dolfin=0.001945s, ode_to_pde=0.001623s, pde_assign_previous_before=0.000889s, pde_set_time=0.001244s, pde_update_matrices=0.010126s, pde_update_rhs=0.036452s, pde_linear_solve=0.135559s, pde_scatter_forward=0.000515s, pde_total_step=0.187518s, pde_step=0.188354s, pde_to_ode=0.001857s, ode_from_dolfin=0.000572s, pde_assign_previous_after=0.000740s, total_step=1.198866s

INFO:beat.telemetry:PDE step timing step=560, t=(13.95000, 14.00000), ksp_iterations=5, ksp_residual_norm=6.424231e-03, ksp_converged_reason=2, ode_function_call=0.936181s, ode_state_update=0.003042s, ode_total_step=0.941191s, ode_step=1.032542s, ode_to_dolfin=0.001999s, ode_to_pde=0.001669s, pde_assign_previous_before=0.000913s, pde_set_time=0.001278s, pde_update_matrices=0.010126s, pde_update_rhs=0.037489s, pde_linear_solve=0.140346s, pde_scatter_forward=0.000529s, pde_total_step=0.193502s, pde_step=0.194361s, pde_to_ode=0.001919s, ode_from_dolfin=0.000591s, pde_assign_previous_after=0.000764s, total_step=1.248298s

INFO:__main__:Solve for t=14.00, v.max() =np.float64(35.6728277838826), v.min() =np.float64(-91.01086871035406)

INFO:__main__:Average solver time per time step: 0.00453 s

INFO:__main__:{'P1': 1.3000000000000005, 'P2': -1.0, 'P3': -1.0, 'P4': -1.0, 'P5': -1.0, 'P6': -1.0, 'P7': -1.0, 'P8': -1.0, 'P9': -1.0}

INFO:beat.telemetry:PDE step timing step=580, t=(14.45000, 14.50000), ksp_iterations=5, ksp_residual_norm=6.392665e-03, ksp_converged_reason=2, ode_function_call=0.973970s, ode_state_update=0.003146s, ode_total_step=0.979170s, ode_step=1.073675s, ode_to_dolfin=0.002075s, ode_to_pde=0.001723s, pde_assign_previous_before=0.000946s, pde_set_time=0.001325s, pde_update_matrices=0.010126s, pde_update_rhs=0.038939s, pde_linear_solve=0.145295s, pde_scatter_forward=0.000551s, pde_total_step=0.200114s, pde_step=0.201008s, pde_to_ode=0.001992s, ode_from_dolfin=0.000613s, pde_assign_previous_after=0.000792s, total_step=1.296879s

INFO:beat.telemetry:PDE step timing step=600, t=(14.95000, 15.00000), ksp_iterations=5, ksp_residual_norm=6.353125e-03, ksp_converged_reason=2, ode_function_call=1.012482s, ode_state_update=0.003250s, ode_total_step=1.017836s, ode_step=1.115405s, ode_to_dolfin=0.002123s, ode_to_pde=0.001765s, pde_assign_previous_before=0.000968s, pde_set_time=0.001357s, pde_update_matrices=0.010126s, pde_update_rhs=0.040000s, pde_linear_solve=0.149966s, pde_scatter_forward=0.000565s, pde_total_step=0.205996s, pde_step=0.206913s, pde_to_ode=0.002047s, ode_from_dolfin=0.000631s, pde_assign_previous_after=0.000816s, total_step=1.345097s

INFO:__main__:Solve for t=15.00, v.max() =np.float64(35.87974315808209), v.min() =np.float64(-91.14592772494588)

INFO:__main__:Average solver time per time step: 0.00455 s

INFO:__main__:{'P1': 1.3000000000000005, 'P2': -1.0, 'P3': -1.0, 'P4': -1.0, 'P5': 14.150000000000066, 'P6': -1.0, 'P7': -1.0, 'P8': -1.0, 'P9': -1.0}

INFO:beat.telemetry:PDE step timing step=620, t=(15.45000, 15.50000), ksp_iterations=5, ksp_residual_norm=6.313597e-03, ksp_converged_reason=2, ode_function_call=1.037584s, ode_state_update=0.003355s, ode_total_step=1.043095s, ode_step=1.143855s, ode_to_dolfin=0.002216s, ode_to_pde=0.001814s, pde_assign_previous_before=0.000994s, pde_set_time=0.001393s, pde_update_matrices=0.010126s, pde_update_rhs=0.040975s, pde_linear_solve=0.154689s, pde_scatter_forward=0.000580s, pde_total_step=0.211851s, pde_step=0.212791s, pde_to_ode=0.002104s, ode_from_dolfin=0.000650s, pde_assign_previous_after=0.000840s, total_step=1.380074s

INFO:beat.telemetry:PDE step timing step=640, t=(15.95000, 16.00000), ksp_iterations=5, ksp_residual_norm=6.281626e-03, ksp_converged_reason=2, ode_function_call=1.062360s, ode_state_update=0.003458s, ode_total_step=1.068029s, ode_step=1.171949s, ode_to_dolfin=0.002268s, ode_to_pde=0.001856s, pde_assign_previous_before=0.001016s, pde_set_time=0.001425s, pde_update_matrices=0.010126s, pde_update_rhs=0.041779s, pde_linear_solve=0.159359s, pde_scatter_forward=0.000592s, pde_total_step=0.217464s, pde_step=0.218424s, pde_to_ode=0.002159s, ode_from_dolfin=0.000668s, pde_assign_previous_after=0.000864s, total_step=1.414355s

INFO:__main__:Solve for t=16.00, v.max() =np.float64(32.71510671586424), v.min() =np.float64(-91.35185801864718)

INFO:__main__:Average solver time per time step: 0.00449 s

INFO:__main__:{'P1': 1.3000000000000005, 'P2': -1.0, 'P3': -1.0, 'P4': -1.0, 'P5': 14.150000000000066, 'P6': -1.0, 'P7': -1.0, 'P8': -1.0, 'P9': -1.0}

INFO:beat.telemetry:PDE step timing step=660, t=(16.45000, 16.50000), ksp_iterations=5, ksp_residual_norm=6.248057e-03, ksp_converged_reason=2, ode_function_call=1.099941s, ode_state_update=0.003569s, ode_total_step=1.105787s, ode_step=1.212977s, ode_to_dolfin=0.002332s, ode_to_pde=0.001901s, pde_assign_previous_before=0.001039s, pde_set_time=0.001460s, pde_update_matrices=0.010126s, pde_update_rhs=0.042717s, pde_linear_solve=0.164088s, pde_scatter_forward=0.000606s, pde_total_step=0.223285s, pde_step=0.224268s, pde_to_ode=0.002216s, ode_from_dolfin=0.000686s, pde_assign_previous_after=0.000889s, total_step=1.461836s

INFO:beat.telemetry:PDE step timing step=680, t=(16.95000, 17.00000), ksp_iterations=5, ksp_residual_norm=6.222907e-03, ksp_converged_reason=2, ode_function_call=1.138832s, ode_state_update=0.003676s, ode_total_step=1.144839s, ode_step=1.255280s, ode_to_dolfin=0.002386s, ode_to_pde=0.001955s, pde_assign_previous_before=0.001060s, pde_set_time=0.001494s, pde_update_matrices=0.010126s, pde_update_rhs=0.043589s, pde_linear_solve=0.168792s, pde_scatter_forward=0.000618s, pde_total_step=0.229009s, pde_step=0.230014s, pde_to_ode=0.002290s, ode_from_dolfin=0.000705s, pde_assign_previous_after=0.000914s, total_step=1.510488s

INFO:__main__:Solve for t=17.00, v.max() =np.float64(31.562363397879523), v.min() =np.float64(-92.23346573815942)

INFO:__main__:Average solver time per time step: 0.00451 s

INFO:__main__:{'P1': 1.3000000000000005, 'P2': -1.0, 'P3': -1.0, 'P4': -1.0, 'P5': 14.150000000000066, 'P6': -1.0, 'P7': -1.0, 'P8': -1.0, 'P9': -1.0}

INFO:beat.telemetry:PDE step timing step=700, t=(17.45000, 17.50000), ksp_iterations=5, ksp_residual_norm=6.223496e-03, ksp_converged_reason=2, ode_function_call=1.176577s, ode_state_update=0.003782s, ode_total_step=1.182744s, ode_step=1.296404s, ode_to_dolfin=0.002454s, ode_to_pde=0.002002s, pde_assign_previous_before=0.001084s, pde_set_time=0.001531s, pde_update_matrices=0.010126s, pde_update_rhs=0.044599s, pde_linear_solve=0.173520s, pde_scatter_forward=0.000633s, pde_total_step=0.234905s, pde_step=0.235933s, pde_to_ode=0.002346s, ode_from_dolfin=0.000723s, pde_assign_previous_after=0.000939s, total_step=1.558156s

INFO:beat.telemetry:PDE step timing step=720, t=(17.95000, 18.00000), ksp_iterations=5, ksp_residual_norm=6.233631e-03, ksp_converged_reason=2, ode_function_call=1.215815s, ode_state_update=0.003888s, ode_total_step=1.222141s, ode_step=1.338979s, ode_to_dolfin=0.002508s, ode_to_pde=0.002046s, pde_assign_previous_before=0.001106s, pde_set_time=0.001564s, pde_update_matrices=0.010126s, pde_update_rhs=0.045446s, pde_linear_solve=0.178324s, pde_scatter_forward=0.000646s, pde_total_step=0.240704s, pde_step=0.241753s, pde_to_ode=0.002404s, ode_from_dolfin=0.000742s, pde_assign_previous_after=0.000965s, total_step=1.607152s

INFO:__main__:Solve for t=18.00, v.max() =np.float64(33.3729493224822), v.min() =np.float64(-92.47996298199891)

INFO:__main__:Average solver time per time step: 0.00453 s

INFO:__main__:{'P1': 1.3000000000000005, 'P2': -1.0, 'P3': -1.0, 'P4': -1.0, 'P5': 14.150000000000066, 'P6': -1.0, 'P7': -1.0, 'P8': -1.0, 'P9': -1.0}

INFO:beat.telemetry:PDE step timing step=740, t=(18.45000, 18.50000), ksp_iterations=5, ksp_residual_norm=6.205766e-03, ksp_converged_reason=2, ode_function_call=1.252629s, ode_state_update=0.004124s, ode_total_step=1.259243s, ode_step=1.379424s, ode_to_dolfin=0.002573s, ode_to_pde=0.002094s, pde_assign_previous_before=0.001130s, pde_set_time=0.001601s, pde_update_matrices=0.010126s, pde_update_rhs=0.046661s, pde_linear_solve=0.183439s, pde_scatter_forward=0.000661s, pde_total_step=0.247194s, pde_step=0.248268s, pde_to_ode=0.002462s, ode_from_dolfin=0.000761s, pde_assign_previous_after=0.000990s, total_step=1.654746s

INFO:beat.telemetry:PDE step timing step=760, t=(18.95000, 19.00000), ksp_iterations=5, ksp_residual_norm=6.164973e-03, ksp_converged_reason=2, ode_function_call=1.290917s, ode_state_update=0.004228s, ode_total_step=1.297702s, ode_step=1.420983s, ode_to_dolfin=0.002632s, ode_to_pde=0.002143s, pde_assign_previous_before=0.001158s, pde_set_time=0.001640s, pde_update_matrices=0.010126s, pde_update_rhs=0.047814s, pde_linear_solve=0.188202s, pde_scatter_forward=0.000678s, pde_total_step=0.253291s, pde_step=0.254391s, pde_to_ode=0.002525s, ode_from_dolfin=0.000781s, pde_assign_previous_after=0.001016s, total_step=1.703097s

INFO:__main__:Solve for t=19.00, v.max() =np.float64(32.413759220069835), v.min() =np.float64(-92.42624050740888)

INFO:__main__:Average solver time per time step: 0.00455 s

INFO:__main__:{'P1': 1.3000000000000005, 'P2': -1.0, 'P3': -1.0, 'P4': -1.0, 'P5': 14.150000000000066, 'P6': -1.0, 'P7': -1.0, 'P8': -1.0, 'P9': -1.0}

INFO:beat.telemetry:PDE step timing step=780, t=(19.45000, 19.50000), ksp_iterations=5, ksp_residual_norm=6.134901e-03, ksp_converged_reason=2, ode_function_call=1.317019s, ode_state_update=0.004339s, ode_total_step=1.323991s, ode_step=1.450928s, ode_to_dolfin=0.002719s, ode_to_pde=0.002204s, pde_assign_previous_before=0.001193s, pde_set_time=0.001696s, pde_update_matrices=0.010126s, pde_update_rhs=0.049462s, pde_linear_solve=0.193334s, pde_scatter_forward=0.000700s, pde_total_step=0.260301s, pde_step=0.261439s, pde_to_ode=0.002603s, ode_from_dolfin=0.000803s, pde_assign_previous_after=0.001044s, total_step=1.740966s

INFO:beat.telemetry:PDE step timing step=800, t=(19.95000, 20.00000), ksp_iterations=5, ksp_residual_norm=6.098557e-03, ksp_converged_reason=2, ode_function_call=1.341481s, ode_state_update=0.004444s, ode_total_step=1.348606s, ode_step=1.478657s, ode_to_dolfin=0.002769s, ode_to_pde=0.002247s, pde_assign_previous_before=0.001215s, pde_set_time=0.001728s, pde_update_matrices=0.010126s, pde_update_rhs=0.050271s, pde_linear_solve=0.198049s, pde_scatter_forward=0.000711s, pde_total_step=0.265966s, pde_step=0.267137s, pde_to_ode=0.002659s, ode_from_dolfin=0.000821s, pde_assign_previous_after=0.001068s, total_step=1.774947s

INFO:__main__:Solve for t=20.00, v.max() =np.float64(31.647177836306682), v.min() =np.float64(-92.50286199984987)

INFO:__main__:Average solver time per time step: 0.00450 s

INFO:__main__:{'P1': 1.3000000000000005, 'P2': -1.0, 'P3': -1.0, 'P4': -1.0, 'P5': 14.150000000000066, 'P6': -1.0, 'P7': -1.0, 'P8': -1.0, 'P9': -1.0}

INFO:beat.telemetry:
==================================================
               PERFORMANCE SUMMARY                
==================================================
Total Steps:           802
KSP Total Iterations:  2005
KSP Max Iterations:    5
--------------------------------------------------
Metric                              |   Time (s)
--------------------------------------------------
total_step                          |     1.7788
ode_step                            |     1.4817
ode_total_step                      |     1.3513
ode_function_call                   |     1.3441
pde_step                            |     0.2679
pde_total_step                      |     0.2667
pde_linear_solve                    |     0.1986
pde_update_rhs                      |     0.0504
pde_update_matrices                 |     0.0101
ode_state_update                    |     0.0045
ode_to_dolfin                       |     0.0028
pde_to_ode                          |     0.0027
ode_to_pde                          |     0.0023
pde_set_time                        |     0.0017
pde_assign_previous_before          |     0.0012
pde_assign_previous_after           |     0.0011
ode_from_dolfin                     |     0.0008
pde_scatter_forward                 |     0.0007
==================================================

INFO:beat.telemetry:Performance summary saved to results-niederer-benchmark/performance_summary.json

_

# 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
[FK98]

Flavio Fenton and Alain Karma. Vortex dynamics in three-dimensional continuous myocardium with fiber rotation: filament instability and fibrillation. Chaos: An Interdisciplinary Journal of Nonlinear Science, 8(1):20–47, 1998. doi:10.1063/1.166311.

[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.

[SLC+07]

Joakim Sundnes, Glenn Terje Lines, Xing Cai, Bjørn Frederik Nielsen, Kent-Andre Mardal, and Aslak Tveito. Computing the electrical activity in the heart. Volume 1. Springer Science & Business Media, 2007.

[tTNNP04]

Kirsten H W J ten Tusscher, Denis Noble, Peter J Noble, and Alexander V Panfilov. A model for human ventricular tissue. American Journal of Physiology-Heart and Circulatory Physiology, 286(4):H1573–H1589, 2004. doi:10.1152/ajpheart.00794.2003.

[tTP06]

Kirsten H W J ten Tusscher and Alexander V Panfilov. Alternans and spiral breakup in a human ventricular tissue model. American Journal of Physiology-Heart and Circulatory Physiology, 291(3):H1088–H1100, 2006. doi:10.1152/ajpheart.00109.2006.

[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.