Scientific computing: SciML and time integration#

semidiscretize(a::BilinearForm, l::LinearForm; kwargs...) -> Semidiscretization

Semidiscretise M u_h' = F(t) - A u_h from the spatial form a and the source l, both posed on the same test space.

Write a as the steady problem is written: assemble(a) is A, so the steady state of the returned system solves A u_h = F.

Keywords

• mass: BilinearForm defining M (default innerₕ(u, v), the discrete L² inner product, which is diagonal).

• dirichlet: constrained labels and, where they carry values, the values – every form assemble accepts, including time-dependent constraints from dirichlet_constraints(Ωₕ, I, :label => (x, t) -> ...), or parameter-dependent ones from (x, t, p) -> ... for a value reached by the residual's own p (default nothing).

• dirichlet_components: leaf components of a composite space the labels bind to (default nothing, all leaves).

• state: VectorElement the current u is copied into before each assembly, for forms whose coefficients read it (default nothing). A u of an element type the state cannot hold (a ForwardDiff.Dual, from automatic differentiation of the residual) goes into a buffer of its own type instead, read by the forms rebound to it, so the user's state never receives it and -(∂A/∂u) u enters the Jacobian. Only the state written in the forms' expressions is rebound. One read inside a Julia function (x -> uₕ(x)), through a component view (uₕ(1)) or by a build-based operator is not, and such a u is refused with an ArgumentError where nothing is rebound; so is any such u beside an update_coefficients! hook, which may read the user's state. A plain number type (a real or complex float or integer) is always copied into the state, converted to its element type.

• update_coefficients!: called with the current t before each assembly, for coefficients that vary in time – t -> Rₕ!(fₕ, x -> f(x, t)) for a time-dependent source, or t -> (α[] = t) for a scalar Ref (default nothing). A two-argument (t, p) -> ... also reaches the residual's own p, the same way a three-argument Dirichlet condition does.

• hook_reads_state: whether update_coefficients! may read the state (default true). With state, a ForwardDiff.Dual u reaches the forms through a buffer of its own type and never the user's state, so a hook reading the state would lose its dependence on u from the derivative, and the pair is refused with an ArgumentError. Give false to assert that the hook reads no state (a source depending on t alone) and the Dual Jacobian is exact; the assertion is not checked, and a hook that does read the state then gets a Jacobian missing that part.

• reassemble: refill A at every step, for an operator whose coefficients change with t or u (default false).

• coefficients: the coefficients of a that depend on u, declared as jacobian_form takes them, (cᵢ => (U -> Lᵢ(U), fᵢ), ...) with cᵢ the grid function a reads (default nothing). Before each assembly of the residual every cᵢ is set to fᵢ.(Lᵢ(u)) from u itself, after update_coefficients! and the state, and A is reassembled (reassemble is then true whatever it was given). jacobian! then writes the exact Jacobian -(A + ∂A) at u and jacobian_prototype carries its pattern, wider than A's where a coefficient reads neighbouring points (using ForwardDiff differentiates each fᵢ). A source l or a mass reading a declared coefficient is refused with an ArgumentError, since jacobian! would miss its derivative. A u of a non-plain element type (a ForwardDiff.Dual) is refused by the residual, semidiscretize_rhs and jacobian! alike, with or without state: automatic differentiation through the coefficients is not supported, and the exact Jacobian is jacobian!'s. Scalar grid spaces only: a on a composite space is refused with an ArgumentError, and the build method and ghost spaces do not take it.

Time and parameter dependence of the Dirichlet values is detected by arity, exactly as dirichlet_constraints validates it: conditions accepting (x, t, p) are evaluated at each step against the residual's own p, conditions accepting (x, t) are evaluated at each step without it, and conditions accepting (x) are not evaluated per step at all. update_coefficients! is detected between its one- and two-argument forms, but prefers the one-argument form: (t, p = default) -> ... is called as f(t) and ignores p. A Base.Fix or composition is classified by the function it wraps.

Examples

Ωₕ = mesh(domain(interval(0.0, 1.0)), 101)
Wₕ = gridspace(Ωₕ)
I = interval(0.0, 1.0)

fₕ = Rₕ(Wₕ, x -> 1.0)
a = form(Wₕ, Wₕ, (u, v) -> inner₊(∇ₕ(u), ∇ₕ(v)))
l = form(Wₕ, v -> innerₕ(fₕ, v))
bcs = dirichlet_constraints(Ωₕ, I, :boundary => (x, t) -> 0.0)

