diff --git a/dev/.documenter-siteinfo.json b/dev/.documenter-siteinfo.json index 715b35a..60f72e6 100644 --- a/dev/.documenter-siteinfo.json +++ b/dev/.documenter-siteinfo.json @@ -1 +1 @@ -{"documenter":{"julia_version":"1.6.7","generation_timestamp":"2025-01-05T16:00:52","documenter_version":"1.8.0"}} \ No newline at end of file +{"documenter":{"julia_version":"1.6.7","generation_timestamp":"2025-01-06T20:39:54","documenter_version":"1.8.0"}} \ No newline at end of file diff --git a/dev/HOWTOshowmodelandoutput/index.html b/dev/HOWTOshowmodelandoutput/index.html index b00401a..137d7f3 100644 --- a/dev/HOWTOshowmodelandoutput/index.html +++ b/dev/HOWTOshowmodelandoutput/index.html @@ -1,2 +1,2 @@ -Displaying model configuration and output from the Julia REPL · PALEOmodel Documentation
GitHub
+Displaying model configuration and output from the Julia REPL · PALEOmodel Documentation
GitHub
diff --git a/dev/HOWTOsmallnegativevalues/index.html b/dev/HOWTOsmallnegativevalues/index.html index fdecc9d..8948edf 100644 --- a/dev/HOWTOsmallnegativevalues/index.html +++ b/dev/HOWTOsmallnegativevalues/index.html @@ -1,2 +1,2 @@ -Managing small and negative values · PALEOmodel Documentation
GitHub

Managing small and negative values

It is common for biogeochemical reservoirs to both (i) be required to be non-negative, and (ii) approach zero (eg oxygen below the oxic layer in a sediment). This requires some explicit management to allow the numerical ODE / DAE solvers to operate stably and efficiently.

PALEO follows the recommended best practice for Sundials CVODE and other adaptive solvers (also including those in MATLAB), which is to allow -ve values within the error tolerance, and set rates to zero when this happens.

This is a FAQ for CVODE, https://computing.llnl.gov/projects/sundials/faq#cvode_negvals, '... Remember that a small negative value in y returned by CVODE, with magnitude comparable to abstol or less, is equivalent to zero as far as the computation is concerned. .... '. See (Shampine et al., 2005) for a detailed discussion.

There are three areas that need to be addressed:

  1. When calculating biogeochemical reaction rates that depend on myvar, use max(myvar, 0.0) or similar everywhere to set rates to zero for -ve values. Linear transport processes (eg diffusion, flux transport) should transport -ve values to maintain conservation properties.
  2. Set the abstol solver option to control errors for near-zero values of state Variables, see DifferentialEquations solvers. The default value will often be too high. In some cases, it may be most efficient to tolerate -ve values, in other cases, it may be most efficient to control errors using a combination of reltol and abstol so that -ve values are not generated. The easiest way to set abstol in PALEO is to use abstol=1e-5*PALEOmodel.get_statevar_norm(modeldata.solver_view_all) to set to a small fraction (here 1e-5, experimentation will be needed) of the state variables' norm_value attributes (these are set in the .yaml configuration file).
  3. Defend against -ve values when using plots with log scales by explicitly setting an axis lower limit eg ylim=(1e-9, Inf) (without this, the autoscaling will fail and possible produce strange-looking plots eg with inverted axes). See https://docs.juliaplots.org/latest/generated/attributes_axis/
+Managing small and negative values · PALEOmodel Documentation
GitHub

Managing small and negative values

It is common for biogeochemical reservoirs to both (i) be required to be non-negative, and (ii) approach zero (eg oxygen below the oxic layer in a sediment). This requires some explicit management to allow the numerical ODE / DAE solvers to operate stably and efficiently.

PALEO follows the recommended best practice for Sundials CVODE and other adaptive solvers (also including those in MATLAB), which is to allow -ve values within the error tolerance, and set rates to zero when this happens.

This is a FAQ for CVODE, https://computing.llnl.gov/projects/sundials/faq#cvode_negvals, '... Remember that a small negative value in y returned by CVODE, with magnitude comparable to abstol or less, is equivalent to zero as far as the computation is concerned. .... '. See (Shampine et al., 2005) for a detailed discussion.

There are three areas that need to be addressed:

  1. When calculating biogeochemical reaction rates that depend on myvar, use max(myvar, 0.0) or similar everywhere to set rates to zero for -ve values. Linear transport processes (eg diffusion, flux transport) should transport -ve values to maintain conservation properties.
  2. Set the abstol solver option to control errors for near-zero values of state Variables, see DifferentialEquations solvers. The default value will often be too high. In some cases, it may be most efficient to tolerate -ve values, in other cases, it may be most efficient to control errors using a combination of reltol and abstol so that -ve values are not generated. The easiest way to set abstol in PALEO is to use abstol=1e-5*PALEOmodel.get_statevar_norm(modeldata.solver_view_all) to set to a small fraction (here 1e-5, experimentation will be needed) of the state variables' norm_value attributes (these are set in the .yaml configuration file).
  3. Defend against -ve values when using plots with log scales by explicitly setting an axis lower limit eg ylim=(1e-9, Inf) (without this, the autoscaling will fail and possible produce strange-looking plots eg with inverted axes). See https://docs.juliaplots.org/latest/generated/attributes_axis/
diff --git a/dev/MathematicalFormulation/index.html b/dev/MathematicalFormulation/index.html index 1c3b72c..97c5156 100644 --- a/dev/MathematicalFormulation/index.html +++ b/dev/MathematicalFormulation/index.html @@ -1,2 +1,2 @@ -Mathematical formulation of the reaction-transport problem · PALEOmodel Documentation
GitHub

Mathematical formulation of the reaction-transport problem

The PALEO models define various special cases of a general DAE problem (these can be combined, providing the number of implicit state variables $S_{impl}$ is equal to the number of algebraic constraints $G$ plus the number of total variables $U$):

Explicit ODE

The time derivative of explicit state variables $S_{explicit}$ (a subset of all state variables $S_{all}$) are an explicit function of time $t$:

\[\frac{dS_{explicit}}{dt} = F(S_{all}, t)\]

where explicit state variables $S_{explicit}$ are identified by PALEO attribute :vfunction = PALEOboxes.VF_StateExplicit and paired time derivatives $F$ by :vfunction = PALEOboxes.VF_Deriv along with the naming convention <statevarname>, <statevarname>_sms.

Algebraic constraints

State variables $S_{impl}$ (a subset of all state variables $S_{all}$) are defined by algebraic constraints $G$:

\[0 = G(S_{all}, t)\]

where implicit state variables $S_{impl}$ are identified by PALEO attribute :vfunction = PALEOboxes.VF_State and algebraic constaints $G$ by :vfunction = PALEOboxes.VF_Constraint (these are not paired).

ODE with variable substitution

