From 9f0d505f1a982b19618f19328549c854590f0d14 Mon Sep 17 00:00:00 2001
From: Brian Groenke <brian.groenke@awi.de>
Date: Tue, 26 Nov 2024 16:17:03 +0100
Subject: [PATCH] Fix more miscellaneous issues
---
examples/heat_freeW_lake_lite_implicit.jl | 7 +++-
examples/heat_freeW_lite_implicit.jl | 40 +++++++++--------------
examples/heat_sfcc_richardseq_samoylov.jl | 2 ++
src/IO/InputOutput.jl | 2 +-
src/IO/forcings/providers.jl | 5 +++
src/Physics/Heat/heat_water.jl | 12 ++++---
src/Tiles/stratigraphy.jl | 4 ++-
src/Tiles/tile.jl | 8 ++++-
8 files changed, 48 insertions(+), 32 deletions(-)
diff --git a/examples/heat_freeW_lake_lite_implicit.jl b/examples/heat_freeW_lake_lite_implicit.jl
index dc4d2da8..2cc4cbaa 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.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),
diff --git a/examples/heat_freeW_lite_implicit.jl b/examples/heat_freeW_lite_implicit.jl
index aa88e557..f818f2f8 100644
--- a/examples/heat_freeW_lite_implicit.jl
+++ b/examples/heat_freeW_lite_implicit.jl
@@ -7,11 +7,17 @@
using CryoGrid
using CryoGrid.LiteImplicit
-# Load forcings and build stratigraphy like before.
-forcings = loadforcings(CryoGrid.Forcings.Samoylov_ERA_MkL3_CCSM4_long_term);
-forcings = Base.rename(forcings, :Ptot => :precip)
-precip_values = forcings.precip.interpolant.coefs
-precip_values .*= 2
+# 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(
@@ -29,7 +35,7 @@ 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),
@@ -40,11 +46,11 @@ 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_sfcc_richardseq_samoylov.jl b/examples/heat_sfcc_richardseq_samoylov.jl
index 877d07f8..a2e57485 100644
--- a/examples/heat_sfcc_richardseq_samoylov.jl
+++ b/examples/heat_sfcc_richardseq_samoylov.jl
@@ -5,6 +5,8 @@
# pressure gradients.
using CryoGrid
+
+
forcings = loadforcings(CryoGrid.Forcings.Samoylov_ERA_obs_fitted_1979_2014_spinup_extended_2044);
tspan = (DateTime(2011,1,1),DateTime(2011,12,31))
T0 = forcings.Tair(tspan[1])
diff --git a/src/IO/InputOutput.jl b/src/IO/InputOutput.jl
index b055d861..aeb237e9 100644
--- a/src/IO/InputOutput.jl
+++ b/src/IO/InputOutput.jl
@@ -67,7 +67,7 @@ export Input, InputProvider, InputFunctionProvider
export inputs
include("input.jl")
-export ForcingFormat, ForcingFormatJSON, ForcingFormatNCD
+export ForcingFormat, ForcingFormatJSON, ForcingFormatNCD, Interpolated1D
export loadforcings
include("forcings/forcings.jl")
diff --git a/src/IO/forcings/providers.jl b/src/IO/forcings/providers.jl
index dfbaf81b..0d3cde62 100644
--- a/src/IO/forcings/providers.jl
+++ b/src/IO/forcings/providers.jl
@@ -41,6 +41,11 @@ function Base.show(io::IO, mimetype::MIME"text/plain", data::Interpolated1D)
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
diff --git a/src/Physics/Heat/heat_water.jl b/src/Physics/Heat/heat_water.jl
index 18a403e5..abdc2cfc 100644
--- a/src/Physics/Heat/heat_water.jl
+++ b/src/Physics/Heat/heat_water.jl
@@ -18,7 +18,6 @@ advectiveflux(jw, Tâ‚, Tâ‚‚, cw, L) = jw*(cw*Tâ‚*(jw > zero(jw)) + cw*Tâ‚‚*(jw
Adds advective energy fluxes for all internal grid cell faces.
"""
function water_energy_advection!(jH, jw, T, cw::Real, L::Real)
- water, heat = ps
@inbounds for i in 2:length(jw)-1
let jw = jw[i],
Tâ‚ = T[i-1],
@@ -86,18 +85,23 @@ function CryoGrid.interact!(
end
end
-# Flux calculation
function CryoGrid.computeprognostic!(sub::SubSurface, ps::Coupled(WaterBalance, HeatBalance), state)
water, heat = ps
- cw = heatcapacitywater(sub, state)
- L = heat.prop.L
CryoGrid.computeprognostic!(sub, water, state)
if heat.advection
+ cw = heatcapacitywater(sub, state)
+ L = heat.prop.L
water_energy_advection!(state.jH, state.jw, state.T, cw, L)
end
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)
if hasproperty(state, :ch_w)
# TODO: this is a temporary solution that works so long as the
diff --git a/src/Tiles/stratigraphy.jl b/src/Tiles/stratigraphy.jl
index e7fab364..2d2d65e3 100644
--- a/src/Tiles/stratigraphy.jl
+++ b/src/Tiles/stratigraphy.jl
@@ -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
diff --git a/src/Tiles/tile.jl b/src/Tiles/tile.jl
index 3b57f17d..96f2cd4f 100644
--- a/src/Tiles/tile.jl
+++ b/src/Tiles/tile.jl
@@ -34,7 +34,13 @@ 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(
--
GitLab