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..2127659debce0107841bfd6d9bc5a06c39ef59e0 100755
--- a/README.md
+++ b/README.md
@@ -49,14 +49,14 @@ 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
+soilprofile, tempprofile = CryoGrid.SamoylovDefault
initT = initializer(:T, tempprofile)
# choose grid with 5cm spacing
-grid = CryoGrid.Presets.DefaultGrid_5cm
+grid = CryoGrid.DefaultGrid_5cm
# build a default Tile with heat conduction on the given soil profile
-tile = CryoGrid.Presets.SoilHeatTile(
+tile = CryoGrid.SoilHeatTile(
TemperatureBC(forcings.Tair),
GeothermalHeatFlux(0.053u"W/m^2"),
soilprofile,
diff --git a/docs/src/manual/architecture.md b/docs/src/manual/architecture.md
index 7bfd05e41342e286e264f579dc36b951979f825e..5cadab0602ac9a2db344c2128fb1e4f708326333 100644
--- a/docs/src/manual/architecture.md
+++ b/docs/src/manual/architecture.md
@@ -54,7 +54,7 @@ The `CryoGrid` module defines three primary methods that can be used to implemen
1. [`computediagnostic!`](@ref) updates all (non-flux) state variables and/or derived quantities based on the current (prognostic) state.
2. [`interact!`](@ref) defines interactions between adjacent layers in the stratigraphy, including fluxes over the layer boundary.
-3. [`computefluxes!`](@ref) computes all internal fluxes (and the divergence thereof) within each layer, after boundary fluxes are taken into account by `interact!`.
+3. [`computeprognostic!`](@ref) computes all internal fluxes (and the divergence thereof) within each layer, after boundary fluxes are taken into account by `interact!`.
Layer and/or process specific implementations of each of these methods can generally assume that the previous methods have already been invoked by the caller (it is the responsibility of the calling code to ensure that this is the case). This is, for example, the order in which these methods will be invoked by `tile(du, u, p t)`.
@@ -65,7 +65,7 @@ using CryoGrid
using CryoGrid.Diagnostics
soil = Ground()
-grid = CryoGrid.Presets.DefaultGrid_5cm
+grid = CryoGrid.DefaultGrid_5cm
state = Diagnostics.build_dummy_state(grid, soil)
@which CryoGrid.computediagnostic!(soil, state)
@@ -82,7 +82,7 @@ In order to facilitate modularity and ease-of-use, CryoGrid.jl provides an autom
`Prognostic`(@ref) state variables fully define the state of the system at any given time `t`. They form what is typically called the "phase space" or "state space" in the mathematics and engineering literature. In order to be compatible with standard ODE solvers (e.g. like those in `OrdinaryDiffEq`), CryoGrid.jl automatically assembles prognostic state variables into a single array `u` (and its corresponding time derivative `du`) which is returned when initializing a `Tile` with the `initialcondition!` method. Note again that this array should always fully define the state of the system.
-`Diagnostic`(@ref) state variables act as caches for intermediate and derived quantities defined by the model. They also may, in some cases, provide a means of coupling between different processes (e.g. the heat and water flux variables `jH` and `jw` might be updated by more than one `Process`). For any model configuration, all diagnostic variables should be fully updated (and thus consistent) with the given prognostic state after invoking `computediagnostic!`, `interact!`, and `computefluxes!`.
+`Diagnostic`(@ref) state variables act as caches for intermediate and derived quantities defined by the model. They also may, in some cases, provide a means of coupling between different processes (e.g. the heat and water flux variables `jH` and `jw` might be updated by more than one `Process`). For any model configuration, all diagnostic variables should be fully updated (and thus consistent) with the given prognostic state after invoking `computediagnostic!`, `interact!`, and `computeprognostic!`.
When a `Tile` is constructed, all variables defined by each layer in the `Stratigraphy` are collected and then intiailized in [`StateVars`](@ref) according to the given `DiscretizationStrategy`.
diff --git a/docs/src/manual/overview.md b/docs/src/manual/overview.md
index bb5cd1ee25704ccceac96d5ed211f58273473b88..414f950a831bd3945b6014483e4099e7e641bd1b 100644
--- a/docs/src/manual/overview.md
+++ b/docs/src/manual/overview.md
@@ -12,12 +12,13 @@ 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"))
);
-grid = CryoGrid.Presets.DefaultGrid_5cm
+grid = CryoGrid.DefaultGrid_5cm
# define initial conditions for temperature using a given profile;
# The default initializer linearly interpolates between profile points.
initT = initializer(:T, tempprofile)
@@ -69,12 +70,12 @@ When the `HeatBalance` process is assigned to a `Soil` layer, `Tile` will invoke
Each variable definition consists of a name (a Julia `Symbol`), a type, and a shape. For variables discretized on the grid, the shape is specified by `OnGrid`, which will generate an array of the appropriate size when the model is compiled. The arguments `Cells` and `Edges` specify whether the variable should be defined on the grid cells or edges respecitvely.
-The real work finally happens in [`computediagnostic!`](@ref) and [`computefluxes!`](@ref), the latter of which should be used to compute the time derivatives (here `dH`). [`interact!`](@ref) defines the behavior at the boundaries and should be used to compute the derivatives (and any other necessary values) at the interface between layers.
+The real work finally happens in [`computediagnostic!`](@ref) and [`computeprognostic!`](@ref), the latter of which should be used to compute the time derivatives (here `dH`). [`interact!`](@ref) defines the behavior at the boundaries and should be used to compute the derivatives (and any other necessary values) at the interface between layers.
-We can take as an example the implementation of `computefluxes!` for enthalpy-based heat conduction (note that `jH` is a diagnostic variable representing the energy flux over each cell edge):
+We can take as an example the implementation of `computeprognostic!` for enthalpy-based heat conduction (note that `jH` is a diagnostic variable representing the energy flux over each cell edge):
```julia
-function CryoGrid.computefluxes!(::SubSurface, ::HeatBalance{<:EnthalpyBased}, state)
+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
diff --git a/docs/src/quickstart.md b/docs/src/quickstart.md
index e3b29ccfe38a21d5abd391795cfa7b1865f5f215..75c28cb83857ba96b4ceb0ec547dd514a05d3b38 100644
--- a/docs/src/quickstart.md
+++ b/docs/src/quickstart.md
@@ -1,6 +1,6 @@
## Quick start
-After [installing](installation.md) `CryoGrid.jl`, you can get started right away with a simple soil heat model. The [`Presets`](@ref) module (aliased `CryoGrid.Presets`) provides pre-specified model configurations that can be obtained with a single function call. It is also possible to provide a custom soil and temperature profile to `SoilHeatTile`; here the values given in `Presets.SamoylovDefault` is used.
+After [installing](installation.md) `CryoGrid.jl`, you can get started right away with a simple soil heat model. The [`Presets`](@ref) module (aliased `CryoGrid.Presets`) provides pre-specified model configurations that can be obtained with a single function call. It is also possible to provide a custom soil and temperature profile to `SoilHeatTile`; here the values given in `SamoylovDefault` is used.
```julia
using CryoGrid
@@ -8,14 +8,14 @@ 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
+soilprofile, tempprofile = CryoGrid.SamoylovDefault
initT = initializer(:T, tempprofile)
# choose grid with 5cm spacing
-grid = CryoGrid.Presets.DefaultGrid_5cm
+grid = CryoGrid.DefaultGrid_5cm
# build Tile from the given soil and temperature profiles
-tile = CryoGrid.Presets.SoilHeatTile(TemperatureBC(forcings.Tair), GeothermalHeatFlux(0.053u"W/m^2"), soilprofile, initT, grid=grid)
+tile = CryoGrid.SoilHeatTile(TemperatureBC(forcings.Tair), GeothermalHeatFlux(0.053u"W/m^2"), soilprofile, initT, grid=grid)
# define time span (1 year)
tspan = (DateTime(2010,11,30),DateTime(2011,11,30))
u0, du0 = initialcondition!(tile, tspan)
diff --git a/examples/cglite_parameter_ensembles.jl b/examples/cglite_parameter_ensembles.jl
index 1eec62a5633f21bb7778954e0b0198a2451eaf43..4c6c40e53418828e36ff15dae11ee1ffe42779f2 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),
@@ -42,7 +42,7 @@ strat = Stratigraphy(
soil_layers,
z_bot => Bottom(GeothermalHeatFlux(0.053u"W/m^2"))
);
-modelgrid = CryoGrid.Presets.DefaultGrid_2cm
+modelgrid = CryoGrid.DefaultGrid_2cm
tile = Tile(strat, modelgrid, 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 10 years.
diff --git a/examples/heat_freeW_bucketW_samoylov.jl b/examples/heat_freeW_bucketW_samoylov.jl
index a93800624378ccf6bf3340fd1bcde6f952ebfb40..355b734e208a175d9dfc6dcf88198eed22599802 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);
-_, tempprofile = CryoGrid.Presets.SamoylovDefault;
+forcings = loadforcings(CryoGrid.Forcings.Samoylov_ERA_obs_fitted_1979_2014_spinup_extended_2044);
+_, tempprofile = CryoGrid.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.
@@ -20,8 +23,8 @@ strat = @Stratigraphy(
2.0u"m" => Ground(SimpleSoil(por=0.10,sat=1.0,org=0.0), heat=HeatBalance(), water=WaterBalance(BucketScheme())),
1000.0u"m" => Bottom(GeothermalHeatFlux(0.053u"W/m^2")),
);
-modelgrid = CryoGrid.Presets.DefaultGrid_2cm
-tile = Tile(strat, modelgrid, initT, initsat);
+modelgrid = CryoGrid.DefaultGrid_2cm
+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..2cc4cbaafb12ce4dda6f7e54e43ab359e24f4f69 100644
--- a/examples/heat_freeW_lake_lite_implicit.jl
+++ b/examples/heat_freeW_lake_lite_implicit.jl
@@ -6,7 +6,12 @@ Plots.plotly()
CryoGrid.debug(true)
-forcings = loadforcings(CryoGrid.Presets.Forcings.Samoylov_ERA_MkL3_CCSM4_long_term);
+raw_forcings = loadforcings(CryoGrid.Forcings.Samoylov_ERA_MkL3_CCSM4_long_term);
+Tair = raw_forcings.data.Tair
+Ptot = uconvert.(u"m/s", raw_forcings.data.Ptot)
+rainfall = Ptot.*(Tair .> 0u"°C")
+snowfall = Ptot.*(Tair .<= 0u"°C")
+forcings = rebuild(raw_forcings; Tair, rainfall, snowfall);
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),
@@ -20,11 +25,11 @@ tempprofile_linear = TemperatureProfile(
10.0u"m" => -10.0u"°C",
1000.0u"m" => 10.2u"°C"
)
-modelgrid = Grid(vcat(-1.0u"m":0.02u"m":-0.02u"m", CryoGrid.Presets.DefaultGrid_2cm))
+modelgrid = Grid(vcat(-1.0u"m":0.02u"m":-0.02u"m", CryoGrid.DefaultGrid_2cm))
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 +40,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..f818f2f8c8382bbe5ec0824c6f38de2b2c7922de 100644
--- a/examples/heat_freeW_lite_implicit.jl
+++ b/examples/heat_freeW_lite_implicit.jl
@@ -7,17 +7,23 @@
using CryoGrid
using CryoGrid.LiteImplicit
-# Load forcings and build stratigraphy like before.
-forcings = loadforcings(CryoGrid.Presets.Forcings.Samoylov_ERA_MkL3_CCSM4_long_term);
-forcings = Base.rename(forcings, :Ptot => :precip)
-precip_values = forcings.precip.interpolant.coefs
-precip_values .*= 2
+# First we prepare the forcing data. This data includes only air temperature and
+# *total* precipitation, so we need to separate this into rainfall and snowfall
+# using air temperature.
+raw_forcings = loadforcings(CryoGrid.Forcings.Samoylov_ERA_MkL3_CCSM4_long_term);
+Tair = raw_forcings.data.Tair
+Ptot = uconvert.(u"m/s", raw_forcings.data.Ptot)
+rainfall = Ptot.*(Tair .> 0u"°C")
+snowfall = Ptot.*(Tair .<= 0u"°C")
+forcings = rebuild(raw_forcings; Tair, rainfall, snowfall);
+
+# Next we define the boundary conditions and stratigraphy.
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(
@@ -29,22 +35,22 @@ soilprofile = SoilProfile(
)
heat = HeatBalance(heatop)
water = WaterBalance()
-soil_layers = map(para -> Ground(para; heat, water), soilprofile)
+soil_layers = map(para -> Ground(para; heat, water), soilprofile);
strat = @Stratigraphy(
z_top => Top(upperbc),
z_top => Snowpack(Snow.LiteGridded(); heat),
soil_layers...,
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);
+modelgrid = Grid(vcat(z_top:0.02u"m":-0.02u"m", CryoGrid.DefaultGrid_2cm))
+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.
-tspan = (DateTime(1990,10,1), DateTime(2020,10,1))
+# Here we specify a time span of 20 years.
+tspan = (DateTime(2000,10,1), DateTime(2020,10,1))
u0, du0 = initialcondition!(tile, tspan);
-prob = CryoGridProblem(tile, u0, tspan, saveat=24*3600.0, savevars=(:T,:θw,:θwi,:kc,:Ïsn))
-sol = @time solve(prob, LiteImplicitEuler())
+prob = CryoGridProblem(tile, u0, tspan, saveat=24*3600, savevars=(:T,:θw))
+sol = @time solve(prob, LiteImplicitEuler(), dt=24*3600)
out = CryoGridOutput(sol)
# Plot the results!
@@ -54,20 +60,6 @@ cg = Plots.cgrad(:copper,rev=true);
Plots.plot(out.T[Z(Near(zs))], color=cg[LinRange(0.0,1.0,length(zs))]', ylabel="Temperature", title="", leg=false, dpi=150)
Plots.plot(out.snowpack.dsn)
-integrator = init(prob, LiteImplicitEuler())
-for i in integrator
- state = getstate(integrator)
- last_above_ground_idx = findlast(<(0), cells(state.grid))
- if state.snowpack.dsn[1] > 0.0 && state.snowpack.swe[last_above_ground_idx] <= 0 && state.snowpack.snowfall[1] <= 0
- println("stopping")
- break
- end
-end
-
-CryoGrid.debug(true)
-
-step!(integrator)
-
# CryoGridLite can also be embedded into integrators from OrdinaryDiffEq.jl via the `NLCGLite` nonlinear solver interface.
# Note that these sovers generally will not be faster (in execution time) but may be more stable in some cases. Adaptive timestepping can be employed by
# removing the `adaptive=false` argument.
diff --git a/examples/heat_freeW_samoylov.jl b/examples/heat_freeW_samoylov.jl
index 8ba3708c22fd68e7d65626a47652c1533405ea8c..a54e08068be210bfbd23febabb54be83d0e00c5b 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)
-tile = CryoGrid.Presets.SoilHeatTile(
+initT = initializer(:T, CryoGrid.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.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
);
@@ -55,4 +59,4 @@ out = CryoGridOutput(sol)
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)
+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)
\ No newline at end of file
diff --git a/examples/heat_freeW_seb_snow_bucketW_samoylov.jl b/examples/heat_freeW_seb_snow_bucketW_samoylov.jl
index 12f7917bd06d99452cda88c1ab1c211e51454260..2f47c9d5f74688b974b886933a78f7d3eb960991 100644
--- a/examples/heat_freeW_seb_snow_bucketW_samoylov.jl
+++ b/examples/heat_freeW_seb_snow_bucketW_samoylov.jl
@@ -8,7 +8,7 @@ using CryoGrid
using OrdinaryDiffEq
# First, load the forcings and construct the Tile.
-modelgrid = CryoGrid.Presets.DefaultGrid_2cm;
+modelgrid = CryoGrid.DefaultGrid_2cm;
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),
@@ -17,13 +17,12 @@ soilprofile = SoilProfile(
10.0u"m" => SimpleSoil(por=0.30,sat=1.0,org=0.0),
);
## mid-winter temperature profile
-tempprofile = CryoGrid.Presets.SamoylovDefault.tempprofile
-forcings = loadforcings(CryoGrid.Presets.Forcings.Samoylov_ERA_obs_fitted_1979_2014_spinup_extended_2044);
-tempprofile = CryoGrid.Presets.SamoylovDefault.tempprofile
+tempprofile = CryoGrid.SamoylovDefault.tempprofile
+forcings = loadforcings(CryoGrid.Forcings.Samoylov_ERA_obs_fitted_1979_2014_spinup_extended_2044);
+tempprofile = CryoGrid.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..2bdc5d3e603f3bf227e5fa890f68378ec9f6c9f9 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),
@@ -15,14 +15,14 @@ soilprofile = SoilProfile(
3.0u"m" => SimpleSoil(por=0.50,sat=1.0,org=0.0),
10.0u"m" => SimpleSoil(por=0.30,sat=1.0,org=0.0),
);
-initT = initializer(:T, CryoGrid.Presets.SamoylovDefault.tempprofile)
+initT = initializer(:T, CryoGrid.SamoylovDefault.tempprofile)
initsat = initializer(:sat, 1.0)
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(),
@@ -39,8 +39,8 @@ strat = @Stratigraphy(
ground_layers...,
z_bot => Bottom(GeothermalHeatFlux(0.053u"J/s/m^2"))
);
-modelgrid = CryoGrid.Presets.DefaultGrid_5cm
-tile = Tile(strat, modelgrid, initT, initsat)
+modelgrid = CryoGrid.DefaultGrid_5cm
+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..6cc3f12c378449ba837f2ee01c268ac4a4ea3288 100644
--- a/examples/heat_sfcc_constantbc.jl
+++ b/examples/heat_sfcc_constantbc.jl
@@ -8,7 +8,7 @@
using CryoGrid
# Select default grid and initial temperature profile.
-grid = CryoGrid.Presets.DefaultGrid_10cm
+grid = CryoGrid.DefaultGrid_10cm
tempprofile = TemperatureProfile(
0.0u"m" => -20.0u"°C",
1000.0u"m" => 10.0u"°C",
@@ -31,12 +31,13 @@ Plots.plot(-2.0u"°C":0.01u"K":0.0u"°C", sfcc)
# specified.
heatop = Heat.Diffusion1D(:H)
initT = initializer(:T, tempprofile)
-tile = CryoGrid.Presets.SoilHeatTile(
+tile = CryoGrid.SoilHeatTile(
heatop,
## 10 W/m^2 in and out, i.e. net zero flux
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..a2e57485e8bc70c946abbbf3865b7c3aa0ca4cf2 100644
--- a/examples/heat_sfcc_richardseq_samoylov.jl
+++ b/examples/heat_sfcc_richardseq_samoylov.jl
@@ -5,9 +5,11 @@
# 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 +26,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.
@@ -38,8 +40,8 @@ strat = @Stratigraphy(
2.0u"m" => Ground(SimpleSoil(por=0.10,sat=1.0,org=0.0,freezecurve=sfcc); heat, water),
1000.0u"m" => Bottom(GeothermalHeatFlux(0.053u"W/m^2"))
);
-grid = CryoGrid.Presets.DefaultGrid_2cm
-tile = Tile(strat, grid, initT);
+grid = CryoGrid.DefaultGrid_2cm
+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_salt_constantbc.jl b/examples/heat_sfcc_salt_constantbc.jl
index bfc22d37b9970456d4b82f8eeef5ffbd46353250..e686b01359d1c4c45c69d2ea091d8a91d3f43c52 100644
--- a/examples/heat_sfcc_salt_constantbc.jl
+++ b/examples/heat_sfcc_salt_constantbc.jl
@@ -14,7 +14,7 @@ strat = @Stratigraphy(
0.0u"m" => SalineGround(SimpleSoil(freezecurve=sfcc), heat=HeatBalance(:T), salt=SaltMassBalance()),
1000.0u"m" => Bottom(GeothermalHeatFlux(0.053u"W/m^2"))
);
-modelgrid = CryoGrid.Presets.DefaultGrid_10cm
+modelgrid = CryoGrid.DefaultGrid_10cm
tile = Tile(strat, modelgrid, initT, initsalt, initpor)
tspan = (DateTime(1990,1,1),DateTime(2000,12,31))
u0, du0 = initialcondition!(tile, tspan)
diff --git a/examples/heat_sfcc_samoylov.jl b/examples/heat_sfcc_samoylov.jl
index 72bf5839cb25cf8cab69af2205e8f1f3d06bbe9c..add6edc0b85956aa87824933299e223342503a08 100644
--- a/examples/heat_sfcc_samoylov.jl
+++ b/examples/heat_sfcc_samoylov.jl
@@ -5,14 +5,15 @@
using CryoGrid
-# First we set up the soil heat model:
-forcings = loadforcings(CryoGrid.Presets.Forcings.Samoylov_ERA5_fitted_daily_1979_2020);
-grid = CryoGrid.Presets.DefaultGrid_5cm
-soilprofile, tempprofile = CryoGrid.Presets.SamoylovDefault
+# First we set up the soil heat model. Note that the default soil profile
+# for Samoylov already includes the appropriate freeze curves.
+forcings = loadforcings(CryoGrid.Forcings.Samoylov_ERA5_fitted_daily_1979_2020);
+grid = CryoGrid.DefaultGrid_5cm
+soilprofile, tempprofile = CryoGrid.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.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..9bc97ad52d69df8ec396d7514cf475d0698b3261 100644
--- a/examples/heat_simple_autodiff_grad.jl
+++ b/examples/heat_simple_autodiff_grad.jl
@@ -3,52 +3,61 @@
# two parameters (summer and winter n-factors) using forward-mode automatic simulation.
#
# TODO: add more detail/background
+using CryoGrid
# 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);
-soilprofile, tempprofile = CryoGrid.Presets.SamoylovDefault
-grid = CryoGrid.Presets.DefaultGrid_5cm
+forcings = loadforcings(CryoGrid.Forcings.Samoylov_ERA_obs_fitted_1979_2014_spinup_extended_2044);
+soilprofile, tempprofile = CryoGrid.SamoylovDefault
+grid = CryoGrid.DefaultGrid_5cm
initT = initializer(:T, tempprofile)
-tile = CryoGrid.Presets.SoilHeatTile(
+tile = CryoGrid.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
)
-tspan = (DateTime(2010,10,1),DateTime(2010,10,2))
+tspan = (DateTime(2010,9,1),DateTime(2011,10,1))
u0, du0 = @time initialcondition!(tile, tspan);
-# Collect model parameters
-p = CryoGrid.parameters(tile)
+# We can retrieve the parameters of the system from `tile`:
+para = CryoGrid.parameters(tile)
# Create the `CryoGridProblem`.
-prob = CryoGridProblem(tile, u0, tspan, p, saveat=3600.0);
+prob = CryoGridProblem(tile, u0, tspan, saveat=3600.0)
# Solve the forward problem with default parameter settings:
sol = @time solve(prob)
out = CryoGridOutput(sol)
-# ForwardDiff provides tools for forward-mode automatic differentiation.
+# Import relevant packages for automatic differentiation.
using ForwardDiff
using SciMLSensitivity
using Zygote
# Define a "loss" function; here we'll just take the mean over the final temperature field.
+using OrdinaryDiffEq
using Statistics
function loss(prob::CryoGridProblem, p)
- # local u0, _ = initialcondition!(tile, tspan, p)
- # local prob = CryoGridProblem(tile, u0, tspan, p, saveat=24*3600.0)
newprob = remake(prob, p=p)
+ # Here we specify the sensitivity algorithm. Note that this is only
+ # necessary for reverse-mode autodiff with Zygote.
+ # autojacvec = true uses ForwardDiff to calculate the jacobian;
+ # enabling checkpointing (theoretically) reduces the memory cost of the backwards pass.
sensealg = InterpolatingAdjoint(autojacvec=true, checkpointing=true)
newsol = solve(newprob, Euler(), dt=300.0, sensealg=sensealg);
newout = CryoGridOutput(newsol)
return mean(ustrip.(newout.T[:,end]))
end
-# Compute gradient with forward diff:
-pvec = vec(p)
-grad = @time ForwardDiff.gradient(páµ¢ -> loss(prob, páµ¢), pvec)
-grad = @time Zygote.gradient(páµ¢ -> loss(prob, páµ¢), pvec)
+# Compute gradient with forward-mode autodiff:
+pvec = vec(prob.p)
+fd_grad = @time ForwardDiff.gradient(páµ¢ -> loss(prob, páµ¢), pvec)
+
+# We can also try with reverse-mode autodiff. This is generally slower for smaller numbers
+# of parmaeters (<100) but could be worthwhile for model configurations with high-dimensional
+# parameterizations.
+zy_grad = @time Zygote.gradient(páµ¢ -> loss(prob, páµ¢), pvec)
+@assert maximum(abs.(fd_grad .- zy_grad)) .< 1e-6 "Forward and reverse gradients don't match!"
diff --git a/src/CryoGrid.jl b/src/CryoGrid.jl
index 97c9a9944a1dbb63027ea5824ee4907d09d247d9..9c227ea0537a09d0314798f675d17755d5d95387 100755
--- a/src/CryoGrid.jl
+++ b/src/CryoGrid.jl
@@ -1,12 +1,34 @@
module CryoGrid
-global CRYOGRID_DEBUG = haskey(ENV,"CG_DEBUG") && ENV["CG_DEBUG"] == "true"
-function debug(debug::Bool)
- global CRYOGRID_DEBUG = debug
+global DEBUG = haskey(ENV,"CG_DEBUG") && ENV["CG_DEBUG"] == "true"
+@deprecate debug(x) debug!(x)
+
+"""
+ debug!(debug::Bool)
+
+Enable or disable global debug mode for CryoGrid. Debug mode disables certain optimizations (e.g. loop vectorizations)
+which sometimes ineterfere with the debugger. It also enables additional numerical stability checks.
+"""
+function debug!(debug::Bool)
+ global DEBUG = debug
# disable loop vectorization in debug mode
Numerics.turbo(!debug)
- CRYOGRID_DEBUG && @warn "Debug mode enabled! Some performance features such as loop vectorization are now turned off by default."
- return CRYOGRID_DEBUG
+ DEBUG && @warn "Debug mode enabled! Some performance features such as loop vectorization are now turned off by default."
+ return DEBUG
+end
+
+global AUTOPARA = haskey(ENV,"CG_AUTOPARA") && ENV["CG_AUTOPARA"] == "true"
+
+"""
+ autoparam!(parameterize::Bool)
+
+Enable or disable automatic parameterization mode. When enabled, model parameters in all layer/process types will be initialized
+as free `Param`s rather than `FixedParam`s.
+"""
+function autoparam!(parameterize::Bool)
+ global AUTOPARA = parameterize
+ AUTOPARA && @warn "Automatic parameterization enabled!" || @warn "Automatic parameterization disabled!"
+ return AUTOPARA
end
using Adapt
@@ -36,6 +58,8 @@ import Interpolations
@reexport using Unitful
@reexport using UnPack
+using SciMLBase: intervals # resolve name collision
+
export Interpolations
# Common types and methods
@@ -54,8 +78,8 @@ include("variables.jl")
export BCKind
include("traits.jl")
-export initialcondition!, computediagnostic!, interact!, interactmaybe!, computefluxes!, resetfluxes!, diagnosticstep!
-export variables, processes, initializers, timestep, isactive, caninteract
+export initialcondition!, computediagnostic!, interact!, interactmaybe!, computeprognostic!, resetfluxes!, diagnosticstep!
+export variables, processes, initializers, timestep, isactive, caninteract, param
export boundaryflux, boundaryvalue, criterion, criterion!, trigger!
include("methods.jl")
@@ -69,6 +93,8 @@ export DiscretizationStrategy, AutoGrid, PresetGrid, LinearSpacing, Grid, cells,
include("Numerics/Numerics.jl")
using .Numerics
+include("default_grids.jl")
+
export initializer
include("initializers.jl")
@@ -77,6 +103,7 @@ include("IO/InputOutput.jl")
include("Tiles/Tiles.jl")
@reexport using .Tiles
+using .Tiles: layers # resolve name collision
parameters = Tiles.parameters
export ConstantBC, PeriodicBC, ConstantValue, PeriodicValue, ConstantFlux, PeriodicFlux
@@ -97,8 +124,22 @@ include("problem.jl")
export CGEuler, CryoGridIntegrator, CryoGridSolution
include("Solvers/Solvers.jl")
-# Presets
-include("Presets/Presets.jl")
+# Built-in model definitions
+export SoilHeatTile, SamoylovDefault
+include("models.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
diff --git a/src/Diagnostics/spinup.jl b/src/Diagnostics/spinup.jl
index 3049132f5c755322bd5c1effcde68d64a10a9255..2da0fbb5ef788870d2b3fc973610e0c8ebb87211 100644
--- a/src/Diagnostics/spinup.jl
+++ b/src/Diagnostics/spinup.jl
@@ -5,9 +5,9 @@ Implements a simple, iterative spin-up procedure.
Runs the model specified by `setup` over `tspan` until the profile mean up to `maxdepth` over the whole time span changes only within the given tolerance `tol`.
Returns the `ODESolution` generated by the final iteration.
"""
-function spinup(setup::Tile, tspan::NTuple{2,DateTime}, p, tol, layername; maxdepth=100u"m", maxiter=1000, saveat=3*3600.0, solver=CGEuler(), dt=60.0, solve_args...)
- u0, du0 = initialcondition!(tile, tspan, p)
- prob = CryoGridProblem(setup, u0, tspan, p)
+function spinup(setup::Tile, tspan::NTuple{2,DateTime}, tol, layername; maxdepth=100u"m", maxiter=1000, saveat=3*3600.0, solver=CGEuler(), dt=60.0, solve_args...)
+ u0, du0 = initialcondition!(tile, tspan)
+ prob = CryoGridProblem(setup, u0, tspan)
@info "Running initial solve ..."
sol = solve(prob, solver, dt=dt, saveat=saveat, solve_args...)
out = CryoGridOutput(sol)
diff --git a/src/IO/InputOutput.jl b/src/IO/InputOutput.jl
index ec01a181d3d6eaa31b66985e688ce6b2b84d0209..aa76df7ba5fc0d3d1ad8c173cc4efd259d9c6074 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
@@ -59,13 +60,12 @@ include("ioutils.jl")
export CryoGridParams
include("params/params.jl")
-export ParamsJSON, ParamsYAML
-include("params/params_loaders.jl")
+export Input, InputProvider, InputFunctionProvider
+export inputs
+include("input.jl")
-export Forcings, Forcing, ConstantForcing, TransformedForcing, TimeVaryingForcing, InterpolatedForcing
-export TemperatureForcing, WindForcing, HumidityForcing, EnergyFluxForcing, PressureForcing, VelocityForcing # aliases
-export ForcingFormat, ForcingFormatJSON, ForcingFormatNCD
-export loadforcings, time_derivative_forcing
+export ForcingFormat, ForcingFormatJSON, ForcingFormatNCD, Interpolated1D
+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..0a19ff877e2597d3dcd632f24d2bc73a06223015 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,11 +37,11 @@ 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)
+function loadforcings(format::ForcingFormatJSON{2}, 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]...)
@@ -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..0d3cde62d4cb57846aac157f945a4f92f9da1a8d
--- /dev/null
+++ b/src/IO/forcings/providers.jl
@@ -0,0 +1,58 @@
+"""
+ 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 DimensionalData.rebuild(data::Interpolated1D, default_interp_mode=Linear(); stripunits=true, interp_modes=(;), datasets...)
+ newdata = DimStack((; datasets...))
+ return Interpolated1D(newdata, default_interp_mode; stripunits, interp_modes...)
+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/IO/params/param_types.jl b/src/IO/params/param_types.jl
new file mode 100644
index 0000000000000000000000000000000000000000..f1fb4e6f12e8da6c109f4d0181b8a0c1d04263d4
--- /dev/null
+++ b/src/IO/params/param_types.jl
@@ -0,0 +1,25 @@
+"""
+ FixedParam(p::NamedTuple)
+ FixedParam(; kw...)
+ FixedParam(val)
+
+Subtype of `AbstractParam` that rerpesents a "fixed" parameter which
+should not be included in the set of free parameters for the model.
+"""
+struct FixedParam{T<:Number} <: AbstractParam{T}
+ parent::NamedTuple
+ function FixedParam{T}(nt::NamedTuple) where {T<:Number}
+ @assert :val ∈ keys(nt)
+ new{T}(nt)
+ end
+end
+
+function FixedParam(nt::NT) where {NT<:NamedTuple}
+ FixedParam{typeof(nt.val)}(nt)
+end
+
+FixedParam(val; kwargs...) = FixedParam((; val, kwargs...))
+
+Base.parent(param::FixedParam) = getfield(param, :parent)
+
+ModelParameters.rebuild(param::FixedParam, nt::NamedTuple) = FixedParam(nt)
diff --git a/src/IO/params/parameterizations.jl b/src/IO/params/parameterizations.jl
index 83878c6cbf4e6b5fe3d2c6f950f134048e78af67..8a0ee5c5ea59f0426ca3ac641c79094e2b7b01ca 100644
--- a/src/IO/params/parameterizations.jl
+++ b/src/IO/params/parameterizations.jl
@@ -127,11 +127,6 @@ function binindex(values::Tuple, st, en, x)
end
end
-CryoGrid.parameters(pw::PiecewiseLinear) = (;
- initialvalue = CryoGrid.parameters(pw.initialvalue),
- knots = CryoGrid.parameters(pw.knots),
-)
-
# Transformed parameterization
"""
diff --git a/src/IO/params/params.jl b/src/IO/params/params.jl
index cd2e93f9ccb44b8f66c8db74b84bebcfd2666597..227c38854a4f7f9080fc8e55a3a1256eea8ffc11 100644
--- a/src/IO/params/params.jl
+++ b/src/IO/params/params.jl
@@ -1,3 +1,6 @@
+export FixedParam
+include("param_types.jl")
+
"""
CryoGridParams{T,TM} <: DenseArray{T,1}
@@ -6,7 +9,7 @@ type directly in math or linear algebra operations but rather to use `Base.value
of parameter values.
"""
struct CryoGridParams{T,TM} <: DenseArray{T,1}
- obj::TM # param obj
+ obj::TM # parameter "model"
CryoGridParams(m::AbstractModel) = new{eltype(m[:val]),typeof(m)}(m)
end
@@ -66,7 +69,7 @@ function Base.show(io::IO, ::MIME"text/plain", ps::CryoGridParams{T}) where T
ModelParameters.printparams(io, ps.obj)
end
-paramname(p::Param, component::Type{T}, fieldname) where {T} = fieldname
+paramname(p::AbstractParam, component::Type{T}, fieldname) where {T} = fieldname
Tables.columns(ps::CryoGridParams) = Tables.columns(ps.obj)
Tables.rows(ps::CryoGridParams) = Tables.rows(ps.obj)
@@ -91,3 +94,9 @@ function _setparafields(m::Model)
newparent = Flatten.reconstruct(parent(m), updated_parameterizations, CryoGrid.Parameterization)
return Model(newparent)
end
+
+export ParamsJSON, ParamsYAML
+include("params_loaders.jl")
+
+export PiecewiseConstant, PiecewiseLinear, LinearTrend, Transformed
+include("parameterizations.jl")
diff --git a/src/Physics/Heat/heat_bc.jl b/src/Physics/Heat/heat_bc.jl
index a09ec57877a632da7b53a989ad231e605ba0ccb3..0f05f1d59a322aadd1e64a5ff811840f51c296aa 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 6a237d616a506920d01035c50434300c7a1b03a6..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.
"""
@@ -104,29 +138,17 @@ function CryoGrid.interact!(sub1::SubSurface, ::HeatBalance, sub2::SubSurface, :
return nothing
end
-function CryoGrid.computefluxes!(::SubSurface, ::HeatBalance{<:EnthalpyBased}, state)
+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.computefluxes!(sub::SubSurface, ::HeatBalance{<:TemperatureBased}, state)
+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 dded8b6fbc456509b4459b8d3c42386fbded221e..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) = (
@@ -82,8 +102,8 @@ function CryoGrid.interact!(sub1::SubSurface, ::HeatBalanceImplicit, sub2::SubSu
return nothing
end
-# do nothing in computefluxes!
-CryoGrid.computefluxes!(::SubSurface, ::HeatBalanceImplicit, state) = nothing
+# do nothing in computeprognostic!
+CryoGrid.computeprognostic!(::SubSurface, ::HeatBalanceImplicit, state) = nothing
function CryoGrid.resetfluxes!(sub::SubSurface, heat::HeatBalanceImplicit, state)
@inbounds for i in 1:length(state.H)
@@ -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 5c513182411145961fb8d90c2f93e4fcab1efe60..abdc2cfc772c5bc5096e4337436b25aedfe60694 100644
--- a/src/Physics/Heat/heat_water.jl
+++ b/src/Physics/Heat/heat_water.jl
@@ -13,18 +13,15 @@ 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)
- 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;
+function water_energy_advection!(jH, jw, T, cw::Real, L::Real)
+ @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
@@ -88,14 +85,21 @@ function CryoGrid.interact!(
end
end
-# Flux calculation
-function CryoGrid.computefluxes!(sub::SubSurface, ps::Coupled(WaterBalance, HeatBalance), state)
+function CryoGrid.computeprognostic!(sub::SubSurface, ps::Coupled(WaterBalance, HeatBalance), state)
water, heat = ps
- CryoGrid.computefluxes!(sub, water, state)
+ CryoGrid.computeprognostic!(sub, water, state)
if heat.advection
- water_energy_advection!(sub, ps, state)
+ cw = heatcapacitywater(sub, state)
+ L = heat.prop.L
+ water_energy_advection!(state.jH, state.jw, state.T, cw, L)
end
- CryoGrid.computefluxes!(sub, heat, state)
+ CryoGrid.computeprognostic!(sub, heat, state)
+end
+
+function CryoGrid.computeprognostic!(sub::SubSurface, ps::Coupled(WaterBalance{NoFlow}, HeatBalance), state)
+ water, heat = ps
+ CryoGrid.computeprognostic!(sub, water, state)
+ CryoGrid.computeprognostic!(sub, heat, state)
end
function heatcapacitywater(sub::SubSurface, state)
diff --git a/src/Physics/Hydrology/water_ET.jl b/src/Physics/Hydrology/water_ET.jl
index 2f5e5e172e23932c1fc44c47c5f74f6e8f5022f9..40c99376c415a287ea5eb8ed23ab89fcb0aaa713 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_balance.jl b/src/Physics/Hydrology/water_balance.jl
index b87ff3ef4ba37a0c7264b86735f6649f214a58a6..0aea1b9ec3c44bfe39a4f73faaab3cf55f5d4eb0 100644
--- a/src/Physics/Hydrology/water_balance.jl
+++ b/src/Physics/Hydrology/water_balance.jl
@@ -219,7 +219,7 @@ function CryoGrid.computediagnostic!(sub::SubSurface, water::WaterBalance, state
evapotranspiration!(sub, water, state)
end
-function CryoGrid.computefluxes!(sub::SubSurface, water::WaterBalance, state)
+function CryoGrid.computeprognostic!(sub::SubSurface, water::WaterBalance, state)
wateradvection!(sub, water, state)
waterdiffusion!(sub, water, state)
balancefluxes!(sub, water, state)
@@ -255,7 +255,7 @@ function CryoGrid.initialcondition!(sub::SubSurface, water::WaterBalance{NoFlow}
end
end
-CryoGrid.computefluxes!(::SubSurface, ::WaterBalance{NoFlow}, state) = nothing
+CryoGrid.computeprognostic!(::SubSurface, ::WaterBalance{NoFlow}, state) = nothing
CryoGrid.resetfluxes!(::SubSurface, ::WaterBalance{NoFlow}, state) = nothing
diff --git a/src/Physics/Hydrology/water_types.jl b/src/Physics/Hydrology/water_types.jl
index b9cb5666d352a5a06b70c972f5ed93df2d09cfac..04c718bd38d36dc86fdd61f02b4d21e3a66fe5d5 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/Lakes/lake_simple.jl b/src/Physics/Lakes/lake_simple.jl
index 643e887edf75f1ce52350092f8f3b3d7ee172a04..7a1ae48dedf522bf98d473bcbbcad49cc2a09e80 100644
--- a/src/Physics/Lakes/lake_simple.jl
+++ b/src/Physics/Lakes/lake_simple.jl
@@ -110,7 +110,7 @@ function CryoGrid.computediagnostic!(sub::Lake, heat::HeatBalance{<:Heat.Diffusi
return nothing
end
-function CryoGrid.computefluxes!(::Lake, ::HeatBalance{<:Heat.Diffusion1D{:H}}, state)
+function CryoGrid.computeprognostic!(::Lake, ::HeatBalance{<:Heat.Diffusion1D{:H}}, state)
Δk = Δ(state.grids.k) # cell sizes
ΔT = Δ(state.grids.T) # midpoint distances
# compute internal fluxes and non-linear diffusion assuming boundary fluxes have been set
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/Salt/salt_diffusion.jl b/src/Physics/Salt/salt_diffusion.jl
index f6175a76e6e2ca37356619d5e28e3d1f6f4f97c8..1401acc4c1109365afef7eba8af88967e5f19f36 100644
--- a/src/Physics/Salt/salt_diffusion.jl
+++ b/src/Physics/Salt/salt_diffusion.jl
@@ -121,7 +121,7 @@ function CryoGrid.timestep(::SalineGround, salt::SaltMassBalance{T,<:CryoGrid.CF
return dtmax
end
-function CryoGrid.computefluxes!(
+function CryoGrid.computeprognostic!(
::SalineGround,
ps::CoupledHeatSalt{THeat},
state
diff --git a/src/Physics/Salt/salt_types.jl b/src/Physics/Salt/salt_types.jl
index 88176b458ab207304489e2368c872170b52871f4..c7df4b8afd48fc99458a54b11f44cdf4b96df6e5 100644
--- a/src/Physics/Salt/salt_types.jl
+++ b/src/Physics/Salt/salt_types.jl
@@ -1,6 +1,6 @@
SaltProperties(
- Ï„ = Param(1.5), # Turtuosity
- dₛ₀ = Param(8.0e-10, units=u"m^2/s"), # salt diffusion coefficient
+ Ï„ = param(1.5), # Turtuosity
+ dₛ₀ = param(8.0e-10, units=u"m^2/s"), # salt diffusion coefficient
) = (; τ, dₛ₀)
abstract type SaltOperator end
diff --git a/src/Physics/Snow/snow_bulk.jl b/src/Physics/Snow/snow_bulk.jl
index 230a7efa4094c219f0d6cfd475454278c202bd65..35a7cbdff37d45f87720db5f60ed9362482e2291 100644
--- a/src/Physics/Snow/snow_bulk.jl
+++ b/src/Physics/Snow/snow_bulk.jl
@@ -167,16 +167,16 @@ function CryoGrid.trigger!(
end
# mass balance fluxes are handled in interact! for bulk snowpack
-CryoGrid.computefluxes!(snow::BulkSnowpack, mass::SnowMassBalance, state) = nothing
+CryoGrid.computeprognostic!(snow::BulkSnowpack, mass::SnowMassBalance, state) = nothing
-# computefluxes! for free water, enthalpy based HeatBalance on bulk snow layer
-function CryoGrid.computefluxes!(
+# computeprognostic! for free water, enthalpy based HeatBalance on bulk snow layer
+function CryoGrid.computeprognostic!(
snow::BulkSnowpack,
ps::CoupledSnowWaterHeat{TM,TW,TH},
state
) where {TM,TW,TH<:HeatBalance{<:EnthalpyBased}}
mass, water, heat = ps
- computefluxes!(snow, mass, state)
+ computeprognostic!(snow, mass, state)
dsn = getscalar(state.dsn)
if dsn < snow.para.thresh
# set fluxes to zero if there is no snow
@@ -185,8 +185,8 @@ function CryoGrid.computefluxes!(
state.dsat .= 0.0
end
else
- # otherwise call computefluxes! for coupled water/heat
- computefluxes!(snow, Coupled(water, heat), state)
+ # otherwise call computeprognostic! for coupled water/heat
+ computeprognostic!(snow, Coupled(water, heat), state)
end
return nothing
end
diff --git a/src/Physics/Snow/snow_lite.jl b/src/Physics/Snow/snow_lite.jl
index 597194d3b6bd3de1fc8c098c2a955cc075cc61e1..1dbbeb700772d977bd16e72c5615a427e32c468b 100644
--- a/src/Physics/Snow/snow_lite.jl
+++ b/src/Physics/Snow/snow_lite.jl
@@ -40,7 +40,7 @@ CryoGrid.variables(::LiteSnowpack, ::SnowMassBalance) = (
CryoGrid.computediagnostic!(::LiteSnowpack, ::SnowMassBalance, state) = nothing
-CryoGrid.computefluxes!(::LiteSnowpack, ::SnowMassBalance, state) = nothing
+CryoGrid.computeprognostic!(::LiteSnowpack, ::SnowMassBalance, state) = nothing
function CryoGrid.interact!(
top::Top,
diff --git a/src/Physics/Snow/snowpack.jl b/src/Physics/Snow/snowpack.jl
index 69d6705044e2c79dd5e91a98c0af3951c40637fc..4644619cf28796433e1be0f01465b73f15c1e5c0 100644
--- a/src/Physics/Snow/snowpack.jl
+++ b/src/Physics/Snow/snowpack.jl
@@ -120,14 +120,14 @@ function CryoGrid.computediagnostic!(
updategrid!(state.grid, z0, state.dsn)
end
-function CryoGrid.computefluxes!(
+function CryoGrid.computeprognostic!(
snow::Snowpack,
state,
)
mass, water, heat = processes(snow)
- computefluxes!(snow, mass, state)
- computefluxes!(snow, water, state)
- computefluxes!(snow, heat, state)
+ computeprognostic!(snow, mass, state)
+ computeprognostic!(snow, water, state)
+ computeprognostic!(snow, heat, state)
end
# Special overrides for heat timestep control on snow layer
diff --git a/src/Physics/Soils/para/simple.jl b/src/Physics/Soils/para/simple.jl
index 2052ae5d4295ed51a6227e1085f058b661ba2068..aaa13fdfcfb4fda0f44c64171e1fe92cddfda822 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, units=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, units=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/Sources/Sources.jl b/src/Physics/Sources/Sources.jl
index 20baf9b080372adf8020d96ce560319ad963dff1..53b58d4dd1917fcf89c087c8b67ad2a7863f30dc 100644
--- a/src/Physics/Sources/Sources.jl
+++ b/src/Physics/Sources/Sources.jl
@@ -1,7 +1,7 @@
module Sources
import CryoGrid: SubSurfaceProcess, SubSurface
-import CryoGrid: computediagnostic!, initialcondition!, interact!, computefluxes!, variables
+import CryoGrid: computediagnostic!, initialcondition!, interact!, computeprognostic!, variables
using ..Heat
using CryoGrid.Numerics
@@ -52,7 +52,7 @@ ConstructionBase.constructorof(::Type{<:Source{P}}) where {P} = (term, aux) -> S
(p::Periodic)(t) = p.amp*sin(2Ï€*p.freq*t - p.shift) + p.level
# HeatBalance sources
-computefluxes!(::SubSurface, s::Source{<:HeatBalance,<:Constant}, state) = @inbounds @. state.dH += s.term.Sâ‚€
-computefluxes!(::SubSurface, src::Source{<:HeatBalance,<:Periodic}, state) = @inbounds @. state.dH += src.term(state.t)
+computeprognostic!(::SubSurface, s::Source{<:HeatBalance,<:Constant}, state) = @inbounds @. state.dH += s.term.Sâ‚€
+computeprognostic!(::SubSurface, src::Source{<:HeatBalance,<:Periodic}, state) = @inbounds @. state.dH += src.term(state.t)
end
\ No newline at end of file
diff --git a/src/Physics/Surface/SEB/seb.jl b/src/Physics/Surface/SEB/seb.jl
index 6bc234d58c68b69a52b20af88e1256c2c176010d..b5245cdb2e20edef8f3b91fd519ebd8f31f49e2f 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..972179688cfe6f8cccbfc129ea896815eff46a8c 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}
+ SEBInputs{TT1,TT2,TR,TP,TQ,TW,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
+Base.@kwdef struct SEBInputs{TT1,TT2,TR,TP,TQ,TW,Tsurf}
+ Ts::TT1 = 0.0u"°C"
+ 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/coupled_sweb.jl b/src/Physics/Surface/SWB/coupled_sweb.jl
index 19753844faa3ab37972afad61e3934e80e2c0bda..53a5de7929f10f3ca873209d6af5eef85fccadd8 100644
--- a/src/Physics/Surface/SWB/coupled_sweb.jl
+++ b/src/Physics/Surface/SWB/coupled_sweb.jl
@@ -13,11 +13,11 @@ CryoGrid.variables(top::Top, bc::SurfaceWaterEnergyBalance) = (
Prognostic(:ET, Scalar, u"m^3"),
)
-function CryoGrid.computefluxes!(top::Top, bc::SurfaceWaterEnergyBalance, state)
+function CryoGrid.computeprognostic!(top::Top, bc::SurfaceWaterEnergyBalance, state)
swb, seb = bc
- computefluxes!(top, swb, state)
+ computeprognostic!(top, swb, state)
ET!(top, swb, state)
- computefluxes!(top, seb, state)
+ computeprognostic!(top, seb, state)
end
function CryoGrid.interact!(top::Top, bc::SurfaceWaterEnergyBalance, sub::SubSurface, ps::Coupled(WaterBalance, HeatBalance), stop, ssub)
diff --git a/src/Physics/Surface/SWB/swb.jl b/src/Physics/Surface/SWB/swb.jl
index c395e12960807a03152d66fed8fd824e21a00391..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,11 +40,11 @@ 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.computefluxes!(top::Top, swb::SurfaceWaterBalance, state)
+function CryoGrid.computeprognostic!(top::Top, swb::SurfaceWaterBalance, state)
runoff!(top, swb, state)
end
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
deleted file mode 100644
index d027a53f4975a512254f23979f75e4813f361de3..0000000000000000000000000000000000000000
--- a/src/Presets/Presets.jl
+++ /dev/null
@@ -1,69 +0,0 @@
-"""
-Pre-built CryoGrid configurations for rapid prototyping.
-"""
-module Presets
-
-using CryoGrid
-using CryoGrid.InputOutput: Resource
-using CryoGrid.Numerics
-using CryoGrid.Utils
-
-# physics modules
-using CryoGrid.Heat
-using CryoGrid.Hydrology
-using CryoGrid.Soils
-
-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...)
- 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...)
-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")
-)
-
-const SamoylovDefault = (
- soilprofile = SoilProfile(
- 0.0u"m" => SimpleSoil(por=0.80,sat=1.0,org=0.75,freezecurve=PainterKarra(swrc=VanGenuchten("organic"))),
- 0.1u"m" => SimpleSoil(por=0.80,sat=1.0,org=0.25,freezecurve=PainterKarra(swrc=VanGenuchten("organic"))),
- 0.4u"m" => SimpleSoil(por=0.55,sat=1.0,org=0.25,freezecurve=PainterKarra(swrc=VanGenuchten("sandy loam"))),
- 3.0u"m" => SimpleSoil(por=0.50,sat=1.0,org=0.0,freezecurve=PainterKarra(swrc=VanGenuchten("sandy loam"))),
- 10.0u"m" => SimpleSoil(por=0.30,sat=1.0,org=0.0,freezecurve=PainterKarra(swrc=VanGenuchten("sandy loam"))),
- 30.0u"m" => SimpleSoil(por=0.05,sat=1.0,org=0.0,freezecurve=PainterKarra(swrc=VanGenuchten("sand"))),
- ),
- tempprofile = TemperatureProfile(
- 0.0u"m" => -1.0u"°C",
- 0.2u"m" => -1.0u"°C",
- 2.0u"m" => -1.0u"°C",
- 5.0u"m" => -3.0u"°C",
- 9.5u"m" => -6.0u"°C",
- 25.0u"m" => -9.0u"°C",
- 100.0u"m" => -9.0u"°C",
- 1000.0u"m" => 10.2u"°C",
- )
-)
-
-end
\ No newline at end of file
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..5e9afe2d97de1aafaaf3df2672a5d26fccf7dd92 100644
--- a/src/Solvers/basic_solvers.jl
+++ b/src/Solvers/basic_solvers.jl
@@ -19,7 +19,7 @@ struct CGEulerCache{Tu} <: SciMLBase.DECache
du::Tu
end
-function DiffEqBase.__init(prob::CryoGridProblem, alg::CGEuler, args...; dt=60.0, saveat=3*3600.0, kwargs...)
+function DiffEqBase.__init(prob::CryoGridProblem, alg::CGEuler, args...; dt=60.0, saveat=nothing, kwargs...)
tile = Tile(prob.f)
u0 = copy(prob.u0)
du0 = zero(u0)
@@ -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))
@@ -42,7 +42,11 @@ function DiffEqBase.__init(prob::CryoGridProblem, alg::CGEuler, args...; dt=60.0
similar(prob.u0),
)
p = isnothing(prob.p) ? prob.p : collect(prob.p)
- saveat = isa(saveat, Number) ? collect(prob.tspan[1]:saveat:prob.tspan[2]) : saveat
+ if isnothing(saveat)
+ saveat = prob.saveat
+ elseif isa(saveat, Number)
+ saveat = collect(prob.tspan[1]:saveat:prob.tspan[2])
+ end
opts = CryoGridIntegratorOptions(; saveat=expandtstep(saveat, prob.tspan), kwargs...)
return CryoGridIntegrator(alg, cache, opts, sol, copy(u0), p, t0*one(eltype(u0)), dt*one(eltype(u0)), 1, 1)
end
@@ -54,7 +58,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/stratigraphy.jl b/src/Tiles/stratigraphy.jl
index a54157c958ba206314aa5741ec19335cd5ef6670..cf1861c34cd8bb94fe91ca5b6ef6c6b42054ca62 100644
--- a/src/Tiles/stratigraphy.jl
+++ b/src/Tiles/stratigraphy.jl
@@ -33,7 +33,7 @@ struct Stratigraphy{N,TLayers<:NamedTuple,TBoundaries}
@nospecialize(sub::Tuple{Vararg{Pair{<:Number}}}),
@nospecialize(bot::Pair{<:Number,<:Bottom})
)
- updateparam(p::Param) = Param(merge(parent(p), (layer=:strat,)))
+ updateparam(p::paraType) where {paraType<:AbstractParam} = paraType(merge(parent(p), (layer=:strat,)))
updateparam(x) = x
# check subsurface layers
@assert length(sub) > 0 "At least one subsurface layer must be specified"
@@ -93,7 +93,9 @@ Base.NamedTuple(strat::Stratigraphy) = layers(strat)
function Base.show(io::IO, ::MIME"text/plain", strat::Stratigraphy)
print(io, "Stratigraphy:\n")
for (k,v,b) in zip(keys(strat), values(strat), boundaries(strat))
- print(repeat(" ", Base.indent_width), "$b: $k :: $(nameof(typeof(v)))\n")
+ procs = isa(processes(v), Tuple) ? processes(v) : (processes(v),)
+ procnames = map(p -> string(nameof(typeof(p))), procs)
+ print(repeat(" ", Base.indent_width), "$(b*u"m"): $k :: $(nameof(typeof(v))) with processes $(join(procnames, ", "))\n")
end
end
@@ -219,9 +221,9 @@ that `getproperty` would return the appropriate state object for the i'th layer
return expr
end
-function CryoGrid.computefluxes!(strat::Stratigraphy, state)
+function CryoGrid.computeprognostic!(strat::Stratigraphy, state)
fastiterate(namedlayers(strat)) do named_layer
- CryoGrid.computefluxes!(named_layer.val, getproperty(state, nameof(named_layer)))
+ CryoGrid.computeprognostic!(named_layer.val, getproperty(state, nameof(named_layer)))
end
end
diff --git a/src/Tiles/tile.jl b/src/Tiles/tile.jl
index 2b8975591a83e5e10c8747b17bb5a777ac4c3ae4..797db01d42502a88a8c9f17f8462c8ef6e6956eb 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,21 +18,29 @@ 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")
+
+function Base.show(io::IO, mime::MIME"text/plain", tile::Tile{TStrat,TGrid,TStates,TInits,TEvents,TInputs,iip}) where {TStrat,TGrid,TStates,TInits,TEvents,TInputs,iip}
+ println(io, "Tile (iip=$iip) with layers $(keys(tile.strat))")
+ print(" "); show(io, mime, tile.grid); println();
+ print(" "); show(io, mime, tile.strat);
+ println(" metadata: $(tile.metadata)")
+end
"""
Tile(
@@ -50,6 +59,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 +92,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 +117,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
"""
- computefluxes!(_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
@@ -117,21 +130,21 @@ in sequence, i.e for each layer 1 <= i <= N:
```julia
computediagnostic!(layer[i], ...)
interact!(layer[i], ..., layer[i+1], ...)
-computefluxes!(layer[i], ...)
+computeprognostic!(layer[i], ...)
```
"""
-function computefluxes!(
- _tile::Tile{TStrat,TGrid,TStates,TInits,TEvents,true},
+function computeprognostic!(
+ _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)
@@ -139,8 +152,8 @@ function computefluxes!(
checkstate!(tile, state, u, du, :computediagnostic!)
CryoGrid.interact!(strat, state)
checkstate!(tile, state, u, du, :interact!)
- CryoGrid.computefluxes!(strat, state)
- checkstate!(tile, state, u, du, :computefluxes!)
+ CryoGrid.computeprognostic!(strat, state)
+ checkstate!(tile, state, u, du, :computeprognostic!)
return nothing
end
@@ -154,7 +167,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,13 +189,13 @@ 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)
CryoGrid.initialcondition!(strat, state, tile.inits...)
# evaluate initial time derivative
- computefluxes!(tile, du, u, p, t0)
+ computeprognostic!(tile, du, u, p, t0)
return u, du
end
@@ -269,28 +282,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,39 +329,49 @@ 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 = Utils.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 = Utils.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
+ if CryoGrid.DEBUG
@inbounds for i in eachindex(u)
if !isfinite(u[i])
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 +402,17 @@ function _initstatevars(@nospecialize(strat::Stratigraphy), @nospecialize(grid::
return states
end
+function _validate_inputs(@nospecialize(tile::Tile), inputprovider::InputProvider)
+ IgnoreTypes = Utils.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 Utils.ignored_types(::Tile{TStrat,TGrid,TStates}) where {TStrat,TGrid,TStates} = Union{TGrid,TStates,TileData,Unitful.Quantity,Numerics.ForwardDiff.Dual}
+
# ===================================================================== #
diff --git a/src/Tiles/tile_base.jl b/src/Tiles/tile_base.jl
index 8340d52affaaa178b0ac422ec65ccd3a925a085c..be02dd9149de9c3c98d1a0b9ee9aa2049bf3bce2 100644
--- a/src/Tiles/tile_base.jl
+++ b/src/Tiles/tile_base.jl
@@ -16,9 +16,9 @@ abstract type AbstractTile{iip} end
(tile::AbstractTile{true})(du,u,p,t,dt=1.0)
(tile::AbstractTile{false})(u,p,t,dt=1.0)
-Invokes `computefluxes!` on `tile` to compute the time derivative du/dt.
+Invokes `computeprognostic!` on `tile` to compute the time derivative du/dt.
"""
-(tile::AbstractTile{true})(du, u, p, t, dt=1.0) = computefluxes!(tile,du,u,p,t,dt)
+(tile::AbstractTile{true})(du, u, p, t, dt=1.0) = computeprognostic!(tile,du,u,p,t,dt)
(tile::AbstractTile{false})(u, p, t, dt=1.0) = computefluxes(tile,u,p,t,dt)
"""
@@ -29,11 +29,11 @@ Returns true if `tile` uses in-place mode, false if out-of-place.
isinplace(::AbstractTile{iip}) where {iip} = iip
"""
- computefluxes!(::T,du,u,p,t) where {T<:AbstractTile}
+ computeprognostic!(::T,du,u,p,t) where {T<:AbstractTile}
In-place update function for tile `T`. Computes du/dt and stores the result in `du`.
"""
-computefluxes!(::T, du, u, p, t, dt=1.0) where {T<:AbstractTile} = error("no implementation of in-place computefluxes! for $T")
+computeprognostic!(::T, du, u, p, t, dt=1.0) where {T<:AbstractTile} = error("no implementation of in-place computeprognostic! for $T")
"""
computefluxes(::T,u,p,t) where {T<:AbstractTile}
diff --git a/src/Utils/Utils.jl b/src/Utils/Utils.jl
index 185c5d7dd86310ce11687cf975c066ed56a0b4cb..4c20a1e76dfd9571f6b8259d85aa4829003f8362 100644
--- a/src/Utils/Utils.jl
+++ b/src/Utils/Utils.jl
@@ -172,22 +172,36 @@ function pstrip(obj; keep_units=false)
end
# TODO: this should be in ModelParameters.jl, not here.
-function Unitful.uconvert(u::Unitful.Units, p::Param)
+function Unitful.uconvert(u::Unitful.Units, p::paraType) where {paraType<:AbstractParam}
nt = parent(p)
@set! nt.val = ustrip(u, stripparams(p))
@set! nt.units = u
- return Param(nt)
+ return paraType(nt)
end
"""
ModelParameters.stripunits(obj)
-Additional override for `stripunits` which reconstructs `obj` with all fields that have unitful quantity
+Additional dispatch for `stripunits` which reconstructs `obj` with all fields that have unitful quantity
types converted to base SI units and then stripped to be unit free.
"""
function ModelParameters.stripunits(obj)
values = Flatten.flatten(obj, Flatten.flattenable, Unitful.AbstractQuantity, Flatten.IGNORE)
return Flatten.reconstruct(obj, map(ustrip ∘ normalize_units, values), Unitful.AbstractQuantity, Flatten.IGNORE)
end
+
+"""
+ ModelParameters.stripparams(::Type{paraType}, obj) where {paraType<:AbstractParam}
+
+Additional dispatch for `stripparams` which allows specification of a specific type.
+"""
+function ModelParameters.stripparams(::Type{paraType}, obj) where {paraType<:AbstractParam}
+ IgnoreTypes = ignored_types(obj)
+ selected_params = Flatten.flatten(obj, Flatten.flattenable, paraType, IgnoreTypes)
+ return Flatten._reconstruct(obj, map(stripparams, selected_params), Flatten.flattenable, paraType, IgnoreTypes, 1)[1]
+end
+
+ignored_types(obj) = ModelParameters.IGNORE
+
# pretty print Param types
Base.show(io::IO, mime::MIME"text/plain", p::Param) = print(io, "Param($(p.val))")
diff --git a/src/coupling.jl b/src/coupling.jl
index 1567f3b0cc0cd7aebb683717b5e49b8980303f1c..785975a68005d150dc9593457ada3a08f6d569f7 100644
--- a/src/coupling.jl
+++ b/src/coupling.jl
@@ -17,7 +17,7 @@ CryoGrid.initialcondition!(layer::Layer, ps::CoupledProcesses, state) = _invoke_
CryoGrid.computediagnostic!(layer::Layer, ps::CoupledProcesses, state) = _invoke_sequential(computediagnostic!, layer, ps, state)
-CryoGrid.computefluxes!(layer::Layer, ps::CoupledProcesses, state) = _invoke_sequential(computefluxes!, layer, ps, state)
+CryoGrid.computeprognostic!(layer::Layer, ps::CoupledProcesses, state) = _invoke_sequential(computeprognostic!, layer, ps, state)
CryoGrid.resetfluxes!(layer::Layer, ps::CoupledProcesses, state) = _invoke_sequential(resetfluxes!, layer, ps, state)
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 a78dce27ee8e3947d0d6f4d8151265bb429caaed..0900440ffe51145db29e8037459503b5f6485288 100644
--- a/src/methods.jl
+++ b/src/methods.jl
@@ -63,15 +63,15 @@ interact!(layer1::Layer, layer2::Layer, state1, state2) = interact!(layer1, proc
interact!(::Layer, ::Process, ::Layer, ::Process, state1, state2) = nothing
"""
- computefluxes!(l::Layer, p::Process, state)
+ computeprognostic!(l::Layer, p::Process, state)
-Calculates all internal fluxes for a given layer. Note that an instance of `computefluxes!` must be provided
+Calculates all internal fluxes for a given layer. Note that an instance of `computeprognostic!` must be provided
for all non-boundary (subsurface) processes/layers.
"""
-computefluxes!(layer::Layer, state) = computefluxes!(layer, processes(layer), state)
-computefluxes!(layer::Layer, proc::Process, state) = error("computefluxes! defined for $(typeof(layer)) with $(typeof(proc))")
-computefluxes!(::Top, ::BoundaryProcess, state) = nothing
-computefluxes!(::Bottom, ::BoundaryProcess, state) = nothing
+computeprognostic!(layer::Layer, state) = computeprognostic!(layer, processes(layer), state)
+computeprognostic!(layer::Layer, proc::Process, state) = error("computeprognostic! defined for $(typeof(layer)) with $(typeof(proc))")
+computeprognostic!(::Top, ::BoundaryProcess, state) = nothing
+computeprognostic!(::Bottom, ::BoundaryProcess, state) = nothing
"""
caninteract(layer1::Layer, layer2::Layer, state1, state2)
@@ -111,7 +111,7 @@ resetfluxes!(layer::Layer, proc::Process, state) = nothing
isactive(::Layer, state)
Returns a boolean whether or not this layer is currently active in the stratigraphy and should interact with other layers.
-Note that `computediagnostic!` and `computefluxes!` are always invoked regardless of the current state of `isactive`.
+Note that `computediagnostic!` and `computeprognostic!` are always invoked regardless of the current state of `isactive`.
The default implementation of `isactive` always returns `true`.
"""
isactive(::Layer, state) = true
@@ -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)
@@ -259,3 +233,18 @@ thickness(::Layer, state, i) = Δ(state.grid)[i]
thickness(l::Layer, state, ::typeof(first)) = thickness(l, state, 1)
thickness(l::Layer, state, ::typeof(last)) = thickness(l, state, lastindex(state.grid)-1)
thickness(::Union{Top,Bottom}, state) = Inf
+
+
+"""
+ param([::Type{paraType}], defval; kwargs...)
+
+Creates a new parameter type from the given default value and keyword properties.
+"""
+param(::Type{paraType}, defval; kwargs...) where {paraType<:AbstractParam} = paraType(defval; kwargs...)
+function param(defval; kwargs...)
+ if AUTOPARA
+ return param(Param, defval; kwargs...)
+ else
+ return param(FixedParam, defval; kwargs...)
+ end
+end
diff --git a/src/models.jl b/src/models.jl
new file mode 100644
index 0000000000000000000000000000000000000000..ce5ebbcef783ebb1a05f06f9a96c275ec1c99135
--- /dev/null
+++ b/src/models.jl
@@ -0,0 +1,44 @@
+"""
+ 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,
+ 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), inputs, inits...; tile_kwargs...)
+end
+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(
+ 0.0u"m" => SimpleSoil(por=0.80,sat=1.0,org=0.75,freezecurve=PainterKarra(swrc=VanGenuchten("organic"))),
+ 0.1u"m" => SimpleSoil(por=0.80,sat=1.0,org=0.25,freezecurve=PainterKarra(swrc=VanGenuchten("organic"))),
+ 0.4u"m" => SimpleSoil(por=0.55,sat=1.0,org=0.25,freezecurve=PainterKarra(swrc=VanGenuchten("sandy loam"))),
+ 3.0u"m" => SimpleSoil(por=0.50,sat=1.0,org=0.0,freezecurve=PainterKarra(swrc=VanGenuchten("sandy loam"))),
+ 10.0u"m" => SimpleSoil(por=0.30,sat=1.0,org=0.0,freezecurve=PainterKarra(swrc=VanGenuchten("sandy loam"))),
+ 30.0u"m" => SimpleSoil(por=0.05,sat=1.0,org=0.0,freezecurve=PainterKarra(swrc=VanGenuchten("sand"))),
+ ),
+ tempprofile = TemperatureProfile(
+ 0.0u"m" => -1.0u"°C",
+ 0.2u"m" => -1.0u"°C",
+ 2.0u"m" => -1.0u"°C",
+ 5.0u"m" => -3.0u"°C",
+ 9.5u"m" => -6.0u"°C",
+ 25.0u"m" => -9.0u"°C",
+ 100.0u"m" => -9.0u"°C",
+ 1000.0u"m" => 10.2u"°C",
+ )
+)
diff --git a/src/problem.jl b/src/problem.jl
index daf876a3f8c7e631ff3a3b826d2fe8ee76f8bff3..6c05881f386d286a86159e3d4c921dab5e90d1f9 100644
--- a/src/problem.jl
+++ b/src/problem.jl
@@ -47,8 +47,7 @@ Constructor for `CryoGridProblem` that automatically generates all necessary cal
function CryoGridProblem(
tile::Tile,
u0::ComponentVector,
- tspan::NTuple{2,Float64},
- p=nothing;
+ tspan::NTuple{2,Float64};
diagnostic_stepsize=3600.0,
saveat=3600.0,
savevars=(),
@@ -66,16 +65,10 @@ function CryoGridProblem(
)
getsavestate(tile::Tile, u, du) = deepcopy(Tiles.getvars(tile.state, Tiles.withaxes(u, tile), Tiles.withaxes(du, tile), savevars...))
savefunc(u, t, integrator) = getsavestate(Tile(integrator), Tiles.withaxes(u, Tile(integrator)), get_du(integrator))
- tile, p = if isnothing(p) && isempty(ModelParameters.params(tile))
- # case 1: no parameters provided
- tile, nothing
- else
- # case 2: parameters are provided; use Model interface to reconstruct Tile with new parameter values
- model_tile = Model(tile)
- p = isnothing(p) ? collect(model_tile[:val]) : p
- model_tile[:val] = p
- parent(model_tile), p
- end
+ # strip all "fixed" parameters
+ tile = stripparams(FixedParam, tile)
+ # retrieve variable parameters
+ p = length(ModelParameters.params(tile)) > 0 ? parameters(tile) : nothing
du0 = zero(u0)
# remove units
tile = stripunits(tile)
@@ -142,7 +135,7 @@ function CryoGrid.odefunction(::DefaultJac, setup::typeof(tile), u0, p, tspan)
# make sure to return an instance of ODEFunction
end
...
-prob = CryoGridProblem(tile, tspan, p)
+prob = CryoGridProblem(tile, u0, tspan)
```
`JacobianStyle` can also be extended to create custom traits which can then be applied to compatible `Tile`s.
diff --git a/test/CryoGridSBIExt/runtests.jl b/test/CryoGridSBIExt/runtests.jl
index 5ddb3e0bbbc20c32b76f87b190f354a6a01131ee..aa2f09df67ee49c4c48760531a26ce590232975a 100644
--- a/test/CryoGridSBIExt/runtests.jl
+++ b/test/CryoGridSBIExt/runtests.jl
@@ -15,14 +15,15 @@ function cryogrid_test_problem(tspan)
);
heatop = Heat.Diffusion1D(:H)
initT = initializer(:T, tempprofile)
- tile = CryoGrid.Presets.SoilHeatTile(
+ tile = CryoGrid.SoilHeatTile(
heatop,
## 10 W/m^2 in and out, i.e. net zero flux
PeriodicBC(HeatBalance, CryoGrid.Dirichlet, 365.0u"d", 15.0u"°C", 0.0, -5.0u"°C"),
ConstantBC(HeatBalance, CryoGrid.Neumann, -0.1u"W/m^2"),
soilprofile,
+ inputs(),
initT;
- grid=CryoGrid.Presets.DefaultGrid_10cm,
+ grid=CryoGrid.DefaultGrid_10cm,
);
u0, _ = initialcondition!(tile, tspan);
prob = CryoGridProblem(tile, u0, tspan)
diff --git a/test/IO/forcing_tests.jl b/test/IO/forcing_tests.jl
index 47be0699ea8110d9f182517f1757e0e4c3732c45..60f2dfdaf90ac19b85edfbbd1cd0abd86c734b22 100644
--- a/test/IO/forcing_tests.jl
+++ b/test/IO/forcing_tests.jl
@@ -2,23 +2,24 @@ using CryoGrid
using Dates
using Test, BenchmarkTools
-@testset "InterpolatedForcing" begin
+@testset "Interpolated1D" begin
ts = DateTime(2020,1,1):Day(1):DateTime(2020,3,1)
ys = randn(length(ts))
- forcing = CryoGrid.InterpolatedForcing(ts, ys, :test)
- @test all(forcing.interpolant.knots[1] .≈ convert_t.(ts))
- @test forcing(Dates.datetime2epochms(ts[1])/1000.0) ≈ ys[1]
- @test forcing(Dates.datetime2epochms(ts[end])/1000.0) ≈ ys[end]
- @test_throws BoundsError forcing(-1.0)
- @test_throws BoundsError forcing(Dates.datetime2epochms(ts[end])/1000.0+1.0)
+ arr = DimStack((test=DimArray(ys, (Ti(ts),)),))
+ interp = InputOutput.Interpolated1D(arr)
+ f = interp.test
+ @test f(Dates.datetime2epochms(ts[1])/1000.0) ≈ ys[1]
+ @test f(Dates.datetime2epochms(ts[end])/1000.0) ≈ ys[end]
+ @test_throws BoundsError f(-1.0)
+ @test_throws BoundsError f(Dates.datetime2epochms(ts[end])/1000.0+1.0)
t1,y1 = ts[1], ys[1]
t2,y2 = ts[2], ys[2]
- @test forcing((Dates.datetime2epochms(t1) + Dates.datetime2epochms(t2))/2000.0) ≈ (y1+y2)/2
+ @test f((Dates.datetime2epochms(t1) + Dates.datetime2epochms(t2))/2000.0) ≈ (y1+y2)/2
t = Dates.datetime2epochms(t1)/1000.0
- benchres = @benchmark $forcing($t)
+ benchres = @benchmark $f($t)
@test benchres.allocs == 0
out = zeros(100)
queries = t .+ (1:100);
- benchres = @benchmark $out .= $forcing.($queries)
+ benchres = @benchmark $out .= $f.($queries)
@test benchres.allocs == 0
end
\ No newline at end of file
diff --git a/test/Physics/Heat/heat_conduction_tests.jl b/test/Physics/Heat/heat_conduction_tests.jl
index 9a2df9563ad30471604a55db7c5e7eea160ee980..78ed483181d7a2d09623112a0fcda6e0722c35df 100644
--- a/test/Physics/Heat/heat_conduction_tests.jl
+++ b/test/Physics/Heat/heat_conduction_tests.jl
@@ -79,19 +79,22 @@ using Test
end
@testset "Boundary conditions" begin
@testset "n-factors" begin
- ts = DateTime(2010,1,1):Hour(1):DateTime(2010,1,1,4)
- forcing = InterpolatedForcing(ts, [1.0,0.5,-0.5,-1.0,0.1]u"°C", :Tair)
- tgrad = TemperatureBC(forcing, NFactor(nf=0.5, nt=1.0))
+ T_pos = TemperatureBC(10.0, NFactor(nf=0.5, nt=0.75))
+ T_neg = TemperatureBC(-10.0, NFactor(nf=0.5, nt=0.75))
heat = HeatBalance()
sub = TestGroundLayer(heat)
zerobc = ConstantBC(HeatBalance, CryoGrid.Dirichlet, 0.0u"°C")
function f1(t)
state = (T_ub=[Inf], nfactor=[Inf], t=t)
- computediagnostic!(Top(zerobc), tgrad, state)
- return boundaryvalue(tgrad, state)
+ computediagnostic!(Top(T_pos), T_pos, state)
+ T_ub = boundaryvalue(T_pos, state)
+ computediagnostic!(Top(T_neg), T_neg, state)
+ T_lb = boundaryvalue(T_neg, state)
+ return T_ub, T_lb
end
- Tres = f1.(Dates.datetime2epochms.(ts)./1000.0)
- @test all(Tres .≈ [1.0,0.5,-0.25,-0.5,0.1])
+ T1, T2 = f1(0.0)
+ @test T1 == 7.5
+ @test T2 == -5.0
end
end
@testset "Fourier solution" begin
@@ -117,7 +120,7 @@ end
state = (T=T,dH=dH,jH=jH,k=k,grid=x,grids=(T=xc,k=x),t=t)
interact!(Top(bc), bc, sub, heat, state, state)
interact!(sub, heat, Bottom(bc), bc, state, state)
- computefluxes!(sub, heat, state)
+ computeprognostic!(sub, heat, state)
# strip units from dH before returning it to the solver;
# note that we do not need to divide by diffusivity since we assume it to be unity
return ustrip.(dH)
diff --git a/test/Physics/Hydrology/water_balance_tests.jl b/test/Physics/Hydrology/water_balance_tests.jl
index 2c10d4f39ea397f6347d2eb8baee7a10a2eec591..b4f0873064469310c2253384728f7208bd17bcd0 100644
--- a/test/Physics/Hydrology/water_balance_tests.jl
+++ b/test/Physics/Hydrology/water_balance_tests.jl
@@ -6,9 +6,9 @@ using CryoGrid.Utils
using Test
@testset "Bucket scheme" begin
- testgrid = CryoGrid.Presets.DefaultGrid_2cm
+ testgrid = CryoGrid.DefaultGrid_2cm
@testset "variables" begin
- sub = TestGroundLayer(WaterBalance(BucketScheme()))
+ sub = TestGroundLayer(WaterBalance(BucketScheme())) |> stripparams
vars = CryoGrid.variables(sub)
prognostic_vars = filter(CryoGrid.isprognostic, vars)
diagnostic_vars = filter(CryoGrid.isdiagnostic, vars)
@@ -23,7 +23,7 @@ using Test
@test isa(proc, WaterBalance)
end
@testset "computediagnostic!" begin
- sub = TestGroundLayer(WaterBalance(BucketScheme()))
+ sub = TestGroundLayer(WaterBalance(BucketScheme())) |> stripparams
state = Diagnostics.build_dummy_state(testgrid, sub, with_units=true)
# initialize fully saturated
state.sat .= 1.0
@@ -32,14 +32,14 @@ using Test
@test allfinite(state.kw)
@test allfinite(state.θw)
end
- @testset "computefluxes!" begin
- sub = TestGroundLayer(WaterBalance(BucketScheme()))
+ @testset "computeprognostic!" begin
+ sub = TestGroundLayer(WaterBalance(BucketScheme())) |> stripparams
state = Diagnostics.build_dummy_state(testgrid, sub, with_units=true)
# initialize fully saturated
state.sat .= 1.0
procs = CryoGrid.processes(sub)
CryoGrid.computediagnostic!(sub, procs, state)
- CryoGrid.computefluxes!(sub, procs, state)
+ CryoGrid.computeprognostic!(sub, procs, state)
@test allfinite(state.kw)
@test allfinite(state.θw)
@test all(iszero.(state.jw))
@@ -47,14 +47,14 @@ using Test
state.sat .= 0.5
procs = CryoGrid.processes(sub)
CryoGrid.computediagnostic!(sub, procs, state)
- CryoGrid.computefluxes!(sub, procs, state)
+ CryoGrid.computeprognostic!(sub, procs, state)
@test allfinite(state.kw)
@test allfinite(state.θw)
@test all(state.jw[2:end-1] .> zero(eltype(state.jw)))
end
@testset "interact!" begin
- sub1 = TestGroundLayer(WaterBalance(BucketScheme()))
- sub2 = TestGroundLayer(WaterBalance(BucketScheme()))
+ sub1 = TestGroundLayer(WaterBalance(BucketScheme())) |> stripparams
+ sub2 = TestGroundLayer(WaterBalance(BucketScheme())) |> stripparams
state1 = Diagnostics.build_dummy_state(testgrid[0.0u"m"..10.0u"m"], sub1, with_units=true)
state2 = Diagnostics.build_dummy_state(testgrid[10.0u"m"..1000.0u"m"], sub2, with_units=true)
# case 1: both saturated
diff --git a/test/Physics/Salt/runtests.jl b/test/Physics/Salt/runtests.jl
index 8363c46cef3cb6044127ca31d8ce226b87910287..70f5154a42c50f92e3e5901d980c07dd764c6dd7 100644
--- a/test/Physics/Salt/runtests.jl
+++ b/test/Physics/Salt/runtests.jl
@@ -5,7 +5,7 @@ using CryoGrid.Utils
using Test
@testset "Salt" begin
- testgrid = CryoGrid.Presets.DefaultGrid_2cm
+ testgrid = CryoGrid.DefaultGrid_2cm
@testset "variables" begin
soil = SalineGround()
vars = CryoGrid.variables(soil)
@@ -63,7 +63,7 @@ using Test
@test state1.jc[end] > zero(eltype(state1.jc))
@test state2.jc[1] == state1.jc[end]
end
- @testset "computefluxes!" begin
+ @testset "computeprognostic!" begin
soil = pstrip(SalineGround())
state = Diagnostics.build_dummy_state(testgrid, soil, with_units=false)
# initialize variables
@@ -72,7 +72,7 @@ using Test
state.por .= porosity(soil)
procs = CryoGrid.processes(soil)
CryoGrid.computediagnostic!(soil, procs, state)
- CryoGrid.computefluxes!(soil, procs, state)
+ CryoGrid.computeprognostic!(soil, procs, state)
# check that all divergences are nonzero for both temperature and salt
@test all(.!iszero.(state.dc))
@test all(.!iszero.(state.dT))
diff --git a/test/Physics/Sources/runtests.jl b/test/Physics/Sources/runtests.jl
index 1338b5c218afc3c60cb6bc71ca1922d3267f2b50..0b72ea448db5b4bf6ced1f54469ebc8fec5284ff 100644
--- a/test/Physics/Sources/runtests.jl
+++ b/test/Physics/Sources/runtests.jl
@@ -6,10 +6,10 @@ using Test
heatsource = Source(HeatBalance, Sources.Constant(1.0u"W/m^3"))
layer = TestGroundLayer(heatsource)
state = (dH=zeros(100)u"W/m^3",)
- computefluxes!(layer, heatsource, state)
+ computeprognostic!(layer, heatsource, state)
@test all(state.dH .≈ 1.0u"W/m^3")
state = (dH=ones(100)u"W/m^3",)
- computefluxes!(layer, heatsource, state)
+ computeprognostic!(layer, heatsource, state)
@test all(state.dH .≈ 2.0u"W/m^3")
end
@testset "Periodic" begin
@@ -17,10 +17,10 @@ using Test
heatsource = Source(HeatBalance, Sources.Periodic(; level, amp, freq, shift))
layer = TestGroundLayer(heatsource)
state = (dH=zeros(100)u"W/m^3", t=1.0u"s")
- computefluxes!(layer, heatsource, state)
+ computeprognostic!(layer, heatsource, state)
@test all(state.dH .≈ level + amp*sin(2π*freq*1.0u"s" - shift))
state = (dH=ones(100)u"W/m^3", t=1.5u"s")
- computefluxes!(layer, heatsource, state)
+ computeprognostic!(layer, heatsource, state)
@test all(state.dH .≈ 1.0u"W/m^3" + level + amp*sin(2π*freq*1.5u"s" - shift))
end
end
diff --git a/test/Physics/Surface/seb_tests.jl b/test/Physics/Surface/seb_tests.jl
index e67ac24eda3ef5e1e428a53eac7d9c38a4bbc068..d4d96beebdb5c9261795ba5a9a4e83848ada29b8 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,
+ seb = SurfaceEnergyBalance(;
+ forcings...,
solscheme = Surface.Iterative(),
)
# TODO: would be good to have SEB work correctly with units to verify correctness;
@@ -35,14 +29,8 @@ Surface.surfaceproperties(seb::SurfaceEnergyBalance, sub::TestGroundLayer) = seb
@test seb_output.Qg > zero(seb_output.Qg)
end
@testset "Analytical" begin
- seb = SurfaceEnergyBalance(
- const_Tair,
- const_pr,
- const_qh,
- const_wind,
- const_Lin,
- const_Sin,
- z,
+ seb = SurfaceEnergyBalance(;
+ forcings...,
solscheme = Surface.Analytical(),
)
seb = pstrip(seb)
@@ -58,16 +46,10 @@ Surface.surfaceproperties(seb::SurfaceEnergyBalance, sub::TestGroundLayer) = seb
@test seb_output.Qg > zero(seb_output.Qg)
end
@testset "Numerical" begin
- seb = SurfaceEnergyBalance(
- const_Tair,
- const_pr,
- const_qh,
- const_wind,
- const_Lin,
- const_Sin,
- z,
- solscheme=Surface.Numerical()
- )
+ seb = SurfaceEnergyBalance(;
+ forcings...,
+ solscheme=Surface.Numerical(),
+ )
seb = pstrip(seb)
grid = Grid([0.0,0.1]u"m")
toplayer = Top(seb)
diff --git a/test/Physics/Surface/swb_tests.jl b/test/Physics/Surface/swb_tests.jl
index c6aabd6b7def090269656a712af194f915830494..96cf4dbd64e2ba98b40dd8b198b6b0001960720b 100644
--- a/test/Physics/Surface/swb_tests.jl
+++ b/test/Physics/Surface/swb_tests.jl
@@ -1,13 +1,8 @@
@testset "Surface Water Balance" begin
- testgrid = CryoGrid.Presets.DefaultGrid_2cm
+ testgrid = CryoGrid.DefaultGrid_2cm
@testset "interact!" begin
- sub = TestGroundLayer(WaterBalance(BucketScheme()))
- top = Top(
- SurfaceWaterBalance(
- ConstantForcing(1e-6u"m/s", :rainfall),
- ConstantForcing(0.0u"m/s", :snowfall),
- )
- )
+ sub = stripparams(TestGroundLayer(WaterBalance(BucketScheme())))
+ top = Top(SurfaceWaterBalance(1e-6u"m/s", 0.0u"m/s"))
stop = Diagnostics.build_dummy_state(testgrid, top, with_units=true)
ssub = Diagnostics.build_dummy_state(testgrid, sub, with_units=true)
# initialize fully saturated
@@ -17,7 +12,7 @@
CryoGrid.computediagnostic!(top, swb, stop)
CryoGrid.computediagnostic!(sub, water, ssub)
CryoGrid.interact!(top, sub, stop, ssub)
- CryoGrid.computefluxes!(top, stop)
+ CryoGrid.computeprognostic!(top, stop)
@test iszero(ssub.jw[1])
@test stop.drunoff[1] > zero(eltype(stop.drunoff))
ssub.sat .= 0.5
@@ -26,7 +21,7 @@
CryoGrid.computediagnostic!(top, swb, stop)
CryoGrid.computediagnostic!(sub, water, ssub)
CryoGrid.interact!(top, sub, stop, ssub)
- CryoGrid.computefluxes!(top, stop)
+ CryoGrid.computeprognostic!(top, stop)
@test ssub.jw[1] > zero(eltype(ssub.jw))
@test iszero(stop.drunoff[1])
end
diff --git a/test/Solvers/cglite_implicit_tests.jl b/test/Solvers/cglite_implicit_tests.jl
index bc91216c520c2d2ce0b19c672dab7edfba6809e7..5c437e86586f25f1c5f4b5030d71c0debac9aa7c 100644
--- a/test/Solvers/cglite_implicit_tests.jl
+++ b/test/Solvers/cglite_implicit_tests.jl
@@ -27,7 +27,7 @@ end
α = upreferred(strat.soil.para.heat.kh_m) / upreferred(strat.soil.para.heat.ch_m)
T_analytic = Heat.heat_conduction_linear_periodic_ub(T₀, A, P, ustrip(α))
initT = initializer(:T, (layer, state) -> state.T .= T_analytic.(cells(state.grid), 0.0))
- modelgrid = CryoGrid.Presets.DefaultGrid_2cm
+ modelgrid = CryoGrid.DefaultGrid_2cm
# modelgrid = Grid(z_top:0.02u"m":z_bot)
# modelgrid = Grid(vcat(0.0:0.02:1.0, 1.05:0.05:5.0, 5.1:0.1:10.0)*u"m")
tile = Tile(strat, modelgrid, initT)
@@ -58,7 +58,7 @@ end
z_bot => Bottom(ConstantFlux(HeatBalance, 0.0))
);
initT = initializer(:T, -1.0)
- modelgrid = CryoGrid.Presets.DefaultGrid_2cm
+ modelgrid = CryoGrid.DefaultGrid_2cm
tile = Tile(strat, modelgrid, initT)
# define time span, 5 years
tspan = (0.0, 5*365*24*3600.0)
@@ -68,7 +68,6 @@ end
sol = solve(prob, LiteImplicitEuler(), dt=24*3600)
out = CryoGridOutput(sol)
- kh_w, kh_i, kh_a, kh_m, kh_o = Heat.thermalconductivities(strat.soil)
θm = mineral(strat.soil)
θo = organic(strat.soil)
θ_s = (θw=0.0, θi=porosity(strat.soil), θa=0.0, θm=θm, θo=θo)
diff --git a/test/benchmarks.jl b/test/benchmarks.jl
index 57680443b97fa31a008ceca7d5a338ed44b79fcd..19e67e87c353b50b959b296f0fa067d7c9410b3a 100644
--- a/test/benchmarks.jl
+++ b/test/benchmarks.jl
@@ -15,7 +15,7 @@ function run_solver_benchmark(solvers, grids, dts, tspan, upperbc;
# high accuracy reference solution, 1 minute time steps
println("Computing reference solution for grid spacing $min_dx with forward Euler @ $dt_ref s time steps ...")
# basic 1-layer heat conduction model (defaults to free water freezing scheme)
- model = Presets.SoilHeat(:H, upperbc, profile; grid=grid, freezecurve=freezecurve)
+ model = SoilHeat(:H, upperbc, profile; grid=grid, freezecurve=freezecurve)
p = copy(model.pproto)
p .= params
prob = CryoGridProblem(model,tspan,p)
@@ -55,7 +55,7 @@ forcings = loadforcings(filename);
# use air temperature as upper boundary forcing
tair = TemperatureBC(forcings.Tair)
solvers = [Euler, DP5, ROCK2, ROCK4, Trapezoid, ROS3P]
-grids = [Presets.DefaultGrid_2cm, Presets.DefaultGrid_5cm, Presets.DefaultGrid_10cm, Presets.DefaultGrid_20cm]
+grids = [DefaultGrid_2cm, DefaultGrid_5cm, DefaultGrid_10cm, DefaultGrid_20cm]
dts = [2*60.0, 10*60.0, 30*60.0, 3600.0]
tspan = (DateTime(2010,10,30),DateTime(2011,10,30))
results_freeW_tair = run_solver_benchmark(solvers, grids, dts, tspan, tair)
@@ -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)
diff --git a/test/test_problems.jl b/test/test_problems.jl
index d12036f201ce6d9fe4bbae0f764efd2c0e35d547..93c81569e42a8b93acceaa2e3f8601472cf7ad07 100644
--- a/test/test_problems.jl
+++ b/test/test_problems.jl
@@ -13,9 +13,9 @@ function test_heat_conduction_freeW_periodic_bc()
initT = initializer(:T, tempprofile)
# Define the simulation time span.
tspan = (0.0,2*24*3600.0)
- discretization = CryoGrid.Presets.DefaultGrid_10cm
+ discretization = CryoGrid.DefaultGrid_10cm
Qbot = -0.01u"W/m^2"
- tile = CryoGrid.Presets.SoilHeatTile(
+ tile = CryoGrid.SoilHeatTile(
heatop,
PeriodicBC(HeatBalance, CryoGrid.Neumann, 24*3600.0, 20.0, 0.0, -5.0),
ConstantBC(HeatBalance, CryoGrid.Neumann, Qbot),
@@ -32,7 +32,7 @@ function test_water_flow_bucket_scheme(topflux=0.0u"m/s")
initsat = initializer(:sat, 0.5)
# Define the simulation time span.
tspan = (0.0, 24*3600.0)
- discretization = CryoGrid.Presets.DefaultGrid_10cm
+ discretization = CryoGrid.DefaultGrid_10cm
strat = Stratigraphy(
0.0u"m" => Top(ConstantBC(WaterBalance, CryoGrid.Neumann, topflux),),
0.0u"m" => Ground(SimpleSoil(por=0.5, sat=0.5), heat=nothing, water=WaterBalance(BucketScheme())),