sd = semidiscretize(a, l; dirichlet = bcs)

Forms posed on a constrained ghost space, whose relation data may depend on time ((x, t) -> ..., see constrained), are taken too, by the method semidiscretize(a::GhostBilinearForm, l::GhostLinearForm; mass, state, update_coefficients!, hook_reads_state, reassemble): the relations eliminate the ghosts and the Dirichlet nodes, so the mass has no zero rows, and the source at t carries the lift of the data at t. That method refuses the dirichlet keyword (give :F => u => g(x, t) in constrained), a space keeping ghosts as unknowns, a mass coupling a degree of freedom to a point whose data depend on time, semidiscretize(build, l) with a ghost load, and semidiscretize_rhs, each with an ArgumentError naming the alternative.

See also ode_function, ode_problem, jacobian!.

semidiscretize(build, l::LinearForm; kwargs...) -> Semidiscretization

Semidiscretise M u_h' = F(t) - A(t) u_h from a spatial operator that genuinely depends on t – built fresh, once per element type t is ever reached at, instead of the single BilinearForm the other method assembles into one fixed matrix.

build(t) is called once for each element type t is ever seen at (the space's element type on a normal step, a ForwardDiff.Dual while a Rosenbrock stepper's tgrad differentiates the residual through t), and must return (a, refill!): the BilinearForm to assemble, built around whatever live coefficient buffer(s) it needs, and a one-argument refill!(t) updating those buffers from the current t – called on every step, cache hit or miss, so a later step at an already-seen type still sees the new t rather than the one build first saw. Same build/refill!/cache discipline as type_cached_assemble!, keyed on t instead of a VectorElement iterate.

The probe at construction calls build at zero(eltype(Wₕ)), Wₕ the space of l, so the time span of the solve should use that same element type: a Float32 space solved over a Float64 time span reaches a second element type and builds once more for it. A step at any element type other than the probe's goes through an untyped cache, so the residual and jacobian! allocate on every call there, not only for a source buffer.

build should be a named function, not a closure literal written inline – the same reason type_cached_assemble!'s own docstring gives: a do ... end block re-literalised on every call allocates a new closure each time.

This is what lets a Rosenbrock method (Rodas5P, Rosenbrock23) capture -Ȧ(t) u_h by differentiating through t directly: the other method's fixed BilinearForm closes over a coefficient buffer of the space's element type, which cannot hold a Dual and throws InexactError under exactly that AD sweep.

Keywords

Same as the BilinearForm method, except reassemble defaults to true – a build-based operator exists specifically to be rebuilt at every step – and coefficients is not taken.

Examples

function build_diffusion_operator(t)
    αₕ = element(Wₕ, typeof(t))
    a = form(Wₕ, Wₕ, (u, v) -> inner₊(αₕ * ∇ₕ(u), ∇ₕ(v)))
    refill!(t) = (Rₕ!(αₕ, x -> α(x, t)); nothing)
    return a, refill!
end

sd = semidiscretize(build_diffusion_operator, l; dirichlet = bcs)

See also ode_function, ode_problem.

semidiscretize(a::GhostBilinearForm, l::GhostLinearForm; mass, dirichlet,
    dirichlet_components, state, update_coefficients!, hook_reads_state, reassemble) ->
    Semidiscretization

Semidiscretise M u' = F(t) − A u for forms posed on one constrained ghost space Cₕ, whose relation data may depend on time. With (E, d(t)) the closure of Cₕ, A = Eᵀ A_pad E, M = Eᵀ M_pad E and F(t) = Eᵀ (F_pad(t) − A_pad d(t)), so at a fixed t the residual is F − A u of assemble(a, l) on the space rebuilt with the data frozen at t.

Keywords

• mass: a form on Cₕ defining M (default innerₕ(u, v), diagonal and positive). A mass that couples a degree of freedom to a point where d depends on time is refused: its term Eᵀ M_pad d′(t) needs the time derivative of the data.

• state: a grid function on Bramble.padded_space(Cₕ) that receives lift(Cₕ, u, t) before each assembly, for forms whose coefficients read it (default nothing). A u it cannot hold (a ForwardDiff.Dual) is lifted into a padded buffer of its own type, read by the forms rebound to it, as on a plain space.

