FEniCSx-Beat

Contents

FEniCSx-Beat#

class beat.BaseMonitor[source]#
class beat.BidomainModel(time: Constant, mesh: Mesh, M_i: Coefficient | float, M_e: Coefficient | float, I_s: Stimulus | Sequence[Stimulus] | Expr | None = None, params: dict[str, Any] | None = None, C_m: float = 1.0, dx: Measure | None = None, **kwargs)[source]#

Solve the bidomain model

\[\begin{split}C_m \frac{\partial v}{\partial t} - \nabla \cdot (M_i \nabla v) - \nabla \cdot (M_i \nabla u_e) - I_{\mathrm{stim}} &= 0 \\ \nabla \cdot (M_i \nabla v) + \nabla \cdot ((M_i + M_e) \nabla u_e) &= 0\end{split}\]

for the transmembrane potential \(v\) and the extracellular potential \(u_e\), with insulating boundary conditions on both.

The two fields live in separate Lagrange spaces and the system is assembled blocked and solved monolithically. Keeping v in a space of its own – rather than making the pair one mixed element – is what lets the ODE solvers consume it unchanged.

\(u_e\) is only determined up to a constant. By default it is normalised to zero mean with a Lagrange multiplier in a real-valued space; supplying bcs that constrain \(u_e\) grounds it instead and the multiplier is dropped, since the two together would over-determine the problem.

Parameters:
property state: BidomainState#

Both solved-for fields, in the order the blocked system holds them.

step(interval) → None[source]#

Perform a single time step.

Parameters:

interval (tuple[float, float]) – The time interval (T0, T)

property u_e: Function#

The extracellular potential.

property v: Function#

The transmembrane potential.

variational_forms(dt: Expr | float) → tuple[Form | Sequence[Sequence[Form | None]], Form | Sequence[Form]][source]#

Blocked theta-rule forms for the coupled system.

Both rows are written in terms of the midpoint pair. The elliptic row carries no time derivative, but it is linear and homogeneous, so imposing it on the midpoint is equivalent to imposing it at the new time given that it held at the old one – and it is the only choice under which the collapse to the monodomain model at equal anisotropy ratio stays exact away from theta = 1.

Parameters:

dt (Expr | float) – The time step

Returns:

The blocked bilinear and linear forms

Return type:

tuple[_BilinearForm, _LinearForm]

class beat.ECGRecovery(v: dolfinx.fem.function.Function, sigma_b: float | dolfinx.fem.function.Constant = 1.0, C_m: float | dolfinx.fem.function.Constant = 1.0, dx: ufl.measure.Measure | None = None, M: float = 1.0, petsc_options: dict[str, typing.Any] = <factory>)[source]#
class beat.ExternalOperatorMultiODESolver(v_ode: Function, v_pde: Function, markers: Function, fun: dict[int, Callable], init_states: dict[int, ndarray], num_states: dict[int, int], v_index: dict[int, int], parameters: dict[int, ndarray] | None = None, monitor: BaseMonitor | None = None)[source]#

Multi-region counterpart to ExternalOperatorODESolver, analogous to how beat.odesolver.DolfinMultiODESolver extends beat.odesolver.DolfinODESolver and beat.irksome_odesolver. IrksomeMultiODESolver extends ExternalOperatorODESolver’s Irksome sibling: a different cell model (fun, init_states, parameters, num_states) can be used in each region of the mesh, as identified by an integer-valued markers function.

As in beat.irksome_odesolver.IrksomeMultiODESolver, one full-mesh external operator is built per marker, and only the degrees of freedom that actually belong to that marker are read back out; the cost of a step therefore scales with the number of regions. Prefer beat.odesolver.DolfinMultiODESolver for problems with many regions or very large meshes.

from_dolfin() → None[source]#

Move the voltage variable from v_ode back into the per-marker state functions.

load_all_states(functions) → None[source]#

Inverse of assign_all_states(): overwrite every ODE state from functions (one per state, on v_ode’s function space), e.g. when restarting from a checkpoint.

to_dolfin() → None[source]#

Move the voltage variable from the per-marker state functions into v_ode.

class beat.ExternalOperatorODESolver(v_ode: Function, v_pde: Function, fun: Callable, init_states: ndarray, num_states: int, v_index: int = 0, parameters: ndarray | None = None, monitor: BaseMonitor | None = None)[source]#

An ODE solver that steps a cell model via dolfinx-external-operator’s FEMExternalOperator, instead of the plain-numpy loop used by beat.odesolver.DolfinODESolver.

