FlashBang.jl

FlashBang

Cardiac tissue electrophysiology on orthogonal, structured grids.

The pipeline mirrors Thunderbolt's by convention — same names, same shapes, no dependency:

model = MonodomainModel(; κ, ion)                       # continuous problem
split = ReactionDiffusionSplit(model)                   # how to split it
f     = semidiscretize(split, FiniteDifferenceDiscretization(), grid)
prob  = OperatorSplittingProblem(f, create_initial_condition(f), tspan)
integrator = init(prob, LieTrotterGodunov((Euler(), Euler())); dt)

Spatial operators come from MatrixFreeOperators.jl and cell kinetics from CytoZoo.jl. Restricting the domain to orthogonal structured grids is what buys the matrix-free operators; unstructured meshes are Thunderbolt's job, not this package's.

See also: MonodomainModel, semidiscretize, create_initial_condition.

AbstractEPModel

Supertype for cardiac electrophysiology tissue models. A model is a declaration of the continuous problem — coefficients, cell model, stimulation — carrying no discretization. semidiscretize turns one into a right-hand side on a given grid.

AbstractStimulationProtocol

Supertype for externally applied stimulation currents.

A protocol is a callable (x, t) -> Iₐₚₚ returning the applied current density at physical position x (an SVector of cell-center coordinates) and time t.

Sign convention

Iₐₚₚ is in the PDE-form convention: positive is depolarizing. It enters the monodomain equation as the χIₛₜᵢₘ term of

χCₘ∂ₜφₘ = ∇⋅κ∇φₘ + χ(Iᵢₒₙ + Iₛₜᵢₘ)

so a positive value drives φₘ up.

Stimulation belongs to the diffusion half of the reaction-diffusion split, uniformly in every dimension. The cell model's own stimulus must therefore be off; MonodomainModel warns when it detects one that is not.

See also: NoStimulationProtocol, TransmembraneStimulationProtocol.

DiffusionFunction(prepared, op, grid, stim, Cₘ, xs)

Diffusion half of the split: ∂ₜφₘ = ∇⋅(κ/(χCₘ))∇φₘ + Iₛₜᵢₘ/Cₘ on the contiguous voltage block of the state-blocked solution vector.

prepared is the MatrixFreeOperators PreparedOperator the right-hand side actually applies; op and grid are kept alongside it because a prepared operator cannot be rebuilt from itself — Adapt needs them to re-prepare on the target device, and tests need them to inspect what was assembled.

Not constructed directly; semidiscretize builds it.

FiniteDifferenceDiscretization(; order = 2)

Finite-difference spatial discretization descriptor.

Deliberately minimal: spacing, boundary conditions, and the device all live on the CartesianGrid, and duplicating them here is exactly the wart this replaces — ArmyHeart's discretization object validated an order and a boundary it then never read. The only thing this type carries is what the grid genuinely does not know, and the only reason it exists at all is to keep the three-argument semidiscretize(split, disc, grid) signature open for a fourth-order stencil later.

order must be 2; MatrixFreeOperators' second-derivative stencils are second order.

See also: semidiscretize.

MonodomainModel(; κ, ion, χ = 1, Cₘ = 1, stim = NoStimulationProtocol(),
                  φ_symbol = :φₘ, state_symbol = :states)

The monodomain equation

χCₘ∂ₜφₘ = ∇⋅κ∇φₘ + χ(Iᵢₒₙ(φₘ, s, t) + Iₛₜᵢₘ(x, t))
                 ∂ₜs = f(φₘ, s, t)

for transmembrane potential φₘ and cell-model states s.

Keyword Arguments

