GitHub

PALEOmodel solvers

Initialization

PALEOmodel.initialize! — Function
initialize!(model::PB.Model; kwargs...) -> (initial_state::Vector, modeldata::PB.ModelData)

Initialize model and return:

  • an initial_state Vector
  • a modeldata struct containing the model arrays.

Initialising the state vector

With default arguments, the model state variables are initialised to the values defined in the .yaml configuration file used to create the model.

The optional pickup_output argument can be used to provide an OutputWriter instance with pickup data to initialise from, using set_statevar_from_output!. This is applied afer the default initialisation, hence can be used to (re)initialise a subset of model state variables.

DataTypes for model arrays

With default arguments, the model arrays use Float64 as the element type. The eltype keyword argument can be used to specify a different Julia DataType, eg for use with automatic differentiation. Per-Variable DataType can be specified by using the :datatype Variable attribute to specify String-valued tags, in combination with the eltypemap keyword argument to provide a Dict of tag names => DataTypes.

Thread safety

A thread-safe model can be created by setting the threadsafe=true global parameter in the YAML config file before calling PB.create_model_from_config (used by Reactions that require eg special handling of global accumulator variables to maintain thread safety), and supplying method_barrier (a thread barrier to add to ReactionMethod dispatch lists between dependency-free groups):

method_barrier = PB.reaction_method_thread_barrier(
    PALEOmodel.ThreadBarriers.ThreadBarrierAtomic("the barrier"),
    PALEOmodel.ThreadBarriers.wait_barrier
)

Keyword summary

  • pickup_output=nothing: OutputWriter with pickup data to initialise from
  • check_units_opt=:no: check linked Variable units are consistent (:no to disable check, :warn to warn and continue, :error to error and stop)
  • eltype::Type=Float64: default data type to use for model arrays
  • eltypemap=Dict{String, DataType}(): Dict of data types to look up Variable :datatype attribute
  • method_barrier=nothing: thread barrier to add to dispatch lists to create a thread-safe model
  • expect_hostdep_varnames=["global.tforce"]: non-state-Variable host-dependent Variable names expected
  • SolverView_all=true: true to create modeldata.solver_view_all
  • create_dispatchlists_all=true: true to create modeldata.dispatchlists_all
  • generated_dispatch=true: true to autogenerate code for modeldata.dispatchlists_all (fast dispatch, slow compile)
source
[deprecated] initialize!(run::Run; kwargs...) -> (initial_state::Vector, modeldata::PB.ModelData)

Call initialize! on run.model.

source
PALEOmodel.set_statevar_from_output! — Function
set_statevar_from_output!(modeldata, output::AbstractOutputWriter) -> initial_state

Initialize model state Variables from last record in output

NB: modeldata must contain solver_view_all

source

DifferentialEquations solvers

Wrappers for the Julia DifferentialEquations package ODE and DAE solvers. These are usually appropriate for smaller biogeochemical models.

NB: see Managing small and negative values for best practices and common issues when using ODE or DAE solvers.

High level wrappers

PALEOmodel.ODE.integrate — Function
integrate(run, initial_state, modeldata, tspan [; kwargs...] ) -> sol::SciMLBase.ODESolution
integrateForwardDiff(run, initial_state, modeldata, tspan [;kwargs...]) -> sol::SciMLBase.ODESolution

Integrate run.model as an ODE or as a DAE with constant mass matrix, and write to outputwriter

Provides a wrapper around the Julia SciML DifferentialEquations package ODE solvers, with PALEO-specific additional setup. Keyword arguments alg and solvekwargs are passed through to the DifferentialEquations solve method.

integrateForwardDiff sets keyword arguments jac_ad=:ForwardDiffSparse, alg=Sundials.CVODE_BDF(linear_solver=:KLU) to use the Julia ForwardDiff package to provide the Jacobian with forward-mode automatic differentiation and automatic sparsity detection.

Implementation

Follows the SciML standard pattern:

  • Create ODEfunction
  • Create SciMLBase.ODEproblem
  • Call SciMLBase.solve

and then

Arguments

  • run::Run: struct with model::PB.Model to integrate and output field
  • initial_state::AbstractVector: initial state vector
  • modeldata::Modeldata: ModelData struct
  • tspan: (tstart, tstop) integration start and stop times

