Interface
- Interface
FVMGeometry
: Defining the meshBoundaryConditions
: Defining boundary conditionsInternalConditions
: Defining internal conditionsFVMProblem
: Defining the PDESteadyFVMProblem
: Making the problem a steady-state problemFVMSystem
: Defining a system of PDEssolve
: Solving the PDE- Custom constraints
- Piecewise linear interpolation
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:
- Use
FVMGeometry
to define the problem's mesh. - Provide boundary conditions using
BoundaryConditions
. - (Optional) Provide internal conditions using
InternalConditions
. - Convert the problem into an
FVMProblem
. - If you want to make the problem steady, use
SteadyFVMProblem
on theFVMProblem
. - If you want a system of equations, construct an
FVMSystem
from multipleFVMProblem
s; if you want this problem to be steady, skip step 5 and only now applySteadyFVMProblem
. - Solve the problem using
solve
. - For a discussion of custom constraints, see the tutorials.
- 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
— TypeFVMGeometry(tri::Triangulation)
This is a constructor for the FVMGeometry
struct, which holds the mesh and associated data for the PDE.
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 underlyingTriangulation
from DelaunayTriangulation.jl.triangulation_statistics
: The statistics of the triangulation.cv_volumes::Vector{Float64}
: AVector
of the volumes of each control volume.triangle_props::Dict{NTuple{3,Int},TriangleProperties}
: ADict
mapping the indices of each triangle to its [TriangleProperties
].
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
— TypeTriangleProperties(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 incv_edge_midpoints
.cv_edge_lengths::NTuple{3,Float64}
: The lengths of the control volume edges, in the same order as incv_edge_midpoints
.
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 n
th vertex, respectively, α = s₁u₁ + s₂u₂ + s₃u₃
, β = s₄u₁ + s₅u₂ + s₆u₃
, and γ = s₇u₁ + s₈u₂ + s₉u₃
.
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
— TypeBoundaryConditions(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 i
th function should correspond to the part of the boundary of the mesh
corresponding to the i
th 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.
There are four types of boundary conditions: Neumann
, Dudt
, Dirichlet
, and Constrained
. These types are defined below.
FiniteVolumeMethod.ConditionType
— TypeConditionType
This is an Enum
-type, with four instances:
This is used for declaring conditions in the PDEs. See the associated docstrings, and also BoundaryConditions
and InternalConditions
.
FiniteVolumeMethod.Neumann
— ConstantNeumann
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.
FiniteVolumeMethod.Dudt
— ConstantDudt
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.
FiniteVolumeMethod.Dirichlet
— ConstantDirichlet
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.
FiniteVolumeMethod.Constrained
— ConstantConstrained
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.
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
— TypeInternalConditions(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.
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.
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
— TypeFVMProblem(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 i
th 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.
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
— Functionconstruct_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.
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
— TypeConditions{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 i
th function should correspond to the part of the boundary of the mesh
corresponding to the i
th boundary index, as defined in DelaunayTriangulation.jl. The i
th function is stored as a ParametrisedFunction
.
Note that the functions in Conditions
get wrapped into a ParametrisedFunction
type, which is public API.
FiniteVolumeMethod.ParametrisedFunction
— TypeParametrisedFunction{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.
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.
FiniteVolumeMethod.SteadyFVMProblem
— TypeSteadyFVMProblem(prob::AbstractFVMProblem)
This is a wrapper for an AbstractFVMProblem
that indicates that the problem is to be solved as a steady-state problem. You can then solve the problem using solve(::SteadyFVMProblem, ::Any; kwargs...)
from NonlinearSolve.jl. Note that you need to have set the final time to Inf
if you want a steady state out at infinity rather than some finite actual time.
See also FVMProblem
and FVMSystem
.
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
— TypeFVMSystem(prob1, prob2, ..., probN)
Constructs a representation for a system of PDEs, where each probi
is a FVMProblem
for the i
th component of the system.
For these FVMProblem
s, 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 Tuple
s and Vector
s will be passed), where uᵢ
is the solution for the i
th 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 Tuple
s so that α[i]*x + β[i]*y + γ[i]
is the approximation to uᵢ
.
Due to this difference in flux functions, and the need to provide α
, β
, and γ
to the flux function, for FVMSystem
s 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 i
th node for the j
th component.
See also FVMProblem
and SteadyFVMProblem
.
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.
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
— Methodsolve(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
.
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 fromjacobian_sparsity
.parallel::Val{<:Bool}=Val(true)
: Whether to use multithreading. UseVal(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 i
th component of sol
refers to the i
th node of the underlying mesh.
In this case, the (j, i)
th component of sol::ODESolution
refers to the i
th node of the underlying mesh for the j
th component of the system.
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
— Methodsolve(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 fromjacobian_sparsity
.parallel::Val{<:Bool}=Val(true)
: Whether to use multithreading. UseVal(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 i
th component of sol
refers to the i
th node of the underlying mesh.
In this case, the (j, i)
th component of sol
refers to the i
th node of the underlying mesh for the j
th component of the system.
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!
— Functionfvm_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
: Theprob <: AbstractFVMProblem
.p.parallel
:Val(false)
.
If the problem is solved in parallel, then the fields are:
p.prob
: Theprob <: AbstractFVMProblem
.p.parallel
:Val(true)
.p.duplicated_du
: AMatrix
of the same size asdu
that is used to store the contributions todu
from each thread.p.solid_triangles
: AVector
of the solid triangles in the triangulation.p.solid_vertices
: AVector
of the solid vertices in the triangulation.p.chunked_solid_triangles
: AVector
of tuples of the form(range, chunk_idx)
whererange
is a range of indices intop.solid_triangles
andchunk_idx
is the index of the chunk.p.boundary_edges
: AVector
of the boundary edges in the triangulation.p.chunked_boundary_edges
: AVector
of tuples of the form(range, chunk_idx)
whererange
is a range of indices intop.boundary_edges
andchunk_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.
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
— Functioncompute_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)
.
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
— Functionpl_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.
Better interpolants are available from NaturalNeighbours.jl - see the this tutorial for some examples.