O, A reservoirs can now be const

Add 'const' Parameter to ReactionAtmOcean_O, ReactionAtmOcean_A
PALEOtoolkit · Sep 21, 2022 · e679d25 · e679d25
1 parent 8921849
commit e679d25
Show file tree

Hide file tree

Showing 9 changed files with 130 additions and 66 deletions.
diff --git a/Project.toml b/Project.toml
@@ -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"
 
@@ -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"

diff --git a/examples/COPSE/COPSE_reloaded_reloaded.jl b/examples/COPSE/COPSE_reloaded_reloaded.jl
@@ -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
@@ -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)
@@ -76,8 +83,6 @@ run = PALEOmodel.Run(model=model, output = PALEOmodel.OutputWriters.OutputMemory
 
 # PB.TestUtils.bench_model(run.model)
 
-
-
 ##############################
 # Plot output
 ###############################

diff --git a/examples/COPSE/COPSE_reloaded_reloaded_cfg.yaml b/examples/COPSE/COPSE_reloaded_reloaded_cfg.yaml
@@ -8,6 +8,8 @@ model1:
         SrIsotope: IsotopeLinear
         CIsotopeReacts: true
         tforcevar: global.tforce
+        const_O: false
+        const_A: false
     domains:
         fluxAtoLand:
 
@@ -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*
@@ -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

diff --git a/examples/COPSE/copse_bergman2004_bergman2004_expts.jl b/examples/COPSE/copse_bergman2004_bergman2004_expts.jl
@@ -1,8 +1,5 @@
 
 
-import PALEOboxes as PB
-import PALEOmodel
-
 function copse_bergman2004_bergman2004_expts(
     expts;   
     modelpars=Dict{}(),

diff --git a/examples/COPSE/copse_reloaded_bergman2004_expts.jl b/examples/COPSE/copse_reloaded_bergman2004_expts.jl
@@ -1,6 +1,4 @@
 
-import PALEOboxes as PB
-import PALEOmodel
 
 function copse_reloaded_bergman2004_expts(
     basemodel, expts;

diff --git a/examples/COPSE/copse_reloaded_reloaded_expts.jl b/examples/COPSE/copse_reloaded_reloaded_expts.jl
@@ -1,9 +1,4 @@
 
-import PALEOboxes as PB
-
-import PALEOmodel
-
-
 function copse_reloaded_reloaded_expts(
     basemodel, expts;
     modelpars=Dict{}(),
@@ -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)
 

diff --git a/examples/Project.toml b/examples/Project.toml
@@ -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"
diff --git a/src/COPSE/MapAtmOceanReservoirs.jl b/src/COPSE/MapAtmOceanReservoirs.jl
@@ -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))")
@@ -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,
@@ -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,
@@ -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))")
@@ -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[]

diff --git a/src/PALEOcopse.jl b/src/PALEOcopse.jl
@@ -1,5 +1,10 @@
 module PALEOcopse
 
+import SnoopPrecompile
+
+import Logging
+import PALEOboxes as PB
+
 function moduledir()
     return dirname(@__DIR__)
 end
@@ -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