CompatHelper: bump compat for DiffEqCallbacks to 3, (keep existing compat)

Project.toml

@@ -49,7 +49,7 @@ ConstructionBase = "1.4"
 DataStructures = "0.18"
 Dates = "1"
 DiffEqBase = "6"
-DiffEqCallbacks = "2"
+DiffEqCallbacks = "2, 3"
 DimensionalData = "0.19,0.20,0.21,0.24, 0.25"
 Downloads = "1"
 Flatten = "0.4"

examples/heat_freeW_lite_implicit.jl

@@ -11,11 +11,11 @@ using CryoGrid.LiteImplicit
 forcings = loadforcings(CryoGrid.Presets.Forcings.Samoylov_ERA_MkL3_CCSM4_long_term);
 forcings = Base.rename(forcings, :Ptot => :precip)
 soilprofile = SoilProfile(
-    0.0u"m" => MineralOrganic(por=0.80,sat=1.0,org=0.75),
-    0.1u"m" => MineralOrganic(por=0.80,sat=1.0,org=0.25),
-    0.4u"m" => MineralOrganic(por=0.55,sat=1.0,org=0.25),
-    3.0u"m" => MineralOrganic(por=0.50,sat=1.0,org=0.0),
-    10.0u"m" => MineralOrganic(por=0.30,sat=1.0,org=0.0),
+    0.0u"m" => MineralOrganic(por=0.80, org=0.75),
+    0.1u"m" => MineralOrganic(por=0.80, org=0.25),
+    0.4u"m" => MineralOrganic(por=0.55, org=0.25),
+    3.0u"m" => MineralOrganic(por=0.50, org=0.0),
+    10.0u"m" => MineralOrganic(por=0.30, org=0.0),
 )
 z_top = -2.0u"m"
 z_bot = 1000.0u"m"

src/Numerics/math.jl

@@ -154,6 +154,13 @@ function harmonicmean!(h::AbstractVector, x::AbstractVector, w::AbstractVector)
     end
 end
 
+"""
+    linearmean(x, w)
+
+Calculates a weighted average of `x` with weights `w`.
+"""
+linearmean(x, w) = sum(map(*, x, w))
+
 """
     tdma_solve!(x, a, b, c, d)
 

src/Physics/Heat/heat_conduction.jl

@@ -12,7 +12,7 @@ function freezethaw!(sub::SubSurface, heat::HeatBalance{FreeWater,<:EnthalpyBase
         H = state.H[i]
         θwi = Hydrology.watercontent(sub, state, i)
         state.θw[i], _ = FreezeCurves.freewater(H, θwi, L)
-        C = state.C[i] = heatcapacity(sub, heat, state, i)
+        C = state.C[i] = heatcapacity(sub, state, i)
         T = state.T[i] = enthalpyinv(sub, heat, H, θwi, C, L)
         # set ∂H∂T (a.k.a ∂H∂T)
         state.∂H∂T[i] = iszero(T) ? 1e8 : C

src/Physics/Heat/heat_methods.jl

@@ -9,18 +9,52 @@ thermalproperties(sub::SubSurface, state) = thermalproperties(sub)
 thermalproperties(sub::SubSurface, state, i) = thermalproperties(sub, state)
 
 """
-    thermalconductivity(::SubSurface, state, i)
+    thermalconductivities(sub::SubSurface)
 
-Computes the thermal conductivity for the given `SubSurface` layer at grid cell `i`.
+Get thermal conductivities for generic `SubSurface` layer.
 """
-thermalconductivity(::SubSurface, state, i) = error("not implemented")
+function thermalconductivities(sub::SubSurface)
+    @unpack kh_w, kh_i, kh_a = thermalproperties(sub)
+    return (kh_w, kh_i, kh_a)
+end
+thermalconductivities(sub::SubSurface, state) = thermalconductivities(sub)
+thermalconductivities(sub::SubSurface, state, i) = thermalconductivities(sub, state)
 
 """
-    heatcapacity(::SubSurface, state, i)
+    heatcapacities(::SubSurface)
 
