GitHub

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:

  1. DiffusionEquations: $\partial_tu = \div[D(\vb x)\grad u]$.
  2. MeanExitTimeProblems: $\div[D(\vb x)\grad T(\vb x)] = -1$.
  3. LinearReactionDiffusionEquations: $\partial_tu = \div[D(\vb x)\grad u] + f(\vb x)u$.
  4. PoissonsEquation: $\div[D(\vb x)\grad u] = f(\vb x)$.
  5. LaplacesEquation: $\div[D(\vb x)\grad u] = 0$.

The docstrings below define the templates for these problems.

FiniteVolumeMethod.AbstractFVMTemplate — Type
abstract type AbstractFVMTemplate <: AbstractFVMProblem

An 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
 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 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.

source
CommonSolve.solve — Method
solve(prob::AbstractFVMTemplate, args...; kwargs...)

Solve a template problem through its canonical SciML problem.

This is equivalent to solve(sciml_problem(prob; kwargs...), args...; kwargs...) and keeps the template families aligned with the same SciML-oriented execution contract used by the rest of the package.

source
FiniteVolumeMethod.DiffusionEquation — Type
DiffusionEquation <: 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.

Warning

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: The FVMGeometry.
  • BCs::BoundaryConditions: The BoundaryConditions. For these boundary conditions, all functions should still be of the form (x, y, t, u, p) -> Number, but the t and u arguments should be unused as they will be replaced with nothing.
  • ICs::InternalConditions=InternalConditions(): The InternalConditions. For these internal conditions, all functions should still be of the form (x, y, t, u, p) -> Number, but the t and u arguments should be unused as they will be replaced with nothing.

Keyword Arguments

  • diffusion_function: The diffusion function. Should be of the form (x, y, p) -> Number, where p = diffusion_parameters below.
  • diffusion_parameters=nothing: The argument p in diffusion_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 the ODEProblem (from DifferentialEquations.jl) that represents the problem.

Fields

The struct has extra fields in addition to the arguments above:

  • A: This is a sparse matrix A so that du/dt = Au + b.
  • b: The b above.
  • Aop: The MatrixOperator that represents the system so that du/dt = Aop*u (with u padded with an extra component since A is now inside Aop).
  • problem: The ODEProblem that represents the problem. This is the problem that is solved when you call solve on the struct.
source
FiniteVolumeMethod.MeanExitTimeProblem — Type
MeanExitTimeProblem

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

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, where p = diffusion_parameters below.
  • diffusion_parameters=nothing: The argument p in diffusion_function.
  • kwargs...: Any other keyword arguments are passed to the LinearProblem (from LinearSolve.jl) that represents the problem.

Fields

The struct has extra fields in addition to the arguments above:

  • A: This is a sparse matrix A so that AT = b.
  • b: The b above.
  • problem: The LinearProblem that represents the problem. This is the problem that is solved when you call solve on the struct.
source
FiniteVolumeMethod.LinearReactionDiffusionEquation — Type
LinearReactionDiffusionEquation

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.

Warning

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: The FVMGeometry.
  • BCs::BoundaryConditions: The BoundaryConditions. For these boundary conditions, all functions should still be of the form (x, y, t, u, p) -> Number, but the t and u arguments should be unused as they will be replaced with nothing.
  • ICs::InternalConditions=InternalConditions(): The InternalConditions. For these internal conditions, all functions should still be of the form (x, y, t, u, p) -> Number, but the t and u arguments should be unused as they will be replaced with nothing.

Keyword Arguments

  • diffusion_function: The diffusion function. Should be of the form (x, y, p) -> Number, where p = diffusion_parameters below.
  • diffusion_parameters=nothing: The argument p in diffusion_function.
  • source_function: The source function. Should be of the form (x, y, p) -> Number, where p = source_parameters below.
  • source_parameters=nothing: The argument p in source_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 the ODEProblem (from DifferentialEquations.jl) that represents the problem.

Fields

The struct has extra fields in addition to the arguments above:

  • A: This is a sparse matrix A so that du/dt = Au + b.
  • b: The b above.
  • Aop: The MatrixOperator that represents the system so that du/dt = Aop*u (with u padded with an extra component since A is now inside Aop).
  • problem: The ODEProblem that represents the problem. This is the problem that is solved when you call solve on the struct.
