Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension


Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 3 additions & 2 deletions .github/workflows/test_package_coverage.yml
Original file line number Diff line number Diff line change
Expand Up @@ -32,8 +32,9 @@ jobs:
- name: Install development versions of dependencies
if: ${{ matrix.container == 'ghcr.io/fenics/dolfinx/dolfinx:nightly' }}
run: |
python3 -m pip install git+https://github.com/jorgensd/adios4dolfinx.git
python3 -m pip install git+https://github.com/scientificcomputing/scifem.git
python3 -m pip install git+https://github.com/scientificcomputing/io4dolfinx.git
python3 -m pip install git+https://github.com/scientificcomputing/scifem.git --no-build-isolation
python3 -m pip install git+https://github.com/ComputationalPhysiology/cardiac-geometriesx.git

- name: Install package
run: |
Expand Down
5 changes: 0 additions & 5 deletions demo/benchmark/problem1.py
Original file line number Diff line number Diff line change
Expand Up @@ -27,7 +27,6 @@
from mpi4py import MPI
import numpy as np
import dolfinx
from dolfinx import log
import pulse

# ## 1. Geometry and Mesh
Expand Down Expand Up @@ -119,8 +118,6 @@ def dirichlet_bc(V: dolfinx.fem.FunctionSpace):

problem = pulse.StaticProblem(model=model, geometry=geo, bcs=bcs)

log.set_log_level(log.LogLevel.INFO)

target_pressure = 0.004
steps = [0.0005, 0.001, 0.002, 0.003, 0.004]

Expand All @@ -129,8 +126,6 @@ def dirichlet_bc(V: dolfinx.fem.FunctionSpace):
pressure.assign(p)
problem.solve()

log.set_log_level(log.LogLevel.WARNING)

# ## 5. Post-processing
#
# We save the final displacement and evaluate the deflection at the specific target point $(10, 0.5, 1)$.
Expand Down
3 changes: 0 additions & 3 deletions demo/benchmark/problem2.py
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,6 @@

from pathlib import Path
from mpi4py import MPI
from dolfinx import log
import logging
import dolfinx
import numpy as np
Expand Down Expand Up @@ -114,8 +113,6 @@
model=model, geometry=geometry, bcs=bcs, parameters={"base_bc": pulse.BaseBC.fixed},
)

log.set_log_level(log.LogLevel.INFO)

# ### Continuation Solver
# We ramp the pressure up to 10 kPa. We use a custom continuation loop to handle the nonlinearity.

Expand Down
3 changes: 0 additions & 3 deletions demo/benchmark/problem3.py
Original file line number Diff line number Diff line change
Expand Up @@ -30,15 +30,12 @@

from pathlib import Path
from mpi4py import MPI
from dolfinx import log
import dolfinx
import numpy as np
import math
import cardiac_geometries
import pulse

# log.set_log_level(log.LogLevel.INFO)

# ## 1. Geometry and Fibers
# We regenerate the mesh with the specific fiber angles required for Problem 3 ($+90^\circ$ to $-90^\circ$).

Expand Down
76 changes: 52 additions & 24 deletions demo/bleeding_biv.py
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,7 @@
import matplotlib.pyplot as plt
import numpy as np
import gotranx
import adios4dolfinx
import io4dolfinx
import cardiac_geometries
import cardiac_geometries.geometry
import pulse
Expand All @@ -30,6 +30,7 @@
def mmHg_to_kPa(x):
return x * 0.133322


def custom_json(obj):
if isinstance(obj, np.float64):
return float(obj)
Expand Down Expand Up @@ -141,8 +142,6 @@ def run_zenker(outdir: Path):
history["before_index"] = before_index
history["after_index"] = after_index



Path(zenker_file).write_text(json.dumps(history, indent=4, default=custom_json))

return json.loads(zenker_file.read_text())
Expand Down Expand Up @@ -202,7 +201,6 @@ def solve_beat(times, states, dt, p, V_index, Ca_index, Vs, Cais, Tas):
monitor = mon(ti, states, p)
Tas[i] = monitor[Ta_index]


