From e679d25292c7487cc44009248426dad055852c1e Mon Sep 17 00:00:00 2001 From: Stuart Daines Date: Wed, 21 Sep 2022 20:23:01 +0100 Subject: [PATCH] O, A reservoirs can now be const Add 'const' Parameter to ReactionAtmOcean_O, ReactionAtmOcean_A --- Project.toml | 5 +- examples/COPSE/COPSE_reloaded_reloaded.jl | 15 ++- .../COPSE/COPSE_reloaded_reloaded_cfg.yaml | 6 +- .../copse_bergman2004_bergman2004_expts.jl | 3 - .../COPSE/copse_reloaded_bergman2004_expts.jl | 2 - .../COPSE/copse_reloaded_reloaded_expts.jl | 14 +-- examples/Project.toml | 1 - src/COPSE/MapAtmOceanReservoirs.jl | 91 +++++++++++++------ src/PALEOcopse.jl | 59 ++++++++---- 9 files changed, 130 insertions(+), 66 deletions(-) diff --git a/Project.toml b/Project.toml index dfc45bc..06d3566 100644 --- 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 index c8dc029..41fdbe2 100644 --- 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 index 88bbaf1..ffb7205 100644 --- 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 index 83fef20..45a67cb 100644 --- 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 index 3875f0c..aacc0c6 100644 --- 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 index c425926..00ce2ed 100644 --- 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, , , , , , , , , 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 index fad9e66..03251f5 100644 --- 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