Forms#

form(Wₕ, f) -> LinearForm

Construct a LinearForm over the test space Wₕ using the linear expression f.

Construction resolves the AST once and runs simplify_ast over it – dropping zero terms, folding nested constant scalings, lifting scalars out of inner products, and distributing sums that mix components – before it is stored; a sum is otherwise kept as written. Grid partitioning for parallel assembly is determined from the resolved AST during assembly (see _term_colour_strides).

Every SourceFunction reachable from the simplified AST – a source term built directly from a plain function, f(x), rather than a VectorElement or a Ref – is then sampled once against its own leaf's space and lowered to a SourceVector because a brand-new closure passed as a source would force a full recompilation of the assembly pipeline on every distinct closure, since Julia gives it its own type; assembling against the fixed SourceVector shape instead means that cost is paid once, here, not on every later assemble!/assemble call. This changes what a raw closure that captures mutable state does: it is evaluated once, now, not re-evaluated on later assemblies. VectorElement and Ref coefficients are unaffected – see the note on _lower_sources in ast/common.jl for why, and for the documented alternative (update_coefficients!) a source meant to keep varying should use instead.

The one exception is an interpolant πₕ(uₕ), which is never sampled: every assemble/assemble! evaluates it on uₕ's current values, and throws an ArgumentError once uₕ's mesh has moved.

A source under restrict_to is also called once here, at the first point of its region that the assembly reads (inside every enclosing region, moved by every enclosing shift), so that outside its region it answers a zero of the type it has inside and its stencil keeps one type at every point. It is never called at a point no assembly reads, with one exception: the point comes from the markers as they stand at form, so after markers! moves the region the source has been called once at a point of the old one. A region with no such point at form keeps the empty stencil. In return:

• a source that throws at a point of its region throws here, not at the first assembly;

• a coefficient that is not finite outside the region gives NaN there when it scales a shift, a difference, an average or a sum holding the restricted source (directly over the restriction, or over restrictions and scales nested in it, it still gives 0);

• assemble_add! into entries holding -0.0 may turn them into +0.0 there.

Examples

using Bramble
Wₕ = gridspace(mesh(domain(interval(0.0, 1.0)), 11))
fₕ = Rₕ(Wₕ, x -> 1.0)
l = form(Wₕ, v -> innerₕ(fₕ, v))    # l(v) = (fₕ, v)ₕ
isapprox(sum(assemble(l)), 1.0; atol = 1.0e-12)

# output
true

form(Wₕ, Vₕ, f) -> BilinearForm

Construct a BilinearForm over the trial space Wₕ and the test space Vₕ from the bilinear expression f (a function of trial and test arguments (u, v)). The AST is resolved once and run through simplify_ast – dropping zero terms, folding nested constant scalings, lifting scalars out of inner products, and distributing sums that mix components – before it is stored; a sum is otherwise kept as written.

An expression may hold one rank-one term, a scalar times a functional of the trial function times a functional of the test function (σ * innerₕ(g, u) * innerₕ(h, v), g and h data). Its matrix σ wₛ wₜᵀ is dense, so the term is kept apart from the stencil: form then returns a RankOneForm, which assemble turns into a RankOneUpdate. A functional of one argument alone, a product of two functionals of the same argument, or two rank-one terms throw an ArgumentError.

Examples

using Bramble: D₋ₓ, inner₊ₓ
# a(u, v) = (∇ₕu, ∇ₕv)₊
a = form(Wₕ, Wₕ, (u, v) -> inner₊ₓ(D₋ₓ(u), D₋ₓ(v)))

# a coupled system, one term per block
a = form(Vₕ, Vₕ, (u, v) -> inner₊ₓ(D₋ₓ(u(1)), D₋ₓ(v(1))) + innerₕ(u(2), v(1)))

# a penalty on the discrete mean, (∇ₕu, ∇ₕv)₊ + σ (1, u)ₕ (1, v)ₕ
a = form(Wₕ, Wₕ, (u, v) -> inner₊(∇ₕ(u), ∇ₕ(v)) + σ * innerₕ(1, u) * innerₕ(1, v))

form(Wₕ::CompositeGridSpace, Vₕ::CompositeGridSpace, f) -> GhostCompositeBilinearForm
form(Vₕ::CompositeGridSpace, f) -> GhostCompositeLinearForm

Forms on composites of ghost spaces (C₁ × C₂, C^Val(n), or constrained(C₁ × C₂, relations)): the form on their padded spaces, reduced through the closure, block diagonal on a product, coupling the components on a composite constrained as a whole wherever a relation reads several of them. f addresses the components as on a plain composite, u(1), v(2). assemble(a), assemble(a, l), assemble! and their dirichlet and dirichlet_components keywords act as on a scalar ghost space, the rows of each selected component replaced through its physical mesh. a(u, v), l(v) and evaluate! give vᴴ A u and Fᵀ v as on a scalar ghost space.

Every component of both spaces must ghost the same facets of meshes of one size; an ArgumentError names the fix otherwise. Components built by polygon_spaces pair with polygon spaces of one mesh: the trial components of a form are read with the reach guard and the conflict moves of a scalar polygon form, and a trial space is refused as a test side. A rank-one term, a scalar space on the other side, semidiscretize, and the matrix-free, Kronecker and multigrid operators are refused with an ArgumentError, as is a(u, v) on a form reading a conflict node of a polygon trial component.

expression(form::Union{LinearForm,BilinearForm}) -> String

Render form's resolved AST (form.ast) in Bramble's own operator notation.

Left untyped on purpose: LinearForm/BilinearForm are defined in assembly/linear.jl and assembly/bilinear.jl, included well after this file, so a Union{LinearForm,BilinearForm} annotation here would need those names to exist at include time. Every LazyOp node has its own, more specific expression method, so this untyped fallback only ever catches forms.

dirac(x0::Union{Real, NTuple{D, Real}, AbstractVector{<:Real}}, strength = 1) -> DiracSource{D}
dirac(points::AbstractVector, strengths = 1) -> DiracSource{D}

Construct a symbolic point (Dirac delta) source term at coordinate x0 with scalar strength, or a collection of point sources at points with corresponding strengths.

Tuples are the type-stable spelling of a point: dirac((0.3, 0.4)) and dirac([(0.3, 0.4), (0.7, 0.2)]) infer a concrete DiracSource{2}. A vector point is accepted, but its dimension is read from its length at run time, so dirac([0.3, 0.4]) and dirac([[0.3, 0.4], [0.7, 0.2]]) are not inferable; they normalise to the same tuples and assemble the same load vector. A vector with a real element type ([0.2, 0.7], Real[0.2, 0.7]), or an untyped Any[0.2, 0.7] holding only reals, is one point, here the 2D point (0.2, 0.7); two 1D points are spelt dirac([(0.2,), (0.7,)]). A list of points may mix tuples and vectors.

The strength differentiates under ForwardDiff and Enzyme. The location is converted to Float64, so ForwardDiff cannot differentiate it and throws, while Enzyme differentiates it in forward and reverse mode.

Mathematical formulation

A point source represents the continuous linear functional:

When assembled into a LinearForm via innerₕ or inner₊:

l = form(Wₕ, v -> innerₕ(dirac(x0, strength), v))

the strength \(S\) is distributed over the \(2^D\) bounding cell corners using multilinear interpolation weights \(w_c\):

• On-grid source (\(x_0 = x_I\)): Exactly 1 nonzero entry \(b_I = S\).

• Off-grid source: Load distributes locally to the surrounding cell vertices, preserving total integral \(\sum b_i = S\).

• Contraction: Contracting the linear form against a discrete grid function \(v_h \in W_h\) evaluates \(\ell(v_h) = b \cdot v_h \approx S \cdot v_h(x_0)\) with \(\mathcal{O}(h^2)\) accuracy.

Arguments

• x0: Point coordinate, a scalar number (1D), NTuple{D, Real} (type-stable), or AbstractVector{<:Real} (one D-dimensional point, not inferable). An empty tuple or vector throws an ArgumentError.