nbeats = 10 if os.environ.get("CI") else 200
times = np.arange(0, T, dt_cell)
all_times = np.arange(0, T * nbeats, dt_cell)
Expand Down Expand Up @@ -261,7 +259,6 @@ def run_3D_model(
mesh_unit: str = "m",
volume2ml: float = 1000.0,
):

geometry = pulse.HeartGeometry.from_cardiac_geometries(
geo,
metadata={"quadrature_degree": 6},
Expand All @@ -274,7 +271,8 @@ def run_3D_model(
# We use an active stress approach with 30% transverse active stress (see {py:meth}`pulse.active_stress.transversely_active_stress`)

Ta = pulse.Variable(
dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(0.0)), "kPa",
dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(0.0)),
"kPa",
)
active_model = pulse.ActiveStress(geo.f0, activation=Ta)

Expand Down Expand Up @@ -312,11 +310,15 @@ def run_3D_model(
robin = [robin_epi, robin_epi_perp, robin_base]

lvv_initial = comm.allreduce(geometry.volume("LV"), op=MPI.SUM)
lv_volume = dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(lvv_initial))
lv_volume = dolfinx.fem.Constant(
geometry.mesh, dolfinx.default_scalar_type(lvv_initial),
)
lv_cavity = pulse.problem.Cavity(marker="LV", volume=lv_volume)

rvv_initial = comm.allreduce(geometry.volume("RV"), op=MPI.SUM)
rv_volume = dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(rvv_initial))
rv_volume = dolfinx.fem.Constant(
geometry.mesh, dolfinx.default_scalar_type(rvv_initial),
)
rv_cavity = pulse.problem.Cavity(marker="RV", volume=rv_volume)

print("Initial volumes", lvv_initial * volume2ml, rvv_initial * volume2ml)
Expand All @@ -325,7 +327,11 @@ def run_3D_model(
parameters = {"base_bc": pulse.problem.BaseBC.free, "mesh_unit": mesh_unit}
bcs = pulse.BoundaryConditions(robin=robin)
problem = pulse.problem.StaticProblem(
model=model, geometry=geometry, bcs=bcs, cavities=cavities, parameters=parameters,
model=model,
geometry=geometry,
bcs=bcs,
cavities=cavities,
parameters=parameters,
)

problem.solve()
Expand All @@ -339,7 +345,7 @@ def run_3D_model(
vtx.write(0.0)

filename = outdir / Path(f"function_checkpoint_{label}.bp")
adios4dolfinx.write_mesh(filename, geometry.mesh)
io4dolfinx.write_mesh(filename, geometry.mesh)
output_file = outdir / f"output_{label}.json"

Ta_history = []
Expand All @@ -348,9 +354,9 @@ def callback(model, i: int, t: float, save=True):
Ta_history.append(get_activation(t))

if save:
adios4dolfinx.write_function(filename, problem.u, time=t, name="displacement")
io4dolfinx.write_function(filename, problem.u, time=t, name="displacement")
vtx.write(t)
out = {k: v[:i+1] for k, v in model.history.items()}
out = {k: v[: i + 1] for k, v in model.history.items()}
out["Ta"] = Ta_history
if comm.rank == 0:
output_file.write_text(json.dumps(out, indent=4, default=custom_json))
Expand Down Expand Up @@ -396,7 +402,9 @@ def callback(model, i: int, t: float, save=True):
def p_BiV_func(V_LV, V_RV, t):
print("Calculating pressure at time", t)
value = get_activation(t)
print(f"Time{t} with activation: {value} and volumes: {V_LV} mL (LV) {V_RV} mL (RV)")
print(
f"Time{t} with activation: {value} and volumes: {V_LV} mL (LV) {V_RV} mL (RV)",
)
old_Ta = Ta.value.value
dTa = value - old_Ta

Expand All @@ -422,8 +430,12 @@ def p_BiV_func(V_LV, V_RV, t):
Ta.assign(value)
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)
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,
)
print(f"Compute pressures: {lv_pendo_mmHg} mmHg (LV) {rv_pendo_mmHg} mmHg (RV)")
return lv_pendo_mmHg, rv_pendo_mmHg

