Add simple impl of tabulated freeze curve

src/Numerics/Numerics.jl

@@ -12,7 +12,7 @@ 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
@@ -35,7 +35,7 @@ struct Cells <: GridSpec end
 abstract type Geometry end
 struct UnitVolume <: Geometry end
 
-export ∇
+export ∇, Tabulated
 include("math.jl")
 
 export Grid, cells, edges, indexmap, subgridinds, Δ, volume, area

src/Numerics/math.jl

@@ -132,6 +132,7 @@ 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)
 
+# Symbolic differentiation
 """
     ∇(f, dvar::Symbol)
 
@@ -170,3 +171,38 @@ function ∇(f, dvar::Symbol; choosefn=first, context_module=Numerics)
     ∇f = @RuntimeGeneratedFunction(context_module, ∇f_expr)
     return ∇f
 end
+
+# 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])
+    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
+    interp = interpolate(Tuple(knots), map(Base.splat(f), arggrid), Gridded(Linear()))
+    return interp
+end
+function ∇(f::AbstractInterpolation)
+    gradient(args...) = Interpolations.gradient(f, args...)
+    return gradient
+end

src/Physics/HeatConduction/soil/sfcc.jl

@@ -163,6 +163,28 @@ function (f::Westermann)(T,Tₘ,θres,θsat,θtot,δ)
         IfElse.ifelse(T<=Tₘ, θres - (θsat-θres)*(δ/(T-δ)), θtot)
     end
 end
+struct SFCCTable{F,I} <: SFCCFunction
+    f::F
+    f_tab::I
+end
+(f::SFCCTable)(args...) = f.f_tab(args...)
+"""
+    Tabulated(f::SFCCFunction, args...)
+
+Produces an `SFCCTable` function which is a tabulation of `f`.
+"""
+Numerics.Tabulated(f::SFCCFunction, args...) = SFCCTable(f, Numerics.tabulate(f, args...))
+"""
+    SFCC(f::SFCCTable, s::SFCCSolver=SFCCNewtonSolver())
+
+Constructs a SFCC from the precomputed `SFCCTable`. The derivative is generated using the
+`gradient` function provided by `Interpolations`.
+"""
+function SFCC(f::SFCCTable, s::SFCCSolver=SFCCNewtonSolver())
+    # we wrap ∇f with Base.splat here to avoid a weird issue with in-place splatting causing allocations
+    # when applied to runtime generated functions.
+    SFCC(f, Base.splat(∇(f.f_tab)), s)
+end
 
 """
 Specialized implementation of Newton's method with backtracking line search for resolving