Unlike beat.irksome_odesolver.IrksomeODESolver, fun keeps the exact same signature already used by beat.odesolver.DolfinODESolver, fun(states, t, parameters, dt) -> new_states (plain NumPy or Numba-jitted arrays) — no UFL-returning variant of the cell model is needed. All num_states state variables are stored in a single vector-valued (blocked) dolfinx.fem.Function, rather than a mixed-element one, since FEMExternalOperator operates on raw (num_cells, num_quadrature_points, num_states) arrays rather than UFL expressions.

from_dolfin() → None[source]#

Move the voltage variable from v_ode back into the ODE state function.

load_all_states(functions) → None[source]#

Inverse of assign_all_states(): overwrite every ODE state from functions (one per state, on v_ode’s function space), e.g. when restarting from a checkpoint.

to_dolfin() → None[source]#

Move the voltage variable from the ODE state function into v_ode.

class beat.IrksomeMonodomainModel(time: Constant, mesh: Mesh, M: Coefficient | float, butcher_tableau, I_s: Stimulus | Sequence[Stimulus] | Expr | None = None, params=None, C_m: float = 1.0, dx: Measure | None = None, **kwargs)[source]#

Solve Monodomain model using Irksome for Runge-Kutta time stepping.

solve(interval, dt=None) → Results[source]#

Solve on the given time interval.

Parameters:
  • interval (tuple[float, float]) – The time interval (T0, T)

  • dt (float, optional) – The time step, by default None

Returns:

The results of the solution

Return type:

Results

property state: Function#

The transmembrane potential.

Kept for compatibility; v says the same thing by name and is preferred.

step(interval)[source]#

Perform a single time step.

Parameters:

interval (tuple[float, float]) – The time interval (T0, T)

class beat.IrksomeMultiODESolver(v_ode: Function, v_pde: Function, markers: Function, fun: dict[int, Callable], init_states: dict[int, ndarray], butcher_tableau, time: Constant, num_states: dict[int, int], v_index: dict[int, int], parameters: dict[int, ndarray] | None = None, monitor: BaseMonitor | None = None)[source]#

Multi-region counterpart to IrksomeODESolver, analogous to how beat.odesolver.DolfinMultiODESolver extends beat.odesolver.DolfinODESolver: a different cell model (fun, init_states, parameters, num_states) can be used in each region of the mesh, as identified by an integer-valued markers function.

All regions share the same butcher_tableau and time/dt stepping. Internally, one full-mesh mixed-element Irksome stepper is built per marker (as in IrksomeODESolver), and only the degrees of freedom that actually belong to that marker are read back out; the values computed at the other degrees of freedom are discarded. This keeps the implementation simple and correct, but means the cost of a step scales with the number of regions (each region’s stepper solves over the whole mesh), unlike beat.odesolver.DolfinMultiODESolver, whose plain-numpy ODE systems are only ever sized to the points that need them. Prefer beat.odesolver.DolfinMultiODESolver for problems with many regions or very large meshes.

from_dolfin() → None[source]#

Move the voltage variable from v_ode back into the per-marker mixed ODE spaces.

load_all_states(functions) → None[source]#

Inverse of assign_all_states(): overwrite every ODE state from functions (one per state, on v_ode’s function space), e.g. when restarting from a checkpoint.

to_dolfin() → None[source]#

Move the voltage variable from the per-marker mixed ODE spaces into v_ode.

class beat.IrksomeODESolver(v_ode: Function, v_pde: Function, fun: Callable, init_states: ndarray, butcher_tableau, time: Constant, num_states: int, v_index: int = 0, parameters: dict | None = None, monitor: BaseMonitor | None = None)[source]#

An ODE Solver mapping the DolfinODESolver interface to an Irksome stepper.

from_dolfin() → None[source]#

Move the voltage variable from v_ode back into the mixed ODE space.

load_all_states(functions) → None[source]#

Inverse of assign_all_states(): overwrite every ODE state from functions (one per state, on v_ode’s function space), e.g. when restarting from a checkpoint.

to_dolfin() → None[source]#

Move the voltage variable from the mixed ODE space into v_ode.

class beat.MonodomainModel(time: Constant, mesh: Mesh, M: Coefficient | float, I_s: Stimulus | Sequence[Stimulus] | Expr | None = None, params=None, C_m: float = 1.0, dx: Measure | None = None, **kwargs)[source]#

Solve

