Sync with PALEOboxes v0.19, PALEOmodel v0.15 setup changes

There is now a new call to PB.dispatch_setup(model, :setup, modeldata) with attribute_name=:setup, hence corresponding callback to setup methods. This commit updates COPSE Reactions to do non-state-variable related initialisation in setup callbacks when attribute_name=:setup (replaces previous (ab)use of :initial_value for this).
PALEOtoolkit · Jun 25, 2022 · 904aa57 · 904aa57
1 parent 13334a2
commit 904aa57
Show file tree

Hide file tree

Showing 8 changed files with 150 additions and 71 deletions.
diff --git a/Project.toml b/Project.toml
@@ -1,7 +1,7 @@
 name = "PALEOcopse"
 uuid = "4a6ed817-0e58-48c6-8452-9e9afc8cb508"
 authors = ["sd336 "]
-version = "0.2.3"
+version = "0.3.0"
 
 [deps]
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
@@ -19,8 +19,8 @@ Documenter = "0.27"
 Infiltrator = "1.0"
 Interpolations = "0.13"
 MAT = "0.10"
-PALEOboxes = "0.18.4"
-PALEOmodel = "0.14"
+PALEOboxes = "0.19.1"
+PALEOmodel = "0.15"
 Plots = "1.0"
 Roots = "1.0, 2.0"
 TestEnv = "1.0"

diff --git a/src/Forcings/COPSEForcings.jl b/src/Forcings/COPSEForcings.jl
@@ -229,101 +229,156 @@ function do_force_CPlandrelbergman2004(m::PB.ReactionMethod, (vars, ), cellrange
     return nothing
 end
 
+const LINEARINTERPOLATION_TEMPLATE = Interpolations.LinearInterpolation(
+    [0.0, 1.0], 
+    [NaN, NaN],
+    extrapolation_bc = Interpolations.Flat()
+)
 
 """
     ReactionForce_spreadsheet
 
-Generic forcing interpolated from values in spreadsheet. Spreadsheet should contain a table of numeric data in Sheet1,
-with a single header row, and time in Ma as the first column.
+Generic forcing interpolated from values in spreadsheet.
 
-Interpolates and optionally applies naive extrapolation of forcings into (constant) Precambrian and future
+Spreadsheet should contain a table of numeric data in sheet `sheetname`,
+with a single header row, and time in `timecolumn` and data in `datacolumn`.
+
+Time from `timecolumn` is multiplied by `timemultiplier` to convert to PALEO model time.
+
+Linearly interpolates in time `tforce` to set Variable `forcename`, and extrapolates into past and future, either to values
+in `extrap_value_past`, `extrap_value_future`, or to the closest spreadsheet value if these Parameters are `NaN`.
+
+# Implementation
+Spreadsheet data is read using the Julia XLSX and DataFrames packages with
+
+    xf = XLSX.readxlsx(forcingfile)
+    df = XLSX.eachtablerow(xf[sheetname]) |> DataFrames.DataFrame
 
 # Parameters
 $(PARS)
 
-# Methods and Variables TODO
-- `tforce` (yr) historical time at which to apply forcing
-- `<forcename>` name from Parameter `forcename`
+# Methods and Variables for default Parameters
+$(METHODS_DO)
 """
 Base.@kwdef mutable struct ReactionForce_spreadsheet{P} <:  PB.AbstractReaction
     base::PB.ReactionBase
 
     pars::P = PB.ParametersTuple(
         PB.ParString("forcename", "FORCE",
             description="name of forcing (this will be a variable name in global Domain)"),
+        PB.ParString("datafolder", PALEOcopse.Forcings.srcdir(),
+            description="folder for spreadsheet with forcing data"),
         PB.ParString("datafile", "", 
-            description="spreadsheet with forcing data (path relative to $(PALEOcopse.Forcings.srcdir()))"),
+            description="spreadsheet with forcing data (path relative to 'datafolder')"),
+        PB.ParString("sheetname", "Sheet1",
+            description="sheet name in spreadsheet"),
         PB.ParInt("timecolumn", 1,
             description="column with time data"),
+        PB.ParDouble("timemultiplier", -1.0e6,
+            description="factor to multiply 'timecolumn' by to convert to yr after present day (so times in past are -ve)"),
         PB.ParInt("datacolumn", 2,
             description="column with forcing data"),
-        PB.ParDouble("extrap_value", 1.0, units="", 
-            description="extrapolate value for out-of-range tforce"),
+        PB.ParDouble("extrap_value_past", 1.0, units="", 
+            description="extrapolate value for tforce before earliest value in spreadsheet (NaN to use earliest value)"),
+        PB.ParDouble("extrap_value_future", 1.0, units="", 
+            description="extrapolate value for tforce after latest value in spreadsheet (NaN to use latest value"),
     )
 
     forcing_data    = DataFrames.DataFrame()  # raw data from spreadsheet
 
