diff --git a/test/test_CreepLaw.jl b/test/test_CreepLaw.jl index c72c9995..284ec753 100644 --- a/test/test_CreepLaw.jl +++ b/test/test_CreepLaw.jl @@ -1,7 +1,6 @@ using Test, Statistics using GeoParams - @testset "CreepLaw" begin #Make sure structs are isbits @@ -65,15 +64,17 @@ using GeoParams τII = 1e0 εII = 1e0 - @test compute_τII(x2, τII) == η0 - @test compute_εII(x2, εII) == 0.36746619407366893 + + (τII / η0)^(1/n) + @test compute_τII(x2, τII) == η0 * εII^n + @test compute_εII(x2, εII) ≈ (τII / η0)^(1/n) @test compute_viscosity_εII(x2, τII, (;)) == 5e0 @test compute_viscosity_τII(x2, εII, (;)) == 1.3606693841876538 x2 = PowerlawViscous() x2 = nondimensionalize(x2, CharUnits_GEO) - @test NumValue(x2.ε0) == 0.001 # powerlaw + @test NumValue(x2.ε0) == 0.001 # powerlaw # Given strainrate compute_τII(x1, 1e-13 / s, args) @@ -160,33 +161,26 @@ using GeoParams @test ε ≈ [2.553516169022911e-6; 5.107032338045822e-6] # vector input # With vector as input - T = Vector(800.0:1400)*K - η_basalt = zeros(size(T))*Pas + T = Vector(800.0:1400)*K + η_basalt = zeros(size(T))*Pas η_rhyolite = zeros(size(T))*Pas - x_basalt = LinearMeltViscosity() x_rhyolite = LinearMeltViscosity(A = -8.1590, B = 2.4050e+04K, T0 = -430.9606K) # Rhyolite for i in eachindex(T) - args_D = (; T = T[i]) - η_basalt[i] = compute_viscosity_τII(x_basalt, 1e6Pa, (; T=T[i])) - η_rhyolite[i] = compute_viscosity_τII(x_rhyolite, 1e6Pa, (; T=T[i])) + args_D = (; T = T[i]) + η_basalt[i] = compute_viscosity_τII(x_basalt, 1e6Pa, (; T=T[i])) + η_rhyolite[i] = compute_viscosity_τII(x_rhyolite, 1e6Pa, (; T=T[i])) end - #using Plots - #plot(ustrip.(T) .- 273.15, log10.(ustrip.(η_basalt)), xlabel="T [C]", ylabel="log10(η [Pas])", label="basalt") - #plot!(ustrip.(T) .- 273.15, log10.(ustrip.(η_rhyolite)), label="rhyolite") - - args_D=(;T=1000K) - η_basalt1 = compute_viscosity_τII(x_basalt, 1e6Pa, args_D) - @test η_basalt1 ≈ 5.2471745814062805e9Pa*s - - η_rhyolite1 = compute_viscosity_τII(x_rhyolite, 1e6Pa, args_D) + args_D = (;T=1000K) + η_basalt1 = compute_viscosity_τII(x_basalt, 1e6Pa, args_D) + @test η_basalt1 ≈ 5.2471745814062805e9Pa*s + η_rhyolite1 = compute_viscosity_τII(x_rhyolite, 1e6Pa, args_D) @test η_rhyolite1 ≈ 4.445205243727607e8Pa*s # ------------------------------------------------------------------- # Test effective viscosity of partially molten rocks ----------------- - x1_D =ViscosityPartialMelt_Costa_etal_2009(η=LinearMeltViscosity(A = -8.1590, B = 2.4050e+04K, T0 = -430.9606K)) # Rhyolite @test isDimensional(x1_D) == true x1_ND = nondimensionalize(x1_D, CharUnits_GEO) # check that we can nondimensionalize all entries within the struct @@ -199,11 +193,11 @@ using GeoParams args_D = (; T = 700K, ϕ=0.5) args_ND = (; T = nondimensionalize(700K, CharUnits_GEO), ϕ=0.5) - ε_D = compute_εII(x1_D, 1e6Pa, args_D) - @test ε_D ≈ 1.5692922423522311e-10 / s # dimensional input - ε_ND = compute_εII(x1_ND, nondimensionalize(1e6Pa, CharUnits_GEO), args_ND) + ε_D = compute_εII(x1_D, 1e6Pa, args_D) + @test ε_D ≈ 1.5692922423522311e-10 / s # dimensional input + ε_ND = compute_εII(x1_ND, nondimensionalize(1e6Pa, CharUnits_GEO), args_ND) @test ε_ND ≈ 156.92922423522444 # non-dimensional - @test ε_ND*CharUnits_GEO.strainrate ≈ ε_D # check that non-dimensiional computations give the same results + @test ε_ND * CharUnits_GEO.strainrate ≈ ε_D # check that non-dimensiional computations give the same results ε = [0.0; 0.0] τ = [1e6; 2e6] @@ -211,59 +205,51 @@ using GeoParams compute_εII!(ε, x1_D, τ, args_ND1) @test ε ≈ [0.0011248074556411618; 0.0022496149112823235] # vector input - # Given strainrate @test compute_τII(x1_D, 1e-13 / s, args_D) ≈ 643.9044043803415Pa # dimensional input - @test compute_τII(x1_ND, 1e-13, args_ND) ≈ 5.64401178083053e-17 # non-dimensional + @test compute_τII(x1_ND, 1e-13, args_ND) ≈ 5.64401178083053e-17 # non-dimensional ε = [0.0; 0.0] compute_τII!(ε, x1_ND, [1e0; 2.0], args_ND1) @test ε ≈ [3.65254588484434e-25, 7.305079089024537e-25] # vector input - - x1_D = ViscosityPartialMelt_Costa_etal_2009() - args_D=(;T=1000K, ϕ=0.5) - ε_D = compute_εII(x1_D, 1e6Pa, args_D) - @test ustrip(ε_D) ≈ 3.4537433628446676e-6 - - η = compute_viscosity_τII(x1_D, 1e6Pa, args_D) - @test ustrip(η) ≈ 1.447704555523709e11 + x1_D = ViscosityPartialMelt_Costa_etal_2009() + args_D = (;T=1000K, ϕ=0.5) + ε_D = compute_εII(x1_D, 1e6Pa, args_D) + @test ustrip(ε_D) ≈ 3.4537433628446676e-6 + η = compute_viscosity_τII(x1_D, 1e6Pa, args_D) + @test ustrip(η) ≈ 1.447704555523709e11 @test dεII_dτII(x1_ND, 1e6, args_ND) ≈ 31884.63277076202 @test dτII_dεII(x1_ND, 1e-15, args_ND) ≈ 0.0005644011780830529 - # ----- # ---- # Create plot - T = Vector(800.0:1400)*K - η_basalt = zeros(size(T))*Pas - η_basalt1 = zeros(size(T))*Pas - η_rhyolite = zeros(size(T))*Pas - ϕ_basalt = zeros(size(T)) - ϕ_rhyolite = zeros(size(T)) - - x_basalt = ViscosityPartialMelt_Costa_etal_2009(η = LinearMeltViscosity()) - x_rhyolite = ViscosityPartialMelt_Costa_etal_2009(η= LinearMeltViscosity(A = -8.1590, B = 2.4050e+04K, T0 = -430.9606K)) # Rhyolite - - mf_basalt = MeltingParam_Smooth3rdOrder() + T = Vector(800.0:1400)*K + η_basalt = zeros(size(T))*Pas + η_basalt1 = zeros(size(T))*Pas + η_rhyolite = zeros(size(T))*Pas + ϕ_basalt = zeros(size(T)) + ϕ_rhyolite = zeros(size(T)) + x_basalt = ViscosityPartialMelt_Costa_etal_2009(η = LinearMeltViscosity()) + x_rhyolite = ViscosityPartialMelt_Costa_etal_2009(η= LinearMeltViscosity(A = -8.1590, B = 2.4050e+04K, T0 = -430.9606K)) # Rhyolite + mf_basalt = MeltingParam_Smooth3rdOrder() mf_rhyolite = MeltingParam_Smooth3rdOrder( a=3043.0, b=-10552.0, c=12204.9, d = -4709.0 ) εII = 1e-5/s; for i in eachindex(T) - args_D = (; T = T[i]) - ϕ_basalt[i] = compute_meltfraction(mf_basalt, (; T=T[i]/K)) - η_basalt[i] = compute_viscosity_εII(x_basalt, εII, (; ϕ=ϕ_basalt[i], T=T[i])) - - η_basalt1[i] = compute_viscosity_εII(x_basalt, 1e-3/s, (; ϕ=ϕ_basalt[i], T=T[i])) - + args_D = (; T = T[i]) + ϕ_basalt[i] = compute_meltfraction(mf_basalt, (; T=T[i]/K)) + η_basalt[i] = compute_viscosity_εII(x_basalt, εII, (; ϕ=ϕ_basalt[i], T=T[i])) + η_basalt1[i] = compute_viscosity_εII(x_basalt, 1e-3/s, (; ϕ=ϕ_basalt[i], T=T[i])) ϕ_rhyolite[i] = compute_meltfraction(mf_rhyolite, (; T=T[i]/K)) - η_rhyolite[i] = compute_viscosity_εII(x_rhyolite, εII, (; ϕ=ϕ_rhyolite[i], T=T[i])) + η_rhyolite[i] = compute_viscosity_εII(x_rhyolite, εII, (; ϕ=ϕ_rhyolite[i], T=T[i])) end - @test mean(η_rhyolite) ≈ 1.26143880075233e19Pas - @test mean(η_basalt) ≈ 5.822731206904198e24Pas - @test mean(η_basalt1) ≈ 8.638313729394237e21Pas + @test mean(η_rhyolite) ≈ 1.262261299272795e19Pas + @test mean(η_basalt) ≈ 5.822731206904198e24Pas + @test mean(η_basalt1) ≈ 8.638313729394237e21Pas #= using Plots