From f2afbb85d0f7d3e941be1f3948c4b78173e29373 Mon Sep 17 00:00:00 2001
From: Brian Groenke <brian.groenke@awi.de>
Date: Wed, 15 Dec 2021 16:24:47 +0100
Subject: [PATCH] Fix freeW init bug and make porosity a state var
---
src/Layers/soil.jl | 2 ++
src/Physics/HeatConduction/heat.jl | 18 ++++++------------
src/Physics/HeatConduction/soil/sfcc.jl | 6 +++---
src/Physics/HeatConduction/soil/soilheat.jl | 18 +++++++++++++++---
test/Physics/HeatConduction/sfcc_tests.jl | 20 ++++++++++----------
5 files changed, 36 insertions(+), 28 deletions(-)
diff --git a/src/Layers/soil.jl b/src/Layers/soil.jl
index cd00f5c5..57e597dd 100644
--- a/src/Layers/soil.jl
+++ b/src/Layers/soil.jl
@@ -67,11 +67,13 @@ end
# Methods
porosity(soil::Soil{T,<:SoilCharacteristicFractions}) where T = soilcomp(Val{:θp}(), soil.para)
variables(::Soil) = (
+ Diagnostic(:θp, Float"1", OnGrid(Cells)), # porosity
Diagnostic(:θw, Float"1", OnGrid(Cells)), # total water content (xice + saturated pore space)
Diagnostic(:θm, Float"1", OnGrid(Cells)), # mineral content
Diagnostic(:θo, Float"1", OnGrid(Cells)), # organic content
)
function initialcondition!(soil::Soil{T,<:SoilCharacteristicFractions}, state) where T
+ state.θp .= soilcomp(Val{:θp}(), soil.para)
state.θw .= soilcomp(Val{:θw}(), soil.para)
state.θm .= soilcomp(Val{:θm}(), soil.para)
state.θo .= soilcomp(Val{:θo}(), soil.para)
diff --git a/src/Physics/HeatConduction/heat.jl b/src/Physics/HeatConduction/heat.jl
index e913691e..6ad85006 100644
--- a/src/Physics/HeatConduction/heat.jl
+++ b/src/Physics/HeatConduction/heat.jl
@@ -75,12 +75,6 @@ variables(heat::Heat{<:FreezeCurve,Temperature}) = (
Diagnostic(:kc, Float"W/m/K", OnGrid(Cells)),
Diagnostic(:θl, Float"1", OnGrid(Cells)),
)
-""" Initial condition for heat conduction (all state configurations) on any subsurface layer. """
-function initialcondition!(layer::SubSurface, heat::Heat, state)
- L = heat.L
- heatcapacity!(layer, heat, state)
- @. state.H = enthalpy(state.T, state.C, L, state.θl)
-end
""" Diagonstic step for heat conduction (all state configurations) on any subsurface layer. """
function diagnosticstep!(layer::SubSurface, heat::Heat, state)
# Reset energy flux to zero; this is redundant when H is the prognostic variable
@@ -88,10 +82,10 @@ function diagnosticstep!(layer::SubSurface, heat::Heat, state)
@. state.dH = zero(eltype(state.dH))
# Evaluate the freeze curve (updates T, C, and θl)
fc! = freezecurve(heat);
- fc!(layer,heat,state)
+ fc!(layer, heat, state)
# Update thermal conductivity
thermalconductivity!(layer, heat, state)
- # Harmonic mean of conductivities
+ # Harmonic mean of inner conductivities
@inbounds let k = (@view state.k[2:end-1]),
Δk = Δ(state.grids.k);
harmonicmean!(k, state.kc, Δk)
@@ -144,7 +138,7 @@ Generic top interaction. Computes flux dH at top cell.
function interact!(top::Top, bc::BoundaryProcess, sub::SubSurface, heat::Heat, stop, ssub)
# thermal conductivity at boundary
# assumes (1) k has already been computed, (2) surface conductivity = cell conductivity
- @inbounds ssub.k[1] = ssub.k[2]
+ @inbounds ssub.k[1] = ssub.kc[1]
# boundary flux
@inbounds ssub.dH[1] += boundaryflux(bc, top, heat, sub, stop, ssub)
return nothing # ensure no allocation
@@ -155,7 +149,7 @@ Generic bottom interaction. Computes flux dH at bottom cell.
function interact!(sub::SubSurface, heat::Heat, bot::Bottom, bc::BoundaryProcess, ssub, sbot)
# thermal conductivity at boundary
# assumes (1) k has already been computed, (2) bottom conductivity = cell conductivity
- @inbounds ssub.k[end] = ssub.k[end-1]
+ @inbounds ssub.k[end] = ssub.kc[end]
# boundary flux
@inbounds ssub.dH[end] += boundaryflux(bc, bot, heat, sub, sbot, ssub)
return nothing # ensure no allocation
@@ -192,7 +186,7 @@ total water content (θw), and liquid water content (θl).
"""
@inline function (fc::FreeWater)(layer::SubSurface, heat::Heat{FreeWater,Enthalpy}, state)
@inline function enthalpyinv(H, C, L, θtot)
- let θtot = max(1.0e-8,θtot),
+ let θtot = max(1.0e-8, θtot),
Lθ = L*θtot,
I_t = H > Lθ,
I_f = H <= 0.0;
@@ -200,7 +194,7 @@ total water content (θw), and liquid water content (θl).
end
end
@inline function freezethaw(H, L, θtot)
- let θtot = max(1.0e-8,θtot),
+ let θtot = max(1.0e-8, θtot),
Lθ = L*θtot,
I_t = H > Lθ,
I_c = (H > 0.0) && (H <= Lθ);
diff --git a/src/Physics/HeatConduction/soil/sfcc.jl b/src/Physics/HeatConduction/soil/sfcc.jl
index 23d4e01f..530905bf 100644
--- a/src/Physics/HeatConduction/soil/sfcc.jl
+++ b/src/Physics/HeatConduction/soil/sfcc.jl
@@ -93,7 +93,7 @@ end
sfccparams(f::DallAmico, soil::Soil, heat::Heat, state) = (
f.Tₘ,
f.θres,
- porosity(soil), # θ saturated = porosity
+ state.θp, # θ saturated = porosity
state.θw, # total water content
heat.L, # specific latent heat of fusion, L
f.α,
@@ -129,7 +129,7 @@ end
sfccparams(f::McKenzie, soil::Soil, heat::Heat, state) = (
f.Tₘ,
f.θres,
- porosity(soil), # θ saturated = porosity
+ state.θp, # θ saturated = porosity
state.θw, # total water content
f.γ,
)
@@ -154,7 +154,7 @@ end
sfccparams(f::Westermann, soil::Soil, heat::Heat, state) = (
f.Tₘ,
f.θres,
- porosity(soil), # θ saturated = porosity
+ state.θp, # θ saturated = porosity
state.θw, # total water content
f.δ,
)
diff --git a/src/Physics/HeatConduction/soil/soilheat.jl b/src/Physics/HeatConduction/soil/soilheat.jl
index 667ebc25..c35eb9b3 100644
--- a/src/Physics/HeatConduction/soil/soilheat.jl
+++ b/src/Physics/HeatConduction/soil/soilheat.jl
@@ -26,11 +26,23 @@ end
include("sfcc.jl")
""" Initial condition for heat conduction (all state configurations) on soil layer. """
-function initialcondition!(soil::Soil, heat::Heat{<:SFCC}, state)
+function initialcondition!(soil::Soil, heat::Heat, state)
initialcondition!(soil, state)
L = heat.L
- sfcc = freezecurve(heat)
- state.θl .= sfcc.f.(state.T, sfccparams(sfcc.f, soil, heat, state)...)
+ fc! = freezecurve(heat)
+ fc!(soil, heat, state)
heatcapacity!(soil, heat, state)
@. state.H = enthalpy(state.T, state.C, L, state.θl)
end
+""" Initial condition for heat conduction (all state configurations) with free water freeze curve on soil layer. """
+function initialcondition!(soil::Soil, heat::Heat{FreeWater}, state)
+ initialcondition!(soil, state)
+ L = heat.L
+ # initialize liquid water content based on temperature
+ @. state.θl = ifelse.(state.T > 0.0, 0.0, state.θp)
+ heatcapacity!(soil, heat, state)
+ @. state.H = enthalpy(state.T, state.C, L, state.θl)
+ # evaluate freeze curve
+ fc! = freezecurve(heat)
+ fc!(soil, heat, state)
+end
diff --git a/test/Physics/HeatConduction/sfcc_tests.jl b/test/Physics/HeatConduction/sfcc_tests.jl
index 76ae35a1..2916e4bc 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,θw=θw.+zero(T),θm=θm.+zero(T),θo=θo.+zero(T))
+ state = (T=T,C=C,Ceff=similar(C),H=H,θl=θl,θw=θw.+zero(T),θm=θm.+zero(T),θo=θo.+zero(T),θp=θp.+zero(T))
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,θw=θw.+zero(T),θm=θm.+zero(T),θo=θo.+zero(T))
+ state = (T=T,C=C,Ceff=similar(C),H=H,θl=θl,θw=θw.+zero(T),θm=θm.+zero(T),θo=θo.+zero(T),θp=θp.+zero(T))
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,θw=θw.+zero(T),θm=θm.+zero(T),θo=θo.+zero(T))
+ state = (T=T,C=C,Ceff=similar(C),H=H,θl=θl,θw=θw.+zero(T),θm=θm.+zero(T),θo=θo.+zero(T),θp=θp.+zero(T))
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,θw=θw.+zero(T),θm=θm.+zero(T),θo=θo.+zero(T))
+ state = (T=T,C=C,Ceff=similar(C),H=H,θl=θl,θw=θw.+zero(T),θm=θm.+zero(T),θo=θo.+zero(T),θp=θp.+zero(T))
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,θw=θw.+zero(T),θm=θm.+zero(T),θo=θo.+zero(T))
+ state = (T=T,C=C,Ceff=similar(C),H=H,θl=θl,θw=θw.+zero(T),θm=θm.+zero(T),θo=θo.+zero(T),θp=θp.+zero(T))
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,θw=θw.+zero(T),θm=θm.+zero(T),θo=θo.+zero(T))
+ state = (T=T,C=C,Ceff=similar(C),H=H,θl=θl,θw=θw.+zero(T),θm=θm.+zero(T),θo=θo.+zero(T),θp=θp.+zero(T))
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,θw=θw.+zero(T),θm=θm.+zero(T),θo=θo.+zero(T))
+ state = (T=T,C=C,Ceff=similar(C),H=H,θl=θl,θw=θw.+zero(T),θm=θm.+zero(T),θo=θo.+zero(T),θp=θp.+zero(T))
@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,θw=θw.+zero(T),θm=θm.+zero(T),θo=θo.+zero(T))
+ state = (T=T,C=C,Ceff=similar(C),H=H,θl=θl,θw=θw.+zero(T),θm=θm.+zero(T),θo=θo.+zero(T),θp=θp.+zero(T))
@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,θw=θw.+zero(T),θm=θm.+zero(T),θo=θo.+zero(T))
+ state = (T=T,C=C,Ceff=similar(C),H=H,θl=θl,θw=θw.+zero(T),θm=θm.+zero(T),θo=θo.+zero(T),θp=θp.+zero(T))
@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,θw=θw.+zero(T),θm=θm.+zero(T),θo=θo.+zero(T))
+ state = (T=T_,C=C,Ceff=similar(C),H=H,θl=θl,θw=θw.+zero(T),θm=θm.+zero(T),θo=θo.+zero(T),θp=θp.+zero(T))
SFCC(McKenzie(γ=Param(p.γ[1])))(soil, heat, state)
state.T[1]
end
--
GitLab