GitHub

Interface

In this section, we describe the basic interface for defining and solving PDEs using this package. This interface will also be made clearer in the tutorials. The basic summary of the discussion below is as follows:

  1. Use FVMGeometry to define the problem's mesh.
  2. Provide boundary conditions using BoundaryConditions.
  3. (Optional) Provide internal conditions using InternalConditions.
  4. Convert the problem into an FVMProblem.
  5. If you want to make the problem steady, use SteadyFVMProblem on the FVMProblem.
  6. If you want a system of equations, construct an FVMSystem from multiple FVMProblems; if you want this problem to be steady, skip step 5 and only now apply SteadyFVMProblem.
  7. Solve the problem using solve.
  8. For a discussion of custom constraints, see the tutorials.
  9. For interpolation, we provide pl_interpolate (but you might prefer NaturalNeighbours.jl - see this tutorial for an example).

FVMGeometry: Defining the mesh

The finite volume method (FVM) requires an underlying triangular mesh, as outlined in the mathematical details section. This triangular mesh is to be defined from DelaunayTriangulation.jl. The FVMGeometry type wraps the resulting Triangulation and computes information about the geometry required for solving the PDEs. The docstring for FVMGeometry is below; the fields of FVMGeometry are public API.

FiniteVolumeMethod.FVMGeometry — Type
FVMGeometry(tri::Triangulation)

This is a constructor for the FVMGeometry struct, which holds the mesh and associated data for the PDE.

Note

It is assumed that all vertices in tri are in the triangulation, meaning v is in tri for each v in DelaunayTriangulation.each_point_index(tri).

Fields

  • triangulation: The underlying Triangulation from DelaunayTriangulation.jl.
  • triangulation_statistics: The statistics of the triangulation.
  • cv_volumes::Vector{Float64}: A Vector of the volumes of each control volume.
  • triangle_props::Dict{NTuple{3,Int},TriangleProperties}: A Dict mapping the indices of each triangle to its [TriangleProperties].
source

The FVMGeometry struct uses TriangleProperties for storing properties of a control volume that intersects a given triangle, defined below. This struct is public API, although it is unlikely you would ever need it.

FiniteVolumeMethod.TriangleProperties — Type
TriangleProperties(shape_function_coefficients, cv_edge_midpoints, cv_edge_normals, cv_edge_lengths)

This is a struct for holding the properties of a control volume's intersection with a triangle.

Fields

  • shape_function_coefficients::NTuple{9,Float64}: The shape function coefficients for the triangle.
  • cv_edge_midpoints::NTuple{3,NTuple{2,Float64}}: The midpoints of the control volume edges. If the triangle is (i, j, k), then the edges are given in the order (i, j), (j, k), and (k, i), where 'edge' refers to the edge joining e.g. the midpoint of the edge (i, j) to the centroid of the triangle.
  • cv_edge_normals::NTuple{3,NTuple{2,Float64}}: The normal vectors to the control volume edges, in the same order as in cv_edge_midpoints.
  • cv_edge_lengths::NTuple{3,Float64}: The lengths of the control volume edges, in the same order as in cv_edge_midpoints.
Notes

The shape function coefficients are defined so that, if s are the coefficients and the triangle is T = (i, j, k), with function values u[i], u[j], and u[k] at the vertices i, j, and k, respectively, then αxₙ + βyₙ + γₙ = u[n] for n = i, j, k, where xₙ and yₙ are the x- and y-coordinates of the nth vertex, respectively, α = s₁u₁ + s₂u₂ + s₃u₃, β = s₄u₁ + s₅u₂ + s₆u₃, and γ = s₇u₁ + s₈u₂ + s₉u₃.

source

BoundaryConditions: Defining boundary conditions

Once a mesh is defined, you need to associate each part of the boundary with a set of boundary nodes. Since you have a Triangulation, the boundary of the mesh already meets the necessary assumptions made by this package about the boundary; these assumptions are simply that they match the specification of a boundary here in DelaunayTriangulation.jl's docs (for example, the boundary points connect, the boundary is positively oriented, etc.).

