Point sources and boundary flux recovery

A well injecting into an aquifer, or a thermal probe in a slab, is a source concentrated at a point rather than spread over a volume. dirac writes one directly in a linear form, and reaction reads back the flux the boundary condition had to supply, which is how the discrete solution can be checked against conservation rather than against a norm alone.

This page is generated from docs/src/examples/point_sources_flux.jl by Literate.jl. The same file runs under test/examples/pages.jl, where the lines marked #src – the assertions the page does not show – execute.

Problem

\[-\nabla \cdot (k \nabla u) = \sum_k Q_k \, \delta(x - x_k) \text{ in } \Omega = (0,1)^2, \qquad u = 0 \text{ on } \partial\Omega\]

with unit conductivity k = 1. Read u as hydraulic head and Q_k as pumping rates.

One well

A single well of strength $Q$ at $x_0$. The Dirichlet condition holds the head at zero on the whole boundary, so all of the injected water has to leave through it.

using Bramble

Q = 3.0
x₀ = (0.35, 0.65)

Ω = domain(interval(0.0, 1.0) × interval(0.0, 1.0),
    :left => :xmin, :right => :xmax, :bottom => :ymin, :top => :ymax)
Ωₕ = mesh(Ω, (61, 61), (true, true))
Wₕ = gridspace(Ωₕ)

a = form(Wₕ, Wₕ, (u, v) -> inner₊(∇ₕ(u), ∇ₕ(v)))
l = form(Wₕ, v -> innerₕ(dirac(x₀, Q), v))

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

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

sum(assemble(l)), maximum(parent(uₕ))

(3.0, 2.3579913160543198)

The assembled load sums to Q exactly. A point evaluation is not weighted by a cell measure the way a density is: dirac puts the strength on the node when the point is one, and splits it across the $2^D$ corners of the containing cell by multilinear weights when it is not, which is what preserves the total.

Checking it against the Green's function

On the unit square with homogeneous Dirichlet data the exact solution is the Green's function, and separating variables gives a series that converges exponentially:

\[G(x, y; x_0, y_0) = \sum_{n \ge 1} \frac{2 \sin(n\pi x)\sin(n \pi x_0) \sinh(n\pi y_<)\sinh(n\pi(1 - y_>))}{n\pi \sinh(n\pi)}\]

with $y_< = \min(y, y_0)$ and $y_> = \max(y, y_0)$. Away from the singularity, where the discrete solution cannot be expected to resolve a logarithm, the two agree:

function green(x, y, x₀, y₀; nterms = 60)
    total = 0.0
    ylo, yhi = minmax(y, y₀)
    for n in 1:nterms
        k = n * π
        total += 2 * sin(k * x) * sin(k * x₀) * sinh(k * ylo) * sinh(k * (1 - yhi)) /
                 (k * sinh(k))
    end
    return total
end

xs, ys = points(Ωₕ)
probes = ((0.15, 0.15), (0.8, 0.2), (0.5, 0.9))

errors = map(probes) do p
    i = argmin(abs.(xs .- p[1]))
    j = argmin(abs.(ys .- p[2]))
    return abs(uₕ[i, j] - Q * green(xs[i], ys[j], x₀...))
end


errors

(3.2249671640671718e-6, 3.071999427165656e-5, 0.00010195851989622251)

What the boundary had to supply

reaction recovers that flux from the unconstrained form pair and the solved uₕ: the constrained rows of A were overwritten by dirichlet_bc!, so the flux they carried is gone from the assembled system, but reassembling without the condition and reading the residual A uₕ - F gives it back. The sign convention is positive for flux leaving the domain.

flux_total = reaction(a, l, uₕ; marker = :boundary)

3.0000000000000586

Every drop injected leaves through the boundary, to round-off, on any mesh: this is a discrete conservation statement, not a discretization error that shrinks under refinement.

flux_total, flux_total - Q

(3.0000000000000586, 5.861977570020827e-14)

Per side, the split follows the geometry. The well sits above centre and left of it, so the top boundary takes the largest share and the bottom the smallest:

sides = (:left, :right, :bottom, :top)
per_side = [reaction(a, l, uₕ; marker = m) for m in sides]


collect(zip(sides, per_side))

4-element Vector{Tuple{Symbol, Float64}}:
 (:left, 1.106994267883872)
 (:right, 0.3930057321161564)
 (:bottom, 0.3930057321161561)
 (:top, 1.106994267883872)

reaction_density is the pointwise counterpart: a grid function, zero away from the marker, carrying the flux density at each marked point rather than its sum. Weighting it by the cell measures recovers the scalar above, which is the check that the two agree, and it is the form to hand an exporter:

rd = reaction_density(a, l, uₕ; marker = :boundary)
w = weights(Wₕ, Bramble.Innerh())


sum(parent(rd) .* w), maximum(parent(rd))

(3.000000000000057, 251.83558609074717)

Two wells

An injector and a producer of equal strength. dirac takes a vector of points and a vector of strengths, which assembles as one term rather than two:

l₂ = form(Wₕ, v -> innerₕ(dirac([(0.3, 0.5), (0.7, 0.5)], [Q, -Q]), v))

A₂, F₂ = assemble(a, l₂; dirichlet = :boundary => x -> 0.0)
u₂ = element(Wₕ)
u₂ .= A₂ \ F₂

net = reaction(a, l₂, u₂; marker = :boundary)

-5.656933255160368e-15

Nothing accumulates and nothing leaves: what one well injects, the other takes back, and the net boundary flux is zero to round-off. The sides do carry flux, in opposite directions, since the two wells are not at the same place.

net, reaction(a, l₂, u₂; marker = :left), reaction(a, l₂, u₂; marker = :right)

(-5.656933255160368e-15, 1.0451997767208245, -1.0451997767208259)