# \[ANN\] FiniteVolumeMethod.jl

**URL:** <https://discourse.julialang.org/t/ann-finitevolumemethod-jl/91395>\
**Category:** Package Announcements\
**Tags:** package, announcement, numerics, differentialequation\
**Created:** [December 7, 2022, 11:10pm UTC](https://discourse.julialang.org/t/ann-finitevolumemethod-jl/91395 "2022-12-07T23:10:31Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![DanielVandH](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielvandh/32/31134_2.png) [@DanielVandH](https://discourse.julialang.org/u/DanielVandH)\
**Post date:** [December 7, 2022, 11:10pm UTC](https://discourse.julialang.org/t/ann-finitevolumemethod-jl/91395/1 "2022-12-07T23:10:31Z")

</div>

Hi all,

I have just released my package [FiniteVolumeMethod.jl](https://github.com/DanielVandH/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.

```julia
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.

---

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [December 8, 2022, 12:57am UTC](https://discourse.julialang.org/t/ann-finitevolumemethod-jl/91395/2 "2022-12-08T00:57:03Z")

</div>

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](https://github.com/JuliaPDE/SurveyofPDEPackages)!

---

<div class="post-metadata">

**Author:** ![DanielVandH](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielvandh/32/31134_2.png) [@DanielVandH](https://discourse.julialang.org/u/DanielVandH)\
**Post date:** [December 8, 2022, 1:47am UTC](https://discourse.julialang.org/t/ann-finitevolumemethod-jl/91395/3 "2022-12-08T01:47:38Z")

</div>

Thanks for the suggestion! Done.

---

<div class="post-metadata">

**Author:** ![kylebeggs](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kylebeggs/32/43348_2.png) [@kylebeggs](https://discourse.julialang.org/u/kylebeggs)\
**Post date:** [December 8, 2022, 2:21am UTC](https://discourse.julialang.org/t/ann-finitevolumemethod-jl/91395/4 "2022-12-08T02:21:21Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![DanielVandH](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielvandh/32/31134_2.png) [@DanielVandH](https://discourse.julialang.org/u/DanielVandH)\
**Post date:** [December 8, 2022, 2:23am UTC](https://discourse.julialang.org/t/ann-finitevolumemethod-jl/91395/5 "2022-12-08T02:23:29Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![DanielVandH](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielvandh/32/31134_2.png) [@DanielVandH](https://discourse.julialang.org/u/DanielVandH)\
**Post date:** [February 27, 2023, 6:36am UTC](https://discourse.julialang.org/t/ann-finitevolumemethod-jl/91395/6 "2023-02-27T06:36:57Z")

</div>

A nice update as of v0.4: The code now supports steady problems, and with recent updates to my other package [DelaunayTriangulation.jl](https://github.com/DanielVandH/DelaunayTriangulation.jl) I have added support for multiply-connected domains (see the docs for DelaunayTriangulation.jl [here](https://danielvandh.github.io/DelaunayTriangulation.jl/stable/interface/interface/) 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:

```julia
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](https://global.discourse-cdn.com/julialang/original/3X/a/7/a774137a209ec57b7bdd815f148448b6a0d7512a.png)

I include plenty of examples in the [documentation](https://danielvandh.github.io/FiniteVolumeMethod.jl/stable/).

---

<div class="post-metadata">

**Author:** ![DanielVandH](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielvandh/32/31134_2.png) [@DanielVandH](https://discourse.julialang.org/u/DanielVandH)\
**Post date:** [April 16, 2023, 2:41pm UTC](https://discourse.julialang.org/t/ann-finitevolumemethod-jl/91395/7 "2023-04-16T14:41:13Z")

</div>

A nice update as of v0.4.3 is that, thanks to upgrades to my other package [DelaunayTriangulation.jl](https://github.com/DanielVandH/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

```julia
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)

```

---

<div class="post-metadata">

**Author:** ![DanielVandH](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielvandh/32/31134_2.png) [@DanielVandH](https://discourse.julialang.org/u/DanielVandH)\
**Post date:** [June 3, 2023, 2:27pm UTC](https://discourse.julialang.org/t/ann-finitevolumemethod-jl/91395/8 "2023-06-03T14:27:16Z")

</div>

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.](https://github.com/DanielVandH/FiniteVolumeMethod1D.jl). 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](https://github.com/DanielVandH/FiniteVolumeMethod.jl/issues/37) 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](https://github.com/DanielVandH/FiniteVolumeMethod1D.jl).