\[\frac{\partial V}{\partial t} - \nabla \cdot \left( M \nabla V \right) - I_{\mathrm{stim}} = 0\]
property state: Function#

The transmembrane potential.

Kept for compatibility; v says the same thing by name and is preferred.

property v: Function#

The transmembrane potential.

variational_forms(dt: Expr | float) → tuple[Form, Form][source]#

Create the variational forms corresponding to the given discretization of the given system of equations.

Parameters:

dt (float) – Time step size

Returns:

The variational form and the precondition

Return type:

tuple[ufl.Form, ufl.Form]

class beat.MonodomainSplittingSolver(pde: beat.monodomain_solver.PDEModel, ode: beat.monodomain_solver.ODESolver, theta: float = 1.0, monitor: beat.telemetry.BaseMonitor = <factory>)[source]#
load_restart(functions: Sequence[tuple[str, Function]], metadata: Mapping[str, Any] | None = None) → None[source]#

Restore the solver from functions as returned by restart_functions().

Overwrites pde.v and every ODE state, restores the PDE’s metadata if given (as returned by restart_metadata()), then sets the PDE’s previous value, exactly as construction does.

restart_functions(state_names: Sequence[str]) → list[tuple[str, Function]][source]#

The Functions that hold everything the splitting scheme carries between steps.

These are ("v", pde.v) followed by ("state_<name>", f) for each ODE state, named by state_names (one per state). Write them to save a checkpoint. On a fresh solver, read into them and pass them to load_restart().

restart_metadata() → dict[str, Any][source]#

State beyond the Functions, to be stored with a checkpoint (JSON-serializable).

Pass it back to load_restart(). It is empty if the PDE model has none.

class beat.NullMonitor[source]#
class beat.PerformanceMonitor(log_frequency: int = 1, comm: ~mpi4py.MPI.Intracomm = <mpi4py.MPI.Intracomm object>)[source]#

A monitor that accumulates timings and KSP stats, logging them periodically, and can save/display a final summary.

display_summary() → None[source]#

Logs a nicely formatted summary of the accumulated timings and metrics.

save_summary(filepath: str | Path) → None[source]#

Saves the performance metrics to a JSON file.

class beat.Stimulus(expr, dZ, marker)[source]#
dZ: Measure#

Alias for field number 1

expr: Expr#

Alias for field number 0

marker: int | None#

Alias for field number 2

class beat.base_model.BaseModel(time: Constant, mesh: Mesh, dx: Measure | None = None, params: dict[str, Any] | None = None, I_s: Stimulus | Sequence[Stimulus] | Expr | None = None, bcs: Callable[[BaseModel], Sequence[DirichletBC]] | None = None, monitor: BaseMonitor | None = None, **kwargs: Any)[source]#

Base class for models.

Parameters:
  • time (dolfinx.fem.Constant) – The current time

  • mesh (dolfinx.mesh.Mesh) – The mesh

  • dx (ufl.Measure, optional) – The measure for the spatial domain, by default None

  • params (dict, optional) – Parameters for the model, by default None

  • I_s (Stimulus | Sequence[Stimulus] | ufl.core.expr.Expr, optional) – The stimulus, by default None

  • bcs (Callable[[BaseModel], Sequence[dolfinx.fem.DirichletBC]], optional) – Factory returning the Dirichlet boundary conditions, by default None (natural boundary conditions). It is called once the state space exists and is handed the model, because a boundary condition is only honoured on the very function space object the model built – one constructed separately on the same mesh and element is silently ignored during assembly.

  • jit_options (dict, optional) – JIT options, by default None

  • form_compiler_options (dict, optional) – Form compiler options, by default None

  • petsc_options (dict, optional) – PETSc options, by default None

load_restart_metadata(data: Mapping[str, Any]) → None[source]#

Restore what restart_metadata() returned, re-assembling the matrix.

property mesh: Mesh#

The mesh the model is discretized on.

restart_metadata() → dict[str, float][source]#

The model’s own state beyond its fields: the time step its matrix was assembled for.

The matrix is kept while dt stays within 1e-12 of this, so a restarted run only reproduces the uninterrupted one bit for bit if it assembles at the same time step.

solve(interval: tuple[float, float], dt: float | None = None) → Results[source]#

Solve on the given time interval.

Parameters:
  • interval (tuple[float, float]) – The time interval (T0, T)

  • dt (float, optional) – The time step, by default None

Returns:

The results of the solution

Return type:

Results

step(interval)[source]#

Perform a single time step.

Parameters:

