# [ANN] FiniteVolumeMethod.jl

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)

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

7 Likes

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!

2 Likes

Thanks for the suggestion! Done.

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.

1 Like

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.