State variables $S_{impl}$ (a subset of all state variables $S_{all}$) are defined the time evolution of total variables $U(S_{all})$ (this case is common in biogeochemistry where the total variables $U$ represent conserved chemical elements, and the state variables eg primary species):

\[\frac{dU(S_{all})}{dt} = F(U(S_{all}), t)\]

where total variables $U$ are identified by PALEO attribute :vfunction = PALEOboxes.VF_Total and paired time derivatives $F$ by :vfunction = PALEOboxes.VF_Deriv along with the naming convention <totalvarname>, <totalvarname>_sms, and implicit state variables $S_{impl}$ are identified by PALEO attribute :vfunction = PALEOboxes.VF_StateTotal.

+Mathematical formulation of the reaction-transport problem · PALEOmodel Documentation
GitHub

Mathematical formulation of the reaction-transport problem

The PALEO models define various special cases of a general DAE problem (these can be combined, providing the number of implicit state variables $S_{impl}$ is equal to the number of algebraic constraints $G$ plus the number of total variables $U$):

Explicit ODE

The time derivative of explicit state variables $S_{explicit}$ (a subset of all state variables $S_{all}$) are an explicit function of time $t$:

\[\frac{dS_{explicit}}{dt} = F(S_{all}, t)\]

where explicit state variables $S_{explicit}$ are identified by PALEO attribute :vfunction = PALEOboxes.VF_StateExplicit and paired time derivatives $F$ by :vfunction = PALEOboxes.VF_Deriv along with the naming convention <statevarname>, <statevarname>_sms.

Algebraic constraints

State variables $S_{impl}$ (a subset of all state variables $S_{all}$) are defined by algebraic constraints $G$:

\[0 = G(S_{all}, t)\]

where implicit state variables $S_{impl}$ are identified by PALEO attribute :vfunction = PALEOboxes.VF_State and algebraic constaints $G$ by :vfunction = PALEOboxes.VF_Constraint (these are not paired).

ODE with variable substitution

State variables $S_{impl}$ (a subset of all state variables $S_{all}$) are defined the time evolution of total variables $U(S_{all})$ (this case is common in biogeochemistry where the total variables $U$ represent conserved chemical elements, and the state variables eg primary species):

\[\frac{dU(S_{all})}{dt} = F(U(S_{all}), t)\]

where total variables $U$ are identified by PALEO attribute :vfunction = PALEOboxes.VF_Total and paired time derivatives $F$ by :vfunction = PALEOboxes.VF_Deriv along with the naming convention <totalvarname>, <totalvarname>_sms, and implicit state variables $S_{impl}$ are identified by PALEO attribute :vfunction = PALEOboxes.VF_StateTotal.

diff --git a/dev/PALEOmodelOutput/index.html b/dev/PALEOmodelOutput/index.html index d3f750a..aeee929 100644 --- a/dev/PALEOmodelOutput/index.html +++ b/dev/PALEOmodelOutput/index.html @@ -1,19 +1,19 @@ -PALEOmodel output · PALEOmodel Documentation
GitHub

PALEOmodel output

Run

PALEOmodel.Run — Type
Run

Container for model and output.

Fields

  • model::Union{Nothing, PB.Model}: The model instance.
  • output::Union{Nothing, AbstractOutputWriter}: model output
source

Output

PALEO model output is accumulated into a container such as an OutputMemory instance that implements the AbstractOutputWriter interface.

AbstractOutputWriter interface

Writing output

PALEOmodel.OutputWriters.initialize! — Function
initialize!(
+PALEOmodel output · PALEOmodel Documentation
GitHub
PALEOmodel output
Run
PALEOmodel.Run — Type
Run
Container for model and output.
Fields
  • model::Union{Nothing, PB.Model}: The model instance.
  • output::Union{Nothing, AbstractOutputWriter}: model output
source
Output
PALEO model output is accumulated into a container such as an OutputMemory instance that implements the AbstractOutputWriter interface.
AbstractOutputWriter interface
Writing output
PALEOmodel.OutputWriters.initialize! — Function
initialize!(
     output::PALEOmodel.AbstractOutputWriter, model, modeldata, [nrecords] 
     [;record_dim_name=:tmodel] [record_coord_units="yr"]
-)
Initialize from a PALEOboxes::Model, optionally reserving memory for an assumed output dataset of nrecords.
The default for record_dim_name is :tmodel, for a sequence of records following the time evolution of the model.
source
PALEOmodel.OutputWriters.add_record! — Function
add_record!(output::PALEOmodel.AbstractOutputWriter, model, modeldata, rec_coord)
Add an output record for current state of model at record coordinate rec_coord. The usual case (set by initialize!) is that the record coordinate is model time tmodel.
source
Modifying output
PALEOboxes.add_field! — Method
add_field!(output::PALEOmodel.AbstractOutputWriter, fr::PALEOmodel.FieldRecord)
Add PALEOmodel.FieldRecord fr to output, with Domain and name defined by fr.attributes[:var_domain] and fr.attributes[:var_name].
source
Querying output
PALEOboxes.get_table — Method
get_table(output::PALEOmodel.AbstractOutputWriter, domainname::AbstractString) -> Table
-get_table(output::PALEOmodel.AbstractOutputWriter, varnames::Vector{<:AbstractString}) -> Table
Return a DataFrame with raw model output data for Domain domainname, or for Variables varnames
source
PALEOboxes.show_variables — Method
show_variables(output::PALEOmodel.AbstractOutputWriter; kwargs...) -> Table
+)
Initialize from a PALEOboxes::Model, optionally reserving memory for an assumed output dataset of nrecords.
The default for record_dim_name is :tmodel, for a sequence of records following the time evolution of the model.
source
PALEOmodel.OutputWriters.add_record! — Function
add_record!(output::PALEOmodel.AbstractOutputWriter, model, modeldata, rec_coord)
Add an output record for current state of model at record coordinate rec_coord. The usual case (set by initialize!) is that the record coordinate is model time tmodel.
source
Modifying output
PALEOboxes.add_field! — Method
add_field!(output::PALEOmodel.AbstractOutputWriter, fr::PALEOmodel.FieldRecord)
Add PALEOmodel.FieldRecord fr to output, with Domain and name defined by fr.attributes[:var_domain] and fr.attributes[:var_name].
source
Querying output
PALEOboxes.get_table — Method
get_table(output::PALEOmodel.AbstractOutputWriter, domainname::AbstractString) -> Table
+get_table(output::PALEOmodel.AbstractOutputWriter, varnames::Vector{<:AbstractString}) -> Table
Return a DataFrame with raw model output data for Domain domainname, or for Variables varnames
source
PALEOboxes.show_variables — Method
show_variables(output::PALEOmodel.AbstractOutputWriter; kwargs...) -> Table
 show_variables(output::PALEOmodel.AbstractOutputWriter, domainname::AbstractString; kwargs...) -> Table
