BiV ellipsoid coupled to a 0D circulatory model and a 0D cell model#
This example is similar to the LV only example.
from pathlib import Path
from mpi4py import MPI
import dolfinx
import logging
import os
from functools import lru_cache
import circulation
from dolfinx import log
import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec
import matplotlib.pyplot as plt
import numpy as np
import gotranx
import adios4dolfinx
import cardiac_geometries
import cardiac_geometries.geometry
import fenicsx_pulse
circulation.log.setup_logging(logging.INFO)
logger = logging.getLogger("pulse")
comm = MPI.COMM_WORLD
geodir = Path("biv_ellipsoid-time-dependent")
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,
)
# If the folder already exist, then we just load the geometry
geo = cardiac_geometries.geometry.Geometry.from_folder(
comm=comm,
folder=geodir,
)
geo.mesh.geometry.x[:] *= 3e-2
[03/05/25 19:48:44] INFO INFO:cardiac_geometries.geometry:Reading geometry from biv_ellipsoid-time-dependent geometry.py:145
geometry = fenicsx_pulse.HeartGeometry.from_cardiac_geometries(geo, metadata={"quadrature_degree": 6})
Next we create the material object, and we will use the transversely isotropic version of the Holzapfel Ogden model
material_params = fenicsx_pulse.HolzapfelOgden.transversely_isotropic_parameters()
# material_params = fenicsx_pulse.HolzapfelOgden.orthotropic_parameters()
material = fenicsx_pulse.HolzapfelOgden(f0=geo.f0, s0=geo.s0, **material_params) # type: ignore
We use an active stress approach with 30% transverse active stress (see fenicsx_pulse.active_stress.transversely_active_stress()
)
Ta = fenicsx_pulse.Variable(dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(0.0)), "kPa")
active_model = fenicsx_pulse.ActiveStress(geo.f0, activation=Ta)
We use an incompressible model
comp_model = fenicsx_pulse.compressibility.Compressible2()
viscoeleastic_model = fenicsx_pulse.viscoelasticity.Viscous()
and assembles the CardiacModel
model = fenicsx_pulse.CardiacModel(
material=material,
active=active_model,
compressibility=comp_model,
viscoelasticity=viscoeleastic_model,
)
alpha_epi = fenicsx_pulse.Variable(
dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(1e8)), "Pa / m",
)
robin_epi = fenicsx_pulse.RobinBC(value=alpha_epi, marker=geometry.markers["EPI"][0])
alpha_base = fenicsx_pulse.Variable(
dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(1e5)), "Pa / m",
)
robin_base = fenicsx_pulse.RobinBC(value=alpha_base, marker=geometry.markers["BASE"][0])
lvv_initial = geometry.volume("ENDO_LV")
lv_volume = dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(lvv_initial))
lv_cavity = fenicsx_pulse.problem.Cavity(marker="ENDO_LV", volume=lv_volume)
rvv_initial = geometry.volume("ENDO_RV")
rv_volume = dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(rvv_initial))
rv_cavity = fenicsx_pulse.problem.Cavity(marker="ENDO_RV", volume=rv_volume)
cavities = [lv_cavity, rv_cavity]
parameters = {"base_bc": fenicsx_pulse.problem.BaseBC.free, "mesh_unit": "m"}
static = True
if static:
outdir = Path("biv_ellipsoid_time_dependent_circulation_static")
bcs = fenicsx_pulse.BoundaryConditions(robin=(robin_epi, robin_base))
problem = fenicsx_pulse.problem.StaticProblem(model=model, geometry=geometry, bcs=bcs, cavities=cavities, parameters=parameters)
else:
outdir = Path("biv_ellipsoid_time_dependent_circulation_dynamic")
beta_epi = fenicsx_pulse.Variable(
dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(5e3)), "Pa s/ m",
)
robin_epi_v = fenicsx_pulse.RobinBC(value=beta_epi, marker=geometry.markers["EPI"][0], damping=True)
beta_base = fenicsx_pulse.Variable(
dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(5e3)), "Pa s/ m",
)
robin_base_v = fenicsx_pulse.RobinBC(value=beta_base, marker=geometry.markers["BASE"][0], damping=True)
bcs = fenicsx_pulse.BoundaryConditions(robin=(robin_epi, robin_epi_v, robin_base, robin_base_v))
problem = fenicsx_pulse.problem.DynamicProblem(model=model, geometry=geometry, bcs=bcs, cavities=cavities, parameters=parameters)
outdir.mkdir(exist_ok=True)
Now we can solve the problem
log.set_log_level(log.LogLevel.INFO)
problem.solve()
INFO INFO:scifem.solvers:Newton iteration 1: r (abs) = 15.92761094832744 (tol=1e-06), r (rel) = 1.0 (tol=1e-10) solvers.py:212
INFO INFO:scifem.solvers:Newton iteration 2: r (abs) = 3.4901575790146033 (tol=1e-06), r (rel) = 0.21912624500544478 (tol=1e-10) solvers.py:212
INFO INFO:scifem.solvers:Newton iteration 3: r (abs) = 0.023725120760079797 (tol=1e-06), r (rel) = 0.00148955928400368 (tol=1e-10) solvers.py:212
INFO INFO:scifem.solvers:Newton iteration 4: r (abs) = 0.00014918274609413277 (tol=1e-06), r (rel) = 9.366297718980792e-06 (tol=1e-10) solvers.py:212
INFO INFO:scifem.solvers:Newton iteration 5: r (abs) = 4.2205476839687285e-10 (tol=1e-06), r (rel) = 2.6498309744387176e-11 (tol=1e-10) solvers.py:212
5
if static:
dt = 0.001
else:
dt = problem.parameters["dt"].to_base_units()
times = np.arange(0.0, 1.0, dt)
ode = gotranx.load_ode("TorOrdLand.ode")
ode = ode.remove_singularities()
code = gotranx.cli.gotran2py.get_code(
ode, scheme=[gotranx.schemes.Scheme.generalized_rush_larsen], shape=gotranx.codegen.base.Shape.single,
)
Path("TorOrdLand.py").write_text(code)
import TorOrdLand
2025-03-05 19:48:44 [info ] Load ode TorOrdLand.ode
2025-03-05 19:48:45 [info ] Num states 52
2025-03-05 19:48:45 [info ] Num parameters 140
2025-03-05 19:48:48 [warning ] Cannot apply black, please install 'black'
TorOrdLand_model = TorOrdLand.__dict__
Ta_index = TorOrdLand_model["monitor_index"]("Ta")
y = TorOrdLand_model["init_state_values"]()
# Get initial parameter values
p = TorOrdLand_model["init_parameter_values"]()
import numba
fgr = numba.njit(TorOrdLand_model["generalized_rush_larsen"])
mon = numba.njit(TorOrdLand_model["monitor_values"])
V_index = TorOrdLand_model["state_index"]("v")
Ca_index = TorOrdLand_model["state_index"]("cai")
# Time in milliseconds
dt_cell = 0.1
state_file = outdir / "state.npy"
if not state_file.is_file():
@numba.jit(nopython=True)
def solve_beat(times, states, dt, p, V_index, Ca_index, Vs, Cais, Tas):
for i, ti in enumerate(times):
states[:] = fgr(states, ti, dt, p)
Vs[i] = states[V_index]
Cais[i] = states[Ca_index]
monitor = mon(ti, states, p)
Tas[i] = monitor[Ta_index]
# Time in milliseconds
nbeats = 200
T = 1000.00
times = np.arange(0, T, dt_cell)
all_times = np.arange(0, T * nbeats, dt_cell)
Vs = np.zeros(len(times) * nbeats)
Cais = np.zeros(len(times) * nbeats)
Tas = np.zeros(len(times) * nbeats)
for beat in range(nbeats):
print(f"Solving beat {beat}")
V_tmp = Vs[beat * len(times) : (beat + 1) * len(times)]
Cai_tmp = Cais[beat * len(times) : (beat + 1) * len(times)]
Ta_tmp = Tas[beat * len(times) : (beat + 1) * len(times)]
solve_beat(times, y, dt_cell, p, V_index, Ca_index, V_tmp, Cai_tmp, Ta_tmp)
fig, ax = plt.subplots(3, 2, sharex="col", sharey="row", figsize=(10, 10))
ax[0, 0].plot(all_times, Vs)
ax[1, 0].plot(all_times, Cais)
ax[2, 0].plot(all_times, Tas)
ax[0, 1].plot(times, Vs[-len(times):])
ax[1, 1].plot(times, Cais[-len(times):])
ax[2, 1].plot(times, Tas[-len(times):])
ax[0, 0].set_ylabel("V")
ax[1, 0].set_ylabel("Cai")
ax[2, 0].set_ylabel("Ta")
ax[2, 0].set_xlabel("Time [ms]")
ax[2, 1].set_xlabel("Time [ms]")
fig.savefig(outdir / "Ta_ORdLand.png")
np.save(state_file, y)
Solving beat 0
Solving beat 1
Solving beat 2
Solving beat 3
Solving beat 4
Solving beat 5
Solving beat 6
Solving beat 7
Solving beat 8
Solving beat 9
Solving beat 10
Solving beat 11
Solving beat 12
Solving beat 13
Solving beat 14
Solving beat 15
Solving beat 16
Solving beat 17
Solving beat 18
Solving beat 19
Solving beat 20
Solving beat 21
Solving beat 22
Solving beat 23
Solving beat 24
Solving beat 25
Solving beat 26
Solving beat 27
Solving beat 28
Solving beat 29
Solving beat 30
Solving beat 31
Solving beat 32
Solving beat 33
Solving beat 34
Solving beat 35
Solving beat 36
Solving beat 37
Solving beat 38
Solving beat 39
Solving beat 40
Solving beat 41
Solving beat 42
Solving beat 43
Solving beat 44
Solving beat 45
Solving beat 46
Solving beat 47
Solving beat 48
Solving beat 49
Solving beat 50
Solving beat 51
Solving beat 52
Solving beat 53
Solving beat 54
Solving beat 55
Solving beat 56
Solving beat 57
Solving beat 58
Solving beat 59
Solving beat 60
Solving beat 61
Solving beat 62
Solving beat 63
Solving beat 64
Solving beat 65
Solving beat 66
Solving beat 67
Solving beat 68
Solving beat 69
Solving beat 70
Solving beat 71
Solving beat 72
Solving beat 73
Solving beat 74
Solving beat 75
Solving beat 76
Solving beat 77
Solving beat 78
Solving beat 79
Solving beat 80
Solving beat 81
Solving beat 82
Solving beat 83
Solving beat 84
Solving beat 85
Solving beat 86
Solving beat 87
Solving beat 88
Solving beat 89
Solving beat 90
Solving beat 91
Solving beat 92
Solving beat 93
Solving beat 94
Solving beat 95
Solving beat 96
Solving beat 97
Solving beat 98
Solving beat 99
Solving beat 100
Solving beat 101
Solving beat 102
Solving beat 103
Solving beat 104
Solving beat 105
Solving beat 106
Solving beat 107
Solving beat 108
Solving beat 109
Solving beat 110
Solving beat 111
Solving beat 112
Solving beat 113
Solving beat 114
Solving beat 115
Solving beat 116
Solving beat 117
Solving beat 118
Solving beat 119
Solving beat 120
Solving beat 121
Solving beat 122
Solving beat 123
Solving beat 124
Solving beat 125
Solving beat 126
Solving beat 127
Solving beat 128
Solving beat 129
Solving beat 130
Solving beat 131
Solving beat 132
Solving beat 133
Solving beat 134
Solving beat 135
Solving beat 136
Solving beat 137
Solving beat 138
Solving beat 139
Solving beat 140
Solving beat 141
Solving beat 142
Solving beat 143
Solving beat 144
Solving beat 145
Solving beat 146
Solving beat 147
Solving beat 148
Solving beat 149
Solving beat 150
Solving beat 151
Solving beat 152
Solving beat 153
Solving beat 154
Solving beat 155
Solving beat 156
Solving beat 157
Solving beat 158
Solving beat 159
Solving beat 160
Solving beat 161
Solving beat 162
Solving beat 163
Solving beat 164
Solving beat 165
Solving beat 166
Solving beat 167
Solving beat 168
Solving beat 169
Solving beat 170
Solving beat 171
Solving beat 172
Solving beat 173
Solving beat 174
Solving beat 175
Solving beat 176
Solving beat 177
Solving beat 178
Solving beat 179
Solving beat 180
Solving beat 181
Solving beat 182
Solving beat 183
Solving beat 184
Solving beat 185
Solving beat 186
Solving beat 187
Solving beat 188
Solving beat 189
Solving beat 190
Solving beat 191
Solving beat 192
Solving beat 193
Solving beat 194
Solving beat 195
Solving beat 196
Solving beat 197
Solving beat 198
Solving beat 199