-Computes the heat capacity for the given `SubSurface` layer at grid cell `i`.
+Get heat capacities for generic `SubSurface` layer.
 """
-heatcapacity(::SubSurface, state, i) = error("not implemented")
+function heatcapacities(sub::SubSurface)
+    @unpack ch_w, ch_i, ch_a = thermalproperties(sub)
+    return ch_w, ch_i, ch_a
+end
+heatcapacities(sub::SubSurface, state) = heatcapacities(sub)
+heatcapacities(sub::SubSurface, state, i) = heatcapacities(sub, state)
+
+"""
+    thermalconductivity(sub::SubSurface, heat::HeatBalance, state, i)
+
+Computes the thermal conductivity as a squared weighted sum over constituent conductivities with volumetric fractions `θfracs`.
+"""
+function thermalconductivity(sub::SubSurface, state, i)
+    θs = volumetricfractions(sub, state, i)
+    ks = thermalconductivities(sub, state, i)
+    f = thermalproperties(sub).thermcond
+    return f(ks, θs)
+end
+
+"""
+    heatcapacity(sub::SubSurface, state, i)
+
+Computes the heat capacity as a weighted average over constituent capacities with volumetric fractions `θfracs`.
+"""
+function heatcapacity(sub::SubSurface, state, i)
+    θs = volumetricfractions(sub, state, i)
+    cs = heatcapacities(sub, state, i)
+    f = thermalproperties(sub).heatcap
+    return f(cs, θs)
+end
 
 """
     freezethaw!(sub::SubSurface, heat::HeatBalance, state)
@@ -61,17 +95,3 @@ 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))
-
-"""
-    thermal_conductivity_function(op::HeatOperator)
-
-Retreives the thermal conductivity function for this heat operator.
-"""
-thermal_conductivity_function(op::HeatOperator) = op.cond
-
-"""
-    heat_capacity_function(op::HeatOperator)
-
-Retreives the volumetric heat capacity function for this heat operator.
-"""
-heat_capacity_function(op::HeatOperator) = op.hc

src/Physics/Heat/heat_types.jl