You can specify boundary conditions using BoundaryConditions, whose docstring is below; the fields of BoundaryConditions are not public API, only this wrapper is.

FiniteVolumeMethod.BoundaryConditions — Type
BoundaryConditions(mesh::FVMGeometry, functions, conditions; parameters=nothing)

This is a constructor for the BoundaryConditions struct, which holds the boundary conditions for the PDE. See also Conditions (which FVMProblem wraps this into), ConditionType, and InternalConditions.

Arguments

  • mesh::FVMGeometry

The mesh on which the PDE is defined.

  • functions::Union{<:Tuple,<:Function}

The functions that define the boundary conditions. The ith function should correspond to the part of the boundary of the mesh corresponding to the ith boundary index, as defined in DelaunayTriangulation.jl.

  • conditions::Union{<:Tuple,<:ConditionType}

The classification for the boundary condition type corresponding to each boundary index as above. See ConditionType for possible conditions - should be one of Neumann, Dudt, Dirichlet, or Constrained.

Keyword Arguments

  • parameters=ntuple(_ -> nothing, length(functions))

The parameters for the functions, with parameters[i] giving the argument p in functions[i].

Outputs

The returned value is the corresponding BoundaryConditions struct.

source

There are four types of boundary conditions: Neumann, Dudt, Dirichlet, and Constrained. These types are defined below.

FiniteVolumeMethod.Neumann — Constant
Neumann

Instance of a ConditionType used for declaring that an edge has a Neumann condition. Neumann conditions take the form

\[\vb q(x, y, t) \vdot \vu n_\sigma = a(x, y, t, u)\]

where $\vb q$ is the flux function and $\vu n_\sigma$ is the outward unit normal vector field on the associated edge (meaning, for example, the normal vector to an edge pq would point to the right of pq).

When providing a Neumann condition, the function you provide takes the form

a(x, y, t, u, p)

where (x, y) is the point, t is the current time, and u is the solution at the point (x, y) at time t, as above, with an extra argument p for additional parameters.

source
FiniteVolumeMethod.Dudt — Constant
Dudt

Instance of a ConditionType used for declaring that an edge, or a point, has a Dudt-type boundary condition. Dudt-type conditions take the form

\[\dv{u(x, y, t)}{t} = a(x, y, t, u).\]

When providing a Dudt condition, the function you provide takes the form

a(x, y, t, u, p)

where (x, y) is the point, t is the current time, and u is the solution at the point (x, y) at time t, as above, with an extra argument p for additional parameters.

source
FiniteVolumeMethod.Dirichlet — Constant
Dirichlet

Instance of a ConditionType used for declaring that an edge, or a point, has a Dirichlet boundary condition. Dirichlet conditions take the form

\[u(x, y, t) = a(x, y, t, u).\]

When providing a Dirichlet condition, the function you provide takes the form

a(x, y, t, u, p)

where (x, y) is the point, t is the current time, and u is the solution at the point (x, y) at time t, as above, with an extra argument p for additional parameters.

source
FiniteVolumeMethod.Constrained — Constant
Constrained

Instance of a ConditionType used for declaring that an edge has a Constrained condition. When an edge has this condition associated with it, it will be treated as any normal edge and no boundary condition will be applied. With this condition, it is assumed that you will later setup your problem as a differential-algebraic equation (DAE) and provide the appropriate constraints. See the docs for some examples.

When you provide a Constrained condition, for certain technical reasons you do still need to provide a function that corresponds to it in the function list provided to BoundaryConditions. For this function, any function will work, e.g. sin - it will not be called. The proper constraint function is to be provided after-the-fact, as explained in the docs.

source

InternalConditions: Defining internal conditions

