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:
DiffusionEquations: $\partial_tu = \div[D(\vb x)\grad u]$.MeanExitTimeProblems: $\div[D(\vb x)\grad T(\vb x)] = -1$.LinearReactionDiffusionEquations: $\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 <: AbstractFVMProblemAn abstract type that defines some specific problems. These problems are those that could be defined directly using FVMProblems, 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
PoissonsEquationThe 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 FVMProblems. 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 <: AbstractFVMTemplateA 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 thetanduarguments 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 thetanduarguments 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_parametersbelow.diffusion_parameters=nothing: The argumentpindiffusion_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 matrixAso thatdu/dt = Au + b.b: Thebabove.Aop: TheMatrixOperatorthat represents the system so thatdu/dt = Aop*u(withupadded with an extra component sinceAis now insideAop).problem: TheODEProblemthat represents the problem. This is the problem that is solved when you callsolveon the struct.
FiniteVolumeMethod.MeanExitTimeProblem — TypeMeanExitTimeProblemA 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_parametersbelow.diffusion_parameters=nothing: The argumentpindiffusion_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 matrixAso thatAT = b.b: Thebabove.problem: TheLinearProblemthat represents the problem. This is the problem that is solved when you callsolveon the struct.
FiniteVolumeMethod.LinearReactionDiffusionEquation — TypeLinearReactionDiffusionEquationA 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 thetanduarguments 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 thetanduarguments 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_parametersbelow.diffusion_parameters=nothing: The argumentpindiffusion_function.source_function: The source function. Should be of the form(x, y, p) -> Number, wherep = source_parametersbelow.source_parameters=nothing: The argumentpinsource_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 matrixAso thatdu/dt = Au + b.b: Thebabove.Aop: TheMatrixOperatorthat represents the system so thatdu/dt = Aop*u(withupadded with an extra component sinceAis now insideAop).problem: TheODEProblemthat represents the problem. This is the problem that is solved when you callsolveon the struct.
FiniteVolumeMethod.PoissonsEquation — TypePoissonsEquationA 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 thetanduarguments 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 thetanduarguments 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_parametersbelow.diffusion_parameters=nothing: The argumentpindiffusion_function.source_function: The source function. Should be of the form(x, y, p) -> Number, wherep = source_parametersbelow.source_parameters=nothing: The argumentpinsource_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 matrixAso thatAu = b.b: Thebabove.problem: TheLinearProblemthat represents the problem. This is the problem that is solved when you callsolveon the struct.
FiniteVolumeMethod.LaplacesEquation — TypeLaplacesEquationA 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 thetanduarguments 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 thetanduarguments 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_parametersbelow.diffusion_parameters=nothing: The argumentpindiffusion_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 matrixAso thatAu = b.b: Thebabove.problem: TheLinearProblemthat represents the problem. This is the problem that is solved when you callsolveon 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_indexth 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 theith vertex and the edge's midpoint.mᵢy: They-coordinate of the midpoint of theith vertex and the edge's midpoint.mⱼx: Thex-coordinate of the midpoint of thejth vertex and the edge's midpoint.mⱼy: They-coordinate of the midpoint of thejth 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: TheTrianglePropertiesforT.
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 ith control volume in mesh.
DelaunayTriangulation.get_point — Methodget_point(mesh, i)Get the ith 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.