Reservoirs and fluxes
These examples illustrate how to code a PALEO Reaction and configure a model, working through a minimal example of first order decay or transformation of a scalar biogeochemical reservoir $A$ into a second reservoir $B$ via a flux $F$, with equations:
\[F = \kappa A\]
\[\frac{dA}{dt} = - F \\\]
\[\frac{dB}{dt} = F \\\]
See Displaying model configuration and output from the Julia REPL for more information on analysing model output.
Example 1 A minimal self-contained PALEO reaction
This is verbose as we have to (re)implement Variable setup and initialisation as well as the biogeochemical reaction of interest, but illustrates the structure of a PALEO model.
This verbose approach is not usually required, it is usually simpler to use predefined PALEO Reservoirs as described below, Example 2 Using a PALEO Reservoir
Additional code files
The code to implement a self-contained PALEO reaction is in file examples/reservoirs/reactions_ex1.jl, and needs to provide three methods for Variable setup, initialization at the start of the main loop, as well as the actual main loop do method. Variables are labelled as state Variables and derivatives by setting the :vfunction attribute to VF_StateExplicit and VF_Deriv.
module Example1
import PALEOboxes as PB
"""
ReactionExample1
Minimal but verbose example with a single state variable, first order decay.
This uses only low-level operations in order to provide a standalone example
that demonstrates all the steps needed for a functioning PALEO model.
See Example 2 for how to implement this functionality with much less boilerplate code.
"""
Base.@kwdef mutable struct ReactionExample1{P} <: PB.AbstractReaction
base::PB.ReactionBase
pars::P = PB.ParametersTuple(
PB.ParDouble("kappa", 1.0, units="yr-1", description="first order decay constant"),
)
end
function PB.register_methods!(rj::ReactionExample1)
# Variables are labelled as state Variables and derivatives by setting the
# :vfunction attribute to VF_StateExplicit and VF_Deriv.
# also need to set :field_data
A = PB.VarStateExplicitScalar("A", "mol", "reservoir for species A",
attributes=(:field_data=>PB.ScalarData,))
# this is equivalent to VarDepScalar with attribute :vfunction=>PB.VF_Deriv
# A = PB.VarDepScalar("A", "mol", "reservoir for species A",
# attributes=(:vfunction=>PB.VF_StateExplicit, :field_data=>PB.ScalarData))
A_sms = PB.VarDerivScalar("A_sms", "mol yr-1", "reservoir A source - sink",
attributes=(:vfunction=>PB.VF_Deriv, :field_data=>PB.ScalarData))
# this is equivalent to VarContribScalar with attribute :vfunction=>PB.VF_Deriv
# A_sms = PB.VarContribScalar("A_sms", "mol yr-1", "reservoir A source - sink",
# attributes=(:vfunction=>PB.VF_Deriv, :field_data=>PB.ScalarData))
# Provide a Property decay_flux as diagnostic output
decay_flux = PB.VarPropScalar("decay_flux", "mol yr-1", "decay flux from reservoir A")
PB.add_method_setup!(rj, setup_example1, (PB.VarList_namedtuple([A]),) )
PB.add_method_initialize!(rj, initialize_example1, (PB.VarList_namedtuple([A_sms]),) )
PB.add_method_do!(rj, do_example1, (PB.VarList_namedtuple([A, A_sms, decay_flux]), ) )
return nothing
end
# setup method, called at model startup
# NB: this is called three times, with attribute_name indicating the action or value that should be set:
# :setup - any non-state-Variable initialisation (not used here)
# :norm_value - Variable normalisation, usually read from the :norm_value Variable attribute
# :initial_value - Variable initial value, usually read from the :initial_value attribute
function setup_example1(m::PB.ReactionMethod, (varsdata, ), cellrange::PB.AbstractCellRange, attribute_name)
attribute_name in (:norm_value, :initial_value) || return
var_A = PB.get_variable(m, "A") # get the Variable as supplied to add_method_setup!
value = PB.get_attribute(var_A, attribute_name) # read attribute
@info "setup_example1: setting A[] to $value read from from attribute $attribute_name"
varsdata.A[] = value
return nothing
end
# initialize method, called at start of each main loop timestep
function initialize_example1(m::PB.ReactionMethod, (varsdata, ), cellrange::PB.AbstractCellRange, _)
varsdata.A_sms[] = 0.0
#varsdata.A[] = 100.0
return nothing
end
# do method, called each main loop timestep
function do_example1(m::PB.ReactionMethod, pars, (varsdata, ), cellrange::PB.AbstractCellRange, deltat)
# mol yr-1 yr-1 mol
varsdata.decay_flux[] = pars.kappa[] * varsdata.A[]
varsdata.A_sms[] -= varsdata.decay_flux[]
return nothing
end
end # moduleyaml configuration file
The model configuration (file examples/reservoirs/config_ex1.yaml) contains just a single Reaction:
############################################
# Example 1
# First order decay of a scalar variable
###########################################
example1:
domains:
global:
# scalar domain
reactions:
Adecay:
class: ReactionExample1
parameters:
kappa: 0.5
variable_attributes:
A:initial_value: 10.0
A:norm_value: 10.0
Run script
The script to run the model (file examples/reservoirs/run_ex1.yaml) contains:
import PALEOboxes as PB
import PALEOmodel
using Plots
include(joinpath(@__DIR__, "reactions_ex1.jl"))
#####################################################
# Create model
#######################################################
model = PB.create_model_from_config(joinpath(@__DIR__, "config_ex1.yaml"), "example1")
#########################################################
# Initialize
##########################################################
initial_state, modeldata = PALEOmodel.initialize!(model)
#####################################################################
# Optional: call ODE function to check derivative
#######################################################################
initial_deriv = similar(initial_state)
PALEOmodel.ODE.ModelODE(modeldata)(initial_deriv, initial_state , nothing, 0.0)
println("initial_state: ", initial_state)
println("initial_deriv: ", initial_deriv)
#################################################################
# Integrate vs time
##################################################################
# create a Run object to hold model and output
paleorun = PALEOmodel.Run(model=model, output = PALEOmodel.OutputWriters.OutputMemory())
println("integrate, ODE")
# first run is slow as it includes JIT time
@time PALEOmodel.ODE.integrate(
paleorun, initial_state, modeldata, (0.0, 10.0),
solvekwargs=(
reltol=1e-5,
# saveat=0.1, # save output every 0.1 yr see https://diffeq.sciml.ai/dev/basics/common_solver_opts/
)
);
########################################
# Plot output
########################################
display(plot(paleorun.output, "global.A"))
display(plot(paleorun.output, "global.decay_flux"))And produces output:
┌ Info:
│ ================================================================================
│ create_model_from_config: configmodel example1
│ config_file: /home/runner/work/PALEOtutorials.jl/PALEOtutorials.jl/examples/reservoirs/config_ex1.yaml
└ ================================================================================
[ Info: Model.parameters:
┌ Info:
│ ================================================================================
│ creating Domains
└ ================================================================================
[ Info: generated Reaction catalog with 74 Reactions
┌ Info:
│ ================================================================================
│ creating domain 'global' ID=1
└ ================================================================================
┌ Info: create_reaction_from_config: global.Adecay classname ReactionExample1
└ set parameters: [config.yaml] kappa =0.5
┌ Info:
│ ================================================================================
│ set_model_geometry
└ ================================================================================
┌ Info:
│ ================================================================================
│ register_reaction_methods!
└ ================================================================================
┌ Info:
│ ================================================================================
│ link_variables: first pass
└ ================================================================================
┌ Info:
│ ================================================================================
│ link_variables: register_reaction_dynamic_methods and configure variables
└ ================================================================================
┌ Info: _configure_variables: ReactionExample1 global.Adecay variable_attributes:
│ set attribute: A :initial_value = 10.0
└ set attribute: A :norm_value = 10.0
┌ Info:
│ ================================================================================
│ link_variables: second pass:
└ ================================================================================
┌ Info:
│ ================================================================================
│ link_variables! unlinked variables:
└ ================================================================================
┌ Info:
│ ================================================================================
│ create_model_from_config: done
└ ================================================================================
┌ Info:
│ ================================================================================
│ PALEOmodel.initialize! start
└ ================================================================================
┌ Info:
│ ================================================================================
│ allocate_variables! (modeldata arrays_idx=1)
└ ================================================================================
[ Info: Domain global data dimensions PALEOboxes.NamedDimension[] allocating 3 variables (hostdep=nothing)
[ Info: set_default_solver_view:
┌ Info: SolverView:
│ host-dependent Variables:
│ ------------------------------------------------------------------------------------------------------------------------------------------------
│ Domain operatorID VF_StateExplicit VF_Total VF_Constraint VF_StateTotal VF_State VF_Undefined
│ global 0 1 0 0 0 0 0
│ ------------------------------------------------------------------------------------------------------------------------------------------------
│ Total - 1 0 0 0 0 0
│ ------------------------------------------------------------------------------------------------------------------------------------------------
│ n_state_vars 1 (stateexplicit 1 + statetotal 0 + state 0)
│ n_equations 1 (stateexplicit 1 + total 0 + constraint 0)
│ including all host-dependent non-state Variables
└ host-dependent non-state Variables (:vfunction PB.VF_Undefined): Any[]
┌ Info:
│ ================================================================================
│ initialize_reactiondata! (modeldata arrays_indices=1:1)
└ ================================================================================
┌ Info:
│ ================================================================================
│ dispatch_setup :setup
└ ================================================================================
┌ Info:
│ ================================================================================
│ dispatch_setup :norm_value
└ ================================================================================
[ Info: setup_example1: setting A[] to 10.0 read from from attribute norm_value
┌ Info:
│ ================================================================================
│ dispatch_setup :initial_value
└ ================================================================================
[ Info: setup_example1: setting A[] to 10.0 read from from attribute initial_value
┌ Info:
│ ================================================================================
│ PALEOmodel.initialize! done
└ ================================================================================
initial_state: [10.0]
initial_deriv: [-5.0]
integrate, ODE
┌ Info:
│ ================================================================================
│ PALEOmodel.ODE.integrate:
│ tspan: (0.0, 10.0)
│ algorithm: Sundials.CVODE_BDF{:Newton, :Dense, Nothing, Nothing}(0, 0, 0, false, 10, 5, 7, 3, 10, nothing, nothing, 0)
│ Jacobian: NoJacobian
│ steadystate: false
│ using 1 BLAS threads
└ ================================================================================
[ Info: ODEfunction: using Jacobian NoJacobian
[ Info: jac_config_ode: jac_ad=NoJacobian
0.580970 seconds (1.13 M allocations: 78.937 MiB, 5.24% gc time, 99.93% compilation time)
┌ Info:
│ ================================================================================
│ print_sol_stats:
│ retcode=Success
│ alg=Sundials.CVODE_BDF{:Newton, :Dense, Nothing, Nothing}(0, 0, 0, false, 10, 5, 7, 3, 10, nothing, nothing, 0)
│ stats=SciMLBase.DEStats
│ Number of function 1 evaluations: 92
│ Number of function 2 evaluations: 0
│ Number of W matrix evaluations: 11
│ Number of linear solves: 0
│ Number of Jacobians created: 2
│ Number of nonlinear solver iterations: 89
│ Number of nonlinear solver convergence failures: 0
│ Number of fixed-point solver iterations: 0
│ Number of fixed-point solver convergence failures: 0
│ Number of rootfind condition calls: 0
│ Number of accepted steps: 72
│ Number of rejected steps: 2
│ length(sol.t) 75
│ size(sol) (1, 75)
└ ================================================================================
[ Info: ODE.calc_output_sol!: 75 records
[ Info: ODE.calc_output_sol!: done
┌ Info:
│ ================================================================================
│ PALEOmodel.ODE.integrate: done
└ ================================================================================
1.211437 seconds (1.77 M allocations: 122.745 MiB, 3.65% gc time, 99.70% compilation time)
Displaying model structure and variables
All metadata for model variables can be displayed with PB.show_variables:
show(PB.show_variables(model), allcols=true, allrows=true) # display in REPL
# vscodedisplay(PB.show_variables(model)) # more convenient when using VS code3×9 DataFrame
Row │ domain name type units vfunction space data_dims field_data description
│ String String DataType String Variable… DataType Tuple{} DataType String
─────┼───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
1 │ global A VariableDomPropDep mol VF_StateExplicit ScalarSpace () ScalarData reservoir for species A
2 │ global A_sms VariableDomContribTarget mol yr-1 VF_Deriv ScalarSpace () ScalarData reservoir A source - sink
3 │ global decay_flux VariableDomPropDep mol yr-1 VF_Undefined ScalarSpace () ScalarData decay flux from reservoir AThe type column shows the two pairings of VariableReactions linked to VariableDomains:
- Reaction Property and Dependency Variables, linked to a
VariableDomPropDep. These are used to represent a quantity calculated in one Reaction (or provided by the numerical solver, in the case of "global.A") that is then used by other Reactions. - Reaction Target and Contributor Variables, linked to a
VariableDomContribTarget. These are used to represent a flux-like quantity, with one Reaction (or the numerical solver, in the case of "global.A_sms") definining the Target and multiple Reactions adding contributions.
Variable links for an individual VariableDomain can be displayed using PB.show_links, where the linked variables are shown as "<domain name>.<reaction name>.<method name>.<local name>":
PB.show_links(model, "global.decay_flux") PALEOboxes.VariableDomPropDep "global.decay_flux" links:
property: "global.Adecay.do_example1.decay_flux"
dependencies: String[](demonstrating that global.decay_flux is set by the do_example1 method of the Reaction named Adecay in the yaml config file)
PB.show_links(model, "global.A") PALEOboxes.VariableDomPropDep "global.A" links:
property: nothing
dependencies: ["global.Adecay.setup_example1.A", "global.Adecay.do_example1.A"](demonstrating that the "global.A" state Variable has both the do_setup_example1 and do_example1 methods of the Reaction named Adecay in the yaml config file, which respectively set the initial value at model start, and then read the value at each timestep)
PB.show_variables with an additional showlinks=true argument will also show variable links:
show(PB.show_variables(model; showlinks=true), allcols=true, allrows=true) # display in REPL
# vscodedisplay(PB.show_variables(model; showlinks=true)) # more convenient when using VS code3×13 DataFrame
Row │ domain name type units vfunction space data_dims field_data description property dependencies target contributors
│ String String DataType String Variable… DataType Tuple{} DataType String Array…? Array…? Missing Array…?
─────┼─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
1 │ global A VariableDomPropDep mol VF_StateExplicit ScalarSpace () ScalarData reservoir for species A missing ["global.Adecay.setup_example1.A… missing missing
2 │ global A_sms VariableDomContribTarget mol yr-1 VF_Deriv ScalarSpace () ScalarData reservoir A source - sink missing missing missing ["global.Adecay.initialize_examp…
3 │ global decay_flux VariableDomPropDep mol yr-1 VF_Undefined ScalarSpace () ScalarData decay flux from reservoir A ["global.Adecay.do_example1.deca… missing missing missingwhere the property, dependencies, target and contributor columns show the VariableReactions (defined by ReactionMethods) that link to each VariableDomain.
Analysing output
Numerical values of variables at each timestep can be displayed with PB.get_table:
show(PB.get_table(paleorun.output, "global"), allcols=true, allrows=true) # display to REPL
# vscodedisplay(PB.get_table(paleorun.output, "global")) # more convenient when using VS code75×4 DataFrame
Row │ tmodel A A_sms decay_flux
│ Float64 Float64 Float64 Float64
─────┼─────────────────────────────────────────────────
1 │ 0.0 10.0 -5.0 5.0
2 │ 0.00449444 9.97758 -4.98879 4.98879
3 │ 0.00898888 9.95521 -4.9776 4.9776
4 │ 0.0201318 9.89995 -4.94997 4.94997
5 │ 0.0312747 9.84496 -4.92248 4.92248
6 │ 0.0424176 9.79027 -4.89513 4.89513
7 │ 0.0593185 9.70789 -4.85394 4.85394
8 │ 0.0890848 9.56447 -4.78223 4.78223
9 │ 0.174343 9.16517 -4.58259 4.58259
10 │ 0.21629 8.97491 -4.48746 4.48746
11 │ 0.258237 8.78861 -4.3943 4.3943
12 │ 0.300184 8.60617 -4.30309 4.30309
13 │ 0.371027 8.30661 -4.15331 4.15331
14 │ 0.441871 8.0175 -4.00875 4.00875
15 │ 0.512715 7.73847 -3.86924 3.86924
16 │ 0.583558 7.46916 -3.73458 3.73458
17 │ 0.743607 6.89479 -3.4474 3.4474
18 │ 0.903655 6.36462 -3.18231 3.18231
19 │ 1.0637 5.87522 -2.93761 2.93761
20 │ 1.22375 5.42345 -2.71172 2.71172
21 │ 1.3838 5.00641 -2.5032 2.5032
22 │ 1.54385 4.62144 -2.31072 2.31072
23 │ 1.7039 4.26607 -2.13304 2.13304
24 │ 1.86395 3.93803 -1.96901 1.96901
25 │ 2.02399 3.63521 -1.81761 1.81761
26 │ 2.18404 3.35568 -1.67784 1.67784
27 │ 2.34409 3.09764 -1.54882 1.54882
28 │ 2.50414 2.85945 -1.42972 1.42972
29 │ 2.66419 2.63957 -1.31978 1.31978
30 │ 2.82424 2.4366 -1.2183 1.2183
31 │ 2.98428 2.24923 -1.12462 1.12462
32 │ 3.14433 2.07628 -1.03814 1.03814
33 │ 3.30438 1.91662 -0.958311 0.958311
34 │ 3.46443 1.76924 -0.884621 0.884621
35 │ 3.62448 1.6332 -0.816598 0.816598
36 │ 3.78453 1.50761 -0.753805 0.753805
37 │ 3.94457 1.39168 -0.695841 0.695841
38 │ 4.10462 1.28467 -0.642334 0.642334
39 │ 4.26467 1.18588 -0.592941 0.592941
40 │ 4.42472 1.09469 -0.547346 0.547346
41 │ 4.58477 1.01052 -0.505258 0.505258
42 │ 4.74482 0.932812 -0.466406 0.466406
43 │ 4.90486 0.861083 -0.430541 0.430541
44 │ 5.06491 0.794869 -0.397435 0.397435
45 │ 5.22496 0.733747 -0.366874 0.366874
46 │ 5.38501 0.677325 -0.338663 0.338663
47 │ 5.54506 0.625242 -0.312621 0.312621
48 │ 5.70511 0.577164 -0.288582 0.288582
49 │ 5.86515 0.532782 -0.266391 0.266391
50 │ 6.0252 0.491814 -0.245907 0.245907
51 │ 6.18525 0.453995 -0.226998 0.226998
52 │ 6.3453 0.419085 -0.209543 0.209543
53 │ 6.50535 0.386859 -0.19343 0.19343
54 │ 6.6654 0.357112 -0.178556 0.178556
55 │ 6.82544 0.329651 -0.164826 0.164826
56 │ 6.98549 0.304303 -0.152151 0.152151
57 │ 7.14554 0.280903 -0.140452 0.140452
58 │ 7.30559 0.259303 -0.129651 0.129651
59 │ 7.46564 0.239364 -0.119682 0.119682
60 │ 7.62569 0.220958 -0.110479 0.110479
61 │ 7.78573 0.203967 -0.101984 0.101984
62 │ 7.94578 0.188283 -0.0941414 0.0941414
63 │ 8.10583 0.173805 -0.0869024 0.0869024
64 │ 8.26588 0.16044 -0.08022 0.08022
65 │ 8.42593 0.148103 -0.0740514 0.0740514
66 │ 8.58598 0.136714 -0.0683572 0.0683572
67 │ 8.74603 0.126202 -0.0631008 0.0631008
68 │ 8.90607 0.116497 -0.0582486 0.0582486
69 │ 9.06612 0.107539 -0.0537696 0.0537696
70 │ 9.22617 0.0992699 -0.0496349 0.0496349
71 │ 9.38622 0.0916365 -0.0458182 0.0458182
72 │ 9.54627 0.08459 -0.042295 0.042295
73 │ 9.70632 0.0780854 -0.0390427 0.0390427
74 │ 9.86636 0.072081 -0.0360405 0.0360405
75 │ 10.0 0.0674227 -0.0337113 0.0337113NB: this displays Variables from the specified model Domain, for this example all Variables are in the "global" Domain.
Output for a single Variable (eg for further analysis) can be retrieved with:
A = PALEOmodel.get_array(paleorun.output, "global.A") # an xarray like object with data and coordinates
A.values # a Vector of values at each model time75-element Vector{Float64}:
10.0
9.977578181715574
9.955206637224666
9.899946959591565
9.844959937226175
9.790266946099884
9.70788605710625
9.564468460248714
9.165173062310783
8.974913766808431
⋮
0.1262016351669181
0.11649728979981157
0.10753916550091484
0.0992698811835523
0.09163646811182256
0.08459003061041624
0.0780854328643422
0.07208100980236412
0.06742268705544686Example 2 Using a PALEO Reservoir
The standard way to configure PALEO models is to use a combination of ReactionReservoirs (from the Reaction catalog provided by the PALEOboxes package) to define model state Variables and provide setup and initialisation, and then link to these Variables from other Reactions.
Additional code files
In this example, the Reaction code now only needs to implement a 'do' method (file examples/reservoirs/reactions_ex2.jl):
module Example2
import PALEOboxes as PB
"""
ReactionExample2
Minimal example, first order decay of a variable.
Use in conjunction with a ReservoirScalar for quantity A.
"""
Base.@kwdef mutable struct ReactionExample2{P} <: PB.AbstractReaction
base::PB.ReactionBase
pars::P = PB.ParametersTuple(
PB.ParDouble("kappa", 1.0, units="yr-1", description="first order decay constant"),
)
end
function PB.register_methods!(rj::ReactionExample2)
vars = [
PB.VarDepScalar("A", "mol", "reservoir for species A"),
PB.VarContribScalar("A_sms", "mol yr-1", "reservoir A source - sink"),
PB.VarPropScalar("decay_flux", "mol yr-1", "decay flux from reservoir A"),
]
PB.add_method_do!(rj, do_example2, (PB.VarList_namedtuple(vars), ) )
return nothing
end
# do method, called each main loop timestep
function do_example2(m::PB.ReactionMethod, pars, (varsdata, ), cellrange::PB.AbstractCellRange, deltat)
# mol yr-1 yr-1 mol
varsdata.decay_flux[] = pars.kappa[] * varsdata.A[]
varsdata.A_sms[] -= varsdata.decay_flux[]
return nothing
end
end # moduleyaml configuration file
The model configuration (file examples/reservoirs/config_ex2.yaml) contains a ReactionReservoirScalar from the generic Reaction catalog provided by the PALEOboxes package:
############################################
# Example 2
# First order decay of a scalar variable
###########################################
example2:
domains:
global:
# scalar domain
reactions:
reservoir_A:
class: ReactionReservoirScalar
variable_links:
R*: A*
variable_attributes:
R:initial_value: 10.0
R:norm_value: 10.0
Adecay:
class: ReactionExample2
parameters:
kappa: 0.5
Run script
The script to run the model (file examples/reservoirs/run_ex2.jl) contains:
import PALEOboxes as PB
import PALEOmodel
using Plots
include(joinpath(@__DIR__, "reactions_ex2.jl"))
#####################################################
# Create model
#######################################################
model = PB.create_model_from_config(joinpath(@__DIR__, "config_ex2.yaml"), "example2")
#########################################################
# Initialize
##########################################################
initial_state, modeldata = PALEOmodel.initialize!(model)
#####################################################################
# Optional: call ODE function to check derivative
#######################################################################
initial_deriv = similar(initial_state)
PALEOmodel.ODE.ModelODE(modeldata)(initial_deriv, initial_state , nothing, 0.0)
println("initial_state: ", initial_state)
println("initial_deriv: ", initial_deriv)
#################################################################
# Integrate vs time
##################################################################
# create a Run object to hold model and output
paleorun = PALEOmodel.Run(model=model, output = PALEOmodel.OutputWriters.OutputMemory())
println("integrate, ODE")
# first run is slow as it includes JIT time
@time PALEOmodel.ODE.integrate(
paleorun, initial_state, modeldata, (0.0, 10.0),
solvekwargs=(
reltol=1e-5,
# saveat=0.1, # save output every 0.1 yr see https://diffeq.sciml.ai/dev/basics/common_solver_opts/
)
);
########################################
# Plot output
########################################
display(plot(paleorun.output, "global.A"))
display(plot(paleorun.output, "global.decay_flux"))
and produces output:
┌ Info:
│ ================================================================================
│ create_model_from_config: configmodel example2
│ config_file: /home/runner/work/PALEOtutorials.jl/PALEOtutorials.jl/examples/reservoirs/config_ex2.yaml
└ ================================================================================
[ Info: Model.parameters:
┌ Info:
│ ================================================================================
│ creating Domains
└ ================================================================================
[ Info: generated Reaction catalog with 74 Reactions
┌ Info:
│ ================================================================================
│ creating domain 'global' ID=1
└ ================================================================================
┌ Info: create_reaction_from_config: global.Adecay classname ReactionExample2
└ set parameters: [config.yaml] kappa =0.5
┌ Info: create_reaction_from_config: global.reservoir_A classname ReactionReservoirScalar
│ set parameters: [Default] field_data =PALEOboxes.ScalarData
│ set parameters: [Default] const =false
└ set parameters: [Default] state_norm =false
┌ Info:
│ ================================================================================
│ set_model_geometry
└ ================================================================================
┌ Info:
│ ================================================================================
│ register_reaction_methods!
└ ================================================================================
┌ Info:
│ ================================================================================
│ link_variables: first pass
└ ================================================================================
┌ Info:
│ ================================================================================
│ link_variables: register_reaction_dynamic_methods and configure variables
└ ================================================================================
┌ Info: _configure_variables: ReactionReservoirScalar global.reservoir_A variable_links:
│ set variable_links: R -> A
│ set variable_links: R_sms -> A_sms
│ set variable_links: R_norm -> A_norm
│ _configure_variables: ReactionReservoirScalar global.reservoir_A variable_attributes:
│ set attribute: R :norm_value = 10.0
└ set attribute: R :initial_value = 10.0
┌ Info:
│ ================================================================================
│ link_variables: second pass:
└ ================================================================================
┌ Info:
│ ================================================================================
│ link_variables! unlinked variables:
└ ================================================================================
┌ Info:
│ ================================================================================
│ create_model_from_config: done
└ ================================================================================
┌ Info:
│ ================================================================================
│ PALEOmodel.initialize! start
└ ================================================================================
┌ Info:
│ ================================================================================
│ allocate_variables! (modeldata arrays_idx=1)
└ ================================================================================
[ Info: Domain global data dimensions PALEOboxes.NamedDimension[] allocating 4 variables (hostdep=nothing)
[ Info: set_default_solver_view:
┌ Info: SolverView:
│ host-dependent Variables:
│ ------------------------------------------------------------------------------------------------------------------------------------------------
│ Domain operatorID VF_StateExplicit VF_Total VF_Constraint VF_StateTotal VF_State VF_Undefined
│ global 0 1 0 0 0 0 0
│ ------------------------------------------------------------------------------------------------------------------------------------------------
│ Total - 1 0 0 0 0 0
│ ------------------------------------------------------------------------------------------------------------------------------------------------
│ n_state_vars 1 (stateexplicit 1 + statetotal 0 + state 0)
│ n_equations 1 (stateexplicit 1 + total 0 + constraint 0)
│ including all host-dependent non-state Variables
└ host-dependent non-state Variables (:vfunction PB.VF_Undefined): Any[]
┌ Info:
│ ================================================================================
│ initialize_reactiondata! (modeldata arrays_indices=1:1)
└ ================================================================================
┌ Info:
│ ================================================================================
│ dispatch_setup :setup
└ ================================================================================
┌ Info:
│ ================================================================================
│ dispatch_setup :norm_value
└ ================================================================================
[ Info: init_values! :norm_value global.A = 10.0
┌ Info:
│ ================================================================================
│ dispatch_setup :initial_value
└ ================================================================================
[ Info: init_values! :initial_value global.A = 10.0
┌ Info:
│ ================================================================================
│ PALEOmodel.initialize! done
└ ================================================================================
initial_state: [10.0]
initial_deriv: [-5.0]
integrate, ODE
┌ Info:
│ ================================================================================
│ PALEOmodel.ODE.integrate:
│ tspan: (0.0, 10.0)
│ algorithm: Sundials.CVODE_BDF{:Newton, :Dense, Nothing, Nothing}(0, 0, 0, false, 10, 5, 7, 3, 10, nothing, nothing, 0)
│ Jacobian: NoJacobian
│ steadystate: false
│ using 1 BLAS threads
└ ================================================================================
[ Info: ODEfunction: using Jacobian NoJacobian
[ Info: jac_config_ode: jac_ad=NoJacobian
0.000167 seconds (920 allocations: 39.797 KiB)
┌ Info:
│ ================================================================================
│ print_sol_stats:
│ retcode=Success
│ alg=Sundials.CVODE_BDF{:Newton, :Dense, Nothing, Nothing}(0, 0, 0, false, 10, 5, 7, 3, 10, nothing, nothing, 0)
│ stats=SciMLBase.DEStats
│ Number of function 1 evaluations: 92
│ Number of function 2 evaluations: 0
│ Number of W matrix evaluations: 11
│ Number of linear solves: 0
│ Number of Jacobians created: 2
│ Number of nonlinear solver iterations: 89
│ Number of nonlinear solver convergence failures: 0
│ Number of fixed-point solver iterations: 0
│ Number of fixed-point solver convergence failures: 0
│ Number of rootfind condition calls: 0
│ Number of accepted steps: 72
│ Number of rejected steps: 2
│ length(sol.t) 75
│ size(sol) (1, 75)
└ ================================================================================
[ Info: ODE.calc_output_sol!: 75 records
[ Info: ODE.calc_output_sol!: done
┌ Info:
│ ================================================================================
│ PALEOmodel.ODE.integrate: done
└ ================================================================================
0.001404 seconds (4.28 k allocations: 1.015 MiB)
Displaying model structure and variables
All metadata for model variables can be displayed with PB.show_variables:
show(PB.show_variables(model), allcols=true, allrows=true) # display in REPL
# vscodedisplay(PB.show_variables(model)) # more convenient when using VS code4×9 DataFrame
Row │ domain name type units vfunction space data_dims field_data description
│ String String DataType String Variable… DataType Tuple{} DataType String
─────┼─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
1 │ global A VariableDomPropDep mol VF_StateExplicit ScalarSpace () ScalarData scalar reservoir
2 │ global A_norm VariableDomPropDep VF_Undefined ScalarSpace () ScalarData scalar reservoir normalized
3 │ global A_sms VariableDomContribTarget mol yr-1 VF_Deriv ScalarSpace () ScalarData scalar reservoir source-sinks
4 │ global decay_flux VariableDomPropDep mol yr-1 VF_Undefined ScalarSpace () ScalarData decay flux from reservoir Athere is one more Variable "global.A_norm", the normalized value of "global.A", provided by the generic ReactionReservoirScalar and not used here.
The linking of the "global.A" variable illustrates the key difference between this example and Example 1:
PB.show_links(model, "global.A") PALEOboxes.VariableDomPropDep "global.A" links:
property: nothing
dependencies: ["global.Adecay.do_example2.A", "global.reservoir_A.setup_reactionreservoirscalar.R", "global.reservoir_A.do_reactionreservoirscalar.R"]The "global.A" state Variable now has three dependencies:
- the
setup_initialvalue_vars_defaultmethod of the reaction namedreservoir_Ain the yaml file (a ReactionReservoirScalar), where the variable has local nameRand has been renamed toAin the yaml file. - the
do_reactionreservoirscalarmethod of the reaction namedreservoir_A(which calculates normalized valueA_norm, not needed here). - the
do_example1method of the Reaction namedAdecayin the yaml config file, where the variable has local nameA, which reads the value.
Example 3 Transfer between two Reservoirs
Generalizing the Reaction to transfer between two reservoirs.
Additional code files
The Reaction code (file examples/reservoirs/reactions_ex3.jl) now produces an output_flux, and has been generalized to operate on a generic input_particle Reservoir:
module Example3
import PALEOboxes as PB
"""
ReactionExample3
Minimal example, first order decay of a variable and transfer of decay flux.
Use config file `variable_links:` to rename `input_particle*` and `output_flux`.
"""
Base.@kwdef mutable struct ReactionExample3{P} <: PB.AbstractReaction
base::PB.ReactionBase
pars::P = PB.ParametersTuple(
PB.ParDouble("kappa", 1.0, units="yr-1", description="first order decay constant"),
)
end
function PB.register_methods!(rj::ReactionExample3)
vars = [
PB.VarDepScalar("input_particle", "mol", "reservoir for input"),
PB.VarContribScalar("input_particle_sms", "mol yr-1", "reservoir input source - sink"),
PB.VarPropScalar("decay_flux", "mol yr-1", "decay flux"),
PB.VarContribScalar("output_flux", "mol yr-1", "output flux"),
]
PB.add_method_do!(rj, do_example3, (PB.VarList_namedtuple(vars), ) )
return nothing
end
# do method, called each main loop timestep
function do_example3(m::PB.ReactionMethod, pars, (varsdata, ), cellrange::PB.AbstractCellRange, deltat)
# mol yr-1 yr-1 mol
varsdata.decay_flux[] = pars.kappa[] * varsdata.input_particle[]
varsdata.input_particle_sms[] -= varsdata.decay_flux[]
varsdata.output_flux[] += varsdata.decay_flux[]
return nothing
end
end # moduleyaml configuration file
The model configuration (file examples/reservoirs/config_ex3.yaml) contains two Reservoirs A and B, and additional configuration for ReactionExample3 to rename the generic input_particle to link to Reservoir A and output_flux to link to Reservoir B:
############################################
# Example 3
# First order decay of a scalar variable
# with transfer of flux to a second variable
###########################################
example3:
domains:
global:
# scalar domain
reactions:
reservoir_A:
class: ReactionReservoirScalar
variable_links:
R*: A*
variable_attributes:
R:initial_value: 10.0
R:norm_value: 10.0
reservoir_B:
class: ReactionReservoirScalar
variable_links:
R*: B*
variable_attributes:
R:initial_value: 0.0
R:norm_value: 10.0
Adecay:
class: ReactionExample3
parameters:
kappa: 0.5
variable_links:
input_particle*: A*
output_flux: B_sms
Run script
The script to run the model (file examples/reservoirs/run_ex3.jl) contains:
import PALEOboxes as PB
import PALEOmodel
using Plots
include(joinpath(@__DIR__, "reactions_ex3.jl"))
#####################################################
# Create model
#######################################################
model = PB.create_model_from_config(joinpath(@__DIR__, "config_ex3.yaml"), "example3")
#########################################################
# Initialize
##########################################################
initial_state, modeldata = PALEOmodel.initialize!(model)
#####################################################################
# Optional: call ODE function to check derivative
#######################################################################
initial_deriv = similar(initial_state)
PALEOmodel.ODE.ModelODE(modeldata)(initial_deriv, initial_state , nothing, 0.0)
println("initial_state: ", initial_state)
println("initial_deriv: ", initial_deriv)
#################################################################
# Integrate vs time
##################################################################
# create a Run object to hold model and output
paleorun = PALEOmodel.Run(model=model, output = PALEOmodel.OutputWriters.OutputMemory())
println("integrate, ODE")
# first run is slow as it includes JIT time
@time PALEOmodel.ODE.integrate(
paleorun, initial_state, modeldata, (0.0, 10.0),
solvekwargs=(
reltol=1e-5,
# saveat=0.1, # save output every 0.1 yr see https://diffeq.sciml.ai/dev/basics/common_solver_opts/
)
);
########################################
# Plot output
########################################
display(plot(paleorun.output, ["global.A", "global.B"]; ylabel="reservoir (mol)"))
display(plot(paleorun.output, "global.decay_flux"))and produces output showing the transfer between two Reservoirs:
┌ Info:
│ ================================================================================
│ create_model_from_config: configmodel example3
│ config_file: /home/runner/work/PALEOtutorials.jl/PALEOtutorials.jl/examples/reservoirs/config_ex3.yaml
└ ================================================================================
[ Info: Model.parameters:
┌ Info:
│ ================================================================================
│ creating Domains
└ ================================================================================
[ Info: generated Reaction catalog with 74 Reactions
┌ Info:
│ ================================================================================
│ creating domain 'global' ID=1
└ ================================================================================
┌ Info: create_reaction_from_config: global.Adecay classname ReactionExample3
└ set parameters: [config.yaml] kappa =0.5
┌ Info: create_reaction_from_config: global.reservoir_B classname ReactionReservoirScalar
│ set parameters: [Default] field_data =PALEOboxes.ScalarData
│ set parameters: [Default] const =false
└ set parameters: [Default] state_norm =false
┌ Info: create_reaction_from_config: global.reservoir_A classname ReactionReservoirScalar
│ set parameters: [Default] field_data =PALEOboxes.ScalarData
│ set parameters: [Default] const =false
└ set parameters: [Default] state_norm =false
┌ Info:
│ ================================================================================
│ set_model_geometry
└ ================================================================================
┌ Info:
│ ================================================================================
│ register_reaction_methods!
└ ================================================================================
┌ Info:
│ ================================================================================
│ link_variables: first pass
└ ================================================================================
┌ Info:
│ ================================================================================
│ link_variables: register_reaction_dynamic_methods and configure variables
└ ================================================================================
┌ Info: _configure_variables: ReactionExample3 global.Adecay variable_links:
│ set variable_links: input_particle -> A
│ set variable_links: input_particle_sms -> A_sms
└ set variable_links: output_flux -> B_sms
┌ Info: _configure_variables: ReactionReservoirScalar global.reservoir_B variable_links:
│ set variable_links: R -> B
│ set variable_links: R_sms -> B_sms
│ set variable_links: R_norm -> B_norm
│ _configure_variables: ReactionReservoirScalar global.reservoir_B variable_attributes:
│ set attribute: R :norm_value = 10.0
└ set attribute: R :initial_value = 0.0
┌ Info: _configure_variables: ReactionReservoirScalar global.reservoir_A variable_links:
│ set variable_links: R -> A
│ set variable_links: R_sms -> A_sms
│ set variable_links: R_norm -> A_norm
│ _configure_variables: ReactionReservoirScalar global.reservoir_A variable_attributes:
│ set attribute: R :norm_value = 10.0
└ set attribute: R :initial_value = 10.0
┌ Info:
│ ================================================================================
│ link_variables: second pass:
└ ================================================================================
┌ Info:
│ ================================================================================
│ link_variables! unlinked variables:
└ ================================================================================
┌ Info:
│ ================================================================================
│ create_model_from_config: done
└ ================================================================================
┌ Info:
│ ================================================================================
│ PALEOmodel.initialize! start
└ ================================================================================
┌ Info:
│ ================================================================================
│ allocate_variables! (modeldata arrays_idx=1)
└ ================================================================================
[ Info: Domain global data dimensions PALEOboxes.NamedDimension[] allocating 7 variables (hostdep=nothing)
[ Info: set_default_solver_view:
┌ Info: SolverView:
│ host-dependent Variables:
│ ------------------------------------------------------------------------------------------------------------------------------------------------
│ Domain operatorID VF_StateExplicit VF_Total VF_Constraint VF_StateTotal VF_State VF_Undefined
│ global 0 2 0 0 0 0 0
│ ------------------------------------------------------------------------------------------------------------------------------------------------
│ Total - 2 0 0 0 0 0
│ ------------------------------------------------------------------------------------------------------------------------------------------------
│ n_state_vars 2 (stateexplicit 2 + statetotal 0 + state 0)
│ n_equations 2 (stateexplicit 2 + total 0 + constraint 0)
│ including all host-dependent non-state Variables
└ host-dependent non-state Variables (:vfunction PB.VF_Undefined): Any[]
┌ Info:
│ ================================================================================
│ initialize_reactiondata! (modeldata arrays_indices=1:1)
└ ================================================================================
┌ Info:
│ ================================================================================
│ dispatch_setup :setup
└ ================================================================================
┌ Info:
│ ================================================================================
│ dispatch_setup :norm_value
└ ================================================================================
[ Info: init_values! :norm_value global.B = 10.0
[ Info: init_values! :norm_value global.A = 10.0
┌ Info:
│ ================================================================================
│ dispatch_setup :initial_value
└ ================================================================================
[ Info: init_values! :initial_value global.B = 0.0
[ Info: init_values! :initial_value global.A = 10.0
┌ Info:
│ ================================================================================
│ PALEOmodel.initialize! done
└ ================================================================================
initial_state: [0.0, 10.0]
initial_deriv: [5.0, -5.0]
integrate, ODE
┌ Info:
│ ================================================================================
│ PALEOmodel.ODE.integrate:
│ tspan: (0.0, 10.0)
│ algorithm: Sundials.CVODE_BDF{:Newton, :Dense, Nothing, Nothing}(0, 0, 0, false, 10, 5, 7, 3, 10, nothing, nothing, 0)
│ Jacobian: NoJacobian
│ steadystate: false
│ using 1 BLAS threads
└ ================================================================================
[ Info: ODEfunction: using Jacobian NoJacobian
[ Info: jac_config_ode: jac_ad=NoJacobian
0.566189 seconds (1.13 M allocations: 78.491 MiB, 5.45% gc time, 99.96% compilation time)
┌ Info:
│ ================================================================================
│ print_sol_stats:
│ retcode=Success
│ alg=Sundials.CVODE_BDF{:Newton, :Dense, Nothing, Nothing}(0, 0, 0, false, 10, 5, 7, 3, 10, nothing, nothing, 0)
│ stats=SciMLBase.DEStats
│ Number of function 1 evaluations: 64
│ Number of function 2 evaluations: 0
│ Number of W matrix evaluations: 16
│ Number of linear solves: 0
│ Number of Jacobians created: 2
│ Number of nonlinear solver iterations: 61
│ Number of nonlinear solver convergence failures: 0
│ Number of fixed-point solver iterations: 0
│ Number of fixed-point solver convergence failures: 0
│ Number of rootfind condition calls: 0
│ Number of accepted steps: 53
│ Number of rejected steps: 1
│ length(sol.t) 55
│ size(sol) (2, 55)
└ ================================================================================
[ Info: ODE.calc_output_sol!: 55 records
[ Info: ODE.calc_output_sol!: done
┌ Info:
│ ================================================================================
│ PALEOmodel.ODE.integrate: done
└ ================================================================================
0.704074 seconds (1.27 M allocations: 89.452 MiB, 4.38% gc time, 99.53% compilation time)
Displaying model structure and variables
All metadata for model variables can be displayed with PB.show_variables:
show(PB.show_variables(model), allcols=true, allrows=true) # display in REPL
# vscodedisplay(PB.show_variables(model)) # more convenient when using VS code7×9 DataFrame
Row │ domain name type units vfunction space data_dims field_data description
│ String String DataType String Variable… DataType Tuple{} DataType String
─────┼─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
1 │ global A VariableDomPropDep mol VF_StateExplicit ScalarSpace () ScalarData scalar reservoir
2 │ global A_norm VariableDomPropDep VF_Undefined ScalarSpace () ScalarData scalar reservoir normalized
3 │ global A_sms VariableDomContribTarget mol yr-1 VF_Deriv ScalarSpace () ScalarData scalar reservoir source-sinks
4 │ global B VariableDomPropDep mol VF_StateExplicit ScalarSpace () ScalarData scalar reservoir
5 │ global B_norm VariableDomPropDep VF_Undefined ScalarSpace () ScalarData scalar reservoir normalized
6 │ global B_sms VariableDomContribTarget mol yr-1 VF_Deriv ScalarSpace () ScalarData scalar reservoir source-sinks
7 │ global decay_flux VariableDomPropDep mol yr-1 VF_Undefined ScalarSpace () ScalarData decay fluxWe now have additional Variables corresponding to the B reservoir.
Example 4 Transfer between Domains
The ReactionExample3 from the previous example can be reused in a different model configuration that includes flux transfer between Reservoirs in two Domains.
Additional code files
This example reuses the PALEO reactions from Example 3 Transfer between two Reservoirs
yaml configuration file
The model configuration (file examples/reservoirs/config_ex4.yaml) contains two Domains Box1 and Box2, and a flux coupler Domain fluxBoxes. The ReactionSum in the global Domain tracks conservation:
############################################
# Example 4
# First order decay of a scalar variable
# with transfer of flux to a second variable
# Transfer between two Domains
###########################################
example4:
domains:
fluxBoxes:
reactions:
flux:
class: ReactionFluxTarget
parameters:
target_prefix: flux_
fluxlist: ["B"]
global:
reactions:
sum_E:
class: ReactionSum
parameters:
vars_to_add: ["Box1.A", "Box2.B"]
variable_links:
sum: E_total
Box1:
# scalar domain
reactions:
reservoir_A:
class: ReactionReservoirScalar
variable_links:
R*: A*
variable_attributes:
R:initial_value: 10.0
R:norm_value: 10.0
Adecay:
class: ReactionExample3
parameters:
kappa: 0.5
variable_links:
input_particle*: A*
output_flux: fluxBoxes.flux_B
Box2:
# scalar domain
reactions:
reservoir_B:
class: ReactionReservoirScalar
variable_links:
R*: B*
variable_attributes:
R:initial_value: 0.0
R:norm_value: 10.0
transfer_fluxBoxes:
class: ReactionFluxTransfer
parameters:
input_fluxes: fluxBoxes.flux_$fluxname$
output_fluxes: $fluxname$_smsRun script
The script to run the model (file examples/reservoirs/run_ex4.jl) contains:
import PALEOboxes as PB
import PALEOmodel
using Plots
include(joinpath(@__DIR__, "reactions_ex3.jl"))
#####################################################
# Create model
#######################################################
model = PB.create_model_from_config(joinpath(@__DIR__, "config_ex4.yaml"), "example4")
#########################################################
# Initialize
##########################################################
initial_state, modeldata = PALEOmodel.initialize!(model)
#####################################################################
# Optional: call ODE function to check derivative
#######################################################################
initial_deriv = similar(initial_state)
PALEOmodel.ODE.ModelODE(modeldata)(initial_deriv, initial_state , nothing, 0.0)
println("initial_state: ", initial_state)
println("initial_deriv: ", initial_deriv)
#################################################################
# Integrate vs time
##################################################################
# create a Run object to hold model and output
paleorun = PALEOmodel.Run(model=model, output = PALEOmodel.OutputWriters.OutputMemory())
println("integrate, ODE")
# first run is slow as it includes JIT time
@time PALEOmodel.ODE.integrate(
paleorun, initial_state, modeldata, (0.0, 10.0),
solvekwargs=(
reltol=1e-5,
# saveat=0.1, # save output every 0.1 yr see https://diffeq.sciml.ai/dev/basics/common_solver_opts/
)
);
############################
# Table of output
###########################
# vscodedisplay(
# PB.get_table(paleorun.output, ["Box1.A", "Box2.B", "global.E_total", "Box1.decay_flux", "fluxBoxes.flux_B"]),
# "Example 4"
# )
########################################
# Plot output
########################################
display(plot(paleorun.output, ["Box1.A", "Box2.B", "global.E_total"]; ylabel="reservoir (mol)"))
display(plot(paleorun.output, ["Box1.decay_flux", "fluxBoxes.flux_B"]; ylabel="flux (mol yr-1)"))
and produces output showing the transfer between two Reservoirs in different Domains via the fluxBoxes flux coupler:
┌ Info:
│ ================================================================================
│ create_model_from_config: configmodel example4
│ config_file: /home/runner/work/PALEOtutorials.jl/PALEOtutorials.jl/examples/reservoirs/config_ex4.yaml
└ ================================================================================
[ Info: Model.parameters:
┌ Info:
│ ================================================================================
│ creating Domains
└ ================================================================================
[ Info: generated Reaction catalog with 74 Reactions
┌ Info:
│ ================================================================================
│ creating domain 'global' ID=1
└ ================================================================================
┌ Info: create_reaction_from_config: global.sum_E classname ReactionSum
│ set parameters: [config.yaml] vars_to_add =["Box1.A", "Box2.B"]
│ set parameters: [Default] vars_prefix =
│ set parameters: [Default] component_to_add =0
└ set parameters: [Default] vectorsum =false
┌ Info:
│ ================================================================================
│ creating domain 'fluxBoxes' ID=2
└ ================================================================================
┌ Info: create_reaction_from_config: fluxBoxes.flux classname ReactionFluxTarget
│ set parameters: [config.yaml] fluxlist =["B"]
│ set parameters: [config.yaml] target_prefix =flux_
│ set parameters: [Default] flux_totals =false
└ set parameters: [Default] const_stub =false
┌ Info:
│ ================================================================================
│ creating domain 'Box2' ID=3
└ ================================================================================
┌ Info: create_reaction_from_config: Box2.reservoir_B classname ReactionReservoirScalar
│ set parameters: [Default] field_data =PALEOboxes.ScalarData
│ set parameters: [Default] const =false
└ set parameters: [Default] state_norm =false
┌ Info: create_reaction_from_config: Box2.transfer_fluxBoxes classname ReactionFluxTransfer
│ set parameters: [config.yaml] input_fluxes =fluxBoxes.flux_$fluxname$
│ set parameters: [config.yaml] output_fluxes =$fluxname$_sms
│ set parameters: [Default] transfer_multiplier =1.0
└ set parameters: [Default] transfer_matrix =Identity
┌ Info:
│ ================================================================================
│ creating domain 'Box1' ID=4
└ ================================================================================
┌ Info: create_reaction_from_config: Box1.Adecay classname ReactionExample3
└ set parameters: [config.yaml] kappa =0.5
┌ Info: create_reaction_from_config: Box1.reservoir_A classname ReactionReservoirScalar
│ set parameters: [Default] field_data =PALEOboxes.ScalarData
│ set parameters: [Default] const =false
└ set parameters: [Default] state_norm =false
┌ Info:
│ ================================================================================
│ set_model_geometry
└ ================================================================================
┌ Info:
│ ================================================================================
│ register_reaction_methods!
└ ================================================================================
┌ Info: register_methods: global.sum_E PALEOboxes.VariableStats.ReactionSum
│ add 1 * Box1.A
└ add 1 * Box2.B
┌ Info:
│ ================================================================================
│ link_variables: first pass
└ ================================================================================
┌ Info:
│ ================================================================================
│ link_variables: register_reaction_dynamic_methods and configure variables
└ ================================================================================
[ Info: Reaction global.sum_E Variable global.E_total adding data=PALEOboxes.ScalarData
┌ Info: _configure_variables: ReactionSum global.sum_E variable_links:
└ set variable_links: sum -> E_total
┌ Info: _configure_variables: ReactionReservoirScalar Box2.reservoir_B variable_links:
│ set variable_links: R -> B
│ set variable_links: R_sms -> B_sms
│ set variable_links: R_norm -> B_norm
│ _configure_variables: ReactionReservoirScalar Box2.reservoir_B variable_attributes:
│ set attribute: R :norm_value = 10.0
└ set attribute: R :initial_value = 0.0
┌ Info: _configure_variables: ReactionExample3 Box1.Adecay variable_links:
│ set variable_links: input_particle -> A
│ set variable_links: input_particle_sms -> A_sms
└ set variable_links: output_flux -> fluxBoxes.flux_B
┌ Info: _configure_variables: ReactionReservoirScalar Box1.reservoir_A variable_links:
│ set variable_links: R -> A
│ set variable_links: R_sms -> A_sms
│ set variable_links: R_norm -> A_norm
│ _configure_variables: ReactionReservoirScalar Box1.reservoir_A variable_attributes:
│ set attribute: R :norm_value = 10.0
└ set attribute: R :initial_value = 10.0
┌ Info:
│ ================================================================================
│ link_variables: second pass:
└ ================================================================================
┌ Info:
│ ================================================================================
│ link_variables! unlinked variables:
└ ================================================================================
┌ Info:
│ ================================================================================
│ create_model_from_config: done
└ ================================================================================
┌ Info:
│ ================================================================================
│ PALEOmodel.initialize! start
└ ================================================================================
┌ Info:
│ ================================================================================
│ allocate_variables! (modeldata arrays_idx=1)
└ ================================================================================
[ Info: Domain global data dimensions PALEOboxes.NamedDimension[] allocating 1 variables (hostdep=nothing)
[ Info: Domain fluxBoxes data dimensions PALEOboxes.NamedDimension[] allocating 1 variables (hostdep=nothing)
[ Info: Domain Box2 data dimensions PALEOboxes.NamedDimension[] allocating 3 variables (hostdep=nothing)
[ Info: Domain Box1 data dimensions PALEOboxes.NamedDimension[] allocating 4 variables (hostdep=nothing)
[ Info: set_default_solver_view:
┌ Info: SolverView:
│ host-dependent Variables:
│ ------------------------------------------------------------------------------------------------------------------------------------------------
│ Domain operatorID VF_StateExplicit VF_Total VF_Constraint VF_StateTotal VF_State VF_Undefined
│ global 0 0 0 0 0 0 0
│ fluxBoxes 0 0 0 0 0 0 0
│ Box2 0 1 0 0 0 0 0
│ Box1 0 1 0 0 0 0 0
│ ------------------------------------------------------------------------------------------------------------------------------------------------
│ Total - 2 0 0 0 0 0
│ ------------------------------------------------------------------------------------------------------------------------------------------------
│ n_state_vars 2 (stateexplicit 2 + statetotal 0 + state 0)
│ n_equations 2 (stateexplicit 2 + total 0 + constraint 0)
│ including all host-dependent non-state Variables
└ host-dependent non-state Variables (:vfunction PB.VF_Undefined): Any[]
┌ Info:
│ ================================================================================
│ initialize_reactiondata! (modeldata arrays_indices=1:1)
└ ================================================================================
┌ Info: prepare_do_transfer: ReactionMethod Box2.transfer_fluxBoxes.do_transfer
│ active fluxes:
│ B fluxBoxes.flux_B -> Box2.B_sms
└ using Identity transfer matrix * 1.0
┌ Info:
│ ================================================================================
│ dispatch_setup :setup
└ ================================================================================
┌ Info:
│ ================================================================================
│ dispatch_setup :norm_value
└ ================================================================================
[ Info: init_values! :norm_value Box2.B = 10.0
[ Info: init_values! :norm_value Box1.A = 10.0
┌ Info:
│ ================================================================================
│ dispatch_setup :initial_value
└ ================================================================================
[ Info: init_values! :initial_value Box2.B = 0.0
[ Info: init_values! :initial_value Box1.A = 10.0
┌ Info:
│ ================================================================================
│ PALEOmodel.initialize! done
└ ================================================================================
initial_state: [0.0, 10.0]
initial_deriv: [5.0, -5.0]
integrate, ODE
┌ Info:
│ ================================================================================
│ PALEOmodel.ODE.integrate:
│ tspan: (0.0, 10.0)
│ algorithm: Sundials.CVODE_BDF{:Newton, :Dense, Nothing, Nothing}(0, 0, 0, false, 10, 5, 7, 3, 10, nothing, nothing, 0)
│ Jacobian: NoJacobian
│ steadystate: false
│ using 1 BLAS threads
└ ================================================================================
[ Info: ODEfunction: using Jacobian NoJacobian
[ Info: jac_config_ode: jac_ad=NoJacobian
0.000153 seconds (684 allocations: 33.984 KiB)
┌ Info:
│ ================================================================================
│ print_sol_stats:
│ retcode=Success
│ alg=Sundials.CVODE_BDF{:Newton, :Dense, Nothing, Nothing}(0, 0, 0, false, 10, 5, 7, 3, 10, nothing, nothing, 0)
│ stats=SciMLBase.DEStats
│ Number of function 1 evaluations: 64
│ Number of function 2 evaluations: 0
│ Number of W matrix evaluations: 16
│ Number of linear solves: 0
│ Number of Jacobians created: 2
│ Number of nonlinear solver iterations: 61
│ Number of nonlinear solver convergence failures: 0
│ Number of fixed-point solver iterations: 0
│ Number of fixed-point solver convergence failures: 0
│ Number of rootfind condition calls: 0
│ Number of accepted steps: 53
│ Number of rejected steps: 1
│ length(sol.t) 55
│ size(sol) (2, 55)
└ ================================================================================
[ Info: ODE.calc_output_sol!: 55 records
[ Info: ODE.calc_output_sol!: done
┌ Info:
│ ================================================================================
│ PALEOmodel.ODE.integrate: done
└ ================================================================================
0.036726 seconds (30.06 k allocations: 2.750 MiB, 94.57% compilation time)
Displaying model structure and variables
All metadata for model variables can be displayed with PB.show_variables:
show(PB.show_variables(model), allcols=true, allrows=true) # display in REPL
# vscodedisplay(PB.show_variables(model)) # more convenient when using VS code9×9 DataFrame
Row │ domain name type units vfunction space data_dims field_data description
│ String String DataType String Variable… DataType Tuple{} DataType String
─────┼────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
1 │ Box1 A VariableDomPropDep mol VF_StateExplicit ScalarSpace () ScalarData scalar reservoir
2 │ Box1 A_norm VariableDomPropDep VF_Undefined ScalarSpace () ScalarData scalar reservoir normalized
3 │ Box1 A_sms VariableDomContribTarget mol yr-1 VF_Deriv ScalarSpace () ScalarData scalar reservoir source-sinks
4 │ Box1 decay_flux VariableDomPropDep mol yr-1 VF_Undefined ScalarSpace () ScalarData decay flux
5 │ Box2 B VariableDomPropDep mol VF_StateExplicit ScalarSpace () ScalarData scalar reservoir
6 │ Box2 B_norm VariableDomPropDep VF_Undefined ScalarSpace () ScalarData scalar reservoir normalized
7 │ Box2 B_sms VariableDomContribTarget mol yr-1 VF_Deriv ScalarSpace () ScalarData scalar reservoir source-sinks
8 │ fluxBoxes flux_B VariableDomContribTarget mol yr-1 VF_Undefined CellSpace () ScalarData flux input
9 │ global E_total VariableDomPropDep unknown VF_Undefined ScalarSpace () ScalarData sum of specified variablesThis shows that we now have four Domains:
Box1containingreservoir_A,Adecay, and associated VariablesBox2containingreservoir_Band associated VariablesfluxBoxescontaining a Variableflux_Bcreated by theReactionFluxTargetglobalcontainingE_totalto check conservation.
An aside on ordering of ReactionMethods
It is rarely necessary or useful to look at this - shown here just to illustrate how PALEO works
PALEO orders ReactionMethods based on the dependency information defined by the linked Variables (a method defining a Property must run before all methods that link to it as Dependencies, a method defining a Target must run after all methods that link to it as Contributors).
The order in which ReactionMethods are called during each timestep is stored in model.sorted_methods_do (a struct PB.MethodSort) and can be displayed using:
model.sorted_methods_doMethodSort
group=BitSet([1, 2, 4, 5])
global.sum_E.do_scalarsum
Box2.reservoir_B.do_reactionreservoirscalar
Box1.Adecay.do_example3
Box1.reservoir_A.do_reactionreservoirscalar
group=BitSet([3])
Box2.transfer_fluxBoxes.do_transfer
Each group of methods has no internal dependencies, but depends on methods in previous groups. Here there is only one dependency, method Box2.transfer_fluxBoxes.do_transfer must run after the decay flux is calculated by Box1.reservoir_A.do_reactionreservoirscalar.
Example 5 Isotopes and Rayleigh fractionation
PALEO Variables can represent isotopes by setting the :field_data Attribute. Currently IsotopeLinear is supported, which represents a linearized approximation to a single isotope using two components total and total x delta. Standard arithmetic (additon/subtraction, multiplication by a scalar) is supported.
Additional code files
The Reaction code (file examples/reservoirs/reactions_ex5.jl) now requires additional Parameters field_data to set the data type, and Delta to set the fractionation assumed to occur during the decay/transfer.
module Example5
import PALEOboxes as PB
"""
ReactionExample5
Minimal example, first order decay of a variable and transfer of decay flux with isotope fractionation
Use config file `variable_links:` to rename `input_particle*` and `output_flux`.
"""
Base.@kwdef mutable struct ReactionExample5{P} <: PB.AbstractReaction
base::PB.ReactionBase
pars::P = PB.ParametersTuple(
PB.ParDouble("kappa", 1.0, units="yr-1",
description="first order decay constant"),
PB.ParType(PB.AbstractData, "field_data", PB.ScalarData,
allowed_values=PB.IsotopeTypes,
description="disable / enable isotopes and specify isotope type"),
PB.ParDouble("Delta", -10.0, units="per mil",
description="isotope fractionation: delta(output_flux) = delta(input_particle) + Delta"),
)
end
function PB.register_methods!(rj::ReactionExample5)
IsotopeType = rj.pars.field_data[]
vars = [
PB.VarDepScalar("input_particle", "mol", "reservoir for input",
attributes=(:field_data=>IsotopeType,)),
PB.VarContribScalar("input_particle_sms", "mol yr-1", "reservoir input source - sink",
attributes=(:field_data=>IsotopeType,)),
PB.VarDepScalar("input_particle_delta", "per mil", "reservoir for input isotope delta"),
PB.VarPropScalar("decay_flux", "mol yr-1", "decay flux",
attributes=(:field_data=>IsotopeType,)),
PB.VarContribScalar("output_flux", "mol yr-1", "output flux",
attributes=(:field_data=>IsotopeType,)),
]
PB.add_method_do!(rj, do_example5, (PB.VarList_namedtuple(vars), ), p=(IsotopeType, ) )
return nothing
end
# do method, called each main loop timestep
function do_example5(m::PB.ReactionMethod, pars, (varsdata, ), cellrange::PB.AbstractCellRange, deltat)
(IsotopeType, ) = m.p
# mol yr-1 yr-1 mol
varsdata.decay_flux[] = pars.kappa[] * @PB.isotope_totaldelta(IsotopeType,
PB.get_total(varsdata.input_particle[]), # total
varsdata.input_particle_delta[] + pars.Delta[]) # delta
varsdata.input_particle_sms[] -= varsdata.decay_flux[]
varsdata.output_flux[] += varsdata.decay_flux[]
return nothing
end
end # moduleyaml configuration file
The model configuration (file examples/reservoirs/config_ex5.yaml) contains a global parameter EIsotope used to set the isotope Type for all Reactions affecting the hypothetical element E.
############################################
# Example 5
# First order decay of a scalar variable
# with transfer of flux to a second variable
# Two domains, with isotopes
###########################################
example5:
parameters:
EIsotope: IsotopeLinear # hypothetical element E
domains:
fluxBoxes:
reactions:
flux:
class: ReactionFluxTarget
parameters:
target_prefix: flux_
fluxlist: ["B::EIsotope"]
global:
reactions:
sum_E:
class: ReactionSum
parameters:
vars_to_add: ["Box1.A", "Box2.B"]
variable_links:
sum: E_total
Box1:
# scalar domain
reactions:
reservoir_A:
class: ReactionReservoirScalar
parameters:
field_data: external%EIsotope
variable_links:
R*: A*
variable_attributes:
R:initial_value: 10.0
R:norm_value: 10.0
Adecay:
class: ReactionExample5
parameters:
field_data: external%EIsotope
kappa: 0.5
variable_links:
input_particle*: A*
output_flux: fluxBoxes.flux_B
Box2:
# scalar domain
reactions:
reservoir_B:
class: ReactionReservoirScalar
parameters:
field_data: external%EIsotope
variable_links:
R*: B*
variable_attributes:
R:initial_value: 0.0
R:norm_value: 10.0
transfer_fluxBoxes:
class: ReactionFluxTransfer
parameters:
input_fluxes: fluxBoxes.flux_$fluxname$
output_fluxes: $fluxname$_smsRun script
The script to run the model (file examples/reservoirs/run_ex5.jl) contains:
import PALEOboxes as PB
import PALEOmodel
using Plots
import Sundials
include(joinpath(@__DIR__, "reactions_ex5.jl"))
#####################################################
# Create model
#######################################################
model = PB.create_model_from_config(joinpath(@__DIR__, "config_ex5.yaml"), "example5")
#########################################################
# Initialize
##########################################################
initial_state, modeldata = PALEOmodel.initialize!(model)
#####################################################################
# Optional: call ODE function to check derivative
#######################################################################
initial_deriv = similar(initial_state)
PALEOmodel.ODE.ModelODE(modeldata)(initial_deriv, initial_state , nothing, 0.0)
println("initial_state: ", initial_state)
println("initial_deriv: ", initial_deriv)
#################################################################
# Integrate vs time
##################################################################
# create a Run object to hold model and output
paleorun = PALEOmodel.Run(model=model, output = PALEOmodel.OutputWriters.OutputMemory())
println("integrate, ODE")
# first run is slow as it includes JIT time
@time PALEOmodel.ODE.integrate(
paleorun, initial_state, modeldata, (0.0, 10.0);
alg=Sundials.CVODE_BDF(), # this is the PALEOmodel default
solvekwargs=(
reltol=1e-5,
# saveat=0.1, # save output every 0.1 yr see https://diffeq.sciml.ai/dev/basics/common_solver_opts/
)
);
############################
# Table of output
###########################
# vscodedisplay(
# PB.get_table(paleorun.output, ["Box1.A", "Box2.B", "global.E_total", "Box1.decay_flux", "fluxBoxes.flux_B"]),
# "Example 5"
# )
########################################
# Plot output
########################################
display(plot(paleorun.output, ["Box1.A", "Box2.B", "global.E_total"]; ylabel="reservoir (mol)"))
display(plot(paleorun.output, ["Box1.decay_flux", "fluxBoxes.flux_B"]; ylabel="flux (mol yr-1)"))
display(plot(paleorun.output, ["Box1.A_delta", "Box1.decay_flux.v_delta", "Box2.B_delta", "global.E_total.v_delta", ]; ylabel="delta (per mil)"))and produces output showing Rayleigh fractionation:
┌ Info:
│ ================================================================================
│ create_model_from_config: configmodel example5
│ config_file: /home/runner/work/PALEOtutorials.jl/PALEOtutorials.jl/examples/reservoirs/config_ex5.yaml
└ ================================================================================
┌ Info: Model.parameters:
└ EIsotope = IsotopeLinear
┌ Info:
│ ================================================================================
│ creating Domains
└ ================================================================================
[ Info: generated Reaction catalog with 74 Reactions
┌ Info:
│ ================================================================================
│ creating domain 'global' ID=1
└ ================================================================================
┌ Info: create_reaction_from_config: global.sum_E classname ReactionSum
│ set parameters: [config.yaml] vars_to_add =["Box1.A", "Box2.B"]
│ set parameters: [Default] vars_prefix =
│ set parameters: [Default] component_to_add =0
└ set parameters: [Default] vectorsum =false
┌ Info:
│ ================================================================================
│ creating domain 'fluxBoxes' ID=2
└ ================================================================================
┌ Info: create_reaction_from_config: fluxBoxes.flux classname ReactionFluxTarget
│ set parameters: [config.yaml] fluxlist =["B::EIsotope"]
│ set parameters: [config.yaml] target_prefix =flux_
│ set parameters: [Default] flux_totals =false
└ set parameters: [Default] const_stub =false
┌ Info:
│ ================================================================================
│ creating domain 'Box2' ID=3
└ ================================================================================
┌ Info: create_reaction_from_config: Box2.reservoir_B classname ReactionReservoirScalar
│ set parameters: [config.yaml] field_data =external%EIsotope
│ expandvalue: external%EIsotope -> IsotopeLinear
│ after substitution field_data=IsotopeLinear
│ set parameters: [Default] const =false
└ set parameters: [Default] state_norm =false
┌ Info: create_reaction_from_config: Box2.transfer_fluxBoxes classname ReactionFluxTransfer
│ set parameters: [config.yaml] input_fluxes =fluxBoxes.flux_$fluxname$
│ set parameters: [config.yaml] output_fluxes =$fluxname$_sms
│ set parameters: [Default] transfer_multiplier =1.0
└ set parameters: [Default] transfer_matrix =Identity
┌ Info:
│ ================================================================================
│ creating domain 'Box1' ID=4
└ ================================================================================
┌ Info: create_reaction_from_config: Box1.Adecay classname ReactionExample5
│ set parameters: [config.yaml] kappa =0.5
│ set parameters: [config.yaml] field_data =external%EIsotope
│ expandvalue: external%EIsotope -> IsotopeLinear
│ after substitution field_data=IsotopeLinear
└ set parameters: [Default] Delta =-10.0
┌ Info: create_reaction_from_config: Box1.reservoir_A classname ReactionReservoirScalar
│ set parameters: [config.yaml] field_data =external%EIsotope
│ expandvalue: external%EIsotope -> IsotopeLinear
│ after substitution field_data=IsotopeLinear
│ set parameters: [Default] const =false
└ set parameters: [Default] state_norm =false
┌ Info:
│ ================================================================================
│ set_model_geometry
└ ================================================================================
┌ Info:
│ ================================================================================
│ register_reaction_methods!
└ ================================================================================
┌ Info: register_methods: global.sum_E PALEOboxes.VariableStats.ReactionSum
│ add 1 * Box1.A
└ add 1 * Box2.B
┌ Info:
│ ================================================================================
│ link_variables: first pass
└ ================================================================================
┌ Info:
│ ================================================================================
│ link_variables: register_reaction_dynamic_methods and configure variables
└ ================================================================================
[ Info: Reaction global.sum_E Variable global.E_total adding data=PALEOboxes.IsotopeLinear
┌ Info: _configure_variables: ReactionSum global.sum_E variable_links:
└ set variable_links: sum -> E_total
┌ Info: _configure_variables: ReactionReservoirScalar Box2.reservoir_B variable_links:
│ set variable_links: R -> B
│ set variable_links: R_sms -> B_sms
│ set variable_links: R_norm -> B_norm
│ set variable_links: R_delta -> B_delta
│ _configure_variables: ReactionReservoirScalar Box2.reservoir_B variable_attributes:
│ set attribute: R :norm_value = 10.0
└ set attribute: R :initial_value = 0.0
┌ Info: _configure_variables: ReactionExample5 Box1.Adecay variable_links:
│ set variable_links: input_particle -> A
│ set variable_links: input_particle_sms -> A_sms
│ set variable_links: input_particle_delta -> A_delta
└ set variable_links: output_flux -> fluxBoxes.flux_B
┌ Info: _configure_variables: ReactionReservoirScalar Box1.reservoir_A variable_links:
│ set variable_links: R -> A
│ set variable_links: R_sms -> A_sms
│ set variable_links: R_norm -> A_norm
│ set variable_links: R_delta -> A_delta
│ _configure_variables: ReactionReservoirScalar Box1.reservoir_A variable_attributes:
│ set attribute: R :norm_value = 10.0
└ set attribute: R :initial_value = 10.0
┌ Info:
│ ================================================================================
│ link_variables: second pass:
└ ================================================================================
┌ Info:
│ ================================================================================
│ link_variables! unlinked variables:
└ ================================================================================
┌ Info:
│ ================================================================================
│ create_model_from_config: done
└ ================================================================================
┌ Info:
│ ================================================================================
│ PALEOmodel.initialize! start
└ ================================================================================
┌ Info:
│ ================================================================================
│ allocate_variables! (modeldata arrays_idx=1)
└ ================================================================================
[ Info: Domain global data dimensions PALEOboxes.NamedDimension[] allocating 1 variables (hostdep=nothing)
[ Info: Domain fluxBoxes data dimensions PALEOboxes.NamedDimension[] allocating 1 variables (hostdep=nothing)
[ Info: Domain Box2 data dimensions PALEOboxes.NamedDimension[] allocating 4 variables (hostdep=nothing)
[ Info: Domain Box1 data dimensions PALEOboxes.NamedDimension[] allocating 5 variables (hostdep=nothing)
[ Info: set_default_solver_view:
┌ Info: SolverView:
│ host-dependent Variables:
│ ------------------------------------------------------------------------------------------------------------------------------------------------
│ Domain operatorID VF_StateExplicit VF_Total VF_Constraint VF_StateTotal VF_State VF_Undefined
│ global 0 0 0 0 0 0 0
│ fluxBoxes 0 0 0 0 0 0 0
│ Box2 0 1 0 0 0 0 0
│ Box1 0 1 0 0 0 0 0
│ ------------------------------------------------------------------------------------------------------------------------------------------------
│ Total - 2 0 0 0 0 0
│ ------------------------------------------------------------------------------------------------------------------------------------------------
│ n_state_vars 2 (stateexplicit 2 + statetotal 0 + state 0)
│ n_equations 2 (stateexplicit 2 + total 0 + constraint 0)
│ including all host-dependent non-state Variables
└ host-dependent non-state Variables (:vfunction PB.VF_Undefined): Any[]
┌ Info:
│ ================================================================================
│ initialize_reactiondata! (modeldata arrays_indices=1:1)
└ ================================================================================
┌ Info: prepare_do_transfer: ReactionMethod Box2.transfer_fluxBoxes.do_transfer
│ active fluxes:
│ B fluxBoxes.flux_B -> Box2.B_sms
└ using Identity transfer matrix * 1.0
┌ Info:
│ ================================================================================
│ dispatch_setup :setup
└ ================================================================================
┌ Info:
│ ================================================================================
│ dispatch_setup :norm_value
└ ================================================================================
[ Info: init_values! :norm_value Box2.B = 10.0, delta=1.0 (fixed delta value to calculate norm)
[ Info: init_values! :norm_value Box1.A = 10.0, delta=1.0 (fixed delta value to calculate norm)
┌ Info:
│ ================================================================================
│ dispatch_setup :initial_value
└ ================================================================================
[ Info: init_values! :initial_value Box2.B = 0.0, delta=0.0
[ Info: init_values! :initial_value Box1.A = 10.0, delta=0.0
┌ Info:
│ ================================================================================
│ PALEOmodel.initialize! done
└ ================================================================================
initial_state: [0.0, 0.0, 10.0, 0.0]
initial_deriv: [5.0, -50.0, -5.0, 50.0]
integrate, ODE
┌ Info:
│ ================================================================================
│ PALEOmodel.ODE.integrate:
│ tspan: (0.0, 10.0)
│ algorithm: Sundials.CVODE_BDF{:Newton, :Dense, Nothing, Nothing}(0, 0, 0, false, 10, 5, 7, 3, 10, nothing, nothing, 0)
│ Jacobian: NoJacobian
│ steadystate: false
│ using 1 BLAS threads
└ ================================================================================
[ Info: ODEfunction: using Jacobian NoJacobian
[ Info: jac_config_ode: jac_ad=NoJacobian
0.573460 seconds (1.12 M allocations: 77.735 MiB, 5.43% gc time, 99.95% compilation time)
┌ Info:
│ ================================================================================
│ print_sol_stats:
│ retcode=Success
│ alg=Sundials.CVODE_BDF{:Newton, :Dense, Nothing, Nothing}(0, 0, 0, false, 10, 5, 7, 3, 10, nothing, nothing, 0)
│ stats=SciMLBase.DEStats
│ Number of function 1 evaluations: 80
│ Number of function 2 evaluations: 0
│ Number of W matrix evaluations: 20
│ Number of linear solves: 0
│ Number of Jacobians created: 2
│ Number of nonlinear solver iterations: 77
│ Number of nonlinear solver convergence failures: 0
│ Number of fixed-point solver iterations: 0
│ Number of fixed-point solver convergence failures: 0
│ Number of rootfind condition calls: 0
│ Number of accepted steps: 66
│ Number of rejected steps: 0
│ length(sol.t) 67
│ size(sol) (4, 67)
└ ================================================================================
[ Info: ODE.calc_output_sol!: 67 records
[ Info: ODE.calc_output_sol!: done
┌ Info:
│ ================================================================================
│ PALEOmodel.ODE.integrate: done
└ ================================================================================
0.789910 seconds (1.32 M allocations: 92.864 MiB, 3.94% gc time, 99.47% compilation time)
Displaying model structure and variables
All metadata for model variables can be displayed with PB.show_variables:
show(PB.show_variables(model), allcols=true, allrows=true) # display in REPL
# vscodedisplay(PB.show_variables(model)) # more convenient when using VS code11×9 DataFrame
Row │ domain name type units vfunction space data_dims field_data description
│ String String DataType String Variable… DataType Tuple{} Type String
─────┼────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
1 │ Box1 A VariableDomPropDep mol VF_StateExplicit ScalarSpace () IsotopeLinear scalar reservoir
2 │ Box1 A_delta VariableDomPropDep per mil VF_Undefined ScalarSpace () ScalarData scalar reservoir isotope delta
3 │ Box1 A_norm VariableDomPropDep VF_Undefined ScalarSpace () ScalarData scalar reservoir normalized
4 │ Box1 A_sms VariableDomContribTarget mol yr-1 VF_Deriv ScalarSpace () IsotopeLinear scalar reservoir source-sinks
5 │ Box1 decay_flux VariableDomPropDep mol yr-1 VF_Undefined ScalarSpace () IsotopeLinear decay flux
6 │ Box2 B VariableDomPropDep mol VF_StateExplicit ScalarSpace () IsotopeLinear scalar reservoir
7 │ Box2 B_delta VariableDomPropDep per mil VF_Undefined ScalarSpace () ScalarData scalar reservoir isotope delta
8 │ Box2 B_norm VariableDomPropDep VF_Undefined ScalarSpace () ScalarData scalar reservoir normalized
9 │ Box2 B_sms VariableDomContribTarget mol yr-1 VF_Deriv ScalarSpace () IsotopeLinear scalar reservoir source-sinks
10 │ fluxBoxes flux_B VariableDomContribTarget mol yr-1 VF_Undefined CellSpace () IsotopeLinear flux input
11 │ global E_total VariableDomPropDep unknown VF_Undefined ScalarSpace () IsotopeLinear sum of specified variablesThe variable names are as in Example 4 Transfer between Domains, however field_data is now IsotopeLinear and not ScalarData for the reservoir and flux variables.