Monolithic 3D-0D coupling: a UK Biobank biventricular mesh in a closed circulation#
This is the biventricular counterpart of Monolithic 3D-0D coupling: an LV in a closed circulation. We couple both ventricles of a clipped UK Biobank atlas mesh to the closed-loop circulation model of Regazzoni et al., solving the displacement, both cavity pressures and all twelve circuit states together in one Newton system. As in the LV demo, the constraint tying each cavity to its chamber is a row of that system, so the two cannot disagree by more than the solver tolerance.
Adding the second ventricle#
Structurally this costs us almost nothing, since the coupling machinery
already takes a list of cavities. Dropping the RV component from the .ode
alongside LV leaves a circuit that expects both p_LV and p_RV, and each
of those is supplied by the Lagrange multiplier its own cavity already
carries. In practice we need two Cavity entries, two ChamberCoupling
entries, and one extra row in the block system.
Why there is no calibration here#
The LV demo spends most of its length measuring the mesh and tuning the circuit’s loading against that measurement. It has to, because an idealized ellipsoid does not behave like a real ventricle: its unloaded cavity holds far more than a person’s, and coupled to the published circuit it fills to nearly 200 mL.
This mesh is a real heart in a real end-diastolic configuration, so we can
get the operating point directly. We seed the circuit with the mesh’s own
end-diastolic volumes, run it alone until it reaches a limit cycle, and then
prestress the mesh to the pressures it settles at. The two then agree at end
diastole by construction, which leaves contractility – TA_SCALE below –
as the only free parameter.
Units#
The UKB mesh comes in millimetres, and we scale it to metres on load, as the other demos built on it do. The chamber coupling converts between the circuit’s millilitres and the mesh’s cubic metres, and it assumes metres, so leaving the mesh in millimetres would couple the circuit to a cavity a billion times the intended size.
import logging
import os
import shutil
from pathlib import Path
import circulation
import dolfinx
import io4dolfinx
import ldrb
import matplotlib.pyplot as plt
import numpy as np
from circulation import bestel, regazzoni2020
from circulation.units import kPa_to_mmHg, mmHg_to_kPa
from matplotlib.gridspec import GridSpec
from scipy.integrate import solve_ivp
# A sibling module in this directory, not a package.
import animation
import cardiac_geometries
import cardiac_geometries.geometry
import pulse
from pulse.circulation import ChamberCoupling, GotranxCirculation, mL, mmHg
circulation.log.setup_logging(logging.INFO)
logging.getLogger("scifem").setLevel(logging.WARNING)
logger = logging.getLogger("pulse")
comm = MPI.COMM_WORLD
_ci = os.getenv("CI", "").strip().lower()
IN_CI = _ci not in ("", "0", "false", "no", "off")
# Inertia works the same way here as in the LV demo. Quasi-static is the
# default, and switching it on also enables the viscous term and the damping
# Robin conditions, without which the cavity pressures ring against their own
# constraints. As there, `PULSE_DYNAMIC=1` sets the flag without editing the
# file.
DYNAMIC = os.getenv("PULSE_DYNAMIC", "0").strip().lower() in ("1", "true", "yes", "on")
ARM = "dynamic" if DYNAMIC else "quasistatic"
# This is the contractility: the Bestel trace sets the shape of the twitch and
# this sets its size. With end diastole pinned by the prestressing, it is what
# ends up deciding the ejection fraction, and we picked the value by running
# the quasi-static arm and reading that ejection fraction off the loop.
TA_SCALE = 1.0
BEAT_LENGTH = 1.0 # s
DT = 0.002 # s
# Two beats are enough to settle the left ventricle: its end-diastolic and
# end-systolic volumes move by half a percent between the first and the second.
# The right ventricle is still drifting by about seven percent, since the
# pulmonary compartment it fills through is more compliant and takes longer to
# settle. Increase this if the right side matters for what you are doing.
NUM_BEATS = 1 if IN_CI else 2
CHAR_LENGTH = 10.0 # mm; the atlas is smooth, so a coarse mesh suffices
outdir = Path("results_monolithic_3d0d_biv")
geodir = Path("ukb-monolithic-3d0d")
outdir.mkdir(exist_ok=True)
Geometry#
We use the mean shape of the atlas (mode=-1, std=0) at end diastole,
clipped at the valve plane so that the mesh has a single BASE surface
rather than four valve annuli, and rotated so that the base normal points
along x. The fibres come from LDRB, with separate angles for the two
ventricles.
if not (geodir / "geometry.bp").exists():
logger.info("Generating the UKB mesh...")
geo = cardiac_geometries.mesh.ukb(
outdir=geodir,
comm=comm,
mode=-1,
std=0,
case="ED",
char_length_max=CHAR_LENGTH,
char_length_min=CHAR_LENGTH,
clipped=True,
)
geo = geo.rotate(target_normal=[1.0, 0.0, 0.0], base_marker="BASE")
system = ldrb.dolfinx_ldrb(
mesh=geo.mesh,
ffun=geo.ffun,
markers=cardiac_geometries.mesh.transform_markers(geo.markers, clipped=True),
alpha_endo_lv=60,
alpha_epi_lv=-60,
alpha_endo_rv=90,
alpha_epi_rv=-25,
beta_endo_lv=-20,
beta_epi_lv=20,
beta_endo_rv=0,
beta_epi_rv=20,
fiber_space="Quadrature_6",
)
if (geodir / "geometry.bp").exists():
shutil.rmtree(geodir / "geometry.bp")
cardiac_geometries.geometry.save_geometry(
path=geodir / "geometry.bp",
mesh=geo.mesh,
ffun=geo.ffun,
markers=geo.markers,
info=geo.info,
f0=system.f0,
s0=system.s0,
n0=system.n0,
)
comm.barrier()
[09/24/26 14:19:35] INFO INFO:pulse:Generating the UKB mesh... 1178313631.py:2
INFO INFO:ukb.surface:Saved ukb-monolithic-3d0d/EPI_ED.stl surface.py:201
INFO INFO:ukb.surface:Saved ukb-monolithic-3d0d/MV_ED.stl surface.py:206
INFO INFO:ukb.surface:Saved ukb-monolithic-3d0d/AV_ED.stl surface.py:206
INFO INFO:ukb.surface:Saved ukb-monolithic-3d0d/TV_ED.stl surface.py:206
INFO INFO:ukb.surface:Saved ukb-monolithic-3d0d/PV_ED.stl surface.py:206
INFO INFO:ukb.surface:Saved ukb-monolithic-3d0d/LV_ED.stl surface.py:214
INFO INFO:ukb.surface:Saved ukb-monolithic-3d0d/RV_ED.stl surface.py:214
INFO INFO:ukb.surface:Saved ukb-monolithic-3d0d/RVFW_ED.stl surface.py:214
Warning: PLY writer doesn't support multidimensional point data yet. Skipping Normals.
Warning: PLY doesn't support 64-bit integers. Casting down to 32-bit.
Warning: PLY writer doesn't support multidimensional point data yet. Skipping Normals.
Warning: PLY doesn't support 64-bit integers. Casting down to 32-bit.
Warning: PLY writer doesn't support multidimensional point data yet. Skipping Normals.
Warning: PLY doesn't support 64-bit integers. Casting down to 32-bit.
0
INFO INFO:ukb.mesh:Creating clipped mesh for ED with char_length_max=10.0, char_length_min=10.0 mesh.py:264
2026-09-24 14:19:36 [debug ] Convert file ukb-monolithic-3d0d/ED_clipped.msh to dolfin
Info : Reading 'ukb-monolithic-3d0d/ED_clipped.msh'...
Info : 11 entities
Info : 696 nodes
Info : 3436 elements
Info : 3 parametrizations
Info : [ 0%] Processing parametrizations
Info : [ 10%] Processing parametrizations
Info : [ 40%] Processing parametrizations
Info : Done reading 'ukb-monolithic-3d0d/ED_clipped.msh'
INFO INFO:ldrb.ldrb:Compute scalar laplacian solutions with the markers: ldrb.py:619 lv: [1] rv: [2] epi: [3] base: [4]
INFO INFO:ldrb.ldrb:alpha: ldrb.py:86 endo_lv: 60 epi_lv: -60 endo_septum: 60 epi_septum: -60 endo_rv: 60 epi_rv: -60
INFO INFO:ldrb.ldrb:beta: ldrb.py:98 endo_lv: 0 epi_lv: 0 endo_septum: 0 epi_septum: 0 endo_rv: 0 epi_rv: 0
2026-09-24 14:19:36 [debug ] Write f0: f0
2026-09-24 14:19:36 [debug ] Write s0: s0
2026-09-24 14:19:36 [debug ] Write n0: n0
2026-09-24 14:19:36 [debug ] Write lv: f
2026-09-24 14:19:36 [debug ] Write rv: f
2026-09-24 14:19:36 [debug ] Write epi: f
2026-09-24 14:19:36 [debug ] Write lv_rv: f
2026-09-24 14:19:36 [debug ] Write apex: f
2026-09-24 14:19:36 [debug ] Write lv_scalar: f
2026-09-24 14:19:36 [debug ] Write rv_scalar: f
2026-09-24 14:19:36 [debug ] Write epi_scalar: f
2026-09-24 14:19:36 [debug ] Write lv_rv_scalar: f
2026-09-24 14:19:36 [debug ] Write lv_gradient: f
2026-09-24 14:19:36 [debug ] Write rv_gradient: f
2026-09-24 14:19:36 [debug ] Write epi_gradient: f
2026-09-24 14:19:36 [debug ] Write apex_gradient: f
2026-09-24 14:19:36 [debug ] Write markers_scalar: f
2026-09-24 14:19:37 [info ] Rotated geometry. Base normal [-0.71608437 0.54439464 0.43687258] aligned to [1.0, 0.0, 0.0]
2026-09-24 14:19:37 [debug ] Rotation matrix:
[[-0.71608437 0.54439464 0.43687258]
[-0.54439464 -0.04385067 -0.83768227]
[-0.43687258 -0.83768227 0.32776631]]
INFO INFO:ldrb.ldrb:Compute scalar laplacian solutions with the markers: ldrb.py:619 lv: [1] rv: [2] epi: [3] base: [4]
geo = cardiac_geometries.geometry.Geometry.from_folder(comm=comm, folder=geodir)
INFO INFO:cardiac_geometries.geometry:Reading geometry from ukb-monolithic-3d0d geometry.py:535
# We rotate here, after loading, rather than relying on the rotation in the
# generation step above. That step does rotate the mesh, but the folder also
# holds the `.msh` it was built from, and that unrotated mesh is what
# `from_folder` gives back. Rotating at this point makes the orientation a
# property of what we actually solve on. It matters because the sliding-base
# condition below constrains a single displacement component: on an unrotated
# mesh it would hold the base in a plane that cuts through the ventricle at an
# angle instead of in the base plane itself.
geo = geo.rotate(target_normal=[1.0, 0.0, 0.0], base_marker="BASE")
geo.mesh.geometry.x[:] *= 1e-3 # mm -> m
geometry = pulse.HeartGeometry.from_cardiac_geometries(geo, metadata={"quadrature_degree": 6})
2026-09-24 14:19:37 [info ] Rotated geometry. Base normal [-1.00000000e+00 -7.77917009e-16 6.58036490e-16] aligned to [1.0, 0.0, 0.0]
2026-09-24 14:19:37 [debug ] Rotation matrix:
[[-1. 0. 0.]
[ 0. -1. 0.]
[ 0. 0. 1.]]
up = animation.base_normal(geometry, "BASE")
if abs(up[0]) < 0.99:
raise RuntimeError(
f"the base normal is {up.round(3)}, not the x axis the sliding-base "
"condition assumes -- the rotation above did not take effect",
)
EDV = {
chamber: comm.allreduce(geometry.volume(chamber), op=MPI.SUM) for chamber in ("LV", "RV")
}
logger.info(f"Mesh end-diastolic volumes: LV {EDV['LV'] / mL:.1f} mL, RV {EDV['RV'] / mL:.1f} mL")
[09/24/26 14:19:38] INFO INFO:pulse:Mesh end-diastolic volumes: LV 108.9 mL, RV 74.9 mL 3936440030.py:4
def build_model(f0, s0, Ta):
material_params = pulse.HolzapfelOgden.transversely_isotropic_parameters()
material = pulse.HolzapfelOgden(f0=f0, s0=s0, **material_params) # type: ignore[arg-type]
# This does nothing without a strain rate, so the static solves below are
# unaffected.
viscoelasticity = (
pulse.viscoelasticity.Viscous() if DYNAMIC else pulse.viscoelasticity.NoneViscoElasticity()
)
return pulse.CardiacModel(
material=material,
active=pulse.ActiveStress(f0, activation=Ta, formulation=pulse.ActiveStressFormulation.stretch),
compressibility=pulse.Compressible(),
viscoelasticity=viscoelasticity,
)
def robin_bcs():
def spring(marker, value, damping=False):
return pulse.RobinBC(
value=pulse.Variable(
dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(value)),
"Pa s/ m" if damping else "Pa / m",
),
marker=geometry.markers[marker][0],
damping=damping,
)
bcs = [spring("EPI", 1.0e5), spring("BASE", 1.0e6)]
if DYNAMIC:
bcs += [spring("EPI", 5.0e3, damping=True), spring("BASE", 5.0e3, damping=True)]
return tuple(bcs)
def sliding_base(V: dolfinx.fem.FunctionSpace):
"""Hold the base in its own plane but let it slide within it.
The mesh is rotated so the base normal is x, which is what makes this one
component rather than a projection.
"""
facets = geometry.facet_tags.find(geometry.markers["BASE"][0])
dofs = dolfinx.fem.locate_dofs_topological(V.sub(0), 2, facets)
return [dolfinx.fem.dirichletbc(0.0, dofs, V.sub(0))]
Activation#
We solve the Bestel twitch once up front, since it depends on time alone and prescribing it therefore introduces no coupling error.
times = np.arange(0.0, BEAT_LENGTH, DT)
activation = solve_ivp(
bestel.BestelActivation(),
[0.0, BEAT_LENGTH],
[0.0],
t_eval=times,
method="Radau",
).y[0]
logger.info(f"Peak activation: {TA_SCALE * activation.max() * 1e-3:.1f} kPa")
INFO INFO:circulation.bestel: bestel.py:72 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 │ └───────────┴──────────┘
INFO INFO:pulse:Peak activation: 118.0 kPa 3671970891.py:9
def activation_at(t: float) -> float:
return TA_SCALE * float(np.interp(t % BEAT_LENGTH, times, activation))
The operating point#
We seed the circuit with this mesh’s own end-diastolic volumes and run it by itself until it reaches a limit cycle. Then we prestress the mesh to the end-diastolic pressures it arrives at, so that afterwards the two agree on both volumes and both pressures at end diastole. This takes the place of the calibration that the LV demo needs a separate module for, and it only works because the geometry is a real one.
state_file = outdir / "circ_state.npy"
if comm.rank == 0 and not state_file.exists():
standalone = regazzoni2020.Regazzoni2020(parameters={"HR": 1.0 / BEAT_LENGTH}, add_units=False)
history_0d = standalone.solve(
num_beats=10,
initial_state={"V_LV": EDV["LV"] / mL, "V_RV": EDV["RV"] / mL},
dt=0.001,
)
np.save(
state_file,
{
"state": dict(zip(standalone.state_names(), standalone.state)),
"p_LV_ED": float(history_0d["p_LV"][-1]),
"p_RV_ED": float(history_0d["p_RV"][-1]),
},
allow_pickle=True,
)
comm.barrier()
INFO INFO:circulation.base: base.py:134 Circulation model parameters (Regazzoni2020) ┏━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓ ┃ Parameter ┃ Value ┃ ┡━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┩ │ HR │ 1.0 │ │ chambers.LA.EA │ 0.07 millimeter_Hg / milliliter │ │ chambers.LA.EB │ 0.18 millimeter_Hg / milliliter │ │ chambers.LA.TC │ 0.17 second │ │ chambers.LA.TR │ 0.17 second │ │ chambers.LA.tC │ 0.9 second │ │ chambers.LA.V0 │ 4.0 milliliter │ │ chambers.LV.EA │ 4.482 millimeter_Hg / milliliter │ │ chambers.LV.EB │ 0.17 millimeter_Hg / milliliter │ │ chambers.LV.TC │ 0.25 second │ │ chambers.LV.TR │ 0.4 second │ │ chambers.LV.tC │ 0.1 second │ │ chambers.LV.V0 │ 42.0 milliliter │ │ chambers.RA.EA │ 0.06 millimeter_Hg / milliliter │ │ chambers.RA.EB │ 0.07 millimeter_Hg / milliliter │ │ chambers.RA.TC │ 0.17 second │ │ chambers.RA.TR │ 0.17 second │ │ chambers.RA.tC │ 0.9 second │ │ chambers.RA.V0 │ 4.0 milliliter │ │ chambers.RV.EA │ 0.2 millimeter_Hg / milliliter │ │ chambers.RV.EB │ 0.029 millimeter_Hg / milliliter │ │ chambers.RV.TC │ 0.25 second │ │ chambers.RV.TR │ 0.4 second │ │ chambers.RV.tC │ 0.1 second │ │ chambers.RV.V0 │ 16.0 milliliter │ │ valves.MV.Rmin │ 0.0075 millimeter_Hg * second / milliliter │ │ valves.MV.Rmax │ 75006.2 millimeter_Hg * second / milliliter │ │ valves.AV.Rmin │ 0.0075 millimeter_Hg * second / milliliter │ │ valves.AV.Rmax │ 75006.2 millimeter_Hg * second / milliliter │ │ valves.TV.Rmin │ 0.0075 millimeter_Hg * second / milliliter │ │ valves.TV.Rmax │ 75006.2 millimeter_Hg * second / milliliter │ │ valves.PV.Rmin │ 0.0075 millimeter_Hg * second / milliliter │ │ valves.PV.Rmax │ 75006.2 millimeter_Hg * second / milliliter │ │ circulation.SYS.R_AR │ 0.733 millimeter_Hg * second / milliliter │ │ circulation.SYS.C_AR │ 1.372 milliliter / millimeter_Hg │ │ circulation.SYS.R_VEN │ 0.32 millimeter_Hg * second / milliliter │ │ circulation.SYS.C_VEN │ 11.363 milliliter / millimeter_Hg │ │ circulation.SYS.L_AR │ 0.005 millimeter_Hg * second ** 2 / milliliter │ │ circulation.SYS.L_VEN │ 0.0005 millimeter_Hg * second ** 2 / milliliter │ │ circulation.PUL.R_AR │ 0.046 millimeter_Hg * second / milliliter │ │ circulation.PUL.C_AR │ 20.0 milliliter / millimeter_Hg │ │ circulation.PUL.R_VEN │ 0.0015 millimeter_Hg * second / milliliter │ │ circulation.PUL.C_VEN │ 16.0 milliliter / millimeter_Hg │ │ circulation.PUL.L_AR │ 0.0005 millimeter_Hg * second ** 2 / milliliter │ │ circulation.PUL.L_VEN │ 0.0005 millimeter_Hg * second ** 2 / milliliter │ │ circulation.external.start_withdrawal │ 0.0 second │ │ circulation.external.end_withdrawal │ 0.0 second │ │ circulation.external.start_infusion │ 0.0 second │ │ circulation.external.end_infusion │ 0.0 second │ │ circulation.external.flow_withdrawal │ 0.0 milliliter / second │ │ circulation.external.flow_infusion │ 0.0 milliliter / second │ └───────────────────────────────────────┴─────────────────────────────────────────────────┘
INFO INFO:circulation.base: base.py:141 Circulation model initial states (Regazzoni2020) ┏━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓ ┃ State ┃ Value ┃ ┡━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┩ │ V_LA │ 87.183 milliliter │ │ V_LV │ 118.52 milliliter │ │ V_RA │ 86.833 milliliter │ │ V_RV │ 166.177 milliliter │ │ p_AR_SYS │ 87.675 millimeter_Hg │ │ p_VEN_SYS │ 35.898 millimeter_Hg │ │ p_AR_PUL │ 19.545 millimeter_Hg │ │ p_VEN_PUL │ 15.004 millimeter_Hg │ │ Q_AR_SYS │ 71.104 milliliter / second │ │ Q_VEN_SYS │ 94.039 milliliter / second │ │ Q_AR_PUL │ 94.084 milliliter / second │ │ Q_VEN_PUL │ 473.279 milliliter / second │ └───────────┴─────────────────────────────┘
cached = np.load(state_file, allow_pickle=True).item()
circ_state = cached["state"]
p_ED = {"LV": mmHg_to_kPa(cached["p_LV_ED"]), "RV": mmHg_to_kPa(cached["p_RV_ED"])}
logger.info(f"End-diastolic pressures from the circuit: "
f"LV {p_ED['LV']:.2f} kPa, RV {p_ED['RV']:.2f} kPa")
[09/24/26 14:19:40] INFO INFO:pulse:End-diastolic pressures from the circuit: LV 2.28 kPa, RV 0.59 kPa 1932125848.py:4
Prestressing#
We recover the unloaded configuration by unloading both cavities together.
Ta = pulse.Variable(dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(0.0)), "Pa")
traction = {
chamber: pulse.Variable(dolfinx.fem.Constant(geometry.mesh, 0.0), "kPa")
for chamber in ("LV", "RV")
}
prestress_fname = outdir / "prestress_biv.bp"
if not prestress_fname.exists():
logger.info("Prestressing to recover the unloaded reference configuration...")
prestress_problem = pulse.unloading.PrestressProblem(
geometry=geometry,
model=build_model(geo.f0, geo.s0, Ta),
bcs=pulse.BoundaryConditions(
robin=robin_bcs(),
dirichlet=(sliding_base,),
neumann=tuple(
pulse.NeumannBC(traction=traction[c], marker=geometry.markers[c][0])
for c in ("LV", "RV")
),
),
parameters={"u_space": "P_2", "mesh_unit": "m"},
targets=[
pulse.unloading.TargetPressure(traction=traction[c], target=p_ED[c], name=c)
for c in ("LV", "RV")
],
ramp_steps=20,
)
u_pre = prestress_problem.unload()
io4dolfinx.write_function_on_input_mesh(prestress_fname, u_pre, time=0.0, name="u_pre")
comm.barrier()
INFO INFO:pulse:Prestressing to recover the unloaded reference configuration... 2427716160.py:3
[09/24/26 14:19:45] INFO INFO:pulse.unloading:Ramping LV traction to 0.0000 unloading.py:499
INFO INFO:pulse.unloading:Ramping RV traction to 0.0000 unloading.py:499
INFO INFO:pulse.unloading:Ramping LV traction to 0.1198 unloading.py:499
INFO INFO:pulse.unloading:Ramping RV traction to 0.0311 unloading.py:499
[09/24/26 14:19:47] INFO INFO:pulse.unloading:Ramping LV traction to 0.2395 unloading.py:499
INFO INFO:pulse.unloading:Ramping RV traction to 0.0621 unloading.py:499
[09/24/26 14:19:49] INFO INFO:pulse.unloading:Ramping LV traction to 0.3593 unloading.py:499
INFO INFO:pulse.unloading:Ramping RV traction to 0.0932 unloading.py:499
[09/24/26 14:19:51] INFO INFO:pulse.unloading:Ramping LV traction to 0.4790 unloading.py:499
INFO INFO:pulse.unloading:Ramping RV traction to 0.1243 unloading.py:499
[09/24/26 14:19:53] INFO INFO:pulse.unloading:Ramping LV traction to 0.5988 unloading.py:499
INFO INFO:pulse.unloading:Ramping RV traction to 0.1553 unloading.py:499
[09/24/26 14:19:55] INFO INFO:pulse.unloading:Ramping LV traction to 0.7186 unloading.py:499
INFO INFO:pulse.unloading:Ramping RV traction to 0.1864 unloading.py:499
[09/24/26 14:19:56] INFO INFO:pulse.unloading:Ramping LV traction to 0.8383 unloading.py:499
INFO INFO:pulse.unloading:Ramping RV traction to 0.2175 unloading.py:499
[09/24/26 14:19:58] INFO INFO:pulse.unloading:Ramping LV traction to 0.9581 unloading.py:499
INFO INFO:pulse.unloading:Ramping RV traction to 0.2486 unloading.py:499
[09/24/26 14:20:00] INFO INFO:pulse.unloading:Ramping LV traction to 1.0778 unloading.py:499
INFO INFO:pulse.unloading:Ramping RV traction to 0.2796 unloading.py:499
[09/24/26 14:20:02] INFO INFO:pulse.unloading:Ramping LV traction to 1.1976 unloading.py:499
INFO INFO:pulse.unloading:Ramping RV traction to 0.3107 unloading.py:499
[09/24/26 14:20:04] INFO INFO:pulse.unloading:Ramping LV traction to 1.3174 unloading.py:499
INFO INFO:pulse.unloading:Ramping RV traction to 0.3418 unloading.py:499
[09/24/26 14:20:06] INFO INFO:pulse.unloading:Ramping LV traction to 1.4371 unloading.py:499
INFO INFO:pulse.unloading:Ramping RV traction to 0.3728 unloading.py:499
[09/24/26 14:20:07] INFO INFO:pulse.unloading:Ramping LV traction to 1.5569 unloading.py:499
INFO INFO:pulse.unloading:Ramping RV traction to 0.4039 unloading.py:499
[09/24/26 14:20:08] INFO INFO:pulse.unloading:Ramping LV traction to 1.6766 unloading.py:499
INFO INFO:pulse.unloading:Ramping RV traction to 0.4350 unloading.py:499
[09/24/26 14:20:10] INFO INFO:pulse.unloading:Ramping LV traction to 1.7964 unloading.py:499
INFO INFO:pulse.unloading:Ramping RV traction to 0.4660 unloading.py:499
[09/24/26 14:20:12] INFO INFO:pulse.unloading:Ramping LV traction to 1.9162 unloading.py:499
INFO INFO:pulse.unloading:Ramping RV traction to 0.4971 unloading.py:499
[09/24/26 14:20:13] INFO INFO:pulse.unloading:Ramping LV traction to 2.0359 unloading.py:499
INFO INFO:pulse.unloading:Ramping RV traction to 0.5282 unloading.py:499
[09/24/26 14:20:14] INFO INFO:pulse.unloading:Ramping LV traction to 2.1557 unloading.py:499
INFO INFO:pulse.unloading:Ramping RV traction to 0.5592 unloading.py:499
[09/24/26 14:20:16] INFO INFO:pulse.unloading:Ramping LV traction to 2.2754 unloading.py:499
INFO INFO:pulse.unloading:Ramping RV traction to 0.5903 unloading.py:499
V_disp = dolfinx.fem.functionspace(geometry.mesh, ("Lagrange", 2, (3,)))
u_pre = dolfinx.fem.Function(V_disp)
io4dolfinx.read_function(prestress_fname, u_pre, time=0.0, name="u_pre")
geometry.deform(u_pre)
f0 = pulse.utils.map_vector_field(f=geo.f0, u=u_pre, normalize=True, name="f0_unloaded")
s0 = pulse.utils.map_vector_field(f=geo.s0, u=u_pre, normalize=True, name="s0_unloaded")
unloaded = {c: comm.allreduce(geometry.volume(c), op=MPI.SUM) for c in ("LV", "RV")}
logger.info(f"Unloaded volumes: LV {unloaded['LV'] / mL:.1f} mL, RV {unloaded['RV'] / mL:.1f} mL")
[09/24/26 14:20:17] INFO INFO:pulse:Unloaded volumes: LV 71.5 mL, RV 59.4 mL 2634031907.py:2
Inflation to end diastole#
Here we go from the unloaded configuration back to the volumes the circuit was seeded with. We prescribe those volumes rather than coupling them, since this is a ramp and not part of the time stepping.
model = build_model(f0, s0, Ta)
bcs = pulse.BoundaryConditions(robin=robin_bcs(), dirichlet=(sliding_base,))
inflation_volume = {
c: dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(unloaded[c]))
for c in ("LV", "RV")
}
inflation = pulse.problem.StaticProblem(
model=model,
geometry=geometry,
bcs=bcs,
cavities=[
pulse.problem.Cavity(marker=c, volume=inflation_volume[c]) for c in ("LV", "RV")
],
parameters={"mesh_unit": "m"},
)
inflation.solve()
True
for fraction in np.linspace(0.0, 1.0, 20)[1:]:
for c in ("LV", "RV"):
inflation_volume[c].value = unloaded[c] + fraction * (EDV[c] - unloaded[c])
if not inflation.solve():
raise RuntimeError(f"inflation failed at {fraction:.2f} of the way to end diastole")
p_inflated = [float(p.x.array[0]) for p in inflation.cavity_pressures]
logger.info(
f"Inflated to LV {comm.allreduce(geometry.volume('LV', u=inflation.u), op=MPI.SUM) / mL:.1f} mL "
f"at {p_inflated[0] / mmHg:.1f} mmHg, "
f"RV {comm.allreduce(geometry.volume('RV', u=inflation.u), op=MPI.SUM) / mL:.1f} mL "
f"at {p_inflated[1] / mmHg:.1f} mmHg",
)
[09/24/26 14:22:06] INFO INFO:pulse:Inflated to LV 108.9 mL at 13.9 mmHg, RV 74.9 mL at 4.1 mmHg 3074010101.py:2
The coupled problem#
Both chamber closures come out of the circuit, and the two cavity pressures take their place.
circulation_model = GotranxCirculation(
ode_file=regazzoni2020.ODE_FILE,
parameters=regazzoni2020.flat_ode_parameters(
circulation.base.remove_units(regazzoni2020.Regazzoni2020.default_parameters())
| {"HR": 1.0 / BEAT_LENGTH},
),
drop_components=("timing", "LV", "RV"),
)
beat_phase = dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(0.0))
2026-09-24 14:22:07 [info ] Load ode /dolfinx-env/lib/python3.12/site-packages/circulation/regazzoni2020.ode
2026-09-24 14:22:07 [info ] Num states 12
2026-09-24 14:22:07 [info ] Num parameters 56
coupled_parameters = {"mesh_unit": "m", "circulation_scheme": "backward_euler"}
if DYNAMIC:
coupled_parameters |= {
"rho": pulse.Variable(1e3, "kg/m^3"),
"dt": pulse.Variable(DT, "s"),
}
Problem = pulse.problem.DynamicProblem if DYNAMIC else pulse.problem.StaticProblem
problem = Problem(
model=model,
geometry=geometry,
bcs=bcs,
cavities=[pulse.problem.Cavity(marker=c, volume=None) for c in ("LV", "RV")],
circulation=circulation_model,
chambers=[
ChamberCoupling(marker="LV", volume_state="V_LV", pressure_missing="p_LV"),
ChamberCoupling(marker="RV", volume_state="V_RV", pressure_missing="p_RV"),
],
circulation_missing={"beat_phase": beat_phase},
parameters=coupled_parameters,
)
# We start from the inflated configuration and from the circuit state that
# matches it.
problem.u.x.array[:] = inflation.u.x.array
problem.u_old.x.array[:] = inflation.u.x.array
for i, pressure in enumerate(p_inflated):
problem.cavity_pressures[i].x.array[:] = pressure
problem.cavity_pressures_old[i].x.array[:] = pressure
if DYNAMIC:
# The inflation gives us a configuration but no motion, so we start from
# rest here as well.
problem.v_old.x.array[:] = 0.0
problem.a_old.x.array[:] = 0.0
names = list(circulation_model.state_names)
for name, state, state_old in zip(
names,
problem.circulation_states,
problem.circulation_states_old,
):
# We take the two chamber volumes from the mesh rather than the circuit.
value = EDV[name[2:]] / mL if name in ("V_LV", "V_RV") else float(circ_state[name])
state.x.array[:] = value
state_old.x.array[:] = value
problem.circulation_dt.value = DT
Stepping#
index = {name: i for i, name in enumerate(names)}
history: dict[str, list[float]] = {
key: [] for key in
(
"time", "V_LV", "V_RV", "p_LV", "p_RV", "Ta",
"iterations", "constraint",
)
}
# We keep the moving geometry every few steps so that `make_animations.py` can
# render it afterwards. Nothing is recorded under CI, where the run is two
# steps rather than a whole beat, so the video on the page comes from a saved
# run instead.
recorder = animation.FrameRecorder(geometry.mesh, every=5, enabled=not IN_CI, up=up)
max_steps = 2 if IN_CI else int(NUM_BEATS * BEAT_LENGTH / DT)
t = 0.0
for step in range(max_steps):
t += DT
problem.circulation_time.value = t
beat_phase.value = t % BEAT_LENGTH
Ta.assign(activation_at(t))
if not problem.solve():
raise RuntimeError(f"Monolithic solve failed at t={t:.4f}")
worst = 0.0
for i, chamber in enumerate(("LV", "RV")):
volume = comm.allreduce(geometry.volume(chamber, u=problem.u), op=MPI.SUM)
state = float(problem.circulation_states[index[f"V_{chamber}"]].x.array[0]) * mL
history[f"V_{chamber}"].append(state / mL)
history[f"p_{chamber}"].append(float(problem.cavity_pressures[i].x.array[0]) / mmHg)
worst = max(worst, abs(volume - state) / state)
history["time"].append(t)
history["Ta"].append(float(Ta.value.value))
history["iterations"].append(int(problem.problem.solver.getIterationNumber()))
history["constraint"].append(worst)
recorder.record(problem.u, t, step)
if step % 50 == 0:
logger.info(
f"t={t:.3f} LV {history['V_LV'][-1]:6.1f} mL {history['p_LV'][-1]:7.1f} mmHg "
f"RV {history['V_RV'][-1]:6.1f} mL {history['p_RV'][-1]:6.1f} mmHg "
f"Ta={history['Ta'][-1] * 1e-3:5.1f} kPa constraint={worst:.1e}",
)
[09/24/26 14:22:29] INFO INFO:pulse:t=0.002 LV 109.3 mL 14.2 mmHg RV 75.5 mL 4.2 mmHg Ta= 0.0 kPa constraint=3.9e-11 197932682.py:27
logger.info(f"Worst constraint violation over the run: {max(history['constraint']):.3e}")
[09/24/26 14:22:30] INFO INFO:pulse:Worst constraint violation over the run: 3.930e-11 4234067209.py:1
saved = recorder.save(outdir / f"frames-{ARM}.npz")
if saved is not None:
logger.info(f"Saved {len(recorder.times)} frames of the moving geometry to {saved}")
for chamber in ("LV", "RV"):
V = np.asarray(history[f"V_{chamber}"])
p = np.asarray(history[f"p_{chamber}"])
if V.size > 10:
EDV_run, ESV_run = float(V.max()), float(V.min())
logger.info(
f"{chamber}: EDV {EDV_run:.1f} mL, ESV {ESV_run:.1f} mL, "
f"SV {EDV_run - ESV_run:.1f} mL, EF {100 * (1 - ESV_run / EDV_run):.1f}%, "
f"peak {p.max():.1f} mmHg",
)
if comm.rank == 0:
np.savez(
outdir / f"traces_biv-{ARM}.npz",
**{k: np.asarray(v) for k, v in history.items()},
)
fig = plt.figure(layout="constrained", figsize=(11, 8))
gs = GridSpec(3, 2, figure=fig)
ax_loop = fig.add_subplot(gs[:, 0])
ax_p = fig.add_subplot(gs[0, 1])
ax_v = fig.add_subplot(gs[1, 1])
ax_ta = fig.add_subplot(gs[2, 1])
for chamber, colour in (("LV", "crimson"), ("RV", "steelblue")):
ax_loop.plot(
history[f"V_{chamber}"], history[f"p_{chamber}"], color=colour,
label=chamber, linewidth=1.1,
)
ax_p.plot(history["time"], history[f"p_{chamber}"], color=colour, label=chamber)
ax_v.plot(history["time"], history[f"V_{chamber}"], color=colour, label=chamber)
ax_loop.set_xlabel("V [mL]")
ax_loop.set_ylabel("p [mmHg]")
ax_loop.set_title(f"Pressure-volume loops ({ARM})")
ax_loop.legend()
ax_p.set_ylabel("p [mmHg]")
ax_p.legend(fontsize="x-small")
ax_v.set_ylabel("V [mL]")
ax_ta.plot(history["time"], np.asarray(history["Ta"]) * 1e-3, color="0.3")
ax_ta.set_ylabel("Ta [kPa]")
ax_ta.set_xlabel("Time [s]")
fig.savefig(outdir / f"monolithic_3d0d_biv-{ARM}.png", dpi=140)
plt.close(fig)
logger.info("Done.")
INFO INFO:pulse:Done. 866737322.py:1
A whole beat#
As in the LV demo, the figure and video come from a full run kept in
_static/ rather than from the two steps this page takes under CI:
python3 monolithic_3d0d_biv.py
python3 make_animations.py monolithic_3d0d_biv
Fig. 6 Both ventricles over two beats. The left ejects 70 mL against a peak of about
100 mmHg, and the right nearly as much against a quarter of that. The left
loop closes on itself, while the right is still drifting, for the reason
given at NUM_BEATS.#