Scientific computing: SciML, AD and solvers

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).
• 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.
• reassemble: refill A at every step, for an operator whose coefficients change with t or u (default: false).

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 the same way, between its one- and two-argument forms.

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)

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 (Float64 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.

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 Float64-typed coefficient buffer, 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.

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.

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, assembled into A.
• 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.
• mass_matrix: M, with zero rows on the constrained degrees of freedom.
• source_vector: the reusable F buffer, refilled in place.
• constraints: how the source's boundary values are reached – one of NoConstraints, LabelsOnly, StaticConstraints, TimeDependentConstraints or TimeParamDependentConstraints.
• 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.
• reassemble: Val(true) to refill A at every step.

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.

Examples

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}

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).

mass_matrix(sd::Semidiscretization) -> AbstractMatrix

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

mass_matrix(sd::SecondOrderSemidiscretization) -> 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) = -A and return it.

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: pass jacobian = nothing to ode_function and let the solver build it by sparse automatic differentiation instead, seeded from jacobian_pattern.

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 – 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.

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)).
• 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::CartesianProduct{1}; kwargs...) -> ODEProblem
ode_problem(a::BilinearForm, l::LinearForm, 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) – pass SciMLBase.FullSpecialize before handing the solved trajectory to Bramble.adjoint_sensitivities: without it, a p-vjp calls the residual with a differently-eltype-p than the forward solve used, which the default specialization cannot dispatch and fails with "No matching function wrapper was found!" rather than differentiating.

The two-form method also 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

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

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.

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. 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).
• 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.

adjoint_sensitivities(sol::ODESolution, alg; kwargs...) -> (du0, dp)

Adjoint sensitivities of a Semidiscretization's solved trajectory, via SciMLSensitivity.adjoint_sensitivities with two Bramble-specific corrections applied. Full documentation lives on BrambleSciMLSensitivityExt's own method, the only one that exists once SciMLSensitivity is loaded – this stub exists so that method has a function to extend, and so calling this without SciMLSensitivity loaded gives a clear error rather than UndefVarError.

Deliberately not exported, unlike pde_solve: SciMLSensitivity itself exports a function of this exact name, so using Bramble, SciMLSensitivity together would collide on the bare name regardless of what Bramble does. Call this one as Bramble.adjoint_sensitivities.

Bramble.adjoint_sensitivities(sol::ODESolution, alg; kwargs...) -> (du0, dp)

Adjoint sensitivities of a Semidiscretization's solved trajectory sol (from Bramble.ode_problem/solve) with respect to its initial condition (du0) and its p (dp), via SciMLSensitivity.adjoint_sensitivities – one backward solve for every parameter at once, the same O(1)-in-parameter-count trade Bramble.pde_solve's own adjoint rule makes for the steady case.

Keywords

Every keyword SciMLSensitivity.adjoint_sensitivities takes, plus these two defaults chosen for a Semidiscretization's index-1 DAE specifically (both overridable):

• sensealg: InterpolatingAdjoint(autojacvec = false) – autojacvec = false uses Bramble.jacobian!'s own exact -A for the u-vjp, so no AD tool ever needs to differentiate through the residual's mutating buffers for that half.
• initializealg: BrownFullBasicInit() – restores the adjoint's own algebraic consistency at t = T rather than merely checking it (the default CheckInit rejects a constrained problem's seed outright); see this file's own module-level comment for why.

Requirements on sol

sol.prob must have been built with ode_problem(sd, u₀, I; p = ..., specialize = SciMLBase.FullSpecialize) – both p and specialize matter: specialize avoids a function-wrapper error the moment a p-vjp is computed (see the note on ode_problem itself), and without a real p there is nothing for dp to be a gradient with respect to.

Examples

using Bramble, SciMLBase, SciMLSensitivity, OrdinaryDiffEqBDF

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], specialize = SciMLBase.FullSpecialize)
sol = solve(prob, FBDF(); saveat = ts)

dgdu!(out, u, p, t, i) = (@. out = 2 * (u - obs[i]); nothing)
du0, dp = Bramble.adjoint_sensitivities(sol, FBDF(); t = ts, dgdu_discrete = dgdu!)

