Skip to content

Commit baba2a0

Browse files
ChrisZYJChrisZYJ
committed
axisym hypo
1 parent e6a7e49 commit baba2a0

File tree

7 files changed

+65
-24
lines changed

7 files changed

+65
-24
lines changed

‎src/common/m_variables_conversion.fpp‎

Lines changed: 6 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -162,10 +162,8 @@ contains
162162
do s = stress_idx%beg, stress_idx%end
163163
if (G > 0) then
164164
E_e = E_e + ((stress/rho)**2._wp)/(4._wp*G)
165-
! Additional terms in 2D and 3D
166-
if ((s == stress_idx%beg + 1) .or. &
167-
(s == stress_idx%beg + 3) .or. &
168-
(s == stress_idx%beg + 4)) then
165+
! Double for shear stresses
166+
if (any(s == shear_indices)) then
169167
E_e = E_e + ((stress/rho)**2._wp)/(4._wp*G)
170168
end if
171169
end if
@@ -1126,10 +1124,8 @@ contains
11261124
if (G_K > verysmall) then
11271125
qK_prim_vf(E_idx)%sf(j, k, l) = qK_prim_vf(E_idx)%sf(j, k, l) - &
11281126
((qK_prim_vf(i)%sf(j, k, l)**2._wp)/(4._wp*G_K))/gamma_K
1129-
! extra terms in 2 and 3D
1130-
if ((i == strxb + 1) .or. &
1131-
(i == strxb + 3) .or. &
1132-
(i == strxb + 4)) then
1127+
! Double for shear stresses
1128+
if (any(i == shear_indices)) then
11331129
qK_prim_vf(E_idx)%sf(j, k, l) = qK_prim_vf(E_idx)%sf(j, k, l) - &
11341130
((qK_prim_vf(i)%sf(j, k, l)**2._wp)/(4._wp*G_K))/gamma_K
11351131
end if
@@ -1396,10 +1392,8 @@ contains
13961392
if (G > verysmall) then
13971393
q_cons_vf(E_idx)%sf(j, k, l) = q_cons_vf(E_idx)%sf(j, k, l) + &
13981394
(q_prim_vf(i)%sf(j, k, l)**2._wp)/(4._wp*G)
1399-
! extra terms in 2 and 3D
1400-
if ((i == stress_idx%beg + 1) .or. &
1401-
(i == stress_idx%beg + 3) .or. &
1402-
(i == stress_idx%beg + 4)) then
1395+
! Double for shear stresses
1396+
if (any(i == shear_indices)) then
14031397
q_cons_vf(E_idx)%sf(j, k, l) = q_cons_vf(E_idx)%sf(j, k, l) + &
14041398
(q_prim_vf(i)%sf(j, k, l)**2._wp)/(4._wp*G)
14051399
end if

‎src/post_process/m_global_parameters.fpp‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -674,7 +674,8 @@ contains
674674
elasticity = .true.
675675
stress_idx%beg = sys_size + 1
676676
stress_idx%end = sys_size + (num_dims*(num_dims + 1))/2
677-
! number of stresses is 1 in 1D, 3 in 2D, 6 in 3D
677+
if (cyl_coord) stress_idx%end = stress_idx%end + 1
678+
! number of stresses is 1 in 1D, 3 in 2D, 4 in 2D-Axisym, 6 in 3D
678679
sys_size = stress_idx%end
679680

680681
! shear stress index is 2 for 2D and 2,4,5 for 3D

‎src/pre_process/m_global_parameters.fpp‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -767,7 +767,8 @@ contains
767767
elasticity = .true.
768768
stress_idx%beg = sys_size + 1
769769
stress_idx%end = sys_size + (num_dims*(num_dims + 1))/2
770-
! number of stresses is 1 in 1D, 3 in 2D, 6 in 3D
770+
if (cyl_coord) stress_idx%end = stress_idx%end + 1
771+
! number of stresses is 1 in 1D, 3 in 2D, 4 in 2D-Axisym, 6 in 3D
771772
sys_size = stress_idx%end
772773