Keywords:
  • attributes=[:units, :vfunction, :space, :field_data, :description]: Variable attributes to include
  • filter = attrb->true: function to filter by Variable attributes. Example: filter=attrb->attrb[:vfunction]!=PB.VF_Undefined to show state Variables and derivatives.
Examples:
julia> vscodedisplay(PB.show_variables(run.output))  # show all output Variables in VS code table viewer
-julia> vscodedisplay(PB.show_variables(run.output, ["atm.pCO2PAL", "fluxOceanBurial.flux_P_total"]))  # show subset of output Variables in VS code table viewer
source
PALEOboxes.has_variable — Method
has_variable(output::PALEOmodel.AbstractOutputWriter, varname::AbstractString)  -> Bool
True if model output contains Variable varname.
source
Accessing output data
PALEOmodel.get_array — Method
get_array(output::PALEOmodel.AbstractOutputWriter, varname::AbstractString [, allselectargs::NamedTuple]; kwargs...) -> FieldArray
-[deprecated] get_array(output::PALEOmodel.AbstractOutputWriter, varname::AbstractString; allselectargs...) -> FieldArray
Return a PALEOmodel.FieldArray containing data values and any attached coordinates.
Equivalent to PALEOmodel.get_array(PB.get_field(output, varname), allselectargs; kwargs), see PALEOmodel.get_array(fr::PALEOmodel.FieldRecord).
source
PALEOboxes.get_data — Method
get_data(output::PALEOmodel.AbstractOutputWriter, varname; records=nothing, kwargs...) -> values
Get Variable varname raw data array, optionally restricting to records.
Equivalent to PB.get_data(PB.get_field(output, varname); records, kwargs...), see PB.get_data(fr::PALEOmodel.FieldRecord).
source
PALEOboxes.get_data — Method
PB.get_data(fr::FieldRecord; records=nothing, squeeze_all_single_dims=true)
Get data records in raw format. Only recommended for variables with scalar data ie one value per record.
records may be nothing to get all records,  an Int to select a single record, or a range to select multiple records.
If squeeze_all_single_dims=true (the default), if each record represents a scalar (eg a PALEO Variable with Space PB.ScalarSpace, or a PB.CellSpace variable in a Domain with a single cell), then data is returned as a Vector. NB: if records is an Int, the single record requested is returned as a length-1 Vector.
Non-scalar data (eg a non-ScalarSpace variable from a Domain with > 1 cell)  is returned in internal format as a Vector-of-Vectors.
source
PALEOboxes.get_mesh — Method
get_mesh(output::PALEOmodel.AbstractOutputWriter, domainname::AbstractString) -> grid::Union{AbstractMesh, Nothing}
Return grid for output Domain domainname.
source
PALEOmodel.FieldRecord — Type
FieldRecord{FieldData <: AbstractData, Space <: AbstractSpace, V, N, Mesh <: AbstractMeshOrNothing, R}
-FieldRecord(dataset, f::PB.Field, attributes; [sizehint=nothing])
A series of records::R each containing values from a PALEOboxes.Field{FieldData, Space, N, V, Mesh}.
Access stored values:
Implementation
Storage in records::R is an internal format that may have:
  • An optimisation to store Fields with single values as a Vector of eltype(Field.values), cf Fields with array values are stored in records as a Vector of arrays.
  • A spatial cartesian grid stored as a linear Vector (eg mesh isa PALEOboxes.Grids.CartesianLinearGrid)
source
OutputMemory
PALEOmodel.OutputWriters.OutputMemory — Type
OutputMemory(; user_data=Dict{String, UserDataTypes}())
In-memory container for model output, organized by model Domains.
Implements the PALEOmodel.AbstractOutputWriter interface, with additional methods save_netcdf and load_netcdf! to save and load data.
Implementation
  • Field domains::Dict{String, OutputMemoryDomain} contains per-Domain model output.
  • Field user_data::Dict{String, UserDataTypes} contains optional user data NB:
    • available types are restricted to those that are compatible with NetCDF attribute types, ie Float64, Int64, String, Vector{Float64}, Vector{Int64}, Vector{String}
    • Vectors with a single element are read back from netcdf as scalars, see https://alexander-barth.github.io/NCDatasets.jl/dev/issues/#Corner-cases
source
PALEOmodel.OutputWriters.save_netcdf — Function
save_netcdf(output::OutputMemory, filename; kwargs...)
Save to filename in netcdf4 format (NB: filename must either have no extension or have extension .nc)
Notes on structure of netcdf output
  • Each PALEO Domain is written to a netcdf4 group. These can be read into a Python xarray using the group=<domainname> argument to open_dataset.
  • Isotope-valued variables (field_data = PB.IsotopeLinear) are written with an extra isotopelinear netCDF dimension, containing the variable total and delta.
  • Any '/' characters in PALEO variables are substited for '%' in the netcdf name.
Keyword arguments
  • check_ext::Bool = true: check that filename ends in ".nc"
source
PALEOmodel.OutputWriters.load_netcdf! — Function
load_netcdf!(output::OutputMemory, filename)
Load from filename in netCDF format, replacing any existing content in output. (NB: filename must either have no extension or have extension .nc).
Example
julia> output = PALEOmodel.OutputWriters.load_netcdf!(PALEOmodel.OutputWriters.OutputMemory(), "savedoutput.nc")
source
Field Array
FieldArray provides a generic array type with named dimensions PALEOboxes.NamedDimension each with optional coordinates for processing of model output.
PALEOmodel.FieldArray — Type
FieldArray
A generic xarray-like or  IRIS-like  n-dimensional Array with named dimensions and optional coordinates.
NB: this aims to be simple and generic, not efficient !!! Intended for representing model output, not for numerically-intensive calculations.
Fields
  • name::String: variable name
  • values::AbstractArray: n-dimensional Array of values
  • dims_coords::Vector{Pair{PALEOboxes.NamedDimension, Vector{PALEOmodel.FieldArray}}}: Names of dimensions with optional attached coordinates
  • attributes::Union{Nothing, Dict{Symbol, Any}}: variable attributes