See also Bramble.ode_problem, Bramble.pde_solve for the steady-state adjoint this mirrors.

semidiscretize_second_order(K::BilinearForm, l::LinearForm; kwargs...) -> SecondOrderSemidiscretization

Semidiscretise M üₕ + C u̇ₕ + K uₕ = F(t) from the spatial stiffness form K and the source l, both posed on the same test space.

Write K as the steady problem is written: assemble(K) is the same matrix semidiscretize calls A.

Keywords

• mass: BilinearForm defining M (default: innerₕ(u, v), the discrete L² inner product, which is diagonal).
• damping: BilinearForm defining C (default: nothing, an undamped system).
• 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) -> ...) (default: nothing).
• dirichlet_components: leaf components of a composite space the labels bind to (default: nothing, all leaves).

Dirichlet conditions constrain u; the consistent velocity at a constrained dof follows from differentiating that constraint in time and is never prescribed independently – see this file's header comment for why, and for the solver-compatibility warning that follows from it.

Examples

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

fₕ = Rₕ(Wₕ, x -> 0.0)
K = form(Wₕ, Wₕ, (u, v) -> inner₊(∇ₕ(u), ∇ₕ(v)))
l = form(Wₕ, v -> innerₕ(fₕ, v))

sd = semidiscretize_second_order(K, l)

See also second_order_ode_function, second_order_ode_problem.

SecondOrderSemidiscretization{...} <: AbstractSemidiscretization

Method-of-lines semidiscretisation of M üₕ + C u̇ₕ + K uₕ = F(t), callable with the (dv, v, u, p, t) signature a second-order time stepper expects.

Built by semidiscretize_second_order; read back with mass_matrix, damping_matrix and stiffness_matrix.

Fields

• stiffness: spatial BilinearForm, assembled into K.
• damping: spatial BilinearForm, assembled into C, or nothing for an undamped system.
• source: source LinearForm, assembled into F(t) at each step.
• space: the test space every form shares.
• stiffness_matrix: K, with eₖ rows on the constrained degrees of freedom.
• damping_matrix: C, with zero rows on the constrained degrees of freedom, or nothing.
• mass_matrix: M, with zero rows on the constrained degrees of freedom.
• source_vector: the reusable F buffer, refilled in place.
• constraints: how the source's boundary values are reached – one of NoConstraints, LabelsOnly, StaticConstraints or TimeDependentConstraints.
• labels: constrained boundary labels.
• components: leaf components the labels bind to, or nothing for all.

damping_matrix(sd::SecondOrderSemidiscretization) -> Union{AbstractMatrix, Nothing}

Return the assembled damping operator C, whose constrained rows are zero, or nothing for an undamped system.

stiffness_matrix(sd::SecondOrderSemidiscretization) -> AbstractMatrix

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

block_mass_matrix(sd::SecondOrderSemidiscretization) -> AbstractMatrix

Return the 2n × 2n block-diagonal blockdiag(mass_matrix(sd), I), the mass matrix second_order_ode_function hands DynamicalODEFunction: the I block leaves the kinematic equation u' = v an ordinary ODE, and the mass_matrix(sd) block carries whatever Dirichlet row-zeroing sd was built with.

second_order_ode_function(sd::SecondOrderSemidiscretization) -> DynamicalODEFunction
second_order_ode_function(K::BilinearForm, l::LinearForm; kwargs...) -> DynamicalODEFunction

Wrap a second-order semidiscretisation as a DynamicalODEFunction carrying its block mass matrix, ready for OrdinaryDiffEq.

The two-form method builds the SecondOrderSemidiscretization first, forwarding every keyword to semidiscretize_second_order.

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

See also second_order_ode_problem, semidiscretize_second_order.

second_order_ode_problem(sd::SecondOrderSemidiscretization, du₀, u₀, I::CartesianProduct{1}) -> SecondOrderODEProblem
second_order_ode_problem(K::BilinearForm, l::LinearForm, du₀, u₀, I; kwargs...) -> SecondOrderODEProblem