interval (tuple[float, float]) – The time interval (T0, T)

abstract variational_forms(dt: Expr | float) → tuple[Form | Sequence[Sequence[Form | None]], Form | Sequence[Form]][source]#

Create the variational forms corresponding to the given discretization of the given system of equations.

Parameters:

dt (Expr | float) – The time step

Returns:

The bilinear and linear form. A model with several coupled fields returns these blocked, as a nested list of forms and a list of forms respectively.

Return type:

tuple[_BilinearForm, _LinearForm]

class beat.base_model.Results(state, status)[source]#
state: Function | Sequence[Function]#

Alias for field number 0

status: Status#

Alias for field number 1

class beat.base_model.Status(value, names=<not given>, *values, module=None, qualname=None, type=None, start=1, boundary=None)[source]#
class beat.monodomain_model.MonodomainModel(time: Constant, mesh: Mesh, M: Coefficient | float, I_s: Stimulus | Sequence[Stimulus] | Expr | None = None, params=None, C_m: float = 1.0, dx: Measure | None = None, **kwargs)[source]#

Solve

\[\frac{\partial V}{\partial t} - \nabla \cdot \left( M \nabla V \right) - I_{\mathrm{stim}} = 0\]
property state: Function#

The transmembrane potential.

Kept for compatibility; v says the same thing by name and is preferred.

property v: Function#

The transmembrane potential.

variational_forms(dt: Expr | float) → tuple[Form, Form][source]#

Create the variational forms corresponding to the given discretization of the given system of equations.

Parameters:

dt (float) – Time step size

Returns:

The variational form and the precondition

Return type:

tuple[ufl.Form, ufl.Form]

class beat.monodomain_solver.MonodomainSplittingSolver(pde: beat.monodomain_solver.PDEModel, ode: beat.monodomain_solver.ODESolver, theta: float = 1.0, monitor: beat.telemetry.BaseMonitor = <factory>)[source]#
load_restart(functions: Sequence[tuple[str, Function]], metadata: Mapping[str, Any] | None = None) → None[source]#

Restore the solver from functions as returned by restart_functions().

Overwrites pde.v and every ODE state, restores the PDE’s metadata if given (as returned by restart_metadata()), then sets the PDE’s previous value, exactly as construction does.

restart_functions(state_names: Sequence[str]) → list[tuple[str, Function]][source]#

The Functions that hold everything the splitting scheme carries between steps.

These are ("v", pde.v) followed by ("state_<name>", f) for each ODE state, named by state_names (one per state). Write them to save a checkpoint. On a fresh solver, read into them and pass them to load_restart().

restart_metadata() → dict[str, Any][source]#

State beyond the Functions, to be stored with a checkpoint (JSON-serializable).

Pass it back to load_restart(). It is empty if the PDE model has none.

class beat.monodomain_solver.ODESolver(*args, **kwargs)[source]#
class beat.monodomain_solver.PDEModel(*args, **kwargs)[source]#

What the splitting driver asks of a PDE model.

It only ever advances the model and asks it to keep its history; the transmembrane potential is shared with the ODE solver directly. So any model offering these is drivable here, whatever else it solves for alongside.

property v: Function#

The transmembrane potential.

class beat.odesolver.BaseDolfinODESolver[source]#
abstract load_all_states(functions: list[Function]) → None[source]#

Inverse of assign_all_states(): overwrite every ODE state from functions (one per state, on v_ode’s function space), e.g. when restarting from a checkpoint.

ode_to_pde() → None[source]#

Projects v_ode (DG0, quadrature space, …) into v_pde (CG1)

pde_to_ode() → None[source]#

Projects v_pde (CG1) into v_ode (DG0, quadrature space, …)

class beat.odesolver.DolfinMultiODESolver(v_ode: 'dolfinx.fem.Function', v_pde: 'dolfinx.fem.Function', markers: 'dolfinx.fem.Function', init_states: 'dict[int, npt.NDArray]', parameters: 'dict[int, npt.NDArray]', fun: 'dict[int, Callable]', num_states: 'dict[int, int]', v_index: 'dict[int, int]', monitor: 'BaseMonitor' = <factory>)[source]#
from_dolfin() → None[source]#

Assign values from dolfinx function to numpy array

load_all_states(functions: list[Function]) → None[source]#

Inverse of assign_all_states(): overwrite every ODE state from functions (one per state, on v_ode’s function space), e.g. when restarting from a checkpoint.

to_dolfin() → None[source]#

Assign values from numpy array to dolfinx function

