Skip to content

Commit 063121f

Browse files
ChrisZYJChrisZYJ
committed
remove clamp
1 parent 5452db0 commit 063121f

2 files changed

Lines changed: 55 additions & 63 deletions

File tree

‎src/common/m_variables_conversion.fpp‎

Lines changed: 13 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -734,15 +734,14 @@ contains
734734
if (cont_damage) G_K = G_K*max((1._wp - qK_cons_vf(eqn_idx%damage)%sf(j, k, l)), 0._wp)
735735
$:GPU_LOOP(parallelism='[seq]')
736736
do i = eqn_idx%stress%beg, eqn_idx%stress%end
737-
! subtracting elastic contribution for pressure calculation
738-
if (G_K > verysmall) then
737+
! subtracting elastic contribution for pressure calculation (unconditional, clamped denominator to match
738+
! Riemann solver)
739+
qK_prim_vf(eqn_idx%E)%sf(j, k, l) = qK_prim_vf(eqn_idx%E)%sf(j, k, l) - ((qK_prim_vf(i)%sf(j, k, &
740+
& l)**2._wp)/max(4._wp*G_K, verysmall))/gamma_K
741+
! Double for shear stresses
742+
if (any(i == shear_indices)) then
739743
qK_prim_vf(eqn_idx%E)%sf(j, k, l) = qK_prim_vf(eqn_idx%E)%sf(j, k, l) - ((qK_prim_vf(i)%sf(j, k, &
740-
& l)**2._wp)/(4._wp*G_K))/gamma_K
741-
! Double for shear stresses
742-
if (any(i == shear_indices)) then
743-
qK_prim_vf(eqn_idx%E)%sf(j, k, l) = qK_prim_vf(eqn_idx%E)%sf(j, k, l) - ((qK_prim_vf(i)%sf(j, &
744-
& k, l)**2._wp)/(4._wp*G_K))/gamma_K
745-
end if
744+
& l)**2._wp)/max(4._wp*G_K, verysmall))/gamma_K
746745
end if
747746
end do
748747
end if
@@ -986,15 +985,13 @@ contains
986985
if (hypoelasticity) then
987986
if (cont_damage) G = G*max((1._wp - q_prim_vf(eqn_idx%damage)%sf(j, k, l)), 0._wp)
988987
do i = eqn_idx%stress%beg, eqn_idx%stress%end
989-
! adding elastic contribution
990-
if (G > verysmall) then
988+
! adding elastic contribution (unconditional, clamped denominator)
989+
q_cons_vf(eqn_idx%E)%sf(j, k, l) = q_cons_vf(eqn_idx%E)%sf(j, k, l) + (q_prim_vf(i)%sf(j, k, &
990+
& l)**2._wp)/max(4._wp*G, verysmall)
991+
! Double for shear stresses
992+
if (any(i == shear_indices)) then
991993
q_cons_vf(eqn_idx%E)%sf(j, k, l) = q_cons_vf(eqn_idx%E)%sf(j, k, l) + (q_prim_vf(i)%sf(j, k, &
992-
& l)**2._wp)/(4._wp*G)
993-
! Double for shear stresses
994-
if (any(i == shear_indices)) then
995-
q_cons_vf(eqn_idx%E)%sf(j, k, l) = q_cons_vf(eqn_idx%E)%sf(j, k, l) + (q_prim_vf(i)%sf(j, k, &
996-
& l)**2._wp)/(4._wp*G)
997-
end if
994+
& l)**2._wp)/max(4._wp*G, verysmall)
998995
end if
999996
end do
1000997
end if

‎src/simulation/m_riemann_solvers.fpp‎