y = np.load(state_file)
class ODEState:
def __init__(self, y, dt_cell, p, t=0.0):
self.y = y
self.dt_cell = dt_cell
self.p = p
self.t = t
def forward(self, t):
for t_cell in np.arange(self.t, t, self.dt_cell):
self.y[:] = fgr(self.y, t_cell, self.dt_cell, self.p)
self.t = t
return self.y[:]
def Ta(self, t):
monitor = mon(t, self.y, p)
return monitor[Ta_index]
ode_state = ODEState(y, dt_cell, p)
@lru_cache
def get_activation(t: float):
# Find index modulo 1000
t_cell_next = t * 1000
ode_state.forward(t_cell_next)
return ode_state.Ta(t_cell_next) * 5.0
vtx = dolfinx.io.VTXWriter(geometry.mesh.comm, f"{outdir}/displacement.bp", [problem.u], engine="BP4")
vtx.write(0.0)
filename = Path("function_checkpoint.bp")
adios4dolfinx.write_mesh(filename, geometry.mesh)
def callback(model, t: float, save=True):
model.results["Ta"].append(get_activation(t))
if save:
adios4dolfinx.write_function(filename, problem.u, time=t, name="displacement")
vtx.write(t)
fig = plt.figure(layout="constrained", figsize=(12, 8))
gs = GridSpec(3, 4, figure=fig)
ax1 = fig.add_subplot(gs[:, 0])
ax2 = fig.add_subplot(gs[:, 1])
ax3 = fig.add_subplot(gs[0, 2])
ax4 = fig.add_subplot(gs[1, 2])
ax5 = fig.add_subplot(gs[0, 3])
ax6 = fig.add_subplot(gs[1, 3])
ax7 = fig.add_subplot(gs[2, 2:])
ax1.plot(model.results["V_LV"], model.results["p_LV"])
ax1.set_xlabel("LVV [mL]")
ax1.set_ylabel("LVP [mmHg]")
ax2.plot(model.results["V_RV"], model.results["p_RV"])
ax2.set_xlabel("RVV [mL]")
ax2.set_ylabel("RVP [mmHg]")
ax3.plot(model.results["time"], model.results["p_LV"])
ax3.set_ylabel("LVP [mmHg]")
ax4.plot(model.results["time"], model.results["V_LV"])
ax4.set_ylabel("LVV [mL]")
ax5.plot(model.results["time"], model.results["p_RV"])
ax5.set_ylabel("RVP [mmHg]")
ax6.plot(model.results["time"], model.results["V_RV"])
ax6.set_ylabel("RVV [mL]")
ax7.plot(model.results["time"], model.results["Ta"])
ax7.set_ylabel("Ta [kPa]")
for axi in [ax3, ax4, ax5, ax6, ax7]:
axi.set_xlabel("Time [s]")
fig.savefig(outdir / "pv_loop_incremental.png")
plt.close(fig)
def p_BiV_func(V_LV, V_RV, t):
print("Calculating pressure at time", t)
value = get_activation(t)
print("Time", t, "Activation", value)
logger.debug(f"Time{t} with activation: {value}")
Ta.assign(value)
lv_volume.value = V_LV * 1e-6
rv_volume.value = V_RV * 1e-6
problem.solve()
lv_pendo_mmHg = circulation.units.kPa_to_mmHg(problem.cavity_pressures[0].x.array[0] * 1e-3)
rv_pendo_mmHg = circulation.units.kPa_to_mmHg(problem.cavity_pressures[1].x.array[0] * 1e-3)
return lv_pendo_mmHg, rv_pendo_mmHg
mL = circulation.units.ureg("mL")
add_units = False
surface_area_lv = geometry.surface_area("ENDO_LV")
lvv_init = geo.mesh.comm.allreduce(geometry.volume("ENDO_LV", u=problem.u), op=MPI.SUM) * 1e6 * 1.0 # Increase the volume by 5%
surface_area_rv = geometry.surface_area("ENDO_RV")
rvv_init = geo.mesh.comm.allreduce(geometry.volume("ENDO_RV", u=problem.u), op=MPI.SUM) * 1e6 * 1.0 # Increase the volume by 5%
logger.info(f"Initial volume (LV): {lvv_init} mL and (RV): {rvv_init} mL")
init_state = {"V_LV": lvv_initial * 1e6 * mL, "V_RV": rvv_initial * 1e6 * mL}
[2025-03-05 19:49:54.000] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-03-05 19:49:54.000] [info] Requesting connectivity (3, 0) - (2, 0)
[2025-03-05 19:49:54.000] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-03-05 19:49:54.000] [info] Requesting connectivity (3, 0) - (2, 0)
[2025-03-05 19:49:54.002] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-03-05 19:49:54.002] [info] Requesting connectivity (3, 0) - (2, 0)
[2025-03-05 19:49:54.002] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-03-05 19:49:54.002] [info] Requesting connectivity (3, 0) - (2, 0)
[2025-03-05 19:49:54.004] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-03-05 19:49:54.004] [info] Requesting connectivity (3, 0) - (2, 0)
[2025-03-05 19:49:54.004] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-03-05 19:49:54.004] [info] Requesting connectivity (3, 0) - (2, 0)
[2025-03-05 19:49:54.006] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-03-05 19:49:54.006] [info] Requesting connectivity (3, 0) - (2, 0)
[2025-03-05 19:49:54.006] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-03-05 19:49:54.006] [info] Requesting connectivity (3, 0) - (2, 0)
[03/05/25 19:49:54] INFO INFO:pulse:Initial volume (LV): 118.5926544370604 mL and (RV): 211.73082566379117 mL 362202336.py:7
circulation_model_3D = circulation.regazzoni2020.Regazzoni2020(
add_units=add_units,
callback=callback,
p_BiV_func=p_BiV_func,
verbose=True,
comm=comm,
outdir=outdir,
initial_state=init_state,
)
# Set end time for early stopping if running in CI
end_time = 2 * dt if os.getenv("CI") else None
circulation_model_3D.solve(num_cycles=5, initial_state=init_state, dt=dt, T=end_time)
circulation_model_3D.print_info()
INFO INFO:circulation.base: base.py:131 Circulation model parameters (Regazzoni2020) ┏━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓ ┃ Parameter ┃ Value ┃ ┡━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┩ │ HR │ 1.0 hertz │ │ chambers.LA.EA │ 0.07 millimeter_Hg / milliliter │ │ chambers.LA.EB │ 0.09 millimeter_Hg / milliliter │ │ chambers.LA.TC │ 0.17 second │ │ chambers.LA.TR │ 0.17 second │ │ chambers.LA.tC │ 0.8 second │ │ chambers.LA.V0 │ 4.0 milliliter │ │ chambers.LV.EA │ 2.75 millimeter_Hg / milliliter │ │ chambers.LV.EB │ 0.08 millimeter_Hg / milliliter │ │ chambers.LV.TC │ 0.34 second │ │ chambers.LV.TR │ 0.17 second │ │ chambers.LV.tC │ 0.0 second │ │ chambers.LV.V0 │ 5.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.8 second │ │ chambers.RA.V0 │ 4.0 milliliter │ │ chambers.RV.EA │ 0.55 millimeter_Hg / milliliter │ │ chambers.RV.EB │ 0.05 millimeter_Hg / milliliter │ │ chambers.RV.TC │ 0.34 second │ │ chambers.RV.TR │ 0.17 second │ │ chambers.RV.tC │ 0.0 second │ │ chambers.RV.V0 │ 10.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.8 millimeter_Hg * second / milliliter │ │ circulation.SYS.C_AR │ 1.2 milliliter / millimeter_Hg │ │ circulation.SYS.R_VEN │ 0.26 millimeter_Hg * second / milliliter │ │ circulation.SYS.C_VEN │ 130.0 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.1625 millimeter_Hg * second / milliliter │ │ circulation.PUL.C_AR │ 10.0 milliliter / millimeter_Hg │ │ circulation.PUL.R_VEN │ 0.1625 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:137 Circulation model initial states (Regazzoni2020) ┏━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓ ┃ State ┃ Value ┃ ┡━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┩ │ V_LA │ 65.0 milliliter │ │ V_LV │ 118.59265443706046 milliliter │ │ V_RA │ 65.0 milliliter │ │ V_RV │ 211.73082566379156 milliliter │ │ p_AR_SYS │ 80.0 millimeter_Hg │ │ p_VEN_SYS │ 30.0 millimeter_Hg │ │ p_AR_PUL │ 35.0 millimeter_Hg │ │ p_VEN_PUL │ 24.0 millimeter_Hg │ │ Q_AR_SYS │ 0.0 milliliter / second │ │ Q_VEN_SYS │ 0.0 milliliter / second │ │ Q_AR_PUL │ 0.0 milliliter / second │ │ Q_VEN_PUL │ 0.0 milliliter / second │ └───────────┴───────────────────────────────┘
Calculating pressure at time 0.0
Time 0.0 Activation 0.2849891894884486
INFO INFO:scifem.solvers:Newton iteration 1: r (abs) = 105.55448698117323 (tol=1e-06), r (rel) = 1.0 (tol=1e-10) solvers.py:212
INFO INFO:scifem.solvers:Newton iteration 2: r (abs) = 1.1148799664042692 (tol=1e-06), r (rel) = 0.010562127658326073 (tol=1e-10) solvers.py:212
INFO INFO:scifem.solvers:Newton iteration 3: r (abs) = 0.1931177904790052 (tol=1e-06), r (rel) = 0.0018295554836380364 (tol=1e-10) solvers.py:212
INFO INFO:scifem.solvers:Newton iteration 4: r (abs) = 0.0017565218871915852 (tol=1e-06), r (rel) = 1.664090212957863e-05 (tol=1e-10) solvers.py:212
INFO INFO:scifem.solvers:Newton iteration 5: r (abs) = 1.0421512345037848e-06 (tol=1e-06), r (rel) = 9.873111644128057e-09 (tol=1e-10) solvers.py:212
INFO INFO:scifem.solvers:Newton iteration 6: r (abs) = 1.5357486418638033e-12 (tol=1e-06), r (rel) = 1.454934494767352e-14 (tol=1e-10) solvers.py:212
Calculating pressure at time 0.0
Time 0.0 Activation 0.2849891894884486
[03/05/25 19:49:55] INFO INFO:scifem.solvers:Newton iteration 1: r (abs) = 1.1778875770831433e-12 (tol=1e-06), r (rel) = 1.1778875770831434e-06 (tol=1e-10) solvers.py:212
INFO INFO:circulation.base: base.py:321 Volumes ┏━━━━━━━━┳━━━━━━━━━┳━━━━━━━━┳━━━━━━━━━┳━━━━━━━━━━┳━━━━━━━━━━━┳━━━━━━━━━━┳━━━━━━━━━━━┳━━━━━━━━━┳━━━━━━━━━━┳━━━━━━━━━┳━━━━━━━━━━┓ ┃ V_LA ┃ V_LV ┃ V_RA ┃ V_RV ┃ V_AR_SYS ┃ V_VEN_SYS ┃ V_AR_PUL ┃ V_VEN_PUL ┃ Heart ┃ SYS ┃ PUL ┃ Total ┃ ┡━━━━━━━━╇━━━━━━━━━╇━━━━━━━━╇━━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━━━╇━━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━╇━━━━━━━━━━┩ │ 63.838 │ 119.755 │ 63.998 │ 212.733 │ 96.000 │ 3900.000 │ 350.000 │ 384.000 │ 460.323 │ 3996.000 │ 734.000 │ 5190.323 │ └────────┴─────────┴────────┴─────────┴──────────┴───────────┴──────────┴───────────┴─────────┴──────────┴─────────┴──────────┘ Pressures ┏━━━━━━━┳━━━━━━━┳━━━━━━━┳━━━━━━━┳━━━━━━━━━━┳━━━━━━━━━━━┳━━━━━━━━━━┳━━━━━━━━━━━┓ ┃ p_LA ┃ p_LV ┃ p_RA ┃ p_RV ┃ p_AR_SYS ┃ p_VEN_SYS ┃ p_AR_PUL ┃ p_VEN_PUL ┃ ┡━━━━━━━╇━━━━━━━╇━━━━━━━╇━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━━━┩ │ 9.440 │ 0.708 │ 7.656 │ 0.126 │ 80.000 │ 30.000 │ 35.000 │ 24.000 │ └───────┴───────┴───────┴───────┴──────────┴───────────┴──────────┴───────────┘ Flows ┏━━━━━━━━━━┳━━━━━━━━┳━━━━━━━━━━┳━━━━━━━━┳━━━━━━━━━━┳━━━━━━━━━━━┳━━━━━━━━━━┳━━━━━━━━━━━┓ ┃ Q_MV ┃ Q_AV ┃ Q_TV ┃ Q_PV ┃ Q_AR_SYS ┃ Q_VEN_SYS ┃ Q_AR_PUL ┃ Q_VEN_PUL ┃ ┡━━━━━━━━━━╇━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━━━┩ │ 1162.166 │ -0.001 │ 1001.810 │ -0.000 │ 10.000 │ 44.688 │ 22.000 │ 29.120 │ └──────────┴────────┴──────────┴────────┴──────────┴───────────┴──────────┴───────────┘
Calculating pressure at time 0.001
Time 0.001 Activation 0.28501827278006825
INFO INFO:scifem.solvers:Newton iteration 1: r (abs) = 41.19892713518743 (tol=1e-06), r (rel) = 1.0 (tol=1e-10) solvers.py:212
INFO INFO:scifem.solvers:Newton iteration 2: r (abs) = 0.2551807803202981 (tol=1e-06), r (rel) = 0.006193869551092066 (tol=1e-10) solvers.py:212
INFO INFO:scifem.solvers:Newton iteration 3: r (abs) = 0.00849688126087356 (tol=1e-06), r (rel) = 0.0002062403526429816 (tol=1e-10) solvers.py:212
[03/05/25 19:49:56] INFO INFO:scifem.solvers:Newton iteration 4: r (abs) = 7.529873299314901e-07 (tol=1e-06), r (rel) = 1.827686744027793e-08 (tol=1e-10) solvers.py:212
INFO INFO:circulation.base: base.py:321 Volumes ┏━━━━━━━━┳━━━━━━━━━┳━━━━━━━━┳━━━━━━━━━┳━━━━━━━━━━┳━━━━━━━━━━━┳━━━━━━━━━━┳━━━━━━━━━━━┳━━━━━━━━━┳━━━━━━━━━━┳━━━━━━━━━┳━━━━━━━━━━┓ ┃ V_LA ┃ V_LV ┃ V_RA ┃ V_RV ┃ V_AR_SYS ┃ V_VEN_SYS ┃ V_AR_PUL ┃ V_VEN_PUL ┃ Heart ┃ SYS ┃ PUL ┃ Total ┃ ┡━━━━━━━━╇━━━━━━━━━╇━━━━━━━━╇━━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━━━╇━━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━╇━━━━━━━━━━┩ │ 62.768 │ 120.854 │ 63.081 │ 213.695 │ 95.990 │ 3899.965 │ 349.978 │ 383.993 │ 460.397 │ 3995.955 │ 733.971 │ 5190.323 │ └────────┴─────────┴────────┴─────────┴──────────┴───────────┴──────────┴───────────┴─────────┴──────────┴─────────┴──────────┘ Pressures ┏━━━━━━━┳━━━━━━━┳━━━━━━━┳━━━━━━━┳━━━━━━━━━━┳━━━━━━━━━━━┳━━━━━━━━━━┳━━━━━━━━━━━┓ ┃ p_LA ┃ p_LV ┃ p_RA ┃ p_RV ┃ p_AR_SYS ┃ p_VEN_SYS ┃ p_AR_PUL ┃ p_VEN_PUL ┃ ┡━━━━━━━╇━━━━━━━╇━━━━━━━╇━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━━━┩ │ 9.240 │ 0.978 │ 7.512 │ 0.279 │ 79.992 │ 30.000 │ 34.998 │ 24.000 │ └───────┴───────┴───────┴───────┴──────────┴───────────┴──────────┴───────────┘ Flows ┏━━━━━━━━━━┳━━━━━━━━┳━━━━━━━━━┳━━━━━━━━┳━━━━━━━━━━┳━━━━━━━━━━━┳━━━━━━━━━━┳━━━━━━━━━━━┓ ┃ Q_MV ┃ Q_AV ┃ Q_TV ┃ Q_PV ┃ Q_AR_SYS ┃ Q_VEN_SYS ┃ Q_AR_PUL ┃ Q_VEN_PUL ┃ ┡━━━━━━━━━━╇━━━━━━━━╇━━━━━━━━━╇━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━━━┩ │ 1099.372 │ -0.001 │ 962.244 │ -0.000 │ 18.400 │ 66.426 │ 36.850 │ 49.176 │ └──────────┴────────┴─────────┴────────┴──────────┴───────────┴──────────┴───────────┘
INFO INFO:circulation.base: base.py:321 Volumes ┏━━━━━━━━┳━━━━━━━━━┳━━━━━━━━┳━━━━━━━━━┳━━━━━━━━━━┳━━━━━━━━━━━┳━━━━━━━━━━┳━━━━━━━━━━━┳━━━━━━━━━┳━━━━━━━━━━┳━━━━━━━━━┳━━━━━━━━━━┓ ┃ V_LA ┃ V_LV ┃ V_RA ┃ V_RV ┃ V_AR_SYS ┃ V_VEN_SYS ┃ V_AR_PUL ┃ V_VEN_PUL ┃ Heart ┃ SYS ┃ PUL ┃ Total ┃ ┡━━━━━━━━╇━━━━━━━━━╇━━━━━━━━╇━━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━━━╇━━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━╇━━━━━━━━━━┩ │ 62.768 │ 120.854 │ 63.081 │ 213.695 │ 95.990 │ 3899.965 │ 349.978 │ 383.993 │ 460.397 │ 3995.955 │ 733.971 │ 5190.323 │ └────────┴─────────┴────────┴─────────┴──────────┴───────────┴──────────┴───────────┴─────────┴──────────┴─────────┴──────────┘ Pressures ┏━━━━━━━┳━━━━━━━┳━━━━━━━┳━━━━━━━┳━━━━━━━━━━┳━━━━━━━━━━━┳━━━━━━━━━━┳━━━━━━━━━━━┓ ┃ p_LA ┃ p_LV ┃ p_RA ┃ p_RV ┃ p_AR_SYS ┃ p_VEN_SYS ┃ p_AR_PUL ┃ p_VEN_PUL ┃ ┡━━━━━━━╇━━━━━━━╇━━━━━━━╇━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━━━┩ │ 9.240 │ 0.978 │ 7.512 │ 0.279 │ 79.992 │ 30.000 │ 34.998 │ 24.000 │ └───────┴───────┴───────┴───────┴──────────┴───────────┴──────────┴───────────┘ Flows ┏━━━━━━━━━━┳━━━━━━━━┳━━━━━━━━━┳━━━━━━━━┳━━━━━━━━━━┳━━━━━━━━━━━┳━━━━━━━━━━┳━━━━━━━━━━━┓ ┃ Q_MV ┃ Q_AV ┃ Q_TV ┃ Q_PV ┃ Q_AR_SYS ┃ Q_VEN_SYS ┃ Q_AR_PUL ┃ Q_VEN_PUL ┃ ┡━━━━━━━━━━╇━━━━━━━━╇━━━━━━━━━╇━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━━━┩ │ 1099.372 │ -0.001 │ 962.244 │ -0.000 │ 18.400 │ 66.426 │ 36.850 │ 49.176 │ └──────────┴────────┴─────────┴────────┴──────────┴───────────┴──────────┴───────────┘

Fig. 2 Pressure volume loop for the BiV.#