class beat.odesolver.DolfinODESolver(v_ode: 'dolfinx.fem.Function', v_pde: 'dolfinx.fem.Function', init_states: 'npt.NDArray', parameters: 'npt.NDArray', fun: 'Callable', num_states: 'int', v_index: 'int' = 0, missing_variables: 'npt.NDArray | None' = None, num_missing_variables: 'int' = 0, monitor: 'BaseMonitor' = <factory>)[source]#
from_dolfin() → None[source]#

Assign values from dolfin function to numpy array

load_all_states(functions: list[Function]) → None[source]#

Inverse of assign_all_states(): overwrite every ODE state from functions (one per state, on v_ode’s function space), e.g. when restarting from a checkpoint.

to_dolfin() → None[source]#

Assign values from numpy array to dolfin function

class beat.odesolver.ODEResults(y, t)[source]#
t: ndarray[tuple[Any, ...], dtype[float64]]#

Alias for field number 1

y: ndarray[tuple[Any, ...], dtype[float64]]#

Alias for field number 0

class beat.odesolver.ODESystemSolver(fun: 'Callable', states: 'npt.NDArray', parameters: 'npt.NDArray', missing_variables: 'npt.NDArray | None' = None, _kwargs: 'dict[str, npt.NDArray]' = <factory>, monitor: 'BaseMonitor' = <factory>)[source]#

Compiled C backend for cell models generated by gotranx

beat advances the cell model through a function fun(states, t, parameters, dt) that acts on a (num_states, num_points) array. This module compiles the C code that gotranx generates for one cell (gotranx ode2c, or gotranx.codegen.CCodeGenerator) together with a small driver that loops over the points, and wraps the library in a CScheme with that signature, so that it can be passed as fun to beat.odesolver.DolfinODESolver and beat.odesolver.DolfinMultiODESolver.

Rank 0 compiles the library into a cache directory, keyed by a hash of the source, the compiler and the flags, and the other ranks load it from there. The cache directory must therefore be on a file system that all ranks can see.

class beat.c_backend.CScheme(library_path: Path, scheme: str, *, num_states: int | None = None, num_parameters: int | None = None, openmp: bool = False, num_threads: int | None = None)[source]#

A compiled cell-model scheme, callable as fun(states, t, parameters, dt)

Parameters:
  • library_path (Path) – The shared library built by compile_scheme()

  • scheme (str) – The name of the per-cell C function

  • num_states (int | None, optional) – The number of states of the model; if given, states is checked against it

  • num_parameters (int | None, optional) – The number of parameters of the model; if given, parameters is checked against it

  • openmp (bool, optional) – Whether the library was compiled with OpenMP

  • num_threads (int | None, optional) – The number of OpenMP threads, by default the OpenMP default

init_parameter_values() → ndarray[tuple[Any, ...], dtype[float64]][source]#

The default parameters of the model

init_state_values() → ndarray[tuple[Any, ...], dtype[float64]][source]#

The default initial states of the model

max_threads() → int[source]#

The number of threads that one call uses (1 without OpenMP)

parameter_index(name: str) → int[source]#

The index of the parameter name

state_index(name: str) → int[source]#

The index of the state name

beat.c_backend.build_library(source: str, *, cc: str, cflags: tuple[str, ...], cache_dir: Path, comm: Intracomm) → Path[source]#

Compile source on rank 0 (unless it is cached) and return the library on all ranks

beat.c_backend.cache_key(source: str, cc: str, cc_version: str, flags: Sequence[str]) → str[source]#

Hash of everything that changes the compiled library

beat.c_backend.compile_scheme(c_code: str, scheme: str = 'generalized_rush_larsen', *, num_states: int | None = None, num_parameters: int | None = None, cache_dir: str | ~pathlib.Path | None = None, cc: str | None = None, cflags: ~collections.abc.Sequence[str] = ('-O3', '-march=native'), openmp: bool = False, num_threads: int | None = None, comm: ~mpi4py.MPI.Intracomm = <mpi4py.MPI.Intracomm object>) → CScheme[source]#

Compile gotranx-generated C code into a fun(states, t, parameters, dt)