773774
! shear stress index is 2 for 2D and 2,4,5 for 3D

‎src/simulation/m_checker.fpp‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -168,6 +168,8 @@ contains
168168
"Only acoustic("//trim(jStr)//")%support = 2, 5, 6, 9, or 10 is allowed for 2D simulations")
169169
@:PROHIBIT(dim == 2 .and. (.not. any(acoustic(j)%support == (/6, 10/))) .and. cyl_coord, &
170170
"Only acoustic("//trim(jStr)//")%support = 6 or 10 is allowed for 2D axisymmetric simulations")
171+
@:PROHIBIT(.not. (dim == 2 .and. cyl_coord) .and. any(acoustic(j)%support == (/6, 10/)), &
172+
"acoustic("//trim(jStr)//")%support = 6 or 10 only works for 2D axisymmetric simulations")
171173
@:PROHIBIT(dim == 2 .and. any(acoustic(j)%support == (/2, 5, 6, 9, 10/)) .and. &
172174
(f_is_default(acoustic(j)%loc(1)) .or. f_is_default(acoustic(j)%loc(2))), &
173175
"acoustic("//trim(jStr)//")%loc(1:2) must be specified for acoustic("//trim(jStr)//")%support = 2")

‎src/simulation/m_global_parameters.fpp‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1039,7 +1039,8 @@ contains
10391039
elasticity = .true.
10401040
stress_idx%beg = sys_size + 1
10411041
stress_idx%end = sys_size + (num_dims*(num_dims + 1))/2
1042-
! number of stresses is 1 in 1D, 3 in 2D, 6 in 3D
1042+
if (cyl_coord) stress_idx%end = stress_idx%end + 1
1043+
! number of stresses is 1 in 1D, 3 in 2D, 4 in 2D-Axisym, 6 in 3D
10431044
sys_size = stress_idx%end
10441045
10451046
! shear stress index is 2 for 2D and 2,4,5 for 3D

‎src/simulation/m_hypoelastic.fpp‎

Lines changed: 33 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -331,6 +331,39 @@ contains
331331
end do
332332
end if
333333

334+
if (cyl_coord) then
335+
336+
!$acc parallel loop collapse(3) gang vector default(present)
337+
do q = 0, p
338+
do l = 0, n
339+
do k = 0, m
340+
! S_xx -= rho * v/r * (tau_xx + 2/3*G)
341+
rhs_vf(strxb)%sf(k, l, q) = rhs_vf(strxb)%sf(k, l, q) - &
342+
rho_K_field(k, l, q)*q_prim_vf(momxb + 1)%sf(k, l, q)/y_cc(l)* &
343+
(q_prim_vf(strxb)%sf(k, l, q) + (2._wp/3._wp)*G_K_field(k, l, q)) ! tau_xx + 2/3*G
344+
345+
! S_xr -= rho * v/r * tau_xr
346+
rhs_vf(strxb + 1)%sf(k, l, q) = rhs_vf(strxb + 1)%sf(k, l, q) - &
347+
rho_K_field(k, l, q)*q_prim_vf(momxb + 1)%sf(k, l, q)/y_cc(l)* &
348+
q_prim_vf(strxb + 1)%sf(k, l, q) ! tau_xx
349+
350+
! S_rr -= rho * v/r * (tau_rr + 2/3*G)
351+
rhs_vf(strxb + 2)%sf(k, l, q) = rhs_vf(strxb + 2)%sf(k, l, q) - &
352+
rho_K_field(k, l, q)*q_prim_vf(momxb + 1)%sf(k, l, q)/y_cc(l)* &
353+
(q_prim_vf(strxb + 2)%sf(k, l, q) + (2._wp/3._wp)*G_K_field(k, l, q)) ! tau_rr + 2/3*G
354+
355+
! S_thetatheta += rho * ( -(tau_thetatheta + 2/3*G)*(du/dx + dv/dr + v/r) + 2*(tau_thetatheta + G)*v/r )
356+
rhs_vf(strxb + 3)%sf(k, l, q) = rhs_vf(strxb + 3)%sf(k, l, q) + &
357+
rho_K_field(k, l, q)*( &
358+
-(q_prim_vf(strxb + 3)%sf(k, l, q) + (2._wp/3._wp)*G_K_field(k, l, q))* &
359+
(du_dx(k, l, q) + dv_dy(k, l, q) + q_prim_vf(momxb + 1)%sf(k, l, q)/y_cc(l)) &
360+
+ 2._wp*(q_prim_vf(strxb + 3)%sf(k, l, q) + G_K_field(k, l, q))*q_prim_vf(momxb + 1)%sf(k, l, q)/y_cc(l))
361+
end do
362+
end do
363+
end do
364+
365+
end if
366+
334367
end subroutine s_compute_hypoelastic_rhs
335368

336369
subroutine s_finalize_hypoelastic_module()

‎src/simulation/m_riemann_solvers.fpp‎

Lines changed: 18 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -637,7 +637,7 @@ contains
637637
E_L = E_L + (tau_e_L(i)*tau_e_L(i))/(4._wp*G_L)
638638
E_R = E_R + (tau_e_R(i)*tau_e_R(i))/(4._wp*G_R)
639639
! Additional terms in 2D and 3D
640-
if ((i == 2) .or. (i == 4) .or. (i == 5)) then
640+
if (any(strxb - 1 + i == shear_indices)) then
641641
E_L = E_L + (tau_e_L(i)*tau_e_L(i))/(4._wp*G_L)
642642
E_R = E_R + (tau_e_R(i)*tau_e_R(i))/(4._wp*G_R)
643643
end if
@@ -1084,17 +1084,26 @@ contains
10841084
flux_gsrc_rs${XYZ}$_vf(j, k, l, i) = flux_rs${XYZ}$_vf(j, k, l, i)
10851085
end do
10861086
! Recalculating the radial momentum geometric source flux
1087-
flux_gsrc_rs${XYZ}$_vf(j, k, l, contxe + dir_idx(1)) = &
1088-
(s_M*(rho_R*vel_R(dir_idx(1)) &
1089-
*vel_R(dir_idx(1))) &
1090-
- s_P*(rho_L*vel_L(dir_idx(1)) &
1091-
*vel_L(dir_idx(1))) &
1092-
+ s_M*s_P*(rho_L*vel_L(dir_idx(1)) &
1093-
- rho_R*vel_R(dir_idx(1)))) &
1094-
/(s_M - s_P)
1087+
flux_gsrc_rs${XYZ}$_vf(j, k, l, contxe + 2) = &
1088+
flux_rs${XYZ}$_vf(j, k, l, contxe + 2) &
1089+
- (s_M*pres_R - s_P*pres_L)/(s_M - s_P)
10951090
! Geometrical source of the void fraction(s) is zero
10961091
!$acc loop seq
10971092
do i = advxb, advxe
1093+
! flux_gsrc_rs${XYZ}$_vf(j, k, l, i) = 0._wp
1094+
flux_gsrc_rs${XYZ}$_vf(j, k, l, i) = flux_rs${XYZ}$_vf(j, k, l, i)
1095+
end do
1096+
end if
1097+
1098+
if (cyl_coord .and. hypoelasticity) then
1099+
! += tau_sigmasigma using HLL
1100+
flux_gsrc_rs${XYZ}$_vf(j, k, l, contxe + 2) = &
1101+
flux_gsrc_rs${XYZ}$_vf(j, k, l, contxe + 2) + &
1102+
(s_M*tau_e_R(4) - s_P*tau_e_L(4)) &
1103+
/(s_M - s_P)
1104+
1105+
!$acc loop seq
1106+
do i = strxb, strxe
10981107
flux_gsrc_rs${XYZ}$_vf(j, k, l, i) = flux_rs${XYZ}$_vf(j, k, l, i)
10991108
end do
11001109
end if

0 commit comments

Comments
 (0)