Custom coupled physics (weak forms + DSL)

Overview

Beyond the built-in physics, any coupled system can be written by hand as a BilinearForm over a MultiFieldFESpace — the extensibility point for new physics. This tutorial builds a two-field coupled diffusion problem (e.g. two inter- converting species) from scratch, in both the functional and DSL styles.

using Rasmah, LinearAlgebra

mesh = tetrahedralize_box(1.0, 1.0, 1.0, 3, 3, 3)
n = size(mesh.nodes, 2)

U1 = FESpace(mesh, ReferenceFE(3, 1))   # species 1 (scalar)
U2 = FESpace(mesh, ReferenceFE(3, 1))   # species 2 (scalar)
V  = MultiFieldFESpace([U1, U2]; style = :block)

MultiFieldFESpace{Tuple{FESpace{TetMesh{Matrix{Float64}, Matrix{Int64}, Nothing}, ReferenceFE, ScalarValue, Float64}, FESpace{TetMesh{Matrix{Float64}, Matrix{Int64}, Nothing}, ReferenceFE, ScalarValue, Float64}}, BlockStyle, Tuple{Matrix{Int64}, Matrix{Int64}}}((FESpace{TetMesh{Matrix{Float64}, Matrix{Int64}, Nothing}, ReferenceFE, ScalarValue, Float64}(TetMesh{Matrix{Float64}, Matrix{Int64}, Nothing}([0.0 0.3333333333333333 … 0.6666666666666666 1.0; 0.0 0.0 … 1.0 1.0; 0.0 0.0 … 1.0 1.0], [1 1 … 43 43; 2 6 … 59 60; 6 5 … 60 44; 22 22 … 64 64], nothing), ReferenceFE(3, :simplex, 1), ScalarValue, [1 1 … 43 43; 2 6 … 59 60; 6 5 … 60 44; 22 22 … 64 64], [1 1 … 43 43; 2 6 … 59 60; 6 5 … 60 44; 22 22 … 64 64], [0.0 0.3333333333333333 … 0.6666666666666666 1.0; 0.0 0.0 … 1.0 1.0; 0.0 0.0 … 1.0 1.0], 64), FESpace{TetMesh{Matrix{Float64}, Matrix{Int64}, Nothing}, ReferenceFE, ScalarValue, Float64}(TetMesh{Matrix{Float64}, Matrix{Int64}, Nothing}([0.0 0.3333333333333333 … 0.6666666666666666 1.0; 0.0 0.0 … 1.0 1.0; 0.0 0.0 … 1.0 1.0], [1 1 … 43 43; 2 6 … 59 60; 6 5 … 60 44; 22 22 … 64 64], nothing), ReferenceFE(3, :simplex, 1), ScalarValue, [1 1 … 43 43; 2 6 … 59 60; 6 5 … 60 44; 22 22 … 64 64], [1 1 … 43 43; 2 6 … 59 60; 6 5 … 60 44; 22 22 … 64 64], [0.0 0.3333333333333333 … 0.6666666666666666 1.0; 0.0 0.0 … 1.0 1.0; 0.0 0.0 … 1.0 1.0], 64)), BlockStyle(), [0, 64, 128], Int64[], 128, ([1 1 … 43 43; 2 6 … 59 60; 6 5 … 60 44; 22 22 … 64 64], [65 65 … 107 107; 66 70 … 123 124; 70 69 … 124 108; 86 86 … 128 128]))

The weak form

For −∇²u₁ + c(u₁ − u₂) = f₁, −∇²u₂ + c(u₂ − u₁) = f₂, the bilinear form has four blocks (the off-diagonal ones are the ±c u v coupling):

c = 1.0
a = BilinearForm(V, V, [
    (1, 1, (x, u, ∇u, v, ∇v) -> dot(∇u, ∇v) + c * u * v),
    (1, 2, (x, u, ∇u, v, ∇v) -> -c * u * v),
    (2, 1, (x, u, ∇u, v, ∇v) -> -c * u * v),
    (2, 2, (x, u, ∇u, v, ∇v) -> dot(∇u, ∇v) + c * u * v),
])
K = assemble_matrix(a)
println(size(K), " block-coupled operator")

(128, 128) block-coupled operator

Each (trial_field, test_field, integrand) triple places one block; the integrand receives the trial value/gradient (u, ∇u) and the test value/gradient (v, ∇v).

The DSL

The @law / ∫ / dΩ sugar compiles to the same integrand closures:

f1 = (x, u, ∇u, v, ∇v) -> dot(∇u, ∇v) + c * u * v
f2 = @law((x, u, ∇u, v, ∇v), dot(∇u, ∇v) + c * u * v)
@assert assemble_matrix(BilinearForm(U1, U1, f1)) ≈
        assemble_matrix(BilinearForm(U1, U1, f2))
println("DSL ≡ functional API")

DSL ≡ functional API

Solve and matrix-free apply

l = LinearForm(V, [(1, (x, v, ∇v) -> x[1] * v)])
u = solve(AffineFEOperator(a, l), Int[])          # all-free solve
println("max |u| = ", round(maximum(abs, u); digits = 4))

v = rand(ndofs(V))
@assert matvec(a, v) ≈ K * v
println("matrix-free matvec matches assembly")

max |u| = 3.953733962284856e13
matrix-free matvec matches assembly

The same pattern extends to any number of fields and to nonlinear operators (FEOperator residual + matrix-free jacobian_v), which is how the built-in physics are implemented.