• strength: Source intensity (default: 1, which takes the space's element type). Accepts any of the following.

  • A constant Number (e.g. 2.5).

  • A dynamic Ref(val), which enables live in-place updates (strength[] = new_val) in time-stepping loops without rebuilding the form and with 0 heap allocations.

  • A zero-argument function thunk, () -> f(t), for time-dependent sources.

• points: Collection of point coordinates for multiple simultaneous sources. The first point sets the dimension D; a later point with a different number of coordinates, or of a kind that cannot have D coordinates, throws an ArgumentError naming its index.

• strengths: Matching collection of intensities, or a single scalar broadcast to all points.

Examples

using Bramble

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

# 1. On-grid point source
l_on = form(Wₕ, v -> innerₕ(dirac(0.3, 2.0), v))
b_on = assemble(l_on)
sum(b_on) ≈ 2.0

# 2. Dynamic time-stepping point source (0 allocations)
s_live = Ref(1.0)
l_live = form(Wₕ, v -> innerₕ(dirac(0.45, s_live), v))
b = assemble(l_live)

for step in 1:10
    s_live[] = sin(0.1 * step)
    assemble!(b, l_live) # 0 bytes allocated
end

See also: DiracSource, innerₕ, inner₊, form, assemble!.

assemble(form::LinearForm; dirichlet = nothing, dirichlet_components = nothing) -> AbstractVector

Assemble the system vector of the LinearForm, applying the boundary values dirichlet carries on the regions it names. dirichlet accepts a label => f Pair, a Tuple of such Pairs, or constraints from dirichlet_constraints – a bare label or Tuple of labels has no values to write and raises an error (see _normalize_dirichlet for every accepted form). dirichlet_components restricts which leaf components of a composite test space they bind to (see dirichlet_bc!).

Runs serially or across threads following test_space(form)'s backend execution_policy: Serial (the default) or Parallel. assemble_parallel! always threads regardless of the backend policy.

assemble(form::BilinearForm; dirichlet = nothing, dirichlet_components = nothing) -> AbstractMatrix

Allocate a matrix with the form's sparsity pattern and assemble into it. The matrix type is matrix_type(backend(test_space(form))) – SparseMatrixCSC{Float64,Int} by default, or whatever backend the space's mesh was built with.

Call this once, then assemble into what it returns. The first call builds the sparsity pattern, the larger part of the work. A repeated call copies the pattern the form recorded and allocates about the matrix it returns. It re-records after the mesh markers change, after another matrix is filled, when a Dirichlet row inserted a missing diagonal, and on every call for a matrix type other than SparseMatrixCSC{T, Int}; a refill allocates nothing. A time loop or Newton iteration benefits from preallocating the pattern once:

A = assemble(a)                        # once: pattern, allocation and initial fill
for step in 1:nsteps
    Rₕ!(cₕ, coefficient_at(step))      # modified in-place
    assemble!(A, a)                    # or assemble_parallel!(A, a)
end

Runs serially or across threads following form.trial_space's backend execution_policy: Serial by default. Optional dirichlet applies boundary conditions to the matrix – a label Symbol, a Tuple of labels, a label => f Pair (values ignored on the matrix side), or constraints from dirichlet_constraints; see _normalize_dirichlet for every accepted form. dirichlet_components restricts which leaf components of a composite test space they bind to (see dirichlet_bc!).

assemble(a::BilinearForm, l::LinearForm; dirichlet = nothing, dirichlet_components = nothing, symmetrize::Bool = false) -> (A, F)

Assemble both the matrix and the vector in one call, applying the same dirichlet constraints to each – the common case of building a system to solve A \ F against. a and l should share the same test space; boundary values are read from dirichlet once rather than once per form.

symmetrize = true additionally restores symmetry in A after dirichlet_bc! (which only clears rows, not the columns) and updates F to match, via symmetrize!. Requires dirichlet to name at least one label – there is nothing to symmetrize against otherwise.

A, F = assemble(a, l; dirichlet = :boundary => sol, symmetrize = true)
u = A \ F

is the one-call equivalent of assembling a and l separately and calling symmetrize! by hand.

assemble(a::RankOneForm; dirichlet = nothing, dirichlet_components = nothing) -> RankOneUpdate

Assemble a form holding a rank-one term into a RankOneUpdate: its other terms as assemble assembles a BilinearForm, with dirichlet applied to their matrix, and the two functionals' vectors through the linear path. wₛ is zero at the Dirichlet rows, so they are identity rows of the result.

assemble(a::RankOneForm, l::LinearForm; dirichlet = nothing, dirichlet_components = nothing, symmetrize::Bool = false) -> (K, F)

The rank-one counterpart of assemble(a::BilinearForm, l::LinearForm): K a RankOneUpdate and F the load, with the same dirichlet constraints applied to both. symmetrize = true runs symmetrize! on them.

assemble(a::RankOneForm, l::GhostLinearForm; dirichlet = nothing, dirichlet_components = nothing, symmetrize::Bool = false) -> (K, F)
assemble!(K::RankOneUpdate, F, a::RankOneForm, l::GhostLinearForm; dirichlet = nothing, dirichlet_components = nothing, symmetrize::Bool = false) -> (K, F)

A form holding a rank-one term and a load on ghost spaces, assembled together as for a ghost form without one: K the RankOneUpdate A + σ wₛ wₜᵀ reduced through the closures (E, d), and the load E_vᵀ(F_pad − A_pad d) − σ wₛ (wₜ_padᵀ d), the rank-one part acting on the padded unknown E u + d. assemble! refills the outputs of an earlier assemble(a, l) in place, allocating nothing. dirichlet replaces the rows of A and zeroes wₛ there, and symmetrize then moves the Dirichlet columns into the load, as on a plain space. Both, and spaces whose relation data depend on time, are refused where the ghost form without the term refuses them.

assemble!(b::AbstractVector, form::LinearForm; dirichlet = nothing,
          dirichlet_components = nothing) -> AbstractVector

Refill b with the assembled form and return it with zero allocations (0 bytes). A product with three or more region-restricted summands on its test side is split by form into products of at most two (gpena/Bramble.jl#727); one it leaves whole, as it would split into more than 16 or its constant scales keep several types, allocates per grid point.

On a test space whose backend lives on a device (a Metal space, say), b is filled on the host and uploaded in one copyto!, as system matrices are: that path allocates a host buffer and host mirrors of the space on every call, so the zero-allocation guarantee is the host path's alone.

assemble! uses the pre-resolved form.ast stored directly inside the form.

Live coefficients

• Grid functions: the stored AST retains references to source VectorElement storage. Mutating values in-place (Rₕ!(uₕ, ...) or parent(uₕ) .= ...) between steps automatically updates the assembled vector without needing to rebuild the form. Nested scales such as uₕ * (wₕ * v) stay live too, because each grid function is read at assembly time and never fused into a copy when the form is built.

• Dynamic scalars: plain numbers work directly for constant scalars. To update a scalar dynamically across loop iterations, wrap it in a Ref(val) (e.g. α = Ref(1.0); l = form(Wₕ, v -> α * innerₕ(uₕ, v))). Mutating α[] = new_val evaluates live during assembly with 0 allocations.

Arguments

• b: Vector to refill, with length ndofs(test_space(form)); any other length raises a DimensionMismatch before anything is written.

• form: Linear form to assemble.

Keywords

• dirichlet: Boundary values to impose after assembly, and the labels to impose them on, together – a label => f Pair, a Tuple of such Pairs, or constraints from dirichlet_constraints (default: nothing; see _normalize_dirichlet for every accepted form).

• dirichlet_components: Restricts which leaf components of a composite test_space(form) the labels bind to (see dirichlet_bc!; default: nothing, targeting all leaves).

Runs serially or across threads following test_space(form)'s backend execution_policy: Serial or Parallel. assemble_parallel! forces threaded execution regardless of backend policy.

See also assemble, assemble_parallel!, and evaluate!.

assemble!(A::AbstractMatrix, form::BilinearForm; dirichlet = nothing, dirichlet_components = nothing) -> AbstractMatrix

Assemble the BilinearForm into the preallocated matrix A, allocating nothing (0 bytes) when A is the default SparseMatrixCSC. A product with three or more region-restricted summands is split by form into products of at most two (gpena/Bramble.jl#727); one it leaves whole, as it would split into more than 16 or its constant scales keep several types, allocates per grid point.

Runs serially or across threads following form.trial_space's backend execution_policy: Serial (the default) or Parallel. assemble_parallel! is a separate, lower-level entry point that always threads, ignoring the backend's policy.

assemble! uses the pre-resolved form.ast stored directly inside the form.

Live coefficients

• Grid functions. The stored AST retains references to source VectorElement storage. Mutating values in-place (Rₕ!(cₕ, ...) or parent(cₕ) .= ...) between steps automatically updates the matrix entries. Nested scales such as uₕ * (wₕ * v) stay live too: each grid function is read at assembly time, never fused into a copy when the form is built. On the host this costs 0 allocations; on a device-backed space each fill copies the coefficient to the host anew, so it stays live but is not allocation-free there.

• Dynamic scalars. Plain numbers work directly for constant scalars. To update a scalar dynamically across loop iterations, wrap it in a Ref(val) (e.g. β = Ref(1.0); a = form(Wₕ, Wₕ, (u, v) -> innerₕ(β * D₋ₓ(u), D₋ₓ(v)))). Mutating β[] = new_val evaluates live during assembly with 0 allocations on the host; the Ref itself adds nothing, but a device-backed space's call still allocates regardless of the scalar.

assemble!(A, F, a::GhostBilinearForm, l::GhostLinearForm;
    dirichlet = nothing, dirichlet_components = nothing, symmetrize = false)

Refill A and F with the matrix and the load that A, F = assemble(a, l) returns, and return A, F. A keeps the pattern and F the element type that assemble(a, l) gave, so call that once and refill its two outputs. The call allocates nothing warm.

Both forms are assembled together because, when the closure of the trial space has an affine part d ≠ 0, the load is E_vᵀ(F − A d) and reads the matrix. assemble!(b, l) alone refuses a test space with d ≠ 0, which is exact when both sides are one space: on a pair whose trial space has data and whose test space has none, it gives E_vᵀ F without the lift. dirichlet, dirichlet_components and symmetrize act as in assemble(a, l).

It throws an ArgumentError when the trial space of a or the test space of l is time-dependent (pose the problem with semidiscretize), for symmetrize = true without dirichlet, for A that is not a SparseMatrixCSC with the reduced pattern of a, for a dirichlet the rows of the space do not allow, and for an F whose element type cannot hold the load's. A DimensionMismatch follows an F of the wrong length and forms whose test spaces differ: not the same space, and meshes with other points or closures (E, d) that are not equal.

assemble!(K::RankOneUpdate, a::RankOneForm; dirichlet = nothing, dirichlet_components = nothing) -> RankOneUpdate

Refill K, from an earlier assemble(a), in place: A, both vectors and σ, allocating nothing on the host. A Ref factor of σ and every coefficient of the functionals are read afresh.

assemble_parallel!(b::AbstractVector, form::LinearForm) -> AbstractVector

Refill b with the assembled form across threads and return it, regardless of test_space(form)'s backend execution policy. Unlike assemble!, does not apply Dirichlet conditions. A stencil reaching neighbours is swept by colour, so the entries can differ from the serial sweep in the last bits; repeated threaded calls agree exactly. On a device test space the threaded sweep runs on a host mirror into a host buffer, uploaded to b in one copyto!, as in assemble!.

assemble_parallel!(A::AbstractMatrix, form::BilinearForm) -> AbstractMatrix

Refill A with the assembled form across threads and return it, regardless of form.trial_space's backend policy. A must already carry the correct sparsity pattern from allocate_system_matrix or a previous assemble. Unlike assemble!, does not apply dirichlet_labels.

Colouring on the test side ensures thread safety when updating stored matrix values concurrently. The band-coloured sweep runs for every matrix type: it reaches storage only through the matrix-type seam (_scatter_position/_scatter_add!), so a dense, banded or CSR backend colours exactly as SparseMatrixCSC does (gpena/Bramble.jl#12, #190).

Shares assemble!'s record/replay cache (gpena/Bramble.jl#338): once the form has recorded where each entry of A lives – A came from assemble or allocate_system_matrix on a serial form, or an earlier fill into this same A recorded it – every sweep writes through those positions instead of searching for them. The first fill into any other matrix object records once, serially, then replays across threads. A device-resident matrix replays too, through positions recorded against its host mirror, so only its recording fill searches (gpena/Bramble.jl#318).

allocate_system_matrix(form::BilinearForm, ast = resolve_form_ast(form)) -> AbstractMatrix

Build the matrix a BilinearForm assembles into: the appropriate size, correct sparsity pattern, and stored zeros throughout, in matrix_type(backend(test_space(form))) – SparseMatrixCSC{Float64,Int} by default, or whatever backend the space's mesh was built with (see _allocate_from_pattern). On a metal_backend space this returns a device-resident BrambleMetalExt.MetalSparseMatrixCSR{Float32,Int64}, not a host SparseMatrixCSC – verified against real triplets (gpena/Bramble.jl#94).

The pattern follows from the stencil rather than coefficient values, remaining invariant while the mesh and expression structure are unchanged. Preallocating the matrix once outside loops allows zero-allocation in-place assembly:

using Bramble: allocate_system_matrix
A = allocate_system_matrix(a)
for step in 1:nsteps
    assemble!(A, a)          # refills values in-place with zero allocations
end

Only the structure is preallocated here; all stored entries are zero until assemble! fills them.

See also assemble and assemble!.

evaluate!(scratch::AbstractVector, form::LinearForm, vₕ::VectorElement) -> Number

Evaluate form at vₕ, assembling into scratch rather than into a newly allocated vector, and return the resulting contracted scalar value.

Useful when both the assembled vector and the scalar value are needed (such as a Newton step requiring the residual vector and its norm). scratch is overwritten and contains the right-hand side upon return.

For the scalar value alone, form(vₕ) fuses the contraction into the assembly sweep in a single pass.

Examples

using Bramble: evaluate!
l = form(Wₕ, v -> innerₕ(fₕ, v))
scratch = zeros(ndofs(Wₕ))
for step in 1:nsteps
    Rₕ!(fₕ, source_at(step))          # modified in-place
    value = evaluate!(scratch, l, uₕ)
end

evaluate!(scratch::AbstractVector, l::GhostLinearForm, vₕ::VectorElement) -> Number
evaluate!(scratch::AbstractVector, l::GhostCompositeLinearForm, vₕ::VectorElement) -> Number

Write the reduced load assemble(l) into scratch, of length ndofs(test_space(l)), and return its contraction with vₕ, the value l(vₕ). It refuses as l(vₕ) does: with an ArgumentError naming assemble(a, l) on a test space whose closure has an affine part d ≠ 0.

RankOneUpdate{T, M, VS, VT} <: AbstractMatrix{T}
RankOneUpdate(A, σ, ws, wt)
RankOneUpdate{T}(A, σ, ws, wt)

The matrix A + σ wₛ wₜᵀ, stored as its four parts and never as a dense block: what assemble returns for a form holding a rank-one term σ * innerₕ(g, u) * innerₕ(h, v). A is the matrix of the form's other terms, wₛ the vector of the test functional innerₕ(h, ·) and wₜ that of the trial functional innerₕ(g, ·), each assembled by the linear path (assemble(form(Vₕ, v -> innerₕ(h, v)))).

Fields

• A::M: the sparse part, an AbstractMatrix of size (length(ws), length(wt)).

• σ::T: the scalar, read when the form is assembled; a Ref factor stays live, so assemble! reads its new value.

• ws::VS: the test-side (row) vector.

• wt::VT: the trial-side (column) vector.

assemble!(K, a) refills A, ws, wt and σ in place, allocating nothing. RankOneUpdate(A, σ, ws, wt) takes T from the four parts' element types; σ is converted to T. Not exported: write Bramble.RankOneUpdate.

Entry (i, j) is A[i, j] + σ * (ws[i] * wt[j]). mul!(y, K, x, α, β) computes α (A x + σ ws (wtᵀ x)) + β y with one product by A and two passes over the vectors, and allocates nothing; so do K * x, K * X and mul! on matrices, apart from their result. Matrix(K) is the dense matrix. K' and transpose(K) are LinearAlgebra's lazy Adjoint and Transpose of K: they read A, σ and the vectors from K, so a refill of K is seen through them, and their products allocate nothing either. similar(K) is a rank-one update with similar parts, and copy(K) one with copied parts.

With dirichlet labels, assemble makes each Dirichlet row of A an identity row, as dirichlet_bc! does, and zeroes ws there, so the row is the identity row of K; the columns keep their entries, as A's do. symmetrize! then moves those columns' entries, of A and of σ wₛ wₜᵀ, into the load and zeroes wt at the Dirichlet indices.

issymmetric(K) is true when A is symmetric and ws == wt (or σ is zero), the case of a symmetric form with one functional on both sides. Any other K reports false without comparing its N² entries, even if they happen to be symmetric (ws a multiple of wt, say). ishermitian(K) is the same test for a real K, and for a complex one asks for a Hermitian A, a real σ and ws == conj.(wt).

Examples

using Bramble, LinearAlgebra
Wₕ = gridspace(mesh(domain(interval(0.0, 1.0)), 6, false))
a = form(Wₕ, Wₕ, (u, v) -> inner₊(∇ₕ(u), ∇ₕ(v)) + 2.0 * innerₕ(1, u) * innerₕ(1, v))
K = assemble(a)
w = assemble(form(Wₕ, v -> innerₕ(1, v)))
A = assemble(form(Wₕ, Wₕ, (u, v) -> inner₊(∇ₕ(u), ∇ₕ(v))))
(K isa Bramble.RankOneUpdate, Matrix(K) ≈ Matrix(A) + 2.0 * w * w')

# output
(true, true)

is_separable(a::BilinearForm) -> Bool

Whether a's resolved AST is a sum of terms each expressible as a sum of Kronecker products of one-dimensional factors over a MeshnD: every addend, its constant (literal or Ref) scalar coefficients stripped, has a projection onto the mesh's axes (_kron_project).

a's trial and test space must both be a ScalarGridSpace sharing one mesh, at least two-dimensional (a 1D mesh has nothing to factor), or composite spaces whose leaves all live on one such mesh; every (test leaf, trial leaf) block of a must then be separable in this sense.

What factors:

• the difference, average, jump and shift families along any axis, and chains of them (D₋ₓ(D₋ᵧ(u))), so mixed derivatives and advection terms, whose factors need not be symmetric;

• innerₕ, inner₊ and inner_Γ weights, the last as one term per face;

• a restriction to :interior, provided the mesh's :interior marker is the geometric one (the product of its axes' own interiors);

• a grid-function coefficient that varies along one axis only;

• a scalar or Ref coefficient around a term or inside a side, so the simplifier merging innerₕ(u, v) + a * innerₕ(D₋ₓ(u), v) into innerₕ(u + a * D₋ₓ(u), v) changes nothing, and a Ref is read at every product.

What does not, and answers false:

• a grid-function coefficient varying along several axes, or living on another mesh;

• a region restriction other than :interior, including Dirichlet rows;

• an interpolation, a 1D mesh, trial and test spaces on different meshes, composite leaves on different meshes;

• the star and cross-weighted differences, and any node without a projection.

A false negative only forgoes the Kronecker fast path, so this never claims separability it cannot back up with factors.

Stale spaces throw. If a's spaces were built before an in-place mutation of their mesh, is_separable, kronecker_operator and fdm_solve(a, F) throw the space's stale-weights ArgumentError, exactly as assemble(a) does, rather than answer.

Examples

Ωₕ = mesh(domain(interval(0.0, 1.0) × interval(0.0, 1.0)), (9, 7), (false, false))
Wₕ = gridspace(Ωₕ)
is_separable(form(Wₕ, Wₕ, (u, v) -> innerₕ(u, v) + inner₊(∇ₕ(u), ∇ₕ(v))))  # true

gₕ = Rₕ(Wₕ, x -> 1 + x[1])
is_separable(form(Wₕ, Wₕ, (u, v) -> innerₕ(gₕ * u, v)))  # true: varies along one axis

fₕ = Rₕ(Wₕ, x -> 1 + x[1] * x[2])
is_separable(form(Wₕ, Wₕ, (u, v) -> innerₕ(fₕ * u, v)))  # false: varies along both axes

See also: kronecker_operator, KroneckerLinearOperator.

kronecker_operator(a::BilinearForm) -> Union{KroneckerLinearOperator, KroneckerBlockOperator}
kronecker_operator(a::GhostBilinearForm) -> KroneckerLinearOperator

Build a matrix-free KroneckerLinearOperator for the separable bilinear form a (see is_separable), without ever assembling the D-dimensional matrix.

Each term's factor on axis d is the 1D form it projects to there, assembled on gridspace(Ωₕ(d)) (_kron_project): a diagonal factor is stored as a Diagonal (the axis's cell measures weights(gridspace(Ωₕ(d)), Innerh()) themselves for a mass factor), and terms sharing a factor on an axis share one matrix. Factors are always built on the host (a device mesh through its host mirror) and then converted to the storage backend(trial_space(a)) uses, so a device-backed form yields device-resident factors. A scalar or Ref coefficient around a term stays live, read at every mul!.

The factors also describe the mesh as it is now: if the mesh is mutated in place after this call, the operator throws instead of applying (see KroneckerLinearOperator).

On a composite trial or test space whose leaves all share one mesh, the result is a Bramble.KroneckerBlockOperator (not exported): one KroneckerLinearOperator per nonzero (test leaf, trial leaf) block, placed at that block's rows and columns, with the same size, getindex, mul!, SparseMatrixCSC and issymmetric interface.

On a form over ghost spaces, form(Cₕ, Cₕ, f) with Cₕ built by constrained, the result is a KroneckerLinearOperator of the reduced matrix E_vᵀ A_pad E_u that assemble(a) returns, on the degrees of freedom of Cₕ. Its closure must be separable: every ghost eliminated by one relation per facet along its normal (u => g, ∂ₙ(u), α * u + β * ∂ₙ(u) with constant α and β, or ∂ₙ(u; K) with diagonal K). The closures E_u and E_v are then Kronecker products of 1D closures, and each term's factor on axis d is E_{v,d}ᵀ F_d E_{u,d}, F_d its factor on the padded axis, so dims counts the points of each axis that no u => g fixes. The relation data enter the load only, not the operator.

Throws

• ArgumentError: on ghost spaces, a closure that is not separable (a conormal K with off-diagonal entries, ghosts left as unknowns), a term that reads the padded points beyond two ghosted facets (a mixed derivative), or a region restriction or grid-function coefficient, each naming assemble(a) as the fix.

• ArgumentError: a form over composite ghost spaces.

• ArgumentError: a is not separable, naming the offending node and, on a composite space, its block (or the dimension, or the spaces) – the same check is_separable runs, made specific; or a is device-backed and a term has more than one non-diagonal factor or a non-symmetric one, which the device kernel cannot apply (_kron_check_device).

Examples

Ωₕ = mesh(domain(interval(0.0, 1.0) × interval(0.0, 1.0)), (9, 7), (false, false))
Wₕ = gridspace(Ωₕ)
a = form(Wₕ, Wₕ, (u, v) -> innerₕ(u, v) + inner₊(∇ₕ(u), ∇ₕ(v)))
K = kronecker_operator(a)
x = rand(ndofs(Wₕ))
K * x ≈ assemble(a) * x

See also: is_separable, KroneckerLinearOperator.

KroneckerLinearOperator{T, D, TermsT <: Tuple, P <: ExecutionPolicy, M}

A matrix-free linear operator for a separable BilinearForm (see is_separable): the sum, over its terms, of a Kronecker product of D one-dimensional factor matrices, applied in one fused pass over the grid (LinearAlgebra.mul!(y, K, x), or the five-argument mul!(y, K, x, α, β) computing α * K * x + β * y) rather than ever materialising the D-dimensional matrix: each entry of y is its own term-weighted combination of x at that grid point and its neighbours along the axes where a term's factor is not diagonal, written once. For a 200^3 mesh the factors together hold O(200) numbers per axis instead of the assembled matrix's O(200^3) stored entries.

A term's factors need be neither diagonal nor symmetric, and a term may have non-diagonal factors on several axes. A mixed derivative D₋ₓ(D₋ᵧ(u)) has non-diagonal factors on two axes, and an advection term has a non-symmetric one. A factor that is not symmetric is applied through the CSC of its transpose, so the product is exact for it. Which forms have such factors is is_separable's list.

Build one with kronecker_operator. Subtypes AbstractMatrix{T} so it plugs into LinearProblem/KrylovJL_CG (LinearSolve.jl) the same way an assembled matrix does, and supports size, eltype, getindex, Base.:*, three- and five-argument mul!, LinearAlgebra.issymmetric, and SparseMatrixCSC(K) (an explicit kron of the factors, for testing and inspection: the very matrix this operator avoids forming). issymmetric(K) is computed from the factors, not assumed: it is false as soon as one factor of one term is not symmetric, even for a form whose matrix happens to be symmetric. mul! refuses vectors that are not 1-based (an OffsetArray, say), which would otherwise give a wrong product silently.

K holds no work buffers: it stores its D one-dimensional factors per term, its execution policy, and the mesh and mesh version it checks for staleness, so it is immutable after construction and safe to share across threads (concurrent mul! calls on one K never race). The fused pass needs no scratch either. A serial host mul! allocates nothing, and so does a CpuPolyester product: each task rebuilds the light structs around plain arrays, so Polyester's argument box stays on the stack. mul! is generic over the element type, so ForwardDiff Duals pass through. The scratch = (b1, b2) keyword that mul!(y, K, x; scratch) and mul!(y, K, x, α, β; scratch) accept is kept so existing callers still work, and is ignored.

P is the execution policy of the space the operator was built on (execution_policy(backend(trial_space(a)))), kept as the policy field. On the host, a CpuThreaded operator runs the grid lines along axis 1 in one contiguous chunk per thread under Threads.@threads :static (serially when nested in a threaded loop), a CpuPolyester one as Polyester.@batch tasks (using Polyester required); every other policy runs them serially. Each line writes its own slice of y, so all give the same result bit for bit. An operator built directly from its terms, KroneckerLinearOperator{T, D, TermsT}(terms, dims, n), is CpuSerial.

On a device-backed form (gpena/Bramble.jl#323) the factors are built on the host and then moved to the space backend's device storage, so mul! with device x/y runs entirely on the device, as one KernelAbstractions kernel with one work item per entry of y (using KernelAbstractions required). The kernel reads at most one non-diagonal factor per term and that factor's columns as its rows, so kronecker_operator refuses, on such a form, a term with non-diagonal factors on several axes and a term with a non-symmetric factor: build those on a host backend. getindex and SparseMatrixCSC(K) stay host-only and throw an ArgumentError on such an operator.

The factors are read from the mesh once, when the operator is built (gpena/Bramble.jl#442): K keeps that mesh (M is its type) and its version then. After an in-place mutation of the mesh (set_points!, change_points!, iterative_refinement!), mul!, getindex, SparseMatrixCSC(K), fdm_solve(K, F) and Kronecker.kronecker(K) throw an ArgumentError, as a space's stale weights do; build the operator again on the mutated mesh. An operator built directly from its terms has no mesh (M is Nothing) and never goes stale.

Dirichlet rows are out of scope: this operator carries no boundary constraint of its own. The Kronecker.jl extension's fast-diagonalisation solve, fdm_solve, imposes homogeneous Dirichlet conditions on the whole boundary through its dirichlet = :boundary keyword.

See also: is_separable, kronecker_operator.

fdm_solve(a::BilinearForm, F::AbstractVector; dirichlet = nothing) -> Vector
fdm_solve(K::KroneckerLinearOperator, F::AbstractVector) -> Vector

Directly solve assemble(a) \ F (or the linear system K represents) for a Laplacian-like BilinearForm, without assembling a's matrix. A form is Laplacian-like when it is separable and every term differs from one mass per axis on at most one axis. Symmetric 1D operators are solved by fast diagonalisation; non-symmetric ones (advection terms) by a complex generalised Schur factorisation per axis and a triangular back substitution. A mixed-derivative form is refused: precondition a Krylov solver with fdm_preconditioner instead. fdm_factorize and fdm_solve! split the factorisation from the solve.

dirichlet, when given, must request homogeneous Dirichlet conditions on the whole mesh boundary; K alone carries no boundary handling, since a KroneckerLinearOperator has no tensor structure of its own to restrict.

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

See also: kronecker_operator, KroneckerLinearOperator, is_separable.

fdm_factorize(a::BilinearForm; dirichlet = nothing)
fdm_factorize(K::KroneckerLinearOperator; dirichlet = nothing)

Factorise, once, the system fdm_solve solves, so that each further right-hand side costs only the solve. The result supports fdm_solve!, ldiv!(x, f, F) (the same solve) and size; its type is internal. fdm_solve(a, F) is fdm_factorize(a) followed by one solve.

a must be a form fdm_solve accepts, and fdm_factorize(a; dirichlet) refuses what fdm_solve(a, F; dirichlet) refuses. K must be Laplacian-like, as for fdm_solve(K, F); unlike that method, fdm_factorize(K) also takes dirichlet = :boundary. With dirichlet = :boundary, each right-hand side must be zero on the boundary and the solution is zero there. The factorisation is built from the mesh and coefficients as they are when fdm_factorize is called: after change_points! or a change of a coefficient (a Ref one included), refill it with fdm_factorize!, or call fdm_factorize again. Construction costs O(n_d^3) per axis; a solve on host vectors costs O(N Σ_d n_d) and allocates nothing. The factorisation is host-only.

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

Throws

• ArgumentError saying fdm_solve does not support the form or operator, with the reason, for one that is not Laplacian-like or is singular.

• ArgumentError: dirichlet is neither nothing nor :boundary.

• ArgumentError naming change_points!: K's mesh was mutated in place after K was built; build the operator again.

Examples

using Bramble, Kronecker
Ωₕ = mesh(domain(interval(0.0, 1.0) × interval(0.0, 1.0)), (25, 19), (false, false))
Wₕ = gridspace(Ωₕ)
a = form(Wₕ, Wₕ, (u, v) -> innerₕ(u, v) + inner₊(∇ₕ(u), ∇ₕ(v)))
f = fdm_factorize(a)
x = zeros(ndofs(Wₕ))
for k in 1:3
    F = rand(ndofs(Wₕ))
    fdm_solve!(x, f, F)   # x ≈ assemble(a) \ F
end

See also: fdm_solve!, fdm_solve, fdm_preconditioner.

fdm_factorize!(f, a::BilinearForm) -> f
fdm_factorize!(f, K::KroneckerLinearOperator) -> f

Refill the factorisation f, which fdm_factorize returned, from a or K, and return f. Use it in place of a new fdm_factorize when the structure is the same and only the numbers changed, as when a coefficient that is a scalar or a Ref took a new value, or change_points! moved the grid points and a is a form on a new gridspace (or K a new kronecker_operator) built after it. f keeps its dirichlet, so there is no keyword here. Afterwards fdm_solve!(x, f, F) solves the new system, and agrees bitwise with the solve of a fresh fdm_factorize, unless a term's coefficient is now zero: that term keeps its place, and the two agree to rounding. Only Float32 and Float64 are supported.

With K, a refill allocates nothing once warm. With a, it allocates nothing once warm for the form f was built (or last refilled) from, when only scalar or Ref coefficients changed. A different form is projected again, and so is one that reads a grid-function coefficient, so that an edit made with Rₕ! is read and never ignored; both allocate. The operator or form must have the structure f was built for: the same dimension, sizes, eltype, number of terms, mass term on each axis, terms left out at build (a zero coefficient, or a factor that is zero on the solved points), and symmetric or non-symmetric (Schur) solve. Anything else needs a new fdm_factorize. A KroneckerLinearOperator backed by a device is refused: the refill is host-only. An FDMPreconditioner is not refilled.

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

Throws

• ArgumentError starting fdm_factorize! cannot refill this factorisation, then the reason: a or K differs from f in dimension, size, eltype, number of terms, mass term, a term left out at build that is present now, or whether the solve is symmetric or Schur, or K is backed by a device. f is left as it was and still solves.

• ArgumentError: the refilled system is one fdm_solve refuses, singular or with a mass that is not symmetric positive definite. f then cannot solve, and fdm_solve! says so, until a later refill succeeds.

• ArgumentError naming change_points!: the mesh of a or K was mutated in place after it was built; build the form on a new gridspace (or the operator again).

• ArgumentError saying fdm_solve does not support the form, with the reason: a is on a composite space, or is not separable.

Examples

using Bramble, Kronecker
Ωₕ = mesh(domain(interval(0.0, 1.0) × interval(0.0, 1.0)), (25, 19), (false, false))
Wₕ = gridspace(Ωₕ)
c = Ref(1.0)
a = form(Wₕ, Wₕ, (u, v) -> innerₕ(u, v) + c * inner₊(∇ₕ(u), ∇ₕ(v)))
f = fdm_factorize(a)
F = rand(ndofs(Wₕ))
x = zeros(ndofs(Wₕ))
for k in 1:3
    c[] = 10.0^k
    fdm_factorize!(f, a)   # no allocation once warm: only a Ref changed
    fdm_solve!(x, f, F)    # x ≈ assemble(a) \ F
end

See also: fdm_factorize, fdm_solve!, fdm_solve.

fdm_solve!(x::AbstractVector, f, F::AbstractVector) -> x

Solve into x, with the factorisation f that fdm_factorize returned, the system it was built for, right-hand side F. With dirichlet = :boundary, x is zero on the boundary. x may alias F. On host vectors it allocates nothing.

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

Throws

• DimensionMismatch: x or F does not have size(f, 1) entries.

• ArgumentError: x or F is not 1-based.

• ArgumentError: F holds a complex value; solve its real and imaginary parts separately.

• ArgumentError naming fdm_factorize!: the last fdm_factorize! of f was refused because the system is singular or its mass is not symmetric positive definite. Refill f before solving again.

See also: fdm_factorize, fdm_factorize!, fdm_solve.

matrix_free_operator(a::BilinearForm; dirichlet = nothing, dirichlet_components = nothing, policy = execution_policy(trial_space(a))) -> MatrixFreeOperator
matrix_free_operator(a::GhostBilinearForm; dirichlet = nothing, dirichlet_components = nothing, policy) -> GhostMatrixFreeOperator

The linear operator of a, applied without assembling its matrix: op * x agrees with assemble(a; dirichlet, dirichlet_components) * x.

A form on ghost spaces (GhostGridSpace) gives an operator that applies its reduced matrix E_vᵀ A_pad E_u: the closure of the trial space without its affine part, the operator of the form on the padded spaces, then the transpose of the test space's closure. The affine part belongs to the load, as in assemble(a, l). Its padded buffers are built once, so a warm product allocates no vector. It applies one product at a time: a product entered on it while another runs throws an ArgumentError, so build one operator per task. dirichlet there takes what assemble(a; dirichlet) takes.

Arguments

• a: The bilinear form. Any form assemble accepts, on scalar or composite spaces or on scalar ghost spaces, except one holding a rank-one term.

Keywords

• dirichlet: The Dirichlet labels whose rows become identity rows, in any form assemble accepts (a symbol, a tuple of labels, a pair with its value, dirichlet_constraints); boundary values are ignored, as they are by the matrix (default: nothing).

• dirichlet_components: The leaves of a composite test space the labels bind to, as in dirichlet_bc!. It is an Int, a Tuple of Ints, or nothing for every leaf (default: nothing).

• policy: The ExecutionPolicy of mul! (default: the trial space's). CpuSerial walks the form on the calling task; any other CPU policy threads it, with each leaf's own mechanism (see MatrixFreeOperator).

Returns

• MatrixFreeOperator: of size (ndofs(test_space(a)), ndofs(trial_space(a))) and the element type of assemble(a). On ghost spaces, a GhostMatrixFreeOperator of that size holding the padded form's MatrixFreeOperator.

Throws

• ArgumentError: policy is a GpuPolicy; device execution is tracked on milestone v4.7.0.

• ArgumentError: dirichlet_components names a leaf the test space does not have.

• ArgumentError: dirichlet on ghost spaces whose unknowns are not the physical nodes, as for assemble(a; dirichlet).

• ArgumentError: a holds a rank-one term.

• ArgumentError: a is a form over composite ghost spaces, which have no matrix-free operator.

• ArgumentError, from mul!: a product on a ghost-space operator while another task's product on it runs.

Examples

using Bramble, LinearAlgebra
Wₕ = gridspace(mesh(domain(interval(0.0, 1.0)), 11, false))
a = form(Wₕ, Wₕ, (u, v) -> innerₕ(u, v) + inner₊(∇ₕ(u), ∇ₕ(v)))
op = matrix_free_operator(a; dirichlet = :boundary)
x = rand(ndofs(Wₕ))
y = similar(x)
mul!(y, op, x)
y ≈ assemble(a; dirichlet = :boundary) * x

# output
true

See also: MatrixFreeOperator, assemble, kronecker_operator.

MatrixFreeOperator{T, Form, DLabels, EP, Plan, Geom} <: AbstractMatrix{T}

A BilinearForm as a linear operator: mul!(y, op, x) computes A * x, and mul!(y, op, x, α, β) computes α * A * x + β * y, for A = assemble(a; dirichlet), without ever storing A. Each product walks the form's stencil once, evaluated afresh, so a coefficient updated in place (Rₕ!, a Ref) is seen by the next product. Build one with matrix_free_operator.

x and y are plain vectors or VectorElements: op * x returns a new vector, op * uₕ a new element of the form's test space. mul! allocates nothing on a serial policy, except on a form with a region restriction (restrict_to), whose stencil evaluation allocates as it does in assemble!.

On a CpuThreaded or CpuPolyester policy the grid is cut into one band per thread along its last axis, and each thread walks every term over its band and adds only into the rows of its band's points, so no two threads add into the same entry of y and the result does not depend on scheduling. A product is one parallel region however many terms the form has; the terms' reach, which widens each band by the points writing into it, is computed when the operator is built. Each row receives its entries in the serial order, so the threaded product is bitwise the serial one (up to rounding when a term runs serially, see below). The threads are those of the leaves' own mechanism: Polyester.@batch on a CpuPolyester leaf, one Threads.@spawned task per band otherwise. A term with a test-side interpolation runs serially. A form whose leaves walk grids of different sizes, or carry different policies, is swept one term at a time instead, each in the colour bands the threaded assemble! uses. A threaded product allocates its task launches, a fixed cost that does not grow with the grid.

Dirichlet rows follow dirichlet_bc!: a row in ΓD is the identity row, so (A * x)[i] = x[i] there and its columns are untouched. A row in ΓD past the last column (a rectangular form) is zero, as in the assembled matrix.

Subtypes AbstractMatrix{T} so it can stand in for the assembled matrix in a Krylov solve (LinearProblem with KrylovJL_CG, say). getindex is supported for inspection and costs one product per entry.

Type parameters

• T: The element type of assemble(a).

• Form: The form's type.

• DLabels: Nothing, or the Tuple of Dirichlet labels.

• EP: The ExecutionPolicy the operator was built with.

• Plan: Nothing, or the type of the threaded sweep's plan (the grid its bands cut and the terms' reach), fixed when the operator is built.

• Geom: Nothing, or the type of the per-axis mesh geometry cached when the operator is built (_MFGeometry, a scalar form's only). It holds no coefficient, and a product on a mesh moved since reads the mesh afresh, as it would without it.

See also: matrix_free_operator, assemble, KroneckerLinearOperator.

assemble_add!(F::AbstractVector, l::LinearForm) -> AbstractVector
assemble_add!(F::AbstractVector, l::LinearForm, α) -> AbstractVector

Add l's contribution to the preallocated F, in place, without first zeroing it – the linear-form counterpart of assemble_add!(A, a). A linear form's own assembly scatters directly into a dense F (no sparsity pattern, no record/replay cache), so this is simply the existing scatter with the initial fill! and the Dirichlet pass left out, same as the bilinear case; see this file's header for the scale factor and the Dirichlet note, both of which apply identically here.

On a device test space (a Metal space, say), F is copied to the host, the contribution added there, and the sum uploaded in one copyto!, so that path allocates. An F whose length is not ndofs(test_space(l)) raises a DimensionMismatch before anything is written.

Examples

F = zeros(ndofs(Wₕ))
assemble_add!(F, l_mass)
assemble_add!(F, l_source, θ)   # F += θ[] * l_source's contribution

See also: assemble!, assemble.

jacobian_pattern(a::BilinearForm, coefficient_dependencies::Function...) -> SparseMatrixCSC{Bool}

Sparsity pattern of the Jacobian of a Newton residual A(u) * u - F, where a is the BilinearForm that assembles A(u) and each function in coefficient_dependencies names, symbolically, the stencil operator one of a's live coefficients was itself computed from – written the same way a form term names an operator, as a function of the trial placeholder. If a's diffusion coefficient was built as αvals = α.(Mₕ(uₕ)), pass U -> Mₕ(U).

Derived entirely from a's AST and each dependency's own stencil reach – no AD tracing, no coefficient values needed. A safe superset of the exact pattern (a pointwise nonlinear α can never narrow the reach its argument already has), suitable for ADTypes.KnownJacobianSparsityDetector:

using Bramble: jacobian_pattern
a = form(Wₕ, Wₕ, (U, V) -> inner₊(αvals * ∇ₕ(U), ∇ₕ(V)))   # αvals = α.(Mₕ(uₕ))
pattern = jacobian_pattern(a, U -> Mₕ(U))
sparse_ad = AutoSparse(AutoForwardDiff();
    sparsity_detector = KnownJacobianSparsityDetector(pattern),
    coloring_algorithm = GreedyColoringAlgorithm())

Composite trial/test spaces

A dependency may also name a different leaf, the same way a form term does – U -> U(2) for a coefficient that is component 2's own value (no stencil op), or U -> Mₕ(U(2)) for one built from a stencil op applied to that other component. nothing named (U -> Mₕ(U), no (k)) means the coefficient depends on this block's own trial leaf, exactly like the non-composite case above. Every dependency still applies to every block the walk visits, whichever leaf it names – a safe superset stays safe however many blocks end up seeing an entry they did not strictly need.

using Bramble: jacobian_pattern
# a = form(Vₕ, Vₕ, (p, q) -> inner₊(∇ₕ(p(1)), ∇ₕ(q(1))) + innerₕ(v_c * p(1), q(1)) +
#                            inner₊(∇ₕ(p(2)), ∇ₕ(q(2))) - innerₕ(u_c * p(2), q(2)))
pattern = jacobian_pattern(a, U -> U(2), U -> U(1))   # block (1,1) reads U(2), (2,2) reads U(1)

Dependencies through πₕ

A dependency like the ones above is read at offsets from the point being walked, so the leaf it names must be discretised on the mesh the term is walked on. When the coefficient lives on that mesh but was computed from a trial function on another one, the dependency says how, through πₕ: U -> πₕ(U) for cvals = α.(πₕ(Wₐ, uₕ)), U -> Mₕ(πₕ(U)) for cvals = α.(Mₕ(πₕ(Wₐ, uₕ))), and U -> πₕ(U(2)) for a coefficient interpolated from component 2 of a composite trial space. At each walked point the pattern then adds, to every row the term reaches there, the source-mesh columns that composition reads at every point the coefficient is read: the reach of the operators enclosing it, so {0, -1} in D₋ₓ(cvals * πₕ(u)) and {-1, 1} in Dcₓ(cvals * Dcₓ(πₕ(u))). Those columns are the corners of the cell that point falls in, and for an operator applied outside πₕ the corners at each neighbour that operator reaches. A dependency naming a leaf on another mesh without πₕ (U -> Mₕ(U)) is refused with an ArgumentError. So is one naming a component the trial space does not have (U -> U(3) on W × W).

using Bramble: jacobian_pattern
# uₕ ∈ Wᵦ, cvals = α.(Mₕ(πₕ(Wₐ, uₕ))) ∈ Wₐ
a = form(Wᵦ, Wₐ, (u, v) -> inner₊(cvals * ∇ₕ(πₕ(u)), ∇ₕ(v)))
pattern = jacobian_pattern(a, U -> Mₕ(πₕ(U)))

ast_sparsity_detector(a::BilinearForm, coefficient_dependencies::Function...)

An ADTypes.AbstractSparsityDetector that supplies jacobian_pattern(a, coefficient_dependencies...) directly as a Newton residual's Jacobian sparsity, in place of one detected by tracing:

using Bramble: ast_sparsity_detector
sparse_ad = AutoSparse(AutoForwardDiff();
    sparsity_detector = ast_sparsity_detector(a, U -> Mₕ(U)),
    coloring_algorithm = GreedyColoringAlgorithm())

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

jacobian_form(a::BilinearForm, coefficients::Pair...; dirichlet = nothing)

The exact Jacobian of the Newton residual R(u) = A(u) u - F as one form, a + ∂a, where a assembles A(u) and reads coefficients that depend on the state. Each is declared as values => (U -> L(U), f): values is the grid function a reads (the same object, by identity), L the stencil operator it is computed from, written as a function of the trial placeholder as jacobian_pattern takes it, and f a scalar function, so that values == f.(L(u)).

using Bramble: jacobian_form, jacobian_prototype, jacobian!
using ForwardDiff                      # differentiates each `f`
αx, αy = element(Wₕ), element(Wₕ)
a = form(Wₕ, Wₕ, (U, V) -> inner₊((αx * ∇ₕ(U)[1], αy * ∇ₕ(U)[2]), ∇ₕ(V)))
jf = jacobian_form(a, αx => (U -> Mₕ(U)[1], α), αy => (U -> Mₕ(U)[2], α);
    dirichlet = :boundary)
J = jacobian_prototype(jf)             # the pattern of a + ∂a, wider than A's in general
jacobian!(J, jf, u)                    # J = ∂R/∂u at u

jacobian_prototype(jf) allocates a matrix with the pattern of a + ∂a, which can be wider than A's. jacobian!(J, jf, u) writes the Jacobian at u into it: it sets every declared coefficient to f.(L(u)) and each derivative buffer to f′.(L(u)) ⊙ T(u), where T is the operator the coefficient multiplies, then assembles a + ∂a once, with the Dirichlet labels dirichlet (as assemble! takes them). The coefficients are written from u itself, so the Jacobian is exact at u whatever state they held before. Each f′ is one ForwardDiff.Dual evaluation of f, so ForwardDiff.jl must be loaded. The refill runs under the mesh's execution policy.

The derivative is built by a walk over a's AST: the product rule where a declared coefficient scales an operator, and every node above it (scales, sums, differences, averages, jumps, shifts, region restrictions, the trial side of an inner product) rebuilt around the derivative of its operand. Every inner product's weight is pointwise, so a boundary term such as a Robin term inner_Γ(c * U, V; markers) needs nothing more, and a coefficient scaling the whole test side, inner₊(∇ₕ(U)[1], c * ∇ₕ(V)[1]), is moved to the trial side. A declared coefficient under a node with no such rule, under an operator on the test side, or in any field a rule does not read, throws MissingJacobianRule naming the node. A declared coefficient the form never reads throws an ArgumentError, and so does a dependency L that reads a declared coefficient (itself included): L may read only the state and constant grid functions. Composite spaces are routed per block.

The trial and test spaces are both scalar grid spaces or both composite ones. On a composite space each declared coefficient scales one component, c * U(i), and its dependency names the component it reads, U -> L(U(j)); the derivative term of c * U(i) in the equation of component k then lands in block (k, j), so a coefficient reading another component couples the two. A coefficient scaling a trial function with no component, or a dependency naming none, throws an ArgumentError, and so does a dependency naming a component on another mesh.

MissingJacobianRule(node::Symbol, reason::String)

Thrown by jacobian_form when a node on a path from the form's root to a declared coefficient has no derivative rule. node names the node's type. The term is refused rather than dropped, so the Jacobian is never silently missing a part.

DirichletConstraint = DomainMarkers

Type alias for Dirichlet boundary constraint storage.

dirichlet_constraints(input, [I::CartesianProduct{1}], pairs::Pair...) -> DomainMarkers

Create Dirichlet boundary constraints.

Each pair is of the form :label => func, where :label identifies the boundary region and func defines the Dirichlet values. If the optional time domain I is provided, func must be a time-dependent function func(x, t), or a time- and parameter-dependent one func(x, t, p) for a value that a semidiscretize residual's own p should reach; this is checked by arity, since nothing about func itself can be evaluated at construction time. The time-dependent pairs given alongside I must agree on which of the two shapes they use, and at least one pair must be time-dependent. A pair whose func takes x alone may stay beside them: it is kept as it is at every t, as semidiscretize and the dirichlet keyword do.

I is used only to require the time-dependent shape of func: it is not stored, and nothing compares it with the time span of an ode_problem. The conditions are evaluated at whatever t the solver passes, so I need not match the time span, and a time span that reaches outside I evaluates func outside I without a warning. Pass the same interval to both so that func is only called on the time span it was written for.

input can be a CartesianProduct, a Domain, an AbstractMeshType, a ScalarGridSpace, or a CompositeGridSpace from which the mesh is extracted. The :label must match a label in the mesh definition.

dirichlet_constraints(X::CartesianProduct, f::Function) -> DomainMarkers

Create a single Dirichlet boundary constraint with function f under the :boundary label.

dirichlet_bc!(A::AbstractMatrix, Ωₕ::AbstractMeshType, labels::Symbol...) -> AbstractMatrix

Apply Dirichlet boundary conditions to matrix A based on marked regions in the mesh Ωₕ.

For each index i associated with the given Dirichlet labels, this function:

1. Sets all elements in the i-th row of A to zero.

2. Sets the diagonal element A[i, i] to one, if A has an i-th column.

dirichlet_bc!(A::AbstractMatrix, space::CompositeGridSpace, labels::Symbol...; components = nothing) -> AbstractMatrix

Apply Dirichlet boundary conditions to matrix A on the regions named by labels, restricted to the leaf components specified in components (1-based positions in leaf_spaces_offsets(space), following the depth-first ordering used by u(1)/u(2) addressing). components = nothing (the default) applies to every leaf component.

This allows coupled systems to constrain selected fields while leaving others unconstrained (for example, prescribing velocity while leaving pressure free in a Stokes problem):

Wₕ = vector_gridspace(Ωₕ, Val(2))   # 1: velocity, 2: pressure
dirichlet_bc!(A, Wₕ, :left, :right; components = 1)   # velocity only

Successive calls with different labels/components pairs compose cleanly.

dirichlet_bc!(v::AbstractVector, Ωₕ::AbstractMeshType, bcs::ConstraintMarkers, labels::Symbol...) -> AbstractVector

Write Dirichlet values into v at the points marked by labels.

bcs may be unevaluated constraints or time-evaluated constraints (see ConstraintMarkers). Only the marked entries are modified, with complexity proportional to the boundary cardinality.

Throws

• ArgumentError: if a label in labels has no condition in bcs.

dirichlet_bc!(v::AbstractVector, Ωₕ::AbstractMeshType, bcs::EvaluatedParametricDomainMarkers, labels::Symbol...) -> AbstractVector

Write Dirichlet values into v at the points marked by labels, from a time- and parameter-evaluated constraint set (see EvaluatedParametricDomainMarkers).

A separate method rather than a third member of ConstraintMarkers: the body is otherwise identical to the bcs::ConstraintMarkers method above, but widening that Union is exactly the shape of change gpena/Bramble.jl#240 removed from the assembly path for Enzyme's sake, and this method sits on the same residual Enzyme differentiates for a p-dependent boundary value.

dirichlet_bc!(A::AbstractMatrix, Vₕ::CompositeGridSpace, labels::Symbol...; components = nothing)
dirichlet_bc!(v::AbstractVector, Vₕ::CompositeGridSpace, bcs, labels::Symbol...; components = nothing)
symmetrize!(A, F, Vₕ::CompositeGridSpace, labels::Symbol...; components = nothing)

Row replacement by hand on a composite of ghost spaces, as the dirichlet keyword of assemble does it: the rows of the nodes of labels (markers or facet names of the physical mesh) in each component components selects (all by default). Each selected component's unknowns must be its physical nodes: every ghost eliminated by a relation and no node fixed by u => g; an ArgumentError names the route otherwise.

dirichlet_bc!(u::AbstractVector, sd::Semidiscretization, t::Number, p = nothing) -> u

Write sd's Dirichlet values at time t (and, for a parameter-dependent condition, parameter p) into u and return it.

An index-1 differential-algebraic system needs its initial condition to satisfy the algebraic rows already: a u disagreeing with g(x, 0) on the boundary is inconsistent, and a stiff solver either rejects it or absorbs it into the first step. ode_problem applies this to a copy of the initial condition it is handed, forwarding whatever p the caller gave the ODEProblem.

symmetrize!(A, F, Wₕ::ScalarGridSpace, labels...; components = nothing): the same on mesh(Wₕ); components may only be nothing or 1

symmetrize!(A::AbstractMatrix, F::AbstractVector, Ωₕ::AbstractMeshType, labels::Symbol...)

Modify the linear system Ax = F to restore symmetry in A after applying Dirichlet conditions. Updates F to account for prescribed boundary data and zeroes the corresponding columns in A.

For each index i with prescribed Dirichlet boundary conditions:

• Calculate dᵢ = cᵢ .* F[i], where cᵢ is the i-th column of A;

• Update F by subtracting dᵢ from F (preserving the i-th component);

• Zero out off-diagonal elements in the i-th column of A.

A non-square A (a form whose test and trial spaces differ in size) has no symmetry to restore and throws an ArgumentError.

symmetrize!(A::AbstractMatrix, F::AbstractVector, Wₕ::CompositeGridSpace, labels::Symbol...; components = nothing)

Symmetrize a coupled linear system leaf space by leaf space.

The counterpart of the composite dirichlet_bc!: each leaf's marker mask is read at that leaf's offset into the global system without allocating full-system masks. leaf_spaces_offsets returns a tuple, enabling loop unrolling and type stability.

Takes the same components keyword as composite dirichlet_bc!, restricting which leaf components labels binds to (1-based positions in leaf_spaces_offsets(Wₕ)). components = nothing (the default) applies to every leaf.

Combines every label into one mask per leaf (_combined_mask, space/innerproduct.jl) rather than visiting each leaf once per label: by the time symmetrize! runs, `dirichletbc!has already zeroed every marked row ofAoff its diagonal (the documented order of use), so for any two marked indicesi,kin the combined set,A[k, i] == 0`. That makes each marked index's elimination read/write entries that no other marked index's elimination touches, so the order they are visited in cannot change the result, and a row marked by several labels is visited once.

symmetrize!(K::RankOneUpdate, F::AbstractVector, Wₕ, labels::Symbol...; components = nothing)

symmetrize! for a RankOneUpdate K = A + σ wₛ wₜᵀ, after the Dirichlet rows are applied (assemble(a; dirichlet)). A is symmetrized as a matrix is, and the rank-one columns of the Dirichlet indices D move into the load: F -= σ wₛ Σⱼ∈D wₜ[j] F[j], then wₜ[D] = 0. As wₛ is zero at D, the Dirichlet rows of F keep their values.

reaction(A::AbstractMatrix, F::AbstractVector, uₕ::VectorElement; marker, components = nothing) -> Real

Net flux through the region(s) named by marker, from the unconstrained operator A and load F (assemble(a)/assemble(l) with no dirichlet keyword – dirichlet_bc! would have already overwritten the rows this needs).

marker is a Symbol or a Tuple of them, exactly as dirichlet accepts one or several labels; overlap between markers is counted once (the same union dirichlet_bc! takes). components restricts a composite space's leaves (1-based positions in leaf_spaces_offsets, nothing meaning every leaf).

Sign and scaling

Positive means flux leaving the domain along the outward normal: for -Δu = f, ∮_Γ reaction = ∫_Ω f to round-off. The residual already carries the cell-measure weighting the form's inner product imposes, so no further rescaling is applied here – see reaction_density for the pointwise quantity that does.

Examples

Ωₕ = mesh(domain(interval(0.0, 1.0), :left => :xmin, :right => :xmax), 41)
Wₕ = gridspace(Ωₕ)
sol(x) = sin(pi * x[1])
src(x) = pi^2 * sin(pi * x[1])

a = form(Wₕ, Wₕ, (u, v) -> inner₊(∇ₕ(u), ∇ₕ(v)))
l = form(Wₕ, v -> innerₕ(Rₕ(Wₕ, src), v))
A, F = assemble(a, l; dirichlet = :boundary => sol)

uₕ = element(Wₕ)
uₕ .= A \ F

reaction(a, l, uₕ; marker = :left)    # ≈ π
reaction(a, l, uₕ; marker = :right)   # ≈ π

See also: reaction_density, assemble, normal_vector.

reaction(a::BilinearForm, l::LinearForm, uₕ::VectorElement; marker, dirichlet_components = nothing) -> Real

Reassembles the unconstrained a/l (assemble with no dirichlet keyword) and calls reaction(A, F, uₕ; marker, components = dirichlet_components). a and l must be posed on the same space as uₕ.

reaction(a::GhostBilinearForm, l::GhostLinearForm, uₕ::VectorElement; marker, dirichlet_components = nothing)

Refuses with an ArgumentError. A reaction reads the rows of the unconstrained plain system at the marked nodes, and a ghost space has replaced those rows by its relations. Pose a and l on the plain space gridspace(Ωₕ), fill an element of it with the physical values of Bramble.lift(C, u), and call reaction there.

reaction!(scratch::AbstractVector, A::AbstractMatrix, F::AbstractVector, uₕ::VectorElement; marker, components = nothing) -> Real

In-place counterpart of reaction: the same net flux, computed with zero allocation using a caller-owned scratch (sized ndofs(space(uₕ))) instead of an internal one – meant for a loop that calls this every step (e.g. a transient solve) and reuses scratch across calls. scratch's ownership/lifetime is the caller's; only the entries marker/components select are written, so reusing scratch across calls with a different marker/components reads stale values outside the new selection – calling with the same selection every time (the intended use) is always safe.

reaction_density(A::AbstractMatrix, F::AbstractVector, uₕ::VectorElement; marker, components = nothing) -> VectorElement

The pointwise counterpart of reaction: a grid function, zero away from marker, equal at each marked point to the boundary flux density q·n there, which converges under refinement. It is reaction's summand divided by the point's surface weight over every boundary face, the weight inner_Γ integrates with, so a marked corner or edge gets the value marker = :boundary gives there, whatever the marker. When every marked corner and edge lies only on marked faces (:boundary, or a tuple of all the sides), inner_Γ(element(space, 1.0), density, faces...) recovers reaction's scalar; a one-face marker's corners carry only part of its flux, so it does not. Meant to be exported and plotted, e.g. export_vtk(filename, Ωₕ, "reaction" => reaction_density(...)).

The faces are read off the marked points, not the label's name: a face counts when every point of it is marked, except possibly its edges and corners. A marked point on no such face (an interior marker, or part of a face) throws ArgumentError; reaction still gives its net flux.

reaction_density(a::BilinearForm, l::LinearForm, uₕ::VectorElement; marker, dirichlet_components = nothing) -> VectorElement

Reassembles the unconstrained a/l and calls reaction_density(A, F, uₕ; marker, components = dirichlet_components).

reaction_density!(dens::AbstractVector, A::AbstractMatrix, F::AbstractVector, uₕ::VectorElement; marker, components = nothing) -> AbstractVector

In-place counterpart of reaction_density: writes the flux density into a caller-owned dens (sized ndofs(space(uₕ))) with zero allocation instead of allocating a fresh one – meant for a loop that calls this every step and reuses dens across calls. As with reaction!, only the entries marker/components select are written, so dens stays exactly zero away from the marker only if it started there (a fresh zeros(...), or a previous call with the same marker/components). Returns dens itself, not a VectorElement; wrap with element(space(uₕ), dens) for a grid function.