source
FiniteVolumeMethod.PoissonsEquation — Type
PoissonsEquation

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: The FVMGeometry.
  • BCs::BoundaryConditions: The BoundaryConditions. For these boundary conditions, all functions should still be of the form (x, y, t, u, p) -> Number, but the t and u arguments should be unused as they will be replaced with nothing.
  • ICs::InternalConditions=InternalConditions(): The InternalConditions. For these internal conditions, all functions should still be of the form (x, y, t, u, p) -> Number, but the t and u arguments should be unused as they will be replaced with nothing.

Keyword Arguments

  • diffusion_function=(x,y,p)->1.0: The diffusion function. Should be of the form (x, y, p) -> Number, where p = diffusion_parameters below.
  • diffusion_parameters=nothing: The argument p in diffusion_function.
  • source_function: The source function. Should be of the form (x, y, p) -> Number, where p = source_parameters below.
  • source_parameters=nothing: The argument p in source_function.
  • kwargs...: Any other keyword arguments are passed to the LinearProblem (from LinearSolve.jl) that represents the problem.

Fields

The struct has extra fields in addition to the arguments above:

  • A: This is a sparse matrix A so that Au = b.
  • b: The b above.
  • problem: The LinearProblem that represents the problem. This is the problem that is solved when you call solve on the struct.
source
FiniteVolumeMethod.LaplacesEquation — Type
LaplacesEquation

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: The FVMGeometry.
  • BCs::BoundaryConditions: The BoundaryConditions. For these boundary conditions, all functions should still be of the form (x, y, t, u, p) -> Number, but the t and u arguments should be unused as they will be replaced with nothing.
  • ICs::InternalConditions=InternalConditions(): The InternalConditions. For these internal conditions, all functions should still be of the form (x, y, t, u, p) -> Number, but the t and u arguments should be unused as they will be replaced with nothing.

Keyword Arguments

  • diffusion_function=(x,y,p)->1.0: The diffusion function. Should be of the form (x, y, p) -> Number, where p = diffusion_parameters below.
  • diffusion_parameters=nothing: The argument p in diffusion_function.
  • kwargs...: Any other keyword arguments are passed to the LinearProblem (from LinearSolve.jl) that represents the problem.

Fields

The struct has extra fields in addition to the arguments above:

  • A: This is a sparse matrix A so that Au = b.
  • b: The b above.
  • problem: The LinearProblem that represents the problem. This is the problem that is solved when you call solve on the struct.
source

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:

  1. The mathematical details involved in implementing each template.
  2. 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_cv_components — Function
get_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: The x-coordinate of the edge's midpoint.
  • y: The y-coordinate of the edge's midpoint.
  • nx: The x-component of the edge's normal vector.
  • ny: The y-component of the edge's normal vector.
  • ℓ: The length of the edge.
source
FiniteVolumeMethod.get_boundary_cv_components — Function
get_boundary_cv_components(tri::Triangulation, i, j)

Get the quantities for both control volume edges lying a boundary edge (i, j).

Outputs

  • nx: The x-component of the edge's normal vector.
  • ny: The y-component of the edge's normal vector.
  • mᵢx: The x-coordinate of the midpoint of the ith vertex and the edge's midpoint.
  • mᵢy: The y-coordinate of the midpoint of the ith vertex and the edge's midpoint.
  • mⱼx: The x-coordinate of the midpoint of the jth vertex and the edge's midpoint.
  • mⱼy: The y-coordinate of the midpoint of the jth 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: The TriangleProperties for T.
source
FiniteVolumeMethod.triangle_contributions! — Function
triangle_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.

source
FiniteVolumeMethod.apply_dirichlet_conditions! — Function
apply_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.

source
FiniteVolumeMethod.apply_dudt_conditions! — Function
apply_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].

source
FiniteVolumeMethod.boundary_edge_contributions! — Function
boundary_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.

source
FiniteVolumeMethod.non_neumann_boundary_edge_contributions! — Function
non_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.

source
FiniteVolumeMethod.neumann_boundary_edge_contributions! — Function
neumann_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.

source
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.

source
FiniteVolumeMethod.create_rhs_b — Function
create_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.

source
FiniteVolumeMethod.apply_steady_dirichlet_conditions! — Function
apply_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.

source
FiniteVolumeMethod.two_point_interpolant — Function
two_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.

source
FiniteVolumeMethod.get_dirichlet_callback — Function
get_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.

source
FiniteVolumeMethod.fix_missing_vertices! — Function
fix_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.

source