Skip to content

Commit a4feb61

Browse files
mrodrig6mrodrig6
committed
RMT now with HLL, need to build test cases
1 parent 49f2dcc commit a4feb61

File tree

3 files changed

+46
-55
lines changed

3 files changed

+46
-55
lines changed

‎src/common/m_derived_types.fpp‎

Lines changed: 0 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -168,7 +168,6 @@ module m_derived_types
168168
!! patch geometries. It is specified through its x-, y-, and z-components
169169
!! respectively.
170170

171-
type(ic_model_parameters) :: model !<
172171
real(wp) :: epsilon, beta !<
173172
!! The isentropic vortex parameters for the amplitude of the disturbance and
174173
!! domain of influence.
@@ -178,9 +177,6 @@ module m_derived_types
178177

179178
logical :: non_axis_sym
180179

181-
real(wp) :: a2, a3, a4, a5, a6, a7, a8, a9, a10, a11, a12 !<
182-
!! The parameters needed for the spherical harmonic patch
183-
184180
real(wp), dimension(3) :: normal !<
185181
!! Normal vector indicating the orientation of the patch. It is specified
186182
!! through its x-, y- and z-components, respectively.

‎src/pre_process/include/3dHardcodedIC.fpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -62,7 +62,7 @@
6262

6363
case (302) ! (3D lung geometry in X direction - axisym, with smoothing)
6464
lam = 200.e-06_wp
65-
amp = patch_icpp(patch_id)%a2
65+
amp = patch_icpp(patch_id)%a(2)
6666
h = 0.125_wp*amp
6767

6868
intH = amp/2._wp*(sin(2._wp*pi*y_cc(j)/lam + pi/2._wp) + sin(2._wp*pi*z_cc(k)/lam + pi/2._wp))

‎src/simulation/m_riemann_solvers.fpp‎

Lines changed: 45 additions & 50 deletions
Original file line numberDiff line numberDiff line change
@@ -549,7 +549,7 @@ contains
549549
tau_e_R(i) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, strxb - 1 + i)
550550
! Elastic contribution to energy if G large enough
551551
!TODO take out if statement if stable without
552-
if ((G_L > 1000) .and. (G_R > 1000)) then
552+
if ((G_L > verysmall) .and. (G_R > verysmall)) then
553553
E_L = E_L + (tau_e_L(i)*tau_e_L(i))/(4._wp*G_L)
554554
E_R = E_R + (tau_e_R(i)*tau_e_R(i))/(4._wp*G_R)
555555
! Additional terms in 2D and 3D
@@ -561,37 +561,31 @@ contains
561561
end do
562562
end if
563563

564-
! elastic energy update
565-
!if ( hyperelasticity ) then
566-
! G_L = 0._wp
567-
! G_R = 0._wp
568-
!
569-
! !$acc loop seq
570-
! do i = 1, num_fluids
571-
! G_L = G_L + alpha_L(i)*Gs(i)
572-
! G_R = G_R + alpha_R(i)*Gs(i)
573-
! end do
574-
! ! Elastic contribution to energy if G large enough
575-
! if ((G_L > 1e-3_wp) .and. (G_R > 1e-3_wp)) then
576-
! E_L = E_L + G_L*qL_prim_rs${XYZ}$_vf(j, k, l, xiend + 1)
577-
! E_R = E_R + G_R*qR_prim_rs${XYZ}$_vf(j + 1, k, l, xiend + 1)
578-
! !$acc loop seq
579-
! do i = 1, b_size-1
580-
! tau_e_L(i) = G_L*qL_prim_rs${XYZ}$_vf(j, k, l, strxb - 1 + i)
581-
! tau_e_R(i) = G_R*qR_prim_rs${XYZ}$_vf(j + 1, k, l, strxb - 1 + i)
582-
! end do
583-
! !$acc loop seq
584-
! do i = 1, b_size-1
585-
! tau_e_L(i) = 0_wp
586-
! tau_e_R(i) = 0_wp
587-
! end do
588-
! !$acc loop seq
589-
! do i = 1, num_dims
590-
! xi_field_L(i) = qL_prim_rs${XYZ}$_vf(j, k, l, xibeg - 1 + i)
591-
! xi_field_R(i) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, xibeg - 1 + i)
592-
! end do
593-
! end if
594-
!end if
564+
! ENERGY ADJUSTMENTS FOR HYPERELASTIC ENERGY
565+
if (hyperelasticity) then
566+
!$acc loop seq
567+
do i = 1, num_dims
568+
xi_field_L(i) = qL_prim_rs${XYZ}$_vf(j, k, l, xibeg - 1 + i)
569+
xi_field_R(i) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, xibeg - 1 + i)
570+
end do
571+
G_L = 0_wp; G_R = 0_wp;
572+
!$acc loop seq
573+
do i = 1, num_fluids
574+
! Mixture left and right shear modulus
575+
G_L = G_L + alpha_L(i)*Gs(i)
576+
G_R = G_R + alpha_R(i)*Gs(i)
577+
end do
578+
! Elastic contribution to energy if G large enough
579+
if (G_L > verysmall .and. G_R > verysmall) then
580+
E_L = E_L + G_L*qL_prim_rs${XYZ}$_vf(j, k, l, xiend + 1)
581+
E_R = E_R + G_R*qR_prim_rs${XYZ}$_vf(j + 1, k, l, xiend + 1)
582+
end if
583+
!$acc loop seq
584+
do i = 1, b_size - 1
585+
tau_e_L(i) = qL_prim_rs${XYZ}$_vf(j, k, l, strxb - 1 + i)
586+
tau_e_R(i) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, strxb - 1 + i)
587+
end do
588+
end if
595589