If you like, you can also put some constraints for nodes away from the boundary. In this case, only Dudt and Dirichlet conditions can be imposed; for Neumann or Constrained conditions, you need to consider differential-algebraic problems as considered in the tutorials. The docstring for InternalConditions is below; the fields of InternalConditions are not public API, only this wrapper is.

FiniteVolumeMethod.InternalConditions — Type
InternalConditions(functions=();
    dirichlet_nodes::Dict{Int,Int}=Dict{Int,Int}(),
    dudt_nodes::Dict{Int,Int}=Dict{Int,Int}(),
    parameters::Tuple=ntuple(_ -> nothing, length(functions)))

This is a constructor for the InternalConditions struct, which holds the internal conditions for the PDE. See also Conditions (which FVMProblem wraps this into), ConditionType, and BoundaryConditions.

Arguments

  • functions::Union{<:Tuple,<:Function}=()

The functions that define the internal conditions. These are the functions referred to in edge_conditions and point_conditions.

Keyword Arguments

  • dirichlet_nodes::Dict{Int,Int}=Dict{Int,Int}()

A Dict that stores all Dirichlet points u as keys, with keys mapping to indices idx that refer to the corresponding condition function and parameters in functions and parameters.

  • dudt_nodes::Dict{Int,Int}=Dict{Int,Int}()

A Dict that stores all Dudt points u as keys, with keys mapping to indices idx that refer to the corresponding condition function and parameters in functions and parameters.

  • parameters::Tuple=ntuple(_ -> nothing, length(functions))

The parameters for the functions, with parameters[i] giving the argument p in functions[i].

Outputs

The returned value is the corresponding InternalConditions struct.

Note

When the internal conditions get merged with the boundary conditions, any internal conditions that are placed onto the boundary will be replaced with the boundary condition at that point on the boundary.

source

FVMProblem: Defining the PDE

Once you have defined the mesh, the boundary conditions, and possibly the internal conditions, you can now construct the PDE itself. This is done using FVMProblem, whose docstring is below; the fields of FVMProblem are public API.

FiniteVolumeMethod.FVMProblem — Type
FVMProblem(mesh, boundary_conditions[, internal_conditions];
    diffusion_function=nothing,
    diffusion_parameters=nothing,
    source_function=nothing,
    source_parameters=nothing,
    flux_function=nothing,
    flux_parameters=nothing,
    initial_condition,
    initial_time=0.0,
    final_time)

Constructs an FVMProblem. See also FVMSystem and SteadyFVMProblem.

Arguments

  • mesh::FVMGeometry

The mesh on which the PDE is defined, given as a FVMGeometry.

  • boundary_conditions::BoundaryConditions

The boundary conditions for the PDE, given as a BoundaryConditions.

  • internal_conditions::InternalConditions=InternalConditions()

The internal conditions for the PDE, given as an InternalConditions. This argument is optional.

Keyword Arguments

  • diffusion_function=nothing

If isnothing(flux_function), then this can be provided to give the diffusion-source formulation. See also construct_flux_function. Should be of the form D(x, y, t, u, p).

  • diffusion_parameters=nothing

The argument p for diffusion_function.

  • source_function=(x, y, t, u, p) -> zero(u)

The source term, given in the form S(x, y, t, u, p).

  • source_parameters=nothing

The argument p for source_function.

  • flux_function=nothing

The flux function, given in the form q(x, y, t, α, β, γ, p) ↦ (qx, qy), where (qx, qy) is the flux vector, (α, β, γ) are the shape function coefficients for estimating u ≈ αx+βy+γ. If isnothing(flux_function), then diffusion_function is used instead to construct the function.

  • flux_parameters=nothing

The argument p for flux_function.

  • initial_condition

The initial condition, with initial_condition[i] the initial value at the ith node of the mesh.

  • initial_time=0.0

The initial time.

  • final_time

The final time.

Outputs

The returned value is the corresponding FVMProblem struct. You can then solve the problem using solve(::Union{FVMProblem,FVMSystem}, ::Any; kwargs...) from DifferentialEquations.jl.