Expand All @@ -432,19 +444,23 @@ def p_BiV_func(V_LV, V_RV, t):
s = circulation.units.ureg("s")
add_units = False
lvv_init = (
geo.mesh.comm.allreduce(geometry.volume("LV", u=problem.u), op=MPI.SUM) * 1e6 * 1.0
geo.mesh.comm.allreduce(geometry.volume("LV", u=problem.u), op=MPI.SUM)
* 1e6
* 1.0
)
rvv_init = (
geo.mesh.comm.allreduce(geometry.volume("RV", u=problem.u), op=MPI.SUM) * 1e6 * 1.0
geo.mesh.comm.allreduce(geometry.volume("RV", u=problem.u), op=MPI.SUM)
* 1e6
* 1.0
)
print(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,
'p_AR_PUL': p_AR_PUL * mmHg,
'p_AR_SYS': p_AR_SYS * mmHg,
'p_VEN_PUL': p_VEN_PUL * mmHg,
'p_VEN_SYS': p_VEN_SYS * mmHg,
"p_AR_PUL": p_AR_PUL * mmHg,
"p_AR_SYS": p_AR_SYS * mmHg,
"p_VEN_PUL": p_VEN_PUL * mmHg,
"p_VEN_SYS": p_VEN_SYS * mmHg,
}

regazzoni_parmeters = circulation.regazzoni2020.Regazzoni2020.default_parameters()
Expand Down Expand Up @@ -475,7 +491,9 @@ def p_BiV_func(V_LV, V_RV, t):
)
# Set end time for early stopping if running in CI
end_time = 2 * dt if os.getenv("CI") else None
regazzoni.solve(num_beats=num_beats, initial_state=init_state, dt=dt, T=end_time) #, checkpoint=RR)
regazzoni.solve(
num_beats=num_beats, initial_state=init_state, dt=dt, T=end_time,
) # , checkpoint=RR)
regazzoni.print_info()


Expand Down Expand Up @@ -534,10 +552,20 @@ def p_BiV_func(V_LV, V_RV, t):
print(f"HR before: {HR_before}, HR after: {HR_after}")

get_activation_before = run_TorOrdLand(
comm, outdir, HR_before, Ta_factor=1, label="before", dt=dt,
comm,
outdir,
HR_before,
Ta_factor=1,
label="before",
dt=dt,
)
get_activation_after = run_TorOrdLand(
comm, outdir, HR_after, Ta_factor=C_PRSW_factor, label="after", dt=dt,
comm,
outdir,
HR_after,
Ta_factor=C_PRSW_factor,
label="after",
dt=dt,
)

run_3D_model(
Expand Down
2 changes: 0 additions & 2 deletions demo/boundary_conditions/lv_ellipsoid_fixed_x.py
Original file line number Diff line number Diff line change
Expand Up @@ -131,8 +131,6 @@ def dirichlet_bc(V: dolfinx.fem.FunctionSpace):
bcs=bcs,
)

log.set_log_level(log.LogLevel.INFO)

# ### Phase 1: Passive Inflation

vtx = dolfinx.io.VTXWriter(geometry.mesh.comm, outdir / "lv_displacement.bp", [problem.u], engine="BP4")
Expand Down
3 changes: 0 additions & 3 deletions demo/boundary_conditions/ukb_bcs.py
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,6 @@
from pathlib import Path
from mpi4py import MPI
import dolfinx
from dolfinx import log
import cardiac_geometries
import cardiac_geometries.geometry
import pulse
Expand Down Expand Up @@ -192,8 +191,6 @@ def dirichlet_bc(V: dolfinx.fem.FunctionSpace):
bcs=bcs,
)

log.set_log_level(log.LogLevel.INFO)

# ### Phase 1: Passive Inflation
# We ramp up the pressure in both ventricles.

Expand Down
5 changes: 0 additions & 5 deletions demo/geometries/biv_ellipsoid.py
Original file line number Diff line number Diff line change
Expand Up @@ -35,16 +35,11 @@
from mpi4py import MPI
import numpy as np
import dolfinx
from dolfinx import log
import ldrb
import cardiac_geometries
import cardiac_geometries.geometry
import pulse

# We enable info logging to track solver progress.

# log.set_log_level(log.LogLevel.INFO)

# ## Geometry and Microstructure
#
# We generate the mesh and fibers if they don't already exist.
Expand Down
Loading