Minor refactoring/cleanup in heat conduction code

src/Numerics/math.jl

@@ -27,7 +27,7 @@ end
 
 
 @propagate_inbounds function _nonlineardiffusion!(∂y, x₁, x₂, x₃, k₁, k₂, Δx₁, Δx₂, Δk)
-    @. ∂y = (k₂*(x₃ - x₂)/Δx₂ - k₁*(x₂-x₁)/Δx₁)/Δk
+    @. ∂y = (k₂*(x₃ - x₂)/Δx₂ - k₁*(x₂ - x₁)/Δx₁)/Δk
 end
 @propagate_inbounds function _nonlineardiffusion!(
     ∂y::AbstractVector{Float64},
@@ -40,7 +40,7 @@ end
     Δx₂::AbstractVector{Float64},
     Δk::AbstractVector{Float64},
 )
-    @turbo @. ∂y = (k₂*(x₃ - x₂)/Δx₂ - k₁*(x₂-x₁)/Δx₁)/Δk
+    @turbo @. ∂y = (k₂*(x₃ - x₂)/Δx₂ - k₁*(x₂ - x₁)/Δx₁)/Δk
 end
 
 """

src/Physics/HeatConduction/HeatConduction.jl

@@ -24,7 +24,6 @@ import Flatten: @flattenable, flattenable
 
 export Heat, TemperatureProfile
 export FreeWater, FreezeCurve, freezecurve
-export enthalpy, heatcapacity, heatcapacity!, thermalconductivity, thermalconductivity!
 
 abstract type FreezeCurve end
 struct FreeWater <: FreezeCurve end
@@ -52,56 +51,8 @@ freezecurve(heat::Heat) = heat.freezecurve
 variables(::SubSurface, ::Heat, ::FreezeCurve) = ()
 # Fallback (error) implementation for freeze curve
 (fc::FreezeCurve)(sub::SubSurface, heat::Heat, state) = error("freeze curve $(typeof(fc)) not implemented for $(typeof(heat)) on layer $(typeof(sub))")
-# Thermal properties (generic)
-"""
-    enthalpy(T, C, L, θ) = T*C + L*θ
 
-Discrete enthalpy function on temperature, heat capacity, specific latent heat of fusion, and liquid water content.
-"""
-@inline enthalpy(T, C, L, θ) = T*C + L*θ
-"""
-    totalwater(sub::SubSurface, heat::Heat, state)
-    totalwater(sub::SubSurface, heat::Heat, state, i)
-
-Retrieves the total water content for the given layer at grid cell `i`, if specified.
-Defaults to using the scalar total water content defined on layer `sub`.
-"""
-@inline totalwater(sub::SubSurface, heat::Heat, state) = totalwater(sub)
-@inline totalwater(sub::SubSurface, heat::Heat, state, i) = Utils.getscalar(totalwater(sub, heat, state), i)
-"""
-    heatcapacity(sub::SubSurface, heat::Heat, state, i)
-
-Computes the heat capacity at grid cell `i` for the given layer from the current state.
-"""
-heatcapacity(sub::SubSurface, heat::Heat, state, i) = error("heatcapacity not defined for $(typeof(heat)) on $(typeof(sub))")
-"""
-    heatcapacity!(sub::SubSurface, heat::Heat, state)
-
-Computes the heat capacity for the given layer from the current state and stores the result in-place in the state variable `C`.
-"""
-@inline function heatcapacity!(sub::SubSurface, heat::Heat, state)
-    @inbounds for i in 1:length(state.T)
-        state.C[i] = heatcapacity(sub, heat, state, i)
-    end
-end
-"""
-    thermalconductivity(sub::SubSurface, heat::Heat, state, i)
-
-Computes the thrmal conductivity at grid cell `i` for the given layer from the current state.
-"""
-thermalconductivity(sub::SubSurface, heat::Heat, state, i) = error("thermalconductivity not defined for $(typeof(heat)) on $(typeof(sub))")
-"""
-    thermalconductivity!(sub::SubSurface, heat::Heat, state)
-
-Computes the thermal conductivity for the given layer from the current state and stores the result in-place in the state variable `C`.
-"""
-@inline function thermalconductivity!(sub::SubSurface, heat::Heat, state)
-    @inbounds for i in 1:length(state.T)
-        state.kc[i] = thermalconductivity(sub, heat, state, i)
-    end
-end
-
-export heatconduction!
+export heatconduction!, enthalpy, heatcapacity, heatcapacity!, thermalconductivity, thermalconductivity!
 include("heat.jl")
 export ConstantTemp, GeothermalHeatFlux, TemperatureGradient, NFactor, Damping
 include("heat_bc.jl")

src/Physics/HeatConduction/heat.jl

+# Thermal properties (generic)
+"""
+    enthalpy(T, C, L, θ) = T*C + L*θ
+
+Discrete enthalpy function on temperature, heat capacity, specific latent heat of fusion, and liquid water content.
+"""
+@inline enthalpy(T, C, L, θ) = T*C + L*θ
+"""
+    totalwater(sub::SubSurface, heat::Heat, state)
+    totalwater(sub::SubSurface, heat::Heat, state, i)
+
+Retrieves the total water content for the given layer at grid cell `i`, if specified.
+Defaults to using the scalar total water content defined on layer `sub`.
+"""
+@inline totalwater(sub::SubSurface, heat::Heat, state) = totalwater(sub)
+@inline totalwater(sub::SubSurface, heat::Heat, state, i) = Utils.getscalar(totalwater(sub, heat, state), i)
+"""
+    liquidwater(sub::SubSurface, heat::Heat, state)
+    liquidwater(sub::SubSurface, heat::Heat, state, i)
+
+Retrieves the liquid water content for the given layer at grid cell `i`, if specified.
+Defaults to using the scalar total water content defined on layer `sub`.
+"""
+@inline liquidwater(sub::SubSurface, heat::Heat, state) = state.θl
+@inline liquidwater(sub::SubSurface, heat::Heat, state, i) = Utils.getscalar(liquidwater(sub, heat, state), i)
+"""
+    heatcapacity(sub::SubSurface, heat::Heat, state, i)
+
+Computes the heat capacity at grid cell `i` for the given layer from the current state.
+"""
+heatcapacity(sub::SubSurface, heat::Heat, state, i) = error("heatcapacity not defined for $(typeof(heat)) on $(typeof(sub))")
+"""
+    heatcapacity!(sub::SubSurface, heat::Heat, state)
+
+Computes the heat capacity for the given layer from the current state and stores the result in-place in the state variable `C`.
+"""
+@inline function heatcapacity!(sub::SubSurface, heat::Heat, state)
+    @inbounds for i in 1:length(state.T)
+        state.C[i] = heatcapacity(sub, heat, state, i)
+    end
+end
+"""
+    thermalconductivity(sub::SubSurface, heat::Heat, state, i)
+
+Computes the thrmal conductivity at grid cell `i` for the given layer from the current state.
+"""
+thermalconductivity(sub::SubSurface, heat::Heat, state, i) = error("thermalconductivity not defined for $(typeof(heat)) on $(typeof(sub))")
+"""
+    thermalconductivity!(sub::SubSurface, heat::Heat, state)
+
+Computes the thermal conductivity for the given layer from the current state and stores the result in-place in the state variable `C`.
+"""
+@inline function thermalconductivity!(sub::SubSurface, heat::Heat, state)
+    @inbounds for i in 1:length(state.T)
+        state.kc[i] = thermalconductivity(sub, heat, state, i)
+    end
+end
 """
     heatconduction!(∂H,T,ΔT,k,Δk)
 
