diff --git a/docs/documentation/case.md b/docs/documentation/case.md index ec82b5bbd..6b7274f24 100644 --- a/docs/documentation/case.md +++ b/docs/documentation/case.md @@ -134,8 +134,6 @@ feature, detecting GPU pointers and performing RDMA accordingly. | `m` | Integer | Number of grid cells in the $x$-coordinate direction | | `n` | Integer | Number of grid cells in the $y$-coordinate direction | | `p` | Integer | Number of grid cells in the $z$-coordinate direction | -| `pref` | Real | Reference pressure | -| `rhoref` | Real | Reference density | The parameters define the boundaries of the spatial and temporal domains, and their discretization that are used in simulation. @@ -489,7 +487,7 @@ See @ref equations "Equations" for the mathematical models these parameters cont | `bc_[x,y,z]%%beg[end]` | Integer | Beginning [ending] boundary condition in the $[x,y,z]$-direction (negative integer, see table [Boundary Conditions](#boundary-conditions)) | | `bc_[x,y,z]%%vb[1,2,3]`‡ | Real | Velocity in the (x,1), (y, 2), (z,3) direction applied to `bc_[x,y,z]%%beg` | | `bc_[x,y,z]%%ve[1,2,3]`‡ | Real | Velocity in the (x,1), (y, 2), (z,3) direction applied to `bc_[x,y,z]%%end` | -| `model_eqns` | Integer | Multicomponent model: [1] \f$\Gamma/\Pi_\infty\f$; [2] 5-equation; [3] 6-equation; [4] 4-equation | +| `model_eqns` | Integer | Multicomponent model: [1] \f$\Gamma/\Pi_\infty\f$; [2] 5-equation; [3] 6-equation | | `alt_soundspeed` * | Logical | Alternate sound speed and \f$K \nabla \cdot u\f$ for 5-equation model | | `adv_n` | Logical | Solving directly for the number density (in the method of classes) and compute void fraction from the number density | | `mpp_lim` | Logical | Mixture physical parameters limits | @@ -873,8 +871,6 @@ Details of the transducer acoustic source model can be found in \cite Maeda17. | `Ca` | Real | Cavitation number | | `Web` | Real | Weber number | | `Re_inv` | Real | Inverse Reynolds number | -| `pref` | Real | Reference pressure for bubble models | -| `rhoref` | Real | Reference density for bubble models | | `fluid_rho` | Real | Reference fluid density | | `bub_pp%%R0ref`*†‡ | Real | Reference bubble radius | | `bub_pp%%p0ref`*†‡ | Real | Reference pressure | diff --git a/docs/documentation/contributing.md b/docs/documentation/contributing.md index ad85b864a..0b261400e 100644 --- a/docs/documentation/contributing.md +++ b/docs/documentation/contributing.md @@ -198,7 +198,7 @@ Both human reviewers and AI code reviewers reference this section. ### Physics and Model Consistency -- **Pressure formula** must match `model_eqns` value. Model 2/3 (multi-fluid), model 4 (bubbles), MHD, and hypoelastic each use different EOS formulations. Wrong formula = wrong physics. +- **Pressure formula** must match `model_eqns` value. Model 2/3 (multi-fluid), MHD, and hypoelastic each use different EOS formulations. Wrong formula = wrong physics. - **Conservative-primitive conversion:** Density recovery, kinetic energy, and pressure each have model-specific paths. Verify the correct branch is taken. - **Volume fractions** must sum to 1. `alpha_rho_K` must be non-negative. Species mass fractions should be clipped to [0,1]. - **Boundary conditions:** Periodic BCs must match at both ends (`bc_x%%beg` and `bc_x%%end`). Cylindrical coordinates have special requirements (`bc_y%%beg = -14` for axis in 3D). diff --git a/docs/documentation/equations.md b/docs/documentation/equations.md index 0c120da36..0d64970e6 100644 --- a/docs/documentation/equations.md +++ b/docs/documentation/equations.md @@ -23,7 +23,7 @@ where: - \f$\mathbf{h}(\mathbf{q})\,\nabla \cdot \mathbf{u}\f$ contains non-conservative terms (volume fraction advection), - \f$\mathbf{s}(\mathbf{q})\f$ is the source vector (bubbles, body forces, chemistry, etc.). -The parameter `model_eqns` (1, 2, 3, or 4) selects the governing equation set. +The parameter `model_eqns` (1, 2, or 3) selects the governing equation set. **Key source files:** `src/simulation/m_rhs.fpp` (RHS evaluation), `src/common/m_variables_conversion.fpp` (EOS and variable conversion). @@ -299,7 +299,6 @@ See Section 8 (Phase Change) below for details. ### 2.3 Other Model Variants - `model_eqns = 1`: **Gamma/pi_inf model** — simplified single-fluid formulation using mixture \f$\gamma\f$ and \f$\pi_\infty\f$ directly without tracking individual volume fractions (\cite Johnsen08). -- `model_eqns = 4`: **Four-equation model** — reduced model from the six-equation system after full pressure-temperature equilibrium relaxation (Tait-like compressible liquid). --- diff --git a/src/common/include/2dHardcodedIC.fpp b/src/common/include/2dHardcodedIC.fpp index 360c49b71..aec33edb3 100644 --- a/src/common/include/2dHardcodedIC.fpp +++ b/src/common/include/2dHardcodedIC.fpp @@ -112,13 +112,13 @@ q_prim_vf(eqn_idx%adv%end)%sf(i, j, 0) = 1._wp - alph q_prim_vf(eqn_idx%cont%beg)%sf(i, j, 0) = alph*rhoH q_prim_vf(eqn_idx%cont%end)%sf(i, j, 0) = (1._wp - alph)*rhoL - q_prim_vf(eqn_idx%E)%sf(i, j, 0) = pref + rhoH*9.81_wp*(1.2_wp - y_cc(j)) + q_prim_vf(eqn_idx%E)%sf(i, j, 0) = pRef + rhoH*9.81_wp*(1.2_wp - y_cc(j)) else q_prim_vf(eqn_idx%adv%beg)%sf(i, j, 0) = alph q_prim_vf(eqn_idx%adv%end)%sf(i, j, 0) = 1._wp - alph q_prim_vf(eqn_idx%cont%beg)%sf(i, j, 0) = alph*rhoH q_prim_vf(eqn_idx%cont%end)%sf(i, j, 0) = (1._wp - alph)*rhoL - pInt = pref + rhoH*9.81_wp*(1.2_wp - intH) + pInt = pRef + rhoH*9.81_wp*(1.2_wp - intH) q_prim_vf(eqn_idx%E)%sf(i, j, 0) = pInt + rhoL*9.81_wp*(intH - y_cc(j)) end if case (205) ! 2D lung wave interaction problem diff --git a/src/common/include/3dHardcodedIC.fpp b/src/common/include/3dHardcodedIC.fpp index 76835a1aa..fd8101e94 100644 --- a/src/common/include/3dHardcodedIC.fpp +++ b/src/common/include/3dHardcodedIC.fpp @@ -155,13 +155,13 @@ q_prim_vf(eqn_idx%adv%end)%sf(i, j, k) = 1._wp - alph q_prim_vf(eqn_idx%cont%beg)%sf(i, j, k) = alph*rhoH q_prim_vf(eqn_idx%cont%end)%sf(i, j, k) = (1._wp - alph)*rhoL - q_prim_vf(eqn_idx%E)%sf(i, j, k) = pref + rhoH*9.81_wp*(1.2_wp - y_cc(j)) + q_prim_vf(eqn_idx%E)%sf(i, j, k) = pRef + rhoH*9.81_wp*(1.2_wp - y_cc(j)) else q_prim_vf(eqn_idx%adv%beg)%sf(i, j, k) = alph q_prim_vf(eqn_idx%adv%end)%sf(i, j, k) = 1._wp - alph q_prim_vf(eqn_idx%cont%beg)%sf(i, j, k) = alph*rhoH q_prim_vf(eqn_idx%cont%end)%sf(i, j, k) = (1._wp - alph)*rhoL - pInt = pref + rhoH*9.81_wp*(1.2_wp - intH) + pInt = pRef + rhoH*9.81_wp*(1.2_wp - intH) q_prim_vf(eqn_idx%E)%sf(i, j, k) = pInt + rhoL*9.81_wp*(intH - y_cc(j)) end if case (301) ! (3D lung geometry in X direction, |sin(*)+sin(*)|) diff --git a/src/common/m_global_parameters_common.fpp b/src/common/m_global_parameters_common.fpp index 76184fe13..28b4f00d9 100644 --- a/src/common/m_global_parameters_common.fpp +++ b/src/common/m_global_parameters_common.fpp @@ -15,8 +15,8 @@ module m_global_parameters_common use m_derived_types use m_thermochem, only: num_species - use m_constants, only: model_eqns_gamma_law, model_eqns_5eq, model_eqns_6eq, model_eqns_4eq, recon_type_weno, & - & recon_type_muscl, name_len, dflt_int, dflt_real + use m_constants, only: model_eqns_gamma_law, model_eqns_5eq, model_eqns_6eq, recon_type_weno, recon_type_muscl, name_len, & + & dflt_int, dflt_real implicit none @@ -190,26 +190,6 @@ contains eqn_idx%int_en%beg = eqn_idx%adv%end + 1 eqn_idx%int_en%end = eqn_idx%adv%end + num_fluids sys_size = eqn_idx%int_en%end - else if (model_eqns == model_eqns_4eq) then - ! 4-equation model with subgrid bubbles - eqn_idx%cont%beg = 1 - eqn_idx%cont%end = 1 - eqn_idx%mom%beg = eqn_idx%cont%end + 1 - eqn_idx%mom%end = eqn_idx%cont%end + num_vels - eqn_idx%E = eqn_idx%mom%end + 1 - eqn_idx%adv%beg = eqn_idx%E + 1 - eqn_idx%adv%end = eqn_idx%adv%beg - eqn_idx%alf = eqn_idx%adv%end - sys_size = eqn_idx%adv%end - - if (bubbles_euler) then - eqn_idx%bub%beg = sys_size + 1 - eqn_idx%bub%end = sys_size + 2*nb_in - if (.not. polytropic) then - eqn_idx%bub%end = sys_size + 4*nb_in - end if - sys_size = eqn_idx%bub%end - end if end if if (model_eqns == model_eqns_5eq .or. model_eqns == model_eqns_6eq) then @@ -374,8 +354,6 @@ contains riemann_solver = dflt_int ! Tait EOS - rhoref = dflt_real - pref = dflt_real ! Bubble modeling flags and parameters R0ref = dflt_real diff --git a/src/common/m_variables_conversion.fpp b/src/common/m_variables_conversion.fpp index f9fb8c73a..78d59323d 100644 --- a/src/common/m_variables_conversion.fpp +++ b/src/common/m_variables_conversion.fpp @@ -14,7 +14,7 @@ module m_variables_conversion use m_helper_basic use m_helper use m_constants, only: riemann_solver_hll, riemann_solver_hlld, model_eqns_gamma_law, model_eqns_5eq, model_eqns_6eq, & - & model_eqns_4eq, avg_state_roe + & avg_state_roe use m_thermochem, only: num_species, get_temperature, get_pressure, gas_constant, get_mixture_molecular_weight, & & get_mixture_energy_mass @@ -94,15 +94,12 @@ contains if (mhd) then ! MHD pressure: subtract magnetic pressure from total energy pres = (energy - dyn_p - pi_inf - qv - pres_mag)/gamma - else if ((model_eqns /= model_eqns_4eq) .and. (bubbles_euler .neqv. .true.)) then + else if (bubbles_euler .neqv. .true.) then ! Gamma/pi_inf model or five-equation model (Allaire et al. JCP 2002): p from mixture EOS pres = (energy - dyn_p - pi_inf - qv)/gamma - else if ((model_eqns /= model_eqns_4eq) .and. bubbles_euler) then + else ! Bubble-augmented pressure with void fraction correction pres = ((energy - dyn_p)/(1._wp - alf) - pi_inf - qv)/gamma - else - ! Four-equation model (Kapila et al. PoF 2001): Tait EOS inversion - pres = (pref + pi_inf)*(energy/(rhoref*(1 - alf)))**(1/gamma + 1) - pi_inf end if if (hypoelasticity .and. present(G)) then @@ -469,26 +466,24 @@ contains call s_compute_species_fraction(qK_cons_vf, j, k, l, alpha_rho_K, alpha_K) - if (model_eqns /= model_eqns_4eq) then #ifdef MFC_GPU - ! Device regions call the device-compiled scalar kernel directly. - if (hypoelasticity) then - call s_convert_species_to_mixture_variables_kernel(rho_K, gamma_K, pi_inf_K, qv_K, alpha_K, & - & alpha_rho_K, Re_K, G_K, Gs_vc) - else - call s_convert_species_to_mixture_variables_kernel(rho_K, gamma_K, pi_inf_K, qv_K, alpha_K, & - & alpha_rho_K, Re_K) - end if + ! Device regions call the device-compiled scalar kernel directly. + if (hypoelasticity) then + call s_convert_species_to_mixture_variables_kernel(rho_K, gamma_K, pi_inf_K, qv_K, alpha_K, alpha_rho_K, & + & Re_K, G_K, Gs_vc) + else + call s_convert_species_to_mixture_variables_kernel(rho_K, gamma_K, pi_inf_K, qv_K, alpha_K, alpha_rho_K, & + & Re_K) + end if #else - ! Host execution uses the wrapper, which also stores requested diagnostics. - if (hypoelasticity) then - call s_convert_to_mixture_variables(qK_cons_vf, j, k, l, rho_K, gamma_K, pi_inf_K, qv_K, Re_K, G_K, & - & fluid_pp(:)%G) - else - call s_convert_to_mixture_variables(qK_cons_vf, j, k, l, rho_K, gamma_K, pi_inf_K, qv_K) - end if -#endif + ! Host execution uses the wrapper, which also stores requested diagnostics. + if (hypoelasticity) then + call s_convert_to_mixture_variables(qK_cons_vf, j, k, l, rho_K, gamma_K, pi_inf_K, qv_K, Re_K, G_K, & + & fluid_pp(:)%G) + else + call s_convert_to_mixture_variables(qK_cons_vf, j, k, l, rho_K, gamma_K, pi_inf_K, qv_K) end if +#endif ! Relativistic MHD primitive variable recovery, Mignone & Bodo A&A (2006) if (relativity) then @@ -597,13 +592,8 @@ contains ! Recover velocity from momentum: u = rho*u / rho, and accumulate dynamic pressure 0.5*rho*|u|^2 $:GPU_LOOP(parallelism='[seq]') do i = eqn_idx%mom%beg, eqn_idx%mom%end - if (model_eqns /= model_eqns_4eq) then - qK_prim_vf(i)%sf(j, k, l) = qK_cons_vf(i)%sf(j, k, l)/rho_K - dyn_pres_K = dyn_pres_K + 5.e-1_wp*qK_cons_vf(i)%sf(j, k, l)*qK_prim_vf(i)%sf(j, k, l) - else - ! Four-equation model (Kapila et al. PoF 2001): divide by total density q_cons(1) - qK_prim_vf(i)%sf(j, k, l) = qK_cons_vf(i)%sf(j, k, l)/qK_cons_vf(1)%sf(j, k, l) - end if + qK_prim_vf(i)%sf(j, k, l) = qK_cons_vf(i)%sf(j, k, l)/rho_K + dyn_pres_K = dyn_pres_K + 5.e-1_wp*qK_cons_vf(i)%sf(j, k, l)*qK_prim_vf(i)%sf(j, k, l) end do if (chemistry) then @@ -868,16 +858,13 @@ contains ! MHD energy includes magnetic pressure contribution q_cons_vf(eqn_idx%E)%sf(j, k, l) = gamma*q_prim_vf(eqn_idx%E)%sf(j, k, & & l) + dyn_pres + pres_mag + pi_inf + qv - else if ((model_eqns /= model_eqns_4eq) .and. (bubbles_euler .neqv. .true.)) then + else if (bubbles_euler .neqv. .true.) then ! Five-equation model (Allaire et al. JCP 2002): E = Gamma*p + 0.5*rho*|u|^2 + pi_inf + qv 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 - else if ((model_eqns /= model_eqns_4eq) .and. (bubbles_euler)) then + else ! Bubble-augmented energy with void fraction correction q_cons_vf(eqn_idx%E)%sf(j, k, l) = dyn_pres + (1._wp - q_prim_vf(eqn_idx%alf)%sf(j, k, & & l))*(gamma*q_prim_vf(eqn_idx%E)%sf(j, k, l) + pi_inf) - else - ! Four-equation model (Kapila et al. PoF 2001): Tait EOS, no conserved energy variable - q_cons_vf(eqn_idx%E)%sf(j, k, l) = 0._wp end if end if @@ -1215,7 +1202,7 @@ contains c = c + adv(q)*gs_min(q)*(pres + pi_infs(q)/(gammas(q) + 1._wp)) end do c = c/rho - else if (((model_eqns == model_eqns_4eq) .or. (model_eqns == model_eqns_5eq .and. bubbles_euler))) then + else if (model_eqns == model_eqns_5eq .and. bubbles_euler) then ! Sound speed for bubble mixture to order O(\alpha) if (mpp_lim .and. (num_fluids > 1)) then diff --git a/src/post_process/m_global_parameters.fpp b/src/post_process/m_global_parameters.fpp index 1423c62fe..a5151f2df 100644 --- a/src/post_process/m_global_parameters.fpp +++ b/src/post_process/m_global_parameters.fpp @@ -386,40 +386,6 @@ contains end if end if - if (model_eqns == model_eqns_4eq .and. bubbles_euler) then - allocate (qbmm_idx%rs(nb), qbmm_idx%vs(nb)) - allocate (qbmm_idx%ps(nb), qbmm_idx%ms(nb)) - allocate (weight(nb), R0(nb)) - - do i = 1, nb - if (polytropic .neqv. .true.) then - fac = 4 - else - fac = 2 - end if - - qbmm_idx%rs(i) = eqn_idx%bub%beg + (i - 1)*fac - qbmm_idx%vs(i) = qbmm_idx%rs(i) + 1 - - if (polytropic .neqv. .true.) then - qbmm_idx%ps(i) = qbmm_idx%vs(i) + 1 - qbmm_idx%ms(i) = qbmm_idx%ps(i) + 1 - end if - end do - - if (nb == 1) then - weight(:) = 1._wp - R0(:) = 1._wp - else if (nb < 1) then - stop 'Invalid value of nb' - end if - - if (polytropic) then - rhoref = 1._wp - pref = 1._wp - end if - end if - if (output_partial_domain) then x_output_idx%beg = 0 x_output_idx%end = 0 diff --git a/src/post_process/m_start_up.fpp b/src/post_process/m_start_up.fpp index 66ed6cfb0..9c697c0ea 100644 --- a/src/post_process/m_start_up.fpp +++ b/src/post_process/m_start_up.fpp @@ -26,7 +26,7 @@ module m_start_up use m_checker use m_thermochem, only: num_species, species_names use m_finite_differences - use m_constants, only: model_eqns_gamma_law, model_eqns_5eq, model_eqns_6eq, model_eqns_4eq, format_silo + use m_constants, only: model_eqns_gamma_law, model_eqns_5eq, model_eqns_6eq, format_silo use m_chemistry #ifdef MFC_MPI @@ -235,14 +235,10 @@ contains call s_compute_finite_difference_coefficients(p, z_cc, fd%fd_coeff_z, buff_size, fd_number, fd_order, offset_z) end if - if ((model_eqns == model_eqns_5eq) .or. (model_eqns == model_eqns_6eq) .or. (model_eqns == model_eqns_4eq)) then + if ((model_eqns == model_eqns_5eq) .or. (model_eqns == model_eqns_6eq)) then do i = 1, num_fluids if (alpha_rho_wrt(i) .or. (cons_vars_wrt .or. prim_vars_wrt)) then - if (model_eqns /= model_eqns_4eq) then - write (varname, '(A,I0)') 'alpha_rho', i - else - write (varname, '(A,I0)') 'rho', i - end if + write (varname, '(A,I0)') 'alpha_rho', i call s_write_field(varname, t_step, q_cons_vf(i), x_beg, x_end, y_beg, y_end, z_beg, z_end) end if end do diff --git a/src/pre_process/m_assign_variables.fpp b/src/pre_process/m_assign_variables.fpp index 20340576a..bf17249c1 100644 --- a/src/pre_process/m_assign_variables.fpp +++ b/src/pre_process/m_assign_variables.fpp @@ -13,7 +13,7 @@ module m_assign_variables use m_variables_conversion use m_helper_basic use m_thermochem, only: num_species, gas_constant, get_mixture_molecular_weight - use m_constants, only: model_eqns_gamma_law, model_eqns_6eq, model_eqns_4eq + use m_constants, only: model_eqns_gamma_law, model_eqns_6eq implicit none @@ -222,7 +222,6 @@ contains ! Density, gamma, and liquid stiffness from current and smoothing patches real(wp) :: rho !< density real(wp) :: gamma - real(wp) :: lit_gamma !< specific heat ratio real(wp) :: pi_inf !< stiffness from SEOS real(wp) :: qv !< reference energy from SEOS real(wp) :: orig_rho @@ -274,20 +273,16 @@ contains end do end if - if (model_eqns /= model_eqns_4eq) then - do i = 1, eqn_idx%cont%end - q_prim_vf(i)%sf(j, k, l) = patch_icpp(patch_id)%alpha_rho(i) - end do - end if + do i = 1, eqn_idx%cont%end + q_prim_vf(i)%sf(j, k, l) = patch_icpp(patch_id)%alpha_rho(i) + end do call s_convert_to_mixture_variables(q_prim_vf, j, k, l, patch_icpp(patch_id)%rho, patch_icpp(patch_id)%gamma, & & patch_icpp(patch_id)%pi_inf, patch_icpp(patch_id)%qv) - if (model_eqns /= model_eqns_4eq) then - do i = 1, eqn_idx%cont%end - q_prim_vf(i)%sf(j, k, l) = patch_icpp(smooth_patch_id)%alpha_rho(i) - end do - end if + do i = 1, eqn_idx%cont%end + q_prim_vf(i)%sf(j, k, l) = patch_icpp(smooth_patch_id)%alpha_rho(i) + end do if (.not. igr .or. num_fluids > 1) then do i = eqn_idx%adv%beg, eqn_idx%adv%end @@ -393,21 +388,10 @@ contains end do end if - if (model_eqns /= model_eqns_4eq) then - ! mixture density is an input - do i = 1, eqn_idx%cont%end - q_prim_vf(i)%sf(j, k, l) = eta*patch_icpp(patch_id)%alpha_rho(i) + (1._wp - eta)*orig_prim_vf(i) - end do - else - ! get mixture density from pressure via Tait EOS - pi_inf = pi_infs(1) - gamma = gammas(1) - lit_gamma = gs_min(1) - - ! \rho = (( p_l + pi_inf)/( p_ref + pi_inf))**(1/little_gam) * rhoref(1-alf) - q_prim_vf(1)%sf(j, k, l) = (((q_prim_vf(eqn_idx%E)%sf(j, k, & - & l) + pi_inf)/(pref + pi_inf))**(1/lit_gamma))*rhoref*(1 - q_prim_vf(eqn_idx%alf)%sf(j, k, l)) - end if + ! mixture density is an input + do i = 1, eqn_idx%cont%end + q_prim_vf(i)%sf(j, k, l) = eta*patch_icpp(patch_id)%alpha_rho(i) + (1._wp - eta)*orig_prim_vf(i) + end do call s_convert_to_mixture_variables(q_prim_vf, j, k, l, rho, gamma, pi_inf, qv) diff --git a/src/pre_process/m_global_parameters.fpp b/src/pre_process/m_global_parameters.fpp index 6a76fc073..e392f467c 100644 --- a/src/pre_process/m_global_parameters.fpp +++ b/src/pre_process/m_global_parameters.fpp @@ -487,40 +487,6 @@ contains end if end if - if (model_eqns == model_eqns_4eq .and. bubbles_euler) then - allocate (qbmm_idx%rs(nb), qbmm_idx%vs(nb)) - allocate (qbmm_idx%ps(nb), qbmm_idx%ms(nb)) - allocate (weight(nb), R0(nb)) - - do i = 1, nb - if (.not. polytropic) then - fac = 4 - else - fac = 2 - end if - - qbmm_idx%rs(i) = eqn_idx%bub%beg + (i - 1)*fac - qbmm_idx%vs(i) = qbmm_idx%rs(i) + 1 - - if (.not. polytropic) then - qbmm_idx%ps(i) = qbmm_idx%vs(i) + 1 - qbmm_idx%ms(i) = qbmm_idx%ps(i) + 1 - end if - end do - - if (nb == 1) then - weight(:) = 1._wp - R0(:) = 1._wp - else if (nb < 1) then - stop 'Invalid value of nb' - end if - - if (polytropic) then - rhoref = 1._wp - pref = 1._wp - end if - end if - if (bubbles_lagrange) fd_number = max(1, fd_order/2) call s_configure_coordinate_bounds(recon_type, weno_polyn, muscl_polyn, igr_order, buff_size, idwint, idwbuff, viscous, & diff --git a/src/pre_process/m_icpp_patches.fpp b/src/pre_process/m_icpp_patches.fpp index 286addf56..4d2116085 100644 --- a/src/pre_process/m_icpp_patches.fpp +++ b/src/pre_process/m_icpp_patches.fpp @@ -16,7 +16,7 @@ module m_icpp_patches use m_model ! Subroutine(s) related to STL files use m_derived_types ! Definitions of the derived types use m_global_parameters - use m_constants, only: max_2d_fourier_modes, max_sph_harm_degree, small_radius, model_eqns_4eq + use m_constants, only: max_2d_fourier_modes, max_sph_harm_degree, small_radius use m_helper_basic use m_helper use m_mpi_common @@ -177,14 +177,10 @@ contains integer :: i, j, k ! Placeholders for the cell boundary values - real(wp) :: pi_inf, gamma, lit_gamma @:HardcodedDimensionsExtrusion() @:Hardcoded1DVariables() - pi_inf = pi_infs(1) - gamma = gammas(1) - lit_gamma = gs_min(1) j = 0 k = 0 @@ -602,16 +598,11 @@ contains integer, dimension(0:m,0:n,0:p), intent(inout) :: patch_id_fp #endif type(scalar_field), dimension(1:sys_size), intent(inout) :: q_prim_vf - integer :: i, j, k !< generic loop iterators - real(wp) :: pi_inf, gamma, lit_gamma !< Equation of state parameters + integer :: i, j, k !< generic loop iterators @:HardcodedDimensionsExtrusion() @:Hardcoded2DVariables() - pi_inf = pi_infs(1) - gamma = gammas(1) - lit_gamma = gs_min(1) - ! Transferring the rectangle's centroid and length information x_centroid = patch_icpp(patch_id)%x_centroid y_centroid = patch_icpp(patch_id)%y_centroid @@ -640,13 +631,6 @@ contains @:Hardcoded2D() end if - if ((q_prim_vf(1)%sf(i, j, 0) < 1.e-10) .and. (model_eqns == model_eqns_4eq)) then - ! zero density, reassign according to Tait EOS - q_prim_vf(1)%sf(i, j, 0) = (((q_prim_vf(eqn_idx%E)%sf(i, j, & - & 0) + pi_inf)/(pref + pi_inf))**(1._wp/lit_gamma))*rhoref*(1._wp & - & - q_prim_vf(eqn_idx%alf)%sf(i, j, 0)) - end if - ! Updating the patch identities bookkeeping variable if (1._wp - eta < sgm_eps) patch_id_fp(i, j, 0) = patch_id end if @@ -728,17 +712,12 @@ contains integer, dimension(0:m,0:n,0:p), intent(inout) :: patch_id_fp #endif type(scalar_field), dimension(1:sys_size), intent(inout) :: q_prim_vf - integer :: i, j, k !< generic loop iterators - real(wp) :: pi_inf, gamma, lit_gamma !< equation of state parameters - real(wp) :: L0, U0 !< Taylor Green Vortex parameters + integer :: i, j, k !< generic loop iterators + real(wp) :: L0, U0 !< Taylor Green Vortex parameters @:HardcodedDimensionsExtrusion() @:Hardcoded2DVariables() - pi_inf = pi_infs(1) - gamma = gammas(1) - lit_gamma = gs_min(1) - ! Transferring the patch's centroid and length information x_centroid = patch_icpp(patch_id)%x_centroid y_centroid = patch_icpp(patch_id)%y_centroid @@ -803,15 +782,10 @@ contains ! Generic loop iterators integer :: i, j, k ! Placeholders for the cell boundary values - real(wp) :: pi_inf, gamma, lit_gamma @:HardcodedDimensionsExtrusion() @:Hardcoded1DVariables() - pi_inf = pi_infs(1) - gamma = gammas(1) - lit_gamma = gs_min(1) - ! Transferring the patch's centroid and length information x_centroid = patch_icpp(patch_id)%x_centroid length_x = patch_icpp(patch_id)%length_x diff --git a/src/simulation/m_acoustic_src.fpp b/src/simulation/m_acoustic_src.fpp index 2ad981fc2..c478405d7 100644 --- a/src/simulation/m_acoustic_src.fpp +++ b/src/simulation/m_acoustic_src.fpp @@ -262,8 +262,7 @@ contains if (dipole(ai)) then ! Double amplitude & No momentum source term (only works for Planar) mass_src(j, k, l) = mass_src(j, k, l) + 2._wp*mom_src_diff/c - if (model_eqns /= model_eqns_4eq) E_src(j, k, l) = E_src(j, k, & - & l) + 2._wp*mom_src_diff*c/(small_gamma - 1._wp) + E_src(j, k, l) = E_src(j, k, l) + 2._wp*mom_src_diff*c/(small_gamma - 1._wp) cycle end if @@ -301,9 +300,7 @@ contains mass_src(j, k, l) = mass_src(j, k, l) + mass_src_diff ! Update energy source term - if (model_eqns /= model_eqns_4eq) then - E_src(j, k, l) = E_src(j, k, l) + mass_src_diff*c**2._wp/(small_gamma - 1._wp) - end if + E_src(j, k, l) = E_src(j, k, l) + mass_src_diff*c**2._wp/(small_gamma - 1._wp) end do $:END_GPU_PARALLEL_LOOP() end if diff --git a/src/simulation/m_bubbles_EL.fpp b/src/simulation/m_bubbles_EL.fpp index fb4ed20c3..b4ed0f831 100644 --- a/src/simulation/m_bubbles_EL.fpp +++ b/src/simulation/m_bubbles_EL.fpp @@ -377,9 +377,6 @@ contains ! Initial particle pressure gas_p(bub_id, 1) = pliq + 2._wp*(1._wp/Web)/bub_R0(bub_id) - if (.not. f_approx_equal((1._wp/Web), 0._wp)) then - pref = gas_p(bub_id, 1) - end if ! Initial particle mass volparticle = 4._wp/3._wp*pi*bub_R0(bub_id)**3._wp ! volume diff --git a/src/simulation/m_data_output.fpp b/src/simulation/m_data_output.fpp index 998ab07b4..e4c32e2c9 100644 --- a/src/simulation/m_data_output.fpp +++ b/src/simulation/m_data_output.fpp @@ -19,7 +19,7 @@ module m_data_output use m_delay_file_access use m_ibm use m_boundary_common - use m_constants, only: model_eqns_5eq, model_eqns_4eq, precision_single + use m_constants, only: model_eqns_5eq, precision_single implicit none @@ -1242,9 +1242,7 @@ contains & dyn_p, pi_inf, gamma, rho, qv, rhoYks, pres, T) end if - if (model_eqns == model_eqns_4eq) then - lit_gamma = gammas(1) - else if (hypoelasticity) then + if (hypoelasticity) then tau_e(1) = q_cons_vf(eqn_idx%stress%end)%sf(j - 2, k, l)/rho end if @@ -1346,9 +1344,7 @@ contains & k - 2, l), dyn_p, pi_inf, gamma, rho, qv, rhoYks, pres, T) end if - if (model_eqns == model_eqns_4eq) then - lit_gamma = gs_min(1) - else if (hypoelasticity) then + if (hypoelasticity) then do s = 1, 3 tau_e(s) = q_cons_vf(s)%sf(j - 2, k - 2, l)/rho end do diff --git a/src/simulation/m_global_parameters.fpp b/src/simulation/m_global_parameters.fpp index 7108ad17d..4af35c426 100644 --- a/src/simulation/m_global_parameters.fpp +++ b/src/simulation/m_global_parameters.fpp @@ -766,27 +766,6 @@ contains end if end if - if (model_eqns == model_eqns_4eq .and. bubbles_euler) then - @:ALLOCATE(qbmm_idx%rs(nb), qbmm_idx%vs(nb)) - @:ALLOCATE(qbmm_idx%ps(nb), qbmm_idx%ms(nb)) - - do i = 1, nb - if (polytropic) then - fac = 2 - else - fac = 4 - end if - - qbmm_idx%rs(i) = eqn_idx%bub%beg + (i - 1)*fac - qbmm_idx%vs(i) = qbmm_idx%rs(i) + 1 - - if (.not. polytropic) then - qbmm_idx%ps(i) = qbmm_idx%vs(i) + 1 - qbmm_idx%ms(i) = qbmm_idx%ps(i) + 1 - end if - end do - end if - ! sim-only: Re_idx (non-gamma-law models only) if (model_eqns /= model_eqns_gamma_law) then ! Count fluids with non-negligible viscous effects (Re > 0) @@ -913,8 +892,8 @@ contains $:GPU_UPDATE(device='[cfl_target, m, n, p]') $:GPU_UPDATE(device='[alt_soundspeed, acoustic_source, num_source]') - $:GPU_UPDATE(device='[dt, sys_size, buff_size, pref, rhoref, eqn_idx, mpp_lim, bubbles_euler, hypoelasticity, & - & alt_soundspeed, avg_state, model_eqns, mixture_err, grid_geometry, cyl_coord, mp_weno, weno_eps, teno_CT, low_Mach]') + $:GPU_UPDATE(device='[dt, sys_size, buff_size, eqn_idx, mpp_lim, bubbles_euler, hypoelasticity, alt_soundspeed, & + & avg_state, model_eqns, mixture_err, grid_geometry, cyl_coord, mp_weno, weno_eps, teno_CT, low_Mach]') $:GPU_UPDATE(device='[Bx0]') diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp index e89c9686a..7d54dba97 100644 --- a/src/simulation/m_ibm.fpp +++ b/src/simulation/m_ibm.fpp @@ -293,14 +293,12 @@ contains end do end if - if (model_eqns /= model_eqns_4eq) then - ! If in simulation, use acc mixture subroutines - if (hypoelasticity) then - call s_convert_species_to_mixture_variables_kernel(rho, gamma, pi_inf, qv_K, alpha_IP, alpha_rho_IP, & - & Re_K, G_K, Gs) - else - call s_convert_species_to_mixture_variables_kernel(rho, gamma, pi_inf, qv_K, alpha_IP, alpha_rho_IP, Re_K) - end if + ! If in simulation, use acc mixture subroutines + if (hypoelasticity) then + call s_convert_species_to_mixture_variables_kernel(rho, gamma, pi_inf, qv_K, alpha_IP, alpha_rho_IP, Re_K, & + & G_K, Gs) + else + call s_convert_species_to_mixture_variables_kernel(rho, gamma, pi_inf, qv_K, alpha_IP, alpha_rho_IP, Re_K) end if if (patch_ib(patch_id)%moving_ibm /= 0) then diff --git a/src/simulation/m_riemann_solver_hllc.fpp b/src/simulation/m_riemann_solver_hllc.fpp index 2f1ab5a54..93a77ac07 100644 --- a/src/simulation/m_riemann_solver_hllc.fpp +++ b/src/simulation/m_riemann_solver_hllc.fpp @@ -14,7 +14,7 @@ module m_riemann_solver_hllc use m_variables_conversion use m_bubbles use m_constants, only: riemann_solver_hll, riemann_solver_hllc, riemann_solver_lax_friedrichs, model_eqns_5eq, & - & model_eqns_6eq, model_eqns_4eq, avg_state_roe, avg_state_arithmetic, wave_speeds_direct, wave_speeds_pressure + & model_eqns_6eq, avg_state_roe, avg_state_arithmetic, wave_speeds_direct, wave_speeds_pressure use m_bubbles_EE use m_surface_tension use m_chemistry @@ -421,215 +421,6 @@ contains end do end do $:END_GPU_PARALLEL_LOOP() - else if (model_eqns == model_eqns_4eq) then - ! 4-equation model (model_eqns=4): single pressure, velocity equilibrium - $:GPU_PARALLEL_LOOP(collapse=3, private='[i, q, alpha_rho_L, alpha_rho_R, vel_L, vel_R, alpha_L, alpha_R, & - & nbub_L, nbub_R, rho_L, rho_R, pres_L, pres_R, E_L, E_R, H_L, H_R, Cp_avg, Cv_avg, & - & T_avg, eps, c_sum_Yi_Phi, T_L, T_R, Y_L, Y_R, MW_L, MW_R, R_gas_L, R_gas_R, Cp_L, Cp_R, & - & Gamm_L, Gamm_R, gamma_L, gamma_R, pi_inf_L, pi_inf_R, qv_L, qv_R, qv_avg, c_L, c_R, & - & rho_avg, H_avg, c_avg, gamma_avg, ptilde_L, ptilde_R, vel_L_rms, vel_R_rms, & - & vel_avg_rms, vel_L_tmp, vel_R_tmp, Ms_L, Ms_R, pres_SL, pres_SR, alpha_L_sum, & - & alpha_R_sum, rho_Star, E_Star, p_Star, p_K_Star, vel_K_star, s_L, s_R, s_M, s_P, s_S, & - & xi_M, xi_P, xi_L, xi_R, xi_L_m1, xi_R_m1, xi_MP, xi_PP, Ys_L, Ys_R, Cp_iL, Cp_iR, Xs_L, & - & Xs_R, Gamma_iL, Gamma_iR, Yi_avg, Phi_avg, h_iL, h_iR, h_avg_2]') - do l = ${Z_BND}$%beg, ${Z_BND}$%end - do k = ${Y_BND}$%beg, ${Y_BND}$%end - do j = ${X_BND}$%beg, ${X_BND}$%end - vel_L_rms = 0._wp; vel_R_rms = 0._wp - rho_L = 0._wp; rho_R = 0._wp - gamma_L = 0._wp; gamma_R = 0._wp - pi_inf_L = 0._wp; pi_inf_R = 0._wp - qv_L = 0._wp; qv_R = 0._wp - - $:GPU_LOOP(parallelism='[seq]') - do i = 1, eqn_idx%cont%end - alpha_rho_L(i) = qL_prim_rsx_vf(${SF('')}$, i) - alpha_rho_R(i) = qR_prim_rsx_vf(${SF(' + 1')}$, i) - end do - - $:GPU_LOOP(parallelism='[seq]') - do i = 1, num_dims - vel_L(i) = qL_prim_rsx_vf(${SF('')}$, eqn_idx%cont%end + i) - vel_R(i) = qR_prim_rsx_vf(${SF(' + 1')}$, eqn_idx%cont%end + i) - vel_L_rms = vel_L_rms + vel_L(i)**2._wp - vel_R_rms = vel_R_rms + vel_R(i)**2._wp - end do - - $:GPU_LOOP(parallelism='[seq]') - do i = 1, num_fluids - alpha_L(i) = qL_prim_rsx_vf(${SF('')}$, eqn_idx%E + i) - alpha_R(i) = qR_prim_rsx_vf(${SF(' + 1')}$, eqn_idx%E + i) - end do - $:GPU_LOOP(parallelism='[seq]') - do i = 1, num_fluids - alpha_L(i) = qL_prim_rsx_vf(${SF('')}$, eqn_idx%E + i) - alpha_R(i) = qR_prim_rsx_vf(${SF(' + 1')}$, eqn_idx%E + i) - end do - - call s_accumulate_mixture_properties(num_fluids, alpha_rho_L, alpha_L, rho_L, gamma_L, pi_inf_L, & - & qv_L) - call s_accumulate_mixture_properties(num_fluids, alpha_rho_R, alpha_R, rho_R, gamma_R, pi_inf_R, & - & qv_R) - - pres_L = qL_prim_rsx_vf(${SF('')}$, eqn_idx%E) - pres_R = qR_prim_rsx_vf(${SF(' + 1')}$, eqn_idx%E) - - E_L = gamma_L*pres_L + pi_inf_L + 5.e-1_wp*rho_L*vel_L_rms + qv_L - E_R = gamma_R*pres_R + pi_inf_R + 5.e-1_wp*rho_R*vel_R_rms + qv_R - - H_L = (E_L + pres_L)/rho_L - H_R = (E_R + pres_R)/rho_R - - @:compute_average_state() - - call s_compute_speed_of_sound(pres_L, rho_L, gamma_L, pi_inf_L, H_L, alpha_L, vel_L_rms, 0._wp, & - & c_L, qv_L) - - call s_compute_speed_of_sound(pres_R, rho_R, gamma_R, pi_inf_R, H_R, alpha_R, vel_R_rms, 0._wp, & - & c_R, qv_R) - - !> The computation of c_avg does not require all the variables, and therefore the non '_avg' - ! variables are placeholders to call the subroutine. - - call s_compute_speed_of_sound(pres_R, rho_avg, gamma_avg, pi_inf_R, H_avg, alpha_R, vel_avg_rms, & - & 0._wp, c_avg, qv_avg) - - if (wave_speeds == wave_speeds_direct) then - s_L = min(vel_L(dir_idx(1)) - c_L, vel_R(dir_idx(1)) - c_R) - s_R = max(vel_R(dir_idx(1)) + c_R, vel_L(dir_idx(1)) + c_L) - - s_S = (pres_R - pres_L + rho_L*vel_L(dir_idx(1))*(s_L - vel_L(dir_idx(1))) & - & - rho_R*vel_R(dir_idx(1))*(s_R - vel_R(dir_idx(1))))/(rho_L*(s_L - vel_L(dir_idx(1))) & - & - rho_R*(s_R - vel_R(dir_idx(1)))) - else if (wave_speeds == wave_speeds_pressure) then - pres_SL = 5.e-1_wp*(pres_L + pres_R + rho_avg*c_avg*(vel_L(dir_idx(1)) - vel_R(dir_idx(1)))) - - pres_SR = pres_SL - - ! Low Mach correction: Thornber et al. JCP (2008) - Ms_L = max(1._wp, & - & sqrt(1._wp + ((5.e-1_wp + gamma_L)/(1._wp + gamma_L))*(pres_SL/pres_L - 1._wp) & - & *pres_L/((pres_L + pi_inf_L/(1._wp + gamma_L))))) - Ms_R = max(1._wp, & - & sqrt(1._wp + ((5.e-1_wp + gamma_R)/(1._wp + gamma_R))*(pres_SR/pres_R - 1._wp) & - & *pres_R/((pres_R + pi_inf_R/(1._wp + gamma_R))))) - - s_L = vel_L(dir_idx(1)) - c_L*Ms_L - s_R = vel_R(dir_idx(1)) + c_R*Ms_R - - s_S = 5.e-1_wp*((vel_L(dir_idx(1)) + vel_R(dir_idx(1))) + (pres_L - pres_R)/(rho_avg*c_avg)) - end if - - ! follows Einfeldt et al. s_M/P = min/max(0.,s_L/R) - s_M = min(0._wp, s_L); s_P = max(0._wp, s_R) - - ! goes with q_star_L/R = xi_L/R * (variable) xi_L/R = ( ( s_L/R - u_L/R )/(s_L/R - s_star) ) - xi_L = (s_L - vel_L(dir_idx(1)))/min(s_L - s_S, -sgm_eps) - xi_R = (s_R - vel_R(dir_idx(1)))/max(s_R - s_S, sgm_eps) - xi_L_m1 = (s_S - vel_L(dir_idx(1)))/min(s_L - s_S, -sgm_eps) - xi_R_m1 = (s_S - vel_R(dir_idx(1)))/max(s_R - s_S, sgm_eps) - - ! goes with numerical velocity in x/y/z directions xi_P/M = 0.5 +/m sgn(0.5,s_star) - xi_M = (5.e-1_wp + sign(5.e-1_wp, s_S)) - xi_P = (5.e-1_wp - sign(5.e-1_wp, s_S)) - - $:GPU_LOOP(parallelism='[seq]') - do i = 1, eqn_idx%cont%end - flux_rsx_vf(${SF('')}$, & - & i) = xi_M*alpha_rho_L(i)*(vel_L(dir_idx(1)) + s_M*xi_L_m1) + xi_P*alpha_rho_R(i) & - & *(vel_R(dir_idx(1)) + s_P*xi_R_m1) - end do - - ! Momentum flux. f = \rho u u + p I, q = \rho u, q_star = \xi * \rho*(s_star, v, w) - $:GPU_LOOP(parallelism='[seq]') - do i = 1, num_dims - flux_rsx_vf(${SF('')}$, & - & eqn_idx%cont%end + dir_idx(i)) = xi_M*(rho_L*(vel_L(dir_idx(1))*vel_L(dir_idx(i) & - & ) + s_M*(xi_L*(dir_flg(dir_idx(i))*s_S + (1._wp - dir_flg(dir_idx(i))) & - & *vel_L(dir_idx(i))) - vel_L(dir_idx(i)))) + dir_flg(dir_idx(i))*pres_L) & - & + xi_P*(rho_R*(vel_R(dir_idx(1))*vel_R(dir_idx(i)) & - & + s_P*(xi_R*(dir_flg(dir_idx(i))*s_S + (1._wp - dir_flg(dir_idx(i))) & - & *vel_R(dir_idx(i))) - vel_R(dir_idx(i)))) + dir_flg(dir_idx(i))*pres_R) - end do - - if (bubbles_euler) then - ! Put p_tilde in - $:GPU_LOOP(parallelism='[seq]') - do i = 1, num_dims - flux_rsx_vf(${SF('')}$, eqn_idx%cont%end + dir_idx(i)) = flux_rsx_vf(${SF('')}$, & - & eqn_idx%cont%end + dir_idx(i)) + xi_M*(dir_flg(dir_idx(i))*(-1._wp*ptilde_L) & - & ) + xi_P*(dir_flg(dir_idx(i))*(-1._wp*ptilde_R)) - end do - end if - - flux_rsx_vf(${SF('')}$, eqn_idx%E) = 0._wp - - $:GPU_LOOP(parallelism='[seq]') - do i = eqn_idx%alf, eqn_idx%alf ! only advect the void fraction - flux_rsx_vf(${SF('')}$, i) = xi_M*qL_prim_rsx_vf(${SF('')}$, & - & i)*(vel_L(dir_idx(1)) + s_M*xi_L_m1) + xi_P*qR_prim_rsx_vf(${SF(' + 1')}$, & - & i)*(vel_R(dir_idx(1)) + s_P*xi_R_m1) - end do - - ! Advection velocity source: interface velocity for volume fraction transport - $:GPU_LOOP(parallelism='[seq]') - do i = 1, num_dims - vel_src_rsx_vf(${SF('')}$, dir_idx(i)) = 0._wp - ! IF ( (model_eqns == 4) .or. (num_fluids==1) ) vel_src_rs_vf(dir_idx(i))%sf(j,k,l) = 0._wp - end do - - flux_src_rsx_vf(${SF('')}$, eqn_idx%adv%beg) = vel_src_rsx_vf(${SF('')}$, dir_idx(1)) - - ! Add advection flux for bubble variables - if (bubbles_euler) then - $:GPU_LOOP(parallelism='[seq]') - do i = eqn_idx%bub%beg, eqn_idx%bub%end - flux_rsx_vf(${SF('')}$, i) = xi_M*nbub_L*qL_prim_rsx_vf(${SF('')}$, & - & i)*(vel_L(dir_idx(1)) + s_M*xi_L_m1) & - & + xi_P*nbub_R*qR_prim_rsx_vf(${SF(' + 1')}$, & - & i)*(vel_R(dir_idx(1)) + s_P*xi_R_m1) - end do - end if - - ! Geometrical source flux for cylindrical coordinates - - #:if (NORM_DIR == 2) - if (cyl_coord) then - ! Substituting the advective flux into the inviscid geometrical source flux - $:GPU_LOOP(parallelism='[seq]') - do i = 1, eqn_idx%E - flux_gsrc_rsx_vf(${SF('')}$, i) = flux_rsx_vf(${SF('')}$, i) - end do - ! Recalculating the radial momentum geometric source flux - flux_gsrc_rsx_vf(${SF('')}$, & - & eqn_idx%cont%end + dir_idx(1)) & - & = f_compute_hllc_star_momentum_flux(rho_L, rho_R, vel_L(dir_idx(1)), & - & vel_R(dir_idx(1)), s_M, s_P, s_S, xi_L, xi_R, xi_M, xi_P, & - & dir_flg(dir_idx(1))) - ! Geometrical source of the void fraction(s) is zero - $:GPU_LOOP(parallelism='[seq]') - do i = eqn_idx%adv%beg, eqn_idx%adv%end - flux_gsrc_rsx_vf(${SF('')}$, i) = 0._wp - end do - end if - #:endif - #:if (NORM_DIR == 3) - if (grid_geometry == 3) then - $:GPU_LOOP(parallelism='[seq]') - do i = 1, sys_size - flux_gsrc_rsx_vf(${SF('')}$, i) = 0._wp - end do - flux_gsrc_rsx_vf(${SF('')}$, & - & eqn_idx%mom%beg + 1) = -f_compute_hllc_star_momentum_flux(rho_L, & - & rho_R, vel_L(dir_idx(1)), vel_R(dir_idx(1)), s_M, s_P, s_S, xi_L, & - & xi_R, xi_M, xi_P, dir_flg(dir_idx(1))) - flux_gsrc_rsx_vf(${SF('')}$, eqn_idx%mom%end) = flux_rsx_vf(${SF('')}$, eqn_idx%mom%beg + 1) - end if - #:endif - end do - end do - end do - $:END_GPU_PARALLEL_LOOP() else if (model_eqns == model_eqns_5eq .and. bubbles_euler) then ! 5-equation model with Euler-Euler bubble dynamics $:GPU_PARALLEL_LOOP(collapse=3, private='[i, q, R0_L, R0_R, V0_L, V0_R, P0_L, P0_R, pbw_L, pbw_R, vel_L, & @@ -945,8 +736,6 @@ contains vel_src_rsx_vf(${SF('')}$, & & dir_idx(i)) = xi_M*(vel_L(dir_idx(i)) + dir_flg(dir_idx(i))*s_M*xi_L_m1) & & + xi_P*(vel_R(dir_idx(i)) + dir_flg(dir_idx(i))*s_P*xi_R_m1) - - ! IF ( (model_eqns == 4) .or. (num_fluids==1) ) vel_src_rs_vf(dir_idx(i))%sf(j,k,l) = 0._wp end do flux_src_rsx_vf(${SF('')}$, eqn_idx%adv%beg) = vel_src_rsx_vf(${SF('')}$, dir_idx(1)) diff --git a/toolchain/mfc/case_validator.py b/toolchain/mfc/case_validator.py index 03d6421b5..e85336a9f 100644 --- a/toolchain/mfc/case_validator.py +++ b/toolchain/mfc/case_validator.py @@ -349,12 +349,11 @@ def check_model_eqns_and_num_fluids(self): cyl_coord = self.get("cyl_coord", "F") == "T" p = self.get("p", 0) - self.prohibit(model_eqns is not None and model_eqns not in [1, 2, 3, 4], "model_eqns must be 1, 2, 3, or 4") + self.prohibit(model_eqns is not None and model_eqns not in [1, 2, 3], "model_eqns must be 1, 2, or 3") self.prohibit(num_fluids is not None and num_fluids < 1, "num_fluids must be positive") self.prohibit(model_eqns == 1 and num_fluids is not None, "num_fluids is not supported for model_eqns = 1") self.prohibit(model_eqns == 2 and num_fluids is None, "5-equation model (model_eqns = 2) requires num_fluids to be set") self.prohibit(model_eqns == 3 and num_fluids is None, "6-equation model (model_eqns = 3) requires num_fluids to be set") - self.prohibit(model_eqns == 4 and num_fluids is None, "4-equation model (model_eqns = 4) requires num_fluids to be set") self.prohibit(model_eqns == 1 and mpp_lim, "model_eqns = 1 does not support mpp_lim") self.prohibit(num_fluids == 1 and mpp_lim, "num_fluids = 1 does not support mpp_lim") self.prohibit(model_eqns == 3 and cyl_coord and p != 0, "6-equation model (model_eqns = 3) does not support cylindrical coordinates (cyl_coord = T and p != 0)") @@ -502,9 +501,6 @@ def check_bubbles_euler(self): thermal = self.get("thermal") model_eqns = self.get("model_eqns") cyl_coord = self.get("cyl_coord", "F") == "T" - rhoref = self.get("rhoref") - pref = self.get("pref") - num_fluids = self.get("num_fluids") self.prohibit(nb is None or nb < 1, "The Ensemble-Averaged Bubble Model requires nb >= 1") self.prohibit(polydisperse and nb == 1, "Polydisperse bubble dynamics requires nb > 1") @@ -512,11 +508,6 @@ def check_bubbles_euler(self): self.prohibit(thermal is not None and thermal > 3, "thermal must be <= 3") self.prohibit(model_eqns == 3, "Bubble models untested with 6-equation model (model_eqns = 3)") self.prohibit(model_eqns == 1, "Bubble models untested with pi-gamma model (model_eqns = 1)") - self.prohibit(model_eqns == 4 and rhoref is None, "rhoref must be set if using bubbles_euler with model_eqns = 4") - self.prohibit(rhoref is not None and rhoref <= 0, "rhoref (reference density) must be positive") - self.prohibit(model_eqns == 4 and pref is None, "pref must be set if using bubbles_euler with model_eqns = 4") - self.prohibit(pref is not None and pref <= 0, "pref (reference pressure) must be positive") - self.prohibit(model_eqns == 4 and num_fluids != 1, "4-equation model (model_eqns = 4) is single-component and requires num_fluids = 1") self.prohibit(cyl_coord, "Bubble models untested in cylindrical coordinates") # BUBBLE PHYSICS PARAMETERS @@ -2016,14 +2007,14 @@ def check_domain_bounds(self): def check_volume_fraction_sum(self): """Warns if volume fractions do not sum to 1 for multi-component models. - For model_eqns in [2, 3, 4], the mixture constraint sum(alpha_j) = 1 + For model_eqns in [2, 3], the mixture constraint sum(alpha_j) = 1 must hold. Skips patches with analytical expressions, alter_patch, hcid, bubbles_euler single-fluid cases (where alpha represents the void fraction, not a partition of unity), and bubbles_lagrange cases (where the Lagrangian phase is not tracked on the Euler grid). """ model_eqns = self.get("model_eqns") - if model_eqns not in [2, 3, 4]: + if model_eqns not in [2, 3]: return num_patches = self.get("num_patches", 0) diff --git a/toolchain/mfc/fp_stability.py b/toolchain/mfc/fp_stability.py index 07ef7aae0..5869d743e 100644 --- a/toolchain/mfc/fp_stability.py +++ b/toolchain/mfc/fp_stability.py @@ -282,8 +282,6 @@ def _merge(*dicts): "polytropic": "T", "polydisperse": "F", "thermal": 3, - "pref": 101325.0, - "rhoref": 1000.0, "patch_icpp(1)%pres": 2.0, "patch_icpp(1)%alpha_rho(1)": 0.96, "patch_icpp(1)%alpha(1)": 0.04, @@ -319,8 +317,6 @@ def _merge(*dicts): "polytropic": "T", "polydisperse": "F", "thermal": 3, - "pref": 101325.0, - "rhoref": 1000.0, "fluid_pp(1)%gamma": 0.16, "fluid_pp(1)%pi_inf": 3515.0, "bub_pp%R0ref": 1.0, diff --git a/toolchain/mfc/params/definitions.py b/toolchain/mfc/params/definitions.py index fd32e8c3a..f5daed04b 100644 --- a/toolchain/mfc/params/definitions.py +++ b/toolchain/mfc/params/definitions.py @@ -337,9 +337,9 @@ def get_value_label(param_name: str, value: int) -> str: }, # Model equations "model_eqns": { - "choices": [1, 2, 3, 4], - "value_labels": {1: "Gamma-law", 2: "5-Equation", 3: "6-Equation", 4: "4-Equation"}, - "names": {"gamma_law": 1, "5eq": 2, "6eq": 3, "4eq": 4}, + "choices": [1, 2, 3], + "value_labels": {1: "Gamma-law", 2: "5-Equation", 3: "6-Equation"}, + "names": {"gamma_law": 1, "5eq": 2, "6eq": 3}, }, # Bubbles "bubble_model": { @@ -392,7 +392,7 @@ def get_value_label(param_name: str, value: int) -> str: "when_true": { "recommends": ["nb", "polytropic"], "requires_value": { - "model_eqns": [2, 4], + "model_eqns": [2], "riemann_solver": [2], "avg_state": [2], }, @@ -402,7 +402,6 @@ def get_value_label(param_name: str, value: int) -> str: "when_value": { 2: {"requires": ["num_fluids"]}, 3: {"requires_value": {"riemann_solver": [2], "avg_state": [2], "wave_speeds": [1]}}, - 4: {"requires": ["rhoref", "pref"], "requires_value": {"num_fluids": [1]}}, } }, "viscous": { @@ -736,9 +735,7 @@ def _load(): "flux_lim", ]: _r(n, INT) - _r("pref", REAL, math=r"\f$p_\text{ref}\f$") _r("poly_sigma", REAL, math=r"\f$\sigma_\text{poly}\f$") - _r("rhoref", REAL, math=r"\f$\rho_\text{ref}\f$") _r("palpha_eps", REAL, math=r"\f$\varepsilon_\alpha\f$") _r("ptgalpha_eps", REAL, math=r"\f$\varepsilon_\alpha\f$") _r("pi_fac", REAL, math=r"\f$\pi\text{-factor}\f$") @@ -1255,8 +1252,6 @@ def _decl(targets: set, *names: str) -> None: "relax_model", "fluid_pp", "bub_pp", - "rhoref", - "pref", "bubbles_euler", "bubbles_lagrange", "R0ref", diff --git a/toolchain/mfc/params/descriptions.py b/toolchain/mfc/params/descriptions.py index bd09bb4f6..020b2160b 100644 --- a/toolchain/mfc/params/descriptions.py +++ b/toolchain/mfc/params/descriptions.py @@ -158,9 +158,7 @@ "thermal": "Thermal model selection", "relax_model": "Relaxation model type", "igr_order": "Implicit gradient reconstruction order", - "pref": "Reference pressure", "poly_sigma": "Polydisperse distribution standard deviation", - "rhoref": "Reference density", "sigma": "Surface tension coefficient", "Bx0": "Background magnetic field in x-direction", "relax": "Enable relaxation terms", diff --git a/toolchain/mfc/params/generators/fortran_gen.py b/toolchain/mfc/params/generators/fortran_gen.py index 8209334e2..958a19b0c 100644 --- a/toolchain/mfc/params/generators/fortran_gen.py +++ b/toolchain/mfc/params/generators/fortran_gen.py @@ -93,14 +93,12 @@ "poly_sigma", "polydisperse", "polytropic", - "pref", "ptgalpha_eps", "qbmm", "recon_type", "relativity", "relax", "relax_model", - "rhoref", "sigma", "surface_tension", "synth_U_inf", diff --git a/toolchain/mfc/params_tests/mutation_tests.py b/toolchain/mfc/params_tests/mutation_tests.py index aecb7a62c..d25eb0a0a 100644 --- a/toolchain/mfc/params_tests/mutation_tests.py +++ b/toolchain/mfc/params_tests/mutation_tests.py @@ -88,9 +88,6 @@ class MutationResult: "bub_pp%mu_g": [-1, -1e-3], # Surface tension must be non-negative "bub_pp%ss": [-1, -0.01], - # Global bubble reference values - "rhoref": [0, -1, -1000], - "pref": [0, -1, -1e5], # PHYSICS: ACOUSTICS # Frequency/wavelength must be positive "acoustic(1)%frequency": [0, -1], diff --git a/toolchain/mfc/params_tests/test_fortran_gen.py b/toolchain/mfc/params_tests/test_fortran_gen.py index a57175514..dffcad4ee 100644 --- a/toolchain/mfc/params_tests/test_fortran_gen.py +++ b/toolchain/mfc/params_tests/test_fortran_gen.py @@ -432,9 +432,9 @@ def test_generate_bcast_fpp_class_a_real_scalars(): assert "call MPI_BCAST(dt, 1, mpi_p, 0, MPI_COMM_WORLD, ierr)" in sim assert "call MPI_BCAST(dt, " not in pre and "call MPI_BCAST(dt, " not in post - # pref is in all three + # poly_sigma is in all three for out, target in [(pre, "pre"), (sim, "sim"), (post, "post")]: - assert "call MPI_BCAST(pref, 1, mpi_p, 0, MPI_COMM_WORLD, ierr)" in out, f"{target}: pref missing" + assert "call MPI_BCAST(poly_sigma, 1, mpi_p, 0, MPI_COMM_WORLD, ierr)" in out, f"{target}: poly_sigma missing" def test_generate_bcast_fpp_class_a_str_scalars(): diff --git a/toolchain/mfc/test/case.py b/toolchain/mfc/test/case.py index 837749791..afd457ec5 100644 --- a/toolchain/mfc/test/case.py +++ b/toolchain/mfc/test/case.py @@ -105,8 +105,6 @@ def get_post_process_mods(case_params: dict) -> dict: "fluid_pp(1)%qv": 0.0, "fluid_pp(1)%qvp": 0.0, "bubbles_euler": "F", - "pref": 101325.0, - "rhoref": 1000.0, "bubble_model": 3, "polytropic": "T", "polydisperse": "F", diff --git a/toolchain/mfc/test_lint_source.py b/toolchain/mfc/test_lint_source.py index 6c788543f..32ff3494e 100644 --- a/toolchain/mfc/test_lint_source.py +++ b/toolchain/mfc/test_lint_source.py @@ -18,11 +18,11 @@ def test_extract_indexed_argument_keeps_root(): def test_extract_fypp_list_with_continuation(): lines = [ " #:for VAR in [ 'weno_eps', 'teno_CT', &", - " & 'pref']", + " & 'poly_sigma']", " call MPI_BCAST(${VAR}$" + BCAST_TAIL, " #:endfor", ] - assert _extract_bcast_roots(lines) == [(3, "weno_eps"), (3, "teno_CT"), (3, "pref")] + assert _extract_bcast_roots(lines) == [(3, "weno_eps"), (3, "teno_CT"), (3, "poly_sigma")] def test_struct_members_and_loop_indices_skipped(): @@ -50,11 +50,11 @@ def _write_proxy(tmp_path, target_dir: str, body: str): def test_manual_broadcast_of_registry_scalar_is_flagged(tmp_path): - _write_proxy(tmp_path, "simulation", f" call MPI_BCAST(rhoref{BCAST_TAIL}\n") + _write_proxy(tmp_path, "simulation", f" call MPI_BCAST(poly_sigma{BCAST_TAIL}\n") errors = check_manual_registry_bcasts(tmp_path) assert len(errors) == 1 - assert "manual MPI_BCAST of registry-bound scalar 'rhoref'" in errors[0] + assert "manual MPI_BCAST of registry-bound scalar 'poly_sigma'" in errors[0] assert "m_mpi_proxy.fpp:1" in errors[0]