GitHub

Solvers

The package comes with a built-in solve method which solves stationary problems, homotopy embedding problems and transient problems via the implicit Euler method. In particular, the transient solver allows to use nonlinear storage terms.

Alternatively, OrdinaryDiffEq.jl based solvers can be used for transient problems.

Built-in solver

This solver and its default parameters are tuned for robustness, possibly at the expense of solution speed. Careful tuning of the parameters, or – in the case of transient problems – the choice of a OrdinaryDiffEq.jl based solver can significantly improve the performance.

Overview:

Solve method

CommonSolve.solve — Method
solve(system; kwargs...)

Built-in solution method for VoronoiFVM.System.

Keyword arguments:

  • General for all solvers

    • inival (default: 0) : Array created via unknowns or number giving the initial value.
    • control (default: nothing): Pass instance of SolverControl
    • All elements of SolverControl can be used as kwargs. Eventually overwrites values given via control
    • params: Parameters (Parameter handling is experimental and may change)

  • Stationary solver: Invoked if neither times nor embed, nor tstep are given as keyword argument.

    • time (default: 0.0): Set time value.

    Returns a DenseSolutionArray or SparseSolutionArray

  • Embedding (homotopy) solver: Invoked if embed kwarg is given. Use homotopy embedding + damped Newton's method to solve a stationary problem or to solve series of parameter dependent problems. Parameter step control is performed according to solver control data. Keyword arguments and default values are:

    • embed (default: nothing ): vector of parameter values to be reached exactly

    In addition, all kwargs of the implicit Euler solver (besides times) are handled, including handle_exceptions (see below). Step size control is performed based on the parameters Δp, Δp_max, Δp_min, Δp_grow, Δp_decrease, Δu_opt, Δu_max_factor. Returns a transient solution object sol containing the stored solution(s), see TransientSolution.

  • Implicit Euler transient solver: Invoked if times kwarg is given. Use implicit Euler method + damped Newton's method to solve time dependent problem. Time step control is performed according to solver control data. kwargs and default values are:

    • times (default: nothing ): vector of time values to be reached exactly
    • pre (default: (sol,t)->nothing ): callback invoked before each time step
    • post (default: (sol,oldsol, t, Δt)->nothing ): callback invoked after each time step
    • sample (default: (sol,t)->nothing ): callback invoked after timestep for all times in times[2:end].
    • delta (default: (system, u,v,t, Δt)->norm(sys,u-v,Inf) ): Value used to control the time step size Δu

    Step size control is performed based on the parameters Δt, Δt_max, Δt_min, Δt_grow, Δt_decrease, Δu_opt, Δu_max_factor. If control.handle_exceptions is true, if time step solution throws an error, stepsize is lowered, and step solution is called again with a smaller time value. If control.Δt<control.Δt_min, solution is aborted with error. Returns a transient solution object sol containing the stored solution, see TransientSolution.

  • Implicit Euler timestep solver. Invoked if tstep kwarg is given. Solve one time step of the implicit Euler method.

    • time (default: 0): Set time value.
    • tstep: time step

    Returns a DenseSolutionArray or SparseSolutionArray

source

Solver control

VoronoiFVM.SolverControl — Type
SolverControl
SolverControl(;kwargs...)

Solver control parameter for time stepping, embedding, Newton method and linear solver control. All field names can be used as keyword arguments for solve(system::VoronoiFVM.AbstractSystem; kwargs...)

