diff --git a/src/Fronts.jl b/src/Fronts.jl index cde98ce0..43a2621f 100644 --- a/src/Fronts.jl +++ b/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 diff --git a/src/PorousModels/PorousModels.jl b/src/PorousModels/PorousModels.jl index 907b95a0..5ba91fe2 100644 --- a/src/PorousModels/PorousModels.jl +++ b/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 diff --git a/src/_Diff.jl b/src/_Diff.jl deleted file mode 100644 index 526c2d6e..00000000 --- a/src/_Diff.jl +++ /dev/null @@ -1,28 +0,0 @@ -module _Diff - -import AbstractDifferentiation -import ForwardDiff - -@inline function derivative(f, x::Real) - return only(AbstractDifferentiation.derivative(AbstractDifferentiation.ForwardDiffBackend(), - f, - x)) -end - -@inline function value_and_derivative(f, x::Real) - a, b = AbstractDifferentiation.value_and_derivative(AbstractDifferentiation.ForwardDiffBackend(), - f, - x) - return a, only(b) -end - -@inline function value_and_derivatives(f, x::Real) - a, b, c = AbstractDifferentiation.value_derivative_and_second_derivative(AbstractDifferentiation.ForwardDiffBackend(), - f, - x) - return a, only(b), only(c) -end - -export derivative, value_and_derivative, value_and_derivatives - -end diff --git a/src/finite.jl b/src/finite.jl index d019e6af..46970b43 100644 --- a/src/finite.jl +++ b/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 diff --git a/src/odes.jl b/src/odes.jl index 424997db..d9771092 100644 --- a/src/odes.jl +++ b/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) diff --git a/test/runtests.jl b/test/runtests.jl index 4159c05c..e6bb189f 100644 --- a/test/runtests.jl +++ b/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") diff --git a/test/test_Diff.jl b/test/test_Diff.jl deleted file mode 100644 index 47652ad7..00000000 --- a/test/test_Diff.jl +++ /dev/null @@ -1,26 +0,0 @@ -@testset "_Diff" 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(D, θ) - @test D_ == D(θ) - @test dD_dθ == derivative(D, θ) == ForwardDiff.derivative(D, θ) - - D_, dD_dθ, d²D_dθ² = value_and_derivatives(D, θ) - @test D_ == D(θ) - @test dD_dθ == derivative(D, θ) == ForwardDiff.derivative(D, θ) - @test d²D_dθ² == derivative(θ -> derivative(D, θ), θ) == - ForwardDiff.derivative(θ -> ForwardDiff.derivative(D, θ), θ) - end -end diff --git a/test/test_PorousModels.jl b/test/test_PorousModels.jl index 76b197ca..33d2a277 100644 --- a/test/test_PorousModels.jl +++ b/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) diff --git a/test/test_differentiation.jl b/test/test_differentiation.jl new file mode 100644 index 00000000..75fa9b1e --- /dev/null +++ b/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