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.065107s, ode_state_update=0.000319s, ode_total_step=0.065654s, ode_step=0.069348s, ode_to_dolfin=0.000274s, ode_to_pde=0.000169s, pde_assign_previous_before=0.000088s, pde_set_time=0.000118s, pde_update_matrices=0.009355s, pde_update_rhs=0.025121s, pde_linear_solve=0.011240s, pde_scatter_forward=0.000046s, pde_total_step=0.046275s, pde_step=0.046340s, pde_to_ode=0.000200s, ode_from_dolfin=0.000059s, pde_assign_previous_after=0.000076s, total_step=0.117593s

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.130046s, ode_state_update=0.000661s, ode_total_step=0.131146s, ode_step=0.139013s, ode_to_dolfin=0.000540s, ode_to_pde=0.000307s, pde_assign_previous_before=0.000169s, pde_set_time=0.000235s, pde_update_matrices=0.009355s, pde_update_rhs=0.050438s, pde_linear_solve=0.020600s, pde_scatter_forward=0.000092s, pde_total_step=0.081395s, pde_step=0.081519s, pde_to_ode=0.000393s, ode_from_dolfin=0.000128s, pde_assign_previous_after=0.000152s, total_step=0.224297s

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

INFO:__main__:Average solver time per time step: 0.01133 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.194780s, ode_state_update=0.000989s, ode_total_step=0.196404s, ode_step=0.208470s, ode_to_dolfin=0.000782s, ode_to_pde=0.000445s, pde_assign_previous_before=0.000275s, pde_set_time=0.000344s, pde_update_matrices=0.009355s, pde_update_rhs=0.075668s, pde_linear_solve=0.029932s, pde_scatter_forward=0.000133s, pde_total_step=0.116398s, pde_step=0.116576s, pde_to_ode=0.000584s, ode_from_dolfin=0.000184s, pde_assign_previous_after=0.000226s, total_step=0.330647s

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.259702s, ode_state_update=0.001324s, ode_total_step=0.261860s, ode_step=0.278189s, ode_to_dolfin=0.001025s, ode_to_pde=0.000580s, pde_assign_previous_before=0.000357s, pde_set_time=0.000459s, pde_update_matrices=0.009355s, pde_update_rhs=0.100891s, pde_linear_solve=0.039311s, pde_scatter_forward=0.000174s, pde_total_step=0.151432s, pde_step=0.151665s, pde_to_ode=0.000787s, ode_from_dolfin=0.000242s, pde_assign_previous_after=0.000301s, total_step=0.437237s

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

INFO:__main__:Average solver time per time step: 0.01104 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.324896s, ode_state_update=0.001668s, ode_total_step=0.327586s, ode_step=0.348252s, ode_to_dolfin=0.001283s, ode_to_pde=0.000718s, pde_assign_previous_before=0.000439s, pde_set_time=0.000578s, pde_update_matrices=0.009355s, pde_update_rhs=0.126419s, pde_linear_solve=0.048761s, pde_scatter_forward=0.000220s, pde_total_step=0.186871s, pde_step=0.187189s, pde_to_ode=0.000983s, ode_from_dolfin=0.000297s, pde_assign_previous_after=0.000379s, total_step=0.544686s

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.391213s, ode_state_update=0.002051s, ode_total_step=0.394481s, ode_step=0.419472s, ode_to_dolfin=0.001528s, ode_to_pde=0.000855s, pde_assign_previous_before=0.000520s, pde_set_time=0.000699s, pde_update_matrices=0.009355s, pde_update_rhs=0.152106s, pde_linear_solve=0.058276s, pde_scatter_forward=0.000268s, pde_total_step=0.222544s, pde_step=0.222926s, pde_to_ode=0.001198s, ode_from_dolfin=0.000356s, pde_assign_previous_after=0.000487s, total_step=0.653494s

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

