Ocean chemistry
These examples illustrate how to implement:
- a minimal model of the marine carbonate system
- air-sea exchange of CO2 between ocean and atmosphere Domains
See PALEOmodel documentation for more information on analysing model output.
Example 1 Minimal Alk-pH model
This example implements a new Reaction_Alk_pH to establish 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.
The carbonate chemistry system is solved as an algebraic constraint, calculating TAlk_calcu as a function of pH and then using a PALEO solver to solve for the pH value that makes TAlk_calcu equal to the required value. In PALEO this is implemented by defining pH as a VarState (attribute :vfunction = PB.VF_State) and the algebraic constraint TAlk_error as a VarConstraint (with attribute :vfunction = PB.VF_Constraint). The combined system of TAlk and DIC reservoirs and pH is then a Differential Algebraic Equation (DAE) with two state variables and one algebraic constraint, and is integrated forward in time by a DAE solver PALEOmodel.ODE.integrateDAE.
Additional code files
The Reaction code (Reaction_Alk_pH in file examples/ocean_chemistry/reactions_Alk_pH.jl) now produces calculation of different carbon species HCO3_conc, CO3_conc, CO2_aq_conc, boron species BOH4_conc, water species H_conc and OH_conc given DIC_conc, TAlk_conc and pH. Difference from required alkalinity TAlk_error is then calculated and labelled as an algebraic constraint. Note that there is loop over ocean cells to support arbitrary ocean models:
module Min_Alk_pH
import PALEOboxes as PB
using PALEOboxes.DocStrings # for $(PARS) and $(METHODS_DO)
"""
Reaction_Alk_pH
Minimal example for aqueous carbonate system.
Solves for carbon, boron and water species given `pH`, calculates difference `TAlk_error` from
required alkalinity.
Use in conjunction with a DAE solver, where this Reaction provides an algebraic constraint `TAlk_error` on `pH`.
# Parameters
$(PARS)
# Methods and Variables
$(METHODS_DO)
"""
Base.@kwdef mutable struct Reaction_Alk_pH{P} <: PB.AbstractReaction
base::PB.ReactionBase
pars::P = PB.ParametersTuple(
PB.ParDouble("K_1", 1.4e-6, units="mol kg-1", description="equilibrium constant of CO2_aq and HCO3-"),
PB.ParDouble("K_2", 1.2e-9, units="mol kg-1", description="equilibrium constant of HCO3- and CO32-"),
PB.ParDouble("K_w", 6.0e-14, units="mol^2 kg-2", description="equilibrium constant of water at S=35, T=25°C"),
PB.ParDouble("K_B", 2.5e-9, units="mol kg-1", description="equilibrium constant of B(OH)4-"),
)
end
function PB.register_methods!(rj::Reaction_Alk_pH)
vars = [
PB.VarDep("DIC_conc", "mol m-3", "DIC concentration"),
PB.VarDep("TAlk_conc", "mol m-3", "TA concentration"),
PB.VarDep("B_conc", "mol m-3", "total Boron concentration"),
PB.VarConstraint("TAlk_error", "mol m-3", "in order to solve TA, we set it"),
PB.VarProp("HCO3_conc", "mol m-3", "HCO3- concentration"),
PB.VarProp("CO3_conc", "mol m-3", "CO32- concentration"),
PB.VarProp("CO2_aq_conc", "mol m-3", "CO2_aq concentration"),
PB.VarProp("BOH3_conc", "mol m-3", "BOH3 concentration"),
PB.VarProp("BOH4_conc", "mol m-3", "BOH4- concentration"),
PB.VarProp("H_conc", "mol m-3", "concentration of H+"),
PB.VarProp("OH_conc", "mol m-3", "concentration of OH-"),
PB.VarState("pH", "", "it is the calcalation for pH"),
PB.VarDep("density", "kg m-3", "ocean density"),
]
PB.add_method_do!(rj, do_Min_Alk_pH, (PB.VarList_namedtuple(vars), ) )
PB.add_method_setup!(rj, setup_carbchem, (PB.VarList_namedtuple(vars), ),)
return nothing
end
# called at model start
# provide an initial value for pH (exact value is not important, just needs to be in a reasonable range)
function setup_carbchem( m::PB.ReactionMethod, pars, (vars, ), cellrange::PB.AbstractCellRange, attribute_name )
attribute_name in (:initial_value, :norm_value) || return
for i in cellrange.indices
vars.pH[i] = 8.0
end
return nothing
end
# do method, called each main loop timestep
function do_Min_Alk_pH(m::PB.ReactionMethod, pars, (varsdata, ), cellrange::PB.AbstractCellRange, deltat)
# rj = m.reaction
for i in cellrange.indices
density = varsdata.density[i]
# mol kg-1 mol m-3 kg m-3
DIC_conc_kg = varsdata.DIC_conc[i] / density
B_total_kg = varsdata.B_conc[i] / density
# mol kg-1 = mol L-1 / kg m-3 * L m-3
H_conc_kg = 10^( -varsdata.pH[i] ) / density * 1000.0
# mol kg-1 mol kg-1 mol kg-1 mol kg-1 mol kg-1 mol kg-1 mol kg-1
CO2_aq_conc_kg = DIC_conc_kg / ( 1 + pars.K_1[]/H_conc_kg + (pars.K_1[]*pars.K_2[])/((H_conc_kg)^2) )
# mol kg-1 mol kg-1 mol kg-1 mol kg-1 mol kg-1 mol kg-1
HCO3_conc_kg = DIC_conc_kg / ( 1 + H_conc_kg/pars.K_1[] + pars.K_2[]/H_conc_kg )
# mol kg-1 mol kg-1 mol kg-1 mol kg-1 mol kg-1 mol kg-1 mol kg-1
CO3_conc_kg = DIC_conc_kg / ( 1 + H_conc_kg/pars.K_2[] + ((H_conc_kg)^2)/(pars.K_1[]*pars.K_2[]) )
# mol kg-1 mol kg-1 mol kg-1 mol kg-1
BOH4_conc_kg = B_total_kg / ( 1 + H_conc_kg/pars.K_B[] )
BOH3_conc_kg = B_total_kg - BOH4_conc_kg
# mol kg-1 mol2 kg-2 mol kg-1
OH_conc_kg = pars.K_w[] / H_conc_kg
# mol kg-1 mol kg-1 mol kg-1 mol kg-1 mol kg-1 mol kg-1
TAlk_conc_kg_calcu = HCO3_conc_kg + 2 * CO3_conc_kg + BOH4_conc_kg + OH_conc_kg - H_conc_kg
# mol m-3 kg m-3 mol kg-1
varsdata.H_conc[i] = density * H_conc_kg
varsdata.OH_conc[i] = density * OH_conc_kg
varsdata.HCO3_conc[i] = density * HCO3_conc_kg
varsdata.CO3_conc[i] = density * CO3_conc_kg
varsdata.CO2_aq_conc[i] = density * CO2_aq_conc_kg
varsdata.BOH4_conc[i] = density * BOH4_conc_kg
varsdata.BOH3_conc[i] = density * BOH3_conc_kg
# mol m-3 mol m-3 mol kg-1 kg m-3
varsdata.TAlk_error[i] = varsdata.TAlk_conc[i] - TAlk_conc_kg_calcu * density
end
return nothing
end
end
# moduleDocumentation (generated by the Julia docstring) reads:
Main.Min_Alk_pH.Reaction_Alk_pH — TypeReaction_Alk_pHMinimal example for aqueous carbonate system.
Solves for carbon, boron and water species given pH, calculates difference TAlk_error from required alkalinity.
Use in conjunction with a DAE solver, where this Reaction provides an algebraic constraint TAlk_error on pH.
Parameters
K_1[Float64]=1.4e-6 (mol kg-1),default_value=1.4e-6,description="equilibrium constant of CO2_aq and HCO3-"K_2[Float64]=1.2e-9 (mol kg-1),default_value=1.2e-9,description="equilibrium constant of HCO3- and CO32-"K_w[Float64]=6.0e-14 (mol^2 kg-2),default_value=6.0e-14,description="equilibrium constant of water at S=35, T=25°C"K_B[Float64]=2.5e-9 (mol kg-1),default_value=2.5e-9,description="equilibrium constant of B(OH)4-"
Methods and Variables
do_Min_Alk_pHDIC_conc(mol m-3),VT_ReactDependency,description="DIC concentration"TAlk_conc(mol m-3),VT_ReactDependency,description="TA concentration"B_conc(mol m-3),VT_ReactDependency,description="total Boron concentration"TAlk_error(mol m-3),VT_ReactContributor,VF_Constraint,description="in order to solve TA, we set it"HCO3_conc(mol m-3),VT_ReactProperty,description="HCO3- concentration"CO3_conc(mol m-3),VT_ReactProperty,description="CO32- concentration"CO2_aq_conc(mol m-3),VT_ReactProperty,description="CO2_aq concentration"BOH3_conc(mol m-3),VT_ReactProperty,description="BOH3 concentration"BOH4_conc(mol m-3),VT_ReactProperty,description="BOH4- concentration"H_conc(mol m-3),VT_ReactProperty,description="concentration of H+"OH_conc(mol m-3),VT_ReactProperty,description="concentration of OH-"pH(),VT_ReactDependency,VF_State,description="it is the calcalation for pH"density(kg m-3),VT_ReactDependency,description="ocean density"
yaml configuration file
The model configuration (file examples/ocean_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
solve_Alk_pH:
class: Reaction_Alk_pH
parameters:
K_1: 1.4e-6
K_2: 1.2e-9
K_w: 6.0e-14
K_B: 2.5e-9
variable_links:
density: rho_ref # OceanNoTransport provides density as rho_ref
oceanfloor:
# unused here, set up by ReactionOceanNoTransport
oceansurface:
# unused here, set up by ReactionOceanNoTransportRun script
The script to run the model (file examples/ocean_chemistry/run_ex1.jl) contains:
import PALEOboxes as PB
import PALEOmodel
import PALEOocean
using Plots
include(joinpath(@__DIR__, "reactions_Alk_pH.jl")) # @__DIR__ so still runs when building docs
#####################################################
# 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, tmodel=(1.0, 1e12)); # omit first point (pH 8 starting condition)
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/ocean_chemistry/config_ex1.yaml
└ ================================================================================
[ Info: Model.parameters:
┌ Info:
│ ================================================================================
│ creating Domains
└ ================================================================================
[ Info: generated Reaction catalog with 73 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.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.constant_B classname ReactionReservoirConst
└ set parameters: [Default] field_data =PALEOboxes.ScalarData
┌ 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.solve_Alk_pH classname Reaction_Alk_pH
│ set parameters: [config.yaml] K_1 =1.4e-6
│ set parameters: [config.yaml] K_2 =1.2e-9
│ set parameters: [config.yaml] K_w =6.0e-14
└ set parameters: [config.yaml] K_B =2.5e-9
┌ 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:
│ ================================================================================
│ 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: 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: 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: Reaction_Alk_pH ocean.solve_Alk_pH variable_links:
└ set variable_links: density -> rho_ref
┌ 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 24 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 1 0 1 0
│ oceansurface 0 0 0 0 0 0 0
│ oceanfloor 0 0 0 0 0 0 0
│ ------------------------------------------------------------------------------------------------------------------------------------------------
│ Total - 2 0 1 0 1 1
│ ------------------------------------------------------------------------------------------------------------------------------------------------
│ n_state_vars 3 (stateexplicit 2 + statetotal 0 + state 1)
│ n_equations 3 (stateexplicit 2 + total 0 + constraint 1)
│ 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.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.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, 8.0]
initial_deriv: [2.66e16, 0.0, -2.299331618163945]
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.726053 seconds (617.42 k allocations: 42.582 MiB, 1.86% gc time, 99.91% 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: 405
│ Number of function 2 evaluations: 0
│ Number of W matrix evaluations: 98
│ Number of linear solves: 0
│ Number of Jacobians created: 98
│ Number of nonlinear solver iterations: 403
│ 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: 261
│ Number of rejected steps: 8
│ length(sol.t) 270
│ size(sol) (3, 270)
└ ================================================================================
[ Info: ODE.calc_output_sol!: 270 records
[ Info: ODE.calc_output_sol!: done
┌ Info:
│ ================================================================================
│ PALEOmodel.ODE.integrateDAE: done
└ ================================================================================
0.933085 seconds (769.02 k allocations: 54.140 MiB, 2.94% gc time, 99.08% 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 code30×8 DataFrame
Row │ domain name type units vfunction space field_data description
│ String String DataType String Variable… DataType 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_Undefined CellSpace ScalarData BOH3 concentration
5 │ ocean BOH4_conc VariableDomPropDep mol m-3 VF_Undefined CellSpace ScalarData BOH4- concentration
6 │ ocean B_conc VariableDomPropDep mol m-3 VF_Undefined CellSpace ScalarData concentration
7 │ ocean CO2_aq_conc VariableDomPropDep mol m-3 VF_Undefined CellSpace ScalarData CO2_aq concentration
8 │ ocean CO3_conc VariableDomPropDep mol m-3 VF_Undefined CellSpace ScalarData CO32- concentration
9 │ ocean DIC VariableDomPropDep mol VF_StateExplicit CellSpace ScalarData vector reservoir
10 │ ocean DIC_conc VariableDomPropDep mol m-3 VF_Undefined CellSpace ScalarData concentration
11 │ ocean DIC_sms VariableDomContribTarget mol yr-1 VF_Deriv CellSpace ScalarData vector reservoir source-sinks
12 │ ocean HCO3_conc VariableDomPropDep mol m-3 VF_Undefined CellSpace ScalarData HCO3- concentration
13 │ ocean H_conc VariableDomPropDep mol m-3 VF_Undefined CellSpace ScalarData concentration of H+
14 │ ocean OH_conc VariableDomPropDep mol m-3 VF_Undefined CellSpace ScalarData concentration of OH-
15 │ ocean TAlk VariableDomPropDep mol VF_StateExplicit CellSpace ScalarData vector reservoir
16 │ ocean TAlk_conc VariableDomPropDep mol m-3 VF_Undefined CellSpace ScalarData concentration
17 │ ocean TAlk_error VariableDomContribTarget mol m-3 VF_Constraint CellSpace ScalarData in order to solve TA, we set it
18 │ ocean TAlk_sms VariableDomContribTarget mol yr-1 VF_Deriv CellSpace ScalarData vector reservoir source-sinks
19 │ ocean pH VariableDomPropDep VF_State CellSpace ScalarData it is the calcalation for pH
20 │ ocean pressure VariableDomPropDep dbar VF_Undefined CellSpace ScalarData Ocean pressure
21 │ ocean rho_ref VariableDomPropDep kg m^-3 VF_Undefined CellSpace ScalarData Ocean transport model density co…
22 │ ocean volume VariableDomPropDep m^3 VF_Undefined CellSpace ScalarData volume of ocean cells
23 │ ocean volume_total VariableDomPropDep m^3 VF_Undefined ScalarSpace ScalarData total volume of ocean cells
24 │ ocean zlower VariableDomPropDep m VF_Undefined CellSpace ScalarData depth of lower surface of box (m)
25 │ ocean zmid VariableDomPropDep m VF_Undefined CellSpace ScalarData mean depth of box
26 │ ocean zupper VariableDomPropDep m VF_Undefined CellSpace ScalarData depth of upper surface of box (m…
27 │ oceanfloor Afloor VariableDomPropDep m^2 VF_Undefined CellSpace ScalarData horizontal area of seafloor at b…
28 │ oceanfloor Afloor_total VariableDomPropDep m^2 VF_Undefined ScalarSpace ScalarData total area of seafloor
29 │ oceanfloor zfloor VariableDomPropDep m VF_Undefined CellSpace ScalarData depth of ocean floor (m, -ve)
30 │ oceansurface Asurf VariableDomPropDep m^2 VF_Undefined CellSpace ScalarData horizontal area of oceansurfaceExample 2 Air-sea exchange
This example adds air-sea exchange of CO2 to Example 1 Minimal Alk-pH model.
This configuration adds an atmosphere Domain atm with a state variable for atmospheric CO2. Following the standard PALEO convention for coupling spatial Domains, air-sea exchange is implemented by a combination of the new reaction Reaction_Min_AirSeaExchange in the oceansurface Domain, a ReactionFluxTarget in a fluxAtmtoOceansurface Domain to store the calculated flux, and a pair of ReactionFluxTransfers to apply the calculated fluxes to the atmosphere CO2 and ocean DIC reservoirs. NB: the ocean.oceansurface subdomain represents the subset of ocean cells adjacent to the surface, and contains the same number of cells as the oceansurface and fluxAtmtoOceansurface Domains, with a 1-1 correspondence.
Additional code files
In order to evaluate the CO2 flux change between air and sea, we add a file (file examples/ocean_chemistry/reactions_AirSeaExchange.jl) that implements a Reaction_Min_AirSeaExchange to calculate air-sea CO2 exchange following Henry's Law. NB: the reaction is implemented for a generic gas X that is then linked to the CO2 and DIC variables using the variable_links: section in the .yaml configuration file.
module Min_AirSeaExchange
import PALEOboxes as PB
using PALEOboxes.DocStrings # for $(PARS) and $(METHODS_DO)
"""
Reaction_Min_AirSeaExchange
Minimal example, just make a easy way to illustrate Air-Sea exchange.
This implements exchange between Air and Sea for a generic gas `X` with fixed Henry law coefficient `K_0`
and piston velocity `vpiston`.
The Reaction-local Variables `X_aq_conc`, `pXatm`, `X_airsea_exchange` should be linked to the appropriate variables using the
`variable_links:` section in the .yaml file.
# Parameters
$(PARS)
# Methods and Variables
$(METHODS_DO)
"""
Base.@kwdef mutable struct Reaction_Min_AirSeaExchange{P} <: PB.AbstractReaction
base::PB.ReactionBase
pars::P = PB.ParametersTuple(
PB.ParDouble("K_0", 3.4e-2*1e3, units="mol m-3 atm-1", description="Henry Law coefficient"),
PB.ParDouble("vpiston", 1138.8, units="m yr-1", description="piston value for a whole year, 365 days"),
)
end
function PB.register_methods!(rj::Reaction_Min_AirSeaExchange)
vars = [
PB.VarDep( "X_aq_conc", "mol m-3", "ocean concentration per cell"),
PB.VarDepScalar( "pXatm", "atm", "atmospheric partial pressure, unit is atm"),
PB.VarContrib( "X_airsea_exchange", "mol yr-1", "calculated airsea exchange flux for gas X"),
PB.VarDep( "area", "m2", "surface area"),
]
PB.add_method_do!(rj, do_Min_AirSeaExchange, (PB.VarList_namedtuple(vars), ) )
return nothing
end
# do method, called each main loop timestep
function do_Min_AirSeaExchange(m::PB.ReactionMethod, pars, (varsdata, ), cellrange::PB.AbstractCellRange, deltat)
# mol m-3 mol m-3 atm-1 atm
X_aq_conc_eqb = pars.K_0[] * varsdata.pXatm[]
for i in cellrange.indices
# mol m-3 mol m-3
X_aq_conc = varsdata.X_aq_conc[i]
# mol yr-1 mol m-3 mol m-3 m yr-1 m2
varsdata.X_airsea_exchange[i] -= (X_aq_conc - X_aq_conc_eqb) * pars.vpiston[] * varsdata.area[i]
end
return nothing
end
end # moduleDocumentation (generated by the Julia docstring) reads:
Main.Min_AirSeaExchange.Reaction_Min_AirSeaExchange — TypeReaction_Min_AirSeaExchangeMinimal example, just make a easy way to illustrate Air-Sea exchange.
This implements exchange between Air and Sea for a generic gas X with fixed Henry law coefficient K_0 and piston velocity vpiston.
The Reaction-local Variables X_aq_conc, pXatm, X_airsea_exchange should be linked to the appropriate variables using the variable_links: section in the .yaml file.
Parameters
K_0[Float64]=34.0 (mol m-3 atm-1),default_value=34.0,description="Henry Law coefficient"vpiston[Float64]=1138.8 (m yr-1),default_value=1138.8,description="piston value for a whole year, 365 days"
Methods and Variables
do_Min_AirSeaExchangeX_aq_conc(mol m-3),VT_ReactDependency,description="ocean concentration per cell"pXatm(atm),VT_ReactDependency,description="atmospheric partial pressure, unit is atm"X_airsea_exchange(mol yr-1),VT_ReactContributor,description="calculated airsea exchange flux for gas X"area(m2),VT_ReactDependency,description="surface area"
yaml configuration file
The model configuration (file examples/ocean_chemistry/config_ex2.yaml) contains three Reservoirs DIC, TAlk and CO2. Following reservoirs Example 4 Transfer between Domains, we use ReactionFluxTarget and ReactionFluxTransfer to transfer CO2_airsea_exchange between DIC reservoir and CO2 reservoir:
Minimal_Alk_pH_AirSea:
domains:
global:
reactions:
sum_C:
class: ReactionSum
parameters:
vars_to_add: ["ocean.DIC_total", "atm.CO2"]
variable_links:
sum: C_total
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: ReactionReservoirTotal
variable_links:
R*: TAlk*
variable_attributes:
R:initial_value: 0.0
R:norm_value: 2.4e0
reservoir_DIC:
class: ReactionReservoirTotal
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
solve_Alk_pH:
class: Reaction_Alk_pH
parameters:
K_1: 1.4e-6
K_2: 1.2e-9
K_w: 6.0e-14
K_B: 2.5e-9
variable_links:
density: rho_ref # OceanNoTransport provides density as rho_ref
# --------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
oceansurface:
# air-sea exchange reactions
reactions:
# calculate air-sea exchange flux
solve_AirSea_Exchange:
class: Reaction_Min_AirSeaExchange
parameters:
K_0: 3.4e1 # (mol m-3 atm-1) = 3.4e-2 (mol L-1 atm-1) * 1e3 (L m-3)
vpiston: 1138.8 # m yr-1
variable_links:
pXatm: atm.pCO2atm
X_aq_conc: ocean.oceansurface.CO2_aq_conc
area: Asurf # ocean surface area
X_airsea_exchange: fluxAtmtoOceansurface.flux_CO2
# apply air-sea flux to ocean surface
transfer_fluxAtmtoOceansurface:
class: ReactionFluxTransfer
parameters:
transfer_matrix: Identity
input_fluxes: fluxAtmtoOceansurface.flux_$fluxname$
output_fluxes: ocean.oceansurface.$fluxname$_sms
variable_links:
output_CO2: ocean.oceansurface.DIC_sms
# --------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
oceanfloor:
# --------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
atm:
reactions:
reservoir_CO2_atm:
class: ReactionReservoirAtm
parameters:
moles1atm: 1.77e20
variable_links:
R*: CO2*
pRatm: pCO2atm
pRnorm: pCO2PAL
variable_attributes:
R:norm_value: 4.956e16 # 280e-6 ppm
R:initial_value: 4.956e16
transfer_AtmtoOceansurface:
class: ReactionFluxTransfer
parameters:
transfer_matrix: Distribute
transfer_multiplier: -1.0
input_fluxes: fluxAtmtoOceansurface.flux_$fluxname$
output_fluxes: $fluxname$_sms
variable_links:
output_CO2: atm.CO2_sms
# --------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
fluxAtmtoOceansurface:
reactions:
fluxtarget:
class: ReactionFluxTarget
parameters:
flux_totals: true
fluxlist: ["CO2"]
Run script
The script to run the model (file examples/ocean_chemistry/run_ex2.jl) contains:
import PALEOboxes as PB
import PALEOmodel
# import PALEOcopse
import PALEOocean
using Plots
include(joinpath(@__DIR__,"reactions_Alk_pH.jl")) # @__DIR__ so still runs when building docs
include(joinpath(@__DIR__,"reactions_AirSeaExchange.jl")) # @__DIR__ so still runs when building docs
include(joinpath(@__DIR__,"reactions_ReservoirAtm.jl")) # @__DIR__ so still runs when building docs
#####################################################
# Create model
model = PB.create_model_from_config(joinpath(@__DIR__, "config_ex2.yaml"), "Minimal_Alk_pH_AirSea")
#########################################################
# 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.DIC_conc", "ocean.HCO3_conc", "ocean.CO3_conc", "ocean.CO2_aq_conc"], (cell=1,);
ylabel="DIC species (mol m-3)", legend=:topleft))
display(plot(paleorun.output, "atm.CO2", ))
display(plot(paleorun.output, "atm.pCO2atm", ))
display(plot(paleorun.output, "atm.CO2_sms", ))
display(plot(paleorun.output, "fluxAtmtoOceansurface.flux_CO2", (cell=1,);
ylabel="air->sea flux (mol yr-1)"))
display(plot(paleorun.output, "global.C_total";
ylabel="total carbon (mol)"))
display(plot(paleorun.output, ["global.C_total", "atm.CO2", "ocean.DIC_total"];
ylabel="atm ocean carbon (mol)", legend=:topleft))
and produces output showing the change, if TAlk_conc increase, how the carbonic acid and pH change in ocean and CO2 change in the air:
display(plot(paleorun.output, ["global.C_total", "atm.CO2", "ocean.DIC_total"] ; ylabel="atm-ocean carbon (mol)"))┌ Info:
│ ================================================================================
│ create_model_from_config: configmodel Minimal_Alk_pH_AirSea
│ config_file: /home/runner/work/PALEOtutorials.jl/PALEOtutorials.jl/examples/ocean_chemistry/config_ex2.yaml
└ ================================================================================
[ Info: Model.parameters:
┌ Info:
│ ================================================================================
│ creating Domains
└ ================================================================================
[ Info: generated Reaction catalog with 73 Reactions
┌ Info:
│ ================================================================================
│ creating domain 'global' ID=1
└ ================================================================================
┌ Info: create_reaction_from_config: global.sum_C classname ReactionSum
│ set parameters: [config.yaml] vars_to_add =["ocean.DIC_total", "atm.CO2"]
│ set parameters: [Default] vars_prefix =
│ set parameters: [Default] component_to_add =0
└ set parameters: [Default] vectorsum =false
┌ 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.reservoir_DIC classname ReactionReservoirTotal
│ set parameters: [Default] field_data =PALEOboxes.ScalarData
│ set parameters: [Default] total =true
│ set parameters: [Default] limit_delta_conc =0.0
└ set parameters: [Default] state_conc =false
┌ Info: create_reaction_from_config: ocean.constant_B classname ReactionReservoirConst
└ set parameters: [Default] field_data =PALEOboxes.ScalarData
┌ Info: create_reaction_from_config: ocean.reservoir_TA classname ReactionReservoirTotal
│ set parameters: [Default] field_data =PALEOboxes.ScalarData
│ set parameters: [Default] total =true
│ 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.solve_Alk_pH classname Reaction_Alk_pH
│ set parameters: [config.yaml] K_1 =1.4e-6
│ set parameters: [config.yaml] K_2 =1.2e-9
│ set parameters: [config.yaml] K_w =6.0e-14
└ set parameters: [config.yaml] K_B =2.5e-9
┌ Info:
│ ================================================================================
│ creating domain 'oceansurface' ID=3
└ ================================================================================
┌ Info: create_reaction_from_config: oceansurface.transfer_fluxAtmtoOceansurface classname ReactionFluxTransfer
│ set parameters: [config.yaml] input_fluxes =fluxAtmtoOceansurface.flux_$fluxname$
│ set parameters: [config.yaml] output_fluxes =ocean.oceansurface.$fluxname$_sms
│ set parameters: [Default] transfer_multiplier =1.0
└ set parameters: [config.yaml] transfer_matrix =Identity
┌ Info: create_reaction_from_config: oceansurface.solve_AirSea_Exchange classname Reaction_Min_AirSeaExchange
│ set parameters: [config.yaml] K_0 =34.0
└ set parameters: [config.yaml] vpiston =1138.8
┌ Info:
│ ================================================================================
│ creating domain 'oceanfloor' ID=4
└ ================================================================================
┌ Info:
│ ================================================================================
│ creating domain 'fluxAtmtoOceansurface' ID=5
└ ================================================================================
┌ Info: create_reaction_from_config: fluxAtmtoOceansurface.fluxtarget classname ReactionFluxTarget
│ set parameters: [config.yaml] fluxlist =["CO2"]
│ set parameters: [Default] target_prefix =flux_
│ set parameters: [config.yaml] flux_totals =true
└ set parameters: [Default] const_stub =false
┌ Info:
│ ================================================================================
│ creating domain 'atm' ID=6
└ ================================================================================
┌ Info: create_reaction_from_config: atm.transfer_AtmtoOceansurface classname ReactionFluxTransfer
│ set parameters: [config.yaml] input_fluxes =fluxAtmtoOceansurface.flux_$fluxname$
│ set parameters: [config.yaml] output_fluxes =$fluxname$_sms
│ set parameters: [config.yaml] transfer_multiplier =-1.0
└ set parameters: [config.yaml] transfer_matrix =Distribute
┌ Info: create_reaction_from_config: atm.reservoir_CO2_atm classname ReactionReservoirAtm
│ set parameters: [Default] field_data =PALEOboxes.ScalarData
│ set parameters: [Default] const =false
└ set parameters: [config.yaml] moles1atm =1.77e20
┌ 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: fluxAtmtoOceansurface Domain size=1 grid=UnstructuredVectorGrid(ncells=1, cellnames=Dict{Symbol, Int64}(), subdomains: ["ocean"])
[ 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: register_methods: global.sum_C PALEOboxes.VariableStats.ReactionSum
│ add 1.0 * ocean.DIC_total
└ add 1.0 * atm.CO2
[ Info: ReactionFluxPerturb.register_methods! global.add_Alk
[ Info: register_methods! ReactionReservoirAtm atm.reservoir_CO2_atm field_data=PALEOboxes.ScalarData
┌ Info:
│ ================================================================================
│ link_variables: first pass
└ ================================================================================
┌ Info:
│ ================================================================================
│ link_variables: register_reaction_dynamic_methods and configure variables
└ ================================================================================
[ Info: Reaction global.sum_C Variable global.C_total adding data=PALEOboxes.ScalarData
┌ Info: _configure_variables: ReactionSum global.sum_C variable_links:
└ set variable_links: sum -> C_total
┌ 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_total -> DIC_total
│ 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: 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: ReactionReservoir ocean.reservoir_TA variable_links:
│ set variable_links: R -> TAlk
│ set variable_links: R_sms -> TAlk_sms
│ set variable_links: R_total -> TAlk_total
│ 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: Reaction_Alk_pH ocean.solve_Alk_pH variable_links:
└ set variable_links: density -> rho_ref
┌ Info: _configure_variables: ReactionFluxTransfer oceansurface.transfer_fluxAtmtoOceansurface variable_links:
└ set variable_links: output_CO2 -> ocean.oceansurface.DIC_sms
┌ Info: _configure_variables: Reaction_Min_AirSeaExchange oceansurface.solve_AirSea_Exchange variable_links:
│ set variable_links: X_airsea_exchange -> fluxAtmtoOceansurface.flux_CO2
│ set variable_links: X_aq_conc -> ocean.oceansurface.CO2_aq_conc
│ set variable_links: area -> Asurf
└ set variable_links: pXatm -> atm.pCO2atm
┌ Info: _configure_variables: ReactionFluxTransfer atm.transfer_AtmtoOceansurface variable_links:
└ set variable_links: output_CO2 -> atm.CO2_sms
┌ Info: _configure_variables: ReactionReservoirAtm atm.reservoir_CO2_atm variable_links:
│ set variable_links: R -> CO2
│ set variable_links: R_sms -> CO2_sms
│ set variable_links: pRatm -> pCO2atm
│ set variable_links: pRnorm -> pCO2PAL
│ _configure_variables: ReactionReservoirAtm atm.reservoir_CO2_atm variable_attributes:
│ set attribute: R :norm_value = 4.956e16
└ set attribute: R :initial_value = 4.956e16
┌ 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 :field_data=PALEOboxes.ScalarData for host-dependent Variable global.tforce
[ Info: Domain ocean data dimensions PALEOboxes.NamedDimension[] allocating 26 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: Domain fluxAtmtoOceansurface data dimensions PALEOboxes.NamedDimension[] allocating 2 variables (hostdep=nothing)
[ Info: Domain atm 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 1
│ ocean 0 2 0 1 0 1 0
│ oceansurface 0 0 0 0 0 0 0
│ oceanfloor 0 0 0 0 0 0 0
│ fluxAtmtoOceansurface 0 0 0 0 0 0 0
│ atm 0 1 0 0 0 0 0
│ ------------------------------------------------------------------------------------------------------------------------------------------------
│ Total - 3 0 1 0 1 1
│ ------------------------------------------------------------------------------------------------------------------------------------------------
│ n_state_vars 4 (stateexplicit 3 + statetotal 0 + state 1)
│ n_equations 4 (stateexplicit 3 + total 0 + constraint 1)
│ 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: prepare_do_transfer: ReactionMethod oceansurface.transfer_fluxAtmtoOceansurface.do_transfer
│ active fluxes:
│ CO2 fluxAtmtoOceansurface.flux_CO2 -> ocean.DIC_sms
└ using Identity transfer matrix * 1.0
┌ Info: prepare_do_transfer: ReactionMethod atm.transfer_AtmtoOceansurface.do_transfer
│ active fluxes:
│ CO2 fluxAtmtoOceansurface.flux_CO2 -> atm.CO2_sms
└ using Distribute transfer matrix size (input, output) (1, 1) * -1.0
┌ 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: atm.reservoir_CO2_atm.setup_initialvalue_vars_default:
[ Info: init_values! :norm_value atm.CO2 = 4.956e16
[ 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: atm.reservoir_CO2_atm.setup_initialvalue_vars_default:
[ Info: init_values! :initial_value atm.CO2 = 4.956e16
[ 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, 4.956e16, 8.0]
initial_deriv: [2.66e16, -1.1428818049576428e15, 1.1428818049576428e15, -2.299331618163945]
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.724461 seconds (619.75 k allocations: 42.694 MiB, 1.92% gc time, 99.87% 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: 595
│ Number of function 2 evaluations: 0
│ Number of W matrix evaluations: 85
│ Number of linear solves: 0
│ Number of Jacobians created: 85
│ Number of nonlinear solver iterations: 593
│ 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: 408
│ Number of rejected steps: 5
│ length(sol.t) 414
│ size(sol) (4, 414)
└ ================================================================================
[ Info: ODE.calc_output_sol!: 414 records
[ Info: ODE.calc_output_sol!: done
┌ Info:
│ ================================================================================
│ PALEOmodel.ODE.integrateDAE: done
└ ================================================================================
0.909246 seconds (778.54 k allocations: 54.642 MiB, 1.53% gc time, 98.61% compilation time)
GKS: Possible loss of precision in routine SET_WINDOW
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 code39×8 DataFrame
Row │ domain name type units vfunction space field_data description
│ String String DataType String Variable… DataType DataType String
─────┼───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
1 │ atm CO2 VariableDomPropDep mol VF_StateExplicit ScalarSpace ScalarData scalar reservoir
2 │ atm CO2_sms VariableDomContribTarget mol yr-1 VF_Deriv ScalarSpace ScalarData scalar reservoir source-sinks
3 │ atm pCO2PAL VariableDomPropDep VF_Undefined ScalarSpace ScalarData scalar atmosphere reservoir norm…
4 │ atm pCO2atm VariableDomPropDep atm VF_Undefined ScalarSpace ScalarData scalar atmosphere reservoir part…
5 │ fluxAtmtoOceansurface flux_CO2 VariableDomContribTarget mol yr-1 VF_Undefined CellSpace ScalarData flux input
6 │ fluxAtmtoOceansurface flux_total_CO2 VariableDomPropDep mol yr-1 VF_Undefined ScalarSpace ScalarData total flux input
7 │ global C_total VariableDomPropDep unknown VF_Undefined ScalarSpace ScalarData sum of specified variables
8 │ global add_Alk/FApplied VariableDomPropDep mol yr-1 VF_Undefined ScalarSpace ScalarData flux perturbation applied, for d…
9 │ global tforce VariableDomPropDep yr VF_Undefined ScalarSpace ScalarData time at which to apply perturbat…
10 │ ocean Abox VariableDomPropDep m^2 VF_Undefined CellSpace ScalarData horizontal area of box
11 │ ocean BOH3_conc VariableDomPropDep mol m-3 VF_Undefined CellSpace ScalarData BOH3 concentration
12 │ ocean BOH4_conc VariableDomPropDep mol m-3 VF_Undefined CellSpace ScalarData BOH4- concentration
13 │ ocean B_conc VariableDomPropDep mol m-3 VF_Undefined CellSpace ScalarData concentration
14 │ ocean CO2_aq_conc VariableDomPropDep mol m-3 VF_Undefined CellSpace ScalarData CO2_aq concentration
15 │ ocean CO3_conc VariableDomPropDep mol m-3 VF_Undefined CellSpace ScalarData CO32- concentration
16 │ ocean DIC VariableDomPropDep mol VF_StateExplicit CellSpace ScalarData vector reservoir
17 │ ocean DIC_conc VariableDomPropDep mol m-3 VF_Undefined CellSpace ScalarData concentration
18 │ ocean DIC_sms VariableDomContribTarget mol yr-1 VF_Deriv CellSpace ScalarData vector reservoir source-sinks
19 │ ocean DIC_total VariableDomPropDep mol VF_Undefined ScalarSpace ScalarData total vector reservoir
20 │ ocean HCO3_conc VariableDomPropDep mol m-3 VF_Undefined CellSpace ScalarData HCO3- concentration
21 │ ocean H_conc VariableDomPropDep mol m-3 VF_Undefined CellSpace ScalarData concentration of H+
22 │ ocean OH_conc VariableDomPropDep mol m-3 VF_Undefined CellSpace ScalarData concentration of OH-
23 │ ocean TAlk VariableDomPropDep mol VF_StateExplicit CellSpace ScalarData vector reservoir
24 │ ocean TAlk_conc VariableDomPropDep mol m-3 VF_Undefined CellSpace ScalarData concentration
25 │ ocean TAlk_error VariableDomContribTarget mol m-3 VF_Constraint CellSpace ScalarData in order to solve TA, we set it
26 │ ocean TAlk_sms VariableDomContribTarget mol yr-1 VF_Deriv CellSpace ScalarData vector reservoir source-sinks
27 │ ocean TAlk_total VariableDomPropDep mol VF_Undefined ScalarSpace ScalarData total vector reservoir
28 │ ocean pH VariableDomPropDep VF_State CellSpace ScalarData it is the calcalation for pH
29 │ ocean pressure VariableDomPropDep dbar VF_Undefined CellSpace ScalarData Ocean pressure
30 │ ocean rho_ref VariableDomPropDep kg m^-3 VF_Undefined CellSpace ScalarData Ocean transport model density co…
31 │ ocean volume VariableDomPropDep m^3 VF_Undefined CellSpace ScalarData volume of ocean cells
32 │ ocean volume_total VariableDomPropDep m^3 VF_Undefined ScalarSpace ScalarData total volume of ocean cells
33 │ ocean zlower VariableDomPropDep m VF_Undefined CellSpace ScalarData depth of lower surface of box (m)
34 │ ocean zmid VariableDomPropDep m VF_Undefined CellSpace ScalarData mean depth of box
35 │ ocean zupper VariableDomPropDep m VF_Undefined CellSpace ScalarData depth of upper surface of box (m…
36 │ oceanfloor Afloor VariableDomPropDep m^2 VF_Undefined CellSpace ScalarData horizontal area of seafloor at b…
37 │ oceanfloor Afloor_total VariableDomPropDep m^2 VF_Undefined ScalarSpace ScalarData total area of seafloor
38 │ oceanfloor zfloor VariableDomPropDep m VF_Undefined CellSpace ScalarData depth of ocean floor (m, -ve)
39 │ oceansurface Asurf VariableDomPropDep m^2 VF_Undefined CellSpace ScalarData horizontal area of oceansurfaceExample 3 Minimal modern earth ocean-air
This example shows how ocean(TAlk-DIC)-air(CO2) exchange in modern state Example 3 Minimal modern earth ocean-air.
This configuration .yaml file is the same as Example 2 Air-sea exchange. What is different is that we set TAlk, DIC and K_0 to be modern state value (from (Sarmiento and Gruber, 2006)) using PB.set_parameter_value! and PB.set_variable_attribute!
NB: this is not an accurate model for the modern system ! See PALEOocean.jl for more complete models, including representations of ocean spatial structure and circulation, and more detailed parameterisations of carbonate chemistry (ReactionCO2SYS) and air-sea exchange (ReactionAirSeaCO2)
Run script
The script to run the model (file examples/ocean_chemistry/run_ex3.jl) contains:
import PALEOboxes as PB
import PALEOmodel
# import PALEOcopse
import PALEOocean
using Plots
include(joinpath(@__DIR__,"reactions_Alk_pH.jl")) # @__DIR__ so still runs when building docs
include(joinpath(@__DIR__,"reactions_AirSeaExchange.jl")) # @__DIR__ so still runs when building docs
include(joinpath(@__DIR__,"reactions_ReservoirAtm.jl")) # @__DIR__ so still runs when building docs
#####################################################
# Create model
model = PB.create_model_from_config(joinpath(@__DIR__, "config_ex2.yaml"), "Minimal_Alk_pH_AirSea")
# set to ~modern (1990s) values from Sarmiento 2006 for global mean ocean surface
# K0 (NB: Table A.3 gives global mean ocean surface temperature 17.88 C)
PB.set_parameter_value!(model, "oceansurface", "solve_AirSea_Exchange", "K_0", 38.3) # mol m-3 atm-1, Table 3.2.3 at S=35. T=15℃
# PB.set_parameter_value!(model, "oceansurface", "solve_AirSea_Exchange", "K_0", 33.11) # mol m-3 atm-1, Table 3.2.3 at S=35. T=20℃
# PB.set_parameter_value!(model, "oceansurface", "solve_AirSea_Exchange", "K_0", 29.6) # mol m-3 atm-1, from Sarmiento 2006, at S=35. T=25℃
PB.set_parameter_value!(model, "global", "add_Alk", "perturb_totals", [0.0, 0.0]) # don't add any Alk
PB.set_variable_attribute!(model, "ocean.TAlk", :initial_value, 2.308) # mol m-3 ~ modern value global mean ocean surface (Table A.3)
PB.set_variable_attribute!(model, "ocean.DIC", :initial_value, 2.026) # mol m-3 ~ modern value global mean ocean surface (Table A.3)
#########################################################
# 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, 100.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.DIC_conc", "ocean.HCO3_conc", "ocean.CO3_conc", "ocean.CO2_aq_conc"], (cell=1,);
ylabel="DIC species (mol m-3)", legend=:topleft))
display(plot(paleorun.output, "atm.CO2", ))
display(plot(paleorun.output, "atm.pCO2atm", ))
display(plot(paleorun.output, "atm.CO2_sms", ))
display(plot(paleorun.output, "fluxAtmtoOceansurface.flux_CO2", (cell=1,);
ylabel="air->sea flux (mol yr-1)"))
display(plot(paleorun.output, ["global.C_total", "atm.CO2", "ocean.DIC_total"];
ylabel="atm ocean carbon (mol)", legend=:topleft))and produces output showing an approximate modern steady state:
display(plot(paleorun.output, ["global.C_total", "atm.CO2", "ocean.DIC_total"] ; ylabel="atm-ocean carbon (mol)"))┌ Info:
│ ================================================================================
│ create_model_from_config: configmodel Minimal_Alk_pH_AirSea
│ config_file: /home/runner/work/PALEOtutorials.jl/PALEOtutorials.jl/examples/ocean_chemistry/config_ex2.yaml
└ ================================================================================
[ Info: Model.parameters:
┌ Info:
│ ================================================================================
│ creating Domains
└ ================================================================================
[ Info: generated Reaction catalog with 73 Reactions
┌ Info:
│ ================================================================================
│ creating domain 'global' ID=1
└ ================================================================================
┌ Info: create_reaction_from_config: global.sum_C classname ReactionSum
│ set parameters: [config.yaml] vars_to_add =["ocean.DIC_total", "atm.CO2"]
│ set parameters: [Default] vars_prefix =
│ set parameters: [Default] component_to_add =0
└ set parameters: [Default] vectorsum =false
┌ 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.reservoir_DIC classname ReactionReservoirTotal
│ set parameters: [Default] field_data =PALEOboxes.ScalarData
│ set parameters: [Default] total =true
│ set parameters: [Default] limit_delta_conc =0.0
└ set parameters: [Default] state_conc =false
┌ Info: create_reaction_from_config: ocean.constant_B classname ReactionReservoirConst
└ set parameters: [Default] field_data =PALEOboxes.ScalarData
┌ Info: create_reaction_from_config: ocean.reservoir_TA classname ReactionReservoirTotal
│ set parameters: [Default] field_data =PALEOboxes.ScalarData
│ set parameters: [Default] total =true
│ 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.solve_Alk_pH classname Reaction_Alk_pH
│ set parameters: [config.yaml] K_1 =1.4e-6
│ set parameters: [config.yaml] K_2 =1.2e-9
│ set parameters: [config.yaml] K_w =6.0e-14
└ set parameters: [config.yaml] K_B =2.5e-9
┌ Info:
│ ================================================================================
│ creating domain 'oceansurface' ID=3
└ ================================================================================
┌ Info: create_reaction_from_config: oceansurface.transfer_fluxAtmtoOceansurface classname ReactionFluxTransfer
│ set parameters: [config.yaml] input_fluxes =fluxAtmtoOceansurface.flux_$fluxname$
│ set parameters: [config.yaml] output_fluxes =ocean.oceansurface.$fluxname$_sms
│ set parameters: [Default] transfer_multiplier =1.0
└ set parameters: [config.yaml] transfer_matrix =Identity
┌ Info: create_reaction_from_config: oceansurface.solve_AirSea_Exchange classname Reaction_Min_AirSeaExchange
│ set parameters: [config.yaml] K_0 =34.0
└ set parameters: [config.yaml] vpiston =1138.8
┌ Info:
│ ================================================================================
│ creating domain 'oceanfloor' ID=4
└ ================================================================================
┌ Info:
│ ================================================================================
│ creating domain 'fluxAtmtoOceansurface' ID=5
└ ================================================================================
┌ Info: create_reaction_from_config: fluxAtmtoOceansurface.fluxtarget classname ReactionFluxTarget
│ set parameters: [config.yaml] fluxlist =["CO2"]
│ set parameters: [Default] target_prefix =flux_
│ set parameters: [config.yaml] flux_totals =true
└ set parameters: [Default] const_stub =false
┌ Info:
│ ================================================================================
│ creating domain 'atm' ID=6
└ ================================================================================
┌ Info: create_reaction_from_config: atm.transfer_AtmtoOceansurface classname ReactionFluxTransfer
│ set parameters: [config.yaml] input_fluxes =fluxAtmtoOceansurface.flux_$fluxname$
│ set parameters: [config.yaml] output_fluxes =$fluxname$_sms
│ set parameters: [config.yaml] transfer_multiplier =-1.0
└ set parameters: [config.yaml] transfer_matrix =Distribute
┌ Info: create_reaction_from_config: atm.reservoir_CO2_atm classname ReactionReservoirAtm
│ set parameters: [Default] field_data =PALEOboxes.ScalarData
│ set parameters: [Default] const =false
└ set parameters: [config.yaml] moles1atm =1.77e20
┌ 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: fluxAtmtoOceansurface Domain size=1 grid=UnstructuredVectorGrid(ncells=1, cellnames=Dict{Symbol, Int64}(), subdomains: ["ocean"])
[ 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: register_methods: global.sum_C PALEOboxes.VariableStats.ReactionSum
│ add 1.0 * ocean.DIC_total
└ add 1.0 * atm.CO2
[ Info: ReactionFluxPerturb.register_methods! global.add_Alk
[ Info: register_methods! ReactionReservoirAtm atm.reservoir_CO2_atm field_data=PALEOboxes.ScalarData
┌ Info:
│ ================================================================================
│ link_variables: first pass
└ ================================================================================
┌ Info:
│ ================================================================================
│ link_variables: register_reaction_dynamic_methods and configure variables
└ ================================================================================
[ Info: Reaction global.sum_C Variable global.C_total adding data=PALEOboxes.ScalarData
┌ Info: _configure_variables: ReactionSum global.sum_C variable_links:
└ set variable_links: sum -> C_total
┌ 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_total -> DIC_total
│ 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: 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: ReactionReservoir ocean.reservoir_TA variable_links:
│ set variable_links: R -> TAlk
│ set variable_links: R_sms -> TAlk_sms
│ set variable_links: R_total -> TAlk_total
│ 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: Reaction_Alk_pH ocean.solve_Alk_pH variable_links:
└ set variable_links: density -> rho_ref
┌ Info: _configure_variables: ReactionFluxTransfer oceansurface.transfer_fluxAtmtoOceansurface variable_links:
└ set variable_links: output_CO2 -> ocean.oceansurface.DIC_sms
┌ Info: _configure_variables: Reaction_Min_AirSeaExchange oceansurface.solve_AirSea_Exchange variable_links:
│ set variable_links: X_airsea_exchange -> fluxAtmtoOceansurface.flux_CO2
│ set variable_links: X_aq_conc -> ocean.oceansurface.CO2_aq_conc
│ set variable_links: area -> Asurf
└ set variable_links: pXatm -> atm.pCO2atm
┌ Info: _configure_variables: ReactionFluxTransfer atm.transfer_AtmtoOceansurface variable_links:
└ set variable_links: output_CO2 -> atm.CO2_sms
┌ Info: _configure_variables: ReactionReservoirAtm atm.reservoir_CO2_atm variable_links:
│ set variable_links: R -> CO2
│ set variable_links: R_sms -> CO2_sms
│ set variable_links: pRatm -> pCO2atm
│ set variable_links: pRnorm -> pCO2PAL
│ _configure_variables: ReactionReservoirAtm atm.reservoir_CO2_atm variable_attributes:
│ set attribute: R :norm_value = 4.956e16
└ set attribute: R :initial_value = 4.956e16
┌ 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 :field_data=PALEOboxes.ScalarData for host-dependent Variable global.tforce
[ Info: Domain ocean data dimensions PALEOboxes.NamedDimension[] allocating 26 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: Domain fluxAtmtoOceansurface data dimensions PALEOboxes.NamedDimension[] allocating 2 variables (hostdep=nothing)
[ Info: Domain atm 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 1
│ ocean 0 2 0 1 0 1 0
│ oceansurface 0 0 0 0 0 0 0
│ oceanfloor 0 0 0 0 0 0 0
│ fluxAtmtoOceansurface 0 0 0 0 0 0 0
│ atm 0 1 0 0 0 0 0
│ ------------------------------------------------------------------------------------------------------------------------------------------------
│ Total - 3 0 1 0 1 1
│ ------------------------------------------------------------------------------------------------------------------------------------------------
│ n_state_vars 4 (stateexplicit 3 + statetotal 0 + state 1)
│ n_equations 4 (stateexplicit 3 + total 0 + constraint 1)
│ 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: prepare_do_transfer: ReactionMethod oceansurface.transfer_fluxAtmtoOceansurface.do_transfer
│ active fluxes:
│ CO2 fluxAtmtoOceansurface.flux_CO2 -> ocean.DIC_sms
└ using Identity transfer matrix * 1.0
┌ Info: prepare_do_transfer: ReactionMethod atm.transfer_AtmtoOceansurface.do_transfer
│ active fluxes:
│ CO2 fluxAtmtoOceansurface.flux_CO2 -> atm.CO2_sms
└ using Distribute transfer matrix size (input, output) (1, 1) * -1.0
┌ 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: atm.reservoir_CO2_atm.setup_initialvalue_vars_default:
[ Info: init_values! :norm_value atm.CO2 = 4.956e16
[ 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: atm.reservoir_CO2_atm.setup_initialvalue_vars_default:
[ Info: init_values! :initial_value atm.CO2 = 4.956e16
[ Info: ocean.reservoir_DIC.setup_initialvalue_vars_default:
[ Info: init_values! :initial_value ocean.DIC = 2.026 * volume
[ Info: ocean.reservoir_TA.setup_initialvalue_vars_default:
[ Info: init_values! :initial_value ocean.TAlk = 2.308 * volume
┌ Info:
│ ================================================================================
│ PALEOmodel.initialize! done
└ ================================================================================
initial_state: [3.06428544e18, 2.6898796799999995e18, 4.956e16, 8.0]
initial_deriv: [0.0, -7.148754361020918e14, 7.148754361020918e14, -0.020006738816293]
integrate, DAE
┌ Info:
│ ================================================================================
│ PALEOmodel.ODE.integrateDAE:
│ tspan: (0.0, 100.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.000229 seconds (1.58 k allocations: 60.672 KiB)
┌ 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: 100
│ Number of function 2 evaluations: 0
│ Number of W matrix evaluations: 18
│ Number of linear solves: 0
│ Number of Jacobians created: 18
│ Number of nonlinear solver iterations: 98
│ 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: 70
│ Number of rejected steps: 2
│ length(sol.t) 73
│ size(sol) (4, 73)
└ ================================================================================
[ Info: ODE.calc_output_sol!: 73 records
[ Info: ODE.calc_output_sol!: done
┌ Info:
│ ================================================================================
│ PALEOmodel.ODE.integrateDAE: done
└ ================================================================================
0.003228 seconds (7.73 k allocations: 1.421 MiB)
For more information and cooperation, please communicate with us!