Mixed softening and new law

src/ConstitutiveRelationships.jl

@@ -62,6 +62,7 @@ export param_info,
     iselastic,
     AbstractSoftening,
     NoSoftening,
+    DecaySoftening,
     LinearSoftening,
     NonLinearSoftening
 

src/GeoParams.jl

@@ -201,6 +201,7 @@
     NoSoftening,
     LinearSoftening,
     NonLinearSoftening,
+    DecaySoftening,
 
     #       Plasticity
     AbstractPlasticity,
@@ -443,6 +444,7 @@
         @inline $(fun)(c::CompositeRheology) = $(fun)(isviscoelastic(c), c)
         @inline $(fun)(::ElasticRheologyTrait, c::CompositeRheology) = mapreduce(x->$(fun)(x), +, c.elements)
         @inline $(fun)(::AbstractCreepLaw) = 0
+        @inline $(fun)(::AbstractPlasticity) = 0
         @inline $(fun)(r::AbstractMaterialParamsStruct) = $(fun)(r.CompositeRheology[1])
         @inline $(fun)(a::NTuple{N, AbstractMaterialParamsStruct}, phase) where N = nphase($(fun), phase, a)
     end

src/Plasticity/DruckerPrager.jl

@@ -54,9 +54,9 @@
 end
 
 # Calculation routines
-function (s::DruckerPrager{_T, U, U1, S, S})(;
+function (s::DruckerPrager{_T, U, U1, S1, S2})(;
     P=0.0, τII=0.0, Pf=0.0, EII=0.0, perturbation_C = 1.0, kwargs...
-) where {_T,U,U1,S<:AbstractSoftening}
+) where {_T,U,U1,S1<:AbstractSoftening,S2<:AbstractSoftening}
     @unpack_val sinϕ, cosϕ, ϕ, C = s
     ϕ = s.softening_ϕ(EII, ϕ)
     C = s.softening_C(EII, C)
@@ -69,9 +69,9 @@
 end
 
 
-function (s::DruckerPrager{_T, U, U1, NoSoftening, S})(;
+function (s::DruckerPrager{_T, U, U1, NoSoftening, AbstractSoftening})(;
     P=0.0, τII=0.0, Pf=0.0, EII=0.0, perturbation_C = 1.0, kwargs...
-) where {_T,U,U1,S}
+) where {_T,U,U1}
     @unpack_val sinϕ, cosϕ, ϕ, C = s
     C = s.softening_C(EII, C)
     C *=  perturbation_C
@@ -80,9 +80,9 @@
     return F
 end
 
-function (s::DruckerPrager{_T, U, U1, S, NoSoftening})(;
+function (s::DruckerPrager{_T, U, U1, AbstractSoftening, NoSoftening})(;
     P=0.0, τII=0.0, Pf=0.0, EII=0.0, perturbation_C = 1.0, kwargs...
-) where {_T,U,U1,S}
+) where {_T,U,U1}
     @unpack_val sinϕ, cosϕ, ϕ, C = s
     ϕ = s.softening_ϕ(EII, ϕ)
     C *=  perturbation_C

src/Plasticity/DruckerPrager_regularised.jl

@@ -55,9 +55,9 @@
 end
 
 # Calculation routines
-function (s::DruckerPrager_regularised{_T})(;
+function (s::DruckerPrager_regularised{T, U, U1, U2, S1, S2})(;
     P=0.0, τII=0.0, Pf=0.0, λ= 0.0, EII=0.0, perturbation_C = 1.0, kwargs...
-) where {_T}
+) where {T, U, U1, U2, S1<:AbstractSoftening, S2<:AbstractSoftening}
     @unpack_val sinϕ, cosϕ, ϕ, C, η_vp = s
     ϕ = s.softening_ϕ(EII, ϕ)
     C = s.softening_C(EII, C)
@@ -67,19 +67,19 @@
     return F
 end
 
-function (s::DruckerPrager_regularised{_T,U,U1,NoSoftening, S})(;
+function (s::DruckerPrager_regularised{_T,U,U1,NoSoftening, AbstractSoftening})(;
     P=0.0, τII=0.0, Pf=0.0, λ= 0.0, EII=0.0, perturbation_C = 1.0, kwargs...
-) where {_T,U,U1,S}
+) where {_T,U,U1,AbstractSoftening}
     @unpack_val sinϕ, cosϕ, ϕ, C, η_vp = s
     C = s.softening_C(EII, C)
     ε̇II_pl = λ*∂Q∂τII(s, τII)  # plastic strainrate
     F = τII - cosϕ * C - sinϕ * (P - Pf)  - 2 * η_vp * ε̇II_pl # with fluid pressure (set to zero by default)
     return F
 end
 
-function (s::DruckerPrager_regularised{_T,U,U1,S, NoSoftening})(;
+function (s::DruckerPrager_regularised{_T,U,U1,AbstractSoftening, NoSoftening})(;
     P=0.0, τII=0.0, Pf=0.0, λ= 0.0, EII=0.0, perturbation_C = 1.0, kwargs...
-) where {_T,U,U1,S}
+) where {_T,U,U1}
     @unpack_val sinϕ, cosϕ, ϕ, C, η_vp = s
     ϕ = s.softening_ϕ(EII, ϕ)
     sinϕ, cosϕ = iszero(EII) ? (sinϕ, cosϕ) : sincosd(ϕ)

