Non linear softening (#92)

* update softening

* bump GeoParams

* update softening miniapp

* format & fma's

Project.toml

@@ -27,7 +27,7 @@ Adapt = "3.7.2"
 CUDA = "4.4.1, 5"
 CellArrays = "0.1 "
 HDF5 = "0.17.1"
-GeoParams = "0.5.4"
+GeoParams = "0.5.5"
 ImplicitGlobalGrid = "0.14.0"
 MPI = "0.20"
 MuladdMacro = "0.2"

miniapps/benchmarks/stokes2D/shear_band/ShearBand2D_softening.jl

@@ -41,7 +41,6 @@ function main(igg; nx=64, ny=64, figdir="model_figs")
     origin       = 0.0, 0.0     # origin coordinates
     grid         = Geometry(ni, li; origin = origin)
     (; xci, xvi) = grid # nodes at the center and vertices of the cells
-    dt           = Inf
 
     # Physical properties using GeoParams ----------------
     τ_y     = 1.6           # yield stress. If do_DP=true, τ_y stand for the cohesion: c*cos(ϕ)
@@ -53,10 +52,12 @@ function main(igg; nx=64, ny=64, figdir="model_figs")
     εbg     = 1.0           # background strain-rate
     η_reg   = 8e-3          # regularisation "viscosity"
     dt      = η0/G0/4.0     # assumes Maxwell time of 4
+    dt /= 5
     el_bg   = ConstantElasticity(; G=G0, Kb=4)
     el_inc  = ConstantElasticity(; G=Gi, Kb=4)
     visc    = LinearViscous(; η=η0)
-    soft_C  = LinearSoftening((C/2, C), (0e0, 1e0))
+    # soft_C  = LinearSoftening((C/2, C), (0e0, 2e0))
+    soft_C  = NonLinearSoftening(;ξ₀ = C, Δ=C/2)
     pl      = DruckerPrager_regularised(;  # non-regularized plasticity
         C    = C,
         ϕ    = ϕ,
@@ -121,7 +122,7 @@ function main(igg; nx=64, ny=64, figdir="model_figs")
 
     # Time loop
     t, it      = 0.0, 0
-    tmax       = 5
+    tmax       = 3.5
     τII        = Float64[]
     sol        = Float64[]
     ttot       = Float64[]
@@ -191,7 +192,7 @@ end
 n      = 128
 nx     = n
 ny     = n
-figdir = "ShearBands2D_softening"
+figdir = "ShearBands2D_lin_softening"
 igg  = if !(JustRelax.MPI.Initialized())
     IGG(init_global_grid(nx, ny, 1; init_MPI = true)...)
 else

src/rheology/StressUpdate.jl

@@ -23,7 +23,8 @@ function _compute_τ_nonlinear!(
 
     # visco-elastic strain rates
     εij_ve = ntuple(Val(N1)) do i
-        @inbounds muladd(0.5 * τij_o[i], _Gdt, εij[i])
+        Base.@_inline_meta
+        fma(0.5 * τij_o[i], _Gdt, εij[i])
     end
     # get plastic parameters (if any...)
     (; is_pl, C, sinϕ, cosϕ, η_reg, volume) = plastic_parameters
@@ -47,11 +48,10 @@ function _compute_τ_nonlinear!(
     end
 
     # fill plastic strain rate tensor
-
     update_plastic_strain_rate!(ε_pl, λdQdτ, idx)
     # update and correct stress
-    τij = τij .+ dτij
-    correct_stress!(τ, τij, idx...)
+    correct_stress!(τ, τij .+ dτij, idx...)
+
     τII[idx...] = τII_ij = second_invariant(τij...)
     η_vep[idx...] = τII_ij * 0.5 * inv(second_invariant(εij_ve...))
 
@@ -69,16 +69,14 @@ end
 # check if plasticity is active
 @inline isyielding(is_pl, τII_trial, τy) = is_pl * (τII_trial > τy)
 
-@inline compute_dτ_r(θ_dτ, ηij, _Gdt) = inv(θ_dτ + muladd(ηij, _Gdt, 1.0))
+@inline compute_dτ_r(θ_dτ, ηij, _Gdt) = inv(θ_dτ + fma(ηij, _Gdt, 1.0))
 
 function compute_stress_increment_and_trial(
     τij::NTuple{N,T}, τij_o::NTuple{N,T}, ηij, εij::NTuple{N,T}, _Gdt, dτ_r
 ) where {N,T}
     dτij = ntuple(Val(N)) do i
         Base.@_inline_meta
-        @inbounds dτ_r * muladd(
-            2.0 * ηij, εij[i], muladd(-((τij[i] - τij_o[i])) * ηij, _Gdt, -τij[i])
-        )
+        dτ_r * fma(2.0 * ηij, εij[i], fma(-((τij[i] - τij_o[i])) * ηij, _Gdt, -τij[i]))
     end
     return dτij, second_invariant((τij .+ dτij)...)
 end
@@ -102,7 +100,7 @@ function compute_dτ_pl(
     dτ_pl = ntuple(Val(N)) do i
         Base.@_inline_meta
         # corrected stress
-        muladd(-dτ_r * 2.0, ηij * λdQdτ[i], dτij[i])
+        fma(-dτ_r * 2.0, ηij * λdQdτ[i], dτij[i])
     end
     return dτ_pl, λ, λdQdτ
 end
@@ -113,7 +111,7 @@ end
 ) where {N1,N2,T}
     quote
         Base.@_inline_meta
-        Base.@nexprs $N1 i -> @inbounds τ[i][idx...] = τij[i]
+        Base.@nexprs $N1 i -> τ[i][idx...] = τij[i]
     end
 end
 
@@ -226,13 +224,13 @@ end
 @inline function soften_cohesion(
     v::DruckerPrager{T,U,U1,S1,S2}, EII::T
 ) where {T,U,U1,S1,S2}
-    return v.softening_C(v.C.val, EII)
+    return v.softening_C(EII, v.C.val)
 end
 
 @inline function soften_cohesion(
     v::DruckerPrager_regularised{T,U,U1,U2,S1,S2}, EII::T
 ) where {T,U,U1,U2,S1,S2}
-    return v.softening_C(v.C.val, EII)
+    return v.softening_C(EII, v.C.val)
 end
 
 @inline function soften_friction_angle(
@@ -250,13 +248,13 @@ end
 @inline function soften_friction_angle(
     v::DruckerPrager{T,U,U1,S1,S2}, EII::T
 ) where {T,U,U1,S1,S2}
-    ϕ = v.softening_ϕ(v.ϕ.val, EII)
-    return cosd(ϕ), sind(ϕ)
+    ϕ = v.softening_ϕ(EII, v.ϕ.val)
+    return sincosd(ϕ)
 end
 
 @inline function soften_friction_angle(
     v::DruckerPrager_regularised{T,U,U1,U2,S1,S2}, EII::T
 ) where {T,U,U1,U2,S1,S2}
-    ϕ = v.softening_ϕ(v.ϕ.val, EII)
-    return cosd(ϕ), sind(ϕ)
+    ϕ = v.softening_ϕ(EII, v.ϕ.val)
+    return sincosd(ϕ)
 end

src/stokes/PressureKernels.jl

@@ -42,7 +42,7 @@ end
 
 function _compute_P!(P, P0, ∇V, η, K, dt, r, θ_dτ)
     _Kdt = inv(K * dt)
-    RP = muladd(-(P - P0), _Kdt, -∇V)
+    RP = fma(-(P - P0), _Kdt, -∇V)
     P += RP / (1.0 / (r / θ_dτ * η) + 1.0 * _Kdt)
     return RP, P
 end

src/stokes/StressRotation.jl

@@ -23,9 +23,9 @@ end
         τ_xy = @cell xy[ip, cell...]
 
         tmp = τ_xy * ω_xy * 2.0
-        @cell xx[ip, cell...] = muladd(dt, cte, τ_xx)
-        @cell yy[ip, cell...] = muladd(dt, cte, τ_yy)
-        @cell xy[ip, cell...] = muladd(dt, (τ_xx - τ_yy) * ω_xy, τ_xy)
+        @cell xx[ip, cell...] = fma(dt, cte, τ_xx)
+        @cell yy[ip, cell...] = fma(dt, cte, τ_yy)
+        @cell xy[ip, cell...] = fma(dt, (τ_xx - τ_yy) * ω_xy, τ_xy)
     end
 
     return nothing
@@ -134,7 +134,7 @@ Base.@propagate_inbounds function rotate_stress!(
 
     ## 3) Update stress
     for k in 1:N
-        τ[k][idx...] = muladd(τij_adv[k] * 0, dt, τr_voigt[k])
+        τ[k][idx...] = fma(τij_adv[k] * 0, dt, τr_voigt[k])
     end
     return nothing
 end