-    interp_FORCE = nothing
+    force_times::Vector{Float64} = Float64[] # times used for interpolation
+    force_values::Vector{Float64} = Float64[] # values used for interpolation
+    interp_FORCE::typeof(LINEARINTERPOLATION_TEMPLATE) = LINEARINTERPOLATION_TEMPLATE
 
 end
 
-function read_xlsx(forcingfile; sheetname="Sheet1")
-
-    @info "read_xlsx: spreadsheet $(forcingfile) sheet $(sheetname)"
-
-    xf = XLSX.readxlsx(forcingfile)
-    # read data assuming first row is column headers
-    # data, column_labels = XLSX.gettable(xf[sheetname])
-    # convert to float64 and create DataFrame
-    # data_float64 = [Float64.(d) for d in data]
-    # df = DataFrames.DataFrame(data_float64, column_labels)
-
-    # Read using built-in iterator - assumes single header row, automatic type conversion
-    df = XLSX.eachtablerow(xf[sheetname]) |> DataFrames.DataFrame
-
-    @info "read_xlsx read $(DataFrames.describe(df))"
-
-    return df
-end
-
 function PB.register_methods!(rj::ReactionForce_spreadsheet)
 
-    forcingfile = joinpath(PALEOcopse.Forcings.srcdir(), rj.pars.datafile.v)
-    @info "ReactionForce_spreadsheet: $(PB.fullname(rj)) loading $(rj.pars.forcename.v) forcing from datafile $(forcingfile)"
-    rj.forcing_data = read_xlsx(forcingfile)
-    @info "ReactionForce_spreadsheet:  $(PB.fullname(rj)) $(rj.pars.forcename.v) from column $(rj.pars.datacolumn.v) ($(names(rj.forcing_data)[rj.pars.datacolumn.v]))"
-    @info "ReactionForce_spreadsheet:  $(PB.fullname(rj)) extrapolating out-of-range tforce to $(rj.pars.forcename.v) = $(rj.pars.extrap_value.v)"
-
-    # create interpolation object
-    rj.interp_FORCE = Interpolations.LinearInterpolation(
-        -1.0e6*rj.forcing_data[:, rj.pars.timecolumn.v],  Float64.(rj.forcing_data[:, rj.pars.datacolumn.v]), 
-        extrapolation_bc = rj.pars.extrap_value.v )  # fill out of range values
-
     var_tforce = PB.VarDepScalar("tforce",     "yr",  
         "historical time at which to apply forcings, present = 0 yr")
     var_FORCE = PB.VarPropScalar(rj.pars.forcename.v,    "",  
-        "forcing interpolated from $(forcingfile) column $(rj.pars.datacolumn.v)")
+        "forcing interpolated from spreadsheet")
+    PB.setfrozen!(rj.pars.forcename)
 
