BiV ellipsoid with time dependent pressure and activation

BiV ellipsoid with time dependent pressure and activation#

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. The reader are referred to the LV example for more details on the models 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 fenicsx_pulse
logging.basicConfig(level=logging.INFO)
comm = MPI.COMM_WORLD
outdir = Path("time-dependent-bestel-biv")
outdir.mkdir(parents=True, exist_ok=True)
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-03-05 19:48:08 [debug    ] Convert file time-dependent-bestel-biv/geometry/biv_ellipsoid.msh to dolfin
Info    : Reading 'time-dependent-bestel-biv/geometry/biv_ellipsoid.msh'...
Info    : 43 entities
Info    : 203 nodes
Info    : 982 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: 
lv: [2]
rv: [3]
epi: [4]
base: [1]
INFO:ldrb.ldrb:  Num vertices: 203
INFO:ldrb.ldrb:  Num cells: 638
INFO:ldrb.ldrb:  Apex coord: (4.00, 0.50, -0.00)
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
0.0 1.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[:] *= 3e-2

We create the geometry object and print the volumes of the LV and RV cavities

geometry = fenicsx_pulse.HeartGeometry.from_cardiac_geometries(geo, metadata={"quadrature_degree": 6})
print(geometry.volume("ENDO_LV") * 1e6, geometry.volume("ENDO_RV") * 1e6)
118.59265443706046 211.73082566379156
material_params = fenicsx_pulse.HolzapfelOgden.orthotropic_parameters()
material = fenicsx_pulse.HolzapfelOgden(f0=geo.f0, s0=geo.s0, **material_params)  # type: ignore
Ta = fenicsx_pulse.Variable(dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(0.0)), "Pa")
active_model = fenicsx_pulse.ActiveStress(geo.f0, activation=Ta)
comp_model = fenicsx_pulse.compressibility.Compressible2()
viscoeleastic_model = fenicsx_pulse.viscoelasticity.Viscous()
model = fenicsx_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 = fenicsx_pulse.Variable(dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(0.0)), "Pa")
neumann_lv = fenicsx_pulse.NeumannBC(traction=traction_lv, marker=geometry.markers["ENDO_LV"][0])
traction_rv = fenicsx_pulse.Variable(dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(0.0)), "Pa")
neumann_rv = fenicsx_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 = fenicsx_pulse.Variable(
    dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(1e8)), "Pa / m",
)
robin_epi_u = fenicsx_pulse.RobinBC(value=alpha_epi, marker=geometry.markers["EPI"][0])
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)
alpha_base = fenicsx_pulse.Variable(
    dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(1e5)), "Pa / m",
)
robin_base_u = fenicsx_pulse.RobinBC(value=alpha_base, marker=geometry.markers["BASE"][0])
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(neumann=(neumann_lv, neumann_rv), robin=(robin_epi_u, robin_epi_v, robin_base_u, robin_base_v))
problem = fenicsx_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.5195246029017035e-05 (tol=1e-06), r (rel) = 1.0 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 2: r (abs) = 1.136043401379393e-07 (tol=1e-06), r (rel) = 0.007476308045358331 (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")
../_images/47754ed7dedb34e3a6dc72eb8c538811ac65f3eef84855adfc40ed2d0ad24adf.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-03-05 19:48:26.558] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-03-05 19:48:26.558] [info] Requesting connectivity (3, 0) - (2, 0)
[2025-03-05 19:48:26.558] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-03-05 19:48:26.558] [info] Requesting connectivity (3, 0) - (2, 0)
[2025-03-05 19:48:26.559] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-03-05 19:48:26.559] [info] Requesting connectivity (3, 0) - (2, 0)
[2025-03-05 19:48:26.559] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-03-05 19:48:26.559] [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.9944426002021674e-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) = 8.8774790379695e-09 (tol=1e-06), r (rel) = 0.00044511078118114963 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 1: r (abs) = 1.9225907390067896e-05 (tol=1e-06), r (rel) = 1.0 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 2: r (abs) = 4.100435169929981e-08 (tol=1e-06), r (rel) = 0.002132765484997741 (tol=1e-10)
Solving for time 0.001, activation 0.0, lv pressure 11.94019949358465 and rv pressure 2.9716796494482693
INFO:scifem.solvers:Newton iteration 1: r (abs) = 2.1490707186081514e-05 (tol=1e-06), r (rel) = 1.0 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 2: r (abs) = 2.8773634311896605e-08 (tol=1e-06), r (rel) = 0.0013388872717288657 (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.0668874395184335e-05 (tol=1e-06), r (rel) = 1.0 (tol=1e-10)
INFO:scifem.solvers:Newton iteration 2: r (abs) = 3.827562049608271e-08 (tol=1e-06), r (rel) = 0.0012480282126719588 (tol=1e-10)
Solving for time 0.003, activation 0.0, lv pressure 35.46523956820262 and rv pressure 8.7481138929159