Skip to content

Commit 1b91f13

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 52e1010 commit 1b91f13

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
@@ -852,11 +852,11 @@ contains
852852
end if
853853

854854
if (hypoelasticity) then
855+
if (cont_damage) G_K = G_K*max((1._wp - qK_cons_vf(damage_idx)%sf(j, k, l)), 0._wp)
855856
$:GPU_LOOP(parallelism='[seq]')
856857
do i = strxb, strxe
857858
! subtracting elastic contribution for pressure calculation
858859
if (G_K > verysmall) then
859-
if (cont_damage) G_K = G_K*max((1._wp - qK_cons_vf(damage_idx)%sf(j, k, l)), 0._wp)
860860
qK_prim_vf(E_idx)%sf(j, k, l) = qK_prim_vf(E_idx)%sf(j, k, l) - &
861861
((qK_prim_vf(i)%sf(j, k, l)**2._wp)/(4._wp*G_K))/gamma_K
862862
! 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)