[ANN] FiniteVolumeMethod.jl

Package Announcements
, , ,
1

Hi all,

I have just released my package FiniteVolumeMethod.jl for solving equations of the form \partial_tu + \nabla \cdot q(x, y, t, u) = R(x, y, t, u) for (x, y) \in \Omega \subset \mathbb R^2 using the finite volume method with a triangular mesh, with an additional easy method for specifying reaction-diffusion equations of type \partial_tu = \nabla \cdot [D(x, y, t, u)\nabla u(x, y, t)] + R(x, y, t, u). The triangular mesh should come from DelaunayTriangulation.jl.

The package uses three main structs for representing information about the PDE. Firstly FVMGeometry which represents the discretised geometry \Omega. Second, BoundaryConditions for the boundary conditions of the PDE. These boundary conditions can be specified on individual parts of the boundary, allowing mixed boundary conditions to be easily specified. Lastly, FVMProblem, which specifies the functions involved, initial condition, etc. Once these are defined, a solve gives the solution in the form of a solution object from DifferentialEquations.jl.

There are six full examples of how these are used in the README of the package, but below is a simple example where I solve

\begin{equation*} \begin{array}{rcll} \dfrac{\partial u(x, y, t)}{\partial t} &=& D\boldsymbol{\nabla} \boldsymbol{\cdot} [u \boldsymbol{\nabla} u] + \lambda u(1-u), & 0 < x < a, 0 < y < b, t > 0, \\ u(x, 0, t) & = & 1 & 0 < x < a, t > 0, \\ u(x, b, t) & = & 0 & 0 < x < a , t > 0, \\ \dfrac{\partial u(0, y, t)}{\partial x} & = & 0 & 0 < y < b, t > 0, \\ \dfrac{\partial u(a, y, t)}{\partial x} & = & 0 & 0 < y < b, t > 0, \\ u(x, y, 0) & = & f(y) & 0 \leq x \leq a, 0 \leq y \leq b, \end{array} \end{equation*}

with f(y)=0, a = 3 and b = 40, and a structured triangulation is used.

using FiniteVolumeMethod
using DelaunayTriangulation
using LinearSolve 
using OrdinaryDiffEq 

## Step 1: Define the mesh 
a, b, c, d, Nx, Ny = 0.0, 3.0, 0.0, 40.0, 60, 80
T, adj, adj2v, DG, points, BN = triangulate_structured(a, b, c, d, Nx, Ny; return_boundary_types=true)
mesh = FVMGeometry(T, adj, adj2v, DG, points, BN)

## Step 2: Define the boundary conditions 
a₁ = ((x, y, t, u::T, p) where {T}) -> one(T)
a₂ = ((x, y, t, u::T, p) where {T}) -> zero(T)
a₃ = ((x, y, t, u::T, p) where {T}) -> zero(T)
a₄ = ((x, y, t, u::T, p) where {T}) -> zero(T)
bc_fncs = (a₁, a₂, a₃, a₄)
types = (:D, :N, :D, :N) # :D is Dirichlet, :N is Neumann
BCs = BoundaryConditions(mesh, bc_fncs, types, BN)

## Step 3: Define the actual PDE  
f = ((x::T, y::T) where {T}) -> zero(T)
diff_fnc = (x, y, t, u, p) -> p * u
reac_fnc = (x, y, t, u, p) -> p * u * (1 - u)
D, λ = 0.9, 0.99
diff_parameters = D
reac_parameters = λ
final_time = 50.0
u₀ = [f(points[:, i]...) for i in axes(points, 2)]
prob = FVMProblem(mesh, BCs; diffusion_function=diff_fnc, reaction_function=reac_fnc,
    diffusion_parameters=diff_parameters, reaction_parameters=reac_parameters,
    initial_condition=u₀, final_time)

## Step 4: Solve
alg = TRBDF2(linsolve=KLUFactorization())
sol = solve(prob, alg; saveat=0.5)

The interface also makes it easy to plot this solution, simply using mesh!. This is shown in the README.

Plenty more detail in the README, along with the mathematical details.

Thoughts and feedback welcome!
Thanks,
Daniel.

17 Likes
2

It would be cool to get a few lines in GitHub - JuliaPDE/SurveyofPDEPackages: Survey of the packages of the Julia ecosystem for solving partial differential equations!

5 Likes
3

Thanks for the suggestion! Done.

4

Exciting contribution to the PDE world in Julia. May I suggest to break up the code into more files instead of one super long file.

2 Likes
5

Thanks! Yeah, I originally had that and then had everything back in a file once things were cleaned up. Will fix that now, thanks for the reminder.

6

A nice update as of v0.4: The code now supports steady problems, and with recent updates to my other package DelaunayTriangulation.jl I have added support for multiply-connected domains (see the docs for DelaunayTriangulation.jl here for more discussion on how these domains are represented). To showcase both of these features, the following code solves the mean exit time problem

\begin{equation*} \begin{array}{rcll} \boldsymbol{\nabla}^2 T & = & -1/D & \boldsymbol x \in \Omega, \\[6pt] T(\boldsymbol x) & = & 0 & \|\boldsymbol{x}\| = R_2(\theta), \\[6pt] \boldsymbol{\nabla} T(\boldsymbol x) \boldsymbol{\cdot} \hat{\boldsymbol{n}} & = & 0 & \|\boldsymbol{x}\| = R_1(\theta), \end{array} \end{equation*}