source
PALEOmodel.FieldArray — Method
FieldArray(
+julia> vscodedisplay(PB.show_variables(run.output, ["atm.pCO2PAL", "fluxOceanBurial.flux_P_total"]))  # show subset of output Variables in VS code table viewer
source
PALEOboxes.has_variable — Method
has_variable(output::PALEOmodel.AbstractOutputWriter, varname::AbstractString)  -> Bool
True if model output contains Variable varname.
source
Accessing output data
PALEOmodel.get_array — Method
get_array(output::PALEOmodel.AbstractOutputWriter, varname::AbstractString [, allselectargs::NamedTuple]; kwargs...) -> FieldArray
+[deprecated] get_array(output::PALEOmodel.AbstractOutputWriter, varname::AbstractString; allselectargs...) -> FieldArray
Return a PALEOmodel.FieldArray containing data values and any attached coordinates.
Equivalent to PALEOmodel.get_array(PB.get_field(output, varname), allselectargs; kwargs), see PALEOmodel.get_array(fr::PALEOmodel.FieldRecord).
source
PALEOboxes.get_data — Method
get_data(output::PALEOmodel.AbstractOutputWriter, varname; records=nothing, kwargs...) -> values
Get Variable varname raw data array, optionally restricting to records.
Equivalent to PB.get_data(PB.get_field(output, varname); records, kwargs...), see PB.get_data(fr::PALEOmodel.FieldRecord).
source
PALEOboxes.get_data — Method
PB.get_data(fr::FieldRecord; records=nothing, squeeze_all_single_dims=true)
Get data records in raw format. Only recommended for variables with scalar data ie one value per record.
records may be nothing to get all records,  an Int to select a single record, or a range to select multiple records.
If squeeze_all_single_dims=true (the default), if each record represents a scalar (eg a PALEO Variable with Space PB.ScalarSpace, or a PB.CellSpace variable in a Domain with a single cell), then data is returned as a Vector. NB: if records is an Int, the single record requested is returned as a length-1 Vector.
Non-scalar data (eg a non-ScalarSpace variable from a Domain with > 1 cell)  is returned in internal format as a Vector-of-Vectors.
source
PALEOboxes.get_mesh — Method
get_mesh(output::PALEOmodel.AbstractOutputWriter, domainname::AbstractString) -> grid::Union{AbstractMesh, Nothing}
Return grid for output Domain domainname.
source
PALEOmodel.FieldRecord — Type
FieldRecord{FieldData <: AbstractData, Space <: AbstractSpace, V, N, Mesh <: AbstractMeshOrNothing, R}
+FieldRecord(dataset, f::PB.Field, attributes; [sizehint=nothing])
A series of records::R each containing values from a PALEOboxes.Field{FieldData, Space, N, V, Mesh}.
Access stored values:
Implementation
Storage in records::R is an internal format that may have:
  • An optimisation to store Fields with single values as a Vector of eltype(Field.values), cf Fields with array values are stored in records as a Vector of arrays.
  • A spatial cartesian grid stored as a linear Vector (eg mesh isa PALEOboxes.Grids.CartesianLinearGrid)
source
OutputMemory
PALEOmodel.OutputWriters.OutputMemory — Type
OutputMemory(; user_data=Dict{String, UserDataTypes}())
In-memory container for model output, organized by model Domains.
Implements the PALEOmodel.AbstractOutputWriter interface, with additional methods save_netcdf and load_netcdf! to save and load data.
Implementation
  • Field domains::Dict{String, OutputMemoryDomain} contains per-Domain model output.
  • Field user_data::Dict{String, UserDataTypes} contains optional user data NB:
    • available types are restricted to those that are compatible with NetCDF attribute types, ie Float64, Int64, String, Vector{Float64}, Vector{Int64}, Vector{String}
    • Vectors with a single element are read back from netcdf as scalars, see https://alexander-barth.github.io/NCDatasets.jl/dev/issues/#Corner-cases
source
PALEOmodel.OutputWriters.save_netcdf — Function
save_netcdf(output::OutputMemory, filename; kwargs...)
Save to filename in netcdf4 format (NB: filename must either have no extension or have extension .nc)
Notes on structure of netcdf output
  • Each PALEO Domain is written to a netcdf4 group. These can be read into a Python xarray using the group=<domainname> argument to open_dataset.
  • Isotope-valued variables (field_data = PB.IsotopeLinear) are written with an extra isotopelinear netCDF dimension, containing the variable total and delta.
  • Any '/' characters in PALEO variables are substited for '%' in the netcdf name.
Keyword arguments
  • check_ext::Bool = true: check that filename ends in ".nc"
source
PALEOmodel.OutputWriters.load_netcdf! — Function
load_netcdf!(output::OutputMemory, filename)
Load from filename in netCDF format, replacing any existing content in output. (NB: filename must either have no extension or have extension .nc).
Example
julia> output = PALEOmodel.OutputWriters.load_netcdf!(PALEOmodel.OutputWriters.OutputMemory(), "savedoutput.nc")
source
Field Array
FieldArray provides a generic array type with named dimensions PALEOboxes.NamedDimension each with optional coordinates for processing of model output.
PALEOmodel.FieldArray — Type
FieldArray
A generic xarray-like or  IRIS-like  n-dimensional Array with named dimensions and optional coordinates.
NB: this aims to be simple and generic, not efficient !!! Intended for representing model output, not for numerically-intensive calculations.
Fields
  • name::String: variable name
  • values::AbstractArray: n-dimensional Array of values
  • dims_coords::Vector{Pair{PALEOboxes.NamedDimension, Vector{PALEOmodel.FieldArray}}}: Names of dimensions with optional attached coordinates
  • attributes::Union{Nothing, Dict{Symbol, Any}}: variable attributes
source
PALEOmodel.FieldArray — Method
FieldArray(
     fr::FieldRecord; 
     expand_cartesian=true, 
     omit_recorddim_if_constant=true,
     lookup_coords=true,
     coords=nothing,
-)
Construct FieldArray from all records in fr::FieldRecord
Provides a view of internal storage format of FieldRecord as an n-dimensional Array.
Keyword arguments
  • expand_cartesian: (spatially resolved output using a CartesianLinearGrid only), true to expand compressed internal data (with spatial dimension cells) to a cartesian array  (with spatial dimensions eg lon, lat, zt)
  • omit_recorddim_if_constant: Specify whether to include multiple (identical) records and record dimension for constant variables  (with attribute is_constant = true). PALEO currently always stores these records in the input fr::FieldRecord;   set false include them in FieldArray output, true to omit duplicate records and record dimension from output.
  • lookup_coords: true to include coordinates, false to omit coordinates.
  • coords: replace the attached coordinates from the input fr::FieldRecord for one or more dimensions. Format is a Vector of Pairs of "<dim_name>"=>("<var_name1>", "<var_name2>", ...),  eg to replace an nD column atmosphere model default pressure coordinate with a z coordinate:
    coords=["cells"=>("atm.zmid", "atm.zlower", "atm.zupper")]
source
PALEOmodel.select — Method
select(fa::FieldArray [, allselectargs::NamedTuple] ; kwargs...) -> fa_out::FieldArray
Select a region of a FieldArray defined by allselectargs.
Selecting records and regions
allselectargs is a NamedTuple of form:
(<dimcoordname> = <filter>, <dimcoordname> = <filter>,  ... [, squeeze_all_single_dims=true])
where <dimcoordname> is of form:
  • <dimname>_isel to select by array indices: <filter> may then be a single Int to select a single index, or a range first:last to select a range of indices.
  • <coordname> to select by coordinate values using the coordinates attached to each dimension: <filter> may then be a single number to select a single index corresponding to the nearest value of the corresponding coordinate, or (first::Float64, last::Float64)  (a Tuple) to select a range starting at the index with the nearest value of fr.coords_record before first and ending at  the nearest index after last.
Available dimensions and coordinates <dimcoordname> depend on the FieldArray dimensions (as returned by get_dimensions, which will be a subset of grid spatial dimensions and Domain data dimensions) and corresponding attached coordinates (as returned by get_coordinates).
Some synonyms are defined for commonly used <dimnamecoordname>, these may require optional record_dim_name and mesh to be supplied in kwargs:
synonymsdimcoordnamecomment
records<record_dim_name>_iselrequires record_dim_name to be suppliied, usually tmodel
cells, cellcells_iselsubstituting named cells requires mesh to be supplied
column=<n>cells_isel = first:lastrequires mesh to be supplied, select range of cells corresponding to column n
NB: Dimensions corresponding to a selection for a single index or coordinate value are always squeezed out from the returned FieldArray. Optional argument squeeze_all_single_dims (default true) controls whether all output dimensions that contain a single index are squeezed out (including eg a selection for a range that results in a dimension with one index, or where the input FieldArray contains a dimension  with a single index).
Selection arguments used are optionally returned as strings in fa_out.attributes[:filter_records] and fa_out.attributes[:filter_region] for use in plot labelling etc.
Keyword arguments
  • record_dim_name=nothing: optionally supply name of record dimension as a String, to substitute for <dimcoordname> records  and to define dimension to use to populate filter_records attribute.
  • mesh=nothing: optionally supply a grid from PALEOboxes.Grids, to use for substituting cell names, looking up column indices.
  • add_attributes=true: true to transfer attributes from input fr::FieldRecord to output FieldArray   adding :filter_region and :filter_records, false to omit.
  • update_name=true: true to update output FieldArray name to add Domain name prefix, and a suffix generated from allselectargs   (NB: requires add_attributes=true), false to use name from input FieldRecord.
Limitations
  • it is only possible to select either a single slice or a contiguous range for each dimension, not a set of slices for a Vector of index or coordinate values.
source
PALEOmodel.get_array — Method
get_array(
+)
Construct FieldArray from all records in fr::FieldRecord
Provides a view of internal storage format of FieldRecord as an n-dimensional Array.
Keyword arguments
  • expand_cartesian: (spatially resolved output using a CartesianLinearGrid only), true to expand compressed internal data (with spatial dimension cells) to a cartesian array  (with spatial dimensions eg lon, lat, zt)
  • omit_recorddim_if_constant: Specify whether to include multiple (identical) records and record dimension for constant variables  (with attribute is_constant = true). PALEO currently always stores these records in the input fr::FieldRecord;   set false include them in FieldArray output, true to omit duplicate records and record dimension from output.
  • lookup_coords: true to include coordinates, false to omit coordinates.
  • coords: replace the attached coordinates from the input fr::FieldRecord for one or more dimensions. Format is a Vector of Pairs of "<dim_name>"=>("<var_name1>", "<var_name2>", ...),  eg to replace an nD column atmosphere model default pressure coordinate with a z coordinate:
    coords=["cells"=>("atm.zmid", "atm.zlower", "atm.zupper")]
source
PALEOmodel.select — Method
select(fa::FieldArray [, allselectargs::NamedTuple] ; kwargs...) -> fa_out::FieldArray
Select a region of a FieldArray defined by allselectargs.
Selecting records and regions
allselectargs is a NamedTuple of form:
(<dimcoordname> = <filter>, <dimcoordname> = <filter>,  ... [, squeeze_all_single_dims=true])
where <dimcoordname> is of form:
  • <dimname>_isel to select by array indices: <filter> may then be a single Int to select a single index, or a range first:last to select a range of indices.
  • <coordname> to select by coordinate values using the coordinates attached to each dimension: <filter> may then be a single number to select a single index corresponding to the nearest value of the corresponding coordinate, or (first::Float64, last::Float64)  (a Tuple) to select a range starting at the index with the nearest value of fr.coords_record before first and ending at  the nearest index after last.
Available dimensions and coordinates <dimcoordname> depend on the FieldArray dimensions (as returned by get_dimensions, which will be a subset of grid spatial dimensions and Domain data dimensions) and corresponding attached coordinates (as returned by get_coordinates).
Some synonyms are defined for commonly used <dimnamecoordname>, these may require optional record_dim_name and mesh to be supplied in kwargs:
synonymsdimcoordnamecomment
records<record_dim_name>_iselrequires record_dim_name to be suppliied, usually tmodel
cells, cellcells_iselsubstituting named cells requires mesh to be supplied
column=<n>cells_isel = first:lastrequires mesh to be supplied, select range of cells corresponding to column n
NB: Dimensions corresponding to a selection for a single index or coordinate value are always squeezed out from the returned FieldArray. Optional argument squeeze_all_single_dims (default true) controls whether all output dimensions that contain a single index are squeezed out (including eg a selection for a range that results in a dimension with one index, or where the input FieldArray contains a dimension  with a single index).
Selection arguments used are optionally returned as strings in fa_out.attributes[:filter_records] and fa_out.attributes[:filter_region] for use in plot labelling etc.
Keyword arguments
  • record_dim_name=nothing: optionally supply name of record dimension as a String, to substitute for <dimcoordname> records  and to define dimension to use to populate filter_records attribute.
  • mesh=nothing: optionally supply a grid from PALEOboxes.Grids, to use for substituting cell names, looking up column indices.
  • add_attributes=true: true to transfer attributes from input fr::FieldRecord to output FieldArray   adding :filter_region and :filter_records, false to omit.
  • update_name=true: true to update output FieldArray name to add Domain name prefix, and a suffix generated from allselectargs   (NB: requires add_attributes=true), false to use name from input FieldRecord.
Limitations
  • it is only possible to select either a single slice or a contiguous range for each dimension, not a set of slices for a Vector of index or coordinate values.
source
PALEOmodel.get_array — Method
get_array(
     fr::FieldRecord [, allselectargs::NamedTuple];
     coords=nothing,
     lookup_coords=true,
@@ -28,19 +28,19 @@
 )
  • select surface 2D array (dimension zt, index 1) from 3D output at nearest available    time to model time 10.0, expanding CartesianLinearGrid:
    get_array(fr, (zt_isel=1, tmodel=10.0, expand_cartesian=true))
  • select section 2D array at nearest index to longitude 200 degrees from 3D output at nearest available  time to model time 10.0, expanding CartesianLinearGrid:
    get_array(fr, (lon=200.0, tmodel=10.0, expand_cartesian=true))
  • get full data cube as used for netcdf output, omitting coordinates and attributes, retaining singleton dimensions, and omitting  record dimension if variable is a constant:
    get_array(
     fr, (expand_cartesian=true, squeeze_all_single_dims=false);
     lookup_coordinates=false, add_attributes=false, omit_recorddim_if_constant=true
-)
    • Limitations
    • it is only possible to select either a single slice or a contiguous range for each dimension, not a set of slices for a Vector of index or coordinate values.
    source
    Plotting output
    Plot recipes
    Plotting using the Julia Plots.jl package is implemented by plot recipes that enable plotting of PALEO data types using the plot command.
    The general principles are that:
    • Data is extracted from model output into FieldArrays with attached coordinates
    • Vector-valued arguments are "broadcast" to allow multiple line plot series to be overlaid in a single plot panel
    RecipesBase.apply_recipe — Method
    plot(output::AbstractOutputWriter, vars::Union{AbstractString, Vector{<:AbstractString}} [, selectargs::NamedTuple] [; coords::AbstractVector])
+)
    Limitations
    • it is only possible to select either a single slice or a contiguous range for each dimension, not a set of slices for a Vector of index or coordinate values.
    source
    Plotting output
    Plot recipes
    Plotting using the Julia Plots.jl package is implemented by plot recipes that enable plotting of PALEO data types using the plot command.
    The general principles are that:
    • Data is extracted from model output into FieldArrays with attached coordinates
    • Vector-valued arguments are "broadcast" to allow multiple line plot series to be overlaid in a single plot panel
    RecipesBase.apply_recipe — Method
    plot(output::AbstractOutputWriter, vars::Union{AbstractString, Vector{<:AbstractString}} [, selectargs::NamedTuple] [; coords::AbstractVector])
 heatmap(output::AbstractOutputWriter, var::AbstractString [, selectargs::NamedTuple] [; coords::AbstractVector])