source

For this problem, you can provide either a diffusion_function or a flux_function. In the former case, the flux_function is constructed from diffusion_function using construct_flux_function, whose docstring is shown below; construct_flux_function is not public API.

FiniteVolumeMethod.construct_flux_function — Function
construct_flux_function(q, D, Dp)

If isnothing(q), then this returns the flux function based on the diffusion function D and diffusion parameters Dp, so that the new function is 

(x, y, t, α, β, γ, _) -> -D(x, y, t, α*x + β*y + γ, Dp)[α, β]

Otherwise, just returns q again.

source

Additionally, FVMProblem merges the provided boundary conditions and internal conditions into a Conditions type, defined below; the documented fields of Conditions are public API.

FiniteVolumeMethod.Conditions — Type
Conditions{F} <: AbstractConditions

This is a struct that holds the boundary and internal conditions for the PDE. The constructor is 

Conditions(mesh::FVMGeometry, bc::BoundaryConditions, ic::InternalConditions=InternalConditions())

The fields are:

Fields

  • neumann_conditions::Dict{NTuple{2,Int},Int}

A Dict that stores all Neumann edges (u, v) as keys, with keys mapping to indices idx that refer to the corresponding condition function and parameters in functions.

  • constrained_conditions::Dict{NTuple{2,Int},Int}

A Dict that stores all Constrained edges (u, v) as keys, with keys mapping to indices idx that refer to the corresponding condition function and parameters in functions.

  • dirichlet_conditions::Dict{Int,Int}

A Dict that stores all Dirichlet points u as keys, with keys mapping to indices idx that refer to the corresponding condition function and parameters in functions.

  • dudt_conditions::Dict{Int,Int}

A Dict that stores all Dudt points u as keys, with keys mapping to indices idx that refer to the corresponding condition function and parameters in functions.

  • functions::F<:Tuple

The functions that define the boundary and internal conditions. The ith function should correspond to the part of the boundary of the mesh corresponding to the ith boundary index, as defined in DelaunayTriangulation.jl. The ith function is stored as a ParametrisedFunction.

source

Note that the functions in Conditions get wrapped into a ParametrisedFunction type, which is public API.

FiniteVolumeMethod.ParametrisedFunction — Type
ParametrisedFunction{F<:Function,P} <: Function

This is a struct that wraps a function f and some parameters p into a single object.

Fields

  • fnc::F

The function that is wrapped.

  • parameters::P

The parameters that are wrapped.

source

SteadyFVMProblem: Making the problem a steady-state problem

To make an FVMProblem a steady-state problem, meaning that you are solving

\[\div \vb q(\vb x, t, u) = S(\vb x, t, u)\]

rather than

\[\pdv{u(\vb x, t)}{t} + \div \vb q(\vb x, t, u) = S(\vb x, t, u),\]

than you need to wrap the FVMProblem inside a SteadyFVMProblem, defined below; the fields of SteadyFVMProblem are public API.

FVMSystem: Defining a system of PDEs

We also allow for systems of PDEs to be defined, where this system should take the form

\[\begin{equation*} \begin{aligned} \pdv{u_1(\vb x, t)}{t} + \div \vb q_1(\vb x, t, u_1, \ldots, u_n) &= S_1(\vb x, t, u_1, \ldots, u_n), \\ \pdv{u_2(\vb x, t)}{t} + \div \vb q_2(\vb x, t, u_1, \ldots, u_n) &= S_2(\vb x, t, u_1, \ldots, u_n), \\ &\vdots \\ \pdv{u_n(\vb x, t)}{t} + \div \vb q_n(\vb x, t, u_1, \ldots, u_n) &= S_n(\vb x, t, u_1, \ldots, u_n). \end{aligned} \end{equation*}\]

To define this system, you need to provide an FVMProblem for each equation, and then construct an FVMSystem from these problems. The docstring for FVMSystem is below; the fields of FVMSystem are public API.

