diff --git a/src/Layers/Layers.jl b/src/Layers/Layers.jl
index 1080bc4b341aecd679c334f86c3dcfa76f9dd237..227bfb1b263a7ca38e74baedf51cb47eb3e37b21 100644
--- a/src/Layers/Layers.jl
+++ b/src/Layers/Layers.jl
@@ -12,7 +12,7 @@ using ModelParameters
using Parameters
using Unitful
-export Soil, SoilParameterization, SoilCharacteristicFractions, SoilProfile, SoilType, Sand, Silt, Clay, soilparameters
+export Soil, SoilParameterization, SoilCharacteristicFractions, SoilProfile, SoilType, Sand, Silt, Clay, soilparameters, soilcomp, porosity
include("soil.jl")
end
\ No newline at end of file
diff --git a/src/Layers/soil.jl b/src/Layers/soil.jl
index 8e8c804fc11010ef836a1e5f752fee6284366516..cd00f5c5d9be97276fd4fa36e794a31891e3a7ef 100644
--- a/src/Layers/soil.jl
+++ b/src/Layers/soil.jl
@@ -29,15 +29,11 @@ function soilparameters(;xic=0.0, por=0.5, sat=1.0, org=0.5)
end
SoilProfile(pairs::Pair{<:DistQuantity,<:SoilParameterization}...) = Profile(pairs...)
# Helper functions for obtaining soil compositions from characteristic fractions.
+soilcomp(::Val{var}, fracs::SoilCharacteristicFractions) where var = soilcomp(Val{var}(), fracs.xic, fracs.por, fracs.sat, fracs.org)
soilcomp(::Val{:θp}, χ, ϕ, θ, ω) = (1-χ)*ϕ
soilcomp(::Val{:θw}, χ, ϕ, θ, ω) = χ + (1-χ)*ϕ*θ
soilcomp(::Val{:θm}, χ, ϕ, θ, ω) = (1-χ)*(1-ϕ)*(1-ω)
soilcomp(::Val{:θo}, χ, ϕ, θ, ω) = (1-χ)*(1-ϕ)*ω
-θx(fracs::SoilCharacteristicFractions) = fracs.xic
-θp(fracs::SoilCharacteristicFractions) = soilcomp(Val{:θp}(), fracs.xic, fracs.por, fracs.sat, fracs.org)
-θw(fracs::SoilCharacteristicFractions) = soilcomp(Val{:θw}(), fracs.xic, fracs.por, fracs.sat, fracs.org)
-θm(fracs::SoilCharacteristicFractions) = soilcomp(Val{:θm}(), fracs.xic, fracs.por, fracs.sat, fracs.org)
-θo(fracs::SoilCharacteristicFractions) = soilcomp(Val{:θo}(), fracs.xic, fracs.por, fracs.sat, fracs.org)
"""
Thermal conductivity constants.
"""
@@ -68,3 +64,15 @@ Basic Soil layer.
hc::SoilHCParams = SoilHCParams()
sp::S = nothing # user-defined specialization
end
+# Methods
+porosity(soil::Soil{T,<:SoilCharacteristicFractions}) where T = soilcomp(Val{:θp}(), soil.para)
+variables(::Soil) = (
+ 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.θw .= soilcomp(Val{:θw}(), soil.para)
+ state.θm .= soilcomp(Val{:θm}(), soil.para)
+ state.θo .= soilcomp(Val{:θo}(), soil.para)
+end
diff --git a/src/Physics/HeatConduction/HeatConduction.jl b/src/Physics/HeatConduction/HeatConduction.jl
index e10bb42f1d908524ffc534557f1c2d8d1ccf3ae4..1db0bf5e24d796ebc35de3b881c01411a2fe747b 100644
--- a/src/Physics/HeatConduction/HeatConduction.jl
+++ b/src/Physics/HeatConduction/HeatConduction.jl
@@ -8,7 +8,7 @@ using CryoGrid.Physics.Boundaries
using CryoGrid.Physics.Water: VanGenuchten
using CryoGrid.Numerics
using CryoGrid.Numerics: nonlineardiffusion!, harmonicmean!, harmonicmean, heaviside
-using CryoGrid.Layers: Soil, θp, θw, θm, θo
+using CryoGrid.Layers: Soil, porosity
using CryoGrid.Utils
using DimensionalData
@@ -39,7 +39,7 @@ struct Temperature <: HeatImpl end
Lsl::Float"J/kg" = 334000.0xu"J/kg" #[J/kg] (latent heat of fusion)
L::Float"J/m^3" = Ï*Lsl #[J/m^3] (specific latent heat of fusion)
freezecurve::F = FreeWater() # freeze curve, defautls to free water fc
- sp::S = Enthalpy()
+ sp::S = Enthalpy() # specialization
end
# convenience constructors for specifying prognostic variable as symbol
Heat(var::Union{Symbol,Tuple{Vararg{Symbol}}}; kwargs...) = Heat(Val{var}(); kwargs...)
diff --git a/src/Physics/HeatConduction/heat.jl b/src/Physics/HeatConduction/heat.jl
index 442c7d7ad6c31e8ce1d1acb380f0969fd0cfede2..e913691e8b799614e273833cc5b0730039556aa4 100644
--- a/src/Physics/HeatConduction/heat.jl
+++ b/src/Physics/HeatConduction/heat.jl
@@ -1,14 +1,3 @@
-# Default implementation of `variables` for freeze curve
-variables(::SubSurface, ::Heat, ::FreezeCurve) = ()
-# Fallback (error) implementation for freeze curve
-(fc::FreezeCurve)(layer::SubSurface, heat::Heat, state) =
- error("freeze curve $(typeof(fc)) not implemented for $(typeof(heat)) on layer $(typeof(layer))")
-
-freezecurve(heat::Heat) = heat.freezecurve
-enthalpy(T::Number"°C", C::Number"J/K/m^3", L::Number"J/m^3", θ::Real) = T*C + L*θ
-totalwater(::SubSurface, ::Heat, state) = state.θw
-heatcapacity!(layer::SubSurface, heat::Heat, state) = error("heatcapacity not defined for $(typeof(heat)) on $(typeof(layer))")
-thermalconductivity!(layer::SubSurface, heat::Heat, state) = error("thermalconductivity not defined for $(typeof(heat)) on $(typeof(layer))")
"""
heatconduction!(∂H,T,ΔT,k,Δk)
@@ -42,8 +31,29 @@ function heatconduction!(∂H,T,ΔT,k,Δk)
end
return nothing
end
-""" Variable definitions for heat conduction (enthalpy) on any subsurface layer. """
-variables(layer::SubSurface, heat::Heat{<:FreezeCurve,Enthalpy}) = (
+
+# Default implementation of `variables` for freeze curve
+variables(::SubSurface, ::Heat, ::FreezeCurve) = ()
+# Fallback (error) implementation for freeze curve
+(fc::FreezeCurve)(layer::SubSurface, heat::Heat, state) =
+ error("freeze curve $(typeof(fc)) not implemented for $(typeof(heat)) on layer $(typeof(layer))")
+freezecurve(heat::Heat) = heat.freezecurve
+enthalpy(T, C, L, θ) = T*C + L*θ
+heatcapacity!(layer::SubSurface, heat::Heat, state) = error("heatcapacity not defined for $(typeof(heat)) on $(typeof(layer))")
+thermalconductivity!(layer::SubSurface, heat::Heat, state) = error("thermalconductivity not defined for $(typeof(heat)) on $(typeof(layer))")
+"""
+Variable definitions for heat conduction on any subsurface layer. Joins variables defined on
+`layer` and `heat` individually as well as variables defined by the freeze curve.
+"""
+variables(layer::SubSurface, heat::Heat) = (
+ variables(layer)..., # layer variables
+ variables(heat)..., # heat variables
+ variables(layer, heat, freezecurve(heat))..., # freeze curve variables
+)
+"""
+Variable definitions for heat conduction (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)),
@@ -51,11 +61,11 @@ variables(layer::SubSurface, heat::Heat{<:FreezeCurve,Enthalpy}) = (
Diagnostic(:k, Float"W/m/K", OnGrid(Edges)),
Diagnostic(:kc, Float"W//m/K", OnGrid(Cells)),
Diagnostic(:θl, Float"1", OnGrid(Cells)),
- # add freeze curve variables (if any are present)
- variables(layer, heat, freezecurve(heat))...,
)
-""" Variable definitions for heat conduction (temperature) on any subsurface layer. """
-variables(layer::SubSurface, heat::Heat{<:FreezeCurve,Temperature}) = (
+"""
+Variable definitions for heat conduction (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)),
@@ -64,8 +74,6 @@ variables(layer::SubSurface, heat::Heat{<:FreezeCurve,Temperature}) = (
Diagnostic(:k, Float"W/m/K", OnGrid(Edges)),
Diagnostic(:kc, Float"W/m/K", OnGrid(Cells)),
Diagnostic(:θl, Float"1", OnGrid(Cells)),
- # add freeze curve variables (if any are present)
- variables(layer, heat, freezecurve(heat))...,
)
""" Initial condition for heat conduction (all state configurations) on any subsurface layer. """
function initialcondition!(layer::SubSurface, heat::Heat, state)
@@ -199,11 +207,8 @@ total water content (θw), and liquid water content (θl).
I_c*(H/Lθ) + I_t
end
end
- let L = heat.L,
- θw = totalwater(layer, heat, state);
- @. state.θl = freezethaw(state.H, L, θw)*θw
- heatcapacity!(layer, heat, state) # update heat capacity, C
- @. state.T = enthalpyinv(state.H, state.C, L, θw)
- end
+ @inbounds @. state.θl = freezethaw(state.H, heat.L, state.θw)*state.θw
+ heatcapacity!(layer, heat, state) # update heat capacity, C
+ @inbounds @. state.T = enthalpyinv(state.H, state.C, heat.L, state.θw)
return nothing
end
diff --git a/src/Physics/HeatConduction/soil/sfcc.jl b/src/Physics/HeatConduction/soil/sfcc.jl
index 65acbfa836e580618905342cd6673e3c5347ec63..df04723e6b92e4f8383c4469d0751275f3c6aa42 100644
--- a/src/Physics/HeatConduction/soil/sfcc.jl
+++ b/src/Physics/HeatConduction/soil/sfcc.jl
@@ -49,47 +49,18 @@ Updates state variables according to the specified SFCC function and solver.
For heat conduction with enthalpy, this is implemented as a simple passthrough to the non-linear solver.
For heat conduction with temperature, we can simply evaluate the freeze curve to get C_eff, θl, and H.
"""
-function (sfcc::SFCC)(soil::Soil, heat::Heat{<:SFCC,Enthalpy}, state)
- sfcc.solver(soil, heat, state, sfcc.f, sfcc.∇f)
-end
-function (sfcc::SFCC)(soil::Soil, heat::Heat{<:SFCC,PartitionedEnthalpy}, state)
- let L = heat.L,
- N = length(state.grids.H),
+(sfcc::SFCC)(soil::Soil, heat::Heat{<:SFCC,Enthalpy}, state) = sfcc.solver(soil, heat, state, sfcc.f, sfcc.∇f)
+function (sfcc::SFCC)(soil::Soil, heat::Heat{<:SFCC,Temperature}, state)
+ @inbounds @fastmath let L = heat.L,
f = sfcc.f,
∇f = sfcc.∇f,
f_args = tuplejoin((state.T,),sfccparams(f,soil,heat,state));
- let θw = θw(soil.para),
- θm = θm(soil.para),
- θo = θo(soil.para);
- @inbounds @fastmath for i in 1:N
- state.H[i] = state.Hâ‚›[i] + state.Hâ‚—[i]
- # It is possible for the integrator to violate physical constraints by integrating
- # Hₗ outside of [0,Lθtot]. Here we clamp the result to physically correct values.
- state.θl[i] = clamp(state.Hₗ[i] / L, 0.0, θw)
- state.C[i] = heatcapacity(soil, θw, state.θl[i], θm, θo)
- state.T[i] = state.Hâ‚›[i] / state.C[i]
- f_argsáµ¢ = Utils.selectat(i, identity, f_args)
- state.dθdT[i] = ∇f(f_argsᵢ)
- state.Ceff[i] = L*state.dΘdT[i] + state.C[i]
- end
- end
- end
-end
-function (sfcc::SFCC)(soil::Soil, heat::Heat{<:SFCC,Temperature}, state)
- let θw = θw(soil.para),
- θm = θm(soil.para),
- θo = θo(soil.para);
- @inbounds @fastmath let L = heat.L,
- f = sfcc.f,
- ∇f = sfcc.∇f,
- f_args = tuplejoin((state.T,),sfccparams(f,soil,heat,state));
- for i in 1:length(state.T)
- f_argsáµ¢ = Utils.selectat(i, identity, f_args)
- 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]
- end
+ for i in 1:length(state.T)
+ f_argsáµ¢ = Utils.selectat(i, identity, f_args)
+ state.θl[i] = f(f_argsᵢ...)
+ state.C[i] = heatcapacity(soil, state.θw[i], state.θl[i], state.θm[i], state.θo[i])
+ state.H[i] = enthalpy(state.T[i], state.C[i], L, state.θl[i])
+ state.Ceff[i] = L*∇f(f_argsᵢ) + state.C[i]
end
end
return nothing
@@ -102,7 +73,7 @@ Retrieves a tuple of values corresponding to each parameter declared by SFCCFunc
Soil layer, Heat process, and model state. The order of parameters *must match* the argument order
of the freeze curve function `f`.
"""
-sfccparams(f::SFCCFunction, soil::Soil, heat::Heat, state) = ()
+sfccparams(::SFCCFunction, ::Soil, ::Heat, state) = ()
# Fallback implementation of variables for SFCCFunction
variables(::Soil, ::Heat, f::SFCCFunction) = ()
variables(::Soil, ::Heat, s::SFCCSolver) = ()
@@ -122,8 +93,8 @@ end
sfccparams(f::DallAmico, soil::Soil, heat::Heat, state) = (
f.Tₘ,
f.θres,
- θp(soil.para), # θ saturated = porosity
- θw(soil.para),
+ porosity(soil), # θ saturated = porosity
+ state.θw, # total water content
heat.L, # specific latent heat of fusion, L
f.α,
f.n,
@@ -158,8 +129,8 @@ end
sfccparams(f::McKenzie, soil::Soil, heat::Heat, state) = (
f.Tₘ,
f.θres,
- θp(soil.para), # θ saturated = porosity
- θw(soil.para),
+ porosity(soil), # θ saturated = porosity
+ state.θw, # total water content
f.γ,
)
function (f::McKenzie)(T,Tₘ,θres,θsat,θtot,γ)
@@ -183,8 +154,8 @@ end
sfccparams(f::Westermann, soil::Soil, heat::Heat, state) = (
f.Tₘ,
f.θres,
- θp(soil.para), # θ saturated = porosity
- θw(soil.para),
+ porosity(soil), # θ saturated = porosity
+ state.θw, # total water content
f.δ,
)
function (f::Westermann)(T,Tₘ,θres,θsat,θtot,δ)
@@ -226,21 +197,21 @@ function (s::SFCCNewtonSolver)(soil::Soil, heat::Heat{<:SFCC,Enthalpy}, state, f
# probably shouldn't incur too much of a performance hit, just a few extra stack pointers!
f_args = sfccparams(f, soil, heat, state)
# iterate over each cell and solve the conservation law: H = TC + Lθ
- @inbounds @fastmath for i in 1:length(state.T)
- itercount = 0
- let Tâ‚€ = state.T[i] |> Utils.adstrip,
+ @threaded for i in 1:length(state.T)
+ @inbounds @fastmath let Tâ‚€ = state.T[i] |> Utils.adstrip,
T = Tâ‚€, # temperature
H = state.H[i] |> Utils.adstrip, # enthalpy
C = state.C[i] |> Utils.adstrip, # heat capacity
θl = state.θl[i] |> Utils.adstrip, # liquid water content
- θw = θw(soil.para), # total water content
- θm = θm(soil.para), # mineral content
- θo = θo(soil.para), # organic content
+ θw = state.θw[i] |> Utils.adstrip, # total water content
+ θm = state.θm[i] |> Utils.adstrip, # mineral content
+ θo = state.θo[i] |> Utils.adstrip, # organic content
L = heat.L, # specific latent heat of fusion
cw = soil.hc.cw, # heat capacity of liquid water
α₀ = s.α₀,
Ï„ = s.Ï„,
- f_argsáµ¢ = Utils.selectat(i, Utils.adstrip, f_args);
+ f_argsáµ¢ = Utils.selectat(i, Utils.adstrip, f_args),
+ itercount = 0;
# compute initial guess T by setting θl according to free water scheme
T = let Lθ = L*θw;
if H < 0
diff --git a/src/Physics/HeatConduction/soil/soilheat.jl b/src/Physics/HeatConduction/soil/soilheat.jl
index 1ebafe563fc3402c0b10b2cf3a5ef270cb196012..667ebc25b6e6f5dba422d35b34415ba742c06b78 100644
--- a/src/Physics/HeatConduction/soil/soilheat.jl
+++ b/src/Physics/HeatConduction/soil/soilheat.jl
@@ -1,5 +1,4 @@
-totalwater(soil::Soil, heat::Heat, state) = θw(soil.para)
-function heatcapacity(soil::Soil, totalwater, liquidwater, mineral, organic)
+@inline function heatcapacity(soil::Soil, totalwater, liquidwater, mineral, organic)
@unpack cw,co,cm,ca,ci = soil.hc
let air = 1.0 - totalwater - mineral - organic,
ice = totalwater - liquidwater,
@@ -7,7 +6,7 @@ function heatcapacity(soil::Soil, totalwater, liquidwater, mineral, organic)
liq*cw + ice*ci + mineral*cm + organic*co + air*ca
end
end
-function thermalconductivity(soil::Soil, totalwater, liquidwater, mineral, organic)
+@inline function thermalconductivity(soil::Soil, totalwater, liquidwater, mineral, organic)
@unpack kw,ko,km,ka,ki = soil.tc
let air = 1.0 - totalwater - mineral - organic,
ice = totalwater - liquidwater,
@@ -16,26 +15,19 @@ function thermalconductivity(soil::Soil, totalwater, liquidwater, mineral, organ
end
end
""" Heat capacity for soil layer """
-function heatcapacity!(soil::Soil, ::Heat, state)
- let w = θw(soil.para),
- m = θm(soil.para),
- o = θo(soil.para);
- @. state.C = heatcapacity(soil, w, state.θl, m, o)
- end
+@inline function heatcapacity!(soil::Soil, ::Heat, state)
+ @inbounds @. state.C = heatcapacity(soil, state.θw, state.θl, state.θm, state.θo)
end
""" Thermal conductivity for soil layer """
-function thermalconductivity!(soil::Soil, ::Heat, state)
- let w = θw(soil.para),
- m = θm(soil.para),
- o = θo(soil.para);
- @. state.kc = thermalconductivity(soil, w, state.θl, m, o)
- end
+@inline function thermalconductivity!(soil::Soil, ::Heat, state)
+ @inbounds @. state.kc = thermalconductivity(soil, state.θw, state.θl, state.θm, state.θo)
end
include("sfcc.jl")
""" Initial condition for heat conduction (all state configurations) on soil layer. """
function initialcondition!(soil::Soil, heat::Heat{<:SFCC}, state)
+ initialcondition!(soil, state)
L = heat.L
sfcc = freezecurve(heat)
state.θl .= sfcc.f.(state.T, sfccparams(sfcc.f, soil, heat, state)...)
diff --git a/src/Physics/coupled.jl b/src/Physics/coupled.jl
index 6a135b613de46a86c395b06a10a3f4fb6e21cbbd..cec35b5bfe58bdbbc5203828e4fd532ae10f3e27 100644
--- a/src/Physics/coupled.jl
+++ b/src/Physics/coupled.jl
@@ -3,6 +3,7 @@ diagnosticstep!(::Top, ::CoupledProcesses, state) = nothing
prognosticstep!(::Top, ::CoupledProcesses, state) = nothing
diagnosticstep!(::Bottom, ::CoupledProcesses, state) = nothing
prognosticstep!(::Bottom, ::CoupledProcesses, state) = nothing
+variables(ps::CoupledProcesses) = tuplejoin((variables(p) for p in ps.processes)...)
variables(layer::Layer, ps::CoupledProcesses) = tuplejoin((variables(layer,p) for p in ps.processes)...)
callbacks(layer::Layer, ps::CoupledProcesses) = tuplejoin((callbacks(layer,p) for p in ps.processes)...)
"""
@@ -54,10 +55,19 @@ Default implementation of `prognosticstep!` for multi-process types. Calls each
return expr
end
"""
- initialcondition!(l::Layer, ps::CoupledProcesses{P}, state) where {P}
+ initialcondition!([l::Layer,] ps::CoupledProcesses{P}, state) where {P}
Default implementation of `initialcondition!` for multi-process types. Calls each process in sequence.
"""
+@generated function initialcondition!(ps::CoupledProcesses{P}, state) where {P}
+ expr = Expr(:block)
+ for i in 1:length(P.parameters)
+ @>> quote
+ initialcondition!(ps[$i],state)
+ end push!(expr.args)
+ end
+ return expr
+end
@generated function initialcondition!(l::Layer, ps::CoupledProcesses{P}, state) where {P}
expr = Expr(:block)
for i in 1:length(P.parameters)
@@ -68,7 +78,7 @@ Default implementation of `initialcondition!` for multi-process types. Calls eac
return expr
end
"""
- initialcondition!(l::Layer, ps::CoupledProcesses{P}, state) where {P}
+ initialcondition!(l1::Layer, ps1::CoupledProcesses{P1}, l2::Layer, ps2::CoupledProcesses{P2}, s1, s2) where {P1,P2}
Default implementation of `initialcondition!` for multi-process types. Calls each process in sequence.
"""
diff --git a/src/Strat/tile.jl b/src/Strat/tile.jl
index 8a9278638e4486c442f85bac0e9f6126e023ecdf..7d4512f6716e69debae6c20869b1390eed499df9 100644
--- a/src/Strat/tile.jl
+++ b/src/Strat/tile.jl
@@ -339,9 +339,6 @@ function _collectvars(@nospecialize(comp::StratComponent))
# check for re-definition of diagnostic variables as prognostic
prog_alg = prog_vars ∪ alg_vars
diag_prog = filter(v -> v ∈ prog_alg, diag_vars)
- if !isempty(diag_prog)
- @warn "Variables $(Tuple(map(varname,diag_prog))) declared as both prognostic/algebraic and diagnostic. In-place modifications outside of callbacks may degrade integration accuracy."
- end
# check for conflicting definitions of differential vars
diff_varnames = map(v -> Symbol(:d, varname(v)), prog_alg)
@assert all((isempty(filter(v -> varname(v) == d, all_vars)) for d in diff_varnames)) "Variable names $(Tuple(diff_varnames)) are reserved for differentials."
diff --git a/src/Utils/Utils.jl b/src/Utils/Utils.jl
index 67c66ff457b426aaee4cfe25fe36d881fd1ec51b..2592e668a956f977c943ca3c70342b8ca0f99357 100644
--- a/src/Utils/Utils.jl
+++ b/src/Utils/Utils.jl
@@ -13,7 +13,7 @@ using Unitful
import CryoGrid
import ForwardDiff
-export @xu_str, @Float_str, @Real_str, @Number_str, @UFloat_str, @UT_str, @setscalar
+export @xu_str, @Float_str, @Real_str, @Number_str, @UFloat_str, @UT_str, @setscalar, @threaded
include("macros.jl")
export DistUnit, DistQuantity, TempUnit, TempQuantity, TimeUnit, TimeQuantity
diff --git a/src/Utils/macros.jl b/src/Utils/macros.jl
index 57635929ab86a446723315dd9e16b0c6c95b53c4..ae0d7e0569307291adeed2d94b8ae9821fd89993 100644
--- a/src/Utils/macros.jl
+++ b/src/Utils/macros.jl
@@ -29,7 +29,6 @@ Similar to Unitful.@u_str (i.e. u"kg") but produces the type of the unit rather
on debug mode.
"""
macro UT_str(unit) :(typeof(@u_str($unit))) end
-
"""
Convenience macro for setting scalar (single-element) arrays/vectors. It turns an expression of the form:
`a.b = ...`
@@ -45,3 +44,19 @@ macro setscalar(expr)
$(esc(refexpr))[1] = $(esc(valexpr))
end
end
+"""
+Prepends `expr` with `Threads.@threads` if and only if `Threads.nthreads() > 1`, thus avoiding the overhead of
+`@threads` when running in single-threaded mode.
+
+Credit to @ranocha (Hendrik Ranocha)
+https://discourse.julialang.org/t/overhead-of-threads-threads/53964/22
+"""
+macro threaded(expr)
+ esc(quote
+ if Threads.nthreads() == 1
+ $(expr)
+ else
+ Threads.@threads $(expr)
+ end
+ end)
+end
\ No newline at end of file
diff --git a/src/methods.jl b/src/methods.jl
index dfa03a22ae2f19f36b876300c7e408d3b697da21..0665f6477577020d1be57dd45cb0b166710b7f05 100644
--- a/src/methods.jl
+++ b/src/methods.jl
@@ -1,22 +1,28 @@
+# Core model functions
"""
- variables(::Layer, ::Process)
+ variables(::Layer)
+ variables(::Process)
+ variables(layer::Layer, process::Process)
-Defines variables for a given Process on the given Layer. Implementations should return a `Tuple` of `Var`s.
+Defines variables for a given Process and/or Layer. Implementations should return a `Tuple` of `Var`s.
"""
-variables(::Layer, ::Process) = ()
-"""
- callbacks(::Layer, ::Process)
-
-Defines callbacks for a given Process on the given Layer. Implementations should return a `Tuple` or `Callback`s.
-"""
-callbacks(::Layer, ::Process) = ()
+variables(::Layer) = ()
+variables(::Process) = ()
+variables(layer::Layer, process::Process) = tuple(variables(layer)..., variables(process)...)
"""
+ initialcondition!(::Layer, state)
+ initialcondition!(::Process, state) = nothing
initialcondition!(::Layer, ::Process, state)
-Defines the initial condition for a given Process on the given Layer. `initialcondition!` should write initial values into all relevant
+Defines the initial condition for a given Process and/or Layer. `initialcondition!` should write initial values into all relevant
state variables in `state`.
"""
-initialcondition!(::Layer, ::Process, state) = nothing
+initialcondition!(::Layer, state) = nothing
+initialcondition!(::Process, state) = nothing
+function initialcondition!(layer::Layer, process::Process, state)
+ initialcondition!(layer, state)
+ initialcondition!(process, state)
+end
"""
initialcondition!(::Layer, ::Process, ::Layer, ::Process, state1, state2)
@@ -45,32 +51,6 @@ follows decreasing depth, i.e. the first layer/process is always on top of the s
and separate dispatches must be provided for interactions in reverse order.
"""
interact!(::Layer, ::Process, ::Layer, ::Process, state1, state2) = nothing
-"""
- boundaryflux(bc::BoundaryProcess, b::Union{Top,Bottom}, p::SubSurfaceProcess, sub::SubSurface, s1, s2)
- boundaryflux(s::BoundaryStyle, bc::BoundaryProcess, b::Union{Top,Bottom}, p::SubSurfaceProcess, sub::SubSurface, s1, s2)
-
-Computes the flux dH/dt at the boundary layer. Calls boundaryflux(BoundaryStyle(B),...) to allow for generic implementations by boundary condition type.
-"""
-boundaryflux(bc::BoundaryProcess, b::Union{Top,Bottom}, p::SubSurfaceProcess, sub::SubSurface, s1, s2) = boundaryflux(BoundaryStyle(bc), bc, b, p, sub, s1, s2)
-boundaryflux(s::BoundaryStyle, bc::BoundaryProcess, b::Union{Top,Bottom}, p::SubSurfaceProcess, sub::SubSurface, s1, s2) = error("missing implementation of $(typeof(b)) $(typeof(s)) boundaryflux for $(typeof(bc)) + $(typeof(p)) on $(typeof(sub))")
-"""
- boundaryvalue(bc::BoundaryProcess, b, p, layer, sbc, ssub)
-
-Computes the value of the boundary condition specified by `bc` for the given layer/process combinations.
-"""
-boundaryvalue(bc::BoundaryProcess, b, p, layer, sbc, ssub) = error("missing implementation of boundaryvalue for $(typeof(b)) $(typeof(bc)) on $(typeof(layer)) with $(typeof(p))")
-"""
- criterion(c::Callback, ::Layer, ::Process, state)
-
-Callback criterion/condition. Should return a `Bool` for discrete callbacks and a real number for continuous callbacks.
-"""
-criterion(c::Callback, ::Layer, ::Process, state) = error("missing implementation of criterion for $(typeof(c))")
-"""
- affect!(c::Callback, ::Layer, ::Process, state)
-
-Callback action executed when `criterion` is met (boolean condition for discrete callbacks, zero for continuous callbacks).
-"""
-affect!(c::Callback, ::Layer, ::Process, state) = error("missing implementation of affect! for $(typeof(c))")
"""
observe(::Val{name}, ::Layer, ::Process, state1)
@@ -91,6 +71,21 @@ setup = Tile(stratigraphy, grid, observed=[:meanT])
```
"""
observe(::Val{name}, ::Layer, ::Process, state) where name = nothing
+# Auxiliary functions for generalized boundary implementations
+"""
+ boundaryflux(bc::BoundaryProcess, b::Union{Top,Bottom}, p::SubSurfaceProcess, sub::SubSurface, s1, s2)
+ boundaryflux(s::BoundaryStyle, bc::BoundaryProcess, b::Union{Top,Bottom}, p::SubSurfaceProcess, sub::SubSurface, s1, s2)
+
+Computes the flux dH/dt at the boundary layer. Calls boundaryflux(BoundaryStyle(B),...) to allow for generic implementations by boundary condition type.
+"""
+boundaryflux(bc::BoundaryProcess, b::Union{Top,Bottom}, p::SubSurfaceProcess, sub::SubSurface, s1, s2) = boundaryflux(BoundaryStyle(bc), bc, b, p, sub, s1, s2)
+boundaryflux(s::BoundaryStyle, bc::BoundaryProcess, b::Union{Top,Bottom}, p::SubSurfaceProcess, sub::SubSurface, s1, s2) = error("missing implementation of $(typeof(b)) $(typeof(s)) boundaryflux for $(typeof(bc)) + $(typeof(p)) on $(typeof(sub))")
+"""
+ boundaryvalue(bc::BoundaryProcess, b, p, layer, sbc, ssub)
+
+Computes the value of the boundary condition specified by `bc` for the given layer/process combinations.
+"""
+boundaryvalue(bc::BoundaryProcess, b, p, layer, sbc, ssub) = error("missing implementation of boundaryvalue for $(typeof(b)) $(typeof(bc)) on $(typeof(layer)) with $(typeof(p))")
"""
BoundaryStyle(::Type{T})
@@ -104,6 +99,31 @@ where `BP` is a `BoundaryProcess` that provides the boundary conditions.
"""
BoundaryStyle(::Type{BP}) where {BP<:BoundaryProcess} = error("No style specified for boundary process $BP")
BoundaryStyle(bc::BoundaryProcess) = BoundaryStyle(typeof(bc))
-# default callback style impl
+# Callbacks
+"""
+ callbacks(::Layer, ::Process)
+
+Defines callbacks for a given Process on the given Layer. Implementations should return a `Tuple` or `Callback`s.
+"""
+callbacks(::Layer, ::Process) = ()
+"""
+ criterion(c::Callback, ::Layer, ::Process, state)
+
+Callback criterion/condition. Should return a `Bool` for discrete callbacks and a real number for continuous callbacks.
+"""
+criterion(c::Callback, ::Layer, ::Process, state) = error("missing implementation of criterion for $(typeof(c))")
+"""
+ affect!(c::Callback, ::Layer, ::Process, state)
+
+Callback action executed when `criterion` is met (boolean condition for discrete callbacks, zero for continuous callbacks).
+"""
+affect!(c::Callback, ::Layer, ::Process, state) = error("missing implementation of affect! for $(typeof(c))")
+"""
+ CallbackStyle(::C)
+ CallbackStyle(::Type{<:Callback})
+
+Trait implementation that defines the "style" or type of the given callback as any subtype of `CallbackStyle`,
+for example `Discrete` or `Continuous`.
+"""
CallbackStyle(::C) where {C<:Callback} = CallbackStyle(C)
CallbackStyle(::Type{<:Callback}) = Discrete()
diff --git a/test/Physics/HeatConduction/sfcc_tests.jl b/test/Physics/HeatConduction/sfcc_tests.jl
index cd994f48193bb9e7a9a4f825d01b5fd5c9ff0913..76ae35a172509fab141cd06fc14212811110129d 100644
--- a/test/Physics/HeatConduction/sfcc_tests.jl
+++ b/test/Physics/HeatConduction/sfcc_tests.jl
@@ -6,11 +6,11 @@ using ComponentArrays
@testset "SFCC" begin
Tₘ = 0.0
θres = 0.0
- soil = Soil()
- θw = CryoGrid.Layers.θw(soil.para)
- θp = CryoGrid.Layers.θp(soil.para)
- θm = CryoGrid.Layers.θm(soil.para)
- θo = CryoGrid.Layers.θo(soil.para)
+ soil = Soil(para=soilparameters())
+ θw = Layers.soilcomp(Val{:θw}(), soil.para)
+ θp = Layers.soilcomp(Val{:θp}(), soil.para)
+ θm = Layers.soilcomp(Val{:θm}(), soil.para)
+ θo = Layers.soilcomp(Val{:θo}(), soil.para)
@testset "McKenzie freeze curve" begin
@testset "Sanity checks" begin
f = McKenzie()
@@ -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,Ceff=similar(C),H=H,θl=θl,θw=θw.+zero(T),θm=θm.+zero(T),θo=θo.+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,)
+ state = (T=T,C=C,Ceff=similar(C),H=H,θl=θl,θw=θw.+zero(T),θm=θm.+zero(T),θo=θo.+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,)
+ state = (T=T,C=C,Ceff=similar(C),H=H,θl=θl,θw=θw.+zero(T),θm=θm.+zero(T),θo=θo.+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,soil)
+ state = (T=T,C=C,Ceff=similar(C),H=H,θl=θl,θw=θw.+zero(T),θm=θm.+zero(T),θo=θo.+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,soil)
+ state = (T=T,C=C,Ceff=similar(C),H=H,θl=θl,θw=θw.+zero(T),θm=θm.+zero(T),θo=θo.+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,soil)
+ state = (T=T,C=C,Ceff=similar(C),H=H,θl=θl,θw=θw.+zero(T),θm=θm.+zero(T),θo=θo.+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,soil)
+ state = (T=T,C=C,Ceff=similar(C),H=H,θl=θl,θw=θw.+zero(T),θm=θm.+zero(T),θo=θo.+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,soil)
+ state = (T=T,C=C,Ceff=similar(C),H=H,θl=θl,θw=θw.+zero(T),θm=θm.+zero(T),θo=θo.+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,soil)
+ state = (T=T,C=C,Ceff=similar(C),H=H,θl=θl,θw=θw.+zero(T),θm=θm.+zero(T),θo=θo.+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)
+ state = (T=T_,C=C,Ceff=similar(C),H=H,θl=θl,θw=θw.+zero(T),θm=θm.+zero(T),θo=θo.+zero(T))
SFCC(McKenzie(γ=Param(p.γ[1])))(soil, heat, state)
state.T[1]
end
@@ -244,10 +244,10 @@ function benchmarksfcc()
L = heat.L
# set up multi-grid-cell state vars
T = [-15.0 for i in 1:10]
- θp = Layers.θp(soilp)
- θw = Layers.θw(soilp)
- θm = Layers.θm(soilp)
- θo = Layers.θo(soilp)
+ θw = Layers.soilcomp(Val{:θw}(), soil.para)
+ θp = Layers.soilcomp(Val{:θp}(), soil.para)
+ θm = Layers.soilcomp(Val{:θm}(), soil.para)
+ θo = Layers.soilcomp(Val{:θo}(), soil.para)
θl = f.(T,Tₘ,θres,θp,θw,L,α,n) # set liquid water content according to freeze curve
C = heatcapacity.(soil,θw,θl,θm,θo)
H = let T = T.+14.999,
@@ -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)
+ state = (T=T,C=C,Ceff=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