src/Softening/Softening.jl

@@ -1,5 +1,6 @@
 
 abstract type AbstractSoftening end
+abstract type AbstractNoSoftening end
 
 struct NoSoftening <: AbstractSoftening end
 
@@ -35,7 +36,7 @@
 # (Duretz et al 2021; https://agupubs.onlinelibrary.wiley.com/doi/pdfdirect/10.1029/2021GC009675)
 using SpecialFunctions
 
-@with_kw struct NonLinearSoftening{T} <: AbstractSoftening
+@with_kw_noshow struct NonLinearSoftening{T} <: AbstractSoftening
     ξ₀::T= 0.0 # maximum value
     Δ::T = 0.0 # amplitude of the softening (i.e. minimum value)
     μ::T = 1.0 # mean of the softening
@@ -49,3 +50,16 @@
 @inline function (softening::NonLinearSoftening)(softening_var::T, args::Vararg{Any, N}) where {T,N}
     return softening.ξ₀ - 0.5 * softening.Δ * erfc(- (softening_var - softening.μ) / softening.σ)
 end
+
+
+# Non linear softening from Taras
+@with_kw_noshow struct DecaySoftening{T} <: AbstractSoftening
+    εref::T = 1e-13
+    n::T = 0.1
+end
+
+@inline (softening::DecaySoftening)(args::Vararg{Any, N}) where N = softening(promote(args...)...)
+
+@inline function (softening::DecaySoftening)(softening_var::T, max_value::T) where T
+    return max_value * inv((softening_var / softening.εref + 1)^softening.n)
+end

test/test_Softening.jl

@@ -30,6 +30,7 @@ using GeoParams, Test
     min_v, max_v = 15e0, 30e0
     lo, hi = 0.0, 1.0
     soft_ϕ = LinearSoftening(min_v, max_v, lo, hi)
+    softening_C = LinearSoftening((0e0, 10e6), (lo, hi))
 
     τII = 20e6
     P = 1e6
@@ -42,11 +43,15 @@ using GeoParams, Test
 
     p1 = DruckerPrager(; softening_ϕ = soft_ϕ)
     p2 = DruckerPrager(; softening_C = LinearSoftening((0e0, 10e6), (lo, hi)))
-    p3 = DruckerPrager(; softening_ϕ = soft_ϕ, softening_C = LinearSoftening((0e0, 10e6), (lo, hi)))
+    p3 = DruckerPrager(; softening_ϕ = soft_ϕ, softening_C = softening_C)
+    p4 = DruckerPrager(; softening_ϕ = DecaySoftening())
+    p5 = DruckerPrager(; softening_C = softening_C, softening_ϕ = DecaySoftening())
 
     @test compute_yieldfunction(p1, args) ≈ 1.0081922692006797e7
     @test compute_yieldfunction(p2, args) ≈ 1.95e7
     @test compute_yieldfunction(p3, args) ≈ 1.974118095489748e7
+    @test compute_yieldfunction(p4, args) ≈ 9.977203951679487e6
+    @test compute_yieldfunction(p5, args) ≈ 1.997376090963723e7
 
     # test regularized Drucker-Prager with softening
     p = DruckerPrager_regularised()
@@ -55,17 +60,17 @@ using GeoParams, Test
 
     args = (P=P, τII=τII, EII=1e0)
 
-    p1 = DruckerPrager(; softening_ϕ = soft_ϕ)
-    p2 = DruckerPrager(; softening_C = LinearSoftening((0e0, 10e6), (lo, hi)))
-    p3 = DruckerPrager(; softening_ϕ = soft_ϕ, softening_C = LinearSoftening((0e0, 10e6), (lo, hi)))
+    p1 = DruckerPrager_regularised(; softening_ϕ = soft_ϕ)
+    p2 = DruckerPrager_regularised(; softening_C = LinearSoftening((0e0, 10e6), (lo, hi)))
+    p3 = DruckerPrager_regularised(; softening_ϕ = soft_ϕ, softening_C = LinearSoftening((0e0, 10e6), (lo, hi)))
 
     @test compute_yieldfunction(p1, args) ≈ 1.0081922692006797e7
     @test compute_yieldfunction(p2, args) ≈ 1.95e7
     @test compute_yieldfunction(p3, args) ≈ 1.974118095489748e7
 
     # non linear softening
-    p4 = DruckerPrager(; softening_C = NonLinearSoftening())
-    p5 = DruckerPrager(; softening_C = NonLinearSoftening(ξ₀ = 30, Δ = 10))
+    p4 = DruckerPrager_regularised(; softening_C = NonLinearSoftening())
+    p5 = DruckerPrager_regularised(; softening_C = NonLinearSoftening(ξ₀ = 30, Δ = 10))
 
     @test compute_yieldfunction(p4, args) ≈ 1.95e7
     @test compute_yieldfunction(p5, (P=P, τII=τII, EII=0e0)) ≈ 1.9499974039493073e7