Coupled systems#

What you will learn. How to write a form over several unknowns, constrain one block and not another, and read a value across leaves that live on different meshes.

What you need first. The form tutorial for assembling and for Dirichlet conditions, and the space tutorial for composite spaces and vector elements.

Where next. The worked examples, which run the whole chain on full problems.

A composite space stacks grid spaces, and a form over one addresses its blocks by component.

Addressing blocks by component#

u[1] (or u(1)) is the trial function of the first block, and v[2] (or v(2)) the test function of the second. Both spellings work on trial and test functions and on compound operators such as (D₋ₓ(u))[i]. components destructures a function into its blocks:

using Bramble
import Bramble: CompositeGridSpace

Ωₕ = mesh(domain(interval(0.0, 1.0)), 33, true)
Wₕ = gridspace(Ωₕ)
Vₕ = Wₕ^Val(2)
ac = form(Vₕ, Vₕ, (u, v) -> begin
    u₁, u₂ = components(u)
    v₁, v₂ = components(v)
    innerₕ(u₁, v₁) + inner₊(∇ₕ(u₂), ∇ₕ(v₂), :x)
end)
Ac = assemble(ac)
size(Ac)

(66, 66)

Sixty-six by sixty-six: two blocks of 33 in one matrix. A term naming u[i] and v[j] lands in block \((j, i)\), so off-diagonal coupling is written the same way. innerₕ(u[1], v[2]) fills the block that couples the first unknown to the second equation.

Indices are checked against the number of blocks when the form is built: u[3] on a two-component space raises an ArgumentError. A term must name both components or neither, and a term that names one and not the other raises when the matrix is assembled:

try
    assemble(form(Vₕ, Vₕ, (u, v) -> innerₕ(u[1], v)))
catch e
    showerror(stdout, e)
end

ArgumentError: a term of this form names its trial component but not its test one. A block of a bilinear form takes a component from each side: write innerₕ(u(i), v(j)) for one block, or innerₕ(u, v) for the same integrand on every diagonal block.

A term naming one component belongs to every equation at once, so it is refused. Naming neither means the diagonal, applied to every block.

Constraining one block, leaving another free#

dirichlet on its own binds to every leaf that shares the named marker. That is right when every block wants the same treatment. A Stokes-style system prescribes velocity and leaves pressure free, so it also needs dirichlet_components: 1-based leaf positions, in the order u(1), u(2) already use.

Ωc = domain(interval(0.0, 1.0), :left => :left, :right => :right)
Ωdc = mesh(Ωc, 21, true)
Vc = gridspace(Ωdc)^Val(2)           # 1: velocity-like, 2: pressure-like
ac2 = form(Vc, Vc, (u, v) -> innerₕ(u(1), v(1)) + innerₕ(u(2), v(2)))
Ac2 = assemble(ac2; dirichlet = (:left, :right), dirichlet_components = 1)

Block 1 (rows 1:21) has its boundary rows pinned. Block 2 is the plain assembled operator, with no row replaced. Leaving dirichlet_components at its default, nothing, applies the labels to every leaf. Call assemble! again with another dirichlet/dirichlet_components pair, or dirichlet_bc! with another labels/components pair, to constrain a different block.

Interpolating between the leaves of a heterogeneous space#

The composite spaces above stack copies of one space, so every leaf shares a mesh. A composite space can also be built from a tuple of leaves over different meshes. A term that couples two such leaves must move a value from one grid to the other. That is the job of πₕ. Applied to a grid function of one leaf, as below, it is a source. Applied to the trial or the test function, it interpolates the unknown and lets a bilinear term couple leaves on different meshes.

πₕ(uₕ) is a source read live: every assemble takes the current values of uₕ. It composes with D₋ₓ, Mₓ and the rest, and sits on the left of innerₕ inside a coupled form:

Ωbig = mesh(domain(box((0.0, 0.0), (1.0, 1.0))), (8, 8), (true, true))
Ωsmall = mesh(domain(box((0.0, 0.0), (1.0, 1.0))), (4, 4), (true, true))
Wbig, Wsmall = gridspace(Ωbig), gridspace(Ωsmall)
Vh = CompositeGridSpace((Wbig, Wsmall))
uv = Rₕ(Vh, (x -> 0.0, x -> x[1] + x[2]))   # only the small leaf (2) carries data

lh = form(Vh, v -> innerₕ(πₕ(uv(2)), v(1)) + innerₕ(∇ₕ[:x](πₕ(uv(2))), ∇ₕ[:x](v(1))))
bh = assemble(lh)

# the differenced term is not a no-op: dropping it changes the answer
b_plain = assemble(form(Vh, v -> innerₕ(πₕ(uv(2)), v(1))))
maximum(abs, bh .- b_plain)

0.14285714285714285

Both terms land in leaf 1, the big mesh, although the source lives on leaf 2's coarser mesh. Leaf 2 can hold one field at the resolution the problem calls for, and a term over leaf 1 can still read it.

The last line compares values: a differenced source whose offsets were discarded assembles to exactly zero, which has the right length and is finite.

innerₕ(D₋ₓ(f), v) is \(\sum_i |\square_i| \, (D_{-x}f)_i \, v_i\): the operator acts on the source, producing another grid function, which is then integrated against the test function. It agrees entry for entry with applying the numeric operator first.

Bilinear coupling across meshes#

A bilinear term coupling two leaves over different meshes has no assembly unless the trial or the test function is wrapped in πₕ. Unwrapped, it is refused:

try
    assemble(form(Vh, Vh, (u, v) -> innerₕ(u(2), v(1))))
catch e
    showerror(stdout, e)
end

ArgumentError: a bilinear term coupling two leaves over different meshes has no assembly: the trial leaf has (4, 4) points and the test leaf (8, 8), so an index on one names no point on the other. Got Bramble.BilinearProduct{2, Bramble.InnerH, Bramble.IndexedTrialFunction{2}, Bramble.IndexedTestFunction{2}}. Couple leaves that share a mesh, or wrap the trial or the test function in an interpolation operator: `πₕ(u)`, `πₕ(v)`.

A coupled block is assembled by walking the test leaf's grid and reading the trial column out of the same index space, so the two leaves must agree on what an index means. Index \((3, 3)\) on an 8-by-8 grid and on a 4-by-4 grid name different points. The error is raised when the matrix is built, because which leaves a term couples is a question about the spaces it is assembled against. allocate_system_matrix refuses it too.

Wrapping the trial function interpolates it onto the test leaf's grid, and the same block assembles:

Ac_pi = assemble(form(Vh, Vh, (u, v) -> innerₕ(πₕ(u(2)), v(1))))
using SparseArrays
nb = ndofs(Wbig)
B = Ac_pi[1:nb, nb+1:end]
nnz(B)

256

Applied to the small leaf's values, the off-diagonal block gives the vector the source form assembled above:

x₂ = parent(Rₕ(Wsmall, x -> x[1] + x[2]))
@assert nnz(B) > 0 && maximum(abs, b_plain[1:nb]) > 0
@assert isapprox(B * x₂, b_plain[1:nb]; atol = 1e-14)
maximum(abs, B * x₂ .- b_plain[1:nb])

6.938893903907228e-18

Leaves that share a mesh are unaffected. That covers every space built by repeating one space (Wₕ^Val(2)), off-diagonal blocks included.