Newton's method solves $F(u)=0$ by the iterative procedure $u_{i+1}=u_{i} - d_i F'(u_i)^{-1}F(u_i)$ starting with some initial value $u_0$, where $d_i$ is a damping parameter.

  • verbose::Union{Bool, String}: Verbosity control. A collection of output categories is given in a string composed of the following letters:

    • a: allocation warnings
    • d: deprecation warnings
    • e: time/parameter evolution log
    • n: newton solver log
    • l: linear solver log

    Alternatively, a Bool value can be given, resulting in

    • true: "neda"
    • false: "da"

    Switch off all output including deprecation warnings via verbose="". In the output, corresponding messages are marked e.g. via '[n]', [a] etc. (besides of '[l]')

    For switching between printing output and logging output, see log_output! and print_output! .

  • abstol::Float64: Tolerance (in terms of norm of Newton update): terminate if $\Delta u_i=||u_{i+1}-u_i||_\infty <$ abstol.
  • reltol::Float64: Tolerance (relative to the size of the first update): terminate if $\Delta u_i/\Delta u_1<$ reltol.
  • maxiters::Int64: Maximum number of newton iterations.
  • tol_round::Float64: Tolerance for roundoff error detection: terminate if $|\;||u_{i+1}||_1 - ||u_{i}||_1\;|/ ||u_{i}||_1<$ tol_round occurred max_round times in a row.
  • tol_mono::Float64: Tolerance for monotonicity test: terminate with error if $\Delta u_i/\Delta u_{i-1}>$ 1/tol_mono.
  • damp_initial::Float64: Initial damping parameter $d_0$. To handle convergence problems, set this to a value less than 1.
  • damp_growth::Float64: Damping parameter growth factor: $d_{i+1}=\min(d_i\cdot$ max_growth $,1)$. It should be larger than 1.
  • max_round::Int64: Maximum number of consecutive iterations within roundoff error tolerance The default effectively disables this criterion.
  • unorm::Function: Calculation of Newton update norm
  • rnorm::Function: Functional for roundoff error calculation

  • method_linear::Union{Nothing, LinearSolve.SciMLLinearSolveAlgorithm}: Solver method for linear systems (see LinearSolve.jl). If given nothing, as default are chosen:

    • UMFPACKFactorization() for sparse matrices with Float64
    • SparspakFactorization() for sparse matrices with general number types.
    • Defaults from LinearSolve.jl for tridiagonal and banded matrices

    Users should experiment with what works best for their problem.

    For an overeview on possible alternatives, see VoronoiFVM.LinearSolverStrategy.

  • reltol_linear::Float64: Relative tolerance of iterative linear solver.
  • abstol_linear::Float64: Absolute tolerance of iterative linear solver.
  • maxiters_linear::Int64: Maximum number of iterations of linear solver

  • precon_linear::Union{Nothing, Function, ExtendableSparse.AbstractFactorization, LinearSolve.SciMLLinearSolveAlgorithm, Type}: Constructor for preconditioner for linear systems. This should work as a function precon_linear(A) which takes an AbstractSparseMatrixCSC (e.g. an ExtendableSparseMatrix) and returns a preconditioner object in the sense of LinearSolve.jl, i.e. which has an ldiv!(u,A,v) method. Useful examples:

    • ExtendableSparse.ILUZero
    • ExtendableSparse.Jacobi

    For easy access to this functionality, see see also VoronoiFVM.LinearSolverStrategy.

  • keepcurrent_linear::Bool: Update preconditioner in each Newton step ? Translates to reuse_precs=!keepcurrent_linear for LinearSolve.
  • Δp::Float64: Initial parameter step for embedding.
  • Δp_max::Float64: Maximal parameter step size.
  • Δp_min::Float64: Minimal parameter step size.
  • Δp_grow::Float64: Maximal parameter step size growth.
  • Δp_decrease::Float64: Parameter step decrease factor upon rejection
  • Δt::Float64: Initial time step size.
  • Δt_max::Float64: Maximal time step size.
  • Δt_min::Float64: Minimal time step size.
  • Δt_grow::Float64: Maximal time step size growth.
  • Δt_decrease::Float64: Time step decrease factor upon rejection
  • Δu_opt::Float64: Optimal size of update for time stepping and embedding. The algorithm tries to keep the difference in norm between "old" and "new" solutions approximately at this value.
  • Δu_max_factor::Float64: Control maximum sice of update Δu for time stepping and embedding relative to Δu_opt. Time steps with Δu > Δu_max_factor*Δu_opt will be rejected.
  • force_first_step::Bool: Force first timestep.
  • num_final_steps::Int64: Number of final steps to adjust at end of time interval in order to prevent breakdown of step size.

  • handle_exceptions::Bool: Handle exceptions during transient solver and parameter embedding. If true, exceptions in Newton solves are caught, the embedding resp. time step is lowered, and solution is retried. Moreover, if embedding or time stepping fails (e.g. due to reaching minimal step size), a warning is issued, and a solution is returned with all steps calculated so far.

    Otherwise (by default) errors are thrown.

  • store_all::Bool: Store all steps of transient/embedding problem:
  • in_memory::Bool: Store transient/embedding solution in memory
  • log::Any: Record history
  • edge_cutoff::Float64: Edge parameter cutoff for rectangular triangles.
  • pre::Function: Function pre(sol,t) called before time/embedding step
  • post::Function: Function post(sol,oldsol,t,Δt) called after successful time/embedding step
  • sample::Function: Function sample(sol,t) to be called for each t in times[2:end]
  • delta::Function: Time step error estimator. A function Δu=delta(system,u,uold,t,Δt) to calculate Δu.

  • tol_absolute::Union{Nothing, Float64}

  • tol_relative::Union{Nothing, Float64}

  • damp::Union{Nothing, Float64}

  • damp_grow::Union{Nothing, Float64}

  • max_iterations::Union{Nothing, Int64}

  • tol_linear::Union{Nothing, Float64}

  • max_lureuse::Union{Nothing, Int64}

  • mynorm::Union{Nothing, Function}

  • myrnorm::Union{Nothing, Function}

