Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • cryogrid/cryogridjulia
  • jnitzbon/cryogridjulia
  • bgroenke/cryogridjulia
3 results
Show changes
Hide whitespace changes
Inline Side-by-side
src/Strat/tile.jl
View file @ 738aa942
......@@ -34,38 +34,38 @@ Out-of-place step function for tile `T`. Computes and returns du/dt as vector wi
step(::T,u,p,t) where {T<:AbstractTile} = error("no implementation of out-of-place step for $T")
"""
Tile{TStrat,TGrid,TStates,TInits,TCallbacks,iip,obsv} <: AbstractTile{iip}
Tile{TStrat,TGrid,TStates,TInits,TEvents,iip,obsv} <: AbstractTile{iip}
Defines the full specification of a single CryoGrid tile; i.e. stratigraphy, grid, and state variables.
"""
struct Tile{TStrat,TGrid,TStates,TInits,TCallbacks,iip,obsv} <: AbstractTile{iip}
struct Tile{TStrat,TGrid,TStates,TInits,TEvents,iip,obsv} <: AbstractTile{iip}
strat::TStrat # stratigraphy
grid::TGrid # grid
state::TStates # state variables
inits::TInits # initializers
callbacks::TCallbacks # callbacks
events::TEvents # events
hist::StateHistory # mutable "history" type for state tracking
function Tile(
strat::TStrat,
grid::TGrid,
state::TStates,
inits::TInits,
callbacks::TCallbacks,
events::TEvents,
hist::StateHistory=StateHistory(),
iip::InPlaceMode=inplace,
observe::Tuple{Vararg{Symbol}}=()) where
{TStrat<:Stratigraphy,TGrid<:Grid{Edges},TStates<:VarStates,TInits<:Tuple,TCallbacks<:NamedTuple}
new{TStrat,TGrid,TStates,TInits,TCallbacks,iip,observe}(strat,grid,state,inits,callbacks,hist)
{TStrat<:Stratigraphy,TGrid<:Grid{Edges},TStates<:VarStates,TInits<:Tuple,TEvents<:NamedTuple}
new{TStrat,TGrid,TStates,TInits,TEvents,iip,observe}(strat,grid,state,inits,events,hist)
end
end
ConstructionBase.constructorof(::Type{Tile{TStrat,TGrid,TStates,TInits,TCallbacks,iip,obsv}}) where {TStrat,TGrid,TStates,TInits,TCallbacks,iip,obsv} =
(strat, grid, state, inits, callbacks, hist) -> Tile(strat, grid, state, inits, callbacks, hist, iip, obsv)
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)
# mark only stratigraphy and initializers fields as flattenable
Flatten.flattenable(::Type{<:Tile}, ::Type{Val{:strat}}) = true
Flatten.flattenable(::Type{<:Tile}, ::Type{Val{:inits}}) = true
Flatten.flattenable(::Type{<:Tile}, ::Type{Val{:callbacks}}) = true
Flatten.flattenable(::Type{<:Tile}, ::Type{Val{:events}}) = true
Flatten.flattenable(::Type{<:Tile}, ::Type{Val{name}}) where name = false
Base.show(io::IO, ::MIME"text/plain", tile::Tile{TStrat,TGrid,TStates,TInits,TCallbacks,iip,obsv}) where {TStrat,TGrid,TStates,TInits,TCallbacks,iip,obsv} = print(io, "Tile ($iip) with layers $(map(componentname, components(tile.strat))), observables=$obsv, $TGrid, $TStrat")
Base.show(io::IO, ::MIME"text/plain", tile::Tile{TStrat,TGrid,TStates,TInits,TEvents,iip,obsv}) where {TStrat,TGrid,TStates,TInits,TEvents,iip,obsv} = print(io, "Tile ($iip) with layers $(map(componentname, components(tile.strat))), observables=$obsv, $TGrid, $TStrat")
"""
Tile(
......@@ -91,7 +91,7 @@ function Tile(
chunksize=nothing,
) where {A<:AbstractArray}
vars = OrderedDict()
callbacks = OrderedDict()
events = OrderedDict()
components = OrderedDict()
for comp in stripunits(strat)
name = componentname(comp)
......@@ -106,15 +106,15 @@ function Tile(
else
components[name] = comp
end
# callbacks
cbs = CryoGrid.callbacks(comp.layer, comp.processes)
# events
cbs = CryoGrid.events(comp.layer, comp.process)
cbparams = ModelParameters.params(cbs)
if length(cbparams) > 0
m_cbs = Model(cbs)
m_cbs[:layer] = repeat([name], length(params))
callbacks[name] = parent(m_cbs)
events[name] = parent(m_cbs)
else
callbacks[name] = cbs
events[name] = cbs
end
end
# set :layer field on initializer parameters (if any)
......@@ -139,95 +139,103 @@ function Tile(
chunksize = isnothing(chunksize) ? length(para) : chunksize
states = VarStates(ntvars, Grid(dustrip(grid), grid.geometry), chunksize, arrayproto)
isempty(inits) && @warn "No initializers provided. State variables without initializers will be set to zero by default."
Tile(strat, grid, states, inits, (;callbacks...), StateHistory(), iip, Tuple(observe))
Tile(strat, grid, states, inits, (;events...), StateHistory(), iip, Tuple(observe))
end
Tile(strat::Stratigraphy, grid::Grid{Cells}; kwargs...) = Tile(strat, edges(grid); kwargs...)
Tile(strat::Stratigraphy, grid::Grid{Edges,<:Numerics.Geometry,T}; kwargs...) where {T} = error("grid must have values with units of length, e.g. try using `Grid((x)u\"m\")` where `x` are your grid points.")
"""
step!(_tile::Tile{TStrat,TGrid,TStates,TInits,TEvents,inplace,obsv}, _du, _u, p, t) where {TStrat,TGrid,TStates,TInits,TEvents,obsv}
Generated step function (i.e. du/dt) for any arbitrary Tile. Specialized code is generated and compiled
on the fly via the @generated macro to ensure type stability. The generated code updates each layer in the stratigraphy
in sequence, i.e for each layer 1 < i < N:
in sequence, i.e for each layer 1 <= i <= N:
diagnosticstep!(layer i, ...)
interact!(layer i-1, ...)
interact!(layer i, ..., layer i+1, ...)
prognosticstep!(layer i, ...)
Note for developers: All sections of code wrapped in quote..end blocks are generated. Code outside of quote blocks
is only executed during compilation and will not appear in the compiled version.
"""
@generated function step!(tile::Tile{TStrat,TGrid,TStates,TInits,TCallbacks,inplace,obsv}, _du, _u, p, t) where {TStrat,TGrid,TStates,TInits,TCallbacks,obsv}
@generated function step!(_tile::Tile{TStrat,TGrid,TStates,TInits,TEvents,inplace,obsv}, _du, _u, p, t) where {TStrat,TGrid,TStates,TInits,TEvents,obsv}
nodetyps = componenttypes(TStrat)
N = length(nodetyps)
expr = Expr(:block)
# Declare variables
@>> quote
quote
_du .= zero(eltype(_du))
du = ComponentArray(_du, getaxes(tile.state.uproto))
u = ComponentArray(_u, getaxes(tile.state.uproto))
tile = updateparams(tile, u, p, t)
du = ComponentArray(_du, getaxes(_tile.state.uproto))
u = ComponentArray(_u, getaxes(_tile.state.uproto))
tile = updateparams(_tile, u, p, t)
strat = tile.strat
state = TileState(tile.state, boundaries(strat), u, du, t, Val{inplace}())
end push!(expr.args)
end |> Base.Fix1(push!, expr.args)
# Generate identifiers for each component
names = map(i -> componentname(nodetyps[i]), 1:N)
comps = map(n -> Symbol(n,:comp), names)
states = map(n -> Symbol(n,:state), names)
layers = map(n -> Symbol(n,:layer), names)
procs = map(n -> Symbol(n,:process), names)
# Initialize variables for all layers
for i in 1:N
n = componentname(nodetyps[i])
nstate = Symbol(n,:state)
ncomp = Symbol(n, :comp)
nlayer = Symbol(n,:layer)
nprocess = Symbol(n,:process)
@>> quote
$nstate = state.$n
$ncomp = components(strat)[$i]
$nlayer = $ncomp.layer
$nprocess = $ncomp.processes
end push!(expr.args)
quote
$(states[i]) = state.$(names[i])
$(comps[i]) = components(strat)[$i]
$(layers[i]) = $(comps[i]).layer
$(procs[i]) = $(comps[i]).process
end |> Base.Fix1(push!, expr.args)
end
# Diagnostic step
for i in 1:N
n = componentname(nodetyps[i])
nstate = Symbol(n,:state)
nlayer = Symbol(n,:layer)
nprocess = Symbol(n,:process)
@>> quote
diagnosticstep!($nlayer, $nstate)
diagnosticstep!($nlayer, $nprocess, $nstate)
end push!(expr.args)
quote
diagnosticstep!($(layers[i]), $(states[i]))
diagnosticstep!($(layers[i]), $(procs[i]), $(states[i]))
end |> Base.Fix1(push!, expr.args)
end
# Interact
can_interact_exprs = map(i -> :(CryoGrid.thickness($(layers[i]), $(states[i])) > 0), 1:N)
quote
can_interact = tuple($(can_interact_exprs...))
end |> Base.Fix1(push!, expr.args)
# We only invoke interact! on pairs of layers for which the following are satisfied:
# 1) Both layers have thickness > 0 (i.e. they occupy non-zero space in the stratigraphy)
# 2) Both layers are adjacent, or more crudely, all layers in between (if any) have zero thickness;
# In order to make this type stable, we pre-generate all possible interact! expressions in the
# downward direction and add a check to each one.
for i in 1:N-1
n1,n2 = componentname(nodetyps[i]), componentname(nodetyps[i+1])
n1state, n2state = Symbol(n1,:state), Symbol(n2,:state)
n1layer, n2layer = Symbol(n1,:layer), Symbol(n2,:layer)
n1process, n2process = Symbol(n1,:process), Symbol(n2,:process)
@>> quote
interact!($n1layer,$n1process,$n2layer,$n2process,$n1state,$n2state)
end push!(expr.args)
for j in i+1:N
if (i == 1 && j == 2) || (i == N-1 && j == N)
# always apply top and bottom interactions
quote
interact!($(layers[i]),$(procs[i]),$(layers[j]),$(procs[j]),$(states[i]),$(states[j]))
end
else
innerchecks = j > i+1 ? map(k -> :(can_interact[$k]), i+1:j-1) : :(false)
quote
if can_interact[$i] && can_interact[$j] && !any(tuple($(innerchecks...)))
interact!($(layers[i]),$(procs[i]),$(layers[j]),$(procs[j]),$(states[i]),$(states[j]))
end
end
end |> Base.Fix1(push!, expr.args)
end
end
# Prognostic step
for i in 1:N
n = componentname(nodetyps[i])
nstate = Symbol(n,:state)
nlayer = Symbol(n,:layer)
nprocess = Symbol(n,:process)
@>> quote
prognosticstep!($nlayer,$nprocess,$nstate)
end push!(expr.args)
quote
prognosticstep!($(layers[i]),$(procs[i]),$(states[i]))
end |> Base.Fix1(push!, expr.args)
end
# Observables
for i in 1:N
n = componentname(nodetyps[i])
nstate = Symbol(n,:state)
nlayer = Symbol(n,:layer)
nprocess = Symbol(n,:process)
for name in obsv
nameval = Val{name}()
@>> quote
observe($nameval,$nlayer,$nprocess,$nstate)
end push!(expr.args)
quote
observe($nameval,$(layers[i]),$(procs[i]),$(states[i]))
end |> Base.Fix1(push!, expr.args)
end
end
# make sure compiled method returns no value
@>> :(return nothing) push!(expr.args)
push!(expr.args, :(return nothing))
# emit generated expression block
return expr
end
......@@ -237,8 +245,8 @@ end
Calls `initialcondition!` on all layers/processes and returns the fully constructed u0 and du0 states.
"""
initialcondition!(tile::Tile, tspan::NTuple{2,DateTime}, p::AbstractVector, args...) = initialcondition!(tile, convert_tspan(tspan), p)
@generated function initialcondition!(tile::Tile{TStrat,TGrid,TStates,TInits,TCallbacks,iip,obsv}, tspan::NTuple{2,Float64}, p::AbstractVector) where {TStrat,TGrid,TStates,TInits,TCallbacks,iip,obsv}
CryoGrid.initialcondition!(tile::Tile, tspan::NTuple{2,DateTime}, p::AbstractVector, args...) = initialcondition!(tile, convert_tspan(tspan), p)
@generated function CryoGrid.initialcondition!(tile::Tile{TStrat,TGrid,TStates,TInits,TEvents,iip,obsv}, tspan::NTuple{2,Float64}, p::AbstractVector) where {TStrat,TGrid,TStates,TInits,TEvents,iip,obsv}
nodetyps = componenttypes(TStrat)
N = length(nodetyps)
expr = Expr(:block)
......@@ -277,9 +285,9 @@ initialcondition!(tile::Tile, tspan::NTuple{2,DateTime}, p::AbstractVector, args
$n1state = state.$n1
$n2state = state.$n2
$n1layer = components(strat)[$i].layer
$n1process = components(strat)[$i].processes
$n1process = components(strat)[$i].process
$n2layer = components(strat)[$(i+1)].layer
$n2process = components(strat)[$(i+1)].processes
$n2process = components(strat)[$(i+1)].process
end push!(expr.args)
# invoke initialcondition! for each layer, then for both (similar to interact!);
# initialcondition! is called for layer1 only on the first iteration to avoid duplicate invocations
......@@ -302,6 +310,68 @@ initialcondition!(tile::Tile, tspan::NTuple{2,DateTime}, p::AbstractVector, args
end push!(expr.args)
return expr
end
@generated function CryoGrid.timestep(_tile::Tile{TStrat,TGrid,TStates,TInits,TEvents,iip,obsv}, _du, _u, p, t) where {TStrat,TGrid,TStates,TInits,TEvents,iip,obsv}
nodetyps = componenttypes(TStrat)
N = length(nodetyps)
expr = Expr(:block)
# Declare variables
quote
du = ComponentArray(_du, getaxes(_tile.state.uproto))
u = ComponentArray(_u, getaxes(_tile.state.uproto))
tile = updateparams(_tile, u, p, t)
strat = tile.strat
state = TileState(tile.state, boundaries(strat), u, du, t, Val{iip}())
end |> Base.Fix1(push!, expr.args)
# Generate identifiers for each component
names = map(i -> componentname(nodetyps[i]), 1:N)
comps = map(n -> Symbol(n,:comp), names)
states = map(n -> Symbol(n,:state), names)
layers = map(n -> Symbol(n,:layer), names)
procs = map(n -> Symbol(n,:process), names)
# Initialize variables for all layers
for i in 1:N
quote
$(states[i]) = state.$(names[i])
$(comps[i]) = components(strat)[$i]
$(layers[i]) = $(comps[i]).layer
$(procs[i]) = $(comps[i]).process
end |> Base.Fix1(push!, expr.args)
end
push!(expr.args, :(dtmax = Inf))
# Retrieve minimum timesteps
for i in 1:N
quote
dtmax = min(dtmax, timestep($(layers[i]), $(procs[i]), $(states[i])))
end |> Base.Fix1(push!, expr.args)
end
push!(expr.args, :(return dtmax))
# emit generated expression block
return expr
end
"""
domain(tile::Tile)
Returns a function `isoutofdomain(u,p,t)` which checks whether any prognostic variable values
in `u` are outside of their respective domains. Only variables which have at least one finite
domain endpoint are checked; variables with unbounded domains are ignored.
"""
function domain(tile::Tile)
# select only variables which have a finite domain
vars = filter(x -> any(map(isfinite, extrema(vardomain(x)))), filter(isprognostic, variables(tile)))
function isoutofdomain(_u,p,t)::Bool
u = withaxes(_u, tile)
for var in vars
domain = vardomain(var)
uvar = getproperty(u, varname(var))
@inbounds for i in 1:length(uvar)
if uvar[i] ∉ domain
return true
end
end
end
return false
end
end
"""
getvar(name::Symbol, tile::Tile, u)
getvar(::Val{name}, tile::Tile, u)
......@@ -318,8 +388,8 @@ Numerics.getvar(::Val{name}, tile::Tile, u) where name = getvar(Val{name}(), til
Constructs a `LayerState` representing the full state of `layername` given `tile`, state vectors `u` and `du`, and the
time step `t`.
"""
getstate(layername::Symbol, tile::Tile, u, du, t) = getstate(Val{layername}(), tile, u, du, t)
function getstate(::Val{layername}, tile::Tile{TStrat,TGrid,<:VarStates{layernames},TInits,TCallbacks,iip}, _u, _du, t) where {layername,TStrat,TGrid,TInits,TCallbacks,iip,layernames}
getstate(layername::Symbol, tile::Tile, u, du, t, dt=nothing) = getstate(Val{layername}(), tile, u, du, t, dt)
function getstate(::Val{layername}, tile::Tile{TStrat,TGrid,<:VarStates{layernames},TInits,TEvents,iip}, _u, _du, t, dt=nothing) where {layername,TStrat,TGrid,TInits,TEvents,iip,layernames}
du = ComponentArray(_du, getaxes(tile.state.uproto))
u = ComponentArray(_u, getaxes(tile.state.uproto))
i = 1
......@@ -330,14 +400,14 @@ function getstate(::Val{layername}, tile::Tile{TStrat,TGrid,<:VarStates{layernam
end
end
z = boundarypairs(map(ustrip, stripparams(boundaries(tile.strat))), ustrip(tile.grid[end]))[i]
return LayerState(tile.state, z, u, du, t, Val{layername}(), Val{iip}())
return LayerState(tile.state, z, u, du, t, dt, Val{layername}(), Val{iip}())
end
"""
variables(tile::Tile)
Returns a tuple of all variables defined in the tile.
"""
variables(tile::Tile) = Tuple(unique(Flatten.flatten(tile.state.vars, Flatten.flattenable, Var)))
CryoGrid.variables(tile::Tile) = Tuple(unique(Flatten.flatten(tile.state.vars, Flatten.flattenable, Var)))
"""
parameters(tile::Tile; kwargs...)
......@@ -352,7 +422,7 @@ as `setup.uproto`.
"""
withaxes(u::AbstractArray, tile::Tile) = ComponentArray(u, getaxes(tile.state.uproto))
withaxes(u::ComponentArray, ::Tile) = u
function getstate(tile::Tile{TStrat,TGrid,TStates,TInits,TCallbacks,iip}, _u, _du, t) where {TStrat,TGrid,TStates,TInits,TCallbacks,iip}
function getstate(tile::Tile{TStrat,TGrid,TStates,TInits,TEvents,iip}, _u, _du, t) where {TStrat,TGrid,TStates,TInits,TEvents,iip}
du = ComponentArray(_du, getaxes(tile.state.uproto))
u = ComponentArray(_u, getaxes(tile.state.uproto))
return TileState(tile.state, map(ustrip ∘ stripparams, boundaries(tile.strat)), u, du, t, Val{iip}())
......@@ -363,20 +433,22 @@ end
Replaces all `ModelParameters.AbstractParam` values in `tile` with their (possibly updated) value from `p`.
Subsequently evaluates and replaces all nested `DynamicParameterization`s.
"""
function updateparams(tile::Tile, u, p, t)
tile_updated = Flatten.reconstruct(tile, p, ModelParameters.AbstractParam, Flatten.IGNORE)
dynamic_ps = Flatten.flatten(tile_updated, Flatten.flattenable, DynamicParameterization, Flatten.IGNORE)
function updateparams(tile::Tile{TStrat,TGrid,TStates}, u, p, t) where {TStrat,TGrid,TStates}
# unfortunately, reconstruct causes allocations due to a mysterious dynamic dispatch when returning the result of _reconstruct;
# I really don't know why, could be a compiler bug, but it doesn't happen if we call the internal _reconstruct directly soooo....
tile_updated = Flatten._reconstruct(tile, p, Flatten.flattenable, ModelParameters.AbstractParam, Union{TGrid,TStates,StateHistory},1)[1]
dynamic_ps = Flatten.flatten(tile_updated, Flatten.flattenable, DynamicParameterization, Union{TGrid,TStates,StateHistory})
# TODO: perhaps should allow dependence on local layer state;
# this would likely require per-layer deconstruction/reconstruction of `StratComponent`s in order to
# build the `LayerState`s and evaluate the dynamic parameters in a fully type stable manner.
dynamic_values = map(d -> d(u, t), dynamic_ps)
return Flatten.reconstruct(tile_updated, dynamic_values, DynamicParameterization, Flatten.IGNORE)
return Flatten._reconstruct(tile_updated, dynamic_values, Flatten.flattenable, DynamicParameterization, Union{TGrid,TStates,StateHistory},1)[1]
end
"""
Collects and validates all declared variables (`Var`s) for the given stratigraphy component.
"""
function _collectvars(@nospecialize(comp::StratComponent))
layer, process = comp.layer, comp.processes
layer, process = comp.layer, comp.process
declared_vars = variables(layer, process)
nested_vars = Flatten.flatten(comp, Flatten.flattenable, Var)
all_vars = tuplejoin(declared_vars, nested_vars)
......
src/Utils/Utils.jl
View file @ 738aa942
......@@ -61,6 +61,7 @@ Converts temperature `x` to Kelvin. If `x` has units, `uconvert` is used. Otherw
"""
normalize_temperature(x) = x + 273.15
normalize_temperature(x::TempQuantity) = uconvert(u"K", x)
normalize_temperature(x::Param) = stripparams(x) |> normalize_temperature
"""
Provides implementation of `Base.iterate` for structs.
......@@ -106,6 +107,10 @@ Concatenates one or more tuples together; should generally be type stable.
Helper method for generalizing between arrays and scalars. Without an index, retrieves
the first element of `x` if `x` is an array, otherwise simply returning `x`. If an index
`i`, is specified, returns the `i`th value of `x` if `x` is an array, or `x` otherwise.
Note that this method is not strictly necessary since Julia allows for scalar quantities
to be accessed at the first index like an array; however, the point is to make it
expliclty clear in scalar-typed code that a state variable is treated as such and is
not a vector valued quantity.
"""
getscalar(x::Number) = x
getscalar(x::Number, i) = x
......@@ -119,7 +124,7 @@ getscalar(x::AbstractArray, i) = x[i]
Convenience method for converting between `Dates.DateTime` and solver time.
"""
convert_t(t::DateTime) = Dates.datetime2epochms(t) / 1000
convert_t(t::Float64) = Dates.epochms2datetime(1000t)
convert_t(t::Float64) = Dates.epochms2datetime(floor(1000t))
"""
convert_tspan(tspan::Tuple{DateTime,DateTime})
convert_tspan(tspan::Tuple{Float64,Float64})
......@@ -139,16 +144,6 @@ 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))]...)))
"""
@generated genmap(f, args::T) where {T<:Tuple}
Generated `map` for `Tuple` types. This function is for use in generated functions where
generators/comprehensions like `map` are not allowed.
"""
@inline @generated function genmap(f, args::T) where {T<:Tuple}
return Expr(:tuple, (:(f(args[$i])) for i in 1:length(T.parameters))...)
end
"""
ffill!(x::AbstractVector{T}) where {E,T<:Union{Missing,E}}
......@@ -206,11 +201,15 @@ end
"""
ModelParameters.stripunits(obj)
Additional override for `stripunits` which reconstructs `obj` with all `AbstractQuantity` fields stripped of units.
Additional override for `stripunits` which reconstructs `obj` with all fields that have unitful quantity
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, values), Unitful.AbstractQuantity, Flatten.IGNORE)
return Flatten.reconstruct(obj, map(ustrip ∘ normalize_units, values), Unitful.AbstractQuantity, Flatten.IGNORE)
end
end
\ No newline at end of file
src/Physics/coupled.jl → src/coupling.jl
View file @ 738aa942
# Default empty implementations of diagnostic_step! and prognostic_step! for boundary layers
diagnosticstep!(::Top, ::BoundaryProcess, state) = nothing
prognosticstep!(::Top, ::BoundaryProcess, state) = nothing
diagnosticstep!(::Bottom, ::BoundaryProcess, state) = nothing
prognosticstep!(::Bottom, ::BoundaryProcess, state) = nothing
variables(ps::CoupledProcesses) = tuplejoin((variables(p) for p in ps.processes)...)
variables(ps::CoupledProcesses) = tuplejoin((variables(p) for p in ps.process)...)
variables(layer::Layer, ps::CoupledProcesses) = tuplejoin((variables(layer,p) for p in ps.processes)...)
callbacks(layer::Layer, ps::CoupledProcesses) = tuplejoin((callbacks(layer,p) for p in ps.processes)...)
events(layer::Layer, ps::CoupledProcesses) = tuplejoin((events(layer,p) for p in ps.processes)...)
"""
interact!(l1::Layer, ps1::CoupledProcesses{P1}, l2::Layer, ps2::CoupledProcesses{P2}, s1, s2) where {P1,P2}
Default implementation of `interact!` for multi-process (CoupledProcesses) types. Generates a specialized implementation that calls
Default implementation of `interact!` for coupled process (CoupledProcesses) types. Generates a specialized implementation that calls
`interact!` on all pairs of processes between the two layers. Since it is a generated function, the process matching
occurs at compile-time and the emitted code will simply be a sequence of `interact!` calls. Pairs of processes which
lack a definition of `interact!` should be automatically omitted by the compiler.
......@@ -17,94 +12,186 @@ lack a definition of `interact!` should be automatically omitted by the compiler
@generated function interact!(l1::Layer, ps1::CoupledProcesses{P1}, l2::Layer, ps2::CoupledProcesses{P2}, s1, s2) where {P1,P2}
p1types = Tuple(P1.parameters)
p2types = Tuple(P2.parameters)
crossprocesses = [(i,j) for (i,p1) in enumerate(p1types) for (j,p2) in enumerate(p2types)]
crossprocesses = [(i,j) for i in 1:length(p1types) for j in 1:length(p2types)]
expr = Expr(:block)
for (i,j) in crossprocesses
@>> quote
interact!(l1,ps1[$i],l2,ps2[$j],s1,s2)
end push!(expr.args)
quote
interact!(l1,ps1[$i],l2,ps2[$j],s1,s2)
end |> Base.Fix1(push!, expr.args)
end
return expr
end
@generated function interact!(l1::Layer, p1::Process, l2::Layer, ps2::CoupledProcesses{P2}, s1, s2) where {P2}
expr = Expr(:block)
for i in 1:length(P2.parameters)
quote
interact!(l1,p1,l2,ps2[$i],s1,s2)
end |> Base.Fix1(push!, expr.args)
end
return expr
end
@generated function interact!(l1::Layer, ps1::CoupledProcesses{P1}, l2::Layer, p2::Process, s1, s2) where {P1}
expr = Expr(:block)
for i in 1:length(P1.parameters)
quote
interact!(l1,ps1[$i],l2,p2,s1,s2)
end |> Base.Fix1(push!, expr.args)
end
return expr
end
"""
diagnosticstep!(l::Layer, ps::CoupledProcesses{P}, state) where {P}
Default implementation of `diagnosticstep!` for multi-process types. Calls each process in sequence.
Default implementation of `diagnosticstep!` for coupled process types. Calls each process in sequence.
"""
@generated function diagnosticstep!(l::Layer, ps::CoupledProcesses{P}, state) where {P}
expr = Expr(:block)
for i in 1:length(P.parameters)
@>> quote
diagnosticstep!(l,ps[$i],state)
end push!(expr.args)
quote
diagnosticstep!(l,ps[$i],state)
end |> Base.Fix1(push!, expr.args)
end
return expr
end
"""
prognosticstep!(l::Layer, ps::CoupledProcesses{P}, state) where {P}
Default implementation of `prognosticstep!` for multi-process types. Calls each process in sequence.
Default implementation of `prognosticstep!` for coupled process types. Calls each process in sequence.
"""
@generated function prognosticstep!(l::Layer, ps::CoupledProcesses{P}, state) where {P}
expr = Expr(:block)
for i in 1:length(P.parameters)
@>> quote
prognosticstep!(l,ps[$i],state)
end push!(expr.args)
quote
prognosticstep!(l,ps[$i],state)
end |> Base.Fix1(push!, expr.args)
end
return expr
end
"""
initialcondition!([l::Layer,] ps::CoupledProcesses{P}, state) where {P}
Default implementation of `initialcondition!` for multi-process types. Calls each process in sequence.
Default implementation of `initialcondition!` for coupled process types. Calls each process in sequence.
"""
@generated function initialcondition!(ps::CoupledProcesses{P}, state) where {P}
expr = Expr(:block)
for i in 1:length(P.parameters)
@>> quote
initialcondition!(ps[$i],state)
end push!(expr.args)
quote
initialcondition!(ps[$i],state)
end |> Base.Fix1(push!, expr.args)
end
return expr
end
@generated function initialcondition!(l::Layer, ps::CoupledProcesses{P}, state) where {P}
expr = Expr(:block)
for i in 1:length(P.parameters)
@>> quote
initialcondition!(l,ps[$i],state)
end push!(expr.args)
quote
initialcondition!(l,ps[$i],state)
end |> Base.Fix1(push!, expr.args)
end
return expr
end
"""
initialcondition!(l1::Layer, ps1::CoupledProcesses{P1}, l2::Layer, ps2::CoupledProcesses{P2}, s1, s2) where {P1,P2}
Default implementation of `initialcondition!` for multi-process types. Calls each process in sequence.
Default implementation of `initialcondition!` for coupled process types. Calls each process in sequence.
"""
@generated function initialcondition!(l1::Layer, ps1::CoupledProcesses{P1}, l2::Layer, ps2::CoupledProcesses{P2}, s1, s2) where {P1,P2}
p1types = Tuple(P1.parameters)
p2types = Tuple(P2.parameters)
crossprocesses = [(i,j) for (i,p1) in enumerate(p1types) for (j,p2) in enumerate(p2types)]
crossprocesses = [(i,j) for i in 1:length(p1types) for j in 1:length(p2types)]
expr = Expr(:block)
for (i,j) in crossprocesses
@>> quote
initialcondition!(l1,ps1[$i],l2,ps2[$j],s1,s2)
end push!(expr.args)
quote
initialcondition!(l1,ps1[$i],l2,ps2[$j],s1,s2)
end |> Base.Fix1(push!, expr.args)
end
return expr
end
@generated function initialcondition!(l1::Layer, p1::Process, l2::Layer, ps2::CoupledProcesses{P2}, s1, s2) where {P2}
expr = Expr(:block)
for i in 1:length(P2.parameters)
quote
initialcondition!(l1,p1,l2,ps2[$i],s1,s2)
end |> Base.Fix1(push!, expr.args)
end
return expr
end
@generated function initialcondition!(l1::Layer, ps1::CoupledProcesses{P1}, l2::Layer, p2::Process, s1, s2) where {P1}
expr = Expr(:block)
for i in 1:length(P1.parameters)
quote
initialcondition!(l1,ps1[$i],l2,p2,s1,s2)
end |> Base.Fix1(push!, expr.args)
end
return expr
end
@generated function criterion(ev::ContinuousEvent, l::Layer, ps::CoupledProcesses{P}, state) where {name,P}
expr = Expr(:block)
push!(expr.args, :(value = 1.0))
for i in 1:length(P.parameters)
quote
value *= criterion(ev,l,ps[$i],state)
end |> Base.Fix1(push!, expr.args)
end
push!(expr.args, :(return value))
return expr
end
@generated function criterion(ev::DiscreteEvent, l::Layer, ps::CoupledProcesses{P}, state) where {name,P}
expr = Expr(:block)
push!(expr.args, :(value = true))
for i in 1:length(P.parameters)
quote
value *= criterion(ev,l,ps[$i],state)
end |> Base.Fix1(push!, expr.args)
end
push!(expr.args, :(return value))
return expr
end
@generated function trigger!(ev::Event, l::Layer, ps::CoupledProcesses{P}, state) where {name,P}
expr = Expr(:block)
for i in 1:length(P.parameters)
quote
trigger!(ev,l,ps[$i],state)
end |> Base.Fix1(push!, expr.args)
end
return expr
end
@generated function trigger!(ev::ContinuousEvent, tr::ContinuousTrigger, l::Layer, ps::CoupledProcesses{P}, state) where {name,P}
expr = Expr(:block)
for i in 1:length(P.parameters)
quote
trigger!(ev,tr,l,ps[$i],state)
end |> Base.Fix1(push!, expr.args)
end
return expr
end
"""
timestep(l::Layer, ps::CoupledProcesses{P}, state) where {P}
Default implementation of `timestep` for coupled process types. Calls each process in sequence.
"""
@generated function timestep(l::Layer, ps::CoupledProcesses{P}, state) where {P}
expr = Expr(:block)
push!(expr.args, :(dtmax = Inf))
for i in 1:length(P.parameters)
quote
dtmax = min(dtmax, timestep(l,ps[$i],state))
end |> Base.Fix1(push!, expr.args)
end
push!(expr.args, :(return dtmax))
return expr
end
"""
observe(::Val{name}, l::Layer, ps::CoupledProcesses{P}, state) where {P}
Default implementation of `observe` for multi-process types. Calls each process in sequence.
Default implementation of `observe` for coupled process types. Calls each process in sequence.
"""
@generated function observe(val::Val{name}, l::Layer, ps::CoupledProcesses{P}, state) where {name,P}
expr = Expr(:block)
for i in 1:length(P.parameters)
@>> quote
observe(val,l,ps[$i],state)
end push!(expr.args)
quote
observe(val,l,ps[$i],state)
end |> Base.Fix1(push!, expr.args)
end
return expr
end
\ No newline at end of file
src/methods.jl
View file @ 738aa942
......@@ -35,6 +35,8 @@ Defines the diagnostic update for a Process on a given Layer.
"""
diagnosticstep!(l::Layer, p::Process, state) = nothing
diagnosticstep!(l::Layer, state) = nothing
diagnosticstep!(::Top, ::BoundaryProcess, state) = nothing
diagnosticstep!(::Bottom, ::BoundaryProcess, state) = nothing
"""
prognosticstep!(l::Layer, p::Process, state)
......@@ -42,6 +44,8 @@ Defines the prognostic update for a Process on a given layer. Note that an insta
for all non-boundary (subsurface) processes/layers.
"""
prognosticstep!(l::Layer, p::Process, state) = error("no prognostic step defined for $(typeof(l)) with $(typeof(p))")
prognosticstep!(::Top, ::BoundaryProcess, state) = nothing
prognosticstep!(::Bottom, ::BoundaryProcess, state) = nothing
"""
interact!(::Layer, ::Process, ::Layer, ::Process, state1, state2)
......@@ -50,6 +54,14 @@ follows decreasing depth, i.e. the first layer/process is always on top of the s
and separate dispatches must be provided for interactions in reverse order.
"""
interact!(::Layer, ::Process, ::Layer, ::Process, state1, state2) = nothing
"""
timestep(::Layer, ::Process, state)
Retrieves the recommended timestep for the given `Process` defined on the given `Layer`.
The default implementation returns `Inf` which indicates no timestep restriction. The
actual chosen timestep will depend on the integrator being used and other user configuration options.
"""
timestep(::Layer, ::Process, state) = Inf
"""
observe(::Val{name}, ::Layer, ::Process, state1)
......@@ -101,44 +113,51 @@ where `BP` is a `BoundaryProcess` that provides the boundary conditions.
"""
BoundaryStyle(::Type{BP}) where {BP<:BoundaryProcess} = error("No style specified for boundary process $BP")
BoundaryStyle(bc::BoundaryProcess) = BoundaryStyle(typeof(bc))
# Callbacks
"""
callbacks(::Layer, ::Process)
Defines callbacks for a given Process on the given Layer. Implementations should return a `Tuple` or `Callback`s.
"""
callbacks(::Layer, ::Process) = ()
# Events
"""
criterion(c::Callback, ::Layer, ::Process, state)
events(::Layer, ::Process)
Callback criterion/condition. Should return a `Bool` for discrete callbacks and a real number for continuous callbacks.
Defines "events" for a given Process on the given Layer. Implementations should return a `Tuple` of `Event`s.
"""
criterion(c::Callback, ::Layer, ::Process, state) = error("missing implementation of criterion for $(typeof(c))")
events(::Layer, ::Process) = ()
"""
affect!(c::Callback, ::Layer, ::Process, state)
criterion(::Event, ::Layer, ::Process, state)
Callback action executed when `criterion` is met (boolean condition for discrete callbacks, zero for continuous callbacks).
Event criterion/condition. Should return a `Bool` for discrete events and a real number for continuous events.
"""
affect!(c::Callback, ::Layer, ::Process, state) = error("missing implementation of affect! for $(typeof(c))")
criterion(ev::Event, ::Layer, ::Process, state) = true
"""
CallbackStyle(::C)
CallbackStyle(::Type{<:Callback})
trigger!(::Event, ::Layer, ::Process, state)
trigger!(ev::ContinuousEvent, ::ContinuousTrigger, ::Layer, ::Process, state) where {name}
Trait implementation that defines the "style" or type of the given callback as any subtype of `CallbackStyle`,
for example `Discrete` or `Continuous`.
Event action executed when `criterion` is met.
"""
CallbackStyle(::C) where {C<:Callback} = CallbackStyle(C)
CallbackStyle(::Type{<:Callback}) = Discrete()
trigger!(ev::Event, ::Layer, ::Process, state) = nothing
trigger!(ev::ContinuousEvent, ::ContinuousTrigger, ::Layer, ::Process, state) = nothing
# Discretization
"""
midpoints(::Layer, state)
midpoint(::Layer, state)
midpoint(::Layer, state, i)
midpoint(::Layer, state, ::typeof(first))
midpoint(::Layer, state, ::typeof(last))
Get midpoint(s) of layer or all grid cells within the layer.
Get midpoint (m) of layer or grid cell `i`.
"""
midpoints(::Layer, state) = cells(state.grid)
midpoint(::Layer, state) = (state.grid[end] - state.grid[1])/2
midpoint(::Layer, state, i) = cells(state.grid)[i]
midpoint(l::Layer, state, ::typeof(first)) = midpoint(l, state, 1)
midpoint(l::Layer, state, ::typeof(last)) = midpoint(l, state, lastindex(state.grid)-1)
midpoint(::Union{Top,Bottom}, state) = Inf
"""
thickness(::Layer, state)
thickness(::Layer, state, i)
thickness(l::Layer, state, ::typeof(first))
thickness(l::Layer, state, ::typeof(last))
Get thickness (m) of layer or all grid cells within the layer.
Get thickness (m) of layer or grid cell `i`.
"""
thickness(::Layer, state) = Δ(state.grid)
thickness(::Layer, state) = state.grid[end] - state.grid[1]
thickness(::Layer, state, i) = Δ(state.grid)[i]
thickness(l::Layer, state, ::typeof(first)) = thickness(l, state, 1)
thickness(l::Layer, state, ::typeof(last)) = thickness(l, state, lastindex(state.grid)-1)
thickness(::Union{Top,Bottom}, state) = Inf
src/types.jl
View file @ 738aa942
......@@ -27,6 +27,13 @@ Base.Broadcast.broadcastable(l::Layer) = Ref(l)
Abstract base type for all dynamical processes.
"""
abstract type Process end
"""
isless(::Process, ::Process)
Defines an ordering between processes, which can be useful for automatically enforcing
order constraints when coupling processes together.
"""
Base.isless(::Process, ::Process) = false
"""
SubSurfaceProcess <: Process
......@@ -51,28 +58,30 @@ of other processes.
"""
struct CoupledProcesses{TProcs} <: Process
processes::TProcs
CoupledProcesses(processes::SubSurfaceProcess...) = CoupledProcesses(processes)
CoupledProcesses(processes::BoundaryProcess...) = CoupledProcesses(processes)
CoupledProcesses(processes::Tuple{Vararg{Process}}) = new{typeof(processes)}(processes)
CoupledProcesses(processes::SubSurfaceProcess...) = new{typeof(processes)}(processes)
CoupledProcesses(processes::BoundaryProcess...) = new{typeof(processes)}(processes)
end
Base.iterate(cp::CoupledProcesses) = Base.iterate(cp.processes)
Base.iterate(cp::CoupledProcesses, state) = Base.iterate(cp.processes, state)
"""
Coupled{P1,P2} = CoupledProcesses{Tuple{T1,T2}} where {T1,T2}
Coupled2{P1,P2} = CoupledProcesses{Tuple{T1,T2}} where {T1,T2}
Represents an explicitly coupled pair processes. Alias for `CoupledProcesses{Tuple{P1,P2}}`.
Represents an explicitly coupled pair of processes. Alias for `CoupledProcesses{Tuple{P1,P2}}`.
`Coupled` provides a simple mechanism for defining new behaviors on multi-processes systems.
"""
const Coupled{P1,P2} = CoupledProcesses{Tuple{T1,T2}} where {T1,T2}
const Coupled2{P1,P2} = CoupledProcesses{Tuple{T1,T2}} where {T1,T2}
const Coupled3{P1,P2,P3} = CoupledProcesses{Tuple{T1,T2,T3}} where {T1,T2,T3}
const Coupled4{P1,P2,P3,P4} = CoupledProcesses{Tuple{T1,T2,T3,T4}} where {T1,T2,T3,T4}
"""
Coupled(p1,p2)
Coupled(ps::Process...)
Alias for `CoupledProcesses(p1,p2)`.
Alias for `CoupledProcesses(ps...)`.
"""
Coupled(p1::P1, p2::P2) where {P1<:Process,P2<:Process} = CoupledProcesses(p1,p2)
Coupled(ps::Process...) = CoupledProcesses(ps...)
# Base methods
Base.show(io::IO, ps::CoupledProcesses{T}) where T = print(io, "$T")
@propagate_inbounds @inline Base.getindex(ps::CoupledProcesses, i) = ps.processes[i]
Base.show(io::IO, ::CoupledProcesses{T}) where T = print(io, "Coupled($(join(T.parameters, " with ")))")
Base.iterate(cp::CoupledProcesses) = Base.iterate(cp.processes)
Base.iterate(cp::CoupledProcesses, state) = Base.iterate(cp.processes, state)
@propagate_inbounds Base.getindex(cp::CoupledProcesses, i) = cp.processes[i]
# allow broadcasting of Process types
Base.Broadcast.broadcastable(p::Process) = Ref(p)
......@@ -92,19 +101,28 @@ struct Dirichlet <: BoundaryStyle end
"""
struct Neumann <: BoundaryStyle end
"""
Callback
Event{name}
Base type for callback implementations.
Base type for integration "events" (i.e. callbacks). `name` should be a `Symbol`
or type which identifies the event for the purposes of dispatch on `criterion` and `trigger!`.
"""
abstract type Callback end
abstract type Event{name} end
struct DiscreteEvent{name} <: Event{name}
DiscreteEvent(name::Symbol) = new{name}()
end
struct ContinuousEvent{name} <: Event{name}
ContinuousEvent(name::Symbol) = new{name}()
end
"""
CallbackStyle
ContinuousTrigger
Trait for callback types.
Base type for continuous trigger flags, `Increasing` and `Decreasing`, which indicate
an upcrossing of the function root (negative to positive) and a downcrossing (positive to negative)
respectively.
"""
abstract type CallbackStyle end
struct Discrete <: CallbackStyle end
struct Continuous <: CallbackStyle end
abstract type ContinuousTrigger end
struct Increasing <: ContinuousTrigger end
struct Decreasing <: ContinuousTrigger end
"""
Parameterization
......
test/Physics/HeatConduction/sfcc_bench.jl
View file @ 738aa942
......@@ -14,32 +14,32 @@ function benchmarksfcc()
sfcc = SFCC(f, solver)
soil = Soil() |> stripparams |> stripunits
heat = Heat(freezecurve=sfcc) |> stripparams |> stripunits
L = heat.L
L = heat.prop.L
Lf = heat.prop.Lf
# set up multi-grid-cell state vars
T = [-5.0 for i in 1:10]
θw = Soils.soilcomponent(Val{:θw}(), soil.para)
θ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)
θl = f.(T,Tₘ,θres,θp,θw,Lf,α,n) # set liquid water content according to freeze curve
C = heatcapacity.(soil,heat,θw,θl,θm,θo)
θ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,
θl = f.(T,Tₘ,θres,θp,θw,Lf,α,n),
C = heatcapacity.(soil,heat,θw,θl,θm,θo);
enthalpy.(T,C,L,θl) # compute enthalpy at "true" temperature
θ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,θl=θl,θw=θw)
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, $θw, $θm, $θo)
# @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, θw, θm, θo, 0.0)
@btime Soils.sfccsolve($solver, $soil, $heat, $sfcc.f, $sfcc.∇f, $params, $H[1], $L, $θw, $θm, $θo, 0.0)
# @btime Soils.residual($soil, $heat, $T[1], $H[1], $L, $sfcc.f, $params, $θw, $θm, $θo)
# res = @btime Soils.sfccsolve($solver, $soil, $heat, $sfcc.f, $sfcc.∇f, $params, $H[1], $L, $θw, $θm, $θo)
# 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.θl[1])")
println("\nsolution: $(state.T[1]), $(state.θw[1])")
end
test/Physics/HeatConduction/sfcc_tests.jl
View file @ 738aa942
......@@ -8,7 +8,7 @@ using Test
Tₘ = 0.0
θres = 0.0
soil = @pstrip Soil(para=Soils.CharacteristicFractions())
θw = Soils.soilcomponent(Val{:θw}(), soil.para)
θ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)
......@@ -20,8 +20,8 @@ using Test
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
θl = f(-0.1u"°C",θsat,θsat,θres,Tₘ,γ)
@test θl > 0.0 && θl < 1.0
θw = f(-0.1u"°C",θsat,θsat,θres,Tₘ,γ)
@test θw > 0.0 && θw < 1.0
end
end
@testset "Newton solver checks" begin
......@@ -31,47 +31,48 @@ using Test
f = @pstrip McKenzie()
sfcc = SFCC(f, SFCCNewtonSolver(abstol=abstol, reltol=reltol))
heat = @pstrip Heat(freezecurve=sfcc)
L = heat.L
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]
θl = f.(T,θp,θw,θres,Tₘ,γ) # set liquid water content according to freeze curve
C = heatcapacity.(soil,heat,θw,θl,θm,θo)
θw = f.(T,θp,θwi,θres,Tₘ,γ) # set liquid water content according to freeze curve
C = heatcap.(θw)
H = let T = T.+1.0,
θl = f.(T,θp,θw,θres,Tₘ,γ)
C = heatcapacity.(soil,heat,θw,θl,θm,θo);
enthalpy.(T,C,L,θl) # compute enthalpy at "true" temperature
θ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,θl=θl,)
state = (T=T,C=C,dHdT=similar(C),H=H,θw=θw,)
freezethaw!(soil, heat, state)
@test all(abs.(T.-(H .- L.*θl)./C) .<= abstol)
@test all(abs.(T.-(H .- L.*θw)./C) .<= abstol)
end
@testset "Right tail" begin
# set up single-grid-cell state vars
T = [5.0]
θl = f.(T,θp,θw,θres,Tₘ,γ) # set liquid water content according to freeze curve
C = heatcapacity.(soil,heat,θw,θl,θm,θo)
θw = f.(T,θp,θwi,θres,Tₘ,γ) # set liquid water content according to freeze curve
C = heatcap.(θw)
H = let T = T.-1.0,
θl = f.(T,θp,θw,θres,Tₘ,γ)
C = heatcapacity.(soil,heat,θw,θl,θm,θo);
enthalpy.(T,C,L,θl) # compute enthalpy at "true" temperature
θ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,θl=θl,)
state = (T=T,C=C,dHdT=similar(C),H=H,θw=θw,)
freezethaw!(soil, heat, state)
@test all(abs.(T.-(H .- L.*θl)./C) .<= abstol)
@test all(abs.(T.-(H .- L.*θw)./C) .<= abstol)
end
@testset "Near zero" begin
# set up single-grid-cell state vars
T = [-0.05]
θl = f.(T,θp,θw,θres,Tₘ,γ) # set liquid water content according to freeze curve
C = heatcapacity.(soil,heat,θw,θl,θm,θo)
θw = f.(T,θp,θwi,θres,Tₘ,γ) # set liquid water content according to freeze curve
C = heatcap.(θw)
H = let T = T.+0.04,
θl = f.(T,θp,θw,θres,Tₘ,γ)
C = heatcapacity.(soil,heat,θw,θl,θm,θo);
enthalpy.(T,C,L,θl) # compute enthalpy at "true" temperature
θ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,θl=θl,)
state = (T=T,C=C,dHdT=similar(C),H=H,θw=θw,)
freezethaw!(soil, heat, state)
@test all(abs.(T.-(H .- L.*θl)./C) .<= abstol)
@test all(abs.(T.-(H .- L.*θw)./C) .<= abstol)
end
end
end
......@@ -85,8 +86,8 @@ using Test
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
θl = f(-0.1u"°C",θtot,θtot,θres,Tₘ,δ)
@test θl > 0.0 && θl < 1.0
θw = f(-0.1u"°C",θtot,θtot,θres,Tₘ,δ)
@test θw > 0.0 && θw < 1.0
end
end
@testset "Newton solver checks" begin
......@@ -96,45 +97,46 @@ using Test
f = @pstrip Westermann()
sfcc = SFCC(f, SFCCNewtonSolver(abstol=abstol, reltol=reltol))
heat = @pstrip Heat(freezecurve=sfcc)
L = heat.L
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]
θl = f.(T,θp,θw,θres,Tₘ,δ) # set liquid water content according to freeze curve
C = heatcapacity.(soil,heat,θw,θl,θm,θo)
θw = f.(T,θp,θwi,θres,Tₘ,δ) # set liquid water content according to freeze curve
C = heatcap.(θw)
H = let T = T.+1.0,
θl = f.(T,θp,θw,θres,Tₘ,δ)
C = heatcapacity.(soil,heat,θw,θl,θm,θo);
enthalpy.(T,C,L,θl) # compute enthalpy at "true" temperature
θ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,θl=θl,)
state = (T=T,C=C,dHdT=similar(C),H=H,θw=θw,)
freezethaw!(soil, heat, state)
@test all(abs.(T.-(H .- L.*θl)./C) .<= abstol)
@test all(abs.(T.-(H .- L.*θw)./C) .<= abstol)
end
@testset "Right tail" begin
T = [5.0]
θl = f.(T,θp,θw,θres,Tₘ,δ) # set liquid water content according to freeze curve
C = heatcapacity.(soil,heat,θw,θl,θm,θo)
θw = f.(T,θp,θwi,θres,Tₘ,δ) # set liquid water content according to freeze curve
C = heatcap.(θw)
H = let T = T.-1,
θl = f.(T,θp,θw,θres,Tₘ,δ)
C = heatcapacity.(soil,heat,θw,θl,θm,θo);
enthalpy.(T,C,L,θl) # compute enthalpy at "true" temperature
θ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,θl=θl,)
state = (T=T,C=C,dHdT=similar(C),H=H,θw=θw,)
freezethaw!(soil, heat, state)
@test all(abs.(T.-(H .- L.*θl)./C) .<= abstol)
@test all(abs.(T.-(H .- L.*θw)./C) .<= abstol)
end
@testset "Near zero" begin
T = [-0.05]
θl = f.(T,θp,θw,θres,Tₘ,δ) # set liquid water content according to freeze curve
C = heatcapacity.(soil,heat,θw,θl,θm,θo)
θw = f.(T,θp,θwi,θres,Tₘ,δ) # set liquid water content according to freeze curve
C = heatcap.(θw)
H = let T = T.+0.04,
θl = f.(T,θp,θw,θres,Tₘ,δ)
C = heatcapacity.(soil,heat,θw,θl,θm,θo);
enthalpy.(T,C,L,θl) # compute enthalpy at "true" temperature
θ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,θl=θl,)
state = (T=T,C=C,dHdT=similar(C),H=H,θw=θw,)
freezethaw!(soil, heat, state)
@test all(abs.(T.-(H .- L.*θl)./C) .<= abstol)
@test all(abs.(T.-(H .- L.*θw)./C) .<= abstol)
end
end
end
......@@ -148,8 +150,8 @@ using Test
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
θl = f(-0.1u"°C",θsat,θsat,Lf,θres,Tₘ,α,n)
@test θl > 0.0 && θl < 1.0
θ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
......@@ -162,46 +164,47 @@ using Test
f = @pstrip DallAmico()
sfcc = SFCC(f, SFCCNewtonSolver(abstol=abstol, reltol=reltol))
heat = @pstrip Heat(freezecurve=sfcc)
L = heat.L
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]
θl = f.(T,θp,θw,Lf,θres,Tₘ,α,n) # set liquid water content according to freeze curve
C = heatcapacity.(soil,heat,θw,θl,θm,θo)
θ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,
θl = f.(T,θp,θw,Lf,θres,Tₘ,α,n)
C = heatcapacity.(soil,heat,θw,θl,θm,θo);
enthalpy.(T,C,L,θl) # compute enthalpy at "true" temperature
θ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,θl=θl,)
state = (T=T,C=C,dHdT=similar(C),H=H,θw=θw,)
@inferred freezethaw!(soil, heat, state)
@test all(abs.(T.-(H .- L.*θl)./C) .<= abstol)
@test all(abs.(T.-(H .- L.*θw)./C) .<= abstol)
end
@testset "Right tail" begin
T = [5.0]
θl = f.(T,θp,θw,Lf,θres,Tₘ,α,n) # set liquid water content according to freeze curve
C = heatcapacity.(soil,heat,θw,θl,θm,θo)
θ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,
θl = f.(T,θp,θw,Lf,θres,Tₘ,α,n)
C = heatcapacity.(soil,heat,θw,θl,θm,θo);
enthalpy.(T,C,L,θl) # compute enthalpy at "true" temperature
θ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,θl=θl,)
state = (T=T,C=C,dHdT=similar(C),H=H,θw=θw,)
@inferred freezethaw!(soil, heat, state)
@test all(abs.(T.-(H .- L.*θl)./C) .<= abstol)
@test all(abs.(T.-(H .- L.*θw)./C) .<= abstol)
end
@testset "Near zero" begin
T = [-0.05]
θl = f.(T,θp,θw,Lf,θres,Tₘ,α,n) # set liquid water content according to freeze curve
C = heatcapacity.(soil,heat,θw,θl,θm,θo)
θ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,
θl = f.(T,θp,θw,Lf,θres,Tₘ,α,n)
C = heatcapacity.(soil,heat,θw,θl,θm,θo);
enthalpy.(T,C,L,θl) # compute enthalpy at "true" temperature
θ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,θl=θl,)
state = (T=T,C=C,dHdT=similar(C),H=H,θw=θw,)
@inferred freezethaw!(soil, heat, state)
@test all(abs.(T.-(H .- L.*θl)./C) .<= abstol)
@test all(abs.(T.-(H .- L.*θw)./C) .<= abstol)
end
end
end
......@@ -214,21 +217,22 @@ using Test
f = @pstrip McKenzie()
sfcc = SFCC(f, SFCCNewtonSolver())
heat = @pstrip Heat(freezecurve=sfcc)
L = heat.L
heatcap = θw -> heatcapacity(soil, heat, θw, θwi - θw, 1-θwi-θm-θo, θm, θo)
L = heat.prop.L
T = [-0.1]
θl = f.(T,θp,θw,θres,Tₘ,γ) # set liquid water content according to freeze curve
C = heatcapacity.(soil,heat,θw,θl,θm,θo)
H = enthalpy.(T.+0.09,C,L,θl) # compute enthalpy at +1 degree
θ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,θw,θres,Tₘ,p.γ)), p)
∂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))
θl = similar(θl,eltype(p))
state = (T=T_,C=C,dHdT=similar(C),H=H,θl=θl,)
θ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]
......
test/Strat/stratigraphy_tests.jl
View file @ 738aa942
......@@ -18,15 +18,15 @@ include("../types.jl")
if i == 1
@test CryoGrid.componentname(component) == :top
@test typeof(component.layer) <: Top
@test typeof(component.processes) <: CoupledProcesses{Tuple{TestBoundary}}
@test typeof(component.process) <: TestBoundary
elseif i == 2
@test CryoGrid.componentname(component) == :testground
@test typeof(component.layer) <: TestGroundLayer
@test typeof(component.processes) <: CoupledProcesses{Tuple{TestGroundProcess}}
@test typeof(component.process) <: TestGroundProcess
elseif i == 3
@test CryoGrid.componentname(component) == :bottom
@test typeof(component.layer) <: Bottom
@test typeof(component.processes) <: CoupledProcesses{Tuple{TestBoundary}}
@test typeof(component.process) <: TestBoundary
end
end
end
......@@ -46,19 +46,51 @@ include("../types.jl")
if i == 1
@test CryoGrid.componentname(component) == :top
@test typeof(component.layer) <: Top
@test typeof(component.processes) <: CoupledProcesses{Tuple{TestBoundary}}
@test typeof(component.process) <: TestBoundary
elseif i == 2
@test CryoGrid.componentname(component) == :testground1
@test typeof(component.layer) <: TestGroundLayer
@test typeof(component.processes) <: CoupledProcesses{Tuple{TestGroundProcess}}
@test typeof(component.process) <: TestGroundProcess
elseif i == 3
@test CryoGrid.componentname(component) == :testground2
@test typeof(component.layer) <: TestGroundLayer
@test typeof(component.processes) <: CoupledProcesses{Tuple{TestGroundProcess}}
@test typeof(component.process) <: TestGroundProcess
elseif i == 4
@test CryoGrid.componentname(component) == :bottom
@test typeof(component.layer) <: Bottom
@test typeof(component.processes) <: CoupledProcesses{Tuple{TestBoundary}}
@test typeof(component.process) <: TestBoundary
end
end
end
@testset "4-layer w/ coupling" begin
bounds = (-1.0u"m",0.0u"m",10.0u"m",100.0u"m")
strat = Stratigraphy(
bounds[1] => top(TestBoundary(), TestBoundary()),(
bounds[2] => subsurface(:testground1, TestGroundLayer(), TestGroundProcess(), TestGroundProcess()),
bounds[3] => subsurface(:testground2, TestGroundLayer(), TestGroundProcess())
),
bounds[4] => bottom(TestBoundary(), TestBoundary())
)
# Check boundaries
@test all([bounds[i] == b for (i,b) in enumerate(bounds)])
# Check iteration and component types
for (i,component) in enumerate(strat)
if i == 1
@test CryoGrid.componentname(component) == :top
@test typeof(component.layer) <: Top
@test typeof(component.process) <: CoupledProcesses{Tuple{TestBoundary,TestBoundary}}
elseif i == 2
@test CryoGrid.componentname(component) == :testground1
@test typeof(component.layer) <: TestGroundLayer
@test typeof(component.process) <: CoupledProcesses{Tuple{TestGroundProcess,TestGroundProcess}}
elseif i == 3
@test CryoGrid.componentname(component) == :testground2
@test typeof(component.layer) <: TestGroundLayer
@test typeof(component.process) <: TestGroundProcess
elseif i == 4
@test CryoGrid.componentname(component) == :bottom
@test typeof(component.layer) <: Bottom
@test typeof(component.process) <: CoupledProcesses{Tuple{TestBoundary,TestBoundary}}
end
end
end
......