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
- Solver control
- Linear solver stragies
- Block preconditioning
- History handling
- Matrix extration
Solve method
CommonSolve.solve
— Methodsolve(system; kwargs...)
Built-in solution method for VoronoiFVM.System
.
Keyword arguments:
General for all solvers
inival
(default: 0) : Array created viaunknowns
or number giving the initial value.control
(default: nothing): Pass instance ofSolverControl
- All elements of
SolverControl
can be used as kwargs. Eventually overwrites values given viacontrol
params
: Parameters (Parameter handling is experimental and may change)
Stationary solver: Invoked if neither
times
norembed
, nortstep
are given as keyword argument.time
(default:0.0
): Set time value.
Returns a
DenseSolutionArray
orSparseSolutionArray
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 objectsol
containing the stored solution(s), seeTransientSolution
.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 exactlypre
(default:(sol,t)->nothing
): callback invoked before each time steppost
(default:(sol,oldsol, t, Δt)->nothing
): callback invoked after each time stepsample
(default:(sol,t)->nothing
): callback invoked after timestep for all times intimes[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. Ifcontrol.Δt<control.Δt_min
, solution is aborted with error. Returns a transient solution objectsol
containing the stored solution, seeTransientSolution
.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
orSparseSolutionArray
Solver control
VoronoiFVM.SolverControl
— TypeSolverControl
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
- true: "neda"
- false: "da"
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
occurredmax_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 givennothing
, as default are chosen (forFloat64
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
.- 1D:
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 functionprecon_linear(A)
which takes an AbstractSparseMatrixCSC (e.g. an ExtendableSparseMatrix) and returns a preconditioner object in the sense ofLinearSolve.jl
, i.e. which has anldiv!(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. Iftrue
, 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
: Functionpre(sol,t)
called before time/embedding step
post::Function
: Functionpost(sol,oldsol,t,Δt)
called after successful time/embedding step
sample::Function
: Functionsample(sol,t)
to be called for eacht 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}
Linear solver strategies
VoronoiFVM.LinearSolverStrategy
— TypeVoronoiFVM.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:
UMFPACKFactorization
(using LinearSolve
)KLUFactorization
(using LinearSolve
)SparspakFactorization
(using LinearSolve
),SparspakLU
(using ExtendableSparse,Sparspak
)MKLPardisoLU
(using ExtendableSparse, Pardiso
), openmp parallelAMGSolver
(using AMGCLWrap
), openmp parallelRLXSolver
(using AMGCLWrap
), openmp parallel
Notable incomplete factorizations/preconditioners
- Incomplete LU factorizations written in Julia:
ILUZeroPreconditioner
ILUTPrecondidtioner
(using ExtendableSparse, IncompleteLU
)
- Algebraic multigrid written in Julia: (
using ExtendableSparse, AlgebraicMultigrid
) - AMGCL based preconditioners (
using ExtendableSparse, AMGCLWrap
), openmp parallel
Blocking strategies are:
VoronoiFVM.DirectSolver
— TypeDirectSolver(;factorization=UMFPACKFactorization())
LU Factorization solver.
VoronoiFVM.CGIteration
— TypeCGIteration(;factorization=UMFPACKFactorization())
CGIteration(factorization)
CG Iteration from Krylov.jl via LinearSolve.jl.
VoronoiFVM.BICGstabIteration
— TypeBICGstabIteration(;factorization=UMFPACKFactorization())
BICGstabIteration(factorization)
BICGstab Iteration from Krylov.jl via LinearSolve.jl.
VoronoiFVM.GMRESIteration
— TypeGMRESIteration(;factorization=ILUZeroFactorization(), memory=20, restart=true)
GMRESIteration(factorization; memory=20, restart=true)
GMRES Iteration from Krylov.jl via LinearSolve.jl.
Block preconditioning
This feature is under development as of 1.6.
VoronoiFVM.BlockStrategy
— TypeVoronoiFVM.BlockStrategy
Abstract supertype for various block preconditioning strategies.
VoronoiFVM.NoBlock
— TypeNoBlock()
No blocking.
VoronoiFVM.EquationBlock
— TypeEquationBlock()
Equation-wise blocking. Can be combined with any preconditioner resp. factorization including direct solvers.
VoronoiFVM.PointBlock
— TypePointBlock()
Point-wise blocking. Currently only together with ILUZeroFactorization. This requires a system with unknown_storage=:dense
.
VoronoiFVM.Equationwise
— Typestruct Equationwise
Equationwise partitioning mode.
VoronoiFVM.partitioning
— Functionpartitioning(system, _)
Calculate partitioning of system unknowns.
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.NewtonSolverHistory
— Typemutable 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 factorizationsnlin::Int64
: number of linear iteration steps / factorization solvestime::Float64
: Elapsed time for solutiontasm::Float64
: Elapsed time for assemblytlinsolve::Float64
: Elapsed time for linear solveupdatenorm::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$
VoronoiFVM.TransientSolverHistory
— Typemutable 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 iterationtimes::Any
: Time valuesupdates::Any
: Update norms used for step control
Base.summary
— Methodsummary(h::NewtonSolverHistory)
Return named tuple summarizing history.
Base.summary
— Methodsummary(h::TransientSolverHistory)
Return named tuple summarizing history.
VoronoiFVM.details
— Functiondetails(h::NewtonSolverHistory)
Return array of named tuples with info on each iteration step
details(h::TransientSolverHistory)
Return array of details of each solver step
VoronoiFVM.history
— Functionhistory(tsol)
Return solver history if log
was set to true. See see NewtonSolverHistory
, TransientSolverHistory
.
history(sys)
Return solver history from last solve
call, if log
was set to true. See see NewtonSolverHistory
, TransientSolverHistory
.
VoronoiFVM.history_details
— Functionhistory_details(tsol)
Return details of solver history from last solve
call, if log
was set to true. See details
.
history_details(sys)
Return details of solver history from last solve
call, if log
was set to true. See details
.
VoronoiFVM.history_summary
— Functionhistory_summary(tsol)
Return summary of solver history from last solve
call, if log
was set to true.
history_summary(sys)
Return summary of solver history from last solve
call, if log
was set to true.
Matrix extraction
For testing and teaching purposes, one can obtain residual and linearization at a given vector of unknowns
VoronoiFVM.evaluate_residual_and_jacobian
— Functionevaluate_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.
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
— MethodODEProblem(system,inival,tspan,callback=SciMLBase.CallbackSet())
Create an ODEProblem in mass matrix form which can be handeled by ODE solvers from DifferentialEquations.jl.
Parameters:
system
: AVoronoiFVM.System
inival
: Initial value. Consider to pass a stationary solution attspan[1]
.tspan
: Time intervalcallback
: (optional) callback for ODE solver
The method returns an ODEProblem which can be solved by solve().
Base.reshape
— Methodreshape(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.
SciMLBase.ODEFunction
— Method 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
: AVoronoiFVM.System
jacval
(optional): Initial value. Default is a zero vector. Consider to pass a stationary solution at timetjac
.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.
Legacy API
During the development of the code, a number of API variants have been developed which are supported for backward compatibility.
CommonSolve.solve
— Method solve(inival, system, times; kwargs...)
Alias for solve(system::VoronoiFVM.AbstractSystem; kwargs...)
with the corresponding keyword arguments.
CommonSolve.solve
— Method solve(inival, system; control=SolverControl(),params, tstep=Inf)
Alias for solve(system::VoronoiFVM.AbstractSystem; kwargs...)
with the corresponding keyword arguments.
Solve stationary problem(if tstep==Inf
) or one step implicit Euler step using Newton's method with inival
as initial value. Returns a solution array.
CommonSolve.solve!
— Methodsolve!(solution, inival, system;
control=SolverControl(),
tstep=Inf)
Mutating version of solve(inival,system)
VoronoiFVM.NewtonControl
— TypeNewtonControl
Legacy name of SolverControl
VoronoiFVM.SolverStrategies
— ModuleSolverStrategies
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())