INFO:__main__:Average solver time per time step: 0.01100 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.454758s, ode_state_update=0.002334s, ode_total_step=0.458480s, ode_step=0.487707s, ode_to_dolfin=0.001742s, ode_to_pde=0.001007s, pde_assign_previous_before=0.000601s, pde_set_time=0.000799s, pde_update_matrices=0.009355s, pde_update_rhs=0.177260s, pde_linear_solve=0.067575s, pde_scatter_forward=0.000305s, pde_total_step=0.257402s, pde_step=0.257837s, pde_to_ode=0.001364s, ode_from_dolfin=0.000418s, pde_assign_previous_after=0.000555s, total_step=0.758326s

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.519861s, ode_state_update=0.002700s, ode_total_step=0.524150s, ode_step=0.557757s, ode_to_dolfin=0.001969s, ode_to_pde=0.001144s, pde_assign_previous_before=0.000684s, pde_set_time=0.000908s, pde_update_matrices=0.009355s, pde_update_rhs=0.202842s, pde_linear_solve=0.077058s, pde_scatter_forward=0.000351s, pde_total_step=0.292914s, pde_step=0.293410s, pde_to_ode=0.001572s, ode_from_dolfin=0.000473s, pde_assign_previous_after=0.000635s, total_step=0.865760s

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

INFO:__main__:Average solver time per time step: 0.01093 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.583440s, ode_state_update=0.002958s, ode_total_step=0.588144s, ode_step=0.625968s, ode_to_dolfin=0.002184s, ode_to_pde=0.001294s, pde_assign_previous_before=0.000764s, pde_set_time=0.001058s, pde_update_matrices=0.009355s, pde_update_rhs=0.227915s, pde_linear_solve=0.086417s, pde_scatter_forward=0.000387s, pde_total_step=0.327819s, pde_step=0.328371s, pde_to_ode=0.001762s, ode_from_dolfin=0.000558s, pde_assign_previous_after=0.000706s, total_step=0.970677s

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.647625s, ode_state_update=0.003251s, ode_total_step=0.652788s, ode_step=0.694854s, ode_to_dolfin=0.002418s, ode_to_pde=0.001427s, pde_assign_previous_before=0.000846s, pde_set_time=0.001163s, pde_update_matrices=0.009355s, pde_update_rhs=0.253413s, pde_linear_solve=0.095776s, pde_scatter_forward=0.000423s, pde_total_step=0.363087s, pde_step=0.363694s, pde_to_ode=0.001941s, ode_from_dolfin=0.000613s, pde_assign_previous_after=0.000783s, total_step=1.076582s

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

INFO:__main__:Average solver time per time step: 0.01087 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.711194s, ode_state_update=0.003536s, ode_total_step=0.716802s, ode_step=0.763135s, ode_to_dolfin=0.002639s, ode_to_pde=0.001557s, pde_assign_previous_before=0.000932s, pde_set_time=0.001264s, pde_update_matrices=0.009355s, pde_update_rhs=0.278518s, pde_linear_solve=0.105056s, pde_scatter_forward=0.000459s, pde_total_step=0.397874s, pde_step=0.398535s, pde_to_ode=0.002122s, ode_from_dolfin=0.000665s, pde_assign_previous_after=0.000853s, total_step=1.181390s

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.775450s, ode_state_update=0.003837s, ode_total_step=0.781540s, ode_step=0.832185s, ode_to_dolfin=0.002869s, ode_to_pde=0.001692s, pde_assign_previous_before=0.001017s, pde_set_time=0.001374s, pde_update_matrices=0.009355s, pde_update_rhs=0.303831s, pde_linear_solve=0.114371s, pde_scatter_forward=0.000510s, pde_total_step=0.432942s, pde_step=0.433698s, pde_to_ode=0.002304s, ode_from_dolfin=0.000723s, pde_assign_previous_after=0.000933s, total_step=1.287339s

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

INFO:__main__:Average solver time per time step: 0.01083 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.839964s, ode_state_update=0.004134s, ode_total_step=0.846524s, ode_step=0.901385s, ode_to_dolfin=0.003103s, ode_to_pde=0.001829s, pde_assign_previous_before=0.001098s, pde_set_time=0.001483s, pde_update_matrices=0.009355s, pde_update_rhs=0.328942s, pde_linear_solve=0.123743s, pde_scatter_forward=0.000550s, pde_total_step=0.467858s, pde_step=0.468671s, pde_to_ode=0.002485s, ode_from_dolfin=0.000776s, pde_assign_previous_after=0.001002s, total_step=1.393239s

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.908803s, ode_state_update=0.004518s, ode_total_step=0.915943s, ode_step=0.975189s, ode_to_dolfin=0.003357s, ode_to_pde=0.001973s, pde_assign_previous_before=0.001182s, pde_set_time=0.001608s, pde_update_matrices=0.009355s, pde_update_rhs=0.354500s, pde_linear_solve=0.133281s, pde_scatter_forward=0.000621s, pde_total_step=0.503447s, pde_step=0.504323s, pde_to_ode=0.002704s, ode_from_dolfin=0.000972s, pde_assign_previous_after=0.001089s, total_step=1.504682s

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

