Automatic Differentiation testing

src/Density/Density.jl

@@ -155,7 +155,7 @@ function param_info(s::Compressible_Density) # info about the struct
     return MaterialParamsInfo(; Equation = L"\rho = \rho_0\exp(\beta*(P-P_0))")
 end
 
-function (s::Compressible_Density{_T})(; P::_T=zero(_T), kwargs...) where {_T}
+function (s::Compressible_Density{_T})(; P=0e0, kwargs...) where {_T}
     if P isa Quantity
         @unpack_units ρ0, β, P0 = s
     else
@@ -199,14 +199,14 @@ function param_info(s::T_Density) # info about the struct
     return MaterialParamsInfo(; Equation = L"\rho = \rho_0*(1 - \alpha*(T-T_0))")
 end
 
-function (s::T_Density{_T})(; T::_T=zero(_T), kwargs...) where {_T}
+function (s::T_Density{_T})(; T = 0e0, kwargs...) where {_T}
     if T isa Quantity
         @unpack_units ρ0, α, T0 = s
     else
         @unpack_val   ρ0, α, T0 = s
     end
 
-    return @muladd ρ0 * (1.0 -  α * (T - T0))
+    return @muladd ρ0 * (1 -  α * (T - T0))
 end
 
 @inline (s::T_Density)(args)                = s(; args...)
@@ -248,7 +248,7 @@ function param_info(s::MeltDependent_Density) # info about the struct
 end
 
 # Calculation routines
-function (rho::MeltDependent_Density{_T})(; ϕ::_T=zero(_T), kwargs...) where {_T}
+function (rho::MeltDependent_Density{_T})(; ϕ=0e0, kwargs...) where {_T}
     ρsolid = compute_density(rho.ρsolid, kwargs)
     ρmelt  = compute_density(rho.ρmelt,  kwargs)
 

src/Energy/Conductivity.jl

@@ -147,7 +147,7 @@
 end
 
 # Calculation routine
-function (s::T_Conductivity_Whittington{_T})(; T::_T=zero(_T), kwargs...) where {_T}
+function (s::T_Conductivity_Whittington{_T})(; T=0e0, kwargs...) where {_T}
     if T isa Quantity
         @unpack_units a0, a1, b0, b1, c0, c1, molmass, Tcutoff, rho, d, e, f, g = s
     else
@@ -161,7 +161,7 @@
     end
 end
 
-function compute_conductivity(s::T_Conductivity_Whittington{_T}; T::_T=zero(_T)) where {_T}
+function compute_conductivity(s::T_Conductivity_Whittington{_T}; T=0e0) where {_T}
     return s(; T=T)
 end
 
@@ -253,7 +253,7 @@
 
 # Calculation routine
 function (s::T_Conductivity_Whittington_parameterised{_T})(;
-    T::_T=zero(_T), kwargs...
+    T=0e0, kwargs...
 ) where {_T}
     if T isa Quantity
         @unpack_units a, b, c, d, Ts = s
@@ -418,7 +418,7 @@
 )
 
 # Calculation routine
-function (s::TP_Conductivity{_T})(; P::_T=zero(_T), T::_T=zero(_T), kwargs...) where {_T}
+function (s::TP_Conductivity{_T})(; P=0e0, T=0e0, kwargs...) where {_T}
     if T isa Quantity
         @unpack_units a, b, c, d = s
     else

test/Project.toml

@@ -1,10 +1,12 @@
 [deps]
+AbstractDifferentiation = "c29ec348-61ec-40c8-8164-b8c60e9d9f3d"
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
+LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Unidecode = "967fb449-e509-55aa-8007-234b4096b967"
-LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"

test/test_AD.jl

