diff --git a/dev/.documenter-siteinfo.json b/dev/.documenter-siteinfo.json index 3aa5b99..1ed255a 100644 --- a/dev/.documenter-siteinfo.json +++ b/dev/.documenter-siteinfo.json @@ -1 +1 @@ -{"documenter":{"julia_version":"1.6.7","generation_timestamp":"2024-02-09T16:54:17","documenter_version":"1.2.1"}} \ No newline at end of file +{"documenter":{"julia_version":"1.6.7","generation_timestamp":"2024-02-09T17:22:43","documenter_version":"1.2.1"}} \ No newline at end of file diff --git a/dev/HOWTOshowmodelandoutput/index.html b/dev/HOWTOshowmodelandoutput/index.html index e360765..5755ea6 100644 --- a/dev/HOWTOshowmodelandoutput/index.html +++ b/dev/HOWTOshowmodelandoutput/index.html @@ -29,4 +29,4 @@ swap_xy=true, xscale=:log10) # plot current value of O2, O3 vs height

Wavelength-dependent data

julia> plot(title="direct transmittance", output, ["atm.direct_trans"], (tmodel=1e12, column=1, cell=[1, 80]),
             ylabel="fraction", labelattribute=:filter_region) # plots vs wavelength

Here tmodel=1e12 selects the last model time output, and column=1, cell=[1, 80] selects the top and bottom cells within the first (only) 1D column. The labelattribute=:filter_region keyword argument is used to generate plot labels from the :filter_region FieldArray attribute, which contains the column and cell values used to select the spatial region.

Examples for a 3D GCM-based model

Visualisation of spatial output from the 3D GENIE transport-matrix example (PALEOdev.jl repository)

Horizontal slices across levels

julia> heatmap(run.output, "ocean.O2_conc", (tmodel=1e12, k=1), swap_xy=true)

Here k=1 selects a horizontal level corresponding to model grid cells with index k=1, which is the ocean surface in the GENIE grid.

Vertical section at constant longitude

julia> heatmap(run.output, "ocean.O2_conc", (tmodel=1e12, i=10), swap_xy=true, mult_y_coord=-1.0)

Here i=10 selects a section at longitude corresponding to model grid cells with index i=10.

Save and load output

Model output can be saved and loaded using the PALEOmodel.OutputWriters.save_netcdf and PALEOmodel.OutputWriters.load_netcdf! methods.

The PALEOmodel.OutputWriters.save_jld2 and PALEOmodel.OutputWriters.load_jld2! methods are still available as these files will not remain compatible with future PALEO versions.

Export output to a CSV file

To write Model output from a single Domain to a CSV file:

julia> import CSV
 julia> CSV.write("copse_land.csv", PB.get_table(run.output, "land")) # all Variables from land Domain
-julia> CSV.write("copse_atm.csv", PB.get_table(run.output, "atm")[!, [:tmodel, :pCO2atm, :pO2atm]]) # subset of Variables from atm Domain
+julia> CSV.write("copse_atm.csv", PB.get_table(run.output, "atm")[!, [:tmodel, :pCO2atm, :pO2atm]]) # subset of Variables from atm Domain diff --git a/dev/HOWTOsmallnegativevalues/index.html b/dev/HOWTOsmallnegativevalues/index.html index d765631..5a2e6ca 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 b34cbe3..e2cf57d 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_State.

+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_State.

diff --git a/dev/PALEOmodelOutput/index.html b/dev/PALEOmodelOutput/index.html index 3e5ebb2..0a8996c 100644 --- a/dev/PALEOmodelOutput/index.html +++ b/dev/PALEOmodelOutput/index.html @@ -1,32 +1,32 @@ -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 
     [;coords_record=:tmodel] [coords_record_units="year"]
