FiniteVolumeMethod.jl v1.0: Solving conservation PDEs in 2D, now with systems and internal conditions

I have just released v1.0 of my FiniteVolumeMethod.jl package, which is a package for solving PDEs of the form

\frac{\partial u(\boldsymbol x, t)}{\partial t} + \boldsymbol\nabla\boldsymbol\cdot\boldsymbol q(\boldsymbol x, t, u) = S(\boldsymbol x, t, u),

with additional support provided for systems of PDEs of the above form, steady state problems where \partial u/\partial t = 0. Previous versions of FiniteVolumeMethod.jl only support scalar PDEs instead of systems. The package uses my other package DelaunayTriangulation.jl for defining the mesh. Additionally, support is now added for inhomogeneous Neumann boundary conditions and for internal conditions. All the new changes are listed in the release. I’ve also rewrote the mathematical details section to be (hopefully) a lot clearer.

There are 15 examples in documentation which cover a range of examples, covering e.g. domains with holes, mixed boundary conditions, interpolation with NaturalNeighbours.jl. Here are two examples showing the support for internal conditions and for systems of PDEs. For an explanation, see the corresponding tutorials in the docs.

Laplace’s Equation with Internal Conditions

As a first example, I’ll solve Laplace’s equation with an internal constraint so that the solution is zero on some line inside the domain. In particular:

\begin{equation*} \begin{aligned} \boldsymbol\nabla^2 u &= 0 & \quad (x, y) \in [0, 1]^2, \\ u(0, y) &= 100y, & 0 \leq y \leq 1, \\ u(1, y) &= 100y, & 0 \leq y \leq 1, \\ u(x, 0) &= 0, & 0 \leq x \leq 1, \\ u(x, 1) &= 100, & 0 \leq x\leq 1, \\ u(1/2, y) &= 0, & 0 \leq y \leq 2/5. \end{aligned} \end{equation*}

This last condition u(1/2, y)=0 for 0 \leq y \leq 2/5 is the interior constraint. Here is how we solve it. The main difference from previous versions is that InternalConditions is needed to enforce constraints at nodes in the interior.

# The mesh. Need to add points on the line to enforce the interior constraints.
using DelaunayTriangulation, FiniteVolumeMethod
tri = triangulate_rectangle(0, 1, 0, 1, 50, 50, single_boundary=false)
new_points = LinRange(0, 2 / 5, 250)
for y in new_points
    add_point!(tri, 1 / 2, y)
end
refine!(tri, max_area=1e-4)
mesh = FVMGeometry(tri)

# The boundary conditions 
bc_bot = (x, y, t, u, p) -> zero(u)
bc_right = (x, y, t, u, p) -> oftype(u, 100y) # helpful to have each bc return the same type
bc_top = (x, y, t, u, p) -> oftype(u, 100)
bc_left = (x, y, t, u, p) -> oftype(u, 100y)
bcs = (bc_bot, bc_right, bc_top, bc_left)
types = (Dirichlet, Dirichlet, Dirichlet, Dirichlet)
BCs = BoundaryConditions(mesh, bcs, types)

# The internal conditions 
function find_all_points_on_line(tri)
    vertices = Int[]
    for i in each_solid_vertex(tri)
        x, y = get_point(tri, i)
        if x == 1 / 2 && 0 ≤ y ≤ 2 / 5
            push!(vertices, i)
        end
    end
    return vertices
end
vertices = find_all_points_on_line(tri)
ICs = InternalConditions((x, y, t, u, p) -> zero(u),
    dirichlet_nodes=Dict(vertices .=> 1))

# The problem 
initial_condition = zeros(num_points(tri)) # initial guess for the steady state
for i in each_solid_vertex(tri)
    x, y = get_point(tri, i)
    initial_condition[i] = ifelse(x == 1 / 2 && 0 ≤ y ≤ 2 / 5, 0, 100y)
end
diffusion_function = (x, y, t, u, p) -> one(u) # ∇²u = ∇⋅[D∇u], D = 1
final_time = Inf
prob = FVMProblem(mesh, BCs, ICs;
    diffusion_function,
    initial_condition,
    final_time)
steady_prob = SteadyFVMProblem(prob)

# The solution 
using SteadyStateDiffEq, LinearSolve, OrdinaryDiffEq
sol = solve(steady_prob, DynamicSS(TRBDF2(linsolve=KLUFactorization())))
using CairoMakie 
fig, ax, sc = tricontourf(tri, sol.u, levels=LinRange(0, 100, 28))
tightlimits!(ax)
fig

Keller-Segel Model of Chemotaxis

Now I’ll solve the following Keller-Segel model of chemotaxis (this example takes quite a long time to run, so I wouldn’t recommend you run it yourself - just has a nice animation). For systems, you need to provide an FVMProblem for each equation.

