Solvers for Specific Problems, and Writing Your Own
The problems solved by this package are quite general, taking the form
\[\pdv{u}{t} + \div\vb q = S.\]
For some problems, though, this is not the most efficient form to implement. For example, the diffusion equation
\[\pdv{u}{t} = D\grad^2 u\]
might be better treated by converting the problem into
\[\dv{\vb u}{t} = \vb A\vb u + \vb b,\]
which is faster to solve than if we were to treat it as a nonlinear problem (which is done by default). For this reason, we define some templates for specific types of problems, namely:
DiffusionEquation
s: $\partial_tu = \div[D(\vb x)\grad u]$.MeanExitTimeProblem
s: $\div[D(\vb x)\grad T(\vb x)] = -1$.LinearReactionDiffusionEquation
s: $\partial_tu = \div[D(\vb x)\grad u] + f(\vb x)u$.PoissonsEquation
: $\div[D(\vb x)\grad u] = f(\vb x)$.LaplacesEquation
: $\div[D(\vb x)\grad u] = 0$.
The docstrings below define the templates for these problems.
FiniteVolumeMethod.AbstractFVMTemplate
— Typeabstract type AbstractFVMTemplate <: AbstractFVMProblem
An abstract type that defines some specific problems. These problems are those that could be defined directly using FVMProblem
s, but are common enough that (1) it is useful to have them defined here, and (2) it is useful to have them defined in a way that is more efficient than with a default implementation (e.g. exploiting linearity). The problems are all defined as subtypes of a common abstract type, namely, AbstractFVMTemplate
(the home of this docstring), which itself is a subtype of AbstractFVMProblem
.
To understand how to make use of these specific problems, either see the docstring for each problem, or see the "Solvers for Specific Problems, and Writing Your Own" section of the docs.
To see the full list of problems, do
julia> using FiniteVolumeMethod
julia> subtypes(FiniteVolumeMethod.AbstractFVMTemplate)
5-element Vector{Any}:
DiffusionEquation
LaplacesEquation
LinearReactionDiffusionEquation
MeanExitTimeProblem
PoissonsEquation
The constructor for each problem is defined in its docstring. Note that all the problems above are exported.
These problems can all be solved using the standard solve
interface from DifferentialEquations.jl, just like for FVMProblem
s. The only exception is for steady state problems, in which case the solve
interface is still used, except the interface is from LinearSolve.jl.
CommonSolve.solve
— Methodsolve(prob::AbstractFVMTemplate, args...; kwargs...)
Solve the problem prob
using the standard solve
interface from DifferentialEquations.jl. For steady state problems, the interface is from LinearSolve.jl.
FiniteVolumeMethod.DiffusionEquation
— TypeDiffusionEquation <: AbstractFVMTemplate
A struct for defining a problem representing a diffusion equation:
\[\pdv{u}{t} = \div\left[D(\vb x)\grad u\right]\]
inside a domain $\Omega$.
You can solve this problem using solve
.
The solution to this problem will have an extra component added to it. The original solution will be inside sol[begin:end-1, :]
, where sol
is the solution returned by solve
.
Constructor
DiffusionEquation(mesh::FVMGeometry,
BCs::BoundaryConditions,
ICs::InternalConditions=InternalConditions();
diffusion_function,
diffusion_parameters=nothing,
initial_condition,
initial_time=0.0,
final_time,
kwargs...)
Arguments
mesh::FVMGeometry
: TheFVMGeometry
.BCs::BoundaryConditions
: TheBoundaryConditions
. For these boundary conditions, all functions should still be of the form(x, y, t, u, p) -> Number
, but thet
andu
arguments should be unused as they will be replaced withnothing
.ICs::InternalConditions=InternalConditions()
: TheInternalConditions
. For these internal conditions, all functions should still be of the form(x, y, t, u, p) -> Number
, but thet
andu
arguments should be unused as they will be replaced withnothing
.
Keyword Arguments
diffusion_function
: The diffusion function. Should be of the form(x, y, p) -> Number
, wherep = diffusion_parameters
below.diffusion_parameters=nothing
: The argumentp
indiffusion_function
.initial_condition
: The initial condition.initial_time=0.0
: The initial time.final_time
: The final time.kwargs...
: Any other keyword arguments are passed to theODEProblem
(from DifferentialEquations.jl) that represents the problem.
Fields
The struct has extra fields in addition to the arguments above:
A
: This is a sparse matrixA
so thatdu/dt = Au + b
.b
: Theb
above.Aop
: TheMatrixOperator
that represents the system so thatdu/dt = Aop*u
(withu
padded with an extra component sinceA
is now insideAop
).problem
: TheODEProblem
that represents the problem. This is the problem that is solved when you callsolve
on the struct.
FiniteVolumeMethod.MeanExitTimeProblem
— TypeMeanExitTimeProblem
A struct for defining a problem representing a mean exit time problem:
\[\div\left[D(\vb x)\grad T\right] =-1\]
inside a domain $\Omega$. This problem is a special case of PoissonsEquation
, but is defined separately since it is common enough to warrant its own definition; MeanExitTimeProblem
is constructed using PoissonsEquation
.
You can solve this problem using solve
.
Constructor
MeanExitTimeProblem(mesh::FVMGeometry,
BCs::BoundaryConditions,
ICs::InternalConditions=InternalConditions();
diffusion_function,
diffusion_parameters=nothing,
kwargs...)
Arguments
mesh::FVMGeometry
: TheFVMGeometry
.BCs::BoundaryConditions
: TheBoundaryConditions
.ICs::InternalConditions=InternalConditions()
: TheInternalConditions
.
The functions for BCs
and ICs
are not used. Whenever a Neumann
condition is encountered, or a Dirichlet
condition, it is assumed that the conditon is homogeneous. If any of the conditions are Dudt
or Constrained
types, then an error is thrown.
Keyword Arguments
diffusion_function
: The diffusion function. Should be of the form(x, y, p) -> Number
, wherep = diffusion_parameters
below.diffusion_parameters=nothing
: The argumentp
indiffusion_function
.kwargs...
: Any other keyword arguments are passed to theLinearProblem
(from LinearSolve.jl) that represents the problem.
Fields
The struct has extra fields in addition to the arguments above:
A
: This is a sparse matrixA
so thatAT = b
.b
: Theb
above.problem
: TheLinearProblem
that represents the problem. This is the problem that is solved when you callsolve
on the struct.
FiniteVolumeMethod.LinearReactionDiffusionEquation
— TypeLinearReactionDiffusionEquation
A struct for defining a problem representing a linear reaction-diffusion equation:
\[\pdv{u}{t} = \div\left[D(\vb x)\grad u\right] + f(\vb x)u\]
inside a domain $\Omega$.
You can solve this problem using solve
.
The solution to this problem will have an extra component added to it. The original solution will be inside sol[begin:end-1, :]
, where sol
is the solution returned by solve
.
Constructor
LinearReactionDiffusionEquation(mesh::FVMGeometry,
BCs::BoundaryConditions,
ICs::InternalConditions=InternalConditions();
diffusion_function,
diffusion_parameters=nothing,
source_function,
source_parameters=nothing,
initial_condition,
initial_time=0.0,
final_time,
kwargs...)
Arguments
mesh::FVMGeometry
: TheFVMGeometry
.BCs::BoundaryConditions
: TheBoundaryConditions
. For these boundary conditions, all functions should still be of the form(x, y, t, u, p) -> Number
, but thet
andu
arguments should be unused as they will be replaced withnothing
.ICs::InternalConditions=InternalConditions()
: TheInternalConditions
. For these internal conditions, all functions should still be of the form(x, y, t, u, p) -> Number
, but thet
andu
arguments should be unused as they will be replaced withnothing
.
Keyword Arguments
diffusion_function
: The diffusion function. Should be of the form(x, y, p) -> Number
, wherep = diffusion_parameters
below.diffusion_parameters=nothing
: The argumentp
indiffusion_function
.source_function
: The source function. Should be of the form(x, y, p) -> Number
, wherep = source_parameters
below.source_parameters=nothing
: The argumentp
insource_function
.initial_condition
: The initial condition.initial_time=0.0
: The initial time.final_time
: The final time.kwargs...
: Any other keyword arguments are passed to theODEProblem
(from DifferentialEquations.jl) that represents the problem.
Fields
The struct has extra fields in addition to the arguments above:
A
: This is a sparse matrixA
so thatdu/dt = Au + b
.b
: Theb
above.Aop
: TheMatrixOperator
that represents the system so thatdu/dt = Aop*u
(withu
padded with an extra component sinceA
is now insideAop
).problem
: TheODEProblem
that represents the problem. This is the problem that is solved when you callsolve
on the struct.
FiniteVolumeMethod.PoissonsEquation
— TypePoissonsEquation
A struct for defining a problem representing a (generalised) Poisson's equation:
\[\div[D(\vb x)\grad u] = f(\vb x)\]
inside a domain $\Omega$. See also LaplacesEquation
, a special case of this problem with $f(\vb x) = 0$.
You can solve this problem using solve
.
Constructor
PoissonsEquation(mesh::FVMGeometry,
BCs::BoundaryConditions,
ICs::InternalConditions=InternalConditions();
diffusion_function=(x,y,p)->1.0,
diffusion_parameters=nothing,
source_function,
source_parameters=nothing,
kwargs...)
Arguments
mesh::FVMGeometry
: TheFVMGeometry
.BCs::BoundaryConditions
: TheBoundaryConditions
. For these boundary conditions, all functions should still be of the form(x, y, t, u, p) -> Number
, but thet
andu
arguments should be unused as they will be replaced withnothing
.ICs::InternalConditions=InternalConditions()
: TheInternalConditions
. For these internal conditions, all functions should still be of the form(x, y, t, u, p) -> Number
, but thet
andu
arguments should be unused as they will be replaced withnothing
.
Keyword Arguments
diffusion_function=(x,y,p)->1.0
: The diffusion function. Should be of the form(x, y, p) -> Number
, wherep = diffusion_parameters
below.diffusion_parameters=nothing
: The argumentp
indiffusion_function
.source_function
: The source function. Should be of the form(x, y, p) -> Number
, wherep = source_parameters
below.source_parameters=nothing
: The argumentp
insource_function
.kwargs...
: Any other keyword arguments are passed to theLinearProblem
(from LinearSolve.jl) that represents the problem.
Fields
The struct has extra fields in addition to the arguments above:
A
: This is a sparse matrixA
so thatAu = b
.b
: Theb
above.problem
: TheLinearProblem
that represents the problem. This is the problem that is solved when you callsolve
on the struct.
FiniteVolumeMethod.LaplacesEquation
— TypeLaplacesEquation
A struct for defining a problem representing a (generalised) Laplace's equation:
\[\div[D(\vb x)\grad u] = 0\]
inside a domain $\Omega$. See also PoissonsEquation
.
You can solve this problem using solve
.
Constructor
LaplacesEquation(mesh::FVMGeometry,
BCs::BoundaryConditions,
ICs::InternalConditions=InternalConditions();
diffusion_function=(x,y,p)->1.0,
diffusion_parameters=nothing,
kwargs...)
Arguments
mesh::FVMGeometry
: TheFVMGeometry
.BCs::BoundaryConditions
: TheBoundaryConditions
. For these boundary conditions, all functions should still be of the form(x, y, t, u, p) -> Number
, but thet
andu
arguments should be unused as they will be replaced withnothing
.ICs::InternalConditions=InternalConditions()
: TheInternalConditions
. For these internal conditions, all functions should still be of the form(x, y, t, u, p) -> Number
, but thet
andu
arguments should be unused as they will be replaced withnothing
.
Keyword Arguments
diffusion_function=(x,y,p)->1.0
: The diffusion function. Should be of the form(x, y, p) -> Number
, wherep = diffusion_parameters
below.diffusion_parameters=nothing
: The argumentp
indiffusion_function
.kwargs...
: Any other keyword arguments are passed to theLinearProblem
(from LinearSolve.jl) that represents the problem.
Fields
The struct has extra fields in addition to the arguments above:
A
: This is a sparse matrixA
so thatAu = b
.b
: Theb
above.problem
: TheLinearProblem
that represents the problem. This is the problem that is solved when you callsolve
on the struct.
Now, again, we note that all these problems can already be implemented using the main interface FVMProblem
. However, the templates we provide are more efficient, and also provide a good starting point for writing your own solver, meaning your own function that evaluates the system of ODEs. In the sections that follow, we will demonstrate two things for each of the problems above:
- The mathematical details involved in implementing each template.
- Examples of using the templates from FiniteVolumeMethod.jl.
With these two steps, you should be able to also know how to write your own solver for any problem you like.
Relevant Docstrings for Writing Your Own Solver
For writing these solvers, there are some specific functions that might be of use to you. Here, we provide the docstrings for these functions. These functions are public API.
FiniteVolumeMethod.get_dirichlet_fidx
— Functionget_dirichlet_fidx(conds, node)
Get the index of the function that corresponds to the Dirichlet
condition at node
.
FiniteVolumeMethod.is_dirichlet_node
— Functionis_dirichlet_node(conds, node)
Check if node
has a Dirichlet
condition.
FiniteVolumeMethod.get_dirichlet_nodes
— Functionget_dirichlet_nodes(conds)
Get all nodes that have a Dirichlet
condition.
FiniteVolumeMethod.has_dirichlet_nodes
— Functionhas_dirichlet_nodes(conds)
Check if any node has a Dirichlet
condition.
FiniteVolumeMethod.get_dudt_fidx
— Functionget_dudt_fidx(conds, node)
Get the index of the function that corresponds to the Dudt
condition at node
.
FiniteVolumeMethod.is_dudt_node
— Functionis_dudt_node(conds, node)
Check if node
has a Dudt
condition.
FiniteVolumeMethod.get_dudt_nodes
— Functionget_dudt_nodes(conds)
Get all nodes that have a Dudt
condition.
FiniteVolumeMethod.has_dudt_nodes
— Functionhas_dudt_nodes(conds)
Check if any node has a Dudt
condition.
FiniteVolumeMethod.get_neumann_fidx
— Functionget_neumann_fidx(conds, i, j)
Get the index of the function that corresponds to the Neumann
condition at the edge (i, j)
.
FiniteVolumeMethod.is_neumann_edge
— Functionis_neumann_edge(conds, i, j)
Check if the edge (i, j)
has a Neumann
condition.
FiniteVolumeMethod.has_neumann_edges
— Functionhas_neumann_edges(conds)
Check if any edge has a Neumann
condition.
FiniteVolumeMethod.get_neumann_edges
— Functionget_neumann_edges(conds)
Get all edges that have a Neumann
condition.
FiniteVolumeMethod.get_constrained_fidx
— Functionget_constrained_fidx(conds, i, j)
Get the index of the function that corresponds to the Constrained
condition at the edge (i, j)
.
FiniteVolumeMethod.is_constrained_edge
— Functionis_constrained_edge(conds, i, j)
Check if the edge (i, j)
has a Constrained
condition.
FiniteVolumeMethod.has_constrained_edges
— Functionhas_constrained_edges(conds)
Check if any edge has a Constrained
condition.
FiniteVolumeMethod.get_constrained_edges
— Functionget_constrained_edges(conds)
Get all edges that have a Constrained
condition.
FiniteVolumeMethod.eval_condition_fnc
— Functioneval_condition_fnc(conds, fidx, x, y, t, u)
Evaluate the function that corresponds to the condition at fidx
at the point (x, y)
at time t
with solution u
.
FiniteVolumeMethod.has_condition
— Functionhas_condition(conds, node)
Check if node
has any condition.
FiniteVolumeMethod.get_cv_components
— Functionget_cv_components(props, edge_index)
Get the quantities for a control volume edge interior to the associated triangulation, relative to the edge_index
th edge of the triangle corresponding to props
.
Outputs
x
: Thex
-coordinate of the edge's midpoint.y
: They
-coordinate of the edge's midpoint.nx
: Thex
-component of the edge's normal vector.ny
: They
-component of the edge's normal vector.ℓ
: The length of the edge.
FiniteVolumeMethod.get_boundary_cv_components
— Functionget_boundary_cv_components(tri::Triangulation, i, j)
Get the quantities for both control volume edges lying a boundary edge (i, j)
.
Outputs
nx
: Thex
-component of the edge's normal vector.ny
: They
-component of the edge's normal vector.mᵢx
: Thex
-coordinate of the midpoint of thei
th vertex and the edge's midpoint.mᵢy
: They
-coordinate of the midpoint of thei
th vertex and the edge's midpoint.mⱼx
: Thex
-coordinate of the midpoint of thej
th vertex and the edge's midpoint.mⱼy
: They
-coordinate of the midpoint of thej
th vertex and the edge's midpoint.ℓᵢ
: Half the length of the boundary edge, which is the length of the control volume edge.T
: The triangle containing the boundary edge.props
: TheTriangleProperties
forT
.
FiniteVolumeMethod.get_triangle_props
— Functionget_triangle_props(mesh, i, j, k)
Get the TriangleProperties
for the triangle (i, j, k)
in mesh
.
FiniteVolumeMethod.get_volume
— Functionget_volume(mesh, i)
Get the volume of the i
th control volume in mesh
.
DelaunayTriangulation.get_point
— Methodget_point(mesh, i)
Get the i
th point in mesh
.
FiniteVolumeMethod.triangle_contributions!
— Functiontriangle_contributions!(A, mesh, conditions, diffusion_function, diffusion_parameters)
Add the contributions from each triangle to the matrix A
, based on the equation
\[\dv{u_i}{t} = \frac{1}{V_i}\sum_{\sigma \in \mathcal E_i} D(\vb x_\sigma)\left[\left(s_{k, 11}n_\sigma^x + s_{k, 21}n_\sigma^y\right)u_{k1} + \left(s_{k, 12}n_\sigma^x + s_{k, 22}n_\sigma^y\right)u_{k2} + \left(s_{k, 13}n_\sigma^x + s_{k, 23}n_\sigma^y\right)u_{k3}\right]L_\sigma + S_i, \]
as explained in the docs. Will not update any rows corresponding to Dirichlet
or Dudt
nodes.
FiniteVolumeMethod.apply_dirichlet_conditions!
— Functionapply_dirichlet_conditions!(initial_condition, mesh, conditions)
Applies the Dirichlet conditions specified in conditions
to the initial_condition
. The boundary conditions are assumed to take the form a(x, y, t, u, p) -> Number
, but t
and u
are passed as nothing
. Note that this assumes that the associated system (A, b)
is such that A[i, :]
is all zero, and b[i]
is zero, where i
is a node with a Dirichlet condition.
FiniteVolumeMethod.apply_dudt_conditions!
— Functionapply_dudt_conditions!(b, mesh, conditions)
Applies the Dudt conditions specified in conditions
to the b
vector. The boundary conditions are assumed to take the form a(x, y, t, u, p) -> Number
, but t
and u
are passed as nothing
. Note that this assumes that the associated system (A, b)
is such that A[i, :]
is all zero, so that replacing b[i]
with the boundary condition will set duᵢ/dt = b[i]
.
FiniteVolumeMethod.boundary_edge_contributions!
— Functionboundary_edge_contributions!(A, b, mesh, conditions, diffusion_function, diffusion_parameters)
Add the contributions from each boundary edge to the matrix A
, based on the equation
\[\dv{u_i}{t} = \frac{1}{V_i}\sum_{\sigma \in \mathcal E_i} D(\vb x_\sigma)\left[\left(s_{k, 11}n_\sigma^x + s_{k, 21}n_\sigma^y\right)u_{k1} + \left(s_{k, 12}n_\sigma^x + s_{k, 22}n_\sigma^y\right)u_{k2} + \left(s_{k, 13}n_\sigma^x + s_{k, 23}n_\sigma^y\right)u_{k3}\right]L_\sigma + S_i, \]
as explained in the docs. Will not update any rows corresponding to Dirichlet
or Dudt
nodes.
FiniteVolumeMethod.non_neumann_boundary_edge_contributions!
— Functionnon_neumann_boundary_edge_contributions!(A, mesh, conditions, diffusion_function, diffusion_parameters)
Add the contributions from each non-Neumann boundary edge to the matrix A
, based on the equation
\[\dv{u_i}{t} = \frac{1}{V_i}\sum_{\sigma \in \mathcal E_i} D(\vb x_\sigma)\left[\left(s_{k, 11}n_\sigma^x + s_{k, 21}n_\sigma^y\right)u_{k1} + \left(s_{k, 12}n_\sigma^x + s_{k, 22}n_\sigma^y\right)u_{k2} + \left(s_{k, 13}n_\sigma^x + s_{k, 23}n_\sigma^y\right)u_{k3}\right]L_\sigma + S_i, \]
as explained in the docs. Will not update any rows corresponding to Dirichlet
or Dudt
nodes.
FiniteVolumeMethod.neumann_boundary_edge_contributions!
— Functionneumann_boundary_edge_contributions!(b, mesh, conditions, diffusion_function, diffusion_parameters)
Add the contributions from each Neumann boundary edge to the vector b
, based on the equation
\[\dv{u_i}{t} = \frac{1}{V_i}\sum_{\sigma \in \mathcal E_i} D(\vb x_\sigma)\left[\grad u(\vb x_\sigma) \vdot \vu n\right]L_\sigma + S_i,\]
as explained in the docs. Will not update any rows corresponding to Dirichlet
or Dudt
nodes. This function will pass nothing
in place of the arguments u
and t
in the boundary condition functions.
neumann_boundary_edge_contributions!(F, mesh, conditions, diffusion_function, diffusion_parameters, u, t)
Add the contributions from each Neumann boundary edge to the vector F
, based on the equation
\[\dv{u_i}{t} = \frac{1}{V_i}\sum_{\sigma \in \mathcal E_i} D(\vb x_\sigma)\left[\grad u(\vb x_\sigma) \vdot \vu n\right]L_\sigma + S_i,\]
as explained in the docs. Will not update any rows corresponding to Dirichlet
or Dudt
nodes.
FiniteVolumeMethod.create_rhs_b
— Functioncreate_rhs_b(mesh, conditions, source_function, source_parameters)
Create the vector b
defined by
b = [source_function(x, y, source_parameters) for (x, y) in DelaunayTriangulation.each_point(mesh.triangulation)],
and b[i] = 0
whenever i
is a Dirichlet node.
FiniteVolumeMethod.apply_steady_dirichlet_conditions!
— Functionapply_steady_dirichlet_conditions!(A, b, mesh, conditions)
Applies the Dirichlet conditions specified in conditions
to the initial_condition
. The boundary conditions are assumed to take the form a(x, y, t, u, p) -> Number
, but t
and u
are passed as nothing
. Note that this assumes that the associated system (A, b)
is such that A[i, :]
is all zero, and b[i]
is zero, where i
is a node with a Dirichlet condition.
For a steady problem Au = b
, applies the Dirichlet boundary conditions specified by conditions
so that A[i, i] = 1
and b[i]
is the condition, where i
is a boundary node. Note that this assumes that all of A[i, :]
is zero before setting A[i, i] = 1
.
FiniteVolumeMethod.two_point_interpolant
— Functiontwo_point_interpolant(mesh, u, i, j, mx, my)
Given a mesh <: FVMGeometry
, a set of function values u
at the nodes of mesh
, and a point (mx, my)
on the line segment between the nodes i
and j
, interpolates the solution at the point (mx, my)
using two-point interpolation.
FiniteVolumeMethod.get_dirichlet_callback
— Functionget_dirichlet_callback(prob[, f=update_dirichlet_nodes!]; kwargs...)
Get the callback for updating Dirichlet
nodes. The kwargs...
argument is ignored, except to detect if a user has already provided a callback, in which case the callback gets merged into a CallbackSet
with the Dirichlet
callback. If the problem prob
has no Dirichlet
nodes, the returned callback does nothing and is never called.
You can provide f
to change the function that updates the Dirichlet
nodes.
FiniteVolumeMethod.jacobian_sparsity
— Functionjacobian_sparsity(prob)
Returns a prototype for the Jacobian of the given prob
.
FiniteVolumeMethod.fix_missing_vertices!
— Functionfix_missing_vertices!(A, b, mesh)
Given a system (A, b)
and a mesh
, sets A[i, i] = 1
and b[i] = 0
for any vertices i
that are a point in mesh
, but not an actual vertex in the mesh
.