Skip to content

Commit 7a78854

Browse files
sbryngelsonclaude
andcommitted
Fix G_K exponential degradation in damage model
The damage factor was applied inside the stress component loop, causing G_K (and G) to be multiplied by the damage factor on every iteration. With N stress components, the effective shear modulus was reduced by damage^N instead of damage^1. Move the damage application before the loop so it is applied exactly once per cell. Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
1 parent b3506d4 commit 7a78854

1 file changed

Lines changed: 2 additions & 3 deletions

File tree

src/common/m_variables_conversion.fpp

Lines changed: 2 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -846,11 +846,11 @@ contains
846846
end if
847847

848848
if (hypoelasticity) then
849+
if (cont_damage) G_K = G_K*max((1._wp - qK_cons_vf(damage_idx)%sf(j, k, l)), 0._wp)
849850
$:GPU_LOOP(parallelism='[seq]')
850851
do i = strxb, strxe
851852
! subtracting elastic contribution for pressure calculation
852853
if (G_K > verysmall) then
853-
if (cont_damage) G_K = G_K*max((1._wp - qK_cons_vf(damage_idx)%sf(j, k, l)), 0._wp)
854854
qK_prim_vf(E_idx)%sf(j, k, l) = qK_prim_vf(E_idx)%sf(j, k, l) - &
855855
((qK_prim_vf(i)%sf(j, k, l)**2._wp)/(4._wp*G_K))/gamma_K
856856
! Double for shear stresses
@@ -1123,11 +1123,10 @@ contains
11231123
end if
11241124

11251125
if (hypoelasticity) then
1126+
if (cont_damage) G = G*max((1._wp - q_prim_vf(damage_idx)%sf(j, k, l)), 0._wp)
11261127
do i = strxb, strxe
11271128
! adding elastic contribution
11281129
if (G > verysmall) then
1129-
if (cont_damage) G = G*max((1._wp - q_prim_vf(damage_idx)%sf(j, k, l)), 0._wp)
1130-
11311130
q_cons_vf(E_idx)%sf(j, k, l) = q_cons_vf(E_idx)%sf(j, k, l) + &
11321131
(q_prim_vf(i)%sf(j, k, l)**2._wp)/(4._wp*G)
11331132
! Double for shear stresses

0 commit comments

Comments
 (0)