From 94720260174c3c36991d6101f0f547e93015b173 Mon Sep 17 00:00:00 2001
From: Brian Groenke <brian.groenke@awi.de>
Date: Mon, 25 Nov 2024 16:54:57 +0100
Subject: [PATCH] Major refactoring of forcing/input handling
Input functions (e.g. forcings) are now defined symbolically in layors
via the `Input` type and later resolved at runtime.
This allows inputs to always be interchangeable with parameters or
dynamic parameterizations. It should also make it easier to swap
out forcing data since it is all sourced from a single place now
rather than spread out among different model components.
Other changes:
- Renamed `computefluxes!` to `computeprognostic!`
- Renamed `resolve` to `materialize` for clarity
Bumped minor version to 0.23
---
Project.toml | 4 +-
README.md | 2 +-
docs/src/manual/overview.md | 3 +-
docs/src/quickstart.md | 2 +-
examples/cglite_parameter_ensembles.jl | 2 +-
examples/heat_freeW_bucketW_samoylov.jl | 9 +-
examples/heat_freeW_lake_lite_implicit.jl | 6 +-
examples/heat_freeW_lite_implicit.jl | 10 +-
examples/heat_freeW_samoylov.jl | 21 ++-
.../heat_freeW_seb_snow_bucketW_samoylov.jl | 32 +++-
examples/heat_freeW_snow_samoylov.jl | 8 +-
examples/heat_sfcc_constantbc.jl | 1 +
examples/heat_sfcc_richardseq_samoylov.jl | 8 +-
examples/heat_sfcc_samoylov.jl | 6 +-
examples/heat_simple_autodiff_grad.jl | 5 +-
src/CryoGrid.jl | 15 ++
src/IO/InputOutput.jl | 9 +-
src/IO/forcings/forcings.jl | 154 +-----------------
src/IO/forcings/forcings_json.jl | 10 +-
src/IO/forcings/forcings_ncd.jl | 4 +-
src/IO/forcings/providers.jl | 53 ++++++
src/IO/input.jl | 85 ++++++++++
src/Physics/Heat/heat_bc.jl | 39 ++---
src/Physics/Heat/heat_conduction.jl | 52 ++++--
src/Physics/Heat/heat_implicit.jl | 42 +++--
src/Physics/Heat/heat_methods.jl | 1 +
src/Physics/Heat/heat_types.jl | 2 -
src/Physics/Heat/heat_water.jl | 18 +-
src/Physics/Hydrology/water_ET.jl | 12 +-
src/Physics/Hydrology/water_types.jl | 14 +-
src/Physics/Salt/salt_bc.jl | 10 +-
src/Physics/Soils/para/simple.jl | 23 +--
src/Physics/Soils/soil_heat.jl | 46 ------
src/Physics/Surface/SEB/seb.jl | 30 ++--
src/Physics/Surface/SEB/seb_heat.jl | 10 +-
src/Physics/Surface/SEB/seb_state.jl | 33 ++--
src/Physics/Surface/SWB/swb.jl | 35 +---
src/Physics/steplimiters.jl | 1 -
src/Presets/Presets.jl | 28 ++--
src/Solvers/LiteImplicit/cglite_step.jl | 2 +-
src/Solvers/LiteImplicit/cglite_types.jl | 2 +-
src/Solvers/basic_solvers.jl | 4 +-
src/Tiles/Tiles.jl | 2 +-
src/Tiles/tile.jl | 130 +++++++++------
.../presetgrids.jl => default_grids.jl} | 0
src/initializers.jl | 3 -
src/methods.jl | 26 ---
test/Physics/Surface/seb_tests.jl | 42 ++---
test/benchmarks.jl | 4 +-
49 files changed, 502 insertions(+), 558 deletions(-)
create mode 100644 src/IO/forcings/providers.jl
create mode 100644 src/IO/input.jl
rename src/{Presets/presetgrids.jl => default_grids.jl} (100%)
diff --git a/Project.toml b/Project.toml
index d831f333..bdc88c13 100755
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
name = "CryoGrid"
uuid = "a535b82e-5f3d-4d97-8b0b-d6483f5bebd5"
authors = ["Brian Groenke <brian.groenke@awi.de>", "Jan Nitzbon <jan.nitzbon@awi.de>", "Moritz Langer <moritz.langer@awi.de>"]
-version = "0.22.3"
+version = "0.23.0"
[deps]
Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
@@ -85,7 +85,7 @@ Tables = "1"
TimeSeries = "0.23, 0.24"
UnPack = "1"
Unitful = "1"
-julia = "1.6, 1.7, 1.8, 1.9, 1.10"
+julia = "1.6, 1.7, 1.8, 1.9, 1.10, 1.11"
[extesnions]
CryoGridSBIExt = ["SimulationBasedInference"]
diff --git a/README.md b/README.md
index d9fb6fb0..7d41667a 100755
--- a/README.md
+++ b/README.md
@@ -49,7 +49,7 @@ using Plots: plot
# load provided forcing data from Samoylov;
# The forcing file will be automatically downloaded to the input/ folder if not already present.
-forcings = loadforcings(CryoGrid.Presets.Forcings.Samoylov_ERA_obs_fitted_1979_2014_spinup_extended_2044);
+forcings = loadforcings(CryoGrid.Forcings.Samoylov_ERA_obs_fitted_1979_2014_spinup_extended_2044);
# get preset soil and initial temperature profile for Samoylov
soilprofile, tempprofile = CryoGrid.Presets.SamoylovDefault
initT = initializer(:T, tempprofile)
diff --git a/docs/src/manual/overview.md b/docs/src/manual/overview.md
index e58d8030..03357710 100644
--- a/docs/src/manual/overview.md
+++ b/docs/src/manual/overview.md
@@ -12,8 +12,9 @@ At the highest level, a model in `CryoGrid.jl` is defined by one or more [`Tile`
```julia
# ... load forcings, set up profiles, etc.
# see examples/heat_vgfc_seb_saoylov_custom.jl for more details
+forcings = (;Tair,pr,q,wind,Lin,Sin,z)
strat = Stratigraphy(
- -2.0u"m" => Top(SurfaceEnergyBalance(Tair,pr,q,wind,Lin,Sin,z)),
+ -2.0u"m" => Top(SurfaceEnergyBalance(forcings)),
0.0u"m" => Ground(soilprofile, HeatBalance(:H; freezecurve=DallAmico())),
1000.0u"m" => Bottom(GeothermalHeatFlux(0.053u"J/s/m^2"))
);
diff --git a/docs/src/quickstart.md b/docs/src/quickstart.md
index e3b29ccf..2d887e00 100644
--- a/docs/src/quickstart.md
+++ b/docs/src/quickstart.md
@@ -8,7 +8,7 @@ using Plots
# load provided forcing data from Samoylov;
# The forcing file will be automatically downloaded to the input/ folder if not already present.
-forcings = loadforcings(CryoGrid.Presets.Forcings.Samoylov_ERA_obs_fitted_1979_2014_spinup_extended_2044, :Tair => u"°C");
+forcings = loadforcings(CryoGrid.Forcings.Samoylov_ERA_obs_fitted_1979_2014_spinup_extended_2044, :Tair => u"°C");
# get preset soil and initial temperature profile for Samoylov
soilprofile, tempprofile = CryoGrid.Presets.SamoylovDefault
initT = initializer(:T, tempprofile)
diff --git a/examples/cglite_parameter_ensembles.jl b/examples/cglite_parameter_ensembles.jl
index 1eec62a5..91bb16c8 100644
--- a/examples/cglite_parameter_ensembles.jl
+++ b/examples/cglite_parameter_ensembles.jl
@@ -15,7 +15,7 @@ if Threads.nthreads() == 1
end
# Load forcings and build stratigraphy like before.
-forcings = loadforcings(CryoGrid.Presets.Forcings.Samoylov_ERA_MkL3_CCSM4_long_term);
+forcings = loadforcings(CryoGrid.Forcings.Samoylov_ERA_MkL3_CCSM4_long_term);
soilprofile = SoilProfile(
0.0u"m" => SimpleSoil(por=Param(0.80, prior=Uniform(0.65,0.95)),sat=1.0,org=0.75),
0.1u"m" => SimpleSoil(por=Param(0.80, prior=Uniform(0.65,0.95)),sat=1.0,org=0.25),
diff --git a/examples/heat_freeW_bucketW_samoylov.jl b/examples/heat_freeW_bucketW_samoylov.jl
index a9380062..96d67625 100644
--- a/examples/heat_freeW_bucketW_samoylov.jl
+++ b/examples/heat_freeW_bucketW_samoylov.jl
@@ -5,11 +5,14 @@
# Frist, load forcings and define boundary conditions.
using CryoGrid
-forcings = loadforcings(CryoGrid.Presets.Forcings.Samoylov_ERA_obs_fitted_1979_2014_spinup_extended_2044);
+forcings = loadforcings(CryoGrid.Forcings.Samoylov_ERA_obs_fitted_1979_2014_spinup_extended_2044);
_, tempprofile = CryoGrid.Presets.SamoylovDefault;
initT = initializer(:T, tempprofile)
initsat = initializer(:sat, (l,state) -> state.sat .= l.para.sat);
-upperbc = WaterHeatBC(SurfaceWaterBalance(forcings), TemperatureBC(forcings.Tair, NFactor(0.5,0.9)));
+upperbc = WaterHeatBC(
+ SurfaceWaterBalance(),
+ TemperatureBC(Input(:Tair), NFactor(0.5,0.9))
+);
# The `@Stratigraphy` macro is just a small convenience that automatically wraps the three subsurface layers in a tuple.
# It would be equivalent to use the `Stratigraphy` constructor directly and wrap these layers in a tuple or list.
@@ -21,7 +24,7 @@ strat = @Stratigraphy(
1000.0u"m" => Bottom(GeothermalHeatFlux(0.053u"W/m^2")),
);
modelgrid = CryoGrid.Presets.DefaultGrid_2cm
-tile = Tile(strat, modelgrid, initT, initsat);
+tile = Tile(strat, modelgrid, forcings, initT, initsat);
# Now we set up the problem and solve using the integrator interface.
tspan = (DateTime(2010,5,30),DateTime(2012,8,30))
diff --git a/examples/heat_freeW_lake_lite_implicit.jl b/examples/heat_freeW_lake_lite_implicit.jl
index 7540a855..c0690472 100644
--- a/examples/heat_freeW_lake_lite_implicit.jl
+++ b/examples/heat_freeW_lake_lite_implicit.jl
@@ -6,7 +6,7 @@ Plots.plotly()
CryoGrid.debug(true)
-forcings = loadforcings(CryoGrid.Presets.Forcings.Samoylov_ERA_MkL3_CCSM4_long_term);
+forcings = loadforcings(CryoGrid.Forcings.Samoylov_ERA_MkL3_CCSM4_long_term);
soilprofile = SoilProfile(
0.0u"m" => SimpleSoil(por=0.80,sat=0.9,org=0.75),
0.1u"m" => SimpleSoil(por=0.80,sat=1.0,org=0.25),
@@ -24,7 +24,7 @@ modelgrid = Grid(vcat(-1.0u"m":0.02u"m":-0.02u"m", CryoGrid.Presets.DefaultGrid_
z_top = -1.0u"m"
z_sub = keys(soilprofile)
z_bot = modelgrid[end]
-upperbc = TemperatureBC(forcings.Tair, NFactor(nf=0.5))
+upperbc = TemperatureBC(Input(:Tair), NFactor(nf=0.5))
initT = initializer(:T, tempprofile_linear)
@info "Building stratigraphy"
heatop = Heat.EnthalpyImplicit()
@@ -35,7 +35,7 @@ strat = @Stratigraphy(
z_bot => Bottom(GeothermalHeatFlux(0.053u"W/m^2"))
);
@info "Building tile"
-tile = @time Tile(strat, modelgrid, initT)
+tile = @time Tile(strat, modelgrid, forcings, initT)
# define time span, 5 years
tspan = (DateTime(2010,12,30), DateTime(2012,12,30))
tspan_sol = convert_tspan(tspan)
diff --git a/examples/heat_freeW_lite_implicit.jl b/examples/heat_freeW_lite_implicit.jl
index 712675d6..eb5b9931 100644
--- a/examples/heat_freeW_lite_implicit.jl
+++ b/examples/heat_freeW_lite_implicit.jl
@@ -8,16 +8,16 @@ using CryoGrid
using CryoGrid.LiteImplicit
# Load forcings and build stratigraphy like before.
-forcings = loadforcings(CryoGrid.Presets.Forcings.Samoylov_ERA_MkL3_CCSM4_long_term);
+forcings = loadforcings(CryoGrid.Forcings.Samoylov_ERA_MkL3_CCSM4_long_term);
forcings = Base.rename(forcings, :Ptot => :precip)
precip_values = forcings.precip.interpolant.coefs
precip_values .*= 2
z_top = -2.0u"m"
z_bot = 1000.0u"m"
upperbc = WaterHeatBC(
- SurfaceWaterBalance(forcings),
- TemperatureBC(forcings.Tair, NFactor()),
-)
+ SurfaceWaterBalance(Input(:rainfall), Input(:snowfall)),
+ TemperatureBC(Input(:Tair), NFactor(0.5,0.9))
+);
ssinit = ThermalSteadyStateInit(T0=-15.0u"°C")
heatop = Heat.EnthalpyImplicit()
soilprofile = SoilProfile(
@@ -37,7 +37,7 @@ strat = @Stratigraphy(
z_bot => Bottom(GeothermalHeatFlux(0.053u"W/m^2"))
);
modelgrid = Grid(vcat(z_top:0.02u"m":-0.02u"m", CryoGrid.Presets.DefaultGrid_2cm))
-tile = Tile(strat, modelgrid, ssinit);
+tile = Tile(strat, modelgrid, forcings, ssinit);
# Since the solver can take daily timesteps, we can easily specify longer simulation time spans at minimal cost.
# Here we specify a time span of 30 years.
diff --git a/examples/heat_freeW_samoylov.jl b/examples/heat_freeW_samoylov.jl
index 8ba3708c..e2eb192b 100644
--- a/examples/heat_freeW_samoylov.jl
+++ b/examples/heat_freeW_samoylov.jl
@@ -9,10 +9,7 @@
# First we load the built-in forcing file from Nitzbon et al. 2020 (CryoGrid 3). Note that this will
# download the forcing files from the AWI NextCloud if they are not already present in the `input/` folder.
using CryoGrid
-forcings = loadforcings(CryoGrid.Presets.Forcings.Samoylov_ERA_obs_fitted_1979_2014_spinup_extended_2044);
-
-# We choose the default grid discretization with 5 cm spacing at the surface.
-grid = CryoGrid.Presets.DefaultGrid_5cm;
+forcings = loadforcings(CryoGrid.Forcings.Samoylov_ERA_obs_fitted_1979_2014_spinup_extended_2044);
# We use a simple 5-layer stratigraphy suitable for Samoylov. This is based on the
# profile provided in `Presets` but uses the default "free water" freezing characteristic
@@ -23,15 +20,22 @@ soilprofile = SoilProfile(
0.4u"m" => SimpleSoil(; por=0.55, org=0.25),
3.0u"m" => SimpleSoil(; por=0.50, org=0.0),
10.0u"m" => SimpleSoil(; por=0.30, org=0.0),
-)
+);
# We construct a state variable initializer for temperature `T` from the temperature profile preset for Samoylov.
-initT = initializer(:T, CryoGrid.Presets.SamoylovDefault.tempprofile)
+initT = initializer(:T, CryoGrid.Presets.SamoylovDefault.tempprofile);
+
+# We choose the default grid discretization with 5 cm spacing at the surface.
+grid = CryoGrid.DefaultGrid_5cm;
+
+# Now we construct the Tile using the built-in model configuration `SoilHeatTile` which defines a
+# standalone soil straigraphy with only heat conduction and no water flow.
tile = CryoGrid.Presets.SoilHeatTile(
:H,
- TemperatureBC(forcings.Tair, NFactor(nf=0.6)),
+ TemperatureBC(Input(:Tair), NFactor(nf=Param(0.6), nt=Param(0.9))),
GeothermalHeatFlux(0.053u"W/m^2"),
soilprofile,
+ forcings,
initT;
grid=grid
);
@@ -56,3 +60,6 @@ import Plots
zs = [1,10,20,30,50,100,200,500,1000]u"cm"
cg = Plots.cgrad(:copper,rev=true);
Plots.plot(out.T[Z(Near(zs))], color=cg[LinRange(0.0,1.0,length(zs))]', ylabel="Temperature", leg=false, size=(800,500), dpi=150)
+
+using BenchmarkTools
+@profview @btime $tile($du0, $u0, $prob.p, $prob.tspan[1])
diff --git a/examples/heat_freeW_seb_snow_bucketW_samoylov.jl b/examples/heat_freeW_seb_snow_bucketW_samoylov.jl
index 12f7917b..4248173e 100644
--- a/examples/heat_freeW_seb_snow_bucketW_samoylov.jl
+++ b/examples/heat_freeW_seb_snow_bucketW_samoylov.jl
@@ -18,12 +18,11 @@ soilprofile = SoilProfile(
);
## mid-winter temperature profile
tempprofile = CryoGrid.Presets.SamoylovDefault.tempprofile
-forcings = loadforcings(CryoGrid.Presets.Forcings.Samoylov_ERA_obs_fitted_1979_2014_spinup_extended_2044);
+forcings = loadforcings(CryoGrid.Forcings.Samoylov_ERA_obs_fitted_1979_2014_spinup_extended_2044);
tempprofile = CryoGrid.Presets.SamoylovDefault.tempprofile
initT = initializer(:T, tempprofile)
-z = 2.0u"m"; # height [m] for which the forcing variables (Temp, humidity, wind, pressure) are provided
-seb = SurfaceEnergyBalance(forcings, z)
-swb = SurfaceWaterBalance(forcings)
+seb = SurfaceEnergyBalance()
+swb = SurfaceWaterBalance()
upperbc = WaterHeatBC(swb, seb)
heat = HeatBalance()
water = WaterBalance(BucketScheme(), DampedET())
@@ -35,7 +34,7 @@ strat = @Stratigraphy(
1000.0u"m" => Bottom(GeothermalHeatFlux(0.053u"J/s/m^2")),
);
## create Tile
-tile = Tile(strat, modelgrid, initT);
+tile = Tile(strat, modelgrid, forcings, initT);
# Set up the problem and solve it!
tspan = (DateTime(2010,10,30), DateTime(2011,10,30))
@@ -87,3 +86,26 @@ Plots.plot(ustrip.(cumsum(out.top.Qg, dims=2)), color=cg[LinRange(0.0,1.0,length
# Integratoed ground latent heat flux:
Plots.plot(ustrip.(cumsum(out.top.Qe, dims=2)), color=cg[LinRange(0.0,1.0,length(zs))]', ylabel="Integrated ground heat flux", leg=false, size=(800,500), dpi=150)
+
+using BenchmarkTools
+using Interpolations
+
+function interpolant(xs::AbstractVector, ts::AbstractVector)
+ f = interpolate((ts,), xs, Gridded(Linear()))
+ return f
+end
+
+function eval(xs::AbstractVector, ts::AbstractVector, t)
+ f = interpolate((ts,), xs, Gridded(Linear()))
+ return f(t)
+end
+
+const xdata = randn(1000)
+const tdata = 1:1000
+@btime eval(xdata, tdata, 10.0)
+
+const itp = interpolant(xdata, tdata)
+@btime itp(10.0)
+
+x = randn(1000,2)
+interpolate((tdata,1:2), x, (Gridded(Linear()), NoInterp()))
\ No newline at end of file
diff --git a/examples/heat_freeW_snow_samoylov.jl b/examples/heat_freeW_snow_samoylov.jl
index 80aaef5e..2a423222 100644
--- a/examples/heat_freeW_snow_samoylov.jl
+++ b/examples/heat_freeW_snow_samoylov.jl
@@ -7,7 +7,7 @@ using CryoGrid
using OrdinaryDiffEq
# First we set up the model:
-forcings = loadforcings(CryoGrid.Presets.Forcings.Samoylov_ERA_obs_fitted_1979_2014_spinup_extended_2044);
+forcings = loadforcings(CryoGrid.Forcings.Samoylov_ERA_obs_fitted_1979_2014_spinup_extended_2044);
soilprofile = SoilProfile(
0.0u"m" => SimpleSoil(por=0.80,sat=1.0,org=0.75),
0.1u"m" => SimpleSoil(por=0.80,sat=1.0,org=0.25),
@@ -21,8 +21,8 @@ z_top = -2.0u"m"
z_sub = keys(soilprofile)
z_bot = 1000.0u"m"
upperbc = WaterHeatBC(
- SurfaceWaterBalance(forcings),
- TemperatureBC(forcings.Tair)
+ SurfaceWaterBalance(),
+ TemperatureBC(Input(:Tair))
)
snowpack = Snowpack(
para=Snow.Bulk(),
@@ -40,7 +40,7 @@ strat = @Stratigraphy(
z_bot => Bottom(GeothermalHeatFlux(0.053u"J/s/m^2"))
);
modelgrid = CryoGrid.Presets.DefaultGrid_5cm
-tile = Tile(strat, modelgrid, initT, initsat)
+tile = Tile(strat, modelgrid, forcings, initT, initsat)
# define time span, 2 years + 3 months
tspan = (DateTime(2010,9,30), DateTime(2012,9,30))
u0, du0 = @time initialcondition!(tile, tspan)
diff --git a/examples/heat_sfcc_constantbc.jl b/examples/heat_sfcc_constantbc.jl
index 2f7bf0fb..11cee116 100644
--- a/examples/heat_sfcc_constantbc.jl
+++ b/examples/heat_sfcc_constantbc.jl
@@ -37,6 +37,7 @@ tile = CryoGrid.Presets.SoilHeatTile(
ConstantBC(HeatBalance, CryoGrid.Neumann, 10.0u"W/m^2"),
ConstantBC(HeatBalance, CryoGrid.Neumann, 10.0u"W/m^2"),
soilprofile,
+ inputs(),
initT;
grid,
);
diff --git a/examples/heat_sfcc_richardseq_samoylov.jl b/examples/heat_sfcc_richardseq_samoylov.jl
index bbdf09cc..36aada34 100644
--- a/examples/heat_sfcc_richardseq_samoylov.jl
+++ b/examples/heat_sfcc_richardseq_samoylov.jl
@@ -5,9 +5,9 @@
# pressure gradients.
using CryoGrid
-forcings = loadforcings(CryoGrid.Presets.Forcings.Samoylov_ERA_obs_fitted_1979_2014_spinup_extended_2044);
+forcings = loadforcings(CryoGrid.Forcings.Samoylov_ERA_obs_fitted_1979_2014_spinup_extended_2044);
tspan = (DateTime(2011,1,1),DateTime(2011,12,31))
-T0 = values(forcings.Tair(tspan[1]))[1]
+T0 = forcings.Tair(tspan[1])
tempprofile = TemperatureProfile(
0.0u"m" => 0.0*u"°C",
1.0u"m" => -8.0u"°C",
@@ -24,7 +24,7 @@ waterflow = RichardsEq(;swrc);
# We use the enthalpy-based heat diffusion with high accuracy Newton-based solver for inverse enthalpy mapping
heatop = Heat.Diffusion1D(:H)
-upperbc = WaterHeatBC(SurfaceWaterBalance(forcings), TemperatureBC(forcings.Tair, NFactor(nf=0.6, nt=0.9)));
+upperbc = WaterHeatBC(SurfaceWaterBalance(), TemperatureBC(Input(:Tair), NFactor(nf=0.6, nt=0.9)));
# We will use a simple stratigraphy with three subsurface soil layers.
# Note that the @Stratigraphy macro lets us list multiple subsurface layers without wrapping them in a tuple.
@@ -39,7 +39,7 @@ strat = @Stratigraphy(
1000.0u"m" => Bottom(GeothermalHeatFlux(0.053u"W/m^2"))
);
grid = CryoGrid.Presets.DefaultGrid_2cm
-tile = Tile(strat, grid, initT);
+tile = Tile(strat, grid, forcings, initT);
u0, du0 = @time initialcondition!(tile, tspan)
prob = CryoGridProblem(tile, u0, tspan, saveat=3*3600, savevars=(:T,:θw,:θwi,:kw));
integrator = init(prob, CGEuler())
diff --git a/examples/heat_sfcc_samoylov.jl b/examples/heat_sfcc_samoylov.jl
index 72bf5839..3b6b8119 100644
--- a/examples/heat_sfcc_samoylov.jl
+++ b/examples/heat_sfcc_samoylov.jl
@@ -6,13 +6,13 @@
using CryoGrid
# First we set up the soil heat model:
-forcings = loadforcings(CryoGrid.Presets.Forcings.Samoylov_ERA5_fitted_daily_1979_2020);
+forcings = loadforcings(CryoGrid.Forcings.Samoylov_ERA5_fitted_daily_1979_2020);
grid = CryoGrid.Presets.DefaultGrid_5cm
soilprofile, tempprofile = CryoGrid.Presets.SamoylovDefault
initT = initializer(:T, tempprofile)
-upperbc = TemperatureBC(forcings.Tair, NFactor(nf=0.6, nt=0.9))
+upperbc = TemperatureBC(Input(:Tair), NFactor(nf=0.6, nt=0.9))
lowerbc = GeothermalHeatFlux(0.053u"W/m^2")
-tile = CryoGrid.Presets.SoilHeatTile(upperbc, lowerbc, soilprofile, initT; grid=grid)
+tile = CryoGrid.Presets.SoilHeatTile(upperbc, lowerbc, soilprofile, forcings, initT; grid=grid)
tspan = (DateTime(2010,10,30),DateTime(2011,10,30))
u0, du0 = @time initialcondition!(tile, tspan)
prob = CryoGridProblem(tile, u0, tspan, saveat=3*3600.0, savevars=(:T,));
diff --git a/examples/heat_simple_autodiff_grad.jl b/examples/heat_simple_autodiff_grad.jl
index 5a1e389e..5b88521e 100644
--- a/examples/heat_simple_autodiff_grad.jl
+++ b/examples/heat_simple_autodiff_grad.jl
@@ -6,15 +6,16 @@
# Set up forcings and boundary conditions similarly to other examples:
using CryoGrid
-forcings = loadforcings(CryoGrid.Presets.Forcings.Samoylov_ERA_obs_fitted_1979_2014_spinup_extended_2044);
+forcings = loadforcings(CryoGrid.Forcings.Samoylov_ERA_obs_fitted_1979_2014_spinup_extended_2044);
soilprofile, tempprofile = CryoGrid.Presets.SamoylovDefault
grid = CryoGrid.Presets.DefaultGrid_5cm
initT = initializer(:T, tempprofile)
tile = CryoGrid.Presets.SoilHeatTile(
:T,
- TemperatureBC(forcings.Tair, NFactor(nf=Param(0.5), nt=Param(0.9))),
+ TemperatureBC(Input(:Tair), NFactor(nf=Param(0.5), nt=Param(0.9))),
GeothermalHeatFlux(0.053u"W/m^2"),
soilprofile,
+ forcings,
initT;
grid=grid
)
diff --git a/src/CryoGrid.jl b/src/CryoGrid.jl
index c3dfbf78..bad4afaa 100755
--- a/src/CryoGrid.jl
+++ b/src/CryoGrid.jl
@@ -69,6 +69,8 @@ export DiscretizationStrategy, AutoGrid, PresetGrid, LinearSpacing, Grid, cells,
include("Numerics/Numerics.jl")
using .Numerics
+include("default_grids.jl")
+
export initializer
include("initializers.jl")
@@ -100,6 +102,19 @@ include("Solvers/Solvers.jl")
# Presets
include("Presets/Presets.jl")
+using .InputOutput: Resource
+
+Forcings = (
+ Samoylov_ERA5_fitted_daily_1979_2020 = Resource("samoylov_era5_fitted_daily_1979-2020", ForcingFormatJSON{2}(), "https://nextcloud.awi.de/s/ScYAoHzeMzAfpjf/download/samoylov_era5_fitted_daily_1979-2020.json"),
+ Samoylov_ERA_obs_fitted_1979_2014_spinup_extended_2044 = Resource("Samoylov_ERA_obs_fitted_1979_2014_spinup_extended_2044", ForcingFormatJSON{1}(), "https://nextcloud.awi.de/s/cbeycGQoQpXi3Ei/download/samoylov_ERA_obs_fitted_1979_2014_spinup_extended2044.json"),
+ Samoylov_ERA_MkL3_CCSM4_long_term = Resource("Samoylov_ERA_MkL3_CCSM4_long_term", ForcingFormatJSON{1}(), "https://nextcloud.awi.de/s/45ax9AsTACxL25Q/download/FORCING_ULC_126_72.json"),
+ Bayelva_ERA5_fitted_daily_1979_2020 = Resource("bayelva_era5_fitted_daily_1979-2020", ForcingFormatJSON{2}(), "https://nextcloud.awi.de/s/5AdbRMYKneCHgx4/download/bayelva_era5_fitted_daily_1979-2020.json")
+)
+Parameters = (
+ # Faroux et al. doi:10.1109/IGARSS.2007.4422971
+ EcoCLimMap_ULC_126_72 = Resource("EcoCLimMap_ULC_126_72", ParamsJSON{1}(), "https://nextcloud.awi.de/s/7F65JET9TzdosMD/download/PARA_ULC_126_72.json")
+)
+
function __init__()
# SimulationBasedInference extension
@require SimulationBasedInference="78927d98-f421-490e-8789-96b006983a5c" begin
diff --git a/src/IO/InputOutput.jl b/src/IO/InputOutput.jl
index ec01a181..b055d861 100644
--- a/src/IO/InputOutput.jl
+++ b/src/IO/InputOutput.jl
@@ -12,6 +12,7 @@ using Dates
using DimensionalData
using Downloads
using Flatten
+using Interpolations
using ModelParameters
using Tables
using Unitful
@@ -62,10 +63,12 @@ include("params/params.jl")
export ParamsJSON, ParamsYAML
include("params/params_loaders.jl")
-export Forcings, Forcing, ConstantForcing, TransformedForcing, TimeVaryingForcing, InterpolatedForcing
-export TemperatureForcing, WindForcing, HumidityForcing, EnergyFluxForcing, PressureForcing, VelocityForcing # aliases
+export Input, InputProvider, InputFunctionProvider
+export inputs
+include("input.jl")
+
export ForcingFormat, ForcingFormatJSON, ForcingFormatNCD
-export loadforcings, time_derivative_forcing
+export loadforcings
include("forcings/forcings.jl")
export CryoGridOutput, write_netcdf!
diff --git a/src/IO/forcings/forcings.jl b/src/IO/forcings/forcings.jl
index 933106ae..cda88094 100644
--- a/src/IO/forcings/forcings.jl
+++ b/src/IO/forcings/forcings.jl
@@ -1,154 +1,9 @@
-"""
- Forcing{unit,T}
-
-Abstract type representing a generic external forcing term.
-"""
-abstract type Forcing{unit,T} end
-Base.nameof(f::Forcing) = f.name
-CryoGrid.parameterize(f::Forcing; ignored...) = f
-@propagate_inbounds (forcing::Forcing)(x::Number) = error("$(typeof(forcing)) not implemented")
-@propagate_inbounds (forcing::Forcing)(t::DateTime) = forcing(convert_t(t))
-
-ConstructionBase.constructorof(::Type{F}) where {unit,F<:Forcing{unit}} = (args...) -> F.name.wrapper(unit, args...)
-
"""
Represents an externally specified format for forcing inputs. IO functions should dispatch on
specific types `T<:ForcingFormat` that they implement.
"""
abstract type ForcingFormat end
-# Aliases for forcing types
-const TemperatureForcing{T} = Forcing{u"°C",T} where {T}
-const VelocityForcing{T} = Forcing{u"m/s",T} where {T}
-const HeightForcing{T} = Forcing{u"m",T} where {T}
-const HumidityForcing{T} = Forcing{u"kg/kg",T} where {T}
-const PressureForcing{T} = Forcing{upreferred(u"Pa"),T} where {T}
-const EnergyFluxForcing{T} = Forcing{upreferred(u"W/m^2"),T} where {T}
-
-"""
- ConstantForcing{unit,T} <: Forcing{unit,T}
-
-Simple `Forcing` type that just returns a constant value. This type is primarily
-intended for writing tests but could still be used in applications if necessary.
-"""
-struct ConstantForcing{unit,T} <: Forcing{unit,T}
- value::T
- name::Symbol
- ConstantForcing(unit::Unitful.Units, value::T, name::Symbol) where {T} = new{unit,T}(value, name)
- ConstantForcing(qty::Unitful.AbstractQuantity, name::Symbol) = new{unit(qty),typeof(qty)}(qty, name)
- ConstantForcing(qty::Number, name::Symbol) = new{Unitful.NoUnits,typeof(qty)}(qty, name)
-end
-
-(f::ConstantForcing)(::Number) = f.value
-
-"""
- TimeVaryingForcing{unit,T} <: Forcing{unit,T}
-
-Simple `Forcing` type that evaluates a given function `f(t)` of where `t` is the current time.
-"""
-struct TimeVaryingForcing{unit,F,T} <: Forcing{unit,T}
- f::F
- name::Symbol
- TimeVaryingForcing(unit, ::Type{T}, f::F, name::Symbol) where {F,T} = new{unit,F,T}(f, name)
- function TimeVaryingForcing(f::F, name::Symbol; initial_value::T=f(0.0)) where {F,T}
- return new{unit(initial_value),F,typeof(ustrip(initial_value))}(f, name)
- end
-end
-
-ConstructionBase.constructorof(::Type{TimeVaryingForcing{unit,F,T}}) where {unit,F,T} = (f, name) -> TimeVaryingForcing(unit, T, f, name)
-
-(f::TimeVaryingForcing)(t::Number) = f.f(t)
-
-"""
- TransformedForcing(transformed_unit,orig_unit,TF,T,TO<:Forcing{orig_unit,T}) <: Forcing{transformed_unit,T}
-
-Wraps another `Forcing` and applies an arbitrary transformation `f` when evaluated. The transformed unit
-can either be automatically determined on construction or specified directly.
-"""
-struct TransformedForcing{transformed_unit,orig_unit,TF,T,TO<:Forcing{orig_unit,T}} <: Forcing{transformed_unit,T}
- orig::TO
- func::TF
- name::Symbol
- function TransformedForcing(transformed_unit::Unitful.Units, orig::TO, func::TF, name::Symbol) where {unit,T,TO<:Forcing{unit,T},TF}
- new{transformed_unit,unit,TF,T,TO}(orig, func, name)
- end
-end
-
-function TransformedForcing(orig::TO, func::TF, name::Symbol) where {unit,T,TO<:Forcing{unit,T},TF}
- val_with_units = one(T)*unit
- transformed_val = func(val_with_units)
- return TransformedForcing(Unitful.unit(transformed_val), orig, func, name)
-end
-
-(f::TransformedForcing)(x::Number) = f.func(f.orig(x))
-
-"""
- InterpolatedForcing{unit,T,TI}
-
-Forcing data provided by a discrete time series of data.
-"""
-struct InterpolatedForcing{unit,T,TI<:Interpolations.AbstractInterpolation{T}} <: Forcing{unit,T}
- interpolant::TI
- name::Symbol
- InterpolatedForcing(unit::Unitful.Units, interpolant::TI, name::Symbol) where {TI} = new{unit,eltype(interpolant),TI}(interpolant, name)
-end
-function InterpolatedForcing(timestamps::AbstractArray{DateTime,1}, values::A, name::Symbol; interpolation_mode=Interpolations.Linear()) where {T,A<:AbstractArray{T,1}}
- ts = convert_t.(timestamps)
- values_converted = Utils.normalize_units.(values)
- interp_values = Numerics.fpzero.(ustrip.(values_converted))
- interp = Interpolations.interpolate((ts,), interp_values, Interpolations.Gridded(interpolation_mode))
- return InterpolatedForcing(unit(eltype(values_converted)), interp, name)
-end
-Flatten.flattenable(::Type{<:InterpolatedForcing}, ::Type) = false
-Base.getindex(f::InterpolatedForcing, i::Integer) = f.interpolant.coefs[i]
-Base.getindex(f::InterpolatedForcing, t) = f(t)
-Base.length(f::InterpolatedForcing) = length(f.interpolant)
-Base.show(io::IO, forcing::InterpolatedForcing{u}) where u = print(io, "InterpolatedForcing $(forcing.name) [$u] of length $(length(forcing)) with time span $(convert_tspan(extrema(forcing.interpolant.knots[1])))")
-# Base.show(io::IO, ::Type{<:InterpolatedForcing}) = print(io, "InterpolatedForcing")
-# Base.show(io::IO, ::Type{<:InterpolatedForcing{unit}}) where {unit} = print(io, "InterpolatedForcing{$unit}")
-# Base.show(io::IO, ::Type{<:InterpolatedForcing{unit,T}}) where {unit,T} = print(io, "InterpolatedForcing{$unit, $T, ...}") # type piracy to truncate type name in stacktraces
-"""
-Get interpolated forcing value at t seconds from t0.
-"""
-@propagate_inbounds (forcing::InterpolatedForcing)(t::Number) = Numerics.fpzero(forcing.interpolant(t))
-
-"""
- time_derivative_forcing(
- f::InterpolatedForcing{unit},
- new_name::Symbol;
- interp=Numerics.Linear()
- ) where {unit}
-
-Computes the finite difference time derivative of the given `InterpolatedForcing`
-time series and returns a new forcing with units `[unit].sâ»Â¹`
-"""
-function time_derivative_forcing(
- f::InterpolatedForcing{unit},
- new_name::Symbol=Symbol(:∂,f.name,:∂t);
- interpolation_mode=Numerics.Constant()
-) where {unit}
- ts = f.interpolant.knots[1]
- ∂f∂t = map(t -> first(Numerics.gradient(f.interpolant, t))*unit/1.0u"s", ts)
- return InterpolatedForcing(convert_t.(ts), ∂f∂t, new_name; interpolation_mode)
-end
-
-"""
- Forcings{names,TF,TMeta}
-
-Generic container for forcing types with optional metadata.
-"""
-struct Forcings{names,TF,TMeta}
- data::NamedTuple{names,TF}
- metadata::TMeta
- Forcings(data::NamedTuple{names,TF}, metadata::NamedTuple=(;)) where {names,TF} = new{names,TF,typeof(metadata)}(data, metadata)
-end
-Flatten.flattenable(::Type{<:Forcings}, ::Val{:metadata}) = false
-CryoGrid.parameterize(f::Forcings; ignored...) = f
-Base.propertynames(::Forcings{names}) where {names} = (:metadata, names...)
-Base.getproperty(fs::Forcings{names}, sym::Symbol) where {names} = sym ∈ names ? getproperty(getfield(fs, :data), sym) : getfield(fs, sym)
-Base.merge(fs::Forcings...) = Forcings(merge(map(f -> getfield(f, :data), fs)...), merge(map(f -> getfield(f, :metadata), fs)...))
-Base.rename(fs::Forcings, names::Pair{Symbol,Symbol}) = Forcings((; map(nm -> nm == names[1] ? names[2] => fs.data[nm] : nm => fs.data[nm], keys(fs.data))...), fs.metadata)
-
forcingunits(::ForcingFormat) = Dict()
detectformat(::Val{x}, filepath) where x = error("unrecognized forcing file suffix $x")
@@ -159,12 +14,12 @@ function detectformat(filepath::String)
end
"""
- loadforcings(filename::String)::Forcings
- loadforcings(resource::Resource; outdir=DEFAULT_FORCINGS_DIR)::Forcings
- loadforcings([format::ForcingFormat], filename::String; outdir=DEFAULT_FORCINGS_DIR)::Forcings
+ loadforcings(filename::String)
+ loadforcings(resource::Resource; outdir=DEFAULT_FORCINGS_DIR)
+ loadforcings([format::ForcingFormat], filename::String; outdir=DEFAULT_FORCINGS_DIR)
Loads forcing data from the given file according to the format specified by `format`. By default, the forcing format
-is automatically detected via `detectformat`. Returns a `Forcings` struct containing all forcing data
+is automatically detected via `detectformat`. Returns a `DimStack` containing all forcing data
and metadata
"""
loadforcings(filename::String) = loadforcings(detectformat(filename), filename)
@@ -174,5 +29,6 @@ loadforcings(f::ForcingFormat, filename::String) = error("loadforcings not imple
_normalize_numeric(x::Number) = convert(Float64, x)
_normalize_numeric(::Union{Missing,Nothing}) = missing
+include("providers.jl")
include("forcings_json.jl")
include("forcings_ncd.jl")
diff --git a/src/IO/forcings/forcings_json.jl b/src/IO/forcings/forcings_json.jl
index ed60224b..3c56b0d8 100644
--- a/src/IO/forcings/forcings_json.jl
+++ b/src/IO/forcings/forcings_json.jl
@@ -18,7 +18,7 @@ forcingunits(::ForcingFormatJSON) = Dict(
:Ptot => u"mm/d",
)
-function loadforcings(format::ForcingFormatJSON{1}, filename::String)
+function loadforcings(format::ForcingFormatJSON{1}, filename::String, kwargs...)
dict = open(filename, "r") do file; JSON3.read(file) end
# convert JSON3 dict for data field to Julia dict
data = Dict(dict[:data]...)
@@ -37,9 +37,9 @@ function loadforcings(format::ForcingFormatJSON{1}, filename::String)
end
end
forcings = map(names, vals_with_units) do name, values
- name => InterpolatedForcing(Array(ts), values, name)
+ name => DimArray(values, (Ti(Array(ts)),); name)
end
- return Forcings((; forcings...))
+ return Interpolated1D(DimStack((; forcings...)); kwargs...)
end
function loadforcings(format::ForcingFormatJSON{2}, filename::String)
dict = open(filename, "r") do file; JSON3.read(file) end
@@ -58,7 +58,7 @@ function loadforcings(format::ForcingFormatJSON{2}, filename::String)
end
end
forcings = map(names, vals_with_units) do name, values
- name => InterpolatedForcing(Array(ts), values, name)
+ name => DimArray(values, (Ti(Array(ts)),); name)
end
- return Forcings((; forcings...))
+ return Interpolated1D(DimStack((; forcings...)); kwargs...)
end
diff --git a/src/IO/forcings/forcings_ncd.jl b/src/IO/forcings/forcings_ncd.jl
index d17509e9..2f3868fe 100644
--- a/src/IO/forcings/forcings_ncd.jl
+++ b/src/IO/forcings/forcings_ncd.jl
@@ -73,10 +73,10 @@ function loadforcings(::ForcingFormatNCD{NCDTopoPyScale}, filepath::String)
# note that this is loading the data entirely into memory;
# may need to consider in the future lazy loading via DiskArray and/or the NetCDF package.
# NetCDF.jl currently seems to have a bug where the data types are not inferred correctly.
- Symbol(name) => InterpolatedForcing(timestamps, var.*units, Symbol(standard_name))
+ Symbol(name) => DimArray(var.*units, (Ti(timestamps),); name=Symbol(standard_name))
end
metadata = (; title, source, creator_name, date_created, latitude, longitude, elevation, calendar)
# note the (;x...) syntax creates a named tuple from a collection of pairs, Symbol => value
- return Forcings((;forcings...), metadata)
+ return DimStack((;forcings...); metadata)
end
end
diff --git a/src/IO/forcings/providers.jl b/src/IO/forcings/providers.jl
new file mode 100644
index 00000000..dfbaf81b
--- /dev/null
+++ b/src/IO/forcings/providers.jl
@@ -0,0 +1,53 @@
+"""
+ ForcingProvider
+
+Base type for forcing "providers", i.e. functional interfaces for
+discrete forcing datasets.
+"""
+abstract type ForcingProvider <: InputProvider end
+
+struct Interpolated1D{names,TVals} <: ForcingProvider
+ data::DimStack
+ interpolants::NamedTuple{names,TVals}
+end
+
+function Interpolated1D(
+ forcingdata::DimStack,
+ default_interp_mode=Linear();
+ stripunits=true,
+ interp_modes...)
+ ts = convert_t.(dims(forcingdata, Ti))
+ interpolants = map(forcingdata) do data
+ @assert length(size(data)) == 1 "Interpolated1D only supports forcing data with a single dimension"
+ data = stripunits ? ustrip.(data) : data
+ interp = get(interp_modes, data.name, default_interp_mode)
+ interpolate((ts,), data, Gridded(interp))
+ end
+ return Interpolated1D(forcingdata, interpolants)
+end
+
+(data::Interpolated1D)(input::Input, t::DateTime) = data(input, convert_t(t))
+function (data::Interpolated1D)(::Input{name}, t::Number) where {name}
+ itp = getproperty(data, name)
+ return itp(t)
+end
+
+# special overrides for prroperty names and getproperty
+Base.propertynames(::Interpolated1D{names}) = (:data, names...)
+Base.getproperty(data::Interpolated1D, name::Symbol) = name == :data ? getfield(data, :data) : getproperty(inputmap(data), name)
+
+function Base.show(io::IO, mimetype::MIME"text/plain", data::Interpolated1D)
+ println("Interpolated1D of")
+ show(io, mimetype, getfield(data, :data))
+end
+
+function inputmap(data::Interpolated1D)
+ return map(_wrap_interpolant, getfield(data, :interpolants))
+end
+
+# wraps the interpolant to support both numeric and DateTime types.
+function _wrap_interpolant(itp::Interpolations.AbstractInterpolation)
+ interpolant(t::Number) = itp(t)
+ interpolant(t::DateTime) = itp(convert_t(t))
+ return interpolant
+end
diff --git a/src/IO/input.jl b/src/IO/input.jl
new file mode 100644
index 00000000..4fbf632f
--- /dev/null
+++ b/src/IO/input.jl
@@ -0,0 +1,85 @@
+"""
+ Input{name,attrsType<:NamedTuple}
+
+The `Input` type represents a placeholder for (typically time-varying)
+inputs that will be replaced with the numeric value at runtime. The
+`name` of the `Input` will be used to match it to the corresponding
+input function.
+"""
+struct Input{name,attrsType<:NamedTuple}
+ attrs::attrsType
+ Input(name::Symbol, attrs::NamedTuple) = new{name,typeof(attrs)}(attrs)
+ """
+ Input(name::Symbol; kwargs...)
+
+ Constructs a `Input` with the given identifier `name` and optional
+ attribtues given as keyword arguments.
+ """
+ function Input(name::Symbol; kwargs...)
+ attrs = (; kwargs...)
+ new{name,typeof(attrs)}(attrs)
+ end
+end
+
+ConstructionBase.constructorof(::Type{<:Input{name}}) where {name} = attrs -> Input(name, attrs)
+
+"""
+ nameof(::Input{name}) where {name}
+
+Extract and return the `name` of the given `Input`.
+"""
+Base.nameof(::Input{name}) where {name} = name
+
+function Base.getproperty(f::Input, name::Symbol)
+ attrs = getfield(f, :attrs)
+ if name ∈ keys(attrs)
+ return getproperty(attrs, name)
+ else
+ return getfield(f, name)
+ end
+end
+
+"""
+ InputProvider
+
+Base type for `Input` providers that realize one or more `Input`s.
+"""
+abstract type InputProvider end
+
+"""
+ inputmap(provider::InputProvider)
+
+Returns a mapping-like object (e.g. a named tuple) with key/value pairs
+corresponding to names and realized inputs.
+"""
+inputmap(::InputProvider) = error("not implemented")
+
+Base.propertynames(provider::InputProvider) = propertynames(inputmap(provider))
+Base.getproperty(provider::InputProvider, name::Symbol) = getproperty(inputmap(provider), name)
+Base.keys(provider::InputProvider) = propertynames(provider)
+Base.values(provider::InputProvider) = values(inputmap(provider))
+
+"""
+ InputFunctionProvider
+
+Generic input provider that simply wraps a `NamedTuple` of functions.
+"""
+struct InputFunctionProvider{TF<:NamedTuple} <: InputProvider
+ inputs::TF
+end
+
+InputFunctionProvider(; input_functions...) = InputFunctionProvider((; input_functions...))
+
+function (ifp::InputFunctionProvider)(::Input{name}, args...)
+ f = ifp.inputs[name]
+ return f(args...)
+end
+
+inputmap(ifp::InputFunctionProvider) = getfield(ifp, :inputs)
+
+"""
+ inputs(; named_inputs...)
+
+Alias for the constructor of `InputFunctionProvider`.
+"""
+inputs(; named_inputs...) = InputFunctionProvider(; named_inputs...)
diff --git a/src/Physics/Heat/heat_bc.jl b/src/Physics/Heat/heat_bc.jl
index a09ec578..0a358042 100644
--- a/src/Physics/Heat/heat_bc.jl
+++ b/src/Physics/Heat/heat_bc.jl
@@ -37,29 +37,37 @@ end
"""
TemperatureBC{E,F} <: BoundaryProcess{HeatBalance}
-Represents a simple, forced Dirichlet temperature boundary condition for `HeatBalance` processes.
+Represents a simple, Dirichlet temperature boundary condition for `HeatBalance` processes.
"""
struct TemperatureBC{E,F} <: BoundaryProcess{HeatBalance}
- T::F # temperature forcing
+ T::F # boundary temperature
effect::E # effect
- TemperatureBC(T::F, effect::E=nothing) where {F<:Forcing{u"°C"},E} = new{E,F}(T, effect)
+ TemperatureBC(T::F, effect::E=nothing) where {F,E} = new{E,F}(T, effect)
end
CryoGrid.BCKind(::Type{<:TemperatureBC}) = Dirichlet()
-CryoGrid.boundaryvalue(bc::TemperatureBC, state) = getscalar(state.T_ub)
+CryoGrid.boundaryvalue(::TemperatureBC, state) = getscalar(state.T_ub)
-CryoGrid.variables(::Union{Top,Bottom}, bc::TemperatureBC) = (
+CryoGrid.variables(::Top, ::TemperatureBC) = (
Diagnostic(:T_ub, Scalar, u"K"),
)
-function CryoGrid.computediagnostic!(::Union{Top,Bottom}, bc::TemperatureBC, state)
- @setscalar state.T_ub = bc.T(state.t)
+CryoGrid.variables(::Bottom, ::TemperatureBC) = (
+ Diagnostic(:T_lb, Scalar, u"K"),
+)
+
+function CryoGrid.computediagnostic!(::Top, bc::TemperatureBC, state)
+ @setscalar state.T_ub = bc.T
+end
+
+function CryoGrid.computediagnostic!(::Bottom, bc::TemperatureBC, state)
+ @setscalar state.T_lb = bc.T
end
Base.@kwdef struct NFactor{W,S} <: CryoGrid.BoundaryEffect
- nf::W = 1.0 # applied when Tair <= 0
- nt::S = 1.0 # applied when Tair > 0
+ nf::W = Param(1.0, domain=0..1) # applied when Tair <= 0
+ nt::S = Param(1.0, domain=0..1) # applied when Tair > 0
end
CryoGrid.variables(::Top, bc::TemperatureBC{<:NFactor}) = (
@@ -67,17 +75,12 @@ CryoGrid.variables(::Top, bc::TemperatureBC{<:NFactor}) = (
Diagnostic(:nfactor, Scalar),
)
-CryoGrid.parameterize(nf::NFactor) = NFactor(
- nf = CryoGrid.parameterize(nf.nf, domain=0..1),
- nt = CryoGrid.parameterize(nf.nt, domain=0..1),
-)
-
nfactor(Tair, nfw, nfs) = (Tair <= zero(Tair))*nfw + (Tair > zero(Tair))*nfs
function CryoGrid.computediagnostic!(::Top, bc::TemperatureBC{<:NFactor}, state)
nfw = bc.effect.nf
nfs = bc.effect.nt
- Tair = bc.T(state.t)
+ Tair = bc.T
@setscalar state.nfactor = nfactor(Tair, nfw, nfs)
@setscalar state.T_ub = getscalar(state.nfactor)*Tair
end
@@ -90,10 +93,10 @@ Represents a simple, forced Neumann heat flux boundary condition for `HeatBalanc
struct GroundHeatFlux{TE,TQ} <: BoundaryProcess{HeatBalance}
Qg::TQ
effect::TE
- GroundHeatFlux(Qg::TQ, effect::TE=nothing) where {TQ<:Forcing{u"W/m^2"},TE} = new{TE,TQ}(Qg, effect)
+ GroundHeatFlux(Qg::TQ, effect::TE=nothing) where {TQ,TE} = new{TE,TQ}(Qg, effect)
end
-CryoGrid.boundaryvalue(bc::GroundHeatFlux, state) = bc.Qg(state.t)
+CryoGrid.boundaryvalue(bc::GroundHeatFlux, state) = bc.Qg
CryoGrid.BCKind(::Type{<:GroundHeatFlux}) = CryoGrid.Neumann()
@@ -109,6 +112,4 @@ end
CryoGrid.boundaryvalue(bc::GeothermalHeatFlux, state) = -bc.Qgeo
-CryoGrid.boundaryvalue(bc::GeothermalHeatFlux{<:Forcing}, state) = -bc.Qgeo(state.t)
-
CryoGrid.BCKind(::Type{<:GeothermalHeatFlux}) = CryoGrid.Neumann()
diff --git a/src/Physics/Heat/heat_conduction.jl b/src/Physics/Heat/heat_conduction.jl
index 785230f7..df28efd5 100644
--- a/src/Physics/Heat/heat_conduction.jl
+++ b/src/Physics/Heat/heat_conduction.jl
@@ -29,6 +29,13 @@ function freezethaw!(fw::FreeWater, sub::SubSurface, heat::HeatBalance{<:Enthalp
return nothing
end
+"""
+ enthalpyinv(::FreeWater, H, θwi, C, L)
+
+Inverse enthalpy function for free water freezing characteristic given
+enthalpy `H`, total water content `θwi`, heat capacity `C`, and latent
+heat of fusion `L`.
+"""
function enthalpyinv(::FreeWater, H, θwi, C, L)
Lθ = L*θwi
T_f = H / C
@@ -49,6 +56,33 @@ function enthalpyinv(::FreeWater, H, θwi, C, L)
return T
end
+"""
+ heatflux!(dH, jH, T, k, ΔT, Δk)
+
+Compute heat fluxes over internal grid cells and divergence for all grid cells.
+Note that fluxes are added to existing values in `dH` and `jH`.
+"""
+function heatflux!(dH, jH, T, k, ΔT, Δk)
+ # compute internal fluxes and non-linear diffusion assuming boundary fluxes have been set
+ # note that we now use the (slightly slower) explicit flux! and divergence! formulation since
+ # the nonlineardiffusion! function is not compatible with additive (existing) fluxes.
+ flux!(jH, T, ΔT, k)
+ divergence!(dH, jH, Δk)
+ return nothing
+end
+
+"""
+ temperatureflux!(dT, dH, ∂H∂T)
+
+Compute the temperature flux `dT = dH / ∂H∂T` as a funtion of the heat flux and effective heat capacity.
+"""
+function temperatureflux!(dT, dH, ∂H∂T)
+ # Compute temperature flux by dividing by ∂H∂T;
+ # ∂H∂T should be computed by the freeze curve.
+ @inbounds @. dT = dH / ∂H∂T
+ return nothing
+end
+
"""
Common variable definitions for all heat implementations.
"""
@@ -107,26 +141,14 @@ end
function CryoGrid.computeprognostic!(::SubSurface, ::HeatBalance{<:EnthalpyBased}, state)
Δk = Δ(state.grid) # cell sizes
ΔT = Δ(cells(state.grid)) # midpoint distances
- # compute internal fluxes and non-linear diffusion assuming boundary fluxes have been set;
- # note that we now use the (slightly slower) explicit flux! and divergence! formulation since
- # the nonlineardiffusion! function is not compatible with additive (existing) fluxes.
- # nonlineardiffusion!(state.dH, state.jH, state.T, ΔT, state.k, Δk)
- flux!(state.jH, state.T, ΔT, state.k)
- divergence!(state.dH, state.jH, Δk)
+ heatflux!(state.dH, state.jH, state.T, state.k, ΔT, Δk)
return nothing
end
function CryoGrid.computeprognostic!(sub::SubSurface, ::HeatBalance{<:TemperatureBased}, state)
Δk = Δ(state.grid) # cell sizes
ΔT = Δ(cells(state.grid)) # midpoint distances
- # compute internal fluxes and non-linear diffusion assuming boundary fluxes have been set
- # note that we now use the (slightly slower) explicit flux! and divergence! formulation since
- # the nonlineardiffusion! function is not compatible with additive (existing) fluxes.
- # nonlineardiffusion!(state.dH, state.jH, state.T, ΔT, state.k, Δk)
- flux!(state.jH, state.T, ΔT, state.k)
- divergence!(state.dH, state.jH, Δk)
- # Compute temperature flux by dividing by ∂H∂T;
- # ∂H∂T should be computed by the freeze curve.
- @inbounds @. state.dT = state.dH / state.∂H∂T
+ heatflux!(state.dH, state.jH, state.T, state.k, ΔT, Δk)
+ temperatureflux!(state.dT, state.dH, state.∂H∂T)
return nothing
end
diff --git a/src/Physics/Heat/heat_implicit.jl b/src/Physics/Heat/heat_implicit.jl
index 254527d8..ad5ceeed 100644
--- a/src/Physics/Heat/heat_implicit.jl
+++ b/src/Physics/Heat/heat_implicit.jl
@@ -3,6 +3,26 @@ Type alias for the implicit enthalpy formulation of HeatBalance.
"""
const HeatBalanceImplicit = HeatBalance{<:EnthalpyImplicit}
+apbc(::Dirichlet, k, Δk) = k / (Δk^2/2)
+apbc(::Dirichlet, k, Δk, Δx) = k / Δx / Δk
+apbc(::Neumann, k, Δk) = 0
+
+function prefactors!(ap, an, as, k, dx, dxp)
+ # loop over grid cells
+ @inbounds for i in eachindex(dxp)
+ if i == 1
+ as[1] = k[2] / dx[1] / dxp[1]
+ elseif i == length(dxp)
+ an[end] = k[end-1] / dx[end] / dxp[end]
+ else
+ an[i] = k[i] / dx[i-1] / dxp[i]
+ as[i] = k[i+1] / dx[i] / dxp[i]
+ end
+ ap[i] = an[i] + as[i]
+ end
+ return nothing
+end
+
# CryoGrid methods
CryoGrid.variables(heat::HeatBalanceImplicit) = (
@@ -95,25 +115,3 @@ function CryoGrid.resetfluxes!(sub::SubSurface, heat::HeatBalanceImplicit, state
state.DT_as[i] = 0.0
end
end
-
-# implementations
-
-apbc(::Dirichlet, k, Δk) = k / (Δk^2/2)
-apbc(::Dirichlet, k, Δk, Δx) = k / Δx / Δk
-apbc(::Neumann, k, Δk) = 0
-
-function prefactors!(ap, an, as, k, dx, dxp)
- # loop over grid cells
- @inbounds for i in eachindex(dxp)
- if i == 1
- as[1] = k[2] / dx[1] / dxp[1]
- elseif i == length(dxp)
- an[end] = k[end-1] / dx[end] / dxp[end]
- else
- an[i] = k[i] / dx[i-1] / dxp[i]
- as[i] = k[i+1] / dx[i] / dxp[i]
- end
- ap[i] = an[i] + as[i]
- end
- return nothing
-end
diff --git a/src/Physics/Heat/heat_methods.jl b/src/Physics/Heat/heat_methods.jl
index e3447bc4..7f198c5e 100644
--- a/src/Physics/Heat/heat_methods.jl
+++ b/src/Physics/Heat/heat_methods.jl
@@ -110,3 +110,4 @@ dHdT(T, C, L, ∂θw∂T, ch_w, ch_i) = C + ∂θw∂T*(L + T*(ch_w - ch_i))
Convenience constructor for `Numerics.Profile` which automatically converts temperature quantities.
"""
TemperatureProfile(pairs::Pair{<:Union{DistQuantity,Param},<:Union{TempQuantity,Param}}...) = Profile(map(p -> uconvert(u"m", p[1]) => uconvert(u"°C", p[2]), pairs))
+TemperatureProfile(pairs::Pair{<:Union{DistQuantity,Param}}...) = Profile(map(p -> uconvert(u"m", p[1]) => p[2], pairs))
diff --git a/src/Physics/Heat/heat_types.jl b/src/Physics/Heat/heat_types.jl
index d4772cc5..66bce2bc 100644
--- a/src/Physics/Heat/heat_types.jl
+++ b/src/Physics/Heat/heat_types.jl
@@ -67,5 +67,3 @@ Utils.@properties HeatBalanceProperties(
Lsl = CryoGrid.Constants.Lsl,
L = Ïw*Lsl,
)
-# do not parameterize heat properties
-CryoGrid.parameterize(prop::HeatBalanceProperties) = prop
diff --git a/src/Physics/Heat/heat_water.jl b/src/Physics/Heat/heat_water.jl
index 8ea58782..18a403e5 100644
--- a/src/Physics/Heat/heat_water.jl
+++ b/src/Physics/Heat/heat_water.jl
@@ -13,18 +13,16 @@ given the heat capacity of water `cw` and latent heat of fusion `L`.
advectiveflux(jw, Tâ‚, Tâ‚‚, cw, L) = jw*(cw*Tâ‚*(jw > zero(jw)) + cw*Tâ‚‚*(jw < zero(jw)) + L)
"""
- water_energy_advection!(::SubSurface, ::Coupled(WaterBalance, HeatBalance), state)
+ water_energy_advection!(jH, jw, T, L::Real)
Adds advective energy fluxes for all internal grid cell faces.
"""
-function water_energy_advection!(sub::SubSurface, ps::Coupled(WaterBalance, HeatBalance), state)
+function water_energy_advection!(jH, jw, T, cw::Real, L::Real)
water, heat = ps
- cw = heatcapacitywater(sub, state)
- @inbounds for i in 2:length(state.jw)-1
- let jw = state.jw[i],
- Tâ‚ = state.T[i-1],
- Tâ‚‚ = state.T[i],
- L = heat.prop.L;
+ @inbounds for i in 2:length(jw)-1
+ let jw = jw[i],
+ Tâ‚ = T[i-1],
+ Tâ‚‚ = T[i];
jH_w = advectiveflux(jw, Tâ‚, Tâ‚‚, cw, L)
state.jH[i] += jH_w
end
@@ -91,9 +89,11 @@ end
# Flux calculation
function CryoGrid.computeprognostic!(sub::SubSurface, ps::Coupled(WaterBalance, HeatBalance), state)
water, heat = ps
+ cw = heatcapacitywater(sub, state)
+ L = heat.prop.L
CryoGrid.computeprognostic!(sub, water, state)
if heat.advection
- water_energy_advection!(sub, ps, state)
+ water_energy_advection!(state.jH, state.jw, state.T, cw, L)
end
CryoGrid.computeprognostic!(sub, heat, state)
end
diff --git a/src/Physics/Hydrology/water_ET.jl b/src/Physics/Hydrology/water_ET.jl
index 2f5e5e17..fe6c3e7e 100644
--- a/src/Physics/Hydrology/water_ET.jl
+++ b/src/Physics/Hydrology/water_ET.jl
@@ -11,15 +11,11 @@ struct EvapTop <: Evapotranspiration end
Corresponds to evapotranspiration scheme 2 described in section 2.2.4 of Westermann et al. (2022).
"""
Base.@kwdef struct DampedET{Tftr,Tdtr,Tdev} <: Evapotranspiration
- f_tr::Tftr = 0.5
- d_tr::Tdtr = 0.5u"m"
- d_ev::Tdev = 0.1u"m"
+ f_tr::Tftr = Param(0.5, domain=0..1, desc="Factor between 0 and 1 weighting transpirative vs. evaporative fluxes.")
+ d_tr::Tdtr = Param(0.5, units=u"m", domain=0..Inf, desc="Damping depth for transpiration.")
+ d_ev::Tdev = Param(0.1, units=u"m", domain=0..Inf, desc="Damping depth for evaporation.")
end
-CryoGrid.parameterize(et::DampedET) = DampedET(
- f_tr = CryoGrid.parameterize(et.f_tr, domain=0..1, desc="Factor between 0 and 1 weighting transpirative vs. evaporative fluxes."),
- d_tr = CryoGrid.parameterize(et.d_tr, domain=0..Inf, desc="Damping depth for transpiration."),
- d_ev = CryoGrid.parameterize(et.d_ev, domain=0..Inf, desc="Damping depth for evaporation."),
-)
+
"""
evapotranspiration!(::SubSurface, ::WaterBalance, state)
diff --git a/src/Physics/Hydrology/water_types.jl b/src/Physics/Hydrology/water_types.jl
index b9cb5666..f5b2c72d 100644
--- a/src/Physics/Hydrology/water_types.jl
+++ b/src/Physics/Hydrology/water_types.jl
@@ -9,27 +9,15 @@ Base.@kwdef struct WaterBalanceProperties{TÏw,TLsg,Trfs}
rf_smoothness::Trfs = 0.3
end
-# do not parameterize water balance constants
-CryoGrid.parameterize(prop::WaterBalanceProperties) = prop
-
"""
HydraulicProperties
Default material hydraulic properties.
"""
Utils.@properties HydraulicProperties(
- kw_sat = 1e-5u"m/s",
+ kw_sat = Param(1e-5, units=u"m/s"),
)
-function CryoGrid.parameterize(prop::HydraulicProperties)
- return HydraulicProperties(
- map(values(prop)) do val
- # this currently assumes that all properties have a strictly positive domain!
- CryoGrid.parameterize(val, domain=StrictlyPositive)
- end
- )
-end
-
"""
WaterFlow
diff --git a/src/Physics/Salt/salt_bc.jl b/src/Physics/Salt/salt_bc.jl
index 87f4efef..cef31f8f 100644
--- a/src/Physics/Salt/salt_bc.jl
+++ b/src/Physics/Salt/salt_bc.jl
@@ -9,18 +9,16 @@ Base.@kwdef struct SaltGradient{TbenthicSalt, TsurfaceState} <: BoundaryProcess{
surfaceState::TsurfaceState
end
-benthicSalt(bc::SaltGradient, t) = bc.benthicSalt
-benthicSalt(bc::SaltGradient{<:Forcing}, t) = bc.benthicSalt(t)
+benthicSalt(bc::SaltGradient) = bc.benthicSalt
-surfaceState(bc::SaltGradient, t) = bc.surfaceState
-surfaceState(bc::SaltGradient{<:Forcing}, t) = bc.surfaceState(t)
+surfaceState(bc::SaltGradient) = bc.surfaceState
CryoGrid.BCKind(::Type{<:SaltGradient}) = CryoGrid.Dirichlet()
function CryoGrid.interact!(::Top, bc::SaltGradient, soil::SalineGround, ::SaltMassBalance, stop, ssed)
# upper boundary
- surfaceState_t = surfaceState(bc, stop.t)
- benthicSalt_t = benthicSalt(bc, stop.t)
+ surfaceState_t = surfaceState(bc)
+ benthicSalt_t = benthicSalt(bc)
# distance to boundary; 1/2 thickness of first grid cell
dz = CryoGrid.thickness(soil, ssed, first) / 2
flux = if (surfaceState_t == 0) # if is inundated
diff --git a/src/Physics/Soils/para/simple.jl b/src/Physics/Soils/para/simple.jl
index 2052ae5d..213b902b 100644
--- a/src/Physics/Soils/para/simple.jl
+++ b/src/Physics/Soils/para/simple.jl
@@ -6,9 +6,9 @@ i.e. natural porosity, saturation, and organic solid fraction. This is the stand
of a discrete soil volume.
"""
Base.@kwdef struct SimpleSoil{Tfc,Tpor,Tsat,Torg,Thp,Twp} <: SoilParameterization
- por::Tpor = 0.5 # natural porosity
- sat::Tsat = 1.0 # saturation
- org::Torg = 0.0 # organic fraction of solid; mineral fraction is 1-org
+ por::Tpor = Param(0.5, domain=0..1, desc="Natural porosity of the soil volume.")
+ sat::Tsat = Param(1.0, domain=0..1, desc="Initial water+ice saturation level of the soil volume.")
+ org::Torg = Param(0.0, domain=0..1, desc="Organic solid fraction of the soil volume.")
freezecurve::Tfc = FreeWater()
heat::Thp = SoilThermalProperties(SimpleSoil)
water::Twp = SoilHydraulicProperties(SimpleSoil, fieldcapacity=0.20)
@@ -34,25 +34,16 @@ SoilThermalProperties(
kh_w = ThermalProperties().kh_w,
kh_i = ThermalProperties().kh_i,
kh_a = ThermalProperties().kh_a,
- kh_o=0.25u"W/m/K", # organic [Hillel (1982)]
- kh_m=3.8u"W/m/K", # mineral [Hillel (1982)]
+ kh_o=Param(0.25, units=u"W/m/K", domain=0..Inf), # organic [Hillel (1982)]
+ kh_m=Param(3.8, untis=u"W/m/K", domain=0..Inf), # mineral [Hillel (1982)]
ch_w = ThermalProperties().ch_w,
ch_i = ThermalProperties().ch_i,
ch_a = ThermalProperties().ch_a,
- ch_o=2.5e6u"J/K/m^3", # heat capacity organic
- ch_m=2.0e6u"J/K/m^3", # heat capacity mineral
+ ch_o=Param(2.5e6, units=u"J/K/m^3", domain=0..Inf), # heat capacity organic
+ ch_m=Param(2.0e6, untis=u"J/K/m^3", domain=0..Inf), # heat capacity mineral
kwargs...,
) = ThermalProperties(; kh_w, kh_i, kh_a, kh_m, kh_o, ch_w, ch_i, ch_a, ch_m, ch_o, kwargs...)
-# CryoGrid methods
-CryoGrid.parameterize(para::SimpleSoil) = SimpleSoil(
- por = CryoGrid.parameterize(para.por, domain=0..1, desc="Natural porosity of the soil volume."),
- sat = CryoGrid.parameterize(para.sat, domain=0..1, desc="Initial water+ice saturation level of the soil volume."),
- org = CryoGrid.parameterize(para.org, domain=0..1, desc="Organic solid fraction of the soil volume."),
- heat = CryoGrid.parameterize(para.heat, domain=0..Inf),
- water = CryoGrid.parameterize(para.water, domain=0..Inf),
-)
-
CryoGrid.variables(soil::Soil{<:SimpleSoil}) = CryoGrid.variables(soil, processes(soil))
CryoGrid.initializers(soil::Soil{<:SimpleSoil,THeat,<:WaterBalance}) where {THeat} = (
diff --git a/src/Physics/Soils/soil_heat.jl b/src/Physics/Soils/soil_heat.jl
index 36790af9..fd53ef4e 100644
--- a/src/Physics/Soils/soil_heat.jl
+++ b/src/Physics/Soils/soil_heat.jl
@@ -157,49 +157,3 @@ function Heat.enthalpyinv(sfcc::SFCC, soil::Soil, heat::HeatBalance{<:EnthalpyBa
return T_sol
end
end
-
-# Freeze curve parameters;
-# Since the freeze curve functions are specified in FreezeCurves.jl, we must (or rather should) provide
-# CryoGrid.parameterize implementations for them here to specify parameter information.
-CryoGrid.parameterize(f::PainterKarra) = PainterKarra(
- f.freezethaw,
- CryoGrid.parameterize(f.β, domain=OpenInterval(0,Inf), desc="Painter-Karra fitting parmaeter which controls the influence of temperature on the matric potential."),
- CryoGrid.parameterize(f.ω, domain=0..(1/f.β), desc="Painter-Karra fitting parameter which controls the depression of the melting point from saturation level."),
- f.g,
- CryoGrid.parameterize(f.swrc),
-)
-CryoGrid.parameterize(f::DallAmico) = DallAmico(
- f.freezethaw,
- f.g,
- CryoGrid.parameterize(f.swrc),
-)
-CryoGrid.parameterize(f::DallAmicoSalt) = DallAmicoSalt(
- f.freezethaw,
- CryoGrid.parameterize(f.saltconc, domain=Interval{:closed,:open}(0,Inf), desc="Assumed salt concentration when not defined as a state variable."), # salt concentration
- f.R,
- f.g,
- CryoGrid.parameterize(f.swrc),
-)
-CryoGrid.parameterize(f::McKenzie) = McKenzie(
- f.freezethaw,
- f.vol,
- CryoGrid.parameterize(f.γ, domain=OpenInterval(0,Inf)),
-)
-CryoGrid.parameterize(f::Westermann) = McKenzie(
- f.freezethaw,
- f.vol,
- CryoGrid.parameterize(f.δ, domain=OpenInterval(0,Inf)),
-)
-CryoGrid.parameterize(f::VanGenuchten) = VanGenuchten(
- f.vol,
- CryoGrid.parameterize(f.α, units=u"1/m", domain=OpenInterval(0,Inf), desc="van Genuchten α parameter which corresponds to the inverse of air entry suction."),
- CryoGrid.parameterize(f.n, domain=OpenInterval(1,Inf), desc="van Genuchten n parameter which controls the shape of the curve. Smaller values generally produce longer tailed curves."),
-)
-CryoGrid.parameterize(f::BrooksCorey) = BrooksCorey(
- f.vol,
- CryoGrid.parameterize(f.ψₛ, domain=OpenInterval(0,Inf), desc="Brooks-Corey suction parameter."),
- CryoGrid.parameterize(f.λ, domain=OpenInterval(0,Inf), desc="Brooks-Corey shape parameter."),
-)
-# do not parameterize default freeze curve properties or SFCC solvers
-CryoGrid.parameterize(prop::FreezeCurves.SoilFreezeThawProperties) = prop
-CryoGrid.parameterize(solver::FreezeCurves.SFCCSolver) = solver
diff --git a/src/Physics/Surface/SEB/seb.jl b/src/Physics/Surface/SEB/seb.jl
index 6bc234d5..4a98aa7b 100644
--- a/src/Physics/Surface/SEB/seb.jl
+++ b/src/Physics/Surface/SEB/seb.jl
@@ -67,36 +67,36 @@ Utils.@properties SEBParams(
βₕ = βₘ/Pr₀,
γₘ = 15.0,
γₕ = 9.0,
+
+ z = 2.0u"m" # height in meters of air and wind inputs
)
"""
- SurfaceEnergyBalance{TSolution,TStabFun,TPara,F} <: BoundaryProcess{HeatBalance}
+ SurfaceEnergyBalance{TSolution,TStabFun,TPara,TInputs} <: BoundaryProcess{HeatBalance}
Surface energy balance upper boundary condition.
"""
-struct SurfaceEnergyBalance{TSolution,TStabFun,TPara,F} <: BoundaryProcess{HeatBalance}
- forcings::F
+struct SurfaceEnergyBalance{TSolution<:SolutionScheme,TStabFun<:StabilityFunctions,TPara,TInputs} <: BoundaryProcess{HeatBalance}
+ inputs::TInputs
para::TPara
# type-dependent parameters
solscheme::TSolution
stabfun::TStabFun
end
-# User facing constructors
-SurfaceEnergyBalance(forcings::Forcings, z=-2.0u"m"; kwargs...) = SurfaceEnergyBalance(forcings.Tair, forcings.pressure, forcings.q, forcings.wind, forcings.Lin, forcings.Sin, z; kwargs...)
-function SurfaceEnergyBalance(
- Tair::TemperatureForcing, # air temperature
- pr::PressureForcing, # air pressure
- qh::HumidityForcing, # specific humidity
- wind::VelocityForcing, # non-directional wind speed
- Lin::EnergyFluxForcing, # long-wave incoming radiation
- Sin::EnergyFluxForcing, # short-wave incoming radiation
- z; # height [m] of air temperature and wind forcing
+# convenience constructor for SEB
+function SurfaceEnergyBalance(;
+ Tair=Input(:Tair),
+ Lin=Input(:Lin),
+ Sin=Input(:Sin),
+ pr=Input(:pr),
+ qh=Input(:qh),
+ wind=Input(:wind),
para::SEBParams = SEBParams(),
solscheme::SolutionScheme = Numerical(),
stabfun::StabilityFunctions = HøgstrømSHEBA(),
)
- forcings = (; Tair, pr, qh, wind, Lin, Sin, z);
- SurfaceEnergyBalance(forcings, para, solscheme, stabfun)
+ inputs = (; Tair, Lin, Sin, pr, qh, wind)
+ SurfaceEnergyBalance(; inputs, para, solscheme, stabfun)
end
"""
diff --git a/src/Physics/Surface/SEB/seb_heat.jl b/src/Physics/Surface/SEB/seb_heat.jl
index c3780d7b..e5af88be 100644
--- a/src/Physics/Surface/SEB/seb_heat.jl
+++ b/src/Physics/Surface/SEB/seb_heat.jl
@@ -14,7 +14,7 @@ function CryoGrid.initialcondition!(top::Top, seb::SurfaceEnergyBalance, state)
resetfluxes!(top, seb, state)
state.Lstar .= -1e5*oneunit(eltype(state.Lstar));
state.ustar .= 10.0*oneunit(eltype(state.ustar));
- state.T_ub .= seb.forcings.Tair(state.t)
+ state.T_ub .= seb.inputs.Tair
end
CryoGrid.BCKind(::Type{<:SurfaceEnergyBalance}) = CryoGrid.Neumann()
@@ -163,7 +163,7 @@ Friction velocity according to Monin-Obukhov theory
function u_star(seb::SurfaceEnergyBalance, state::SEBState)
let κ = seb.para.κ,
uz = state.inputs.wind, # wind speed at height z
- z = state.inputs.z, # height z of wind forcing
+ z = seb.para.z, # height z of wind forcing
zâ‚€ = state.inputs.surf.zâ‚€, # aerodynamic roughness length [m]
Lstar = state.Lstar;
κ * uz ./ (log(z / z₀) - Ψ_M(seb, z / Lstar, z₀ / Lstar)) # Eq. (7) in Westermann et al. (2016)
@@ -209,7 +209,7 @@ function L_star(seb::SurfaceEnergyBalance{Analytical}, state::SEBState)
p = state.inputs.pr, # atmospheric pressure at surface (height z)
pâ‚€ = state.inputs.pr, # normal pressure (for now assumed to be equal to p)
uz = state.inputs.wind, # wind speed at height z
- z = state.inputs.z, # height z of wind forcing
+ z = seb.para.z, # height z of wind forcing
zâ‚€ = state.inputs.surf.zâ‚€, # aerodynamic roughness length [m]
Prâ‚€ = seb.para.Prâ‚€, # turbulent Prandtl number
γₕ = seb.para.γₕ,
@@ -260,7 +260,7 @@ function Q_H(seb::SurfaceEnergyBalance, state::SEBState)
Tâ‚• = state.inputs.Tair, # air temperature
Tâ‚€ = state.inputs.Ts, # surface temperature
câ‚š = seb.para.câ‚ / seb.para.Ïâ‚, # specific heat capacity of air at constant pressure
- z = state.inputs.z, # height at which forcing data are provided
+ z = seb.para.z, # height at which forcing data are provided
Lstar = state.Lstar,
ustar = state.ustar,
pr = state.inputs.pr,
@@ -287,7 +287,7 @@ function Q_E(seb::SurfaceEnergyBalance, state::SEBState)
Tâ‚€ = state.inputs.Ts, # surface temperature
p = state.inputs.pr, # atmospheric pressure at surface
qâ‚• = state.inputs.qh, # specific humidity at height h over surface
- z = state.inputs.z, # height at which forcing data are provided
+ z = seb.para.z, # height at which forcing data are provided
râ‚› = state.inputs.surf.râ‚›, # surface resistance against evapotranspiration / sublimation [1/m]
Lstar = state.Lstar,
ustar = state.ustar,
diff --git a/src/Physics/Surface/SEB/seb_state.jl b/src/Physics/Surface/SEB/seb_state.jl
index 39cffc44..d2703a48 100644
--- a/src/Physics/Surface/SEB/seb_state.jl
+++ b/src/Physics/Surface/SEB/seb_state.jl
@@ -1,31 +1,30 @@
"""
SEBInputs{TT1,TT2,TR,TP,TQ,TW,TZ,Tsurf}
-Non-prognostic input variables to the SEB that are assumed to be known a priori.
+Non-prognostic input variables to the SEB that are assumed to be given.
"""
Base.@kwdef struct SEBInputs{TT1,TT2,TR,TP,TQ,TW,TZ,Tsurf}
- Ts::TT1
- Tair::TT2
- Lin::TR
- Sin::TR
- pr::TP
- qh::TQ
- wind::TW
- z::TZ
+ Ts::TT1 = nothing
+ Tair::TT2 = 0.0u"°C"
+ Lin::TR = 0.0u"W/m^2"
+ Sin::TR = 0.0u"W/m^2"
+ pr::TP = 0.0u"m/s"
+ qh::TQ = 0.0
+ wind::TW = 1e-8u"m/s"
# surface properties
surf::Tsurf = SurfaceProperties()
end
+
function SEBInputs(seb::SurfaceEnergyBalance, sub::SubSurface, stop, ssub)
Ts = ssub.T[1] # soil temperature at surface
- Tair = seb.forcings.Tair(stop.t)
- Lin = seb.forcings.Lin(stop.t)
- Sin = seb.forcings.Sin(stop.t)
- pr = seb.forcings.pr(stop.t)
- qh = seb.forcings.qh(stop.t)
- wind = seb.forcings.wind(stop.t)
- z = seb.forcings.z
+ Tair = seb.inputs.Tair
+ Lin = seb.inputs.Lin
+ Sin = seb.inputs.Sin
+ pr = seb.inputs.pr
+ qh = seb.inputs.qh
+ wind = seb.inputs.wind
surf = surfaceproperties(seb, sub)
- return SEBInputs(; Ts, Tair, Lin, Sin, pr, qh, wind, z, surf)
+ return SEBInputs(; Ts, Tair, Lin, Sin, pr, qh, wind, surf)
end
"""
diff --git a/src/Physics/Surface/SWB/swb.jl b/src/Physics/Surface/SWB/swb.jl
index 3003d3df..1c05ae3b 100644
--- a/src/Physics/Surface/SWB/swb.jl
+++ b/src/Physics/Surface/SWB/swb.jl
@@ -8,35 +8,8 @@ a Neumann-type upper boundary condition for snow and water fluxes as well as an
water mass balance at the surface.
"""
Base.@kwdef struct SurfaceWaterBalance{TR,TS} <: BoundaryProcess{Union{WaterBalance, SnowMassBalance}}
- rainfall::TR = ConstantForcing(0.0u"m/s")
- snowfall::TS = ConstantForcing(0.0u"m/s")
-end
-
-"""
- SurfaceWaterBalance(precip::Forcing{u"m/s",T1}, Tair::Forcing{u"°C",T2}) where {T1,T2}
-
-SWB constructor which derives rainfall and snowfall forcings from total precipitation and air temperature.
-"""
-function SurfaceWaterBalance(precip::Forcing{u"m/s",T1}, Tair::Forcing{u"°C",T2}) where {T1,T2}
- return SurfaceWaterBalance(
- TimeVaryingForcing(u"m/s", T1, t -> Tair(t) > 0.0 ? precip(t) : 0.0, :rainfall),
- TimeVaryingForcing(u"m/s", T1, t -> Tair(t) <= 0.0 ? precip(t) : 0.0, :snowfall),
- )
-end
-
-"""
- SurfaceWaterBalance(forcings::Forcings)
-
-Automatically attempts to extract precipitation forcings from `forcings`.
-"""
-function SurfaceWaterBalance(forcings::Forcings)
- if hasproperty(forcings, :rainfall) && hasproperty(forcings, :snowfall)
- return SurfaceWaterBalance(forcings.rainfall, forcings.snowfall)
- elseif hasproperty(forcings, :precip) && hasproperty(forcings, :Tair)
- return SurfaceWaterBalance(forcings.precip, forcings.Tair)
- else
- error("forcings must provide either 'rainfall' and 'snowfall' or 'precip' and 'Tair'")
- end
+ rainfall::TR = Input(:rainfall, units=u"m/s")
+ snowfall::TS = Input(:snowfall, units=u"m/s")
end
function infiltrate!(top::Top, swb::SurfaceWaterBalance, sub::SubSurface, water::WaterBalance, stop, ssub)
@@ -67,8 +40,8 @@ CryoGrid.variables(::Top, ::SurfaceWaterBalance) = (
)
function CryoGrid.computediagnostic!(::Top, swb::SurfaceWaterBalance, stop)
- @setscalar stop.jw_snow = swb.snowfall(stop.t)
- @setscalar stop.jw_rain = swb.rainfall(stop.t)
+ @setscalar stop.jw_snow = swb.snowfall
+ @setscalar stop.jw_rain = swb.rainfall
end
function CryoGrid.computeprognostic!(top::Top, swb::SurfaceWaterBalance, state)
diff --git a/src/Physics/steplimiters.jl b/src/Physics/steplimiters.jl
index 11bd26e3..eb033344 100644
--- a/src/Physics/steplimiters.jl
+++ b/src/Physics/steplimiters.jl
@@ -1,6 +1,5 @@
# Generic step limiter types
abstract type StepLimiter end
-CryoGrid.parameterize(limiter::StepLimiter) = limiter
"""
MaxDelta{T}
diff --git a/src/Presets/Presets.jl b/src/Presets/Presets.jl
index d027a53f..2978ec7e 100644
--- a/src/Presets/Presets.jl
+++ b/src/Presets/Presets.jl
@@ -17,33 +17,29 @@ using Statistics
export SoilHeatTile, SamoylovDefault
-include("presetgrids.jl")
-
"""
SoilHeatTile([heatop=:H], upperbc::BoundaryProcess, soilprofile::Profile, init::CryoGrid.Initializer...; grid::Grid=DefaultGrid_10cm, tile_kwargs...) where {F<:FreezeCurve}
Builds a simple one-layer soil/heat-conduction model with the given grid and configuration.
"""
-function SoilHeatTile(heatop, upperbc::BoundaryProcess, lowerbc::BoundaryProcess, soilprofile::Profile, inits::CryoGrid.Initializer...; grid::Grid=DefaultGrid_5cm, tile_kwargs...)
+function SoilHeatTile(
+ heatop,
+ upperbc::BoundaryProcess,
+ lowerbc::BoundaryProcess,
+ soilprofile::Profile,
+ inputs::InputProvider,
+ inits::CryoGrid.Initializer...;
+ grid::Grid=DefaultGrid_5cm,
+ tile_kwargs...
+)
strat = Stratigraphy(
grid[1] => Top(upperbc),
Tuple(d => Ground(para, heat=HeatBalance(heatop), water=nothing) for (i,(d,para)) in enumerate(soilprofile)),
grid[end] => Bottom(lowerbc)
)
- return Tile(strat, PresetGrid(grid), inits...; tile_kwargs...)
+ return Tile(strat, PresetGrid(grid), inputs, inits...; tile_kwargs...)
end
-SoilHeatTile(upperbc::BoundaryProcess, lowerbc::BoundaryProcess, soilprofile::Profile, inits::CryoGrid.Initializer...; kwargs...) = SoilHeatTile(:H, upperbc, lowerbc, soilprofile, inits...; kwargs...)
-
-Forcings = (
- Samoylov_ERA5_fitted_daily_1979_2020 = Resource("samoylov_era5_fitted_daily_1979-2020", ForcingFormatJSON{2}(), "https://nextcloud.awi.de/s/ScYAoHzeMzAfpjf/download/samoylov_era5_fitted_daily_1979-2020.json"),
- Samoylov_ERA_obs_fitted_1979_2014_spinup_extended_2044 = Resource("Samoylov_ERA_obs_fitted_1979_2014_spinup_extended_2044", ForcingFormatJSON{1}(), "https://nextcloud.awi.de/s/cbeycGQoQpXi3Ei/download/samoylov_ERA_obs_fitted_1979_2014_spinup_extended2044.json"),
- Samoylov_ERA_MkL3_CCSM4_long_term = Resource("Samoylov_ERA_MkL3_CCSM4_long_term", ForcingFormatJSON{1}(), "https://nextcloud.awi.de/s/45ax9AsTACxL25Q/download/FORCING_ULC_126_72.json"),
- Bayelva_ERA5_fitted_daily_1979_2020 = Resource("bayelva_era5_fitted_daily_1979-2020", ForcingFormatJSON{2}(), "https://nextcloud.awi.de/s/5AdbRMYKneCHgx4/download/bayelva_era5_fitted_daily_1979-2020.json")
-)
-Parameters = (
- # Faroux et al. doi:10.1109/IGARSS.2007.4422971
- EcoCLimMap_ULC_126_72 = Resource("EcoCLimMap_ULC_126_72", ParamsJSON{1}(), "https://nextcloud.awi.de/s/7F65JET9TzdosMD/download/PARA_ULC_126_72.json")
-)
+SoilHeatTile(upperbc::BoundaryProcess, lowerbc::BoundaryProcess, soilprofile::Profile, inputs::InputProvider, inits::CryoGrid.Initializer...; kwargs...) = SoilHeatTile(:H, upperbc, lowerbc, soilprofile, inputs, inits...; kwargs...)
const SamoylovDefault = (
soilprofile = SoilProfile(
diff --git a/src/Solvers/LiteImplicit/cglite_step.jl b/src/Solvers/LiteImplicit/cglite_step.jl
index 974d34eb..0605863a 100644
--- a/src/Solvers/LiteImplicit/cglite_step.jl
+++ b/src/Solvers/LiteImplicit/cglite_step.jl
@@ -7,7 +7,7 @@ function CryoGrid.perform_step!(integrator::CGLiteIntegrator)
p = integrator.p
dt = integrator.dt
t = tâ‚€ + dt
- tile = Tiles.resolve(Tile(integrator.sol.prob.f), p, t)
+ tile = Tiles.materialize(Tile(integrator.sol.prob.f), p, t)
# explicit update, if necessary
explicit_step!(integrator, tile, du, u, p, t)
# implicit update for energy state
diff --git a/src/Solvers/LiteImplicit/cglite_types.jl b/src/Solvers/LiteImplicit/cglite_types.jl
index 6c7872d6..47c993f9 100644
--- a/src/Solvers/LiteImplicit/cglite_types.jl
+++ b/src/Solvers/LiteImplicit/cglite_types.jl
@@ -32,7 +32,7 @@ function DiffEqBase.__init(prob::CryoGridProblem, alg::LiteImplicitEuler, args..
u_storage = [copy(u0)]
t_storage = [prob.tspan[1]]
# evaluate tile at initial condition
- tile = Tiles.resolve(Tile(prob.f), prob.p, t0)
+ tile = Tiles.materialize(Tile(prob.f), prob.p, t0)
tile(zero(u0), u0, prob.p, t0, dt)
# reset SavedValues on tile.data
initialsave = prob.savefunc(tile, u0, similar(u0))
diff --git a/src/Solvers/basic_solvers.jl b/src/Solvers/basic_solvers.jl
index 83422b84..0a6c3a53 100644
--- a/src/Solvers/basic_solvers.jl
+++ b/src/Solvers/basic_solvers.jl
@@ -28,7 +28,7 @@ function DiffEqBase.__init(prob::CryoGridProblem, alg::CGEuler, args...; dt=60.0
u_storage = [copy(u0)]
t_storage = [prob.tspan[1]*one(eltype(u0))]
# evaluate tile at initial condition
- tile = Tiles.resolve(Tile(prob.f), prob.p, t0)
+ tile = Tiles.materialize(Tile(prob.f), prob.p, t0)
tile(du0, u0, prob.p, t0, dt)
# reset SavedValues on tile.data
initialsave = prob.savefunc(tile, u0, similar(u0))
@@ -54,7 +54,7 @@ function perform_step!(integrator::CryoGridIntegrator{CGEuler})
tâ‚€ = integrator.t
p = integrator.p
initial_tile = Tile(integrator.sol.prob.f)
- tile = Tiles.resolve(initial_tile, p, tâ‚€)
+ tile = Tiles.materialize(initial_tile, p, tâ‚€)
# compute time derivative du
tile(du, u, p, tâ‚€)
dt = integrator.dt
diff --git a/src/Tiles/Tiles.jl b/src/Tiles/Tiles.jl
index 1cf7f226..72bf64c9 100644
--- a/src/Tiles/Tiles.jl
+++ b/src/Tiles/Tiles.jl
@@ -2,7 +2,7 @@ module Tiles
using CryoGrid
using CryoGrid: Parameterization, DynamicParameterization, DVar
-using CryoGrid.InputOutput: Forcing, CryoGridParams
+using CryoGrid.InputOutput: CryoGridParams
using CryoGrid.Numerics
using CryoGrid.Utils
diff --git a/src/Tiles/tile.jl b/src/Tiles/tile.jl
index ed4d2c29..3b57f17d 100644
--- a/src/Tiles/tile.jl
+++ b/src/Tiles/tile.jl
@@ -1,14 +1,15 @@
"""
- Tile{TStrat,TGrid,TStates,TInits,TEvents,iip} <: AbstractTile{iip}
+ Tile{TStrat,TGrid,TStates,TInits,TEvents,TInputs,iip} <: AbstractTile{iip}
Defines the full specification of a single CryoGrid tile; i.e. stratigraphy, grid, and state variables.
"""
-struct Tile{TStrat,TGrid,TStates,TInits,TEvents,iip} <: AbstractTile{iip}
+struct Tile{TStrat,TGrid,TStates,TInits,TEvents,TInputs,iip} <: AbstractTile{iip}
strat::TStrat # stratigraphy
grid::TGrid # grid
state::TStates # state variables
inits::TInits # initializers
events::TEvents # events
+ inputs::TInputs # inputs
data::TileData # output data
metadata::Dict # metadata
function Tile(
@@ -17,19 +18,21 @@ struct Tile{TStrat,TGrid,TStates,TInits,TEvents,iip} <: AbstractTile{iip}
state::TStates,
inits::TInits,
events::TEvents,
+ inputs::TInputs,
data::TileData=TileData(),
metadata::Dict=Dict(),
iip::Bool=true) where
- {TStrat<:Stratigraphy,TGrid<:Grid{Edges},TStates<:StateVars,TInits<:Tuple,TEvents<:NamedTuple}
- new{TStrat,TGrid,TStates,TInits,TEvents,iip}(strat,grid,state,inits,events,data,metadata)
+ {TStrat<:Stratigraphy,TGrid<:Grid{Edges},TStates<:StateVars,TInits<:Tuple,TEvents<:NamedTuple,TInputs<:InputProvider}
+ new{TStrat,TGrid,TStates,TInits,TEvents,TInputs,iip}(strat, grid, state, inits, events, inputs, data, metadata)
end
end
-ConstructionBase.constructorof(::Type{Tile{TStrat,TGrid,TStates,TInits,TEvents,iip}}) where {TStrat,TGrid,TStates,TInits,TEvents,iip} =
- (strat, grid, state, inits, events, data, metadata) -> Tile(strat, grid, state, inits, events, data, metadata, iip)
+ConstructionBase.constructorof(::Type{Tile{TStrat,TGrid,TStates,TInits,TEvents,TInputs,iip}}) where {TStrat,TGrid,TStates,TInits,TEvents,TInputs,iip} =
+ (strat, grid, state, inits, events, inputs, data, metadata) -> Tile(strat, grid, state, inits, events, inputs, data, metadata, iip)
# mark only stratigraphy and initializers fields as flattenable
Flatten.flattenable(::Type{<:Tile}, ::Type{Val{:strat}}) = true
Flatten.flattenable(::Type{<:Tile}, ::Type{Val{:inits}}) = true
Flatten.flattenable(::Type{<:Tile}, ::Type{Val{:events}}) = true
+# all other fields are non-flattenable by default
Flatten.flattenable(::Type{<:Tile}, ::Type{Val{name}}) where name = false
Base.show(io::IO, ::MIME"text/plain", tile::Tile{TStrat,TGrid,TStates,TInits,TEvents,iip}) where {TStrat,TGrid,TStates,TInits,TEvents,iip} = print(io, "Tile (iip=$iip) with layers $(keys(tile.strat)), $TGrid, $TStrat")
@@ -50,6 +53,7 @@ Constructs a `Tile` from the given stratigraphy and discretization strategy. `ar
function Tile(
@nospecialize(strat::Stratigraphy),
@nospecialize(discretization_strategy::DiscretizationStrategy),
+ @nospecialize(inputs::InputProvider),
@nospecialize(inits::CryoGrid.Initializer...);
metadata::Dict=Dict(),
cachetype::Type{T}=DiffCache,
@@ -82,11 +86,14 @@ function Tile(
_addlayerfield(init, Symbol(:init))
end
end
- return Tile(strat, grid, states, inits, (;events...), TileData(), metadata, iip)
+ tile = Tile(strat, grid, states, inits, (;events...), inputs, TileData(), metadata, iip)
+ _validate_inputs(tile, inputs)
+ return tile
end
Tile(strat::Stratigraphy, grid::Grid{Cells}, inits...; kwargs...) = Tile(strat, edges(grid), inits...; kwargs...)
Tile(strat::Stratigraphy, grid::Grid{Edges}, inits...; kwargs...) = Tile(strat, PresetGrid(grid), inits...; kwargs...)
-Tile(strat::Stratigraphy, inits...; discretization_strategy=AutoGrid(), kwargs...) = Tile(strat, discretization_strategy, inits...; kwargs...)
+Tile(strat::Stratigraphy, disc::DiscretizationStrategy, inits::CryoGrid.Initializer...; kwargs...) = Tile(strat, disc, InputFunctionProvider(), inits...; kwargs...)
+Tile(strat::Stratigraphy, args...; discretization_strategy=AutoGrid(), kwargs...) = Tile(strat, discretization_strategy, args...; kwargs...)
# convenience function to unwrap Tile from ODEFunction
function Tile(f::ODEFunction)
extract_f(tile::Tile) = tile
@@ -104,11 +111,11 @@ Constructs a `Tile` from a `SciMLBase` DEIntegrator.
function Tiles.Tile(integrator::SciMLBase.DEIntegrator)
tile = Tiles.Tile(integrator.sol.prob.f)
u = integrator.u
- return Tiles.resolve(tile, integrator.p, integrator.t)
+ return Tiles.materialize(tile, integrator.p, integrator.t)
end
"""
- computeprognostic!(_tile::Tile{TStrat,TGrid,TStates,TInits,TEvents,true}, _du, _u, p, t) where {TStrat,TGrid,TStates,TInits,TEvents}
+ computeprognostic!(_tile::Tile{TStrat,TGrid,TStates,TInits,TEvents,TInputs,true}, _du, _u, p, t) where {TStrat,TGrid,TStates,TInits,TEvents}
Time derivative step function (i.e. du/dt) for any arbitrary `Tile`. Specialized code is generated and compiled
on the fly via the @generated macro to ensure type stability. The generated code updates each layer in the stratigraphy
@@ -121,17 +128,17 @@ computeprognostic!(layer[i], ...)
```
"""
function computeprognostic!(
- _tile::Tile{TStrat,TGrid,TStates,TInits,TEvents,true},
+ _tile::Tile{TStrat,TGrid,TStates,TInits,TEvents,TInputs,true},
_du::AbstractVector,
_u::AbstractVector,
p::Union{Nothing,AbstractVector},
t::Number,
dt=1.0,
-) where {N,TStrat<:Stratigraphy{N},TGrid,TStates,TInits,TEvents}
+) where {N,TStrat<:Stratigraphy{N},TGrid,TStates,TInits,TEvents,TInputs}
_du .= zero(eltype(_du))
du = withaxes(_du, _tile)
u = withaxes(_u, _tile)
- tile = resolve(_tile, p, t)
+ tile = materialize(_tile, p, t)
strat = tile.strat
state = TileState(tile.strat, tile.grid, tile.state, du, u, t, dt)
CryoGrid.resetfluxes!(strat, state)
@@ -154,7 +161,7 @@ Computes the maximum permissible forward timestep for this `Tile` given the curr
function CryoGrid.timestep(_tile::Tile, _du, _u, p, t)
du = withaxes(_du, _tile)
u = withaxes(_u, _tile)
- tile = resolve(_tile, p, t)
+ tile = materialize(_tile, p, t)
strat = tile.strat
state = TileState(tile.strat, tile.grid, tile.state, du, u, t, 1.0)
dtmax = CryoGrid.timestep(strat::Stratigraphy, state)
@@ -176,7 +183,7 @@ function CryoGrid.initialcondition!(tile::Tile, tspan::NTuple{2,Float64}, p::Abs
utype = isempty(p) ? eltype(tile.state.uproto) : eltype(p)
du = zero(similar(tile.state.uproto, utype))
u = zero(similar(tile.state.uproto, utype))
- tile = resolve(tile, p, t0)
+ tile = materialize(tile, p, t0)
strat = tile.strat
state = TileState(tile.strat, tile.grid, tile.state, du, u, t0, 1.0)
CryoGrid.initialcondition!(tile.grid, state)
@@ -269,28 +276,28 @@ getstate(integrator::SciMLBase.DEIntegrator) = Tiles.getstate(Tile(integrator),
"""
Numerics.getvar(var::Symbol, integrator::SciMLBase.DEIntegrator; interp=true) = Numerics.getvar(Val{var}(), Tile(integrator), integrator.u; interp)
-"""
- parameterize(tile::Tile)
-
-Adds parameter information to all nested types in `tile` by recursively calling `parameterize`.
-"""
-function CryoGrid.parameterize(tile::Tile)
- ctor = ConstructionBase.constructorof(typeof(tile))
- new_layers = map(namedlayers(tile.strat)) do named_layer
- name = nameof(named_layer)
- layer = CryoGrid.parameterize(named_layer.val)
- name => _addlayerfield(layer, name)
- end
- new_inits = map(tile.inits) do init
- _addlayerfield(CryoGrid.parameterize(init), :init)
- end
- new_events = map(keys(tile.events)) do name
- evs = map(CryoGrid.parameterize, getproperty(tile.events, name))
- name => _addlayerfield(evs, name)
- end
- new_strat = Stratigraphy(boundaries(tile.strat), (;new_layers...))
- return ctor(new_strat, tile.grid, tile.state, new_inits, (;new_events...), tile.data, tile.metadata)
-end
+# """
+# parameterize(tile::Tile)
+
+# Adds parameter information to all nested types in `tile` by recursively calling `parameterize`.
+# """
+# function CryoGrid.parameterize(tile::Tile)
+# ctor = ConstructionBase.constructorof(typeof(tile))
+# new_layers = map(namedlayers(tile.strat)) do named_layer
+# name = nameof(named_layer)
+# layer = CryoGrid.parameterize(named_layer.val)
+# name => _addlayerfield(layer, name)
+# end
+# new_inits = map(tile.inits) do init
+# _addlayerfield(CryoGrid.parameterize(init), :init)
+# end
+# new_events = map(keys(tile.events)) do name
+# evs = map(CryoGrid.parameterize, getproperty(tile.events, name))
+# name => _addlayerfield(evs, name)
+# end
+# new_strat = Stratigraphy(boundaries(tile.strat), (;new_layers...))
+# return ctor(new_strat, tile.grid, tile.state, new_inits, (;new_events...), tile.data, tile.metadata)
+# end
"""
variables(tile::Tile)
@@ -316,30 +323,40 @@ withaxes(u::AbstractArray, tile::Tile) = ComponentArray(u, getaxes(tile.state.up
withaxes(u::ComponentArray, ::Tile) = u
"""
- resolve(tile::Tile, p, t)
+ materialize(tile::Tile, p, t)
-Resolves/updates the given `tile` by:
-(1) Replacing all `ModelParameters.AbstractParam` values in `tile` with their (possibly updated) value from `p`.
-(2) Resolving the boundary depths of the `Stratigraphy` layers by invoking `resolveboundaries`.
-(3) Replacing all instances of `DynamicParameterization` with their resolved values given the current state.
+Materializes the given `tile` by:
+- Replacing all `ModelParameters.AbstractParam` values in `tile` with their (possibly updated) value from `p`.
+- Evaluating and replacing all `DynamicParameterization`s given the time `t`.
+- Evaluating and replacing all `Input`s given the time `t`.
Returns the reconstructed `Tile` instance.
"""
-function resolve(tile::Tile{TStrat,TGrid,TStates}, p::AbstractVector, t::Number) where {TStrat,TGrid,TStates}
- IgnoreTypes = Union{TGrid,TStates,TileData,Unitful.Quantity,Numerics.ForwardDiff.Dual}
+function materialize(tile::Tile, p::AbstractVector, t::Number)
+ IgnoreTypes = _ignored_types(tile)
+ # ==== Update parameter values ==== #
# unfortunately, reconstruct causes allocations due to a mysterious dynamic dispatch when returning the result of _reconstruct;
# I really don't know why, could be a compiler bug, but it doesn't happen if we call the internal _reconstruct method directly...
# so that's what we do here. The last integer argument denotes the index of the first parameter.
- tile_updated = Flatten._reconstruct(tile, p, Flatten.flattenable, ModelParameters.AbstractParam, IgnoreTypes, 1)[1]
- dynamic_ps = Flatten.flatten(tile_updated, Flatten.flattenable, DynamicParameterization, IgnoreTypes)
+ parameterized_tile = Flatten._reconstruct(tile, p, Flatten.flattenable, ModelParameters.AbstractParam, IgnoreTypes, 1)[1]
+ # call materialize on parameterized tile
+ return materialize(parameterized_tile, nothing, t)
+end
+function materialize(tile::Tile, ::Nothing, t::Number)
+ IgnoreTypes = _ignored_types(tile)
+ # ==== Compute dynamic parameter values ==== #
# TODO: perhaps should allow dependence on local layer state;
# this would likely require deconstruction/reconstruction of layers in order to
# build the `LayerState`s and evaluate the dynamic parameters in a fully type stable manner.
+ dynamic_ps = Flatten.flatten(tile, Flatten.flattenable, DynamicParameterization, IgnoreTypes)
dynamic_values = map(d -> d(t), dynamic_ps)
- reconstructed_tile = Flatten._reconstruct(tile_updated, dynamic_values, Flatten.flattenable, DynamicParameterization, IgnoreTypes, 1)[1]
- return reconstructed_tile
+ tile2 = Flatten._reconstruct(tile, dynamic_values, Flatten.flattenable, DynamicParameterization, IgnoreTypes, 1)[1]
+ # ==== Compute input values ==== #
+ inputs = Flatten.flatten(tile2, Flatten.flattenable, Input, IgnoreTypes)
+ input_values = map(input -> tile2.inputs(input, t), inputs)
+ concrete_tile = Flatten._reconstruct(tile2, input_values, Flatten.flattenable, Input, IgnoreTypes, 1)[1]
+ return concrete_tile
end
-resolve(tile::Tile, p::Nothing, t::Number) = tile
function checkstate!(tile::Tile, state::TileState, u, du, label::Symbol)
if CryoGrid.CRYOGRID_DEBUG
@@ -348,7 +365,7 @@ function checkstate!(tile::Tile, state::TileState, u, du, label::Symbol)
debughook!(tile, state, AssertionError("[$label] Found NaN/Inf value in current state vector at index $i"))
end
if !isfinite(du[i])
- debughook!(tile, state, AssertionError("[$label] Found NaN/Inf value in computed time derivatives at index $i"))
+ debughook!(tile, state, AssertionError("[$label] Found NaN/Inf value in current time derivatives at index $i"))
end
end
end
@@ -379,4 +396,17 @@ function _initstatevars(@nospecialize(strat::Stratigraphy), @nospecialize(grid::
return states
end
+function _validate_inputs(@nospecialize(tile::Tile), inputprovider::InputProvider)
+ IgnoreTypes = _ignored_types(tile)
+ inputs = Flatten.flatten(tile, Flatten.flattenable, Input, IgnoreTypes)
+ names = keys(inputprovider)
+ for input in inputs
+ name = nameof(input)
+ @assert name ∈ names "input $name does not exist in the given input provider with keys $names; check for typos or mismatched names"
+ end
+end
+
+# helper method that returns a Union type of all types that should be ignored by Flatten.flatten
+@inline _ignored_types(::Tile{TStrat,TGrid,TStates}) where {TStrat,TGrid,TStates} = Union{TGrid,TStates,TileData,Unitful.Quantity,Numerics.ForwardDiff.Dual}
+
# ===================================================================== #
diff --git a/src/Presets/presetgrids.jl b/src/default_grids.jl
similarity index 100%
rename from src/Presets/presetgrids.jl
rename to src/default_grids.jl
diff --git a/src/initializers.jl b/src/initializers.jl
index 5167f029..6a7b03c6 100644
--- a/src/initializers.jl
+++ b/src/initializers.jl
@@ -3,9 +3,6 @@ Base.isless(init1::Initializer, init2::Initializer) = false
# add varname dispatch for initializer types
CryoGrid.varname(::VarInitializer{varname}) where {varname} = varname
-# default behavior is to not automatically parameterize initializers
-CryoGrid.parameterize(init::VarInitializer) = init
-
# default to invoking initializer
CryoGrid.initialcondition!(init::VarInitializer, layer::Layer, state) = init(layer, state)
diff --git a/src/methods.jl b/src/methods.jl
index b9511218..0017f588 100644
--- a/src/methods.jl
+++ b/src/methods.jl
@@ -134,32 +134,6 @@ if the prognostic state was modified and `false` otherwise. Defaults to returnin
"""
diagnosticstep!(layer::Layer, state) = false
-"""
- parameterize(x::T) where {T}
- parameterize(x::Unitful.AbstractQuantity; props...)
- parameterize(p::AbstractParam; ignored...)
-
-Recursively wraps `x` or nested numeric quantities in `x` with `Param` to mark them as parameters.
-If `x` is already a `Param` type, `x` will be returned as-is.
-If `x` is a numeric type, `x` will be wrapped in `Param` with associated properties `props`.
-If `x` is a struct type, `x` will be recursively unpacked and `parameterize` called on each field.
-"""
-parameterize(x::Number; type=Param, props...) = type(x; props...)
-parameterize(x::Unitful.AbstractQuantity; type=Param, props...) = type(ustrip(x); untis=unit(x), props...)
-parameterize(p::AbstractParam; ignored...) = p
-parameterize(f::Function; ignored...) = f
-function parameterize(x::T; props...) where {T}
- # get field names of T, if available
- T_fieldnames = isabstracttype(T) ? () : fieldnames(T)
- # invoke parameterize on all fields
- new_fields = map(T_fieldnames) do fieldname
- fieldvalue = getfield(x, fieldname)
- parameterize(fieldvalue)
- end
- ctor = ConstructionBase.constructorof(T)
- return ctor(new_fields...)
-end
-
"""
initializers(::Layer)
initializers(::Layer, ::Process)
diff --git a/test/Physics/Surface/seb_tests.jl b/test/Physics/Surface/seb_tests.jl
index e67ac24e..c7762670 100644
--- a/test/Physics/Surface/seb_tests.jl
+++ b/test/Physics/Surface/seb_tests.jl
@@ -1,22 +1,16 @@
Surface.surfaceproperties(seb::SurfaceEnergyBalance, sub::TestGroundLayer) = seb.para.soil
@testset "Surface Energy Balance" begin
- const_Tair = ConstantForcing(0.0u"°C", :Tair)
- const_pr = ConstantForcing(upreferred(1.0u"atm"), :pr)
- const_qh = ConstantForcing(0.001u"kg/kg", :qh)
- const_wind = ConstantForcing(1.0u"m/s", :wind)
- const_Lin = ConstantForcing(upreferred(300.0u"W/m^2"), :Lin)
- const_Sin = ConstantForcing(upreferred(200.0u"W/m^2"), :Sin)
- z = 2.0u"m"
+ Tair = 0.0u"°C"
+ pr = upreferred(1.0u"atm")
+ qh = 0.001u"kg/kg"
+ wind = 1.0u"m/s"
+ Lin = upreferred(300.0u"W/m^2")
+ Sin = upreferred(200.0u"W/m^2")
+ forcings = (; Tair, pr, qh, wind, Lin, Sin)
@testset "Iterative" begin
seb = SurfaceEnergyBalance(
- const_Tair,
- const_pr,
- const_qh,
- const_wind,
- const_Lin,
- const_Sin,
- z,
+ forcings;
solscheme = Surface.Iterative(),
)
# TODO: would be good to have SEB work correctly with units to verify correctness;
@@ -36,13 +30,7 @@ Surface.surfaceproperties(seb::SurfaceEnergyBalance, sub::TestGroundLayer) = seb
end
@testset "Analytical" begin
seb = SurfaceEnergyBalance(
- const_Tair,
- const_pr,
- const_qh,
- const_wind,
- const_Lin,
- const_Sin,
- z,
+ forcings;
solscheme = Surface.Analytical(),
)
seb = pstrip(seb)
@@ -59,15 +47,9 @@ Surface.surfaceproperties(seb::SurfaceEnergyBalance, sub::TestGroundLayer) = seb
end
@testset "Numerical" begin
seb = SurfaceEnergyBalance(
- const_Tair,
- const_pr,
- const_qh,
- const_wind,
- const_Lin,
- const_Sin,
- z,
- solscheme=Surface.Numerical()
- )
+ forcings;
+ solscheme=Surface.Numerical(),
+ )
seb = pstrip(seb)
grid = Grid([0.0,0.1]u"m")
toplayer = Top(seb)
diff --git a/test/benchmarks.jl b/test/benchmarks.jl
index 57680443..87c01f4f 100644
--- a/test/benchmarks.jl
+++ b/test/benchmarks.jl
@@ -78,8 +78,8 @@ StatsPlots.groupedbar(xnames, results_freeW_tair.error[:,:,1], yerror=results_fr
savefig("solver_benchmark_freeW_tair_error.png")
# Test 2: Energy + SEB upper bc + free water freeze curve
-z = 2.; # height [m] for which the forcing variables (Temp, humidity, wind, pressure) are provided
-seb = SurfaceEnergyBalance(Tair,pr,q,wind,Lin,Sin,z)
+forcings = (; Tair,pr,q,wind,Lin,Sin)
+seb = SurfaceEnergyBalance()
solvers = [Euler, DP5, ROCK2, ROCK4]
results_freeW_seb = run_solver_benchmark(solvers, grids, dts, tspan, seb)
serialize("solver_benchmark_freeW_seb.ser", results_freeW_seb)
--
GitLab