diff --git a/src/Density/Density.jl b/src/Density/Density.jl index bec89599..1b3b6500 100644 --- a/src/Density/Density.jl +++ b/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) diff --git a/src/Energy/Conductivity.jl b/src/Energy/Conductivity.jl index 42b04cbb..2b9847b5 100644 --- a/src/Energy/Conductivity.jl +++ b/src/Energy/Conductivity.jl @@ -147,7 +147,7 @@ function param_info(s::T_Conductivity_Whittington) # info about the structwhere 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 @@ function (s::T_Conductivity_Whittington{_T})(; T::_T=zero(_T), kwargs...) where 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 @@ end # 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 @@ TP_Conductivity_info = Dict( ) # 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 diff --git a/test/Project.toml b/test/Project.toml index 4432c3a8..c23eaddf 100644 --- a/test/Project.toml +++ b/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" diff --git a/test/test_AD.jl b/test/test_AD.jl new file mode 100644 index 00000000..4c3381b1 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/test/test_Energy.jl b/test/test_Energy.jl index 470366a4..07d9bbbc 100644 --- a/test/test_Energy.jl +++ b/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) diff --git a/test/test_TensorAlgebra.jl b/test/test_TensorAlgebra.jl index c67431f7..d6dff0f9 100644 --- a/test/test_TensorAlgebra.jl +++ b/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