From 4a6f64f1f01dba424b98cd94cfc7c64412bef571 Mon Sep 17 00:00:00 2001 From: Brian Groenke Date: Fri, 22 Sep 2023 13:12:57 +0200 Subject: [PATCH 1/8] Rework layer name handling in stratigraphy Also adds support for autogenerated subsurface layer names (finally). --- examples/04_heat_freeW_snow_samoylov.jl | 10 +- examples/05_heat_freeW_bucketW_samoylov.jl | 8 +- ...06_heat_freeW_seb_snow_bucketW_samoylov.jl | 12 +- examples/07_heat_freeW_lite_implicit.jl | 10 +- examples/08_heat_sfcc_richardseq_samoylov.jl | 6 +- examples/09_heat_sfcc_salt_constantbc.jl | 2 +- src/Presets/Presets.jl | 2 +- src/Tiles/Tiles.jl | 2 +- src/Tiles/stratigraphy.jl | 137 +++++++++++------- src/Tiles/tile.jl | 10 +- src/Utils/types.jl | 1 + src/problem.jl | 12 +- test/Physics/Surface/swb_tests.jl | 7 +- test/Tiles/stratigraphy_tests.jl | 60 ++++---- test/Tiles/tile_tests.jl | 10 +- 15 files changed, 161 insertions(+), 128 deletions(-) diff --git a/examples/04_heat_freeW_snow_samoylov.jl b/examples/04_heat_freeW_snow_samoylov.jl index 3d46284f..8888fd5b 100644 --- a/examples/04_heat_freeW_snow_samoylov.jl +++ b/examples/04_heat_freeW_snow_samoylov.jl @@ -29,11 +29,11 @@ snowpack = Snowpack( strat = @Stratigraphy( z_top => Top(upperbc), z_top => :snowpack => snowpack, - z_sub[1] => :topsoil1 => Ground(soilprofile[1].value, heat=HeatBalance(), water=WaterBalance(BucketScheme())), - z_sub[2] => :topsoil2 => Ground(soilprofile[2].value, heat=HeatBalance(), water=WaterBalance(BucketScheme())), - z_sub[3] => :sediment1 => Ground(soilprofile[3].value, heat=HeatBalance(), water=WaterBalance(BucketScheme())), - z_sub[4] => :sediment2 => Ground(soilprofile[4].value, heat=HeatBalance(), water=WaterBalance(BucketScheme())), - z_sub[5] => :sediment3 => Ground(soilprofile[5].value, heat=HeatBalance(), water=WaterBalance(BucketScheme())), + z_sub[1] => Ground(soilprofile[1].value, heat=HeatBalance(), water=WaterBalance(BucketScheme())), + z_sub[2] => Ground(soilprofile[2].value, heat=HeatBalance(), water=WaterBalance(BucketScheme())), + z_sub[3] => Ground(soilprofile[3].value, heat=HeatBalance(), water=WaterBalance(BucketScheme())), + z_sub[4] => Ground(soilprofile[4].value, heat=HeatBalance(), water=WaterBalance(BucketScheme())), + z_sub[5] => Ground(soilprofile[5].value, heat=HeatBalance(), water=WaterBalance(BucketScheme())), z_bot => Bottom(GeothermalHeatFlux(0.053u"J/s/m^2")) ); modelgrid = CryoGrid.Presets.DefaultGrid_5cm diff --git a/examples/05_heat_freeW_bucketW_samoylov.jl b/examples/05_heat_freeW_bucketW_samoylov.jl index cc158d9c..78ec304d 100644 --- a/examples/05_heat_freeW_bucketW_samoylov.jl +++ b/examples/05_heat_freeW_bucketW_samoylov.jl @@ -9,15 +9,15 @@ forcings = loadforcings(CryoGrid.Presets.Forcings.Samoylov_ERA_obs_fitted_1979_2 _, tempprofile = CryoGrid.Presets.SamoylovDefault; initT = initializer(:T, tempprofile) initsat = initializer(:sat, (l,state) -> state.sat .= l.para.sat); -upperbc = WaterHeatBC(SurfaceWaterBalance(rainfall=forcings.rainfall), TemperatureGradient(forcings.Tair, NFactor(0.5,0.9))); +upperbc = WaterHeatBC(SurfaceWaterBalance(forcings), TemperatureGradient(forcings.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. strat = @Stratigraphy( 0.0u"m" => Top(upperbc), - 0.0u"m" => :topsoil => Ground(MineralOrganic(por=0.80,sat=0.5,org=0.75), heat=HeatBalance(), water=WaterBalance(BucketScheme())), - 0.2u"m" => :subsoil => Ground(MineralOrganic(por=0.40,sat=0.75,org=0.10), heat=HeatBalance(), water=WaterBalance(BucketScheme())), - 2.0u"m" => :substrat => Ground(MineralOrganic(por=0.10,sat=1.0,org=0.0), heat=HeatBalance(), water=WaterBalance(BucketScheme())), + 0.0u"m" => Ground(MineralOrganic(por=0.80,sat=0.5,org=0.75), heat=HeatBalance(), water=WaterBalance(BucketScheme())), + 0.2u"m" => :middle => Ground(MineralOrganic(por=0.40,sat=0.75,org=0.10), heat=HeatBalance(), water=WaterBalance(BucketScheme())), + 2.0u"m" => Ground(MineralOrganic(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 diff --git a/examples/06_heat_freeW_seb_snow_bucketW_samoylov.jl b/examples/06_heat_freeW_seb_snow_bucketW_samoylov.jl index 035a5456..632dfb64 100644 --- a/examples/06_heat_freeW_seb_snow_bucketW_samoylov.jl +++ b/examples/06_heat_freeW_seb_snow_bucketW_samoylov.jl @@ -30,12 +30,12 @@ water = WaterBalance(BucketScheme(), DampedET()) ## build stratigraphy strat = @Stratigraphy( -z => Top(upperbc), - -z => :snowpack => Snowpack(heat=HeatBalance(), water=water), - soilprofile[1].depth => :soil1 => Ground(soilprofile[1].value; heat, water), - soilprofile[2].depth => :soil2 => Ground(soilprofile[2].value; heat, water), - soilprofile[3].depth => :soil3 => Ground(soilprofile[3].value; heat, water), - soilprofile[4].depth => :soil4 => Ground(soilprofile[4].value; heat, water), - soilprofile[5].depth => :soil5 => Ground(soilprofile[5].value; heat, water), + -z => Snowpack(heat=HeatBalance(), water=water), + soilprofile[1].depth => Ground(soilprofile[1].value; heat, water), + soilprofile[2].depth => Ground(soilprofile[2].value; heat, water), + soilprofile[3].depth => Ground(soilprofile[3].value; heat, water), + soilprofile[4].depth => Ground(soilprofile[4].value; heat, water), + soilprofile[5].depth => Ground(soilprofile[5].value; heat, water), 1000.0u"m" => Bottom(GeothermalHeatFlux(0.053u"J/s/m^2")), ); ## create Tile diff --git a/examples/07_heat_freeW_lite_implicit.jl b/examples/07_heat_freeW_lite_implicit.jl index 6e6e6c6a..d7ca206c 100644 --- a/examples/07_heat_freeW_lite_implicit.jl +++ b/examples/07_heat_freeW_lite_implicit.jl @@ -23,11 +23,11 @@ freezecurve = FreeWater() # freezecurve = PainterKarra(swrc=VanGenuchten(α=0.1, n=1.8)) strat = @Stratigraphy( z_top => Top(upperbc), - 0.0u"m" => :topsoil1 => Ground(MineralOrganic(por=0.80,sat=1.0,org=0.75), heat=HeatBalance(heatop; freezecurve)), - 0.1u"m" => :topsoil2 => Ground(MineralOrganic(por=0.80,sat=1.0,org=0.25), heat=HeatBalance(heatop; freezecurve)), - 0.4u"m" => :sediment1 => Ground(MineralOrganic(por=0.55,sat=1.0,org=0.25), heat=HeatBalance(heatop; freezecurve)), - 3.0u"m" => :sediment2 => Ground(MineralOrganic(por=0.50,sat=1.0,org=0.0), heat=HeatBalance(heatop; freezecurve)), - 10.0u"m" => :sediment3 => Ground(MineralOrganic(por=0.30,sat=1.0,org=0.0), heat=HeatBalance(heatop; freezecurve)), + 0.0u"m" => Ground(MineralOrganic(por=0.80,sat=1.0,org=0.75), heat=HeatBalance(heatop; freezecurve)), + 0.1u"m" => Ground(MineralOrganic(por=0.80,sat=1.0,org=0.25), heat=HeatBalance(heatop; freezecurve)), + 0.4u"m" => Ground(MineralOrganic(por=0.55,sat=1.0,org=0.25), heat=HeatBalance(heatop; freezecurve)), + 3.0u"m" => Ground(MineralOrganic(por=0.50,sat=1.0,org=0.0), heat=HeatBalance(heatop; freezecurve)), + 10.0u"m" => Ground(MineralOrganic(por=0.30,sat=1.0,org=0.0), heat=HeatBalance(heatop; freezecurve)), z_bot => Bottom(GeothermalHeatFlux(0.053u"W/m^2")) ); modelgrid = CryoGrid.Presets.DefaultGrid_2cm diff --git a/examples/08_heat_sfcc_richardseq_samoylov.jl b/examples/08_heat_sfcc_richardseq_samoylov.jl index d59ededd..4be164af 100644 --- a/examples/08_heat_sfcc_richardseq_samoylov.jl +++ b/examples/08_heat_sfcc_richardseq_samoylov.jl @@ -30,9 +30,9 @@ upperbc = WaterHeatBC(SurfaceWaterBalance(rainfall=forcings.rainfall), Temperatu # Note that the @Stratigraphy macro lets us list multiple subsurface layers without wrapping them in a tuple. strat = @Stratigraphy( -2.0u"m" => Top(upperbc), - 0.0u"m" => :topsoil => Ground(MineralOrganic(por=0.80,sat=0.7,org=0.75), heat=HeatBalance(heatop, freezecurve=sfcc), water=WaterBalance(RichardsEq(;swrc))), - 0.2u"m" => :subsoil => Ground(MineralOrganic(por=0.40,sat=0.8,org=0.10), heat=HeatBalance(heatop, freezecurve=sfcc), water=WaterBalance(RichardsEq(;swrc))), - 2.0u"m" => :substrat => Ground(MineralOrganic(por=0.10,sat=1.0,org=0.0), heat=HeatBalance(heatop, freezecurve=sfcc), water=WaterBalance(RichardsEq(;swrc))), + 0.0u"m" => Ground(MineralOrganic(por=0.80,sat=0.7,org=0.75), heat=HeatBalance(heatop, freezecurve=sfcc), water=WaterBalance(RichardsEq(;swrc))), + 0.2u"m" => Ground(MineralOrganic(por=0.40,sat=0.8,org=0.10), heat=HeatBalance(heatop, freezecurve=sfcc), water=WaterBalance(RichardsEq(;swrc))), + 2.0u"m" => Ground(MineralOrganic(por=0.10,sat=1.0,org=0.0), heat=HeatBalance(heatop, freezecurve=sfcc), water=WaterBalance(RichardsEq(;swrc))), 1000.0u"m" => Bottom(GeothermalHeatFlux(0.053u"W/m^2")) ); grid = CryoGrid.Presets.DefaultGrid_2cm diff --git a/examples/09_heat_sfcc_salt_constantbc.jl b/examples/09_heat_sfcc_salt_constantbc.jl index 1c0d98b9..eac2e23b 100644 --- a/examples/09_heat_sfcc_salt_constantbc.jl +++ b/examples/09_heat_sfcc_salt_constantbc.jl @@ -11,7 +11,7 @@ sfcc = DallAmicoSalt(swrc=VanGenuchten(α=4.06, n=2.03)) upperbc = SaltHeatBC(SaltGradient(benthicSalt=900.0, surfaceState=0), ConstantTemperature(0.0)) strat = @Stratigraphy( 0.0u"m" => Top(upperbc), - 0.0u"m" => :sediment => SalineGround(heat=HeatBalance(:T, freezecurve=sfcc), salt=SaltMassBalance()), + 0.0u"m" => SalineGround(heat=HeatBalance(:T, freezecurve=sfcc), salt=SaltMassBalance()), 1000.0u"m" => Bottom(GeothermalHeatFlux(0.053u"W/m^2")) ); modelgrid = CryoGrid.Presets.DefaultGrid_10cm diff --git a/src/Presets/Presets.jl b/src/Presets/Presets.jl index 201f34f1..32188187 100644 --- a/src/Presets/Presets.jl +++ b/src/Presets/Presets.jl @@ -28,7 +28,7 @@ but this can be changed via the `freezecurve` parameter. For example, to use the function SoilHeatTile(heatop, upperbc::BoundaryProcess, lowerbc::BoundaryProcess, soilprofile::Profile, init::VarInitializer; grid::Grid=DefaultGrid_5cm, freezecurve::F=FreeWater(), tile_kwargs...) where {F<:FreezeCurve} strat = Stratigraphy( grid[1] => Top(upperbc), - Tuple(knot.depth => Symbol(:soil,i) => Ground(knot.value, heat=HeatBalance(heatop, freezecurve=freezecurve)) for (i,knot) in enumerate(soilprofile)), + Tuple(knot.depth => Ground(knot.value, heat=HeatBalance(heatop, freezecurve=freezecurve)) for (i,knot) in enumerate(soilprofile)), grid[end] => Bottom(lowerbc) ) return Tile(strat, PresetGrid(grid), init; tile_kwargs...) diff --git a/src/Tiles/Tiles.jl b/src/Tiles/Tiles.jl index 19e5c11c..312457fc 100644 --- a/src/Tiles/Tiles.jl +++ b/src/Tiles/Tiles.jl @@ -32,7 +32,7 @@ export TileState, LayerState include("state.jl") export Stratigraphy, @Stratigraphy, NamedLayer -export layernames, layertypes, layers, boundaries, layername +export layernames, layertypes, layers, boundaries, layername, namedlayers export top, bottom, subsurface include("stratigraphy.jl") diff --git a/src/Tiles/stratigraphy.jl b/src/Tiles/stratigraphy.jl index 5a681b2c..d8842426 100644 --- a/src/Tiles/stratigraphy.jl +++ b/src/Tiles/stratigraphy.jl @@ -2,51 +2,51 @@ const RESERVED_LAYER_NAMES = (:top, :bottom, :strat, :init, :event) # NamedLayer type (alias of Named for Layer types) const NamedLayer{name,TLayer} = Named{name,TLayer} where {name,TLayer<:Layer} -layertype(layer::NamedLayer) = layertype(typeof(layer)) -layertype(::Type{<:NamedLayer{name,TLayer}}) where {name,TLayer} = TLayer -layername(layer::NamedLayer) = layername(typeof(layer)) -layername(::Type{<:NamedLayer{name}}) where {name} = name """ - Stratigraphy{N,TLayers,TBoundaries} + Stratigraphy{N,TLayers<:NamedTuple,TBoundaries} Defines a 1-dimensional stratigraphy by connecting a top and bottom layer to one or more subsurface layers. """ -struct Stratigraphy{N,TLayers,TBoundaries} +struct Stratigraphy{N,TLayers<:NamedTuple,TBoundaries} boundaries::TBoundaries layers::TLayers - Stratigraphy(boundaries::NTuple{N,Any}, layers::NTuple{N,NamedLayer}) where {N} = new{N,typeof(layers),typeof(boundaries)}(boundaries, layers) + Stratigraphy(boundaries::NTuple{N,Any}, layers::NamedTuple) where {N} = new{N,typeof(layers),typeof(boundaries)}(boundaries, layers) Stratigraphy( top::Pair{<:DistQuantity,<:Top}, - sub::Pair{<:DistQuantity,<:Pair{Symbol,<:SubSurface}}, + sub::Pair{<:DistQuantity}, bot::Pair{<:DistQuantity,<:Bottom} ) = Stratigraphy(top,(sub,),bot) Stratigraphy( top::Pair{<:DistQuantity,<:Top}, - sub::AbstractVector{<:Pair{<:DistQuantity,<:Pair{Symbol,<:SubSurface}}}, + sub::AbstractVector{<:Pair{<:DistQuantity}}, bot::Pair{<:DistQuantity,<:Bottom} ) = Stratigraphy(top, Tuple(sub), bot) function Stratigraphy( # use @nospecialize to (hopefully) reduce compilation overhead @nospecialize(top::Pair{<:DistQuantity,<:Top}), - @nospecialize(sub::Tuple{Vararg{Pair{<:DistQuantity,<:Pair{Symbol,<:SubSurface}}}}), + @nospecialize(sub::Tuple{Vararg{Pair{<:DistQuantity}}}), @nospecialize(bot::Pair{<:DistQuantity,<:Bottom}) ) + # check subsurface layers @assert length(sub) > 0 "At least one subsurface layer must be specified" length(sub) > 18 && @warn "Stratigraphies with more than 20 layers will result in very long compile times and are not recommended. Consider creating heterogeneous layer types." top = top[1] => Named(:top => top[2]) bot = bot[1] => Named(:bottom => bot[2]) - sub = map(x -> x[1] => Named(x[2]), sub) - names = map(nameof, map(last, sub)) - @assert length(unique(names)) == length(names) "All layer names in Stratigraphy must be unique" + sub = _withnames(sub) + sub_names = map(nameof, map(last, sub)) + @assert all(map(∉(RESERVED_LAYER_NAMES), sub_names)) "Subsurface layer names may not be one of: $RESERVED_LAYER_NAMES" boundaries = Tuple(map(first, (top, sub..., bot))) + # check boundary depths @assert issorted(boundaries) "Stratigraphy boundary locations must be in strictly increasing order." # wrap layers - layers = tuple(last(top), map(last, sub)..., last(bot)) + named_layers = tuple(last(top), map(last, sub)..., last(bot)) + layers = NamedTuple(named_layers) # construct type new{length(layers),typeof(layers),typeof(boundaries)}(boundaries, layers) end end + """ Convenience macro for defining stratigraphies with multiple subsurface layers. """ @@ -58,18 +58,20 @@ macro Stratigraphy(args...) :(Stratigraphy($(esc(args[1])), tuple($(esc.(args[2:end-1])...)), $(esc(args[end])))) end end + layers(strat::Stratigraphy) = getfield(strat, :layers) boundaries(strat::Stratigraphy) = getfield(strat, :boundaries) boundarypairs(strat::Stratigraphy) = boundarypairs(boundaries(strat)) boundarypairs(bounds::NTuple) = tuple(map(tuple, bounds[1:end-1], bounds[2:end])..., (bounds[end], bounds[end])) -layernames(strat::Stratigraphy) = map(layername, layers(strat)) +layernames(strat::Stratigraphy) = keys(layers(strat)) layertypes(::Type{<:Stratigraphy{N,TLayers}}) where {N,TLayers} = map(layertype, TLayers.parameters) +namedlayers(strat::Stratigraphy) = Tuple(map(Named, layernames(strat), layers(strat))) + Base.keys(strat::Stratigraphy) = layernames(strat) -Base.values(strat::Stratigraphy) = layers(strat) +Base.values(strat::Stratigraphy) = values(layers(strat)) Base.propertynames(strat::Stratigraphy) = Base.keys(strat) -Base.getproperty(strat::Stratigraphy, sym::Symbol) = strat[Val{sym}()].val -Base.getindex(strat::Stratigraphy, sym::Symbol) = strat[Val{sym}()].val -@generated Base.getindex(strat::Stratigraphy{N,TC}, ::Val{sym}) where {N,TC,sym} = :(layers(strat)[$(findfirst(T -> layername(T) == sym, TC.parameters))]) +Base.getproperty(strat::Stratigraphy, name::Symbol) = getproperty(layers(strat), name) +Base.getindex(strat::Stratigraphy, name::Symbol) = getproperty(strat, name) # Array and iteration overrides Base.size(strat::Stratigraphy) = size(layers(strat)) Base.length(strat::Stratigraphy) = length(layers(strat)) @@ -79,23 +81,20 @@ Base.lastindex(strat::Stratigraphy) = length(strat) Base.iterate(strat::Stratigraphy) = (layers(strat)[1],layers(strat)[2:end]) Base.iterate(strat::Stratigraphy, itrstate::Tuple) = (itrstate[1],itrstate[2:end]) Base.iterate(strat::Stratigraphy, itrstate::Tuple{}) = nothing -Base.show(io::IO, strat::Stratigraphy) = print(io, "Stratigraphy($(prod("$b => $n" for (n,b) in zip(layers(strat), boundaries(strat))))") -Base.NamedTuple(strat::Stratigraphy) = (; map(named_layer -> nameof(named_layer) => named_layer, layers(strat))...) -# ConstructionBase -ConstructionBase.getproperties(strat::Stratigraphy) = (;map(Pair, Base.keys(strat), Base.values(strat))...) -function ConstructionBase.setproperties(strat::Stratigraphy, patch::NamedTuple) - layers_patched = map(layers(strat)) do layer - Named(layername(layer), get(patch, layername(layer), layer.val)) +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") end - return Stratigraphy(boundaries(strat), layers_patched) end function Numerics.makegrid(strat::Stratigraphy, strategy::DiscretizationStrategy) strat_grid = nothing - for (bounds, named_layer) in zip(boundarypairs(strat)[2:end-1], layers(strat)[2:end-1]) + for (bounds, named_layer) in zip(boundarypairs(strat)[2:end-1], namedlayers(strat)[2:end-1]) @assert bounds[2] - bounds[1] > zero(bounds[1]) "Subsurface layers must have thickness greater than zero in the stratigraphy. The initial thickness can be later updated in `initialcondition!`." - layer = named_layer.val - layer_grid = Numerics.makegrid(layer, strategy, bounds) + layer = named_layer + layer_grid = Numerics.makegrid(layer.val, strategy, bounds) if !isnothing(strat_grid) # check that grid edges line up at layer boundary @assert strat_grid[end] == layer_grid[1] "Upper boundary of layer $(nameof(named_layer)) does not match the previous layer." @@ -113,8 +112,8 @@ function CryoGrid.initialcondition!(strat::Stratigraphy, state, inits) # we can just loop over everything. # first invoke initialconditon! with initializers for i in 1:length(strat) - layerᵢ = strat[i].val - stateᵢ = getproperty(state, layername(strat[i])) + layerᵢ = strat[i] + stateᵢ = getproperty(state, layernames(strat)[i]) for init in tuplejoin(CryoGrid.initializers(layerᵢ), inits) if haskey(stateᵢ.states, varname(init)) CryoGrid.initialcondition!(init, layerᵢ, stateᵢ) @@ -123,10 +122,10 @@ function CryoGrid.initialcondition!(strat::Stratigraphy, state, inits) end # then invoke non-specific layer inits for i in 1:length(strat)-1 - layerᵢ = strat[i].val - layerᵢ₊₁ = strat[i+1].val - stateᵢ = getproperty(state, layername(strat[i])) - stateᵢ₊₁ = getproperty(state, layername(strat[i+1])) + layerᵢ = strat[i] + layerᵢ₊₁ = strat[i+1] + stateᵢ = getproperty(state, layernames(strat)[i]) + stateᵢ₊₁ = getproperty(state, layernames(strat)[i+1]) if i == 1 CryoGrid.initialcondition!(layerᵢ, stateᵢ) end @@ -136,14 +135,14 @@ function CryoGrid.initialcondition!(strat::Stratigraphy, state, inits) end function CryoGrid.resetfluxes!(strat::Stratigraphy, state) - fastiterate(layers(strat)) do named_layer - CryoGrid.resetfluxes!(named_layer.val, getproperty(state, layername(named_layer))) + fastiterate(namedlayers(strat)) do named_layer + CryoGrid.resetfluxes!(named_layer, getproperty(state, nameof(named_layer))) end end function CryoGrid.updatestate!(strat::Stratigraphy, state) fastiterate(layers(strat)) do named_layer - CryoGrid.updatestate!(named_layer.val, getproperty(state, layername(named_layer))) + CryoGrid.updatestate!(named_layer, getproperty(state, nameof(named_layer))) end end @@ -152,12 +151,12 @@ end Special implementation of `interact!` that iterates over each pair of layers in the stratigraphy which are adjacent and "active" based on the current `state`. `state` must have properties defined corresponding to the name of each layer such -that `getproperty(state, layername(strat[i]))` would return the appropriate state object for the i'th layer in the stratigraphy. +that `getproperty` would return the appropriate state object for the i'th layer in the stratigraphy. """ @generated function CryoGrid.interact!(strat::Stratigraphy{N}, state) where {N} expr = Expr(:block) # build expressions for checking whether each layer is active - is_active_exprs = map(i -> :(CryoGrid.isactive(strat[$i].val, getproperty(state, layername(strat[$i])))), tuple(1:N...)) + is_active_exprs = map(i -> :(CryoGrid.isactive(strat[$i], getproperty(state, layernames(strat)[$i]))), tuple(1:N...)) # header code; get layer names and evaluate `isactive` for each layer push!( expr.args, @@ -177,8 +176,8 @@ that `getproperty(state, layername(strat[i]))` would return the appropriate stat if j == i + 1 # always try to invoke interact! for adjacent layers quote - if caninteract(strat[$i].val, strat[$j].val, getproperty(state, names[$i]), getproperty(state, names[$j])) - interact!(strat[$i].val, strat[$j].val, getproperty(state, names[$i]), getproperty(state, names[$j])) + if caninteract(strat[$i], strat[$j], getproperty(state, names[$i]), getproperty(state, names[$j])) + interact!(strat[$i], strat[$j], getproperty(state, names[$i]), getproperty(state, names[$j])) end end else @@ -192,9 +191,9 @@ that `getproperty(state, layername(strat[i]))` would return the appropriate stat # if layers i and j are both active, and all imbetween layers are inactive, invoke f!; # note the splat syntax here: $(inactive_imbetween_layers...) simply expands the tuple of # expressions 'inactive_imbetween_layers' as arguments to `tuple`. - can_interact = caninteract(strat[$i].val, strat[$j].val, getproperty(state, names[$i]), getproperty(state, names[$j])) + can_interact = caninteract(strat[$i], strat[$j], getproperty(state, names[$i]), getproperty(state, names[$j])) if is_active[$i] && is_active[$j] && all(tuple($(inactive_imbetween_layers...))) && can_interact - interact!(strat[$i].val, strat[$j].val, getproperty(state, names[$i]), getproperty(state, names[$j])) + interact!(strat[$i], strat[$j], getproperty(state, names[$i]), getproperty(state, names[$j])) end end end @@ -206,20 +205,20 @@ that `getproperty(state, layername(strat[i]))` would return the appropriate stat end function CryoGrid.computefluxes!(strat::Stratigraphy, state) - fastiterate(layers(strat)) do named_layer - CryoGrid.computefluxes!(named_layer.val, getproperty(state, layername(named_layer))) + fastiterate(namedlayers(strat)) do named_layer + CryoGrid.computefluxes!(named_layer.val, getproperty(state, nameof(named_layer))) end end function CryoGrid.timestep(strat::Stratigraphy, state) - max_dts = fastmap(layers(strat)) do named_layer - CryoGrid.timestep(named_layer.val, getproperty(state, layername(named_layer))) + max_dts = fastmap(namedlayers(strat)) do named_layer + CryoGrid.timestep(named_layer.val, getproperty(state, nameof(named_layer))) end return minimum(max_dts) end # collecting/grouping components -CryoGrid.events(strat::Stratigraphy) = map(named_layer -> _addlayerfield(CryoGrid.events(named_layer.val), nameof(named_layer)), NamedTuple(strat)) +CryoGrid.events(strat::Stratigraphy) = (; map(named_layer -> nameof(named_layer) => _addlayerfield(CryoGrid.events(named_layer.val), nameof(named_layer)), namedlayers(strat))...) CryoGrid.variables(::Union{FixedVolume,DiagnosticVolume}) = ( Diagnostic(:Δz, Scalar, u"m", domain=0..Inf), @@ -233,18 +232,16 @@ CryoGrid.variables(::PrognosticVolume) = ( ) # collects all variables in the stratgriphy, returning a NamedTuple of variable sets. function CryoGrid.variables(strat::Stratigraphy) - strat_nt = NamedTuple(strat) - layervars = map(CryoGrid.variables, strat_nt) - return map(layervars, strat_nt) do vars, named_layer + layervars = map(_collectvars, layers(strat)) + return map(layervars, layers(strat)) do vars, layer ( vars..., - CryoGrid.variables(CryoGrid.Volume(named_layer.val))..., + CryoGrid.variables(CryoGrid.Volume(layer))..., ) end end # collects and validates all variables in a given layer -function CryoGrid.variables(@nospecialize(named_layer::NamedLayer)) - layer = named_layer.val +function _collectvars(@nospecialize(layer::Layer)) declared_vars = variables(layer) nested_vars = Flatten.flatten(layer, Flatten.flattenable, Var) all_vars = vcat(collect(declared_vars), collect(nested_vars)) @@ -292,3 +289,31 @@ function _addlayerfield(@nospecialize(obj), name::Symbol) return obj end end + +function _withnames(pairs::Tuple{Vararg{Pair}}) + # case 1: user specified name as :name => layer or Named + withname(l::Pair{Symbol,<:Layer}) = Named(l) + withname(l::NamedLayer) = l + # case 2: no name specified, autogenerate based on type name + withname(l::Layer) = Named(Symbol(lowercase(string(nameof(typeof(l))))), l) + # extract depths and layers + depths = map(first, pairs) + layers = map(last, pairs) + # get named layers + named_layers = map(withname, layers) + # group by name + grouped_layers = Utils.groupby(nameof, named_layers) + # iterate over each name group + for name in keys(grouped_layers) + if length(grouped_layers[name]) > 1 + # if there is more than one layer with this name, append an integer to deduplicate them. + grouped_layers[name] = map(enumerate(grouped_layers[name])) do (i, named_layer) + # rename layer with integer id to deduplicate + Named(Symbol(nameof(named_layer), i), named_layer.val) + end + end + end + return map(depths, named_layers) do depth, named_layer + depth => popfirst!(grouped_layers[nameof(named_layer)]) + end +end diff --git a/src/Tiles/tile.jl b/src/Tiles/tile.jl index ce6bba22..e3a428df 100644 --- a/src/Tiles/tile.jl +++ b/src/Tiles/tile.jl @@ -103,11 +103,11 @@ function Tile( strat = stripunits(strat) events = CryoGrid.events(strat) vars = CryoGrid.variables(strat) - layers = map(NamedTuple(strat)) do named_layer - _addlayerfield(named_layer, nameof(named_layer)) + layers = map(namedlayers(strat)) do named_layer + nameof(named_layer) => _addlayerfield(named_layer, nameof(named_layer)) end # rebuild stratigraphy with updated parameters - strat = Stratigraphy(boundaries(strat), Tuple(values(layers))) + strat = Stratigraphy(boundaries(strat), (;layers...)) # construct state variables states = _initstatevars(strat, grid, vars, cachetype, arraytype; chunk_size) if isempty(inits) @@ -378,8 +378,8 @@ Adds parameter information to all nested types in `tile` by recursively calling """ function CryoGrid.parameterize(tile::Tile) ctor = ConstructionBase.constructorof(typeof(tile)) - new_layers = map(tile.strat) do named_layer - name = layername(named_layer) + new_layers = map(namedlayers(tile.strat)) do named_layer + name = nameof(named_layer) layer = CryoGrid.parameterize(named_layer.val) Named(name, _addlayerfield(layer, name)) end diff --git a/src/Utils/types.jl b/src/Utils/types.jl index e552d9ba..37d7152f 100644 --- a/src/Utils/types.jl +++ b/src/Utils/types.jl @@ -25,6 +25,7 @@ struct Named{name,T} end Named(values::Pair{Symbol,T}) where T = Named(values[1], values[2]) Base.nameof(::Named{name}) where name = name +Base.NamedTuple(ns::Tuple{Vararg{<:Named}}) = (; map(n -> nameof(n) => n.val, ns)...) ConstructionBase.constructorof(::Type{<:Named{name}}) where name = val -> Named(name, val) """ NamedTupleWrapper diff --git a/src/problem.jl b/src/problem.jl index 9d52bf29..87357fde 100644 --- a/src/problem.jl +++ b/src/problem.jl @@ -143,9 +143,9 @@ function _makecallbacks(tile::Tile) isgridevent(::GridContinuousEvent) = true isgridevent(::Event) = false callbacks = [] - for (i,named_layer) in enumerate(tile.strat) + for (i,named_layer) in enumerate(namedlayers(tile.strat)) events = CryoGrid.events(named_layer.val) - layer_name = CryoGrid.layername(named_layer) + layer_name = nameof(named_layer) for ev in events # if ev is a GridContinuousEvent, and was already added in a different layer, skip it. # GridContinuousEvents are defined on the whole grid/domain and so do not need to be duplicated @@ -167,7 +167,7 @@ function _criterionfunc(::Val{layername}, ev::Event, i_layer::Int) where layerna du = Tiles.withaxes(get_du(integrator), tile), t = t, state = Tiles.getstate(Val{layername}(), tile, u, du, t, integrator.dt); - return criterion(ev, layer.val, state) + return criterion(ev, layer, state) end end end @@ -179,7 +179,7 @@ function _gridcriterionfunc(::Val{layername}, ev::Event) where layername du = Tiles.withaxes(get_du(integrator), tile), t = t, state = Tiles.getstate(Val{layername}(), tile, u, du, t, integrator.dt); - criterion!(view(out, Numerics.bounds(state.grid)), ev, layer.val, process, state) + criterion!(view(out, Numerics.bounds(state.grid)), ev, layer, process, state) end end end @@ -194,7 +194,7 @@ function _triggerfunc(::Val{layername}, ev::Event, trig::Union{Nothing,T}, i_lay du = Tiles.withaxes(get_du(integrator), tile), t = integrator.t, state = Tiles.getstate(Val{layername}(), tile, u, du, t, integrator.dt); - _invoke_trigger!(ev, trig, layer.val, state) + _invoke_trigger!(ev, trig, layer, state) end end end @@ -209,7 +209,7 @@ function _gridtriggerfunc(::Val{layername}, ev::GridContinuousEvent, grid::Grid, t = integrator.t state = Tiles.getstate(Val{layername}(), tile, u, du, t, integrator.dt) if event_idx ∈ Numerics.bounds(state.grid) - _invoke_trigger!(ev, T(nothing), layer.val, state) + _invoke_trigger!(ev, T(nothing), layer, state) break end end diff --git a/test/Physics/Surface/swb_tests.jl b/test/Physics/Surface/swb_tests.jl index 6effaf9b..53113663 100644 --- a/test/Physics/Surface/swb_tests.jl +++ b/test/Physics/Surface/swb_tests.jl @@ -2,7 +2,12 @@ testgrid = CryoGrid.Presets.DefaultGrid_2cm @testset "interact!" begin sub = TestGroundLayer(WaterBalance(BucketScheme())) - top = Top(SurfaceWaterBalance(rainfall=ConstantForcing(1e-6u"m/s", :rainfall))) + top = Top( + SurfaceWaterBalance( + ConstantForcing(1e-6u"m/s", :rainfall), + ConstantForcing(0.0u"m/s", :snowfall), + ) + ) stop = Diagnostics.build_dummy_state(testgrid, top, with_units=true) ssub = Diagnostics.build_dummy_state(testgrid, sub, with_units=true) # initialize fully saturated diff --git a/test/Tiles/stratigraphy_tests.jl b/test/Tiles/stratigraphy_tests.jl index 6414cd93..959b1674 100644 --- a/test/Tiles/stratigraphy_tests.jl +++ b/test/Tiles/stratigraphy_tests.jl @@ -8,23 +8,23 @@ include("../types.jl") bounds = (-1.0u"m",0.0u"m",10.0u"m") strat = Stratigraphy( bounds[1] => Top(TestBoundary()), - bounds[2] => :testground => TestGroundLayer(TestGroundProcess()), + bounds[2] => TestGroundLayer(TestGroundProcess()), bounds[3] => Bottom(TestBoundary()) ) # Check boundaries @test all([bounds[i] == b for (i,b) in enumerate(bounds)]) # Check iteration and layer types - for (i,named_layer) in enumerate(strat) + for (i,named_layer) in enumerate(namedlayers(strat)) if i == 1 - @test layername(named_layer) == :top + @test nameof(named_layer) == :top @test typeof(named_layer.val) <: Top @test typeof(named_layer.val.proc) <: TestBoundary elseif i == 2 - @test layername(named_layer) == :testground + @test nameof(named_layer) == :testgroundlayer @test typeof(named_layer.val) <: TestGroundLayer @test typeof(named_layer.val.proc) <: TestGroundProcess elseif i == 3 - @test layername(named_layer) == :bottom + @test nameof(named_layer) == :bottom @test typeof(named_layer.val) <: Bottom @test typeof(named_layer.val.proc) <: TestBoundary end @@ -34,29 +34,29 @@ include("../types.jl") bounds = (-1.0u"m",0.0u"m",10.0u"m",100.0u"m") strat = Stratigraphy( bounds[1] => Top(TestBoundary()),( - bounds[2] => :testground1 => TestGroundLayer(TestGroundProcess()), - bounds[3] => :testground2 => TestGroundLayer(TestGroundProcess()), + bounds[2] => TestGroundLayer(TestGroundProcess()), + bounds[3] => TestGroundLayer(TestGroundProcess()), ), bounds[4] => Bottom(TestBoundary()) ) # Check boundaries @test all([bounds[i] == b for (i,b) in enumerate(bounds)]) # Check iteration and layer types - for (i,named_layer) in enumerate(strat) + for (i,named_layer) in enumerate(namedlayers(strat)) if i == 1 - @test layername(named_layer) == :top + @test nameof(named_layer) == :top @test typeof(named_layer.val) <: Top @test typeof(named_layer.val.proc) <: TestBoundary elseif i == 2 - @test layername(named_layer) == :testground1 + @test nameof(named_layer) == :testgroundlayer1 @test typeof(named_layer.val) <: TestGroundLayer @test typeof(named_layer.val.proc) <: TestGroundProcess elseif i == 3 - @test layername(named_layer) == :testground2 + @test nameof(named_layer) == :testgroundlayer2 @test typeof(named_layer.val) <: TestGroundLayer @test typeof(named_layer.val.proc) <: TestGroundProcess elseif i == 4 - @test layername(named_layer) == :bottom + @test nameof(named_layer) == :bottom @test typeof(named_layer.val) <: Bottom @test typeof(named_layer.val.proc) <: TestBoundary end @@ -66,62 +66,64 @@ include("../types.jl") bounds = (-1.0u"m",0.0u"m",10.0u"m",100.0u"m") strat = Stratigraphy( bounds[1] => Top(Coupled(TestBoundary(), TestBoundary())),( - bounds[2] => :testground1 => TestGroundLayer(Coupled(TestGroundProcess(), TestGroundProcess())), - bounds[3] => :testground2 => TestGroundLayer(TestGroundProcess()), + bounds[2] => TestGroundLayer(Coupled(TestGroundProcess(), TestGroundProcess())), + bounds[3] => TestGroundLayer(TestGroundProcess()), ), bounds[4] => Bottom(Coupled(TestBoundary(), TestBoundary())) ) # Check boundaries @test all([bounds[i] == b for (i,b) in enumerate(bounds)]) # Check iteration and component types - for (i,named_layer) in enumerate(strat) + for (i,named_layer) in enumerate(namedlayers(strat)) if i == 1 - @test layername(named_layer) == :top + @test nameof(named_layer) == :top @test typeof(named_layer.val) <: Top @test typeof(named_layer.val.proc) <: CoupledProcesses{Tuple{TestBoundary,TestBoundary}} elseif i == 2 - @test layername(named_layer) == :testground1 + @test nameof(named_layer) == :testgroundlayer1 @test typeof(named_layer.val) <: TestGroundLayer @test typeof(named_layer.val.proc) <: CoupledProcesses{Tuple{TestGroundProcess,TestGroundProcess}} elseif i == 3 - @test layername(named_layer) == :testground2 + @test nameof(named_layer) == :testgroundlayer2 @test typeof(named_layer.val) <: TestGroundLayer @test typeof(named_layer.val.proc) <: TestGroundProcess elseif i == 4 - @test layername(named_layer) == :bottom + @test nameof(named_layer) == :bottom @test typeof(named_layer.val) <: Bottom @test typeof(named_layer.val.proc) <: CoupledProcesses{Tuple{TestBoundary,TestBoundary}} end end end bounds = (-1.0u"m",0.0u"m",10.0u"m",100.0u"m") - # Check that duplicated layer names throws an error - @test_throws AssertionError Stratigraphy( + # Check that layer names are automatically deduplicated + strat = Stratigraphy( bounds[1] => Top(TestBoundary()), ( - bounds[2] => :duplicatedname => TestGroundLayer(TestGroundProcess()), - bounds[3] => :duplicatedname => TestGroundLayer(TestGroundProcess()), + bounds[2] => Named(:duplicatedname, TestGroundLayer(TestGroundProcess())), + bounds[3] => Named(:duplicatedname, TestGroundLayer(TestGroundProcess())), ), bounds[4] => Bottom(TestBoundary()) ) + @test hasproperty(strat, :duplicatedname1) + @test hasproperty(strat, :duplicatedname2) # Check that mis-specified stratigraphies throw method errors @test_throws MethodError Stratigraphy( bounds[1] => Bottom(TestBoundary()), - bounds[2] => :testground => TestGroundLayer(TestGroundProcess()), + bounds[2] => TestGroundLayer(TestGroundProcess()), bounds[3] => Top(TestBoundary()) ) @test_throws MethodError Stratigraphy( bounds[1] => Top(TestBoundary()), - bounds[2] => :testground => TestGroundLayer(TestGroundProcess()), + bounds[2] => TestGroundLayer(TestGroundProcess()), ) @test_throws MethodError Stratigraphy( - bounds[1] => :testground => TestGroundLayer(TestGroundProcess()), + bounds[1] => TestGroundLayer(TestGroundProcess()), bounds[2] => Bottom(TestBoundary()) ) # Check that out-of-order bounds throw an error @test_throws AssertionError Stratigraphy( bounds[2] => Top(TestBoundary()), ( - bounds[3] => :duplicatedname => TestGroundLayer(TestGroundProcess()), - bounds[1] => :duplicatedname => TestGroundLayer(TestGroundProcess()), + bounds[3] => TestGroundLayer(TestGroundProcess()), + bounds[1] => TestGroundLayer(TestGroundProcess()), ), bounds[4] => Bottom(TestBoundary()) ) @@ -129,7 +131,7 @@ include("../types.jl") bounds = (0.0u"m",0.0u"m",10.0u"m") @test_nowarn strat = Stratigraphy( bounds[1] => Top(TestBoundary()), - bounds[2] => :testground => TestGroundLayer(TestGroundProcess()), + bounds[2] => TestGroundLayer(TestGroundProcess()), bounds[3] => Bottom(TestBoundary()) ) end diff --git a/test/Tiles/tile_tests.jl b/test/Tiles/tile_tests.jl index 642f07d3..d7abd695 100644 --- a/test/Tiles/tile_tests.jl +++ b/test/Tiles/tile_tests.jl @@ -11,7 +11,7 @@ include("../types.jl") grid = Grid(Vector(0.0:10.0:1000.0)u"m") strat = Stratigraphy( -1.0u"m" => Top(TestBoundary()), - 0.0u"m" => :testground => TestGroundLayer(TestGroundProcess()), + 0.0u"m" => Named(:testground, TestGroundLayer(TestGroundProcess())), 1000.0u"m" => Bottom(TestBoundary()) ) function checkfields(model) @@ -89,8 +89,8 @@ include("../types.jl") grid = Grid(Vector(0.0:10.0:1000.0)u"m") strat = Stratigraphy( -1.0u"m" => Top(TestBoundary()), ( - 0.0u"m" => :testground1 => TestGroundLayer(TestGroundProcess()), - 100.0u"m" => :testground2 => TestGroundLayer(TestGroundProcess()), + 0.0u"m" => Named(:testground1, TestGroundLayer(TestGroundProcess())), + 100.0u"m" => Named(:testground2, TestGroundLayer(TestGroundProcess())), ), 1000.0u"m" => Bottom(TestBoundary()) ) @@ -115,8 +115,8 @@ include("../types.jl") grid = Grid(Vector(0.0:10.0:1000.0)u"m") strat = Stratigraphy( -1.0u"m" => Top(TestBoundary()), ( - 0.0u"m" => :testground1 => TestGroundLayer(TestGroundProcess()), - 100.0u"m" => :testground2 => TestGroundLayer(TestGroundProcess()), + 0.0u"m" => Named(:testgroundlayer1, TestGroundLayer(TestGroundProcess())), + 100.0u"m" => Named(:testgroundlayer2, TestGroundLayer(TestGroundProcess())), ), 1000.0u"m" => Bottom(TestBoundary()) ) From be1642e10d2f98d472e4b86d2c184445f5a63f97 Mon Sep 17 00:00:00 2001 From: Brian Groenke Date: Sat, 30 Sep 2023 00:39:28 +0200 Subject: [PATCH 2/8] Fix lingering issues from refactoring --- examples/07_heat_freeW_lite_implicit.jl | 1 + src/Tiles/stratigraphy.jl | 10 +++++----- src/Tiles/tile.jl | 6 +++--- 3 files changed, 9 insertions(+), 8 deletions(-) diff --git a/examples/07_heat_freeW_lite_implicit.jl b/examples/07_heat_freeW_lite_implicit.jl index d7ca206c..fa02a43a 100644 --- a/examples/07_heat_freeW_lite_implicit.jl +++ b/examples/07_heat_freeW_lite_implicit.jl @@ -51,6 +51,7 @@ Plots.plot(out.T[Z(Near(zs))], color=cg[LinRange(0.0,1.0,length(zs))]', ylabel=" # 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. +using OrdinaryDiffEq sol2 = @time solve(prob, ImplicitEuler(nlsolve=NLCGLite(max_iter=1000)), adaptive=false, dt=24*3600.0) ## 3rd order additive scheme from Kennedy and Alan 2001 sol3 = @time solve(prob, KenCarp3(nlsolve=NLCGLite(max_iter=1000)), adaptive=false, dt=24*3600.0) diff --git a/src/Tiles/stratigraphy.jl b/src/Tiles/stratigraphy.jl index d8842426..02fa921a 100644 --- a/src/Tiles/stratigraphy.jl +++ b/src/Tiles/stratigraphy.jl @@ -64,8 +64,8 @@ boundaries(strat::Stratigraphy) = getfield(strat, :boundaries) boundarypairs(strat::Stratigraphy) = boundarypairs(boundaries(strat)) boundarypairs(bounds::NTuple) = tuple(map(tuple, bounds[1:end-1], bounds[2:end])..., (bounds[end], bounds[end])) layernames(strat::Stratigraphy) = keys(layers(strat)) -layertypes(::Type{<:Stratigraphy{N,TLayers}}) where {N,TLayers} = map(layertype, TLayers.parameters) -namedlayers(strat::Stratigraphy) = Tuple(map(Named, layernames(strat), layers(strat))) +layertypes(::Type{<:Stratigraphy{N,NamedTuple{names,TLayers}}}) where {N,TLayers,names} = Tuple(TLayers.parameters) +namedlayers(strat::Stratigraphy) = map(Named, layernames(strat), values(layers(strat))) Base.keys(strat::Stratigraphy) = layernames(strat) Base.values(strat::Stratigraphy) = values(layers(strat)) @@ -136,13 +136,13 @@ end function CryoGrid.resetfluxes!(strat::Stratigraphy, state) fastiterate(namedlayers(strat)) do named_layer - CryoGrid.resetfluxes!(named_layer, getproperty(state, nameof(named_layer))) + CryoGrid.resetfluxes!(named_layer.val, getproperty(state, nameof(named_layer))) end end function CryoGrid.updatestate!(strat::Stratigraphy, state) - fastiterate(layers(strat)) do named_layer - CryoGrid.updatestate!(named_layer, getproperty(state, nameof(named_layer))) + fastiterate(namedlayers(strat)) do named_layer + CryoGrid.updatestate!(named_layer.val, getproperty(state, nameof(named_layer))) end end diff --git a/src/Tiles/tile.jl b/src/Tiles/tile.jl index e3a428df..79947e21 100644 --- a/src/Tiles/tile.jl +++ b/src/Tiles/tile.jl @@ -104,7 +104,7 @@ function Tile( events = CryoGrid.events(strat) vars = CryoGrid.variables(strat) layers = map(namedlayers(strat)) do named_layer - nameof(named_layer) => _addlayerfield(named_layer, nameof(named_layer)) + nameof(named_layer) => _addlayerfield(named_layer.val, nameof(named_layer)) end # rebuild stratigraphy with updated parameters strat = Stratigraphy(boundaries(strat), (;layers...)) @@ -289,7 +289,7 @@ end function initboundaries!(tile::Tile{TStrat}, u) where {TStrat} bounds = boundarypairs(tile.strat) - map(bounds, layers(tile.strat)) do (z1, z2), named_layer + map(bounds, namedlayers(tile.strat)) do (z1, z2), named_layer name = nameof(named_layer) diag_layer = getproperty(tile.state.diag, name) z = retrieve(diag_layer.z, u) @@ -303,7 +303,7 @@ end function update_layer_boundaries(tile::Tile, u) # calculate grid boundaries starting from the bottom moving up to the surface zbot = tile.state.grid[end] - return accumulate(reverse(layers(tile.strat)); init=zbot) do z_acc, named_layer + return accumulate(reverse(namedlayers(tile.strat)); init=zbot) do z_acc, named_layer name = nameof(named_layer) diag_layer = getproperty(tile.state.diag, name) z_state = retrieve(diag_layer.z, u) From 7f87a20a966333427a96f1ce39c9a9934c23a30d Mon Sep 17 00:00:00 2001 From: Brian Groenke Date: Sat, 30 Sep 2023 11:43:48 +0200 Subject: [PATCH 3/8] Rename heat operator types --- docs/src/manual/overview.md | 4 +-- examples/02_heat_sfcc_constantbc.jl | 2 +- examples/03_heat_sfcc_samoylov.jl | 2 +- examples/07_heat_freeW_lite_implicit.jl | 4 +-- examples/08_heat_sfcc_richardseq_samoylov.jl | 2 +- src/Physics/Heat/heat_conduction.jl | 24 ++++++++-------- src/Physics/Heat/heat_methods.jl | 2 +- src/Physics/Heat/heat_types.jl | 30 ++++++++++---------- src/Physics/Salt/Salt.jl | 2 +- src/Physics/Salt/salt_diffusion.jl | 6 ++-- src/Physics/Snow/snow_bulk.jl | 12 ++++---- src/Physics/Soils/Soils.jl | 4 +-- src/Physics/Soils/soil_heat.jl | 12 ++++---- src/Physics/Soils/soil_water_heat.jl | 6 ++-- src/Solvers/DiffEq/DiffEq.jl | 2 ++ src/Solvers/DiffEq/ode_solvers.jl | 1 - src/Solvers/Solvers.jl | 1 + test/Solvers/cglite_implicit_tests.jl | 4 +-- test/test_problems.jl | 2 +- 19 files changed, 62 insertions(+), 60 deletions(-) diff --git a/docs/src/manual/overview.md b/docs/src/manual/overview.md index 493f6e01..37a36698 100644 --- a/docs/src/manual/overview.md +++ b/docs/src/manual/overview.md @@ -55,7 +55,7 @@ julia> out.T Notice that, in the example above, it is types such as `Ground`, `HeatBalance`, `DallAmico`, etc. that specify which components the model should use. These components are defined by adding method dispatches to the [CryoGrid interface](@ref toplevel) methods. State variables are declared via the [`variables`](@ref) method, e.g: ```julia -variables(soil::Soil, heat::HeatBalance{<:Enthalpy}) = ( +variables(soil::Soil, heat::HeatBalance{<:EnthalpyBased}) = ( Prognostic(:H, OnGrid(Cells), u"J/m^3"), Diagnostic(:T, OnGrid(Cells), u"°C"), Diagnostic(:C, OnGrid(Cells), u"J//K*/m^3"), @@ -74,7 +74,7 @@ The real work finally happens in [`updatestate!`](@ref) and [`computefluxes!`](@ 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): ```julia -function CryoGrid.computefluxes!(::SubSurface, ::HeatBalance{<:FreezeCurve,<:Enthalpy}, state) +function CryoGrid.computefluxes!(::SubSurface, ::HeatBalance{<:FreezeCurve,<: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/examples/02_heat_sfcc_constantbc.jl b/examples/02_heat_sfcc_constantbc.jl index d72cd660..918e45b5 100644 --- a/examples/02_heat_sfcc_constantbc.jl +++ b/examples/02_heat_sfcc_constantbc.jl @@ -30,7 +30,7 @@ Plots.plot(-2.0u"°C":0.01u"K":0.0u"°C", sfcc) # Enthalpy form of the heat transfer operator (i.e. prognostic :H). In this case, this is equivalent to # the shorthand `SoilHeatTile(:H, ...)`. However, it's worth demonstrating how the operator can be explicitly # specified. -heatop = Heat.EnthalpyForm(SFCCPreSolver()) +heatop = Heat.MOLEnthalpy(SFCCPreSolver()) initT = initializer(:T, tempprofile) tile = CryoGrid.Presets.SoilHeatTile( heatop, diff --git a/examples/03_heat_sfcc_samoylov.jl b/examples/03_heat_sfcc_samoylov.jl index f6fad1eb..4cd459cf 100644 --- a/examples/03_heat_sfcc_samoylov.jl +++ b/examples/03_heat_sfcc_samoylov.jl @@ -4,7 +4,7 @@ # forced using air temperatures from Samoylov Island. We use the # SFCC formulation of Painter and Karra (2014). For the purpose # of demonstration, we use the apparent heat capacity form of the -# heat equation in this example (i.e. [`Heat.TemperatureForm`](@ref)). +# heat equation in this example (i.e. [`Heat.MOLTemperature`](@ref)). using CryoGrid diff --git a/examples/07_heat_freeW_lite_implicit.jl b/examples/07_heat_freeW_lite_implicit.jl index fa02a43a..da2db198 100644 --- a/examples/07_heat_freeW_lite_implicit.jl +++ b/examples/07_heat_freeW_lite_implicit.jl @@ -35,11 +35,11 @@ tile = Tile(strat, modelgrid, initT); # 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 5 years. -tspan = (DateTime(2010,12,30), DateTime(2011,12,30)) +tspan = (DateTime(2010,1,1), DateTime(2015,1,1)) tspan_sol = convert_tspan(tspan) u0, du0 = initialcondition!(tile, tspan); prob = CryoGridProblem(tile, u0, tspan, saveat=24*3600.0, savevars=(:T,), step_limiter=nothing) -sol = @time solve(prob, LiteImplicitEuler(), dt=3*3600.0) +sol = @time solve(prob, LiteImplicitEuler(), dt=24*3600.0) out = CryoGridOutput(sol) # Plot the results! diff --git a/examples/08_heat_sfcc_richardseq_samoylov.jl b/examples/08_heat_sfcc_richardseq_samoylov.jl index 4be164af..226ff474 100644 --- a/examples/08_heat_sfcc_richardseq_samoylov.jl +++ b/examples/08_heat_sfcc_richardseq_samoylov.jl @@ -23,7 +23,7 @@ sfcc = PainterKarra(ω=0.0, swrc=swrc) waterflow = RichardsEq(swrc=swrc); # We use the enthalpy-based heat diffusion with high accuracy Newton-based solver for inverse enthalpy mapping -heatop = Heat.EnthalpyForm(SFCCNewtonSolver()) +heatop = Heat.MOLEnthalpy(SFCCNewtonSolver()) upperbc = WaterHeatBC(SurfaceWaterBalance(rainfall=forcings.rainfall), TemperatureGradient(forcings.Tair, NFactor(nf=0.6, nt=0.9))); # We will use a simple stratigraphy with three subsurface soil layers. diff --git a/src/Physics/Heat/heat_conduction.jl b/src/Physics/Heat/heat_conduction.jl index 8906a3e9..e978a345 100644 --- a/src/Physics/Heat/heat_conduction.jl +++ b/src/Physics/Heat/heat_conduction.jl @@ -9,7 +9,7 @@ function heatflux(T₁, T₂, k₁, k₂, Δ₁, Δ₂, z₁, z₂) end # Free water freeze curve -@propagate_inbounds function enthalpyinv(sub::SubSurface, heat::HeatBalance{FreeWater,<:Enthalpy}, state, i) +@propagate_inbounds function enthalpyinv(sub::SubSurface, heat::HeatBalance{FreeWater,<:EnthalpyBased}, state, i) θwi = Hydrology.watercontent(sub, state, i) H = state.H[i] L = heat.prop.L @@ -35,14 +35,14 @@ end end """ - freezethaw!(sub::SubSurface, heat::HeatBalance{FreeWater,<:Enthalpy}, state) + freezethaw!(sub::SubSurface, heat::HeatBalance{FreeWater,<:EnthalpyBased}, state) Implementation of "free water" freezing characteristic for any subsurface layer. Assumes that `state` contains at least temperature (T), enthalpy (H), heat capacity (C), total water content (θwi), and liquid water content (θw). """ freezethaw!(sub, state) = freezethaw!(sub, processes(sub), state) -function freezethaw!(sub::SubSurface, heat::HeatBalance{FreeWater,<:Enthalpy}, state) +function freezethaw!(sub::SubSurface, heat::HeatBalance{FreeWater,<:EnthalpyBased}, state) @inbounds for i in 1:length(state.H) # update T, θw, C state.T[i], state.θw[i], state.C[i] = enthalpyinv(sub, heat, state, i) @@ -68,7 +68,7 @@ heatvariables(::HeatBalance) = ( """ Variable definitions for heat conduction (enthalpy) on any SubSurface layer. """ -CryoGrid.variables(heat::HeatBalance{<:FreezeCurve,<:Enthalpy}) = ( +CryoGrid.variables(heat::HeatBalance{<:FreezeCurve,<:EnthalpyBased}) = ( Prognostic(:H, OnGrid(Cells), u"J/m^3"), Diagnostic(:T, OnGrid(Cells), u"°C"), heatvariables(heat)..., @@ -76,7 +76,7 @@ CryoGrid.variables(heat::HeatBalance{<:FreezeCurve,<:Enthalpy}) = ( """ Variable definitions for heat conduction (temperature). """ -CryoGrid.variables(heat::HeatBalance{<:FreezeCurve,<:Temperature}) = ( +CryoGrid.variables(heat::HeatBalance{<:FreezeCurve,<:TemperatureBased}) = ( Prognostic(:T, OnGrid(Cells), u"°C"), Diagnostic(:H, OnGrid(Cells), u"J/m^3"), Diagnostic(:∂H∂t, OnGrid(Cells), u"W/m^3"), @@ -107,7 +107,7 @@ function CryoGrid.interact!(sub1::SubSurface, ::HeatBalance, sub2::SubSurface, : return nothing end -function CryoGrid.computefluxes!(::SubSurface, ::HeatBalance{<:FreezeCurve,<:Enthalpy}, state) +function CryoGrid.computefluxes!(::SubSurface, ::HeatBalance{<:FreezeCurve,<: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; @@ -118,7 +118,7 @@ function CryoGrid.computefluxes!(::SubSurface, ::HeatBalance{<:FreezeCurve,<:Ent divergence!(state.∂H∂t, state.jH, Δk) return nothing end -function CryoGrid.computefluxes!(sub::SubSurface, ::HeatBalance{<:FreezeCurve,<:Temperature}, state) +function CryoGrid.computefluxes!(sub::SubSurface, ::HeatBalance{<:FreezeCurve,<: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 @@ -134,7 +134,7 @@ function CryoGrid.computefluxes!(sub::SubSurface, ::HeatBalance{<:FreezeCurve,<: end # CFL not defined for free-water freeze curve -CryoGrid.timestep(::SubSurface, heat::HeatBalance{FreeWater,<:Enthalpy,<:CryoGrid.CFL}, state) = Inf +CryoGrid.timestep(::SubSurface, heat::HeatBalance{FreeWater,<:EnthalpyBased,<:CryoGrid.CFL}, state) = Inf """ timestep(::SubSurface, ::HeatBalance{Tfc,THeatOp,CFL}, state) where {Tfc,THeatOp} @@ -143,10 +143,10 @@ defined as: Δt_max = u*Δx^2, where`u` is the "characteristic velocity" which h is taken to be the diffusivity: `∂H∂T / kc`. """ function CryoGrid.timestep(::SubSurface, heat::HeatBalance{Tfc,THeatOp,<:CryoGrid.CFL}, state) where {Tfc,THeatOp} - derivative(::Enthalpy, state) = state.∂H∂t - derivative(::Temperature, state) = state.∂T∂t - prognostic(::Enthalpy, state) = state.H - prognostic(::Temperature, state) = state.T + derivative(::EnthalpyBased, state) = state.∂H∂t + derivative(::TemperatureBased, state) = state.∂T∂t + prognostic(::EnthalpyBased, state) = state.H + prognostic(::TemperatureBased, state) = state.T Δx = Δ(state.grid) dtmax = Inf @inbounds for i in eachindex(Δx) diff --git a/src/Physics/Heat/heat_methods.jl b/src/Physics/Heat/heat_methods.jl index 98c41802..702e84a5 100644 --- a/src/Physics/Heat/heat_methods.jl +++ b/src/Physics/Heat/heat_methods.jl @@ -76,6 +76,6 @@ heatcapacity(op::HeatOperator) = op.hc Retreives the nonlinear solver for the enthalpy/temperature constitutive relation, if defined. """ fcsolver(::HeatOperator) = nothing -fcsolver(op::Enthalpy) = op.fcsolver +fcsolver(op::EnthalpyBased) = op.fcsolver fcsolver(heat::HeatBalance) = fcsolver(heat.op) fcsolver(::Nothing) = nothing diff --git a/src/Physics/Heat/heat_types.jl b/src/Physics/Heat/heat_types.jl index 96b98a78..c8dbfa3b 100644 --- a/src/Physics/Heat/heat_types.jl +++ b/src/Physics/Heat/heat_types.jl @@ -7,15 +7,15 @@ abstract type HeatOperator{progvar} end """ Type alias for `HeatOperator{:H}`, i.e. enthalpy-based heat transfer operators. """ -const Enthalpy = HeatOperator{:H} +const EnthalpyBased = HeatOperator{:H} """ Type alias for `HeatOperator{:T}`, i.e. temperature-based heat transfer operators. """ -const Temperature = HeatOperator{:T} +const TemperatureBased = HeatOperator{:T} # default step limiters -default_dtlim(::Temperature) = CryoGrid.CFL(maxdelta=CryoGrid.MaxDelta(Inf)) -default_dtlim(::Enthalpy) = CryoGrid.MaxDelta(50u"kJ") +default_dtlim(::TemperatureBased) = CryoGrid.CFL(maxdelta=CryoGrid.MaxDelta(Inf)) +default_dtlim(::EnthalpyBased) = CryoGrid.MaxDelta(50u"kJ") default_dtlim(::HeatOperator) = nothing # default freezecurve solvers @@ -31,7 +31,7 @@ the `HeatOperator`, `op`. """ Base.@kwdef struct HeatBalance{Tfc<:FreezeCurve,THeatOp<:HeatOperator,Tdt,Tprop} <: SubSurfaceProcess freezecurve::Tfc = FreeWater() - op::THeatOp = EnthalpyForm(default_fcsolver(freezecurve)) + op::THeatOp = MOLEnthalpy(default_fcsolver(freezecurve)) prop::Tprop = HeatBalanceProperties() dtlim::Tdt = default_dtlim(op) # timestep limiter function HeatBalance(freezecurve, op, prop, dtlim) @@ -42,39 +42,39 @@ Base.@kwdef struct HeatBalance{Tfc<:FreezeCurve,THeatOp<:HeatOperator,Tdt,Tprop} end # convenience constructors for HeatBalance HeatBalance(var::Symbol; kwargs...) = HeatBalance(Val{var}(); kwargs...) -HeatBalance(::Val{:H}; freezecurve::FreezeCurve=FreeWater(), fcsolver=default_fcsolver(freezecurve), kwargs...) = HeatBalance(; op=EnthalpyForm(fcsolver), freezecurve, kwargs...) -HeatBalance(::Val{:T}; freezecurve::FreezeCurve, kwargs...) = HeatBalance(; op=TemperatureForm(), freezecurve, kwargs...) +HeatBalance(::Val{:H}; freezecurve::FreezeCurve=FreeWater(), fcsolver=default_fcsolver(freezecurve), kwargs...) = HeatBalance(; op=MOLEnthalpy(fcsolver), freezecurve, kwargs...) +HeatBalance(::Val{:T}; freezecurve::FreezeCurve, kwargs...) = HeatBalance(; op=MOLTemperature(), freezecurve, kwargs...) HeatBalance(op::HeatOperator; kwargs...) = HeatBalance(; op, kwargs...) # validation of HeatBalance freezecurve/operator configuration _validate_heat_config(::FreezeCurve, ::HeatOperator) = nothing # do nothing when valid -_validate_heat_config(::FreeWater, ::Temperature) = error("Invalid heat balance configuration; temperature formulations of the heat operator are not compatible with the free water freeze curve.") +_validate_heat_config(::FreeWater, ::TemperatureBased) = error("Invalid heat balance configuration; temperature formulations of the heat operator are not compatible with the free water freeze curve.") # Heat operators """ - TemperatureForm{Tcond,Thc} <: HeatOperator{:T} + MOLTemperature{Tcond,Thc} <: HeatOperator{:T} Represents a standard method-of-lines (MOL) forward diffusion operator for heat conduction with temperature `T` as the prognostic variable. The time derivative is scaled by the reciprocal of the apparent heat capacity `dH/dT` to account for latent heat effects due to phase change. """ -struct TemperatureForm{Tcond,Thc} <: HeatOperator{:T} +struct MOLTemperature{Tcond,Thc} <: HeatOperator{:T} cond::Tcond hc::Thc - TemperatureForm(cond=quadratic_parallel_conductivity, hc=weighted_average_heatcapacity) = new{typeof(cond),typeof(hc)}(cond, hc) + MOLTemperature(cond=quadratic_parallel_conductivity, hc=weighted_average_heatcapacity) = new{typeof(cond),typeof(hc)}(cond, hc) end """ - EnthalpyForm{Tsolver,Tcond,Thc} <: HeatOperator{:H} + MOLEnthalpy{Tsolver,Tcond,Thc} <: HeatOperator{:H} Represents a standard method-of-lines (MOL) forward diffusion operator for heat conduction with enthalpy `H` as the prognostic variable and a nonlinear solver for resolving the inverse enthalpy -> temperature mapping when applicable. This formulation should generally be preferred -over `TemperatureForm` since it is energy-conserving and embeds the latent heat storage directly +over `MOLTemperature` since it is energy-conserving and embeds the latent heat storage directly in the prognostic state. """ -struct EnthalpyForm{Tsolver,Tcond,Thc} <: HeatOperator{:H} +struct MOLEnthalpy{Tsolver,Tcond,Thc} <: HeatOperator{:H} fcsolver::Tsolver cond::Tcond hc::Thc - EnthalpyForm(fcsolver=nothing, cond=quadratic_parallel_conductivity, hc=weighted_average_heatcapacity) = new{typeof(fcsolver),typeof(cond),typeof(hc)}(fcsolver, cond, hc) + MOLEnthalpy(fcsolver=nothing, cond=quadratic_parallel_conductivity, hc=weighted_average_heatcapacity) = new{typeof(fcsolver),typeof(cond),typeof(hc)}(fcsolver, cond, hc) end """ EnthalpyImplicit <: HeatOperator{:H} diff --git a/src/Physics/Salt/Salt.jl b/src/Physics/Salt/Salt.jl index b71a17a0..643b475c 100644 --- a/src/Physics/Salt/Salt.jl +++ b/src/Physics/Salt/Salt.jl @@ -2,7 +2,7 @@ module Salt using CryoGrid using CryoGrid.Heat -using CryoGrid.Heat: Temperature +using CryoGrid.Heat: TemperatureBased using CryoGrid.Hydrology using CryoGrid.Numerics using CryoGrid.Soils diff --git a/src/Physics/Salt/salt_diffusion.jl b/src/Physics/Salt/salt_diffusion.jl index f9364103..a600699d 100644 --- a/src/Physics/Salt/salt_diffusion.jl +++ b/src/Physics/Salt/salt_diffusion.jl @@ -4,7 +4,7 @@ function Heat.freezethaw!( soil::SalineGround, ps::CoupledHeatSalt{THeat}, state -) where {THeat<:HeatBalance{<:DallAmicoSalt,<:Temperature}} +) where {THeat<:HeatBalance{<:DallAmicoSalt,<:TemperatureBased}} salt, heat = ps sfcc = heat.freezecurve thermalprops = Heat.thermalproperties(soil) @@ -54,7 +54,7 @@ function CryoGrid.updatestate!( soil::SalineGround, ps::CoupledHeatSalt{THeat}, state -) where {THeat<:HeatBalance{<:SFCC,<:Temperature}} +) where {THeat<:HeatBalance{<:SFCC,<:TemperatureBased}} # Reset energy flux to zero; this is redundant when H is the prognostic variable # but necessary when it is not. salt, heat = ps @@ -130,7 +130,7 @@ function CryoGrid.computefluxes!( ::SalineGround, ps::CoupledHeatSalt{THeat}, state -) where {THeat<:HeatBalance{<:SFCC,<:Temperature}} +) where {THeat<:HeatBalance{<:SFCC,<:TemperatureBased}} salt, heat = ps #read current state diff --git a/src/Physics/Snow/snow_bulk.jl b/src/Physics/Snow/snow_bulk.jl index a8f5de0d..e1a057dd 100644 --- a/src/Physics/Snow/snow_bulk.jl +++ b/src/Physics/Snow/snow_bulk.jl @@ -16,8 +16,8 @@ end Type alias for Snowpack with `Bulk` parameterization. """ const BulkSnowpack{T} = Snowpack{<:Bulk{T}} where {T} -# Local alias for Heat Enthalpy type -const Enthalpy = Heat.Enthalpy +# Local alias for Heat EnthalpyBased type +const EnthalpyBased = Heat.EnthalpyBased threshold(snow::BulkSnowpack) = snow.para.thresh @@ -159,7 +159,7 @@ end # computefluxes! for free water, enthalpy based HeatBalance on bulk snow layer function CryoGrid.computefluxes!( snow::BulkSnowpack, - ps::Coupled(SnowMassBalance, HeatBalance{FreeWater,<:Enthalpy}), + ps::Coupled(SnowMassBalance, HeatBalance{FreeWater,<:EnthalpyBased}), state ) mass, heat = ps @@ -175,7 +175,7 @@ function CryoGrid.computefluxes!( end function CryoGrid.computefluxes!( snow::BulkSnowpack, - ps::Coupled(SnowMassBalance, WaterBalance, HeatBalance{FreeWater,<:Enthalpy}), + ps::Coupled(SnowMassBalance, WaterBalance, HeatBalance{FreeWater,<:EnthalpyBased}), state ) mass, water, heat = ps @@ -233,7 +233,7 @@ function CryoGrid.updatestate!( procs::Coupled( DynamicSnowMassBalance{TAcc,<:DegreeDayMelt}, WaterBalance, - HeatBalance{FreeWater,<:Enthalpy} + HeatBalance{FreeWater,<:EnthalpyBased} ), state ) where {TAcc} @@ -304,7 +304,7 @@ function CryoGrid.trigger!( end function CryoGrid.updatestate!( snow::BulkSnowpack, - procs::Coupled(PrescribedSnowMassBalance,HeatBalance{FreeWater,<:Enthalpy}), + procs::Coupled(PrescribedSnowMassBalance,HeatBalance{FreeWater,<:EnthalpyBased}), state ) mass, heat = procs diff --git a/src/Physics/Soils/Soils.jl b/src/Physics/Soils/Soils.jl index 751817a0..5e2aa337 100644 --- a/src/Physics/Soils/Soils.jl +++ b/src/Physics/Soils/Soils.jl @@ -22,8 +22,8 @@ using UnPack import ConstructionBase # aliases for heat formulations in Heat module -const Temperature = Heat.Temperature -const Enthalpy = Heat.Enthalpy +const TemperatureBased = Heat.TemperatureBased +const EnthalpyBased = Heat.EnthalpyBased const EnthalpyImplicit = Heat.EnthalpyImplicit export Ground, AbstractGround diff --git a/src/Physics/Soils/soil_heat.jl b/src/Physics/Soils/soil_heat.jl index fe4494bb..2edc460c 100644 --- a/src/Physics/Soils/soil_heat.jl +++ b/src/Physics/Soils/soil_heat.jl @@ -35,7 +35,7 @@ function sfccsolver!(soil::Soil, heat::HeatBalance{<:SFCC,TOp}, state) where {TO sat = saturation(soil, state) θsat = porosity(soil, state) hc = sfccheatcap(soil, heat, state, 1) - if TOp <: Enthalpy + if TOp <: EnthalpyBased solver = Heat.fcsolver(heat) @assert !isnothing(solver) "SFCC solver must be provided in HeatBalance operator. Check the model configuration." FreezeCurves.initialize!(solver, fc, hc; sat, θsat) @@ -111,14 +111,14 @@ function CryoGrid.initialcondition!(soil::Soil, heat::HeatBalance{FreeWater}, st end """ - freezethaw!(soil::Soil, heat::HeatBalance{<:SFCC,<:Temperature}, state) - freezethaw!(soil::Soil, heat::HeatBalance{<:SFCC,<:Enthalpy}, state) + freezethaw!(soil::Soil, heat::HeatBalance{<:SFCC,<:TemperatureBased}, state) + freezethaw!(soil::Soil, heat::HeatBalance{<:SFCC,<:EnthalpyBased}, state) Updates state variables according to the specified SFCC function and solver. For heat conduction with enthalpy, evaluation of the inverse enthalpy function is performed using the given solver. For heat conduction with temperature, we can simply evaluate the freeze curve to get dHdT, θw, and H. """ -function Heat.freezethaw!(soil::Soil, heat::HeatBalance{<:SFCC,<:Temperature}, state) +function Heat.freezethaw!(soil::Soil, heat::HeatBalance{<:SFCC,<:TemperatureBased}, state) sfcc = heat.freezecurve L = heat.prop.L @unpack ch_w, ch_i = Heat.thermalproperties(soil) @@ -135,7 +135,7 @@ function Heat.freezethaw!(soil::Soil, heat::HeatBalance{<:SFCC,<:Temperature}, s end # freezethaw! implementation for enthalpy and implicit enthalpy formulations -function Heat.freezethaw!(soil::Soil, heat::HeatBalance{<:SFCC,<:Enthalpy}, state) +function Heat.freezethaw!(soil::Soil, heat::HeatBalance{<:SFCC,<:EnthalpyBased}, state) sfcc = heat.freezecurve solver = Heat.fcsolver(heat) @inbounds for i in 1:length(state.H) @@ -162,7 +162,7 @@ function Heat.freezethaw!(soil::Soil, heat::HeatBalance{<:SFCC,<:Enthalpy}, stat end end -function Heat.enthalpyinv(soil::Soil, heat::HeatBalance{<:SFCC,<:Enthalpy}, state, i) +function Heat.enthalpyinv(soil::Soil, heat::HeatBalance{<:SFCC,<:EnthalpyBased}, state, i) sfcc = heat.freezecurve solver = Heat.fcsolver(heat) @inbounds let H = state.H[i], # enthalpy diff --git a/src/Physics/Soils/soil_water_heat.jl b/src/Physics/Soils/soil_water_heat.jl index 38e10a7c..9079117d 100644 --- a/src/Physics/Soils/soil_water_heat.jl +++ b/src/Physics/Soils/soil_water_heat.jl @@ -11,8 +11,8 @@ function Heat.freezethaw!( sfcc = heat.freezecurve swrc = FreezeCurves.swrc(sfcc) # helper function for computing temperature (inverse enthalpy, if necessary) - _get_temperature(::Type{<:Temperature}, i) = state.T[i] - _get_temperature(::Type{<:Enthalpy}, i) = enthalpyinv(soil, heat, state, i) + _get_temperature(::Type{<:TemperatureBased}, i) = state.T[i] + _get_temperature(::Type{<:EnthalpyBased}, i) = enthalpyinv(soil, heat, state, i) L = heat.prop.L @unpack ch_w, ch_i = thermalproperties(soil) @inbounds @fastmath for i in 1:length(state.T) @@ -31,7 +31,7 @@ function Heat.freezethaw!( state.C[i] = C state.∂H∂T[i] = ∂H∂T state.∂θw∂ψ[i] = ∂θw∂ψ - if THeatForm == Temperature + if THeatForm == TemperatureBased state.H[i] = Heat.enthalpy(T, C, L, θw) end end diff --git a/src/Solvers/DiffEq/DiffEq.jl b/src/Solvers/DiffEq/DiffEq.jl index fd1c92af..1da723a3 100644 --- a/src/Solvers/DiffEq/DiffEq.jl +++ b/src/Solvers/DiffEq/DiffEq.jl @@ -45,6 +45,8 @@ function __init__() export Euler, Heun, DP5, Tsit5, ROCK2, ROCK4, SSPRK22, SSPRK33, SSPRK43 # implicit methods export ImplicitEuler, ImplicitMidpoint, Trapezoid, SSPSDIRK2 + # nonlinear solvers + export NLNewton, NLCGLite include("ode_solvers.jl") end end diff --git a/src/Solvers/DiffEq/ode_solvers.jl b/src/Solvers/DiffEq/ode_solvers.jl index 1303615a..6b7a6668 100644 --- a/src/Solvers/DiffEq/ode_solvers.jl +++ b/src/Solvers/DiffEq/ode_solvers.jl @@ -9,5 +9,4 @@ function DiffEqBase.__init(prob::CryoGridProblem, alg::Union{OrdinaryDiffEqAlgor end # custom nonlinear solvers -export NLCGLite include("nlsolve/nlsolve.jl") diff --git a/src/Solvers/Solvers.jl b/src/Solvers/Solvers.jl index 03f6e9aa..bc11f318 100644 --- a/src/Solvers/Solvers.jl +++ b/src/Solvers/Solvers.jl @@ -14,5 +14,6 @@ include("basic_solvers.jl") export LiteImplicit include("LiteImplicit/LiteImplicit.jl") +export NLCGLite include("DiffEq/DiffEq.jl") @reexport using .DiffEq diff --git a/test/Solvers/cglite_implicit_tests.jl b/test/Solvers/cglite_implicit_tests.jl index 8c9550ca..33b5f91c 100644 --- a/test/Solvers/cglite_implicit_tests.jl +++ b/test/Solvers/cglite_implicit_tests.jl @@ -18,7 +18,7 @@ end z_top = 0.0u"m" z_bot = 1000.0u"m" heatop = Heat.EnthalpyImplicit() - # heatop = Heat.EnthalpyForm() + # heatop = Heat.MOLEnthalpy() soil = Ground(MineralOrganic(por=0.0, org=0.0), heat=HeatBalance(heatop)) strat = @Stratigraphy( z_top => Top(PeriodicBC(HeatBalance, CryoGrid.Dirichlet, P, 1.0, 0.0, T₀)), @@ -51,7 +51,7 @@ end z_top = 0.0u"m" z_bot = 1000.0u"m" heatop = Heat.EnthalpyImplicit() - # heatop = Heat.EnthalpyForm(SFCCPreSolver()) + # heatop = Heat.MOLEnthalpy(SFCCPreSolver()) soil = Ground(MineralOrganic(por=0.3, sat=1.0, org=0.0), heat=HeatBalance(heatop)) strat = @Stratigraphy( z_top => Top(ConstantTemperature(1.0u"°C")), diff --git a/test/test_problems.jl b/test/test_problems.jl index 745683a0..88ff4c84 100644 --- a/test/test_problems.jl +++ b/test/test_problems.jl @@ -9,7 +9,7 @@ function test_heat_conduction_freeW_periodic_bc() soilprofile = SoilProfile( 0.0u"m" => MineralOrganic(por=0.50), ); - heatop = Heat.EnthalpyForm(SFCCPreSolver()) + heatop = Heat.MOLEnthalpy(SFCCPreSolver()) initT = initializer(:T, tempprofile) # Define the simulation time span. tspan = (0.0,2*24*3600.0) From e4f437191182dac6f98a1f9a1e2fa6b3b8f2215f Mon Sep 17 00:00:00 2001 From: Brian Groenke Date: Sat, 30 Sep 2023 13:27:29 +0200 Subject: [PATCH 4/8] Fix SWB cumulative ET not being updated --- .../06_heat_freeW_seb_snow_bucketW_samoylov.jl | 9 +++++++++ examples/08_heat_sfcc_richardseq_samoylov.jl | 15 ++++++++------- src/Physics/Heat/heat_water.jl | 10 +++------- src/Physics/Surface/SWB/coupled_sweb.jl | 11 +++++++++++ 4 files changed, 31 insertions(+), 14 deletions(-) diff --git a/examples/06_heat_freeW_seb_snow_bucketW_samoylov.jl b/examples/06_heat_freeW_seb_snow_bucketW_samoylov.jl index 632dfb64..0248ba4e 100644 --- a/examples/06_heat_freeW_seb_snow_bucketW_samoylov.jl +++ b/examples/06_heat_freeW_seb_snow_bucketW_samoylov.jl @@ -72,8 +72,17 @@ Plots.plot(ustrip.(out.T[Z(Near(zs))]), color=cg[LinRange(0.0,1.0,length(zs))]', # Saturation: Plots.plot(ustrip.(out.sat[Z(Near(zs))]), color=cg[LinRange(0.0,1.0,length(zs))]', ylabel="Soil saturation", leg=false, size=(800,500), dpi=150) +# Runoff +Plots.plot(ustrip.(out.top.runoff), color=cg[LinRange(0.0,1.0,length(zs))]', ylabel="Total runoff", leg=false, size=(800,500), dpi=150) + +# Evapotranspiration +Plots.plot(ustrip.(out.top.ET), color=cg[LinRange(0.0,1.0,length(zs))]', ylabel="Total ET", leg=false, size=(800,500), dpi=150) + # Snow depth: Plots.plot(ustrip.(out.snowpack.dsn), color=cg[LinRange(0.0,1.0,length(zs))]', ylabel="Snow depth", leg=false, size=(800,500), dpi=150) # Integrated ground heat flux: Plots.plot(ustrip.(cumsum(out.top.Qg, dims=2)), color=cg[LinRange(0.0,1.0,length(zs))]', ylabel="Integrated ground heat flux", leg=false, size=(800,500), dpi=150) + +# 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) diff --git a/examples/08_heat_sfcc_richardseq_samoylov.jl b/examples/08_heat_sfcc_richardseq_samoylov.jl index 226ff474..3e8e56e3 100644 --- a/examples/08_heat_sfcc_richardseq_samoylov.jl +++ b/examples/08_heat_sfcc_richardseq_samoylov.jl @@ -24,7 +24,7 @@ waterflow = RichardsEq(swrc=swrc); # We use the enthalpy-based heat diffusion with high accuracy Newton-based solver for inverse enthalpy mapping heatop = Heat.MOLEnthalpy(SFCCNewtonSolver()) -upperbc = WaterHeatBC(SurfaceWaterBalance(rainfall=forcings.rainfall), TemperatureGradient(forcings.Tair, NFactor(nf=0.6, nt=0.9))); +upperbc = WaterHeatBC(SurfaceWaterBalance(forcings), TemperatureGradient(forcings.Tair, NFactor(nf=0.6, nt=0.9))); # We will use a simple stratigraphy with three subsurface soil layers. # Note that the @Stratigraphy macro lets us list multiple subsurface layers without wrapping them in a tuple. @@ -39,10 +39,6 @@ grid = CryoGrid.Presets.DefaultGrid_2cm tile = Tile(strat, grid, initT); u0, du0 = initialcondition!(tile, tspan) prob = CryoGridProblem(tile, u0, tspan, saveat=3*3600, savevars=(:T,:θw,:θwi,:kw)); - -# This is currently somewhat slow since the integrator must take very small time steps during the thawed season; -# expect it to take about 3-5 minutes per year on a typical workstation/laptop. -# Minor speed-ups might be possible by tweaking the dt limiters or by using the `SFCCPreSolver` for the freeze curve. integrator = init(prob, CGEuler()) # Here we take just one step to check if it's working. @@ -52,9 +48,12 @@ step!(integrator) # We can use the `getstate` function to construct the current Tile state from the integrator. # We then check that all water fluxes are near zero since we're starting in frozen conditions. state = getstate(integrator); -@assert all(isapprox.(0.0, state.topsoil.jw, atol=1e-14)) +@assert all(isapprox.(0.0, state.ground1.jw, atol=1e-14)) -# Run the integrator forward in time until the end of the tspan: +# Run the integrator forward in time until the end of the tspan; +# Note that this is currently somewhat slow since the integrator must take very small time steps during the thawed season; +# expect it to take about 3-5 minutes per year on a typical workstation/laptop. +# Minor speed-ups might be possible by tweaking the dt limiters or by using the `SFCCPreSolver` for the freeze curve. @time while integrator.t < prob.tspan[end] @assert all(isfinite.(integrator.u)) @assert all(0 .<= integrator.u.sat .<= 1) @@ -65,6 +64,8 @@ state = getstate(integrator); end; out = CryoGridOutput(integrator.sol) +@run step!(integrator) + # Check mass conservation... water_added = values(sum(upreferred.(forcings.rainfall.(tspan[1]:Hour(3):tspan[2]).*u"m/s".*3u"hr")))[1] water_mass = Diagnostics.integrate(out.θwi, tile.grid) diff --git a/src/Physics/Heat/heat_water.jl b/src/Physics/Heat/heat_water.jl index 07c6d28a..72e0a8d9 100644 --- a/src/Physics/Heat/heat_water.jl +++ b/src/Physics/Heat/heat_water.jl @@ -5,10 +5,10 @@ const WaterHeatBC{TWater,THeat} = Coupled2{TWater,THeat} where {TWater<:WaterBC, WaterHeatBC(waterbc::WaterBC, heatbc::HeatBC) = Coupled(waterbc, heatbc) """ - advectiveflux(jw, T, cw, L) + advectiveflux(jw, T₁, T₂, cw, L) -Computes the advective energy flux from a grid cell with temperature `T`, water heat capacity `cw`, -and latent heat of fusion `L`. +Computes the advective energy flux between grid cells with temperatures `T₁` and `T₂` +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) @@ -24,7 +24,6 @@ function water_energy_advection!(sub::SubSurface, ps::Coupled(WaterBalance, Heat let jw = state.jw[i], T₁ = state.T[i-1], T₂ = state.T[i], - # T = (jw > zero(jw))*T₁ + (jw < zero(jw))*T₂, L = heat.prop.L; jH_w = advectiveflux(jw, T₁, T₂, cw, L) state.jH[i] += jH_w @@ -45,7 +44,6 @@ function CryoGrid.interact!(top::Top, bc::WaterBC, sub::SubSurface, heat::HeatBa cw = heatcapacity_water(sub) L = heat.prop.L jw = ssub.jw[1] - # T = (jw > zero(jw))*T_ub + (jw < zero(jw))*Ts jH_w = advectiveflux(jw, T_ub, Ts, cw, L) ssub.jH[1] += jH_w return nothing @@ -56,7 +54,6 @@ function CryoGrid.interact!(sub::SubSurface, heat::HeatBalance, bot::Bottom, bc: cw = heatcapacity_water(sub) L = heat.prop.L jw = ssub.jw[end] - # T = (jw > zero(jw))*Ts + (jw < zero(jw))*T_lb jH_w = advectiveflux(jw, Ts, T_lb, cw, L) ssub.jH[end] += jH_w return nothing @@ -81,7 +78,6 @@ function CryoGrid.interact!( cw1 = heatcapacity_water(sub1) cw2 = heatcapacity_water(sub2) cw = (jw >= zero(jw))*cw1 + (jw < zero(jw))*cw2 - # T = (jw > zero(jw))*T₁ + (jw < zero(jw))*T₂ jH_w = advectiveflux(jw, T₁, T₂, cw, L) state1.jH[end] = state2.jH[1] += jH_w end diff --git a/src/Physics/Surface/SWB/coupled_sweb.jl b/src/Physics/Surface/SWB/coupled_sweb.jl index 90d2f5c0..83b3fcc7 100644 --- a/src/Physics/Surface/SWB/coupled_sweb.jl +++ b/src/Physics/Surface/SWB/coupled_sweb.jl @@ -9,6 +9,17 @@ CryoGrid.variables(top::Top, bc::SurfaceWaterEnergyBalance) = ( Prognostic(:ET, Scalar, u"m^3"), ) +function ET!(::Top, ::SurfaceWaterBalance, state) + @setscalar state.∂ET∂t += state.jw_ET[1]*area(state.grid) +end + +function CryoGrid.computefluxes!(top::Top, bc::SurfaceWaterEnergyBalance, state) + swb, seb = bc + computefluxes!(top, swb, state) + ET!(top, swb, state) + computefluxes!(top, seb, state) +end + function CryoGrid.interact!(top::Top, bc::SurfaceWaterEnergyBalance, sub::SubSurface, ps::Coupled(WaterBalance, HeatBalance), stop, ssub) swb, seb = bc water, heat = ps From 9bead18fa01f3acd74cec7b3da5da1b4d57bcb80 Mon Sep 17 00:00:00 2001 From: Brian Groenke Date: Sat, 30 Sep 2023 15:06:31 +0200 Subject: [PATCH 5/8] Add option to disable heat advection --- examples/08_heat_sfcc_richardseq_samoylov.jl | 15 +++-- src/Physics/Heat/heat_types.jl | 5 +- src/Physics/Heat/heat_water.jl | 58 +++++++++++--------- 3 files changed, 46 insertions(+), 32 deletions(-) diff --git a/examples/08_heat_sfcc_richardseq_samoylov.jl b/examples/08_heat_sfcc_richardseq_samoylov.jl index 3e8e56e3..4d82bf00 100644 --- a/examples/08_heat_sfcc_richardseq_samoylov.jl +++ b/examples/08_heat_sfcc_richardseq_samoylov.jl @@ -19,7 +19,7 @@ initT = initializer(:T, tempprofile); # Here we conigure the water retention curve and freeze curve. The van Genuchten parameters coorespond to that # which would be reasonable for a silty soil. swrc = VanGenuchten(α=0.1, n=1.8) -sfcc = PainterKarra(ω=0.0, swrc=swrc) +sfcc = PainterKarra(ω=0.0; swrc) waterflow = RichardsEq(swrc=swrc); # We use the enthalpy-based heat diffusion with high accuracy Newton-based solver for inverse enthalpy mapping @@ -28,18 +28,23 @@ upperbc = WaterHeatBC(SurfaceWaterBalance(forcings), TemperatureGradient(forcing # 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. +heat = HeatBalance(heatop, freezecurve=sfcc, advection=false) +water = WaterBalance(RichardsEq(; swrc)) strat = @Stratigraphy( -2.0u"m" => Top(upperbc), - 0.0u"m" => Ground(MineralOrganic(por=0.80,sat=0.7,org=0.75), heat=HeatBalance(heatop, freezecurve=sfcc), water=WaterBalance(RichardsEq(;swrc))), - 0.2u"m" => Ground(MineralOrganic(por=0.40,sat=0.8,org=0.10), heat=HeatBalance(heatop, freezecurve=sfcc), water=WaterBalance(RichardsEq(;swrc))), - 2.0u"m" => Ground(MineralOrganic(por=0.10,sat=1.0,org=0.0), heat=HeatBalance(heatop, freezecurve=sfcc), water=WaterBalance(RichardsEq(;swrc))), + 0.0u"m" => Ground(MineralOrganic(por=0.80,sat=0.7,org=0.75); heat, water), + 0.2u"m" => Ground(MineralOrganic(por=0.40,sat=0.8,org=0.10); heat, water), + 2.0u"m" => Ground(MineralOrganic(por=0.10,sat=1.0,org=0.0); heat, water), 1000.0u"m" => Bottom(GeothermalHeatFlux(0.053u"W/m^2")) ); grid = CryoGrid.Presets.DefaultGrid_2cm tile = Tile(strat, grid, initT); u0, du0 = initialcondition!(tile, tspan) prob = CryoGridProblem(tile, u0, tspan, saveat=3*3600, savevars=(:T,:θw,:θwi,:kw)); -integrator = init(prob, CGEuler()) +integrator = init(prob, Euler(), dt=60.0) + +using BenchmarkTools +@btime $tile($du0, $u0, $prob.p, $prob.tspan[1]) # Here we take just one step to check if it's working. step!(integrator) diff --git a/src/Physics/Heat/heat_types.jl b/src/Physics/Heat/heat_types.jl index c8dbfa3b..072693dc 100644 --- a/src/Physics/Heat/heat_types.jl +++ b/src/Physics/Heat/heat_types.jl @@ -34,10 +34,11 @@ Base.@kwdef struct HeatBalance{Tfc<:FreezeCurve,THeatOp<:HeatOperator,Tdt,Tprop} op::THeatOp = MOLEnthalpy(default_fcsolver(freezecurve)) prop::Tprop = HeatBalanceProperties() dtlim::Tdt = default_dtlim(op) # timestep limiter - function HeatBalance(freezecurve, op, prop, dtlim) + advection::Bool = true # whether or not to include advective fluxes when coupled with WaterBalance + function HeatBalance(freezecurve, op, prop, dtlim, advection) # check that heat configuration is valid _validate_heat_config(freezecurve, op) - return new{typeof(freezecurve),typeof(op),typeof(dtlim),typeof(prop)}(freezecurve, op, prop, dtlim) + return new{typeof(freezecurve),typeof(op),typeof(dtlim),typeof(prop)}(freezecurve, op, prop, dtlim, advection) end end # convenience constructors for HeatBalance diff --git a/src/Physics/Heat/heat_water.jl b/src/Physics/Heat/heat_water.jl index 72e0a8d9..b60dc245 100644 --- a/src/Physics/Heat/heat_water.jl +++ b/src/Physics/Heat/heat_water.jl @@ -39,23 +39,27 @@ end # Water/heat interactions function CryoGrid.interact!(top::Top, bc::WaterBC, sub::SubSurface, heat::HeatBalance, stop, ssub) - T_ub = getscalar(stop.T_ub) - Ts = ssub.T[1] - cw = heatcapacity_water(sub) - L = heat.prop.L - jw = ssub.jw[1] - jH_w = advectiveflux(jw, T_ub, Ts, cw, L) - ssub.jH[1] += jH_w + if heat.advection + T_ub = getscalar(stop.T_ub) + Ts = ssub.T[1] + cw = heatcapacity_water(sub) + L = heat.prop.L + jw = ssub.jw[1] + jH_w = advectiveflux(jw, T_ub, Ts, cw, L) + ssub.jH[1] += jH_w + end return nothing end function CryoGrid.interact!(sub::SubSurface, heat::HeatBalance, bot::Bottom, bc::WaterBC, ssub, sbot) - Ts = ssub.T[end] - T_lb = getscalar(sbot.T_ub) - cw = heatcapacity_water(sub) - L = heat.prop.L - jw = ssub.jw[end] - jH_w = advectiveflux(jw, Ts, T_lb, cw, L) - ssub.jH[end] += jH_w + if heat.advection + Ts = ssub.T[end] + T_lb = getscalar(sbot.T_ub) + cw = heatcapacity_water(sub) + L = heat.prop.L + jw = ssub.jw[end] + jH_w = advectiveflux(jw, Ts, T_lb, cw, L) + ssub.jH[end] += jH_w + end return nothing end function CryoGrid.interact!( @@ -70,23 +74,27 @@ function CryoGrid.interact!( interact!(sub1, p1[1], sub2, p2[1], state1, state2) # heat interaction interact!(sub1, p1[2], sub2, p2[2], state1, state2) - # advective energy flux - T₁ = state1.T[end] - T₂ = state2.T[1] - jw = state1.jw[end] - L = p1[2].prop.L - cw1 = heatcapacity_water(sub1) - cw2 = heatcapacity_water(sub2) - cw = (jw >= zero(jw))*cw1 + (jw < zero(jw))*cw2 - jH_w = advectiveflux(jw, T₁, T₂, cw, L) - state1.jH[end] = state2.jH[1] += jH_w + if p1[2].advection && p2[2].advection + # advective energy flux + T₁ = state1.T[end] + T₂ = state2.T[1] + jw = state1.jw[end] + L = p1[2].prop.L + cw1 = heatcapacity_water(sub1) + cw2 = heatcapacity_water(sub2) + cw = (jw >= zero(jw))*cw1 + (jw < zero(jw))*cw2 + jH_w = advectiveflux(jw, T₁, T₂, cw, L) + state1.jH[end] = state2.jH[1] += jH_w + end end # Flux calculation function CryoGrid.computefluxes!(sub::SubSurface, ps::Coupled(WaterBalance, HeatBalance), state) water, heat = ps CryoGrid.computefluxes!(sub, water, state) - water_energy_advection!(sub, ps, state) + if heat.advection + water_energy_advection!(sub, ps, state) + end CryoGrid.computefluxes!(sub, heat, state) end From c4a091fbfb07c155a6ba03b38e8dc2db2106680e Mon Sep 17 00:00:00 2001 From: Brian Groenke Date: Sat, 30 Sep 2023 17:31:29 +0200 Subject: [PATCH 6/8] Rename file with Ground type definitions --- src/Physics/Soils/Soils.jl | 2 +- src/Physics/Soils/{ground_types.jl => ground.jl} | 0 2 files changed, 1 insertion(+), 1 deletion(-) rename src/Physics/Soils/{ground_types.jl => ground.jl} (100%) diff --git a/src/Physics/Soils/Soils.jl b/src/Physics/Soils/Soils.jl index 5e2aa337..3fac8310 100644 --- a/src/Physics/Soils/Soils.jl +++ b/src/Physics/Soils/Soils.jl @@ -27,7 +27,7 @@ const EnthalpyBased = Heat.EnthalpyBased const EnthalpyImplicit = Heat.EnthalpyImplicit export Ground, AbstractGround -include("ground_types.jl") +include("ground.jl") export Soil, SoilParameterization, Heterogeneous include("soil_types.jl") diff --git a/src/Physics/Soils/ground_types.jl b/src/Physics/Soils/ground.jl similarity index 100% rename from src/Physics/Soils/ground_types.jl rename to src/Physics/Soils/ground.jl From 58617b0ceaddd02849ebf6878180b6b096e8e69d Mon Sep 17 00:00:00 2001 From: Brian Groenke Date: Sat, 30 Sep 2023 17:39:01 +0200 Subject: [PATCH 7/8] Ignore richards eq example in doc build --- docs/make.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/docs/make.jl b/docs/make.jl index a3881c14..49861618 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -36,7 +36,8 @@ rm(examples_output_dir, recursive=true, force=true) # recreate directory mkpath(examples_output_dir) # generate example docs from scripts -example_docfiles = map(filter(∉(["Manifest.toml", "Project.toml"]), readdir(examples_dir))) do f +ignored_files = ["Manifest.toml", "Project.toml", "08_heat_sfcc_richardseq_samoylov.jl"] +example_docfiles = map(filter(∉(ignored_files), readdir(examples_dir))) do f infile = joinpath(examples_dir, f) @info "Generating docpage for example script $infile and writing to directory $examples_output_dir" Literate.markdown(infile, examples_output_dir, execute=false, documenter=true) From 8bf74fcff5f54193a1b78863210086e71a7e2570 Mon Sep 17 00:00:00 2001 From: Brian Groenke Date: Sat, 30 Sep 2023 17:39:36 +0200 Subject: [PATCH 8/8] Bump revision --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 43fd625c..9da003c9 100755 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "CryoGrid" uuid = "a535b82e-5f3d-4d97-8b0b-d6483f5bebd5" authors = ["Brian Groenke ", "Jan Nitzbon ", "Moritz Langer "] -version = "0.19.2" +version = "0.19.3" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"