• update_coefficients!, hook_reads_state, reassemble: as for plain forms.

The dirichlet and dirichlet_components keywords are refused: give a Dirichlet condition as :F => u => g(x, t) in constrained. So are a space that keeps ghosts as unknowns (their mass rows vanish), semidiscretize(build, l) with a ghost load (use update_coefficients! and reassemble = true), and semidiscretize_rhs.

Semidiscretization{...}

Method-of-lines semidiscretisation of M u_h' = F(t) - A u_h, callable with the (du, u, p, t) signature a time stepper expects.

Built by semidiscretize; read back with mass_matrix and operator_matrix.

Fields

• operator: spatial BilinearForm (or a ghost-space form), assembled into A; for a time-dependent operator, a TypeCachedOperator wrapping the form builder.

• source: source LinearForm, assembled into F(t) at each step.

• space: the test space both forms share.

• operator_matrix: A, with eₖ rows on the constrained degrees of freedom (none on a ghost space, whose relations eliminate them).

• mass_matrix: M, with zero rows on the constrained degrees of freedom (none on a ghost space, whose relations eliminate them).

• source_vector: the reusable F buffer, refilled in place.

• constraints: how the source's boundary values are reached – one of NoConstraints, LabelsOnly, StaticConstraints, TimeDependentConstraints, TimeParamDependentConstraints or, on a ghost space, GhostConstraints.

• labels: constrained boundary labels.

• components: leaf components the labels bind to, or nothing for all.

• state: VectorElement receiving u before assembly, or nothing.

