From 33023df8277687a587a0d3549a1b22de8bacf7c1 Mon Sep 17 00:00:00 2001 From: Marcel Thielmann Date: Thu, 12 Feb 2026 15:19:03 +0100 Subject: [PATCH 01/10] Added a variant of the Herschel Bulkley rheology based on an earlier version of Albert --- src/CreepLaw/CreepLaw.jl | 3 +- src/CreepLaw/HerschelBulkley.jl | 84 +++++++++++++++++++++++++++++++++ src/GeoParams.jl | 1 + src/Viscosity/Viscosity.jl | 19 ++++++++ 4 files changed, 106 insertions(+), 1 deletion(-) create mode 100644 src/CreepLaw/HerschelBulkley.jl diff --git a/src/CreepLaw/CreepLaw.jl b/src/CreepLaw/CreepLaw.jl index 644aada13..850d96857 100644 --- a/src/CreepLaw/CreepLaw.jl +++ b/src/CreepLaw/CreepLaw.jl @@ -16,7 +16,7 @@ abstract type AbstractCreepLaw{T} <: AbstractConstitutiveLaw{T} end export isvolumetric, - LinearViscous, PowerlawViscous, CorrectionFactor, AbstractCreepLaw, ArrheniusType + LinearViscous, PowerlawViscous, CorrectionFactor, AbstractCreepLaw, ArrheniusType, HerschelBulkley isvolumetric(a::AbstractCreepLaw) = false @@ -75,6 +75,7 @@ include("PeierlsCreep.jl") include("NonLinearPeierlsCreep.jl") include("CustomRheology.jl") include("MeltViscosity.jl") +include("HerschelBulkley.jl") # Linear viscous rheology ------------------------------------------------ """ diff --git a/src/CreepLaw/HerschelBulkley.jl b/src/CreepLaw/HerschelBulkley.jl new file mode 100644 index 000000000..f7cb7ce7b --- /dev/null +++ b/src/CreepLaw/HerschelBulkley.jl @@ -0,0 +1,84 @@ + + +struct HerschelBulkley{T, U1, U2, U3} <: AbstractCreepLaw{T} + n::T # shear thinning exponent + η0::GeoUnit{T, U1} # "rigid" viscosity + τ0::GeoUnit{T, U2} # critical stress + ηr::GeoUnit{T, U1} # reference viscosity at the critical strain rate, which is given by 0.5*τ0/η0 and the critical temperature + Q::GeoUnit{T, U3} # temperature dependence of ηr, activation energy divided by R, unit is K + Tr::GeoUnit{T,U3} # reference temperature + function HerschelBulkley(; + n = 3, + η0 = 1e24Pa*s, + τ0 = 100e6Pa, + ηr = 1e20Pa*s, + Q = 0K, + Tr = 1600K, + ) + # Convert to GeoUnits + η0U = convert(GeoUnit, η0) + τ0U = convert(GeoUnit, τ0) + ηrU = convert(GeoUnit, ηr) + QU = convert(GeoUnit,Q) + Tr = convert(GeoUnit,Tr) + # Extract struct types + T = typeof(η0U).types[1] + U1 = typeof(η0U).types[2] + U2 = typeof(σ0U).types[2] + U3 = typeof(Tr).types[2] + # Create struct + return new{T, U1, U2, U3}( + n, η0U,τ0U,ηrU,QU, Tr + ) + end + + function HerschelBulkley(n, η0,τ0,ηr,Q,Tr) + return HerschelBulkley(;n = n, η0 = η0,τ0 = τ0,ηr = ηr,Q = Q,Tr = Tr) + end +end + + +# function when eII is given +function HB_rheology(εII,η0,τ0,εr,ηr,n,ηl) + # compute viscosity using the Souza-Mendes approach + η = (1.0 - exp(-2.0*η0*εII/τ0)) * (0.5*τ0/εII + ηr*(εII./εr)^(1/n - 1)) + ηl + return η +end + + + +# @inline function HerschelBulkley( +# a::DislocationCreep, +# EpsII; +# T = one(precision(a)), +# args..., +# ) +# n, η0, σ0, Kv, A, B= if EpsII isa Quantity +# @unpack_units n, η0, σ0, Kv, A, B = a +# n, η0, σ0, Kv, A, B +# else +# @unpack_val n, η0, σ0, Kv, A, B = a +# n, η0, σ0, Kv, A, B +# end + +# KvT = Kv * A * exp(B * (T - 273.15)) +# τ = @pow (1 - exp(-η0 * EpsII / σ0)) * (σ0 + KvT * (εII^n)) + +# end + +# @inline function compute_viscosity( +# a::HerschelBulkley; T = 0.0, kwargs... +# ) +# n, η0, σ0, Kv, A, B= if EpsII isa Quantity +# @unpack_units n, η0, σ0, Kv, A, B = a +# n, η0, σ0, Kv, A, B +# else +# @unpack_val n, η0, σ0, Kv, A, B = a +# n, η0, σ0, Kv, A, B +# end + +# KvT = Kv * A * exp(B * (T - 273.15)) + +# η = (1 - exp(-η0 *εII / σ0)) * (σ0 + KvT * (εII^n)) +# return η +# end \ No newline at end of file diff --git a/src/GeoParams.jl b/src/GeoParams.jl index 598db8002..e46d8c0be 100644 --- a/src/GeoParams.jl +++ b/src/GeoParams.jl @@ -174,6 +174,7 @@ module GeoParams GiordanoMeltViscosity, PowerlawViscous, ArrheniusType, + HerschelBulkley, CustomRheology, DislocationCreep, SetDislocationCreep, diff --git a/src/Viscosity/Viscosity.jl b/src/Viscosity/Viscosity.jl index f7aa12209..acca3363a 100644 --- a/src/Viscosity/Viscosity.jl +++ b/src/Viscosity/Viscosity.jl @@ -18,6 +18,25 @@ end @inline compute_viscosity_εII(v::LinearViscous, εII, args) = v.η.val @inline compute_viscosity_εII(v::ConstantElasticity, εII, args) = v.G * args.dt +@inline compute_viscosity_εII(v::HerschelBulkley, εII, args) = compute_viscosity_εII(v, εII; args...) + +@inline function compute_viscosity_εII(v::HerschelBulkley, εII; T = 0.0, kwargs...) + η0, τ0, ηr, Q, Tr = if εII isa Quantity + @unpack_units η0, τ0, ηr, Q, Tr = v + η0, τ0, ηr, Q, Tr + else + @unpack_val η0, τ0, ηr, Q, Tr = v + η0, τ0, ηr, Q, Tr + end + (; n) = v + + ηT = ηr * exp(B * (Q * (1/T-1/Tr))) # temperature dependence + εr = 0.5 * τ0/η0 # strain rate at which the Bingham yield stress is reached, this is defined as the reference strain rate + η = (1.0 - exp(-2.0*η0*εII/τ0)) * (0.5*τ0/εII + ηT*(εII./εr)^(1/n - 1)) + return η +end + + # compute effective "creep" viscosity from deviatoric stress tensor """ From 3169be3ae7f2c399d502439adb68caabf16f1c13 Mon Sep 17 00:00:00 2001 From: Marcel Thielmann Date: Fri, 13 Feb 2026 09:58:36 +0100 Subject: [PATCH 02/10] removed a function where it did not belong, --- src/CreepLaw/HerschelBulkley.jl | 9 --------- 1 file changed, 9 deletions(-) diff --git a/src/CreepLaw/HerschelBulkley.jl b/src/CreepLaw/HerschelBulkley.jl index f7cb7ce7b..8de26bfd1 100644 --- a/src/CreepLaw/HerschelBulkley.jl +++ b/src/CreepLaw/HerschelBulkley.jl @@ -38,15 +38,6 @@ struct HerschelBulkley{T, U1, U2, U3} <: AbstractCreepLaw{T} end -# function when eII is given -function HB_rheology(εII,η0,τ0,εr,ηr,n,ηl) - # compute viscosity using the Souza-Mendes approach - η = (1.0 - exp(-2.0*η0*εII/τ0)) * (0.5*τ0/εII + ηr*(εII./εr)^(1/n - 1)) + ηl - return η -end - - - # @inline function HerschelBulkley( # a::DislocationCreep, # EpsII; From 2da15ac75b640e32e92b78df930e28f9c2e5debd Mon Sep 17 00:00:00 2001 From: Marcel Thielmann Date: Wed, 18 Feb 2026 16:42:02 +0100 Subject: [PATCH 03/10] Modified according to Albert's suggestions --- src/CreepLaw/HerschelBulkley.jl | 37 --------------------------------- src/Viscosity/Viscosity.jl | 2 +- 2 files changed, 1 insertion(+), 38 deletions(-) diff --git a/src/CreepLaw/HerschelBulkley.jl b/src/CreepLaw/HerschelBulkley.jl index 8de26bfd1..d47e80b97 100644 --- a/src/CreepLaw/HerschelBulkley.jl +++ b/src/CreepLaw/HerschelBulkley.jl @@ -36,40 +36,3 @@ struct HerschelBulkley{T, U1, U2, U3} <: AbstractCreepLaw{T} return HerschelBulkley(;n = n, η0 = η0,τ0 = τ0,ηr = ηr,Q = Q,Tr = Tr) end end - - -# @inline function HerschelBulkley( -# a::DislocationCreep, -# EpsII; -# T = one(precision(a)), -# args..., -# ) -# n, η0, σ0, Kv, A, B= if EpsII isa Quantity -# @unpack_units n, η0, σ0, Kv, A, B = a -# n, η0, σ0, Kv, A, B -# else -# @unpack_val n, η0, σ0, Kv, A, B = a -# n, η0, σ0, Kv, A, B -# end - -# KvT = Kv * A * exp(B * (T - 273.15)) -# τ = @pow (1 - exp(-η0 * EpsII / σ0)) * (σ0 + KvT * (εII^n)) - -# end - -# @inline function compute_viscosity( -# a::HerschelBulkley; T = 0.0, kwargs... -# ) -# n, η0, σ0, Kv, A, B= if EpsII isa Quantity -# @unpack_units n, η0, σ0, Kv, A, B = a -# n, η0, σ0, Kv, A, B -# else -# @unpack_val n, η0, σ0, Kv, A, B = a -# n, η0, σ0, Kv, A, B -# end - -# KvT = Kv * A * exp(B * (T - 273.15)) - -# η = (1 - exp(-η0 *εII / σ0)) * (σ0 + KvT * (εII^n)) -# return η -# end \ No newline at end of file diff --git a/src/Viscosity/Viscosity.jl b/src/Viscosity/Viscosity.jl index acca3363a..ec5926689 100644 --- a/src/Viscosity/Viscosity.jl +++ b/src/Viscosity/Viscosity.jl @@ -32,7 +32,7 @@ end ηT = ηr * exp(B * (Q * (1/T-1/Tr))) # temperature dependence εr = 0.5 * τ0/η0 # strain rate at which the Bingham yield stress is reached, this is defined as the reference strain rate - η = (1.0 - exp(-2.0*η0*εII/τ0)) * (0.5*τ0/εII + ηT*(εII./εr)^(1/n - 1)) + η = @fastpow (1.0 - exp(-2.0*η0*εII/τ0)) * (0.5*τ0/εII + ηT*(εII./εr)^(1/n - 1)) return η end From 724391d9570b7279b0490cef917b213edba921bb Mon Sep 17 00:00:00 2001 From: Marcel Thielmann Date: Tue, 17 Mar 2026 17:44:24 +0100 Subject: [PATCH 04/10] changed fastpow to pow in Viscosity.jl --- src/Viscosity/Viscosity.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Viscosity/Viscosity.jl b/src/Viscosity/Viscosity.jl index ec5926689..1b1ef0031 100644 --- a/src/Viscosity/Viscosity.jl +++ b/src/Viscosity/Viscosity.jl @@ -32,7 +32,7 @@ end ηT = ηr * exp(B * (Q * (1/T-1/Tr))) # temperature dependence εr = 0.5 * τ0/η0 # strain rate at which the Bingham yield stress is reached, this is defined as the reference strain rate - η = @fastpow (1.0 - exp(-2.0*η0*εII/τ0)) * (0.5*τ0/εII + ηT*(εII./εr)^(1/n - 1)) + η = @pow (1.0 - exp(-2.0*η0*εII/τ0)) * (0.5*τ0/εII + ηT*(εII./εr)^(1/n - 1)) return η end From 1831933e2c1157bf938adfb6d24099e0d57fb76c Mon Sep 17 00:00:00 2001 From: Marcel Thielmann Date: Thu, 19 Mar 2026 11:08:17 +0100 Subject: [PATCH 05/10] added functions to compute viscosity, stress and strain rate --- src/CreepLaw/HerschelBulkley.jl | 41 +++++++++++++++++++++--- src/Viscosity/Viscosity.jl | 56 +++++++++++++++++++++++++++++++-- 2 files changed, 90 insertions(+), 7 deletions(-) diff --git a/src/CreepLaw/HerschelBulkley.jl b/src/CreepLaw/HerschelBulkley.jl index d47e80b97..1aec754ba 100644 --- a/src/CreepLaw/HerschelBulkley.jl +++ b/src/CreepLaw/HerschelBulkley.jl @@ -8,12 +8,12 @@ struct HerschelBulkley{T, U1, U2, U3} <: AbstractCreepLaw{T} Q::GeoUnit{T, U3} # temperature dependence of ηr, activation energy divided by R, unit is K Tr::GeoUnit{T,U3} # reference temperature function HerschelBulkley(; - n = 3, + n = 3.0, η0 = 1e24Pa*s, τ0 = 100e6Pa, ηr = 1e20Pa*s, - Q = 0K, - Tr = 1600K, + Q = 0.0K, + Tr = 1273K, ) # Convert to GeoUnits η0U = convert(GeoUnit, η0) @@ -24,7 +24,7 @@ struct HerschelBulkley{T, U1, U2, U3} <: AbstractCreepLaw{T} # Extract struct types T = typeof(η0U).types[1] U1 = typeof(η0U).types[2] - U2 = typeof(σ0U).types[2] + U2 = typeof(τ0U).types[2] U3 = typeof(Tr).types[2] # Create struct return new{T, U1, U2, U3}( @@ -36,3 +36,36 @@ struct HerschelBulkley{T, U1, U2, U3} <: AbstractCreepLaw{T} return HerschelBulkley(;n = n, η0 = η0,τ0 = τ0,ηr = ηr,Q = Q,Tr = Tr) end end + +function compute_εII(a::HershelBulkley, TauII; T = one(precision(a)), kwargs...) + η = compute_viscosity_τII(a, TauII; T = T) + EpsII = 0.5 * TauII / η + return EpsII +end + +function compute_εII(a::HershelBulkley, TauII; T = 1K, kwargs...) + η = compute_viscosity_τII(a, TauII; T = T) + EpsII = 0.5 * TauII / η + return EpsII +end + + +function compute_τII(a::HershelBulkley, EpsII; T = one(precision(a)), kwargs...) + η = compute_viscosity_εII(a, EpsII; T = T) + TauII = 2 * EpsII * η + return TauII +end + +function compute_τII(a::HershelBulkley, EpsII; T = 1K, kwargs...) + η = compute_viscosity_εII(a, EpsII; T = T) + TauII = 2 * EpsII * η + return TauII +end + +# print info +function show(io::IO, g::HershelBulkley) + return print( + io, + "Hershel Bulkley viscosity: η0=$(Value(g.η0)), τ0=$(Value(g.τ0)), ηr=$(Value(g.ηr)), n=$(Value(g.n)), Q=$(Value(g.Q)), Tr=$(Value(g.Tr))", + ) +end \ No newline at end of file diff --git a/src/Viscosity/Viscosity.jl b/src/Viscosity/Viscosity.jl index 1b1ef0031..4f879be37 100644 --- a/src/Viscosity/Viscosity.jl +++ b/src/Viscosity/Viscosity.jl @@ -20,7 +20,7 @@ end @inline compute_viscosity_εII(v::ConstantElasticity, εII, args) = v.G * args.dt @inline compute_viscosity_εII(v::HerschelBulkley, εII, args) = compute_viscosity_εII(v, εII; args...) -@inline function compute_viscosity_εII(v::HerschelBulkley, εII; T = 0.0, kwargs...) +@inline function compute_viscosity_εII(v::HerschelBulkley, εII; T = 1.0, kwargs...) η0, τ0, ηr, Q, Tr = if εII isa Quantity @unpack_units η0, τ0, ηr, Q, Tr = v η0, τ0, ηr, Q, Tr @@ -30,12 +30,62 @@ end end (; n) = v - ηT = ηr * exp(B * (Q * (1/T-1/Tr))) # temperature dependence + ηT = ηr * exp(Q * (1/T-1/Tr)) # temperature dependence εr = 0.5 * τ0/η0 # strain rate at which the Bingham yield stress is reached, this is defined as the reference strain rate - η = @pow (1.0 - exp(-2.0*η0*εII/τ0)) * (0.5*τ0/εII + ηT*(εII./εr)^(1/n - 1)) + η = @pow (1.0 - exp(-2.0*η0*εII/τ0)) * (0.5*τ0/εII + ηT*(εII/εr)^(1/n - 1)) return η end +@inline function compute_viscosity_τII(v::HerschelBulkley, τII; T = one(precision(a)), kwargs...) + + η0, τ0, ηr, Q, Tr = if τII isa Quantity + @unpack_units η0, τ0, ηr, Q, Tr = v + η0, τ0, ηr, Q, Tr + else + @unpack_val η0, τ0, ηr, Q, Tr = v + η0, τ0, ηr, Q, Tr + end + (; n) = v + + ηT = ηr * exp(Q * (1/T-1/Tr)) # temperature dependence + + # define residual function --> could we also call compute_viscosity_εII here? + function fres(εII,τII,η0,τ0,ηT,n) + η = @pow (1.0 - exp(-2.0*η0*εII/τ0)) * (0.5*τ0/εII + ηT*(εII/εr)^(1/n - 1)) + return 2.0 * η * εII - τII + end + + # define an initial guess for η + if τII < τ0 + η = η0 + εII = 0.5*τII/η0 + elseif τII > τ0 + εII = ((τII - τ0) / (2 * ηr * εr^(1 - one(n)/n)))^n + η = 0.5*τII/εII + end + + # Newton-Raphson method to get a solution for η + it_max = 100 + tol = 1e-8 + + # compute the residual of the initial guess + res = fres(εII,τII,η0,τ0,ηT,n) + + if abs(res) > tol + for _ in 1:it_max + f, dfdε = value_and_partial(εII -> fres(εII,τII,η0,τ0,ηT,n), εII) # value and dF/dεII via ForwardDiff + Δε = f / dfdε + εII -= Δε # adapt εII + abs(Δε) / (abs(εII) + eps(typeof(εII))) < tol && break # stop if we are below the tolerance + end + η = @pow (1.0 - exp(-2.0*η0*εII/τ0)) * (0.5*τ0/εII + ηT*(εII/εr)^(1/n - 1)) # --> could we also call compute_viscosity_εII here? + end + return η +end + + + + # compute effective "creep" viscosity from deviatoric stress tensor From a09cc51d906c611b1fae2741079d3fac6f78ee8b Mon Sep 17 00:00:00 2001 From: Marcel Thielmann Date: Thu, 19 Mar 2026 11:39:15 +0100 Subject: [PATCH 06/10] fixed compute_viscosity_tauII --- src/Viscosity/Viscosity.jl | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/Viscosity/Viscosity.jl b/src/Viscosity/Viscosity.jl index 4f879be37..f09ea6907 100644 --- a/src/Viscosity/Viscosity.jl +++ b/src/Viscosity/Viscosity.jl @@ -48,9 +48,10 @@ end (; n) = v ηT = ηr * exp(Q * (1/T-1/Tr)) # temperature dependence - + εr = 0.5 * τ0/η0 + # define residual function --> could we also call compute_viscosity_εII here? - function fres(εII,τII,η0,τ0,ηT,n) + function fres(εII,τII,η0,τ0,ηT,n,εr) η = @pow (1.0 - exp(-2.0*η0*εII/τ0)) * (0.5*τ0/εII + ηT*(εII/εr)^(1/n - 1)) return 2.0 * η * εII - τII end @@ -69,11 +70,11 @@ end tol = 1e-8 # compute the residual of the initial guess - res = fres(εII,τII,η0,τ0,ηT,n) + res = fres(εII,τII,η0,τ0,ηT,n,εr) if abs(res) > tol for _ in 1:it_max - f, dfdε = value_and_partial(εII -> fres(εII,τII,η0,τ0,ηT,n), εII) # value and dF/dεII via ForwardDiff + f, dfdε = value_and_partial(εII -> fres(εII,τII,η0,τ0,ηT,n,εr), εII) # value and dF/dεII via ForwardDiff Δε = f / dfdε εII -= Δε # adapt εII abs(Δε) / (abs(εII) + eps(typeof(εII))) < tol && break # stop if we are below the tolerance From ba80216be4e3a11f4e4ee3e3683f759d50fcf009 Mon Sep 17 00:00:00 2001 From: Pascal Aellig Date: Thu, 19 Mar 2026 12:01:40 +0100 Subject: [PATCH 07/10] progress, fix a few bugs, missing import? --- src/ConstitutiveRelationships.jl | 2 +- src/CreepLaw/HerschelBulkley.jl | 59 ++++++++++++++++++++++++++++---- test/test_CreepLaw.jl | 45 ++++++++++++++++++++++++ 3 files changed, 98 insertions(+), 8 deletions(-) diff --git a/src/ConstitutiveRelationships.jl b/src/ConstitutiveRelationships.jl index c984f0663..5c1985fa7 100644 --- a/src/ConstitutiveRelationships.jl +++ b/src/ConstitutiveRelationships.jl @@ -70,7 +70,7 @@ export param_info, for myType in ( :LinearViscous, :LinearMeltViscosity, :ViscosityPartialMelt_Costa_etal_2009, :GiordanoMeltViscosity, :DiffusionCreep, :DislocationCreep, :ConstantElasticity, :DruckerPrager, :ArrheniusType, - :GrainBoundarySliding, :PeierlsCreep, :NonLinearPeierlsCreep, :PowerlawViscous, + :GrainBoundarySliding, :PeierlsCreep, :NonLinearPeierlsCreep, :PowerlawViscous, :HerschelBulkley ) @eval begin @inline compute_εII(a::$(myType), TauII, args) = compute_εII(a, TauII; args...) diff --git a/src/CreepLaw/HerschelBulkley.jl b/src/CreepLaw/HerschelBulkley.jl index 1aec754ba..539d12840 100644 --- a/src/CreepLaw/HerschelBulkley.jl +++ b/src/CreepLaw/HerschelBulkley.jl @@ -1,4 +1,8 @@ - +export HerschelBulkley, + compute_εII!, + compute_εII, + compute_τII!, + compute_τII struct HerschelBulkley{T, U1, U2, U3} <: AbstractCreepLaw{T} n::T # shear thinning exponent @@ -37,35 +41,76 @@ struct HerschelBulkley{T, U1, U2, U3} <: AbstractCreepLaw{T} end end -function compute_εII(a::HershelBulkley, TauII; T = one(precision(a)), kwargs...) +function compute_εII(a::HerschelBulkley, TauII; T = one(precision(a)), kwargs...) η = compute_viscosity_τII(a, TauII; T = T) EpsII = 0.5 * TauII / η return EpsII end -function compute_εII(a::HershelBulkley, TauII; T = 1K, kwargs...) +function compute_εII(a::HerschelBulkley, TauII::Quantity; T = 1K, kwargs...) η = compute_viscosity_τII(a, TauII; T = T) EpsII = 0.5 * TauII / η return EpsII end +""" + compute_εII!(EpsII::AbstractArray{T, N}, a::HerschelBulkley, TauII::AbstractArray{T, N}; T = ones(size(TauII)), kwargs...) + +In-place function for the second invariant of the strain rate for Herschel-Bulkley rheology. +""" +function compute_εII!( + EpsII::AbstractArray{_T, N}, + a::HerschelBulkley, + TauII::AbstractArray{_T, N}; + T = ones(size(TauII))::AbstractArray{_T, N}, + kwargs..., +) where {_T, N} + @inbounds for i in eachindex(EpsII) + EpsII[i] = compute_εII(a, TauII[i]; T = T[i]) + end + + return nothing +end + +""" + compute_τII(a::HerschelBulkley, EpsII; T = one(precision(a)), kwargs...) -function compute_τII(a::HershelBulkley, EpsII; T = one(precision(a)), kwargs...) +""" +function compute_τII(a::HerschelBulkley, EpsII; T = one(precision(a)), kwargs...) η = compute_viscosity_εII(a, EpsII; T = T) TauII = 2 * EpsII * η return TauII end -function compute_τII(a::HershelBulkley, EpsII; T = 1K, kwargs...) +function compute_τII(a::HerschelBulkley, EpsII::Quantity; T = 1K, kwargs...) η = compute_viscosity_εII(a, EpsII; T = T) TauII = 2 * EpsII * η return TauII end +""" + compute_τII!(TauII::AbstractArray{T, N}, a::HerschelBulkley, EpsII::AbstractArray{T, N}; T = ones(size(EpsII)), kwargs...) + +In-place function for the second invariant of the stress for Herschel-Bulkley rheology. +""" +function compute_τII!( + TauII::AbstractArray{_T, N}, + a::HerschelBulkley, + EpsII::AbstractArray{_T, N}; + T = ones(size(EpsII)), + kwargs..., +) where {_T, N} + @inbounds for i in eachindex(TauII) + TauII[i] = compute_τII(a, EpsII[i]; T = T[i]) + end + + return nothing +end + # print info -function show(io::IO, g::HershelBulkley) +function show(io::IO, g::HerschelBulkley) return print( io, "Hershel Bulkley viscosity: η0=$(Value(g.η0)), τ0=$(Value(g.τ0)), ηr=$(Value(g.ηr)), n=$(Value(g.n)), Q=$(Value(g.Q)), Tr=$(Value(g.Tr))", ) -end \ No newline at end of file +end diff --git a/test/test_CreepLaw.jl b/test/test_CreepLaw.jl index 1ed916467..b74ee859f 100644 --- a/test/test_CreepLaw.jl +++ b/test/test_CreepLaw.jl @@ -22,6 +22,9 @@ using GeoParams x = GiordanoMeltViscosity() @test isbits(x) + + x = HerschelBulkley() + @test isbits(x) # This tests the MaterialParameters structure CharUnits_GEO = GEO_units(; viscosity = 1.0e19, length = 1000km) @@ -343,6 +346,48 @@ using GeoParams end T_plots = 10000 ./ ustrip.(T) + + # Herschel-Bulkley viscosity (Regularized Bingham rheology) + x1 = HerschelBulkley() + @test x1.n == 3 + @test x1.η0.val == 1e24 + @test x1.τ0.val == 100e6 + @test x1.ηr.val == 1e20 + @test x1.Q.val == 0 + @test x1.Tr.val == 1273 + + x1_D = HerschelBulkley() + @test isDimensional(x1_D) == true + x1_ND = nondimensionalize(x1_D, CharUnits_GEO) + @test isDimensional(x1_ND) == false + @test x1_ND.η0 * 1.0 == 1e5 + @test x1_ND.τ0 * 1.0 ≈ 10.0 + @test x1_ND.ηr * 1.0 == 1e1 + @test x1_ND.Q * 1.0 == 0.0 + @test x1_ND.Tr * 1.0 ≈ 1.0 rtol = 1.0e-3 + @test dimensionalize(x1_ND, CharUnits_GEO) == x1_D + + args = (; T = 1273K,) + Tau_II = 1.0e6Pa + compute_εII(x1_D, Tau_II, args) + + + # and strain rate + ε1 = 0.5 + tau = compute_τII(x1, ε1, args) + @test tau ≈ 1.0 + # using a vector input ------------------ + T2 = [1.0; 0.9; 0.8] + # and stress + ε21 = [0.0; 0.0; 0.0] + τ21 = [1.0; 0.8; 0.9] + args = (T = T2,) + compute_εII!(ε21, x1, τ21, args) + # and strain rate + τ22 = [0.0; 0.0; 0.0] + ε22 = [0.5; 1.0; 0.2] + compute_τII!(τ22, x1, ε22, args) + # fig = Figure(size = (800, 600)) # ax = Axis(fig[1, 1], xlabel = "T [C]", ylabel = "log10(η [Pas])") # lines!(ax, ustrip.(T) .- 273.15, log10.(ustrip.(η_MORB)), label = "MORB, εII=1e-5/s") From 1af64268bc93b69899152386438613fb71d3e813 Mon Sep 17 00:00:00 2001 From: Marcel Thielmann Date: Mon, 23 Mar 2026 13:18:07 +0100 Subject: [PATCH 08/10] modified HerschelBulkley.jl and test_CreepLaw.jl --- src/CreepLaw/HerschelBulkley.jl | 105 +++++++++++++++++++-- src/Viscosity/Viscosity.jl | 6 +- test/test_CreepLaw.jl | 157 ++++++++++++++++++++++++++------ 3 files changed, 231 insertions(+), 37 deletions(-) diff --git a/src/CreepLaw/HerschelBulkley.jl b/src/CreepLaw/HerschelBulkley.jl index 539d12840..8d6cf3e60 100644 --- a/src/CreepLaw/HerschelBulkley.jl +++ b/src/CreepLaw/HerschelBulkley.jl @@ -2,7 +2,9 @@ export HerschelBulkley, compute_εII!, compute_εII, compute_τII!, - compute_τII + compute_τII, + compute_hb_viscosity_εII, + compute_hb_viscosity_τII struct HerschelBulkley{T, U1, U2, U3} <: AbstractCreepLaw{T} n::T # shear thinning exponent @@ -42,13 +44,13 @@ struct HerschelBulkley{T, U1, U2, U3} <: AbstractCreepLaw{T} end function compute_εII(a::HerschelBulkley, TauII; T = one(precision(a)), kwargs...) - η = compute_viscosity_τII(a, TauII; T = T) + η = compute_hb_viscosity_τII(a, TauII; T = T) EpsII = 0.5 * TauII / η return EpsII end function compute_εII(a::HerschelBulkley, TauII::Quantity; T = 1K, kwargs...) - η = compute_viscosity_τII(a, TauII; T = T) + η = compute_hb_viscosity_τII(a, TauII; T = T) EpsII = 0.5 * TauII / η return EpsII end @@ -77,13 +79,13 @@ end """ function compute_τII(a::HerschelBulkley, EpsII; T = one(precision(a)), kwargs...) - η = compute_viscosity_εII(a, EpsII; T = T) + η = compute_hb_viscosity_εII(a, EpsII; T = T) TauII = 2 * EpsII * η return TauII end function compute_τII(a::HerschelBulkley, EpsII::Quantity; T = 1K, kwargs...) - η = compute_viscosity_εII(a, EpsII; T = T) + η = compute_hb_viscosity_εII(a, EpsII; T = T) TauII = 2 * EpsII * η return TauII end @@ -107,10 +109,101 @@ function compute_τII!( return nothing end + +""" + compute_hb_viscosity_εII(a::HerschelBulkley, EpsII; T = one(precision(a)), kwargs...) + +function to compute the viscosity if EpsII is given +""" +@inline function compute_hb_viscosity_εII(v::HerschelBulkley, εII; T = 1.0, kwargs...) + η0, τ0, ηr, Q, Tr = if εII isa Quantity + @unpack_units η0, τ0, ηr, Q, Tr = v + η0, τ0, ηr, Q, Tr + else + @unpack_val η0, τ0, ηr, Q, Tr = v + η0, τ0, ηr, Q, Tr + end + (; n) = v + + ηT = ηr * exp(Q * (1/T-1/Tr)) # temperature dependence + εr = 0.5 * τ0/η0 # strain rate at which the Bingham yield stress is reached, this is defined as the reference strain rate + η = @pow (1.0 - exp(-2.0*η0*εII/τ0)) * (0.5*τ0/εII + ηT*(εII/εr)^(one(n)/n - 1)) + return η +end + + +""" +compute_hb_viscosity_τII(a::HerschelBulkley, EpsII; T = one(precision(a)), kwargs...) + +function to compute the viscosity if TauII is given +""" + +@inline function compute_hb_viscosity_τII(v::HerschelBulkley, τII; T = one(precision(v)), kwargs...) + + η0, τ0, ηr, Q, Tr = if τII isa Quantity + @unpack_units η0, τ0, ηr, Q, Tr = v + η0, τ0, ηr, Q, Tr + else + @unpack_val η0, τ0, ηr, Q, Tr = v + T = ustrip(T) + η0, τ0, ηr, Q, Tr + end + (; n) = v + + ηT = ηr * exp(Q * (1/T - 1/Tr)) + εr = 0.5 * τ0 / η0 + + # initial guess + η = if τII < τ0 + η0 + elseif τII == τ0 + (1 - exp(-one(η0))) * (η0 + ηT) + else + @pow (ηT / (1 - τ0/τII))^n * (τII / (2*εr))^(1-n) + end + + εII = 0.5 * τII / η + εII_unit = τII isa Quantity ? unit(εII) : one(εII) + + # strip ALL quantities to plain floats before Newton iteration + # so that ForwardDiff never sees Quantity{Dual} types + τII_s = ustrip(τII) + η0_s = ustrip(η0) + τ0_s = ustrip(τ0) + ηT_s = ustrip(ηT) + εr_s = ustrip(εr) + εII_s = ustrip(εII) + + tol = 1e-10 + it_max = 100 + + for _ in 1:it_max + f, dfdε = value_and_partial( + ε -> fres_hb(ε, τII_s, η0_s, τ0_s, ηT_s, n, εr_s), # all plain floats + εII_s + ) + Δε = f / dfdε + εII_s -= Δε + abs(Δε) / (abs(εII_s) + eps(typeof(εII_s))) < tol && break + end + + εII = εII_s * εII_unit + η = @pow (1.0 - exp(-2.0*η0*εII/τ0)) * (0.5*τ0/εII + ηT*(εII/εr)^(one(n)/n - 1)) + + return η +end + +# non-dimensional residual +function fres_hb(εII::Number, τII::Number, η0::Number, τ0::Number, ηT::Number, n::Number, εr::Number) + η = @pow (1.0 - exp(-2.0*η0*εII/τ0)) * (0.5*τ0/εII + ηT*(εII/εr)^(one(n)/n - 1)) + return 2.0 * η * εII - τII +end + + # print info function show(io::IO, g::HerschelBulkley) return print( io, - "Hershel Bulkley viscosity: η0=$(Value(g.η0)), τ0=$(Value(g.τ0)), ηr=$(Value(g.ηr)), n=$(Value(g.n)), Q=$(Value(g.Q)), Tr=$(Value(g.Tr))", + "Hershel Bulkley viscosity: η0=$(Value(g.η0)), τ0=$(Value(g.τ0)), ηr=$(Value(g.ηr)), n=$(g.n), Q=$(Value(g.Q)), Tr=$(Value(g.Tr))", ) end diff --git a/src/Viscosity/Viscosity.jl b/src/Viscosity/Viscosity.jl index f09ea6907..6110308b8 100644 --- a/src/Viscosity/Viscosity.jl +++ b/src/Viscosity/Viscosity.jl @@ -18,8 +18,9 @@ end @inline compute_viscosity_εII(v::LinearViscous, εII, args) = v.η.val @inline compute_viscosity_εII(v::ConstantElasticity, εII, args) = v.G * args.dt -@inline compute_viscosity_εII(v::HerschelBulkley, εII, args) = compute_viscosity_εII(v, εII; args...) +@inline compute_viscosity_εII(v::HerschelBulkley, εII, args) = compute_hb_viscosity_εII(v, εII; args...) +#= @inline function compute_viscosity_εII(v::HerschelBulkley, εII; T = 1.0, kwargs...) η0, τ0, ηr, Q, Tr = if εII isa Quantity @unpack_units η0, τ0, ηr, Q, Tr = v @@ -84,7 +85,7 @@ end return η end - +=# @@ -103,6 +104,7 @@ end @inline compute_viscosity_τII(v::LinearViscous, τII, args) = v.η.val @inline compute_viscosity_τII(v::ConstantElasticity, τII, args) = v.G * args.dt +@inline compute_viscosity_τII(v::HerschelBulkley, εII, args) = compute_hb_viscosity_τII(v, εII; args...) for fn in (:compute_viscosity_εII, :compute_viscosity_τII) @eval begin diff --git a/test/test_CreepLaw.jl b/test/test_CreepLaw.jl index b74ee859f..6b9a8aaf8 100644 --- a/test/test_CreepLaw.jl +++ b/test/test_CreepLaw.jl @@ -348,16 +348,21 @@ using GeoParams # Herschel-Bulkley viscosity (Regularized Bingham rheology) - x1 = HerschelBulkley() - @test x1.n == 3 - @test x1.η0.val == 1e24 - @test x1.τ0.val == 100e6 - @test x1.ηr.val == 1e20 - @test x1.Q.val == 0 - @test x1.Tr.val == 1273 - - x1_D = HerschelBulkley() + x1_D = HerschelBulkley() @test isDimensional(x1_D) == true + @test x1_D.n == 3 + @test x1_D.η0.val == 1e24 + @test x1_D.τ0.val == 100e6 + @test x1_D.ηr.val == 1e20 + @test x1_D.Q.val == 0 + @test x1_D.Tr.val == 1273 + + @test x1_D.η0.unit == Pas + @test x1_D.τ0.unit == Pa + @test x1_D.ηr.unit == Pas + @test x1_D.Q.unit == K + @test x1_D.Tr.unit == K + x1_ND = nondimensionalize(x1_D, CharUnits_GEO) @test isDimensional(x1_ND) == false @test x1_ND.η0 * 1.0 == 1e5 @@ -367,26 +372,120 @@ using GeoParams @test x1_ND.Tr * 1.0 ≈ 1.0 rtol = 1.0e-3 @test dimensionalize(x1_ND, CharUnits_GEO) == x1_D - args = (; T = 1273K,) - Tau_II = 1.0e6Pa - compute_εII(x1_D, Tau_II, args) - - - # and strain rate - ε1 = 0.5 - tau = compute_τII(x1, ε1, args) - @test tau ≈ 1.0 - # using a vector input ------------------ - T2 = [1.0; 0.9; 0.8] - # and stress - ε21 = [0.0; 0.0; 0.0] - τ21 = [1.0; 0.8; 0.9] - args = (T = T2,) - compute_εII!(ε21, x1, τ21, args) - # and strain rate - τ22 = [0.0; 0.0; 0.0] - ε22 = [0.5; 1.0; 0.2] - compute_τII!(τ22, x1, ε22, args) + # below yield case - dimensional + εt = 0.5e-20 * s^-1 + τ1273 = compute_τII(x1_D, εt, (; T = 1273.0K)) + + @test ustrip(τ1273) ≈ 2*x1_D.η0.val*ustrip(εt) rtol=1e-4 + ε1273 = compute_εII(x1_D, τ1273, (; T = 1273.0K)) + @test εt ≈ ε1273 rtol = 1e-6 # see if we can recover the strain rate + + # below yield case - nondimensional + εr = 0.5 * x1_ND.τ0/x1_ND.η0 # yield strain rate + εt = 1e-4 * εr # well below the yield + + τ1 = compute_τII(x1_ND, εt, (; T = 1.0)) + @test τ1 ≈ 2*x1_ND.η0.val*ustrip(εt) rtol=1e-4 + ε1 = compute_εII(x1_ND, τ1, (; T = 1.0)) + @test εt ≈ ε1 rtol = 1e-6 # see if we can revocer the strain rate + + εt = [1e-4* εr; 1e-4 * εr] + τ1 = [1e3, 1e3] # initialize τ1 + compute_τII!(τ1,x1_ND, εt, (; T = ones(size(εt)))) + @test τ1 ≈ [2*x1_ND.η0.val*ustrip(εt[1]),2*x1_ND.η0.val*ustrip(εt[2])] rtol=1e-4 + + # above yield case, dimensional + x1_D = HerschelBulkley( + n = 3.0, + η0 = 1e24Pa*s, + τ0 = 100e6Pa, + ηr = 1e20Pa*s, + Q = 0.0K, + Tr = 1273K, + ) + + εt = 0.5e-9 * s^-1 + εr = 0.5 * x1_D.τ0.val/x1_D.η0.val # yield strain rate + + η_approx = 0.5*x1_D.τ0.val/ustrip(εt) + x1_D.ηr.val*(ustrip(εt)/εr)^(one(x1_D.n)/x1_D.n - 1) + τ1273 = compute_τII(x1_D, εt ; T = 1273.0K) + ε1273 = compute_εII(x1_D,τ1273; T = 1273.0K) + + @test ustrip(τ1273) ≈ 2*η_approx*ustrip(εt) rtol=1e-6 + @test ε1273 ≈ εt rtol = 1e-6 +# +# + x1_ND = nondimensionalize(x1_D, CharUnits_GEO) + εt_ND = nondimensionalize(εt, CharUnits_GEO) +# + τ1_ND = compute_τII(x1_ND, εt_ND ; T = 1.0) + ε1_ND = compute_εII(x1_ND,τ1_ND; T = 1.0) + + @test ε1_ND ≈ εt_ND rtol = 1e-6 + + + # temperature dependence test + x1_Q = HerschelBulkley( + n = 3.0, + η0 = 1e24Pa*s, + τ0 = 100e6Pa, + ηr = 1e20Pa*s, + Q = 63772.0K, # E/R for olivine, E ≈ 530 kJ/mol + Tr = 1273K, + ) + + # compute ηT at both temperatures analytically + Q_val = 63772.0 + Tr_val = 1273.0 + T1_val = 1273.0 # reference temperature → ηT = ηr + T2_val = 1473.0 # higher temperature → ηT < ηr + + ηT1 = x1_Q.ηr.val * exp(Q_val * (1/T1_val - 1/Tr_val)) # = 1e20 + ηT2 = x1_Q.ηr.val * exp(Q_val * (1/T2_val - 1/Tr_val)) # << 1e20 + + εr_Q = 0.5 * x1_Q.τ0.val / x1_Q.η0.val # reference strain rate + + # --- above yield, dimensional --- + εt_Q = 0.5e-9 * s^-1 # well above yield + + # at T = 1273K (reference): ηT = ηr, same as Q=0 case + τ_T1 = compute_τII(x1_Q, εt_Q; T = T1_val*K) + ε_T1 = compute_εII(x1_Q, τ_T1; T = T1_val*K) + η_approx_T1 = 0.5*x1_Q.τ0.val/ustrip(εt_Q) + ηT1*(ustrip(εt_Q)/εr_Q)^(one(x1_Q.n)/x1_Q.n - 1) + @test ustrip(τ_T1) ≈ 2*η_approx_T1*ustrip(εt_Q) rtol = 1e-6 + @test ε_T1 ≈ εt_Q rtol = 1e-6 + + # at T = 1473K: ηT < ηr → lower viscosity → lower stress at same strain rate + τ_T2 = compute_τII(x1_Q, εt_Q; T = T2_val*K) + ε_T2 = compute_εII(x1_Q, τ_T2; T = T2_val*K) + η_approx_T2 = 0.5*x1_Q.τ0.val/ustrip(εt_Q) + ηT2*(ustrip(εt_Q)/εr_Q)^(one(x1_Q.n)/x1_Q.n - 1) + @test ustrip(τ_T2) ≈ 2*η_approx_T2*ustrip(εt_Q) rtol = 1e-6 + @test ε_T2 ≈ εt_Q rtol = 1e-6 + + # monotonicity: higher T → lower τII at same εII (thermal softening) + @test ustrip(τ_T2) < ustrip(τ_T1) + + # monotonicity: higher T → higher εII at same τII (thermal softening) + τ_fixed = compute_τII(x1_Q, εt_Q; T = T1_val*K) + ε_at_T1 = compute_εII(x1_Q, τ_fixed; T = T1_val*K) + ε_at_T2 = compute_εII(x1_Q, τ_fixed; T = T2_val*K) + @test ustrip(ε_at_T2) > ustrip(ε_at_T1) + + # nondimensional temperature dependence + x1_Q_ND = nondimensionalize(x1_Q, CharUnits_GEO) + εt_ND = nondimensionalize(εt_Q, CharUnits_GEO) + T1_ND = nondimensionalize(T1_val*K, CharUnits_GEO) + T2_ND = nondimensionalize(T2_val*K, CharUnits_GEO) + + τ_T1_ND = compute_τII(x1_Q_ND, εt_ND; T = T1_ND) + τ_T2_ND = compute_τII(x1_Q_ND, εt_ND; T = T2_ND) + + # ND results should be consistent with dimensional results after rescaling + @test ustrip(nondimensionalize(τ_T1, CharUnits_GEO)) ≈ τ_T1_ND rtol = 1e-6 + @test ustrip(nondimensionalize(τ_T2, CharUnits_GEO)) ≈ τ_T2_ND rtol = 1e-6 + + # thermal softening should hold in ND as well + @test τ_T2_ND < τ_T1_ND # fig = Figure(size = (800, 600)) # ax = Axis(fig[1, 1], xlabel = "T [C]", ylabel = "log10(η [Pas])") From 738f7b4f9e00b211ef0747dccb65f9b055887dde Mon Sep 17 00:00:00 2001 From: Pascal Aellig Date: Mon, 23 Mar 2026 13:36:14 +0100 Subject: [PATCH 09/10] format --- src/ConstitutiveRelationships.jl | 2 +- src/CreepLaw/HerschelBulkley.jl | 82 ++++++++++---------- src/Viscosity/Viscosity.jl | 2 - test/test_CreepLaw.jl | 126 +++++++++++++++---------------- 4 files changed, 105 insertions(+), 107 deletions(-) diff --git a/src/ConstitutiveRelationships.jl b/src/ConstitutiveRelationships.jl index 5c1985fa7..bd849e359 100644 --- a/src/ConstitutiveRelationships.jl +++ b/src/ConstitutiveRelationships.jl @@ -70,7 +70,7 @@ export param_info, for myType in ( :LinearViscous, :LinearMeltViscosity, :ViscosityPartialMelt_Costa_etal_2009, :GiordanoMeltViscosity, :DiffusionCreep, :DislocationCreep, :ConstantElasticity, :DruckerPrager, :ArrheniusType, - :GrainBoundarySliding, :PeierlsCreep, :NonLinearPeierlsCreep, :PowerlawViscous, :HerschelBulkley + :GrainBoundarySliding, :PeierlsCreep, :NonLinearPeierlsCreep, :PowerlawViscous, :HerschelBulkley, ) @eval begin @inline compute_εII(a::$(myType), TauII, args) = compute_εII(a, TauII; args...) diff --git a/src/CreepLaw/HerschelBulkley.jl b/src/CreepLaw/HerschelBulkley.jl index 8d6cf3e60..501604611 100644 --- a/src/CreepLaw/HerschelBulkley.jl +++ b/src/CreepLaw/HerschelBulkley.jl @@ -12,21 +12,21 @@ struct HerschelBulkley{T, U1, U2, U3} <: AbstractCreepLaw{T} τ0::GeoUnit{T, U2} # critical stress ηr::GeoUnit{T, U1} # reference viscosity at the critical strain rate, which is given by 0.5*τ0/η0 and the critical temperature Q::GeoUnit{T, U3} # temperature dependence of ηr, activation energy divided by R, unit is K - Tr::GeoUnit{T,U3} # reference temperature + Tr::GeoUnit{T, U3} # reference temperature function HerschelBulkley(; - n = 3.0, - η0 = 1e24Pa*s, - τ0 = 100e6Pa, - ηr = 1e20Pa*s, - Q = 0.0K, - Tr = 1273K, + n = 3.0, + η0 = 1.0e24Pa * s, + τ0 = 100.0e6Pa, + ηr = 1.0e20Pa * s, + Q = 0.0K, + Tr = 1273K, ) # Convert to GeoUnits η0U = convert(GeoUnit, η0) τ0U = convert(GeoUnit, τ0) ηrU = convert(GeoUnit, ηr) - QU = convert(GeoUnit,Q) - Tr = convert(GeoUnit,Tr) + QU = convert(GeoUnit, Q) + Tr = convert(GeoUnit, Tr) # Extract struct types T = typeof(η0U).types[1] U1 = typeof(η0U).types[2] @@ -34,12 +34,12 @@ struct HerschelBulkley{T, U1, U2, U3} <: AbstractCreepLaw{T} U3 = typeof(Tr).types[2] # Create struct return new{T, U1, U2, U3}( - n, η0U,τ0U,ηrU,QU, Tr + n, η0U, τ0U, ηrU, QU, Tr ) end - function HerschelBulkley(n, η0,τ0,ηr,Q,Tr) - return HerschelBulkley(;n = n, η0 = η0,τ0 = τ0,ηr = ηr,Q = Q,Tr = Tr) + function HerschelBulkley(n, η0, τ0, ηr, Q, Tr) + return HerschelBulkley(; n = n, η0 = η0, τ0 = τ0, ηr = ηr, Q = Q, Tr = Tr) end end @@ -61,12 +61,12 @@ end In-place function for the second invariant of the strain rate for Herschel-Bulkley rheology. """ function compute_εII!( - EpsII::AbstractArray{_T, N}, - a::HerschelBulkley, - TauII::AbstractArray{_T, N}; - T = ones(size(TauII))::AbstractArray{_T, N}, - kwargs..., -) where {_T, N} + EpsII::AbstractArray{_T, N}, + a::HerschelBulkley, + TauII::AbstractArray{_T, N}; + T = ones(size(TauII))::AbstractArray{_T, N}, + kwargs..., + ) where {_T, N} @inbounds for i in eachindex(EpsII) EpsII[i] = compute_εII(a, TauII[i]; T = T[i]) end @@ -96,12 +96,12 @@ end In-place function for the second invariant of the stress for Herschel-Bulkley rheology. """ function compute_τII!( - TauII::AbstractArray{_T, N}, - a::HerschelBulkley, - EpsII::AbstractArray{_T, N}; - T = ones(size(EpsII)), - kwargs..., -) where {_T, N} + TauII::AbstractArray{_T, N}, + a::HerschelBulkley, + EpsII::AbstractArray{_T, N}; + T = ones(size(EpsII)), + kwargs..., + ) where {_T, N} @inbounds for i in eachindex(TauII) TauII[i] = compute_τII(a, EpsII[i]; T = T[i]) end @@ -118,17 +118,17 @@ function to compute the viscosity if EpsII is given @inline function compute_hb_viscosity_εII(v::HerschelBulkley, εII; T = 1.0, kwargs...) η0, τ0, ηr, Q, Tr = if εII isa Quantity @unpack_units η0, τ0, ηr, Q, Tr = v - η0, τ0, ηr, Q, Tr + η0, τ0, ηr, Q, Tr else @unpack_val η0, τ0, ηr, Q, Tr = v η0, τ0, ηr, Q, Tr end (; n) = v - ηT = ηr * exp(Q * (1/T-1/Tr)) # temperature dependence - εr = 0.5 * τ0/η0 # strain rate at which the Bingham yield stress is reached, this is defined as the reference strain rate - η = @pow (1.0 - exp(-2.0*η0*εII/τ0)) * (0.5*τ0/εII + ηT*(εII/εr)^(one(n)/n - 1)) - return η + ηT = ηr * exp(Q * (1 / T - 1 / Tr)) # temperature dependence + εr = 0.5 * τ0 / η0 # strain rate at which the Bingham yield stress is reached, this is defined as the reference strain rate + η = @pow (1.0 - exp(-2.0 * η0 * εII / τ0)) * (0.5 * τ0 / εII + ηT * (εII / εr)^(one(n) / n - 1)) + return η end @@ -150,8 +150,8 @@ function to compute the viscosity if TauII is given end (; n) = v - ηT = ηr * exp(Q * (1/T - 1/Tr)) - εr = 0.5 * τ0 / η0 + ηT = ηr * exp(Q * (1 / T - 1 / Tr)) + εr = 0.5 * τ0 / η0 # initial guess η = if τII < τ0 @@ -159,22 +159,22 @@ function to compute the viscosity if TauII is given elseif τII == τ0 (1 - exp(-one(η0))) * (η0 + ηT) else - @pow (ηT / (1 - τ0/τII))^n * (τII / (2*εr))^(1-n) + @pow (ηT / (1 - τ0 / τII))^n * (τII / (2 * εr))^(1 - n) end - εII = 0.5 * τII / η + εII = 0.5 * τII / η εII_unit = τII isa Quantity ? unit(εII) : one(εII) # strip ALL quantities to plain floats before Newton iteration # so that ForwardDiff never sees Quantity{Dual} types τII_s = ustrip(τII) - η0_s = ustrip(η0) - τ0_s = ustrip(τ0) - ηT_s = ustrip(ηT) - εr_s = ustrip(εr) + η0_s = ustrip(η0) + τ0_s = ustrip(τ0) + ηT_s = ustrip(ηT) + εr_s = ustrip(εr) εII_s = ustrip(εII) - tol = 1e-10 + tol = 1.0e-10 it_max = 100 for _ in 1:it_max @@ -182,20 +182,20 @@ function to compute the viscosity if TauII is given ε -> fres_hb(ε, τII_s, η0_s, τ0_s, ηT_s, n, εr_s), # all plain floats εII_s ) - Δε = f / dfdε + Δε = f / dfdε εII_s -= Δε abs(Δε) / (abs(εII_s) + eps(typeof(εII_s))) < tol && break end εII = εII_s * εII_unit - η = @pow (1.0 - exp(-2.0*η0*εII/τ0)) * (0.5*τ0/εII + ηT*(εII/εr)^(one(n)/n - 1)) + η = @pow (1.0 - exp(-2.0 * η0 * εII / τ0)) * (0.5 * τ0 / εII + ηT * (εII / εr)^(one(n) / n - 1)) return η end # non-dimensional residual function fres_hb(εII::Number, τII::Number, η0::Number, τ0::Number, ηT::Number, n::Number, εr::Number) - η = @pow (1.0 - exp(-2.0*η0*εII/τ0)) * (0.5*τ0/εII + ηT*(εII/εr)^(one(n)/n - 1)) + η = @pow (1.0 - exp(-2.0 * η0 * εII / τ0)) * (0.5 * τ0 / εII + ηT * (εII / εr)^(one(n) / n - 1)) return 2.0 * η * εII - τII end diff --git a/src/Viscosity/Viscosity.jl b/src/Viscosity/Viscosity.jl index 6110308b8..ceb993a62 100644 --- a/src/Viscosity/Viscosity.jl +++ b/src/Viscosity/Viscosity.jl @@ -88,8 +88,6 @@ end =# - - # compute effective "creep" viscosity from deviatoric stress tensor """ compute_viscosity_τII(s::AbstractConstitutiveLaw, τII, kwargs...) diff --git a/test/test_CreepLaw.jl b/test/test_CreepLaw.jl index 6b9a8aaf8..d0f64831d 100644 --- a/test/test_CreepLaw.jl +++ b/test/test_CreepLaw.jl @@ -348,12 +348,12 @@ using GeoParams # Herschel-Bulkley viscosity (Regularized Bingham rheology) - x1_D = HerschelBulkley() + x1_D = HerschelBulkley() @test isDimensional(x1_D) == true @test x1_D.n == 3 - @test x1_D.η0.val == 1e24 - @test x1_D.τ0.val == 100e6 - @test x1_D.ηr.val == 1e20 + @test x1_D.η0.val == 1.0e24 + @test x1_D.τ0.val == 100.0e6 + @test x1_D.ηr.val == 1.0e20 @test x1_D.Q.val == 0 @test x1_D.Tr.val == 1273 @@ -365,9 +365,9 @@ using GeoParams x1_ND = nondimensionalize(x1_D, CharUnits_GEO) @test isDimensional(x1_ND) == false - @test x1_ND.η0 * 1.0 == 1e5 + @test x1_ND.η0 * 1.0 == 1.0e5 @test x1_ND.τ0 * 1.0 ≈ 10.0 - @test x1_ND.ηr * 1.0 == 1e1 + @test x1_ND.ηr * 1.0 == 1.0e1 @test x1_ND.Q * 1.0 == 0.0 @test x1_ND.Tr * 1.0 ≈ 1.0 rtol = 1.0e-3 @test dimensionalize(x1_ND, CharUnits_GEO) == x1_D @@ -375,73 +375,73 @@ using GeoParams # below yield case - dimensional εt = 0.5e-20 * s^-1 τ1273 = compute_τII(x1_D, εt, (; T = 1273.0K)) - - @test ustrip(τ1273) ≈ 2*x1_D.η0.val*ustrip(εt) rtol=1e-4 + + @test ustrip(τ1273) ≈ 2 * x1_D.η0.val * ustrip(εt) rtol = 1.0e-4 ε1273 = compute_εII(x1_D, τ1273, (; T = 1273.0K)) - @test εt ≈ ε1273 rtol = 1e-6 # see if we can recover the strain rate - + @test εt ≈ ε1273 rtol = 1.0e-6 # see if we can recover the strain rate + # below yield case - nondimensional - εr = 0.5 * x1_ND.τ0/x1_ND.η0 # yield strain rate - εt = 1e-4 * εr # well below the yield + εr = 0.5 * x1_ND.τ0 / x1_ND.η0 # yield strain rate + εt = 1.0e-4 * εr # well below the yield τ1 = compute_τII(x1_ND, εt, (; T = 1.0)) - @test τ1 ≈ 2*x1_ND.η0.val*ustrip(εt) rtol=1e-4 + @test τ1 ≈ 2 * x1_ND.η0.val * ustrip(εt) rtol = 1.0e-4 ε1 = compute_εII(x1_ND, τ1, (; T = 1.0)) - @test εt ≈ ε1 rtol = 1e-6 # see if we can revocer the strain rate + @test εt ≈ ε1 rtol = 1.0e-6 # see if we can revocer the strain rate - εt = [1e-4* εr; 1e-4 * εr] - τ1 = [1e3, 1e3] # initialize τ1 - compute_τII!(τ1,x1_ND, εt, (; T = ones(size(εt)))) - @test τ1 ≈ [2*x1_ND.η0.val*ustrip(εt[1]),2*x1_ND.η0.val*ustrip(εt[2])] rtol=1e-4 + εt = [1.0e-4 * εr; 1.0e-4 * εr] + τ1 = [1.0e3, 1.0e3] # initialize τ1 + compute_τII!(τ1, x1_ND, εt, (; T = ones(size(εt)))) + @test τ1 ≈ [2 * x1_ND.η0.val * ustrip(εt[1]), 2 * x1_ND.η0.val * ustrip(εt[2])] rtol = 1.0e-4 # above yield case, dimensional x1_D = HerschelBulkley( - n = 3.0, - η0 = 1e24Pa*s, - τ0 = 100e6Pa, - ηr = 1e20Pa*s, - Q = 0.0K, - Tr = 1273K, + n = 3.0, + η0 = 1.0e24Pa * s, + τ0 = 100.0e6Pa, + ηr = 1.0e20Pa * s, + Q = 0.0K, + Tr = 1273K, ) - εt = 0.5e-9 * s^-1 - εr = 0.5 * x1_D.τ0.val/x1_D.η0.val # yield strain rate + εt = 0.5e-9 * s^-1 + εr = 0.5 * x1_D.τ0.val / x1_D.η0.val # yield strain rate - η_approx = 0.5*x1_D.τ0.val/ustrip(εt) + x1_D.ηr.val*(ustrip(εt)/εr)^(one(x1_D.n)/x1_D.n - 1) - τ1273 = compute_τII(x1_D, εt ; T = 1273.0K) - ε1273 = compute_εII(x1_D,τ1273; T = 1273.0K) + η_approx = 0.5 * x1_D.τ0.val / ustrip(εt) + x1_D.ηr.val * (ustrip(εt) / εr)^(one(x1_D.n) / x1_D.n - 1) + τ1273 = compute_τII(x1_D, εt; T = 1273.0K) + ε1273 = compute_εII(x1_D, τ1273; T = 1273.0K) - @test ustrip(τ1273) ≈ 2*η_approx*ustrip(εt) rtol=1e-6 - @test ε1273 ≈ εt rtol = 1e-6 -# -# + @test ustrip(τ1273) ≈ 2 * η_approx * ustrip(εt) rtol = 1.0e-6 + @test ε1273 ≈ εt rtol = 1.0e-6 + # + # x1_ND = nondimensionalize(x1_D, CharUnits_GEO) εt_ND = nondimensionalize(εt, CharUnits_GEO) -# - τ1_ND = compute_τII(x1_ND, εt_ND ; T = 1.0) - ε1_ND = compute_εII(x1_ND,τ1_ND; T = 1.0) + # + τ1_ND = compute_τII(x1_ND, εt_ND; T = 1.0) + ε1_ND = compute_εII(x1_ND, τ1_ND; T = 1.0) - @test ε1_ND ≈ εt_ND rtol = 1e-6 + @test ε1_ND ≈ εt_ND rtol = 1.0e-6 # temperature dependence test x1_Q = HerschelBulkley( - n = 3.0, - η0 = 1e24Pa*s, - τ0 = 100e6Pa, - ηr = 1e20Pa*s, - Q = 63772.0K, # E/R for olivine, E ≈ 530 kJ/mol + n = 3.0, + η0 = 1.0e24Pa * s, + τ0 = 100.0e6Pa, + ηr = 1.0e20Pa * s, + Q = 63772.0K, # E/R for olivine, E ≈ 530 kJ/mol Tr = 1273K, ) # compute ηT at both temperatures analytically - Q_val = 63772.0 + Q_val = 63772.0 Tr_val = 1273.0 T1_val = 1273.0 # reference temperature → ηT = ηr T2_val = 1473.0 # higher temperature → ηT < ηr - ηT1 = x1_Q.ηr.val * exp(Q_val * (1/T1_val - 1/Tr_val)) # = 1e20 - ηT2 = x1_Q.ηr.val * exp(Q_val * (1/T2_val - 1/Tr_val)) # << 1e20 + ηT1 = x1_Q.ηr.val * exp(Q_val * (1 / T1_val - 1 / Tr_val)) # = 1e20 + ηT2 = x1_Q.ηr.val * exp(Q_val * (1 / T2_val - 1 / Tr_val)) # << 1e20 εr_Q = 0.5 * x1_Q.τ0.val / x1_Q.η0.val # reference strain rate @@ -449,40 +449,40 @@ using GeoParams εt_Q = 0.5e-9 * s^-1 # well above yield # at T = 1273K (reference): ηT = ηr, same as Q=0 case - τ_T1 = compute_τII(x1_Q, εt_Q; T = T1_val*K) - ε_T1 = compute_εII(x1_Q, τ_T1; T = T1_val*K) - η_approx_T1 = 0.5*x1_Q.τ0.val/ustrip(εt_Q) + ηT1*(ustrip(εt_Q)/εr_Q)^(one(x1_Q.n)/x1_Q.n - 1) - @test ustrip(τ_T1) ≈ 2*η_approx_T1*ustrip(εt_Q) rtol = 1e-6 - @test ε_T1 ≈ εt_Q rtol = 1e-6 + τ_T1 = compute_τII(x1_Q, εt_Q; T = T1_val * K) + ε_T1 = compute_εII(x1_Q, τ_T1; T = T1_val * K) + η_approx_T1 = 0.5 * x1_Q.τ0.val / ustrip(εt_Q) + ηT1 * (ustrip(εt_Q) / εr_Q)^(one(x1_Q.n) / x1_Q.n - 1) + @test ustrip(τ_T1) ≈ 2 * η_approx_T1 * ustrip(εt_Q) rtol = 1.0e-6 + @test ε_T1 ≈ εt_Q rtol = 1.0e-6 # at T = 1473K: ηT < ηr → lower viscosity → lower stress at same strain rate - τ_T2 = compute_τII(x1_Q, εt_Q; T = T2_val*K) - ε_T2 = compute_εII(x1_Q, τ_T2; T = T2_val*K) - η_approx_T2 = 0.5*x1_Q.τ0.val/ustrip(εt_Q) + ηT2*(ustrip(εt_Q)/εr_Q)^(one(x1_Q.n)/x1_Q.n - 1) - @test ustrip(τ_T2) ≈ 2*η_approx_T2*ustrip(εt_Q) rtol = 1e-6 - @test ε_T2 ≈ εt_Q rtol = 1e-6 + τ_T2 = compute_τII(x1_Q, εt_Q; T = T2_val * K) + ε_T2 = compute_εII(x1_Q, τ_T2; T = T2_val * K) + η_approx_T2 = 0.5 * x1_Q.τ0.val / ustrip(εt_Q) + ηT2 * (ustrip(εt_Q) / εr_Q)^(one(x1_Q.n) / x1_Q.n - 1) + @test ustrip(τ_T2) ≈ 2 * η_approx_T2 * ustrip(εt_Q) rtol = 1.0e-6 + @test ε_T2 ≈ εt_Q rtol = 1.0e-6 # monotonicity: higher T → lower τII at same εII (thermal softening) @test ustrip(τ_T2) < ustrip(τ_T1) # monotonicity: higher T → higher εII at same τII (thermal softening) - τ_fixed = compute_τII(x1_Q, εt_Q; T = T1_val*K) - ε_at_T1 = compute_εII(x1_Q, τ_fixed; T = T1_val*K) - ε_at_T2 = compute_εII(x1_Q, τ_fixed; T = T2_val*K) + τ_fixed = compute_τII(x1_Q, εt_Q; T = T1_val * K) + ε_at_T1 = compute_εII(x1_Q, τ_fixed; T = T1_val * K) + ε_at_T2 = compute_εII(x1_Q, τ_fixed; T = T2_val * K) @test ustrip(ε_at_T2) > ustrip(ε_at_T1) # nondimensional temperature dependence x1_Q_ND = nondimensionalize(x1_Q, CharUnits_GEO) - εt_ND = nondimensionalize(εt_Q, CharUnits_GEO) - T1_ND = nondimensionalize(T1_val*K, CharUnits_GEO) - T2_ND = nondimensionalize(T2_val*K, CharUnits_GEO) + εt_ND = nondimensionalize(εt_Q, CharUnits_GEO) + T1_ND = nondimensionalize(T1_val * K, CharUnits_GEO) + T2_ND = nondimensionalize(T2_val * K, CharUnits_GEO) τ_T1_ND = compute_τII(x1_Q_ND, εt_ND; T = T1_ND) τ_T2_ND = compute_τII(x1_Q_ND, εt_ND; T = T2_ND) # ND results should be consistent with dimensional results after rescaling - @test ustrip(nondimensionalize(τ_T1, CharUnits_GEO)) ≈ τ_T1_ND rtol = 1e-6 - @test ustrip(nondimensionalize(τ_T2, CharUnits_GEO)) ≈ τ_T2_ND rtol = 1e-6 + @test ustrip(nondimensionalize(τ_T1, CharUnits_GEO)) ≈ τ_T1_ND rtol = 1.0e-6 + @test ustrip(nondimensionalize(τ_T2, CharUnits_GEO)) ≈ τ_T2_ND rtol = 1.0e-6 # thermal softening should hold in ND as well @test τ_T2_ND < τ_T1_ND From 269e0cc5051eadbc4c4e600f57d3f3c2edc9b231 Mon Sep 17 00:00:00 2001 From: Marcel Thielmann Date: Fri, 27 Mar 2026 13:37:14 +0100 Subject: [PATCH 10/10] removed dead code --- src/Viscosity/Viscosity.jl | 66 -------------------------------------- 1 file changed, 66 deletions(-) diff --git a/src/Viscosity/Viscosity.jl b/src/Viscosity/Viscosity.jl index ceb993a62..5fb72067c 100644 --- a/src/Viscosity/Viscosity.jl +++ b/src/Viscosity/Viscosity.jl @@ -20,72 +20,6 @@ end @inline compute_viscosity_εII(v::ConstantElasticity, εII, args) = v.G * args.dt @inline compute_viscosity_εII(v::HerschelBulkley, εII, args) = compute_hb_viscosity_εII(v, εII; args...) -#= -@inline function compute_viscosity_εII(v::HerschelBulkley, εII; T = 1.0, kwargs...) - η0, τ0, ηr, Q, Tr = if εII isa Quantity - @unpack_units η0, τ0, ηr, Q, Tr = v - η0, τ0, ηr, Q, Tr - else - @unpack_val η0, τ0, ηr, Q, Tr = v - η0, τ0, ηr, Q, Tr - end - (; n) = v - - ηT = ηr * exp(Q * (1/T-1/Tr)) # temperature dependence - εr = 0.5 * τ0/η0 # strain rate at which the Bingham yield stress is reached, this is defined as the reference strain rate - η = @pow (1.0 - exp(-2.0*η0*εII/τ0)) * (0.5*τ0/εII + ηT*(εII/εr)^(1/n - 1)) - return η -end - -@inline function compute_viscosity_τII(v::HerschelBulkley, τII; T = one(precision(a)), kwargs...) - - η0, τ0, ηr, Q, Tr = if τII isa Quantity - @unpack_units η0, τ0, ηr, Q, Tr = v - η0, τ0, ηr, Q, Tr - else - @unpack_val η0, τ0, ηr, Q, Tr = v - η0, τ0, ηr, Q, Tr - end - (; n) = v - - ηT = ηr * exp(Q * (1/T-1/Tr)) # temperature dependence - εr = 0.5 * τ0/η0 - - # define residual function --> could we also call compute_viscosity_εII here? - function fres(εII,τII,η0,τ0,ηT,n,εr) - η = @pow (1.0 - exp(-2.0*η0*εII/τ0)) * (0.5*τ0/εII + ηT*(εII/εr)^(1/n - 1)) - return 2.0 * η * εII - τII - end - - # define an initial guess for η - if τII < τ0 - η = η0 - εII = 0.5*τII/η0 - elseif τII > τ0 - εII = ((τII - τ0) / (2 * ηr * εr^(1 - one(n)/n)))^n - η = 0.5*τII/εII - end - - # Newton-Raphson method to get a solution for η - it_max = 100 - tol = 1e-8 - - # compute the residual of the initial guess - res = fres(εII,τII,η0,τ0,ηT,n,εr) - - if abs(res) > tol - for _ in 1:it_max - f, dfdε = value_and_partial(εII -> fres(εII,τII,η0,τ0,ηT,n,εr), εII) # value and dF/dεII via ForwardDiff - Δε = f / dfdε - εII -= Δε # adapt εII - abs(Δε) / (abs(εII) + eps(typeof(εII))) < tol && break # stop if we are below the tolerance - end - η = @pow (1.0 - exp(-2.0*η0*εII/τ0)) * (0.5*τ0/εII + ηT*(εII/εr)^(1/n - 1)) # --> could we also call compute_viscosity_εII here? - end - return η -end - -=# # compute effective "creep" viscosity from deviatoric stress tensor