-)
Initialize from a PALEOboxes::Model, reserving memory for an assumed output dataset of nrecords.
The default for coords_record 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, reserving memory for an assumed output dataset of nrecords.
The default for coords_record 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] [; coords::AbstractVector]) -> FieldArray
-get_array(output::PALEOmodel.AbstractOutputWriter, varname::AbstractString; allselectargs...) -> FieldArray
Return a PALEOmodel.FieldArray containing data values and any attached coordinates, for the  PALEOmodel.FieldRecord for varname, and records and spatial region defined by selectargs
If coords is not supplied, equivalent to PALEOmodel.get_array(PB.get_field(output, varname), allselectargs).
Optional argument coords can be used to supply plot coordinates from Variable in output, to replace any default coordinates. 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")]
NB: the coordinates will be generated by applying selectargs, so the supplied coordinate Variables must have the same dimensionality as vars.
source
PALEOboxes.get_data — Method
get_data(output::PALEOmodel.AbstractOutputWriter, varname; records=nothing) -> values
Get Variable varname raw data array, optionally restricting to records
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{D <: AbstractData, S <: AbstractSpace, V, N, M, R}
+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] [; coords::AbstractVector]) -> FieldArray
+get_array(output::PALEOmodel.AbstractOutputWriter, varname::AbstractString; allselectargs...) -> FieldArray
Return a PALEOmodel.FieldArray containing data values and any attached coordinates, for the  PALEOmodel.FieldRecord for varname, and records and spatial region defined by selectargs
If coords is not supplied, equivalent to PALEOmodel.get_array(PB.get_field(output, varname), allselectargs).
Optional argument coords can be used to supply plot coordinates from Variable in output, to replace any default coordinates. 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")]
NB: the coordinates will be generated by applying selectargs, so the supplied coordinate Variables must have the same dimensionality as vars.
source
PALEOboxes.get_data — Method
get_data(output::PALEOmodel.AbstractOutputWriter, varname; records=nothing) -> values
Get Variable varname raw data array, optionally restricting to records
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{D <: AbstractData, S <: AbstractSpace, V, N, M, R}
 FieldRecord(
     f::PB.Field{D, S, V, N, M}, attributes; 
     coords_record, 
     sizehint::Union{Nothing, Int}=nothing
-) -> fr
A series of records::R each containing the values from a PALEOboxes.Field{D, S, N, V, M}.
A coords_record may be attached to provide a coordinate (eg model time) corresponding to records.
Implementation
Fields with array values are stored in records as a Vector of arrays. Fields with single values (field_single_element true) are stored as a Vector of eltype(Field.values). 
source
OutputMemory
PALEOmodel.OutputWriters.OutputMemoryDomain — Type
OutputMemoryDomain
In-memory model output, for one model Domain.
Includes an additional coords_record (usually :tmodel, when storing output vs time).
Implementation
data::DataFrame contains columns of same type as FieldRecord.records for each Variable.
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"
  • add_coordinates::Bool = false: true to attempt to add CF convention coords to variables (experimental, doesn't look that useful)
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
PALEOmodel.OutputWriters.load_jld2! — Function
load_jld2!(output::OutputMemory, filename)
Load from filename in JLD2 format, replacing any existing content in output. (NB: filename must either have no extension or have extension .jld2).
Example
julia> output = PALEOmodel.OutputWriters.load_jld2!(PALEOmodel.OutputWriters.OutputMemory(), "savedoutput.jld2")
source
Field Array
FieldArray provides a generic array type with named dimensions PALEOboxes.NamedDimension and optional coordinates PALEOboxes.FixedCoord for processing of model output.
PALEOmodel.FieldArray — Type
FieldArray
A generic xarray-like or  IRIS-like  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.
source
PALEOmodel.get_array — Method
get_array(f::Field [, selectargs::NamedTuple]; [attributes=nothing]) -> FieldArray
Return a FieldArray containing f::Field data values and any attached coordinates, for the spatial region defined by selectargs.
Available selectargs depend on the grid f.mesh, and  are passed to PB.Grids.get_region.
attributes (if present) are added to FieldArray
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])
+) -> fr
A series of records::R each containing the values from a PALEOboxes.Field{D, S, N, V, M}.
A coords_record may be attached to provide a coordinate (eg model time) corresponding to records.
Implementation
Fields with array values are stored in records as a Vector of arrays. Fields with single values (field_single_element true) are stored as a Vector of eltype(Field.values). 
source
OutputMemory
PALEOmodel.OutputWriters.OutputMemoryDomain — Type
OutputMemoryDomain
In-memory model output, for one model Domain.
Includes an additional coords_record (usually :tmodel, when storing output vs time).
Implementation
data::DataFrame contains columns of same type as FieldRecord.records for each Variable.
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"
  • add_coordinates::Bool = false: true to attempt to add CF convention coords to variables (experimental, doesn't look that useful)
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
PALEOmodel.OutputWriters.load_jld2! — Function
load_jld2!(output::OutputMemory, filename)
Load from filename in JLD2 format, replacing any existing content in output. (NB: filename must either have no extension or have extension .jld2).
Example
julia> output = PALEOmodel.OutputWriters.load_jld2!(PALEOmodel.OutputWriters.OutputMemory(), "savedoutput.jld2")
source
Field Array
FieldArray provides a generic array type with named dimensions PALEOboxes.NamedDimension and optional coordinates PALEOboxes.FixedCoord for processing of model output.
PALEOmodel.FieldArray — Type
FieldArray
A generic xarray-like or  IRIS-like  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.
source
PALEOmodel.get_array — Method
get_array(f::Field [, selectargs::NamedTuple]; [attributes=nothing]) -> FieldArray
Return a FieldArray containing f::Field data values and any attached coordinates, for the spatial region defined by selectargs.
Available selectargs depend on the grid f.mesh, and  are passed to PB.Grids.get_region.
attributes (if present) are added to FieldArray
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(modeldata::AbstractModelData, vars::Union{AbstractString, Vector{<:AbstractString}} [, selectargs::NamedTuple] [; coords::AbstractVector])
-heatmap(modeldata::AbstractModelData, var::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")]
NB: the coordinates will be generated by applying selectargs, so the supplied coordinate Variables must have the same dimensionality as vars.
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_records can be used to supply plot coordinates from FieldRecords. Format is a Vector of Pairs of "coordname"=>(cr1::FieldRecord, cr2::FieldRecord, ...) Example:     coordsrecords=["z"=>(zmid::FieldRecord, zlower::FieldRecord, zupper::FieldRecord)] to replace a 1D column default pressure coordinate with a z coordinate.   NB: the coordinates will be generated by applying selectargs, so the supplied coordinate FieldRecords must have the same dimensionality as fr.
source
RecipesBase.apply_recipe — Method
plot(fa::FieldArray; kwargs...)
+heatmap(modeldata::AbstractModelData, var::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")]
NB: the coordinates will be generated by applying selectargs, so the supplied coordinate Variables must have the same dimensionality as vars.
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_records can be used to supply plot coordinates from FieldRecords. Format is a Vector of Pairs of "coordname"=>(cr1::FieldRecord, cr2::FieldRecord, ...) Example:     coordsrecords=["z"=>(zmid::FieldRecord, zlower::FieldRecord, zupper::FieldRecord)] to replace a 1D column default pressure coordinate with a z coordinate.   NB: the coordinates will be generated by applying selectargs, so the supplied coordinate FieldRecords must have the same dimensionality as fr.
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(model, 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(model, output, "atm")
+PALEOmodel.ReactionNetwork.show_ratesummaries(ratesummaries)
source
diff --git a/dev/PALEOmodelSolvers/index.html b/dev/PALEOmodelSolvers/index.html index 81411e2..d187e64 100644 --- a/dev/PALEOmodelSolvers/index.html +++ b/dev/PALEOmodelSolvers/index.html @@ -1,23 +1,23 @@ -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 with threadsafe=true (to create Atomic Variables for those Variables with attribute :atomic==true), and supplying method_barrier (a thread barrier to add to ReactionMethod dispatch lists between dependency-free groups)

Keyword summary

  • pickup_output=nothing: OutputWriter with pickup data to initialise from
  • 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
  • threadsafe=false: true to create thread safe Atomic Variables where Variable attribute :atomic==true
  • method_barrier=nothing: thread barrier to add to dispatch lists if threadsafe==true
  • 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.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:

and then

Keywords

As integrate, with defaults:

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

Low level 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.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(
+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 with threadsafe=true (to create Atomic Variables for those Variables with attribute :atomic==true), and supplying method_barrier (a thread barrier to add to ReactionMethod dispatch lists between dependency-free groups)
Keyword summary
  • pickup_output=nothing: OutputWriter with pickup data to initialise from
  • 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
  • threadsafe=false: true to create thread safe Atomic Variables where Variable attribute :atomic==true
  • method_barrier=nothing: thread barrier to add to dispatch lists if threadsafe==true
  • 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.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:
and then 
Keywords
As integrate, with defaults:
  • alg=Sundials.IDA() (integrateDAE)
  • alg=Sundials.IDA(linear_solver=:KLU) (integrateDAEForwardDiff)
source
Low level 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.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
Create (inconsistent) initial_deriv for a DAE problem: ODE variables are consistent, DAE variables set to zero  ie rely on DAE solver to find them
source
PALEOmodel.ODE.print_sol_stats — Function
print_sol_stats(sol::SciMLBase.ODESolution)
+) -> initial_deriv
Create (inconsistent) initial_deriv for a DAE problem: ODE variables are consistent, DAE variables set to zero  ie rely on DAE solver to find them
source
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 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.NonlinearSolution)
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.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
  • tss_output=Float64[]: Vector of model 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
  • max_iter=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: - 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
  • tss_output=Float64[]: Vector of model 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
  • max_iter=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: - 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
@@ -31,7 +31,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[],
@@ -41,7 +41,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[],
@@ -50,51 +50,51 @@
 ) -> 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] [,tss_output] [,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   (Domains, indices within spatial Domains) and Variables, with:     - ODE paired stateexplicit (S) and stateexplicit_deriv (dS/dt), where dS/dt = F(S).     - Implicit-ODE paired total (T) and total_deriv (dT/dt), where dT(S)/dt = F(T(S)) with         total a function of explicit and implicit state Variables statexplicit and state (S).     - Algebraic constraints (C), where C(S) = 0 with C a function of explicit and implicit state Variables statexplicit and state (S),
The number of total + number of constraint Variables must equal the number of implicit state Variables.
Optional access methods provide an ODE/DAE solver view with composite statevar and statevar_sms,where:
- `statevar` is a concatenation of `stateexplicit` 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; 
     solver_view=modeldata.solver_view_all,
     dispatchlists=modeldata.dispatchlists_all
-) -> f::ModelODE
Function object to calculate model time derivative and adapt to SciML ODE solver interface
Call as f(du,u, p, t)
source
PALEOmodel.SolverFunctions.JacODEForwardDiffDense — Type
JacODEForwardDiffDense(
     modeldata; 
     solver_view=modeldata.solver_view_all,
     dispatchlists=modeldata.dispatchlists_all,
     du_template, 
     jacconf,
-) -> jac::JacODEForwardDiffDense
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
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

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
diff --git a/dev/References/index.html b/dev/References/index.html index 0f29aee..7994b20 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 b5da889..4901ce0 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 396d150..d3dc72b 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