Skip to content

Commit

Permalink
Merge pull request #18 from PALEOtoolkit/const_O_A
Browse files Browse the repository at this point in the history 
O, A reservoirs can now be const
  • Loading branch information
@sjdaines
sjdaines authored Sep 21, 2022
2 parents 8921849 + e679d25 commit 093272f
Show file tree
Hide file tree
Showing 9 changed files with 130 additions and 66 deletions.
5 changes: 4 additions & 1 deletion Project.toml
View file
Original file line number Diff line number Diff line change
@@ -1,15 +1,17 @@
name = "PALEOcopse"
uuid = "4a6ed817-0e58-48c6-8452-9e9afc8cb508"
authors = ["sd336 "]
version = "0.4.3"
version = "0.4.4"

[deps]
DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
Infiltrator = "5903a43b-9cc3-4c30-8d17-598619ec4e9b"
Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
MAT = "23992714-dd62-5051-b70f-ba57cb901cac"
PALEOboxes = "804b410e-d900-4b2a-9ecd-f5a06d4c1fd4"
Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
SnoopPrecompile = "66db9d55-30c0-4569-8b51-7e840670fc0c"
TestEnv = "1e6cf692-eddd-4d53-88a5-2d735e33781b"
XLSX = "fdbf4ff8-1666-58a4-91e7-1b58723a45e0"

Expand All @@ -23,6 +25,7 @@ PALEOboxes = "0.20.4, 0.21"
PALEOmodel = "0.15.8"
Plots = "1.0"
Roots = "1.0, 2.0"
SnoopPrecompile = "1.0"
TestEnv = "1.0"
XLSX = "0.7, 0.8"
julia = "1.6"
Expand Down
15 changes: 10 additions & 5 deletions examples/COPSE/COPSE_reloaded_reloaded.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,6 @@ global_logger(ConsoleLogger(stderr,Logging.Info))
# load archived model output
# comparemodel=nothing


# Baseline Phanerozoic configuration
# comparemodel = CompareOutput.copse_output_load("reloaded","reloaded_baseline")
comparemodel = nothing
Expand All @@ -32,15 +31,23 @@ model = copse_reloaded_reloaded_expts(
)
tspan=(-1000e6, 0)

# Fix oxygen
# comparemodel = nothing
# model = copse_reloaded_reloaded_expts(
# "reloaded",
# ["baseline", ("set_initial_value", "atmocean", "O", 3.7e19)],
# modelpars=Dict("const_O"=>true),
# )
# tspan=(-1000e6, 0)

