Configurable chemistry using PALEOaqchem
These examples illustrate how to use the yaml configuration file to configure generic chemistry Reactions from PALEOaqchem to implement the same minimal model of the marine carbonate system shown in the Ocean chemistry examples.
See PALEOaqchem documentation for more detail on the generic aqueous chemistry.
See PALEOmodel documentation for more information on analysing model output.
Example 1 Configuring a minimal Alk-pH model
This example configures generic chemistry reactions from PALEOaqchem to define a minimal Alk-pH model for aqueous carbonate chemistry. Carbonate system equations are from (Zeebe and Wolf-Gladrow, 2001).
The ocean domain contains a ReactionNoTransport from the PALEOocean Julia package. This defines standard ocean variables including cell volume, and is configured here to provide one ocean cell.
There are two state variables for DIC and TAlk, implemented using a ReactionReservoir from the PALEOboxes Julia package. This Reaction also provides concentrations in mol m-3. Constant Boron concentration B_conc is defined using a ReactionReservoirConst.
The carbonate chemistry system includes the speciation of water, boron and DIC. B, DIC, and TAlk are the three total variables with corresponding primary species concentrations BOH3_conc, CO2_aq_conc, and H_conc (defined by pH). Primary species are defined using the PALEOaqchem ReactionConstraintReservoir, which defines a state variable (attribute :vfunction = PB.VF_State) eg BOH3_conc, an algebraic constraint for the corresponding total eg B_constraint_conc (with attribute :vfunction = PB.VF_Constraint), and a variable eg B_calc to accumulate primary and secondary species contributions. Secondary species eg BOH4 are defined using the PALEOaqchem ReactionAqEqb, and their contribution added to the calculated total eg B_calc. The PALEO DAE solver then solves for the primary species concentration BOH3_conc that makes B_calc equal to the required value B. Note that the ReactionConstraintReservoir for TAlk defines pH as the state variable (from which H_conc is calculated), and that many species contribute to both TAlk_calc and their own total.
The combined system is then a Differential Algebraic Equation (DAE) system with two ODE state variables and corresponding time derivatives for DIC and TAlk, three algebraic constraints for B, DIC, and TAlk, and three additional state variables for the primary species BOH3_conc, CO2_aq_conc, and pH, and is integrated forward in time by a DAE solver PALEOmodel.ODE.integrateDAE.
yaml configuration file
The model configuration (file examples/configurable_chemistry/config_ex1.yaml) contains two Reservoirs DIC, TAlk.
A ReactionFluxPerturb from the PALEOboxes.jl reaction catalog is used to add a constant TAlk flux.
Minimal_Alk_pH:
domains:
global:
reactions:
add_Alk:
# from PALEOboxes reaction catalog: constant input flux of Alk
class: ReactionFluxPerturb
parameters:
# linear interpolation for input flux - set to constant
perturb_times: [-1e30, 1e30]
# ocean volume = 3.6e14 * 3688.0 = 1.33e18 m^3
# DIC in ocean = 2 * 1.33e18 = 2.66e18 mol
# set Alk flux so Alk = DIC after 100 yr-1 = 2.66e18/100 = 2.66e16 mol yr-1
perturb_totals: [2.66e16, 2.66e16]
variable_links:
F: ocean.TAlk_sms # add flux to TAlk reservoir NB: works for single-box ocean only !!
# --------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
ocean:
reactions:
oceantrasport:
class: ReactionOceanNoTransport
parameters:
area: [3.6e14]
depth: [3688.0]
reservoir_TA:
class: ReactionReservoir
variable_links:
R*: TAlk*
variable_attributes:
R:initial_value: 0.0
R:norm_value: 2.4e0
reservoir_DIC:
class: ReactionReservoir
variable_links:
R*: DIC*
variable_attributes:
R:initial_value: 2.0e0
R:norm_value: 2.0e0
constant_B:
class: ReactionReservoirConst
variable_links:
R*: B
variable_attributes:
R_conc:initial_value: 0.427 # mol m-3 contemporary value
##########################################################################
# TAlk, primary species H_conc (as pH)
########################################################################
# [H+] primary species (as pH) for TAlk total
H_primary_species: {class: ReactionConstraintReservoir,
variable_links: {Primary_pconc: pH, Primary_conc: H_conc, R*: TAlk*},
parameters: {primary_total_stoich: -1.0, primary_variable: p_concentration, constraint_variable: concentration},
variable_attributes: {Primary_pconc%initial_value: 7.0, Primary_pconc%norm_value: 1.0, R_constraint_conc%norm_value: 1e-3}}
# Define OH_conc and add to TAlk
# OH- + H+ <--> H2O, [OH] = K_w / [H]
OH_conc: { class: ReactionAqEqb, parameters: { K_power: 1.0, # we need K_eqb in numerator
K_eqb: 6.0e-14, # K_w 6.0e-14 (mol^2 kg-2) equilibrium constant of water at S=35, T=25°C,
K_density_power: 2.0, # K_eqb*rho_ref^2 to convert to (mol^2 m-6)
Reactants: ["OH_conc", "H_conc"], Products: [],
N_components: ["TAlk_calc"],
}
}
##########################################################################
# Boron speciation
##########################################################################
# BOH3 primary species for B total
BOH3_primary_species: {class: ReactionConstraintReservoir,
variable_links: {Primary*: BOH3*, R*: B*},
parameters: {primary_total_stoich: 1.0, primary_variable: concentration, constraint_variable: concentration},
variable_attributes: {Primary_conc%initial_value: 1.0, Primary_conc%norm_value: 1e-3, R_constraint_conc%norm_value: 1e-3}}
# Define BOH4_conc and add to B, TAlk
# B(OH)_4^- + H+ <--> B(OH)3 + H2O, [B(OH)_4^-] = K_B [B(OH)3] / [H+]
BOH4_conc: { class: ReactionAqEqb, parameters: { K_power: 1.0, # we need K_eqb in numerator
K_eqb: 2.5e-9, # (mol kg-1) K_B equilibrium constant of B(OH)4-
K_density_power: 1.0, # K_eqb*rho_ref to convert to (mol m-3)
Reactants: ["BOH4_conc", "H_conc"], Products: ["BOH3_conc"],
N_components: ["B_calc", "TAlk_calc"],
}
}
##########################################################################
# Carbon speciation
##########################################################################
# CO2_aq primary species for DIC total
CO2_aq_primary_species: {class: ReactionConstraintReservoir,
variable_links: {Primary*: CO2_aq*, R*: DIC*},
parameters: {primary_total_stoich: 1.0, primary_variable: concentration, constraint_variable: concentration},
variable_attributes: {Primary_conc%initial_value: 1.0, Primary_conc%norm_value: 1e-3, R_constraint_conc%norm_value: 1.0}}
# Define HCO3_conc and add to DIC, TAlk
# HCO3- + H+ <--> CO2 + H2O, [HCO3-] = K_1 [CO2] / [H+]
HCO3_conc: { class: ReactionAqEqb, parameters: { K_power: 1.0, # we need K_eqb in numerator
K_eqb: 1.4e-6, # (mol kg-1) K_1 equilibrium constant of CO2_aq and HCO3-
K_density_power: 1.0, # K_eqb*rho_ref to convert to (mol m-3)
Reactants: ["HCO3_conc", "H_conc"], Products: ["CO2_aq_conc"],
N_components: ["DIC_calc", "TAlk_calc"],
}
}
# Define CO3_conc and add to DIC, TAlk
# CO3-- + H+ <--> HCO3- + H2O, [CO3--] = K_2 [HCO3-] / [H+]
CO3_conc: { class: ReactionAqEqb, parameters: { K_power: 1.0, # we need K_eqb in numerator
K_eqb: 1.2e-9, # (mol kg-1) K_2 equilibrium constant of HCO3- and CO32-
K_density_power: 1.0, # K_eqb*rho_ref to convert to (mol m-3)
Reactants: ["CO3_conc", "H_conc"], Products: ["HCO3_conc"],
N_components: ["DIC_calc", "2*TAlk_calc"],
}
}
oceanfloor:
# unused here, set up by ReactionOceanNoTransport
oceansurface:
# unused here, set up by ReactionOceanNoTransportRun script
The script to run the model (file examples/configurable_chemistry/run_ex1.jl) contains:
import PALEOboxes as PB
import PALEOmodel
import PALEOocean
import PALEOaqchem
using Plots
#####################################################
# Create model
model = PB.create_model_from_config(joinpath(@__DIR__, "config_ex1.yaml"), "Minimal_Alk_pH")
#########################################################
# 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, DAE")
# first run is slow as it includes JIT time
@time PALEOmodel.ODE.integrateDAE(
paleorun, initial_state, modeldata, (0.0, 200.0),
solvekwargs=(
reltol=1e-6,
# 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, ["ocean.TAlk_conc", "ocean.DIC_conc"], (cell=1,);
ylabel="TAlk, DIC conc (mol m-3)"))
display(plot(paleorun.output, ["ocean.TAlk","ocean.TAlk_sms"], (cell=1,)))
display(plot(paleorun.output, "ocean.pH", (cell=1,)))
display(plot(paleorun.output, ["ocean.H_conc", "ocean.OH_conc"], (cell=1,);
ylabel="H2O species (mol m-3)", ylim=(0, Inf)))
display(plot(paleorun.output, ["ocean.B_conc", "ocean.BOH4_conc", "ocean.BOH3_conc"], (cell=1,);
ylabel="B species (mol m-3)", ylim=(0, Inf)))
display(plot(paleorun.output, ["ocean.DIC_conc", "ocean.HCO3_conc", "ocean.CO3_conc", "ocean.CO2_aq_conc"], (cell=1,);
ylabel="DIC species (mol m-3)", ylim=(0, Inf)))
display(
plot(
paleorun.output,
["ocean.HCO3_conc", "ocean.CO3_conc", "ocean.CO2_aq_conc", "ocean.BOH4_conc", "ocean.BOH3_conc", "ocean.H_conc", "ocean.OH_conc",],
(cell=1, );
coords=["tmodel"=>("ocean.pH",),], # plot against pH instead of tmodel
ylabel="H2O, B, DIC species (mol m-3)", ylim=(0.5e-3, 0.5e1), yscale=:log10,
legend_background_color=nothing,
legend=:bottom,
)
)
and produces output showing the change, if TAlk_conc increase, how the carbonic acid and pH change in ocean:
[ Info: OceanTransportTMM adding default TMMDir => /home/runner/.julia/packages/PALEOocean/TMM to [PALEOocean.Ocean.OceanTransportTMM] LocalPreferences.toml (modify for your local setup)
[ Info: OceanTransportTMM adding default TMMDir => /home/runner/.julia/packages/PALEOocean/TMM to [PALEOboxes] LocalPreferences.toml (modify for your local setup)
┌ Info:
│ ================================================================================
│ create_model_from_config: configmodel Minimal_Alk_pH
│ config_file: /home/runner/work/PALEOtutorials.jl/PALEOtutorials.jl/examples/configurable_chemistry/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.add_Alk classname ReactionFluxPerturb
│ set parameters: [Default] field_data =PALEOboxes.ScalarData
│ set parameters: [config.yaml] perturb_times =[-1.0e30, 1.0e30]
│ set parameters: [config.yaml] perturb_totals =[2.66e16, 2.66e16]
│ set parameters: [Default] perturb_deltas =[0.0, 0.0]
│ set parameters: [Default] extrapolate =throw
│ set parameters: [Default] extrapolate_before =NaN
│ set parameters: [Default] extrapolate_before_delta=NaN
│ set parameters: [Default] extrapolate_after =NaN
└ set parameters: [Default] extrapolate_after_delta=NaN
┌ Info:
│ ================================================================================
│ creating domain 'ocean' ID=2
└ ================================================================================
┌ Info: create_reaction_from_config: ocean.CO3_conc classname ReactionAqEqb
│ set parameters: [config.yaml] Reactants =["CO3_conc", "H_conc"]
│ set parameters: [config.yaml] Products =["HCO3_conc"]
│ set parameters: [config.yaml] K_eqb =1.2e-9
│ set parameters: [config.yaml] K_density_power =1.0
│ set parameters: [config.yaml] K_power =1.0
└ set parameters: [config.yaml] N_components =["DIC_calc", "2*TAlk_calc"]
┌ Info: create_reaction_from_config: ocean.reservoir_DIC classname ReactionReservoir
│ set parameters: [Default] field_data =PALEOboxes.ScalarData
│ set parameters: [Default] total =false
│ set parameters: [Default] limit_delta_conc =0.0
└ set parameters: [Default] state_conc =false
┌ Info: create_reaction_from_config: ocean.reservoir_TA classname ReactionReservoir
│ set parameters: [Default] field_data =PALEOboxes.ScalarData
│ set parameters: [Default] total =false
│ set parameters: [Default] limit_delta_conc =0.0
└ set parameters: [Default] state_conc =false
┌ Info: create_reaction_from_config: ocean.oceantrasport classname ReactionOceanNoTransport
│ set parameters: [config.yaml] area =[3.6e14]
└ set parameters: [config.yaml] depth =[3688.0]
┌ Info: create_reaction_from_config: ocean.HCO3_conc classname ReactionAqEqb
│ set parameters: [config.yaml] Reactants =["HCO3_conc", "H_conc"]
│ set parameters: [config.yaml] Products =["CO2_aq_conc"]
│ set parameters: [config.yaml] K_eqb =1.4e-6
│ set parameters: [config.yaml] K_density_power =1.0
│ set parameters: [config.yaml] K_power =1.0
└ set parameters: [config.yaml] N_components =["DIC_calc", "TAlk_calc"]
┌ Info: create_reaction_from_config: ocean.BOH3_primary_species classname ReactionConstraintReservoir
│ set parameters: [config.yaml] primary_total_stoich=1.0
│ set parameters: [Default] primary_other_components=String[]
│ set parameters: [config.yaml] primary_variable =concentration
└ set parameters: [config.yaml] constraint_variable =concentration
┌ Info: create_reaction_from_config: ocean.constant_B classname ReactionReservoirConst
└ set parameters: [Default] field_data =PALEOboxes.ScalarData
┌ Info: create_reaction_from_config: ocean.BOH4_conc classname ReactionAqEqb
│ set parameters: [config.yaml] Reactants =["BOH4_conc", "H_conc"]
│ set parameters: [config.yaml] Products =["BOH3_conc"]
│ set parameters: [config.yaml] K_eqb =2.5e-9
│ set parameters: [config.yaml] K_density_power =1.0
│ set parameters: [config.yaml] K_power =1.0
└ set parameters: [config.yaml] N_components =["B_calc", "TAlk_calc"]
┌ Info: create_reaction_from_config: ocean.CO2_aq_primary_species classname ReactionConstraintReservoir
│ set parameters: [config.yaml] primary_total_stoich=1.0
│ set parameters: [Default] primary_other_components=String[]
│ set parameters: [config.yaml] primary_variable =concentration
└ set parameters: [config.yaml] constraint_variable =concentration
┌ Info: create_reaction_from_config: ocean.OH_conc classname ReactionAqEqb
│ set parameters: [config.yaml] Reactants =["OH_conc", "H_conc"]
│ set parameters: [config.yaml] Products =Any[]
│ set parameters: [config.yaml] K_eqb =6.0e-14
│ set parameters: [config.yaml] K_density_power =2.0
│ set parameters: [config.yaml] K_power =1.0
└ set parameters: [config.yaml] N_components =["TAlk_calc"]
┌ Info: create_reaction_from_config: ocean.H_primary_species classname ReactionConstraintReservoir
│ set parameters: [config.yaml] primary_total_stoich=-1.0
│ set parameters: [Default] primary_other_components=String[]
│ set parameters: [config.yaml] primary_variable =p_concentration
└ set parameters: [config.yaml] constraint_variable =concentration
┌ Info:
│ ================================================================================
│ creating domain 'oceansurface' ID=3
└ ================================================================================
┌ Info:
│ ================================================================================
│ creating domain 'oceanfloor' ID=4
└ ================================================================================
┌ Info:
│ ================================================================================
│ set_model_geometry
└ ================================================================================
[ Info: set ocean.oceansurface Subdomain size=1
[ Info: set ocean.oceanfloor Subdomain size=1
[ Info: Ocean.set_model_domains:
[ Info: ocean Domain size=1 grid=UnstructuredVectorGrid(ncells=1, cellnames=Dict{Symbol, Int64}(), subdomains: ["oceansurface", "oceanfloor"])
[ Info: oceansurface Domain size=1 grid=UnstructuredVectorGrid(ncells=1, cellnames=Dict{Symbol, Int64}(), subdomains: ["ocean"])
[ Info: oceanfloor Domain size=1 grid=UnstructuredVectorGrid(ncells=1, cellnames=Dict{Symbol, Int64}(), subdomains: ["ocean"])
┌ Info:
│ ================================================================================
│ register_reaction_methods!
└ ================================================================================
[ Info: ReactionFluxPerturb.register_methods! global.add_Alk
[ Info: register_methods! ocean.CO3_conc
[ Info: register_methods! ocean.HCO3_conc
[ Info: register_methods! ocean.BOH4_conc
[ Info: register_methods! ocean.OH_conc
┌ Info:
│ ================================================================================
│ link_variables: first pass
└ ================================================================================
┌ Info:
│ ================================================================================
│ link_variables: register_reaction_dynamic_methods and configure variables
└ ================================================================================
┌ Info: _configure_variables: ReactionFluxPerturb global.add_Alk variable_links:
└ set variable_links: F -> ocean.TAlk_sms
┌ Info: _configure_variables: ReactionReservoir ocean.reservoir_DIC variable_links:
│ set variable_links: R -> DIC
│ set variable_links: R_sms -> DIC_sms
│ set variable_links: R_conc -> DIC_conc
│ _configure_variables: ReactionReservoir ocean.reservoir_DIC variable_attributes:
│ set attribute: R :norm_value = 2.0
└ set attribute: R :initial_value = 2.0
┌ Info: _configure_variables: ReactionReservoir ocean.reservoir_TA variable_links:
│ set variable_links: R -> TAlk
│ set variable_links: R_sms -> TAlk_sms
│ set variable_links: R_conc -> TAlk_conc
│ _configure_variables: ReactionReservoir ocean.reservoir_TA variable_attributes:
│ set attribute: R :norm_value = 2.4
└ set attribute: R :initial_value = 0.0
┌ Info: _configure_variables: ReactionConstraintReservoir ocean.BOH3_primary_species variable_links:
│ set variable_links: Primary_conc -> BOH3_conc
│ set variable_links: R_constraint_conc -> B_constraint_conc
│ set variable_links: R_calc -> B_calc
│ set variable_links: R_conc -> B_conc
│ _configure_variables: ReactionConstraintReservoir ocean.BOH3_primary_species variable_attributes:
│ set attribute: Primary_conc :initial_value = 1.0
│ set attribute: Primary_conc :norm_value = 0.001
└ set attribute: R_constraint_conc :norm_value = 0.001
┌ Info: _configure_variables: ReactionReservoirConst ocean.constant_B variable_links:
│ set variable_links: R_conc -> B_conc
│ _configure_variables: ReactionReservoirConst ocean.constant_B variable_attributes:
└ set attribute: R_conc :initial_value = 0.427
┌ Info: _configure_variables: ReactionConstraintReservoir ocean.CO2_aq_primary_species variable_links:
│ set variable_links: Primary_conc -> CO2_aq_conc
│ set variable_links: R_constraint_conc -> DIC_constraint_conc
│ set variable_links: R_calc -> DIC_calc
│ set variable_links: R_conc -> DIC_conc
│ _configure_variables: ReactionConstraintReservoir ocean.CO2_aq_primary_species variable_attributes:
│ set attribute: Primary_conc :initial_value = 1.0
│ set attribute: Primary_conc :norm_value = 0.001
└ set attribute: R_constraint_conc :norm_value = 1.0
┌ Info: _configure_variables: ReactionConstraintReservoir ocean.H_primary_species variable_links:
│ set variable_links: R_constraint_conc -> TAlk_constraint_conc
│ set variable_links: R_calc -> TAlk_calc
│ set variable_links: R_conc -> TAlk_conc
│ set variable_links: Primary_conc -> H_conc
│ set variable_links: Primary_pconc -> pH
│ _configure_variables: ReactionConstraintReservoir ocean.H_primary_species variable_attributes:
│ set attribute: Primary_pconc :initial_value = 7.0
│ set attribute: R_constraint_conc :norm_value = 0.001
└ set attribute: Primary_pconc :norm_value = 1.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 2 variables (hostdep=nothing)
[ Info: set :field_data=PALEOboxes.ScalarData for host-dependent Variable global.tforce
[ Info: Domain ocean data dimensions PALEOboxes.NamedDimension[] allocating 29 variables (hostdep=nothing)
[ Info: Domain oceansurface data dimensions PALEOboxes.NamedDimension[] allocating 1 variables (hostdep=nothing)
[ Info: Domain oceanfloor 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 0 0 0 0 0 1
│ ocean 0 2 0 3 0 3 0
│ oceansurface 0 0 0 0 0 0 0
│ oceanfloor 0 0 0 0 0 0 0
│ ------------------------------------------------------------------------------------------------------------------------------------------------
│ Total - 2 0 3 0 3 1
│ ------------------------------------------------------------------------------------------------------------------------------------------------
│ n_state_vars 5 (stateexplicit 2 + statetotal 0 + state 3)
│ n_equations 5 (stateexplicit 2 + total 0 + constraint 3)
│ including all host-dependent non-state Variables
└ host-dependent non-state Variables (:vfunction PB.VF_Undefined): ["global.tforce"]
┌ Info:
│ ================================================================================
│ initialize_reactiondata! (modeldata arrays_indices=1:1)
└ ================================================================================
┌ Info:
│ ================================================================================
│ dispatch_setup :setup
└ ================================================================================
[ Info: ocean.constant_B.setup_initialvalue_vars_default:
[ Info: init_values! :initial_value ocean.B_conc = 0.427
┌ Info:
│ ================================================================================
│ dispatch_setup :norm_value
└ ================================================================================
[ Info: ocean.BOH3_primary_species.setup_initialvalue_vars_default:
[ Info: init_values! :norm_value ocean.BOH3_conc = 0.001
[ Info: init_values! :norm_value ocean.B_constraint_conc = 0.001
[ Info: ocean.CO2_aq_primary_species.setup_initialvalue_vars_default:
[ Info: init_values! :norm_value ocean.CO2_aq_conc = 0.001
[ Info: init_values! :norm_value ocean.DIC_constraint_conc = 1.0
[ Info: ocean.H_primary_species.setup_initialvalue_vars_default:
[ Info: init_values! :norm_value ocean.pH = 1.0
[ Info: init_values! :norm_value ocean.TAlk_constraint_conc = 0.001
[ Info: ocean.reservoir_DIC.setup_initialvalue_vars_default:
[ Info: init_values! :norm_value ocean.DIC = 2.0 * volume
[ Info: ocean.reservoir_TA.setup_initialvalue_vars_default:
[ Info: init_values! :norm_value ocean.TAlk = 2.4 * volume
┌ Info:
│ ================================================================================
│ dispatch_setup :initial_value
└ ================================================================================
[ Info: ocean.BOH3_primary_species.setup_initialvalue_vars_default:
[ Info: init_values! :initial_value ocean.BOH3_conc = 1.0
[ Info: init_values! :initial_value ocean.B_constraint_conc = 0.0
[ Info: ocean.CO2_aq_primary_species.setup_initialvalue_vars_default:
[ Info: init_values! :initial_value ocean.CO2_aq_conc = 1.0
[ Info: init_values! :initial_value ocean.DIC_constraint_conc = 0.0
[ Info: ocean.H_primary_species.setup_initialvalue_vars_default:
[ Info: init_values! :initial_value ocean.pH = 7.0
[ Info: init_values! :initial_value ocean.TAlk_constraint_conc = 0.0
[ Info: ocean.reservoir_DIC.setup_initialvalue_vars_default:
[ Info: init_values! :initial_value ocean.DIC = 2.0 * volume
[ Info: ocean.reservoir_TA.setup_initialvalue_vars_default:
[ Info: init_values! :initial_value ocean.TAlk = 0.0 * volume
┌ Info:
│ ================================================================================
│ PALEOmodel.initialize! done
└ ================================================================================
initial_state: [0.0, 2.65536e18, 1.0, 1.0, 7.0]
initial_deriv: [2.66e16, 0.0, -0.5979999999999999, -13.168, -14.361513500000001]
integrate, DAE
┌ Info:
│ ================================================================================
│ PALEOmodel.ODE.integrateDAE:
│ tspan: (0.0, 200.0)
│ algorithm: Sundials.IDA{:Dense, Nothing, Nothing}(0, 0, 0, 5, 7, 0.33, 3, 10, 0.0033, 5, 4, 10, 100, true, false, nothing, nothing)
│ Jacobian: NoJacobian
│ using 1 BLAS threads
└ ================================================================================
[ Info: DAEfunction: using Jacobian NoJacobian
[ Info: jac_config_dae: jac_ad=NoJacobian
[ Info: calling get_inconsistent_initial_deriv
0.940230 seconds (1.08 M allocations: 73.833 MiB, 3.89% gc time, 99.85% compilation time)
┌ Info:
│ ================================================================================
│ print_sol_stats:
│ retcode=Success
│ alg=Sundials.IDA{:Dense, Nothing, Nothing}(0, 0, 0, 5, 7, 0.33, 3, 10, 0.0033, 5, 4, 10, 100, true, false, nothing, nothing)
│ stats=SciMLBase.DEStats
│ Number of function 1 evaluations: 547
│ Number of function 2 evaluations: 0
│ Number of W matrix evaluations: 109
│ Number of linear solves: 0
│ Number of Jacobians created: 109
│ Number of nonlinear solver iterations: 545
│ 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: 292
│ Number of rejected steps: 14
│ length(sol.t) 307
│ size(sol) (5, 307)
└ ================================================================================
[ Info: ODE.calc_output_sol!: 307 records
[ Info: ODE.calc_output_sol!: done
┌ Info:
│ ================================================================================
│ PALEOmodel.ODE.integrateDAE: done
└ ================================================================================
4.682510 seconds (6.13 M allocations: 421.506 MiB, 3.76% gc time, 99.71% 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 code35×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 add_Alk/FApplied VariableDomPropDep mol yr-1 VF_Undefined ScalarSpace () ScalarData flux perturbation applied, for d…
2 │ global tforce VariableDomPropDep yr VF_Undefined ScalarSpace () ScalarData time at which to apply perturbat…
3 │ ocean Abox VariableDomPropDep m^2 VF_Undefined CellSpace () ScalarData horizontal area of box
4 │ ocean BOH3_conc VariableDomPropDep mol m-3 VF_State CellSpace () ScalarData concentration of primary species
5 │ ocean BOH4_conc VariableDomPropDep mol m-3 VF_Undefined CellSpace () ScalarData aqueous concentration or activity
6 │ ocean B_calc VariableDomContribTarget mol VF_Undefined CellSpace () ScalarData contributions to total R_calc_co…
7 │ ocean B_conc VariableDomPropDep mol m-3 VF_Undefined CellSpace () ScalarData concentration
8 │ ocean B_constraint_conc VariableDomContribTarget mol m-3 VF_Constraint CellSpace () ScalarData algebraic constraint on R_conc (…
9 │ ocean CO2_aq_conc VariableDomPropDep mol m-3 VF_State CellSpace () ScalarData concentration of primary species
10 │ ocean CO3_conc VariableDomPropDep mol m-3 VF_Undefined CellSpace () ScalarData aqueous concentration or activity
11 │ ocean DIC VariableDomPropDep mol VF_StateExplicit CellSpace () ScalarData vector reservoir
12 │ ocean DIC_calc VariableDomContribTarget mol VF_Undefined CellSpace () ScalarData contributions to total R_calc_co…
13 │ ocean DIC_conc VariableDomPropDep mol m-3 VF_Undefined CellSpace () ScalarData concentration
14 │ ocean DIC_constraint_conc VariableDomContribTarget mol m-3 VF_Constraint CellSpace () ScalarData algebraic constraint on R_conc (…
15 │ ocean DIC_sms VariableDomContribTarget mol yr-1 VF_Deriv CellSpace () ScalarData vector reservoir source-sinks
16 │ ocean HCO3_conc VariableDomPropDep mol m-3 VF_Undefined CellSpace () ScalarData aqueous concentration or activity
17 │ ocean H_conc VariableDomPropDep mol m-3 VF_Undefined CellSpace () ScalarData concentration of primary species
18 │ ocean OH_conc VariableDomPropDep mol m-3 VF_Undefined CellSpace () ScalarData aqueous concentration or activity
19 │ ocean TAlk VariableDomPropDep mol VF_StateExplicit CellSpace () ScalarData vector reservoir
20 │ ocean TAlk_calc VariableDomContribTarget mol VF_Undefined CellSpace () ScalarData contributions to total R_calc_co…
21 │ ocean TAlk_conc VariableDomPropDep mol m-3 VF_Undefined CellSpace () ScalarData concentration
22 │ ocean TAlk_constraint_conc VariableDomContribTarget mol m-3 VF_Constraint CellSpace () ScalarData algebraic constraint on R_conc (…
23 │ ocean TAlk_sms VariableDomContribTarget mol yr-1 VF_Deriv CellSpace () ScalarData vector reservoir source-sinks
24 │ ocean pH VariableDomPropDep VF_State CellSpace () ScalarData -log10(concentration of primary …
25 │ ocean pressure VariableDomPropDep dbar VF_Undefined CellSpace () ScalarData Ocean pressure
26 │ ocean rho_ref VariableDomPropDep kg m^-3 VF_Undefined CellSpace () ScalarData Ocean transport model density co…
27 │ ocean volume VariableDomPropDep m^3 VF_Undefined CellSpace () ScalarData volume of ocean cells
28 │ ocean volume_total VariableDomPropDep m^3 VF_Undefined ScalarSpace () ScalarData total volume of ocean cells
29 │ ocean zlower VariableDomPropDep m VF_Undefined CellSpace () ScalarData depth of lower surface of box (m)
30 │ ocean zmid VariableDomPropDep m VF_Undefined CellSpace () ScalarData mean depth of box
31 │ ocean zupper VariableDomPropDep m VF_Undefined CellSpace () ScalarData depth of upper surface of box (m…
32 │ oceanfloor Afloor VariableDomPropDep m^2 VF_Undefined CellSpace () ScalarData horizontal area of seafloor at b…
33 │ oceanfloor Afloor_total VariableDomPropDep m^2 VF_Undefined ScalarSpace () ScalarData total area of seafloor
34 │ oceanfloor zfloor VariableDomPropDep m VF_Undefined CellSpace () ScalarData depth of ocean floor (m, -ve)
35 │ oceansurface Asurf VariableDomPropDep m^2 VF_Undefined CellSpace () ScalarData horizontal area of oceansurfaceExample 2 Minimal Alk-pH model with implicit variables
This example contains the same chemistry as Example 1 Configuring a minimal Alk-pH model, using the Direct Substitution Approach (DSA) to eliminate the two algebraic equations corresponding to the mass action equations for DIC and TAlk.
The state variables for DIC_conc and TAlk_conc are now implicit total variables defined by PALEOaqchem ReactionImplicitReservoirs. These define a PALEO Total variable (attribute :vfunction = PB.VF_Total) eg DIC_conc which is a function of a primary species concentration implemented as a PALEO StateTotal variable (attribute :vfunction = PB.VF_StateTotal) eg CO_2_aq_conc, and a time derivative eg DIC_conc_sms implemented as a PALEO Deriv variable (attribute :vfunction = PB.VF_Deriv).
The combined system is then a Differential Algebraic Equation (DAE) system with three implicit total variables for DIC_conc, B_conc, and TAlk_conc, with corresponding time derivatives DIC_conc_sms, B_conc_sms, and TAlk_conc_sms and state variables for the primary species CO2_aq_conc, BOH3_conc and pH, and is integrated forward in time by a DAE solver PALEOmodel.ODE.integrateDAE.
Notes:
- The combination of constant
BOH3_concand aReactionConstraintReservoirdefining primary speciesBOH3_conchas been replaced with an implicit state variable forB_concand primary speciesBOH3_conc: this is required to use the IDA DAE solver, which cannot handle initialisation of that combination of implicit and constraint variables (seePALEOmodeldocumentation). - Initial values are now specified for the primary species (eg
CO2_conc,BOH3_conc,pH), not the corresponding totals (egDIC_conc). This is less convenient here (where we are want to setDIC_concandB_concthat remain constant), but may not be an issue in a more complete model that includes additional processes that define time evolution that is insensitive to initial conditions after a spinup. - This configuration provides concentrations, instead of mol (per cell) to the solver. This is not required here, but is good practice as it helps the numerical solver by providing variables with better scaling.
yaml configuration file
The model configuration (file examples/configurable_chemistry/config_ex2.yaml) contains:
Minimal_Alk_pH:
domains:
global:
reactions:
add_Alk:
# from PALEOboxes reaction catalog: constant input flux of Alk
class: ReactionFluxPerturb
parameters:
# linear interpolation for input flux - set to constant
perturb_times: [-1e30, 1e30]
# ocean volume = 3.6e14 * 3688.0 = 1.33e18 m^3
# DIC in ocean = 2 * 1.33e18 = 2.66e18 mol
# set Alk flux so Alk = DIC after 100 yr-1 = 2.66e18/100 = 2.66e16 mol yr-1
perturb_totals: [2.66e16, 2.66e16]
variable_links:
F: ocean.TAlk_sms # add flux to TAlk reservoir NB: works for single-box ocean only !!
# --------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
ocean:
reactions:
oceantransport:
class: ReactionOceanNoTransport
parameters:
area: [3.6e14]
depth: [3688.0]
# TA as total and H (as pH) as primary species
TA_total_H_primary:
class: ReactionImplicitReservoir
variable_links:
Primary_pconc: pH
Primary_conc: H_conc
R*: TAlk*
parameters:
primary_total_stoich: -1.0
primary_variable: p_concentration # provide solver with -log10(H_conc)
total_variable: concentration # provide solver with TAlk_conc
variable_attributes:
Primary_pconc%initial_value: 4.28 # value tuned to give TAlk ~ 0.0 initially
Primary_pconc%norm_value: 1.0
R_conc%norm_value: 1.0
# DIC as total and CO2_aq as primary species
DIC_total_CO2_aq_primary:
class: ReactionImplicitReservoir
variable_links:
R*: DIC*
Primary*: CO2_aq*
parameters:
primary_total_stoich: 1.0 # add CO2 to DIC, no contribution to TAlk
primary_variable: concentration # provide solver with CO2_aq_conc (mol m-3)
total_variable: concentration # provide solver with DIC_conc (mol m-3)
variable_attributes:
Primary_conc%initial_value: 1.94 # This value has been tuned to produce initial DIC_conc = 2.0
Primary_conc%norm_value: 1.0
R_conc%norm_value: 1.0
# B as total and BOH3 as primary species
B_total_BOH3_primary:
class: ReactionImplicitReservoir
variable_links:
R*: B*
Primary*: BOH3*
parameters:
primary_total_stoich: 1.0 # add BOH3 to B, no contribution to TAlk
primary_variable: concentration # provide solver with BOH3_conc (mol m-3)
total_variable: concentration # provide solver with B_conc (mol m-3)
variable_attributes:
Primary_conc%initial_value: 0.427 # This produces initial B_conc ~ 0.427 (as initial BOH4_conc is very small)
Primary_conc%norm_value: 1.0
R_conc%norm_value: 1.0
##########################################################################
# H2O speciation
########################################################################
# Define OH_conc and add to TAlk
# OH- + H+ <--> H2O, [OH] = K_w / [H]
OH_conc: { class: ReactionAqEqb, parameters: { K_power: 1.0, # we need K_eqb in numerator
K_eqb: 6.0e-14, # K_w 6.0e-14 (mol^2 kg-2) equilibrium constant of water at S=35, T=25°C,
K_density_power: 2.0, # K_eqb*rho_ref^2 to convert to (mol^2 m-6)
Reactants: ["OH_conc", "H_conc"], Products: [],
N_components: ["TAlk_calc"],
}
}
##########################################################################
# Boron speciation
##########################################################################
# Define BOH4_conc and add to B, TAlk
# B(OH)_4^- + H+ <--> B(OH)3 + H2O, [B(OH)_4^-] = K_B [B(OH)3] / [H+]
BOH4_conc: { class: ReactionAqEqb, parameters: { K_power: 1.0, # we need K_eqb in numerator
K_eqb: 2.5e-9, # (mol kg-1) K_B equilibrium constant of B(OH)4-
K_density_power: 1.0, # K_eqb*rho_ref to convert to (mol m-3)
Reactants: ["BOH4_conc", "H_conc"], Products: ["BOH3_conc"],
N_components: ["B_calc", "TAlk_calc"],
}
}
##########################################################################
# Carbon speciation
##########################################################################
# Define HCO3_conc and add to DIC, TAlk
# HCO3- + H+ <--> CO2 + H2O, [HCO3-] = K_1 [CO2] / [H+]
HCO3_conc: { class: ReactionAqEqb, parameters: { K_power: 1.0, # we need K_eqb in numerator
K_eqb: 1.4e-6, # (mol kg-1) K_1 equilibrium constant of CO2_aq and HCO3-
K_density_power: 1.0, # K_eqb*rho_ref to convert to (mol m-3)
Reactants: ["HCO3_conc", "H_conc"], Products: ["CO2_aq_conc"],
N_components: ["DIC_calc", "TAlk_calc"],
}
}
# Define CO3_conc and add to DIC, TAlk
# CO3-- + H+ <--> HCO3- + H2O, [CO3--] = K_2 [HCO3-] / [H+]
CO3_conc: { class: ReactionAqEqb, parameters: { K_power: 1.0, # we need K_eqb in numerator
K_eqb: 1.2e-9, # (mol kg-1) K_2 equilibrium constant of HCO3- and CO32-
K_density_power: 1.0, # K_eqb*rho_ref to convert to (mol m-3)
Reactants: ["CO3_conc", "H_conc"], Products: ["HCO3_conc"],
N_components: ["DIC_calc", "2*TAlk_calc"],
}
}
oceanfloor:
# unused here, set up by ReactionOceanNoTransport
oceansurface:
# unused here, set up by ReactionOceanNoTransportRun script
The script to run the model (file examples/configurable_chemistry/run_ex2.jl) contains:
import PALEOboxes as PB
import PALEOmodel
import PALEOocean
import PALEOaqchem
using Plots
#####################################################
# Create model
model = PB.create_model_from_config(joinpath(@__DIR__, "config_ex2.yaml"), "Minimal_Alk_pH")
#########################################################
# 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, DAE")
# first run is slow as it includes JIT time
# NB: Jacobian is required to use implicit variables, so use integrateDAEForwardDiff instead of integrateDAE
@time PALEOmodel.ODE.integrateDAEForwardDiff(
paleorun, initial_state, modeldata, (0.0, 200.0),
solvekwargs=(
reltol=1e-6,
# 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, ["ocean.TAlk_conc", "ocean.DIC_conc"], (cell=1,);
ylabel="TAlk, DIC conc (mol m-3)"))
display(plot(paleorun.output, ["ocean.TAlk_conc","ocean.TAlk_conc_sms"], (cell=1,)))
display(plot(paleorun.output, "ocean.pH", (cell=1,)))
display(plot(paleorun.output, ["ocean.H_conc", "ocean.OH_conc"], (cell=1,);
ylabel="H2O species (mol m-3)", ylim=(0, Inf)))
display(plot(paleorun.output, ["ocean.B_conc", "ocean.BOH4_conc", "ocean.BOH3_conc"], (cell=1,);
ylabel="B species (mol m-3)", ylim=(0, Inf)))
display(plot(paleorun.output, ["ocean.DIC_conc", "ocean.HCO3_conc", "ocean.CO3_conc", "ocean.CO2_aq_conc"], (cell=1,);
ylabel="DIC species (mol m-3)", ylim=(0, Inf)))
display(
plot(
paleorun.output,
["ocean.HCO3_conc", "ocean.CO3_conc", "ocean.CO2_aq_conc", "ocean.BOH4_conc", "ocean.BOH3_conc", "ocean.H_conc", "ocean.OH_conc",],
(cell=1, );
coords=["tmodel"=>("ocean.pH",),], # plot against pH instead of tmodel
ylabel="H2O, B, DIC species (mol m-3)", ylim=(0.5e-3, 0.5e1), yscale=:log10,
legend_background_color=nothing,
legend=:bottom,
)
)
and produces output showing the change, if TAlk_conc increase, how the carbonic acid and pH change in ocean:
┌ Info:
│ ================================================================================
│ create_model_from_config: configmodel Minimal_Alk_pH
│ config_file: /home/runner/work/PALEOtutorials.jl/PALEOtutorials.jl/examples/configurable_chemistry/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.add_Alk classname ReactionFluxPerturb
│ set parameters: [Default] field_data =PALEOboxes.ScalarData
│ set parameters: [config.yaml] perturb_times =[-1.0e30, 1.0e30]
│ set parameters: [config.yaml] perturb_totals =[2.66e16, 2.66e16]
│ set parameters: [Default] perturb_deltas =[0.0, 0.0]
│ set parameters: [Default] extrapolate =throw
│ set parameters: [Default] extrapolate_before =NaN
│ set parameters: [Default] extrapolate_before_delta=NaN
│ set parameters: [Default] extrapolate_after =NaN
└ set parameters: [Default] extrapolate_after_delta=NaN
┌ Info:
│ ================================================================================
│ creating domain 'ocean' ID=2
└ ================================================================================
┌ Info: create_reaction_from_config: ocean.OH_conc classname ReactionAqEqb
│ set parameters: [config.yaml] Reactants =["OH_conc", "H_conc"]
│ set parameters: [config.yaml] Products =Any[]
│ set parameters: [config.yaml] K_eqb =6.0e-14
│ set parameters: [config.yaml] K_density_power =2.0
│ set parameters: [config.yaml] K_power =1.0
└ set parameters: [config.yaml] N_components =["TAlk_calc"]
┌ Info: create_reaction_from_config: ocean.oceantransport classname ReactionOceanNoTransport
│ set parameters: [config.yaml] area =[3.6e14]
└ set parameters: [config.yaml] depth =[3688.0]
┌ Info: create_reaction_from_config: ocean.DIC_total_CO2_aq_primary classname ReactionImplicitReservoir
│ set parameters: [config.yaml] primary_total_stoich=1.0
│ set parameters: [Default] primary_other_components=String[]
│ set parameters: [config.yaml] primary_variable =concentration
│ set parameters: [config.yaml] total_variable =concentration
└ set parameters: [Default] total =false
┌ Info: create_reaction_from_config: ocean.BOH4_conc classname ReactionAqEqb
│ set parameters: [config.yaml] Reactants =["BOH4_conc", "H_conc"]
│ set parameters: [config.yaml] Products =["BOH3_conc"]
│ set parameters: [config.yaml] K_eqb =2.5e-9
│ set parameters: [config.yaml] K_density_power =1.0
│ set parameters: [config.yaml] K_power =1.0
└ set parameters: [config.yaml] N_components =["B_calc", "TAlk_calc"]
┌ Info: create_reaction_from_config: ocean.TA_total_H_primary classname ReactionImplicitReservoir
│ set parameters: [config.yaml] primary_total_stoich=-1.0
│ set parameters: [Default] primary_other_components=String[]
│ set parameters: [config.yaml] primary_variable =p_concentration
│ set parameters: [config.yaml] total_variable =concentration
└ set parameters: [Default] total =false
┌ Info: create_reaction_from_config: ocean.B_total_BOH3_primary classname ReactionImplicitReservoir
│ set parameters: [config.yaml] primary_total_stoich=1.0
│ set parameters: [Default] primary_other_components=String[]
│ set parameters: [config.yaml] primary_variable =concentration
│ set parameters: [config.yaml] total_variable =concentration
└ set parameters: [Default] total =false
┌ Info: create_reaction_from_config: ocean.HCO3_conc classname ReactionAqEqb
│ set parameters: [config.yaml] Reactants =["HCO3_conc", "H_conc"]
│ set parameters: [config.yaml] Products =["CO2_aq_conc"]
│ set parameters: [config.yaml] K_eqb =1.4e-6
│ set parameters: [config.yaml] K_density_power =1.0
│ set parameters: [config.yaml] K_power =1.0
└ set parameters: [config.yaml] N_components =["DIC_calc", "TAlk_calc"]
┌ Info: create_reaction_from_config: ocean.CO3_conc classname ReactionAqEqb
│ set parameters: [config.yaml] Reactants =["CO3_conc", "H_conc"]
│ set parameters: [config.yaml] Products =["HCO3_conc"]
│ set parameters: [config.yaml] K_eqb =1.2e-9
│ set parameters: [config.yaml] K_density_power =1.0
│ set parameters: [config.yaml] K_power =1.0
└ set parameters: [config.yaml] N_components =["DIC_calc", "2*TAlk_calc"]
┌ Info:
│ ================================================================================
│ creating domain 'oceansurface' ID=3
└ ================================================================================
┌ Info:
│ ================================================================================
│ creating domain 'oceanfloor' ID=4
└ ================================================================================
┌ Info:
│ ================================================================================
│ set_model_geometry
└ ================================================================================
[ Info: set ocean.oceansurface Subdomain size=1
[ Info: set ocean.oceanfloor Subdomain size=1
[ Info: Ocean.set_model_domains:
[ Info: ocean Domain size=1 grid=UnstructuredVectorGrid(ncells=1, cellnames=Dict{Symbol, Int64}(), subdomains: ["oceansurface", "oceanfloor"])
[ Info: oceansurface Domain size=1 grid=UnstructuredVectorGrid(ncells=1, cellnames=Dict{Symbol, Int64}(), subdomains: ["ocean"])
[ Info: oceanfloor Domain size=1 grid=UnstructuredVectorGrid(ncells=1, cellnames=Dict{Symbol, Int64}(), subdomains: ["ocean"])
┌ Info:
│ ================================================================================
│ register_reaction_methods!
└ ================================================================================
[ Info: ReactionFluxPerturb.register_methods! global.add_Alk
[ Info: register_methods! ocean.OH_conc
[ Info: register_methods! ocean.BOH4_conc
[ Info: register_methods! ocean.HCO3_conc
[ Info: register_methods! ocean.CO3_conc
┌ Info:
│ ================================================================================
│ link_variables: first pass
└ ================================================================================
┌ Info:
│ ================================================================================
│ link_variables: register_reaction_dynamic_methods and configure variables
└ ================================================================================
┌ Info: _configure_variables: ReactionFluxPerturb global.add_Alk variable_links:
└ set variable_links: F -> ocean.TAlk_sms
┌ Info: _configure_variables: ReactionImplicitReservoir ocean.DIC_total_CO2_aq_primary variable_links:
│ set variable_links: Primary_conc -> CO2_aq_conc
│ set variable_links: R_conc -> DIC_conc
│ set variable_links: R_sms -> DIC_sms
│ set variable_links: R_conc_sms -> DIC_conc_sms
│ set variable_links: R_calc -> DIC_calc
│ set variable_links: R -> DIC
│ _configure_variables: ReactionImplicitReservoir ocean.DIC_total_CO2_aq_primary variable_attributes:
│ set attribute: Primary_conc :initial_value = 1.94
│ set attribute: Primary_conc :norm_value = 1.0
└ set attribute: R_conc :norm_value = 1.0
┌ Info: _configure_variables: ReactionImplicitReservoir ocean.TA_total_H_primary variable_links:
│ set variable_links: R_conc -> TAlk_conc
│ set variable_links: R_sms -> TAlk_sms
│ set variable_links: R_conc_sms -> TAlk_conc_sms
│ set variable_links: R_calc -> TAlk_calc
│ set variable_links: R -> TAlk
│ set variable_links: Primary_conc -> H_conc
│ set variable_links: Primary_pconc -> pH
│ _configure_variables: ReactionImplicitReservoir ocean.TA_total_H_primary variable_attributes:
│ set attribute: Primary_pconc :initial_value = 4.28
│ set attribute: R_conc :norm_value = 1.0
└ set attribute: Primary_pconc :norm_value = 1.0
┌ Info: _configure_variables: ReactionImplicitReservoir ocean.B_total_BOH3_primary variable_links:
│ set variable_links: Primary_conc -> BOH3_conc
│ set variable_links: R_conc -> B_conc
│ set variable_links: R_sms -> B_sms
│ set variable_links: R_conc_sms -> B_conc_sms
│ set variable_links: R_calc -> B_calc
│ set variable_links: R -> B
│ _configure_variables: ReactionImplicitReservoir ocean.B_total_BOH3_primary variable_attributes:
│ set attribute: Primary_conc :initial_value = 0.427
│ set attribute: Primary_conc :norm_value = 1.0
└ set attribute: R_conc :norm_value = 1.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 2 variables (hostdep=nothing)
[ Info: set :field_data=PALEOboxes.ScalarData for host-dependent Variable global.tforce
[ Info: Domain ocean data dimensions PALEOboxes.NamedDimension[] allocating 31 variables (hostdep=nothing)
[ Info: Domain oceansurface data dimensions PALEOboxes.NamedDimension[] allocating 1 variables (hostdep=nothing)
[ Info: Domain oceanfloor 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 0 0 0 0 0 1
│ ocean 0 0 3 0 3 0 0
│ oceansurface 0 0 0 0 0 0 0
│ oceanfloor 0 0 0 0 0 0 0
│ ------------------------------------------------------------------------------------------------------------------------------------------------
│ Total - 0 3 0 3 0 1
│ ------------------------------------------------------------------------------------------------------------------------------------------------
│ n_state_vars 3 (stateexplicit 0 + statetotal 3 + state 0)
│ n_equations 3 (stateexplicit 0 + total 3 + constraint 0)
│ including all host-dependent non-state Variables
└ host-dependent non-state Variables (:vfunction PB.VF_Undefined): ["global.tforce"]
┌ Info:
│ ================================================================================
│ initialize_reactiondata! (modeldata arrays_indices=1:1)
└ ================================================================================
┌ Info:
│ ================================================================================
│ dispatch_setup :setup
└ ================================================================================
┌ Info:
│ ================================================================================
│ dispatch_setup :norm_value
└ ================================================================================
[ Info: ocean.DIC_total_CO2_aq_primary.setup_initialvalue_vars_default:
[ Info: init_values! :norm_value ocean.CO2_aq_conc = 1.0
[ Info: init_values! :norm_value ocean.DIC_conc = 1.0
[ Info: ocean.TA_total_H_primary.setup_initialvalue_vars_default:
[ Info: init_values! :norm_value ocean.pH = 1.0
[ Info: init_values! :norm_value ocean.TAlk_conc = 1.0
[ Info: ocean.B_total_BOH3_primary.setup_initialvalue_vars_default:
[ Info: init_values! :norm_value ocean.BOH3_conc = 1.0
[ Info: init_values! :norm_value ocean.B_conc = 1.0
┌ Info:
│ ================================================================================
│ dispatch_setup :initial_value
└ ================================================================================
[ Info: ocean.DIC_total_CO2_aq_primary.setup_initialvalue_vars_default:
[ Info: init_values! :initial_value ocean.CO2_aq_conc = 1.94
[ Info: init_values! :initial_value ocean.DIC_conc = 0.0
[ Info: ocean.TA_total_H_primary.setup_initialvalue_vars_default:
[ Info: init_values! :initial_value ocean.pH = 4.28
[ Info: init_values! :initial_value ocean.TAlk_conc = 0.0
[ Info: ocean.B_total_BOH3_primary.setup_initialvalue_vars_default:
[ Info: init_values! :initial_value ocean.BOH3_conc = 0.427
[ Info: init_values! :initial_value ocean.B_conc = 0.0
┌ Info:
│ ================================================================================
│ PALEOmodel.initialize! done
└ ================================================================================
initial_state: [0.427, 1.94, 4.28]
initial_deriv: [0.0, 0.0, 0.020034948180284406]
integrate, DAE
┌ Info:
│ ================================================================================
│ PALEOmodel.ODE.integrateDAE:
│ tspan: (0.0, 200.0)
│ algorithm: Sundials.IDA{:KLU, Nothing, Nothing}(0, 0, 0, 5, 7, 0.33, 3, 10, 0.0033, 5, 4, 10, 100, true, false, nothing, nothing)
│ Jacobian: ForwardDiffSparse
│ using 1 BLAS threads
└ ================================================================================
[ Info: DAEfunction: using Jacobian ForwardDiffSparse
[ Info: jac_config_dae: jac_ad=ForwardDiffSparse
[ Info: using ForwardDiff sparse Jacobian with sparsity calculated at t=0.0
[ Info: calcJacobianSparsitySparsityTracing!
[ Info: add_arrays_data! (arrays_eltype=SparsityTracing.ADval{Float64}, arrays_tagname=sparsity_tracing)
[ Info: calcJacobianSparsitySparsityTracing! do_deriv
[ Info: calcJacobianSparsitySparsityTracing! initial_jac size=(3, 3) nnz=0 non-zero=0 at time=0.0
[ Info: fill_sparse_jac: initial_jac nnz 0 (0.0%) non-zero 0
[ Info: fill_sparse_jac: after filling diagonal nnz 3
[ Info: jac_prototype nnz=3 num colors=1
[ Info: calculating dTdS for 3 Total Variables
[ Info: add_arrays_data! (arrays_eltype=SparsityTracing.ADval{Float64}, arrays_tagname=sparsity_tracing_implicit)
[ Info: initial_dTdS size=(3, 3) nnz=7 non-zero=7 at time=0.0
[ Info: fill_sparse_jac: initial_jac nnz 7 (77.8%) non-zero 7
[ Info: implicit_prototype nnz=7 num colors=3
[ Info: combined jac_prototype nnz=9 num colors=3
[ Info: add_arrays_data! (arrays_eltype=ForwardDiff.Dual{Nothing, Float64, 3}, arrays_tagname=jac_ad)
[ Info: calling get_inconsistent_initial_deriv
3.091375 seconds (4.09 M allocations: 278.691 MiB, 11.04% gc time, 99.95% compilation time)
┌ Info:
│ ================================================================================
│ print_sol_stats:
│ retcode=Success
│ alg=Sundials.IDA{:KLU, Nothing, Nothing}(0, 0, 0, 5, 7, 0.33, 3, 10, 0.0033, 5, 4, 10, 100, true, false, nothing, nothing)
│ stats=SciMLBase.DEStats
│ Number of function 1 evaluations: 370
│ Number of function 2 evaluations: 0
│ Number of W matrix evaluations: 48
│ Number of linear solves: 0
│ Number of Jacobians created: 48
│ Number of nonlinear solver iterations: 368
│ Number of nonlinear solver convergence failures: 1
│ 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: 191
│ Number of rejected steps: 5
│ length(sol.t) 197
│ size(sol) (3, 197)
└ ================================================================================
[ Info: ODE.calc_output_sol!: 197 records
[ Info: ODE.calc_output_sol!: done
┌ Info:
│ ================================================================================
│ PALEOmodel.ODE.integrateDAE: done
└ ================================================================================
6.718112 seconds (7.85 M allocations: 539.950 MiB, 6.33% gc time, 99.68% 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 code37×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 add_Alk/FApplied VariableDomPropDep mol yr-1 VF_Undefined ScalarSpace () ScalarData flux perturbation applied, for d…
2 │ global tforce VariableDomPropDep yr VF_Undefined ScalarSpace () ScalarData time at which to apply perturbat…
3 │ ocean Abox VariableDomPropDep m^2 VF_Undefined CellSpace () ScalarData horizontal area of box
4 │ ocean B VariableDomPropDep mol VF_Undefined CellSpace () ScalarData total or component R
5 │ ocean BOH3_conc VariableDomPropDep mol m-3 VF_StateTotal CellSpace () ScalarData concentration of primary species
6 │ ocean BOH4_conc VariableDomPropDep mol m-3 VF_Undefined CellSpace () ScalarData aqueous concentration or activity
7 │ ocean B_calc VariableDomContribTarget mol VF_Undefined CellSpace () ScalarData contributions to total R_calc_co…
8 │ ocean B_conc VariableDomContribTarget mol m-3 VF_Total CellSpace () ScalarData total or component R_conc
9 │ ocean B_conc_sms VariableDomContribTarget mol m-3 yr-1 VF_Deriv CellSpace () ScalarData total or component R_conc source…
10 │ ocean B_sms VariableDomContribTarget mol yr-1 VF_Undefined CellSpace () ScalarData total or component R source - si…
11 │ ocean CO2_aq_conc VariableDomPropDep mol m-3 VF_StateTotal CellSpace () ScalarData concentration of primary species
12 │ ocean CO3_conc VariableDomPropDep mol m-3 VF_Undefined CellSpace () ScalarData aqueous concentration or activity
13 │ ocean DIC VariableDomPropDep mol VF_Undefined CellSpace () ScalarData total or component R
14 │ ocean DIC_calc VariableDomContribTarget mol VF_Undefined CellSpace () ScalarData contributions to total R_calc_co…
15 │ ocean DIC_conc VariableDomContribTarget mol m-3 VF_Total CellSpace () ScalarData total or component R_conc
16 │ ocean DIC_conc_sms VariableDomContribTarget mol m-3 yr-1 VF_Deriv CellSpace () ScalarData total or component R_conc source…
17 │ ocean DIC_sms VariableDomContribTarget mol yr-1 VF_Undefined CellSpace () ScalarData total or component R source - si…
18 │ ocean HCO3_conc VariableDomPropDep mol m-3 VF_Undefined CellSpace () ScalarData aqueous concentration or activity
19 │ ocean H_conc VariableDomPropDep mol m-3 VF_Undefined CellSpace () ScalarData concentration of primary species
20 │ ocean OH_conc VariableDomPropDep mol m-3 VF_Undefined CellSpace () ScalarData aqueous concentration or activity
21 │ ocean TAlk VariableDomPropDep mol VF_Undefined CellSpace () ScalarData total or component R
22 │ ocean TAlk_calc VariableDomContribTarget mol VF_Undefined CellSpace () ScalarData contributions to total R_calc_co…
23 │ ocean TAlk_conc VariableDomContribTarget mol m-3 VF_Total CellSpace () ScalarData total or component R_conc
24 │ ocean TAlk_conc_sms VariableDomContribTarget mol m-3 yr-1 VF_Deriv CellSpace () ScalarData total or component R_conc source…
25 │ ocean TAlk_sms VariableDomContribTarget mol yr-1 VF_Undefined CellSpace () ScalarData total or component R source - si…
26 │ ocean pH VariableDomPropDep VF_StateTotal CellSpace () ScalarData -log10(concentration of primary …
27 │ ocean pressure VariableDomPropDep dbar VF_Undefined CellSpace () ScalarData Ocean pressure
28 │ ocean rho_ref VariableDomPropDep kg m^-3 VF_Undefined CellSpace () ScalarData Ocean transport model density co…
29 │ ocean volume VariableDomPropDep m^3 VF_Undefined CellSpace () ScalarData volume of ocean cells
30 │ ocean volume_total VariableDomPropDep m^3 VF_Undefined ScalarSpace () ScalarData total volume of ocean cells
31 │ ocean zlower VariableDomPropDep m VF_Undefined CellSpace () ScalarData depth of lower surface of box (m)
32 │ ocean zmid VariableDomPropDep m VF_Undefined CellSpace () ScalarData mean depth of box
33 │ ocean zupper VariableDomPropDep m VF_Undefined CellSpace () ScalarData depth of upper surface of box (m…
34 │ oceanfloor Afloor VariableDomPropDep m^2 VF_Undefined CellSpace () ScalarData horizontal area of seafloor at b…
35 │ oceanfloor Afloor_total VariableDomPropDep m^2 VF_Undefined ScalarSpace () ScalarData total area of seafloor
36 │ oceanfloor zfloor VariableDomPropDep m VF_Undefined CellSpace () ScalarData depth of ocean floor (m, -ve)
37 │ oceansurface Asurf VariableDomPropDep m^2 VF_Undefined CellSpace () ScalarData horizontal area of oceansurface