INFO:__main__:Average solver time per time step: 0.01085 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.973757s, ode_state_update=0.004838s, ode_total_step=0.981403s, ode_step=1.045018s, ode_to_dolfin=0.003602s, ode_to_pde=0.002111s, pde_assign_previous_before=0.001264s, pde_set_time=0.001728s, pde_update_matrices=0.009355s, pde_update_rhs=0.379787s, pde_linear_solve=0.142632s, pde_scatter_forward=0.000662s, pde_total_step=0.538543s, pde_step=0.539477s, pde_to_ode=0.002900s, ode_from_dolfin=0.001026s, pde_assign_previous_after=0.001166s, total_step=1.611406s

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=1.038269s, ode_state_update=0.005167s, ode_total_step=1.046448s, ode_step=1.114426s, ode_to_dolfin=0.003840s, ode_to_pde=0.002250s, pde_assign_previous_before=0.001345s, pde_set_time=0.001878s, pde_update_matrices=0.009355s, pde_update_rhs=0.405422s, pde_linear_solve=0.151990s, pde_scatter_forward=0.000701s, pde_total_step=0.574047s, pde_step=0.575037s, pde_to_ode=0.003074s, ode_from_dolfin=0.001078s, pde_assign_previous_after=0.001237s, total_step=1.718081s

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

INFO:__main__:Average solver time per time step: 0.01084 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=1.102370s, ode_state_update=0.005447s, ode_total_step=1.111012s, ode_step=1.183370s, ode_to_dolfin=0.004071s, ode_to_pde=0.002383s, pde_assign_previous_before=0.001426s, pde_set_time=0.001987s, pde_update_matrices=0.009355s, pde_update_rhs=0.430638s, pde_linear_solve=0.161339s, pde_scatter_forward=0.000739s, pde_total_step=0.609052s, pde_step=0.610096s, pde_to_ode=0.003256s, ode_from_dolfin=0.001169s, pde_assign_previous_after=0.001308s, total_step=1.823815s

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=1.167097s, ode_state_update=0.005785s, ode_total_step=1.176281s, ode_step=1.253085s, ode_to_dolfin=0.004334s, ode_to_pde=0.002519s, pde_assign_previous_before=0.001509s, pde_set_time=0.002098s, pde_update_matrices=0.009355s, pde_update_rhs=0.456000s, pde_linear_solve=0.170698s, pde_scatter_forward=0.000781s, pde_total_step=0.644200s, pde_step=0.645302s, pde_to_ode=0.003455s, ode_from_dolfin=0.001230s, pde_assign_previous_after=0.001385s, total_step=1.930585s

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

INFO:__main__:Average solver time per time step: 0.01083 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=1.231556s, ode_state_update=0.006095s, ode_total_step=1.241231s, ode_step=1.322506s, ode_to_dolfin=0.004594s, ode_to_pde=0.002660s, pde_assign_previous_before=0.001590s, pde_set_time=0.002211s, pde_update_matrices=0.009355s, pde_update_rhs=0.481445s, pde_linear_solve=0.180014s, pde_scatter_forward=0.000821s, pde_total_step=0.679409s, pde_step=0.680568s, pde_to_ode=0.003647s, ode_from_dolfin=0.001288s, pde_assign_previous_after=0.001462s, total_step=2.037114s

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=1.296447s, ode_state_update=0.006419s, ode_total_step=1.306645s, ode_step=1.392469s, ode_to_dolfin=0.004853s, ode_to_pde=0.002797s, pde_assign_previous_before=0.001671s, pde_set_time=0.002325s, pde_update_matrices=0.009355s, pde_update_rhs=0.506902s, pde_linear_solve=0.189347s, pde_scatter_forward=0.000861s, pde_total_step=0.714636s, pde_step=0.715852s, pde_to_ode=0.003838s, ode_from_dolfin=0.001345s, pde_assign_previous_after=0.001538s, total_step=2.144176s

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