# OOE oscillations with VCI and linear mocb
# comparemodel = CompareOutput.copse_output_load("reloaded","reloaded_baseline")
# comparemodel = nothing
# model = copse_reloaded_reloaded_expts(
# "reloaded",
# [
# "VCI",
# "mocbProdLinear",
# ("setpar", "global", "force_LIPs", "co2releasefield", "CO2max"),
# ("set_par", "global", "force_LIPs", "co2releasefield", "CO2max"),
# ],
# )
# tspan=(-1000e6, 1e6)
Expand Down Expand Up @@ -76,8 +83,6 @@ run = PALEOmodel.Run(model=model, output = PALEOmodel.OutputWriters.OutputMemory

# PB.TestUtils.bench_model(run.model)



##############################
# Plot output
###############################
Expand Down
6 changes: 5 additions & 1 deletion examples/COPSE/COPSE_reloaded_reloaded_cfg.yaml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,8 @@ model1:
SrIsotope: IsotopeLinear
CIsotopeReacts: true
tforcevar: global.tforce
const_O: false
const_A: false
domains:
fluxAtoLand:

Expand Down Expand Up @@ -269,6 +271,7 @@ model1:
A:initial_delta: 0.0 # per mil
parameters:
f_atfrac: quadratic
const: external%const_A
variable_links:
TEMP: global.TEMP
pCO2*: atm.pCO2*
Expand All @@ -279,7 +282,8 @@ model1:

reservoir_O:
class: ReactionAtmOcean_O

parameters:
const: external%const_O
variable_attributes:
O:norm_value: 3.7e19
O:initial_value: 3.7e19
Expand Down
3 changes: 0 additions & 3 deletions examples/COPSE/copse_bergman2004_bergman2004_expts.jl
View file
Original file line number Diff line number Diff line change
@@ -1,8 +1,5 @@


import PALEOboxes as PB
import PALEOmodel

function copse_bergman2004_bergman2004_expts(
expts;
modelpars=Dict{}(),
Expand Down
2 changes: 0 additions & 2 deletions examples/COPSE/copse_reloaded_bergman2004_expts.jl
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,4 @@

import PALEOboxes as PB
import PALEOmodel

function copse_reloaded_bergman2004_expts(
basemodel, expts;
Expand Down
14 changes: 7 additions & 7 deletions examples/COPSE/copse_reloaded_reloaded_expts.jl
View file
Original file line number Diff line number Diff line change
@@ -1,9 +1,4 @@

import PALEOboxes as PB

import PALEOmodel


function copse_reloaded_reloaded_expts(
basemodel, expts;
modelpars=Dict{}(),
Expand Down Expand Up @@ -66,11 +61,16 @@ function copse_reloaded_reloaded_expts(
PB.set_parameter_value!(model, "ocean", "marinebiota_copse", "f_ncycle", false)
PB.set_parameter_value!(model, "ocean", "oceanburial_copse", "f_ncycle", false)

elseif length(expt) == 5 && expt[1] == "setpar"
# generic parameter set (setpar, <domain>, <reaction>, <parname>, <parvalue)
elseif length(expt) == 5 && expt[1] == "set_par"
# generic parameter set (set_par, <domain>, <reaction>, <parname>, <parvalue)
_, domname, reactname, parname, parvalue = expt
PB.set_parameter_value!(model, domname, reactname, parname, parvalue)

elseif length(expt) == 4 && expt[1] == "set_initial_value"
# generic :initial_value set (set_initial_value, <domain>, <varname>, <initial_value)
_, domname, varname, initial_value = expt
PB.set_variable_attribute!(model, domname, varname, :initial_value, initial_value)

elseif length(expt)==4 && expt[1] == "CO2pulse"
_, size, pstart, duration = expt # (_, mol C, yr, yr)

Expand Down
1 change: 0 additions & 1 deletion examples/Project.toml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,5 +14,4 @@ Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228"

[compat]
PALEOboxes = "0.20.4"
PALEOmodel = "0.15.8"
91 changes: 64 additions & 27 deletions src/COPSE/MapAtmOceanReservoirs.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,22 +22,19 @@ $(PARS)
# Methods and Variables
$(METHODS_DO)
"""
Base.@kwdef mutable struct ReactionAtmOcean_O <: PB.AbstractReaction
Base.@kwdef mutable struct ReactionAtmOcean_O{P} <: PB.AbstractReaction
base::PB.ReactionBase

pars::P = PB.ParametersTuple(
PB.ParBool("const", false,
description="true to provide constant value, ignoring fluxes from _sms Variable"),
)

norm_value::Float64 = NaN
end

function PB.register_methods!(rj::ReactionAtmOcean_O)

vars = [
PB.VarStateExplicitScalar( "O", "mol", "atm-ocean oxygen"),
PB.VarDerivScalar( "O_sms", "mol yr-1", "atm-ocean oxygen source-sinks"),
PB.VarPropScalar( "O_norm", "", "atm-ocean normalized"),
PB.VarPropScalar( "pO2atm", "atm", "atmospheric pO2"),
PB.VarPropScalar( "pO2PAL", "", "atmospheric pO2 normalized to present day"),
]

# callback function to store Variable norm during setup
function setup_callback(m, attribute_value, v, vdata)
v.localname == "O" || error("setup_callback unexpected Variable $(PB.fullname(v))")
Expand All @@ -47,7 +44,31 @@ function PB.register_methods!(rj::ReactionAtmOcean_O)
return nothing
end

PB.add_method_setup_initialvalue_vars_default!(rj, vars, setup_callback=setup_callback) # initialise state variables
if rj.pars.const[]
O = PB.VarPropScalar( "O", "mol", "atm-ocean oxygen")
O_sms = PB.VarTargetScalar( "O_sms", "mol yr-1", "atm-ocean oxygen source-sinks")
PB.add_method_setup_initialvalue_vars_default!(
rj, [O],
filterfn = v->true, # force setup even though O is not a state Variable
force_initial_norm_value=true, # force initialize as a state variable including :norm_value, even though O is not a state Variable
setup_callback=setup_callback
)
else
O = PB.VarStateExplicitScalar( "O", "mol", "atm-ocean oxygen")
O_sms = PB.VarDerivScalar( "O_sms", "mol yr-1", "atm-ocean oxygen source-sinks")
PB.add_method_setup_initialvalue_vars_default!(rj, [O], setup_callback=setup_callback) # initialise state variables
end
PB.setfrozen!(rj.pars.const) # prevent modification

# sms variable not used by us, but must appear in a method to be linked and created
PB.add_method_do_nothing!(rj, [O_sms])

vars = [
O,
PB.VarPropScalar( "O_norm", "", "atm-ocean normalized"),
PB.VarPropScalar( "pO2atm", "atm", "atmospheric pO2"),
PB.VarPropScalar( "pO2PAL", "", "atmospheric pO2 normalized to present day"),
]

PB.add_method_do!(
rj,
Expand Down Expand Up @@ -117,6 +138,8 @@ Base.@kwdef mutable struct ReactionAtmOcean_A{P} <: PB.AbstractReaction
description="calculate d13CO2, d13DIC relative to A"),
PB.ParBool("fix_cisotopefrac_T", false,
description="remove temperature dependence of d13CO2, d13DIC relative to A"),
PB.ParBool("const", false,
description="true to provide constant value, ignoring fluxes from _sms Variable"),
PB.ParType(PB.AbstractData, "CIsotope", PB.ScalarData,
external=true,
allowed_values=PB.IsotopeTypes,
Expand All @@ -130,17 +153,6 @@ function PB.register_methods!(rj::ReactionAtmOcean_A)

CIsotopeType = rj.pars.CIsotope[]

vars = [
PB.VarStateExplicitScalar("A", "mol", "atm-ocean inorganic carbon (CO2 + DIC)",
attributes=(:field_data=>CIsotopeType,)),
PB.VarDerivScalar( "A_sms", "mol yr-1", "atm-ocean inorganic carbon (CO2 + DIC) source-sinks",
attributes=(:field_data=>CIsotopeType,)),
PB.VarPropScalar( "A_norm", "", "atm-ocean inorganic carbon (CO2 + DIC) normalized to present day"),
PB.VarPropScalar( "pCO2atm","atm", "atmospheric pCO2"),
PB.VarPropScalar( "pCO2PAL","", "atmospheric pCO2 normalized to present day"),
PB.VarPropScalar( "phi", "", "atmospheric pCO2 fraction"),
]

# callback function to store Variable norm during setup
function setup_callback(m, attribute_value, v, vdata)
v.localname == "A" || error("setup_callback unexpected Variable $(PB.fullname(v))")
Expand All @@ -150,21 +162,46 @@ function PB.register_methods!(rj::ReactionAtmOcean_A)
return nothing
end

PB.add_method_setup_initialvalue_vars_default!(rj, vars, setup_callback=setup_callback) # initialise state variables
if rj.pars.const[]
A = PB.VarPropScalar( "A", "mol", "atm-ocean inorganic carbon (CO2 + DIC)",
attributes=(:field_data=>CIsotopeType,))
A_sms = PB.VarTargetScalar( "A_sms", "mol yr-1", "atm-ocean inorganic carbon (CO2 + DIC) source-sinks",
attributes=(:field_data=>CIsotopeType,))
PB.add_method_setup_initialvalue_vars_default!(
rj, [A],
filterfn = v->true, # force setup even though A is not a state Variable
force_initial_norm_value=true, # force initialize as a state variable including :norm_value, even though A is not a state Variable
setup_callback=setup_callback
)
else
A = PB.VarStateExplicitScalar( "A", "mol", "atm-ocean inorganic carbon (CO2 + DIC)",
attributes=(:field_data=>CIsotopeType,))
A_sms = PB.VarDerivScalar( "A_sms", "mol yr-1", "atm-ocean inorganic carbon (CO2 + DIC) source-sinks",
attributes=(:field_data=>CIsotopeType,))
PB.add_method_setup_initialvalue_vars_default!(rj, [A], setup_callback=setup_callback) # initialise state variables
end
PB.setfrozen!(rj.pars.const) # prevent modification

# sms variable not used by us, but must appear in a method to be linked and created
PB.add_method_do_nothing!(rj, [A_sms])

vars = [
A,
PB.VarPropScalar( "A_norm", "", "atm-ocean inorganic carbon (CO2 + DIC) normalized to present day"),
PB.VarPropScalar( "pCO2atm","atm", "atmospheric pCO2"),
PB.VarPropScalar( "pCO2PAL","", "atmospheric pCO2 normalized to present day"),
PB.VarPropScalar( "phi", "", "atmospheric pCO2 fraction"),
]

norm_value = NaN # will be updated in prepare_
PB.add_method_do!(
rj,
do_AtmOcean_A,
(PB.VarList_namedtuple(vars),),
)


if CIsotopeType <: PB.AbstractIsotopeScalar
vars_isotope = [
PB.VarDepScalar("A", "mol", "atm-ocean inorganic carbon (CO2 + DIC)",
attributes=(:field_data=>CIsotopeType,)),

PB.VarDep(A),
PB.VarPropScalar("A_delta", "per mil", "atm-ocean inorganic carbon (CO2 + DIC) delta13C"),
]
if rj.pars.delta_atm_ocean[]
Expand Down
59 changes: 40 additions & 19 deletions src/PALEOcopse.jl
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,10 @@
module PALEOcopse

import SnoopPrecompile

import Logging
import PALEOboxes as PB

function moduledir()
return dirname(@__DIR__)
end
Expand All @@ -15,24 +20,40 @@ include("sedcrust/SedCrust.jl")
include("biogeochem/BioGeoChem.jl")


# Negligible difference Julia 1.7.3
# function precompile_reactions()
# rdict = PB.find_all_reactions()

# reactionlist = ["ReactionSrSed", "ReactionLandWeatheringFluxes", "ReactionSrLand", "ReactionLandArea", "ReactionForce_LIPs",
# "ReactionForce_CK_Solar", "ReactionLandWeatheringRates", "ReactionModelBergman2004",
# "ReactionCIsotopes", "ReactionOceanCOPSE", "ReactionAtmOcean_A", "ReactionSedCrustCOPSE", "ReactionGlobalTemperatureBerner",
# "ReactionSrOceanfloor", "ReactionForce_CPlandrelbergman2004", "ReactionForce_UDWEbergman2004", "ReactionSeafloorWeathering",
# "ReactionSrMantleCrust", "ReactionLandBergman2004", "ReactionAtmOcean_O", "ReactionForce_spreadsheet",
# "ReactionGlobalTemperatureCK1992", "ReactionLandBiota", "ReactionForce_Bbergman2004",
# ]
# for r in reactionlist
# PB.precompile_reaction(rdict, r)
# end

# return nothing
# end

# precompile_reactions()
# TODO gains here with Julia 1.8.1 are marginal (~20s)
if false # VERSION >= v"1.8.0" # negligible benefit from precompile prior to Julia 1.8.0
@SnoopPrecompile.precompile_setup begin
# create Reactions and register methods to precompile this code

# Putting some things in `setup` can reduce the size of the
# precompile file and potentially make loading faster.


rdict = PB.find_all_reactions()
reactionlist = ["ReactionSrSed", "ReactionLandWeatheringFluxes", "ReactionSrLand", "ReactionLandArea", "ReactionForce_LIPs",
"ReactionForce_CK_Solar", "ReactionLandWeatheringRates", "ReactionModelBergman2004",
"ReactionCIsotopes", "ReactionMarineBiotaCOPSE", "ReactionOceanBurialCOPSE", "ReactionAtmOcean_A", "ReactionSedCrustCOPSE", "ReactionGlobalTemperatureBerner",
"ReactionSrOceanfloor", "ReactionForce_CPlandrelbergman2004", "ReactionForce_UDWEbergman2004", "ReactionSeafloorWeathering",
"ReactionSrMantleCrust", "ReactionLandBergman2004", "ReactionAtmOcean_O", "ReactionForce_spreadsheet",
"ReactionGlobalTemperatureCK1992", "ReactionLandBiota", "ReactionForce_Bbergman2004",
]

@SnoopPrecompile.precompile_all_calls begin
# all calls in this block will be precompiled, regardless of whether
# they belong to your package or not (on Julia 1.8 and higher)

Logging.with_logger(Logging.NullLogger()) do
for r in reactionlist
PB.precompile_reaction(rdict, r)
end

#
PB.run_model(joinpath(@__DIR__, "../examples/COPSE/COPSE_bergman2004_bergman2004_cfg.yaml"), "Bergman2004")
PB.run_model(joinpath(@__DIR__, "../examples/COPSE/COPSE_reloaded_bergman2004_cfg.yaml"), "model1")
PB.run_model(joinpath(@__DIR__, "../examples/COPSE/COPSE_reloaded_reloaded_cfg.yaml"), "model1")
end
end
end
end

end

2 comments on commit 093272f

@sjdaines
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator register

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Registration pull request created: JuliaRegistries/General/68725

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via:

git tag -a v0.4.4 -m "<description of version>" 093272f20816713e9c4611e916fd8a0726cb2604
git push origin v0.4.4

Please sign in to comment.