-    PB.setfrozen!.(PB.get_parameters(rj)) # no modification after spreadsheet read
+    PB.add_method_setup!(
+        rj, 
+        setup_force_spreadsheet,
+        (),
+    )
 
     PB.add_method_do!(
         rj, 
         do_force_spreadsheet,
         (PB.VarList_single(var_tforce), PB.VarList_single(var_FORCE)),
-        p=rj.interp_FORCE
+    )
+
+    return nothing
+end
+
+function setup_force_spreadsheet(m::PB.ReactionMethod, (), cellrange::PB.AbstractCellRange, attribute_name)
+    rj = m.reaction
+
+    attribute_name == :setup || return nothing
+
+    forcingfile = joinpath(rj.pars.datafolder.v, rj.pars.datafile.v)
+    @info "setup_force_spreadsheet ReactionForce_spreadsheet $(PB.fullname(rj)): "*
+        "loading $(rj.pars.forcename.v) forcing from 'datafolder/datafile'='$(forcingfile)'"
+    rj.forcing_data = _read_xlsx(forcingfile, sheetname=rj.pars.sheetname.v)
+
+    sp_times = rj.pars.timemultiplier.v*rj.forcing_data[:, rj.pars.timecolumn.v]
+    @info "    'tforce' from $(rj.pars.timemultiplier.v) * column $(rj.pars.timecolumn.v) ($(names(rj.forcing_data)[rj.pars.timecolumn.v]))"
+
+    sp_values = Float64.(rj.forcing_data[:, rj.pars.datacolumn.v])
+    @info "    '$(rj.pars.forcename.v)' from column $(rj.pars.datacolumn.v) ($(names(rj.forcing_data)[rj.pars.datacolumn.v]))"
+
+    # sort in ascending time order
+    sp_perm = sortperm(sp_times)
+    rj.force_times = sp_times[sp_perm]
+    rj.force_values = sp_values[sp_perm]
+
+    extrap_past = isnan(rj.pars.extrap_value_past.v) ? "earlist value in spreadsheet = $(first(rj.force_values))" : "'extrap_value_past' = $(rj.pars.extrap_value_past.v)"
+    @info "    extrapolating out-of-range tforce < $(first(rj.force_times)) (yr) to $extrap_past"
+    extrap_future = isnan(rj.pars.extrap_value_future.v) ? "latest value in spreadsheet = $(last(rj.force_values))" : "'extrap_value_future' = $(rj.pars.extrap_value_future.v)"
+    @info "    extrapolating out-of-range tforce > $(last(rj.force_times)) (yr) to $extrap_future"
+
+    # create interpolation object
+    rj.interp_FORCE = Interpolations.LinearInterpolation(
+        rj.force_times, rj.force_values, 
+        extrapolation_bc = Interpolations.Flat() # only used for extrap_value_past, future == NaN
     )
 
     return nothing
 end
 
 function do_force_spreadsheet(m::PB.ReactionMethod, (var_tforce, var_FORCE), cellrange::PB.AbstractCellRange, deltat)
-    interp_FORCE = m.p
+    rj = m.reaction
 
     tforce = var_tforce[]
-    var_FORCE[] =  interp_FORCE(tforce)
+
+    if tforce < first(rj.force_times) && !isnan(rj.pars.extrap_value_past.v)
+        var_FORCE[] = rj.pars.extrap_value_past.v
+    elseif tforce > last(rj.force_times) && !isnan(rj.pars.extrap_value_future.v)
+        var_FORCE[] = rj.pars.extrap_value_future.v
+    else
+        # extrapolation_bc = Flat() will extrapolate to first/last constant value
+        var_FORCE[] =  rj.interp_FORCE(tforce)
+    end
 
     return nothing
 end
 
+function _read_xlsx(forcingfile; sheetname="Sheet1")
+
+    @info "    read_xlsx: spreadsheet $(forcingfile) sheet $(sheetname)"
+
+    xf = XLSX.readxlsx(forcingfile)
+
+    # Read using built-in iterator - assumes single header row, automatic type conversion
+    df = XLSX.eachtablerow(xf[sheetname]) |> DataFrames.DataFrame
+
+    @info "    read_xlsx read $(DataFrames.describe(df))"
+
+    return df
+end
 
 end # module