INFO:__main__:Average solver time per time step: 0.01083 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=1.360141s, ode_state_update=0.006695s, ode_total_step=1.370774s, ode_step=1.461204s, ode_to_dolfin=0.005070s, ode_to_pde=0.002953s, pde_assign_previous_before=0.001754s, pde_set_time=0.002427s, pde_update_matrices=0.009355s, pde_update_rhs=0.531909s, pde_linear_solve=0.198705s, pde_scatter_forward=0.000899s, pde_total_step=0.749433s, pde_step=0.750702s, pde_to_ode=0.004015s, ode_from_dolfin=0.001396s, pde_assign_previous_after=0.001623s, total_step=2.249467s

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=1.424056s, ode_state_update=0.006968s, ode_total_step=1.435123s, ode_step=1.530321s, ode_to_dolfin=0.005283s, ode_to_pde=0.003088s, pde_assign_previous_before=0.001836s, pde_set_time=0.002532s, pde_update_matrices=0.009355s, pde_update_rhs=0.557100s, pde_linear_solve=0.208052s, pde_scatter_forward=0.000936s, pde_total_step=0.784383s, pde_step=0.785708s, pde_to_ode=0.004187s, ode_from_dolfin=0.001450s, pde_assign_previous_after=0.001694s, total_step=2.355279s

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

INFO:__main__:Average solver time per time step: 0.01081 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=1.488395s, ode_state_update=0.007236s, ode_total_step=1.499895s, ode_step=1.599835s, ode_to_dolfin=0.005596s, ode_to_pde=0.003223s, pde_assign_previous_before=0.001919s, pde_set_time=0.002635s, pde_update_matrices=0.009355s, pde_update_rhs=0.582829s, pde_linear_solve=0.217425s, pde_scatter_forward=0.000974s, pde_total_step=0.819902s, pde_step=0.821304s, pde_to_ode=0.004365s, ode_from_dolfin=0.001507s, pde_assign_previous_after=0.001785s, total_step=2.462176s

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=1.553003s, ode_state_update=0.007529s, ode_total_step=1.564972s, ode_step=1.669652s, ode_to_dolfin=0.005833s, ode_to_pde=0.003358s, pde_assign_previous_before=0.002000s, pde_set_time=0.002741s, pde_update_matrices=0.009355s, pde_update_rhs=0.608130s, pde_linear_solve=0.226836s, pde_scatter_forward=0.001013s, pde_total_step=0.855038s, pde_step=0.856498s, pde_to_ode=0.004546s, ode_from_dolfin=0.001565s, pde_assign_previous_after=0.001879s, total_step=2.568922s

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

INFO:__main__:Average solver time per time step: 0.01081 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=1.617337s, ode_state_update=0.007840s, ode_total_step=1.629802s, ode_step=1.739177s, ode_to_dolfin=0.006085s, ode_to_pde=0.003495s, pde_assign_previous_before=0.002106s, pde_set_time=0.002854s, pde_update_matrices=0.009355s, pde_update_rhs=0.633570s, pde_linear_solve=0.236301s, pde_scatter_forward=0.001053s, pde_total_step=0.890382s, pde_step=0.891898s, pde_to_ode=0.004747s, ode_from_dolfin=0.001655s, pde_assign_previous_after=0.001953s, total_step=2.675691s

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=1.682050s, ode_state_update=0.008145s, ode_total_step=1.695011s, ode_step=1.809101s, ode_to_dolfin=0.006333s, ode_to_pde=0.003631s, pde_assign_previous_before=0.002187s, pde_set_time=0.002967s, pde_update_matrices=0.009355s, pde_update_rhs=0.658968s, pde_linear_solve=0.245656s, pde_scatter_forward=0.001095s, pde_total_step=0.925576s, pde_step=0.927150s, pde_to_ode=0.004956s, ode_from_dolfin=0.001715s, pde_assign_previous_after=0.002032s, total_step=2.782700s

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

INFO:__main__:Average solver time per time step: 0.01081 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=1.746559s, ode_state_update=0.008463s, ode_total_step=1.760031s, ode_step=1.878850s, ode_to_dolfin=0.006587s, ode_to_pde=0.003770s, pde_assign_previous_before=0.002270s, pde_set_time=0.003078s, pde_update_matrices=0.009355s, pde_update_rhs=0.684202s, pde_linear_solve=0.255038s, pde_scatter_forward=0.001150s, pde_total_step=0.960637s, pde_step=0.962270s, pde_to_ode=0.005151s, ode_from_dolfin=0.001775s, pde_assign_previous_after=0.002111s, total_step=2.889376s

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=1.811562s, ode_state_update=0.008794s, ode_total_step=1.825566s, ode_step=1.949084s, ode_to_dolfin=0.006902s, ode_to_pde=0.003907s, pde_assign_previous_before=0.002354s, pde_set_time=0.003194s, pde_update_matrices=0.009355s, pde_update_rhs=0.709694s, pde_linear_solve=0.264445s, pde_scatter_forward=0.001193s, pde_total_step=0.996021s, pde_step=0.997712s, pde_to_ode=0.005355s, ode_from_dolfin=0.001841s, pde_assign_previous_after=0.002194s, total_step=2.997018s

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

