System

The computational grid required is assumed to correspond to a domain $\Omega=\cup_{r=1}^{n_\Omega} \Omega_r$

Grids for VoronoiFVM are managed by the packages ExtendableGrids.jl and SimplexGridFactory.jl

with boundary $\partial\Omega=\Gamma=\cup_{b=1}^{n_\Gamma} \Gamma_b$.

The subdomains $\Omega_r$ are called "regions" and the boundary subdomains $\Gamma_b$ are called "boundary regions".

On this complex of domains "lives" a number of species which are either attached to a number of regions or to a number of boundary regions.

All these data, the matrix for the linear system and other things are hold together by a struct VoronoiFVM.System. This type is not exported to avoid name clashes.

System constructors

VoronoiFVM.SystemMethod
System(grid; kwargs...)

Create structure of type VoronoiFVM.System{Tv,Ti, Tm, TSpecMat<:AbstractMatrix, TSolArray<:AbstractMatrix} holding data for finite volume system solution.

Parameters:

  • grid::ExtendableGrid: 1, 2 or 3D computational grid

Keyword arguments:

  • species: vector of integer species indices. Added to all grid regions, avoiding the need to call enable_species! for this default case. If it is kept empty, species have be added to the system after creation via enable_species!.
  • unknown_storage: string or symbol. Information on species distribution is kept in sparse or dense matrices matrices and, correspondingly, the solution array is of type SparseSolutionArray or matrix, respectively. In the case of sparse unknown storage, the system matrix handles exactly those degrees of freedom which correspond to unknowns. However, handling of the sparse matrix structures for the bookkeeping of the unknowns creates overhead.
    • :dense : solution vector is an nspecies x nnodes dense matrix
    • :sparse : solution vector is an nspecies x nnodes sparse matrix
  • matrixindextype: Integer type. Index type for sparse matrices created in the system.
  • is_linear: whether the system is linear or not. If it is linear, only one Newton step is used to solve it.
  • assembly: either :cellwise (default) or :edgewise. Determine, how the assembly loop is organized. :cellwise means that the outer loop goes over grid cells (triangles, tetrahedra), and contributions to edge fluxes and node reactions are calculated for each cell. As a consequence, e.g. im 2D for all interior edges, flux functions are callled twice, once for each adjacent cell. Especially in 3D, this becomes a significant overhead. With :edgewise, geometry factors of these edges are pre-assembled, and the outer assembly loops go over all grid edges resp. nodes, still with separate calls if neigboring cells belong to different regions.
Note

It is planned to make :edgewise the default in a later version.

Physics keyword arguments:

  • flux: Function. Flux between neighboring control volumes: flux(f,u,edge) or flux(f,u,edge,data) should return in f[i] the flux of species i along the edge joining circumcenters of neighboring control volumes. For species i,u[i,1] and u[i,2] contain the unknown values at the corresponding ends of the edge.

  • storage: Function. Storage term (term under time derivative): storage(f,u,node) or storage(f,u,node,data) It should return in f[i] the storage term for the i-th equation. u[i] contains the value of the i-th unknown.

  • reaction: Function. Reaction term: reaction(f,u,node) or reaction(f,u,node,data) It should return in f[i] the reaction term for the i-th equation. u[i] contains the value of the i-th unknown.

  • edgereaction: Function. Edge reeaction term: edgereaction(f,u,edge) or edgereaction(f,u,edge,data) It should return in f[i] the reaction term for the i-th equation. For species i,u[i,1] and u[i,2] contain the unknown values at the corresponding ends of the edge.

  • source: Function. Source term: source(f,node) or source(f,node,data). It should return the in f[i] the value of the source term for the i-th equation.

  • bflux: Function. Flux between neighboring control volumes on the boundary

  • breaction Function. Boundary reaction term: breaction(f,u,node) or breaction(f,u,node,data) Similar to reaction, but restricted to the inner or outer boundaries.

  • bcondition Function. Alias for breaction.

  • bsource: Function. Boundary source term: bsource(f,node) or bsource(f,node,data). It should return in f[i] the value of the source term for the i-th equation.

  • bstorage: Function. Boundary storage term: bstorage(f,u,node) or bstorage(f,u,node,data) Similar to storage, but restricted to the inner or outer boundaries.

  • generic_operator: Function. Generic operator generic_operator(f,u,sys). This operator acts on the full solution u of a system. Sparsity is detected automatically unless generic_operator_sparsity is given.

  • generic_operator_sparsity: Function defining the sparsity structure of the generic operator. This should return the sparsity pattern of the generic_operator.

  • nparams: number of parameters the system is depending on, and with respect to which the derivatives need to be obtained

  • data: User data (parameters). This allows to pass various parameters to the callback functions. If data is given, all callback functions should accept a last data argument. Otherwise, no data are passed explicitly, and constitutive callbacks can take parameters from the closure where the function is defined.

  • matrixtype: :default, :sparse, :tridiagonal, :banded