• update_coefficients: callable invoked with t (or, wrapped in ParametricUpdate, with t and the residual's own p) before assembly, or nothing. Wrapped in Bramble._StateFreeUpdate when the user declared hook_reads_state = false.

• reassemble: Val(true) to refill A at every step.

• twins: by element type of u, the semidiscretisations a u the state cannot hold (a ForwardDiff.Dual) is evaluated by: a state buffer of that type, the forms rebound to it and their own operator matrix, built at the first call there.

• jacobian_form: the form a + ∂a of jacobian_form built from the coefficients keyword, whose refills set each declared coefficient from u before every assembly, or nothing.

semidiscretize_rhs(sd::Semidiscretization) -> SemidiscretizeRHS

Build the matrix-free right-hand side of uₕ' = M⁻¹(F(t) - A uₕ) from sd, folding M's diagonal into a precomputed scaling instead of leaving M for a solver to factorise at every step – valid only because M is diagonal and, here, invertible everywhere.

Requires sd to have been built with dirichlet = nothing: any Dirichlet label makes the corresponding row of M zero by construction (see semidiscretize), which is an algebraic constraint on uₕ, not an equation for uₕ' – there is no uₕ' value that divides it away. Also requires mass_matrix(sd) to actually be diagonal: the default mass (the discrete L² inner product) always assembles diagonally, but a caller-supplied mass keyword to semidiscretize need not.

When the space's execution policy is CpuPolyester (using Polyester required), sd was built with reassemble = false and A is a SparseMatrixCSC, the product A u runs one row per Polyester.@batch iteration, each row writing its own entry of du. That needs A's pattern in compressed-row order, which is built here once. It takes two Int arrays of length nnz(A) (the column of each entry, and its position in nonzeros(A)) and one of length size(A, 1) + 1, and costs about three products to build. The values are not copied. They are read from A itself at every call, so editing nonzeros(operator_matrix(sd)) afterwards is seen. A change to A's stored pattern is not seen, and if it changes nnz(A) the call falls back to the serial product. The row sums run in a different order than the serial column product, so the result agrees with CpuSerial's to rounding, not bit for bit. With reassemble = true, A changes at every step and a fresh copy would cost more than the threads save. That case, like every other policy, keeps the serial product.

Examples

using Bramble: semidiscretize_rhs
sd = semidiscretize(a, l)  # no `dirichlet` keyword: NoConstraints
rhs = semidiscretize_rhs(sd)
prob = ODEProblem(rhs, parent(u₀), (0.0, 1.0))  # no mass_matrix to factorise
sol = solve(prob, Tsit5())

See also semidiscretize, ode_problem.

SemidiscretizeRHS{S,V,C}

The right-hand side of uₕ' = M⁻¹(F(t) - A uₕ), M folded into a precomputed diagonal scaling rather than solved for – for a Semidiscretization with no Dirichlet constraints, where that fold is valid. Built by semidiscretize_rhs; callable with the (du, u, p, t) signature a plain (non-mass-matrix) SciMLBase.ODEProblem expects.

Named SemidiscretizeRHS/semidiscretize_rhs rather than the free function semidiscretize_rhs!(du, u, p, t) gpena/Bramble.jl#163 proposes: that signature is exactly SciMLBase.ODEFunction's own, which leaves no argument to pass a Semidiscretization or the precomputed scaling through – both have to be closed over somehow, and a callable struct is what every other stateful callable in this package does instead of a closure (Semidiscretization itself, TypeCachedOperator).

The csr field holds the compressed-row pattern of A that a CpuPolyester space with reassemble = false multiplies by (see semidiscretize_rhs), and is nothing otherwise. The twins field caches, per element type of a u the state cannot hold (a ForwardDiff.Dual), the right-hand side of the residual's twin at that type.

mass_matrix(sd::Semidiscretization) -> AbstractMatrix

Return the constant mass matrix M, whose constrained rows are zero.

operator_matrix(sd::Semidiscretization) -> AbstractMatrix

Return the assembled spatial operator A, whose constrained rows are eₖ.

jacobian!(J, sd::Semidiscretization, u, p, t) -> J

Fill J with ∂/∂u (F(t) - A(u) u) and return it.

Without coefficients this is -A, exact whenever the operator's coefficients do not read u – that is, whenever sd was built without state. With a state the operator also varies with u, the term -(∂A/∂u) u is missing, and this is a Picard linearisation rather than a Jacobian: declare the coefficients that read u with the coefficients keyword of semidiscretize, or pass jacobian = nothing to ode_function and let the solver build it by sparse automatic differentiation instead, seeded from jacobian_pattern. That one is exact: a ForwardDiff.Dual u goes into a state buffer of its own type, read by the forms rebound to it (see the state keyword of semidiscretize).

With coefficients it is the exact Jacobian -(A + ∂A) at u, assembled from the form a + ∂a of jacobian_form after every declared coefficient is set from u, so J must carry jacobian_prototype(sd)'s pattern. using ForwardDiff is needed.

See also jacobian_prototype.

jacobian_prototype(sd::Semidiscretization) -> AbstractMatrix

Return a matrix carrying the sparsity of ∂/∂u (F(t) - A u), for a solver to use as its Jacobian cache.

This is the pattern of the assembled operator, Dirichlet rows included, which is what the system's Jacobian has when the operator does not read u. With coefficients it is the pattern of a + ∂a (jacobian_form), wider than A's where a coefficient reads neighbouring points. jacobian_pattern answers the different question of what a Newton residual's Jacobian looks like when the form's coefficients depend on the solution.

See also jacobian!.

ode_function(sd::Semidiscretization; kwargs...) -> ODEFunction
ode_function(a::BilinearForm, l::LinearForm; kwargs...) -> ODEFunction

Wrap a semidiscretisation as an ODEFunction carrying its mass matrix, Jacobian and sparsity, ready for OrdinaryDiffEq. On a csr_backend() the solver receives SparseMatrixCSC copies of the mass matrix and of the Jacobian prototype, since no solver factorizes a SparseMatrixCSR.

The two-form method builds the Semidiscretization first, forwarding every keyword to semidiscretize.

Keywords

• jacobian: jacobian! (the default) to hand the solver the exact -A, or nothing to let it build one by automatic differentiation from jac_prototype.

• jac_prototype: sparsity for the solver's Jacobian cache (default: jacobian_prototype(sd)). A SparseMatrixCSR, given or defaulted, reaches the solver as a SparseMatrixCSC copy.

• tgrad: analytical ∂f/∂t, for a Rosenbrock method to use instead of differentiating through t (default: nothing, AD). Either SciMLBase's own (dT, u, p, t) -> ... signature or the Bramble-aware (dT, sd, u, p, t) -> ... one, told apart by arity.

The resulting system is a differential-algebraic one whenever any Dirichlet label is constrained, since those rows of the mass matrix are zero. Solve it with a method that admits a singular mass matrix – FBDF, QNDF, Rodas5P, RadauIIA5 – not an explicit one.

Requires SciMLBase.jl; call using SciMLBase (or any package that loads it, such as OrdinaryDiffEq) before calling this function.

See also ode_problem, linear_problem, nonlinear_problem.

ode_problem(sd::Semidiscretization, u₀, I; kwargs...) -> ODEProblem
ode_problem(a::BilinearForm, l::LinearForm, u₀, I; kwargs...) -> ODEProblem
ode_problem(a::Bramble.GhostBilinearForm, l::Bramble.GhostLinearForm, u₀, I; kwargs...) -> ODEProblem

Build the ODEProblem stepping sd over the time domain I, from the initial condition u₀ – a VectorElement or a plain vector. I may equally be a (t₀, t₁) tuple.

u₀ is copied, never mutated, and the copy is made consistent with the Dirichlet rows at t₀ (see dirichlet_bc!), against p when one is given.

Keywords are those of ode_function, plus:

• p (default SciMLBase.NullParameters()) for a residual whose update_coefficients! or Dirichlet conditions were given a parameter-dependent, (t, p)/(x, t, p) form (see semidiscretize's own keywords).

• specialize (default nothing, ODEProblem's own choice untouched).

The two-form method, for plain forms and for forms on one constrained ghost space, forwards its other keywords to semidiscretize.

Requires SciMLBase.jl.

Examples

sd = semidiscretize(a, l; dirichlet = bcs)
prob = ode_problem(sd, Rₕ(Wₕ, x -> sinpi(x[1])), interval(0.0, 1.0))
sol = solve(prob, FBDF())

bcs = dirichlet_constraints(Ωₕ, I, :boundary => (x, t, p) -> p[1] * t)
sd = semidiscretize(a, l; dirichlet = bcs)
prob = ode_problem(sd, u₀, I; p = [0.7])

See also ode_function, semidiscretize.

ode_problem(rhs::SemidiscretizeRHS, u₀, I; kwargs...) -> ODEProblem

Build the plain (non-mass-matrix) ODEProblem stepping rhs over the time domain I, from the initial condition u₀ – a VectorElement or a plain vector, copied, never mutated.

Unlike ode_problem(sd::Semidiscretization, ...), there is no Dirichlet consistency step: semidiscretize_rhs only ever builds rhs from a Semidiscretization with dirichlet = nothing, so there are no boundary rows to make consistent.

Requires SciMLBase.jl.

Examples

using Bramble: semidiscretize_rhs
sd = semidiscretize(a, l)
rhs = semidiscretize_rhs(sd)
prob = ode_problem(rhs, Rₕ(Wₕ, x -> sinpi(x[1])), interval(0.0, 1.0))
sol = solve(prob, Tsit5())

See also semidiscretize_rhs, ode_problem(::Semidiscretization, ...).

linear_problem(a::BilinearForm, l::LinearForm; kwargs...) -> LinearProblem
linear_problem(a::GhostBilinearForm, l::GhostLinearForm; kwargs...) -> LinearProblem

Assemble a and l into the LinearProblem that LinearSolve.solve takes, so the steady system reaches the factorisations, iterative solvers and preconditioners in that stack without being assembled by hand first.

On ghost spaces (a constrained space, or the pair polygon_spaces returns), the problem is the reduced system A, F = assemble(a, l), whose unknowns are the ndofs of the trial space, exactly as pde_solve(A, F) would take it.

Keywords

• dirichlet, dirichlet_components: as assemble takes them.

• symmetrize: restore symmetry after imposing the conditions (default: false; see symmetrize!).

Requires SciMLBase.jl, which defines LinearProblem; LinearSolve itself is needed only to solve the result.

Examples

using LinearSolve, IncompleteLU

prob = linear_problem(a, l; dirichlet = bcs)
sol = solve(prob, KrylovJL_GMRES())

See also assemble, ode_problem, nonlinear_problem.

nonlinear_problem(residual, u0; jacobian = nothing, jac_prototype = nothing, kwargs...) -> NonlinearProblem

Wrap residual – the residual of a steady nonlinear discretisation F(u) = 0, as F(u, p) or in-place F!(res, u, p) – into the NonlinearProblem that NonlinearSolve.solve takes.

u0, the initial guess, is a VectorElement or a plain vector; copied, never mutated. Every Bramble worked example builds residual from an assembled BilinearForm, A(u) * u - F evaluated at the current iterate, so Dirichlet rows already carry their own identity/value pair and a Newton step corrects them from any initial guess – no separate consistency step, unlike ode_problem's differential-algebraic system.

Keywords

• jacobian: jac(J, u, p) if residual is in-place, or J = jac(u, p) if it is not – the two conventions cannot be mixed, the same requirement NonlinearFunction itself enforces. Or a jacobian_form jf, whose exact Jacobian jacobian!(J, jf, u) is then used in that same convention, with jac_prototype defaulting to jacobian_prototype(jf). Default nothing: the solver differentiates residual itself.

• jac_prototype: sparsity for the solver's Jacobian cache, e.g. built from jacobian_pattern. Default nothing (dense), or jacobian_prototype(jf) when jacobian is a jacobian_form.

• Every other keyword forwards to NonlinearProblem.

Requires SciMLBase.jl; call using SciMLBase (or any package that loads it, such as NonlinearSolve) before calling this function.

Examples

function residual!(res, u_vec::AbstractVector{T}, p) where {T}
    uₕ = element(Wₕ, T)
    uₕ .= u_vec
    A = assemble(diffusion_form(uₕ); dirichlet = :boundary)
    mul!(res, A, u_vec)
    return res .-= F
end

prob = nonlinear_problem(residual!, zeros(ndofs(Wₕ)))
sol = solve(prob, NewtonRaphson())

See also ode_problem, linear_problem, jacobian_pattern.

type_cached_assemble!(build, cache::AbstractDict, uₕ::VectorElement;
    dirichlet = nothing, dirichlet_components = nothing) -> AbstractMatrix

Assembles a coefficient-dependent BilinearForm into a matrix whose sparsity pattern is built once per distinct element type uₕ is ever passed at, rather than on every call, as diffusion_matrix in the nonlinear worked example does without it: assemble/allocate_system_matrix rebuild the whole matrix, pattern included, every time otherwise.

build(uₕ) is called once for each element type uₕ is ever seen at, and must return (a, refill!): the BilinearForm to assemble, built around whatever live coefficient buffer(s) it needs (see the forms tutorial), and a one-argument function refill!(uₕ) updating those buffers from the current uₕ – called on every invocation, cache hit or miss, so a later call at an already-seen type still sees the new guess rather than the one build first saw.

build itself should be a named function defined once, not a closure literal written inside whatever function calls type_cached_assemble!, the same reason the Picard loop in poisson_nonlinear.jl builds its own form once, outside the loop, rather than on every iteration: a do ... end block re-literalized on every call allocates a new closure each time, which is exactly the cost this function exists to avoid paying more than once.

using Bramble: Mₓ!, type_cached_assemble!
function build_diffusion(uₕ)
    Mu = element(Wₕ, eltype(uₕ))     # scratch for Mₓ!'s own output
    αvals = element(Wₕ, eltype(uₕ))
    a = form(Wₕ, Wₕ, (U, V) -> inner₊(αvals * ∇ₕ(U), ∇ₕ(V)))
    refill!(uₕ) = begin
        Mₓ!(Mu, uₕ)          # in place: `Mₓ(uₕ)` alone would allocate a fresh result
        αvals .= α.(Mu)
    end
    return a, refill!
end

cache = Dict()
diffusion_matrix(uₕ) = type_cached_assemble!(
    build_diffusion, cache, uₕ; dirichlet = :boundary)

refill! reaches for Mₓ! rather than the non-mutating Mₓ/Mₕ deliberately: the latter allocates a fresh result every call (the same similar-based cost every allocating stencil operator has), which would silently reintroduce an O(n) allocation this function's whole point is to stop paying repeatedly. Mₓ! alone covers the 1D case above; a D-dimensional coefficient needs one scratch buffer and one Mₓ!/Mᵧ!/M₂! call per direction, the same way poisson_nonlinear.jl's own nonlinear_series builds a D-dimensional coefficient tuple.

cache is shared across an entire Newton (or Picard) loop, one Dict per residual: the first call at a given type pays build's own cost plus allocate_system_matrix's; every later call at that same type pays only assemble!'s refill plus a small, fixed dictionary/dynamic-dispatch overhead fetching the cached entry back out (independent of ndofs) – not the O(ndofs) pattern rebuild a cache miss (or no cache at all) pays every time.

Not thread-safe. cache is a bare Dict with no lock, one per residual. Assembling from more than one task needs a lock or one cache per task.