Remove _Diff module in favor of AbstractDifferentiation

src/Fronts.jl

@@ -1,8 +1,5 @@
 module Fronts
 
-include("_Diff.jl")
-using ._Diff: derivative, value_and_derivative, value_and_derivatives
-
 include("_Rootfinding.jl")
 using ._Rootfinding: bracket_bisect
 
@@ -11,6 +8,9 @@ using ._Chebyshev: chebdif
 
 using LinearAlgebra: Diagonal, Tridiagonal, SingularException
 
+using ForwardDiff: derivative
+using AbstractDifferentiation: value_and_derivative,
+                               value_derivative_and_second_derivative, ForwardDiffBackend
 using ArgCheck: @argcheck
 using StaticArrays: @SVector, @SMatrix
 using RecursiveArrayTools: ArrayPartition

src/PorousModels/PorousModels.jl

@@ -1,8 +1,8 @@
 module PorousModels
 
-using .._Diff: derivative
 import ..Fronts: DiffusionEquation
 
+using ForwardDiff: derivative
 using NaNMath: pow
 using ArgCheck: @argcheck
 

src/finite.jl

@@ -457,6 +457,6 @@ d_dr(sol::FiniteSolution, r, t) = derivative(r -> sol(r, t), r)
 d_dt(sol::FiniteSolution, r, t) = derivative(t -> sol(r, t), t)
 
 function flux(sol::FiniteSolution, r, t)
-    val, d_dr = value_and_derivative(r -> sol(r, t), r)
+    val, (d_dr,) = value_and_derivative(ForwardDiffBackend(), r -> sol(r, t), r)
     return -conductivity(sol.prob.eq, val) * d_dr
 end

src/odes.jl

