Skip to content

Commit ab8ae6c

Browse files
sbryngelsonsbryngelson
committed
revert: remove reduced-energy (E-tilde) scheme — moving to separate PR
1 parent 7e80db0 commit ab8ae6c

5 files changed

Lines changed: 30 additions & 158 deletions

File tree

‎src/common/m_variables_conversion.fpp‎

Lines changed: 2 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -104,12 +104,8 @@ contains
104104
if (mhd) then
105105
! MHD pressure: subtract magnetic pressure from total energy
106106
pres = (energy - dyn_p - pi_inf - qv - pres_mag)/gamma
107-
else if (model_eqns == 2 .and. (bubbles_euler .neqv. .true.)) then
108-
! Allaire et al. (JCP 2002): store reduced energy Ẽ = E - pi_inf_mix so that p = (Ẽ - KE - qv)/gamma avoids
109-
! cancellation when pi_inf_mix >> p.
110-
pres = (energy - dyn_p - qv)/gamma
111107
else if ((model_eqns /= 4) .and. (bubbles_euler .neqv. .true.)) then
112-
! model_eqns=1 or 3: physical E stored, p = (E - pi_inf - KE - qv)/gamma
108+
! Gamma/pi_inf model or five-equation model (Allaire et al. JCP 2002): p from mixture EOS
113109
pres = (energy - dyn_p - pi_inf - qv)/gamma
114110
else if ((model_eqns /= 4) .and. bubbles_euler) then
115111
! Bubble-augmented pressure with void fraction correction
@@ -923,11 +919,8 @@ contains
923919
! MHD energy includes magnetic pressure contribution
924920
q_cons_vf(eqn_idx%E)%sf(j, k, l) = gamma*q_prim_vf(eqn_idx%E)%sf(j, k, &
925921
& l) + dyn_pres + pres_mag + pi_inf + qv
926-
else if (model_eqns == 2 .and. (bubbles_euler .neqv. .true.)) then
927-
! Store reduced energy Ẽ = gamma*p + KE + qv (no pi_inf) for numerical stability.
928-
q_cons_vf(eqn_idx%E)%sf(j, k, l) = gamma*q_prim_vf(eqn_idx%E)%sf(j, k, l) + dyn_pres + qv
929922
else if ((model_eqns /= 4) .and. (bubbles_euler .neqv. .true.)) then
930-
! model_eqns=1 or 3: store physical E = gamma*p + pi_inf + KE + qv
923+
! Five-equation model (Allaire et al. JCP 2002): E = Gamma*p + 0.5*rho*|u|^2 + pi_inf + qv
931924
q_cons_vf(eqn_idx%E)%sf(j, k, l) = gamma*q_prim_vf(eqn_idx%E)%sf(j, k, l) + dyn_pres + pi_inf + qv
932925
else if ((model_eqns /= 4) .and. (bubbles_euler)) then
933926
! Bubble-augmented energy with void fraction correction

‎src/simulation/m_cbc.fpp‎

Lines changed: 3 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -658,11 +658,10 @@ contains
658658
gamma = 1.0_wp/(Cp/Cv - 1.0_wp)
659659
end if
660660
else
661-
! Reduced energy Ẽ = E - pi_inf; add pi_inf back for physical H.
662-
E = gamma*pres + 5.e-1_wp*rho*vel_K_sum
661+
E = gamma*pres + pi_inf + 5.e-1_wp*rho*vel_K_sum
663662
end if
664663
665-
H = (E + pi_inf + pres)/rho
664+
H = (E + pres)/rho
666665
667666
! Compute mixture sound speed
668667
call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, adv_local, vel_K_sum, 0._wp, c, qv)
@@ -891,7 +890,7 @@ contains
891890
end do
892891
else
893892
flux_rs${XYZ}$_vf_l(-1, k, r, eqn_idx%E) = flux_rs${XYZ}$_vf_l(0, k, r, &
894-
& eqn_idx%E) + ds(0)*(pres*dgamma_dt + gamma*dpres_dt + dqv_dt &
893+
& eqn_idx%E) + ds(0)*(pres*dgamma_dt + gamma*dpres_dt + dpi_inf_dt + dqv_dt &
895894
& + rho*vel_dv_dt_sum + 5.e-1_wp*drho_dt*vel_K_sum)
896895
end if
897896