where R_1(\theta) = 2(1 + \varepsilon g_1(\theta)), R_2(\theta) = 3(1 + \varepsilon g_2(\theta)), \varepsilon = 1/20, g_1(\theta) = \sin(3\theta) + \cos(5\theta), g_2(\theta) = \cos(3\theta), and \Omega is the region R_1(\theta) < r < R_2(\theta), where \theta is the polar angle of \boldsymbol x in each case. The code is essentially as it would be for an unsteady problem, except now initial_condition is interpreted as an initial guess for the nonlinear solvers compatible with NonlinearSolve.jl, and final_time is just \infty. The two separate boundaries are represented by placing the set of outer boundary segments first in counter-clockwise order, and then the inner boundary segments clockwise (so that the entire boundary is positively oriented). You just set steady = true in FVMProblem and pass a nonlinear solver algorithm into solve and it should just work:

using FiniteVolumeMethod
using DelaunayTriangulation
using CairoMakie
using OrdinaryDiffEq
using SteadyStateDiffEq
using LinearSolve
using Krylov

## Define the multiply-connected geometry
R₁ = 2.0
R₂ = 3.0
g₁ = θ -> sin(3θ) + cos(5θ)
g₂ = θ -> cos(3θ)
θ = LinRange(0, 2π, 250)
ε = 1/20
inner_R = R₁ .* (1 .+ ε .* g₁.(θ))
outer_R = R₂ .* ( 1 .+ ε .* g₂.(θ))
x = [
    [outer_R .* cos.(θ)],
    [reverse(inner_R .* cos.(θ))] # inner boundaries are clockwise
]
y = [
    [outer_R .* sin.(θ)],
    [reverse(inner_R .* sin.(θ))] # inner boundaries are clockwise
]
GMSH_PATH = "./gmsh-4.9.4-Windows64/gmsh.exe"

tri = generate_mesh(x, y, 0.1; gmsh_path=GMSH_PATH)

## Define the FVMProblem
mesh = FVMGeometry(tri)
outer_bc = (x, y, t, u, p) -> zero(u)
inner_bc = (x, y, t, u, p) -> zero(u)
type = [:D, :N] # absorbing outer, reflecting inner
BCs = BoundaryConditions(mesh, [outer_bc, inner_bc], type)
D = 6.25e-4
initial_guess = zeros(num_points(tri)) # start the initial guess at the zero vector
diffusion_function = (x, y, t, u, D) -> D
diffusion_parameters = D
reaction = (x, y, t, u, p) -> one(u)
final_time = Inf # doesn't actually get used, but we still need to provide it
prob = FVMProblem(mesh, BCs;
    diffusion_function,
    diffusion_parameters,
    reaction_function=reaction,
    final_time=final_time,
    initial_condition=initial_guess,
    steady=true)

## Now solve - we use a SteadyState solver with a Krylov method
sol = solve(prob, DynamicSS(TRBDF2(linsolve=KrylovJL_GMRES())), parallel = true)

## Visualise
fig = Figure(fontsize=38)
ax = Axis(fig[1, 1], xlabel=L"x", ylabel=L"y", width=600, height=300)
pt_mat = Matrix(get_points(tri)')
T_mat = [T[j] for T in each_triangle(tri), j in 1:3]
msh = mesh!(ax, pt_mat, T_mat, color=sol, colorrange=(0, 900))
Colorbar(fig[1, 2], msh)
resize_to_layout!(fig)

perturbed_annulus_mean_exit_time
perturbed_annulus_mean_exit_time832×425 83.5 KB

I include plenty of examples in the documentation.

5 Likes
7

A nice update as of v0.4.3 is that, thanks to upgrades to my other package DelaunayTriangulation.jl, is that there is no longer any need to use Gmsh so it has been completely removed from all examples. Hopefully that is less prohibitive. For example, the annulus example’s mesh above could be generated with

R₁ = 2.0
R₂ = 3.0
g₁ = θ -> sin(3θ) + cos(5θ)
g₂ = θ -> cos(3θ)
θ = (collect ∘ LinRange)(0, 2π, 250)
θ[end] = θ[begin] # close the boundary properly
ε = 1 / 20
inner_R = R₁ .* (1 .+ ε .* g₁.(θ))
outer_R = R₂ .* (1 .+ ε .* g₂.(θ))
x = [
    [outer_R .* cos.(θ)],
    [reverse(inner_R .* cos.(θ))] # inner boundaries are clockwise
]
y = [
    [outer_R .* sin.(θ)],
    [reverse(inner_R .* sin.(θ))] # inner boundaries are clockwise
]
boundary_nodes, points = convert_boundary_points_to_indices(x, y)
tri = triangulate(points; boundary_nodes)
A = get_total_area(tri)
refine!(tri; max_area=1e-4A)
6 Likes
8

In case it is useful to anyone, I also have made public a one-dimensional version of this package at GitHub - DanielVandH/FiniteVolumeMethod1D.jl: Implementation of the finite volume method in 1D.. It’s pretty much for the same type of equations as in FiniteVolumeMethod.jl, although I support Robin boundary conditions also. It’s completely documented and tested and only has a few dependencies (SparseArrays.jl, SciMLBase.jl, CommonSolve.jl).

I probably won’t register FiniteVolumeMethod1D.jl, so just posting this in case anyone wants it. Not registering since maybe it’s better one day merged into FiniteVolumeMethod.jl (less inclined to do this, though - maybe something for 1.0 if I ever do decide), and there’s probably other packages that handle these types of PDEs more effectively / efficiently, but if there’s any interest feel free to let me know.

Edit: I’ve decided, and I’ve now registered FiniteVolumeMethod1D.

2 Likes