-plot(outputs::Vector{<:AbstractOutputWriter}, vars::Union{AbstractString, Vector{<:AbstractString}} [, selectargs::NamedTuple] [; coords::AbstractVector])
    Plot recipe that calls PB.get_field(output, var), and passes on to plot(fr::FieldRecord, selectargs) (see RecipesBase.apply_recipe(::Dict{Symbol, Any}, fr::FieldRecord, selectargs::NamedTuple))
    Vector-valued outputs or vars are "broadcast" to create a plot series for each element. A labelprefix (index in outputs Vector) is added to identify each plot series produced.
    If var is of form <domain>.<name>.<structfield>, then sets the structfield keyword argument to take a single field from a struct Variable.
    Optional argument coords can be used to supply plot coordinates from Variable in output. Format is a Vector of Pairs of "coordname"=>("varname1", "var_name2", ...)
    Example: to replace a 1D column default pressure coordinate with a z coordinate:
    coords=["z"=>("atm.zmid", "atm.zlower", "atm.zupper")]
    source
    RecipesBase.apply_recipe — Method
    plot(fr::FieldRecord, selectargs::NamedTuple)
-heatmap(fr::FieldRecord, selectargs::NamedTuple)
    Plot recipe to plot a single FieldRecord
    Calls get_array(fr, selectargs) and passes on to plot(fa::FieldArray) (see RecipesBase.apply_recipe(::Dict{Symbol, Any}, fa::FieldArray)).
    Vector-valued fields in selectargs are "broadcasted" (generating a separate plot series for each combination)
    Optional argument coords can be used to supply plot coordinates from Variable in output. Format is a Vector of Pairs of "coordname"=>("varname1", "var_name2", ...)
    Example: to replace a 1D column default pressure coordinate with a z coordinate:
    coords=["z"=>("atm.zmid", "atm.zlower", "atm.zupper")]
    source
    RecipesBase.apply_recipe — Method
    plot(fa::FieldArray; kwargs...)
