MFlowCode · sbryngelson · Aug 9, 2026 · Aug 9, 2026
@@ -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 |

@@ -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).

@@ -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).
 
 ---
 

@@ -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

@@ -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(*)|)

@@ -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

@@ -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

@@ -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

@@ -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