Coupled nonlinear reaction-diffusion system

Two species, coupled through a quadratic reaction term, solved with Newton's method the same way the nonlinear Poisson example does — except now the Jacobian differentiates through a composite space's assembly, not a scalar one. Every number and every plot below was produced by the code shown.

Problem

\[\begin{aligned} -\Delta u + u + uv &= f_1 \\ -\Delta v + v - uv &= f_2 \end{aligned} \qquad \text{in } \Omega = (0,1)^2, \qquad u = v = 0 \text{ on } \partial\Omega,\]

a predator-prey-shaped coupling without the time derivative: u grows through the interaction term, v is depleted by it. The manufactured solutions vanish on the boundary already, so homogeneous Dirichlet data is all that is needed:

using Bramble
using Random

u_ex(x) = sin(π * x[1]) * sin(π * x[2])
v_ex(x) = sin(2π * x[1]) * sin(2π * x[2])
f1(x) = 2π^2 * u_ex(x) + u_ex(x) + u_ex(x) * v_ex(x)
f2(x) = 8π^2 * v_ex(x) + v_ex(x) - u_ex(x) * v_ex(x)

Random.seed!(20260903)
Ω = domain(interval(0.0, 1.0) × interval(0.0, 1.0))
Ωₕ = mesh(Ω, (24, 24), (false, false))
Wₕ = gridspace(Ωₕ)
Vₕ = Wₕ^Val(2)

bcs = dirichlet_constraints(Ω, :boundary => x -> 0.0)
f1ₕ = element(Wₕ)
avgₕ!(f1ₕ, f1)
f2ₕ = element(Wₕ)
avgₕ!(f2ₕ, f2)
l = form(Vₕ, q -> innerₕ(f1ₕ, q(1)) + innerₕ(f2ₕ, q(2)))
F = assemble(l; dirichlet = bcs)

The seed is what makes the numbers below reproducible: (false, false) draws the interior points from the global RNG, so without it the mesh – and every figure on this page – would differ from build to build, and the suite could not assert what the page prints.

Newton's method on a composite residual

uv is quadratic in the unknowns, so it cannot sit inside a matrix independent of w = (u, v) the way the linear terms can — but it can sit inside a matrix that depends on the current guess, the same trick the nonlinear Poisson example uses for a single species, extended to a second one. Writing the coupling as v_current * u(1) in u's own equation and -u_current * u(2) in v's reproduces uv and -uv exactly once the trial function is evaluated at the current w — which is all A(w) needs to do. Nothing here works out ∂(uv)/∂u and ∂(uv)/∂v by hand; ForwardDiff differentiates through how A itself depends on w automatically:

function coupled_matrix(wₕ)
    u_c, v_c = components(wₕ)
    a = form(Vₕ,
        Vₕ,
        (p, q) -> inner₊(∇₋ₕ(p(1)), ∇₋ₕ(q(1))) + innerₕ(p(1), q(1)) + innerₕ(v_c * p(1), q(1)) +
                  inner₊(∇₋ₕ(p(2)), ∇₋ₕ(q(2))) + innerₕ(p(2), q(2)) - innerₕ(u_c * p(2), q(2)))
    return assemble(a; dirichlet = :boundary)
end

Its Jacobian is sparse for the same reason the nonlinear Poisson example's is — the reaction term couples u and v only pointwise, so it adds nothing to the diffusion stencil's own reach — so the same sparse AD setup applies unchanged, just over twice as many unknowns:

using ForwardDiff, DifferentiationInterface
import SparseConnectivityTracer, SparseMatrixColorings

const sparse_ad = AutoSparse(AutoForwardDiff();
    sparsity_detector = SparseConnectivityTracer.TracerSparsityDetector(),
    coloring_algorithm = SparseMatrixColorings.GreedyColoringAlgorithm())

function residual(w::AbstractVector{T}) where {T}
    wₕ = element(Vₕ, T)
    wₕ .= w
    A = coupled_matrix(wₕ)
    return A * w .- F
end

w = zeros(ndofs(Vₕ))
prep = prepare_jacobian(residual, sparse_ad, w)
J = DifferentiationInterface.jacobian(residual, prep, sparse_ad, w)  # once, for its sparse structure
newton_residuals = Float64[]
for it in 1:20
    r = residual(w)
    push!(newton_residuals, sqrt(sum(abs2, r)))
    newton_residuals[end] < 1e-10 && break
    DifferentiationInterface.jacobian!(residual, J, prep, sparse_ad, w)
    w .-= J \ r
end
length(newton_residuals), newton_residuals

(4, [2.5762825732242143, 0.02192029918716173, 1.2892755987507398e-6, 4.9583421918330345e-14])

Quadratic convergence, same as the single-species case — the composite space changes what the Jacobian differentiates through, not how well Newton converges once it has a correct one.

Skipping the tracer here too

v_c scaling a term routed into block (1,1) is a different leaf's component reaching into this one — jacobian_pattern reads that the same way a form term names a component, U -> U(2) rather than a stencil op, since v_c is read directly rather than averaged first. Block (2,2)'s own u_c dependency is named the same way, U -> U(1):

using ADTypes

u_c0, v_c0 = components(element(Vₕ, 0.0))
a_for_pattern = form(Vₕ,
    Vₕ,
    (p, q) -> inner₊(∇₋ₕ(p(1)), ∇₋ₕ(q(1))) + innerₕ(p(1), q(1)) + innerₕ(v_c0 * p(1), q(1)) +
              inner₊(∇₋ₕ(p(2)), ∇₋ₕ(q(2))) + innerₕ(p(2), q(2)) - innerₕ(u_c0 * p(2), q(2)))

