diff --git a/Project.toml b/Project.toml
index b94f4b7a509784e63acc344b9b9b89ea319f7700..a5c8b3b9aa6f3651532ebc5b3fbdddc46c326ee7 100755
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
name = "CryoGrid"
uuid = "a535b82e-5f3d-4d97-8b0b-d6483f5bebd5"
authors = ["Brian Groenke <brian.groenke@awi.de>", "Moritz Langer <moritz.langer@awi.de>"]
-version = "0.10.4"
+version = "0.11.0"
[deps]
ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66"
@@ -14,6 +14,7 @@ DimensionalData = "0703355e-b756-11e9-17c0-8b28908087d0"
Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
Flatten = "4c728ea3-d9ee-5c9a-9642-b6f7d7dc04fa"
ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+FreezeCurves = "71e4ad71-e4f2-45a3-aa0a-91ffaa9676be"
IfElse = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953"
@@ -23,7 +24,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890"
ModelParameters = "4744a3fa-6c31-4707-899e-a3298e4618ad"
-NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
+NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
PreallocationTools = "d236fae5-4411-538c-8e31-a6e3d9e00b46"
RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
@@ -31,6 +32,7 @@ Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"
SimulationLogs = "e3c6736c-93d9-479d-93f3-96ff5e8836d5"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
StructTypes = "856f2bd8-1eba-4b0a-8007-ebc267875bd4"
Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
@@ -38,8 +40,10 @@ TimeSeries = "9e3dc215-6440-5c97-bce1-76c03772f85e"
Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
[compat]
-julia = ">= 1.6"
-PreallocationTools = ">= 0.3.1"
+FreezeCurves = "0.2"
+OrdinaryDiffEq = "6"
+PreallocationTools = "0.4"
+julia = "1.6,1.7"
[extras]
BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
diff --git a/src/IO/InputOutput.jl b/src/IO/InputOutput.jl
index 8a4ef38d228fc055e43397eed454e6b36f33689c..20ab7081be3ec91d4da5c4fbb967b7ba219708c2 100644
--- a/src/IO/InputOutput.jl
+++ b/src/IO/InputOutput.jl
@@ -1,10 +1,13 @@
module InputOutput
import CryoGrid
+
using CryoGrid.Numerics
+using CryoGrid.Utils
using Base: @propagate_inbounds
using ComponentArrays
+using ConstructionBase
using DataStructures: DefaultDict
using Dates
using DimensionalData
@@ -23,10 +26,10 @@ import DimensionalData: stack
export loadforcings
include("ioutils.jl")
+export CryoGridParams, parameterize
+include("params.jl")
export Forcing, TimeSeriesForcing
include("forcings.jl")
-export CryoGridParams
-include("params.jl")
export CryoGridOutput
include("output.jl")
diff --git a/src/IO/forcings.jl b/src/IO/forcings.jl
index f262035188cb1ebf77d713dac60e1ba0db6a9028..8226fa0ec3c378c326ff57308ef4962a06719aaa 100644
--- a/src/IO/forcings.jl
+++ b/src/IO/forcings.jl
@@ -8,6 +8,7 @@ abstract type Forcing{T,N} end
@inline @propagate_inbounds (forcing::Forcing)(t::DateTime) = forcing(ustrip(u"s", float(Dates.datetime2epochms(t))u"ms"))
# disable flattening for all fields of forcing types by default
Flatten.flattenable(::Type{<:Forcing}, ::Type) = false
+InputOutput.parameterize(f::Forcing; fields...) = f
"""
TimeSeriesForcing{T,A,I}
diff --git a/src/IO/params.jl b/src/IO/params.jl
index aeeb5a7283934b62e4d135f8839f08d1ec1a6fce..2f181bab0a396dc7d999021dcef1f6ebd2bb005e 100644
--- a/src/IO/params.jl
+++ b/src/IO/params.jl
@@ -40,14 +40,6 @@ function Base.show(io::IO, ::MIME"text/plain", ps::CryoGridParams{T}) where T
end
Tables.columns(ps::CryoGridParams) = Tables.columns(ps.table)
Tables.rows(ps::CryoGridParams) = Tables.rows(ps.table)
-function CryoGridParams(obj; full_metadata=false)
- m = Model(obj)
- if full_metadata
- m[:idx] = 1:length(m)
- m = _setparafields(m)
- end
- return CryoGridParams(m)
-end
function _setparafields(m::Model)
function _setparafield(name, type::Type, para::CryoGrid.Parameterization)
if length(ModelParameters.params(para)) > 0
@@ -68,3 +60,37 @@ function _setparafields(m::Model)
newparent = Flatten.reconstruct(parent(m), updated_parameterizations, CryoGrid.Parameterization)
return Model(newparent)
end
+function CryoGridParams(obj; full_metadata=false)
+ m = Model(obj)
+ if full_metadata
+ m[:idx] = 1:length(m)
+ m = _setparafields(m)
+ end
+ return CryoGridParams(m)
+end
+"""
+ parameterize(x::T) where {T}
+ parameterize(x::Unitful.AbstractQuantity; fields...)
+ parameterize(x::Number; fields...)
+ parameterize(p::Param; ignored...)
+
+If `x` is a numeric type, `x` will be wrapped in a `ModelParameters.Param`, including a `units` field if `x` has units.
+If `x` is a `Param` type, `x` will be returned as-is.
+If `x` is some other type, `x` will be recursively unpacked and `parameterize` called on each field. It may be necessary
+or desirable for some types to override `parameterize` to define custom behavior, e.g. if only some fields should be
+parameterized.
+"""
+function parameterize(x::Unitful.AbstractQuantity; fields...)
+ let x = normalize_units(x);
+ Param(ustrip(x); untis=unit(x), fields...)
+ end
+end
+function parameterize(x::T; fields...) where {T}
+ _parameterize(x) = parameterize(x; fields...)
+ new_fields = map(_parameterize, getfields(x))
+ ctor = ConstructionBase.constructorof(T)
+ return ctor(new_fields...)
+end
+parameterize(x::Number; fields...) = Param(x; fields...)
+parameterize(x::Param; ignored...) = x
+parameterize(x::AbstractArray; ignored...) = x
diff --git a/src/Numerics/init.jl b/src/Numerics/init.jl
index dada3b3fddeaca70e008d558f2c161a58dba533d..6ee788520ab1a54037cf3ea49ba38bb5997654ee 100644
--- a/src/Numerics/init.jl
+++ b/src/Numerics/init.jl
@@ -25,7 +25,7 @@ struct InterpInitializer{varname,P,I,E} <: VarInitializer{varname}
extrap::E
InterpInitializer(varname::Symbol, profile::P, interp::I=Linear(), extrap::E=Flat()) where {P<:Profile,I,E} = new{varname,P,I,E}(profile, interp, extrap)
end
-function (init::InterpInitializer{var})(state) where var
+function (init::InterpInitializer{var})(layer, state) where var
profile, interp, extrap = init.profile, init.interp, init.extrap
depths = collect(map(knot -> dustrip(knot.depth), profile.knots))
u = getproperty(state, var)
@@ -51,6 +51,6 @@ Base.getindex(init::InterpInitializer{var}, itrv::Interval) where var = InterpIn
Convenience constructor for `VarInitializer` that selects the appropriate initializer type based on the arguments.
"""
-initializer(varname::Symbol, x::Number) = FunctionInitializer(varname, state -> getproperty(state, varname) .= x)
+initializer(varname::Symbol, x::Number) = FunctionInitializer(varname, (layer,state) -> getproperty(state, varname) .= x)
initializer(varname::Symbol, f::Function) = FunctionInitializer(varname, f)
initializer(varname::Symbol, profile::Profile, interp=Linear(), extrap=Flat()) = InterpInitializer(varname, profile, interp, extrap)
diff --git a/src/Numerics/math.jl b/src/Numerics/math.jl
index cb562acd79d33054efa886bedac9cd8135040e1a..e179282c1cdcfc3c7a68f4a7161eb908ea7cfbbe 100644
--- a/src/Numerics/math.jl
+++ b/src/Numerics/math.jl
@@ -1,77 +1,87 @@
-∇(f, x) = ∇(typeof(x), f, x)
-function ∇(::Type{T}, f, x) where {T}
- res = ForwardDiff.derivative!(ForwardDiff.DiffResult(zero(T), x), f, x)
+function ∇(f::F, x::Number) where {F}
+ res = ForwardDiff.derivative!(ForwardDiff.DiffResult(zero(x), zero(x)), f, x)
return res.value, res.derivs[1]
end
-
+function ∇(f::F, x::AbstractArray) where {F}
+ res = ForwardDiff.gradient!(ForwardDiff.DiffResult(eltype(x), zero(x)), f, x)
+ return res.value, res.derivs
+end
+# Flux calculations
+@propagate_inbounds @inline _flux_kernel(xâ‚, xâ‚‚, Δx, k) = -k*(xâ‚‚ - xâ‚)/ Δx
+@propagate_inbounds @inline function _flux!(j::AbstractVector, xâ‚::AbstractVector, xâ‚‚::AbstractVector, Δx::AbstractVector, k::AbstractVector, ::Val{use_turbo}) where {use_turbo}
+ @. j = _flux_kernel(xâ‚, xâ‚‚, Δx, k)
+end
+@propagate_inbounds @inline function _flux!(j::AbstractVector{Float64}, xâ‚::AbstractVector{Float64}, xâ‚‚::AbstractVector{Float64}, Δx::AbstractVector{Float64}, k::AbstractVector{Float64}, ::Val{true})
+ @turbo @. j = _flux_kernel(xâ‚, xâ‚‚, Δx, k)
+end
"""
- finitediff!(∂x::AbstractVector, x::AbstractVector, Δ::AbstractVector)
+ flux!(j::AbstractVector, x::AbstractVector, Δx::AbstractVector, k::AbstractVector)
-First order forward finite difference operator.
+Calculates the first-order, non-linear spatial flux over a discretized variable `x` with conductivity `k`.
+`x` is assumed to have shape `(N,)`, `Δx` shape `(N-1,)`, and `j` and `k` shape `(N+1,)` such that `j[2:end-1]` represents
+the fluxes over the inner grid cell faces. Fluxes are added to existing values in `j`.
"""
-function finitediff!(∂x::AbstractVector, x::AbstractVector, Δ::AbstractVector)
- @inbounds let xâ‚ = (@view x[1:end-1]),
- xâ‚‚ = (@view x[2:end]);
- @. ∂x = (xâ‚‚ - xâ‚) / Δ
+function flux!(j::AbstractVector, x::AbstractVector, Δx::AbstractVector, k::AbstractVector)
+ @inbounds let xâ‚ = @view(x[1:end-1]),
+ xâ‚‚ = @view(x[2:end]),
+ k = @view(k[2:end-1]),
+ j = @view(j[2:end-1]);
+ _flux!(j, xâ‚, xâ‚‚, Δx, k, Val{USE_TURBO}())
end
end
+# Divergence
+@propagate_inbounds @inline _div_kernel(jâ‚, jâ‚‚, Δj) = (jâ‚ - jâ‚‚) / Δj
+@propagate_inbounds @inline function _div!(dx::AbstractVector, jâ‚::AbstractVector, jâ‚‚::AbstractVector, Δj::AbstractVector, ::Val{use_turbo}) where {use_turbo}
+ @. dx = _div_kernel(jâ‚, jâ‚‚, Δj)
+end
+@propagate_inbounds @inline function _div!(dx::AbstractVector{Float64}, jâ‚::AbstractVector{Float64}, jâ‚‚::AbstractVector{Float64}, Δj::AbstractVector{Float64}, ::Val{true})
+ @turbo @. dx = _div_kernel(jâ‚, jâ‚‚, Δj)
+end
"""
- lineardiffusion!(∂y::AbstractVector, x::AbstractVector, Δ::AbstractVector, k::Number)
+ divergence!(dx::AbstractVector, j::AbstractVector, Δj::AbstractVector)
-Second order Laplacian with constant diffusion k.
+Calculates the first-order divergence over a 1D flux vector field `j` and grid cell lengths `Δj`. Divergences are added to existing values in `dx`.
"""
-function lineardiffusion!(∂y::AbstractVector, x::AbstractVector, Δ::AbstractVector, k::Number)
- @inbounds let xâ‚ = (@view x[1:end-2]),
- xâ‚‚ = (@view x[2:end-1]),
- x₃ = (@view x[3:end]),
- Δ₠= (@view Δ[1:end-1]),
- Δ₂ = (@view Δ[2:end]);
- @. ∂y = k*((x₃ - xâ‚‚)/Δ₂ - (xâ‚‚-xâ‚)/Δâ‚)/Δâ‚
+function divergence!(dx::AbstractVector, j::AbstractVector, Δj::AbstractVector)
+ @inbounds let jâ‚ = @view(j[1:end-1]),
+ jâ‚‚ = @view(j[2:end]);
+ _div!(dx, jâ‚, jâ‚‚, Δj, Val{USE_TURBO}())
end
end
-"""
- nonlineardiffusion!(∂y, x, Δx, k, Δk)
-
-Second order Laplacian with non-linear diffusion operator, `k`. Accelerated using `LoopVectorization.@turbo` for `Float64` vectors.
-"""
-function nonlineardiffusion!(∂y::AbstractVector, x::AbstractVector, Δx::AbstractVector, k::AbstractVector, Δk::AbstractVector)
- @inbounds let xâ‚ = (@view x[1:end-2]),
- xâ‚‚ = (@view x[2:end-1]),
- x₃ = (@view x[3:end]),
- kâ‚ = (@view k[1:end-1]),
- kâ‚‚ = (@view k[2:end]),
- Δx₠= (@view Δx[1:end-1]),
- Δx₂ = (@view Δx[2:end]);
- nonlineardiffusion!(∂y, xâ‚, xâ‚‚, x₃, kâ‚, kâ‚‚, Δxâ‚, Δxâ‚‚, Δk)
+# non-linear diffusion
+@propagate_inbounds @inline function _nonlineardiffusion_kernel(xâ‚, xâ‚‚, Δx, jâ‚, k, Δj)
+ jâ‚‚ = _flux_kernel(xâ‚, xâ‚‚, Δx, k)
+ div = _div_kernel(jâ‚, jâ‚‚, Δj)
+ return jâ‚‚, div
+end
+@propagate_inbounds @inline function _nonlineardiffusion!(dx::AbstractVector, j::AbstractVector, x::AbstractVector, Δx::AbstractVector, k::AbstractVector, Δk::AbstractVector, i::Integer)
+ let xâ‚ = x[i-1],
+ xâ‚‚ = x[i],
+ jâ‚ = j[i-1],
+ k = k[i],
+ δx = Δx[i-1],
+ δj = Δk[i-1];
+ jâ‚‚, divxâ‚ = _nonlineardiffusion_kernel(xâ‚, xâ‚‚, δx, jâ‚, k, δj)
+ j[i] += jâ‚‚
+ dx[i-1] += divxâ‚
end
end
"""
- nonlineardiffusion(xâ‚, xâ‚‚, x₃, kâ‚, kâ‚‚, Δxâ‚, Δxâ‚‚, Δk)
+ nonlineardiffusion!(dx::AbstractVector, j::AbstractVector, x::AbstractVector, Δx::AbstractVector, k::AbstractVector, Δk::AbstractVector)
-Scalar second order Laplacian with non-linear diffusion operator, `k`.
+Fast alternative to `flux!` and `divergence!` which computes fluxes and divergences (via `_flux_kernel` and `_div_kernel`) in a single pass. Note, however, that
+loop vectorization with `@turbo` is not possible because of necessary loop-carried dependencies. Fluxes and divergences are added to the existing values stored
+in `j` and `dx`.
"""
-@inline nonlineardiffusion(xâ‚, xâ‚‚, x₃, kâ‚, kâ‚‚, Δxâ‚, Δxâ‚‚, Δk) = (kâ‚‚*(x₃ - xâ‚‚)/Δxâ‚‚ - kâ‚*(xâ‚‚ - xâ‚)/Δxâ‚)/Δk
-@propagate_inbounds function nonlineardiffusion!(∂y, xâ‚, xâ‚‚, x₃, kâ‚, kâ‚‚, Δxâ‚, Δxâ‚‚, Δk)
- @. ∂y = nonlineardiffusion(xâ‚, xâ‚‚, x₃, kâ‚, kâ‚‚, Δxâ‚, Δxâ‚‚, Δk)
-end
-@propagate_inbounds function nonlineardiffusion!(
- ∂y::AbstractVector{Float64},
- xâ‚::AbstractVector{Float64},
- xâ‚‚::AbstractVector{Float64},
- x₃::AbstractVector{Float64},
- kâ‚::AbstractVector{Float64},
- kâ‚‚::AbstractVector{Float64},
- Δxâ‚::AbstractVector{Float64},
- Δx₂::AbstractVector{Float64},
- Δk::AbstractVector{Float64},
-)
- if USE_TURBO
- @turbo @. ∂y = nonlineardiffusion(xâ‚, xâ‚‚, x₃, kâ‚, kâ‚‚, Δxâ‚, Δxâ‚‚, Δk)
- else
- @. ∂y = nonlineardiffusion(xâ‚, xâ‚‚, x₃, kâ‚, kâ‚‚, Δxâ‚, Δxâ‚‚, Δk)
+function nonlineardiffusion!(dx::AbstractVector, j::AbstractVector, x::AbstractVector, Δx::AbstractVector, k::AbstractVector, Δk::AbstractVector)
+ @inbounds for i in 2:length(j)-1
+ _nonlineardiffusion!(dx, j, x, Δx, k, Δk, i)
end
+ # compute divergence for last grid cell
+ @inbounds dx[end] += _div_kernel(j[end-1], j[end], Δk[end])
+ return nothing
end
-
+# other helper functions
"""
harmonicmean(xâ‚, xâ‚‚, wâ‚, wâ‚‚)
@@ -169,38 +179,3 @@ softplusinv(x) = let x = clamp(x, eps(), Inf); IfElse.ifelse(x > 34, x, log(exp(
# sqrt ∘ plusone == x -> sqrt(x. + 1.0)
minusone(x) = x .- one.(x)
plusone(x) = x .+ one.(x)
-
-# Function tabulation
-"""
- Tabulated(f, argknots...)
-
-Alias for `tabulate` intended for function types.
-"""
-Tabulated(f, argknots...) = tabulate(f, argknots...)
-"""
- tabulate(f, argknots::Pair{Symbol,<:Union{Number,AbstractArray}}...)
-
-Tabulates the given function `f` using a linear, multi-dimensional interpolant.
-Knots should be given as pairs `:arg => A` where `A` is a `StepRange` or `Vector`
-of input values (knots) at which to evaluate the function. `A` may also be a
-`Number`, in which case a pseudo-point interpolant will be used (i.e valid on
-`[A,A+ϵ]`). No extrapolation is provided by default but can be configured via
-`Interpolations.extrapolate`.
-"""
-function tabulate(f, argknots::Pair{Symbol,<:Union{Number,AbstractArray}}...)
- initknots(a::AbstractArray) = Interpolations.deduplicate_knots!(a)
- initknots(x::Number) = initknots([x,x])
- interp(::AbstractArray) = Gridded(Linear())
- interp(::Number) = Gridded(Constant())
- extrap(::AbstractArray) = Flat()
- extrap(::Number) = Throw()
- names = map(first, argknots)
- # get knots for each argument, duplicating if only one value is provided
- knots = map(initknots, map(last, argknots))
- f_argnames = Utils.argnames(f)
- @assert all(map(name -> name ∈ names, f_argnames)) "Missing one or more arguments $f_argnames in $f"
- arggrid = Iterators.product(knots...)
- # evaluate function construct interpolant
- f = extrapolate(interpolate(Tuple(knots), map(Base.splat(f), arggrid), map(interp ∘ last, argknots)), map(extrap ∘ last, argknots))
- return f
-end
diff --git a/src/Numerics/variables.jl b/src/Numerics/variables.jl
index 880cf7900e807e29d7f0142185c3d8f8f8d36aa9..8a4a068ca87ea1b78de3623004a0f50078fb6d23 100644
--- a/src/Numerics/variables.jl
+++ b/src/Numerics/variables.jl
@@ -37,16 +37,25 @@ end
"""
Delta{dname,name,S,T,units,domain} <: Var{dname,S,T,units,domain}
-Defines a "delta" term `du` for variable `u`, which is the time-derivative or flux for prognostic variables and
+Defines a "delta" term `du` for variable `u`, which is the time-derivative/divergence for prognostic variables and
the residual for algebraic variables.
"""
struct Delta{dname,name,S,T,units,domain} <: Var{dname,S,T,units,domain}
dim::S
- Delta(::Symbol, dims::OnGrid, args...) = error("Off-cell prognostic/algebraic spatial variables are not currently supported.")
Delta(dname::Symbol, name::Symbol, dims::Union{<:Shape,OnGrid{Cells,typeof(identity)}}, units=NoUnits, ::Type{T}=Float64; domain=-Inf..Inf) where {T} = new{dname,name,typeof(dims),T,units,domain}(dims)
Delta(var::Prognostic{name,S,T,units,domain}) where {name,S,T,units,domain} = let dims=vardims(var); new{Symbol(:d,name),name,typeof(dims),T,upreferred(units)/u"s",domain}(dims) end
Delta(var::Algebraic{name,S,T,units,domain}) where {name,S,T,units,domain} = let dims=vardims(var); new{Symbol(:d,name),name,typeof(dims),T,units,domain}(dims) end
end
+"""
+ Flux{jname,name,S,T,units,domain} <: Var{jname,S,T,units,domain}
+
+Defines a "flux" term `ju` for prognostic variable `u`, which is typically a function of the gradient w.r.t space.
+"""
+struct Flux{jname,name,S,T,units,domain} <: Var{jname,S,T,units,domain}
+ dim::S
+ Flux(jname::Symbol, name::Symbol, dims::Union{<:Shape,OnGrid{Cells,typeof(identity)}}, units=NoUnits, ::Type{T}=Float64; domain=-Inf..Inf) where {T} = new{jname,name,typeof(dims),T,units,domain}(dims)
+ Flux(var::Prognostic{name,OnGrid{Cells,typeof(identity)},T,units,domain}) where {name,T,units,domain} = let dims=OnGrid(Edges, vardims(var).f); new{Symbol(:j,name),name,typeof(dims),T,upreferred(units)/u"s",domain}(dims) end
+end
"""
Diagnostic{name,S,T,units,domain} <: Var{name,S,T,units,domain}
@@ -61,6 +70,7 @@ ConstructionBase.constructorof(::Type{Prognostic{name,S,T,units,domain}}) where
ConstructionBase.constructorof(::Type{Diagnostic{name,S,T,units,domain}}) where {name,S,T,units,domain} = s -> Diagnostic(name, s, units, T; domain)
ConstructionBase.constructorof(::Type{Algebraic{name,S,T,units,domain}}) where {name,S,T,units,domain} = s -> Algebraic(name, s, units, T; domain)
ConstructionBase.constructorof(::Type{Delta{dname,name,S,T,units,domain}}) where {dname,name,S,T,units,domain} = s -> Delta(dname, name, s, units, T; domain)
+ConstructionBase.constructorof(::Type{Flux{jname,name,S,T,units,domain}}) where {jname,name,S,T,units,domain} = s -> Flux(jname, name, s, units, T; domain)
==(::Var{N1,S1,T1,u1,d1},::Var{N2,S2,T2,u2,d2}) where {N1,N2,S1,S2,T1,T2,u1,u2,d1,d2} = (N1 == N2) && (S1 == S2) && (T1 == T2) && (u1 == u2) && (d1 == d2)
varname(::Var{name}) where {name} = name
varname(::Type{<:Var{name}}) where {name} = name
diff --git a/src/Numerics/varstates.jl b/src/Numerics/varstates.jl
index f15779c861c80a641e3705c0a774c18667fffbfb..dcb0b574d8c2f2f9735282e291e0ed5b07059a5a 100644
--- a/src/Numerics/varstates.jl
+++ b/src/Numerics/varstates.jl
@@ -97,8 +97,10 @@ function VarStates(vars::GroupedVars, D::Numerics.AbstractDiscretization, chunks
diagvars = map(group -> filter(isdiagnostic, group), vars)
progvars = map(group -> filter(isprognostic, group), vars)
algvars = map(group -> filter(isalgebraic, group), vars)
- # create variables for gradients/fluxes
+ # create variables for time delta variables (divergence/residual)
dpvars = map(group -> map(Delta, filter(var -> isalgebraic(var) || isprognostic(var), group)), vars)
+ # create variables for fluxes (spatial gradients)
+ jpvars = map(group -> map(Flux, filter(var -> isprognostic(var) && isongrid(var), group)), vars)
allprogvars = tuplejoin(_flatten(progvars), _flatten(algvars)) # flattened prognostic/algebraic variable group
# TODO: not a currently necessary use case, but could be supported by partitioning the state array further by layer/group name;
# this would require passing the name to discretize(...)
@@ -111,9 +113,9 @@ function VarStates(vars::GroupedVars, D::Numerics.AbstractDiscretization, chunks
freediagvars = map(group -> filter(!isongrid, group), diagvars)
freediagstate = map(group -> (;map(v -> varname(v) => DiffCache(varname(v), discretize(A, D, v), chunksize), group)...), freediagvars)
# build gridded diagnostic state vectors
- griddiagvars = filter(isongrid, _flatten(diagvars))
+ griddiagvars = Tuple(unique(filter(isongrid, _flatten(map(tuplejoin, diagvars, jpvars)))))
griddiagstate = map(v -> varname(v) => DiffCache(varname(v), discretize(A, D, v), chunksize), griddiagvars)
- # join prognostic variables with flux variables, then build nested named tuples in each group with varnames as keys
- allvars = map(vars -> NamedTuple{map(varname, vars)}(vars), map(tuplejoin, vars, dpvars))
+ # join prognostic variables with delta and flux variables, then build nested named tuples in each group with varnames as keys
+ allvars = map(vars -> NamedTuple{map(varname, vars)}(vars), map(tuplejoin, vars, dpvars, jpvars))
VarStates(uproto, D, allvars, (;freediagstate...), (;griddiagstate...))
end
diff --git a/src/Physics/HeatConduction/HeatConduction.jl b/src/Physics/HeatConduction/HeatConduction.jl
index fa37f1d6e7cc49345101694e764d13f268bd0e83..cdcde0c97aa4d581782608b817eae95e99aa56e6 100644
--- a/src/Physics/HeatConduction/HeatConduction.jl
+++ b/src/Physics/HeatConduction/HeatConduction.jl
@@ -10,6 +10,7 @@ using CryoGrid.Utils
using Base: @propagate_inbounds, @kwdef
using IfElse
+using FreezeCurves: FreezeCurves, FreezeCurve, FreeWater
using ModelParameters
using Unitful
@@ -19,9 +20,6 @@ import CryoGrid.Physics
export Heat, TemperatureProfile
export FreeWater, FreezeCurve, freezecurve
-abstract type FreezeCurve end
-struct FreeWater <: FreezeCurve end
-
"""
TemperatureProfile(pairs::Pair{<:Union{DistQuantity,Param},<:Union{TempQuantity,Param}}...)
@@ -38,24 +36,17 @@ abstract type HeatParameterization end
struct Enthalpy <: HeatParameterization end
struct Temperature <: HeatParameterization end
-abstract type StepLimiter end
-@kwdef struct CFL <: StepLimiter
- fallback_dt::Float64 = 60.0 # fallback dt [s]
+@Base.kwdef struct ThermalProperties{Tconsts,TL,Tkw,Tki,Tka,Tcw,Tci,Tca}
+ consts::Tconsts = Physics.Constants()
+ L::TL = consts.Ïw*consts.Lsl
+ kw::Tkw = 0.57u"W/m/K" # thermal conductivity of water [Hillel(1982)]
+ ki::Tki = 2.2u"W/m/K" # thermal conductivity of ice [Hillel(1982)]
+ ka::Tka = 0.025u"W/m/K" # air [Hillel(1982)]
+ cw::Tcw = 4.2e6u"J/K/m^3" # heat capacity of water
+ ci::Tci = 1.9e6u"J/K/m^3" # heat capacity of ice
+ ca::Tca = 0.00125e6u"J/K/m^3" # heat capacity of air
end
-ThermalProperties(
- consts=Physics.Constants();
- Ïw = consts.Ïw,
- Lf = consts.Lf,
- L = consts.Ïw*consts.Lf,
- kw = Param(0.57, units=u"W/m/K"), # thermal conductivity of water [Hillel(1982)]
- ki = Param(2.2, units=u"W/m/K"), # thermal conductivity of ice [Hillel(1982)]
- ka = Param(0.025, units=u"W/m/K"), # air [Hillel(1982)]
- cw = Param(4.2e6, units=u"J/K/m^3"), # heat capacity of water
- ci = Param(1.9e6, units=u"J/K/m^3"), # heat capacity of ice
- ca = Param(0.00125e6, units=u"J/K/m^3"), # heat capacity of air
-) = (; Ïw, Lf, L, kw, ki, ka, cw, ci, ca)
-
struct Heat{Tfc<:FreezeCurve,TPara<:HeatParameterization,Tdt,Tinit,TProp} <: SubSurfaceProcess
para::TPara
prop::TProp
@@ -68,7 +59,7 @@ Heat(var::Symbol=:H; kwargs...) = Heat(Val{var}(); kwargs...)
Heat(::Val{:H}; kwargs...) = Heat(Enthalpy(); kwargs...)
Heat(::Val{:T}; kwargs...) = Heat(Temperature(); kwargs...)
Heat(para::Enthalpy; freezecurve=FreeWater(), prop=ThermalProperties(), dtlim=nothing, init=nothing) = Heat(para, prop, deepcopy(freezecurve), dtlim, init)
-Heat(para::Temperature; freezecurve, prop=ThermalProperties(), dtlim=CFL(), init=nothing) = Heat(para, prop, deepcopy(freezecurve), dtlim, init)
+Heat(para::Temperature; freezecurve, prop=ThermalProperties(), dtlim=Physics.CFL(), init=nothing) = Heat(para, prop, deepcopy(freezecurve), dtlim, init)
# getter functions
thermalproperties(heat::Heat) = heat.prop
diff --git a/src/Physics/HeatConduction/heat.jl b/src/Physics/HeatConduction/heat.jl
index 5692abfb61075ca1b574f4eeba3a224c25b370b8..2cae8271eab3cb51a25bb879196da7882453af55 100644
--- a/src/Physics/HeatConduction/heat.jl
+++ b/src/Physics/HeatConduction/heat.jl
@@ -127,56 +127,23 @@ CryoGrid.variables(heat::Heat{<:FreezeCurve,Temperature}) = (
Common variable definitions for all heat implementations.
"""
heatvariables(::Heat) = (
- Diagnostic(:dH_upper, Scalar, u"J/K/m^2"),
- Diagnostic(:dH_lower, Scalar, u"J/K/m^2"),
Diagnostic(:dHdT, OnGrid(Cells), u"J/K/m^3", domain=0..Inf),
Diagnostic(:C, OnGrid(Cells), u"J/K/m^3"),
Diagnostic(:k, OnGrid(Edges), u"W/m/K"),
Diagnostic(:kc, OnGrid(Cells), u"W/m/K"),
Diagnostic(:θw, OnGrid(Cells), domain=0..1),
)
-"""
- heatconduction!(∂H,T,ΔT,k,Δk)
-
-1-D heat conduction/diffusion given T, k, and their deltas. Resulting enthalpy gradient is stored in ∂H.
-Note that this function does not perform bounds checking. It is up to the user to ensure that all variables are
-arrays of the correct length.
-"""
-function heatconduction!(∂H,T,ΔT,k,Δk)
- # upper boundary
- @inbounds ∂H[1] += let T₂=T[2],
- Tâ‚=T[1],
- k=k[2],
- ϵ=ΔT[1],
- δ=Δk[1];
- k*(Tâ‚‚-Tâ‚)/ϵ/δ
- end
- # diffusion on non-boundary cells
- @inbounds let T = T,
- k = (@view k[2:end-1]),
- Δk = (@view Δk[2:end-1]),
- ∂H = (@view ∂H[2:end-1]);
- nonlineardiffusion!(∂H, T, ΔT, k, Δk)
- end
- # lower boundary
- @inbounds ∂H[end] += let T₂=T[end],
- Tâ‚=T[end-1],
- k=k[end-1],
- ϵ=ΔT[end],
- δ=Δk[end];
- -k*(Tâ‚‚-Tâ‚)/ϵ/δ
- end
- return nothing
+function resetfluxes!(sub::SubSurface, heat::Heat, state)
+ # Reset energy fluxes to zero; this is redundant when H is the prognostic variable
+ # but necessary when it is not.
+ @. state.dH = zero(eltype(state.dH))
+ @. state.jH = zero(eltype(state.jH))
end
"""
Diagonstic step for heat conduction (all state configurations) on any subsurface layer.
"""
function CryoGrid.diagnosticstep!(sub::SubSurface, heat::Heat, state)
- # Reset energy flux to zero; this is redundant when H is the prognostic variable
- # but necessary when it is not.
- @. state.dH = zero(eltype(state.dH))
- @. state.dH_upper = zero(eltype(state.dH_upper))
- @. state.dH_lower = zero(eltype(state.dH_lower))
+ resetfluxes!(sub, heat, state)
# Evaluate freeze/thaw processes
freezethaw!(sub, heat, state)
# Update thermal conductivity
@@ -196,20 +163,16 @@ end
Generic top interaction. Computes flux dH at top cell.
"""
function CryoGrid.interact!(top::Top, bc::HeatBC, sub::SubSurface, heat::Heat, stop, ssub)
- Δk = CryoGrid.thickness(sub, ssub, first) # using `thickness` allows for generic layer implementations
# boundary flux
- @setscalar ssub.dH_upper = boundaryflux(bc, top, heat, sub, stop, ssub)
- @inbounds ssub.dH[1] += getscalar(ssub.dH_upper) / Δk
+ ssub.jH[1] += boundaryflux(bc, top, heat, sub, stop, ssub)
return nothing # ensure no allocation
end
"""
Generic bottom interaction. Computes flux dH at bottom cell.
"""
function CryoGrid.interact!(sub::SubSurface, heat::Heat, bot::Bottom, bc::HeatBC, ssub, sbot)
- Δk = CryoGrid.thickness(sub, ssub, last) # using `thickness` allows for generic layer implementations
- # boundary flux
- @setscalar ssub.dH_lower = boundaryflux(bc, bot, heat, sub, sbot, ssub)
- @inbounds ssub.dH[end] += getscalar(ssub.dH_lower) / Δk
+ # boundary flux; here we flip the sign since a positive flux is by convention downward
+ ssub.jH[end] -= boundaryflux(bc, bot, heat, sub, sbot, ssub)
return nothing # ensure no allocation
end
"""
@@ -226,18 +189,17 @@ function CryoGrid.interact!(sub1::SubSurface, ::Heat, sub2::SubSurface, ::Heat,
Δ₂ = Δk₂[1];
harmonicmean(kâ‚, kâ‚‚, Δâ‚, Δ₂)
end
- # calculate heat flux between cells
+ # calculate heat flux between cells (positive downward)
Qáµ¢ = @inbounds let zâ‚ = CryoGrid.midpoint(sub1, s1, last),
zâ‚‚ = CryoGrid.midpoint(sub2, s2, first),
δ = zâ‚‚ - zâ‚;
- k*(s2.T[1] - s1.T[end]) / δ
+ -k*(s2.T[1] - s1.T[end]) / δ
end
- # diagnostics
- @setscalar s1.dH_lower = Qáµ¢
- @setscalar s2.dH_upper = -Qáµ¢
- # add fluxes scaled by grid cell size
- @inbounds s1.dH[end] += Qáµ¢ / Δkâ‚[end]
- @inbounds s2.dH[1] += -Qᵢ / Δk₂[1]
+ # set inner boundary flux
+ s1.jH[end] += Qáµ¢
+ # these should already be equal since jH[1] on s2 and jH[end] on s1 should point to the same array index;
+ # but we set them explicitly just to be safe.
+ s2.jH[1] = s1.jH[end]
return nothing # ensure no allocation
end
"""
@@ -245,18 +207,19 @@ Prognostic step for heat conduction (enthalpy) on subsurface layer.
"""
function CryoGrid.prognosticstep!(::SubSurface, ::Heat{<:FreezeCurve,Enthalpy}, state)
Δk = Δ(state.grids.k) # cell sizes
- ΔT = Δ(state.grids.T)
- # Diffusion on non-boundary cells
- heatconduction!(state.dH, state.T, ΔT, state.k, Δk)
+ ΔT = Δ(state.grids.T) # midpoint distances
+ # compute internal fluxes and non-linear diffusion assuming boundary fluxes have been set
+ nonlineardiffusion!(state.dH, state.jH, state.T, ΔT, state.k, Δk)
+ return nothing
end
"""
Prognostic step for heat conduction (temperature) on subsurface layer.
"""
function CryoGrid.prognosticstep!(sub::SubSurface, ::Heat{<:FreezeCurve,Temperature}, state)
Δk = Δ(state.grids.k) # cell sizes
- ΔT = Δ(state.grids.T)
- # Diffusion on non-boundary cells
- heatconduction!(state.dH, state.T, ΔT, state.k, Δk)
+ ΔT = Δ(state.grids.T) # midpoint distances
+ # compute internal fluxes and non-linear diffusion assuming boundary fluxes have been set
+ nonlineardiffusion!(state.dH, state.jH, state.T, ΔT, state.k, Δk)
# Compute temperature flux by dividing by dHdT;
# dHdT should be computed by the freeze curve.
@inbounds @. state.dT = state.dH / state.dHdT
@@ -280,11 +243,12 @@ end
Tsub=ssub.T[end],
k=ssub.k[end],
δ=Δk/2; # distance to boundary
- -k*(Tsub-Tlower)/δ
+ # note again the inverted sign; positive here means *upward from* the bottom boundary
+ k*(Tlower-Tsub)/δ
end
end
# CFL not defined for free-water freeze curve
-CryoGrid.timestep(::SubSurface, heat::Heat{FreeWater,Enthalpy,CFL}, state) = Inf
+CryoGrid.timestep(::SubSurface, heat::Heat{FreeWater,Enthalpy,Physics.CFL}, state) = Inf
"""
timestep(::SubSurface, ::Heat{Tfc,TPara,CFL}, state) where {TPara}
@@ -292,7 +256,7 @@ Implementation of `timestep` for `Heat` using the Courant-Fredrichs-Lewy conditi
defined as: Δt_max = u*Δx^2, where`u` is the "characteristic velocity" which here
is taken to be the diffusivity: `dHdT / kc`.
"""
-@inline function CryoGrid.timestep(::SubSurface, heat::Heat{Tfc,TPara,CFL}, state) where {Tfc,TPara}
+@inline function CryoGrid.timestep(::SubSurface, heat::Heat{Tfc,TPara,Physics.CFL}, state) where {Tfc,TPara}
Δx = Δ(state.grid)
dtmax = Inf
@inbounds for i in 1:length(Δx)
@@ -307,21 +271,14 @@ is taken to be the diffusivity: `dHdT / kc`.
end
# Free water freeze curve
@inline function enthalpyinv(sub::SubSurface, heat::Heat{FreeWater,Enthalpy}, state, i)
- f_hc = partial(heatcapacity, liquidwater, sub, heat, state, i)
- return enthalpyinv(heat.freezecurve, f_hc, state.H[i], heat.prop.L, f_hc.θwi)
+ hc = partial(heatcapacity, liquidwater, sub, heat, state, i)
+ return enthalpyinv(heat.freezecurve, hc, state.H[i], hc.θwi, heat.prop.L)
end
-@inline function enthalpyinv(::FreeWater, f_hc::Function, H, L, θtot)
- let θtot = max(1e-8, θtot),
- Lθ = L*θtot,
- I_t = H > Lθ,
- I_f = H <= 0.0,
- I_c = (H > 0.0) && (H <= Lθ);
- # compute liquid water content -> heat capacity -> temperature
- θw = (I_c*(H/Lθ) + I_t)θtot
- C = f_hc(θw)
- T = (I_t*(H-Lθ) + I_f*H)/C
- return T, θw, C
- end
+@inline function enthalpyinv(::FreeWater, hc::F, H, θwi, L) where {F}
+ θw, I_t, I_f, I_c, Lθ = FreezeCurves.freewater(H, θwi, L)
+ C = hc(θw)
+ T = (I_t*(H-Lθ) + I_f*H)/C
+ return T, θw, C
end
"""
freezethaw!(sub::SubSurface, heat::Heat{FreeWater,Enthalpy}, state)
diff --git a/src/Physics/HeatConduction/heat_bc.jl b/src/Physics/HeatConduction/heat_bc.jl
index 9690d4bf666c60239e9facad09c33d791794896f..af84a0a663bde6de247490c4a4c2a2a429cdb4ce 100644
--- a/src/Physics/HeatConduction/heat_bc.jl
+++ b/src/Physics/HeatConduction/heat_bc.jl
@@ -24,9 +24,10 @@ function CryoGrid.diagnosticstep!(::Top, bc::TemperatureGradient, state)
end
Base.@kwdef struct NFactor{W,S} <: BoundaryEffect
- nf::W = Param(1.0, domain=0..1) # applied when Tair <= 0
- nt::S = Param(1.0, domain=0..1) # applied when Tair > 0
+ nf::W = 1.0 # applied when Tair <= 0
+ nt::S = 1.0 # applied when Tair > 0
end
+InputOutput.parameterize(nf::NFactor; fields...) = NFactor(Param(nf.nf; domain=0..1, fields...), Param(nf.nt; domain=0..1, fields...))
CryoGrid.variables(::Top, bc::TemperatureGradient{<:NFactor}) = (
Diagnostic(:T_ub, Scalar, u"K"),
Diagnostic(:nfactor, Scalar),
diff --git a/src/Physics/Hydrology/Hydrology.jl b/src/Physics/Hydrology/Hydrology.jl
index badff06f60125b0924b9a7c278ae85367c0ef4f3..76c41f60eccedeab9820fa89217d79e1e3fc79ff 100644
--- a/src/Physics/Hydrology/Hydrology.jl
+++ b/src/Physics/Hydrology/Hydrology.jl
@@ -5,7 +5,7 @@ import CryoGrid
using CryoGrid
using CryoGrid.Physics
using CryoGrid.Numerics
-using CryoGrid.Numerics: nonlineardiffusion!
+using CryoGrid.Numerics: flux!, divergence!
using IfElse
using ModelParameters
diff --git a/src/Physics/Physics.jl b/src/Physics/Physics.jl
index d6040936da429298bd2c97e80f44e1cd1e81b4ff..177c89580844edc097584c8fd2a0f6880d2730f3 100644
--- a/src/Physics/Physics.jl
+++ b/src/Physics/Physics.jl
@@ -1,11 +1,11 @@
module Physics
+import CryoGrid
+
using CryoGrid
using CryoGrid.Numerics
using CryoGrid.Utils
-import CryoGrid
-
using Reexport
using Unitful
diff --git a/src/Physics/Snow/Snow.jl b/src/Physics/Snow/Snow.jl
index 94afde13508c3f616e388ce4b59029da01b97b65..9c6aa4fc63caf163263b2370dd5578ea0e4d9fa0 100644
--- a/src/Physics/Snow/Snow.jl
+++ b/src/Physics/Snow/Snow.jl
@@ -8,6 +8,7 @@ using CryoGrid.Numerics
using CryoGrid.Utils
import CryoGrid
+import CryoGrid.InputOutput
import CryoGrid.Physics
import CryoGrid.Physics.HeatConduction
@@ -20,8 +21,8 @@ export Snowpack, SnowProperties, SnowMassBalance, Snowfall
SnowProperties(
consts=Physics.Constants();
Ïw = consts.Ïw,
- Ïsn_new = Param(250.0, units=u"kg/m^3"),
- Ïsn_old = Param(500.0, units=u"kg/m^3"),
+ Ïsn_new = 250.0u"kg/m^3",
+ Ïsn_old = 500.0u"kg/m^3",
) = (; Ïw, Ïsn_new, Ïsn_old)
"""
@@ -50,8 +51,9 @@ end
abstract type SnowAccumulationScheme end
Base.@kwdef struct LinearAccumulation{S} <: SnowAccumulationScheme
- rate_scale::S = Param(1.0, domain=0..Inf) # scaling factor for snowfall rate
+ rate_scale::S = 1.0 # scaling factor for snowfall rate
end
+InputOutput.parameterize(acc::LinearAccumulation{<:Real}; fields...) = LinearAccumulation(Param(acc.rate_sacle; domain=0..Inf, fields...))
abstract type SnowDensityScheme end
# constant density (using Snowpack properties)
diff --git a/src/Physics/Snow/snow_bulk.jl b/src/Physics/Snow/snow_bulk.jl
index 83e246af103beb88f80cf57287664556794752fe..88caac402a30c06062830c3ee96cbfda6b079596 100644
--- a/src/Physics/Snow/snow_bulk.jl
+++ b/src/Physics/Snow/snow_bulk.jl
@@ -75,14 +75,12 @@ end
# heat upper boundary (for all bulk implementations)
CryoGrid.interact!(top::Top, bc::HeatBC, snow::BulkSnowpack, heat::Heat, stop, ssnow) = CryoGrid.interact!(CryoGrid.BoundaryStyle(bc), top, bc, snow, heat, stop, ssnow)
function CryoGrid.interact!(::CryoGrid.Dirichlet, top::Top, bc::HeatBC, snow::BulkSnowpack, heat::Heat, stop, ssnow)
- Δk = CryoGrid.thickness(snow, ssnow) # using `thickness` allows for generic layer implementations
@setscalar ssnow.T_ub = CryoGrid.boundaryvalue(bc, top, heat, snow, stop, ssnow)
if getscalar(ssnow.dsn) < threshold(snow)
@setscalar ssnow.T = getscalar(ssnow.T_ub)
end
# boundary flux
- @setscalar ssnow.dH_upper = CryoGrid.boundaryflux(bc, top, heat, snow, stop, ssnow)
- @inbounds ssnow.dH[1] += getscalar(ssnow.dH_upper) / Δk[1]
+ ssnow.jH[1] += CryoGrid.boundaryflux(bc, top, heat, snow, stop, ssnow)
return nothing # ensure no allocation
end
@@ -101,6 +99,7 @@ function CryoGrid.diagnosticstep!(
) where {TAcc,TAbl<:DegreeDayMelt,TDen<:ConstantDensity}
smb, heat = procs
Ïsn = snow.prop.Ïsn_new
+ HeatConduction.resetfluxes!(snow, heat, state)
@setscalar state.θwi = θwi = Ïsn / snow.prop.Ïw
@setscalar state.Ïsn = Ïsn
dsn = getscalar(state.swe) / θwi
@@ -112,7 +111,7 @@ function CryoGrid.diagnosticstep!(
# but capping the liquid fraction according to the 'max_unfrozen' parameter.
max_unfrozen = ablation(smb).max_unfrozen
θwi_cap = θwi*max_unfrozen
- T, θw, C = HeatConduction.enthalpyinv(heat.freezecurve, f_hc, getscalar(state.H), heat.prop.L, θwi_cap)
+ T, θw, C = HeatConduction.enthalpyinv(heat.freezecurve, f_hc, getscalar(state.H), θwi_cap, heat.prop.L)
# do not allow temperature to exceed 0°C
@. state.T = min(T, zero(T))
@. state.θw = θw
@@ -120,6 +119,7 @@ function CryoGrid.diagnosticstep!(
# compute thermal conductivity
HeatConduction.thermalconductivity!(snow, heat, state)
@. state.k = state.kc
+ return nothing
end
# snowfall upper boundary
function CryoGrid.interact!(
@@ -147,15 +147,15 @@ function CryoGrid.prognosticstep!(
end
if getscalar(state.swe) > 0.0 && getscalar(state.T_ub) > 0.0
ddf = ablation(smb).factor # [m/K/s]
- dH_upper = getscalar(state.dH_upper) # [J/m^3]
+ jH_upper = state.jH[1] # [J/m^3]
T_ub = getscalar(state.T_ub) # upper boundary temperature
Tref = 0.0*unit(T_ub) # just in case T_ub has units
# calculate the melt rate per second via the degree day model
- dmelt = max(ddf*(T_ub-Tref), zero(dH_upper)) # [m/s]
+ dmelt = max(ddf*(T_ub-Tref), zero(eltype(state.dswe))) # [m/s]
@. state.dswe += -dmelt
- # remove upper heat flux from dH if dmelt > 0;
- # this is due to the energy being "consumed" to melt the snow
- @. state.dH += -dH_upper*(dmelt > zero(dmelt))
+ # set upper heat flux to zero if dmelt > 0;
+ # this is due to the energy being (theoretically) "consumed" to melt the snow
+ state.jH[1] *= 1 - (dmelt > zero(dmelt))
end
return nothing
end
@@ -212,10 +212,11 @@ function CryoGrid.diagnosticstep!(
)
smb, heat = procs
Ïw = snow.prop.Ïw
+ HeatConduction.resetfluxes!(snow, heat, state)
new_swe = swe(snow, smb, state)
new_Ïsn = snowdensity(snow, smb, state)
new_dsn = new_swe*Ïw/new_Ïsn
- Ïw = heat.prop.Ïw
+ Ïw = heat.prop.consts.Ïw
if new_dsn > threshold(snow)
# if new snow depth is above threshold, set state variables
@setscalar state.swe = new_swe
@@ -238,14 +239,17 @@ function CryoGrid.diagnosticstep!(
@setscalar state.kc = heat.prop.ka
end
end
-# override prognosticstep! for incompatible heat types to prevent incorrect usage;
-# this will actually result in an ambiguous dispatch error before this error is thrown.
-CryoGrid.prognosticstep!(snow::BulkSnowpack, heat::Heat, state) = error("prognosticstep! not implemented for $(typeof(heat)) on $(typeof(snow))")
# prognosticstep! for free water, enthalpy based Heat on snow layer
-function CryoGrid.prognosticstep!(snow::BulkSnowpack, ::Heat{FreeWater,Enthalpy}, state)
+function CryoGrid.prognosticstep!(snow::BulkSnowpack, ps::Coupled2{<:SnowMassBalance,<:Heat{FreeWater,Enthalpy}}, state)
+ smb, heat = ps
+ prognosticstep!(snow, smb, state)
dsn = getscalar(state.dsn)
if dsn < snow.para.thresh
- # set energy flux to zero if there is no snow
+ # set divergence to zero if there is no snow
@. state.dH = zero(eltype(state.H))
+ else
+ # otherwise call prognosticstep! for heat
+ prognosticstep!(snow, heat, state)
end
end
+``
\ No newline at end of file
diff --git a/src/Physics/Soils/Soils.jl b/src/Physics/Soils/Soils.jl
index 8388962b7011f6f3843c59b79abf22e6a4c96e47..32416176e8f4d7518ee08893819db52412684726 100644
--- a/src/Physics/Soils/Soils.jl
+++ b/src/Physics/Soils/Soils.jl
@@ -13,10 +13,13 @@ using IfElse
using Interpolations: Interpolations
using IntervalSets
using ForwardDiff
+using FreezeCurves
using ModelParameters
+using Setfield
using Unitful
import CryoGrid
+import CryoGrid.InputOutput
import CryoGrid.Physics
import CryoGrid.Physics.HeatConduction
@@ -25,6 +28,9 @@ import Flatten: @flattenable, flattenable
export Soil, SoilParameterization, CharacteristicFractions, SoilProfile
export soilparameters, soilcomponent, porosity, mineral, organic
+# from FreezeCurves
+export SFCC, DallAmico, DallAmicoSalt, Westermann, McKenzie, VanGenuchten
+
const Enthalpy = HeatConduction.Enthalpy
const Temperature = HeatConduction.Temperature
@@ -48,12 +54,19 @@ abstract type SoilParameterization end
Represents uniform composition of a soil volume in terms of fractions: excess ice, natural porosity, saturation, and organic solid fraction.
"""
Base.@kwdef struct CharacteristicFractions{P1,P2,P3,P4} <: SoilParameterization
- xic::P1 = Param(0.0, domain=0..1) # excess ice fraction
- por::P2 = Param(0.5, domain=0..1) # natural porosity
- sat::P3 = Param(1.0, domain=0..1) # saturation
- org::P4 = Param(0.0, domain=0..1) # organic fraction of solid; mineral fraction is 1-org
+ xic::P1 = 0.0 # excess ice fraction
+ por::P2 = 0.5 # natural porosity
+ sat::P3 = 1.0 # saturation
+ org::P4 = 0.0 # organic fraction of solid; mineral fraction is 1-org
+end
+function soilparameters(::Type{CharacteristicFractions}=CharacteristicFractions; xic, por, sat, org)
+ CharacteristicFractions(
+ xic=Param(xic, domain=0..1),
+ por=Param(por, domain=0..1),
+ sat=Param(sat, domain=0..1),
+ org=Param(org, domain=0..1)
+ )
end
-soilparameters(::Type{CharacteristicFractions}=CharacteristicFractions; xic, por, sat, org) = CharacteristicFractions(xic=Param(xic, domain=0..1), por=Param(por, domain=0..1), sat=Param(sat, domain=0..1), org=Param(org, domain=0..1))
# Type alias for CharacteristicFractions with all scalar/numeric constituents
const HomogeneousCharacteristicFractions = CharacteristicFractions{<:Number,<:Number,<:Number,<:Number}
SoilProfile(pairs::Pair{<:DistQuantity,<:SoilParameterization}...) = Profile(pairs...)
@@ -66,12 +79,12 @@ soilcomponent(::Val{:θo}, χ, ϕ, θ, ω) = (1-χ)*(1-ϕ)*ω
"""
Soil thermal properties.
"""
-SoilThermalProperties(;
- ko = Param(0.25, units=u"W/m/K"), # organic [Hillel(1982)]
- km = Param(3.8, units=u"W/m/K"), # mineral [Hillel(1982)]
- co = Param(2.5e6, units=u"J/K/m^3"), # heat capacity organic
- cm = Param(2.0e6, units=u"J/K/m^3"), # heat capacity mineral
-) = (; ko, km, co, cm)
+Base.@kwdef struct SoilThermalProperties{Tko,Tkm,Tco,Tcm}
+ ko::Tko = 0.25u"W/m/K" # organic [Hillel(1982)]
+ km::Tkm = 3.8u"W/m/K" # mineral [Hillel(1982)]
+ co::Tco = 2.5e6u"J/K/m^3" # heat capacity organic
+ cm::Tcm = 2.0e6u"J/K/m^3" # heat capacity mineral
+end
"""
Basic Soil layer.
"""
@@ -80,6 +93,7 @@ Basic Soil layer.
prop::TProp = SoilThermalProperties()
sp::TSp = nothing # user-defined specialization
end
+InputOutput.parameterize(soil::Soil; fields...) = Soil(para=InputOutput.parameterize(soil.para; fields...), prop=soil.prop, sp=soil.sp)
HeatConduction.thermalproperties(soil::Soil) = soil.prop
# SoilComposition trait impl
SoilComposition(soil::Soil) = SoilComposition(typeof(soil))
@@ -153,10 +167,6 @@ Defaults to using the scalar porosity defined on `soil`.
"""
@inline porosity(soil::Soil, state, i) = Utils.getscalar(porosity(soil, state), i)
-export SWRC, VanGenuchten
-include("sfcc/swrc.jl")
-export SFCC, DallAmico, Westermann, McKenzie, SFCCNewtonSolver, SFCCPreSolver
-include("sfcc/sfcc.jl")
include("soilheat.jl")
end
diff --git a/src/Physics/Soils/sfcc/sfcc.jl b/src/Physics/Soils/sfcc/sfcc.jl
deleted file mode 100644
index d191ebea338f3916cdddcc6871fcb4844ad7f7bb..0000000000000000000000000000000000000000
--- a/src/Physics/Soils/sfcc/sfcc.jl
+++ /dev/null
@@ -1,143 +0,0 @@
-"""
-Abstract representation of a soil freeze characteristic curve (SFCC) function.
-Subtypes should be callable structs that implement the freeze curve and contain
-any necessary additional constants or configuration options. User-specified parameters
-can either be supplied in the struct or declared as model parameters via the `variables`
-method.
-"""
-abstract type SFCCFunction end
-"""
-Abstract type for SFCC H <--> T solvers.
-"""
-abstract type SFCCSolver end
-"""
- SFCC{F,S} <: FreezeCurve
-
-Generic representation of the soil freeze characteristic curve. The shape and parameters
-of the curve are determined by the implementation of SFCCFunction `f`. Also requires
-an implementation of SFCCSolver which provides the solution to the non-linear mapping H <--> T.
-"""
-struct SFCC{F,S} <: FreezeCurve
- f::F # freeze curve function f: (T,...) -> θ
- solver::S # solver for H -> T or T -> H
- SFCC(f::F, s::S=SFCCPreSolver()) where {F<:SFCCFunction,S<:SFCCSolver} = new{F,S}(f,s)
-end
-
-# Join the declared state variables of the SFCC function and the solver
-CryoGrid.variables(soil::Soil, heat::Heat, sfcc::SFCC) = tuplejoin(CryoGrid.variables(soil, heat, sfcc.f), CryoGrid.variables(soil, heat, sfcc.solver))
-# Default SFCC initialization
-function CryoGrid.initialcondition!(soil::Soil, heat::Heat, sfcc::SFCC, state)
- HeatConduction.freezethaw!(soil, heat, state)
-end
-
-"""
- sfccargs(f::SFCCFunction, soil::Soil, heat::Heat, state)
-
-Retrieves a tuple of values corresponding to each parameter declared by SFCCFunction `f` given the
-Soil layer, Heat process, and model state. The order of parameters *must match* the argument order
-of the freeze curve function `f`.
-"""
-sfccargs(::SFCCFunction, ::Soil, ::Heat, state) = ()
-# Fallback implementation of variables for SFCCFunction
-CryoGrid.variables(::Soil, ::Heat, f::SFCCFunction) = ()
-CryoGrid.variables(::Soil, ::Heat, s::SFCCSolver) = ()
-"""
- DallAmico <: SFCCFunction
-
-Dall'Amico M, 2010. Coupled water and heat transfer in permafrost modeling. Ph.D. Thesis, University of Trento, pp. 43.
-"""
-Base.@kwdef struct DallAmico{T,Θ,G,Tvg<:VanGenuchten} <: SFCCFunction
- Tₘ::T = Param(0.0, units=u"°C")
- θres::Θ = Param(0.0, domain=0..1)
- g::G = 9.80665u"m/s^2" # acceleration due to gravity
- swrc::Tvg = VanGenuchten()
-end
-sfccargs(f::DallAmico, soil::Soil, heat::Heat, state) = (
- porosity(soil, state), # θ saturated = porosity
- waterice(soil, heat, state), # total water content
- heat.prop.Lf, # specific latent heat of fusion, L
- f.θres,
- f.Tₘ,
- f.swrc.α,
- f.swrc.n,
-)
-# pressure head at T
-@inline ψ(T,Tstar,ψ₀,Lf,g) = ψ₀ + Lf/(g*Tstar)*(T-Tstar)*heaviside(Tstar-T)
-@inline function (f::DallAmico)(T, θsat, θtot, Lf, θres=f.θres, Tₘ=f.Tₘ, α=f.swrc.α, n=f.swrc.n)
- let θsat = max(θtot, θsat),
- g = f.g,
- Tₘ = normalize_temperature(Tₘ),
- ψ₀ = f.swrc(inv, θtot, θres, θsat, α, n),
- Tstar = Tₘ + g*Tₘ/Lf*ψ₀,
- T = normalize_temperature(T),
- ψ = ψ(T, Tstar, ψ₀, Lf, g);
- f.swrc(ψ, θres, θsat, α, n)
- end
-end
-
-"""
- McKenzie <: SFCCFunction
-
-McKenzie JM, Voss CI, Siegel DI, 2007. Groundwater flow with energy transport and water-ice phase change:
- numerical simulations, benchmarks, and application to freezing in peat bogs. Advances in Water Resources,
- 30(4): 966–983. DOI: 10.1016/j.advwatres.2006.08.008.
-"""
-Base.@kwdef struct McKenzie{T,Θ,Γ} <: SFCCFunction
- Tₘ::T = Param(0.0, units=u"°C")
- θres::Θ = Param(0.0, domain=0..1)
- γ::Γ = Param(0.1, domain=0..Inf, units=u"K")
-end
-sfccargs(f::McKenzie, soil::Soil, heat::Heat, state) = (
- porosity(soil, state), # θ saturated = porosity
- waterice(soil, heat, state), # total water content
- f.θres,
- f.Tₘ,
- f.γ,
-)
-function (f::McKenzie)(T, θsat, θtot, θres=f.θres, Tₘ=f.Tₘ, γ=f.γ)
- let T = normalize_temperature(T),
- Tₘ = normalize_temperature(Tₘ),
- θsat = max(θtot, θsat);
- IfElse.ifelse(T<=Tₘ, θres + (θsat-θres)*exp(-((T-Tₘ)/γ)^2), θtot)
- end
-end
-
-"""
- Westermann <: SFCCFunction
-
-Westermann, S., Boike, J., Langer, M., Schuler, T. V., and Etzelmüller, B.: Modeling the impact of
- wintertime rain events on the thermal regime of permafrost, The Cryosphere, 5, 945–959,
- https://doi.org/10.5194/tc-5-945-2011, 2011.
-"""
-Base.@kwdef struct Westermann{T,Θ,Δ} <: SFCCFunction
- Tₘ::T = Param(0.0, units=u"°C")
- θres::Θ = Param(0.0, domain=0..1)
- δ::Δ = Param(0.1, domain=0..Inf, units=u"K")
-end
-sfccargs(f::Westermann, soil::Soil, heat::Heat, state) = (
- porosity(soil, state), # θ saturated = porosity
- waterice(soil, heat, state), # total water content
- f.θres,
- f.Tₘ,
- f.δ,
-)
-function (f::Westermann)(T,θsat,θtot,θres=f.θres,Tₘ=f.Tₘ,δ=f.δ)
- let T = normalize_temperature(T),
- Tₘ = normalize_temperature(Tₘ),
- θsat = max(θtot, θsat);
- IfElse.ifelse(T<=Tₘ, θres - (θsat-θres)*(δ/(T-Tₘ-δ)), θtot)
- end
-end
-struct SFCCTable{F,I} <: SFCCFunction
- f::F
- f_tab::I
-end
-(f::SFCCTable)(args...) = f.f_tab(args...)
-"""
- Tabulated(f::SFCCFunction, args...)
-
-Produces an `SFCCTable` function which is a tabulation of `f`.
-"""
-Numerics.Tabulated(f::SFCCFunction, args...; kwargs...) = SFCCTable(f, Numerics.tabulate(f, args...; kwargs...))
-
-include("sfcc_solvers.jl")
diff --git a/src/Physics/Soils/sfcc/sfcc_solvers.jl b/src/Physics/Soils/sfcc/sfcc_solvers.jl
deleted file mode 100644
index 178c6f5ad4f159751eee51a5aa5d31b006bade49..0000000000000000000000000000000000000000
--- a/src/Physics/Soils/sfcc/sfcc_solvers.jl
+++ /dev/null
@@ -1,215 +0,0 @@
-using NLsolve
-
-"""
- SFCCNewtonSolver <: SFCCSolver
-
-Specialized implementation of Newton's method with backtracking line search for resolving
-the energy conservation law, H = TC + Lθ. Attempts to find the root of the corresponding
-temperature residual: ϵ = T - (H - Lθ(T)) / C(θ(T)) and uses backtracking to avoid
-jumping over the solution. This prevents convergence issues that arise due to
-discontinuities and strong non-linearity in most common soil freeze curves.
-"""
-Base.@kwdef struct SFCCNewtonSolver <: SFCCSolver
- maxiter::Int = 100 # maximum number of iterations
- abstol::Float64 = 1e-2 # absolute tolerance for convergence
- reltol::Float64 = 1e-4 # relative tolerance for convergence
- α₀::Float64 = 1.0 # initial step size multiplier
- Ï„::Float64 = 0.7 # step size decay for backtracking
-end
-# Helper function for updating θw, C, and the residual.
-@inline function sfccresidual(soil::Soil, heat::Heat, f::F, f_args::Fargs, f_hc, T, H) where {F,Fargs}
- L = heat.prop.L
- θw = f(T, f_args...)
- C = f_hc(θw)
- Tres = T - (H - θw*L) / C
- return Tres, θw, C
-end
-# Newton solver implementation
-function sfccsolve(solver::SFCCNewtonSolver, soil::Soil, heat::Heat, f, f_args, f_hc, H, Tâ‚€::Nothing=nothing)
- Tâ‚€ = IfElse.ifelse(H < zero(H), H / f_hc(0.0), zero(H))
- return sfccsolve(solver, soil, heat, f, f_args, f_hc, H, Tâ‚€)
-end
-using ForwardDiff
-function sfccsolve(solver::SFCCNewtonSolver, soil::Soil, heat::Heat, f::F, f_args, f_hc, H, Tâ‚€) where {F}
- T = Tâ‚€
- L = heat.prop.L
- cw = heat.prop.cw # heat capacity of liquid water
- ci = heat.prop.ci # heat capacity of ice
- α₀ = solver.α₀
- Ï„ = solver.Ï„
- # compute initial residual
- Tres, θw, C = sfccresidual(soil, heat, f, f_args, f_hc, T, H)
- itercount = 0
- T_converged = false
- while abs(Tres) > solver.abstol && abs(Tres) / abs(T) > solver.reltol
- if itercount >= solver.maxiter
- return (;T, Tres, θw, itercount, T_converged)
- end
- # derivative of freeze curve
- ∂θ∂T = ForwardDiff.derivative(T -> f(T, f_args...), T)
- # derivative of residual by quotient rule;
- # note that this assumes heatcapacity to be a simple weighted average!
- # in the future, it might be a good idea to compute an automatic derivative
- # of heatcapacity in addition to the freeze curve function.
- ∂Tres∂T = 1.0 - ∂θ∂T*(-L*C - (H - θw*L)*(cw-ci))/C^2
- α = α₀ / ∂Tres∂T
- T̂ = T - α*Tres
- # do first residual check outside of loop;
- # this way, we don't decrease α unless we have to.
- T̂res, θw, C = sfccresidual(soil, heat, f, f_args, f_hc, T̂, H)
- inneritercount = 0
- # simple backtracking line search to avoid jumping over the solution
- while sign(TÌ‚res) != sign(Tres)
- if inneritercount > 100
- @warn "Backtracking failed; this should not happen. Current state: α=$α, T=$T, T̂=$T̂, residual $(T̂res), initial residual: $(Tres)"
- return (;T, Tres, θw, itercount, T_converged)
- end
- α = α*τ # decrease step size by τ
- T̂ = T - α*Tres # new guess for T
- T̂res, θw, C = sfccresidual(soil, heat, f, f_args, f_hc, T̂, H)
- inneritercount += 1
- end
- T = TÌ‚ # update T
- Tres = TÌ‚res # update residual
- itercount += 1
- end
- T_converged = true
- return (;T, Tres, θw, itercount, T_converged)
-end
-function HeatConduction.enthalpyinv(soil::Soil, heat::Heat{<:SFCC{F,SFCCNewtonSolver},Enthalpy}, state, i) where {F}
- sfcc = freezecurve(heat)
- f = sfcc.f
- # get f arguments; note that this does create some redundancy in the arguments
- # eventually passed to the `residual` function; this is less than ideal but
- # probably shouldn't incur too much of a performance hit, just a few extra stack pointers!
- f_args = sfccargs(f, soil, heat, state)
- solver = sfcc.solver
- @inbounds let Tâ‚€ = i > 1 ? state.T[i-1] : nothing,
- H = state.H[i] |> Utils.adstrip, # enthalpy
- f_hc = partial(heatcapacity, liquidwater, soil, heat, state, i),
- f_argsáµ¢ = Utils.selectat(i, Utils.adstrip, f_args);
- T, _, _, _, _ = sfccsolve(solver, soil, heat, f, f_argsáµ¢, f_hc, H, Tâ‚€)
- return T
- end
-end
-"""
- SFCCPreSolver{TCache} <: SFCCSolver
-
-A fast SFCC "solver" which pre-builds an interpolant for the freeze curve in terms of enthalpy, θ(H).
-Note that this solver is **only valid when all freeze curve parameters are held constant** and will
-produce incorrect results otherwise.
-"""
-@flattenable struct SFCCPreSolver{TCache} <: SFCCSolver
- cache::TCache | false
- Tmin::Float64 | false
- errtol::Float64 | false
- SFCCPreSolver(cache, Tmin, errtol) = new{typeof(cache)}(cache, Tmin, errtol)
- """
- SFCCPreSolver(Tmin=-60.0, errtol=1e-4)
-
- Constructs a new `SFCCPreSolver` with minimum temperature `Tmin` and integration step `dH`.
- Enthalpy values below `H(Tmin)` under the given freeze curve will be extrapolated with a
- constant/flat function. `errtol` determines the permitted local error in the interpolant.
- """
- function SFCCPreSolver(;Tmin=-60.0, errtol=1e-4)
- cache = SFCCPreSolverCache()
- new{typeof(cache)}(cache, Tmin, errtol)
- end
-end
-mutable struct SFCCPreSolverCache{TFH,T∇FH,T∇FT}
- f_H::TFH # θw(H) interpolant
- ∇f_H::T∇FH # derivative of θw(H) interpolant
- ∇f_T::T∇FT # ∂θ∂T interpolant
- function SFCCPreSolverCache()
- # initialize with dummy functions to get type information;
- # this is just to make sure that the compiler can still infer all of the types.
- x = -3e8:1e6:3e8
- dummy_f = _build_interpolant(x, zeros(length(x)))
- dummy_∇f = first ∘ _interpgrad(dummy_f)
- return new{typeof(dummy_f),typeof(dummy_∇f),typeof(dummy_f)}(dummy_f, dummy_∇f, dummy_f)
- end
-end
-_interpgrad(f) = (args...) -> Interpolations.gradient(f, args...)
-function _build_interpolant(Hs, θs)
- return Interpolations.extrapolate(
- Interpolations.interpolate((Vector(Hs),), θs, Interpolations.Gridded(Interpolations.Linear())),
- Interpolations.Flat()
- )
-end
-function CryoGrid.initialcondition!(soil::Soil{<:HomogeneousCharacteristicFractions}, heat::Heat, sfcc::SFCC{F,<:SFCCPreSolver}, state) where {F}
- L = heat.prop.L
- args = sfccargs(sfcc.f, soil, heat, state)
- state.θw .= sfcc.f.(state.T, args...)
- heatcapacity!(soil, heat, state)
- @. state.H = enthalpy(state.T, state.C, L, state.θw)
- # pre-solve freeze curve;
- # note that this is only valid given that the following assumptions hold:
- # 1) none of the freeze curve parameters (e.g. soil properties) change
- # 2) soil properties are uniform in `soil`
- let θsat = args[1],
- θtot = args[2],
- tail_args = args[3:end],
- θm = mineral(soil, state),
- θo = organic(soil, state),
- θa = 1-θm-θo-θtot,
- L = heat.prop.L,
- Tmin = sfcc.solver.Tmin,
- Tmax = 0.0,
- f(T) = sfcc.f(T, θsat, θtot, tail_args...),
- C(θw) = heatcapacity(soil, heat, θw, θtot-θw, θa, θm, θo),
- Hmin = enthalpy(Tmin, C(f(Tmin)), L, f(Tmin)),
- Hmax = enthalpy(Tmax, C(f(Tmax)), L, f(Tmax));
- # residual as a function of T and H
- resid(T,H) = sfccresidual(soil, heat, sfcc.f, args, C, T, H)
- function solve(H,Tâ‚€)
- local T = nlsolve(T -> resid(T[1],H)[1], [Tâ‚€]).zero[1]
- θw = f(T)
- return (; T, θw)
- end
- function deriv(T) # implicit partial derivative w.r.t H as a function of T
- θw, ∂θ∂T = ∇(f, T)
- # get C_eff, i.e. dHdT
- ∂H∂T = HeatConduction.C_eff(T, C(θw), heat.prop.L, ∂θ∂T, heat.prop.cw, heat.prop.ci)
- # valid by chain rule and inverse function theorem
- return ∂θ∂T / ∂H∂T, ∂θ∂T
- end
- function step(ΔH, H, θ, ∂θ∂H, T₀)
- # get first order estimate
- θest = θ + ΔH*∂θ∂H
- # get true θ at H + ΔH
- θsol = solve(H + ΔH, T₀).θw
- err = abs(θsol - θest)
- # return residual of error with target error
- return err
- end
- T = [Tmin]
- H = [Hmin]
- θ = [f(T[1])]
- dθdT₀, dθdH₀ = deriv(T[1])
- dθdT = [dθdT₀]
- dθdH = [dθdH₀]
- @assert isfinite(H[1]) && isfinite(θ[1]) "H=$H, θ=$θ"
- while H[end] < Hmax
- # find the optimal step size
- ϵ = Inf
- ΔH = heat.prop.L*θtot/10 # initially set to large value
- while abs(ϵ) > sfcc.solver.errtol
- ϵ = step(ΔH, H[end], θ[end], dθdH[end], T[end])
- # iteratively halve the step size until error tolerance is satisfied
- ΔH *= 0.5
- end
- Hnew = H[end] + ΔH
- @assert isfinite(Hnew) "isfinite(ΔH) failed; H=$(H[end]), T=$(T[end]), ΔH=$ΔH"
- opt = solve(Hnew, T[end])
- dθdHᵢ, dθdTᵢ = deriv(opt.T)
- push!(H, Hnew)
- push!(θ, opt.θw)
- push!(T, opt.T)
- push!(dθdT, dθdTᵢ)
- push!(dθdH, dθdHᵢ)
- end
- sfcc.solver.cache.f_H = _build_interpolant(H, θ)
- sfcc.solver.cache.∇f_H = first ∘ _interpgrad(sfcc.solver.cache.f_H)
- sfcc.solver.cache.∇f_T = _build_interpolant(T, dθdT)
- end
-end
diff --git a/src/Physics/Soils/sfcc/swrc.jl b/src/Physics/Soils/sfcc/swrc.jl
deleted file mode 100644
index 21d89e6976f57c3fd262c94d2cef1862857004e9..0000000000000000000000000000000000000000
--- a/src/Physics/Soils/sfcc/swrc.jl
+++ /dev/null
@@ -1,36 +0,0 @@
-"""
- SWRCFunction
-
-Base type for soil water retention curve SWRC function implementations.
-"""
-abstract type SWRCFunction end
-Base.inv(f::SWRCFunction) = (args...) -> f(inv, args...)
-"""
- SWRC{F}
-
-Soil water retention curve with function type `F`.
-"""
-struct SWRC{F}
- f::F # soil water retention curve function f: (ψ,...) -> θ
- SWRC(f::F) where {F<:SWRCFunction} = new{F}(f)
-end
-"""
- VanGenuchten <: SWRCFunction
-
-van Genuchten MT, 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils.
- Soil Science Society of America Journal, 44(5): 892–898. DOI: 10.2136/sssaj 1980.03615995004400050002x.
-"""
-Base.@kwdef struct VanGenuchten{Tα,Tn} <: SWRCFunction
- α::Tα = Param(1.0, units=u"1/m")
- n::Tn = Param(2.0)
-end
-function (f::VanGenuchten)(ψ, θres, θsat, α=f.α, n=f.n)
- let m = 1-1/n;
- IfElse.ifelse(ψ <= zero(ψ), θres + (θsat - θres)*(1 + abs(-α*ψ)^n)^(-m), θsat)
- end
-end
-function (f::VanGenuchten)(::typeof(inv), θ, θres, θsat, α=f.α, n=f.n)
- let m = 1-1/n;
- IfElse.ifelse(θ < θsat, -1/α*(((θ-θres)/(θsat-θres))^(-1/m)-1.0)^(1/n), zero(1/α))
- end
-end
diff --git a/src/Physics/Soils/soilheat.jl b/src/Physics/Soils/soilheat.jl
index 4d730f985da3c814934a630b138abe70353d648b..2e9b30eea4ea425c1c2590b929154c13ac4a6424 100644
--- a/src/Physics/Soils/soilheat.jl
+++ b/src/Physics/Soils/soilheat.jl
@@ -17,11 +17,44 @@
(θw, θi, θa, θm, θo)
end
end
+# Default implementations of variables for SFCCFunction
+CryoGrid.variables(::Soil, ::Heat, f::SFCCFunction) = ()
+CryoGrid.variables(::Soil, ::Heat, s::SFCCSolver) = ()
+# Join the declared state variables of the SFCC function and the solver
+CryoGrid.variables(soil::Soil, heat::Heat, sfcc::SFCC) = tuplejoin(CryoGrid.variables(soil, heat, sfcc.f), CryoGrid.variables(soil, heat, sfcc.solver))
+"""
+ sfcckwargs(f::SFCCFunction, soil::Soil, heat::Heat, state, i)
+Builds a named tuple of values corresponding to each keyword arguments of the SFCCFunction `f`
+which should be set according to the layer/process properties or state. The default implementation
+sets only the total water content, θtot = θwi, and the saturated water content, θsat = θp.
+"""
+sfcckwargs(::SFCCFunction, soil::Soil, heat::Heat, state, i) = (
+ θtot = waterice(soil, heat, state, i), # total water content
+ θsat = porosity(soil, state, i), # θ saturated = porosity
+)
"""
Initial condition for heat conduction (all state configurations) on soil layer w/ SFCC.
"""
-CryoGrid.initialcondition!(soil::Soil, heat::Heat{<:SFCC}, state) = HeatConduction.initialcondition!(soil, heat, freezecurve(heat), state)
+function CryoGrid.initialcondition!(soil::Soil, heat::Heat{<:SFCC}, state)
+ fc = freezecurve(heat)
+ L = heat.prop.L
+ @inbounds for i in 1:length(state.T)
+ fc_kwargsáµ¢ = sfcckwargs(fc.f, soil, heat, state, i)
+ hc = partial(heatcapacity, liquidwater, soil, heat, state, i)
+ if i == 1
+ # TODO: this is currently only relevant for the pre-solver scheme and assumes that
+ # the total water content is uniform throughout the layer and does not change over time.
+ FreezeCurves.Solvers.initialize!(fc.solver, fc.f, hc; fc_kwargsáµ¢...)
+ end
+ T = state.T[i]
+ θw, dθdT = ∇(T -> fc(T; fc_kwargsᵢ...), T)
+ state.θw[i] = θw
+ state.C[i] = hc(θw)
+ state.H[i] = enthalpy(state.T[i], state.C[i], L, state.θw[i])
+ state.dHdT[i] = HeatConduction.C_eff(T, state.C[i], L, dθdT, heat.prop.cw, heat.prop.ci)
+ end
+end
"""
Initial condition for heat conduction (all state configurations) on soil layer w/ free water freeze curve.
"""
@@ -46,12 +79,11 @@ For heat conduction with temperature, we can simply evaluate the freeze curve to
function HeatConduction.freezethaw!(soil::Soil, heat::Heat{<:SFCC,Temperature}, state)
sfcc = freezecurve(heat)
let L = heat.prop.L,
- f = sfcc.f,
- f_args = sfccargs(f,soil,heat,state);
+ f = sfcc.f;
@inbounds @fastmath for i in 1:length(state.T)
T = state.T[i]
- f_argsáµ¢ = Utils.selectat(i, identity, f_args)
- θw, dθdT = ∇(T -> f(T, f_argsᵢ...), T)
+ f_argsáµ¢ = sfcckwargs(f, soil, heat, state, i)
+ θw, dθdT = ∇(T -> f(T; f_argsᵢ...), T)
state.θw[i] = θw
state.dθdT[i] = dθdT
state.C[i] = C = heatcapacity(soil, heat, volumetricfractions(soil, heat, state, i)...)
@@ -60,42 +92,40 @@ function HeatConduction.freezethaw!(soil::Soil, heat::Heat{<:SFCC,Temperature},
end
end
end
-function HeatConduction.freezethaw!(soil::Soil, heat::Heat{<:SFCC{F,SFCCNewtonSolver},Enthalpy}, state) where {F}
+function HeatConduction.freezethaw!(soil::Soil, heat::Heat{<:SFCC,Enthalpy}, state)
sfcc = freezecurve(heat)
- f_args = sfccargs(sfcc.f, soil, heat, state)
@inbounds for i in 1:length(state.H)
- let f_argsáµ¢ = Utils.selectat(i, identity, f_args),
- f = sfcc.f;
- # Evaluate inverse enthalpy function
- T = enthalpyinv(soil, heat, state, i)
- # Here we apply the recovered temperature to the state variables;
- # Since we perform iteration on untracked variables, we need to
- # recompute θw, C, and T here with the tracked variables.
- # Note that this results in one additional freeze curve function evaluation.
- state.θw[i], dθdT = ∇(T -> f(T, f_argsᵢ...), T)
- let H = state.H[i],
- L = heat.prop.L,
- cw = heat.prop.cw,
- ci = heat.prop.ci,
- θw = state.θw[i],
- (_, θi, θa, θm, θo) = volumetricfractions(soil, heat, state, i);
- state.C[i] = heatcapacity(soil, heat, θw, θi, θa, θm, θo)
- state.T[i] = (H - L*θw) / state.C[i]
- state.dHdT[i] = HeatConduction.C_eff(state.T[i], state.C[i], L, dθdT, cw, ci)
- end
+ @inbounds let H = state.H[i], # enthalpy
+ L = heat.prop.L,
+ cw = heat.prop.cw,
+ ci = heat.prop.ci,
+ θwi = waterice(soil, heat, state, i), # total water content
+ T₀ = i > 1 ? state.T[i-1] : FreezeCurves.freewater(H, θwi, L), # initial guess for T
+ hc = partial(heatcapacity, liquidwater, soil, heat, state, i),
+ f = sfcc.f,
+ f_kwargsáµ¢ = sfcckwargs(f, soil, heat, state, i),
+ obj = FreezeCurves.SFCCInverseEnthalpyObjective(f, f_kwargsáµ¢, hc, L, H);
+ res = FreezeCurves.sfccsolve(obj, sfcc.solver, Tâ‚€, Val{true}())
+ state.T[i] = res.T
+ state.θw[i] = res.θw
+ state.C[i] = res.C
+ state.dHdT[i] = HeatConduction.C_eff(state.T[i], state.C[i], L, res.dθwdT, cw, ci)
end
end
end
-function HeatConduction.freezethaw!(soil::Soil{<:HomogeneousCharacteristicFractions}, heat::Heat{<:SFCC{F,<:SFCCPreSolver},Enthalpy}, state) where {F}
- solver = freezecurve(heat).solver
- f_H = solver.cache.f_H
- ∇f_T = solver.cache.∇f_T
- L, cw, ci = heat.prop.L, heat.prop.cw, heat.prop.ci
- @inbounds for i in 1:length(state.H)
- H = state.H[i]
- state.θw[i] = θw = f_H(H)
- state.C[i] = C = heatcapacity(soil, heat, volumetricfractions(soil, heat, state, i)...)
- state.T[i] = T = enthalpyinv(H, C, heat.prop.L, θw)
- state.dHdT[i] = HeatConduction.C_eff(T, C, L, ∇f_T(T), cw, ci)
+function HeatConduction.enthalpyinv(soil::Soil, heat::Heat{<:SFCC,Enthalpy}, state, i)
+ sfcc = freezecurve(heat)
+ @inbounds let H = state.H[i], # enthalpy
+ L = heat.prop.L, # latent heat of fusion of water
+ θwi = waterice(soil, heat, state, i), # total water content
+ T₀ = i > 1 ? state.T[i-1] : freewater(H, θwi, L),
+ hc = partial(heatcapacity, liquidwater, soil, heat, state, i),
+ f = sfcc.f,
+ f_kwargsáµ¢ = sfcckwargs(f, soil, heat, state, i),
+ obj = FreezeCurves.SFCCInverseEnthalpyObjective(f, f_kwargsáµ¢, hc, L, H);
+ T_sol = FreezeCurves.sfccsolve(obj, sfcc.solver, Tâ‚€, Val{false}())
+ return T_sol
end
end
+# Parameterization
+InputOutput.parameterize(solver::FreezeCurves.SFCCSolver; fields...) = solver
diff --git a/src/Physics/Sources/Sources.jl b/src/Physics/Sources/Sources.jl
index 6e60ed49416a92f779af426b69b89505350e091c..a9fc83e7059a3f422698ffcb2f1863dde06f4e9e 100644
--- a/src/Physics/Sources/Sources.jl
+++ b/src/Physics/Sources/Sources.jl
@@ -23,7 +23,7 @@ abstract type SourceTerm end
Parametric source term that is constant through time and space.
"""
Base.@kwdef struct Constant{S} <: SourceTerm
- Sâ‚€::S = Param(0.0)
+ Sâ‚€::S = 0.0
end
"""
Periodic <: SourceTerm
@@ -31,10 +31,10 @@ end
Parametric source term that is periodic through time and constant through space.
"""
Base.@kwdef struct Periodic{A,F,S,L} <: SourceTerm
- amp::A = Param(1.0)
- freq::F = Param(1.0/(3600*24))
- shift::S = Param(0.0)
- level::L = Param(0.0)
+ amp::A = 1.0
+ freq::F = 1.0/(3600*24)
+ shift::S = 0.0
+ level::L = 0.0
end
"""
Source{P,T,S} <: SubSurfaceProcess
diff --git a/src/Physics/common.jl b/src/Physics/common.jl
index 2d09ead7ee0e66f22844d1f407b9a8f233fc471a..4f15637dea071bb0337a9c5f3e2e031b00093231 100644
--- a/src/Physics/common.jl
+++ b/src/Physics/common.jl
@@ -1,9 +1,14 @@
Constants() = (
Ïw = 1000.0u"kg/m^3", # density of water at standard conditions
- Lf = 3.34e5u"J/kg", # specific latent heat of fusion of water [J/kg]
+ Lsl = 3.34e5u"J/kg", # specific latent heat of fusion of water [J/kg]
g = 9.80665u"m/s^2", # gravitational constant
)
-# Composition
+# Generic step limiter types
+abstract type StepLimiter end
+Base.@kwdef struct CFL <: StepLimiter
+ fallback_dt::Float64 = 60.0 # fallback dt [s]
+end
+# Volume material composition
"""
volumetricfractions(::SubSurface, ::SubSurfaceProcess, state)
volumetricfractions(::SubSurface, ::SubSurfaceProcess, state, i)
diff --git a/src/Strat/Strat.jl b/src/Strat/Strat.jl
index 26edb14f8fb0aab9bedbf14b906991d6b5c5c201..3a3fcfb623d944d2a7ea9f6020469b11cb4969c1 100644
--- a/src/Strat/Strat.jl
+++ b/src/Strat/Strat.jl
@@ -1,6 +1,7 @@
module Strat
import CryoGrid
+import CryoGrid.InputOutput
using CryoGrid
using CryoGrid: Parameterization, DynamicParameterization
diff --git a/src/Strat/state.jl b/src/Strat/state.jl
index b3ba3fecaf1df5e12bb44458a9dda75d110cf753..7bcaf1606863a78af1f0d91f1a7a548b31852281 100644
--- a/src/Strat/state.jl
+++ b/src/Strat/state.jl
@@ -1,4 +1,4 @@
-using CryoGrid.Numerics: Delta
+using CryoGrid.Numerics: Delta, Flux
"""
Enumeration for in-place vs out-of-place mode.
"""
@@ -86,11 +86,12 @@ end
@inline _makegrid(::Var{name,<:OnGrid{Edges}}, vs::VarStates, z_inds) where {name} = vs.grid[z_inds]
@inline _makegrid(var::Var, vs::VarStates, z_inds) = 1:dimlength(vardims(var), vs.grid)
@inline _makestate(::Val, ::Prognostic{name,<:OnGrid{Cells}}, vs::VarStates, z_inds, u, du, t) where {name} = view(view(u, Val{name}()), infimum(z_inds):supremum(z_inds)-1)
-@inline _makestate(::Val, ::Prognostic{name}, vs::VarStates, z_inds, u, du, t) where {name} = view(u, Val{name}())
+@inline _makestate(::Val, ::Prognostic{name,<:Shape}, vs::VarStates, z_inds, u, du, t) where {name} = view(u, Val{name}())
@inline _makestate(::Val, ::Algebraic{name,<:OnGrid{Cells}}, vs::VarStates, z_inds, u, du, t) where {name} = view(view(u, Val{name}()), infimum(z_inds):supremum(z_inds)-1)
-@inline _makestate(::Val, ::Algebraic{name}, vs::VarStates, z_inds, u, du, t) where {name} = view(u, Val{name}())
+@inline _makestate(::Val, ::Algebraic{name,<:Shape}, vs::VarStates, z_inds, u, du, t) where {name} = view(u, Val{name}())
@inline _makestate(::Val, ::Delta{dname,name,<:OnGrid{Cells}}, vs::VarStates, z_inds, u, du, t) where {dname,name} = view(view(du, Val{name}()), infimum(z_inds):supremum(z_inds)-1)
-@inline _makestate(::Val, ::Delta{dname,name}, vs::VarStates, z_inds, u, du, t) where {dname,name} = view(du, Val{name}())
+@inline _makestate(::Val, ::Delta{dname,name,<:Shape}, vs::VarStates, z_inds, u, du, t) where {dname,name} = view(du, Val{name}())
+@inline _makestate(::Val, ::Flux{jname,name,<:OnGrid{Edges}}, vs::VarStates, z_inds, u, du, t) where {jname,name} = view(retrieve(vs.griddiag[jname], u, t), infimum(z_inds):supremum(z_inds))
@inline _makestate(::Val, ::Diagnostic{name,<:OnGrid{Cells}}, vs::VarStates, z_inds, u, du, t) where {name} = view(retrieve(vs.griddiag[name], u, t), infimum(z_inds):supremum(z_inds)-1)
@inline _makestate(::Val, ::Diagnostic{name,<:OnGrid{Edges}}, vs::VarStates, z_inds, u, du, t) where {name} = view(retrieve(vs.griddiag[name], u, t), infimum(z_inds):supremum(z_inds))
@inline _makestate(::Val{layername}, ::Diagnostic{name}, vs::VarStates, z_inds, u, du, t) where {name,layername} = retrieve(vs.diag[layername][name], u, t)
diff --git a/src/Strat/stratigraphy.jl b/src/Strat/stratigraphy.jl
index d35ce6c7e76c3a4da46d53e48555648aa1e1a0d4..97d1a219061f2065aeae53e79129442dd0c0cc09 100644
--- a/src/Strat/stratigraphy.jl
+++ b/src/Strat/stratigraphy.jl
@@ -1,4 +1,4 @@
-const RESERVED_COMPONENT_NAMES = (:top, :bottom, :strat, :init)
+const RESERVED_COMPONENT_NAMES = (:top, :bottom, :strat, :init, :event)
"""
StratComponent{TLayer,TProc,name}
@@ -97,6 +97,7 @@ boundarypairs(strat::Stratigraphy, z_bottom) = boundarypairs(boundaries(strat),
boundarypairs(bounds::NTuple, z_bottom) = tuplejoin(map(tuple, bounds[1:end-1], bounds[2:end]), ((bounds[end], z_bottom),))
componentnames(strat::Stratigraphy) = map(componentname, components(strat))
componenttypes(::Type{<:Stratigraphy{N,TComponents}}) where {N,TComponents} = Tuple(TComponents.parameters)
+InputOutput.parameterize(strat::Stratigraphy; fields...) = Stratigraphy(getfield(strat, :boundaries), map(comp -> InputOutput.parameterize(comp; layer=componentname(comp), fields...), getfield(strat, :components)))
Base.keys(strat::Stratigraphy) = componentnames(strat)
Base.values(strat::Stratigraphy) = components(strat)
@inline Base.propertynames(strat::Stratigraphy) = Base.keys(strat)
diff --git a/src/Strat/tile.jl b/src/Strat/tile.jl
index c1a5b6d542218ab648bb09a98c000e0f6fd07617..dedcf99a31696a420a9c64fcc83e83dd6c628c09 100644
--- a/src/Strat/tile.jl
+++ b/src/Strat/tile.jl
@@ -60,6 +60,14 @@ struct Tile{TStrat,TGrid,TStates,TInits,TEvents,iip,obsv} <: AbstractTile{iip}
end
ConstructionBase.constructorof(::Type{Tile{TStrat,TGrid,TStates,TInits,TEvents,iip,obsv}}) where {TStrat,TGrid,TStates,TInits,TEvents,iip,obsv} =
(strat, grid, state, inits, events, hist) -> Tile(strat, grid, state, inits, events, hist, iip, obsv)
+InputOutput.parameterize(tile::T) where {T<:Tile} = ConstructionBase.constructorof(T)(
+ InputOutput.parameterize(tile.strat),
+ tile.grid,
+ tile.state,
+ InputOutput.parameterize(tile.inits, layer=:init),
+ InputOutput.parameterize(tile.events, layer=:event),
+ tile.hist,
+)
# mark only stratigraphy and initializers fields as flattenable
Flatten.flattenable(::Type{<:Tile}, ::Type{Val{:strat}}) = true
Flatten.flattenable(::Type{<:Tile}, ::Type{Val{:inits}}) = true
@@ -263,10 +271,11 @@ CryoGrid.initialcondition!(tile::Tile, tspan::NTuple{2,DateTime}, p::AbstractVec
for i in 1:N
for j in 1:length(TInits.parameters)
@>> quote
- let layerstate = state[$i],
+ let layer = strat[$i].layer,
+ layerstate = state[$i],
init! = tile.inits[$j];
if haskey(layerstate.states, varname(init!))
- init!(layerstate)
+ init!(layer, layerstate)
end
end
end push!(expr.args)
diff --git a/src/Utils/Utils.jl b/src/Utils/Utils.jl
index 760346d431f94b535e4042c58c0f019e5c6e6ec8..af4295669a02ef636a38a7f0ed11bd1167bcc34e 100644
--- a/src/Utils/Utils.jl
+++ b/src/Utils/Utils.jl
@@ -14,15 +14,17 @@ using Unitful
import CryoGrid
import ForwardDiff
+import FreezeCurves
import Unitful
+import FreezeCurves: normalize_temperature
import ModelParameters: stripunits
export @xu_str, @Float_str, @Real_str, @Number_str, @UFloat_str, @UT_str, @setscalar, @threaded, @sym_str, @pstrip
include("macros.jl")
export DistUnit, DistQuantity, TempUnit, TempQuantity, TimeUnit, TimeQuantity
-export dustrip, duconvert, applyunit, normalize_temperature, pstrip
+export dustrip, duconvert, applyunits, normalize_units, normalize_temperature, pstrip
export structiterate, getscalar, tuplejoin, convert_t, convert_tspan, haskeys
export IterableStruct
@@ -40,12 +42,12 @@ StructTypes.lowertype(value::Type{<:Quantity}) = String
StructTypes.construct(::Type{Q}, value::String) where {Q<:Quantity} = uconvert(Q, uparse(replace(value, " " => "")))
"""
- applyunit(u::Unitful.Units, x::Number)
+ applyunits(u::Unitful.Units, x::Number)
Conditionally applies unit `u` to `x` if and only if `x` is a unit-free quantity.
If `x` is a unitful quantity, asserts that the unit matches `u`.
"""
-function applyunit(u::Unitful.Units, x::Number)
+function applyunits(u::Unitful.Units, x::Number)
if typeof(unit(x)) <: Unitful.FreeUnits{(), NoDims, nothing}
return x*u
else
@@ -54,14 +56,11 @@ function applyunit(u::Unitful.Units, x::Number)
end
end
-"""
- normalize_temperature(x)
-
-Converts temperature `x` to Kelvin. If `x` has units, `uconvert` is used. Otherwise, if `x` a general numeric type, it is assumed that `x` is in celsius.
-"""
-normalize_temperature(x) = x + 273.15
-normalize_temperature(x::TempQuantity) = uconvert(u"K", x)
-normalize_temperature(x::Param) = stripparams(x) |> normalize_temperature
+# special case: make sure temperatures are in °C
+normalize_units(x::Unitful.AbstractQuantity{T,Unitful.ðš¯}) where T = uconvert(u"°C", x)
+normalize_units(x::Unitful.AbstractQuantity) = upreferred(x)
+# Add method dispatch for normalize_temperature in FreezeCurves.jl
+normalize_temperature(x::Param) = normalize_temperature(stripparams(x))
"""
Provides implementation of `Base.iterate` for structs.
@@ -73,7 +72,6 @@ function structiterate(obj::A) where {A}
(val,state) = iterate(gen)
(val, (gen,state))
end
-
function structiterate(obj, state)
gen, genstate = state
nextitr = iterate(gen,genstate)
@@ -134,16 +132,6 @@ Convenience method for converting between `Dates.DateTime` and solver time.
convert_tspan(tspan::NTuple{2,DateTime}) = convert_t.(tspan)
convert_tspan(tspan::NTuple{2,Float64}) = convert_t.(tspan)
-"""
- @generated selectat(i::Int, f, args::T) where {T<:Tuple}
-
-Helper function for handling mixed vector/scalar arguments to runtime generated functions.
-select calls getindex(i) for all array-typed arguments leaves non-array arguments as-is.
-We use a generated function to expand the arguments into an explicitly defined tuple to preserve type-stability (i.e. it's an optmization);
-function `f` is then applied to each element.
-"""
-@generated selectat(i::Int, f, args::T) where {T<:Tuple} = :(tuple($([typ <: AbstractArray ? :(f(args[$k][i])) : :(f(args[$k])) for (k,typ) in enumerate(Tuple(T.parameters))]...)))
-
"""
ffill!(x::AbstractVector{T}) where {E,T<:Union{Missing,E}}
@@ -158,19 +146,6 @@ function ffill!(x::AbstractVector{T}) where {E,T<:Union{Missing,E}}
return x
end
-"""
- adstrip(x::Number)
- adstrip(x::ForwardDiff.Dual)
- adstrip(x::ReverseDiff.TrackedReal)
-
-adstrip extracts the underlying numeric value from `x` if `x` is a tracked value from
-an autodiff library (e.g. ForwardDiff or ReverseDiff). If `x` is not an AD type, then
-`adstrip` simply returns `x` as-is.
-"""
-adstrip(x::Number) = x
-adstrip(x::ForwardDiff.Dual) = ForwardDiff.value(x) |> adstrip
-adstrip(x::Param{T}) where {T<:ForwardDiff.Dual} = Param(NamedTuple{Tuple(keys(x))}((adstrip(x.val), Base.tail(parent(x))...)))
-
"""
Debug ustrip. Remove units if and only if debug mode is NOT enabled.
"""
@@ -178,7 +153,9 @@ dustrip(x::Number) = CryoGrid.CRYOGRID_DEBUG ? x : ustrip(x)
dustrip(x::AbstractVector{<:Quantity{T}}) where {T} = CryoGrid.CRYOGRID_DEBUG ? x : reinterpret(T, x)
dustrip(u::Unitful.Units, x::Number) = CryoGrid.CRYOGRID_DEBUG ? x : ustrip(u,x)
dustrip(u::Unitful.Units, x::AbstractVector{<:Quantity{T}}) where {T} = CryoGrid.CRYOGRID_DEBUG ? x : reinterpret(T, uconvert.(u, x))
-
+"""
+Debug uconvert.
+"""
duconvert(u::Unitful.Units, x::Number) = CryoGrid.CRYOGRID_DEBUG ? x : uconvert(u, x)
"""
@@ -205,9 +182,6 @@ Additional override for `stripunits` which reconstructs `obj` with all fields th
types converted to base SI units and then stripped to be unit free.
"""
function ModelParameters.stripunits(obj)
- # special case: make sure temperatures are in °C
- normalize_units(x::Unitful.AbstractQuantity{T,Unitful.ðš¯}) where T = uconvert(u"°C", x)
- normalize_units(x::Unitful.AbstractQuantity) = upreferred(x)
values = Flatten.flatten(obj, Flatten.flattenable, Unitful.AbstractQuantity, Flatten.IGNORE)
return Flatten.reconstruct(obj, map(ustrip ∘ normalize_units, values), Unitful.AbstractQuantity, Flatten.IGNORE)
end
diff --git a/test/Numerics/math_tests.jl b/test/Numerics/math_tests.jl
index 1fea8ad43c1264b3696a168668f5e94d9d108c7a..288c33ca0c31f97eeca23763f22b7ef0e55379f8 100644
--- a/test/Numerics/math_tests.jl
+++ b/test/Numerics/math_tests.jl
@@ -1,36 +1,27 @@
-using CryoGrid.Numerics: finitediff!, lineardiffusion!, nonlineardiffusion!
+using CryoGrid.Numerics
+using CryoGrid.Numerics: flux!, divergence!, nonlineardiffusion!
using Test
include("../testutils.jl")
-@testset "finitediff!" begin
- f(x) = (1/2)x^2 # function to differntiate
- df(x) = x # analytical 2nd derivative
- x = exp.(0.0:0.001:0.05)
- y = f.(x)
- dy = df.(x[1:end-1])
- δ = x[2:end] .- x[1:end-1]
- ∂y = zeros(length(y)-1)
- finitediff!(∂y,y,δ)
- @test allfinite(∂y)
- @test allequal(∂y,dy,atol=0.01)
-end
-
-@testset "lineardiffusion!" begin
- f(x) = (1/6)x^3 # function to differntiate
- d2f(x) = x # analytical 2nd derivative
- x = exp.(0.0:0.001:0.05)
- y = f.(x)
- d2y = d2f.(x[2:end-1])
- k = 2.0 # constant diffusion
- δ = x[2:end] .- x[1:end-1]
- ∂²y = zeros(length(y)-2)
- lineardiffusion!(∂²y,y,δ,k)
- @test allfinite(∂²y)
- @test allequal(∂²y,k*d2y,atol=0.01)
+@inline function test_flux!_and_divergence!()
+ x = 0.0:0.25:1.0
+ xc = (x[1:end-1] .+ x[2:end])./2
+ Δx = x[2:end] .- x[1:end-1]
+ Δxc = xc[2:end] .- xc[1:end-1]
+ y = [0.0,1.0,0.0]
+ k = ones(length(x))
+ jy = zeros(length(x))
+ flux!(jy, y, Δxc, k)
+ @test jy[1] == jy[end] == zero(eltype(jy))
+ @test jy[2] ≈ -1/0.25
+ @test jy[3] ≈ 1/0.25
+ div = zeros(length(y))
+ divergence!(div, jy, Δx)
+ @test allequal(div, [1/0.25^2,-2/0.25^2,1/0.25^2])
end
-@testset "nonlineardiffusion!" begin
+@inline function test_nonlineardiffusion!()
f(x) = (1/6)x^3 # function to differntiate
df(x) = (1/2)x^2 # analytical 1st derivative
d2f(x) = x # analytical 2nd derivative
@@ -40,12 +31,25 @@ end
x = exp.(0.0:0.001:0.05)
xc = (x[1:end-1].+x[2:end])/2
y = f.(xc)
- kâ‚“ = k.(x)[2:end-1]
- d2y = d2kf.(xc[2:end-1]) # analytical solution
- δx = x[2:end] .- x[1:end-1]
- δxc = xc[2:end] .- xc[1:end-1]
- out = zeros(length(y)-2)
- nonlineardiffusion!(out,y,δxc,kₓ,δx[2:end-1])
- @test allfinite(out)
- @test allequal(out,d2y,atol=0.01)
+ kâ‚“ = k.(x)
+ d2y = d2kf.(xc) # analytical solution
+ Δx = x[2:end] .- x[1:end-1]
+ jy = zeros(length(x))
+ # manually set boundary fluxes
+ jy[1] = -kâ‚“[1]*df(x[1])
+ jy[end] = -kâ‚“[end]*df(x[end])
+ Δxc = xc[2:end] .- xc[1:end-1]
+ div = zeros(length(y))
+ nonlineardiffusion!(div, jy, y, Δxc, kₓ, Δx)
+ @test allfinite(div)
+ @test allequal(div, d2y, atol=0.01)
+end
+
+@testset "Math" begin
+ @testset "flux! and divergence!" begin
+ test_flux!_and_divergence!()
+ end
+ @testset "Non-linear diffusion" begin
+ test_nonlineardiffusion!()
+ end
end
diff --git a/test/Physics/HeatConduction/heat_conduction_tests.jl b/test/Physics/HeatConduction/heat_conduction_tests.jl
index 2d41dc55868986b19478569d660df578c825fd82..65eb79cee7201dbf8c999b5fd614a2e961dc765d 100644
--- a/test/Physics/HeatConduction/heat_conduction_tests.jl
+++ b/test/Physics/HeatConduction/heat_conduction_tests.jl
@@ -1,4 +1,5 @@
using CryoGrid
+using CryoGrid.Numerics: flux!, nonlineardiffusion!
using Dates
using DiffEqBase
using OrdinaryDiffEq
@@ -16,12 +17,12 @@ include("../../types.jl")
k = collect(LinRange(0.5,5.0,length(x)))u"W/m/K"
ΔT = Δ(xc)
Δk = Δ(x)
- ∂H = zeros(length(T₀))u"W/m^3"
- @inferred heatconduction!(∂H,T₀,ΔT,k,Δk)
+ jH = zeros(length(x))u"W/m^2"
+ dH = zeros(length(Tâ‚€))u"W/m^3"
+ @inferred nonlineardiffusion!(dH, jH, T₀, ΔT, k, Δk)
# conditions based on initial temperature gradient
- @test ∂H[1] < 0.0u"W/m^3"
- @test ∂H[end] > 0.0u"W/m^3"
- @test sum(∂H) <= 0.0u"W/m^3"
+ @test dH[1] < 0.0u"W/m^3"
+ @test dH[end] > 0.0u"W/m^3"
end
# check boundary fluxes
@testset "Boundary fluxes" begin
@@ -35,37 +36,42 @@ include("../../types.jl")
bc = ConstantBC(Heat, CryoGrid.Dirichlet, 0.0u"°C")
@testset "top: +, bot: -" begin
T₀ = Vector(LinRange(-23,27,length(xc)))u"°C"
- ∂H = zeros(length(T₀))u"W/m^3"
- state = (T=T₀,k=k,dH=∂H,grid=x,grids=(T=xc,k=x),t=0.0)
+ jH = zeros(length(x))u"W/m^2"
+ dH = zeros(length(Tâ‚€))u"W/m^3"
+ state = (T=Tâ‚€,k=k,dH=dH,jH=jH,grid=x,grids=(T=xc,k=x),t=0.0)
@test boundaryflux(bc,Top(),heat,sub,state,state) > 0.0u"W/m^2"
@test boundaryflux(bc,Bottom(),heat,sub,state,state) < 0.0u"W/m^2"
end
@testset "top: -, bot: +" begin
T₀ = Vector(LinRange(27,-23,length(xc)))u"°C"
- ∂H = zeros(length(T₀))u"W/m^3"
- state = (T=T₀,k=k,dH=∂H,grid=x,grids=(T=xc,k=x),t=0.0)
+ jH = zeros(length(x))u"W/m^2"
+ dH = zeros(length(Tâ‚€))u"W/m^3"
+ state = (T=Tâ‚€,k=k,dH=dH,grid=x,grids=(T=xc,k=x),t=0.0)
@test boundaryflux(bc,Top(),heat,sub,state,state) < 0.0u"W/m^2"
@test boundaryflux(bc,Bottom(),heat,sub,state,state) > 0.0u"W/m^2"
end
@testset "inner edge boundary (positive)" begin
T₀ = Vector(sin.(ustrip.(xc).*π))u"°C"
- ∂H = zeros(length(T₀))u"W/m^3"
- heatconduction!(∂H,T₀,ΔT,k,Δk)
- @test ∂H[1] > 0.0u"W/m^3"
- @test ∂H[end] > 0.0u"W/m^3"
+ jH = zeros(length(x))u"W/m^2"
+ dH = zeros(length(Tâ‚€))u"W/m^3"
+ nonlineardiffusion!(dH,jH,T₀,ΔT,k,Δk)
+ @test dH[1] > 0.0u"W/m^3"
+ @test dH[end] > 0.0u"W/m^3"
end
@testset "inner edge boundary (negative)" begin
T₀ = Vector(-sin.(ustrip.(xc).*π))u"°C"
- ∂H = zeros(length(T₀))u"W/m^3"
- heatconduction!(∂H,T₀,ΔT,k,Δk)
- @test ∂H[1] < 0.0u"W/m^3"
- @test ∂H[end] < 0.0u"W/m^3"
+ jH = zeros(length(x))u"W/m^2"
+ dH = zeros(length(Tâ‚€))u"W/m^3"
+ nonlineardiffusion!(dH,jH,T₀,ΔT,k,Δk)
+ @test dH[1] < 0.0u"W/m^3"
+ @test dH[end] < 0.0u"W/m^3"
end
@testset "Neumann boundary" begin
bc = ConstantBC(Heat, CryoGrid.Neumann, -1.0u"W/m^2")
T₀ = Vector(LinRange(-23,27,length(xc)))u"°C"
- ∂H = zeros(length(T₀))u"W/m^3"
- state = (T=T₀,k=k,dH=∂H,grid=x,grids=(T=xc,k=x),t=0.0)
+ jH = zeros(length(x))u"W/m^2"
+ dH = zeros(length(Tâ‚€))u"W/m^3"
+ state = (T=Tâ‚€,k=k,dH=dH,jH=jH,grid=x,grids=(T=xc,k=x),t=0.0)
@test boundaryflux(bc,Top(),heat,sub,state,state) == -1.0u"W/m^2"
end
end
@@ -98,17 +104,20 @@ end
sub = TestGroundLayer()
heat = Heat()
bc = ConstantBC(Heat, CryoGrid.Dirichlet, 0.0u"°C")
- function dTdt(T,p,t)
- dH = similar(T)u"W/m^3"
+ function dTdt(u,p,t)
+ dH = similar(u)u"W/m^3"
dH .= zero(eltype(dH))
- T_K = (T)u"°C"
- heatconduction!(dH,T_K,ΔT,k,Δk)
+ jH = similar(u, length(u)+1)u"W/m^2"
+ jH .= zero(eltype(jH))
+ T = (u)u"°C"
# compute boundary fluxes;
# while not correct in general, for this test case we can just re-use state for both layers.
- state = (T=T_K,k=k,dH=dH,grid=x,grids=(T=xc,k=x),t=t)
- dH[1] += boundaryflux(bc,Top(),heat,sub,state,state)/Δk[1]
- dH[end] += boundaryflux(bc,Bottom(),heat,sub,state,state)/Δk[end]
- # strip units from dH before returning it to the solver
+ state = (T=T,dH=dH,jH=jH,k=k,grid=x,grids=(T=xc,k=x),t=t)
+ interact!(Top(), bc, sub, heat, state, state)
+ interact!(sub, heat, Bottom(), bc, state, state)
+ prognosticstep!(sub, heat, state)
+ # strip units from dH before returning it to the solver;
+ # note that we do not need to divide by diffusivity since we assume it to be unity
return ustrip.(dH)
end
tspan = (0.0,0.5)
diff --git a/test/Physics/HeatConduction/runtests.jl b/test/Physics/HeatConduction/runtests.jl
index 989eb9a6316b90028fc882dba178de731f5371bc..3591176108ec5be3ce5a95a197ed3579b97b9e6d 100644
--- a/test/Physics/HeatConduction/runtests.jl
+++ b/test/Physics/HeatConduction/runtests.jl
@@ -1,4 +1,3 @@
@testset "Heat Conduction" begin
include("heat_conduction_tests.jl")
- include("sfcc_tests.jl")
end
diff --git a/test/Physics/HeatConduction/sfcc_bench.jl b/test/Physics/HeatConduction/sfcc_bench.jl
deleted file mode 100644
index 2c1b66a0505482318268875428ff45a82a8bf393..0000000000000000000000000000000000000000
--- a/test/Physics/HeatConduction/sfcc_bench.jl
+++ /dev/null
@@ -1,45 +0,0 @@
-using BenchmarkTools
-using CryoGrid
-using CryoGrid.Utils
-using Profile
-using Test
-
-function benchmarksfcc()
- θres = 0.05
- α = 4.0
- n = 2.0
- Tₘ = 0.0
- f = DallAmico(;θres,α,n) |> stripparams |> stripunits
- solver = SFCCNewtonSolver(α₀=1.0, τ=0.75, onfail=:error)
- sfcc = SFCC(f, solver)
- soil = Soil() |> stripparams |> stripunits
- heat = Heat(freezecurve=sfcc) |> stripparams |> stripunits
- L = heat.prop.L
- Lf = heat.prop.Lf
- # set up multi-grid-cell state vars
- T = [-5.0 for i in 1:10]
- θwi = Soils.soilcomponent(Val{:θwi}(), soil.para)
- θp = Soils.soilcomponent(Val{:θp}(), soil.para)
- θm = Soils.soilcomponent(Val{:θm}(), soil.para)
- θo = Soils.soilcomponent(Val{:θo}(), soil.para)
- θw = f.(T,Tₘ,θres,θp,θwi,Lf,α,n) # set liquid water content according to freeze curve
- C = heatcapacity.(soil,heat,θwi,θw,θm,θo)
- H = let T = T .+ 4.99,
- θw = f.(T,Tₘ,θres,θp,θwi,Lf,α,n),
- C = heatcapacity.(soil,heat,θwi,θw,θm,θo);
- enthalpy.(T,C,L,θw) # compute enthalpy at "true" temperature
- end
- state = (T=T,C=C,dHdT=similar(C),H=H,θw=θw,θwi=θwi)
- params = Utils.selectat(1, identity, Soils.sfccargs(f, soil, heat, state))
- f_params = tuplejoin((T[1],), params)
- # @btime $sfcc.f($f_params...)
- # @btime $sfcc.∇f($f_params...)
- # @btime Soils.residual($soil, $heat, $T[1], $H[1], $L, $sfcc.f, $params, $θwi, $θm, $θo)
- # @code_warntype HeatConduction.enthalpyinv(soil, heat, state, 1)
- # Soils.sfccsolve(solver, soil, heat, sfcc.f, sfcc.∇f, params, H[1], L, θwi, θm, θo, 0.0)
- @btime Soils.sfccsolve($solver, $soil, $heat, $sfcc.f, $sfcc.∇f, $params, $H[1], $L, $θwi, $θm, $θo, 0.0)
- # @btime Soils.residual($soil, $heat, $T[1], $H[1], $L, $sfcc.f, $params, $θwi, $θm, $θo)
- # res = @btime Soils.sfccsolve($solver, $soil, $heat, $sfcc.f, $sfcc.∇f, $params, $H[1], $L, $θwi, $θm, $θo)
- @btime freezethaw!($soil, $heat, $state)
- println("\nsolution: $(state.T[1]), $(state.θw[1])")
-end
diff --git a/test/Physics/HeatConduction/sfcc_tests.jl b/test/Physics/HeatConduction/sfcc_tests.jl
deleted file mode 100644
index 15c7577cfd7b40acfe9124959fd3cbb6a8043c79..0000000000000000000000000000000000000000
--- a/test/Physics/HeatConduction/sfcc_tests.jl
+++ /dev/null
@@ -1,244 +0,0 @@
-using CryoGrid
-using ComponentArrays
-using ForwardDiff
-using Setfield
-using Test
-
-@testset "SFCC" begin
- Tₘ = 0.0
- θres = 0.0
- soil = @pstrip Soil(para=Soils.CharacteristicFractions())
- θwi = Soils.soilcomponent(Val{:θwi}(), soil.para)
- θp = Soils.soilcomponent(Val{:θp}(), soil.para)
- θm = Soils.soilcomponent(Val{:θm}(), soil.para)
- θo = Soils.soilcomponent(Val{:θo}(), soil.para)
- @testset "McKenzie freeze curve" begin
- @testset "Sanity checks" begin
- f = @pstrip McKenzie() keep_units=true
- let θsat = 0.8,
- γ = 1.0u"K",
- Tₘ = 0.0u"°C";
- @test isapprox(f(-10.0u"°C",θsat,θsat,θres,Tₘ,γ), 0.0, atol=1e-6)
- @test f(0.0u"°C",θsat,θsat,θres,Tₘ,γ) ≈ θsat
- θw = f(-0.1u"°C",θsat,θsat,θres,Tₘ,γ)
- @test θw > 0.0 && θw < 1.0
- end
- end
- @testset "Newton solver checks" begin
- abstol = 1e-2
- reltol = 1e-4
- γ = 0.1
- f = @pstrip McKenzie()
- sfcc = SFCC(f, SFCCNewtonSolver(abstol=abstol, reltol=reltol))
- heat = @pstrip Heat(freezecurve=sfcc)
- heatcap = θw -> heatcapacity(soil, heat, θw, θwi - θw, 1-θwi-θm-θo, θm, θo)
- L = heat.prop.L
- @testset "Left tail" begin
- T = [-5.0]
- θw = f.(T,θp,θwi,θres,Tₘ,γ) # set liquid water content according to freeze curve
- C = heatcap.(θw)
- H = let T = T.+1.0,
- θw = f.(T,θp,θwi,θres,Tₘ,γ)
- C = heatcap.(θw);
- enthalpy.(T,C,L,θw) # compute enthalpy at "true" temperature
- end
- state = (T=T,C=C,dHdT=similar(C),H=H,θw=θw,)
- freezethaw!(soil, heat, state)
- @test all(abs.(T.-(H .- L.*θw)./C) .<= abstol)
- end
- @testset "Right tail" begin
- # set up single-grid-cell state vars
- T = [5.0]
- θw = f.(T,θp,θwi,θres,Tₘ,γ) # set liquid water content according to freeze curve
- C = heatcap.(θw)
- H = let T = T.-1.0,
- θw = f.(T,θp,θwi,θres,Tₘ,γ)
- C = heatcap.(θw);
- enthalpy.(T,C,L,θw) # compute enthalpy at "true" temperature
- end
- state = (T=T,C=C,dHdT=similar(C),H=H,θw=θw,)
- freezethaw!(soil, heat, state)
- @test all(abs.(T.-(H .- L.*θw)./C) .<= abstol)
- end
- @testset "Near zero" begin
- # set up single-grid-cell state vars
- T = [-0.05]
- θw = f.(T,θp,θwi,θres,Tₘ,γ) # set liquid water content according to freeze curve
- C = heatcap.(θw)
- H = let T = T.+0.04,
- θw = f.(T,θp,θwi,θres,Tₘ,γ)
- C = heatcap.(θw);
- enthalpy.(T,C,L,θw) # compute enthalpy at "true" temperature
- end
- state = (T=T,C=C,dHdT=similar(C),H=H,θw=θw,)
- freezethaw!(soil, heat, state)
- @test all(abs.(T.-(H .- L.*θw)./C) .<= abstol)
- end
- end
- end
- # TODO: DRY violation; a lot of this code is redundant and could possibly be
- # shared between different freeze curve tests.
- @testset "Westermann freeze curve" begin
- @testset "Sanity checks" begin
- f = @pstrip Westermann() keep_units=true
- let θtot = 0.8,
- δ = 0.1u"K",
- Tₘ = 0.0u"°C";
- @test isapprox(f(-10.0u"°C",θtot,θtot,θres,Tₘ,δ), 0.0, atol=1e-2)
- @test f(0.0u"°C",θtot,θtot,θres,Tₘ,δ) ≈ θtot
- θw = f(-0.1u"°C",θtot,θtot,θres,Tₘ,δ)
- @test θw > 0.0 && θw < 1.0
- end
- end
- @testset "Newton solver checks" begin
- abstol = 1e-2
- reltol = 1e-4
- δ = 0.1
- f = @pstrip Westermann()
- sfcc = SFCC(f, SFCCNewtonSolver(abstol=abstol, reltol=reltol))
- heat = @pstrip Heat(freezecurve=sfcc)
- heatcap = θw -> heatcapacity(soil, heat, θw, θwi - θw, 1-θwi-θm-θo, θm, θo)
- L = heat.prop.L
- @testset "Left tail" begin
- T = [-5.0]
- θw = f.(T,θp,θwi,θres,Tₘ,δ) # set liquid water content according to freeze curve
- C = heatcap.(θw)
- H = let T = T.+1.0,
- θw = f.(T,θp,θwi,θres,Tₘ,δ)
- C = heatcap.(θw);
- enthalpy.(T,C,L,θw) # compute enthalpy at "true" temperature
- end
- state = (T=T,C=C,dHdT=similar(C),H=H,θw=θw,)
- freezethaw!(soil, heat, state)
- @test all(abs.(T.-(H .- L.*θw)./C) .<= abstol)
- end
- @testset "Right tail" begin
- T = [5.0]
- θw = f.(T,θp,θwi,θres,Tₘ,δ) # set liquid water content according to freeze curve
- C = heatcap.(θw)
- H = let T = T.-1,
- θw = f.(T,θp,θwi,θres,Tₘ,δ)
- C = heatcap.(θw);
- enthalpy.(T,C,L,θw) # compute enthalpy at "true" temperature
- end
- state = (T=T,C=C,dHdT=similar(C),H=H,θw=θw,)
- freezethaw!(soil, heat, state)
- @test all(abs.(T.-(H .- L.*θw)./C) .<= abstol)
- end
- @testset "Near zero" begin
- T = [-0.05]
- θw = f.(T,θp,θwi,θres,Tₘ,δ) # set liquid water content according to freeze curve
- C = heatcap.(θw)
- H = let T = T.+0.04,
- θw = f.(T,θp,θwi,θres,Tₘ,δ)
- C = heatcap.(θw);
- enthalpy.(T,C,L,θw) # compute enthalpy at "true" temperature
- end
- state = (T=T,C=C,dHdT=similar(C),H=H,θw=θw,)
- freezethaw!(soil, heat, state)
- @test all(abs.(T.-(H .- L.*θw)./C) .<= abstol)
- end
- end
- end
- @testset "Dall'Amico freeze curve" begin
- @testset "Sanity checks" begin
- f = @pstrip DallAmico() keep_units=true
- let θsat = 0.8,
- α = 4.0u"1/m",
- n = 2.0,
- Tₘ = 0.0u"°C";
- Lf = stripparams(Heat()).prop.Lf
- @test isapprox(f(-10.0u"°C",θsat,θsat,Lf,θres,Tₘ,α,n), 0.0, atol=1e-3)
- @test f(0.0u"°C",θsat,θsat,Lf,θres,Tₘ,α,n) ≈ θsat
- θw = f(-0.1u"°C",θsat,θsat,Lf,θres,Tₘ,α,n)
- @test θw > 0.0 && θw < 1.0
- end
- end
- @testset "Newton solver checks" begin
- abstol = 1e-2
- reltol = 1e-4
- θsat = 0.8
- α = 4.0
- n = 2.0
- Tₘ = 0.0
- f = @pstrip DallAmico()
- sfcc = SFCC(f, SFCCNewtonSolver(abstol=abstol, reltol=reltol))
- heat = @pstrip Heat(freezecurve=sfcc)
- heatcap = θw -> heatcapacity(soil, heat, θw, θwi - θw, 1-θwi-θm-θo, θm, θo)
- L = heat.prop.L
- Lf = heat.prop.Lf
- @testset "Left tail" begin
- T = [-5.0]
- θw = f.(T,θp,θwi,Lf,θres,Tₘ,α,n) # set liquid water content according to freeze curve
- C = heatcap.(θw)
- H = let T = T.+1.0,
- θw = f.(T,θp,θwi,Lf,θres,Tₘ,α,n)
- C = heatcap.(θw);
- enthalpy.(T,C,L,θw) # compute enthalpy at "true" temperature
- end
- state = (T=T,C=C,dHdT=similar(C),H=H,θw=θw,)
- @inferred freezethaw!(soil, heat, state)
- @test all(abs.(T.-(H .- L.*θw)./C) .<= abstol)
- end
- @testset "Right tail" begin
- T = [5.0]
- θw = f.(T,θp,θwi,Lf,θres,Tₘ,α,n) # set liquid water content according to freeze curve
- C = heatcap.(θw)
- H = let T = T.-1.0,
- θw = f.(T,θp,θwi,Lf,θres,Tₘ,α,n)
- C = heatcap.(θw);
- enthalpy.(T,C,L,θw) # compute enthalpy at "true" temperature
- end
- state = (T=T,C=C,dHdT=similar(C),H=H,θw=θw,)
- @inferred freezethaw!(soil, heat, state)
- @test all(abs.(T.-(H .- L.*θw)./C) .<= abstol)
- end
- @testset "Near zero" begin
- T = [-0.05]
- θw = f.(T,θp,θwi,Lf,θres,Tₘ,α,n) # set liquid water content according to freeze curve
- C = heatcap.(θw)
- H = let T = T.+0.04,
- θw = f.(T,θp,θwi,Lf,θres,Tₘ,α,n)
- C = heatcap.(θw);
- enthalpy.(T,C,L,θw) # compute enthalpy at "true" temperature
- end
- state = (T=T,C=C,dHdT=similar(C),H=H,θw=θw,)
- @inferred freezethaw!(soil, heat, state)
- @test all(abs.(T.-(H .- L.*θw)./C) .<= abstol)
- end
- end
- end
- @testset "Newton solver autodiff" begin
- # set up
- θsat = 0.8
- θres = 0.05
- Tₘ = 0.0
- γ = 0.1
- f = @pstrip McKenzie()
- sfcc = SFCC(f, SFCCNewtonSolver())
- heat = @pstrip Heat(freezecurve=sfcc)
- heatcap = θw -> heatcapacity(soil, heat, θw, θwi - θw, 1-θwi-θm-θo, θm, θo)
- L = heat.prop.L
- T = [-0.1]
- θw = f.(T,θp,θwi,θres,Tₘ,γ) # set liquid water content according to freeze curve
- C = heatcap.(θw)
- H = enthalpy.(T.+0.09,C,L,θw) # compute enthalpy at +1 degree
- # test gradients
- p = ComponentArray(γ=γ)
- ∂f∂p = ForwardDiff.gradient(p -> sum(f.(T,θp,θwi,θres,Tₘ,p.γ)), p)
- @test all(isfinite.(∂f∂p))
- function F(p)
- T_ = similar(T,eltype(p))
- T_ .= T
- C = similar(C,eltype(p))
- θw = similar(θw,eltype(p))
- state = (T=T_,C=C,dHdT=similar(C),H=H,θw=θw,)
- @set! heat.freezecurve.f.γ = p.γ
- freezethaw!(soil, heat, state)
- state.T[1]
- end
- p = ComponentArray(γ=[γ])
- ∂F∂p = ForwardDiff.gradient(F, p)
- @test all(isfinite.(∂F∂p))
- end
-end;
diff --git a/test/Physics/Sources/runtests.jl b/test/Physics/Sources/runtests.jl
index 87e4afab7c45772eb1df55991e46f39f106273e6..500be9d0bdda8042328807fcb83eaa485b62dd05 100644
--- a/test/Physics/Sources/runtests.jl
+++ b/test/Physics/Sources/runtests.jl
@@ -6,7 +6,7 @@ include("../../types.jl")
@testset "Sources" begin
layer = TestGroundLayer()
@testset "Constant" begin
- heatsource = Source(Heat, Sources.Constant())
+ heatsource = parameterize(Source(Heat, Sources.Constant()))
heatsource = CryoGrid.update(heatsource, [1.0u"W/m^3"])
state = (dH=zeros(100)u"W/m^3",)
prognosticstep!(layer, heatsource, state)
@@ -17,7 +17,7 @@ include("../../types.jl")
end
@testset "Periodic" begin
level, amp, freq, shift = 0.0u"W/m^3", 2.0u"W/m^3", 0.5u"Hz", π/2
- heatsource = Source(Heat, Sources.Periodic())
+ heatsource = parameterize(Source(Heat, Sources.Periodic()))
heatsource = CryoGrid.update(heatsource, [amp, freq, shift, level])
state = (dH=zeros(100)u"W/m^3", t=1.0u"s")
prognosticstep!(layer, heatsource, state)