source

Adding species by species numbers

VoronoiFVM.enable_species!Method
enable_species!(system,ispec,regions)

Add species ispec to a list of bulk regions. Species numbers for bulk and boundary species have to be distinct. Once a species has been added, it cannot be removed.

source
VoronoiFVM.enable_species!Method
enable_species!(system; kwargs...)

Keyword arguments:

  • species: Integer or vector of integers. Species to be added to the system.
  • regions: Vector of integers. Regions, where these species shall be added.If nothing, they are added to all species.

Once a species has been added, it cannot be removed.

source
VoronoiFVM.enable_boundary_species!Function
enable_boundary_species!(system,ispec,regions)

Add species ispec to a list of boundary regions. Species numbers for bulk and boundary species have to be distinct. Once a species has been added, it cannot be removed.

source

Handling boundary conditions

Boundary conditions are handled in the bcondition callback passed to the system constructor. For being called in this callback, the following functions are available

VoronoiFVM.boundary_dirichlet!Method
 boundary_dirichlet!(y,u,bnode, args...; kwargs...)

Keyword argument version:

  • species: species number. Default: 1
  • region: boundary region number. By default, all boundary regions.
  • value: value
source
VoronoiFVM.boundary_neumann!Method
 boundary_neumann!(y,u,bnode, args...; kwargs...)

Keyword argument version:

  • species: species number. Default: 1
  • region: boundary region number. By default, all boundary regions.
  • value: value
source
VoronoiFVM.boundary_robin!Method
 boundary_robin!(y,u,bnode, args...; kwargs...)

Keyword argument version:

  • species: species number. Default: 1
  • region: boundary region number. By default, all boundary regions.
  • factor: factor
  • value: value
source
VoronoiFVM.rampFunction
   ramp(t; kwargs...)

Ramp function for specifying time dependent boundary conditions

Keyword arguments:

  • dt: Tuple: start and end time of ramp. Default: (0,0.1)
  • du: Tuple: values at start and end time. Default: (0,0)
source

Outflow boundary conditions

These are characterized by the boutflow physics callback and and the outflowboundaries keyword argument in the system resp. physics constructor. See also the corresponding notebook

VoronoiFVM.hasoutflownodeFunction
hasoutflownode(edge)

Check if one node of the edge is situated on a boundary region listed in outflowboundaries, see [struct Physics].

source
VoronoiFVM.isoutflownodeFunction
isoutflownode(edge,inode)

Check if inode (1 or 2) is an outflow node.

source
isoutflownode(edge,inode,irefgion)

Check if inode (1 or 2) is an outflow node on boundary region iregion.

source

Allocation warnings

The code checks for allocations in the assembly loop. Care has been taken to ensure that allocations in the assembly loop don't emerge from VoronoiFVM.jl code.

If allocations occur in the assembly loop, they happen in the physics callbacks. The corresponding warnings can bee switched off by passing a verbosity strings without 'a' to the solver. If no data are allocated in the physics callbacks, these allocations are probably due to type instabilities in physics callbacks, see the the discussion here. Type instabilities can be debugged via the @time macro applied to expressions in a physics callback.

The following cases provide some ideas where to look for reasons of the problem and possible remedies:

Case 1: a parameter changes its value, and Julia is not sure about the type.

eps=1.0

flux(f,u,edge)
    f[1]=eps*(u[1,1]-[1,2])
end
... solve etc ...
eps=2.0

Remedy: use a type annotation eps::Float64=... to signalize your intent to Julia. This behaviour is explained in the Julia documentation.

