API updates

Source

Here we describe some updates for the API of VoronoiFVM.jl. These have been implemented mostly on top of the existing API, whose functionality is not affected.

This is a breaking release. Implementations using default solver settings should continue to work (albeit possibly with deprecation and allocation warnings). Really breaking is control of iterative linear solvers and allocation checks.

Solve now a method of CommonSolve.solve

As a consequence, all VoronoiFVM.solve methods with signatures others than solve(system; kwargs...) are now deprecated

Deprecated call:

begin inival = unknowns(sys0; inival = 0.1) sol00 = unknowns(sys0) solve!(sol00, inival, sys0) end

Replace this by:

Docstring of solve

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

Docstring of SolverControl

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]')

• 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 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}

Rely on LinearSolve.jl for linear system solution

This provides easy access to a large variety of linear solvers:

LU factorization from UMFPACK

LU factorization from Sparspak.jl

Iterative solvers

BICGstab from Krylov.jl with diagonal (Jacobi) preconditioner

The Jacobi preconditioner is defined in ExtendableSparse.jl.

BICGstab from Krylov.jl with delayed factorization preconditioner

BICGstab from Krylov.jl with ilu0 preconditioner

ILUZeroPreconditioner is exported from ExtendableSparse and wraps the predonditioner defined in ILUZero.jl .

New verbosity handling

• verbose can now be a Bool or a String of flag characters, allowing for control of different output categories. I would love to do this via logging, but there is still a long way to go IMHO

• Allocation check is active by default with warnings which can be muted by passing a verbose string without 'a'. This is now the only control in this respect. All check_allocs methods/kwargs, control via environment variables have been removed.

• Deprecation warnings can be switched off by passing a verbose string without 'd'.

• Improve iteration logging etc., allow for logging of linear iterations ('l' flag character)

The following example gives some information in this respect:

If we find these warnings annoying, we can switch them off:

Or we get some more logging:

But we can also look for the reasons of the allocations. Here, global values should be declared as constants.

VoronoiFVM.System constructor

Implicit creation of physics

The VoronoiFVM.Physics struct almost never was used outside of the constructor of VoronoiFVM.System. Now it is possible to specify the flux functions directly in the system constructor. By default, it is also possible to set a list of species which are attached to all interior and boundary regions of the grid.

Boundary conditions as part of physics

This makes the API more consistent and opens an easy possibility to have space and time dependent boundary conditions. One can specify them either in breaction or the synonymous bcondition.

time:

4.99

Implicit creation of grid

By passing data for grid creation (one to three abstract vectors) instead a grid, a tensor product grid is implicitly created. This example also demonstrates position dependent boundary values.

GridVisualize API extended to System

Instead of a grid, a system can be passed to gridplot and scalarplot.

Parameters of solve

The solve API has been simplified and made more Julian. All entries of VoronoiFVM.NewtonControl can be now passed as keyword arguments to solve.

Another new keyword argument is inival which allows to pass an initial value which by default is initialized to zero. Therefore we now can write solve(system) as we already have seen above.