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.solveMethod
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 stationary problem or to solve series of parameter dependent problems. Parameter step control is performed according to solver control data. kwargs 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. 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

    If control.handle_error 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.SolverControlType
SolverControl
SolverControl(;kwargs...)
SolverControl(linear_solver_strategy, sys; 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.

For linear solver strategies, see VoronoiFVM.LinearSolverStrategy.

  • 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]')
  • 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 (for Float64 calculations):

    • 1D: KLUFactorization()
    • 2D: SparspakFactorization()
    • 3D: UMFPACKFactorization()

    SparspakFactorization() is the default choice for general number types. Users should experiment with what works best for their problem.

    For easy access to this functionality, see see also 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 ?
  • Δ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 embeding 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

Linear solver strategies

VoronoiFVM.LinearSolverStrategyType
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.CGIterationType
CGIteration(;factorization=UMFPACKFactorization())
CGIteration(factorization)

CG Iteration from Krylov.jl via LinearSolve.jl.

source
VoronoiFVM.BICGstabIterationType
BICGstabIteration(;factorization=UMFPACKFactorization())
BICGstabIteration(factorization)

BICGstab Iteration from Krylov.jl via LinearSolve.jl.

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

GMRES Iteration from Krylov.jl via LinearSolve.jl.

source

Block preconditioning

This feature is under development as of 1.6.

VoronoiFVM.EquationBlockType
EquationBlock()

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

source
VoronoiFVM.PointBlockType
PointBlock()

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

source

History handling

If log is set to true in solve, the history of newton iterations and time/embedding steps is recorded and. For the respective previous solution step it can be obtained via history(system).

VoronoiFVM.NewtonSolverHistoryType
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.TransientSolverHistoryType
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
Base.summaryMethod
summary(h::NewtonSolverHistory)

Return named tuple summarizing history.

source
Base.summaryMethod
summary(h::TransientSolverHistory)

Return named tuple summarizing history.

source
VoronoiFVM.detailsFunction
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
VoronoiFVM.history_detailsFunction
history_details(tsol)

Return details of solver history from last solve call, if log was set to true. See details.

source
history_details(sys)

Return details of solver history from last solve call, if log was set to true. See details.

source
VoronoiFVM.history_summaryFunction
history_summary(tsol)

Return summary of solver history from last solve call, if log was set to true.

source
history_summary(sys)

Return summary of solver history from last solve call, if log was set to true.

source

Matrix extraction

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

VoronoiFVM.evaluate_residual_and_jacobianFunction
evaluate_residual_and_jacobian(system,u;
                               t=0.0, tstep=Inf,embed=0.0)

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.

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.

Base.reshapeMethod
reshape(ode_solution, system, times=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 ordee) interpolated solution at the given moments of time will be returned.

source
SciMLBase.ODEFunctionMethod
 ODEFunction(system,inival=unknowns(system,inival=0),t0=0)

Create an ODEPFunction in mass matrix form to be handeled 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.

source

Legacy API

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

VoronoiFVM.SolverStrategiesModule
SolverStrategies
Only available in 1.5

Please replace this functionality by the new strategy API in 1.6 as follows:

direct_umfpack() = DirectSolver(UMFPACKFactorization())                            
gmres_umfpack() = GMRESIteration(UMFPACKFactorization())                           
gmres_eqnblock_umfpack() = GMRESIteration(UMFPACKFactorization(), EquationBlock()) 
gmres_iluzero() = GMRESIteration(ILUZeroPreconditioner())                          
gmres_eqnblock_iluzero() = GMRESIteration(ILUZeroPreconditioner(), EquationBlock())
gmres_pointblock_iluzero() = GMRESIteration(ILUZeroPreconditioner(), PointBlock()) 
source