Skip to content

Commit

Permalink
fix test
Browse files Browse the repository at this point in the history
  • Loading branch information
@albert-de-montserrat
albert-de-montserrat committed Oct 19, 2024
1 parent aa97968 commit 6dfd983
Showing 1 changed file with 42 additions and 56 deletions.
98 changes: 42 additions & 56 deletions test/test_CreepLaw.jl
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,6 @@
using Test, Statistics
using GeoParams


@testset "CreepLaw" begin

#Make sure structs are isbits
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand All @@ -199,71 +193,63 @@ 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]
args_ND1 = (; T = 1000*ones(size(τ)), ϕ=0.5*ones(size(τ)))
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
Expand Down

0 comments on commit 6dfd983

Please sign in to comment.