source
VoronoiFVM.fixed_timesteps! — Function
fixed_timesteps!(control,Δt; grow=1.0)

Modify control data such that the time steps are fixed to a geometric sequence such that Δtnew=Δtold*grow

source

System state

VoronoiFVM.SystemState — Type
mutable struct SystemState{Tv, Tp, TMatrix<:AbstractArray{Tv, 2}, TGenericMatrix, TSolArray<:AbstractArray{Tv, 2}, TData}

Structure holding state information for finite volume system.

Type parameters:

  • Tv: element type of solution vectors and matrix
  • TMatrix: matrix type
  • TSolArray: type of solution vector: (Matrix or SparseMatrixCSC)
  • TData: type of user data

Type fields:

  • system::VoronoiFVM.System: Related finite volume system
  • data::Any: User data
  • solution::AbstractMatrix: Solution vector
  • matrix::AbstractMatrix: Jacobi matrix for nonlinear problem
  • generic_matrix::Any: Sparse matrix for generic operator handling
  • dudp::Vector{TSolArray} where {Tv, TSolArray<:AbstractMatrix{Tv}}: Parameter derivative (vector of solution arrays)
  • update::AbstractMatrix: Vector holding Newton update
  • residual::AbstractMatrix: Vector holding Newton residual
  • linear_cache::Any: Linear solver cache
  • params::Vector: Parameter vector
  • uhash::UInt64: Hash value of latest unknowns vector the assembly was called with (used by differential equation interface)
  • history::Union{Nothing, VoronoiFVM.DiffEqHistory}: History record for solution process (used by differential equation interface)
source
VoronoiFVM.SystemState — Method
SystemState(Tv, system; data=system.physics.data, matrixtype=system.matrixtype)

Create state information for finite volume system.

Arguments:

  • Tv: value type of unknowns, matrix
  • system: Finite volume system

Keyword arguments:

  • data: User data. Default: data(system)
  • matrixtype. Default: system.matrixtype
source
VoronoiFVM.SystemState — Method
SystemState(Tv, system; data=system.physics.data, matrixtype=system.matrixtype)

Create state information for finite volume system.

Arguments:

  • Tv: value type of unknowns, matrix
  • system: Finite volume system

Keyword arguments:

  • data: User data. Default: data(system)
  • matrixtype. Default: system.matrixtype
source
SystemState(system; kwargs...)

Shortcut for creating state with value type defined by Tv type parameter of system

source
Base.similar — Method
similar(state; data=state.data)

Create a new state of with the same system, different work arrays, and possibly different data. The matrix of the new state initially shares the sparsity structure with state.

source

Linear solver control