+plot(outputs::Vector{<:AbstractOutputWriter}, vars::Union{AbstractString, Vector{<:AbstractString}} [, selectargs::NamedTuple] [; coords::AbstractVector])
    Plot recipe that calls PB.get_field(output, var), and passes on to plot(fr::FieldRecord, selectargs) (see RecipesBase.apply_recipe(::Dict{Symbol, Any}, fr::FieldRecord, selectargs::NamedTuple))
    Vector-valued outputs or vars are "broadcast" to create a plot series for each element. A labelprefix (index in outputs Vector) is added to identify each plot series produced.
    If var is of form <domain>.<name>.<structfield>, then sets the structfield keyword argument to take a single field from a struct Variable.
    Optional argument coords can be used to supply plot coordinates from Variable in output. Format is a Vector of Pairs of "coordname"=>("varname1", "var_name2", ...)
    Example: to replace a 1D column default pressure coordinate with a z coordinate:
    coords=["z"=>("atm.zmid", "atm.zlower", "atm.zupper")]
    source
    RecipesBase.apply_recipe — Method
    plot(fr::FieldRecord, selectargs::NamedTuple)
+heatmap(fr::FieldRecord, selectargs::NamedTuple)
    Plot recipe to plot a single FieldRecord
    Calls get_array(fr, selectargs) and passes on to plot(fa::FieldArray) (see RecipesBase.apply_recipe(::Dict{Symbol, Any}, fa::FieldArray)).
    Vector-valued fields in selectargs are "broadcasted" (generating a separate plot series for each combination)
    Optional argument coords can be used to supply plot coordinates from Variable in output. Format is a Vector of Pairs of "coordname"=>("varname1", "var_name2", ...)
    Example: to replace a 1D column default pressure coordinate with a z coordinate:
    coords=["z"=>("atm.zmid", "atm.zlower", "atm.zupper")]
    source
    RecipesBase.apply_recipe — Method
    plot(fa::FieldArray; kwargs...)
 heatmap(fa::FieldArray; kwargs...)
-plot(fas::Vector{<:FieldArray}; kwargs...)
    Plot recipe that plots a single [FieldArray] or Vector of [FieldArray]. 
    If fa has a single dimension, this is suitable for a line-like plot, if two dimensions, a heatmap.
    If fas::Vector is supplied, this is "broadcast" generating one plot series for each element,  with the Vector index prepended to labelprefix to identify the plot series (unless overridden by labellist or labelattribute)
    Keywords
    • swap_xy::Bool=false: true to swap x and y axes 
    • mult_y_coord=1.0: workaround for bugs in Plots.jl heatmap yflip - multiply y coordinate by constant factor.
    • structfield::Union{Symbol, Nothing}=nothing: use field structfield from a struct-valued array.
    • map_values=PB.get_total: function to apply to y (for a 1D series) or z (for a 2D heatmap etc) before plotting 
    • labelprefix="": prefix for plot label.
    • labellist=[]: list of labels to override defaults
    • labelattribute=nothing: FieldArray attribute to use as label
    source
    Assembling multi-plot panels
    PALEOmodel.PlotPager — Type
    PlotPager(layout [, kwargs=NamedTuple()][; displayfunc=(plot, nplot)->display(plot)])
    Accumulate plots into subplots.
    layout is supplied to Plots.jl layout keyword, may be an Int or a Tuple (ny, nx), see https://docs.juliaplots.org/latest/.
    Optional argument kwargs::NamedTuple provides additional keyword arguments passed through to plot (eg (legend_background_color=nothing, ) to set all subplot legends to transparent backgrounds)
    Optional keyword argument displayfunc::(plot, nplot)->display(plot) provides the function used to display a screen of plots, where plot is a Plot object and nplot::Integer is the sequential number of this screen.
    Usage
    julia> pp = PlotPager((2,2))  # 4 panels per screen (2 down, 2 across)