Parameters:
  • c_code (str) – C code with the per-cell function scheme, as generated by gotranx

  • scheme (str, optional) – The name of the per-cell function, by default “generalized_rush_larsen”

  • num_states (int | None, optional) – The number of states, used to check the arrays that the scheme is called with. The generated C code has no state or parameter counts, so without num_states and num_parameters the arrays are not checked against the model, and arrays that are too small make the C code read out of bounds. from_ode() fills them in.

  • num_parameters (int | None, optional) – The number of parameters, used to check the arrays that the scheme is called with (see num_states)

  • cache_dir (str | Path | None, optional) – Where to put the library, by default default_cache_dir(). Must be visible to all ranks

  • cc (str | None, optional) – The C compiler, by default default_compiler()

  • cflags (Sequence[str], optional) – The compiler flags, by default ("-O3", "-march=native")

  • openmp (bool, optional) – Compile with -fopenmp and run the point loop in parallel, by default False

  • num_threads (int | None, optional) – The number of OpenMP threads, by default OMP_NUM_THREADS or the OpenMP default

  • comm (MPI.Intracomm, optional) – The communicator; rank 0 compiles, by default MPI.COMM_WORLD

Returns:

The compiled scheme

Return type:

CScheme

beat.c_backend.compiler_version(cc: str) → str[source]#

The first line of cc --version

beat.c_backend.default_cache_dir() → Path[source]#

$BEAT_C_CACHE if it is set, else ~/.cache/beat/c_backend

beat.c_backend.default_compiler() → str[source]#

$CC if it is set, else cc

beat.c_backend.from_ode(ode: ~typing.Any, scheme: str = 'generalized_rush_larsen', *, codegen_kwargs: dict[str, ~typing.Any] | None = None, comm: ~mpi4py.MPI.Intracomm = <mpi4py.MPI.Intracomm object>, **compile_kwargs: ~typing.Any) → CScheme[source]#

Generate C for a gotranx ODE and compile it with compile_scheme()

The code is generated on rank 0 only. compile_kwargs are passed on to compile_scheme() (cache_dir, cc, cflags, …).

Parameters:
  • ode (gotranx.ode.ODE) – The cell model

  • scheme (str, optional) – The name of a gotranx scheme, by default “generalized_rush_larsen”

  • codegen_kwargs (dict[str, Any] | None, optional) – Extra keyword arguments for gotranx.codegen.CCodeGenerator.scheme()

  • comm (MPI.Intracomm, optional) – The communicator, by default MPI.COMM_WORLD

beat.c_backend.full_source(c_code: str, scheme: str) → str[source]#

The generated C code followed by the driver that loops over the points

beat.c_backend.generate_c_code(ode: Any, scheme: str = 'generalized_rush_larsen', codegen_kwargs: dict[str, Any] | None = None) → str[source]#

The C code for one cell, with index and initial value helpers and scheme

Parameters:
  • ode (gotranx.ode.ODE) – The cell model

  • scheme (str, optional) – The name of a gotranx scheme, by default “generalized_rush_larsen”

  • codegen_kwargs (dict[str, Any] | None, optional) – Extra keyword arguments for gotranx.codegen.CCodeGenerator.scheme()

beat.c_backend.openmp_available(cc: str | None = None) → bool[source]#

Whether the C compiler can build an OpenMP shared library

beat.c_backend.target_signature(cc: str, flags: Sequence[str]) → str[source]#

The compiler’s predefined macros for flags if they target the native CPU, else “”

With -march=native (or -mcpu=native) the library depends on the features of the CPU it is compiled on, which platform.machine() does not tell apart (e.g. with and without AVX-512), so the macros that the compiler defines for these flags go into the cache key.

class beat.irksome_model.IrksomeMonodomainModel(time: Constant, mesh: Mesh, M: Coefficient | float, butcher_tableau, I_s: Stimulus | Sequence[Stimulus] | Expr | None = None, params=None, C_m: float = 1.0, dx: Measure | None = None, **kwargs)[source]#

Solve Monodomain model using Irksome for Runge-Kutta time stepping.

solve(interval, dt=None) → Results[source]#

Solve on the given time interval.

Parameters:
  • interval (tuple[float, float]) – The time interval (T0, T)

  • dt (float, optional) – The time step, by default None

Returns:

The results of the solution

Return type:

Results

property state: Function#

The transmembrane potential.

Kept for compatibility; v says the same thing by name and is preferred.

step(interval)[source]#

Perform a single time step.

Parameters:

interval (tuple[float, float]) – The time interval (T0, T)

class beat.irksome_odesolver.IrksomeMultiODESolver(v_ode: Function, v_pde: Function, markers: Function, fun: dict[int, Callable], init_states: dict[int, ndarray], butcher_tableau, time: Constant, num_states: dict[int, int], v_index: dict[int, int], parameters: dict[int, ndarray] | None = None, monitor: BaseMonitor | None = None)[source]#

