diff --git a/Project.toml b/Project.toml
index d831f333b38578798f449a8f23e1848f9a2f42e3..bdc88c13834039a21c4a1af20ae675a82e5ec0aa 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 d9fb6fb091b32f6adf2361048f48d3d687cec269..7d41667a47a7faf284043a246fdfef29421a265b 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 e58d8030fc029780948a000c2996cfdc20e9425e..03357710a8154c4c8cd6f447ced8154f8f156534 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 e3b29ccfe38a21d5abd391795cfa7b1865f5f215..2d887e003a9e1d0a1ee98299b674b802e6694fa9 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 1eec62a5633f21bb7778954e0b0198a2451eaf43..91bb16c8143935a083c52e3a29916cba06831d0c 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 a93800624378ccf6bf3340fd1bcde6f952ebfb40..96d6762546aca339005bf94d9726ef2b6b8ba548 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 7540a8555279ece68793d9e1b61ac7271e6bfdc7..c06904729a386e354d93503f19d70dcf37d8d2d5 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 712675d6d9fcea9840764d7c3b4c709127f6d104..eb5b993126799fab3fe36b4c50aa12d5338be341 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 8ba3708c22fd68e7d65626a47652c1533405ea8c..e2eb192b1713ba630d7b8b914b3dbe0bd683b936 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 12f7917bd06d99452cda88c1ab1c211e51454260..4248173e186576545ba3f53daf17cc8f10d13d80 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 80aaef5e8b5889330b5e73f46132fa9bc8fb8202..2a42322224ff9297e48c922294b4bd661c5632d1 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 2f7bf0fbcb7ccf38d8f390479590d7102dd7efa4..11cee1160fc93ed5e5008d37d4444574b8f78c9a 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 bbdf09cc9493808fb5917ec3916cb1b0cd362eea..36aada3402318d30617c5748d24f5b211be30c5f 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 72bf5839cb25cf8cab69af2205e8f1f3d06bbe9c..3b6b81192cb4a7f9177701624e1851fcc295dd73 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 5a1e389e908b31436ce9e7359054a41ef279ebbb..5b88521e0d56bae9f67a88819fa8cc649ab7f689 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 c3dfbf7814a95f2e865d111a4c5c0fc24719c1cb..bad4afaa280860cd99f62a43c8eed3a23314875d 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 ec01a181d3d6eaa31b66985e688ce6b2b84d0209..b055d861aab31c6a0201e7e8ed4275193ee83493 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 933106aedbc1a03ea37fa5adb5cba9c225c66ed2..cda880945d3c75ac19deddd0fd92fffe923ef85b 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 ed60224b2e1c56bd6dd074c46ea4958689890388..3c56b0d8593fab94e1629d3aa9cd032b73992de3 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 d17509e95b6876e78e2e07bf7fdc725de7958103..2f3868fe9e1ef9b7431b6c8551d293232dfecfee 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 0000000000000000000000000000000000000000..dfbaf81bb8169a32e7bd08c55a6b7185f17728ac
--- /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 0000000000000000000000000000000000000000..4fbf632fa7d9533503aa67cda33f449736ce2e44
--- /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 a09ec57877a632da7b53a989ad231e605ba0ccb3..0a35804213784ec05fb83c0602c829f543e30de5 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 785230f74da2649b4968391a4659bec16745eee8..df28efd59d403e85ef154f443dbe6d90da86db52 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 254527d8d0e65c4c9aa831ed941375896a9e3d60..ad5ceeeddd50f74f81f2ce6c773a99464a890123 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 e3447bc4fb1cd819c24c2af4237a87bd85e7000b..7f198c5e65d280ac69b577bb4ca14d2a8197ba13 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 d4772cc5d00af02d2fad5aafd5ced6ab358bd46d..66bce2bc4327293f623b67a5eabf32d4505870e1 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 8ea5878298ba8d9961c3055864e05acdef3b11d2..18a403e5cab4a767e1fb9d32293db568a791e043 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 2f5e5e172e23932c1fc44c47c5f74f6e8f5022f9..fe6c3e7e3fcf95a8746f37c9fabea97e079942b0 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 b9cb5666d352a5a06b70c972f5ed93df2d09cfac..f5b2c72dd812d10a2907ecbc3bdd046a681963c2 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 87f4efefe6773c92be6118f5e275f680f3fe7c78..cef31f8fb3a04a365f67be4b5da573bca8967c41 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 2052ae5d4295ed51a6227e1085f058b661ba2068..213b902b0a277d1675f925a17a2afb111d3f6ce5 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 36790af9899e8780335e0d46fb98fc51f72d608f..fd53ef4e66933fe0c436e59f62063e90fbcf0296 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 6bc234d58c68b69a52b20af88e1256c2c176010d..4a98aa7b9b0fa9e8fb37fd1b1c9fd9f003f20524 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 c3780d7bf5b2254eb1d9dc528ef34da186231a64..e5af88bed2d815aa167684e2613132c3ad70f77c 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 39cffc441c473a69b34d8d762a8d5c2cbfcfb466..d2703a4887ed8bcfdb12040756ea06976a66aa20 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 3003d3df44e0e77701f34c0d2454533d49d77c17..1c05ae3b23840b6f240d2414e5e7e99b85e2dc29 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 11bd26e3b49e2b007705d08f1a1b39050ec2dbe7..eb0333440763a0f394e08cb151f80c178caa7f80 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 d027a53f4975a512254f23979f75e4813f361de3..2978ec7e9cf19b78341bab86a69e68df7660dc65 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 974d34eb020b956dbbdbb90d463065e01db3d47f..0605863a8b19a42f60bb017ef97c983b348e2f1b 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 6c7872d6df35ce37d9a46515eb813584a993bc76..47c993f91d840839ecd1e45eca97049199fb59df 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 83422b84da625705ad8907012b4c0aa6e3d8220a..0a6c3a53bb89a3daa1f148a4be6ec5f05e38633c 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 1cf7f226550cd266945dd34c4e8443be9f469171..72bf64c900be2c04f56c81cf48b4e567623225d4 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 ed4d2c295b07729492f45188d08f55e24c2bd703..3b57f17d884aac90853f5506faaa0803f17d91b4 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 5167f029059f0ff82dca5aa72ef3570a10dc8a33..6a7b03c627e950cf31e8b20c7b5107bfbddabd36 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 b9511218b2698fe8532055d60e2de9327cf7252c..0017f58884e8f5b110ab1d850cf2b1f0958728da 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 e67ac24eda3ef5e1e428a53eac7d9c38a4bbc068..c77626706b4246b3a4f4a2a961fbac46dd1446a5 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 57680443b97fa31a008ceca7d5a338ed44b79fcd..87c01f4f0c7a6146d96635ec8a73988dcb0077d3 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)