+plot(fas::Vector{<:FieldArray}; kwargs...)
    Plot recipe that plots a single [FieldArray] or Vector of [FieldArray]. 
    If fa has a single dimension, this is suitable for a line-like plot, if two dimensions, a heatmap.
    If fas::Vector is supplied, this is "broadcast" generating one plot series for each element,  with the Vector index prepended to labelprefix to identify the plot series (unless overridden by labellist or labelattribute)
    Keywords
    • swap_xy::Bool=false: true to swap x and y axes 
    • mult_y_coord=1.0: workaround for bugs in Plots.jl heatmap yflip - multiply y coordinate by constant factor.
    • structfield::Union{Symbol, Nothing}=nothing: use field structfield from a struct-valued array.
    • map_values=PB.get_total: function to apply to y (for a 1D series) or z (for a 2D heatmap etc) before plotting 
    • labelprefix="": prefix for plot label.
    • labellist=[]: list of labels to override defaults
    • labelattribute=nothing: FieldArray attribute to use as label
    source
    Assembling multi-plot panels
    PALEOmodel.PlotPager — Type
    PlotPager(layout [, kwargs=NamedTuple()][; displayfunc=(plot, nplot)->display(plot)])
    Accumulate plots into subplots.
    layout is supplied to Plots.jl layout keyword, may be an Int or a Tuple (ny, nx), see https://docs.juliaplots.org/latest/.
    Optional argument kwargs::NamedTuple provides additional keyword arguments passed through to plot (eg (legend_background_color=nothing, ) to set all subplot legends to transparent backgrounds)
    Optional keyword argument displayfunc::(plot, nplot)->display(plot) provides the function used to display a screen of plots, where plot is a Plot object and nplot::Integer is the sequential number of this screen.
    Usage
    julia> pp = PlotPager((2,2))  # 4 panels per screen (2 down, 2 across)
 julia> pp(plot(1:3))  # Accumulate
 julia> pp(:skip, plot(1:4), plot(1:5), plot(1:6))  # add multiple panels in one command
 julia> pp(:newpage) # flush any partial screen and start new page (NB: always add this at end of a sequence!)
 
-julia> pp = PlotPager((2, 2); displayfunc=(plot, nplot)->savefig(plot, "plot_$nplot.png")) # save to file instead of default display
    Commands
    • pp(p): accumulate plot p
    • pp(:skip): leave a blank panel
    • pp(:newpage): fill with blank panels and start new page
    • pp(p1, p2, ...): multiple plots/commands in one call 
    source
    Analyze reaction network
    PALEOmodel.ReactionNetwork — Module
    ReactionNetwork
    Functions to analyze a PALEOboxes.Model or PALEOmodel output that contains a reaction network
    Compiles reaction stoichiometry and rate information from attributes attached to reaction rate variables:
    • rate_processname::String: a process name (eg "photolysis", "reaction", ...)
    • rate_species::Vector{String} names of reactant and product species
    • rate_stoichiometry::Vector{Float64} stoichiometry of reactant and product species
    source
    PALEOmodel.ReactionNetwork.get_ratetable — Function
    get_ratetable(model, domainname) -> DataFrame
-get_ratetable(output, domainname) -> DataFrame
    Get table of rate Variables and reactions
    Returns a DataFrame with columns :name, :rate_processname, :rate_species, :rate_stoichiometry
    source
    PALEOmodel.ReactionNetwork.get_all_species_ratevars — Function
    get_all_species_ratevars(model, domainname) -> OrderedDict(speciesname => [(stoich, ratevarname, processname), ...])
+julia> pp = PlotPager((2, 2); displayfunc=(plot, nplot)->savefig(plot, "plot_$nplot.png")) # save to file instead of default display
    Commands
    • pp(p): accumulate plot p
    • pp(:skip): leave a blank panel
    • pp(:newpage): fill with blank panels and start new page
    • pp(p1, p2, ...): multiple plots/commands in one call 
    source
    Analyze reaction network
    PALEOmodel.ReactionNetwork — Module
    ReactionNetwork
    Functions to analyze a PALEOboxes.Model or PALEOmodel output that contains a reaction network
    Compiles reaction stoichiometry and rate information from attributes attached to reaction rate variables:
    • rate_processname::String: a process name (eg "photolysis", "reaction", ...)
    • rate_species::Vector{String} names of reactant and product species
    • rate_stoichiometry::Vector{Float64} stoichiometry of reactant and product species
    source
    PALEOmodel.ReactionNetwork.get_ratetable — Function
    get_ratetable(model, domainname) -> DataFrame
+get_ratetable(output, domainname) -> DataFrame
    Get table of rate Variables and reactions
    Returns a DataFrame with columns :name, :rate_processname, :rate_species, :rate_stoichiometry
    source
    PALEOmodel.ReactionNetwork.get_all_species_ratevars — Function
    get_all_species_ratevars(model, domainname) -> OrderedDict(speciesname => [(stoich, ratevarname, processname), ...])
 get_all_species_ratevars(output, domainname) -> OrderedDict(speciesname => [(stoich, ratevarname, processname), ...])
-get_all_species_ratevars(ratetable::DataFrame) -> OrderedDict(speciesname => [(stoich, ratevarname, processname), ...])
    Get all species and contributing reaction rate Variable names as Dict of Tuples (stoich, ratevarname, processname) where ratevarname is the name of an output Variable with a reaction rate, stoich is the stoichiometry of that rate applied to species, and processname identifies the nature of the reaction.
    source
    PALEOmodel.ReactionNetwork.get_rates — Function
    get_rates(output, domainname [, outputrec] [, indices] [, scalefac] [, add_equations] [, ratetable_source]) -> DataFrame
    Get all reaction rates as column rate_total for domainname from output record outputrec (defaults to last time record), for subset of cells in indices (defaults to whole domain).
    Set optional ratetable_source = model to use with older output that doesn't include rate variable attributes.
    source
    PALEOmodel.ReactionNetwork.get_all_species_ratesummaries — Function
    get_all_species_ratesummaries(output, domainname [, outputrec] [, indices] [, scalefac] [, ratetable_source]) 
