# 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
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)

# The solution
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.

22 Likes
10 Likes

Great work

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.

7 Likes

Now that DelaunayTriangulation v1.0 has been released, support for curve-bounded domains now exists. The documentation of FiniteVolumeMethod.jl has now been updated, where applicable, to use these new features. This allows for us to specify, for example, a circle as a function rather than manually discretising it. For example, instead of specifying the circle as

using DelaunayTriangulation, CairoMakie
R = 2.0
θ = LinRange(0, 2π, 250)
x = @. R * cos(θ)
y = @. R * sin(θ)
x[end] = x[1]
y[end] = y[1]
boundary_nodes, points = convert_boundary_points_to_indices(x, y)
tri = triangulate(points; boundary_nodes)
refine!(tri; max_area=1e-3 * π * R^2)


you can do

circle = CircularArc((0.0, R), (0.0, R), (0.0, 0.0))
points = NTuple{2,Float64}[] # still need to provide some points, even if empty
tri2 = triangulate(points; boundary_nodes = [circle])
refine!(tri2; max_area=1e-3 * π * R^2)


which can significantly improve the quality of the mesh and, thus, the solution of your PDE - just see how the quality of the mesh around the boundary changes here:

fig = Figure()
ax=Axis(fig[1,1],width=400,height=400,title="Manual discretisation")
triplot!(ax, tri, color=:black)
ax=Axis(fig[1,2],width=400,height=400,title="Automatic discretisation")
triplot!(ax, tri2, color=:black)
resize_to_layout!(fig)
fig


Please see the docs for more examples, and the docs of DelaunayTriangulation.jl for more examples of specifying curve-bounded domains.

5 Likes