‎src/simulation/m_rhs.fpp‎

Lines changed: 2 additions & 116 deletions
Original file line numberDiff line numberDiff line change
@@ -1118,9 +1118,9 @@ contains
11181118
type(vector_field), intent(in) :: flux_src_n_vf_arg
11191119
! CORRECTED DECLARATION FOR Kterm_arg:
11201120
real(wp), allocatable, dimension(:,:,:), intent(in) :: Kterm_arg
1121-
integer :: j_adv, k_idx, l_idx, q_idx, i_fluid
1121+
integer :: j_adv, k_idx, l_idx, q_idx
11221122
real(wp) :: local_inv_ds, local_term_coeff, local_flux1, local_flux2
1123-
real(wp) :: local_q_cons_val, local_k_term_val, pi_inf_src
1123+
real(wp) :: local_q_cons_val, local_k_term_val
11241124
logical :: use_standard_riemann
11251125
11261126
select case (current_idir)
@@ -1191,44 +1191,6 @@ contains
11911191
$:END_GPU_PARALLEL_LOOP()
11921192
end if
11931193
end if
1194-
! Ẽ non-conservative source for model_eqns=2: S_Ẽ = -\Sigma_i \pi_inf_i * rhs(\alpha_i)
1195-
if (model_eqns == 2 .and. (bubbles_euler .neqv. .true.) .and. .not. mhd) then
1196-
if (use_standard_riemann) then
1197-
$:GPU_PARALLEL_LOOP(collapse=3, private='[k_idx, l_idx, q_idx, local_inv_ds, local_term_coeff, &
1198-
& pi_inf_src, i_fluid]')
1199-
do q_idx = 0, p; do l_idx = 0, n; do k_idx = 0, m
1200-
local_inv_ds = 1._wp/dx(k_idx)
1201-
local_term_coeff = q_prim_vf_arg%vf(eqn_idx%cont%end + 1)%sf(k_idx, l_idx, q_idx)
1202-
pi_inf_src = 0._wp
1203-
do i_fluid = 1, num_fluids
1204-
pi_inf_src = pi_inf_src + pi_infs(i_fluid)*(flux_src_n_vf_arg%vf(eqn_idx%adv%beg + i_fluid - 1) &
1205-
& %sf(k_idx - 1, l_idx, &
1206-
& q_idx) - flux_src_n_vf_arg%vf(eqn_idx%adv%beg + i_fluid - 1) &
1207-
& %sf(k_idx, l_idx, q_idx))
1208-
end do
1209-
rhs_vf_arg(eqn_idx%E)%sf(k_idx, l_idx, q_idx) = rhs_vf_arg(eqn_idx%E)%sf(k_idx, l_idx, &
1210-
& q_idx) - local_inv_ds*local_term_coeff*pi_inf_src
1211-
end do; end do; end do
1212-
$:END_GPU_PARALLEL_LOOP()
1213-
else
1214-
$:GPU_PARALLEL_LOOP(collapse=3, private='[k_idx, l_idx, q_idx, local_inv_ds, pi_inf_src, i_fluid]')
1215-
do q_idx = 0, p; do l_idx = 0, n; do k_idx = 0, m
1216-
local_inv_ds = 1._wp/dx(k_idx)
1217-
pi_inf_src = 0._wp
1218-
do i_fluid = 1, num_fluids
1219-
pi_inf_src = pi_inf_src + pi_infs(i_fluid)*q_cons_vf_arg%vf(eqn_idx%adv%beg + i_fluid - 1) &
1220-
& %sf(k_idx, l_idx, &
1221-
& q_idx)*(flux_src_n_vf_arg%vf(eqn_idx%adv%beg + i_fluid - 1) &
1222-
& %sf(k_idx, l_idx, &
1223-
& q_idx) - flux_src_n_vf_arg%vf(eqn_idx%adv%beg + i_fluid - 1) &
1224-
& %sf(k_idx - 1, l_idx, q_idx))
1225-
end do
1226-
rhs_vf_arg(eqn_idx%E)%sf(k_idx, l_idx, q_idx) = rhs_vf_arg(eqn_idx%E)%sf(k_idx, l_idx, &
1227-
& q_idx) - local_inv_ds*pi_inf_src
1228-
end do; end do; end do
1229-
$:END_GPU_PARALLEL_LOOP()
1230-
end if
1231-
end if
12321194
case (2)
12331195
! y-direction: loops q_idx (x), k_idx (y), l_idx (z); sf(q_idx, k_idx, l_idx); dy(k_idx); Kterm(q_idx,k_idx,l_idx)
12341196
use_standard_riemann = (riemann_solver == 1 .or. riemann_solver == 4)
@@ -1305,44 +1267,6 @@ contains
13051267
$:END_GPU_PARALLEL_LOOP()
13061268
end if
13071269
end if
1308-
! Ẽ non-conservative source for model_eqns=2: S_Ẽ = -\Sigma_i \pi_inf_i * rhs(\alpha_i)
1309-
if (model_eqns == 2 .and. (bubbles_euler .neqv. .true.) .and. .not. mhd) then
1310-
if (use_standard_riemann) then
1311-
$:GPU_PARALLEL_LOOP(collapse=3, private='[k_idx, l_idx, q_idx, local_inv_ds, local_term_coeff, &
1312-
& pi_inf_src, i_fluid]')
1313-
do l_idx = 0, p; do k_idx = 0, n; do q_idx = 0, m
1314-
local_inv_ds = 1._wp/dy(k_idx)
1315-
local_term_coeff = q_prim_vf_arg%vf(eqn_idx%cont%end + 2)%sf(q_idx, k_idx, l_idx)
1316-
pi_inf_src = 0._wp
1317-
do i_fluid = 1, num_fluids
1318-
pi_inf_src = pi_inf_src + pi_infs(i_fluid)*(flux_src_n_vf_arg%vf(eqn_idx%adv%beg + i_fluid - 1) &
1319-
& %sf(q_idx, k_idx - 1, &
1320-
& l_idx) - flux_src_n_vf_arg%vf(eqn_idx%adv%beg + i_fluid - 1) &
1321-
& %sf(q_idx, k_idx, l_idx))
1322-
end do
1323-
rhs_vf_arg(eqn_idx%E)%sf(q_idx, k_idx, l_idx) = rhs_vf_arg(eqn_idx%E)%sf(q_idx, k_idx, &
1324-
& l_idx) - local_inv_ds*local_term_coeff*pi_inf_src
1325-
end do; end do; end do
1326-
$:END_GPU_PARALLEL_LOOP()
1327-
else
1328-
$:GPU_PARALLEL_LOOP(collapse=3, private='[k_idx, l_idx, q_idx, local_inv_ds, pi_inf_src, i_fluid]')
1329-
do l_idx = 0, p; do k_idx = 0, n; do q_idx = 0, m
1330-
local_inv_ds = 1._wp/dy(k_idx)
1331-
pi_inf_src = 0._wp
1332-
do i_fluid = 1, num_fluids
1333-
pi_inf_src = pi_inf_src + pi_infs(i_fluid)*q_cons_vf_arg%vf(eqn_idx%adv%beg + i_fluid - 1) &
1334-
& %sf(q_idx, k_idx, &
1335-
& l_idx)*(flux_src_n_vf_arg%vf(eqn_idx%adv%beg + i_fluid - 1) &
1336-
& %sf(q_idx, k_idx, &
1337-
& l_idx) - flux_src_n_vf_arg%vf(eqn_idx%adv%beg + i_fluid - 1) &
1338-
& %sf(q_idx, k_idx - 1, l_idx))
1339-
end do
1340-
rhs_vf_arg(eqn_idx%E)%sf(q_idx, k_idx, l_idx) = rhs_vf_arg(eqn_idx%E)%sf(q_idx, k_idx, &
1341-
& l_idx) - local_inv_ds*pi_inf_src
1342-
end do; end do; end do
1343-
$:END_GPU_PARALLEL_LOOP()
1344-
end if
1345-
end if
13461270
case (3)
13471271
! z-direction: loops l_idx (x), q_idx (y), k_idx (z); sf(l_idx, q_idx, k_idx); dz(k_idx); Kterm(l_idx,q_idx,k_idx)
13481272
if (grid_geometry == 3) then
@@ -1416,44 +1340,6 @@ contains
14161340
$:END_GPU_PARALLEL_LOOP()
14171341
end if
14181342
end if
1419-
! Ẽ non-conservative source for model_eqns=2: S_Ẽ = -\Sigma_i \pi_inf_i * rhs(\alpha_i)
1420-
if (model_eqns == 2 .and. (bubbles_euler .neqv. .true.) .and. .not. mhd) then
1421-
if (use_standard_riemann) then
1422-
$:GPU_PARALLEL_LOOP(collapse=3, private='[k_idx, l_idx, q_idx, local_inv_ds, local_term_coeff, &
1423-
& pi_inf_src, i_fluid]')
1424-
do k_idx = 0, p; do q_idx = 0, n; do l_idx = 0, m
1425-
local_inv_ds = 1._wp/dz(k_idx)
1426-
local_term_coeff = q_prim_vf_arg%vf(eqn_idx%cont%end + 3)%sf(l_idx, q_idx, k_idx)
1427-
pi_inf_src = 0._wp
1428-
do i_fluid = 1, num_fluids
1429-
pi_inf_src = pi_inf_src + pi_infs(i_fluid)*(flux_src_n_vf_arg%vf(eqn_idx%adv%beg + i_fluid - 1) &
1430-
& %sf(l_idx, q_idx, &
1431-
& k_idx - 1) - flux_src_n_vf_arg%vf(eqn_idx%adv%beg + i_fluid &
1432-
& - 1)%sf(l_idx, q_idx, k_idx))
1433-
end do
1434-
rhs_vf_arg(eqn_idx%E)%sf(l_idx, q_idx, k_idx) = rhs_vf_arg(eqn_idx%E)%sf(l_idx, q_idx, &
1435-
& k_idx) - local_inv_ds*local_term_coeff*pi_inf_src
1436-
end do; end do; end do
1437-
$:END_GPU_PARALLEL_LOOP()
1438-
else
1439-
$:GPU_PARALLEL_LOOP(collapse=3, private='[k_idx, l_idx, q_idx, local_inv_ds, pi_inf_src, i_fluid]')
1440-
do k_idx = 0, p; do q_idx = 0, n; do l_idx = 0, m
1441-
local_inv_ds = 1._wp/dz(k_idx)
1442-
pi_inf_src = 0._wp
1443-
do i_fluid = 1, num_fluids
1444-
pi_inf_src = pi_inf_src + pi_infs(i_fluid)*q_cons_vf_arg%vf(eqn_idx%adv%beg + i_fluid - 1) &
1445-
& %sf(l_idx, q_idx, &
1446-
& k_idx)*(flux_src_n_vf_arg%vf(eqn_idx%adv%beg + i_fluid - 1) &
1447-
& %sf(l_idx, q_idx, &
1448-
& k_idx) - flux_src_n_vf_arg%vf(eqn_idx%adv%beg + i_fluid - 1) &
1449-
& %sf(l_idx, q_idx, k_idx - 1))
1450-
end do
1451-
rhs_vf_arg(eqn_idx%E)%sf(l_idx, q_idx, k_idx) = rhs_vf_arg(eqn_idx%E)%sf(l_idx, q_idx, &
1452-
& k_idx) - local_inv_ds*pi_inf_src
1453-
end do; end do; end do
1454-
$:END_GPU_PARALLEL_LOOP()
1455-
end if
1456-
end if
14571343
end select
14581344
14591345
end subroutine s_add_directional_advection_source_terms

0 commit comments

Comments
 (0)