avoid plastic correction if iszero(tau_II)

src/stokes/Stokes2D.jl

@@ -571,6 +571,8 @@ function _solve!(
 
             @parallel (@idx ni.+1) update_stresses_center_vertex_ps!(
                 @strain(stokes),
+                @tensor_center(stokes.ε_pl),
+                stokes.EII_pl,
                 @tensor_center(stokes.τ),
                 (stokes.τ.xy,),
                 @tensor_center(stokes.τ_o),

src/stokes/StressKernels.jl

@@ -499,13 +499,13 @@ function update_stress!(
     return nothing
 end
 
-
-
 #####
 @parallel_indices (I...) function update_stresses_center_vertex_ps!(
     ε::NTuple{N, T},          # normal components @ centers; shear components @ vertices
+    ε_pl::NTuple{N, T},       # whole Voigt tensor @ centers
+    EII,                      # accumulated plastic strain rate @ centers
     τ::NTuple{N, T},          # whole Voigt tensor @ centers
-    τshear_v::NTuple{N2, T},   # shear tensor components @ vertices
+    τshear_v::NTuple{N2, T},  # shear tensor components @ vertices
     τ_o::NTuple{N, T},
     τshear_ov::NTuple{N2, T}, # shear tensor components @ vertices
     Pr,
@@ -523,22 +523,10 @@ end
     phase_vertex,
 ) where {N, N2, T}
 
-
     τxyv     = τshear_v[1]
     τxyv_old = τshear_ov[1]
     ni       = size(Pr)
-    # for j in axes(Pr,2), i in axes(Pr,1)
-    #     I   = i,j
-    # I   = i,j
-    Ic  = clamped_indices(ni, I...)
-    ## vertex
-    phase  = @inbounds phase_vertex[I...]
-    _Gvdt  = inv(fn_ratio(get_shear_modulus, rheology, phase) * dt)
-    is_pl, Cv, sinϕv, cosϕv, sinψv, η_regv = plastic_params_phase(rheology, 0e0, phase)
-    Kv     = fn_ratio(get_bulk_modulus, rheology, phase)
-    volumev= isinf(Kv) ? 0.0 : Kv * dt * sinϕv * sinψv # plastic volumetric change K * dt * sinϕ * sinψ
-    ηv_ij  = av_clamped_indices(η, Ic...)
-    dτ_rv  = 1.0/(θ_dτ + ηv_ij * _Gvdt + 1.0)
+    Ic       = clamped_indices(ni, I...)
 
     # interpolate to ith vertex
     Pv_ij       = av_clamped_indices(Pr, Ic...)
@@ -548,6 +536,16 @@ end
     τyyv_ij     = av_clamped_indices(τ[2], Ic...)
     τxxv_old_ij = av_clamped_indices(τ_o[1], Ic...)
     τyyv_old_ij = av_clamped_indices(τ_o[2], Ic...)
+    EIIv_ij     = av_clamped_indices(EII, Ic...)
+
+    ## vertex
+    phase  = @inbounds phase_vertex[I...]
+    is_pl, Cv, sinϕv, cosϕv, sinψv, η_regv = plastic_params_phase(rheology, EIIv_ij, phase)
+    _Gvdt  = inv(fn_ratio(get_shear_modulus, rheology, phase) * dt)
+    Kv     = fn_ratio(get_bulk_modulus, rheology, phase)
+    volumev= isinf(Kv) ? 0.0 : Kv * dt * sinϕv * sinψv # plastic volumetric change K * dt * sinϕ * sinψ
+    ηv_ij  = av_clamped_indices(η, Ic...)
+    dτ_rv  = 1.0/(θ_dτ + ηv_ij * _Gvdt + 1.0)
 
     # stress increments @ vertex
     dτxxv   = (-(τxxv_ij - τxxv_old_ij) * ηv_ij * _Gvdt - τxxv_ij + 2.0 * ηv_ij * εxxv_ij) * dτ_rv
@@ -557,26 +555,23 @@ end
 
     # yield function @ center
     Fv = τIIv_ij - Cv - Pv_ij *sinϕv
-    if is_pl #&& F > 0
+    if is_pl && !iszero(τIIv_ij)
         # stress correction @ vertex
         λv[I...]    = (1.0 - relλ) * λv[I...] + relλ*(max(Fv, 0.0) / (ηv_ij * dτ_rv + η_regv + volumev))
         dQdτxy      = 0.5 * (τxyv[I...] + dτxyv) / τIIv_ij
         τxyv[I...] += dτxyv - 2.0 * ηv_ij  * 0.5 * λv[I...] * dQdτxy * dτ_rv
         # τxyv[I...] += (-(τxyv[I...] - τxyv_old[I...] ) * ηv_ij * _Gvdt - τxyv[I...]  + 2.0 * ηv_ij  *(ε[3][I...] - 0.5 * λv[I...] * dQdτxy)) * dτ_rv
-
-
     else
         # stress correction @ vertex
         τxyv[I...] += dτxyv
     end
-    # end
 
     ## center
     if all(I .≤ ni)
         # Material properties
         phase = @inbounds phase_center[I...]
         _Gdt  = inv(fn_ratio(get_shear_modulus, rheology, phase) * dt)
-        is_pl, C, sinϕ, cosϕ, sinψ, η_reg = plastic_params_phase(rheology, 0e0, phase)
+        is_pl, C, sinϕ, cosϕ, sinψ, η_reg = plastic_params_phase(rheology, EII[I...], phase)
         K      = fn_ratio(get_bulk_modulus, rheology, phase)
         volume = isinf(K) ? 0.0 : K * dt * sinϕ * sinψ # plastic volumetric change K * dt * sinϕ * sinψ
         ηij    = η[I...]
@@ -594,21 +589,22 @@ end
         # yield function @ center
         F  = τII_ij - C - Pr[I...] *sinϕ
 
-        if is_pl #&& F > 0
+        if is_pl && !iszero(τII_ij)
             # stress correction @ center
             λ[I...]     = (1.0 - relλ)*λ[I...]  + relλ.*(max(F,0.0)/(η[I...] *dτ_r + η_reg + volume))
             dQdτij      = @. 0.5 * (τij + dτij) / τII_ij
             # dτij        = @. (-(τij - τij_o) * ηij * _Gdt - τij .+ 2.0 * ηij * (εij  - λ[I...] *dQdτij )) * dτ_r
-            dτij        = @. dτij - 2.0 * ηij *λ[I...] * dQdτij * dτ_r
+            εij_pl      =  λ[I...] .* dQdτij 
+            dτij        = @. dτij - 2.0 * ηij * εij_pl * dτ_r
             τij         = dτij .+ τij
             setindex!.(τ, τij, I...)
+            setindex!.(ε_pl, εij_pl, I...)
             τII[I...]   = GeoParams.second_invariant(τij)
-            Pr_c[I...]  = Pr[I...] + K*dt*λ[I...] * sinψ
-            η_vep[I...] = τII_ij / 2.0 / εII_ve
+            Pr_c[I...]  = Pr[I...] + K * dt * λ[I...] * sinψ
+            η_vep[I...] = 0.5 * τII_ij / εII_ve
         else
             # stress correction @ center
             setindex!.(τ, dτij .+ τij, I...)
-            Fchk[I...]  = 0e0
             η_vep[I...] = ηij
             τII[I...]   = τII_ij
         end