• κ: conductivity. A Number is isotropic; an NTuple{N} or SVector{N} is axis-aligned diagonal anisotropy, one entry per grid dimension (fibres along an axis). Spatially varying conductivity fields are not supported in v0.
• ion: any CytoZoo.AbstractCellModel. Supplies Iᵢₒₙ and the state kinetics; it is the reaction half of the split.
• χ: surface-to-volume ratio.
• Cₘ: membrane capacitance per unit area.
• stim: an AbstractStimulationProtocol. Applied in the diffusion half, as a source term, uniformly in every dimension.
• φ_symbol, state_symbol: names the solution helpers (getvariable, setvariable!) accept for the voltage block and for the whole non-voltage state block. state_symbol is deliberately plural: a cell model may itself name a state :s (CytoZoo's FHNModel does), and a singular default would shadow it.

The constructor is keyword-only on purpose. Thunderbolt's tutorial calls MonodomainModel(Cₘ, χ, …) while its fields are ordered (χ, Cₘ, …); the two are dimensionally interchangeable in the equation above, so a swapped pair is silent. Naming them removes the trap.

Units

Nothing here fixes a unit system — κ, χ, Cₘ, the stimulus, and the cell model must simply agree. What the code does fix is that κ enters as the diffusivity κ/(χCₘ) and the stimulus as Iₛₜᵢₘ/Cₘ.

Examples

using CytoZoo: FHNModel

model = MonodomainModel(; κ = 1.0e-3, ion = FHNModel())

# fibres along x, three-dimensional anisotropy
model = MonodomainModel(;
    κ = (0.133, 0.0176, 0.0176),
    χ = 140.0,
    Cₘ = 0.01,
    ion = FHNModel(),
)

See also: ReactionDiffusionSplit, semidiscretize.

NoStimulationProtocol()

The absence of an applied stimulus. Dispatching on this type removes the stimulus evaluation from the diffusion right-hand side entirely rather than multiplying by zero.

PointwiseODEFunction(ion, xs, overrides, nnodes, nstates, φ_symbol, state_symbol)

Reaction half of the split: the cell model ion evaluated independently at every node.

The solution vector is state-blocked (structure of arrays):

u = [φₘ(1:N); s₁(1:N); s₂(1:N); … ; s_{M-1}(1:N)]

so node i's full state is the strided slice u[i:N:N*M] and each state is contiguous. That layout is what makes the split cheap: the diffusion half acts on φₘ alone, which is the contiguous leading block 1:N, so its dof range needs no synchronizer object and no gather/scatter. It is also the coalesced layout on a GPU — neighbouring threads read neighbouring nodes of the same state.

Fields beyond the kinetics — nnodes, nstates, and the two symbols — are metadata for the solution helpers in solution.jl; the right-hand side itself uses only the first three plus the two counts.

Not constructed directly; semidiscretize builds it.

ReactionDiffusionSplit(model)
ReactionDiffusionSplit(model, overrides)

Annotate model for a reaction-diffusion operator split: the ionic kinetics are advanced pointwise, the diffusion operator globally. This is a declaration only — no work happens until semidiscretize.

overrides is the NamedTuple of spatial parameter overrides handed to the cell model through a CytoZoo.SpatialContext at every node. This is how spatial heterogeneity enters: cell type, pH, hypoxia, or any other parameter the cell model resolves spatially. Values may be numbers, callables (x, t) -> v, or isbits CytoZoo.SpatialFunction functors (Constant, SpatialStep, SpatialGradient) — only the last kind is GPU-safe.

With no overrides the reaction half calls the cell model with p = nothing, which is the fast non-spatial dispatch every CytoZoo model compiles down to.

Examples

using CytoZoo: SpatialStep

split = ReactionDiffusionSplit(model)

# endocardial on the left half of x, epicardial on the right
split = ReactionDiffusionSplit(model, (celltype = SpatialStep(1, 0.5, 0.0, 1.0),))

TransmembraneStimulationProtocol(f; nonzero_intervals = nothing)

Applied stimulus current f(x, t) -> Iₐₚₚ, in the positive-is-depolarizing convention documented on AbstractStimulationProtocol.

nonzero_intervals is an optional collection of (t₀, t₁) windows outside which f is known to vanish. When supplied, the diffusion right-hand side skips the whole stimulus evaluation for times outside every window — worth having when a single 2 ms pulse rides on a 500 ms simulation. It is stored as a Tuple, not a Vector, so the protocol stays isbits and can be captured by value in a GPU broadcast.

The protocol is isbits (and therefore GPU-safe) exactly when f is: a plain function or a small immutable functor is fine, a closure over a boxed variable is not.

Examples

# 2 ms depolarizing pulse over the left 1 mm of a cable
stim = TransmembraneStimulationProtocol(
    (x, t) -> (t <= 2.0 && x[1] <= 1.0) ? 50.0 : 0.0;
    nonzero_intervals = ((0.0, 2.0),),
)

apply_stimulus!(du, stim, xs, Cₘ, t) -> du

Add Iₛₜᵢₘ(x, t)/Cₘ to the voltage derivative at every node. A NoStimulationProtocol dispatches to a no-op — the evaluation is removed, not multiplied by zero — and a protocol with declared nonzero_intervals skips the broadcast outside them.

create_initial_condition(f::GenericSplitFunction, [T]) -> AbstractVector{T}

Allocate the solution vector and fill every node with the cell model's default_initial_state, on the same device the semidiscretization lives on.

T defaults to the grid's coordinate element type (Float64 for a Float64 grid), so the initial condition matches the operator rather than silently mixing precisions.

This is the spatially uniform starting point. Impose structure on top of it with setvariable! — an initial voltage profile, a heterogeneous gating state, and so on.

Examples

u₀ = create_initial_condition(f)
setvariable!(u₀, f, :φₘ) do x
    x[1] < 1.0 ? 1.0 : 0.0
end

diffusion_operator(model::MonodomainModel, grid::CartesianGrid) -> AbstractOperator

The matrix-free operator ∇⋅(κ/(χCₘ))∇, folded once at setup so the right-hand side carries no coefficient arithmetic.

A scalar conductivity gives κ_eff∇². An axis-aligned diagonal κ gives

κ_min∇² + Σ_d (κ_eff,d − κ_min)∂²/∂x_d²

which is exact for a diagonal tensor and cheaper than N separate second derivatives: the Laplacian carries the isotropic part in one halo exchange and only the genuinely anisotropic axes pay for a correction term. Axes equal to κ_min contribute nothing and are dropped, so an isotropic tuple collapses back to the single-Laplacian form.

getvariable(u, f::GenericSplitFunction, name::Symbol) -> SubArray

View of one named block of the solution vector.

name is a cell-model state name, the model's φ_symbol (:φₘ by default, aliasing whichever state the cell model reports as the transmembrane potential), or its state_symbol (:states by default) for the whole non-voltage block.

Being a view, writing through it writes the solution.

Examples

φ = getvariable(u, f, :φₘ)      # length N
s = getvariable(u, f, :states)  # length N*(M-1), all non-voltage states
w = getvariable(u, f, :CaMKt)   # length N, one named state

is_active(protocol, t) -> Bool

Whether protocol can be nonzero at time t. Conservative: a protocol without declared nonzero_intervals is always reported active. Used by the diffusion right-hand side to skip stimulus evaluation, never to decide the value.

node_coordinates(grid::CartesianGrid{N}) -> AbstractVector{SVector{N}}
node_coordinates(f::GenericSplitFunction) -> AbstractVector{SVector{N}}

Cell-center coordinates of every interior cell, flattened in the same order MatrixFreeOperators' flatten uses, on the grid's device. Given a semidiscretization it returns the very vector the right-hand sides use, so plotting φₘ against position needs no second construction.

num_nodes(f::GenericSplitFunction) -> Int

Number of spatial nodes N in a semidiscretized problem.

semidiscretize(split::ReactionDiffusionSplit, disc::FiniteDifferenceDiscretization,
               grid::CartesianGrid) -> GenericSplitFunction

Turn a continuous ReactionDiffusionSplit into a right-hand side on grid, ready for OperatorSplittingProblem.

The result is a two-operator GenericSplitFunction:

with N nodes and M states. No synchronizer objects are needed: under the state-blocked layout the voltage block the diffusion half acts on is the leading contiguous range, so the two dof ranges overlap exactly where the physics couples them.

Only zero-flux (homogeneous Neumann) boundaries are supported in v0; anything else raises an ArgumentError naming the offending face.

Examples

using CytoZoo: FHNModel

grid = CartesianGrid(((0.0, 20.0),), (200,); bc = ((Neumann(), Neumann()),))
model = MonodomainModel(; κ = 1.0e-3, ion = FHNModel())
f = semidiscretize(ReactionDiffusionSplit(model), FiniteDifferenceDiscretization(), grid)
u₀ = create_initial_condition(f)
prob = OperatorSplittingProblem(f, u₀, (0.0, 100.0))

See also: create_initial_condition, solution_size.

setvariable!(fun, u, f::GenericSplitFunction, name::Symbol) -> u
setvariable!(u, f::GenericSplitFunction, name::Symbol, value::Number) -> u

Fill one named state block by evaluating fun(x) at every node's cell-center coordinate, or with a constant value.

Reads naturally with do-block syntax, since the function is the first argument.

name must name a single state (a cell-model state name or the model's φ_symbol); the whole-state-block symbol is rejected, because one function of position cannot say what to put in each of several states.

Examples

setvariable!(u, f, :φₘ) do x
    cos(π * x[1] / L)
end
setvariable!(u, f, :φₘ, 0.0)

solution_size(f::GenericSplitFunction) -> Int

Length of the solution vector: N * M for N nodes and M cell-model states per node.

variable_range(f::GenericSplitFunction, name::Symbol) -> UnitRange{Int}

Index range of a named block in the solution vector — the raw form of getvariable, useful for building views of an integrator's u without holding one.