\begin{equation*} \begin{aligned} \frac{\partial u}{\partial t} &= \boldsymbol\nabla^2 u - \boldsymbol\nabla\boldsymbol\cdot\left(\frac{cu}{1+u^2}\boldsymbol\nabla v\right) + u(1-u), \\ \frac{\partial v}{\partial t} &= D\boldsymbol\nabla^2 v + u - av, \\ \end{aligned} \end{equation*}

where (x, y) \in [0, 100]^2 and we have homogeneous Neumann boundary conditions. Note that I’ve had to implement this problem by rewriting in the conservation form.

# Defining the problem
using FiniteVolumeMethod, DelaunayTriangulation
tri = triangulate_rectangle(0, 100, 0, 100, 250, 250, single_boundary=true)
mesh = FVMGeometry(tri)
bc_u = (x, y, t, (u, v), p) -> zero(u)
bc_v = (x, y, t, (u, v), p) -> zero(v)
BCs_u = BoundaryConditions(mesh, bc_u, Neumann)
BCs_v = BoundaryConditions(mesh, bc_v, Neumann)
q_u = (x, y, t, (αu, αv), (βu, βv), (γu, γv), p) -> begin
    u = αu * x + βu * y + γu
    ∇u = (αu, βu)
    ∇v = (αv, βv)
    χu = p.c * u / (1 + u^2)
    _q = χu .* ∇v .- ∇u
    return _q
end
q_v = (x, y, t, (αu, αv), (βu, βv), (γu, γv), p) -> begin
    ∇v = (αv, βv)
    _q = -p.D .* ∇v
    return _q
end
S_u = (x, y, t, (u, v), p) -> begin
    return u * (1 - u)
end
S_v = (x, y, t, (u, v), p) -> begin
    return u - p.a * v
end
q_u_parameters = (c=4.0,)
q_v_parameters = (D=1.0,)
S_v_parameters = (a=0.1,)
u_initial_condition = 0.01rand(num_points(tri))
v_initial_condition = zeros(num_points(tri))
final_time = 1000.0
u_prob = FVMProblem(mesh, BCs_u;
    flux_function=q_u, flux_parameters=q_u_parameters,
    source_function=S_u,
    initial_condition=u_initial_condition, final_time=final_time)
v_prob = FVMProblem(mesh, BCs_v;
    flux_function=q_v, flux_parameters=q_v_parameters,
    source_function=S_v, source_parameters=S_v_parameters,
    initial_condition=v_initial_condition, final_time=final_time)
prob = FVMSystem(u_prob, v_prob)

# Solving and visualising
using OrdinaryDiffEq, LinearSolve, CairoMakie
sol = solve(prob, TRBDF2(linsolve=KLUFactorization()), saveat=1.0) # very slow
fig = Figure(fontsize=44)
x = LinRange(0, 100, 250)
y = LinRange(0, 100, 250)
i = Observable(1)
axu = Axis(fig[1, 1], width=600, height=600,
    title=map(i -> L"u(x,~ y,~ %$(sol.t[i]))", i), xlabel=L"x", ylabel=L"y")
axv = Axis(fig[1, 2], width=600, height=600,
    title=map(i -> L"v(x,~ y,~ %$(sol.t[i]))", i), xlabel=L"x", ylabel=L"y")
u = map(i -> reshape(sol.u[i][1, :], 250, 250), i)
v = map(i -> reshape(sol.u[i][2, :], 250, 250), i)
heatmap!(axu, x, y, u, colorrange=(0.0, 2.5), colormap=:turbo)
heatmap!(axv, x, y, v, colorrange=(0.0, 10.0), colormap=:turbo)
resize_to_layout!(fig)
record(fig, "keller_segel_chemotaxis.mp4", eachindex(sol);
    framerate=60) do _i
    i[] = _i
end

The complete list of examples that I give in the docs are below.

18 Likes

Also just another fun example: Solving mazes with Laplace’s equation:

9 Likes

Great work :clap:

2 Likes

cool

1 Like

I’ve just updated the package so that I also provide some specialised solvers for specific problems, namely:

  1. DiffusionEquations: \partial_tu = \boldsymbol\nabla\boldsymbol\cdot[D(\boldsymbol x)\boldsymbol\nabla u].
  2. MeanExitTimeProblems: \boldsymbol\nabla\boldsymbol\cdot[D(\boldsymbol x)\boldsymbol\nabla T(\boldsymbol x)] = -1.
  3. LinearReactionDiffusionEquations: \partial_tu = \boldsymbol\nabla\boldsymbol\cdot[D(\boldsymbol x)\boldsymbol\nabla u] + f(\boldsymbol x)u.
  4. PoissonsEquation: \boldsymbol\nabla\boldsymbol\cdot[D(\boldsymbol x)\boldsymbol\nabla u] = f(\boldsymbol x).
  5. LaplacesEquation: \boldsymbol\nabla\boldsymbol\cdot[D(\boldsymbol x)\boldsymbol\nabla u] = 0.

This list will hopefully include semilinear systems in the future.

