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

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.