BiV ellipsoid with time dependent pressure and activation - Benchmark#
This example is very similar to the LV example, only that we now will use a biventricular geometry, meaning that there are two different pressure boundary conditions, one for the left ventricle and one for the right ventricle. We will also use different parameters for the Bestel pressure model for the left and right ventricle.
Background#
This example implements the Bi-Ventricular (BiV) version of the benchmark described in the paper [ArosticaNB+25].
The formulation and numerical methods (generalized \(\alpha\)-method, viscoelasticity, active stress, boundary conditions) are identical to the LV case. The primary difference lies in the geometry and the application of distinct pressure loads on the LV and RV endocardial surfaces.
The reader is referred to the LV example for a detailed description of the mathematical models and equations used.
from pathlib import Path
import logging
import os
import numpy as np
import matplotlib.pyplot as plt
from mpi4py import MPI
import dolfinx
from dolfinx import log
from scipy.integrate import solve_ivp
import circulation.bestel
import cardiac_geometries
import cardiac_geometries.geometry
import pulse
logging.basicConfig(level=logging.INFO)
comm = MPI.COMM_WORLD
geodir = outdir / "geometry"
if not geodir.exists():
comm.barrier()
cardiac_geometries.mesh.biv_ellipsoid(
outdir=geodir,
create_fibers=True,
fiber_space="Quadrature_6",
comm=comm,
char_length=1.0,
fiber_angle_epi=-60,
fiber_angle_endo=60,
)
2025-12-01 21:46:46 [debug ] Convert file time-dependent-bestel-biv/geometry/biv_ellipsoid.msh to dolfin
Info : Reading 'time-dependent-bestel-biv/geometry/biv_ellipsoid.msh'...
Info : 64 entities
Info : 439 nodes
Info : 2188 elements
Info : Done reading 'time-dependent-bestel-biv/geometry/biv_ellipsoid.msh'
INFO:ldrb.ldrb:Calculating scalar fields
INFO:ldrb.ldrb:Compute scalar laplacian solutions with the markers:
base: [3]
lv: [4, 7]
rv: [5, 6]
epi: [1, 2]
INFO:ldrb.ldrb: Num vertices: 439
INFO:ldrb.ldrb: Num cells: 1314
INFO:ldrb.ldrb: Apex coord: (0.00, 0.00, -3.75)
INFO:ldrb.ldrb:
Calculating gradients
INFO:ldrb.ldrb:Compute fiber-sheet system
INFO:ldrb.ldrb:Angles:
INFO:ldrb.ldrb:alpha:
endo_lv: 60
epi_lv: -60
endo_septum: 60
epi_septum: -60
endo_rv: 60
epi_rv: -60
INFO:ldrb.ldrb:beta:
endo_lv: 0
epi_lv: 0
endo_septum: 0
epi_septum: 0
endo_rv: 0
epi_rv: 0
INFO:cardiac_geometries.geometry:Reading geometry from time-dependent-bestel-biv/geometry
geo = cardiac_geometries.geometry.Geometry.from_folder(
comm=comm,
folder=geodir,
)
INFO:cardiac_geometries.geometry:Reading geometry from time-dependent-bestel-biv/geometry
In this case we scale the geometry to be in meters
geo.mesh.geometry.x[:] *= 1.5e-2
Now we need to redefine the markers to have so that facets on the endo- and epicardium combine both free wall and the septum.
markers = {"ENDO_LV": [1, 2], "ENDO_RV": [2, 2], "BASE": [3, 2], "EPI": [4, 2]}
marker_values = geo.ffun.values.copy()
marker_values[np.isin(geo.ffun.indices, geo.ffun.find(geo.markers["LV_ENDO_FW"][0]))] = markers["ENDO_LV"][0]
marker_values[np.isin(geo.ffun.indices, geo.ffun.find(geo.markers["LV_SEPTUM"][0]))] = markers["ENDO_LV"][0]
marker_values[np.isin(geo.ffun.indices, geo.ffun.find(geo.markers["RV_ENDO_FW"][0]))] = markers["ENDO_RV"][0]
marker_values[np.isin(geo.ffun.indices, geo.ffun.find(geo.markers["RV_SEPTUM"][0]))] = markers["ENDO_RV"][0]
marker_values[np.isin(geo.ffun.indices, geo.ffun.find(geo.markers["BASE"][0]))] = markers["BASE"][0]
marker_values[np.isin(geo.ffun.indices, geo.ffun.find(geo.markers["LV_EPI_FW"][0]))] = markers["EPI"][0]
marker_values[np.isin(geo.ffun.indices, geo.ffun.find(geo.markers["RV_EPI_FW"][0]))] = markers["EPI"][0]
geo.markers = markers
ffun = dolfinx.mesh.meshtags(
geo.mesh,
geo.ffun.dim,
geo.ffun.indices,
marker_values,
)
geo.ffun = ffun
We create the geometry object and print the volumes of the LV and RV cavities
geometry = pulse.HeartGeometry.from_cardiac_geometries(geo, metadata={"quadrature_degree": 6})
print(geometry.volume("ENDO_LV") * 1e6, geometry.volume("ENDO_RV") * 1e6)
139.6685201557857 121.83581009633141
material_params = pulse.HolzapfelOgden.orthotropic_parameters()
material = pulse.HolzapfelOgden(f0=geo.f0, s0=geo.s0, **material_params) # type: ignore
Ta = pulse.Variable(dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(0.0)), "Pa")
active_model = pulse.ActiveStress(geo.f0, activation=Ta)
comp_model = pulse.compressibility.Compressible2()
viscoeleastic_model = pulse.viscoelasticity.Viscous()
model = pulse.CardiacModel(
material=material,
active=active_model,
compressibility=comp_model,
viscoelasticity=viscoeleastic_model,
)
One difference with the LV example is that we now have two different pressure boundary conditions, one for the LV and one for the RV
traction_lv = pulse.Variable(dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(0.0)), "Pa")
neumann_lv = pulse.NeumannBC(traction=traction_lv, marker=geometry.markers["ENDO_LV"][0])
traction_rv = pulse.Variable(dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(0.0)), "Pa")
neumann_rv = pulse.NeumannBC(traction=traction_rv, marker=geometry.markers["ENDO_RV"][0])
Otherwize we have the same Robin boundary conditions as in the LV example
alpha_epi = pulse.Variable(
dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(1e8)), "Pa / m",
)
robin_epi_u = pulse.RobinBC(value=alpha_epi, marker=geometry.markers["EPI"][0])
beta_epi = pulse.Variable(
dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(5e3)), "Pa s/ m",
)
robin_epi_v = pulse.RobinBC(value=beta_epi, marker=geometry.markers["EPI"][0], damping=True)
alpha_base = pulse.Variable(
dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(1e5)), "Pa / m",
)
robin_base_u = pulse.RobinBC(value=alpha_base, marker=geometry.markers["BASE"][0])
beta_base = pulse.Variable(
dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(5e3)), "Pa s/ m",
)
robin_base_v = pulse.RobinBC(value=beta_base, marker=geometry.markers["BASE"][0], damping=True)
bcs = pulse.BoundaryConditions(neumann=(neumann_lv, neumann_rv), robin=(robin_epi_u, robin_epi_v, robin_base_u, robin_base_v))
problem = pulse.problem.DynamicProblem(model=model, geometry=geometry, bcs=bcs)
log.set_log_level(log.LogLevel.INFO)
problem.solve()
INFO:scifem.solvers:Newton iteration 1: r (abs) = 1.5444104306679715e-05 (tol=1e-06), r (rel) = 1.0 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 2: r (abs) = 2.044647598642954e-07 (tol=1e-06), r (rel) = 0.013239017025795568 (tol=1e-10)
2
dt = problem.parameters["dt"].to_base_units()
times = np.arange(0.0, 1.0, dt)
We now solve the Bestel pressure model using different parameters for the LV and RV
lv_pressure_model = circulation.bestel.BestelPressure(
parameters=dict(
t_sys_pre=0.17,
t_dias_pre=0.484,
gamma=0.005,
a_max=5.0,
a_min=-30.0,
alpha_pre=5.0,
alpha_mid=15.0,
sigma_pre=12000.0,
sigma_mid=16000.0,
),
)
res_lv = solve_ivp(
lv_pressure_model,
[0.0, 1.0],
[0.0],
t_eval=times,
method="Radau",
)
lv_pressure = res_lv.y[0]
INFO:circulation.bestel:
Bestel pressure model
parameters
┏━━━━━━━━━━━━┳━━━━━━━━━┓
┃ Parameter ┃ Value ┃
┡━━━━━━━━━━━━╇━━━━━━━━━┩
│ t_sys_pre │ 0.17 │
│ t_dias_pre │ 0.484 │
│ gamma │ 0.005 │
│ a_max │ 5.0 │
│ a_min │ -30.0 │
│ alpha_pre │ 5.0 │
│ alpha_mid │ 15.0 │
│ sigma_pre │ 12000.0 │
│ sigma_mid │ 16000.0 │
└────────────┴─────────┘
rv_pressure_model = circulation.bestel.BestelPressure(
parameters=dict(
t_sys_pre=0.17,
t_dias_pre=0.484,
gamma=0.005,
a_max=5.0,
a_min=-30.0,
alpha_pre=1.0,
alpha_mid=10.0,
sigma_pre=3000.0,
sigma_mid=4000.0,
),
)
res_rv = solve_ivp(
rv_pressure_model,
[0.0, 1.0],
[0.0],
t_eval=times,
method="Radau",
)
rv_pressure = res_rv.y[0]
INFO:circulation.bestel:
Bestel pressure model
parameters
┏━━━━━━━━━━━━┳━━━━━━━━┓
┃ Parameter ┃ Value ┃
┡━━━━━━━━━━━━╇━━━━━━━━┩
│ t_sys_pre │ 0.17 │
│ t_dias_pre │ 0.484 │
│ gamma │ 0.005 │
│ a_max │ 5.0 │
│ a_min │ -30.0 │
│ alpha_pre │ 1.0 │
│ alpha_mid │ 10.0 │
│ sigma_pre │ 3000.0 │
│ sigma_mid │ 4000.0 │
└────────────┴────────┘
activation_model = circulation.bestel.BestelActivation()
res = solve_ivp(
activation_model,
[0.0, 1.0],
[0.0],
t_eval=times,
method="Radau",
)
# Convert the pressure from Pa to kPa
activation = res.y[0]
INFO:circulation.bestel:
Bestel activation model
parameters
┏━━━━━━━━━━━┳━━━━━━━━━━┓
┃ Parameter ┃ Value ┃
┡━━━━━━━━━━━╇━━━━━━━━━━┩
│ t_sys │ 0.16 │
│ t_dias │ 0.484 │
│ gamma │ 0.005 │
│ a_max │ 5.0 │
│ a_min │ -30.0 │
│ sigma_0 │ 150000.0 │
└───────────┴──────────┘
fig, ax = plt.subplots(2, 1, sharex=True, figsize=(10, 10))
ax[0].plot(times, lv_pressure, label="LV")
ax[0].plot(times, rv_pressure, label="RV")
ax[0].set_title("Pressure")
ax[0].legend()
ax[1].plot(times, activation)
ax[1].set_title("Activation")
fig.savefig(outdir / "pressure_activation.png")
vtx = dolfinx.io.VTXWriter(geometry.mesh.comm, outdir / "displacement.bp", [problem.u], engine="BP4")
vtx.write(0.0)
volume_form = geometry.volume_form(u=problem.u)
lv_volume_form = dolfinx.fem.form(volume_form * geometry.ds(geometry.markers["ENDO_LV"][0]))
rv_volume_form = dolfinx.fem.form(volume_form * geometry.ds(geometry.markers["ENDO_RV"][0]))
[2025-12-01 21:47:07.223] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-12-01 21:47:07.223] [info] Requesting connectivity (3, 0) - (2, 0)
[2025-12-01 21:47:07.223] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-12-01 21:47:07.223] [info] Requesting connectivity (3, 0) - (2, 0)
[2025-12-01 21:47:07.672] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-12-01 21:47:07.672] [info] Requesting connectivity (3, 0) - (2, 0)
[2025-12-01 21:47:07.672] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-12-01 21:47:07.672] [info] Requesting connectivity (3, 0) - (2, 0)
lv_volumes = []
rv_volumes = []
for i, (tai, plv, prv, ti) in enumerate(zip(activation, lv_pressure, rv_pressure, times)):
print(f"Solving for time {ti}, activation {tai}, lv pressure {plv} and rv pressure {prv}")
traction_lv.assign(plv)
traction_rv.assign(prv)
Ta.assign(tai)
problem.solve()
vtx.write((i + 1) * dt)
lv_volumes.append(dolfinx.fem.assemble_scalar(lv_volume_form))
rv_volumes.append(dolfinx.fem.assemble_scalar(rv_volume_form))
if geo.mesh.comm.rank == 0:
fig, ax = plt.subplots(4, 1, figsize=(10, 10))
ax[0].plot(times[:i + 1], lv_pressure[:i + 1], label="LV")
ax[0].plot(times[:i + 1], rv_pressure[:i + 1], label="RV")
ax[0].set_title("Pressure")
ax[1].plot(times[:i + 1], activation[:i + 1], label="Activation")
ax[1].set_title("Activation")
ax[2].plot(times[:i + 1], lv_volumes, label="LV")
ax[2].plot(times[:i + 1], rv_volumes, label="RV")
ax[2].set_title("Volume")
ax[3].plot(lv_volumes, lv_pressure[:i + 1], label="LV")
ax[3].plot(rv_volumes, rv_pressure[:i + 1], label="RV")
for a in ax:
a.legend()
fig.savefig(outdir / "biv_ellipsoid_time_dependent_bestel.png")
plt.close(fig)
if os.getenv("CI") and i > 2:
# Early stopping for CI
break
INFO:scifem.solvers:Newton iteration 1: r (abs) = 1.175915601068858e-05 (tol=1e-06), r (rel) = 1.0 (tol=1e-10)
Solving for time 0.0, activation 0.0, lv pressure 0.0 and rv pressure 0.0
INFO:scifem.solvers:Newton iteration 2: r (abs) = 6.54898426806494e-09 (tol=1e-06), r (rel) = 0.0005569263867332135 (tol=1e-10)
Solving for time 0.001, activation 0.0, lv pressure 11.940199493584645 and rv pressure 2.9716796494482707
INFO:scifem.solvers:Newton iteration 1: r (abs) = 1.0547528142226222e-05 (tol=1e-06), r (rel) = 1.0 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 2: r (abs) = 1.8638807870237084e-08 (tol=1e-06), r (rel) = 0.0017671256828050586 (tol=1e-10)
Solving for time 0.002, activation 0.0, lv pressure 23.761484968094862 and rv pressure 5.887258837637258
INFO:scifem.solvers:Newton iteration 1: r (abs) = 3.8307606197928834e-05 (tol=1e-06), r (rel) = 1.0 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 2: r (abs) = 2.7149119015022267e-08 (tol=1e-06), r (rel) = 0.0007087135352375563 (tol=1e-10)
Solving for time 0.003, activation 0.0, lv pressure 35.46523956820262 and rv pressure 8.748113892915901
INFO:scifem.solvers:Newton iteration 1: r (abs) = 8.106357283748186e-05 (tol=1e-06), r (rel) = 1.0 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 2: r (abs) = 3.961092321288332e-08 (tol=1e-06), r (rel) = 0.0004886402341566689 (tol=1e-10)
References#
[ArosticaNB+25]
Reidmen Aróstica, David Nolte, Aaron Brown, Amadeus Gebauer, Elias Karabelas, Javiera Jilberto, Matteo Salvador, Michele Bucelli, Roberto Piersanti, Kasra Osouli, and others. A software benchmark for cardiac elastodynamics. Computer Methods in Applied Mechanics and Engineering, 435:117485, 2025.