Multi-region counterpart to IrksomeODESolver, analogous to how beat.odesolver.DolfinMultiODESolver extends beat.odesolver.DolfinODESolver: a different cell model (fun, init_states, parameters, num_states) can be used in each region of the mesh, as identified by an integer-valued markers function.

All regions share the same butcher_tableau and time/dt stepping. Internally, one full-mesh mixed-element Irksome stepper is built per marker (as in IrksomeODESolver), and only the degrees of freedom that actually belong to that marker are read back out; the values computed at the other degrees of freedom are discarded. This keeps the implementation simple and correct, but means the cost of a step scales with the number of regions (each region’s stepper solves over the whole mesh), unlike beat.odesolver.DolfinMultiODESolver, whose plain-numpy ODE systems are only ever sized to the points that need them. Prefer beat.odesolver.DolfinMultiODESolver for problems with many regions or very large meshes.

from_dolfin() → None[source]#

Move the voltage variable from v_ode back into the per-marker mixed ODE spaces.

load_all_states(functions) → None[source]#

Inverse of assign_all_states(): overwrite every ODE state from functions (one per state, on v_ode’s function space), e.g. when restarting from a checkpoint.

to_dolfin() → None[source]#

Move the voltage variable from the per-marker mixed ODE spaces into v_ode.

class beat.irksome_odesolver.IrksomeODESolver(v_ode: Function, v_pde: Function, fun: Callable, init_states: ndarray, butcher_tableau, time: Constant, num_states: int, v_index: int = 0, parameters: dict | None = None, monitor: BaseMonitor | None = None)[source]#

An ODE Solver mapping the DolfinODESolver interface to an Irksome stepper.

from_dolfin() → None[source]#

Move the voltage variable from v_ode back into the mixed ODE space.

load_all_states(functions) → None[source]#

Inverse of assign_all_states(): overwrite every ODE state from functions (one per state, on v_ode’s function space), e.g. when restarting from a checkpoint.

to_dolfin() → None[source]#

Move the voltage variable from the mixed ODE space into v_ode.

class beat.external_operator_odesolver.ExternalOperatorMultiODESolver(v_ode: Function, v_pde: Function, markers: Function, fun: dict[int, Callable], init_states: dict[int, ndarray], num_states: dict[int, int], v_index: dict[int, int], parameters: dict[int, ndarray] | None = None, monitor: BaseMonitor | None = None)[source]#

Multi-region counterpart to ExternalOperatorODESolver, analogous to how beat.odesolver.DolfinMultiODESolver extends beat.odesolver.DolfinODESolver and beat.irksome_odesolver. IrksomeMultiODESolver extends ExternalOperatorODESolver’s Irksome sibling: a different cell model (fun, init_states, parameters, num_states) can be used in each region of the mesh, as identified by an integer-valued markers function.

As in beat.irksome_odesolver.IrksomeMultiODESolver, one full-mesh external operator is built per marker, and only the degrees of freedom that actually belong to that marker are read back out; the cost of a step therefore scales with the number of regions. Prefer beat.odesolver.DolfinMultiODESolver for problems with many regions or very large meshes.

from_dolfin() → None[source]#

Move the voltage variable from v_ode back into the per-marker state functions.

load_all_states(functions) → None[source]#

Inverse of assign_all_states(): overwrite every ODE state from functions (one per state, on v_ode’s function space), e.g. when restarting from a checkpoint.

to_dolfin() → None[source]#

Move the voltage variable from the per-marker state functions into v_ode.

class beat.external_operator_odesolver.ExternalOperatorODESolver(v_ode: Function, v_pde: Function, fun: Callable, init_states: ndarray, num_states: int, v_index: int = 0, parameters: ndarray | None = None, monitor: BaseMonitor | None = None)[source]#

An ODE solver that steps a cell model via dolfinx-external-operator’s FEMExternalOperator, instead of the plain-numpy loop used by beat.odesolver.DolfinODESolver.

Unlike beat.irksome_odesolver.IrksomeODESolver, fun keeps the exact same signature already used by beat.odesolver.DolfinODESolver, fun(states, t, parameters, dt) -> new_states (plain NumPy or Numba-jitted arrays) — no UFL-returning variant of the cell model is needed. All num_states state variables are stored in a single vector-valued (blocked) dolfinx.fem.Function, rather than a mixed-element one, since FEMExternalOperator operates on raw (num_cells, num_quadrature_points, num_states) arrays rather than UFL expressions.

