From 778163796a0c8202139824c26a6db5ece88e3d21 Mon Sep 17 00:00:00 2001 From: Albert de Montserrat Date: Fri, 20 Sep 2024 16:57:18 +0200 Subject: [PATCH] avoid plastic correction if `iszero(tau_II)` --- src/stokes/Stokes2D.jl | 2 ++ src/stokes/StressKernels.jl | 48 +++++++++++++++++-------------------- 2 files changed, 24 insertions(+), 26 deletions(-) diff --git a/src/stokes/Stokes2D.jl b/src/stokes/Stokes2D.jl index 2d1fc7a4..5cff3d81 100644 --- a/src/stokes/Stokes2D.jl +++ b/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), diff --git a/src/stokes/StressKernels.jl b/src/stokes/StressKernels.jl index 190824de..1017245b 100644 --- a/src/stokes/StressKernels.jl +++ b/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