I have just released v1.0 of my FiniteVolumeMethod.jl package, which is a package for solving PDEs of the form
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:
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.
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.
- Diffusion Equation on a Square Plate
- Diffusion Equation in a Wedge with Mixed Boundary Conditions
- Reaction Diffusion Equation with a Time-dependent Dirichlet Boundary Condition on a Disk
- Porous-Medium Equation
- Porous-Fisher Equation and Travelling Waves
- Piecewise Linear and Natural Neighbour Interpolation for an Advection-Diffusion Equation
- Helmholtz Equation with Inhomogeneous Boundary Conditions
- Laplace’s Equation with Internal Dirichlet Conditions
- Equilibrium Temperature Distribution with Mixed Boundary Conditions and using EnsembleProblems
- A Reaction-Diffusion Brusselator System of PDEs
- Gray-Scott Model: Turing Patterns from a Coupled Reaction-Diffusion System
- Diffusion Equation on an Annulus
- Mean Exit Time
- Solving Mazes with Laplace’s Equation
- Keller-Segel Model of Chemotaxis