Keywords

  • alg=Sundials.CVODE_BDF(): ODE algorithm to use, passed through to DifferentialEquations.jl solve method. The default is appropriate for a stiff system of equations (common in biogeochemical models), see https://diffeq.sciml.ai/dev/solvers/ode_solve/ for other options.
  • solvekwargs=NamedTuple(): NamedTuple of keyword arguments passed through to DifferentialEquations.jl solve (eg to set abstol, reltol, saveat, see https://diffeq.sciml.ai/dev/basics/common_solver_opts/)
  • outputwriter=run.output: PALEOmodel.AbstractOutputWriter instance to hold output
  • jac_ad=:NoJacobian: Jacobian to generate and use (:NoJacobian, :ForwardDiffSparse, :ForwardDiff)
  • jac_ad_t_sparsity=tspan[1]: model time at which to calculate Jacobian sparsity pattern
  • steadystate=false: true to use DifferentialEquations.jl SteadyStateProblem (not recommended, see PALEOmodel.SteadyState.steadystate).
  • BLAS_num_threads=1: number of LinearAlgebra.BLAS threads to use
  • generated_dispatch=true: true to autogenerate code (fast solve, slow compile)
source
PALEOmodel.ODE.integrateDAE — Function
integrateDAE(run, initial_state, modeldata, tspan [; kwargs...]) -> sol::SciMLBase.DAESolution
integrateDAEForwardDiff(run, initial_state, modeldata, tspan [; kwargs...]) -> sol::SciMLBase.DAESolution

Integrate run.model as a DAE and copy output to outputwriter.

Provides a wrapper around the Julia SciML DifferentialEquations package DAE solvers, with PALEO-specific additional setup. Keyword arguments alg and solvekwargs are passed through to the DifferentialEquations solve method.

integrateDAEForwardDiff sets keyword arguments jac_ad=:ForwardDiffSparse, alg=Sundials.IDA(linear_solver=:KLU) to use the Julia ForwardDiff package to provide the Jacobian with forward-mode automatic differentiation and automatic sparsity detection.

Limitations

  • arbitrary combinations of implicit total (T) and algebraic constraint (C) variables are not supported by Sundials IDA

as used here, as the IDA solver option used to find consistent initial conditions requires a partioning into differential and algebraic variables (see SolverView documentation).

Implementation

Follows the SciML standard pattern:

and then

Keywords

As integrate, with defaults:

  • alg=Sundials.IDA() (integrateDAE)
  • alg=Sundials.IDA(linear_solver=:KLU) (integrateDAEForwardDiff)
source

Low level functions

ODE functions:

PALEOmodel.ODE.ODEfunction — Function
ODEfunction(model::PB.Model, modeldata [; kwargs]) -> SciMLBase.ODEFunction

Contruct SciML ODEfunction https://diffeq.sciml.ai/latest/ with PALEO-specific setup

Keyword arguments are required to generate a Jacobian function (using automatic differentation).

Keywords

  • jac_ad=:NoJacobian: Jacobian to generate and use (:NoJacobian, :ForwardDiffSparse, :ForwardDiff)
  • initial_state::AbstractVector: initial state vector
  • jac_ad_t_sparsity::Float64: model time at which to calculate Jacobian sparsity pattern
source
PALEOmodel.JacobianAD.jac_config_ode — Function
jac_config_ode(
    jac_ad, model, initial_state, modeldata, jac_ad_t_sparsity;
    kwargs...
)-> (jac, jac_prototype::Union{Nothing, SparseMatrixCSC})

Create and return jac (ODE Jacobian function object), and jac_prototype (sparsity pattern as SparseMatrixCSC, or nothing for dense Jacobian)

jac_ad defines Jacobian type (:ForwardDiffSparse, :ForwardDiff)

Adds an array set to modeldata with appropriate datatypes for ForwardDiff AD Dual numbers, sets up cache for ForwardDiff, calculates Jacobian sparsity (if required) at time jac_ad_t_sparsity.

NB: there is a profusion of different Julia APIs here:

  • ForwardDiff Sparse and dense Jacobian use different APIs and have different cache setup requirements.
  • ForwardDiff requires f!(du, u) hence a closure or function object, DifferentialEquations allows context objects to be passed around.

Keyword arguments

  • parameter_aggregator::Union{Nothing, PB.ParameterAggregator}=nothing: provide Jacobian that supports parameters p as defined by parameter_aggregator
  • jac_cellranges=modeldata.cellranges_all: restrict Jacobian to this subset of Domains and Reactions (via operatorID).
  • request_adchunksize=ForwardDiff.DEFAULT_CHUNK_THRESHOLD: chunk size for ForwardDiff automatic differentiation
  • fill_jac_diagonal=true: (jac=:ForwardDiffSparse only) true to fill diagonal of jac_prototype
  • generated_dispatch=true: true to autogenerate code for dispatch (fast dispatch, slow compile)
  • throw_on_nan=false: true to error if NaN detected in Jacobian
  • use_base_vars=String[]: additional Variable full names not calculated by Jacobian, which instead use arrays from modeldata base arrays (arrays_idx=1) instead of allocating new AD Variables
source
PALEOmodel.JacobianAD.paramjac_config_ode — Function
paramjac_config_ode(model::PB.Model, pa::PB.ParameterAggregator, modeldata::PB.ModelData; kwargs...)
    -> paramjac::ParamJacODEForwardDiffDense

Create parameter Jacobian function object paramjac to calculate (dense) parameter Jacobian using ForwardDiff.

A set of model arrays of appropriate Dual number type is added to modeldata.

Keyword arguments

  • filterT=nothing: enable optimisation for the case where only a small number of parameter indices contribute to paramjac at model time t - supply a function filterT(t) -> (pi_1, pi_2, ...) to specify the parameter indices.
  • paramjac_cellranges=modeldata.cellranges_all: restrict Jacobian to this subset of Domains and Reactions (via operatorID).
  • request_adchunksize=ForwardDiff.DEFAULT_CHUNK_THRESHOLD: chunk size for ForwardDiff automatic differentiation
  • generated_dispatch=true: true to autogenerate code for dispatch (fast dispatch, slow compile)
  • use_base_vars=String[]: additional Variable full names not calculated by Jacobian, which instead use arrays from modeldata base arrays (arrays_idx=1) instead of allocating new AD Variables
source

DAE functions:

PALEOmodel.ODE.DAEfunction — Function
DAEfunction(model::PB.Model, modeldata [; kwargs]) -> SciMLBase.DAEFunction

Contruct SciML DAEfunction https://diffeq.sciml.ai/latest/ with PALEO-specific setup

Keyword arguments are required to generate a Jacobian function (using automatic differentation).

Keywords

  • jac_ad=:NoJacobian: Jacobian to generate and use (:NoJacobian, :ForwardDiffSparse, :ForwardDiff)
  • initial_state::AbstractVector: initial state vector
  • jac_ad_t_sparsity::Float64: model time at which to calculate Jacobian sparsity pattern
source
PALEOmodel.ODE.get_inconsistent_initial_deriv — Function
get_inconsistent_initial_deriv(
    initial_state, modeldata, initial_t, differential_vars, modeldae::SolverFunctions.ModelDAE
) -> initial_deriv

NB: IDA initialisation seems not fully understood: with Julia Sundials.IDA(initall=false) (the Sundials.jl default) corresponding to the IDA option `IDAYAYDPINITtoIDACalcIC(), this should "direct IDACalcIC() to compute the algebraic components of y and differential components of ydot, given the differential components of y." But it seems to have some sensitivity toinitial_deriv` (ydot), which shouldn't be used according to the above ?

This function takes a guess at what is needed for initial_deriv:

  • initial_deriv for ODE variables will now be consistent
  • set initial_deriv (constraint for algebraic variables) = 0, this will not be satisfied (ie will not be consistent)
  • initial_deriv for implicit variables will not be consistent
source

Model solution:

PALEOmodel.ODE.print_sol_stats — Function
print_sol_stats(sol::SciMLBase.ODESolution)
print_sol_stats(sol::SciMLBase.DAESolution)
print_sol_stats(sol::SciMLBase.NonlinearSolution)
print_sol_stats(sol::SciMLBase.RODESolution)

Print solution statistics

source
PALEOmodel.ODE.calc_output_sol! — Function
calc_output_sol!(outputwriter, model::PB.Model, sol::SciMLBase.ODESolution, tspan, initial_state, modeldata)
calc_output_sol!(outputwriter, model::PB.Model, sol::SciMLBase.DAESolution, tspan, initial_state, modeldata)
calc_output_sol!(outputwriter, model::PB.Model, sol::SciMLBase.RODESolution, tspan, initial_state, modeldata)
calc_output_sol!(outputwriter, model::PB.Model, sol::SciMLBase.NonlinearSolution, tspan, initial_state, modeldata)
calc_output_sol!(outputwriter, model::PB.Model, tsoln::AbstractVector, soln::AbstractVector,  modeldata)

Iterate through solution and recalculate model fields (functions of state variables and time) and store in outputwriter.

Arguments

  • outputwriter::PALEOmodel.AbstractOutputWriter: container for output
  • model::PB.Model used to calculate solution
  • sol: SciML solution object
  • tspan: (tstart, tstop) integration start and stop times
  • initial_state::AbstractVector: initial state vector
  • tsoln::AbstractVector: solution times
  • soln::AbstractVector: solution state variables
  • modeldata::PB.Modeldata: ModelData struct
source

Steady-state solvers (Julia NLsolve based)

PALEOmodel.SteadyState.steadystate — Function
steadystate(run, initial_state, modeldata, tss [; kwargs...] )
steadystateForwardDiff(run, initial_state, modeldata, tss [; kwargs...] )

Find steady-state solution (using NLsolve.jl package) and write to outputwriter (two records are written, for initial_state and the steady-state solution).

steadystateForwardDiff has default keyword argument jac_ad=:ForwardDiffSparse to use automatic differentiation for sparse Jacobian.

Arguments

  • run::Run: struct with model::PB.Model to integrate and output field
  • initial_state::AbstractVector: initial state vector
  • modeldata::Modeldata: ModelData struct with appropriate element type for forward model
  • tss: (yr) model tforce time for steady state solution

Optional Keywords

  • outputwriter::PALEOmodel.AbstractOutputWriter=run.output: container to hold output
  • initial_time=-1.0: tmodel to write for first output record
  • solvekwargs=NamedTuple(): NamedTuple of keyword arguments passed through to NLsolve.jl (eg to set method, ftol, iteration, show_trace, store_trace).
  • jac_ad: :NoJacobian, :ForwardDiffSparse, :ForwardDiff
  • BLAS_num_threads=1: number of LinearAlgebra.BLAS threads to use
  • generated_dispatch=true: true to use autogenerated code (fast solve, slow compile)
  • [Deprecated: use_norm=false: (must be false)]
source
PALEOmodel.SteadyState.steadystate_ptc — Function
steadystate_ptc(run, initial_state, modeldata, tspan, deltat_initial; kwargs...) 
steadystate_ptcForwardDiff(run, initial_state, modeldata, tspan, deltat_initial; kwargs...)

Find steady-state solution and write to outputwriter, using solve_ptc naive pseudo-transient-continuation with first order implicit Euler pseudo-timesteps.

steadystate_ptcForwardDiff sets keyword argument default jac_ad=:ForwardDiffSparse to use automatic differentiation for sparse Jacobian.

Arguments

  • run::Run: struct with model::PB.Model to integrate and output field
  • initial_state::AbstractVector: initial state vector
  • modeldata::Modeldata: ModelData struct with appropriate element type for forward model
  • tspan: Vector or Tuple with (initial_time, final_time)
  • deltat_initial: initial pseudo-timestep

Keywords in addition to those supported by solve_ptc:

  • jac_ad: AD Jacobian to use (steadystate_ptc default :NoJacobian, steadystate_ptcForwardDiff default :ForwardDiff)
  • request_adchunksize=10: ForwardDiff chunk size to request.
  • jac_cellranges=modeldata.cellranges_all: CellRanges to use for Jacobian calculation (eg to restrict to an approximate Jacobian by using a cellrange with a non-default operatorID: in this case, Variables that are not calculated but needed for the Jacobian should set the transfer_jacobian attribute so that they will be copied)
  • generated_dispatch=true: true for fast solve, slow compile.
source
PALEOmodel.SteadyState.steadystate_ptc_splitdae — Function
steadystate_ptc_splitdae(run, initial_state, modeldata, tspan, deltat_initial; kwargs...)

Find steady-state solution and write to outputwriter, using solve_ptc naive pseudo-transient-continuation with first order implicit Euler pseudo-timesteps, with inner Newton solve for per-cell DAE algebraic constraints (eg quasi-steady-state reaction intermediates)

Arguments

  • run::Run: struct with model::PB.Model to integrate and output field
  • initial_state::AbstractVector: initial state vector
  • modeldata::Modeldata: ModelData struct with appropriate element type for forward model
  • tspan: Vector or Tuple with (initial_time, final_time)
  • deltat_initial: initial pseudo-timestep

Keywords in addition to those supported by solve_ptc:

  • request_adchunksize=10: ForwardDiff chunk size to request for outer Newton solve.
  • jac_cellranges=modeldata.cellranges_all: CellRanges to use for Jacobian calculation (eg to restrict to an approximate Jacobian by using a cellrange with a non-default operatorID: in this case, Variables that are not calculated but needed for the Jacobian should set the transfer_jacobian attribute so that they will be copied)
  • generated_dispatch=true: true for fast solve, slow compile.
  • operatorID_inner=3: operatorID for Reactions to run for inner solve (typically all reservoirs and chemical reactions)
  • transfer_inner_vars=["tmid", "volume", "ntotal", "Abox"]: Variables not calculated by operatorID_inner that need to be copied for inner solve (additional to those with transfer_jacobian set).
  • inner_jac_ad::Symbol=:ForwardDiff: form of automatic differentiation to use for Jacobian for inner NonlinearNewton.solve solver (options :ForwardDiff, :ForwardDiffSparse)
  • inner_start::Symbol=:current: start value for inner solve (options :initial, :current, :zero)
  • inner_kwargs::NamedTuple=(verbose=0, miniters=2, reltol=1e-12, jac_constant=true, project_region=identity): keywords for inner NonlinearNewton.solve solver.
source
PALEOmodel.SteadyState.solve_ptc — Function
solve_ptc(run, initial_state, nlsolveF, tspan, deltat_initial::Float64; kwargs...)

Pseudo-transient continuation using NLsolve.jl

Not called directly: see steadystate_ptc, steadystate_ptc_splitdae

Find the approximate time evolution or accurate steady-state solution for state variables S satisfying the ODE dS/dt = dSdt(S, t) , using naive pseudo-transient-continuation with first order implicit Euler pseudo-timesteps from tspan[1] to tspan[2] and NLsolve.jl as the non-linear solver.

Each pseudo-timestep from time t to t+Δt uses a variant of Newton's method to solve the nonlinear system S_new = S(t) + Δt dSdt(S_new, t+Δt) for S_new = S(t+Δt).

The initial pseudo-timestep Δt is deltat_initial, this is multiplied by deltat_fac for the next iteration until pseudo-time tss_max is reached. If an iteration fails, Δt is divided by deltat_fac and the iteration retried. NB: this is a very naive initial implementation, there is currently no reliable error control to adapt pseudo-timesteps to the rate of convergence, so requires some trial-and-error to set an appropiate deltat_fac for each problem.

Solver options to pass through to the outer NLsolve nonlinear solver are set by solvekwargs

Arguments

  • run::Run: struct with model::PB.Model to integrate and output field
  • initial_state::AbstractVector: initial state vector
  • nlsolveF = (ssFJ!, nldf): function objects created by high-level wrapper.
  • tspan: Vector or Tuple with (initial_time, final_time)
  • deltat_initial: initial pseudo-timestep

Keyword Arguments

  • deltat_fac=2.0: factor to increase pseudo-timestep on success
  • 'dtmax=Inf': maximum pseudo-timestep
  • saveat=Float64[]: Vector of model pseudo-times at which to save output (empty Vector to save all output timesteps)
  • outputwriter=run.output: container for output
  • solvekwargs=NamedTuple(): arguments to pass through to NLsolve
  • maxiters=1000: maximum number of PTC iterations
  • max_failed_iter=20: maximum number of iterations that make no progress (either fail, or no change in solution) before exiting
  • check_callbacks=[]: vector of callback functions to check validity of solution [check_callback(sol, tss, deltat, model, modeldata)]
  • step_callbacks=[]: vector of callback functions to follow a successful step [step_callback(sol, tss, deltat, model, modeldata)], eg to project solution onto invariant manifold (Shampine, 1999)
  • verbose=false: true for detailed output
  • BLAS_num_threads=1: restrict threads used by Julia BLAS (likely irrelevant if using sparse Jacobian?)
  • [Deprecated]: tss_output=Float64[]: renamed to 'saveat'
  • [Deprecated]: - max_iter: renamed to maxiters
  • [Deprecated: - enforce_noneg=nothing: true to fail pseudo-timesteps that generate negative values for state variables.]
  • [Deprecated: use_norm=false: not supported (must be false)]
  • [Deprecated: sol_min: now has no effect, replace with solve_kwargs=(apply_step! = PALEOmodel.SolverFunctions.StepClampAll!(sol_min, sol_max), )]

Function objects nlsolveF = (ssFJ!, nldf)

The nldf function object should adapt ssFJ! to implement the NLsolve interface to provide the residual function f(S_new) = S_new - (S(t) + Δt dSdt(S_new, t+Δt) (and its jacobian) where the equation f(S_new) = 0 defines the new state variable S_new = S(t+Δt) at new time t+Δt after implicit timestep Δt.

The ssFJ! struct should additionally provide a function object modelode to calculate the ODE time derivative dSdt(S, p, t) (including a solve for any DAE algebraic variables), NB: parameters p are not used.

  • ssFJ!.modelode(dSdt, S, p, t)

and should implement:

  • get_state(ssFJ!) -> inner_state: record any state additional to that defined by ODE state variables (eg state variables for short lived species)
  • set_step!(ssFJ!, t+Δt, Δt, S(t), inner_state): reset ssFJ! and hence nldf to calculate the residual f(S_new) for an implicit timestep

Δt to t+Δt from previous state S(t), optionally using inner_state to eg configure an internal solver for DAE algebraic variables.

source

Sparse linear solvers adapted to NLsolve interface:

Function objects to project Newton iterations into valid regions:

PALEOmodel.SolverFunctions.StepClampMultAll! — Type
StepClampMultAll!(minvalue, maxvalue, minmult, maxmult) -> scma!
StepClampMultAll!(minvalue, maxvalue, maxratio) = StepClampMultAll!(minvalue, maxvalue, 1.0/maxratio, maxratio)
scma!(x, x_old, newton_step)

Function object to take Newton step x .= x_old .+ newton_step and then clamp all values in Vector x to specified range using:

  • clamp!(x, x_old*minmult, x_old*maxmult)
  • clamp!(x, minvalue, maxvalue)
source
PALEOmodel.SolverFunctions.StepClampAll! — Type
StepClampAll!(minvalue, maxvalue) -> sca!
sca!(x, x_old, newton_step)

Function object to take Newton step x .= x_old .+ newton_step and then clamp all values in Vector x to specified range using clamp!(x, minvalue, maxvalue)

source
PALEOmodel.SolverFunctions.ClampAll — Type
ClampAll(minvalue, maxvalue) -> ca
ca(v) -> v

Function object to clamp all values in Vector v to specified range using clamp.(v, minvalue, maxvalue) (out-of-place version)

source
PALEOmodel.SolverFunctions.ClampAll! — Type
ClampAll!(minvalue, maxvalue) -> ca!
ca!(v)

Function object to clamp all values in Vector v to specified range using clamp!(v, minvalue, maxvalue) (in-place, mutating version)

source

Function objects to manage outer nonlinear steps:

PALEOmodel.SteadyState.ConservationCallback — Type
ConservationCallback(
    tmodel_start::Float64 # earliest model time to apply correction
    content_name::String # variable with a total of X
    flux_name::String  # variable with a corresponding boundary flux of X
    reservoir_total_name::String  # total for reservoir to apply correction to
    reservoir_statevar_name::String # state variable for reservoir to apply correction to
    reservoir_fac::Float64  # stoichiometric factor (1 / moles of quantity X per reservoir molecule)
) -> ccb

Provides a callback function with signature

ccb(state, tmodel, deltat, model, modeldata)

that modifies modeldata arrays and state to enforce budget conservation.

General principle is to project the solution onto an invariant manifold (modifying state variables), as discussed in (Shampine, 1999)

Example

conservation_callback_H = PALEOmodel.SteadyState.ConservationCallback(
    tmodel_start=1e5,
    content_name="global.content_H_atmocean",
    flux_name="global.total_H",
    reservoir_total_name="atm.CH4_total", # "atm.H2_total",
    reservoir_statevar_name="atm.CH4_mr", # "atm.H2_mr",
    reservoir_fac=0.25 # 0.5, # H per reservoir molecule
)

then add to eg steadystate_ptc_splitdae with step_callbacks keyword argument:

step_callbacks = [conservation_callback_H]
source
PALEOmodel.SteadyState.RestartSmallValuesCallback — Type
RestartSmallValuesCallback(
    stateexplicit::PB.VariableAggregator; 
    include_variables::Vector{String} = String[v.domain.name*"."*v.name for v in stateexplicit.vars],
    exclude_variables::Vector{String} = String[],
    modify_threshold=1e-80,
    modify_val=1e-30
) -> rsvc

Provides a callback function with signature

rsvc(state, tmodel, deltat, model, modeldata)

that modifies state to reset variables with values < modify_threshold to value modify_val, hence restarting the Newton solve from modify_val on the next iteration.

Example

restart_small_values_callback = PALEOmodel.SteadyState.RestartSmallValuesCallback(
    modeldata.solver_view_all.stateexplicit;
    modify_threshold = 1e-80,
    modify_val = 1e-30,
)

then add to eg steadystate_ptc_splitdae with step_callbacks keyword argument:

step_callbacks = [restart_small_values_callback]
source
PALEOmodel.SteadyState.CheckValuesCallback — Type
CheckValuesCallback(
    stateexplicit::PB.VariableAggregator; 
    include_variables::Vector{String} = String[v.domain.name*"."*v.name for v in stateexplicit.vars],
    exclude_variables::Vector{String} = String[],
    check_min_threshold=-Inf,
    check_max_threshold=Inf,
) -> cvc

Provides a callback function with signature

cvc(state, tmodel, deltat, model, modeldata)::Bool

that checks variables have check_min_threshold < values < check_max_threshold.

Example

check_min_values_callback = PALEOmodel.SteadyState.CheckValuesCallback(
    modeldata.solver_view_all.stateexplicit;
    check_min_value = 1e-80,
)

then add to eg steadystate_ptc_splitdae with check_callbacks keyword argument:

check_callbacks = [check_min_values_callback]
source

Steady-state solvers (Sundials Kinsol based):

PALEOmodel.SteadyStateKinsol.steadystate_ptc — Function
steadystate_ptc(run, initial_state, modeldata, tspan, deltat_initial; 
    [,deltat_fac=2.0] [,saveat] [,outputwriter] [,createkwargs] [,solvekwargs]
    [, use_jac_preconditioner] [,jac_cellranges] [, use_directional_ad] [, directional_ad_eltypestomap]
    [,verbose] [,  BLAS_num_threads]
)

Find steady-state solution and write to outputwriter, using naive pseudo-transient-continuation with first order implicit Euler pseudo-timesteps and PALEOmodel.Kinsol as the non-linear solver.

Each pseudo-timestep solves the nonlinear system S(t+Δt) = S(t) + Δt dS/dt(t+Δt) for S(t+Δt), using a variant of Newton's method (preconditioned Newton-Krylov, with the Jacobian as preconditioner)

Initial pseudo-timestep Δt is deltat_initial, this is multiplied by deltat_fac for the next iteration until pseudo-time tss_max is reached. If an iteration fails, Δt is divided by deltat_fac and the iteration retried.

NB: this is a very naive initial implementation, there is currently no reliable error control to adapt pseudo-timesteps to the rate of convergence, so requires some trial-and-error to set an appropiate deltat_fac for each problem.

Solver PALEOmodel.Kinsol options are set by arguments createkwargs (passed through to PALEOmodel.Kinsol.kin_create) and solvekwargs (passed through to PALEOmodel.Kinsol.kin_solve).

If use_jac_ad_preconditioner is true, the AD Jacobian is used as preconditioner. The preconditioner (Jacobian) calculation can be modified by jac_cellranges, to specify a operatorIDs so use only a subset of Reactions in order to calculate an approximate Jacobian to use as the preconditioner.

If use_directional_ad is true, the Jacobian-vector product will be calculated using automatic differentiation (instead of the default finite difference approximation). directional_ad_eltypestomap can be used to specify Variable :datatype tags (strings) that should be mapped to the AD directional derivative datatype hence included in the AD directional derivative.

source
PALEOmodel.Kinsol.kin_create — Function
kin_create(f, y0 [; kwargs...]) -> kin

Create and return a kinsol solver context kin, which can then be passed to kin_solve

Arguments

  • f: Function of form f(fy::Vector{Float64}, y::Vector{Float64}, userdata)
  • y0::Vector template Vector of initial values (used only to define problem dimension)

Keywords

  • userdata: optional user data, passed through to f etc.
  • linear_solver: linear solver to use (only partially implemented, supports :Dense, :Band, :FGMRES)
  • psolvefun: optional preconditioner solver function (for :FGMRES)
  • psetupfun: optional preconditioner setup function
  • jvfun: optional Jacobian*vector function (for :FGMRES)
source
PALEOmodel.Kinsol.kin_solve — Function
kin_solve(
    kin, y0::Vector;
    [strategy] [, fnormtol] [, mxiter] [, print_level] [,y_scale] [, f_scale] [, noInitSetup]
) -> (y, kin_stats)

Solve nonlinear system using kinsol solver context kin (created by kin_create) and initial conditions y0. Returns solution y and solver statistics kinstats. kinstats.returnflag indicates success/failure.

source

Fixed timestep solvers

PALEO native fixed-timestep, first-order Euler integrators, with split-explicit and multi-threaded options. These are usually appropriate for larger biogeochemical models (eg ocean models using GCM transport matrices).

The low-level timestepping is provided by integrateFixed and integrateFixedthreads, with higher-level wrappers for common options provided by integrateEuler etc.

High-level wrappers

PALEOmodel.ODEfixed.integrateEuler — Function
integrateEuler(run, initial_state, modeldata, tspan, Δt [; kwargs])

Integrate run.model from initial_state using first-order Euler with fixed timestep.

Calls integrateFixed

Arguments

  • run::Run: struct with model::PB.Model to integrate and output field
  • initial_state::AbstractVector: initial state vector
  • modeldata::Modeldata: ModelData struct with appropriate element type for forward model
  • tspan: (tstart, toutput1, toutput2, ..., tstop) integration start, output, stop times
  • Δt: (yr) fixed timestep

Keywords

  • outputwriter::PALEOmodel.AbstractOutputWriter=run.output: container to write model output to
  • report_interval=1000: number of timesteps between progress update to console
source
PALEOmodel.ODEfixed.integrateSplitEuler — Function
integrateSplitEuler(run, initial_state, modeldata, tspan, Δt_outer, n_inner;
                        cellranges_outer, cellranges_inner,
                        [,outputwriter])

Integrate run.model representing:

\[\frac{dS}{dt} = f_{outer}(t, S) + f_{inner}(t, S)\]

using split first-order Euler with Δt_outer for f_outer and a shorter timestep Δt_outer/n_inner for f_inner.

f_outer is defined by calling PALEOboxes.do_deriv with cellranges_outer hence corresponds to those Reactions with operatorID of cellranges_outer. f_inner is defined by calling PALEOboxes.do_deriv with cellranges_inner hence corresponds to those Reactions with operatorID of cellranges_inner.

NB: the combined time derivative is written to outputwriter.

Arguments

  • run::Run: struct with model::PB.Model to integrate and output field
  • initial_state::AbstractVector: initial state vector
  • modeldata::Modeldata: ModelData struct with appropriate element type for forward model
  • tspan: (tstart, toutput1, toutput2, ..., tstop) integration start, output, stop times
  • Δt_outer: (yr) fixed outer timestep
  • n_inner: number of inner timesteps per outer timestep

Keywords

  • cellranges_outer: Vector of CellRange with operatorID defining f_outer.
  • cellranges_inner: Vector of CellRange with operatorID defining f_inner.
  • outputwriter::PALEOmodel.AbstractOutputWriter=run.output: container to write model output to
  • report_interval=1000: number of outer timesteps between progress update to console
source
PALEOmodel.ODEfixed.integrateEulerthreads — Function
integrateEulerthreads(run, initial_state, modeldata, cellranges, tspan, Δt;
    outputwriter=run.output, report_interval=1000)

Integrate run.model using first-order Euler with fixed timestep Δt, with tiling over multiple threads.

Arguments

  • run::Run: struct with model::PB.Model to integrate and output field
  • initial_state::AbstractVector: initial state vector
  • modeldata::Modeldata: ModelData struct with appropriate element type for forward model
  • cellranges::Vector{Vector{AbstractCellRange}}: Vector of Vector-of-cellranges, one per thread (so length(cellranges) == Threads.nthreads).
  • tspan: (tstart, toutput1, toutput2, ..., tstop) integration start, output, stop times
  • Δt: (yr) fixed outer timestep

Keywords

  • outputwriter::PALEOmodel.AbstractOutputWriter=run.output: container to write model output to
  • report_interval=1000: number of outer timesteps between progress update to console
source
PALEOmodel.ODEfixed.integrateSplitEulerthreads — Function
integrateSplitEulerthreads(run, initial_state, modeldata, tspan, Δt_outer, n_inner::Int; 
                            cellranges_outer, cellranges_inner, [,outputwriter] [, report_interval])

Integrate run.model using split first-order Euler with Δt_outer for f_outer and a shorter timestep Δt_outer/n_inner for f_inner.

f_outer is defined by calling PALEOboxes.do_deriv with cellrange_outer hence corresponds to those Reactions with operatorID of cellrange_outer. f_inner is defined by calling PALEOboxes.do_deriv with cellrange_inner hence corresponds to those Reactions with operatorID of cellrange_inner.

Uses Threads.nthreads threads and tiling described by cellranges_inner and cellranges_outer (each a Vector of Vector{AbstractCellRange}, one per thread).

Arguments

  • run::Run: struct with model::PB.Model to integrate and output field
  • initial_state::AbstractVector: initial state vector
  • modeldata::Modeldata: ModelData struct with appropriate element type for forward model
  • tspan: (tstart, toutput1, toutput2, ..., tstop) integration start, output, stop times
  • Δt_outer: (yr) fixed outer timestep
  • n_inner: number of inner timesteps per outer timestep (0 for non-split solver)

Keywords

  • cellranges_outer::Vector{Vector{AbstractCellRange}}: Vector of list-of-cellranges, one list per thread (so length(cellranges) == Threads.nthreads), with operatorID defining f_outer.
  • cellranges_inner::Vector{Vector{AbstractCellRange}}: As cellranges_outer, with operatorID defining f_inner.
  • outputwriter::PALEOmodel.AbstractOutputWriter=run.output: container to write model output to
  • report_interval=1000: number of outer timesteps between progress update to console
source
PALEOmodel.ODELocalIMEX.integrateLocalIMEXEuler — Function
integrateLocalIMEXEuler(run, initial_state, modeldata, tspan, Δt_outer [; kwargs...])

Integrate run.model representing:

\[\frac{dS}{dt} = f_{outer}(t, S) + f_{inner}(t, S)\]

using first-order Euler with Δt_outer for f_outer and implicit first-order Euler for f_inner, where f_inner is local (within-cell, ie no transport), for a single Domain, and uses only StateExplicit and Deriv variables.

f_outer is defined by calling PALEOboxes.do_deriv with cellranges_outer hence corresponds to those Reactions with operatorID of cellranges_outer. f_inner is defined by calling PALEOboxes.do_deriv with cellrange_inner hence corresponds to those Reactions with operatorID of cellrange_inner.

NB: the combined time derivative is written to outputwriter.

Arguments

  • run::Run: struct with model::PB.Model to integrate and output field
  • initial_state::AbstractVector: initial state vector
  • modeldata::Modeldata: ModelData struct with appropriate element type for forward model
  • tspan: (tstart, toutput1, toutput2, ..., tstop) integration start, output, stop times
  • Δt_outer: (yr) fixed timestep

Keywords

  • cellranges_outer: Vector of CellRange with operatorID defining f_outer.
  • cellrange_inner: A single CellRange with operatorID defining f_inner.
  • exclude_var_nameroots: State variables that are modified by Reactions in cellrange_inner, but not needed to find implicit solution (ie reaction rates etc don't depend on them).
  • [outputwriter=run.output: PALEOmodel.AbstractOutputWriter instance to write model output to]
  • [report_interval=1000: number of outer timesteps between progress update to console]
  • [Lnorm_inf_max=1e-3: normalized error tolerance for implicit solution]
  • [niter_max=10]: maximum number of Newton iterations
  • [request_adchunksize=4]: request ForwardDiff AD chunk size (will be restricted to an upper limit)
source

Low-level timesteppers

PALEOmodel.ODEfixed.integrateFixed — Function
integrateFixed(run, initial_state, modeldata, tspan, Δt_outer;
            timesteppers, outputwriter=run.output, report_interval=1000)

Fixed timestep integration, with time step implemented by timesteppers,

`timesteppers = [ [(timestep_function, cellranges, timestep_ctxt), ...], ...]`

Where timestep_function(model, modeldata, cellranges, timestep_ctxt, touter, Δt, barrier)

source
PALEOmodel.ODEfixed.integrateFixedthreads — Function
integrateFixedthreads(run, initial_state, modeldata, tspan, Δt_outer;
            timesteppers, outputwriter=run.output, report_interval=1000)

Fixed timestep integration using Threads.nthreads() threads, with time step implemented by timesteppers,

`timesteppers = [ [(timestep_function, cellranges, timestep_ctxt), ... (1 per thread)], ...]`

Where timestep_function(model, modeldata, cellranges, timestep_ctxt, touter, Δt, barrier).

source

Thread barriers

PALEOmodel.ThreadBarriers.ThreadBarrierAtomic — Type
ThreadBarrierAtomic

Thread synchronisation barrier for Threads.nthreads() Julia Threads. Uses a pair of Atomic variables to avoid the need for locks. Resets so can be reused.

Example:

barrier = ThreadBarrierAtomic("my barrier")    
Threads.@threads for t in 1:Threads.nthreads()  
    # do stuff
    
    wait_barrier(barrier)  # blocks until all threads reach this point

    # do more stuff
end
source
PALEOmodel.ThreadBarriers.ThreadBarrierCond — Type
ThreadBarrierCond

Thread synchronisation barrier for Threads.nthreads() Julia Threads using Condition variable and associated lock. Resets so can be reused.

NB: testing on Julia 1.6 suggests this is slow.

Example:

barrier = ThreadBarrierCond("my barrier")    
Threads.@threads for t in 1:Threads.nthreads()  
    # do stuff
    
    wait_barrier(barrier)  # blocks until all threads reach this point

    # do more stuff
end

Implementation:

Uses a condition variable (with associated lock) and a counter. See eg http://web.eecs.utk.edu/~huangj/cs360/360/notes/CondVar/lecture.html

source

Variable aggregation

A SolverView uses a collection of PALEOboxes.VariableAggregators to assemble model state Variables and associated time derivatives into contiguous Vectors, for the convenience of standard numerical ODE / DAE solvers. See Mathematical formulation of the reaction-transport problem.

PALEOmodel.SolverView — Type
SolverView(model, modeldata, arrays_idx; [verbose=true]) 
SolverView(model, modeldata, cellranges; [verbose=false], [indices_from_cellranges=true])

Provides a view on the whole or some part of the Model for a numerical solver.

Contains PALEOboxes.VariableAggregators for a subset of spatial locations, defining state variables S_all composed of subsets stateexplicit (S_e), statetotal (S_t) and state (S_a).

Differential and algebraic constraints are provided by:

  • ODE paired stateexplicit (S_e) and stateexplicit_deriv (dS_e/dt), where dS_e/dt = F(S_all).
  • Implicit-ODE paired total (T) and total_deriv (dT/dt), where dT(S_all)/dt = F(S_all), This implicitly determines the time derivative dS_t/dt, as dT(S_all)/dt = dT(S_all)/dS_all * dS_all/dt. The number of total (T) variables must equal the number of statetotal variables S_t.
  • Algebraic constraints (C), where C(S_all) = 0 and the number of constraint Variables C must equal the number of state variables S_a.

If there are either no total (T) variables, or they are functions T(S_e, S_t) only of S_e and S_t (not of S_a), then the state variables can be partitioned into subsets of differential variables S_d (consisting of stateexplicit (S_e) and statetotal (S_t)), and algebraic variables state (S_a).

The available choice of numerical solvers depends on the types of state variables in the model configuration:

  • A system with only stateexplicit (S_e) variables can be solved by an ODE solver (eg as a SciML ODEProblem by PALEOmodel.ODE.integrate).
  • Some ODE solvers (those that support a constant mass matrix on the LHS) can also solve systems with constraints C (eg as a SciML ODEProblem by PALEOmodel.ODE.integrate).
  • A DAE solver such as Sundials IDA is required to solve systems with implicit total (T) variables (eg as a SciML DAEProblem by PALEOmodel.ODE.integrateDAE). NB: arbitrary combinations of total (T) and constraint (C) variables are not supported by Sundials IDA via PALEOmodel.ODE.integrateDAE, as the IDA solver option used to find consistent initial conditions requires a partioning into differential and algebraic variables as described above.

Optional access methods provide an ODE/DAE solver view with composite statevar and statevar_sms,where:

- `statevar` is a concatenation of `stateexplicit`, `statetotal`, and `state` ([`set_statevar!`](@ref))
- `statevar_sms` is a concatenation of `stateexplicit_deriv`, `total_deriv`, `constraints` ([`get_statevar_sms!`](@ref))

Constructors create a SolverView for the entire model from modeldata array set arrays_idx, or for a subset of Model Variables defined by the Domains and operatorIDs of cellranges.

Keywords

  • indices_from_cellranges=true: true to restrict to the index ranges from cellranges, false to just use cellranges to define Domains

and take the whole of each Domain.

  • hostdep_all=true: true to include host dependent not-state Variables from all Domains
  • reallocate_hostdep_eltype=Float64: a data type to reallocate hostdep Variables eg to replace any

AD types.

source
PALEOmodel.set_default_solver_view! — Function
set_default_solver_view!(model, modeldata)

(Optional, used to set modeldata.solver_view_all to a SolverView) for the whole model, and set modeldata.hostdep_data to any non-state-variable host dependent Variables)

reallocate_hostdep_eltype a data type to reallocate hostdep_data eg to replace any AD types.

source

Function objects

Function objects are callable structs with function signatures required by DifferentialEquations or other solvers to calculate model time derivative, Jacobian, etc. They combine variable aggregation (using PALEOboxes.VariableAggregators or PALEOmodel.SolverView) with corresponding Reaction dispatch lists.

ODE function objects

PALEOmodel.SolverFunctions.ModelODE — Type
ModelODE(
    modeldata, [parameter_aggregator]; 
    solver_view=modeldata.solver_view_all,
    dispatchlists=modeldata.dispatchlists_all
) -> f::ModelODE

Function object to calculate model time derivative and adapt to SciML ODE solver interface

Call as f(du,u, p, t)

source
PALEOmodel.SolverFunctions.JacODEForwardDiffDense — Type
JacODEForwardDiffDense(
    modeldata; 
    solver_view=modeldata.solver_view_all,
    dispatchlists=modeldata.dispatchlists_all,
    du_template, 
    jacconf,
) -> jac::JacODEForwardDiffDense

JacODEForwardDiffDense_p(
    modeldata, pa::PB.ParameterAggregator; 
    solver_view=modeldata.solver_view_all,
    dispatchlists=modeldata.dispatchlists_all,
    du_template, 
    jacconf,
) -> jac::JacODEForwardDiffDense_p

Function object to calculate dense Jacobian in form required for SciML ODE solver.

solver_view, dispatchlists should correspond to modeldata, which should have the appropriate element type for ForwardDiff Dual numbers.

Call as jac(J, u, p, t)

source
PALEOmodel.SolverFunctions.JacODEForwardDiffSparse — Type
JacODEForwardDiffSparse(
    modeldata; 
    solver_view=modeldata.solver_view_all,
    dispatchlists=modeldata.dispatchlists_all,
    du_template,
    throw_on_nan, 
    jac_cache,
) -> jac::JacODEForwardDiffSparse

Function object to calculate sparse Jacobian in form required for SciML ODE solver.

solver_view, dispatchlists should correspond to modeldata, which should have the appropriate element type for ForwardDiff Dual numbers.

Call as jac(J, u, p, t)

source

DAE function objects

PALEOmodel.SolverFunctions.ModelDAE — Type
ModelDAE

Function object to calculate model residual G and adapt to SciML DAE solver interface.

If using Total variables, odeimplicit should be an ImplicitForwardDiffDense or ImplicitForwardDiffSparse, otherwise nothing.

Provides function signature:

(fdae::ModelDAE)(G, dsdt, s, p, t)

where residual G(dsdt,s,p,t) is:

  • -dsdt + F(s) (for ODE-like state Variables s with time derivative F given explicitly in terms of s)
  • F(s) (for algebraic constraints)
  • duds*dsdt + F(s, u(s)) (for Total variables u that depend implicitly on state Variables s)
source