-    -> OrderedDict(speciesname => (source, sink, net, source_rxs, sink_rxs))
    Get source, sink, net rates and rates of source_rxs and sink_rxs  for all species in domainname from output record outputrec (defaults to last record),  cells in indices (defaults to whole domain),
    Optional scalefac to convert units, eg scalefac=1.90834e12 to convert mol m-2 yr-1 to molecule cm-2 s-1
    Set optional ratetable_source = model to use with older output that doesn't include rate variable attributes.
    source
    PALEOmodel.ReactionNetwork.show_ratesummaries — Function
    show_ratesummaries([io::IO = stdout], ratesummaries [; select_species=[]])
    Print per-species reaction rates from ratesummaries to output stream io (defaults to stdout),  optionally selecting species to print.
    Example
    ratesummaries = PALEOmodel.ReactionNetwork.get_all_species_ratesummaries(output, "atm")
-PALEOmodel.ReactionNetwork.show_ratesummaries(ratesummaries)
    source
    
+get_all_species_ratevars(ratetable::DataFrame) -> OrderedDict(speciesname => [(stoich, ratevarname, processname), ...])

    Get all species and contributing reaction rate Variable names as Dict of Tuples (stoich, ratevarname, processname) where ratevarname is the name of an output Variable with a reaction rate, stoich is the stoichiometry of that rate applied to species, and processname identifies the nature of the reaction.

    source
    PALEOmodel.ReactionNetwork.get_rates — Function
    get_rates(output, domainname [, outputrec] [, indices] [, scalefac] [, add_equations] [, ratetable_source]) -> DataFrame

    Get all reaction rates as column rate_total for domainname from output record outputrec (defaults to last time record), for subset of cells in indices (defaults to whole domain).

    Set optional ratetable_source = model to use with older output that doesn't include rate variable attributes.

    source
    PALEOmodel.ReactionNetwork.get_all_species_ratesummaries — Function
    get_all_species_ratesummaries(output, domainname [, outputrec] [, indices] [, scalefac] [, ratetable_source]) 
+    -> OrderedDict(speciesname => (source, sink, net, source_rxs, sink_rxs))

    Get source, sink, net rates and rates of source_rxs and sink_rxs for all species in domainname from output record outputrec (defaults to last record), cells in indices (defaults to whole domain),

    Optional scalefac to convert units, eg scalefac=1.90834e12 to convert mol m-2 yr-1 to molecule cm-2 s-1

    Set optional ratetable_source = model to use with older output that doesn't include rate variable attributes.

    source
    PALEOmodel.ReactionNetwork.show_ratesummaries — Function
    show_ratesummaries([io::IO = stdout], ratesummaries [; select_species=[]])

    Print per-species reaction rates from ratesummaries to output stream io (defaults to stdout), optionally selecting species to print.

    Example

    ratesummaries = PALEOmodel.ReactionNetwork.get_all_species_ratesummaries(output, "atm")
+PALEOmodel.ReactionNetwork.show_ratesummaries(ratesummaries)
    source
    diff --git a/dev/PALEOmodelSolvers/index.html b/dev/PALEOmodelSolvers/index.html index 326c9db..4f61380 100644 --- a/dev/PALEOmodelSolvers/index.html +++ b/dev/PALEOmodelSolvers/index.html @@ -2,34 +2,34 @@ PALEOmodel solvers · PALEOmodel Documentation
    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(
+)

    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.JacobianAD.jac_config_dae — Function
    jac_config_dae(
+)-> (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)
+) -> 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)
+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!
+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
    ca = ClampAll(minvalue, maxvalue)
-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(
+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
    ca = ClampAll(minvalue, maxvalue)
+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
@@ -43,7 +43,7 @@
     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[],
@@ -53,7 +53,7 @@
     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[],
@@ -62,40 +62,40 @@
 ) -> 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; 
+)

    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(
+)

    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;
+) -> (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")    
+                        [,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")    
+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.

    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(
+- `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,
@@ -109,12 +109,12 @@
     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(
+) -> 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
    +) -> 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
    PALEOmodel.SolverFunctions.JacODE_at_t — Type
    JacODE_at_t

    Function object to calculate ODE model Jacobian at t, eg to adapt to NLsolve interface

    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
    PALEOmodel.SolverFunctions.JacDAE — Type
    JacDAE

    Function object to calculate Jacobian in form required for SciML DAE solver

    odejac should be a JacODEForwardDiffDense or JacODEForwardDiffSparse

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

    Provides function signature:

    (jdae::JacDAE)(J, dsdt, s, p, γ, t)

    Calculates Jacobian J in the form γ*dG/d(dsdt) + dG/ds where γ is given by the solver

    source
    PALEOmodel.SolverFunctions.TotalForwardDiff — Type
    TotalForwardDiff

    Calculate Total variables, with function signature required by ForwardDiff

    Calling:

    set_t!(tfd::TotalForwardDiff, t)
+tfd(T, u)
    source
    PALEOmodel.SolverFunctions.ImplicitForwardDiffDense — Type
    ImplicitForwardDiffDense

    Calculate dT/dS required for a model containing implicit Total variables, using ForwardDiff and dense AD

    source
    PALEOmodel.SolverFunctions.ImplicitForwardDiffSparse — Type
    ImplicitForwardDiffSparse

    Calculate dT/dS required for a model containing implicit Total variables, using ForwardDiff and sparse AD with SparseDiffTools.forwarddiff_color_jacobian!

    source
    diff --git a/dev/References/index.html b/dev/References/index.html index 6a3d751..68f137c 100644 --- a/dev/References/index.html +++ b/dev/References/index.html @@ -1,2 +1,2 @@ -References · PALEOmodel Documentation
    GitHub
    +References · PALEOmodel Documentation
    GitHub
    diff --git a/dev/index.html b/dev/index.html index d639c45..ef8fdcf 100644 --- a/dev/index.html +++ b/dev/index.html @@ -1,2 +1,2 @@ -PALEOmodel.jl · PALEOmodel Documentation
    GitHub

    PALEOmodel.jl

    The PALEOmodel Julia package provides modules to create and solve a standalone PALEOboxes.Model, and to analyse output interactively from the Julia REPL. It implements:

    For examples and further information see the documentation for other PALEOtoolkit components and repositories:

    PALEO documentation follows the recommendations from https://documentation.divio.com/

    +PALEOmodel.jl · PALEOmodel Documentation
    GitHub

    PALEOmodel.jl

    The PALEOmodel Julia package provides modules to create and solve a standalone PALEOboxes.Model, and to analyse output interactively from the Julia REPL. It implements:

    For examples and further information see the documentation for other PALEOtoolkit components and repositories:

    PALEO documentation follows the recommendations from https://documentation.divio.com/

    diff --git a/dev/indexpage/index.html b/dev/indexpage/index.html index 28ac6b6..4399b7e 100644 --- a/dev/indexpage/index.html +++ b/dev/indexpage/index.html @@ -1,2 +1,2 @@ -Index · PALEOmodel Documentation
    GitHub

    Index

    +Index · PALEOmodel Documentation
    GitHub

    Index