Lines changed: 42 additions & 47 deletions
Original file line numberDiff line numberDiff line change
@@ -2050,15 +2050,13 @@ contains
20502050
end do
20512051
$:GPU_LOOP(parallelism='[seq]')
20522052
do i = 1, eqn_idx%stress%end - eqn_idx%stress%beg + 1
2053-
! Elastic contribution to energy if G large enough
2054-
if ((G_L > verysmall) .and. (G_R > verysmall)) then
2055-
E_L = E_L + (tau_e_L(i)*tau_e_L(i))/(4._wp*G_L)
2056-
E_R = E_R + (tau_e_R(i)*tau_e_R(i))/(4._wp*G_R)
2057-
! Additional terms in 2D and 3D
2058-
if ((i == 2) .or. (i == 4) .or. (i == 5)) then
2059-
E_L = E_L + (tau_e_L(i)*tau_e_L(i))/(4._wp*G_L)
2060-
E_R = E_R + (tau_e_R(i)*tau_e_R(i))/(4._wp*G_R)
2061-
end if
2053+
! Elastic contribution to energy (unconditional, clamped denominator)
2054+
E_L = E_L + (tau_e_L(i)*tau_e_L(i))/max(4._wp*G_L, verysmall)
2055+
E_R = E_R + (tau_e_R(i)*tau_e_R(i))/max(4._wp*G_R, verysmall)
2056+
! Additional terms in 2D and 3D
2057+
if ((i == 2) .or. (i == 4) .or. (i == 5)) then
2058+
E_L = E_L + (tau_e_L(i)*tau_e_L(i))/max(4._wp*G_L, verysmall)
2059+
E_R = E_R + (tau_e_R(i)*tau_e_R(i))/max(4._wp*G_R, verysmall)
20622060
end if
20632061
end do
20642062
end if
@@ -2121,14 +2119,14 @@ contains
21212119
if (wave_speeds == 1) then
21222120
if (elasticity) then
21232121
! Elastic wave speed, Rodriguez et al. JCP (2019)
2124-
s_L = min(vel_L(dir_idx(1)) - sqrt(c_L*c_L + (((4._wp*G_L)/3._wp) + tau_e_L(dir_idx_tau(1) &
2125-
& ))/rho_L), &
2126-
& vel_R(dir_idx(1)) - sqrt(c_R*c_R + (((4._wp*G_R)/3._wp) &
2127-
& + tau_e_R(dir_idx_tau(1)))/rho_R))
2128-
s_R = max(vel_R(dir_idx(1)) + sqrt(c_R*c_R + (((4._wp*G_R)/3._wp) + tau_e_R(dir_idx_tau(1) &
2129-
& ))/rho_R), &
2130-
& vel_L(dir_idx(1)) + sqrt(c_L*c_L + (((4._wp*G_L)/3._wp) &
2131-
& + tau_e_L(dir_idx_tau(1)))/rho_L))
2122+
s_L = min(vel_L(dir_idx(1)) - sqrt(max(verysmall, &
2123+
& c_L*c_L + (((4._wp*G_L)/3._wp) + tau_e_L(dir_idx_tau(1)))/rho_L)), &
2124+
& vel_R(dir_idx(1)) - sqrt(max(verysmall, &
2125+
& c_R*c_R + (((4._wp*G_R)/3._wp) + tau_e_R(dir_idx_tau(1)))/rho_R)))
2126+
s_R = max(vel_R(dir_idx(1)) + sqrt(max(verysmall, &
2127+
& c_R*c_R + (((4._wp*G_R)/3._wp) + tau_e_R(dir_idx_tau(1)))/rho_R)), &
2128+
& vel_L(dir_idx(1)) + sqrt(max(verysmall, &
2129+
& c_L*c_L + (((4._wp*G_L)/3._wp) + tau_e_L(dir_idx_tau(1)))/rho_L)))
21322130
s_S = (pres_R - tau_e_R(dir_idx_tau(1)) - pres_L + tau_e_L(dir_idx_tau(1)) &
21332131
& + rho_L*vel_L(dir_idx(1))*(s_L - vel_L(dir_idx(1))) - rho_R*vel_R(dir_idx(1)) &
21342132
& *(s_R - vel_R(dir_idx(1))))/(rho_L*(s_L - vel_L(dir_idx(1))) - rho_R*(s_R &
@@ -3211,15 +3209,13 @@ contains
32113209
end do
32123210
$:GPU_LOOP(parallelism='[seq]')
32133211
do i = 1, eqn_idx%stress%end - eqn_idx%stress%beg + 1
3214-
! Elastic contribution to energy if G large enough
3215-
if ((G_L > verysmall) .and. (G_R > verysmall)) then
3216-
E_L = E_L + (tau_e_L(i)*tau_e_L(i))/(4._wp*G_L)
3217-
E_R = E_R + (tau_e_R(i)*tau_e_R(i))/(4._wp*G_R)
3218-
! Additional terms in 2D and 3D
3219-
if ((i == 2) .or. (i == 4) .or. (i == 5)) then
3220-
E_L = E_L + (tau_e_L(i)*tau_e_L(i))/(4._wp*G_L)
3221-
E_R = E_R + (tau_e_R(i)*tau_e_R(i))/(4._wp*G_R)
3222-
end if
3212+
! Elastic contribution to energy (unconditional, clamped denominator)
3213+
E_L = E_L + (tau_e_L(i)*tau_e_L(i))/max(4._wp*G_L, verysmall)
3214+
E_R = E_R + (tau_e_R(i)*tau_e_R(i))/max(4._wp*G_R, verysmall)
3215+
! Additional terms in 2D and 3D
3216+
if ((i == 2) .or. (i == 4) .or. (i == 5)) then
3217+
E_L = E_L + (tau_e_L(i)*tau_e_L(i))/max(4._wp*G_L, verysmall)
3218+
E_R = E_R + (tau_e_R(i)*tau_e_R(i))/max(4._wp*G_R, verysmall)
32233219
end if
32243220
end do
32253221
end if
@@ -3285,14 +3281,14 @@ contains
32853281
if (wave_speeds == 1) then
32863282
if (elasticity) then
32873283
! Elastic wave speed, Rodriguez et al. JCP (2019)
3288-
s_L = min(vel_L(dir_idx(1)) - sqrt(c_L*c_L + (((4._wp*G_L)/3._wp) + tau_e_L(dir_idx_tau(1) &
3289-
& ))/rho_L), &
3290-
& vel_R(dir_idx(1)) - sqrt(c_R*c_R + (((4._wp*G_R)/3._wp) &
3291-
& + tau_e_R(dir_idx_tau(1)))/rho_R))
3292-
s_R = max(vel_R(dir_idx(1)) + sqrt(c_R*c_R + (((4._wp*G_R)/3._wp) + tau_e_R(dir_idx_tau(1) &
3293-
& ))/rho_R), &
3294-
& vel_L(dir_idx(1)) + sqrt(c_L*c_L + (((4._wp*G_L)/3._wp) &
3295-
& + tau_e_L(dir_idx_tau(1)))/rho_L))
3284+
s_L = min(vel_L(dir_idx(1)) - sqrt(max(verysmall, &
3285+
& c_L*c_L + (((4._wp*G_L)/3._wp) + tau_e_L(dir_idx_tau(1)))/rho_L)), &
3286+
& vel_R(dir_idx(1)) - sqrt(max(verysmall, &
3287+
& c_R*c_R + (((4._wp*G_R)/3._wp) + tau_e_R(dir_idx_tau(1)))/rho_R)))
3288+
s_R = max(vel_R(dir_idx(1)) + sqrt(max(verysmall, &
3289+
& c_R*c_R + (((4._wp*G_R)/3._wp) + tau_e_R(dir_idx_tau(1)))/rho_R)), &
3290+
& vel_L(dir_idx(1)) + sqrt(max(verysmall, &
3291+
& c_L*c_L + (((4._wp*G_L)/3._wp) + tau_e_L(dir_idx_tau(1)))/rho_L)))
32963292
s_S = (pres_R - tau_e_R(dir_idx_tau(1)) - pres_L + tau_e_L(dir_idx_tau(1)) &
32973293
& + rho_L*vel_L(dir_idx(1))*(s_L - vel_L(dir_idx(1))) - rho_R*vel_R(dir_idx(1)) &
32983294
& *(s_R - vel_R(dir_idx(1))))/(rho_L*(s_L - vel_L(dir_idx(1))) - rho_R*(s_R &
@@ -4593,15 +4589,14 @@ contains
45934589
45944590
$:GPU_LOOP(parallelism='[seq]')
45954591
do i = 1, eqn_idx%stress%end - eqn_idx%stress%beg + 1
4596-
! Elastic contribution to energy if G large enough
4597-
if ((G_L > verysmall) .and. (G_R > verysmall)) then
4598-
E%L = E%L + (tau_e_L(i)*tau_e_L(i))/(4._wp*G_L)
4599-
E%R = E%R + (tau_e_R(i)*tau_e_R(i))/(4._wp*G_R)
4600-
! Shear terms doubled: 2D/2D-axisym i==2 only; 3D i==2,4,5
4601-
if ((n > 0 .and. p == 0 .and. i == 2) .or. (p > 0 .and. (i == 2 .or. i == 4 .or. i == 5))) then
4602-
E%L = E%L + (tau_e_L(i)*tau_e_L(i))/(4._wp*G_L)
4603-
E%R = E%R + (tau_e_R(i)*tau_e_R(i))/(4._wp*G_R)
4604-
end if
4592+
! Elastic contribution to energy if G large enough Elastic contribution to energy (unconditional,
4593+
! clamped denominator)
4594+
E%L = E%L + (tau_e_L(i)*tau_e_L(i))/max(4._wp*G_L, verysmall)
4595+
E%R = E%R + (tau_e_R(i)*tau_e_R(i))/max(4._wp*G_R, verysmall)
4596+
! Shear terms doubled: 2D/2D-axisym i==2 only; 3D i==2,4,5
4597+
if ((n > 0 .and. p == 0 .and. i == 2) .or. (p > 0 .and. (i == 2 .or. i == 4 .or. i == 5))) then
4598+
E%L = E%L + (tau_e_L(i)*tau_e_L(i))/max(4._wp*G_L, verysmall)
4599+
E%R = E%R + (tau_e_R(i)*tau_e_R(i))/max(4._wp*G_R, verysmall)
46054600
end if
46064601
end do
46074602
@@ -4615,10 +4610,10 @@ contains
46154610
call s_compute_speed_of_sound(pres%R, rho%R, gamma%R, pi_inf%R, H%R, alpha_R, vel_rms%R, 0._wp, c%R, &
46164611
& qv%R)
46174612
4618-
S_L = min(u_n_L - sqrt(c%L*c%L + ((4._wp/3._wp)*G_L + tau_nn_L)/rho%L), &
4619-
& u_n_R - sqrt(c%R*c%R + ((4._wp/3._wp)*G_R + tau_nn_R)/rho%R))
4620-
S_R = max(u_n_R + sqrt(c%R*c%R + ((4._wp/3._wp)*G_R + tau_nn_R)/rho%R), &
4621-
& u_n_L + sqrt(c%L*c%L + ((4._wp/3._wp)*G_L + tau_nn_L)/rho%L))
4613+
S_L = min(u_n_L - sqrt(max(verysmall, c%L*c%L + ((4._wp/3._wp)*G_L + tau_nn_L)/rho%L)), &
4614+
& u_n_R - sqrt(max(verysmall, c%R*c%R + ((4._wp/3._wp)*G_R + tau_nn_R)/rho%R)))
4615+
S_R = max(u_n_R + sqrt(max(verysmall, c%R*c%R + ((4._wp/3._wp)*G_R + tau_nn_R)/rho%R)), &
4616+
& u_n_L + sqrt(max(verysmall, c%L*c%L + ((4._wp/3._wp)*G_L + tau_nn_L)/rho%L)))
46224617
46234618
! Two-component 2D only (enforced by checker restrictions)
46244619

0 commit comments

Comments
 (0)