INFO:__main__:Average solver time per time step: 0.01081 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=1.876477s, ode_state_update=0.009089s, ode_total_step=1.890980s, ode_step=2.019268s, ode_to_dolfin=0.007161s, ode_to_pde=0.004048s, pde_assign_previous_before=0.002437s, pde_set_time=0.003338s, pde_update_matrices=0.009355s, pde_update_rhs=0.735103s, pde_linear_solve=0.273822s, pde_scatter_forward=0.001234s, pde_total_step=1.031285s, pde_step=1.033034s, pde_to_ode=0.005566s, ode_from_dolfin=0.001915s, pde_assign_previous_after=0.002271s, total_step=3.104376s

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.942907s, ode_state_update=0.009469s, ode_total_step=1.958004s, ode_step=2.090969s, ode_to_dolfin=0.007436s, ode_to_pde=0.004187s, pde_assign_previous_before=0.002521s, pde_set_time=0.003486s, pde_update_matrices=0.009355s, pde_update_rhs=0.760593s, pde_linear_solve=0.283299s, pde_scatter_forward=0.001282s, pde_total_step=1.066765s, pde_step=1.068575s, pde_to_ode=0.005785s, ode_from_dolfin=0.001972s, pde_assign_previous_after=0.002348s, total_step=3.213492s

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

INFO:__main__:Average solver time per time step: 0.01082 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=2.009045s, ode_state_update=0.009830s, ode_total_step=2.024721s, ode_step=2.162382s, ode_to_dolfin=0.007740s, ode_to_pde=0.004351s, pde_assign_previous_before=0.002603s, pde_set_time=0.003606s, pde_update_matrices=0.009355s, pde_update_rhs=0.785980s, pde_linear_solve=0.292694s, pde_scatter_forward=0.001325s, pde_total_step=1.101998s, pde_step=1.103870s, pde_to_ode=0.005982s, ode_from_dolfin=0.002029s, pde_assign_previous_after=0.002436s, total_step=3.322124s

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=2.074247s, ode_state_update=0.010176s, ode_total_step=2.090508s, ode_step=2.232866s, ode_to_dolfin=0.008042s, ode_to_pde=0.004493s, pde_assign_previous_before=0.002685s, pde_set_time=0.003725s, pde_update_matrices=0.009355s, pde_update_rhs=0.811276s, pde_linear_solve=0.302067s, pde_scatter_forward=0.001368s, pde_total_step=1.137122s, pde_step=1.139057s, pde_to_ode=0.006179s, ode_from_dolfin=0.002094s, pde_assign_previous_after=0.002515s, total_step=3.429672s

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

INFO:__main__:Average solver time per time step: 0.01082 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=2.140663s, ode_state_update=0.010503s, ode_total_step=2.157456s, ode_step=2.304553s, ode_to_dolfin=0.008321s, ode_to_pde=0.004633s, pde_assign_previous_before=0.002768s, pde_set_time=0.003842s, pde_update_matrices=0.009355s, pde_update_rhs=0.836884s, pde_linear_solve=0.312065s, pde_scatter_forward=0.001410s, pde_total_step=1.173177s, pde_step=1.175173s, pde_to_ode=0.006382s, ode_from_dolfin=0.002155s, pde_assign_previous_after=0.002595s, total_step=3.539322s

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=2.205901s, ode_state_update=0.010819s, ode_total_step=2.223223s, ode_step=2.375020s, ode_to_dolfin=0.008580s, ode_to_pde=0.004775s, pde_assign_previous_before=0.002850s, pde_set_time=0.003983s, pde_update_matrices=0.009355s, pde_update_rhs=0.862206s, pde_linear_solve=0.321420s, pde_scatter_forward=0.001450s, pde_total_step=1.208320s, pde_step=1.210386s, pde_to_ode=0.006578s, ode_from_dolfin=0.002218s, pde_assign_previous_after=0.002677s, total_step=3.646826s

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