FiniteVolumeMethod.FVMSystem — Type
FVMSystem(prob1, prob2, ..., probN)

Constructs a representation for a system of PDEs, where each probi is a FVMProblem for the ith component of the system.

For these FVMProblems, the functions involved, such as the condition functions, should all be defined so that the u argument assumes the form u = (u₁, u₂, ..., uN) (both Tuples and Vectors will be passed), where uᵢ is the solution for the ith component of the system. For the flux functions, which for a FVMProblem takes the form

q(x, y, t, α, β, γ, p) ↦ (qx, qy),

the same form is used, except α, β, γ are all Tuples so that α[i]*x + β[i]*y + γ[i] is the approximation to uᵢ.

Providing default flux functions

Due to this difference in flux functions, and the need to provide α, β, and γ to the flux function, for FVMSystems you need to provide a flux function rather than a diffusion function. If you do provide a diffusion function, it will error when you try to solve the problem.

This problem is solved in the same way as a FVMProblem, except the problem is defined such that the solution returns a matrix at each time, where the (j, i)th component corresponds to the solution at the ith node for the jth component.

See also FVMProblem and SteadyFVMProblem.

Note

To construct a steady-state problem for an FVMSystem, you need to apply SteadyFVMProblem to the system rather than first applying it to each individual FVMProblem in the system.

source

If you want to make a steady-state problem for an FVMSystem, you should apply SteadyFVMProblem to FVMSystem rather than to each FVMProblem individually.

solve: Solving the PDE

You can use solve from the SciMLBase ecosystem to solve these PDEs. This allows you to use any of the available algorithms from DifferentialEquations.jl for solving these problems. For non-steady problems, the relevant function is (which is public API)

CommonSolve.solve — Method
solve(prob::Union{FVMProblem,FVMSystem}, alg; 
    specialization=SciMLBase.AutoSpecialize, 
    jac_prototype=jacobian_sparsity(prob), 
    parallel::Val{<:Bool}=Val(true),
    kwargs...)

Solves the given FVMProblem or FVMSystem prob with the algorithm alg.

Missing vertices

When the underlying triangulation, tri, has points in get_points(tri) that are not vertices of the triangulation itself, the associated values of the solution at these points will not be updated by the solver, and will remain at their initial values.

Arguments

  • prob: The problem to be solved.
  • alg: The algorithm to be used to solve the problem. This can be any of the algorithms in DifferentialEquations.jl.

Keyword Arguments

  • specialization=SciMLBase.AutoSpecialize: The type of specialization to be used. See https://docs.sciml.ai/DiffEqDocs/stable/features/low_dep/#Controlling-Function-Specialization-and-Precompilation.
  • jac_prototype=jacobian_sparsity(prob): The prototype for the Jacobian matrix, constructed by default from jacobian_sparsity.
  • parallel::Val{<:Bool}=Val(true): Whether to use multithreading. Use Val(false) to disable multithreading.
  • kwargs...: Any other keyword arguments to be passed to the solver.

Outputs

The returned value sol depends on the type of the problem.

In this case, sol::ODESolution is such that the ith component of sol refers to the ith node of the underlying mesh.

In this case, the (j, i)th component of sol::ODESolution refers to the ith node of the underlying mesh for the jth component of the system.

source

For steady-state problems, the algorithms to use are those from NonlinearSolve.jl. The relevant function is still solve and is public API:

CommonSolve.solve — Method
solve(prob::SteadyFVMProblem, alg; 
    specialization=SciMLBase.AutoSpecialize, 
    jac_prototype=jacobian_sparsity(prob),
    parallel::Val{<:Bool}=Val(true),
    kwargs...)

Solves the given FVMProblem or FVMSystem prob with the algorithm alg.

Arguments

  • prob: The problem to be solved.
  • alg: The algorithm to be used to solve the problem. This can be any of the algorithms in NonlinearSolve.jl.