Build the SecondOrderODEProblem stepping sd over the time domain I, from the initial velocity du₀ and initial displacement u₀ – each 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!). du₀ is passed through unchanged: velocity at a constrained dof is not an independent algebraic unknown here (see this file's header comment).

Keywords are those of second_order_ode_function; the two-form method also forwards to semidiscretize_second_order.

Requires SciMLBase.jl.

Examples

sd = semidiscretize_second_order(K, l; dirichlet = bcs)
prob = second_order_ode_problem(sd, Rₕ(Wₕ, x -> 0.0), Rₕ(Wₕ, x -> sinpi(x[1])), interval(0.0, 1.0))
sol = solve(prob, Rodas5P(); reltol = 1e-11, abstol = 1e-13)

See also second_order_ode_function, semidiscretize_second_order.

pde_solve(A::SparseMatrixCSC, F::AbstractVector; solver = :default, sym = :auto, kwargs...) -> Vector
pde_solve(fact::MUMPSFactorization, F::AbstractVector) -> Vector

Solve A u = F (or fact u = F) and return u.

With no keyword arguments this is A \ F, unless AppleAccelerate.jl is loaded on macOS and the system is symmetric (see :default below), in which case it is accelerate_solve(A, F; sym). The name exists so that reverse-mode AD tools have one function to attach an adjoint rule to, which is what makes a whole θ -> assemble -> pde_solve -> J(u) chain differentiable.

Keywords

• solver: :default – A \ F on Linux, Windows, macOS without AppleAccelerate.jl loaded, and macOS on an unsymmetric A; on macOS with using AppleAccelerate in effect and A symmetric, dispatches to accelerate_solve instead (same as passing solver = :accelerate explicitly). Symmetry is issymmetric(A) under sym = :auto (the default), or trusted outright from an explicit sym = :spd/:definite/:symmetric (and conversely :unsymmetric skips straight to A \ F) – the caller has already asserted the property, so it is not checked again. The narrowing exists because Accelerate is only a win on the symmetric factorisations it reaches (SPD/Cholesky, LDLᵀ): measured on this host against forms assembled on real 2D grid spaces, :default was a 1.2-1.3x win on a symmetric Poisson-plus-mass system and a 2.3-3.6x loss on an unsymmetric convection-diffusion one (gpena/Bramble.jl#246). This only ever narrows which solver runs on macOS; Linux and Windows are never affected, and a matrix Accelerate cannot factor (e.g. non-square) throws from accelerate_solve exactly as solver = :accelerate would – :default never silently falls back to \ after picking Accelerate. :suitesparse (CHOLMOD/UMFPACK), :spqr (sparse QR, for a least-squares or rectangular A), :accelerate (Apple libSparse, needs AppleAccelerate.jl, honoured unconditionally regardless of symmetry), :mumps (needs MUMPS.jl), :sparspak (pure Julia, needs Sparspak.jl).
• sym: symmetry hint for :suitesparse, :accelerate and :mumps (and, on macOS, for :default's own choice of solver – see above). :auto (default) detects it; :spd/:definite/1, :symmetric/2 and :unsymmetric/0 state it.

Returns

• Vector: the solution, of the promoted element type of A and F.

Reverse-mode differentiation

With ChainRulesCore.jl loaded, BrambleChainRulesExt's rrule solves the adjoint system Aᵀ λ = ∂J/∂u once, reusing the forward solve's own factorisation, and returns ∂J/∂A = -λ uᵀ restricted to A's sparsity (never densified) and ∂J/∂F = λ. using Enzyme is enough for the same adjoint through BrambleEnzymeExt's native EnzymeRules rule; do not call Enzyme.@import_rrule, whose bridge drops the cotangent's explicit zeros from nzval and returns a wrong gradient without any error. Mooncake cannot represent a SparseMatrixCSC cotangent at all and is unsupported.

A closure capturing a grid space or a form needs Enzyme.Const(f) and Enzyme.set_runtime_activity(Enzyme.Reverse), as the automatic differentiation tutorial describes for every other path. A gradient with respect to a runtime scalar coefficient of the form should be a Float64 or a Ref: an Integer coefficient known only at run time costs form its inferred return type, which Enzyme's type analysis rejects.

Examples

A, F = assemble(a, l; dirichlet = :boundary => x -> 0.0)
u = pde_solve(A, F)

using SuiteSparse
u_spd = pde_solve(A, F; solver = :suitesparse, sym = :spd)

See also assemble, sparse_factorize, suitesparse_solve, suitesparse_qr_solve, accelerate_solve, mumps_solve, sparspak_solve, linear_problem.

amg_preconditioner(A::AbstractMatrix; method = :smoothed_aggregation, kwargs...) -> MultiLevel
amg_preconditioner(a::BilinearForm; method = :smoothed_aggregation, dirichlet = nothing,
                    dirichlet_components = nothing, kwargs...) -> MultiLevel

Build an algebraic multigrid hierarchy for the symmetric positive-definite matrix A – or for a assembled with the given Dirichlet conditions – ready to turn into a preconditioner with AlgebraicMultigrid.aspreconditioner.

The elliptic SBP forms Bramble assembles (Poisson, variable-coefficient diffusion, Helmholtz) have condition numbers scaling as O(h^-2): a direct sparse solve suffers severe fill-in past a few hundred thousand degrees of freedom, and an unpreconditioned Krylov method needs O(h^-1) iterations. Algebraic multigrid builds its coarse-grid hierarchy from the graph of A alone, giving a preconditioned Krylov method grid-independent, O(1) iteration counts instead.

Keywords

• method: :smoothed_aggregation (default) or :ruge_stuben, AlgebraicMultigrid's two hierarchy constructions.
• Every other keyword forwards to the chosen AlgebraicMultigrid constructor.

The BilinearForm method matches assemble(a::BilinearForm; ...)'s own default and does not symmetrize – Dirichlet rows become eₖ, but the matching columns are left alone, so A is not exactly symmetric even though the underlying operator is. AMG still builds a usable hierarchy from it, but for the SPD matrix the theory assumes, assemble through assemble(a, l; symmetrize = true) and pass that matrix to the AbstractMatrix method instead.

Requires AlgebraicMultigrid.jl; call using AlgebraicMultigrid before calling this function.

Examples

using AlgebraicMultigrid

A = assemble(a; dirichlet = :boundary)
ml = amg_preconditioner(A)
P = aspreconditioner(ml)
sol = solve(LinearProblem(A, F), KrylovJL_CG(); Pl = P)

# Or, directly from the form:
uₕ = solve(a, l; dirichlet = bcs, preconditioner = :amg, solver = KrylovJL_CG())

See also linear_problem, assemble.

ilu_preconditioner(A::AbstractMatrix; kwargs...) -> ILUZero.ILU0Precon
ilu_preconditioner(a::BilinearForm; dirichlet = nothing,
                    dirichlet_components = nothing, kwargs...) -> ILUZero.ILU0Precon

Build a zero-fill incomplete LU (ILU(0)) preconditioner for the matrix A – or for a assembled with the given Dirichlet conditions – ready to use directly as LinearSolve's Pl.

Convection-dominated forms (advection large relative to diffusion) assemble unsymmetric matrices far from an M-matrix, where algebraic multigrid (amg_preconditioner) does not degrade gracefully: measured on a 2D convection-diffusion system with diffusion 1e-2 against unit advection, ruge_stuben AMG failed to converge in 2000 GMRES iterations while ILU(0) converged in 18 (gpena/Bramble.jl#244). ILU(0) reuses A's own sparsity pattern for its factors, so it has no fill-in parameter to tune and is cheap to build, at the cost of a weaker preconditioner than a tuned incomplete factorization on harder systems.

Keywords

Forwarded to ILUZero.ilu0, which currently takes none beyond A itself.

The BilinearForm method matches assemble(a::BilinearForm; ...)'s own default and does not symmetrize; ILU(0) does not assume symmetry, so this is the natural entry point for the unsymmetric, convection-dominated forms this preconditioner targets.

Requires ILUZero.jl; call using ILUZero before calling this function.

Examples

using ILUZero

A = assemble(a; dirichlet = :boundary)
P = ilu_preconditioner(A)
sol = solve(LinearProblem(A, F), KrylovJL_GMRES(); Pl = P)

# Or, directly from the form:
uₕ = solve(a, l; dirichlet = bcs, preconditioner = :ilu0, solver = KrylovJL_GMRES())

See also amg_preconditioner, linear_problem, assemble.

sparse_factorize(A::SparseMatrixCSC; solver::Symbol = :default, sym = :auto, kwargs...) -> Factorization
sparse_factorize(a::BilinearForm; solver::Symbol = :default, sym = :auto,
                 dirichlet = nothing, dirichlet_components = nothing, kwargs...) -> Factorization

Compute the sparse direct factorization of A (or the assembled matrix of a) using the requested solver backend.

Solvers (solver)

• :default or :suitesparse: SuiteSparse (CHOLMOD for SPD/symmetric, UMFPACK for unsymmetric).
• :accelerate: Apple Accelerate native libSparse on macOS (Cholesky, LDLᵀ, LUTPP, QR).
• :mumps: MUMPS multifrontal parallel direct solver.
• :sparspak: pure-Julia sparse direct LU, zero binary dependencies.

Symmetry options

• :auto: automatically detect matrix symmetry (and diagonal positivity).
• :spd, :definite, or 1: symmetric positive definite.
• :symmetric or 2: general symmetric.
• :unsymmetric or 0: general unsymmetric.

Ignored by :sparspak, which always factors as general unsymmetric LU.

Examples

fact = sparse_factorize(A; solver = :suitesparse, sym = :spd)
u = fact \ F

See also refactor!, pde_solve, suitesparse_factorize, suitesparse_qr_factorize, accelerate_factorize, mumps_factorize, sparspak_factorize.

refactor!(fact::Factorization, A::SparseMatrixCSC) -> Factorization
refactor!(fact::Factorization, a::BilinearForm; dirichlet = nothing, dirichlet_components = nothing) -> Factorization

Recompute the numerical values of fact for updated matrix A (or assembled bilinear form a) reusing the existing symbolic factorization (fill-reducing ordering and elimination tree). A must have the exact same sparsity pattern as the matrix originally factored.

Dispatches automatically via multiple dispatch to the appropriate backend:

• SuiteSparseFactorization: updates CHOLMOD or UMFPACK numeric values.
• AccelerateFactorization: updates Apple Accelerate libSparse numeric values.
• MUMPSFactorization: updates MUMPS multifrontal numerical factorization (job = 2).
• SparspakFactorization: updates Sparspak's numeric LU values.

See also sparse_factorize, pde_solve.

SuiteSparseFactorization{T} <: Factorization{T}

Wrapper type representing a factorized SuiteSparse linear system (CHOLMOD Cholesky or UMFPACK LU), supporting in-place solves via LinearAlgebra.ldiv!, back-substitution via \, and non-allocating numeric refactoring via suitesparse_refactor!.

suitesparse_factorize(A::AbstractMatrix; sym = :auto, kwargs...) -> SuiteSparseFactorization
suitesparse_factorize(a::BilinearForm; dirichlet = nothing, dirichlet_components = nothing,
                      symmetrize = false, sym = :auto, kwargs...) -> SuiteSparseFactorization

Compute the sparse direct factorization of A (or the assembled matrix of a) using SuiteSparse (CHOLMOD for symmetric positive-definite systems, UMFPACK for unsymmetric).

Symmetry options

• :auto (default): automatically detects symmetry. If A is symmetric and positive definite (or symmetrize = true), uses CHOLMOD sparse Cholesky. Otherwise, uses UMFPACK sparse LU.
• :spd, :definite, or 1: symmetric positive definite (CHOLMOD Cholesky).
• :symmetric or 2: symmetric factorization.
• :unsymmetric or 0: general unsymmetric (UMFPACK LU).

Requires SuiteSparse.jl; call using SuiteSparse before calling this function.

Examples

using SuiteSparse

A, F = assemble(a, l; dirichlet = :boundary => x -> 0.0, symmetrize = true)
fact = suitesparse_factorize(A; sym = :spd)
u = fact \ F

See also suitesparse_solve, suitesparse_refactor!, pde_solve, assemble.

suitesparse_solve(A::SparseMatrixCSC, F::AbstractVector; sym = :auto, kwargs...) -> Vector
suitesparse_solve(a::BilinearForm, l::LinearForm; dirichlet = nothing, dirichlet_components = nothing,
                  symmetrize = false, sym = :auto, kwargs...) -> VectorElement

Directly solve A u = F (or assemble(a, l) system) using SuiteSparse factorization.

Requires SuiteSparse.jl; call using SuiteSparse before calling this function.

See also suitesparse_factorize, refactor!, pde_solve.

suitesparse_refactor!(fact::SuiteSparseFactorization, A::SparseMatrixCSC) -> SuiteSparseFactorization

Recompute the numeric factorization of A inside fact reusing the existing symbolic factorization (fill-reducing ordering and elimination tree). A must have the exact same sparsity pattern as the matrix originally factored.

See also refactor!, suitesparse_factorize.

suitesparse_qr_factorize(A::AbstractMatrix; tol = ..., ordering = ..., kwargs...) -> QRSparse
suitesparse_qr_factorize(a::BilinearForm; dirichlet = nothing, dirichlet_components = nothing,
                         kwargs...) -> QRSparse

Compute the sparse direct QR factorization of A (or the assembled matrix of a) using SuiteSparse's SPQR, suited to overdetermined least-squares systems and the rectangular blocks of a constrained saddle-point form. Unlike suitesparse_factorize, A need not be square, and this needs only SparseArrays – already a dependency of Bramble, no using SuiteSparse required.

kwargs (tol, ordering, ...) are forwarded to SparseArrays.SPQR.qr.

Examples

A, F = assemble(a, l; dirichlet = :boundary => x -> 0.0)
fact = suitesparse_qr_factorize(A)
u = fact \ F

See also suitesparse_qr_solve, suitesparse_factorize, pde_solve.

suitesparse_qr_solve(A::SparseMatrixCSC, F::AbstractVector; tol = ..., ordering = ..., kwargs...) -> Vector
suitesparse_qr_solve(a::BilinearForm, l::LinearForm; dirichlet = nothing, dirichlet_components = nothing,
                    kwargs...) -> VectorElement

Directly solve A u = F (exactly if square, least-squares if overdetermined) using SuiteSparse's SPQR sparse QR factorization. Needs only SparseArrays; no using SuiteSparse required.

See also suitesparse_qr_factorize, pde_solve.

AccelerateFactorization{T} <: Factorization{T}

Wrapper type representing a factorized Apple Accelerate (libSparse) linear system, supporting in-place solves via LinearAlgebra.ldiv!, back-substitution via \, and non-allocating symbolic reuse via accelerate_refactor!. Available only on macOS.

accelerate_factorize(A::SparseMatrixCSC; sym = :auto, kind = :auto, kwargs...) -> AccelerateFactorization
accelerate_factorize(a::BilinearForm; dirichlet = nothing, dirichlet_components = nothing,
                     symmetrize = false, sym = :auto, kind = :auto, kwargs...) -> AccelerateFactorization

Compute the sparse direct factorization of A (or the assembled matrix of a) using Apple Accelerate's native libSparse on macOS.

Symmetry and factorization options

• sym = :auto (default): automatically detects symmetry. If A is symmetric with strictly positive diagonal, uses Cholesky (SparseFactorizationCholesky). If symmetric, uses $LDL^T$ (SparseFactorizationLDLT). Otherwise uses threshold partial pivoting LU (SparseFactorizationLUTPP).
• sym = :spd, :definite, or 1 (or kind = :cholesky): symmetric positive definite Cholesky.
• sym = :symmetric or 2 (or kind = :ldlt): symmetric indefinite $LDL^T$.
• sym = :unsymmetric or 0 (or kind = :lu / :lutpp): general unsymmetric LU with threshold partial pivoting.
• kind = :qr: sparse QR factorization.

Requires macOS and AppleAccelerate.jl; call using AppleAccelerate before calling this function.

Examples

using AppleAccelerate

A, F = assemble(a, l; dirichlet = :boundary => x -> 0.0, symmetrize = true)
fact = accelerate_factorize(A; sym = :spd)
u = fact \ F

See also accelerate_solve, accelerate_refactor!, pde_solve, assemble.

accelerate_factorize(A::AbstractMatrix; sym = :auto, kind = :auto, kwargs...) -> LinearAlgebra.Factorization

Dense factorization of A, via ordinary LinearAlgebra.lu, LinearAlgebra.cholesky, or LinearAlgebra.qr – not Apple Accelerate's libSparse used by the SparseMatrixCSC method above, and not gated on macOS.

AppleAccelerate.jl (checked against v0.7.0) defines no dense lu, cholesky, getrf, potrf, gemm, or other LAPACK/BLAS bindings for Bramble to call. What its __init__ does is register Accelerate with libblastrampoline (BLAS.lbt_forward), so plain LinearAlgebra.lu/cholesky/qr are already running on Accelerate's BLAS/LAPACK the moment using AppleAccelerate has been evaluated, on every call site in Bramble or anywhere else, with no Bramble dispatch involved. This method adds no acceleration beyond that; it exists so dense callers can spell accelerate_factorize with the same sym/kind vocabulary as the sparse methods above. kwargs... is accepted but unused.

Symmetry and factorization options

• sym = :auto (default): LinearAlgebra.cholesky if A is symmetric positive definite, otherwise LinearAlgebra.lu.
• sym = :spd, :definite, or 1 (or kind = :cholesky): LinearAlgebra.cholesky.
• sym = :unsymmetric or 0 (or kind = :lu): LinearAlgebra.lu.
• kind = :qr: LinearAlgebra.qr.
• sym = :symmetric, 2, or kind = :ldlt: not implemented here. The dense analogue of a sparse LDLᵀ factorization is LinearAlgebra.bunchkaufman, which this method does not wrap; call it directly.

See also accelerate_solve.

accelerate_solve(A::SparseMatrixCSC, F::AbstractVector; sym = :auto, kind = :auto, kwargs...) -> Vector
accelerate_solve(a::BilinearForm, l::LinearForm; dirichlet = nothing, dirichlet_components = nothing,
                 symmetrize = false, sym = :auto, kind = :auto, kwargs...) -> VectorElement

Directly solve A u = F (or assemble(a, l) system) using Apple Accelerate direct factorization.

Requires macOS and AppleAccelerate.jl; call using AppleAccelerate before calling this function.

See also accelerate_factorize, refactor!, pde_solve.

accelerate_refactor!(fact::AccelerateFactorization, A::SparseMatrixCSC) -> AccelerateFactorization

Recompute the numeric factorization of A stored in fact reusing the existing symbolic factorization (fill-reducing ordering and sparsity analysis). A must have the exact same sparsity pattern as the matrix originally factored.

See also refactor!, accelerate_factorize.

MUMPSFactorization{T} <: Factorization{T}

Wrapper type representing a factorized MUMPS linear system, supporting in-place solves via LinearAlgebra.ldiv!, back-substitution via \, and non-allocating reuse for transient PDE time stepping.

mumps_factorize(A::AbstractMatrix; sym = :auto, icntl = nothing, cntl = nothing, kwargs...) -> MUMPSFactorization
mumps_factorize(a::BilinearForm; dirichlet = nothing, dirichlet_components = nothing,
                symmetrize = false, sym = :auto, kwargs...) -> MUMPSFactorization

Compute the sparse direct multifrontal factorization of A (or the assembled matrix of a) using MUMPS.jl.

Symmetry options

• :auto (default): automatically detects symmetry. If A is symmetric and positive definite (or symmetrize = true), uses symmetric positive-definite factorization (sym = 1). If symmetric, uses general symmetric $LDL^T$ (sym = 2). Otherwise, uses unsymmetric LU (sym = 0).
• :spd, :definite, or 1: symmetric positive definite (Cholesky / $LL^T$).
• :symmetric or 2: general symmetric ($LDL^T$ with Bunch-Kaufman pivoting).
• :unsymmetric or 0: general unsymmetric LU.

Control parameters

• icntl: Optional dictionary or collection of pairs of integer control parameters (e.g., 7 => 1 for user/METIS ordering, 14 => 30 for memory relaxation).
• cntl: Optional dictionary or collection of pairs of real control parameters (e.g., 1 => 0.01 for numerical pivoting threshold).

Requires MUMPS.jl; call using MUMPS before calling this function.

Examples

using MUMPS

A, F = assemble(a, l; dirichlet = :boundary => x -> 0.0, symmetrize = true)
fact = mumps_factorize(A; sym = :spd)
u = fact \ F

See also pde_solve, mumps_solve, assemble.

mumps_solve(A::SparseMatrixCSC, F::AbstractVector; sym = :auto, kwargs...) -> Vector
mumps_solve(a::BilinearForm, l::LinearForm; dirichlet = nothing, dirichlet_components = nothing,
            symmetrize = false, sym = :auto, kwargs...) -> VectorElement

Directly solve A u = F (or assemble(a, l) system) using MUMPS direct factorization.

Symmetry options

• :auto (default): automatically detects symmetry.
• :spd, :definite, or 1: symmetric positive definite.
• :symmetric or 2: general symmetric.
• :unsymmetric or 0: general unsymmetric.

Requires MUMPS.jl; call using MUMPS before calling this function.

Examples

using MUMPS

A, F = assemble(a, l; dirichlet = :boundary => x -> 0.0)
u = mumps_solve(A, F)

See also mumps_factorize, refactor!, pde_solve.

mumps_refactor!(fact::MUMPSFactorization, A::SparseMatrixCSC) -> MUMPSFactorization

Recompute the numeric factorization of A inside fact reusing the existing symbolic factorization (fill-reducing analysis and ordering). A must have the exact same sparsity pattern as the matrix originally factored.

See also refactor!, mumps_factorize.

SparspakFactorization{T} <: Factorization{T}

Wrapper type representing a factorized Sparspak.jl sparse LU system, supporting in-place solves via LinearAlgebra.ldiv!, back-substitution via \, and non-allocating reuse for transient PDE time stepping.

sparspak_factorize(A::AbstractMatrix) -> SparspakFactorization
sparspak_factorize(a::BilinearForm; dirichlet = nothing, dirichlet_components = nothing) -> SparspakFactorization

Compute the pure-Julia sparse direct LU factorization of A (or the assembled matrix of a) using Sparspak.jl.

Sparspak has no binary dependency, so it factors matrices whose entries are not Float64/ComplexF64 – Float32, BigFloat, or a ForwardDiff.Dual – where SuiteSparse and MUMPS cannot.

Requires Sparspak.jl; call using Sparspak before calling this function.

Examples

using Sparspak

A, F = assemble(a, l; dirichlet = :boundary => x -> 0.0)
fact = sparspak_factorize(A)
u = fact \ F

See also sparspak_solve, sparspak_refactor!, pde_solve, assemble.

sparspak_solve(A::SparseMatrixCSC, F::AbstractVector) -> Vector
sparspak_solve(a::BilinearForm, l::LinearForm; dirichlet = nothing, dirichlet_components = nothing,
               symmetrize = false) -> VectorElement

Directly solve A u = F (or assemble(a, l) system) using Sparspak's pure-Julia sparse direct LU factorization.

Requires Sparspak.jl; call using Sparspak before calling this function.

Examples

using Sparspak

A, F = assemble(a, l; dirichlet = :boundary => x -> 0.0)
u = sparspak_solve(A, F)

See also sparspak_factorize, sparspak_refactor!, pde_solve.

sparspak_refactor!(fact::SparspakFactorization, A::SparseMatrixCSC) -> SparspakFactorization

Recompute the numeric factorization of A inside fact reusing the existing symbolic factorization (fill-reducing ordering). A must have the exact same sparsity pattern as the matrix originally factored.

See also refactor!, sparspak_factorize.

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 – the fix diffusion_matrix-style Newton residuals in the nonlinear worked examples name and deliberately defer, since 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.

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 (a few KB, 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 plain, unlocked Dict, sized for the one-cache-per-residual usage above. A form assembled from more than one task needs a lock or a per-task cache, the same as any other shared mutable Dict.