native_ad = AutoSparse(AutoForwardDiff();
    sparsity_detector = ast_sparsity_detector(a_for_pattern, U -> U(2), U -> U(1)),
    coloring_algorithm = SparseMatrixColorings.GreedyColoringAlgorithm())

w_native = zeros(ndofs(Vₕ))
prep_native = prepare_jacobian(residual, native_ad, w_native)
J_native = DifferentiationInterface.jacobian(residual, prep_native, native_ad, w_native)
newton_residuals_native = Float64[]
for it in 1:20
    r = residual(w_native)
    push!(newton_residuals_native, sqrt(sum(abs2, r)))
    newton_residuals_native[end] < 1e-10 && break
    DifferentiationInterface.jacobian!(residual, J_native, prep_native, native_ad, w_native)
    w_native .-= J_native \ r
end
newton_residuals_native

4-element Vector{Float64}:
 2.5762825732242143
 0.02192029918716133
 1.2892755987293357e-6
 4.083592084640994e-14

Same quadratic convergence, no tracing pass paid for it: v_c0/u_c0 are read at w = 0 only to build some concrete BilinearForm — the pattern is a property of a's AST, not of those values. SparseConnectivityTracer's tracer still works here regardless of how coupled_matrix was built, composite space and all — the case to reach for it is a residual whose matrix does not come from a BilinearForm in the first place, which is not this one.

wₕ = element(Vₕ)
wₕ .= w
uₕ, vₕ = components(wₕ)
uexact, vexact = Rₕ(Wₕ, u_ex), Rₕ(Wₕ, v_ex)
norm₁ₕ(uₕ .- uexact), norm₁ₕ(vₕ .- vexact)

(0.019779214731516854, 0.10747520241691955)

Visualizing the solution

Each species is its own 2D scalar field — components(wₕ) gives a view directly onto it, no new solve or copy needed:

Checking the answer

The same nested-random-mesh pattern as every other example, checking each species' own error separately — a routing mistake would show up as one converging correctly while the other silently used the wrong block, which a single combined error could hide. A dense ForwardDiff.jacobian over two coupled species would cost (2n)^2 against the scalar examples' n^2, and was what forced this example to stay at three small refinement levels before switching to sparse AD; with it, this reaches five levels — the same order of tens of thousands of degrees of freedom the linear coupled example reaches at six — in about a second per level:

Random.seed!(20260903)

function coupled_series(; n0::Int = 5, levels::Int)
    Ωc = mesh(Ω, (n0, n0), (false, false))
    hs = Float64[]
    erru, errv = Float64[], Float64[]
    for level in 1:levels
        Wc = gridspace(Ωc)
        Vc = Wc^Val(2)
        bcs_c = dirichlet_constraints(Ω, :boundary => x -> 0.0)
        f1_c = element(Wc)
        avgₕ!(f1_c, f1)
        f2_c = element(Wc)
        avgₕ!(f2_c, f2)
        l_c = form(Vc, q -> innerₕ(f1_c, q(1)) + innerₕ(f2_c, q(2)))
        F_c = assemble(l_c; dirichlet = bcs_c)

        Ac(wₕ) = begin
            u_c, v_c = components(wₕ)
            assemble(
                form(Vc,
                    Vc,
                    (p, q) -> inner₊(∇₋ₕ(p(1)), ∇₋ₕ(q(1))) + innerₕ(p(1), q(1)) +
                              innerₕ(v_c * p(1), q(1)) +
                              inner₊(∇₋ₕ(p(2)), ∇₋ₕ(q(2))) + innerₕ(p(2), q(2)) -
                              innerₕ(u_c * p(2), q(2)));
                dirichlet = :boundary)
        end
        rc(w::AbstractVector{T}) where {T} = begin
            wₕ = element(Vc, T)
            wₕ .= w
            Ac(wₕ) * w .- F_c
        end

        w_c = zeros(ndofs(Vc))
        prep_c = prepare_jacobian(rc, sparse_ad, w_c)
        J_c = DifferentiationInterface.jacobian(rc, prep_c, sparse_ad, w_c)
        for it in 1:20
            r = rc(w_c)
            sqrt(sum(abs2, r)) < 1e-10 && break
            DifferentiationInterface.jacobian!(rc, J_c, prep_c, sparse_ad, w_c)
            w_c .-= J_c \ r
        end
        w_ch = element(Vc)
        w_ch .= w_c
        u_ch, v_ch = components(w_ch)
        uexact_c, vexact_c = Rₕ(Wc, u_ex), Rₕ(Wc, v_ex)

        push!(hs, hₘₐₓ(Ωc))
        push!(erru, norm₁ₕ(u_ch .- uexact_c))
        push!(errv, norm₁ₕ(v_ch .- vexact_c))
        level < levels && iterative_refinement!(Ωc)
    end
    return hs, erru, errv
end

hs, erru, errv = coupled_series(; levels = 5)
order_u = log(erru[end - 1] / erru[end]) / log(hs[end - 1] / hs[end])
order_v = log(errv[end - 1] / errv[end]) / log(hs[end - 1] / hs[end])
(order_u, order_v)

(1.998070809315084, 1.9904359603737625)

order_u > 1.9 && order_v > 1.9

true

Second order for both species, same rate as every other example — the composite space and the quadratic coupling change how the residual and its Jacobian are built, not the discretization's own accuracy once Newton has converged to it.

coupled_series above uses sparse_ad, the tracer, at every level — native_ad's substitution (ast_sparsity_detector(a, U -> U(2), U -> U(1)) in place of sparse_ad's sparsity_detector) works here unchanged too. Not re-run a second time here, the same reason the nonlinear Poisson example does not re-run its own convergence sweep a second time either.