Reaction API
API calls used by AbstractReaction
implementations.
Implementations should create a subclass of AbstractReaction
containing a ReactionBase
and Parameter
s, optionally provide a create_reaction
, and implement callbacks to define VariableReaction
s and register ReactionMethod
s.
Reaction struct
PALEOboxes.AbstractReaction
— TypeAbstractReaction
Abstract base Type for Reactions.
Implementation
Derived types should include a field base::
ReactionBase
, and usually a ParametersTuple
, eg
Base.@kwdef mutable struct ReactionHello{P} <: PB.AbstractReaction
base::PB.ReactionBase
pars::P = PB.ParametersTuple(
PB.ParDouble("my_par", 42.0, units="yr",
description="an example of a Float64-valued scalar parameter called 'my_par'"),
)
some_additional_field::Float64 # additional fields to eg cache data read from disk etc
end
Derived types should implement register_methods!
, and may optionally implement create_reaction
, set_model_geometry
, check_configuration
, register_dynamic_methods!
.
Methods should be registered using add_method_setup!
, add_method_initialize!
, add_method_do!
.
Any parameters not included in pars
should be added explicitly with add_par
(this is rarely needed).
PALEOboxes.ReactionBase
— TypeReactionBase
Base Type for a biogeochemical Reaction.
Implementation
Include as field base
in types derived from AbstractReaction
Parameters
PALEOboxes.AbstractParameter
— TypePALEOboxes.Parameter
— TypeParameter{T, ParseFromString}
A reaction parameter of type T
.
Create using short names ParDouble
, ParInt
, ParBool
, ParString
.
Read value as <par>[]
, set with setvalue!
.
Parameters with external=true
may be set from the Model-level Parameters list, if name
is present in that list.
ParseFromString
should usually be Nothing
: a value of Type T
is then required when calling setvalue!
. If ParseFromString
is not Nothing
, then setvalue!
will accept an AbstractString
and call Base.parse(ParseFromString, strvalue)
. This allows eg an enum-valued Parameter to be defined by Parameter{EnumType, EnumType} and implementing parse(EnumType, rawvalue::AbstractString)
PALEOboxes.VecParameter
— TypeVecParameter{T, ParseFromString}
A reaction parameter of type Vector{T}
.
Create using short names ParDoubleVec
, ParStringVec
.
Read values as <par>[i]
, access raw Vector{T}
as <par>.v
.
Set with setvalue!
, in config file use standard yaml syntax for a vector eg [1, 2, 3]
See Parameter
for additional documentation.
Implementation
Implements part of the AbstractVector
interface, sufficient to access elements as <par>[i]
, and to support iteration eg for v in <par>; ...; end
.
PALEOboxes.VecVecParameter
— TypeVecVecParameter{T, ParseFromString}
A reaction parameter of type Vector{Vector{T}}.
Create using short names ParDoubleVecVec
.
Read values as <par>.v::Vector{Vector{T}}
.
Set using standard yaml syntax for a vector of vectors eg [[1, 2, 3], [4, 5, 6]]
See Parameter
for additional documentation.
PALEOboxes.ParametersTuple
— FunctionParametersTuple(parameters::AbstractParameter...) -> NamedTuple
ParametersTuple(parameters) -> NamedTuple
Create a NamedTuple of Parameters.
PALEOboxes.add_par
— Functionadd_par(reaction::AbstractReaction, par::AbstractParameter)
add_par(reaction::AbstractReaction, objectwithpars)
Add a single parameter or parameters from fields of objectwithpars
to a new Reaction.
Not usually needed: Parameters in pars::ParametersTuple
will be added automatically, only needed if there are additional Parameters that are not members of
pars`.
PALEOboxes.setvalue!
— Functionsetvalue!(par::Parameter, value)
Set Parameter to value
.
Optionally (if Parameter
has Type parameter ParseFromString != Nothing
) parse value
from a String.
Registering Reactions with the PALEOboxes framework
All subtypes of AbstractReaction
that are present in loaded modules are available to the PALEO framework. Available Reactions can be listed with find_reaction
and find_all_reactions
. The default create_reaction
is called to create Reactions when the model is created (this method can be overridden if needed).
PALEOboxes.create_reaction
— Methodcreate_reaction(ReactionType::Type{<:AbstractReaction}, base::ReactionBase) -> reaction::AbstractReaction
Default method to create a ReactionType
and set base
field.
A reaction implementation may optionally implement a custom method eg to set additional fields
PALEOboxes.find_reaction
— Functionfind_reaction(class::AbstractString) -> ReactionType
Look up "class" in list of Reactions available in currently loaded modules (using find_all_reactions
), and return fully-qualified Reaction Type (including module prefixes).
PALEOboxes.find_all_reactions
— Functionfind_all_reactions() -> Dict{String, Type}
Use InteractiveUtils.subtypes(AbstractReaction)
to find all currently loaded subtypes off AbstractReaction, and create a Dict
with last part of the name of the Type as key (ie without the module prefix) and Type as value.
Any Types that generate non-unique keys (eg Module1.MyReactionType and Module2.MyReactionType) will generate a warning, and no entry will be added to the Dict (so if this Reaction is present in a config file, it will not be found and will error).
Defining Domain Grids and array sizes
PALEOboxes.set_model_geometry
— Functionset_model_geometry(reaction, model)
Optional: define Domain
grid
, data_dims
.
One Reaction per Domain may create a grid (an AbstractMesh
subtype) and set the domain.grid
field.
Multiple Reactions per domain may call set_data_dimension!
to define (different) named data dimensions (eg wavelength grids).
PALEOboxes.set_data_dimension!
— Functionset_data_dimension!(domain::Domain, dim::NamedDimension [, coordinates], ; allow_exists=false)
Define a Domain data dimension as a NamedDimension
, optionally attaching a Vector of coordinate names.
Variables may then specify data dimensions as a list of names using the :data_dims
Variable Attribute.
Creating and registering Reaction methods
All Reactions should implement register_methods!
, and may optionally implement register_dynamic_methods!
.
PALEOboxes.register_methods!
— Functionregister_methods!(reaction)
register_methods!(reaction, model)
Add ReactionMethod
s, using add_method_setup!
, add_method_initialize!
add_method_do!
. See also register_dynamic_methods!
.
PALEOboxes.register_dynamic_methods!
— Functionregister_dynamic_methods!(reaction)
register_dynamic_methods!(reaction, model)
Optional: called after first variable link pass, to add ReactionMethod
s that depend on Variables generated by other Reactions (see register_methods!
).
These methods should then define one or more ReactionMethod
s, which requires:
- Defining collections of
VariableReaction
s (see Defining VariableReactions, Defining collections of VariableReactions). - Implementing a function to iterate over model cells and calculate the required fluxes etc (see Implementing method functions)
- Adding methods (see Adding ReactionMethods).
In addition it is possible to add Predefined ReactionMethods for some common operations (Variable initialisation, calculating totals, etc).
Defining VariableReactions
PALEOboxes.VariableReaction
— TypeVariableReaction(VT, localname => link_namestr, units, description; attributes=Tuple()) -> VariableReaction{VT}
VariableReaction(VT, linklocal_namestr, units, description; attributes=Tuple()) -> VariableReaction{VT}
VarProp, VarPropScalar, VarPropStateIndep, VarPropScalarStateIndep -> VariableReaction{VT_ReactProperty}
VarDep, VarDepColumn, VarDepScalar, VarDepStateIndep, VarDepColumnStateIndep, VarDepScalarStateIndep -> VariableReaction{VT_ReactDependency}
VarTarget, VarTargetScalar -> VariableReaction{VT_ReactTarget}
VarContrib, VarContribColumn, VarContribScalar -> VariableReaction{VT_ReactContributor}
VarStateExplicit, VarStateExplicitScalar -> VariableReaction{VT_ReactDependency}
VarDeriv, VarDerivScalar -> VariableReaction{VT_ReactContributor}
VarState, VarStateScalar -> VariableReaction{VT_ReactDependency}
VarConstraint, VarConstraintScalar -> VariableReaction{VT_ReactDependency}
[deprecated] VariableReaction(VT, localname, units, description; link_namestr, attributes=Tuple()) -> VariableReaction{VT}
Reaction view on a model variable.
Reactions define AbstractVarList
s of VariableReaction
s when creating a ReactionMethod
. Within a ReactionMethod
, a VariableReaction
is referred to by localname
. When the model is initialised, VariableDomain
variables are created that link together VariableReaction
s with the same link_namestr
, and data Arrays are allocated. views
on the VariableDomain
data Arrays are then passed to the ReactionMethod
function at each timestep.
Subtypes
The Type parameter VT
is one of VT_ReactProperty
, VT_ReactDependency
, VT_ReactContributor
, VT_ReactTarget
, where short names are defined for convenience:
const VarPropT = VariableReaction{VT_ReactProperty}
const VarDepT = VariableReaction{VT_ReactDependency}
const VarTargetT = VariableReaction{VT_ReactTarget}
const VarContribT = VariableReaction{VT_ReactContributor}
There are two pairings of VariableReaction
s with VariableDomain
s:
- Reaction Property and Dependency Variables, linked to a
VariableDomPropDep
. These are used to represent a quantity calculated in one Reaction that is then used by other Reactions. - Reaction Target and Contributor Variables, linked to a
VariableDomContribTarget
. These are used to represent a flux-like quantity, with one Reaction definining the Target and multiple Reactions adding contributions.
Variable Attributes and constructor convenience functions
Variable attributes are used to define Variable :space
AbstractSpace
(scalar, per-cell, etc) and data content :field_data
AbstractData
, and to label state Variables for use by numerical solvers.
VariableReaction
is usually not called directly, instead convenience functions are defined that provide commonly-used combinations of VT
and attributes
:
short name | VT | attributes | |||||
---|---|---|---|---|---|---|---|
:space | :field_data | :vfunction | :initialize_to_zero | :datatype | :is_constant | ||
VarProp | VT_ReactProperty | CellSpace | ScalarData | VF_Undefined | - | - | false |
VarPropScalar | VT_ReactProperty | ScalarSpace | ScalarData | VF_Undefined | - | - | false |
VarPropStateIndep | VT_ReactProperty | CellSpace | ScalarData | VF_Undefined | - | Float64 | true |
VarPropScalarStateIndep | VT_ReactProperty | ScalarSpace | ScalarData | VF_Undefined | - | Float64 | true |
VarDep | VT_ReactDependency | CellSpace | UndefinedData | VF_Undefined | - | - | false |
VarDepColumn | VT_ReactDependency | ColumnSpace | UndefinedData | VF_Undefined | - | - | false |
VarDepScalar | VT_ReactDependency | ScalarSpace | UndefinedData | VF_Undefined | - | - | false |
VarDepStateIndep | VT_ReactDependency | CellSpace | UndefinedData | VF_Undefined | - | Float64 | true |
VarDepColumnStateIndep | VT_ReactDependency | ColumnSpace | UndefinedData | VF_Undefined | - | Float64 | true |
VarDepScalarStateIndep | VT_ReactDependency | ScalarSpace | UndefinedData | VF_Undefined | - | Float64 | true |
VarTarget | VT_ReactTarget | CellSpace | ScalarData | VF_Undefined | true | - | false |
VarTargetScalar | VT_ReactTarget | ScalarSpace | ScalarData | VF_Undefined | true | - | false |
VarContrib | VT_ReactContributor | CellSpace | UndefinedData | VF_Undefined | - | - | false |
VarContribColumn | VT_ReactContributor | ColumnSpace | UndefinedData | VF_Undefined | - | - | false |
VarContribScalar | VT_ReactContributor | ScalarSpace | UndefinedData | VF_Undefined | - | - | false |
VarStateExplicit | VT_ReactDependency | CellSpace | ScalarData | VF_StateExplicit | - | - | false |
VarStateExplicitScalar | VT_ReactDependency | ScalarSpace | ScalarData | VF_StateExplicit | - | - | false |
VarDeriv | VT_ReactContributor | CellSpace | ScalarData | VF_Deriv | true | - | false |
VarDerivScalar | VT_ReactContributor | ScalarSpace | ScalarData | VF_Deriv | true | - | false |
VarState | VT_ReactDependency | CellSpace | ScalarData | VF_State | - | - | false |
VarStateScalar | VT_ReactDependency | ScalarSpace | ScalarData | VF_State | - | - | false |
VarConstraint | VT_ReactContributor | CellSpace | ScalarData | VF_Constraint | true | - | false |
VarConstraintScalar | VT_ReactContributor | ScalarSpace | ScalarData | VF_Constraint | true | - | false |
This illustrates some general principles for the use of attributes:
- All Variables must define the
:space
VariableAttribute (a subtype ofAbstractSpace
) to specify whether they are Domain scalars, per-cell quantities, etc. This is used to define array dimensions, and to check that dimensions match when variables are linked. - The
:field_data
attribute (a subtype ofAbstractData
) specifies the quantity that Property and Target Variables represent. This defaults toScalarData
to represent a scalar value. To eg represent a single isotope the:field_data
attribute should be set toIsotopeLinear
. Dependency and Contributor Variables with:field_data = UndefinedData
then acquire this value when they are linked, or may specify:field_data
to constrain allowed links. - The
:initialize_to_zero
attribute is set for Target variables, this is than used (by the ReactionMethod created byadd_method_initialize_zero_vars_default!
) to identify variables that should be initialised to zero at the start of each timestep. - The
:vfunction
attribute is used to label state Variables and corresponding time derivatives, for access by a numerical solver.- An ODE-like combination of a state variable and time derivative are defined by a paired
VarStateExplicit
andVarDeriv
. Note that that these are justVarDep
andVarContrib
with the:vfunction
attribute set, and that there is noVarProp
andVarTarget
defined in the model (these are effectively provided by the numerical solver). The pairing is defined by the naming convention ofvarname
andvarname_sms
. - An algebraic constraint (for a DAE) is defined by a
VarState
andVarConstraint
. Note that that these are justVarDep
andVarContrib
with the:vfunction
attribute set, and that there is noVarProp
andVarTarget
defined in the model (these are effectively provided by the numerical solver). These variables are not paired. - The :initializetozero attribute is also set for Contributor variables VarDeriv and VarConstraint (as there is no corresponding Target variable in the model).
- The :field_data attribute should be set on labelled state etc Variables (as there are no corresponding Property or Target variables in the model to define this).
- An ODE-like combination of a state variable and time derivative are defined by a paired
- The
:is_constant
attribute is used to identify constant Property Variables (not modified after initialisation). A Dependency Variable with :is_constant = true can only link to a constant Property Variable. - The
:datatype
attribute is used both to provide a concrete datatype for constant Variables, and to exclude a non-constant Variable from automatic differentiation (TODO document that usage).
Additional attributes can be specified to provide model-specific information, with defaults defined in the Reaction .jl code that can often then be overridden in the .yaml configuration file, see StandardAttributes
. Examples include:
:initial_value
,:norm_value
,:initial_delta
for state variables or constant variables. NB: the Reaction creating these variables is responsible for implementing a setup method to read the attributes and set the variable data array appropriately at model initialisation.- An
:advect
attribute is used to label tracer variables to indicate that they should have advective transport applied by a transport Reaction included in the model.
NB: after Variables are linked to Domain Variables, the attributes used are those from the master
Variable (either a Property or Target variable, or a labelled state variable Dependency or Contributor with no corresponding Property or Target). Additional attributes must therefore be set on this master Variable.
Specifying links
localname
identifies the VariableReaction
within the Reaction
, and can be used to set variable_attributes:
and variable_links:
in the .yaml configuration file.
linkreq_domain.linkreq_subdomain.linkreq_name
defines the Domain, Subdomain and name for run-time linking to VariableDomain
variables.
Arguments
VT::VariableType
: one ofVT_ReactProperty
,VT_ReactDependency
,VT_ReactContributor
,VT_ReactTarget
localname::AbstractString
: Reaction-local Variable namelink_namestr::AbstractString
:<linkreq_domain>.[linkreq_subdomain.]linkreq_name
. Parsed byparse_variablereaction_namestr
to define the requested linking toDomain
Variable.linklocal_namestr::AbstractString
:<linkreq_domain>.[linkreq_subdomain.]localname
. Convenience form to define bothlocalname
and requested linking to Domain Variable, for the common case wherelinkreq_name == localname
.units::AbstractString
: units ("" if not applicable)description::AbstractString
: text describing the variable
Keywords
attributes::Tuple(:attrb1name=>attrb1value, :attrb2name=>attrb2value, ...)
: variable attributes, seeStandardAttributes
,set_attribute!
,get_attribute
PALEOboxes.parse_variablereaction_namestr
— Functionparse_variablereaction_namestr(linkstr)
-> (linkreq_domain, linkreq_subdomain, linkreq_name, link_optional)
Parse a linkstr into component parts.
linkstr
is of format: [(][<linkreq_domain>.][<linkreq_subdomain>.]<linkreq_name>[)]
- Optional brackets
( ... )
setlink_optional=true
linkreq_name
may contain%reaction%
which will later be substituted with<Reaction name>/
Examples:
julia> PALEOboxes.parse_variablereaction_namestr("foo") # Common case eg for a property that should be public
("", "", "foo", false)
julia> PALEOboxes.parse_variablereaction_namestr("%reaction%foo") # Reaction-private by default
("", "", "%reaction%foo", false)
julia> PALEOboxes.parse_variablereaction_namestr("ocean.foo") # Request link to variable of same name in ocean Domain
("ocean", "", "foo", false)
julia> PALEOboxes.parse_variablereaction_namestr("(ocean.oceansurface.goo)") # Full syntax
("ocean", "oceansurface", "goo", true)
PALEOboxes.set_attribute!
— Functionset_attribute!(var::VariableBase, name::Symbol, value, allow_create=false) -> var
Set Variable attribute.
Defining collections of VariableReactions
PALEOboxes.AbstractVarList
— TypeAbstractVarList
Variables required by a ReactionMethod
methodfn
are specified by a Tuple of VarList_xxx <: AbstractVarList, each containing a collection of VariableReaction
.
These are then converted (by the create_accessors
method) to a corresponding Tuple of collections of views on Domain data arrays , which are then be passed to the ReactionMethod
methodfn
.
NB: creates and uses a copy of supplied Variables including metadata, so set / modify Variable attributes before creating a VarList.
Implementation
Subtypes of AbstractVarList
should implement:
- a constructor that takes a collection of VariableReactions
create_accessors
, returning the views on Domain data arrays in a subtype-specific collection.get_variables
, returning the collection of VariableReactions (as a flat list).
PALEOboxes.VarList_single
— TypeVarList_single(var; components=false) -> VarList_single
Create a VarList_single
describing a single VariableReaction
, create_accessors
will then return a single accessor.
PALEOboxes.VarList_namedtuple
— TypeVarList_namedtuple(varcollection; components=false) -> VarList_namedtuple
Create a VarList_namedtuple
describing a collection of VariableReaction
s, create_accessors
will then return a NamedTuple with field names = VariableReaction.localname
and field values = corresponding data arrays.
If components = true
, each NamedTuple field will be a Vector of data array components.
PALEOboxes.VarList_tuple
— TypeVarList_tuple(varcollection; components=false) -> VarList_tuple
Create a VarList_tuple
describing a collection of VariableReaction
s, create_accessors
will then return a Tuple of data arrays.
An entry of nothing
in varcollection
will produce nothing
in the corresponding entry in the Tuple of arrays supplied to the Reaction method.
If components = true
, each Tuple field will be a Vector of data array components.
PALEOboxes.VarList_ttuple
— TypeVarList_ttuple(varcollection) -> VarList_ttuple
Create a VarList_ttuple
describing a collection of collections of VariableReaction
s, create_accessors
will then return a Tuple of Tuples of data arrays.
PALEOboxes.VarList_vector
— TypeVarList_vector(varcollection; components=false, forceview=false) -> VarList_vector
Create a VarList_vector
describing a collection of VariableReaction
s, create_accessors
will then return a Vector of data arrays.
If components = true
, each Vector element will be a Vector of data array components.
If forceview = true
, each accessor will be a 1-D view
to help type stability, even if this is redundant (ie no view
required, v::Vector -> view(v, 1:length(v)))
PALEOboxes.VarList_vvector
— TypeVarList_vvector(Vector{Vector{VariableReaction}}::vars; components=false) -> VarList_vvector
Create a VarList_vvector
describing a Vector of Vectors of VariableReaction
s, create_accessors
will then return a Vector of Vectors of data arrays.
If components = true
, each Vector of Vectors element will be a Vector of data array components.
PALEOboxes.VarList_nothing
— TypeVarList_nothing() -> VarList_nothing
Create a placeholder for a missing/unavailable VariableReaction. create_accessors
will then return nothing
.
Implementing method functions
Reaction method functions should iterate over the cells in the Domain supplied by cellrange
argument and calculate appropriate biogeochemical fluxes etc (which may include the model time derivative and any intermediate or diagnostic output).
Iterating over cells
The simplest case is a method function that iterates over individual cells, with skeleton form:
function do_something_cellwise(m::PB.AbstractReactionMethod, pars, (vars, ), cellrange::PB.AbstractCellRange, deltat)
@inbounds for i in cellrange.indices
vars.A[i] = something*vars.B[i]*vars.C[i] # in general A is some function of B, C, etc
# etc
end
return nothing
end
Iterating over cells in columns
If necessary (eg to calculate vertical transport), provided the model grid and cellrange allow, it is possible to iterate over columns and then cells within columns (in order from top to bottom):
function do_something_columnwise(m::PB.AbstractReactionMethod, pars, (vars, ), cellrange::PB.AbstractCellRange, deltat)
@inbounds for (icol, colindices) in cellrange.columns
accum = zero(vars.A[first(colindices)]) # accumulator of appropriate type
for i in colindices # order is top to bottom
accum += vars.A[i]
vars.C[i] = accum # C = sum of A in cells above
# etc
end
vars.floor_C[icol] = vars.C[last(colindices)] # assumes model has a floor domain with one floor cell per column in the interior domain
end
return nothing
end
Iteration from bottom to top within a column can be implemented using Iterators.reverse
, eg
function do_something_columnwise(m::PB.AbstractReactionMethod, pars, (vars, ), cellrange::PB.AbstractCellRange, deltat)
@inbounds for (icol, colindices) in cellrange.columns
colreverse = Iterators.reverse(colindices)
for i in colreverse # order is bottom to top
# etc
end
end
return nothing
end
The method function shouldn't make any assumptions about colindices
other than that it is a list of indices ordered from top to bottom in a column. Depending on the grid in use, the indices may not be contiguous, and may not be integers.
The example above made the additional assumption that a floor
domain had been defined (containing Variable floor_C
) with one floor cell per column. This is determined by the model configuration, and is not true in general.
In rare cases where it is necessary to operate on a Vector representing a quantity for the whole column (rather than just iterate through it), this can be implemented using view
, eg
function do_something_columnwise(m::PB.AbstractReactionMethod, pars, (vars, ), cellrange::PB.AbstractCellRange, deltat)
@inbounds for (icol, colindices) in cellrange.columns
A_col = view(vars.A, colindices) # A_col is an AbstractVector with contiguous indices 1:length(colindices)
B_col = view(vars.B, colindices) # B_col is an AbstractVector with contiguous indices 1:length(colindices)
# do something that needs a vector of cells for a whole column
end
return nothing
end
Writing automatic-differentiation-compatible code
Stiff systems of equations (with a wide range of timescales) will require a solver that uses Jacobian information (either calculated numerically or by automatic differentiation).
To be numerically efficient, this means that the Jacobian needs to exist mathematically, so any reaction rates should change continuously. Any discontinuous steps (eg a rate that suddenly increases to a finite value at a threshold) should be smoothed out, eg using the smoothstepcubic
function or ideally by reformulating the physical model in such a way that discontinuities are avoided.
For large models, a sparse Jacobian may be required, where the sparsity pattern is generated by tracing the code with the model variables set eg to their initial values. Conditional logic (if-then-else blocks or elseif) may then cause this tracing to fail (omit some terms from the Jacobian) if the conditional branch taken sets some rates to zero losing dependency information. Workaround is to use zero_ad
to generate a zero value that retains variable dependency information. NB: the max
and min
Julia functions are safe to use with sparsity tracing.
It is sometimes desirable to omit terms from the Jacobian, eg an insignificantly small term that would greatly reduce the Jacobian sparsity. This can be accomplished by using the value_ad
function to drop automatic-differentiation information and retain only the value.
PALEOboxes.smoothstepcubic
— Functionsmoothstepcubic(x, xedge, xwidth) -> y
Smoothed step function over width xwidth
at location xedge
.
Provides a way of smoothing a discontinuous step function so that the derivative exists, to help numerical solution methods that require a Jacobian (eg stiff ODE solvers).
Uses a cubic function in range xedge +/- xwidth/2
, so first derivative is continuous, higher derivatives are not.
Uses zero_ad
to retain AD dependency information for tracing Jacobian sparsity pattern.
Returns:
- 0.0 for x < (xedge - xwidth/2)
- 1.0 for x > (xedge + xwidth/2)
- a smoothed step for xedge-xwidth/2 < x < xedge+xwidth/2
PALEOboxes.zero_ad
— Functionzero_ad(x...) -> 0.0*x
Provide a zero of type of x
(or type of x1*x2*...
if given multiple arguments), retaining AD dependency information.
Workaround to enable use of conditional logic whilst retaining dependency information for tracing Jacobian sparsity pattern.
PALEOboxes.value_ad
— Functionvalue_ad(x::ADT) -> x
Get scalar value from variable x
(discarding any AD derivatives).
This can be used to exclude x
from automatic differentation and hence a Jacobian calculation, eg where x
is a small contribution but would make the Jacobian much denser.
Model code should implement this for any AD types used, eg
value_ad(x::SparsityTracing.ADval) = SparsityTracing.value(x)
value_ad(x::ForwardDiff.Dual) = ForwardDiff.value(x)
Writing thread-safe code
Reactions that accumulate per-cell quantities into a scalar variable (eg to calculate Domain totals) and are intended for use in large models that use a multithreaded solver with tiled Domains should use atomic_add!
:
PALEOboxes.atomic_add!
— Functionatomic_add!(values::AbstractArray, v)
Add v
to values[]
using Thread-safe atomic operation.
This is a generic fallback for ScalarData containing values
of any data_type
.
Optimising loops over cells using explicit SIMD instructions
Reactions with simple loops over cellindices
that implement time-consuming per-cell calculations may be optimised by using explicit SIMD instructions.
PALEOboxes.SIMDutils.SIMDIter
— TypeSIMDIter(baseiter, Val{N})
SIMDIter(baseiter, ::Type{SIMD.Vec{N, U}})
SIMDIter(baseiter, ::Val{1}) # scalar fallback
SIMDIter(baseiter, ::Type{U}) where {U <: Real} # scalar fallback
Iterator that takes up to N
SIMD elements at a time from baseiter
(which should represent indices into a Vector). See Julia package SIMD.jl
If baseiter
contained 1 or more but less then N
elements, then indices
is filled with repeats of the last available element.
Returns Tuple of indices (length N
).
Examples
v_a = [1.0, 2.0, 3.0, 4.0, 5.0]
v_b = similar(v_a)
iter = eachindex(v_a) # iter should represent indices into a Vector
# simplest version - Float64 x 4, ie type of v_a x 4
for i in SIMDIter(iter, Val(4))
x = v_a[i] # x is a packed SIMD vector
v_b[i] = x
end
# with type conversion - Float32 x 8, ie explicitly change Type of SIMD vector
ST = SIMD.Vec{8, Float32}}
for i in SIMDIter(iter, ST)
# v = vec[i] <--> vgatherind(ST, vec, i)
# vec[i] = v <--> vscatterind!(v, vec, i)
# vec[i] += v <--> vaddind!(v, vec, i)
x = vgatherind(ST, v_a, i) # x is a packed SIMD vector with type conversion to Float32
vscatterind!(x, v_b, i)
end
Adding ReactionMethods
PALEOboxes.ReactionMethod
— TypeReactionMethod(
methodfn::Function,
reaction::AbstractReaction,
name::String,
varlists::Tuple{Vararg{AbstractVarList}},
p,
operatorID::Vector{Int64},
domain::AbstractDomain;
preparefn = (m, vardata) -> vardata
) -> m::ReactionMethod
Defines a callback function methodfn
with Variables varlists
, to be called from the Model framework either during setup or as part of the main loop.
Fields
methodfn
: callback from Model frameworkreaction
: the Reaction that created this ReactionMethodname
: a descriptive name, eg generated from the name of methodfnvarlists
: Tuple of VarLists, each representing a list ofVariableReaction
s. Corresponding Variable accessorsvardata
(views on Arrays) will be provided to themethodfn
callback. NB: not concretely typed to reduce compile time, as not performance-criticalp
: optional context field (of arbitrary type) to store data needed by methodfn.operatorID
domain
preparefn
: preparefn(m::ReactionMethod, vardata::Tuple) -> modified_vardata::Tuple optionally modifyvardata
to eg add buffers. NB: not concretely typed as not performance-critical
methodfn
The methodfn
callback is:
methodfn(m::ReactionMethod, pars, vardata::Tuple, cellrange::AbstractCellRange, modelctxt)
or (if Parameters are not required):
methodfn(m::ReactionMethod, vardata::Tuple, cellrange::AbstractCellRange, modelctxt)
With arguments:
m::ReactionMethod
: context is available asm.reaction::AbstractReaction
(the Reaction that defined theReactionMethod
), andm.p
(an arbitrary extra context field supplied whenReactionMethod
created).pars
: a struct with Parameters as fields (current just the ParametersTuple defined asreaction.pars
)vardata
: A Tuple of collections of views on Domain data arrays corresponding toVariableReaction
s defined byvarlists
cellrange::AbstractCellRange
: range of cells to calculate.modelctxt
:- for a setup method,
:setup
,:initial_value
or:norm_value
defining the type of setup requested - for a main loop method
deltat
providing timestep information eg for rate throttling.
- for a setup method,
preparefn
An optional preparefn
callback can be supplied eg to allocate buffers that require knowledge of the data types of vardata
or to cache expensive calculations:
preparefn(m::ReactionMethod, vardata::Tuple) -> modified_vardata::Tuple
This is called after model arrays are allocated, and prior to setup.
PALEOboxes.add_method_setup!
— Functionadd_method_setup!(reaction::AbstractReaction, method::AbstractReactionMethod)
add_method_setup!(reaction::AbstractReaction, methodfn::Function, vars::Tuple{Vararg{AbstractVarList}}; kwargs...) -> ReactionMethod
Add or create-and-add a setup method (called before main loop) eg to set persistent data or initialize state variables. methodfn
, vars
, kwargs
are passed to ReactionMethod
.
PALEOboxes.add_method_initialize!
— Functionadd_method_initialize!(reaction::AbstractReaction, method::AbstractReactionMethod)
add_method_initialize!(reaction::AbstractReaction, methodfn::Function, vars::Tuple{Vararg{AbstractVarList}}; kwargs...) -> ReactionMethod
Add or create-and-add an initialize method (called at start of each main loop iteration) eg to zero out accumulator Variables. methodfn
, vars
, kwargs
are passed to ReactionMethod
.
PALEOboxes.add_method_do!
— Functionadd_method_do!(reaction::AbstractReaction, method::AbstractReactionMethod)
add_method_do!(reaction::AbstractReaction, methodfn::Function, vars::Tuple{Vararg{AbstractVarList}}; kwargs...) -> ReactionMethod
Add or create and add a main loop method. methodfn
, vars
, kwargs
are passed to ReactionMethod
.
Predefined ReactionMethods
Setup and initialization of Variables
PALEOboxes.add_method_setup_initialvalue_vars_default!
— Functionadd_method_setup_initialvalue_vars_default!(react::AbstractReaction, variables [; kwargs...])
Create and add a default method to initialize Variables matching filterfn
(defaults to state Variables) at beginning of integration.
Setup callbacks used
- State Variables and similar (
:vfunction != VF_Undefined
) are initialized in a setup callback withattribute_name in (:initial_value, :norm_value)
, with values from those Variable attributes. - If
force_state_norm_value=false
, other Variables (with:vfunction == VF_Undefined
) are initialized in a setup callback withattribute_name=:setup
, with values from the:initial_value
Variable attribute. NB:filterfn
must be set to include these Variables. - If
force_initial_norm_value=true
, all Variables (including those with:vfunction == VF_Undefined
) are initialised as state Variables
Keywords
filterfn
: set to f(var)::Bool to override the default selection for state variables only (Variables with:vfunction in (VF_StateExplicit, VF_State, VF_Total, VF_StateTotal, VF_Constraint)
)force_initial_norm_value=false
:true
to always use:initial_value
,:norm_value
, even for variables with:vfunction=VF_Undefined
transfer_attribute_vars=[]
: Set to a list of the same length asvariables
to initialisevariables
from attributes oftransfer_attribute_vars
.setup_callback=(method, attribute_name, var, vardata) -> nothing
: Set to a function that is called after each Variable initialisation eg to store norm values.convertvars=[]
convertfn = (convertvars_tuple, i) -> 1.0
convertinfo = ""
Including volume etc conversion
Set convertvars
to a Vector of Variables (eg for cell volume) and supply convertfn
and convertinfo
to initialize to :initial_value*convertfn(convertvars_tuple, i)
where the argument of convertfn
is the Tuple generated by VarList_tuple(convertvars)
.
Example: To interpret :initial_value
as a concentration-like quantity:
convertvars = [volume],
convertfn = ((volume, ), i) -> volume[i],
convertinfo = " * volume"
PALEOboxes.add_method_initialize_zero_vars_default!
— Functionadd_method_initialize_zero_vars_default!(react::AbstractReaction, variables=PB.get_variables(react))
Create and add a default method to initialize Variables to zero at beginning of each timestep. Defaults to adding all Variables from react
with :initialize_to_zero
attribute true
.
NB: TODO variables are converted to VarDep (no dependency checking or sorting needed, could define a VarInit or similar?)
Adding totals for Variables
PALEOboxes.add_method_do_totals_default!
— Functionadd_method_do_totals_default!(react::AbstractReaction, total_candidates=PB.get_variables(react);
[filterfn] [, methodname] [, total_localnames] [, operatorID])
Create and add a method to add total variables (Scalar Properties), for Variables in collection total_candidates
that match filterfn
(defaults to those that are Array Variables and have attribute `:calc_total == true
).
NB: total Variables will require initialization to zero using add_method_initialize_zero_vars_default!
Chemical reactions
PALEOboxes.RateStoich
— TypeRateStoich(
ratevartemplate, stoich_statevarname;
deltavarname_eta=nothing, prcessname="", sms_prefix="", sms_suffix="_sms"
) -> RateStoich
Calculate fluxes for a biogeochemical reaction given rate, stoichiometry, and optionally isotope eta.
Add to a Reaction using create_ratestoich_method
and add_method_do!
.
A Property Variable should be set to provide the reaction rate (often this is implemented by another method of the same Reaction). This method will then link to that (using the local and link names supplied by ratevartemplate
) and calculate the appropriate product rates, omitting products that are not present (VariableReaction
not linked) in the Model
configuration. Metadata for use when analysing model output should be added to the rate variable using add_rate_stoichiometry!
, in the usual case where this Variable is supplied as ratevartemplate
this will happen automatically.
Arguments:
ratevartemplate::Union{VarPropT, VarDepT}
: used to define the rate variable local and link names.stoich_statevarname
: collection of Tuple(stoichiometry, name) eg ((-2.0, "O2"), (-1.0,"H2S::Isotope"), (1.0, "SO4::Isotope"))deltavarname_eta
: optional tuple of variable delta + eta ("SO4_delta", -30.0) or ("SO4_delta", rj.pars.delta). If a Parameter is supplied, this is read indo_react_ratestoich
to allow modification.processname::String
: optional tag to identify the type of reaction in model outputadd_rate_stoichiometry=true
:true
to add calladd_rate_stoichiometry!
to add metadata toratevartemplate
.
Examples:
Create a RateStoich
representing the reaction 2 O2 + H2S -> H2SO4
julia> myratevar = PALEOboxes.VarProp("myrate", "mol yr-1", "a rate");
julia> rs = PALEOboxes.RateStoich(myratevar, ((-2.0, "O2"), (-1.0,"H2S"), (1.0, "SO4")));
julia> rs.stoich_statevarname
((-2.0, "O2"), (-1.0, "H2S"), (1.0, "SO4"))
PALEOboxes.create_ratestoich_method
— Functioncreate_ratestoich_method(reaction::AbstractReaction, ratestoich::RateStoich; isotope_data=ScalarData)
-> ReactionMethod
Create method (see RateStoich
).
PALEOboxes.add_rate_stoichiometry!
— Functionadd_rate_stoichiometry!(ratevar::VarPropT, ratestoich::RateStoich)
Add metadata to rate variable ratevar
for use when analysing model output.
Only needs to be called explicitly if RateStoich
was supplied with a VarDep that links to the rate variable, not the rate variable itself.
Adds Variable attributes:
rate_processname::String = ratestoich.processname
rate_species::Vector{String}
reactants + products fromratestoich.stoich_statevarname
rate_stoichiometry::Vector{Float64}
reaction stoichiometry fromratestoich.stoich_statevarname
Internal details of Variable arrays accessor generation
VariableReaction
s in the AbstractVarList
s for a ReactionMethod
are processed by create_accessor
to supply view
s on arrays as corresponding arguments to the ReactionMethod
function.
PALEOboxes.create_accessor
— Function create_accessor(var::VariableReaction, modeldata, arrays_idx, components, [,forceview=false])
-> accessor or (accessor, subdomain_indices)
Creates a view
on a (single) VariableDomain
data array linked by var::
VariableReaction
. Called by an AbstractVarList
create_accessors
implementation to generate a collection of views
for multiple VariableReaction
s.
Returns:
- if
var
is linked, anaccessor
or Tuple(accessor, subdomain_indices)
that provides aview
on variable data. - if
var
is not linked,nothing
ifvar
is optional, or errors and doesn't return ifvar
is non-optional.
Mapping of indices for Subdomains <–> Domains:
- if no Subdomain, returns unmodified indices (if
forceview=false
), or an equivalent view (ifforceview=true
, this is to help type stability) - if Variable is a Domain interior and Subdomain is a Domain boundary,
accessor
is a view with a subset of Domain indices. - if Variable is a Domain boundary and Subdomain is the Domain interior, returns a Tuple with
subdomain_indices
, length(Subdomain size), withmissing
for interior points.
Mapping of multi-component (Isotope) Variables:
- If
components=false
:- map multi-component Variable to
accessor::IsotopeArray
- return a single-component Variable as a
accessor::AbstractArray
.
- map multi-component Variable to
- If
components=true
:- variable data as a
accessor::Vector{Array}
, length=number of components
- variable data as a
PALEOboxes.create_accessors
— Functioncreate_accessors(varlist::AbstractVarList, modeldata::AbstractModelData, arrays_idx::Int) -> vardata
Return a collection vardata
of views on Domain data arrays for VariableReactions in varlist
. Collection and view are determined by varlist
Type.