src/Drivers/DiffEq/solvers.jl

+import LinearSolve
+
+struct TDMASolver <: SciMLBase.AbstractLinearAlgorithm end
+
+LinearSolve.init_cacheval(alg::TDMASolver, A, b, u, Pl, Pr, maxiters, abstol, reltol, verbose) = Tridiagonal(similar(A))
+LinearSolve.needs_concrete_A(alg::TDMASolver) = true
+
+function SciMLBase.solve(cache::LinearSolve.LinearCache, alg::TDMASolver; kwargs...)
+    factorized = !cache.isfresh
+    if !factorized
+        copyto!(cache.cacheval, cache.A)
+    end
+    a = diag(cache.cacheval,-1)
+    b = diag(cache.cacheval,0)
+    c = diag(cache.cacheval,1)
+    d = cache.b
+    x = cache.u
+    Numerics.tdma_solve!(x, a, b, c, d, factorized)
+    SciMLBase.build_linear_solution(alg,x,nothing,cache)
+end

src/Drivers/DiffEq/steplimiters.jl

+struct CFLStepLimiter{TTile,A}
+    tile::TTile
+    Δt::A
+    default_dt::Float64
+end
+function (cfl::CFLStepLimiter{<:Tile,<:AbstractArray})(u,p,t)
+    let Δt = cfl.Δt,
+        cfl! = Drivers.cfl!(Heat),
+        Δx = dustrip(Δ(cfl.tile.grid)),
+        dHdT = getvar(:dHdT, cfl.tile, u), # apparent heat capacity
+        kc = getvar(:kc, cfl.tile, u); # thermal cond. at grid centers
+        cfl!(Δt, Δx, dHdT, kc)
+        Δt_min = minimum(Δt)
+        IfElse.ifelse(isfinite(Δt_min) && Δt_min > 0, Δt_min, cfl.default_dt)
+    end
+end
+function (cfl::CFLStepLimiter{<:Tile,Nothing})(u,p,t)
+    let cfl = Drivers.cfl(Heat),
+        Δx = dustrip(Δ(cfl.tile.grid)),
+        dHdT = getvar(:dHdT, cfl.tile, u), # apparent heat capacity
+        kc = getvar(:kc, cfl.tile, u); # thermal cond. at grid centers
+        Δt = cfl(Δt, Δx, dHdT, kc)
+        Δt_min = minimum(Δt)
+        IfElse.ifelse(isfinite(Δt_min) && Δt_min > 0, Δt_min, cfl.default_dt)
+    end
+end
+function CFLStepLimiter(tile::HeatOnlyTile; courant_number=1//2, default_dt=60.0, iip=true)
+    cfl = iip ? CFLStepLimiter(tile, zero(dustrip(Δ(tile.grid))), default_dt) : CFLStepLimiter(tile, nothing, default_dt)
+    StepsizeLimiter(cfl; safety_factor=courant_number, max_step=true)
+end

src/Drivers/Drivers.jl

 """
-Module containing types and functions for enabling time-integration of a CryoGrid land model.
+Module containing types and functions for enabling time-integration of a CryoGrid `Tile` model.
 """
 module Drivers
 
-using CryoGrid: Land, SubSurface, CompoundProcess, variables
-using CryoGrid.InputOutput
-using CryoGrid.Numerics
-using CryoGrid.Physics: Heat
+using CryoGrid
 using CryoGrid.Utils
 
-using ComponentArrays
-using Dates
-using DimensionalData
-using Flatten
-using ModelParameters
-using LinearAlgebra
 using Reexport
-using Unitful
 
-import CryoGrid.Land: LandModel, Stratigraphy, StratComponent
-
-export DefaultJac, TridiagJac, JacobianStyle, HeatOnlyLandModel
+export JacobianStyle, DefaultJac, TridiagJac, HeatOnlyTile
 
 """
     JacobianStyle
@@ -31,19 +19,31 @@ abstract type JacobianStyle end
 struct DefaultJac <: JacobianStyle end
 struct TridiagJac <: JacobianStyle end
 """
-    JacobianStyle(::Type{<:LandModel})
+    JacobianStyle(::Type{<:Tile})
 
-Can be overriden/extended to specify Jacobian structure for specific `LandModel`s.
+Can be overriden/extended to specify Jacobian structure for specific `Tile`s.
 """
-JacobianStyle(::Type{<:LandModel}) = DefaultJac()
+JacobianStyle(::Type{<:Tile}) = DefaultJac()
 # Auto-detect Jacobian sparsity for problems with one or more heat-only layers.
 # Note: This assumes that the processes/forcings on the boundary layers do not violate the tridiagonal structure!
 # Unfortunately, the Stratigraphy type signature is a bit nasty to work with :(
-const HeatOnlyLandModel = LandModel{<:Stratigraphy{N,<:Tuple{TTop,Vararg{<:Union{<:StratComponent{<:SubSurface, <:CompoundProcess{<:Tuple{<:Heat}}},TBot}}}}} where {N,TTop,TBot}
-JacobianStyle(::Type{<:HeatOnlyLandModel}) = TridiagJac()
+const HeatOnlyTile = Tile{<:Stratigraphy{N,<:Tuple{TTop,Vararg{<:Union{<:StratComponent{<:SubSurface, <:CoupledProcesses{<:Tuple{<:Heat}}},TBot}}}}} where {N,TTop,TBot}
+JacobianStyle(::Type{<:HeatOnlyTile}) = TridiagJac()
+
+# CFL conditions
+"""
+    cfl(::Type{<:SubSurfaceProcess})
+    cfl!(::Type{<:SubSurfaceProcess})
+
+Returns a function of the form (Δx, args...) -> Δt (or in-place, (Δt, Δx, args...) -> Δt)
+which comptues the CFL condition with process-specific parameters `args`.
+"""
+cfl(::Type{<:SubSurfaceProcess}) = error("not implemented")
+cfl(::Type{<:Heat}) = (Δx, dHdT, kc) -> Utils.adstrip(Δx^2 * dHdT / kc)
+cfl!(::Type{T}) where {T<:SubSurfaceProcess} = (Δt, Δx, dHdT, kc) -> @. Δt = cfl(T)(Δx, dHdT, kc)
 
-# DiffEq/SciML driver
-export CryoGridProblem
-include("diffeq.jl")
+# DiffEq/SciML driver (possibly should be a soft dependency with Requires.jl)
+include("DiffEq/DiffEq.jl")
+@reexport using .DiffEq
 
 end

src/IO/InputOutput.jl

 module InputOutput
 
+import CryoGrid
 using CryoGrid.Numerics
-using CryoGrid.Land
 
+using Base: @propagate_inbounds
+using ComponentArrays
 using DataStructures: DefaultDict
 using Dates
 using DimensionalData
 using Downloads
+using Flatten
+using Interpolations
 using JSON3
-using Lazy: @>>
+using Lazy: @>>, groupby
+using ModelParameters
+using Tables
+using TimeSeries
 using Unitful
 
+import DimensionalData: stack
+
 export loadforcings
 
 include("ioutils.jl")
-
+export Forcing, TimeSeriesForcing
+include("forcings.jl")
+export CryoGridParams
+include("params.jl")
 export CryoGridOutput
 include("output.jl")
 

src/Physics/Boundaries/forcing.jl → src/IO/forcings.jl

@@ -4,8 +4,8 @@
 Abstract type representing a generic external boundary condition (i.e. "forcing").
 """
 abstract type Forcing{T,N} end
-(forcing::Forcing)(x::Number) = error("$(typeof(forcing)) not implemented")
-(forcing::Forcing)(t::DateTime) = forcing(ustrip(u"s", float(Dates.datetime2epochms(t))u"ms"))
+@inline @propagate_inbounds (forcing::Forcing)(x::Number) = error("$(typeof(forcing)) not implemented")
+@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
 
@@ -36,9 +36,9 @@ Base.show(io::IO, forcing::TimeSeriesForcing{T}) where T = print(io, "TimeSeries
 """
 Get interpolated forcing value at t seconds from t0.
 """
-(forcing::TimeSeriesForcing)(t::Number) = forcing.interp(t) # extract interpolation and evaluate
+@inline @propagate_inbounds (forcing::TimeSeriesForcing)(t::Number) = forcing.interp(t) # extract interpolation and evaluate
 
-Base.getindex(f::TimeSeriesForcing, i) = forcing.tarray[i]
+@inline @propagate_inbounds Base.getindex(forcing::TimeSeriesForcing, i) = forcing.tarray[i]
 function Base.getindex(f::TimeSeriesForcing{T,A,I}, range::StepRange{DateTime,TStep}) where {T,A,I,TStep}
       order(::Interpolations.GriddedInterpolation{T1,N,T2,Gridded{Torder}}) where {T1,N,T2,Torder} = Torder()
       subseries = f.tarray[range]

src/IO/output.jl

 """
-    CryoGridOutput{TRes}
+    CryoGridOutput{TSol}
 
 Helper type that stores the raw output from a CryoGrid run along with `DimArray` views of all
 logged variables. `CryoGridOutput` overrides `Base.getproperty` to allow for direct dot-syntax
 access of state variables. For example, if your model has a grid variable named `T`, `out.T` returns a `DimArray`
-with indexed time and depth axes. For OrdinaryDiffEq.jl outputs, the `ODESolution` can be accessed via `out.res`,
+with indexed time and depth axes. For OrdinaryDiffEq.jl outputs, the `ODESolution` can be accessed via `out.sol`,
 or for convenience, the continuous solution at time `t` can be computed via `out(t)` which is equivalent to
-`withaxes(out.res(t))`.
+`withaxes(out.sol(t))`.
 """
-struct CryoGridOutput{TRes}
+struct CryoGridOutput{TSol}
     ts::Vector{DateTime}
-    res::TRes
+    sol::TSol
     data::NamedTuple
-    CryoGridOutput(ts::Vector{DateTime}, res::TRes, data::NamedTuple) where TRes = new{TRes}(ts, res, data)
+    CryoGridOutput(ts::Vector{DateTime}, sol::TSol, data::NamedTuple) where TSol = new{TSol}(ts, sol, data)
 end
 # Overrides from Base
 function Base.show(io::IO, out::CryoGridOutput)
@@ -28,10 +28,14 @@ function Base.show(io::IO, out::CryoGridOutput)
         println(io, "   $r")
     end
 end
+stack(out::CryoGridOutput, var::Symbol, vars::Symbol...) = DimStack((;map(n -> n => getproperty(out, n), tuple(var, vars...))...))
+Base.keys(out::CryoGridOutput) = Base.propertynames(out.data)
+Base.propertynames(out::CryoGridOutput) = tuple(fieldnames(typeof(out))..., propertynames(out.data)...)
 function Base.getproperty(out::CryoGridOutput, sym::Symbol)
-    if sym in (:res,:ts,:data)
+    if sym in (:sol,:ts,:data)
         getfield(out, sym)
     else
         out.data[sym]
     end
 end
+Base.Dict(out::CryoGridOutput) = Dict(map(k -> string(k) => getproperty(out, k), keys(out))...)

src/IO/params.jl

+# Parameter array
+"""
+    CryoGridParams{T} <: DenseArray{T,1}
+
+Wraps a `ModelParameters.Model` parameter table for CryoGrid types. It is recommended *not* to use this
+type directly in math or linear algebra operations but rather to use `Base.values` to obtain a normal array
+of parameter values.
+"""
+struct CryoGridParams{T} <: DenseArray{T,1}
+    table::Model # param table
+    CryoGridParams(table::AbstractModel) where {T} = new{eltype(table[:val])}(table)
+end
+Base.values(ps::CryoGridParams) = ps.table[:val]
+Base.axes(ps::CryoGridParams) = axes(collect(values(ps)))
+Base.LinearIndices(ps::CryoGridParams) = LinearIndices(collect(values(ps)))
+Base.IndexStyle(::Type{<:CryoGridParams}) = Base.IndexLinear()
+Base.similar(ps::CryoGridParams) = CryoGridParams(Model(ps.table))
+Base.similar(ps::CryoGridParams, ::Type{T}) where T = CryoGridParams(Model(parent(ps.table)))
+Base.length(ps::CryoGridParams) = length(ps.table)
+Base.size(ps::CryoGridParams) = size(ps.table)
+Base.keys(ps::CryoGridParams) = keys(ps.table)
+Base.getindex(ps::CryoGridParams, i::Int) = values(ps)[i]
+Base.getindex(ps::CryoGridParams, col::Symbol) = ps.table[col]
+Base.setindex!(ps::CryoGridParams, val, i::Int) = ps[:val] = map(enumerate(values(ps))) do (j,p)
+    j == i ? val : p
+end
+function Base.setindex!(ps::CryoGridParams, vals, col::Symbol; kwargs...)
+    # TODO: replace this implementation when ModelParameters supports table row assignment
+    inds = findall(1:length(ps)) do i
+        all(ps[first(kw)][i] == last(kw) for kw in kwargs)
+    end
+    ps.table[col] = map(1:length(ps)) do i
+        r = searchsorted(inds, i)
+        length(r) == 1 ?  vals[first(r)] : ps[col][i]
+    end
+end
+function Base.show(io::IO, ::MIME"text/plain", ps::CryoGridParams{T}) where T
+    println(io, "CryoGridParams{$T} with $(length(ps)) parameters")
+    ModelParameters.printparams(io, ps.table)
+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
+            mpara = Model(para)
+            mpara[:paratype] = repeat([type], length(mpara))
+            mpara[:parafield] = repeat([name], length(mpara))
+            return parent(mpara)
+        else
+            return para
+        end
+    end
+    parameterizations = Flatten.flatten(parent(m), Flatten.flattenable, CryoGrid.Parameterization)
+    parameterization_fields = Flatten.fieldnameflatten(parent(m), Flatten.flattenable, CryoGrid.Parameterization)
+    parameterization_types = Flatten.metaflatten(parent(m), ModelParameters._fieldparentbasetype, CryoGrid.Parameterization)
+    updated_parameterizations = map(_setparafield, parameterization_fields, parameterization_types, parameterizations)
+    # reconstruct parent of `m` using updated parameterization parameters and then normalize parameters
+    # by rebuilding the parameter Model `m`.
+    newparent = Flatten.reconstruct(parent(m), updated_parameterizations, CryoGrid.Parameterization)
+    return Model(newparent)
+end

src/Land/params.jl

-abstract type ParamTransform end
-
-struct ParamMapping{T,name,layer}
-    transform::T
-    ParamMapping(transform::T, name::Symbol, layer::Symbol) where {T<:ParamTransform} = new{T,name,layer}(transform)
-end
-
-struct ParameterVector{T,TV,P,M} <: DenseArray{T,1}
-    vals::TV # input/reparameterized param vector
-    params::P # parameters grouped by layer and name
-    mappings::M # mapping metadata
-    ParameterVector(vals::TV,params::P,mappings::ParamMapping...) where {T,TV<:ComponentVector{T},P<:NamedTuple} = new{T,TV,P,typeof(mappings)}(vals,params,mappings)
-end
-vals(rv::ParameterVector) = getfield(rv, :vals)
-mappings(rv::ParameterVector) = getfield(rv, :mappings)
-Base.axes(rv::ParameterVector) = axes(getfield(rv, :vals))
-Base.LinearIndices(rv::ParameterVector) = LinearIndices(getfield(rv, :vals))
-Base.IndexStyle(::Type{<:ParameterVector}) = Base.IndexLinear()
-Base.similar(rv::ParameterVector) = ParameterVector(similar(vals(rv)), similar(pout(rv)), mappings(rv))
-Base.similar(rv::ParameterVector, ::Type{T}) where T = ParameterVector(similar(vals(rv), T), similar(pout(rv), T), mappings(rv))
-Base.length(rv::ParameterVector) = length(getfield(rv, :vals))
-Base.size(rv::ParameterVector) = size(getfield(rv, :vals))
-Base.getproperty(rv::ParameterVector, sym::Symbol) = getproperty(getfield(rv, :vals), sym)
-Base.getindex(rv::ParameterVector, i) = getfield(rv, :vals)[i]
-Base.setproperty!(rv::ParameterVector, val, i) = setproperty!(getfield(rv, :vals), val, sym)
-Base.setindex!(rv::ParameterVector, val, i) = setindex!(getfield(rv, :vals), val, i)
-Base.show(io, rv::ParameterVector) = show(io, getfield(rv, :vals))
-ComponentArrays.ComponentArray(rv::ParameterVector) = getfield(rv, :vals)
-
-_paramval(p::Param) = ustrip(p.val) # extracts value from Param type and strips units
-function parameters(model::LandModel, transforms::Pair{Symbol,<:Pair{Symbol,<:ParamTransform}}...)
-    function getparam(p)
-        # currently, we assume only one variable of each name in each layer;
-        # this could be relaxed in the future but will need to be appropriately handled
-        @assert length(p) == 1 "Found more than one parameter with name $var in $layer; this is not currently supported."
-        return p[1]
-    end
-    m = Model(model)
-    nestedparams = mapflat(getparam, groupparams(m, :layer, :fieldname); maptype=NamedTuple)
-    mappedparams = nestedparams
-    mappings = ParamMapping[]
-    for (layer,(var,transform)) in transforms
-        @set! mappedparams[layer][var] = mapflat(getparam, groupparams(Model(transform), :fieldname); maptype=NamedTuple)
-        push!(mappings, ParamMapping(transform, var, layer))
-    end
-    mappedarr = ComponentArray(mapflat(_paramval, mappedparams))
-    return ParameterVector(mappedarr, mappedparams, mappings...)
-end
-@inline @generated function updateparams!(v::AbstractVector, model::LandModel, u, du, t)
-    quote
-        p = ModelParameters.update(ModelParameters.params(model), v)
-        return Utils.genmap(_paramval, p)
-    end
-end
-@inline @generated function updateparams!(rv::ParameterVector{T,TV,P,M}, model::LandModel, u, du, t) where {T,TV,P,M}
-    expr = quote
-        pvals = vals(rv)
-        pmodel = ModelParameters.update(getfield(rv, :params), pvals)
-    end
-    # apply parameter transforms
-    for i in 1:length(M.parameters)
-        push!(expr.args, :(pmodel = _updateparam(pmodel, mappings(rv)[$i], model, u, du, t)))
-    end
-    # flatten parameters and strip Param types
-    push!(expr.args, :(return Utils.genmap(_paramval, ModelParameters.params(pmodel))))
-    return expr
-end
-@inline @generated function _updateparam(pmodel, mapping::ParamMapping{T,name,layer}, model::LandModel, u, du, t) where {T,name,layer}
-    quote
-        state = getstate(Val($(QuoteNode(layer))), model, u, du, t)
-        p = pmodel.$layer.$name
-        # reconstruct transform with new parameter values
-        op = ModelParameters.update(mapping.transform, ModelParameters.params(p))
-        # apply transform and replace parameter in named tuple
-        @set! pmodel.$layer.$name = Param(transform(state, op))
-        return pmodel
-    end
-end
-
-# Transform implementations
-@with_kw struct LinearTrend{P} <: ParamTransform
-    slope::P = Param(0.0)
-    intercept::P = Param(0.0)
-    tstart::Float64 = 0.0
-    tstop::Float64 = Inf; @assert tstop > tstart
-    minval::Float64 = -Inf
-    maxval::Float64 = Inf
-end
-function transform(state, trend::LinearTrend)
-    let t = min(state.t - trend.tstart, trend.tstop),
-        β = trend.slope,
-        α = trend.intercept;
-        min(max(β*t + α, trend.minval), trend.maxval)
-    end
-end

src/Layers/Layers.jl

-module Layers
-
-import CryoGrid: SubSurface
-import CryoGrid: initialcondition!, variables
-
-using CryoGrid.Numerics
-using CryoGrid.Utils
-
-using DimensionalData
-using IntervalSets
-using ModelParameters
-using Parameters
-using Unitful
-
-export Soil, SoilParameterization, SoilCharacteristicFractions, SoilProfile, SoilType, Sand, Silt, Clay, soilparameters
-include("soil.jl")
-
-end
\ No newline at end of file

src/Layers/soil.jl

-"""
-Represents the texture classification of the soil. Sand, Silt, and Clay are provided by default.
-"""
-abstract type SoilTexture end
-struct Sand <: SoilTexture end
-struct Silt <: SoilTexture end
-struct Clay <: SoilTexture end
-"""
-    SoilParameterization
-
-Abstract base type for parameterizations of soil properties.
-"""
-abstract type SoilParameterization end
-"""
-    SoilCharacteristicFractions{P1,P2,P3,P4} <: SoilParameterization
-
-Represents the composition of the soil in terms of fractions: excess ice, natural porosity, saturation, and organic/(mineral + organic).
-"""
-struct SoilCharacteristicFractions{P1,P2,P3,P4} <: SoilParameterization
-    xic::P1 # excess ice fraction
-    por::P2 # natural porosity
-    sat::P3 # saturation
-    org::P4 # organic fraction of solid; mineral fraction is 1-org
-    SoilCharacteristicFractions(xic::P1, por::P2, sat::P3, org::P4) where {P1,P2,P3,P4} = new{P1,P2,P3,P4}(xic,por,sat,org)
-end
-function soilparameters(;xic=0.0, por=0.5, sat=1.0, org=0.5)
-    params = Tuple(Param(p, bounds=(0.0,1.0)) for p in [xic,por,sat,org])
-    SoilCharacteristicFractions(params...)
-end
-SoilProfile(pairs::Pair{<:DistQuantity,<:SoilParameterization}...) = Profile(pairs...)
-# Helper functions for obtaining soil compositions from characteristic fractions.
-soilcomp(::Val{:θp}, χ, ϕ, θ, ω) = (1-χ)*ϕ
-soilcomp(::Val{:θw}, χ, ϕ, θ, ω) = χ + (1-χ)*ϕ*θ
-soilcomp(::Val{:θm}, χ, ϕ, θ, ω) = (1-χ)*(1-ϕ)*(1-ω)
-soilcomp(::Val{:θo}, χ, ϕ, θ, ω) = (1-χ)*(1-ϕ)*ω
-θx(fracs::SoilCharacteristicFractions) = fracs.xic
-θp(fracs::SoilCharacteristicFractions) = soilcomp(Val{:θp}(), fracs.xic, fracs.por, fracs.sat, fracs.org)
-θw(fracs::SoilCharacteristicFractions) = soilcomp(Val{:θw}(), fracs.xic, fracs.por, fracs.sat, fracs.org)
-θm(fracs::SoilCharacteristicFractions) = soilcomp(Val{:θm}(), fracs.xic, fracs.por, fracs.sat, fracs.org)
-θo(fracs::SoilCharacteristicFractions) = soilcomp(Val{:θo}(), fracs.xic, fracs.por, fracs.sat, fracs.org)
-"""
-Thermal conductivity constants.
-"""
-@with_kw struct SoilTCParams @deftype Float"W/(m*K)"
-    kw = 0.57xu"W/(m*K)" #water [Hillel(1982)]
-    ko = 0.25xu"W/(m*K)" #organic [Hillel(1982)]
-    km = 3.8xu"W/(m*K)" #mineral [Hillel(1982)]
-    ka = 0.025xu"W/(m*K)" #air [Hillel(1982)]
-    ki = 2.2xu"W/(m*K)" #ice [Hillel(1982)]
-end
-"""
-Heat capacity constants.
-"""
-@with_kw struct SoilHCParams @deftype Float"J/(K*m^3)"
-    cw = 4.2e6xu"J/(K*m^3)" #[J/m^3K] heat capacity water
-    co = 2.5e6xu"J/(K*m^3)" #[J/m^3K]  heat capacity organic
-    cm = 2e6xu"J/(K*m^3)" #[J/m^3K]  heat capacity mineral
-    ca = 0.00125e6xu"J/(K*m^3)" #[J/m^3K]  heat capacity pore space
-    ci = 1.9e6xu"J/(K*m^3)" #[J/m^3K]  heat capacity ice
-end
-"""
-Basic Soil layer.
-"""
-@with_kw struct Soil{T,P<:SoilParameterization,S} <: SubSurface
-    texture::T = Sand()
-    para::P = soilparameters()
-    tc::SoilTCParams = SoilTCParams()
-    hc::SoilHCParams = SoilHCParams()
-    sp::S = nothing # user-defined specialization
-end

src/Numerics/Numerics.jl

 module Numerics
 
 import Base.==
-import ExprTools
+import ConstructionBase
 import ForwardDiff
+import ModelParameters
 import PreallocationTools as Prealloc
 
 using CryoGrid.Utils
 
 using Base: @inbounds, @propagate_inbounds
-using ConstructionBase
 using ComponentArrays
 using DimensionalData: AbstractDimArray, DimArray, Dim, At, dims, Z
 using Flatten
 using IfElse
-using Interpolations: Interpolations, Gridded, Linear, Flat, Line, interpolate, extrapolate
+using Interpolations
 using IntervalSets
 using LinearAlgebra
 using LoopVectorization
 using RuntimeGeneratedFunctions
 using Unitful
 using StructTypes
-using Symbolics
-using SymbolicUtils
 
 RuntimeGeneratedFunctions.init(@__MODULE__)
 
 export GridSpec, Edges, Cells
+export Profile, ProfileKnot
+
+USE_TURBO = true
+function turbo(value::Bool)
+    global USE_TURBO = value
+end
 
 abstract type AbstractDiscretization{Q,N} <: DenseArray{Q,N} end
 
@@ -36,17 +40,39 @@ struct Cells <: GridSpec end
 abstract type Geometry end
 struct UnitVolume <: Geometry end
 
-export ∇
+struct ProfileKnot{D,T}
+    depth::D
+    value::T
+end
+Base.show(io::IO, knot::ProfileKnot) = print(io, "$(knot.depth): $(knot.value)")
+struct Profile{N,TKnots}
+    knots::TKnots
+    Profile(::Tuple{}) = new{0,Tuple{}}(())
+    Profile(knots::Tuple{Vararg{<:ProfileKnot,N}}) where N = new{N,typeof(knots)}(knots)
+    Profile(pairs::Tuple{Vararg{<:Pair}}) where D = Profile(map(Base.splat(ProfileKnot), pairs))
+    Profile(pairs::Pair...) = Profile(pairs)
+end
+Flatten.flattenable(::Type{<:ProfileKnot}, ::Type{Val{:depth}}) = false
+Flatten.flattenable(::Type{ProfileKnot{D,T}}, ::Type{Val{:value}}) where {D,T<:Number} = false
+Flatten.flattenable(::Type{ProfileKnot{D,T}}, ::Type{Val{:depth}}) where {D<:ModelParameters.Param,T} = true
+Flatten.flattenable(::Type{ProfileKnot{D,T}}, ::Type{Val{:value}}) where {D,T<:ModelParameters.Param} = true
+Base.length(::Profile{N}) where N = N
+Base.iterate(profile::Profile) = iterate(profile.knots)
+Base.iterate(profile::Profile, state) = iterate(profile.knots, state)
+Base.getindex(profile::Profile, itrv::Interval) = Profile(Tuple(knot for knot in profile.knots if knot.depth ∈ itrv))
+Base.getindex(profile::Profile, i::Int) = profile.knots[i]
+Base.getindex(profile::Profile, i) = Profile(profile.knots[i])
+Base.lastindex(profile::Profile) = lastindex(profile.knots)
+StructTypes.StructType(::Type{<:Profile}) = StructTypes.UnorderedStruct()
+
+export Tabulated, ∇
 include("math.jl")
 
-export Grid, cells, edges, indexmap, subgridinds, Δ, volume, area
+export Grid, cells, edges, subgridinds, Δ, volume, area, updategrid!
 include("grid.jl")
 
-export Profile, profile2array, array2profile
-include("profile.jl")
-
 export Var, Prognostic, Algebraic, Diagnostic, VarDim, OnGrid, Shape, Scalar
-export varname, vartype, vardims, isprognostic, isalgebraic, isflux, isdiagnostic, isongrid, dimlength
+export varname, vartype, vardims, varunits, isprognostic, isalgebraic, isflux, isdiagnostic, isongrid, dimlength
 include("variables.jl")
 
 export VarStates, DiffCache, retrieve, getvar, getvars
@@ -54,7 +80,7 @@ include("varstates.jl")
 
 include("discretize.jl")
 
-export initializer, init!
+export ConstantInitializer, InterpInitializer, initializer, init!
 include("init.jl")
 
 end

src/Numerics/grid.jl

-const GridValues{Q,A} = NamedTuple{(:edges,:cells),NTuple{2,SubArray{Q,1,A,Tuple{StepRange{Int64,Int64}},true}}} where {Q,A<:AbstractArray}
+const GridValues{A} = NamedTuple{(:edges,:cells),NTuple{2,A}} where {A<:AbstractVector}
 
 """
     struct Grid{S,G,Q,A} <: AbstractDiscretization{Q,1}
@@ -11,72 +11,87 @@ array type, and `Q` is the numerical type (e.g. `Float64` or a `Unitful.Quantity
 """
 struct Grid{S,G,Q,A} <: AbstractDiscretization{Q,1}
     geometry::G
-    values::GridValues{Q,A}
-    deltas::GridValues{Q,A}
-    Grid(::Type{S}, values::GridValues{Q,A}, deltas::GridValues{Q,A}, geom::G) where {S<:GridSpec,Q,A,G<:Geometry} = new{S,G,Q,A}(geom,values,deltas)
+    values::GridValues{A}
+    deltas::GridValues{A}
+    bounds::UnitRange{Int}
+    Grid(::Type{S}, values::GridValues{A}, deltas::GridValues{A}, geom::G, bounds::UnitRange{Int}=1:length(values)) where {S<:GridSpec,Q,A<:AbstractVector{Q},G<:Geometry} = new{S,G,Q,A}(geom,values,deltas,bounds)
     function Grid(vals::AbstractVector{Q}, geometry::G=UnitVolume()) where {G<:Geometry,Q<:Number}
         @assert issorted(vals) "grid values should be in ascending order"
         nedges = length(vals)
-        ncells = nedges - 1 
+        ncells = nedges - 1
         edges = copyto!(similar(vals), vals)
         cells = similar(edges, ncells)
         @. cells = (edges[1:end-1] + edges[2:end]) / (2*one(Q))
-        # note that the distinction between edges and cells here is kind of redundant since the distance between
-        # midpoints should always be the same as distance between boundaries...
         Δedges = edges[2:end] .- edges[1:end-1]
         Δcells = cells[2:end] .- cells[1:end-1]
-        new{Edges,G,Q,typeof(edges)}(geometry,
-            (edges=view(edges,1:1:nedges),cells=view(cells,1:1:ncells)),
-            (edges=view(Δedges,1:1:nedges-1),cells=view(Δcells,1:1:ncells-1)),
-        )
+        new{Edges,G,Q,typeof(edges)}(geometry, (edges=edges, cells=cells), (edges=Δedges, cells=Δcells), 1:length(edges))
     end
-    function Grid(grid::Grid{Edges,G,Q,A}, interval::ClosedInterval{Int}) where {G,Q,A}
-        start, stop = interval.left, interval.right
-        edges = @view grid.values.edges[start:stop]
-        cells = @view grid.values.cells[start:stop-1]
-        Δedges = @view grid.deltas.edges[start:stop-1]
-        Δcells = @view grid.deltas.cells[start:stop-2]
-        new{Edges,G,Q,A}(grid.geometry, (edges=edges,cells=cells), (edges=Δedges,cells=Δcells))
+    function Grid(grid::Grid{S,G,Q,A}, interval::ClosedInterval{Int}) where {S,G,Q,A}
+        new{S,G,Q,A}(grid.geometry, grid.values, grid.deltas, leftendpoint(interval):rightendpoint(interval))
     end
-    function Grid(grid::Grid{Cells,G,Q,A}, interval::ClosedInterval{Int}) where {G,Q,A}
-        start, stop = interval.left, interval.right
-        edges = @view grid.values.edges[start:stop+1]
-        cells = @view grid.values.cells[start:stop]
-        Δedges = @view grid.deltas.edges[start:stop]
-        Δcells = @view grid.deltas.cells[start:stop-1]
-        new{Cells,G,Q,A}(grid.geometry, (edges=edges,cells=cells), (edges=Δedges,cells=Δcells))
+    Grid(::Type{Cells}, grid::Grid{Edges,G,Q,A}) where {G,Q,A} = new{Cells,G,Q,A}(grid.geometry, grid.values, grid.deltas, grid.bounds)
+    Grid(::Type{Edges}, grid::Grid{Cells,G,Q,A}) where {G,Q,A} = new{Edges,G,Q,A}(grid.geometry, grid.values, grid.deltas, grid.bounds)
+end
+ConstructionBase.constructorof(::Type{Grid{S,G,Q,A}}) where {S,G,Q,A} = (geom,values,deltas,bounds) -> Grid(S,values,deltas,geom,bounds)
+Base.show(io::IO, ::MIME"text/plain", grid::Grid) = show(io, grid)
+function Base.show(io::IO, grid::Grid{S,G}) where {S,G}
+    if length(grid) == length(grid.values.edges)
+        print(io, "Grid{$S}($(grid[1])..$(grid[end])) of length $(length(grid)) with geometry $G")
+    else
+        print(io, "Grid{$S}($(grid[1])..$(grid[end])) of length $(length(grid)) (child of Grid{$S}$(parent(grid)[1])..$(parent(grid)[end]) of length $(length(parent(grid)))) with geometry $G")
     end
-    Grid(::Type{Cells}, grid::Grid{Edges,G,Q,A}) where {G,Q,A} = new{Cells,G,Q,A}(grid.geometry,grid.values,grid.deltas)
-    Grid(::Type{Edges}, grid::Grid{Cells,G,Q,A}) where {G,Q,A} = new{Edges,G,Q,A}(grid.geometry,grid.values,grid.deltas)
 end
-ConstructionBase.constructorof(::Type{Grid{S,G,Q,A}}) where {S,G,Q,A} = (geom,values,deltas) -> Grid(S,values,deltas,geom)
-Base.show(io::IO, grid::Grid{S,G}) where {S,G} = print(io, "Grid{$S}($(grid[1])..$(grid[end])) of length $(length(grid)) with geometry $G")
-Base.show(io::IO, ::MIME{Symbol("text/plain")}, grid::Grid) = show(io, grid)
 
-function subgridinds(grid::Grid{S,G,Q,A}, interval::Interval{L,R,Q}) where {S,G,Q,A,L,R}
+function subgridinds(grid::Grid, interval::Interval{L,R}) where {L,R}
     @assert interval.left <= interval.right "Invalid interval: $interval"
-    vals = values(grid)
     # Determine indices which lie in the given interval
-    l_ind = searchsortedfirst(vals, interval.left)
-    r_ind = searchsortedlast(vals, interval.right)
+    l_ind = searchsortedfirst(grid, interval.left)
+    r_ind = searchsortedlast(grid, interval.right)
     return (L == :closed ? l_ind : l_ind + 1)..(R == :closed ? r_ind : r_ind - 1)
 end
-@inline Δ(grid::Grid{Edges}) = grid.deltas.edges
-@inline Δ(grid::Grid{Cells}) = grid.deltas.cells
+@inline bounds(grid::Grid{Edges}) = grid.bounds
+@inline bounds(grid::Grid{Cells}) = first(grid.bounds):last(grid.bounds)-1
+@inline Δbounds(grid::Grid{Edges}) = first(grid.bounds):last(grid.bounds)-1
+@inline Δbounds(grid::Grid{Cells}) = first(grid.bounds):last(grid.bounds)-2
+@inline Δ(grid::Grid{Edges}) = view(grid.deltas.edges, Δbounds(grid))
+@inline Δ(grid::Grid{Cells}) = view(grid.deltas.cells, Δbounds(grid))
 @inline cells(grid::Grid{Edges}) = Grid(Cells, grid)
 @inline cells(grid::Grid{Cells}) = grid
 @inline edges(grid::Grid{Edges}) = grid
 @inline edges(grid::Grid{Cells}) = Grid(Edges, grid)
-@inline Base.values(grid::Grid{Edges}) = grid.values.edges
-@inline Base.values(grid::Grid{Cells}) = grid.values.cells
+@inline Base.parent(grid::Grid{Edges}) = Grid(grid, 1..length(grid.values.edges))
+@inline Base.parent(grid::Grid{Cells}) = Grid(grid, 1..length(grid.values.cells))
+@inline Base.collect(grid::Grid) = collect(values(grid))
+@inline Base.values(grid::Grid{Edges}) = view(grid.values.edges, bounds(grid))
+@inline Base.values(grid::Grid{Cells}) = view(grid.values.cells, bounds(grid))
 @inline Base.similar(grid::Grid{Edges}) = Grid(copy(grid.values.edges))
 @inline Base.similar(grid::Grid{Cells}) = similar(edges(grid)) |> cells
-@inline Base.size(grid::Grid) = size(values(grid))
-@inline Base.length(grid::Grid) = length(values(grid))
+@inline Base.size(grid::Grid) = (length(grid),)
+@inline Base.length(grid::Grid) = last(bounds(grid)) - first(bounds(grid)) + 1
+@inline Base.firstindex(grid::Grid) = 1
+@inline Base.lastindex(grid::Grid) = length(grid)
 @propagate_inbounds Base.getindex(grid::Grid, i::Int) = values(grid)[i]
-@propagate_inbounds Base.getindex(grid::Grid{S,G,Q,A}, interval::Interval{L,R,Q}) where {S,G,Q,A,L,R} = Grid(grid, subgridinds(grid,interval))
-@propagate_inbounds Base.getindex(grid::Grid, interval::Interval{L,R,Int}) where {L,R} = Grid(grid, interval)
-Base.setindex!(grid::Grid, args...) = error("setindex! is not allowed for Grid types")
+@propagate_inbounds Base.getindex(grid::Grid, i::AbstractRange) = grid[grid[first(i)]..grid[last(i)]]
+@propagate_inbounds Base.getindex(grid::Grid{S,G,Q,A}, interval::Interval{L,R,Q}) where {S,G,Q,A,L,R} = grid[subgridinds(grid, interval)]
+@propagate_inbounds Base.getindex(grid::Grid, interval::Interval{L,R,Int}) where {L,R} = Grid(grid, first(grid.bounds)+interval.left-1..first(grid.bounds)+interval.right-1)
+Base.setindex!(::Grid, args...) = error("setindex! is not allowed for Grid types")
+"""
+    updategrid!(grid::Grid{Edges,G,Q}, vals::Q) where {G,Q}
+
+Overwrites `grid` edges with `vals`, and recomputes grid centers/deltas to be consistent with the new grid.
+"""
+function updategrid!(grid::Grid{Edges,G,Q}, vals::AbstractVector{Q}) where {G,Q}
+    z_edges = values(grid)
+    z_cells = values(cells(grid))
+    Δz_edges = Δ(grid)
+    Δz_cells = Δ(cells(grid))
+    z_edges .= vals
+    z_cells .= (z_edges[1:end-1] .+ z_edges[2:end]) ./ (2*one(Q))
+    Δz_edges .= z_edges[2:end] .- z_edges[1:end-1]
+    Δz_cells .= z_cells[2:end] .- z_cells[1:end-1]
+    @assert issorted(parent(grid)) "updated grid values are invalid; grid edges must be strictly non-decreasing"
+    return grid
+end
 
 # unit volume
 @inline volume(grid::Grid{Cells,UnitVolume,Q}) where Q = Δ(edges(grid)).*oneunit(Q)^2

src/Numerics/init.jl

-abstract type VarInit{varname} end
-Numerics.varname(::VarInit{varname}) where {varname} = varname
+abstract type VarInitializer{varname} end
+Numerics.varname(::VarInitializer{varname}) where {varname} = varname
+ConstructionBase.constructorof(::Type{T}) where {varname,T<:VarInitializer{varname}} = (args...) -> T.name.wrapper(varname, args...)
 """
-    init!(x::AbstractVector, init::VarInit, args...)
+    FunctionInitializer{varname,F} <: VarInitializer{varname}
 
-Initializes state variable `x` using initializer `init` and implementation-specific additional
-arguments.
+Initializes a scalar or vector-valued state variable using an arbitrary method/function.
 """
-init!(x::AbstractVector, init::VarInit, args...) = error("not implemented")
-"""
-    ConstantInit{varname,T} <: VarInit
-
-Initializes a scalar or vector-valued state variable with a pre-specified constant value.
-"""
-struct ConstantInit{varname,T} <: VarInit{varname}
-    value::T
-    ConstantInit(varname::Symbol, value::T) where {T} = new{varname,T}(value)
+struct FunctionInitializer{varname,F} <: VarInitializer{varname}
+    f::F
+    FunctionInitializer(varname::Symbol, f::F) where {F} = new{varname,F}(f)
 end
-init!(x::AbstractVector, init::ConstantInit) = @. x = init.value
+(init::FunctionInitializer)(args...) = f(args...)
+Base.getindex(init::FunctionInitializer, itrv::Interval) = init
 """
-    InterpInit{varname,P,I,E} <: VarInit
+    InterpInitializer{varname,P,I,E} <: VarInitializer{varname}
 
-Init for on-grid variables that takes a `Profile` as initial values
-and interpolates along the model grid. The interpolation mode is linear by default,
-but can also be quadratic, cubic, or any other `Gridded` interpolation mode supported
+Initializer for on-grid variables that takes a `Profile` as initial values and interpolates along the model grid.
+The interpolation mode is linear by default, but can also be any other `Gridded` interpolation mode supported
 by `Interpolations.jl`.
 """
-struct InterpInit{varname,P,I,E} <: VarInit{varname}
+struct InterpInitializer{varname,P,I,E} <: VarInitializer{varname}
     profile::P
     interp::I
     extrap::E
-    InterpInit(varname::Symbol, profile::P, interp::I=Linear(), extrap::E=Flat()) where {P<:Profile,I,E} = new{varname,P,I,E}(profile, interp)
+    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
+    profile, interp, extrap = init.profile, init.interp, init.extrap
+    depths = collect(map(knot -> dustrip(knot.depth), profile.knots))
+    u = getproperty(state, var)
+    z = getproperty(state.grids, var)
+    if length(depths) > 1
+        f = extrapolate(interpolate((depths,), collect(map(knot -> dustrip(knot.value), profile.knots)), Gridded(interp)), extrap)
+        @. u = f(z)
+        return u
+    else
+        # if only one knot is defined, set to this value over all z;
+        # this is necessary because interpolate does not support interpolation grids of length 1
+        y = dustrip(profile.knots[1].value)
+        @. u = y
+        return u
+    end
 end
-function init!(x::AbstractVector, init::InterpInit, grid::AbstractVector)
-    z = collect(map(dustrip, init.profile.depths))
-    f = extrapolate(interpolate((z,), collect(map(dustrip, init.profile.values)), Gridded(init.interp)), init.extrap)
-    @. x = f(grid)
+# automatic partitioning of profile based on interval
+Base.getindex(init::InterpInitializer{var}, itrv::Interval) where var = InterpInitializer(var, init.profile[itrv], init.interp, init.extrap)
+# constructor for constant initializer function
+function _constantinit(::Val{varname}, x) where {varname}
+    return state -> state[varname] .= x
 end
 """
-    initializer(varname::Symbol, value::T) where {T} => ConstantInit
-    initializer(varname::Symbol, profile::Profile, interp=Linear(), extrap=Flat()) => InterpInit
+    initializer(varname::Symbol, x::Number) => FunctionInitializer w/ constant
+    initializer(varname::Symbol, f::Function) => FunctionInitializer
+    initializer(varname::Symbol, profile::Profile, interp=Linear(), extrap=Flat()) => InterpInitializer
 
-Convenience constructor for `VarInit`s that selects the appropriate initializer type based on the arguments.
+Convenience constructor for `VarInitializer` that selects the appropriate initializer type based on the arguments.
 """
-initializer(varname::Symbol, value::T) where {T} = ConstantInit(varname, value)
-initializer(varname::Symbol, profile::Profile, interp=Linear(), extrap=Flat()) = InterpInit(varname, profile, interp, extrap)
+initializer(varname::Symbol, x::Number) = FunctionInitializer(varname, _constantinit(Val{varname}(), x))
+initializer(varname::Symbol, f::Function) = FunctionInitializer(varname, f)
+initializer(varname::Symbol, profile::Profile, interp=Linear(), extrap=Flat()) = InterpInitializer(varname, profile, interp, extrap)

src/Numerics/math.jl

+∇(f, x) = ∇(typeof(x), f, x)
+function ∇(::Type{T}, f, x) where {T}
+    res = ForwardDiff.derivative!(ForwardDiff.DiffResult(zero(T), x), f, x)
+    return res.value, res.derivs[1]
+end
+
 """
     finitediff!(∂x::AbstractVector, x::AbstractVector, Δ::AbstractVector)
 
@@ -9,7 +15,6 @@ function finitediff!(∂x::AbstractVector, x::AbstractVector, Δ::AbstractVector
         @. ∂x = (x₂ - x₁) / Δ
     end
 end
-
 """
     lineardiffusion!(∂y::AbstractVector, x::AbstractVector, Δ::AbstractVector, k::Number)
 
@@ -24,12 +29,32 @@ function lineardiffusion!(∂y::AbstractVector, x::AbstractVector, Δ::AbstractV
         @. ∂y = k*((x₃ - x₂)/Δ₂ - (x₂-x₁)/Δ₁)/Δ₁
     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)
+    end
+end
+"""
+    nonlineardiffusion(x₁, x₂, x₃, k₁, k₂, Δx₁, Δx₂, Δk)
 
-@propagate_inbounds function _nonlineardiffusion!(∂y, x₁, x₂, x₃, k₁, k₂, Δx₁, Δx₂, Δk)
-    @. ∂y = (k₂*(x₃ - x₂)/Δx₂ - k₁*(x₂-x₁)/Δx₁)/Δk
+Scalar second order Laplacian with non-linear diffusion operator, `k`.
+"""
+@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!(
+@propagate_inbounds function nonlineardiffusion!(
     ∂y::AbstractVector{Float64},
     x₁::AbstractVector{Float64},
     x₂::AbstractVector{Float64},
@@ -40,29 +65,10 @@ end
     Δx₂::AbstractVector{Float64},
     Δk::AbstractVector{Float64},
 )
-    @turbo @. ∂y = (k₂*(x₃ - x₂)/Δx₂ - k₁*(x₂-x₁)/Δx₁)/Δk
-end
-
-"""
-    nonlineardiffusion!(
-        ∂y::AbstractVector{Float64},
-        x::AbstractVector{Float64}, 
-        Δx::AbstractVector{Float64},
-        k::AbstractVector{Float64},
-        Δk::AbstractVector{Float64}
-    )
-
-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)
+    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)
     end
 end
 
@@ -87,6 +93,38 @@ function harmonicmean!(h::AbstractVector, x::AbstractVector, w::AbstractVector)
     end
 end
 
+"""
+    tdma_solve!(x, a, b, c, d)
+
+Tridiagonal matrix solver; borrowed from CryoGridLite. Modifies all input vectors in-place.
+"""
+function tdma_solve!(x, a, b, c, d, factorized=false)
+    #a, b, c are the column vectors for the compressed tridiagonal matrix, d is the right vector
+    n = length(b); # n is the number of rows
+    #x = zeros(n,1);
+    if !factorized
+        # Modify the first-row coefficients
+        c[1] = c[1] / b[1];    # Division by zero risk.
+        d[1] = d[1] / b[1];    # Division by zero would imply a singular matrix.
+
+        @inbounds @fastmath for i = 2:n-1
+            temp = b[i] - a[i] * c[i-1];
+            c[i] = c[i] / temp;
+            d[i] = (d[i] - a[i] * d[i-1]) / temp;
+        end
+
+        d[n] = (d[n] - a[n] * d[n-1])/( b[n] - a[n] * c[n-1]);
+    end
+
+    # Now back substitute.
+    x[n] = d[n];
+    @inbounds @fastmath for i = n-1:-1:1
+        x[i] = d[i] - c[i] * x[i + 1];
+    end
+
+    return x
+end
+
 """
     heaviside(x)
 
@@ -94,7 +132,7 @@ Differentiable implementation of heaviside step function, i.e:
 
 ``h(x) = \\begin{cases} 1 & x ≥ 0 \\\\ 0 & x < 0 \\end{cases}``
 """
-heaviside(x) = IfElse.ifelse(x >= 0.0, 1.0, 0.0)
+heaviside(x) = IfElse.ifelse(x >= zero(x), one(x), zero(x))
 """
     logistic(x)
 
@@ -102,7 +140,7 @@ Numerically stable logistic function.
 
 ``σ(x) = \\begin{cases} \\frac{1}{1+\\exp(-x)} & x ≥ 0 \\\\ \\frac{\\exp(x)}{1+\\exp(x)} & x < 0 \\end{cases}``
 """
-logistic(x) = IfElse.ifelse(x >= 0, 1 / (1 + exp(-x)), exp(x) / (1 + exp(x)))
+logistic(x) = IfElse.ifelse(x >= zero(x), 1 / (1 + exp(-x)), exp(x) / (1 + exp(x)))
 """
     logit(x)
 
@@ -132,45 +170,37 @@ softplusinv(x) = let x = clamp(x, eps(), Inf); IfElse.ifelse(x > 34, x, log(exp(
 minusone(x) = x .- one.(x)
 plusone(x) = x .+ one.(x)
 
-"""
-    ∇(f, dvar::Symbol)
-
-Automatically generates an analytical partial derivative of `f` w.r.t `dvar` using Symbolics.jl.
-To avoid symbolic tracing issues, the function should 1) be pure (no side effects or non-mathematical behavior) and 2) avoid
-indeterminate control flow such as if-else or while blocks (technically should work but sometimes doesn't...). Conditional
-logic can be included using `IfElse.ifelse`. Additional argument names are extracted automatically from the method signature
-of `f`. Keyword arg `choosefn` should be a function which selects from available methods of `f` (returned by `methods`); defaults to `first`.
-Note that `∇` uses `RuntimeGeneratedFunction` to produce a fully specialized and compiled Julia function; it may be slow on the first call
-(due to compilation), but should be just as fast as handwriting it on subsequent calls.
-
-Example:
-
-For ``f(x,y) = 2x + xy``, ``\\frac{\\partial f}{\\partial x} = 2 + y``. Using `∇`, we can obtain this automagically:
-
-```jldoctest
-f(x,y) = 2*x + x*y
-∇f_x = ∇(f,:x)
-∇f_x(2.0,3.0)
-
-# output
-
-5.0
-```
-"""
-function ∇(f, dvar::Symbol; choosefn=first, context_module=Numerics)
-    # Parse function parameter names using ExprTools
-    fms = ExprTools.methods(f)
-    symbol(arg::Symbol) = arg
-    symbol(expr::Expr) = expr.args[1]
-    argnames = map(symbol, ExprTools.signature(choosefn(fms))[:args])
-    @assert dvar in argnames "function must have $dvar as an argument"
-    dind = findfirst(s -> s == dvar, argnames)
-    # Convert to symbols
-    argsyms = map(s -> Symbolics.Num(SymbolicUtils.Sym{Real}(s)), argnames)
-    # Generate analytical derivative of f
-    x = argsyms[dind]
-    ∂x = Differential(x)
-    ∇f_expr = build_function(∂x(f(argsyms...)) |> expand_derivatives,argsyms...)
-    ∇f = @RuntimeGeneratedFunction(context_module, ∇f_expr)
-    return ∇f
+# 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

src/Numerics/profile.jl

-struct Profile{N,D<:DistQuantity,T}
-    depths::NTuple{N,D}
-    values::NTuple{N,T}
-    Profile(depths::NTuple{N,D}, values::NTuple{N,T}) where {N,D,T} = new{N,D,T}(depths,values)
-    Profile(pairs::NTuple{N,Pair{D,T}}) where {N,D,T} = new{N,D,T}(map(first,pairs), map(last,pairs))
-    Profile(pairs::Pair{D,T}...) where {D,T} = Profile(pairs)
-end
-Flatten.flattenable(::Type{<:Profile}, ::Type{Val{:depths}}) = false
-Base.length(::Profile{N}) where N = N
-Base.iterate(profile::Profile) = iterate(zip(profile.depths,profile.values))
-Base.iterate(profile::Profile, state) = iterate(zip(profile.depths,profile.values),state)
-StructTypes.StructType(::Type{<:Profile}) = StructTypes.UnorderedStruct()
-"""
-    profile2array(profile::Profile{N,D,T};names) where {N,D,T}
-
-Constructs a DimArray from the given Profile, i.e. pairs Q => (x1,...,xn) where x1...xn are the values defined at Q.
-Column names for the resulting DimArray can be set via the names parameter which accepts an NTuple of symbols,
-where N must match the number of parameters given (i.e. n).
-"""
-function profile2array(profile::Profile{N,D,T};names::Union{Nothing,NTuple{M,Symbol}}=nothing) where {M,N,D,T}
-    depths, vals = zip(profile.values...)
-    params = hcat(collect.(vals)...)'
-    names = isnothing(names) ? [Symbol(:x,:($i)) for i in 1:N] : collect(names)
-    DimArray(params, (Z(collect(depths)), Dim{:var}(names)))
-end
-array2profile(arr::AbstractDimArray) = Profile(collect(dims(profile,Z)), mapslices(collect, arr; dims=2))

src/Numerics/variables.jl

@@ -12,41 +12,65 @@ dimlength(::Shape{dims}, grid::Grid) where dims = prod(dims)
 dimlength(d::OnGrid{Cells}, grid::Grid) = d.f(length(cells(grid)))
 dimlength(d::OnGrid{Edges}, grid::Grid) = d.f(length(edges(grid)))
 
-abstract type Var{name,T,S<:VarDim} end
-struct Prognostic{name,T,S} <: Var{name,T,S}
+abstract type Var{name,S<:VarDim,T,units} end
+"""
+    Prognostic{name,S,T,units} <: Var{name,S,T,units}
+
+Defines a prognostic (time-integrated) state variable.
+"""
+struct Prognostic{name,S,T,units} <: Var{name,S,T,units}
     dim::S
-    Prognostic(name::Symbol, ::Type{T}, dims::Union{<:Shape,OnGrid{Cells,typeof(identity)}}) where {T} = new{name, T, typeof(dims)}(dims)
-    Prognostic(::Symbol, ::Type{T}, dims::OnGrid) where {T} = error("Off-cell prognostic/algebraic spatial variables are not currently supported.")
+    Prognostic(name::Symbol, dims::Union{<:Shape,OnGrid{Cells,typeof(identity)}}, units=NoUnits, ::Type{T}=Float64) where {T} = new{name,typeof(dims),T,units}(dims)
+    Prognostic(::Symbol, dims::OnGrid, args...) where {T} = error("Off-cell prognostic/algebraic spatial variables are not currently supported.")
 end
-struct Algebraic{name,T,S} <: Var{name,T,S}
+"""
+    Algebraic{name,S,T,units} <: Var{name,S,T,units}
+
+Defines an algebraic (implicit) state variable.
+"""
+struct Algebraic{name,S,T,units} <: Var{name,S,T,units}
     dim::S
     # maybe a mass matrix init function?
-    Algebraic(name::Symbol, ::Type{T}, dims::Union{<:Shape,OnGrid{Cells,typeof(identity)}}) where {T} = new{name, T, typeof(dims)}(dims)
-    Algebraic(::Symbol, ::Type{T}, dims::OnGrid) where {T} = error("Off-cell prognostic/algebraic spatial variables are not currently supported.")
+    Algebraic(name::Symbol, dims::Union{<:Shape,OnGrid{Cells,typeof(identity)}}, units=NoUnits, ::Type{T}=Float64) where {T} = new{name,typeof(dims),T,units}(dims)
+    Algebraic(::Symbol, dims::OnGrid, args...) where {T} = error("Off-cell prognostic/algebraic spatial variables are not currently supported.")
 end
-struct Flux{dname,name,T,S} <: Var{dname,T,S}
+"""
+    Delta{dname,name,S,T,units} <: Var{dname,S,T,units}
+
+Defines a "delta" term `du` for variable `u`, which is the time-derivative or flux for prognostic variables and
+the residual for algebraic variables.
+"""
+struct Delta{dname,name,S,T,units} <: Var{dname,S,T,units}
     dim::S
-    Flux(::Symbol, ::Type{T}, dims::OnGrid) where {T} = error("Off-cell prognostic/algebraic spatial variables are not currently supported.")
-    Flux(dname::Symbol, name::Symbol, ::Type{T}, dims::Union{<:Shape,OnGrid{Cells,typeof(identity)}}) where {T} = new{dname, name, T, typeof(dims)}(dims)
-    Flux(var::Prognostic{name,T,S}) where {name,T,S} = let dims=vardims(var); new{Symbol(:d,name), name, T, typeof(dims)}(dims) end
-    Flux(var::Algebraic{name,T,S}) where {name,T,S} = let dims=vardims(var); new{Symbol(:d,name), name, T, typeof(dims)}(dims) end
+    Delta(::Symbol, dims::OnGrid, args...) where {T} = 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) where {T} = new{dname,name,typeof(dims),T,units}(dims)
+    Delta(var::Prognostic{name,S,T,units}) where {name,S,T,units} = let dims=vardims(var); new{Symbol(:d,name),name,typeof(dims),T,upreferred(units)/u"s"}(dims) end
+    Delta(var::Algebraic{name,S,T,units}) where {name,S,T,units} = let dims=vardims(var); new{Symbol(:d,name),name,typeof(dims),T,units}(dims) end
 end
-struct Diagnostic{name,T,S} <: Var{name,T,S}
+"""
+    Diagnostic{name,S,T,units} <: Var{name,S,T,units}
+
+Defines a diagnostic variable which is allocated and cached per timestep but not integrated over time.
+"""
+struct Diagnostic{name,S,T,units} <: Var{name,S,T,units}
     dim::S
-    Diagnostic(name::Symbol, ::Type{T}, dims::VarDim) where {T} = new{name, T, typeof(dims)}(dims)
+    Diagnostic(name::Symbol, dims::VarDim, units=NoUnits, ::Type{T}=Float64) where {T} = new{name,typeof(dims),T,units}(dims)
 end
-ConstructionBase.constructorof(::Type{Prognostic{name,T,S}}) where {name,T,S} = s -> Prognostic(name, T, s)
-ConstructionBase.constructorof(::Type{Diagnostic{name,T,S}}) where {name,T,S} = s -> Diagnostic(name, T, s)
-ConstructionBase.constructorof(::Type{Algebraic{name,T,S}}) where {name,T,S} = s -> Algebraic(name, T, s)
-ConstructionBase.constructorof(::Type{Flux{dname,name,T,S}}) where {dname,name,T,S} = s -> Flux(dname, name, T, s)
-==(var1::Var{N1,T1,D1},var2::Var{N2,T2,D2}) where {N1,N2,T1,T2,D1,D2} = (N1==N2) && (T1==T2) && (D1==D2)
+# constructors for Flatten.reconstruct
+ConstructionBase.constructorof(::Type{Prognostic{name,S,T,units}}) where {name,S,T,units} = s -> Prognostic(name, s, T, units)
+ConstructionBase.constructorof(::Type{Diagnostic{name,S,T,units}}) where {name,S,T,units} = s -> Diagnostic(name, s, T, units)
+ConstructionBase.constructorof(::Type{Algebraic{name,S,T,units}}) where {name,S,T,units} = s -> Algebraic(name, s, T, units)
+ConstructionBase.constructorof(::Type{Delta{dname,name,S,T,units}}) where {dname,name,S,T,units} = s -> Delta(dname, name, s, T, units)
+==(var1::Var{N1,S1,T1,u1},var2::Var{N2,S2,T2,u2}) where {N1,N2,S1,S2,T1,T2,u1,u2} = (N1==N2) && (S1==S2) && (T1==T2) && (u1 == u2)
 varname(::Var{name}) where {name} = name
 varname(::Type{<:Var{name}}) where {name} = name
-vartype(::Var{name,T}) where {name,T} = T
-vartype(::Type{<:Var{name,T}}) where {name,T} = T
-vardims(var::Var{name,T,S}) where {name,T,S} = var.dim
+vartype(::Var{name,S,T}) where {name,S,T} = T
+vartype(::Type{<:Var{name,S,T}}) where {name,S,T} = T
+varunits(::Var{name,S,T,units}) where {name,S,T,units} = units
+varunits(::Type{<:Var{name,S,T,units}}) where {name,S,T,units} = units
+vardims(var::Var) = var.dim
 isprognostic(::T) where {T<:Var} = T <: Prognostic
 isalgebraic(::T) where {T<:Var} = T <: Algebraic
-isflux(::T) where {T<:Var} = T <: Flux
+isflux(::T) where {T<:Var} = T <: Delta
 isdiagnostic(::T) where {T<:Var} = T <: Diagnostic
-isongrid(::Var{name,T,S}) where {name,T,S} = S <: OnGrid
+isongrid(::Var{name,S}) where {name,S} = S <: OnGrid

src/Numerics/varstates.jl

@@ -6,14 +6,13 @@ prognostic state vector (`uproto`).
 """
 struct VarStates{names,griddvars,TU,TD,TV,DF,DG}
     uproto::TU # prognostic state vector prototype
-    grid::TD
+    grid::TD # model grid/discretization
     vars::NamedTuple{names,TV} # variable metadata
     diag::NamedTuple{names,DF} # off-grid non-prognostic variables
     griddiag::NamedTuple{griddvars,DG} # on-grid non-prognostic variables
-    gridcache::Dict{ClosedInterval{Int},TD} # grid cache; indices -> subgrid
 end
-@generated function getvar(::Val{name}, vs::VarStates{layers,griddvars,TU,TD,TV}, u::TU, du::Union{Nothing,TU}=nothing) where
-    {name,layers,griddvars,T,A,pax,TU<:ComponentVector{T,A,Tuple{Axis{pax}}},TD,TV}
+@generated function getvar(::Val{name}, vs::VarStates{layers,griddvars}, u, du=nothing) where {name,layers,griddvars}
+    pax = ComponentArrays.indexmap(first(ComponentArrays.getaxes(u)))
     dnames = map(n -> Symbol(:d,n), keys(pax))
     if name ∈ griddvars
         quote
@@ -91,7 +90,7 @@ function VarStates(vars::GroupedVars, D::Numerics.AbstractDiscretization, chunks
     progvars = map(group -> filter(isprognostic, group), vars)
     algvars = map(group -> filter(isalgebraic, group), vars)
     # create variables for gradients/fluxes
-    dpvars = map(group -> map(Flux, filter(var -> isalgebraic(var) || isprognostic(var), group)), vars)
+    dpvars = map(group -> map(Delta, filter(var -> isalgebraic(var) || isprognostic(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(...)
@@ -108,5 +107,5 @@ function VarStates(vars::GroupedVars, D::Numerics.AbstractDiscretization, chunks
     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))
-    VarStates(uproto, D, allvars, (;freediagstate...), (;griddiagstate...), Dict(1..length(D) => D))
+    VarStates(uproto, D, allvars, (;freediagstate...), (;griddiagstate...))
 end

src/Physics/Boundaries/Boundaries.jl

 module Boundaries
 
-import CryoGrid: BoundaryProcess, BoundaryStyle, Dirichlet, Neumann, Top
+import CryoGrid: SubSurfaceProcess, BoundaryProcess, BoundaryStyle, Dirichlet, Neumann, Top
 import CryoGrid: variables, boundaryvalue
 
 using CryoGrid.Numerics
@@ -8,83 +8,21 @@ using CryoGrid.Utils
 
 using ConstructionBase
 using Dates
-using Flatten
-using Interpolations
-using ModelParameters
-using Parameters
-using TimeSeries
 using Unitful
 
 import Flatten: flattenable
 
-export Constant, Periodic, Bias
-export BoundaryEffect, Damping
-include("effects.jl")
-export Forcing, TimeSeriesForcing, ForcingData
-include("forcing.jl")
-
-"""
-    struct Constant{S,T} <: BoundaryProcess
-
-Constant boundary condition (of any type/unit) specified by `value`.
-"""
-struct Constant{S,T} <: BoundaryProcess
-    value::T
-    Constant(::Type{S}, value::T) where {S<:BoundaryStyle,T} = new{S,T}(value)
-end
-ConstructionBase.constructorof(::Type{<:Constant{S}}) where {S} = value -> Constant(S,value)
-boundaryvalue(bc::Constant{S,T},l1,p2,l2,s1,s2) where {S,T} = bc.value
-
-BoundaryStyle(::Type{<:Constant{S}}) where {S} = S()
-
 """
-    struct Periodic{S,T} <: BoundaryProcess
+    BoundaryEffect
 
-Periodic boundary condition (of any type/unit) specified by `period`, `amplitude`, and `phaseshift`.
+Base type for boundary "effects" which modify boundary conditions based on some
+given parameterization.
 """
-struct Periodic{S,T} <: BoundaryProcess
-    period::Float"s"
-    amplitude::T
-    phaseshift::T
-    Periodic{S}(period::Q, amplitude::T=one(T), phaseshift::T=one(T)) where
-        {S<:BoundaryStyle,Q<:TimeQuantity,T} =
-        new{S,T}(uconvert(u"s",period) |> dustrip, amplitude, phaseshift)
-end
+abstract type BoundaryEffect end
 
-@inline boundaryvalue(bc::Periodic,l1,p2,l2,s1,s2) = bc.amplitude*sin(π*(1/bc.period)*t + bc.phaseshift)
+export BoundaryEffect
 
-BoundaryStyle(::Type{<:Periodic{S}}) where {S} = S()
-
-@with_kw struct Bias{P} <: BoundaryProcess
-    bias::P = Param(0.0)
-end
-@inline boundaryvalue(bc::Bias,l1,p2,l2,s1,s2) = bc.bias
-
-BoundaryStyle(::Type{<:Bias}) = Dirichlet()
-
-"""
-    struct CombinedBoundaryProcess{B1,B2,F,S} <: BoundaryProcess
-
-Represents a composition of two boundary processes, `B1` and `B2`, via an operator `F`.
-A typical use case is combining `Constant` with a forcing-driven boundary process to
-scale or shift the forcing.
-"""
-struct CombinedBoundaryProcess{B1,B2,F,S} <: BoundaryProcess
-    op::F
-    bc1::B1
-    bc2::B2
-    function CombinedBoundaryProcess(op::F, bc1::B1, bc2::B2) where {F,B1<:BoundaryProcess,B2<:BoundaryProcess}
-        @assert BoundaryStyle(bc1) == BoundaryStyle(bc2) "boundary condition styles (e.g. Dirichlet vs Neumann) must match"
-        new{B1,B2,F,typeof(BoundaryStyle(bc1))}(op,bc1,bc2)
-    end
-end
-@inline boundaryvalue(cbc::CombinedBoundaryProcess,l1,p2,l2,s1,s2) = cbc.op(cbc.bc1(l1,l2,p2,s1,s2), cbc.bc2(l1,l2,p2,s1,s2))
-variables(top::Top, cbc::CombinedBoundaryProcess) = tuplejoin(variables(top, cbc.bc1), variables(top, cbc.bc2))
-BoundaryStyle(::Type{CombinedBoundaryProcess{B1,B2,F,S}}) where {F,B1,B2,S} = S()
-# Overload arithmetic operators on boundary processes.
-Base.:+(bc1::BoundaryProcess, bc2::BoundaryProcess) = CombinedBoundaryProcess(+, bc1, bc2)
-Base.:-(bc1::BoundaryProcess, bc2::BoundaryProcess) = CombinedBoundaryProcess(-, bc1, bc2)
-Base.:*(bc1::BoundaryProcess, bc2::BoundaryProcess) = CombinedBoundaryProcess(*, bc1, bc2)
-Base.:/(bc1::BoundaryProcess, bc2::BoundaryProcess) = CombinedBoundaryProcess(/, bc1, bc2)
+export ConstantBC, PeriodicBC, Bias
+include("bc.jl")
 
 end

src/Physics/Boundaries/bc.jl

+"""
+    ConstantBC{P,S,T} <: BoundaryProcess{P}
+
+Constant boundary condition (of any type/unit) specified by `value`.
+"""
+struct ConstantBC{P,S,T} <: BoundaryProcess{P}
+    value::T
+    ConstantBC(::Type{P}, ::Type{S}, value::T) where {P<:SubSurfaceProcess,S<:BoundaryStyle,T} = new{P,S,T}(value)
+end
+ConstructionBase.constructorof(::Type{<:ConstantBC{P,S}}) where {P,S} = value -> ConstantBC(P, S, value)
+boundaryvalue(bc::ConstantBC,l1,p2,l2,s1,s2) = bc.value
+
+BoundaryStyle(::Type{<:ConstantBC{P,S}}) where {P,S} = S()
+
+"""
+    PeriodicBC{P,S,T1,T2,T3} <: BoundaryProcess{P}
+
+Periodic boundary condition (of any type/unit) specified by `period`, `amplitude`, and `phaseshift`.
+"""
+struct PeriodicBC{P,S,T1,T2,T3} <: BoundaryProcess{P}
+    period::T1
+    amplitude::T2
+    phaseshift::T3
+    PeriodicBC(::Type{P}, ::Type{S}, period::T1, amplitude::T2=one(T), phaseshift::T3=one(T)) where
+        {P<:SubSurfaceProcess,S<:BoundaryStyle,T1<:TimeQuantity,T2,T3} =
+        new{P,S,T1,T2,T3}(period, amplitude, phaseshift)
+end
+ConstructionBase.constructorof(::Type{<:PeriodicBC{P,S}}) where {P,S} = (args...) -> PeriodicBC(P, S, args...)
+@inline boundaryvalue(bc::PeriodicBC,l1,p2,l2,s1,s2) = bc.amplitude*sin(π*(1/bc.period)*t + bc.phaseshift)
+
+BoundaryStyle(::Type{<:PeriodicBC{P,S}}) where {P,S} = S()
+
+"""
+    Bias{P,Tb} <: BoundaryProcess{P}
+
+Boundary process which adds a constant shift/offset to the boundary condition.
+"""
+Base.@kwdef struct Bias{P,Tb} <: BoundaryProcess{P}
+    bias::Tb = Param(0.0)
+    Bias(::Type{P}, bias::Tb) where {P<:SubSurfaceProcess,Tb} = new{P,Tb}(bias)
+end
+@inline boundaryvalue(bc::Bias,l1,p2,l2,s1,s2) = bc.bias
+
+BoundaryStyle(::Type{<:Bias}) = Dirichlet()
+
+"""
+    CombinedBoundaryProcess{B1,B2,F,S} <: BoundaryProcess
+
+Represents a composition of two boundary processes, `B1` and `B2`, via an operator `F`.
+A typical use case is combining `ConstantBC` with a forcing-driven boundary process to
+scale or shift the forcing.
+"""
+struct CombinedBoundaryProcess{B1,B2,P,F,S} <: BoundaryProcess{P}
+    op::F
+    bc1::B1
+    bc2::B2
+    function CombinedBoundaryProcess(op::F, bc1::B1, bc2::B2) where {F,P,B1<:BoundaryProcess{P},B2<:BoundaryProcess{P}}
+        @assert BoundaryStyle(bc1) == BoundaryStyle(bc2) "boundary condition styles (e.g. Dirichlet vs Neumann) must match"
+        new{B1,B2,P,F,typeof(BoundaryStyle(bc1))}(op,bc1,bc2)
+    end
+end
+@inline boundaryvalue(cbc::CombinedBoundaryProcess{B1,B2},l1,p2,l2,s1,s2) where {B1<:BoundaryProcess,B2<:BoundaryProcess} = cbc.op(boundaryvalue(cbc.bc1,l1,p2,l2,s1,s2), boundaryvalue(cbc.bc2,l1,p2,l2,s1,s2))
+@inline boundaryvalue(cbc::CombinedBoundaryProcess{B1,B2},l1,p2,l2,s1,s2) where {B1,B2<:BoundaryProcess} = cbc.op(cbc.bc1, boundaryvalue(cbc.bc2,l1,p2,l2,s1,s2))
+@inline boundaryvalue(cbc::CombinedBoundaryProcess{B1,B2},l1,p2,l2,s1,s2) where {B1<:BoundaryProcess,B2} = cbc.op(boundaryvalue(cbc.bc1,l1,p2,l2,s1,s2), cbc.bc2)
+variables(top::Top, cbc::CombinedBoundaryProcess) = tuplejoin(variables(top, cbc.bc1), variables(top, cbc.bc2))
+BoundaryStyle(::Type{CombinedBoundaryProcess{B1,B2,P,F,S}}) where {F,P,B1,B2,S} = S()
+# Overload arithmetic operators on boundary processes.
+Base.:+(bc1::BoundaryProcess, bc2::BoundaryProcess) = CombinedBoundaryProcess(+, bc1, bc2)
+Base.:-(bc1::BoundaryProcess, bc2::BoundaryProcess) = CombinedBoundaryProcess(-, bc1, bc2)
+Base.:*(bc1::BoundaryProcess, bc2::BoundaryProcess) = CombinedBoundaryProcess(*, bc1, bc2)
+Base.:/(bc1::BoundaryProcess, bc2::BoundaryProcess) = CombinedBoundaryProcess(/, bc1, bc2)

src/Physics/Boundaries/effects.jl

-abstract type BoundaryEffect end
-"""
-    Damping{D,K} <: BoundaryEffect
-
-Generic implementation of bulk conductive damping at the boundary.
-"""
-@with_kw struct Damping{D,K} <: BoundaryEffect
-    depth::D = t -> 0.0 # function of t -> damping depth; defaults to zero function
-    k::K = Param(1.0, bounds=(0.0,Inf)) # conductivity of medium
-end
-function (damp::Damping)(u_top, u_sub, k_sub, Δsub, t)
-    let d_med = damp.depth(t), # length of medium
-        d_sub = Δsub / 2, # half length of upper grid cell
-        d = d_med + d_sub, # total flux distance (from grid center to "surface")
-        k_med = damp.k, # conductivity of damping medium (assumed uniform)
-        k_sub = k_sub, # conductivity at grid center
-        k = Numerics.harmonicmean(k_med, k_sub, d_med, d_sub);
-        -k*(u_sub - u_top)/d/Δsub # flux per unit volume
-    end
-end