diff --git a/src/Forcings/LipForcing.jl b/src/Forcings/LipForcing.jl
@@ -88,13 +88,13 @@ read LIP data from spreadsheet
 """
 function read_lips_xlsx(lipfile)
 
-    @info "read_lips_xlsx: spreadsheet $(lipfile)"
+    @info "    read_lips_xlsx: spreadsheet $(lipfile)"
 
     xf = XLSX.readxlsx(lipfile)
     sh = xf["Sheet1"]
 
     column_names = sh["A1:I1"]
-    @info "read_lips_xlsx: column_names: $(column_names)"
+    @info "    read_lips_xlsx: column_names: $(column_names)"
     expected_column_names = ["Age" "Name" "Type" "Continental areas" "All volumes" "CO2 release (min)" "CO2 release (max)" "Degassing duration" "Present day area"]
     column_names == expected_column_names || error("spreadsheet column names don't match expected_column_names=", expected_column_names)
 
@@ -229,16 +229,17 @@ If Parameter `co2releasefield` != `NoCO2`, a CO2 pulse (duration 2e5yr, magnitud
 # Parameters
 $(PARS)
 
-# Methods and Variables
+# Methods and Variables for default Parameters
 $(METHODS_DO)
 """
 Base.@kwdef mutable struct ReactionForce_LIPs{P} <:  PB.AbstractReaction
     base::PB.ReactionBase
 
     pars::P = PB.ParametersTuple(
+        PB.ParString("datafolder", PALEOcopse.Forcings.srcdir(),
+            description="folder for forcing data 'datafile'"),
         PB.ParString("datafile", "CR_force_LIP_table.xlsx", 
-            description="spreadsheet with forcing data (path relative to $(PALEOcopse.Forcings.srcdir()))"),
-
+            description="spreadsheet with forcing data (path relative to 'datafolder')"),    
         PB.ParDouble("present_day_CFB_area", 4.8e6, units="km^2", 
             description="present day continental flood basalt area"),
         PB.ParBool("smoothempl", false, 
@@ -260,14 +261,10 @@ end
 
 
 function PB.register_methods!(rj::ReactionForce_LIPs)
-
-    lipfile = joinpath(PALEOcopse.Forcings.srcdir(), rj.pars.datafile.v)
-    @info "register_methods! ReactionForce_LIPs: $(PB.fullname(rj)) loading LIP forcing from datafile $(lipfile)"
-
-    rj.LIP_data = read_lips_xlsx(lipfile)
-
+
     CIsotopeType = rj.pars.CIsotope.v
-    @info "    CIsotopeType=$(CIsotopeType)"
+    PB.setfrozen!(rj.pars.CIsotope)
+    @info "register_methods! ReactionForce_LIPs: $(PB.fullname(rj)) CIsotopeType=$(CIsotopeType)"
 
     vars = [
         PB.VarDepScalar("tforce",     "yr",  "historical time at which to apply forcings, present = 0 yr"),
@@ -277,7 +274,6 @@ function PB.register_methods!(rj::ReactionForce_LIPs)
         PB.VarContribScalar("flux_C"=>"fluxSedCrusttoAOcean.flux_C",  "mol yr-1",  "LIP CO2 release", 
             attributes=(:field_data=>CIsotopeType,)),
     ]
-
 
     PB.add_method_setup!(
         rj,
@@ -295,16 +291,20 @@ function PB.register_methods!(rj::ReactionForce_LIPs)
     return nothing
 end
 
-function setup_force_LIPs(m::PB.ReactionMethod, (), cellrange::PB.AbstractCellRange, attribute_value)
+function setup_force_LIPs(m::PB.ReactionMethod, (), cellrange::PB.AbstractCellRange, attribute_name)
     rj = m.reaction
 
-    attribute_value == :initial_value || return nothing
+    attribute_name == :setup || return nothing
+
+    lipfile = joinpath(rj.pars.datafolder.v, rj.pars.datafile.v)
+    @info "setup_force_LIPs! ReactionForce_LIPs: $(PB.fullname(rj)) loading LIP forcing from 'datafolder/datafile'='$(lipfile)'"
+
+    rj.LIP_data = read_lips_xlsx(lipfile)
 
     lipkwargs = (smoothempl=rj.pars.smoothempl.v, )
 
     rj.default_lambda = find_default_lambda(rj.pars.present_day_CFB_area.v, rj.LIP_data, lipkwargs )
 
-    @info "ReactionForce_LIPs.setup_force_LIPs $(PB.fullname(rj))"
     @info "    found default_lambda=$(rj.default_lambda) "*
         "to match present_day_CFB_area = $(rj.pars.present_day_CFB_area.v) km^2"
 

diff --git a/src/PALEOcopse.jl b/src/PALEOcopse.jl
@@ -14,4 +14,25 @@ include("oceanfloor/Oceanfloor.jl")
 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()
+
 end
diff --git a/src/oceanfloor/Weathering.jl b/src/oceanfloor/Weathering.jl
@@ -60,7 +60,7 @@ end
 function PB.register_methods!(rj::ReactionSeafloorWeathering)
 
     CIsotopeType = rj.pars.CIsotope.v
-
+    PB.setfrozen!(rj.pars.CIsotope)
     vars = [
         PB.VarDep("global.RHOSFW",          "", "seafloor weathering-specific additional forcing (usually 1.0)"),
         PB.VarDep("(global.TEMP)",          "K", "global mean temperature"),
@@ -208,10 +208,12 @@ function set_distributionfromdepths(
     m::PB.ReactionMethod, 
     (vars, ), 
     cellrange::PB.AbstractCellRange, 
-    attribute_value
+    attribute_name
 )
     rj = m.reaction
-    attribute_value == :initial_value || return nothing
+    attribute_name == :setup || return nothing
+
+    @info "$(fullname(m)) ReactionSeafloorWeathering:"
 
     PB.get_length(rj.domain) == length(cellrange.indices) || 
         error("ReactionSeafloorWeathering $(PB.fullname(rj)) set_distributionfromdepths cellrange does not cover whole Domain")
@@ -230,8 +232,7 @@ function set_distributionfromdepths(
 
     normsum = sum(sfw_distribution)
     Atot = sum(vars.Afloor)
-    @info "ReactionSeafloorWeathering $(PB.fullname(rj)) "*
-        "sfw distributed among $nsfw of $(length(sfw_distribution)) boxes, area $normsum m^2 of $Atot m^2"
+    @info "    sfw distributed among $nsfw of $(length(sfw_distribution)) boxes, area $normsum m^2 of $Atot m^2"
     PB.setvalue!(rj.pars.sfw_distribution, sfw_distribution./normsum)
 
     return nothing

diff --git a/test/forcings/PlotForcings.ipynb b/test/forcings/PlotForcings.ipynb
@@ -98,7 +98,7 @@
     "# Initialize model Reactions (will read forcing data from files etc)\n",
     "PB.initialize_reactiondata!(model, modeldata)\n",
     "\n",
-    "PB.dispatch_setup(model, :initial_value, modeldata)\n",
+    "PB.dispatch_setup(model, :setup, modeldata)\n",
     "\n",
     "dispatchlists = modeldata.dispatchlists_all\n"
    ]

diff --git a/test/forcings/runforcingtests.jl b/test/forcings/runforcingtests.jl
@@ -34,7 +34,7 @@ import PALEOcopse
     PB.initialize_reactiondata!(model, modeldata)    
 
     @info "dispatch_setup"
-    PB.dispatch_setup(model, :initial_value, modeldata)
+    PB.dispatch_setup(model, :setup, modeldata)
 
     dispatchlists = modeldata.dispatchlists_all
 

diff --git a/test/global/runtemperaturetests.jl b/test/global/runtemperaturetests.jl
@@ -33,6 +33,8 @@ import PALEOcopse
 
     # check Reaction configuration
     PB.check_configuration(model)
+    # Initialise Reactions and non-state variables
+    PB.dispatch_setup(model, :setup, modeldata)
     # Initialise state variables to norm_value
     PB.dispatch_setup(model, :norm_value, modeldata)
     PB.copy_norm!(modeldata.solver_view_all)