Add analytical derivatives for SWRCs

Author: bgroenks96

Project.toml

@@ -23,6 +23,7 @@ ForwardDiff = "0.10, 1"
 IfElse = "0.1"
 Interpolations = "0.14, 0.15, 0.16"
 IntervalSets = "0.7"
+QuadGK = "2"
 NLsolve = "4"
 NonlinearSolve = "1.1, 2, 3, 4"
 RecipesBase = "1"
@@ -36,9 +37,10 @@ julia = "1.6"
 [extras]
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
+QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Turing = "fce5fe82-541a-59a6-adf8-730c64b5f9a0"
 
 [targets]
-test = ["Test", "BenchmarkTools", "NonlinearSolve", "Optim", "Turing"]
+test = ["Test", "BenchmarkTools", "NonlinearSolve", "QuadGK", "Optim", "Turing"]

src/swrc.jl

@@ -18,11 +18,17 @@ abstract type SWRC end
 """
     Base.inv(f::SWRC)
 
-Returns a function `(args...) -> f(inv, args...)` which, when implemented, evaluates the inverse of the
-soil water retention curve.
+Returns a function that evaluates the inverse of the soil water retention curve, `f`.
 """
 Base.inv(f::SWRC) = (args...; kwargs...) -> f(inv, args...; kwargs...)
 
+"""
+    derivative(f::SWRC)
+
+Returns a function that evaluates the derivative `∂f∂ψ` of the soil water retention curve, `f`.
+"""
+derivative(f::SWRC) = (args...; kwargs...) -> f(derivative, args..., kwargs...)
+
 """
     SoilWaterVolume{Tρw,Tθres,Tθsat}
 
@@ -78,9 +84,10 @@ function (f::VanGenuchten)(ψ; θsat=f.vol.θsat, θres=f.vol.θres, α=f.α, n=
     let m = 1-1/n,
         x = 1 + (-α*ψ)^n;
         # @assert isnan(x) || x > zero(x) "van Genuchten violated constraint $x > 0 with input $ψ and parameters θsat=$θsat, θres=$θres, α=$α, n=$n"
-        IfElse.ifelse(ψ <= zero(ψ), θres + (θsat - θres)*x^(-m), θsat*one(ψ))
+        IfElse.ifelse(ψ < zero(ψ), θres + (θsat - θres)*x^(-m), θsat*one(ψ))
     end
 end
+
 function (f::VanGenuchten)(
     ::typeof(inv),
     θ;
@@ -94,6 +101,21 @@ function (f::VanGenuchten)(
     end
 end
 
+function (f::VanGenuchten)(
+    ::typeof(derivative),
+    ψ;
+    θsat=f.vol.θsat,
+    θres=f.vol.θres,
+    α=f.α,
+    n=f.n
+)
+    let m = 1-1/n,
+        x = 1 + (-α*ψ)^n;
+        # @assert isnan(x) || x > zero(x) "van Genuchten violated constraint $x > 0 with input $ψ and parameters θsat=$θsat, θres=$θres, α=$α, n=$n"
+        IfElse.ifelse(ψ < zero(ψ), α*n*m*(θsat - θres)*x^(-m-1)*(-α*ψ)^(n-1), 0)
+    end
+end
+
 inflectionpoint(f::VanGenuchten; α=f.α, n=f.n) = -1/α*((n - 1) / n)^(1/n)
 
 # Preset Van Genucthen curves from soil types;
@@ -121,6 +143,7 @@ Base.@kwdef struct BrooksCorey{Tvol,Tψₛ,Tλ} <: SWRC
     ψₛ::Tψₛ = 0.01u"m"
     λ::Tλ = 0.2 # domain: (0,Inf)
 end
+
 function (f::BrooksCorey)(
     ψ;
     θsat=f.vol.θsat,
@@ -130,6 +153,7 @@ function (f::BrooksCorey)(
 )
     IfElse.ifelse(ψ < -ψₛ, θres + (θsat - θres)*(-ψₛ / ψ)^λ, θsat*one(ψ))
 end
+
 function (f::BrooksCorey)(
     ::typeof(inv),
     θ;
@@ -140,4 +164,16 @@ function (f::BrooksCorey)(
 )
     IfElse.ifelse(θ < θsat, -ψₛ*((θ - θres)/(θsat - θres))^(-1/λ), -ψₛ*one(θ))
 end
+
+function (f::BrooksCorey)(
+    ::typeof(derivative),
+    ψ;
+    θsat=f.vol.θsat,
+    θres=f.vol.θres,
+    ψₛ=f.ψₛ,
+    λ=f.λ
+)
+    IfElse.ifelse(ψ < -ψₛ, (θsat - θres)*λ*(-ψₛ / ψ)^(λ-1)*(ψₛ*ψ^-2), θsat*one(ψ))
+end
+
 inflectionpoint(f::BrooksCorey; ψₛ=f.ψₛ) = ψₛ

test/env/Project.toml

@@ -4,6 +4,7 @@ ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 FreezeCurves = "71e4ad71-e4f2-45a3-aa0a-91ffaa9676be"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
+QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Turing = "fce5fe82-541a-59a6-adf8-730c64b5f9a0"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"

test/swrc_tests.jl

@@ -1,4 +1,5 @@
 using FreezeCurves
+using QuadGK
 using Test
 using Unitful
 
@@ -16,6 +17,12 @@ using Unitful
             # check inverse
             @test isapprox(Base.inv(f)(θw; θsat, θres, α, n), -0.1u"m", atol=1e-6u"m")
             @test isapprox(Base.inv(f)(θsat; θsat, θres, α, n), 0.0u"m", atol=1e-6u"m")
+            # check derivative using fundamental theorem of calculus
+            a = -10.0u"m"
+            b = -0.01u"m"
+            ∂f∂ψ = FreezeCurves.derivative(f)
+            f̃, _ = quadgk(∂f∂ψ, a, b)
+            @test f̃ ≈ f(b) - f(a)
         end
     end
     @testset "VanGenuchten (units)" begin
@@ -37,6 +44,12 @@ using Unitful
             # check inverse
             @test isapprox(Base.inv(f)(θw; θsat, θres, ψₛ, λ), -0.5u"m", atol=1e-6u"m")
             @test isapprox(Base.inv(f)(θsat; θsat, θres, ψₛ, λ), -ψₛ, atol=1e-6u"m")
+            # check derivative using fundamental theorem of calculus
+            a = -10.0u"m"
+            b = -0.01u"m"
+            ∂f∂ψ = FreezeCurves.derivative(f)
+            f̃, _ = quadgk(∂f∂ψ, a, b)
+            @test f̃ ≈ f(b) - f(a)
         end
     end
 end