@@ -11,8 +68,8 @@ function heatconduction!(∂H,T,ΔT,k,Δk)
         T₁=T[1],
         k=k[2],
         δ=ΔT[1],
-        a=Δk[1];
-        k*(T₂-T₁)/δ/a
+        d=Δk[1];
+        k*(T₂-T₁)/δ/d
     end
     # diffusion on non-boundary cells
     @inbounds let T = T,
@@ -26,8 +83,8 @@ function heatconduction!(∂H,T,ΔT,k,Δk)
         T₁=T[end-1],
         k=k[end-1],
         δ=ΔT[end],
-        a=Δk[end];
-        -k*(T₂-T₁)/δ/a
+        d=Δk[end];
+        -k*(T₂-T₁)/δ/d
     end
     return nothing
 end
@@ -192,13 +249,18 @@ total water content (θw), and liquid water content (θl).
         end
     end
     @inbounds for i in 1:length(state.H)
-        let θw = totalwater(sub, heat, state, i),
-            H = state.H[i],
-            L = heat.L;
-            state.θl[i] = θw*freezethaw(H, L, θw)
-            state.C[i] = heatcapacity(sub, heat, state, i) # update heat capacity, C
-            state.T[i] = enthalpyinv(H, state.C[i], L, θw)
-        end
+        θw = totalwater(sub, heat, state, i)
+        H = state.H[i]
+        L = heat.L
+        # liquid water content = (total water content) * (liquid fraction)
+        liqfrac = freezethaw(H, L, θw)
+        state.θl[i] = θw*liqfrac
+        # update heat capacity
+        state.C[i] = heatcapacity(sub, heat, state, i)
+        # enthalpy inverse function
+        state.T[i] = enthalpyinv(H, state.C[i], L, θw)
+        # set dHdT (a.k.a Ceff)
+        state.Ceff[i] = liqfrac > 0.0 ? 1e8 : 1/state.C[i]
     end
     return nothing
 end

src/Physics/HeatConduction/soil/soilheat.jl

@@ -34,12 +34,12 @@ Defaults to using the scalar porosity defined on `soil`.
 @inline heatcapacity(soil::Soil, heat::Heat, state, i) = heatcapacity(
     soil,
     totalwater(soil, heat, state, i),
-    state.θl[i],
+    liquidwater(soil, heat, state, i),
     mineral(soil, heat, state, i),
     organic(soil, heat, state, i),
 )
 @inline function heatcapacity(soil::Soil, totalwater, liquidwater, mineral, organic)
-    @unpack cw,co,cm,ca,ci = soil.hc
+    @unpack cw, co, cm, ca, ci = soil.hc
     let air = 1.0 - totalwater - mineral - organic,
         ice = totalwater - liquidwater,
         liq = liquidwater;
@@ -49,12 +49,12 @@ end
 @inline thermalconductivity(soil::Soil, heat::Heat, state, i) = thermalconductivity(
     soil,
     totalwater(soil, heat, state, i),
-    state.θl[i],
+    liquidwater(soil, heat, state, i),
     mineral(soil, heat, state, i),
     organic(soil, heat, state, i),
 )
 @inline function thermalconductivity(soil::Soil, totalwater, liquidwater, mineral, organic)
-    @unpack kw,ko,km,ka,ki = soil.tc
+    @unpack kw, ko, km, ka, ki = soil.tc
     let air = 1.0 - totalwater - mineral - organic,
         ice = totalwater - liquidwater,
         liq = liquidwater;