Keyword Arguments

  • specialization=SciMLBase.AutoSpecialize: The type of specialization to be used. See https://docs.sciml.ai/DiffEqDocs/stable/features/low_dep/#Controlling-Function-Specialization-and-Precompilation.
  • jac_prototype=jacobian_sparsity(prob): The prototype for the Jacobian matrix, constructed by default from jacobian_sparsity.
  • parallel::Val{<:Bool}=Val(true): Whether to use multithreading. Use Val(false) to disable multithreading.
  • kwargs...: Any other keyword arguments to be passed to the solver.

Outputs

The returned value sol depends on whether the underlying problem is a FVMProblem or an FVMSystem, but in each case it is an ODESolution type that can be accessed like the solutions in DifferentialEquations.jl:

In this case, sol is such that the ith component of sol refers to the ith node of the underlying mesh.

In this case, the (j, i)th component of sol refers to the ith node of the underlying mesh for the jth component of the system.

source

These solve functions rely on fvm_eqs! for evaluating the equations. You should never need to use fvm_eqs! directly, unless you are using a differential-algebraic equation. The docstring for fvm_eqs! is below; this function is public API.

FiniteVolumeMethod.fvm_eqs! — Function
fvm_eqs!(du, u, p, t)

Computes the finite volume method equations for the current time t and solution u. This function is public API.

The argument p depends on whether the problem is being solved in parallel or not. If it is solved serially, than the fields are:

  • p.prob: The prob <: AbstractFVMProblem.
  • p.parallel: Val(false).

If the problem is solved in parallel, then the fields are:

  • p.prob: The prob <: AbstractFVMProblem.
  • p.parallel: Val(true).
  • p.duplicated_du: A Matrix of the same size as du that is used to store the contributions to du from each thread.
  • p.solid_triangles: A Vector of the solid triangles in the triangulation.
  • p.solid_vertices: A Vector of the solid vertices in the triangulation.
  • p.chunked_solid_triangles: A Vector of tuples of the form (range, chunk_idx) where range is a range of indices into p.solid_triangles and chunk_idx is the index of the chunk.
  • p.boundary_edges: A Vector of the boundary edges in the triangulation.
  • p.chunked_boundary_edges: A Vector of tuples of the form (range, chunk_idx) where range is a range of indices into p.boundary_edges and chunk_idx is the index of the chunk.

These fields are public API, although note that they are not intended to be modified by the user, and we may freely add in new fields over new versions.

source

Custom constraints

You can also provide custom constraints. Rather than outlining this precisely here, it is best explained in the tutorials, namely this tutorial. We note that one useful function for this is compute_flux, which allows you to compute the flux across a given edge. The docstring for compute_flux is below, and this function is public API.

FiniteVolumeMethod.compute_flux — Function
compute_flux(prob::AbstractFVMProblem, i, j, u, t)

Given an edge with indices (i, j), a solution vector u, a current time t, and a problem prob, computes the flux ∇⋅q(x, y, t, α, β, γ, p) ⋅ n, where n is the normal vector to the edge, q is the flux function from prob, (x, y) is the midpoint of the edge, (α, β, γ) are the shape function coefficients, and p are the flux parameters from prob. If prob is an FVMSystem, the returned value is a Tuple for each individual flux. The normal vector n is a clockwise rotation of the edge, meaning pointing right of the edge (i, j).

source

Piecewise linear interpolation

You can evaluate the piecewise linear interpolation corresponding to a solution using pl_interpolate, defined below; this function is public API.

FiniteVolumeMethod.pl_interpolate — Function
pl_interpolate(prob, T, u, x, y)

Given a prob <: AbstractFVMProblem, a triangle T containing a point (x, y), and a set of function values u at the corresponding nodes of prob, interpolates the solution at the point (x, y) using piecewise linear interpolation.

source

Better interpolants are available from NaturalNeighbours.jl - see the this tutorial for some examples.