import inspect
import logging
import typing
from collections.abc import Mapping
from dataclasses import dataclass, field
from enum import Enum
from mpi4py import MPI
import basix
import dolfinx
import dolfinx.fem.petsc
import numpy as np
import scifem
import ufl
from packaging.version import Version
from .boundary_conditions import BoundaryConditions
from .cardiac_model import CardiacModel
from .circulation import ChamberCoupling, CirculationModel, mL, mmHg
from .geometry import HeartGeometry
from .telemetry import BaseMonitor, NullMonitor
from .units import Variable, mesh_factor
T = typing.TypeVar("T", dolfinx.fem.Function, np.ndarray)
_dolfinx_version = Version(dolfinx.__version__)
logger = logging.getLogger(__name__)
#: Coefficients of (y, y_old, y_prev) in ``dt * dy/dt``; see
#: :meth:`StaticProblem._circulation_form`.
BACKWARD_EULER_STENCIL = (1.0, -1.0, 0.0)
BDF2_STENCIL = (1.5, -2.0, 0.5)
CIRCULATION_SCHEMES = {"backward_euler", "bdf2"}
[docs]
def interpolate(x0: T, x1: T, alpha: float):
r"""Interpolate between :math:`x_0` and :math:`x_1`
to find `math:`x_{1-\alpha}`
Parameters
----------
x0 : T
First point
x1 : T
Second point
alpha : float
Amount of interpolate
Returns
-------
T
`math:`x_{1-\alpha}`
"""
return alpha * x0 + (1 - alpha) * x1
[docs]
class Geometry(typing.Protocol):
"""Protocol for geometry objects used in mechanics problems."""
dx: ufl.Measure
ds: ufl.Measure
mesh: dolfinx.mesh.Mesh
facet_tags: dolfinx.mesh.MeshTags | None
markers: dict[str, tuple[int, int]]
@property
def facet_normal(self) -> ufl.FacetNormal: ...
def surface_area(self, marker: str) -> float: ...
[docs]
class CavityControl:
"""Which constraint a cavity's pressure unknown satisfies, set at run time.
A controlled cavity has one pressure unknown ``p``, as a prescribed-volume
cavity does, but its equation is chosen by the constants held here:
- volume mode (``mode == 1``): ``V(u) = V_target``;
- pressure mode (``mode == 0``): ``p = A + B V(u)``. ``B = 0`` prescribes
the pressure, and ``B != 0`` ties the pressure to the volume inside
Newton, as a Windkessel does during ejection.
Every one of these is a `Constant` read when the form is assembled, so
switching between them does not rebuild the problem. A new control starts
in pressure mode at zero pressure. Values are in SI units: ``V_target`` in
m^3, ``A`` in Pa and ``B`` in Pa/m^3 -- which only means what it says when
the problem's own ``V(u)`` is in cubic metres, so a controlled cavity
requires ``parameters["mesh_unit"] == "m"``; `StaticProblem` refuses one
otherwise.
Do not also put a Neumann pressure on the cavity's marker: the load on the
wall comes from ``p``.
"""
def __init__(self, mesh: dolfinx.mesh.Mesh):
def constant(value: float) -> dolfinx.fem.Constant:
return dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type(value))
self.mode = constant(0.0)
self.V_target = constant(0.0)
self.A = constant(0.0)
self.B = constant(0.0)
[docs]
def set_volume(self, V: float) -> None:
"""Hold the cavity volume at `V` (m^3)."""
_assign(self.mode, 1.0)
_assign(self.V_target, V)
[docs]
def set_pressure(self, P: float) -> None:
"""Hold the cavity pressure at `P` (Pa)."""
self.set_affine_pressure(P, 0.0)
[docs]
def set_affine_pressure(self, A: float, B: float) -> None:
"""Make the cavity pressure ``A + B V(u)``, with `A` in Pa and `B` in Pa/m^3."""
_assign(self.mode, 0.0)
_assign(self.A, A)
_assign(self.B, B)
def _assign(constant: dolfinx.fem.Constant, value: float) -> None:
# `Constant.value`'s setter is typed for an array, not a scalar.
constant.value = np.asarray(value)
[docs]
def volume_scale(mesh_unit: str) -> float:
"""The factor turning a volume in ``mesh_unit``^3 into mL.
Every cavity's volume row is multiplied by it, so the row's residual is in
mL whatever the mesh unit. In m^3 a heart's volume changes by less than
`snes_atol` (1e-6, i.e. 1 mL) in a typical step, so Newton could stop
before the wall had moved at all.
"""
return mesh_factor(mesh_unit) ** 3 / mL
#: The rows of a controlled cavity measure the volume in mL and the pressure in
#: kPa, so that both modes have residuals of order one for a heart and Newton's
#: tolerance means the same thing whichever mode is active.
CONTROLLED_VOLUME_SCALE = volume_scale("m")
#: 1 kPa in pascals -- the pressure analogue of `mL` (1 mL in cubic metres)
#: above. There is no shared `kPa` constant to import for this, unlike `mL`,
#: so it is defined right here, next to the one row that uses it.
kPa = 1e3
CONTROLLED_PRESSURE_SCALE = 1 / kPa
[docs]
class Cavity(typing.NamedTuple):
"""A chamber whose pressure is an unknown of the problem, and its constraint.
Exactly one of these constrains it:
- `volume`: the volume is held at that value, with the pressure as its
Lagrange multiplier. It is usually a `Constant` you set each step.
- `control`: a :class:`CavityControl`, whose constraint (volume, pressure,
or pressure affine in the volume) can be switched at run time.
- neither, when the marker is coupled to a chamber of a 0D circulation
model (see :mod:`pulse.circulation`): the problem then points `volume` at
that chamber's volume state, so the volume is itself an unknown and the
constraint couples the two.
"""
marker: str
volume: dolfinx.fem.Constant | dolfinx.fem.Function | ufl.core.expr.Expr | None = None
control: CavityControl | None = None
[docs]
class BaseBC(str, Enum):
"""Base boundary condition"""
fixed = "fixed"
free = "free"
[docs]
@dataclass
class StaticProblem:
model: CardiacModel
geometry: Geometry
parameters: dict[str, typing.Any] = field(default_factory=dict)
bcs: BoundaryConditions = field(default_factory=BoundaryConditions)
cavities: list[Cavity] = field(default_factory=list)
circulation: CirculationModel | None = None
chambers: list[ChamberCoupling] = field(default_factory=list)
circulation_missing: dict[str, typing.Any] = field(default_factory=dict)
NonlinearProblem: typing.Type[dolfinx.fem.petsc.NonlinearProblem] = (
dolfinx.fem.petsc.NonlinearProblem
)
Function: typing.Type[dolfinx.fem.Function] = dolfinx.fem.Function
monitor: BaseMonitor = field(default_factory=NullMonitor, repr=False)
def __post_init__(self):
parameters = type(self).default_parameters()
parameters.update(self.parameters)
self.parameters = parameters
self._check_cavities()
self._init_spaces()
self._init_forms()
logger.debug("Initialized StaticProblem with parameters:")
for key, value in self.parameters.items():
logger.debug(f" {key}: {value}")
logger.debug(f"Number of cavities: {len(self.cavities)}")
if self.circulation is not None:
logger.debug(
f"Circulation states: {list(self.circulation.state_names)}",
)
logger.debug(f"Boundary conditions: {self.bcs}")
def _check_cavities(self):
"""Refuse a cavity that is not constrained exactly once, or a control the mesh_unit breaks.
Checked before anything is built, because a cavity with no constraint
would otherwise surface only as a row that fails to compile, far from
the cavity that caused it -- and a coupled cavity that also carries an
explicit volume would otherwise have that volume silently replaced by
the circulation rewrite (`_init_circulation_spaces`), rather than
refused.
"""
coupled = set()
if self.circulation is not None:
coupled = {chamber.marker for chamber in self.chambers}
for cavity in self.cavities:
is_coupled = cavity.marker in coupled
if cavity.control is not None:
if cavity.volume is not None:
raise ValueError(
f"Cavity {cavity.marker!r} has both a volume and a control. "
"Give it one: a control can hold the volume itself.",
)
if is_coupled:
raise ValueError(
f"Cavity {cavity.marker!r} has a control and is also coupled to "
"a circulation chamber, which constrains its volume already.",
)
if str(self.parameters["mesh_unit"]) != "m":
raise ValueError(
f"Cavity {cavity.marker!r} has a control, whose V_target/A/B are "
"in SI units (m^3, Pa, Pa/m^3), so it needs mesh_unit == 'm'; this "
f"problem's mesh_unit is {self.parameters['mesh_unit']!r}.",
)
elif is_coupled:
if cavity.volume is not None:
raise ValueError(
f"Cavity {cavity.marker!r} has a volume and is also coupled to a "
"circulation chamber, which would silently replace it with the "
"chamber's own volume state. Give it one: drop the volume, or "
"uncouple the chamber.",
)
elif cavity.volume is None:
raise ValueError(
f"Cavity {cavity.marker!r} has no constraint: give it a volume or a "
"control, or couple it to a chamber of a circulation model.",
)
def _init_spaces(self):
"""Initialize function spaces"""
logger.debug("Initializing function spaces...")
self._init_u_space()
self._init_p_space()
self._init_cavity_pressure_spaces()
self._init_circulation_spaces()
self._init_rigid_body()
self.update_fields()
def _init_p_space(self):
logger.debug("Initializing pressure function space...")
if self.is_incompressible:
logger.debug(
"Model is incompressible, initializing pressure space with Lagrange multiplier",
)
# Need lagrange multiplier for incompressible model
p_family, p_degree = self.parameters["p_space"].split("_")
p_element = basix.ufl.element(
family=p_family,
cell=self.geometry.mesh.basix_cell(),
degree=int(p_degree),
)
self.p_space = dolfinx.fem.functionspace(self.geometry.mesh, p_element)
self.p = self.Function(self.p_space)
self.p_old = self.Function(self.p_space)
self.p_test = ufl.TestFunction(self.p_space)
self.dp = ufl.TrialFunction(self.p_space)
else:
logger.debug("Model is compressible, no pressure space needed")
self.p_space = None
self.p_old = None
self.p = None
self.p_test = None
self.dp = None
self.model.compressibility.register(self.p)
def _init_u_space(self):
logger.debug("Initializing displacement function space...")
u_family, u_degree = self.parameters["u_space"].split("_")
logger.debug(f"Displacement space: family={u_family}, degree={u_degree}")
u_element = basix.ufl.element(
family=u_family,
cell=self.geometry.mesh.basix_cell(),
degree=int(u_degree),
shape=(self.geometry.mesh.topology.dim,),
)
self.u_space = dolfinx.fem.functionspace(self.geometry.mesh, u_element)
self.u = self.Function(self.u_space, name="u")
self.u_old = self.Function(self.u_space, name="u_old")
self.u_test = ufl.TestFunction(self.u_space)
self.du = ufl.TrialFunction(self.u_space)
self.u_full = dolfinx.fem.Function(self.u_space, name="u_full")
self.model.active.register(self.u)
@property
def is_incompressible(self):
return not self.model.compressibility.is_compressible()
@staticmethod
def default_parameters():
return {
"u_space": "P_2",
"p_space": "P_1",
"base_bc": BaseBC.free,
"rigid_body_constraint": False,
"mesh_unit": "m",
"base_marker": "BASE",
# How the circulation states are stepped; see `_circulation_form`.
"circulation_scheme": "backward_euler",
# Whether `solve` raises on non-convergence instead of returning
# False. False keeps the documented contract; set it True to get
# the older behaviour back everywhere at once, or pass
# `raise_on_failure` to a single `solve` call.
"raise_on_failure": False,
"petsc_options": {
"ksp_type": "preonly",
"pc_type": "lu",
"pc_factor_mat_solver_type": "mumps",
# Set to match `raise_on_failure` above, but not read from
# here: `solve` sets both on the solver itself every call, so
# these two are the value the solver is built with and nothing
# more. Change `raise_on_failure`, not these.
"snes_error_if_not_converged": False,
"ksp_error_if_not_converged": False,
# "snes_monitor": None,
# "ksp_monitor": None,
# "snes_linesearch_monitor": None,
"snes_type": "newtonls",
"snes_atol": 1e-6,
"snes_rtol": 1e-10,
"snes_stol": 1e-8,
"snes_max_it": 50,
# "snes_linesearch_maxstep": 50,
# "snes_view": None,
# "snes_type": "newtontr",
# "snes_type": "vinewtonrsls",
# "snes_linesearch_type": "none",
"snes_linesearch_type": "l2",
# "mat_mumps_icntl_24": 1, # Zero pivot detection
# "mat_mumps_icntl_25": 0, # Which nullspace to extract
# "mat_mumps_icntl_4": 1, # Verbosity
# "mat_mumps_icntl_2": 1, # std out
# "mat_mumps_cntl_3": 1e-6, # Threshold factor
},
}
@property
def num_cavity_pressure_states(self):
return len(self.cavities)
@property
def incompressibility_index(self) -> int:
"""Index of the incompressibility row in the block system.
The order is (u, cavity pressures, rigid body, p, circulation states),
so `p` is only the last row when there is no circulation model. Naming
the row rather than counting back from the end keeps the constraint
where it belongs when unknowns are added after it.
"""
return 1 + self.num_cavity_pressure_states + int(self.parameters["rigid_body_constraint"])
@property
def num_circulation_states(self):
if self.circulation is None:
return 0
return len(self.circulation.state_names)
def _init_circulation_spaces(self):
"""Give every circulation state its own unknown, and tie the chambers.
Each state gets one degree of freedom on the real space, the same space
the cavity pressures already use, since a circuit state is a single
global number rather than a field.
Coupling a chamber then needs no new machinery. The cavity constraint
row already reads ``pendo * (volume / area - V(u)) * ds``, which fixes
the deformed cavity volume to whatever ``volume`` says; pointing it at
the chamber's volume state instead of a prescribed constant turns that
row into the coupling, and `ufl.derivative` picks up the cross term.
"""
self.circulation_states: list[dolfinx.fem.Function] = []
self.circulation_states_old: list[dolfinx.fem.Function] = []
self.circulation_states_prev: list[dolfinx.fem.Function] = []
self.circulation_states_test: list[ufl.Argument] = []
self.circulation_states_trial: list[ufl.Argument] = []
if self.circulation is None:
return
scheme = self.parameters["circulation_scheme"]
if scheme not in CIRCULATION_SCHEMES:
raise ValueError(
f"Unknown circulation_scheme {scheme!r}, expected one of "
f"{sorted(CIRCULATION_SCHEMES)}",
)
logger.debug(f"Number of circulation states: {self.num_circulation_states}")
if getattr(self, "real_space", None) is None:
self.real_space = scifem.create_real_functionspace(self.geometry.mesh)
initial = np.asarray(self.circulation.initial_states, dtype=np.float64)
for name, value in zip(self.circulation.state_names, initial):
state = self.Function(self.real_space, name=name)
state_old = self.Function(self.real_space, name=f"{name}_old")
state_prev = self.Function(self.real_space, name=f"{name}_prev")
state.x.array[:] = value
state_old.x.array[:] = value
state_prev.x.array[:] = value
self.circulation_states.append(state)
self.circulation_states_old.append(state_old)
self.circulation_states_prev.append(state_prev)
self.circulation_states_test.append(ufl.TestFunction(self.real_space))
self.circulation_states_trial.append(ufl.TrialFunction(self.real_space))
# Coefficients of (y, y_old, y_prev) in dt * dy/dt. Constants rather
# than literals so BDF2 can take its first step as backward Euler
# without rebuilding the form.
self._circulation_stencil = tuple(
dolfinx.fem.Constant(self.geometry.mesh, dolfinx.default_scalar_type(value))
for value in BACKWARD_EULER_STENCIL
)
self._circulation_steps = -1
# The circuit carries its own time, since it is the only part of a
# static problem with a time derivative in it.
self.circulation_time = dolfinx.fem.Constant(self.geometry.mesh, 0.0)
self.circulation_dt = dolfinx.fem.Constant(self.geometry.mesh, 1.0)
# A real-space test function is constant, so integrating a residual
# against it over the mesh multiplies it by the mesh volume. Divide that
# back out, so each row is the ODE residual itself and can be compared
# against a reference implementation without carrying a stray factor.
one = dolfinx.fem.form(
dolfinx.fem.Constant(self.geometry.mesh, 1.0) * ufl.dx(domain=self.geometry.mesh),
)
volume = self.geometry.mesh.comm.allreduce(
dolfinx.fem.assemble_scalar(one),
op=MPI.SUM,
)
self._circulation_scale = dolfinx.fem.Constant(self.geometry.mesh, 1.0 / volume)
# Point each coupled cavity at its chamber's volume state. The circuit
# works in milliliters and the mechanics in cubic metres, so this is one
# of the two places a unit conversion belongs.
index = {name: i for i, name in enumerate(self.circulation.state_names)}
by_marker = {chamber.marker: chamber for chamber in self.chambers}
for i, cavity in enumerate(self.cavities):
chamber = by_marker.get(cavity.marker)
if chamber is None:
continue
if chamber.volume_state not in index:
raise KeyError(
f"Chamber {chamber.marker!r} refers to volume state "
f"{chamber.volume_state!r}, which the circulation model does not "
f"have. Known states: {list(self.circulation.state_names)}",
)
volume_state = self.circulation_states[index[chamber.volume_state]]
self.cavities[i] = Cavity(marker=cavity.marker, volume=volume_state * mL)
logger.debug(
f"Coupled cavity {cavity.marker!r} to circulation state {chamber.volume_state!r}",
)
def _circulation_missing_values(self):
"""Assemble the values the circuit expects to be supplied.
Chamber pressures come from the cavity-pressure unknowns the mechanics
problem already carries, converted from pascals to the millimeters of
mercury the circuit is written in. Anything else, such as an activation
phase, has to be provided by the caller through `circulation_missing`.
"""
assert self.circulation is not None
supplied = dict(self.circulation_missing)
cavity_of = {cavity.marker: i for i, cavity in enumerate(self.cavities)}
for chamber in self.chambers:
if chamber.marker not in cavity_of:
raise KeyError(
f"Chamber {chamber.marker!r} has no matching cavity. "
f"Known cavities: {sorted(cavity_of)}",
)
pressure = self.cavity_pressures[cavity_of[chamber.marker]]
supplied[chamber.pressure_missing] = pressure / mmHg
missing = []
for name in self.circulation.missing_names:
if name not in supplied:
raise KeyError(
f"The circulation model needs {name!r}, which is neither a coupled "
f"chamber pressure nor given in `circulation_missing`. It needs "
f"{list(self.circulation.missing_names)}.",
)
missing.append(supplied[name])
return missing
def _circulation_form(self):
"""One residual row per circuit state.
The model supplies only ``f`` in ``dy/dt = f(t, y, m)``; the scheme is
chosen here, so the same circuit can be advanced differently without
touching the model. Both schemes evaluate ``f`` at the end of the step:
``backward_euler``
``(y - y_old)/dt = f(t, y, m)``. First order, L-stable.
``bdf2``
``(3y - 4y_old + y_prev)/(2dt) = f(t, y, m)``. Second order and
still L-stable, at the same cost. The first step falls back to
backward Euler, since only one past level exists then.
Midpoint and trapezoidal rules are not offered. They are second order
on the circuit alone, but the cavity constraint ties the chamber volume
to the deformed cavity at the end of the step, so a circuit evaluated
at the midpoint is coupled half a step away from it. On an LV coupled
to Regazzoni's circuit, midpoint came out five times less accurate than
backward Euler at the same step size.
"""
if self.circulation is None:
return self._empty_form()
rhs = self.circulation.rhs(
self.circulation_time,
self.circulation_states,
self._circulation_missing_values(),
)
if len(rhs) != self.num_circulation_states:
raise ValueError(
f"Circulation model returned {len(rhs)} equations for "
f"{self.num_circulation_states} states",
)
c_new, c_old, c_prev = self._circulation_stencil
form = self._empty_form()
offset = self.num_states - self.num_circulation_states
for i, (state, state_old, state_prev, test, f) in enumerate(
zip(
self.circulation_states,
self.circulation_states_old,
self.circulation_states_prev,
self.circulation_states_test,
rhs,
),
):
derivative = (
c_new * state + c_old * state_old + c_prev * state_prev
) / self.circulation_dt
residual = derivative - f
form[offset + i] = (
self._circulation_scale * residual * test * ufl.dx(domain=self.geometry.mesh)
)
return form
def _init_cavity_pressure_spaces(self):
logger.debug("Initializing cavity pressure function spaces...")
self.cavity_pressures = []
self.cavity_pressures_test = []
self.cavity_pressures_trial = []
self.cavity_pressures_old = []
if self.num_cavity_pressure_states > 0:
logger.debug(f"Number of cavity pressure states: {self.num_cavity_pressure_states}")
self.real_space = scifem.create_real_functionspace(self.geometry.mesh)
for _ in range(self.num_cavity_pressure_states):
cavity_pressure = self.Function(self.real_space)
cavity_pressure_old = self.Function(self.real_space)
cavity_pressure_test = ufl.TestFunction(self.real_space)
cavity_pressure_trial = ufl.TrialFunction(self.real_space)
self.cavity_pressures.append(cavity_pressure)
self.cavity_pressures_old.append(cavity_pressure_old)
self.cavity_pressures_test.append(cavity_pressure_test)
self.cavity_pressures_trial.append(cavity_pressure_trial)
else:
logger.debug("No cavity pressure states needed")
def _init_rigid_body(self):
logger.debug("Initializing rigid body function space...")
if self.parameters["rigid_body_constraint"]:
logger.debug("Rigid body constraint enabled")
self.rigid_space = scifem.create_real_functionspace(
self.geometry.mesh,
value_shape=(6,),
)
self.r = self.Function(self.rigid_space)
self.r_old = self.Function(self.rigid_space)
self.dr = ufl.TrialFunction(self.rigid_space)
self.q = ufl.TestFunction(self.rigid_space)
else:
logger.debug("No rigid body constraint needed")
self.rigid_space = None
self.r = None
self.r_old = None
self.dr = None
self.q = None
def _create_residual_form(self, form: dolfinx.fem.Form) -> list[dolfinx.fem.Form]:
return [
ufl.derivative(form, f, f_test) for f, f_test in zip(self.states, self.test_functions)
]
def _empty_form(self):
return [ufl.as_ufl(0.0) for _ in range(self.num_states)]
def _rigid_body_form(self, u: dolfinx.fem.Function) -> list[dolfinx.fem.Form]:
if not self.parameters["rigid_body_constraint"]:
return self._empty_form()
logger.debug("Creating rigid body constraint form...")
X = ufl.SpatialCoordinate(self.geometry.mesh)
RM = [
ufl.as_vector((1, 0, 0)),
ufl.as_vector((0, 1, 0)),
ufl.as_vector((0, 0, 1)),
ufl.cross(X, ufl.as_vector((1, 0, 0))),
ufl.cross(X, ufl.as_vector((0, 1, 0))),
ufl.cross(X, ufl.as_vector((0, 0, 1))),
]
form = sum(ufl.inner(u, zi) * self.r[i] * ufl.dx for i, zi in enumerate(RM))
forms = self._create_residual_form(form)
return forms
# def _material_form(self, u: dolfinx.fem.Function, p: dolfinx.fem.Function):
# logger.debug("Creating material form...")
# F = ufl.grad(u) + ufl.Identity(3)
# J = ufl.det(F)
# var_C = ufl.grad(self.u_test).T * F + F.T * ufl.grad(self.u_test)
# C = ufl.variable(F.T * F)
# forms = self._empty_form()
# forms[0] += ufl.inner(self.model.S(C), 0.5 * var_C) * self.geometry.dx
# if self.is_incompressible:
# forms[-1] += (J - 1.0) * self.p_test * self.geometry.dx
# return forms
def _material_form(self, u: dolfinx.fem.Function, p: dolfinx.fem.Function):
logger.debug("Creating material form...")
I = ufl.Identity(3)
F = I + ufl.grad(u)
C = ufl.variable(F.T * F)
J = ufl.det(F)
# Automatic differentiation for the variation of C
var_C = ufl.derivative(C, u, self.u_test)
forms = self._empty_form()
# Integrate over reference configuration dX = J_map * dx
forms[0] += ufl.inner(self.model.S(C), 0.5 * var_C) * self.geometry.dx
if self.is_incompressible:
forms[self.incompressibility_index] += (J - 1.0) * self.p_test * self.geometry.dx
return forms
def _robin_form(
self,
u: dolfinx.fem.Function,
v: dolfinx.fem.Function | None = None,
) -> list[dolfinx.fem.Form]:
forms = self._empty_form()
if not self.bcs.robin:
return forms
logger.debug("Creating Robin boundary condition form...")
form = ufl.as_ufl(0.0)
N = self.geometry.facet_normal
F = ufl.grad(u) + ufl.Identity(3)
for robin in self.bcs.robin:
if robin.damping:
# Should be applied to the velocity
continue
k = robin.value.to_base_units() * mesh_factor(str(self.parameters["mesh_unit"]))
Q, ratio = robin.projection(N, F)
form += k * ufl.dot(Q * u, self.u_test) * ratio * self.geometry.ds(robin.marker)
forms[0] += form
return forms
def _neumann_form(self, u: dolfinx.fem.Function) -> list[dolfinx.fem.Form]:
forms = self._empty_form()
if not self.bcs.neumann:
return forms
logger.debug("Creating Neumann boundary condition form...")
I = ufl.Identity(3)
F = I + ufl.grad(u)
N = self.geometry.facet_normal
ds = self.geometry.ds
form = ufl.as_ufl(0.0)
for neumann in self.bcs.neumann:
t = neumann.traction.to_base_units()
n = t * ufl.det(F) * ufl.inv(F).T * N
form += ufl.inner(self.u_test, n) * ds(neumann.marker)
forms[0] += form
return forms
def _body_force_form(self, u: dolfinx.fem.Function) -> list[dolfinx.fem.Form]:
forms = self._empty_form()
if not self.bcs.body_force:
return forms
logger.debug("Creating body force form...")
form = ufl.as_ufl(0.0)
for body_force in self.bcs.body_force:
form += -ufl.derivative(ufl.inner(body_force, u) * self.geometry.dx, u, self.u_test)
forms[0] += form
return forms
def _cavity_pressure_form(
self,
u: dolfinx.fem.Function,
cavity_pressures: list[dolfinx.fem.Function] | None = None,
):
logger.debug("Creating cavity pressure form...")
if self.num_cavity_pressure_states == 0:
return self._empty_form()
if not isinstance(self.geometry, HeartGeometry):
raise RuntimeError("Cavity pressures are only supported for HeartGeometry")
V_u = self.geometry.volume_form(u)
scale = volume_scale(str(self.parameters["mesh_unit"]))
# The rows are written out rather than derived from the Lagrangian
# p (V - V(u)): the volume row is scaled to mL (see `volume_scale`)
# while the pressure stays in Pa, and with B = 0 a controlled cavity's
# pressure mode is not the stationary point of one. Differentiating
# that Lagrangian would also put a `pendo` term on the row of a
# circulation volume state, which is the chamber's own differential
# equation and must not have one. The displacement row is the
# Lagrangian's; the cavity pressures sit right after u in the block
# order.
residual = self._empty_form()
assert cavity_pressures is not None
assert len(self.cavities) == self.num_cavity_pressure_states
for i, cavity in enumerate(self.cavities):
area = self.geometry.surface_area(cavity.marker)
pendo = cavity_pressures[i]
ds = self.geometry.ds(self.geometry.markers[cavity.marker][0])
residual[0] += ufl.derivative(-pendo * V_u * ds, self.u, self.u_test)
control = cavity.control
if control is None:
# `_check_cavities` refused a cavity with neither, and a
# circuit-coupled one has had its volume set by now.
assert cavity.volume is not None
row = (cavity.volume / area - V_u) * scale
else:
# Blend the two constraints by `mode`: (V_target - V(u)) in mL
# or (A + B V(u) - p) in kPa.
volume_row = (control.V_target / area - V_u) * CONTROLLED_VOLUME_SCALE
pressure_row = ((control.A - pendo) / area + control.B * V_u) * (
CONTROLLED_PRESSURE_SCALE
)
row = control.mode * volume_row + (1.0 - control.mode) * pressure_row
residual[1 + i] += row * self.cavity_pressures_test[i] * ds
return residual
@property
def base_dirichlet(self):
bcs = []
# First add boundary conditions for the base
if self.parameters["base_bc"] == BaseBC.fixed:
marker = self.geometry.markers[self.parameters["base_marker"]][0]
base_facets = self.geometry.facet_tags.find(marker)
dofs_base = dolfinx.fem.locate_dofs_topological(self.u_space, 2, base_facets)
u_bc_base = self.Function(self.u_space)
u_bc_base.x.array[:] = 0
bcs.append(dolfinx.fem.dirichletbc(u_bc_base, dofs_base))
# Now check if we have any for bcs
for dirichlet_bc in self.bcs.dirichlet:
if callable(dirichlet_bc):
bcs += dirichlet_bc(self.u_space)
return bcs
@property
def num_states(self):
return (
1
+ self.num_cavity_pressure_states
+ int(self.parameters["rigid_body_constraint"])
+ int(self.is_incompressible)
+ self.num_circulation_states
)
@property
def R(self):
# Order is always (u, cavity pressures, rigid body, p)
R = self._empty_form()
R_material = self._material_form(self.u, p=self.p)
R_cavity = self._cavity_pressure_form(self.u, self.cavity_pressures)
R_robin = self._robin_form(self.u)
R_neumann = self._neumann_form(self.u)
R_rigid = self._rigid_body_form(self.u)
R_body_force = self._body_force_form(self.u)
R_circulation = self._circulation_form()
for i in range(self.num_states):
R[i] += R_material[i]
R[i] += R_cavity[i]
R[i] += R_robin[i]
R[i] += R_neumann[i]
R[i] += R_rigid[i]
R[i] += R_body_force[i]
R[i] += R_circulation[i]
return R
@property
def states(self):
u = [self.u]
if self.num_cavity_pressure_states > 0:
u += self.cavity_pressures
if self.parameters["rigid_body_constraint"]:
u.append(self.r)
if self.is_incompressible:
u.append(self.p)
u += self.circulation_states
return u
@property
def old_states(self):
u = [self.u_old]
if self.num_cavity_pressure_states > 0:
u += self.cavity_pressures_old
if self.parameters["rigid_body_constraint"]:
u.append(self.r_old)
if self.is_incompressible:
u.append(self.p_old)
u += self.circulation_states_old
return u
[docs]
def update_old_states(self):
"""Update old states to current values"""
logger.debug("Updating old states to current values...")
logger.debug("Updating old displacement state to current value...")
self.u_old.x.array[:] = self.u.x.array.copy()
for i in range(self.num_cavity_pressure_states):
logger.debug(f"Updating old cavity pressure state {i} to current value...")
self.cavity_pressures_old[i].x.array[:] = self.cavity_pressures[i].x.array.copy()
if self.parameters["rigid_body_constraint"]:
logger.debug("Updating old rigid body constraint state to current value...")
assert self.r is not None
assert self.r_old is not None
self.r_old.x.array[:] = self.r.x.array.copy()
if self.is_incompressible:
assert self.p is not None
assert self.p_old is not None
logger.debug("Updating old pressure state to current value...")
self.p_old.x.array[:] = self.p.x.array.copy()
for state, state_old in zip(self.circulation_states, self.circulation_states_old):
state_old.x.array[:] = state.x.array.copy()
[docs]
def reset_states(self):
"""Reset states to old values"""
logger.debug("Resetting states to old values...")
logger.debug("Resetting displacement state to old value...")
self.u.x.array[:] = self.u_old.x.array.copy()
for i in range(self.num_cavity_pressure_states):
logger.debug(f"Resetting cavity pressure state {i} to old value...")
self.cavity_pressures[i].x.array[:] = self.cavity_pressures_old[i].x.array.copy()
if self.parameters["rigid_body_constraint"]:
logger.debug("Resetting rigid body constraint state to old value...")
assert self.r is not None
assert self.r_old is not None
self.r.x.array[:] = self.r_old.x.array.copy()
if self.is_incompressible:
assert self.p is not None
assert self.p_old is not None
logger.debug("Resetting pressure state to old value...")
self.p.x.array[:] = self.p_old.x.array.copy()
for state, state_old in zip(self.circulation_states, self.circulation_states_old):
state.x.array[:] = state_old.x.array.copy()
[docs]
def restart_functions(self) -> list[tuple[str, dolfinx.fem.Function]]:
"""Every Function a restart needs, under ``mechanics_*`` names, in a fixed order.
These are the problem's own Functions, not copies: write a checkpoint
from them, or restore one by writing into them in place, then call
:meth:`load_restart_metadata`. The order is ``u``; ``p`` if
incompressible; each cavity's pressure, in ``cavities`` order; the
rigid-body multiplier ``r``; then each circulation state. Each is
followed by its ``_old`` copy, and a circulation state also by its
``_prev`` (BDF2's second level). `DynamicProblem` appends ``v_old``
and ``a_old``. With no cavity, rigid body or circulation, this is the
list the CLI wrote before it called this method, so its earlier
checkpoints still restore.
A solve overwrites the ``_old`` copies before it reads them, unless
called with ``update_old_states=False``.
"""
out = [("mechanics_u", self.u), ("mechanics_u_old", self.u_old)]
if self.is_incompressible:
out += [("mechanics_p", self.p), ("mechanics_p_old", self.p_old)]
for cavity, pressure, pressure_old in zip(
self.cavities,
self.cavity_pressures,
self.cavity_pressures_old,
):
name = f"mechanics_cavity_pressure_{cavity.marker}"
out += [(name, pressure), (f"{name}_old", pressure_old)]
if self.parameters["rigid_body_constraint"]:
out += [("mechanics_r", self.r), ("mechanics_r_old", self.r_old)]
if self.circulation is not None:
for name, state, state_old, state_prev in zip(
self.circulation.state_names,
self.circulation_states,
self.circulation_states_old,
self.circulation_states_prev,
):
key = f"mechanics_circulation_{name}"
out += [(key, state), (f"{key}_old", state_old), (f"{key}_prev", state_prev)]
return out
@property
def test_functions(self):
u = [self.u_test]
if self.num_cavity_pressure_states > 0:
u += self.cavity_pressures_test
if self.parameters["rigid_body_constraint"]:
u.append(self.q)
if self.is_incompressible:
u.append(self.p_test)
u += self.circulation_states_test
return u
@property
def trial_functions(self):
u = [self.du]
if self.num_cavity_pressure_states > 0:
u += self.cavity_pressures_trial
if self.parameters["rigid_body_constraint"]:
u.append(self.dr)
if self.is_incompressible:
u.append(self.dp)
u += self.circulation_states_trial
return u
def K(self, R):
K = []
for i in range(self.num_states):
K_row = [
ufl.derivative(R[i], f, df) for f, df in zip(self.states, self.trial_functions)
]
K.append(K_row)
return K
def _init_forms(self) -> None:
"""Initialize ufl forms"""
logger.debug("Initializing ufl forms...")
# Markers
if self.geometry.markers is None:
raise RuntimeError("Missing markers in geometry")
R = self.R
u = self.states
K = self.K(R)
assert len(R) == self.num_states
assert len(K) == self.num_states
assert all(len(Ki) == self.num_states for Ki in K)
bcs = self.base_dirichlet
logger.debug("Creating Newton solver...")
# Hack to pass petsc options to dolfinx_adjoint NonlinearProblem
kwargs = {}
if "adjoint_petsc_options" in inspect.getfullargspec(self.NonlinearProblem).kwonlyargs:
kwargs["adjoint_petsc_options"] = self.parameters["petsc_options"]
if _dolfinx_version >= Version("0.10"):
# Until we have implemented NonlinearProblem for blocked systems
# in dolfinx_adjoint, we only support single state problems here.
# Therefore we extract the first block.
if self.num_states == 1:
R = R[0]
u = u[0]
K = K[0][0]
kwargs["petsc_options_prefix"] = "pulse_problem_"
kwargs["petsc_options"] = self.parameters["petsc_options"]
self.problem = self.NonlinearProblem(
F=R,
J=K,
u=u,
bcs=bcs,
**kwargs,
)
if kwargs["petsc_options"].get("pc_type") == "bddc":
# Create the unassembled MATIS matrix required by BDDC
# Note: problem.a is the compiled bilinear form of the Jacobian
from petsc4py import PETSc
J_form = dolfinx.fem.form(K)
A_is = dolfinx.fem.petsc.create_matrix(J_form, kind=PETSc.Mat.Type.IS)
# Define a custom Python callback to assemble the
# MATIS matrix at every Newton step
def compute_jacobian(snes, x, J, P):
J.zeroEntries()
dolfinx.fem.petsc.assemble_matrix(J, J_form, bcs=bcs)
J.assemble()
# 4. Retrieve the SNES solver and swap both the matrix AND the callback
snes = self.problem.solver
snes.setJacobian(compute_jacobian, J=A_is, P=A_is)
# Apply options so PETSc registers the MATIS matrix before setup
snes.setFromOptions()
else:
petsc_options = self.parameters["petsc_options"]
# Pop options that are not supported in older dolfinx versions
petsc_options.pop("snes_error_if_not_converged", None)
petsc_options.pop("snes_type", None)
petsc_options.pop("snes_linesearch_type", None)
petsc_options.pop("snes_atol", None)
petsc_options.pop("snes_rtol", None)
petsc_options.pop("snes_stol", None)
petsc_options.pop("snes_max_it", None)
# Keep old behavior for older dolfinx versions
self._solver = scifem.NewtonSolver(
R,
K,
u,
bcs=bcs,
max_iterations=25,
petsc_options=petsc_options,
)
[docs]
def update_fields(self):
"""Shift the circuit's second history level, for BDF2.
Called only after a converged solve. `update_old_states` runs before
Newton and again on every retry, so shifting there would consume a
history level per attempt instead of per step.
`circulation_states_old` still holds the level `update_old_states` is
about to overwrite, so copying it into `_prev` here lines up
(y_{n+1}, y_n, y_{n-1}) for the next step.
"""
if self.circulation is None:
return
for prev, old in zip(self.circulation_states_prev, self.circulation_states_old):
prev.x.array[:] = old.x.array
# `_init_spaces` calls this once before the first solve, hence the
# count starting below zero: the first step has only one past level and
# must be backward Euler whichever scheme is chosen.
self._circulation_steps += 1
self._select_circulation_stencil()
def _select_circulation_stencil(self) -> None:
"""BDF2 once the circuit has a second past level, backward Euler otherwise."""
if self.parameters["circulation_scheme"] == "bdf2" and self._circulation_steps >= 1:
self._set_circulation_stencil(BDF2_STENCIL)
else:
self._set_circulation_stencil(BACKWARD_EULER_STENCIL)
def _set_circulation_stencil(self, stencil):
for constant, value in zip(self._circulation_stencil, stencil):
constant.value = value
def reset(self):
logger.debug("Resetting problem...")
self.reset_states()
self._init_forms()
[docs]
def solve(self, update_old_states: bool = True, raise_on_failure: bool | None = None) -> bool:
"""Solve nonlinear problem with Newton solver
Parameters
----------
update_old_states : bool, optional
Whether to update old states before solving, by default True
raise_on_failure : bool | None, optional
Whether a Newton solve that does not converge should raise rather
than return ``False``. ``None``, the default, defers to
``parameters["raise_on_failure"]``, which is itself ``False``.
Returns
-------
bool
True if converged, False otherwise. When this is ``False`` the
state is whatever the failed solve left behind and should not be
used; :meth:`reset_states` puts it back.
Notes
-----
This overrides ``snes_error_if_not_converged`` (and the corresponding
KSP option) on every call, so setting them through ``petsc_options`` has
no effect -- ``raise_on_failure`` is the single control.
Non-convergence is logged at warning level whichever way it is
reported, because a caller that ignores the return value would
otherwise carry on against a failed solve in silence.
"""
logger.debug("Solving the system...")
if raise_on_failure is None:
raise_on_failure = bool(self.parameters["raise_on_failure"])
if update_old_states:
self.update_old_states()
with self.monitor.track_time("newton_solve"):
if _dolfinx_version >= Version("0.10"):
solver = self.problem.solver
solver.setErrorIfNotConverged(raise_on_failure)
solver.getKSP().setErrorIfNotConverged(raise_on_failure)
try:
self.problem.solve()
finally:
# Runs even when raise_on_failure made solve() raise on
# divergence, so a failure that propagates as an exception
# is still counted rather than silently dropped.
self.monitor.record_snes(solver)
reason = typing.cast(int, solver.getConvergedReason())
converged = reason > 0
iters = solver.getIterationNumber()
logger.debug(f"Solved in {iters} iterations, converged: {converged}")
if not converged:
logger.warning(
f"Newton did not converge after {iters} iterations "
f"(SNES converged reason {reason})",
)
else:
# scifem's Newton solver returns the iteration count, not a flag,
# and raises when it gives up -- so the old code here assigned an
# int to `converged` and could only ever report success.
try:
self._solver.solve(rtol=1e-10, atol=1e-6)
except RuntimeError:
if raise_on_failure:
raise
logger.warning("Newton did not converge")
converged = False
else:
converged = True
# Derived fields are only meaningful for a solution that exists. A
# `DynamicProblem` in particular reconstructs velocity and acceleration
# from the displacement, so updating them from a failed iterate would
# write the failure into the history and outlast the rollback.
if converged:
with self.monitor.track_time("update_fields"):
self.update_fields()
return converged
[docs]
class DynamicProblem(StaticProblem):
r"""Second-order elastodynamics, via the generalized-:math:`\alpha` method.
Every term of the residual is evaluated at one of three points:
- at the :math:`\alpha_f` point (``interpolate(u_old, u, alpha_f)``,
built in :attr:`R`): the material, compressibility and viscous stress,
the Robin and Neumann loads, and the body force -- the genuine
second-order dynamics, for which alpha_f-interpolation is the
consistent choice.
- at the :math:`\alpha_m` point: the inertia term.
- at the end of the step, i.e. the true current ``self.u``: the
cavity-volume constraint, :math:`J - 1`, and the circuit. These are
algebraic constraints (Lagrange multipliers), not part of the
differential dynamics an alpha-blend applies to -- enforcing them
against the alpha_f-interpolated configuration instead would leave
``self.u``'s actual volume/incompressibility/circuit coupling
unconstrained. The active stress joins this group precisely when
:attr:`~pulse.active_model.ActiveModel.evaluate_at_end_of_step` is set
on the active model: such a model's stress depends on state advanced
over the step (a stretch rate, or ODE states), and the alpha_f point
would feed it a blended stretch and a rate scaled by
:math:`1 - \alpha_f` rather than the ones it actually advanced with.
"""
def __post_init__(self):
super().__post_init__()
# Just make sure we have units on rho and dt
for key, default_unit in [("rho", "kg/m^3"), ("dt", "s")]:
if not isinstance(self.parameters[key], Variable):
self.parameters[key] = Variable(self.parameters[key], default_unit)
def _init_u_space(self):
super()._init_u_space()
self.u_old = self.Function(self.u_space)
self.v_old = self.Function(self.u_space)
self.a_old = self.Function(self.u_space)
def _init_p_space(self):
super()._init_p_space()
if self.is_incompressible:
self.p_old = self.Function(self.p_space)
else:
self.p_old = None
def _init_cavity_pressure_spaces(self):
super()._init_cavity_pressure_spaces()
self.cavity_pressures_old = []
for _ in range(self.num_cavity_pressure_states):
cavity_pressure_old = self.Function(self.real_space)
self.cavity_pressures_old.append(cavity_pressure_old)
def _material_form(self, u, v, p):
F = ufl.grad(u) + ufl.Identity(3)
C = ufl.variable(F.T * F)
F_dot = ufl.grad(v)
C_dot = ufl.variable(F_dot.T * F + F.T * F_dot)
var_C = ufl.grad(self.u_test).T * F + F.T * ufl.grad(self.u_test)
forms = self._empty_form()
if self.model.active.evaluate_at_end_of_step:
# The active model is stateful/rate-dependent: assemble
# everything else at the alpha_f point as usual, but the active
# stress at the true end-of-step displacement self.u, built the
# same way the J - 1 row below builds F_true.
forms[0] += (
ufl.inner(self.model.S(C, C_dot=C_dot, active=False), 0.5 * var_C)
* self.geometry.dx
)
F1 = ufl.grad(self.u) + ufl.Identity(3)
C1 = ufl.variable(F1.T * F1)
var_C1 = ufl.grad(self.u_test).T * F1 + F1.T * ufl.grad(self.u_test)
forms[0] += ufl.inner(self.model.active.S(C1), 0.5 * var_C1) * self.geometry.dx
else:
forms[0] += ufl.inner(self.model.S(C, C_dot=C_dot), 0.5 * var_C) * self.geometry.dx
if self.is_incompressible:
# Incompressibility is an algebraic constraint (like the
# cavity-volume Lagrange multiplier in _cavity_pressure_form):
# enforce it against the true current self.u, not the
# alpha_f-interpolated u, or self.u's actual volume is left
# under-constrained and p picks up spurious high-frequency
# oscillation even under smooth forcing.
F_true = ufl.grad(self.u) + ufl.Identity(3)
J_true = ufl.det(F_true)
forms[self.incompressibility_index] += (J_true - 1.0) * self.p_test * self.geometry.dx
return forms
def _acceleration_form(self, a: dolfinx.fem.Function):
rho = self.parameters["rho"].to_base_units()
forms = self._empty_form()
forms[0] += ufl.inner(rho * a, self.u_test) * self.geometry.dx
return forms
def _robin_form(self, u: dolfinx.fem.Function, v: dolfinx.fem.Function | None = None):
forms = super()._robin_form(u)
N = self.geometry.facet_normal
F = ufl.grad(u) + ufl.Identity(3)
for robin in self.bcs.robin:
if not robin.damping:
# Should be applied to the velocity
continue
assert v is not None
c = robin.value.to_base_units() * mesh_factor(str(self.parameters["mesh_unit"]))
Q, ratio = robin.projection(N, F)
forms[0] += c * ufl.dot(Q * v, self.u_test) * ratio * self.geometry.ds(robin.marker)
return forms
@property
def R(self):
# Order is always (u, cavity pressures, rigid body, p)
alpha_m = self.parameters["alpha_m"]
alpha_f = self.parameters["alpha_f"]
a_new = self.a(u=self.u, u_old=self.u_old, v_old=self.v_old, a_old=self.a_old)
v_new = self.v(a=a_new, v_old=self.v_old, a_old=self.a_old)
u = interpolate(self.u_old, self.u, alpha_f)
v = interpolate(self.v_old, v_new, alpha_f)
a = interpolate(self.a_old, a_new, alpha_m)
R = self._empty_form()
R_material = self._material_form(u=u, v=v, p=self.p)
# The cavity-volume constraint is algebraic (a Lagrange multiplier),
# not part of the generalized-alpha differential dynamics: enforcing
# it against the alpha_f-interpolated u instead of the true current
# self.u leaves self.u's actual volume unconstrained and excites a
# spurious high-frequency oscillation in the multiplier (cavity
# pressure) even under smooth volume forcing.
R_cavity = self._cavity_pressure_form(self.u, self.cavity_pressures)
R_robin = self._robin_form(u=u, v=v)
R_neumann = self._neumann_form(u)
R_rigid = self._rigid_body_form(u)
R_body_force = self._body_force_form(u)
R_acceleration = self._acceleration_form(a)
# Evaluated at the end of the step, like the cavity constraint, rather
# than at the alpha_f point. Without these rows every circuit state
# gets a zero row, which does not compile into a form.
R_circulation = self._circulation_form()
for i in range(self.num_states):
R[i] += R_material[i]
R[i] += R_cavity[i]
R[i] += R_robin[i]
R[i] += R_neumann[i]
R[i] += R_rigid[i]
R[i] += R_body_force[i]
R[i] += R_acceleration[i]
R[i] += R_circulation[i]
return R
@staticmethod
def default_parameters():
parameters = StaticProblem.default_parameters()
parameters.update(
{
"dt": Variable(1e-3, "s"),
"rho": Variable(1e3, "kg/m^3"),
"alpha_m": 0.2,
"alpha_f": 0.4,
},
)
return parameters
def _dt_float(self) -> float:
"""The current time step in seconds as a plain float (for numpy updates)."""
dt = self.parameters["dt"]
value = dt.value
if isinstance(value, dolfinx.fem.Constant):
value = float(value.value)
return float(value) * dt.factor
[docs]
def v(
self,
a: T,
v_old: T,
a_old: T,
dt: typing.Any = None,
) -> T:
r"""
Velocity computed using the generalized
:math:`alpha`-method
.. math::
v_{i+1} = v_i + (1-\gamma) \Delta t a_i + \gamma \Delta t a_{i+1}
Parameters
----------
a : T
Current acceleration
v_old : T
Previous velocity
a_old: T
Previous acceleration
dt: optional
Time step; defaults to ``parameters["dt"]`` in base units (a UFL
expression when ``dt`` wraps a Constant). Pass a float for numpy input.
Returns
-------
T
The current velocity
"""
if dt is None:
dt = self.parameters["dt"].to_base_units()
return v_old + (1 - self._gamma) * dt * a_old + self._gamma * dt * a
[docs]
def a(
self,
u: T,
u_old: T,
v_old: T,
a_old: T,
dt: typing.Any = None,
) -> T:
r"""
Acceleration computed using the generalized
:math:`alpha`-method
.. math::
a_{i+1} = \frac{u_{i+1} - (u_i + \Delta t v_i +
(0.5 - \beta) \Delta t^2 a_i)}{\beta \Delta t^2}
Parameters
----------
u : T
Current displacement
u_old : T
Previous displacement
v_old : T
Previous velocity
a_old: T
Previous acceleration
dt: optional
Time step; see :meth:`v`.
Returns
-------
T
The current acceleration
"""
if dt is None:
dt = self.parameters["dt"].to_base_units()
dt2 = dt**2
beta = self._beta
return (u - (u_old + dt * v_old + (0.5 - beta) * dt2 * a_old)) / (beta * dt2)
[docs]
def update_fields(self) -> None:
"""Update old values of displacement, velocity
and acceleration
"""
super().update_fields()
dt = self._dt_float()
u = self.u.x.array.copy()
u_old = self.u_old.x.array.copy()
v_old = self.v_old.x.array.copy()
a_old = self.a_old.x.array.copy()
a = self.a(u=u, u_old=u_old, v_old=v_old, a_old=a_old, dt=dt)
v = self.v(a=a, v_old=v_old, a_old=a_old, dt=dt)
self.a_old.x.array[:] = a
self.v_old.x.array[:] = v
self.u_old.x.array[:] = u
# self.update_base_values()
for i in range(self.num_cavity_pressure_states):
self.cavity_pressures_old[i].x.array[:] = self.cavity_pressures[i].x.array.copy()
[docs]
def restart_functions(self) -> list[tuple[str, dolfinx.fem.Function]]:
"""`StaticProblem.restart_functions`, then the velocity and acceleration history."""
return [
*super().restart_functions(),
("mechanics_v_old", self.v_old),
("mechanics_a_old", self.a_old),
]
@property
def _gamma(self) -> float:
"""Parameter in the generalized alpha-method"""
alpha_m = self.parameters["alpha_m"]
alpha_f = self.parameters["alpha_f"]
assert isinstance(alpha_m, float)
assert isinstance(alpha_f, float)
return 0.5 + alpha_f - alpha_m
@property
def _beta(self) -> float:
"""Parameter in the generalized alpha-method"""
return (self._gamma + 0.5) ** 2 / 4.0