Linear systems are solved using LinearSolve.jl. Linear solve compatible solver strategies (factorizations, iterative solvers) can be specified wia method_linear keyword argument to LinearSolve (equivalent to the method_linear entry of SolverControl.

Currently supported possibilities are documented in the documentation of ExtendableSparse.jl.

VoronoiFVM.jl provides partitioning methods for block preconditioners.

History handling

If log is set to true in solve, the history of newton iterations and time/embedding steps is recorded and returned as history(solution)

VoronoiFVM.NewtonSolverHistory — Type
mutable struct NewtonSolverHistory <: AbstractVector{Float64}

History information for one Newton solve of a nonlinear system. As an abstract vector it provides the history of the update norms. See summary and details for other ways to extract information.

  • nlu::Int64: number of Jacobi matrix factorizations

  • nlin::Int64: number of linear iteration steps / factorization solves

  • time::Float64: Elapsed time for solution

  • tasm::Float64: Elapsed time for assembly

  • tlinsolve::Float64: Elapsed time for linear solve

  • updatenorm::Any: History of norms of $||u_{i+1}-u_i||$

  • l1normdiff::Any: History of norms of $|\;||u_{i+1}||_1 - ||u_{i}||_1\;|/ ||u_{i}||_1$

source
VoronoiFVM.TransientSolverHistory — Type
mutable struct TransientSolverHistory <: AbstractVector{NewtonSolverHistory}

History information for transient solution/parameter embedding

As an abstract vector it provides the histories of each implicit Euler/embedding step. See summary and details for other ways to extract information.

  • histories::Any: Histories of each implicit Euler Newton iteration

  • times::Any: Time values

  • updates::Any: Update norms used for step control

source
VoronoiFVM.DiffEqHistory — Type
mutable struct DiffEqHistory

History information for DiffEqHistory

  • nd::Any: number of combined jacobi/rhs evaluations

  • njac::Any: number of combined jacobi evaluations

  • nf::Any: number of rhs evaluations

source
Base.summary — Method
summary(h::NewtonSolverHistory)

Return named tuple summarizing history.

source
Base.summary — Method
summary(h::TransientSolverHistory)

Return named tuple summarizing history.

source
VoronoiFVM.details — Function
details(h::NewtonSolverHistory)

Return array of named tuples with info on each iteration step

source
details(h::TransientSolverHistory)

Return array of details of each solver step

source

Matrix extraction

For testing and teaching purposes, one can obtain residual and linearization at a given vector of unknowns

VoronoiFVM.evaluate_residual_and_jacobian — Function
evaluate_residual_and_jacobian(system,u;
                               t=0.0, tstep=Inf, embed=0.0, init_dirichlet=false)

Evaluate residual and jacobian at solution value u. Returns a solution vector containing a copy of residual, and an ExendableSparseMatrix containing a copy of the linearization at u. The flag init_dirichlet controls whether u should be adjusted to satisfy the Dirichlet boundary conditions specified by the system.

source

OrdinaryDiffEq.jl transient solver

For transient problems, as an alternative to the built-in implicit Euler method, (stiff) ODE solvers from OrdinaryDiffEq.jl can be used.

The interface just provides two methods: creation of an ODEProblem from a VoronoiFVM.System and a reshape method which turns the output of the ode solver into a TransientSolution.

The basic usage pattern is as follows: use OrdinaryDiffEq.jl and replace the call to the built-in solver

sol=solve(sys; times=(t0,t1), inival=inival)

by

problem = ODEProblem(sys,inival,(t0,t1))
odesol = solve(problem, solver)
sol=reshape(odesol,sys)

Here, solver is some ODE/DAE solver from OrdinaryDiffEq.jl. It is preferable to choose methods able to handle stiff problems. Moreover, often, discretized PDE systems (e.g. containing elliptic equations) are differential agebraic equation (DAE) systems which should be solved by DAE solvers. Some choices to start with are Rosenbrock methods like Rosenbrock23 and multistep methods like QNDF and FBDF.

If the DifferentialEquations.jl package is loaded, the solver parameter can be omitted, and some default is chosen.

The solution odesol returned by solve conforms to the ArrayInterface but "forgot" the VoronoiFVM species structure. Using

Accessing odesol(t) will return an interpolated solution vector giving the value of the solution at moment t. Using reshape(::AbstractVector, ::VoronoiFVM.AbstractSystem) on (odesol(t),system) it can be turned into into a sparse or dense array reflecting the species structure of system. The order of the interpolation depends on the ODE solver.

Using reshape(::AbstractDiffEqArray,::VoronoiFVM.AbstractSystem) on (odesol, system) returns a TransientSolution knowing the species structure.

SciMLBase.ODEProblem — Type
ODEProblem(state,inival,tspan,callback=SciMLBase.CallbackSet())

Create an ODEProblem from a system state. See SciMLBase.ODEProblem(sys::VoronoiFVM.System, inival, tspan;kwargs...) for more documentation.

Defined in VoronoiFVM.jl.

source
ODEProblem(system,inival,tspan,callback=SciMLBase.CallbackSet())

Create an ODEProblem in mass matrix form which can be handled by ODE solvers from DifferentialEquations.jl.

Parameters:

  • system: A VoronoiFVM.System
  • inival: Initial value. Consider to pass a stationary solution at tspan[1].
  • tspan: Time interval
  • callback : (optional) callback for ODE solver

The method returns an ODEProblem which can be solved by solve().

Defined in VoronoiFVM.jl.

source
Base.reshape — Method
reshape(ode_solution, system; times=nothing, state=nothing)

Create a TransientSolution from the output of the ode solver which reflects the species structure of the system ignored by the ODE solver. Howvever the interpolation behind reshaped_sol(t) will be linear and ignores the possibility of higher order interpolations with ode_sol.

If times is specified, the (possibly higher order) interpolated solution at the given moments of time will be returned.

Defined in VoronoiFVM.jl.

source
SciMLBase.ODEFunction — Type
 ODEFunction(state; jacval = unknowns(sys, 0), tjac = 0 )

Create an ODEFunction. For more documentation, see SciMLBase.ODEFunction(state::VoronoiFVM.SystemState; kwargs...) jacval and tjac are passed to prepare_diffeq! and used there to calculate the Jacobian prototype.

Defined in VoronoiFVM.jl.

source
 ODEFunction(system; jacval=unknowns(system,inival=0),tjac=0)

Create an ODEFunction in mass matrix form to be handled by ODE solvers from DifferentialEquations.jl.

Parameters:

  • system: A VoronoiFVM.System
  • jacval (optional): Initial value. Default is a zero vector. Consider to pass a stationary solution at time tjac.
  • tjac (optional): tjac, Default: 0

The jacval and tjac are passed for a first evaluation of the Jacobian, allowing to detect the sparsity pattern which is passed to the solver.

Defined in VoronoiFVM.jl.

source

Legacy API

During the development of the code, a number of API variants have been developed which are supported for backward compatibility.

Deprecated API

The methods and struct in this section are deprecated as of version 2.4 and will be removed in version 3.

Linear solver strategies

VoronoiFVM.LinearSolverStrategy — Type
VoronoiFVM.LinearSolverStrategy

An linear solver strategy provides the possibility to construct SolverControl objects as follows:

    SolverControl(strategy,sys;kwargs...)

, e.g.

    SolverControl(GMRESIteration(UMFPackFactorization(), EquationBlock()),sys; kwargs...)

A linear solver strategy combines a Krylov method with a preconditioner which by default is calculated from the linearization of the initial value of the Newton iteration. For coupled systems, a blocking strategy can be chosen. The EquationBlock strategy calculates preconditioners or LU factorization separately for each species equation and combines them to a block Jacobi preconditioner. The PointBlock strategy treats the linear system as consisting of nspecies x nspecies blocks.

Which is the best strategy, depends on the space dimension. The following is a rule of thumb for choosing strategies

  • For 1D problems use direct solvers
  • For 2D stationary problems, use direct solvers, for transient problems consider iterative solvers which can take advantage of the diagonal dominance of the implicit timestep problem
  • For 3D problems avoid direct solvers

Currently available strategies are:

Notable LU Factorizations/direct solvers are:

Notable incomplete factorizations/preconditioners

Blocking strategies are:

source
VoronoiFVM.CGIteration — Type
CGIteration(;factorization=UMFPACKFactorization())
CGIteration(factorization)

CG Iteration from Krylov.jl via LinearSolve.jl.

source
VoronoiFVM.BICGstabIteration — Type
BICGstabIteration(;factorization=UMFPACKFactorization())
BICGstabIteration(factorization)

BICGstab Iteration from Krylov.jl via LinearSolve.jl.

source
VoronoiFVM.GMRESIteration — Type
GMRESIteration(;factorization=ILUZeroFactorization(), memory=20, restart=true)
GMRESIteration(factorization; memory=20, restart=true)

GMRES Iteration from Krylov.jl via LinearSolve.jl.

source

Block preconditioning

VoronoiFVM.EquationBlock — Type
EquationBlock()

Equation-wise blocking. Can be combined with any preconditioner resp. factorization including direct solvers.

source
VoronoiFVM.PointBlock — Type
PointBlock()

Point-wise blocking. Currently only together with ILUZeroFactorization. This requires a system with unknown_storage=:dense.

source