596590
! Enthalpy with elastic energy
597591
H_L = (E_L + pres_L)/rho_L
@@ -648,6 +642,7 @@ contains
648642
(s_R - vel_R(dir_idx(1)))) &
649643
/(rho_L*(s_L - vel_L(dir_idx(1))) - &
650644
rho_R*(s_R - vel_R(dir_idx(1))))
645+
651646
elseif (wave_speeds == 2) then
652647
pres_SL = 5e-1_wp*(pres_L + pres_R + rho_avg*c_avg* &
653648
(vel_L(dir_idx(1)) - &
@@ -705,7 +700,7 @@ contains
705700
- rho_R*vel_R(dir_idx(i)))) &
706701
/(s_M - s_P)
707702
end do
708-
else if (hypoelasticity) then
703+
else if (elasticity) then
709704
!$acc loop seq
710705
do i = 1, num_dims
711706
flux_rs${XYZ}$_vf(j, k, l, contxe + dir_idx(i)) = &
@@ -737,14 +732,14 @@ contains
737732
end do
738733
end if
739734

740-
! Energy
735+
! ENERGY UPDATE
741736
if (bubbles) then
742737
flux_rs${XYZ}$_vf(j, k, l, E_idx) = &
743738
(s_M*vel_R(dir_idx(1))*(E_R + pres_R - ptilde_R) &
744739
- s_P*vel_L(dir_idx(1))*(E_L + pres_L - ptilde_L) &
745740
+ s_M*s_P*(E_L - E_R)) &
746741
/(s_M - s_P)
747-
else if (hypoelasticity) then
742+
else if (elasticity) then
748743
!TODO: simplify this so it's not split into 3
749744
if (num_dims == 1) then
750745
flux_rs${XYZ}$_vf(j, k, l, E_idx) = &
@@ -785,9 +780,9 @@ contains
785780
/(s_M - s_P)
786781
end if
787782
788-
! Elastic Stresses
783+
! ELASTIC STRESSES FLUX.
789784
if (hypoelasticity) then
790-
do i = 1, strxe - strxb + 1 !TODO: this indexing may be slow
785+
do i = 1, strxe - strxb + 1
791786
flux_rs${XYZ}$_vf(j, k, l, strxb - 1 + i) = &
792787
(s_M*(rho_R*vel_R(dir_idx(1)) &
793788
*tau_e_R(i)) &
@@ -799,7 +794,19 @@ contains
799794
end do
800795
end if
801796
802-
! Advection
797+
! REFERENCE MAP FLUX.
798+
if (hyperelasticity) then
799+
do i = 1, num_dims
800+
flux_rs${XYZ}$_vf(j, k, l, xibeg - 1 + i) = &
801+
(s_M*rho_R*vel_R(dir_idx(1))*xi_field_R(i) &
802+
- s_P*rho_L*vel_L(dir_idx(1))*xi_field_L(i) &
803+
+ s_M*s_P*(rho_L*xi_field_L(i) &
804+
- rho_R*xi_field_R(i))) &
805+
/(s_M - s_P)
806+
end do
807+
end if
808+
809+
! ADVECTION FLUX.
803810
!$acc loop seq
804811
do i = advxb, advxe
805812
flux_rs${XYZ}$_vf(j, k, l, i) = &
@@ -812,18 +819,6 @@ contains
812819
/(s_M - s_P)
813820
end do
814821
815-
! Xi field
816-
!if ( hyperelasticity ) then
817-
! do i = 1, num_dims
818-
! flux_rs${XYZ}$_vf(j, k, l, xibeg - 1 + i) = &
819-
! (s_M*rho_R*vel_R(dir_idx(1))*xi_field_R(i) &
820-
! - s_P*rho_L*vel_L(dir_idx(1))*xi_field_L(i) &
821-
! + s_M*s_P*(rho_L*xi_field_L(i) &
822-
! - rho_R*xi_field_R(i))) &
823-
! /(s_M - s_P)
824-
! end do
825-
!end if
826-
827822
! Div(U)?
828823
!$acc loop seq
829824
do i = 1, num_dims

0 commit comments

Comments
 (0)