@@ -49,30 +49,24 @@ _validate_heat_config(::FreeWater, ::TemperatureBased) = error("Invalid heat bal
 
 # Heat operators
 """
-    Diffusion1D{progvar,Tcond,Thc} <: HeatOperator{progvar}
+    Diffusion1D{progvar} <: HeatOperator{progvar}
 
 Represents a standard method-of-lines (MOL) forward diffusion operator in 1 dimension.
 """
-struct Diffusion1D{progvar,Tcond,Thc} <: HeatOperator{progvar}
-    cond::Tcond
-    hc::Thc
-    Diffusion1D(progvar::Symbol=:H, cond=quadratic_parallel_conductivity, hc=weighted_average_heatcapacity) = new{progvar,typeof(cond),typeof(hc)}(cond, hc)
+struct Diffusion1D{progvar} <: HeatOperator{progvar}
+    Diffusion1D(progvar::Symbol=:H) = new{progvar}()
 end
-ConstructionBase.constructorof(::Type{<:Diffusion1D{progvar}}) where {progvar} = (cond, hc) -> Diffusion1D(progvar, cond, hc)
+ConstructionBase.constructorof(::Type{<:Diffusion1D{progvar}}) where {progvar} = () -> Diffusion1D(progvar)
 
 """
-    EnthalpyImplicit{Tcond,Thc} <: HeatOperator{:H}
+    EnthalpyImplicit <: HeatOperator{:H}
 
 Implicit enthalpy formulation of Swaminathan and Voller (1992) and Langer et al. (2022). Note that this
 heat operator formulation does not compute a divergence `dH` but only computes the necessary diffusion
 coefficients for use by an appropriate solver. See the `LiteImplicit` module for the appropriate
 solver algorithms.
 """
-struct EnthalpyImplicit{Tcond,Thc} <: HeatOperator{:H}
-    cond::Tcond
-    hc::Thc
-    EnthalpyImplicit(cond=quadratic_parallel_conductivity, hc=weighted_average_heatcapacity) = new{typeof(cond),typeof(hc)}(cond, hc)
-end
+struct EnthalpyImplicit <: HeatOperator{:H} end
 
 """
 Numerical constants for pararameterizing heat processes.

src/Physics/Heat/thermal_properties.jl

@@ -8,6 +8,8 @@ this includes the thermal properties of water, ice, and air. This can be extende
 by passing additional properties into the constructor.
 """
 Utils.@properties ThermalProperties(
+    thermcond = quadratic_parallel_conductivity,
+    heatcap = Numerics.linearmean,
     kh_w = 0.57u"J/s/m/K", # thermal conductivity of water [Hillel (1982)]
     kh_i = 2.2u"J/s/m/K", # thermal conductivity of ice [Hillel (1982)]
     kh_a = 0.025u"J/s/m/K", # thermal conductivity of air [Hillel (1982)]
@@ -18,18 +20,6 @@ Utils.@properties ThermalProperties(
 
 # Thermal conductivity
 
-"""
-    thermalconductivities(sub::SubSurface)
-
-Get thermal conductivities for generic `SubSurface` layer.
-"""
-function thermalconductivities(sub::SubSurface)
-    @unpack kh_w, kh_i, kh_a = thermalproperties(sub)
-    return (kh_w, kh_i, kh_a)
-end
-thermalconductivities(sub::SubSurface, state) = thermalconductivities(sub)
-thermalconductivities(sub::SubSurface, state, i) = thermalconductivities(sub, state)
-
 """
     quadratic_parallel_conductivity(ks, θs)
 
@@ -55,26 +45,14 @@ Woodside, W. & Messmer, J.H. 1961. Thermal conductivity of porous media. I. Unco
 """
 geometric_conductivity(ks, θs) = prod(map((k,θ) -> k^θ, ks, θs))
 
-"""
-    thermalconductivity(sub::SubSurface, heat::HeatBalance, state, i)
-
-Computes the thermal conductivity as a squared weighted sum over constituent conductivities with volumetric fractions `θfracs`.
-"""
-function thermalconductivity(sub::SubSurface, heat::HeatBalance, state, i)
-    θs = volumetricfractions(sub, state, i)
-    ks = thermalconductivities(sub, state, i)
-    f = thermal_conductivity_function(heat.op)
-    return f(ks, θs)
-end
-
 """
     thermalconductivity!(sub::SubSurface, heat::HeatBalance, state)
 
 Computes the thermal conductivity for the given layer from the current state and stores the result in-place in the state variable `k`.
 """
-function thermalconductivity!(sub::SubSurface, heat::HeatBalance, state)
+function thermalconductivity!(sub::SubSurface, ::HeatBalance, state)
     @inbounds for i in 1:length(state.T)
-        state.kc[i] = thermalconductivity(sub, heat, state, i)
+        state.kc[i] = thermalconductivity(sub, state, i)
         if i > 1
             Δk₁ = CryoGrid.thickness(sub, state, i-1)
             Δk₂ = CryoGrid.thickness(sub, state, i)
@@ -90,44 +68,12 @@ end
 # Heat capacity
 
 """
-    heatcapacities(::SubSurface)
-
-Get heat capacities for generic `SubSurface` layer.
-"""
-function heatcapacities(sub::SubSurface)
-    @unpack ch_w, ch_i, ch_a = thermalproperties(sub)
-    return ch_w, ch_i, ch_a
-end
-heatcapacities(sub::SubSurface, state) = heatcapacities(sub)
-heatcapacities(sub::SubSurface, state, i) = heatcapacities(sub, state)
-
-"""
-    weighted_average_heatcapacity(cs, θs)
-
-Represents a simple composite heat capacity that is the sum of each constituent heat capacity weighted
-by the volumetric fraction.
-"""
-weighted_average_heatcapacity(cs, θs) = sum(map(*, cs, θs))
-
-"""
-    heatcapacity(sub::SubSurface, heat::HeatBalance, state, i)
-
-Computes the heat capacity as a weighted average over constituent capacities with volumetric fractions `θfracs`.
-"""
-function heatcapacity(sub::SubSurface, heat::HeatBalance, state, i)
-    θs = volumetricfractions(sub, state, i)
-    cs = heatcapacities(sub, state, i)
-    f = heat_capacity_function(heat.op)
-    return f(cs, θs)
-end
-
-"""
-    heatcapacity!(sub::SubSurface, heat::HeatBalance, state)
+    heatcapacity!(sub::SubSurface, state)
 
 Computes the heat capacity for the given layer from the current state and stores the result in-place in the state variable `C`.
 """
-function heatcapacity!(sub::SubSurface, heat::HeatBalance, state)
+function heatcapacity!(sub::SubSurface, state)
     @inbounds for i in 1:length(state.T)
-        state.C[i] = heatcapacity(sub, heat, state, i)
+        state.C[i] = heatcapacity(sub, state, i)
     end
 end

src/Physics/Salt/salt_diffusion.jl

@@ -20,7 +20,7 @@ function Heat.freezethaw!(
             state.θw[i] = θw = ForwardDiff.value(res_dual.θw)
             state.∂θw∂T[i] = ∂θw∂T = ForwardDiff.partials(res_dual.θw)[1]
             state.∂θw∂c[i] = ForwardDiff.partials(res_dual.θw)[2]
-            state.C[i] = C = heatcapacity(soil, heat, state, i)
+            state.C[i] = C = heatcapacity(soil, state, i)
             state.∂H∂T[i] =  C + L*∂θw∂T
             state.H[i] = enthalpy(T, C, L, θw)
             state.dₛ_mid[i] = salt.prop.dₛ₀ * θw / salt.prop.τ

src/Physics/Snow/snow_heat.jl

@@ -27,7 +27,7 @@ function Heat.thermalconductivity(sq::SturmQuadratic, ρsn::Number)
 end
 
 # extract thermal conductivity scheme from thermal properties struct and invoke special dispatches defined above.
-Heat.thermalconductivity(snow::Snowpack, heat::HeatBalance, state, i) = thermalconductivity(snow.para.heat.cond, state.ρsn[i])
+Heat.thermalconductivity(snow::Snowpack, state, i) = thermalconductivity(snow.para.heat.cond, state.ρsn[i])
 
 function Heat.enthalpyinv(::Snowpack, heat::HeatBalance{FreeWater,<:EnthalpyBased}, H, θwi, C, L)
     T_f = H / C

src/Physics/Soils/Soils.jl

@@ -38,11 +38,9 @@ include("soil_methods.jl")
 export SoilTexture
 include("soil_texture.jl")
 
-export Heterogeneous, MineralOrganic, soilcomponent
-include("soil_para.jl")
-
-export SURFEX
-include("soil_surfex.jl")
+export Heterogeneous, MineralOrganic, SURFEX
+export soilcomponent
+include("para/soil_para.jl")
 
 export RichardsEq
 include("soil_water.jl")

src/Physics/Soils/ground.jl

@@ -1,21 +1,21 @@
 abstract type GroundParameterization end
 
 """
-    AbstractGround{Tpara<:GroundParameterization,Theat<:Optional{HeatBalance},Twater<:Optional{WaterBalance}} <: SubSurface
+    AbstractGround{Tpara<:GroundParameterization,Theat<:HeatBalance,Twater<:WaterBalance} <: SubSurface
 
 Base type for all ground layers defining heat and water balance schemes.
 """
 abstract type AbstractGround{Tpara<:GroundParameterization,Theat<:Optional{HeatBalance},Twater<:Optional{WaterBalance}} <: SubSurface end
 
 """
-    Ground{Tpara,Theat<:Optional{HeatBalance},Twater<:Optional{WaterBalance},Taux} <: Soil{Tpara,Theat,Twater}
+    Ground{Tpara,Theat<:HeatBalance,Twater<:WaterBalance,Taux} <: Soil{Tpara,Theat,Twater}
 
 Generic representation of a `Ground` layer with material parameterization `para`.
 """
 Base.@kwdef struct Ground{Tpara,Theat<:Optional{HeatBalance},Twater<:Optional{WaterBalance},Tsolver,Taux} <: AbstractGround{Tpara,Theat,Twater}
     para::Tpara = MineralOrganic() # ground parameterization
     heat::Theat = HeatBalance() # heat conduction
-    water::Twater = nothing # water balance
+    water::Twater = WaterBalance() # water balance
     fcsolver::Tsolver = default_fcsolver(para, heat, water)
     aux::Taux = nothing # user-defined specialization
 end
@@ -27,3 +27,9 @@ default_fcsolver(::GroundParameterization, ::HeatBalance{<:SFCC}, ::Nothing) = S
 default_fcsolver(::GroundParameterization, ::HeatBalance{<:SFCC}, ::WaterBalance) = SFCCPreSolver(FreezeCurves.SFCCPreSolverCacheND())
 
 fcsolver(ground::Ground) = ground.fcsolver
+
+# CryoGrid methods
+
+CryoGrid.processes(g::Ground{<:GroundParameterization,<:HeatBalance,Nothing}) = g.heat
+CryoGrid.processes(g::Ground{<:GroundParameterization,Nothing,<:WaterBalance}) = g.water
+CryoGrid.processes(g::Ground{<:GroundParameterization,<:HeatBalance,<:WaterBalance}) = Coupled(g.water, g.heat)

src/Physics/Soils/soil_para.jl → src/Physics/Soils/para/mineral_organic.jl

@@ -1,5 +1,5 @@
 """
-    MineralOrganic{Tpor,Tsat,Torg} <: SoilParameterization
+    MineralOrganic{Tpor,Tsat,Torg,Thp,Twp} <: SoilParameterization
 
 Represents a simple organic/mineral soil mixutre in terms of its characteristic fractions:
 i.e. natural porosity, saturation, and organic solid fraction. This is the standard CryoGrid representation
@@ -40,7 +40,8 @@ SoilThermalProperties(
     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
-) = ThermalProperties(; kh_w, kh_i, kh_a, kh_m, kh_o, ch_w, ch_i, ch_a, ch_m, ch_o)
+    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::MineralOrganic) = MineralOrganic(
@@ -58,6 +59,8 @@ CryoGrid.initializers(soil::Soil{<:MineralOrganic,THeat,<:WaterBalance}) where {
     initializers(soil, processes(soil))...,
 )
 
+# Heat
+
 function Heat.thermalconductivities(soil::Soil{<:MineralOrganic})
     @unpack kh_w, kh_i, kh_a, kh_m, kh_o = thermalproperties(soil)
     return kh_w, kh_i, kh_a, kh_m, kh_o
@@ -73,6 +76,8 @@ Gets the `ThermalProperties` for the given soil layer.
 """
 Heat.thermalproperties(soil::Soil{<:MineralOrganic}) = soil.para.heat
 
+# Hydrology
+
 """
 Gets the `HydraulicProperties` for the given soil layer.
 """
@@ -88,26 +93,7 @@ Hydrology.maxwater(soil::Soil{<:MineralOrganic}, water::WaterBalance) = porosity
 Hydrology.watercontent(soil::Soil{<:MineralOrganic,THeat,Nothing}, state) where {THeat} = soilcomponent(Val{:θwi}(), soil.para)
 Hydrology.watercontent(soil::Soil{<:MineralOrganic,THeat,<:WaterBalance}, state) where {THeat} = state.θwi
 
-"""
-    Heterogeneous{V,N,D,Taux} <: SoilParameterization
-
-Special `SoilParameterization` which wraps a `Profile` of another soil parameterization type
-to indicate that it should be heterogeneous with over depth.
-"""
-Base.@kwdef struct Heterogeneous{V,N,IT,VT,Taux} <: SoilParameterization
-    profile::SoilProfile{N,IT,VT}
-    aux::Taux = nothing
-    Heterogeneous(::SoilProfile) = error("SoilProfile for heterogeneous layer must have uniform parameterization types (but the parameters may vary).")
-    Heterogeneous(profile::SoilProfile{N,IT,VT}, aux=nothing) where {N,V<:SoilParameterization,IT<:NTuple{N},VT<:NTuple{N,V}} = new{V,N,IT,VT,typeof(aux)}(profile, aux)
-end
-
-"""
-    Ground(soilprofile::SoilProfile; kwargs...)
-"""
-Ground(soilprofile::SoilProfile; kwargs...) = Ground(Heterogeneous(soilprofile); kwargs...)
-
-# add dispatch for default_fcsolver that selects the ND presolver
-default_fcsolver(::Heterogeneous, ::HeatBalance{<:SFCC}, ::Nothing) = SFCCPreSolver(FreezeCurves.SFCCPreSolverCacheND())
+# Heterogeneous
 
 saturation(::Soil{<:Heterogeneous{<:MineralOrganic}}, state) = state.sat
 porosity(::Soil{<:Heterogeneous{<:MineralOrganic}}, state) = state.por