INFO:__main__:Average solver time per time step: 0.01083 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=2.272769s, ode_state_update=0.011244s, ode_total_step=2.290741s, ode_step=2.447346s, ode_to_dolfin=0.008840s, ode_to_pde=0.004915s, pde_assign_previous_before=0.002933s, pde_set_time=0.004100s, pde_update_matrices=0.009355s, pde_update_rhs=0.887956s, pde_linear_solve=0.331125s, pde_scatter_forward=0.001502s, pde_total_step=1.244264s, pde_step=1.246390s, pde_to_ode=0.006795s, ode_from_dolfin=0.002286s, pde_assign_previous_after=0.002759s, total_step=3.757062s

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=2.338191s, ode_state_update=0.011596s, ode_total_step=2.356709s, ode_step=2.518175s, ode_to_dolfin=0.009103s, ode_to_pde=0.005055s, pde_assign_previous_before=0.003016s, pde_set_time=0.004216s, pde_update_matrices=0.009355s, pde_update_rhs=0.913407s, pde_linear_solve=0.340556s, pde_scatter_forward=0.001547s, pde_total_step=1.279599s, pde_step=1.281784s, pde_to_ode=0.006996s, ode_from_dolfin=0.002345s, pde_assign_previous_after=0.002849s, total_step=3.865146s

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

INFO:__main__:Average solver time per time step: 0.01084 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=2.403238s, ode_state_update=0.011939s, ode_total_step=2.422288s, ode_step=2.588500s, ode_to_dolfin=0.009349s, ode_to_pde=0.005192s, pde_assign_previous_before=0.003095s, pde_set_time=0.004329s, pde_update_matrices=0.009355s, pde_update_rhs=0.938773s, pde_linear_solve=0.349938s, pde_scatter_forward=0.001591s, pde_total_step=1.314792s, pde_step=1.317035s, pde_to_ode=0.007212s, ode_from_dolfin=0.002404s, pde_assign_previous_after=0.002927s, total_step=3.972507s

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=2.468525s, ode_state_update=0.012284s, ode_total_step=2.488110s, ode_step=2.659063s, ode_to_dolfin=0.009583s, ode_to_pde=0.005328s, pde_assign_previous_before=0.003177s, pde_set_time=0.004462s, pde_update_matrices=0.009355s, pde_update_rhs=0.964057s, pde_linear_solve=0.359412s, pde_scatter_forward=0.001637s, pde_total_step=1.350014s, pde_step=1.352315s, pde_to_ode=0.007403s, ode_from_dolfin=0.002464s, pde_assign_previous_after=0.003008s, total_step=4.080092s

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

INFO:__main__:Average solver time per time step: 0.01084 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=2.533506s, ode_state_update=0.012607s, ode_total_step=2.553614s, ode_step=2.729295s, ode_to_dolfin=0.009849s, ode_to_pde=0.005469s, pde_assign_previous_before=0.003258s, pde_set_time=0.004575s, pde_update_matrices=0.009355s, pde_update_rhs=0.989468s, pde_linear_solve=0.368772s, pde_scatter_forward=0.001679s, pde_total_step=1.385224s, pde_step=1.387583s, pde_to_ode=0.007600s, ode_from_dolfin=0.002523s, pde_assign_previous_after=0.003085s, total_step=4.187441s

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=2.598405s, ode_state_update=0.012943s, ode_total_step=2.619037s, ode_step=2.799431s, ode_to_dolfin=0.010098s, ode_to_pde=0.005634s, pde_assign_previous_before=0.003339s, pde_set_time=0.004690s, pde_update_matrices=0.009355s, pde_update_rhs=1.014907s, pde_linear_solve=0.378211s, pde_scatter_forward=0.001738s, pde_total_step=1.420564s, pde_step=1.422981s, pde_to_ode=0.007806s, ode_from_dolfin=0.002583s, pde_assign_previous_after=0.003164s, total_step=4.294850s

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

INFO:__main__:Average solver time per time step: 0.01084 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                          |     4.3056
ode_step                            |     2.8064
ode_total_step                      |     2.6255
ode_function_call                   |     2.6049
pde_step                            |     1.4266
pde_total_step                      |     1.4241
pde_update_rhs                      |     1.0175
pde_linear_solve                    |     0.3791
ode_state_update                    |     0.0130
ode_to_dolfin                       |     0.0101
pde_update_matrices                 |     0.0094
pde_to_ode                          |     0.0078
ode_to_pde                          |     0.0056
pde_set_time                        |     0.0047
pde_assign_previous_before          |     0.0033
pde_assign_previous_after           |     0.0032
ode_from_dolfin                     |     0.0026
pde_scatter_forward                 |     0.0017
==================================================

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.