diff --git a/src/CreepLaw/DislocationCreep.jl b/src/CreepLaw/DislocationCreep.jl index 9d28cccb5..15335b16c 100644 --- a/src/CreepLaw/DislocationCreep.jl +++ b/src/CreepLaw/DislocationCreep.jl @@ -70,6 +70,14 @@ function computeCreepLaw_EpsII(TauII, a::DislocationCreep, p::CreepLawVariables) @unpack_val n,r,A,E,V,R = a @unpack_val P,T,f = p + FT, FE = CorrectionFactor(a) + + return A*(TauII*FT)^n*f^r*exp(-(E + P*V)/(R*T))/FE +end + +function computeCreepLaw_EpsII(TauII, a::DislocationCreep, P::_R, T::_R, f::_R) where _R<:Real + @unpack_val n,r,A,E,V,R = a + FT, FE = CorrectionFactor(a); return A*(TauII*FT)^n*f^r*exp(-(E + P*V)/(R*T))/FE; @@ -81,11 +89,23 @@ function computeCreepLaw_TauII(EpsII, a::DislocationCreep, p::CreepLawVariables) @unpack_val n,r,A,E,V,R = a @unpack_val P,T,f = p - FT, FE = CorrectionFactor(a); + FT, FE = CorrectionFactor(a) return A^(-1/n)*(EpsII*FE)^(1/n)*f^(-r/n)*exp((E + P*V)/(n * R*T))/FT; end + +# EpsII .= A.*(TauII.*FT).^n.*f.^r.*exp.(-(E.+P.*V)./(R.*T))./FE; Once we have a +# All inputs must be non-dimensionalized (or converted to consistent units) GeoUnits +function computeCreepLaw_TauII(EpsII, a::DislocationCreep, P::_R, T::_R, f::_R) where _R<:Real + @unpack_val n,r,A,E,V,R = a + + FT, FE = CorrectionFactor(a); + + return A^(-1/n)*(EpsII*FE)^(1/n)*f^(-r/n)*exp((E + P*V)/(n * R*T))/FT +end + + # Print info function show(io::IO, g::DislocationCreep) print(io, "DislocationCreep: Name = $(String(collect(g.Name))), n=$(g.n.val), r=$(g.r.val), A=$(g.A.val), E=$(g.E.val), V=$(g.V.val), Apparatus=$(g.Apparatus)" ) diff --git a/test/test_DislocationCreep.jl b/test/test_DislocationCreep.jl index fb9add3eb..a0ec98e6d 100644 --- a/test/test_DislocationCreep.jl +++ b/test/test_DislocationCreep.jl @@ -3,51 +3,53 @@ using GeoParams @testset "DislocationCreepLaws" begin -# This tests the MaterialParameters structure -CharUnits_GEO = GEO_units(viscosity=1e19, length=1000km); - -# Define a linear viscous creep law --------------------------------- -x1 = DislocationCreep() -@test isbits(x1) -@test x1.n.val == 1.0 -@test x1.A.val == 1.5MPa^-1*s^-1 + # This tests the MaterialParameters structure + CharUnits_GEO = GEO_units(; viscosity=1e19, length=1000km) -x2 = DislocationCreep(n=3) -@test x2.A.val == 1.5MPa^-3*s^-1 + # Define a linear viscous creep law --------------------------------- + x1 = DislocationCreep() + @test isbits(x1) + @test x1.n.val == 1.0 + @test x1.A.val == 1.5 + x2 = DislocationCreep(; n=3) + @test x2.A.val == 1.5 - -# perform a computation with the dislocation creep laws + # perform a computation with the dislocation creep laws # Calculate EpsII, using a set of pre-defined values -CharDim = GEO_units() -EpsII = GeoUnit(0s^-1) - Nondimensionalize!(EpsII,CharDim) -TauII = GeoUnit(100MPa) - Nondimensionalize!(TauII,CharDim) -P = GeoUnit(100MPa) - Nondimensionalize!(P,CharDim) -T = GeoUnit(500C) - Nondimensionalize!(T,CharDim) -f = GeoUnit(50MPa) - Nondimensionalize!(f,CharDim) -p = CreepLawVariables(P=P,T=T,f=f) -Phase = SetMaterialParams(Name="Viscous Matrix", Phase=2, - Density = ConstantDensity(), - CreepLaws = DislocationCreep(n=3NoUnits, r=1NoUnits), CharDim = CharDim) -EpsII.val = computeCreepLaw_EpsII(TauII,Phase.CreepLaws[1],p) -@test EpsII.val ≈ 2.1263214994323903e-11 rtol = 1e-8 + CharDim = GEO_units() + EpsII = GeoUnit(0s^-1) + EpsII = nondimensionalize(EpsII, CharDim) + TauII = GeoUnit(100MPa) + TauII = nondimensionalize(TauII, CharDim) + P = GeoUnit(100MPa) + P = nondimensionalize(P, CharDim) + T = GeoUnit(500C) + T = nondimensionalize(T, CharDim) + f = GeoUnit(50MPa) + f = nondimensionalize(f, CharDim) + p = CreepLawVariables(; P=P, T=T, f=f) + Phase = SetMaterialParams(; + Name="Viscous Matrix", + Phase=2, + Density=ConstantDensity(), + CreepLaws=DislocationCreep(; n=3NoUnits, r=1NoUnits), + CharDim=CharDim, + ) + εII = computeCreepLaw_EpsII(TauII, Phase.CreepLaws[1], p) + @test εII ≈ 2.1263214994323903e-11 rtol = 1e-8 + + @test εII == computeCreepLaw_EpsII(TauII, Phase.CreepLaws[1], P.val, T.val, f.val) # Check that once inverted, we get back the TauII that we used to calculate EpsII -NewTau = computeCreepLaw_EpsII(EpsII,Phase.CreepLaws[1],p) -@test NewTau ≈ TauII.val - - -# Given stress -#@test computeCreepLaw_EpsII(1e6Pa, x1, CreepLawParams())==5e-13/s # dimensional input - -# Given strainrate -#@test computeCreepLaw_EpsII(1e-13/s, x1, CreepLawParams())==1e18*2*1e-13Pa # dimensional input -# ------------------------------------------------------------------- + NewTau = computeCreepLaw_TauII(εII, Phase.CreepLaws[1], p) + @test NewTau ≈ TauII.val + @test NewTau == computeCreepLaw_TauII(εII, Phase.CreepLaws[1], P.val, T.val, f.val) + # Given stress + #@test computeCreepLaw_EpsII(1e6Pa, x1, CreepLawParams())==5e-13/s # dimensional input -end \ No newline at end of file + # Given strainrate + #@test computeCreepLaw_EpsII(1e-13/s, x1, CreepLawParams())==1e18*2*1e-13Pa # dimensional input + # ------------------------------------------------------------------- +end