@@ -12,15 +12,18 @@ See also: [`DifferentialEquations`](https://diffeq.sciml.ai/stable/), [`StaticAr
 function boltzmann(eq::DiffusionEquation{1})
     let K = u -> conductivity(eq, u), C = u -> capacity(eq, u)
         function f((u, du_do), ::SciMLBase.NullParameters, o)
-            K_, dK_du = value_and_derivative(K, u)
+            K_, (dK_du,) = value_and_derivative(ForwardDiffBackend(), K, u)
 
             d²u_do² = -((C(u) * o / 2 + dK_du * du_do) / K_) * du_do
 
             return @SVector [du_do, d²u_do²]
         end
         function jac((u, du_do), ::SciMLBase.NullParameters, o)
-            K_, dK_du, d²K_du² = value_and_derivatives(K, u)
-            C_, dC_du = value_and_derivative(C, u)
+            K_, (dK_du,), (d²K_du²,) = value_derivative_and_second_derivative(
+                ForwardDiffBackend(),
+                K,
+                u)
+            C_, (dC_du,) = value_and_derivative(ForwardDiffBackend(), C, u)
 
             j21 = -du_do * (K_ * (2 * d²K_du² * du_do + dC_du * o) -
                    dK_du * (C_ * o + 2 * dK_du * du_do)) / (2K_^2)
@@ -37,15 +40,18 @@ function boltzmann(eq::DiffusionEquation{m}) where {m}
     @assert m in 2:3
     let K = u -> conductivity(eq, u), C = u -> capacity(eq, u), k = m - 1
         function f((u, du_do), ::SciMLBase.NullParameters, o)
-            K_, dK_du = value_and_derivative(K, u)
+            K_, (dK_du,) = value_and_derivative(ForwardDiffBackend(), K, u)
 
             d²u_do² = -((C(u) * o / 2 + dK_du * du_do) / K_ + k / o) * du_do
 
             return @SVector [du_do, d²u_do²]
         end
         function jac((u, du_do), ::SciMLBase.NullParameters, o)
-            K_, dK_du, d²K_du² = value_and_derivatives(K, u)
-            C_, dC_du = value_and_derivative(C, u)
+            K_, (dK_du,), (d²K_du²,) = value_derivative_and_second_derivative(
+                ForwardDiffBackend(),
+                K,
+                u)
+            C_, (dC_du,) = value_and_derivative(ForwardDiffBackend(), C, u)
 
             j21 = -du_do * (K_ * (2 * d²K_du² * du_do + dC_du * o) -
                    dK_du * (C_ * o + 2 * dK_du * du_do)) / (2K_^2)

test/runtests.jl

@@ -1,11 +1,12 @@
 using Fronts
-using Fronts._Diff
 using Fronts.PorousModels
 using Fronts.ParamEstim
 using Fronts.SciMLBase: NullParameters
 using Test
 
 import ForwardDiff
+using AbstractDifferentiation: derivative, second_derivative, value_and_derivative,
+                               value_derivative_and_second_derivative, ForwardDiffBackend
 import NaNMath
 using NumericalIntegration
 using OrdinaryDiffEq: ODEFunction, ODEProblem
@@ -14,8 +15,8 @@ using StaticArrays: @SVector, SVector
 using Plots: plot
 
 @testset "Fronts.jl" begin
-    include("test_Diff.jl")
     include("test_PorousModels.jl")
+    include("test_differentiation.jl")
     include("test_boltzmann.jl")
     include("test_dirichlet.jl")
     include("test_neumann.jl")

test/test_PorousModels.jl

@@ -14,7 +14,7 @@
         htest = -10.0
 
         @test (@inferred hθ(bc, @inferred θh(bc, htest))) ≈ htest
-        @test (@inferred Ch(bc, htest)) ≈ derivative(h -> θh(bc, h), htest)
+        @test (@inferred Ch(bc, htest)) ≈ ForwardDiff.derivative(h -> θh(bc, h), htest)
         @test (@inferred Cθ(bc, @inferred θh(bc, htest))) ≈ @inferred Ch(bc, htest)
         @test (@inferred Kθ(bc, @inferred θh(bc, htest))) ≈ @inferred Kh(bc, htest)
 
@@ -47,7 +47,7 @@
         @test Kh(vg, htest) ≈ 9.430485870291618e-9
 
         @test (@inferred hθ(vg, @inferred θh(vg, htest))) ≈ htest
-        @test (@inferred Ch(vg, htest)) ≈ derivative(h -> θh(vg, h), htest)
+        @test (@inferred Ch(vg, htest)) ≈ ForwardDiff.derivative(h -> θh(vg, h), htest)
         @test (@inferred Cθ(vg, @inferred θh(vg, htest))) ≈ @inferred Ch(vg, htest)
         @test (@inferred Kθ(vg, @inferred θh(vg, htest))) ≈ @inferred Kh(vg, htest)
 

test/test_differentiation.jl

@@ -0,0 +1,31 @@
+@testset "differentiation" begin
+    # HF135 membrane, Van Genuchten model
+    # Data from Buser (PhD thesis, 2016)
+    # http://hdl.handle.net/1773/38064
+    θr = 0.0473
+    θs = 0.945
+    k = 5.50e-13  # m^2
+    α = 0.2555  # 1/m
+    n = 2.3521
+
+    model = VanGenuchten(n = n, α = α, k = k, θr = θr, θs = θs)
+
+    D(θ) = Dθ(model, θ)
+
+    for θ in range(θr + eps(θr), θs - eps(θs), length = 10)
+        D_, (dD_dθ,) = value_and_derivative(ForwardDiffBackend(), D, θ)
+        @test D_ == D(θ)
+        @test dD_dθ == only(derivative(ForwardDiffBackend(), D, θ)) ==
+              ForwardDiff.derivative(D, θ)
+
+        D_, (dD_dθ,), (d²D_dθ²,) = value_derivative_and_second_derivative(
+            ForwardDiffBackend(),
+            D,
+            θ)
+        @test D_ == D(θ)
+        @test dD_dθ == only(derivative(ForwardDiffBackend(), D, θ)) ==
+              ForwardDiff.derivative(D, θ)
+        @test d²D_dθ² == only(second_derivative(ForwardDiffBackend(), D, θ)) ==
+              ForwardDiff.derivative(θ -> ForwardDiff.derivative(D, θ), θ)
+    end
+end