Case 2: variables in the closure have the same name as a variable introduced in a callback.

flux(f,u,edge)
    x=(u[1,1]-[1,2])
    f[1]=x
end

... create etc ...

x=solve(...)

Remedy: rename e.g. x=solve() to sol=solve()

Various tools

VoronoiFVM.num_speciesFunction
num_species(edge::VoronoiFVM.AbstractEdge) -> Any

Return number of species for edge

source
num_species(system)

Number of species in system

source
num_species(a)

Number of species (size of first dimension) of solution array.

source
VoronoiFVM.unknownsMethod
unknowns(system; inival, inifunc)

Create a solution vector for system. If inival is not specified, the entries of the returned vector are undefined.

source
unknowns(impedance_system)

Create a vector of unknowns of the impedance system

source
VoronoiFVM.unknownsMethod
unknowns(Tu, system; inival, inifunc)

Create a solution vector for system with elements of type Tu. If inival is not specified, the entries of the returned vector are undefined.

source
Base.mapFunction
map(inifunc, sys)

Create a solution vector for system using the callback inifunc which has the same signature as a source term.

source
map(inival, sys)

Create a solution vector for system using a constant initial value

source
Base.map!Function
map!(inifunc, U, system)

Map inifunc onto solution array U

source
Base.reshapeMethod
reshape(v, system)

Reshape vector to fit as solution to system.

source

Types

VoronoiFVM.AbstractSystemType
abstract type AbstractSystem{Tv<:Number, Tc<:Number, Ti<:Integer, Tm<:Integer}

Abstract type for finite volume system structure.

source
VoronoiFVM.SystemType
mutable struct System{Tv, Tc, Ti, Tm, TSpecMat<:(AbstractMatrix), TSolArray<:(AbstractMatrix)} <: VoronoiFVM.AbstractSystem{Tv, Tc, Ti, Tm}

Structure holding data for finite volume system.

source

Legacy API

VoronoiFVM.boundary_dirichlet!Method
boundary_dirichlet!(system, ispec, ibc, v)

Set Dirichlet boundary condition for species ispec at boundary ibc:

$u_{ispec}=v$ on $\Gamma_{ibc}$

Info

Starting with version 0.14, it is preferable to define boundary conditions within the bcondition physics callback

source
VoronoiFVM.boundary_dirichlet!Method
  boundary_dirichlet!(system; kwargs...)

Keyword argument version:

  • species: species number
  • region: region number
  • value: value
Info

Starting with version 0.14, it is preferable to define boundary conditions within the bcondition physics callback

source
VoronoiFVM.boundary_neumann!Method
boundary_neumann!(system, ispec, ibc, v)

Set Neumann boundary condition for species ispec at boundary ibc:

$\mathrm{flux}_{ispec}\cdot \vec n=v$ on $\Gamma_{ibc}$

Info

Starting with version 0.14, it is preferable to define boundary conditions within the bcondition physics callback

source
VoronoiFVM.boundary_neumann!Method
  boundary_neumann!(system; kwargs...)

Keyword argument version:

  • species: species number
  • region: region number
  • value: value
Info

Starting with version 0.14, it is preferable to define boundary conditions within the bcondition physics callback

source
VoronoiFVM.boundary_robin!Method
boundary_robin!(system, ispec, ibc, α, v)

Set Robin boundary condition for species ispec at boundary ibc:

$\mathrm{flux}_{ispec}\cdot \vec n + \alpha u_{ispec}=v$ on $\Gamma_{ibc}$

Info

Starting with version 0.14, it is preferable to define boundary conditions within the bcondition physics callback

source
VoronoiFVM.boundary_robin!Method
  boundary_robin!(system; kwargs...)

Keyword argument version:

  • species: species number
  • region: region number
  • factor: factor
  • value: value
Info

Starting with version 0.14, it is preferable to define boundary conditions within the bcondition physics callback

source
VoronoiFVM.viewKFunction
viewK(
    edge::VoronoiFVM.AbstractEdge,
    u
) -> VoronoiFVM.VectorUnknowns

Solution view on first edge node

source
VoronoiFVM.viewLFunction
viewL(
    edge::VoronoiFVM.AbstractEdge,
    u
) -> VoronoiFVM.VectorUnknowns

Solution view on second edge node

source