I also derive, in the linked section, the details for all of these methods, so that it’s also clear if anyone ever wants to make their own specialised solver for their problem. These are significantly faster than if you use the more general FVMProblem approach, since the linearity in the diffusion term can be easily exploited. I also fixed a design issue that was accidentally causing FVMSystem problems to take 50x longer than they should… so now the Keller-Segel model problem above can be solved relatively quickly (if you also use CVODE_BDF(linear_solver=:GMRES) from Sundials.jl.

Here’s an example of one of these template problems (see here), solving Poisson’s equation describing a problem with two parallel capacitor plates with a space-varying dielectric function and charge density so that \boldsymbol\nabla\boldsymbol\cdot[\epsilon(\boldsymbol x)\boldsymbol\nabla V(\boldsymbol x)] = -\rho(\boldsymbol x)/\epsilon_0. Here is the setup of the problem. It’s a bit long since \epsilon and \rho have complicated definitions and we need the indices of the two plates for the internal Dirichlet conditions; none of this is new in 1.1.

Setting up the problem
using DelaunayTriangulation, FiniteVolumeMethod, CairoMakie, LinearAlgebra, SpecialFunctions
## The mesh
a, b, c, d = 0.0, 10.0, 0.0, 10.0
e, f = (2.0, 3.0), (8.0, 3.0)
g, h = (2.0, 7.0), (8.0, 7.0)
points = [(a, c), (b, c), (b, d), (a, d), e, f, g, h]
boundary_nodes = [[1, 2], [2, 3], [3, 4], [4, 1]]
edges = Set(((5, 6), (7, 8))) # the two capacitor plates
tri = triangulate(points; boundary_nodes, edges, delete_ghosts=false)
refine!(tri; max_area=1e-4get_total_area(tri))
mesh = FVMGeometry(tri)
## The boundary and internal conditions
zero_f = (x, y, t, u, p) -> zero(x)
bc_f = (zero_f, zero_f, zero_f, zero_f)
bc_types = (Dirichlet, Neumann, Dirichlet, Neumann) # bot, right, top, left
BCs = BoundaryConditions(mesh, bc_f, bc_types)
function find_vertices_on_plates(tri)
    lower_plate = Int[]
    upper_plate = Int[]
    for i in each_solid_vertex(tri)
        x, y = get_point(tri, i)
        in_range = 2 ≤ x ≤ 8
        if in_range
            y == 3 && push!(lower_plate, i)
            y == 7 && push!(upper_plate, i)
        end
    end
    return lower_plate, upper_plate
end
lower_plate, upper_plate = find_vertices_on_plates(tri)
dirichlet_nodes = Dict(
    (lower_plate .=> 1)...,
    (upper_plate .=> 2)...
)
internal_f = ((x, y, t, u, p) -> -one(x), (x, y, t, u, p) -> one(x))
ICs = InternalConditions(internal_f, dirichlet_nodes=dirichlet_nodes)
## The dielectric function and charge density
function dist_to_line(p, a, b)
    ℓ² = norm(a .- b)^2
    t = max(0.0, min(1.0, dot(p .- a, b .- a) / ℓ²))
    pv = a .+ t .* (b .- a)
    return norm(p .- pv)
end
function dist_to_plates(x, y)
    p = (x, y)
    a1, b1 = (2.0, 3.0), (8.0, 3.0)
    a2, b2 = (2.0, 7.0), (8.0, 7.0)
    d1 = dist_to_line(p, a1, b1)
    d2 = dist_to_line(p, a2, b2)
    return min(d1, d2)
end
function dielectric_function(x, y, p)
    r = dist_to_plates(x, y)
    return 1 + (p.ϵ₀ - 1) * (1 + erf(r / p.Δ)) / 2
end
function charge_density(x, y, p)
    r = dist_to_plates(x, y)
    return p.Q / (2π) * exp(-r^2 / 2)
end
function plate_source_function(x, y, p)
    ρ = charge_density(x, y, p)
    return -ρ / p.ϵ₀
end
diffusion_parameters = (ϵ₀=8.8541878128e-12, Δ=4.0)
source_parameters = (ϵ₀=8.8541878128e-12, Q=1e-6)

Here is the new part. Rather than using FVMProblem, you can use PoissonsEquation:

diffusion_parameters = (ϵ₀=8.8541878128e-12, Δ=4.0)
source_parameters = (ϵ₀=8.8541878128e-12, Q=1e-6)
prob = PoissonsEquation(mesh, BCs, ICs;
    diffusion_function=dielectric_function,
    diffusion_parameters=diffusion_parameters,
    source_function=plate_source_function,
    source_parameters=source_parameters)
sol = solve(prob, KLUFactorization()) # steady state problems are a linear solve

This takes about 9.68 ms, while the corresponding FVMProblem takes 352 ms - a great speedup! In the docs, I also show how the voltage and electric fields can be visualised, using NaturalNeighbours.jl to compute the gradients, giving the figure below.

6 Likes