@@ -0,0 +1,104 @@
+using Test, GeoParams
+import AbstractDifferentiation as AD, ReverseDiff, ForwardDiff
+
+T        = 1e3
+P        = 1e9
+backends = AD.ReverseDiffBackend(), AD.ForwardDiffBackend()
+pkg      = "ReverseDiff", "ForwardDiff"
+
+for (name, backend) in zip(pkg, backends)
+    @testset "$name compatibility with: density" begin
+        ρ = ConstantDensity()
+        @test AD.derivative(backend, x -> compute_density(ρ, (; T=x, P=1e9)), T)[1] == 0
+        @test AD.derivative(backend, x -> compute_density(ρ, (; T=1e3, P=x)), P)[1] == 0
+
+        ρ = PT_Density()
+        @test AD.derivative(backend, x -> compute_density(ρ, (; T=x, P=1e9)), T)[1] ≈ -0.087
+        @test AD.derivative(backend, x -> compute_density(ρ, (; T=1e3, P=x)), P)[1] == 2.9e-6
+
+        ρ = Compressible_Density()
+        @test AD.derivative(backend, x -> compute_density(ρ, (; T=1e3, P=x)), P)[1] ≈ 7.88301730e-6
+
+        ρ = MeltDependent_Density()
+        @test AD.derivative(backend, x -> compute_density(ρ, (; ϕ =x)), 0.2)[1] == -700
+
+        ρ = T_Density()
+        @test AD.derivative(backend, x -> compute_density(ρ, (; T=x, P=1e9)), T)[1] ≈ -0.087
+    end
+
+    @testset "$name compatibility with: heat capacity" begin
+        Cp = T_HeatCapacity_Whittington()
+        @test AD.derivative(backend, x -> compute_heatcapacity(Cp, (; T=x)), T)[1] ≈ 0.145639823
+
+        Cp = Latent_HeatCapacity(Q_L=500e3)
+        @test AD.derivative(backend, x -> compute_heatcapacity(Cp, (; T=x)), T)[1] == 0 
+    end
+
+    @testset "$name compatibility with: conductivity" begin
+        K = ConstantConductivity()
+        @test AD.derivative(backend, x -> compute_conductivity(K, (; T=x)), T)[1] == 0
+
+        K = T_Conductivity_Whittington()
+        @test AD.derivative(backend, x -> compute_conductivity(K, (; T=x)), T)[1] ≈ -0.00019522103
+
+        K = T_Conductivity_Whittington_parameterised()
+        @test AD.derivative(backend, x -> compute_conductivity(K, (; T=x)), T)[1] ≈ -0.00064766553
+
+        K = TP_Conductivity()
+        @test AD.derivative(backend, x -> compute_conductivity(K, (; T=x)), T)[1] ≈ -0.00040864570
+    end
+
+    @testset "$name compatibility with: Diffusion" begin
+        import GeoParams.Diffusion
+        # Define a linear viscous creep law ---------------------------------
+        diffusion_law = Diffusion.dry_anorthite_Rybacki_2006
+        p             = SetDiffusionCreep(diffusion_law; n = 1NoUnits)
+        args          = (; T=T, P=P)
+        TauII         = 1e6
+        @test AD.derivative(backend, x -> compute_εII(p, TauII, (; T=x, P=P)), T)[1] ≈ 5.7562812856e-33
+        @test AD.derivative(backend, x -> compute_εII(p, TauII, (; T=T, P=x)), P)[1] ≈ -2.8543543565e-40
+
+        εII           = compute_εII(p, TauII, args)
+        @test AD.derivative(backend, x -> compute_τII(p, εII, (; T=x, P=P)), T)[1] ≈ -58211.55812135427
+        @test AD.derivative(backend, x -> compute_τII(p, εII, (; T=T, P=x)), P)[1] ≈ 0.0028865235432076496
+    end
+
+    @testset "$name compatibility with: Dislocation" begin
+        import GeoParams.Dislocation
+        # Define a linear viscous creep law ---------------------------------
+        diffusion_law = Dislocation.dry_olivine_Hirth_2003
+        p             = SetDislocationCreep(diffusion_law; n = 1NoUnits)
+        args          = (; T=T, P=P)
+        TauII         = 1e6
+        @test AD.derivative(backend, x -> compute_εII(p, TauII, (; T=x, P=P)), T)[1] ≈ 4.1522654949e-40
+        @test AD.derivative(backend, x -> compute_εII(p, TauII, (; T=T, P=x)), P)[1] ≈ -1.0685977376e-47
+
+        εII           = compute_εII(p, TauII, args)
+        @test AD.derivative(backend, x -> compute_τII(p, εII, (; T=x, P=P)), T)[1] ≈ -65427.866979
+        @test AD.derivative(backend, x -> compute_τII(p, εII, (; T=T, P=x)), P)[1] ≈ 0.00168380540
+    end
+
+    @testset "$name compatibility with: CompositeRheology" begin
+        # Define a range of rheological components
+        v1       = SetDiffusionCreep(Diffusion.dry_anorthite_Rybacki_2006)
+        v2       = SetDislocationCreep(Dislocation.dry_anorthite_Rybacki_2006)
+        v3       = LinearViscous()
+        el       = ConstantElasticity()
+        # composite rheology
+        c1       = CompositeRheology(v1, v2, el)
+        c2       = CompositeRheology(v3, el)       # composite rheology
+        # arguments
+        args     = (T=900.0, d=100e-6, τII_old=1e6, dt=1e8)
+        εII, τII = 1e-12, 2e6
+        # test non-linear rheology
+        @test AD.derivative(backend, x -> compute_τII(c1, εII, (;T=x, P=P, d=100e-6, τII_old=1e6, dt=1e8)), T)[1] ≈ -1.074788789
+        @test AD.derivative(backend, x -> compute_τII(c1, εII, (;T=T, P=x, d=100e-6, τII_old=1e6, dt=1e8)), P)[1] ≈  4.73352298e-8
+        @test AD.derivative(backend, x -> compute_εII(c1, τII, (;T=x, P=P, d=100e-6, τII_old=1e6, dt=1e8)), T)[1] ≈ 1.17780761876e-20
+        @test AD.derivative(backend, x -> compute_εII(c1, τII, (;T=T, P=x, d=100e-6, τII_old=1e6, dt=1e8)), P)[1] ≈ -5.8045331761e-28
+        # test linear rheology
+        @test iszero(AD.derivative(backend, x -> compute_τII(c2, εII, (;T=x, P=P, d=100e-6, τII_old=1e6, dt=1e8)), T)[1])
+        @test iszero(AD.derivative(backend, x -> compute_τII(c2, εII, (;T=T, P=x, d=100e-6, τII_old=1e6, dt=1e8)), P)[1])
+        @test iszero(AD.derivative(backend, x -> compute_εII(c2, τII, (;T=x, P=P, d=100e-6, τII_old=1e6, dt=1e8)), T)[1])
+        @test iszero(AD.derivative(backend, x -> compute_εII(c2, τII, (;T=T, P=x, d=100e-6, τII_old=1e6, dt=1e8)), P)[1])
+    end
+end

test/test_Energy.jl

@@ -199,7 +199,7 @@ using StaticArrays
     @test NumValue(cond.k) ≈ 3.8194500000000007
 
     @test compute_conductivity(cond; T=100.0) ≈ 3.8194500000000007 # compute
-
+    
     # Temperature-dependent conductivity
     # dimensional
     T = collect(250e0:100:1250e0)

test/test_TensorAlgebra.jl

@@ -70,6 +70,7 @@ end
         3 4 5
     ]
 end
+
 @testset "Tensor algebra" begin
     # J2 TENSOR INVARIANT: 2D TESTS 
     τ_xx, τ_yy, τ_xy = 1.0, 2.0, 3.0