from_dolfin() → None[source]#

Move the voltage variable from v_ode back into the ODE state function.

load_all_states(functions) → None[source]#

Inverse of assign_all_states(): overwrite every ODE state from functions (one per state, on v_ode’s function space), e.g. when restarting from a checkpoint.

to_dolfin() → None[source]#

Move the voltage variable from the ODE state function into v_ode.

beat.utils.cpp_function_space(V)[source]#

The C++ function space behind V, whichever wrapper V arrives in.

Some dolfinx APIs return the Python FunctionSpace and others the C++ space it wraps; DirichletBC.function_space has returned each of them in different releases. Comparisons such as contains are only defined between two C++ spaces, so bring both sides down to that one shape rather than assuming either.

Parameters:

V (dolfinx.fem.FunctionSpace | dolfinx.cpp.fem.FunctionSpace) – The function space, in either wrapper

Returns:

The C++ function space

Return type:

dolfinx.cpp.fem.FunctionSpace

beat.utils.expand_layer(V: FunctionSpace, ft: MeshTags, endo_marker: int, epi_marker: int, endo_size: float, epi_size: float, output_mid_marker: int = 0, output_endo_marker: int = 1, output_epi_marker: int = 2) → Function[source]#

Expand the endo and epi markers to the rest of the mesh with a given size

Parameters:
  • V (dolfinx.fem.FunctionSpace) – The function space for your return function

  • ft (dolfinx.mesh.MeshTags) – The facet tags

  • endo_marker (int) – The endocardium marker

  • epi_marker (int) – The epicardium marker

  • endo_size (float) – The endocardium size

  • epi_size (float) – The epicardium size

  • output_mid_marker (int, optional) – The marker to be set for the midwall, by default 0

  • output_endo_marker (int, optional) – The marker to be set for the endocardium, by default 1

  • output_epi_marker (int, optional) – The marker to be set for the epicardium, by default 2

Returns:

The expanded markers as a function in V

Return type:

dolfinx.fem.Function

beat.utils.expand_layer_biv(V: FunctionSpace, ft: MeshTags, endo_lv_marker: int, endo_rv_marker: int, epi_marker: int, endo_size: float, epi_size: float, output_mid_marker: int = 0, output_endo_marker: int = 1, output_epi_marker: int = 2) → Function[source]#

Expand the endo and epi markers to the rest of the mesh with a given size

Parameters:
  • V (dolfinx.fem.FunctionSpace) – The function space for your return function

  • ft (dolfinx.mesh.MeshTags) – The facet tags

  • endo_lv_marker (int) – The LV endocardium marker

  • endo_rv_marker (int) – The LV endocardium marker

  • epi_marker (int) – The epicardium marker

  • endo_size (float) – The endocardium size

  • epi_size (float) – The epicardium size

  • output_mid_marker (int, optional) – The marker to be set for the midwall, by default 0

  • output_endo_marker (int, optional) – The marker to be set for the endocardium, by default 1

  • output_epi_marker (int, optional) – The marker to be set for the epicardium, by default 2

Returns:

The expanded markers as a function in V

Return type:

dolfinx.fem.Function

beat.utils.local_project(v: Function, V: FunctionSpace, u: Function | None = None) → Function | None[source]#

Element-wise projection using LocalSolver

Parameters:
Returns:

The projected function

Return type:

dolfinx.fem.Function | None

beat.utils.parse_element(space_string: str, mesh: Mesh, dim: int) → _ElementBase[source]#

Parse a string representation of a basix element family

beat.utils.real_function_space(mesh: Mesh) → FunctionSpace[source]#

A scalar space holding a single global degree of freedom.

dolfinx supports real spaces natively from 0.11; before that, basix.ufl.real_element exists but cannot be compiled into a function space, so fall back to scifem’s implementation.

Parameters:

mesh (dolfinx.mesh.Mesh) – The mesh the space is defined on

Returns:

The real function space

Return type:

dolfinx.fem.FunctionSpace

beat.utils.space_from_string(space_string: str, mesh: Mesh, dim: int = 1) → FunctionSpace[source]#

Constructed a finite elements space from a string representation of the space

Parameters:
  • space_string (str) – A string on the form {family}_{degree} which determines the space. Example ‘Lagrange_1’.

  • mesh (df.Mesh) – The mesh

  • dim (int) – 1 for scalar space, 3 for vector space.

Returns:

The function space

Return type:

df.FunctionSpace