src/Numerics/init.jl

+abstract type VarInit{varname} end
+Numerics.varname(::VarInit{varname}) where {varname} = varname
+"""
+    init!(x::AbstractVector, init::VarInit, args...)
+
+Initializes state variable `x` using initializer `init` and implementation-specific additional
+arguments.
+"""
+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)
+end
+init!(x::AbstractVector, init::ConstantInit) = @. x = init.value
+"""
+    InterpInit{varname,P,I,E} <: VarInit
+
+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
+by `Interpolations.jl`.
+"""
+struct InterpInit{varname,P,I,E} <: VarInit{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)
+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)
+end
+"""
+    initializer(varname::Symbol, value::T) where {T} => ConstantInit
+    initializer(varname::Symbol, profile::Profile, interp=Linear(), extrap=Flat()) => InterpInit
+
+Convenience constructor for `VarInit`s 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)

src/Common/Numerics/math.jl → src/Numerics/math.jl

@@ -25,21 +25,22 @@ function lineardiffusion!(∂y::AbstractVector, x::AbstractVector, Δ::AbstractV
     end
 end
 
-"""
-    nonlineardiffusion!(∂y::AbstractVector, x::AbstractVector, Δx::AbstractVector, k::AbstractVector, Δk::AbstractArray)
 
-Second order Laplacian with non-linear diffusion function, k.
-"""
-function nonlineardiffusion!(∂y::AbstractVector, x::AbstractVector, Δx::AbstractVector, k::AbstractVector, Δk::AbstractArray)
-    @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]);
-        @. ∂y = (k₂*(x₃ - x₂)/Δx₂ - k₁*(x₂-x₁)/Δx₁)/Δk
-    end
+@propagate_inbounds function _nonlineardiffusion!(∂y, x₁, x₂, x₃, k₁, k₂, Δx₁, Δx₂, Δk)
+    @. ∂y = (k₂*(x₃ - x₂)/Δx₂ - k₁*(x₂-x₁)/Δx₁)/Δk
+end
+@propagate_inbounds function _nonlineardiffusion!(
+    ∂y::AbstractVector{Float64},
+    x₁::AbstractVector{Float64},
+    x₂::AbstractVector{Float64},
+    x₃::AbstractVector{Float64},
+    k₁::AbstractVector{Float64},
+    k₂::AbstractVector{Float64},
+    Δx₁::AbstractVector{Float64},
+    Δx₂::AbstractVector{Float64},
+    Δk::AbstractVector{Float64},
+)
+    @turbo @. ∂y = (k₂*(x₃ - x₂)/Δx₂ - k₁*(x₂-x₁)/Δx₁)/Δk
 end
 
 """
@@ -51,15 +52,9 @@ end
         Δk::AbstractVector{Float64}
     )
 
-Second order Laplacian with non-linear diffusion function, k. Accelerated using `LoopVectorization.@turbo` for `Float64` vectors.
+Second order Laplacian with non-linear diffusion operator, `k`. Accelerated using `LoopVectorization.@turbo` for `Float64` vectors.
 """
-function nonlineardiffusion!(
-    ∂y::AbstractVector{Float64},
-    x::AbstractVector{Float64}, 
-    Δx::AbstractVector{Float64},
-    k::AbstractVector{Float64},
-    Δk::AbstractVector{Float64}
-)
+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]),
@@ -67,7 +62,28 @@ function nonlineardiffusion!(
         k₂ = (@view k[2:end]),
         Δx₁ = (@view Δx[1:end-1]),
         Δx₂ = (@view Δx[2:end]);
-        @turbo @. ∂y = (k₂*(x₃ - x₂)/Δx₂ - k₁*(x₂-x₁)/Δx₁)/Δk
+        _nonlineardiffusion!(∂y, x₁, x₂, x₃, k₁, k₂, Δx₁, Δx₂, Δk)
+    end
+end
+
+"""
+    harmonicmean(x₁, x₂, w₁, w₂)
+
+Simple weighted harmonic mean of two values, x₁ and x₂.
+"""
+harmonicmean(x₁, x₂, w₁, w₂) = (w₁ + w₂) / (w₁*x₁^-1 + w₂*x₂^-1)
+"""
+    harmonicmean!(h::AbstractVector, x::AbstractVector, w::AbstractVector)
+
+Vectorized harmonic mean of elements in `x` with weights `w`. Output is stored in `h`,
+which should have size `length(x)-1`.
+"""
+function harmonicmean!(h::AbstractVector, x::AbstractVector, w::AbstractVector)
+    @inbounds let x₁ = (@view x[1:end-1]),
+        x₂ = (@view x[2:end]),
+        w₁ = (@view w[1:end-1]),
+        w₂ = (@view w[2:end]);
+        @. h = harmonicmean(x₁,x₂,w₁,w₂)
     end
 end
 
@@ -158,5 +174,3 @@ function ∇(f, dvar::Symbol; choosefn=first, context_module=Numerics)
     ∇f = @RuntimeGeneratedFunction(context_module, ∇f_expr)
     return ∇f
 end
-
-export ∇

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/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

src/Processes/SEB/SEB.jl → src/Physics/SEB/SEB.jl

 module SEB
 
 using ..HeatConduction: Heat
-using ..Processes
-using CryoGrid.Forcings
-using CryoGrid.Interface
+using ..Physics
+using ..Boundaries
 using CryoGrid.Layers: Soil
 using CryoGrid.Numerics
 using CryoGrid.Utils
 
 using Parameters
 
-import CryoGrid.Interface: BoundaryStyle, initialcondition!, variables
+import CryoGrid: BoundaryProcess, BoundaryStyle, Neumann, Top
+import CryoGrid: initialcondition!, variables
 
 @with_kw struct SEBParams <: Params
     # surface properties --> should be associated with the Stratigraphy and maybe made state variables
@@ -31,7 +31,7 @@ import CryoGrid.Interface: BoundaryStyle, initialcondition!, variables
     cₐ::Float"J/(m^3*K)" = 1005.7xu"J/(kg*K)" * ρₐ  # volumetric heat capacity of dry air at standard pressure and 0°C [J/(m^3*K)]
 end
 
-struct SurfaceEnergyBalance{F} <: BoundaryProcess{Heat}
+struct SurfaceEnergyBalance{F} <: BoundaryProcess
     forcing::F
     sebparams::SEBParams
     function SurfaceEnergyBalance(