diff --git a/docs/src/manual/overview.md b/docs/src/manual/overview.md
index 02a51354026a33c44dfb647b0a71442b70f814a7..155977f0cac455a4c30fb82a7530678ed69b876a 100644
--- a/docs/src/manual/overview.md
+++ b/docs/src/manual/overview.md
@@ -53,7 +53,7 @@ variables(soil::Soil, heat::Heat{:H}) = (
Prognostic(:H, Float"J/m^3", OnGrid(Cells)),
Diagnostic(:T, Float"°C", OnGrid(Cells)),
Diagnostic(:C, Float"J//K*/m^3", OnGrid(Cells)),
- Diagnostic(:Ceff, Float"J/K/m^3", OnGrid(Cells)),
+ Diagnostic(:dHdT, Float"J/K/m^3", OnGrid(Cells)),
Diagnostic(:k, Float"W/m/K", OnGrid(Edges)),
Diagnostic(:kc, Float"W//m/K", OnGrid(Cells)),
# this last line just appends any state variables or parameters
diff --git a/src/Diagnostics/Diagnostics.jl b/src/Diagnostics/Diagnostics.jl
index f7b5835781259f0a8db29a3c1874ab04b9023f6d..f59f85547fd83f2d56598e9d9f540be9e5874f16 100644
--- a/src/Diagnostics/Diagnostics.jl
+++ b/src/Diagnostics/Diagnostics.jl
@@ -27,8 +27,8 @@ function zero_annual_amplitude(T::AbstractDimArray{<:TempQuantity}; threshold=0.
Tarr = TimeArray(dims(T, Ti).val, permutedims(T, (Ti,Z)))
annual_min = collapse(Tarr, year, first, minimum)
annual_max = collapse(Tarr, year, first, maximum)
- z_inds = argmax(abs.(values(annual_max) .- values(annual_min)) .< threshold, dims=2)
- rebuild(T, dims(T, Z)[[i[2] for i in z_inds]], (Ti(timestamp(annual_max)),))[:,1]
+ z_inds = vec(argmax(abs.(values(annual_max) .- values(annual_min)) .< threshold, dims=2))
+ rebuild(T, collect(dims(T, Z)[[i[2] for i in z_inds]]), (Ti(timestamp(annual_max)),))
end
"""
permafrosttable(T::AbstractDimArray{<:TempQuantity})
@@ -42,8 +42,8 @@ function permafrosttable(T::AbstractDimArray{<:TempQuantity})
Tarr = TimeArray(dims(T, Ti).val, permutedims(T, (Ti,Z)))
annual_max = collapse(Tarr, year, first, maximum)
# find argmax (first one/"true" bit) of annual_max < 0
- z_inds = argmax(values(annual_max) .< 0u"°C", dims=2)
- rebuild(T, dims(T, Z)[[i[2] for i in z_inds]], (Ti(timestamp(annual_max)),))[:,1]
+ z_inds = vec(argmax(values(annual_max) .< 0u"°C", dims=2))
+ rebuild(T, collect(dims(T, Z)[[i[2] for i in z_inds]]), (Ti(timestamp(annual_max)),))
end
"""
permafrostbase(T::AbstractDimArray{<:TempQuantity})
@@ -57,8 +57,8 @@ function permafrostbase(T::AbstractDimArray{<:TempQuantity})
Tarr = TimeArray(dims(T, Ti).val, permutedims(T, (Ti,Z)))
annual_max = collapse(Tarr, year, first, maximum)
# find argmax (first one/"true" bit) of annual_max < 0
- z_inds = size(Tarr,2) .- [i[2] for i in argmax(reverse(values(annual_max) .< 0u"°C", dims=2), dims=2)]
- rebuild(T, dims(T, Z)[[i for i in z_inds]], (Ti(timestamp(annual_max)),))[:,1]
+ z_inds = vec(size(Tarr,2) .- [i[2] for i in argmax(reverse(values(annual_max) .< 0u"°C", dims=2), dims=2)])
+ rebuild(T, collect(dims(T, Z)[[i for i in z_inds]]), (Ti(timestamp(annual_max)),))
end
"""
thawdepth(T::AbstractDimArray{<:TempQuantity})
@@ -69,8 +69,8 @@ and `Z` (depth) in any order.
function thawdepth(T::AbstractDimArray{<:TempQuantity})
_check_arr_dims(T)
T = permutedims(T, (Z,Ti))
- z_inds = argmax((T .<= 0u"°C").data, dims=1)
- rebuild(T, dims(T, Z)[[i[1] for i in z_inds]], (dims(T,Ti),))[1,:]
+ z_inds = vec(argmax((T .<= 0u"°C").data, dims=1))
+ rebuild(T, collect(dims(T, Z)[[i[1] for i in z_inds]]), (dims(T,Ti),))
end
"""
active_layer_thickness(T::AbstractDimArray{<:TempQuantity})
diff --git a/src/Drivers/courant_step.jl b/src/Drivers/courant_step.jl
index 70fd96850293cae5b9135b3fdbdddf28c6b21993..863a943a0d6da910e17f4ed589205e719a30461a 100644
--- a/src/Drivers/courant_step.jl
+++ b/src/Drivers/courant_step.jl
@@ -6,9 +6,9 @@ end
function (cfl::CFLStepLimiter{<:Tile,<:AbstractArray})(u,p,t)
let Δt = cfl.Δt,
Δx = dustrip(Δ(cfl.tile.grid)),
- Ceff = getvar(:Ceff, cfl.tile, u), # apparent heat capacity
+ dHdT = getvar(:dHdT, cfl.tile, u), # apparent heat capacity
kc = getvar(:kc, cfl.tile, u); # thermal cond. at grid centers
- @. Δt = Utils.adstrip(Δx^2 * Ceff / kc)
+ @. Δt = Utils.adstrip(Δx^2 * dHdT / kc)
Δt_min = minimum(Δt)
IfElse.ifelse(isfinite(Δt_min) && Δt_min > 0, Δt_min, cfl.default_dt)
end
@@ -16,9 +16,9 @@ end
function (cfl::CFLStepLimiter{<:Tile,Nothing})(u,p,t)
let Δt = cfl.Δt,
Δx = dustrip(Δ(cfl.tile.grid)),
- Ceff = getvar(:Ceff, cfl.tile, u), # apparent heat capacity
+ dHdT = getvar(:dHdT, cfl.tile, u), # apparent heat capacity
kc = getvar(:kc, cfl.tile, u); # thermal cond. at grid centers
- Δt = Utils.adstrip(Δx^2 * Ceff / kc)
+ Δt = Utils.adstrip(Δx^2 * dHdT / kc)
Δt_min = minimum(Δt)
IfElse.ifelse(isfinite(Δt_min) && Δt_min > 0, Δt_min, cfl.default_dt)
end
diff --git a/src/Numerics/math.jl b/src/Numerics/math.jl
index 8b0d7d262db735e2e32f90d01b46548d87b42a93..f6fd90c8c1109be037d76bc0802250d7943d7007 100644
--- a/src/Numerics/math.jl
+++ b/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
"""
diff --git a/src/Physics/HeatConduction/HeatConduction.jl b/src/Physics/HeatConduction/HeatConduction.jl
index b688de7a7dfba1b9ffed2a1fa01da611c7ad2b31..8d3c16b68b6e47956b2afec09f54fbb737431691 100644
--- a/src/Physics/HeatConduction/HeatConduction.jl
+++ b/src/Physics/HeatConduction/HeatConduction.jl
@@ -18,13 +18,13 @@ using Interpolations: Linear, Flat
using IntervalSets
using ModelParameters
using Parameters
+using SimulationLogs
using Unitful
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 +52,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")
diff --git a/src/Physics/HeatConduction/heat.jl b/src/Physics/HeatConduction/heat.jl
index 9949f11c583b2d45858e1054c791a9cb2c523bfd..94d07ef147d5e12826c2bcc23e4def90f3cb8713 100644
--- a/src/Physics/HeatConduction/heat.jl
+++ b/src/Physics/HeatConduction/heat.jl
@@ -1,3 +1,60 @@
+# 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
@@ -43,24 +100,26 @@ variables(sub::SubSurface, heat::Heat) = (
"""
Variable definitions for heat conduction (enthalpy).
"""
-variables(::Heat{<:FreezeCurve,Enthalpy}) = (
+variables(heat::Heat{<:FreezeCurve,Enthalpy}) = (
Prognostic(:H, Float"J/m^3", OnGrid(Cells)),
Diagnostic(:T, Float"°C", OnGrid(Cells)),
- Diagnostic(:C, Float"J//K/m^3", OnGrid(Cells)),
- Diagnostic(:Ceff, Float"J/K/m^3", OnGrid(Cells)),
- Diagnostic(:k, Float"W/m/K", OnGrid(Edges)),
- Diagnostic(:kc, Float"W//m/K", OnGrid(Cells)),
- Diagnostic(:θl, Float"1", OnGrid(Cells)),
+ _variables(heat)...,
)
"""
Variable definitions for heat conduction (temperature).
"""
-variables(::Heat{<:FreezeCurve,Temperature}) = (
+variables(heat::Heat{<:FreezeCurve,Temperature}) = (
Prognostic(:T, Float"°C", OnGrid(Cells)),
Diagnostic(:H, Float"J/m^3", OnGrid(Cells)),
- Diagnostic(:dH, Float"J/s/m^3", OnGrid(Cells)),
+ Diagnostic(:dH, Float"W/m^3", OnGrid(Cells)),
+ _variables(heat)...,
+)
+"""
+Common variable definitions for all heat implementations.
+"""
+_variables(::Heat) = (
+ Diagnostic(:dHdT, Float"J/K/m^3", OnGrid(Cells)),
Diagnostic(:C, Float"J/K/m^3", OnGrid(Cells)),
- Diagnostic(:Ceff, Float"J/K/m^3", OnGrid(Cells)),
Diagnostic(:k, Float"W/m/K", OnGrid(Edges)),
Diagnostic(:kc, Float"W/m/K", OnGrid(Cells)),
Diagnostic(:θl, Float"1", OnGrid(Cells)),
@@ -97,7 +156,7 @@ function prognosticstep!(::SubSurface, ::Heat{<:FreezeCurve,Temperature}, state)
heatconduction!(state.dH, state.T, ΔT, state.k, Δk)
# Compute temperature flux by dividing by C_eff;
# C_eff should be computed by the freeze curve.
- @inbounds @. state.dT = state.dH / state.Ceff
+ @inbounds @. state.dT = state.dH / state.dHdT
return nothing
end
@inline boundaryflux(::Neumann, bc::BoundaryProcess, top::Top, heat::Heat, sub::SubSurface, stop, ssub) = boundaryvalue(bc,top,heat,sub,stop,ssub)
@@ -130,7 +189,8 @@ function interact!(top::Top, bc::BoundaryProcess, sub::SubSurface, heat::Heat, s
# assumes (1) k has already been computed, (2) surface conductivity = cell conductivity
@inbounds ssub.k[1] = ssub.kc[1]
# boundary flux
- @inbounds ssub.dH[1] += boundaryflux(bc, top, heat, sub, stop, ssub)
+ @log dH_upper = boundaryflux(bc, top, heat, sub, stop, ssub)
+ @inbounds ssub.dH[1] += dH_upper
return nothing # ensure no allocation
end
"""
@@ -141,7 +201,8 @@ function interact!(sub::SubSurface, heat::Heat, bot::Bottom, bc::BoundaryProcess
# assumes (1) k has already been computed, (2) bottom conductivity = cell conductivity
@inbounds ssub.k[end] = ssub.kc[end]
# boundary flux
- @inbounds ssub.dH[end] += boundaryflux(bc, bot, heat, sub, sbot, ssub)
+ @log dH_lower = boundaryflux(bc, bot, heat, sub, sbot, ssub)
+ @inbounds ssub.dH[end] += dH_lower
return nothing # ensure no allocation
end
"""
@@ -192,13 +253,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 dHdT)
+ state.dHdT[i] = liqfrac > 0.0 ? 1e8 : 1/state.C[i]
end
return nothing
end
diff --git a/src/Physics/HeatConduction/heat_bc.jl b/src/Physics/HeatConduction/heat_bc.jl
index e7ea2fdf7135e2c204e09bd0f34872c4b2863eaa..4b793ee31eb64b76d02ba55a3a6c9faa4089bd7a 100644
--- a/src/Physics/HeatConduction/heat_bc.jl
+++ b/src/Physics/HeatConduction/heat_bc.jl
@@ -15,12 +15,12 @@ BoundaryStyle(::Type{<:TemperatureGradient{<:Damping}}) = Neumann()
@inline boundaryvalue(bc::TemperatureGradient, l1, p2, l2, s1, s2) = bc.T(s1.t)
@with_kw struct NFactor{W,S} <: BoundaryEffect
- winterfactor::W = Param(1.0, bounds=(0.0,1.0)) # applied when Tair <= 0
- summerfactor::S = Param(1.0, bounds=(0.0,1.0)) # applied when Tair > 0
+ nf::W = Param(1.0, bounds=(0.0,1.0)) # applied when Tair <= 0
+ nt::S = Param(1.0, bounds=(0.0,1.0)) # applied when Tair > 0
end
@inline function boundaryvalue(bc::TemperatureGradient{<:NFactor}, l1, ::Heat, l2, s1, s2)
- let nfw = bc.effect.winterfactor,
- nfs = bc.effect.summerfactor,
+ let nfw = bc.effect.nf,
+ nfs = bc.effect.nt,
Tair = bc.T(s1.t),
nf = (Tair <= zero(Tair))*nfw + (Tair > zero(Tair))*nfs;
nf*Tair # apply damping factor to air temperature
diff --git a/src/Physics/HeatConduction/soil/sfcc.jl b/src/Physics/HeatConduction/soil/sfcc.jl
index 75fdc97c0341010f117e13e69c970642ff3005b1..1038467ddba8551a4448d3f100017485c1caaa8a 100644
--- a/src/Physics/HeatConduction/soil/sfcc.jl
+++ b/src/Physics/HeatConduction/soil/sfcc.jl
@@ -63,7 +63,7 @@ function (sfcc::SFCC)(soil::Soil, heat::Heat{<:SFCC,Temperature}, state)
state.θl[i] = f(f_argsᵢ...)
state.C[i] = heatcapacity(soil, θw, state.θl[i], θm, θo)
state.H[i] = enthalpy(state.T[i], state.C[i], L, state.θl[i])
- state.Ceff[i] = L*∇f(f_argsᵢ) + state.C[i]
+ state.dHdT[i] = L*∇f(f_argsᵢ) + state.C[i]
end
end
return nothing
@@ -295,7 +295,7 @@ function (s::SFCCNewtonSolver)(soil::Soil, heat::Heat{<:SFCC,Enthalpy}, state, f
let θl = state.θl[i],
H = state.H[i];
state.C[i] = heatcapacity(soil,θw,θl,θm,θo)
- state.Ceff[i] = state.C[i] + dθdT
+ state.dHdT[i] = state.C[i] + dθdT
state.T[i] = (H - L*θl) / state.C[i]
end
end
diff --git a/src/Physics/HeatConduction/soil/soilheat.jl b/src/Physics/HeatConduction/soil/soilheat.jl
index ed11bac9a31cbe84a8a7357269e27a009be3674b..be3d56d5862d7d8723a54b528467fa517c1c1dba 100644
--- a/src/Physics/HeatConduction/soil/soilheat.jl
+++ b/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;
diff --git a/test/Physics/HeatConduction/heat_conduction_tests.jl b/test/Physics/HeatConduction/heat_conduction_tests.jl
index 50cffb0650f03f2dd541e521522b5e54c7a08c7f..91b11c129b81f8e483ea6f88bc9567607873da3f 100644
--- a/test/Physics/HeatConduction/heat_conduction_tests.jl
+++ b/test/Physics/HeatConduction/heat_conduction_tests.jl
@@ -74,7 +74,7 @@ end
@testset "n-factors" begin
ts = DateTime(2010,1,1):Hour(1):DateTime(2010,1,1,4)
forcing = TimeSeriesForcing([1.0,0.5,-0.5,-1.0,0.1], ts, :Tair)
- tgrad = TemperatureGradient(forcing, NFactor(winterfactor=0.5, summerfactor=1.0))
+ tgrad = TemperatureGradient(forcing, NFactor(nf=0.5, nt=1.0))
sub = TestGroundLayer()
heat = Heat()
f1(t) = boundaryvalue(tgrad,Top(),heat,sub,(t=t,),(t=t,))
diff --git a/test/Physics/HeatConduction/sfcc_tests.jl b/test/Physics/HeatConduction/sfcc_tests.jl
index d1c6748a95cd5713e601335c700be10de6f2bcd0..b17b2ff82d6cea4ec56ec3b8540a91ff6e9ce312 100644
--- a/test/Physics/HeatConduction/sfcc_tests.jl
+++ b/test/Physics/HeatConduction/sfcc_tests.jl
@@ -36,7 +36,7 @@ using ComponentArrays
C = heatcapacity.(soil,θw,θl,θm,θo);
enthalpy.(T,C,L,θl) # compute enthalpy at "true" temperature
end
- state = (T=T,C=C,Ceff=similar(C),H=H,θl=θl,)
+ state = (T=T,C=C,dHdT=similar(C),H=H,θl=θl,)
sfcc(soil, heat, state)
@test all(abs.(T.-(H .- L.*θl)./C) .<= tol)
end
@@ -51,7 +51,7 @@ using ComponentArrays
C = heatcapacity.(soil,θw,θl,θm,θo);
enthalpy.(T,C,L,θl) # compute enthalpy at "true" temperature
end
- state = (T=T,C=C,Ceff=similar(C),H=H,θl=θl,)
+ state = (T=T,C=C,dHdT=similar(C),H=H,θl=θl,)
sfcc(soil, heat, state)
@test all(abs.(T.-(H .- L.*θl)./C) .<= tol)
end
@@ -66,7 +66,7 @@ using ComponentArrays
C = heatcapacity.(soil,θw,θl,θm,θo);
enthalpy.(T,C,L,θl) # compute enthalpy at "true" temperature
end
- state = (T=T,C=C,Ceff=similar(C),H=H,θl=θl,)
+ state = (T=T,C=C,dHdT=similar(C),H=H,θl=θl,)
sfcc(soil, heat, state)
@test all(abs.(T.-(H .- L.*θl)./C) .<= tol)
end
@@ -100,7 +100,7 @@ using ComponentArrays
C = heatcapacity.(soil,θw,θl,θm,θo);
enthalpy.(T,C,L,θl) # compute enthalpy at "true" temperature
end
- state = (T=T,C=C,Ceff=similar(C),H=H,θl=θl,)
+ state = (T=T,C=C,dHdT=similar(C),H=H,θl=θl,)
sfcc(soil, heat, state)
@test all(abs.(T.-(H .- L.*θl)./C) .<= tol)
end
@@ -113,7 +113,7 @@ using ComponentArrays
C = heatcapacity.(soil,θw,θl,θm,θo);
enthalpy.(T,C,L,θl) # compute enthalpy at "true" temperature
end
- state = (T=T,C=C,Ceff=similar(C),H=H,θl=θl,)
+ state = (T=T,C=C,dHdT=similar(C),H=H,θl=θl,)
sfcc(soil, heat, state)
@test all(abs.(T.-(H .- L.*θl)./C) .<= tol)
end
@@ -126,7 +126,7 @@ using ComponentArrays
C = heatcapacity.(soil,θw,θl,θm,θo);
enthalpy.(T,C,L,θl) # compute enthalpy at "true" temperature
end
- state = (T=T,C=C,Ceff=similar(C),H=H,θl=θl,)
+ state = (T=T,C=C,dHdT=similar(C),H=H,θl=θl,)
sfcc(soil, heat, state)
@test all(abs.(T.-(H .- L.*θl)./C) .<= tol)
end
@@ -163,7 +163,7 @@ using ComponentArrays
C = heatcapacity.(soil,θw,θl,θm,θo);
enthalpy.(T,C,L,θl) # compute enthalpy at "true" temperature
end
- state = (T=T,C=C,Ceff=similar(C),H=H,θl=θl,)
+ state = (T=T,C=C,dHdT=similar(C),H=H,θl=θl,)
@inferred sfcc(soil, heat, state)
@test all(abs.(T.-(H .- L.*θl)./C) .<= tol)
end
@@ -176,7 +176,7 @@ using ComponentArrays
C = heatcapacity.(soil,θw,θl,θm,θo);
enthalpy.(T,C,L,θl) # compute enthalpy at "true" temperature
end
- state = (T=T,C=C,Ceff=similar(C),H=H,θl=θl,)
+ state = (T=T,C=C,dHdT=similar(C),H=H,θl=θl,)
@inferred sfcc(soil, heat, state)
@test all(abs.(T.-(H .- L.*θl)./C) .<= tol)
end
@@ -189,7 +189,7 @@ using ComponentArrays
C = heatcapacity.(soil,θw,θl,θm,θo);
enthalpy.(T,C,L,θl) # compute enthalpy at "true" temperature
end
- state = (T=T,C=C,Ceff=similar(C),H=H,θl=θl,)
+ state = (T=T,C=C,dHdT=similar(C),H=H,θl=θl,)
@inferred sfcc(soil, heat, state)
@test all(abs.(T.-(H .- L.*θl)./C) .<= tol)
end
@@ -219,7 +219,7 @@ using ComponentArrays
T_ .= T
C = similar(C,eltype(p))
θl = similar(θl,eltype(p))
- state = (T=T_,C=C,Ceff=similar(C),H=H,θl=θl,)
+ state = (T=T_,C=C,dHdT=similar(C),H=H,θl=θl,)
SFCC(McKenzie(γ=Param(p.γ[1])))(soil, heat, state)
state.T[1]
end
@@ -255,7 +255,7 @@ function benchmarksfcc()
C = heatcapacity.(soil,θw,θl,θm,θo);
enthalpy.(T,C,L,θl) # compute enthalpy at "true" temperature
end
- state = (T=T,C=C,Ceff=similar(C),H=H,θl=θl,θw=θw.+zero(T))
+ state = (T=T,C=C,dHdT=similar(C),H=H,θl=θl,θw=θw.+zero(T))
# sfcc.solver(soil, heat, state, sfcc.f, sfcc.∇f)
# @time begin
# state.T .= -0.05