Make fluid variables uniform through the code (#1063)

Co-authored-by: qodo-merge-pro[bot] <151058649+qodo-merge-pro[bot]@users.noreply.github.com>
Co-authored-by: cubic-dev-ai[bot] <191113872+cubic-dev-ai[bot]@users.noreply.github.com>

Author: 3 people

src/common/m_variables_conversion.fpp

@@ -1306,8 +1306,7 @@ contains
                         elseif ((model_eqns /= 4) .and. (bubbles_euler .neqv. .true.)) then
                             ! E = Gamma*P + \rho u u /2 + \pi_inf + (\alpha\rho qv)
                             q_cons_vf(E_idx)%sf(j, k, l) = &
-                                gamma*q_prim_vf(E_idx)%sf(j, k, l) + dyn_pres + pi_inf &
-                                + qv
+                                gamma*q_prim_vf(E_idx)%sf(j, k, l) + dyn_pres + pi_inf + qv
                         else if ((model_eqns /= 4) .and. (bubbles_euler)) then
                             ! \tilde{E} = dyn_pres + (1-\alf)(\Gamma p_l + \Pi_inf)
                             q_cons_vf(E_idx)%sf(j, k, l) = dyn_pres + &
@@ -1323,11 +1322,9 @@ contains
                     if (model_eqns == 3) then
                         do i = 1, num_fluids
                             ! internal energy calculation for each of the fluids
-                            q_cons_vf(i + internalEnergies_idx%beg - 1)%sf(j, k, l) = &
-                                q_cons_vf(i + adv_idx%beg - 1)%sf(j, k, l)* &
-                                (fluid_pp(i)%gamma*q_prim_vf(E_idx)%sf(j, k, l) + &
-                                 fluid_pp(i)%pi_inf) + &
-                                q_cons_vf(i + cont_idx%beg - 1)%sf(j, k, l)*fluid_pp(i)%qv
+                            q_cons_vf(i + intxb - 1)%sf(j, k, l) = q_cons_vf(i + advxb - 1)%sf(j, k, l)* &
+                                                                   (gammas(i)*q_prim_vf(E_idx)%sf(j, k, l) + pi_infs(i)) + &
+                                                                   q_cons_vf(i + contxb - 1)%sf(j, k, l)*qvs(i)
                         end do
                     end if
 
@@ -1637,7 +1634,7 @@ contains
                 c = 0._wp
                 $:GPU_LOOP(parallelism='[seq]')
                 do q = 1, num_fluids
-                    c = c + adv(q)*(1._wp/gammas(q) + 1._wp)* &
+                    c = c + adv(q)*gs_min(q)* &
                         (pres + pi_infs(q)/(gammas(q) + 1._wp))
                 end do
                 c = c/rho

src/post_process/m_data_output.fpp

@@ -21,6 +21,8 @@ module m_data_output
 
     use m_helper
 
+    use m_variables_conversion
+
     implicit none
 
     private; public :: s_initialize_data_output_module, &
@@ -1612,7 +1614,7 @@ contains
                     pi_inf = 0._wp
                     qv = 0._wp
                     pres = q_prim_vf(E_idx)%sf(i, j, k)
-                    Egint = Egint + q_prim_vf(E_idx + 2)%sf(i, j, k)*(fluid_pp(2)%gamma*pres)*dV
+                    Egint = Egint + q_prim_vf(E_idx + 2)%sf(i, j, k)*(gammas(2)*pres)*dV
                     do s = 1, num_vels
                         vel(s) = q_prim_vf(num_fluids + s)%sf(i, j, k)
                         Egk = Egk + 0.5_wp*q_prim_vf(E_idx + 2)%sf(i, j, k)*q_prim_vf(2)%sf(i, j, k)*vel(s)*vel(s)*dV
@@ -1623,10 +1625,10 @@ contains
                     end do
                     do l = 1, adv_idx%end - E_idx
                         adv(l) = q_prim_vf(E_idx + l)%sf(i, j, k)
-                        gamma = gamma + adv(l)*fluid_pp(l)%gamma
-                        pi_inf = pi_inf + adv(l)*fluid_pp(l)%pi_inf
+                        gamma = gamma + adv(l)*gammas(l)
+                        pi_inf = pi_inf + adv(l)*pi_infs(l)
                         rho = rho + adv(l)*q_prim_vf(l)%sf(i, j, k)
-                        qv = qv + adv(l)*q_prim_vf(l)%sf(i, j, k)*fluid_pp(l)%qv
+                        qv = qv + adv(l)*q_prim_vf(l)%sf(i, j, k)*qvs(l)
                     end do
 
                     H = ((gamma + 1._wp)*pres + pi_inf + qv)/rho

src/post_process/m_derived_variables.fpp

@@ -204,10 +204,10 @@ contains
                                           pi_inf_sf(i, j, k))/(gamma_sf(i, j, k)* &
                                                                rho_sf(i, j, k)))
                     else
-                        blkmod1 = ((fluid_pp(1)%gamma + 1._wp)*q_prim_vf(E_idx)%sf(i, j, k) + &
-                                   fluid_pp(1)%pi_inf)/fluid_pp(1)%gamma
-                        blkmod2 = ((fluid_pp(2)%gamma + 1._wp)*q_prim_vf(E_idx)%sf(i, j, k) + &
-                                   fluid_pp(2)%pi_inf)/fluid_pp(2)%gamma
+                        blkmod1 = ((gammas(1) + 1._wp)*q_prim_vf(E_idx)%sf(i, j, k) + &
+                                   pi_infs(1))/gammas(1)
+                        blkmod2 = ((gammas(2) + 1._wp)*q_prim_vf(E_idx)%sf(i, j, k) + &
+                                   pi_infs(2))/gammas(2)
                         q_sf(i, j, k) = (1._wp/(rho_sf(i, j, k)*(q_prim_vf(adv_idx%beg)%sf(i, j, k)/blkmod1 + &
                                                                  (1._wp - q_prim_vf(adv_idx%beg)%sf(i, j, k))/blkmod2)))
                     end if

src/pre_process/m_assign_variables.fpp

@@ -213,11 +213,8 @@ contains
         p0 = 101325._wp
         pres_mag = 1.e-1_wp
         loc = x_cc(177)
-        n_tait = fluid_pp(1)%gamma
-        B_tait = fluid_pp(1)%pi_inf
-
-        n_tait = 1._wp/n_tait + 1._wp
-        B_tait = B_tait*(n_tait - 1._wp)/n_tait
+        n_tait = gs_min(1)
+        B_tait = ps_inf(1)
 
         if (j < 177) then
             q_prim_vf(E_idx)%sf(j, k, l) = 0.5_wp*q_prim_vf(E_idx)%sf(j, k, l)
@@ -563,9 +560,9 @@ contains
             end do
         else
             !get mixture density from pressure via Tait EOS
-            pi_inf = fluid_pp(1)%pi_inf
-            gamma = fluid_pp(1)%gamma
-            lit_gamma = (1._wp + gamma)/gamma
+            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) = &

src/pre_process/m_data_output.fpp

@@ -273,10 +273,10 @@ contains
             end do
         end if
 
-        gamma = fluid_pp(1)%gamma
-        lit_gamma = 1._wp/fluid_pp(1)%gamma + 1._wp
-        pi_inf = fluid_pp(1)%pi_inf
-        qv = fluid_pp(1)%qv
+        gamma = gammas(1)
+        lit_gamma = gs_min(1)
+        pi_inf = pi_infs(1)
+        qv = qvs(1)
 
         if (precision == 1) then
             FMT = "(2F30.3)"

src/pre_process/m_icpp_patches.fpp

@@ -31,6 +31,8 @@ module m_icpp_patches
 
     use m_ib_patches
 
+    use m_variables_conversion
+
     implicit none
 
     private; public :: s_apply_icpp_patches
@@ -210,9 +212,9 @@ contains
         @:HardcodedDimensionsExtrusion()
         @:Hardcoded1DVariables()
 
-        pi_inf = fluid_pp(1)%pi_inf
-        gamma = fluid_pp(1)%gamma
-        lit_gamma = (1._wp + gamma)/gamma
+        pi_inf = pi_infs(1)
+        gamma = gammas(1)
+        lit_gamma = gs_min(1)
         j = 0
         k = 0
 
@@ -748,9 +750,9 @@ contains
         @:HardcodedDimensionsExtrusion()
         @:Hardcoded2DVariables()
 
-        pi_inf = fluid_pp(1)%pi_inf
-        gamma = fluid_pp(1)%gamma
-        lit_gamma = (1._wp + gamma)/gamma
+        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
@@ -911,9 +913,9 @@ contains
         @:HardcodedDimensionsExtrusion()
         @:Hardcoded2DVariables()
 
-        pi_inf = fluid_pp(1)%pi_inf
-        gamma = fluid_pp(1)%gamma
-        lit_gamma = (1._wp + gamma)/gamma
+        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
@@ -997,9 +999,9 @@ contains
         @:HardcodedDimensionsExtrusion()
         @:Hardcoded1DVariables()
 
-        pi_inf = fluid_pp(1)%pi_inf
-        gamma = fluid_pp(1)%gamma
-        lit_gamma = (1._wp + gamma)/gamma
+        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

src/simulation/m_data_output.fpp

@@ -507,10 +507,10 @@ contains
             write (2) ib_markers%sf(0:m, 0:n, 0:p); close (2)
         end if
 
-        gamma = fluid_pp(1)%gamma
-        lit_gamma = 1._wp/fluid_pp(1)%gamma + 1._wp
-        pi_inf = fluid_pp(1)%pi_inf
-        qv = fluid_pp(1)%qv
+        gamma = gammas(1)
+        lit_gamma = gs_min(1)
+        pi_inf = pi_infs(1)
+        qv = qvs(1)
 
         if (precision == 1) then
             FMT = "(2F30.3)"
@@ -1239,7 +1239,7 @@ contains
                     end if
 
                     if (model_eqns == 4) then
-                        lit_gamma = 1._wp/fluid_pp(1)%gamma + 1._wp
+                        lit_gamma = gammas(1)
                     else if (elasticity) then
                         tau_e(1) = q_cons_vf(stress_idx%end)%sf(j - 2, k, l)/rho
                     end if
@@ -1353,7 +1353,7 @@ contains
                         end if
 
                         if (model_eqns == 4) then
-                            lit_gamma = 1._wp/fluid_pp(1)%gamma + 1._wp
+                            lit_gamma = gs_min(1)
                         else if (elasticity) then
                             do s = 1, 3
                                 tau_e(s) = q_cons_vf(s)%sf(j - 2, k - 2, l)/rho

src/simulation/m_global_parameters.fpp

@@ -943,8 +943,6 @@ contains
                         end if
                     end if
                     sys_size = bub_idx%end
-                    ! print*, 'alf idx', alf_idx
-                    ! print*, 'bub -idx beg end', bub_idx%beg, bub_idx%end
 
                     if (adv_n) then
                         n_idx = bub_idx%end + 1

src/simulation/m_pressure_relaxation.fpp

@@ -20,9 +20,6 @@ module m_pressure_relaxation
  s_initialize_pressure_relaxation_module, &
  s_finalize_pressure_relaxation_module
 
-    real(wp), allocatable, dimension(:) :: gamma_min, pres_inf
-    $:GPU_DECLARE(create='[gamma_min, pres_inf]')
-
     real(wp), allocatable, dimension(:, :) :: Res_pr
     $:GPU_DECLARE(create='[Res_pr]')
 
@@ -33,14 +30,6 @@ contains
 
         integer :: i, j
 
-        @:ALLOCATE(gamma_min(1:num_fluids), pres_inf(1:num_fluids))
-
-        do i = 1, num_fluids
-            gamma_min(i) = 1._wp/fluid_pp(i)%gamma + 1._wp
-            pres_inf(i) = fluid_pp(i)%pi_inf/(1._wp + fluid_pp(i)%gamma)
-        end do
-        $:GPU_UPDATE(device='[gamma_min, pres_inf]')
-
         if (viscous) then
             @:ALLOCATE(Res_pr(1:2, 1:Re_size_max))
             do i = 1, 2
@@ -56,7 +45,6 @@ contains
     !> Finalize the pressure relaxation module
     impure subroutine s_finalize_pressure_relaxation_module
 
-        @:DEALLOCATE(gamma_min, pres_inf)
         if (viscous) then
             @:DEALLOCATE(Res_pr)
         end if
@@ -170,8 +158,8 @@ contains
             if (q_cons_vf(i + advxb - 1)%sf(j, k, l) > sgm_eps) then
                 pres_K_init(i) = (q_cons_vf(i + intxb - 1)%sf(j, k, l)/ &
                                   q_cons_vf(i + advxb - 1)%sf(j, k, l) - pi_infs(i))/gammas(i)
-                if (pres_K_init(i) <= -(1._wp - 1.e-8_wp)*pres_inf(i) + 1.e-8_wp) &
-                    pres_K_init(i) = -(1._wp - 1.e-8_wp)*pres_inf(i) + 1.e-8_wp
+                if (pres_K_init(i) <= -(1._wp - 1.e-8_wp)*ps_inf(i) + 1.e-8_wp) &
+                    pres_K_init(i) = -(1._wp - 1.e-8_wp)*ps_inf(i) + 1.e-8_wp
             else
                 pres_K_init(i) = 0._wp
             end if
@@ -188,8 +176,8 @@ contains
 
                 ! Enforce pressure bounds
                 do i = 1, num_fluids
-                    if (pres_relax <= -(1._wp - 1.e-8_wp)*pres_inf(i) + 1.e-8_wp) &
-                        pres_relax = -(1._wp - 1.e-8_wp)*pres_inf(i) + 1._wp
+                    if (pres_relax <= -(1._wp - 1.e-8_wp)*ps_inf(i) + 1.e-8_wp) &
+                        pres_relax = -(1._wp - 1.e-8_wp)*ps_inf(i) + 1.e-8_wp
                 end do
 
                 ! Newton-Raphson step
@@ -200,11 +188,11 @@ contains
                     if (q_cons_vf(i + advxb - 1)%sf(j, k, l) > sgm_eps) then
                         rho_K_s(i) = q_cons_vf(i + contxb - 1)%sf(j, k, l)/ &
                                      max(q_cons_vf(i + advxb - 1)%sf(j, k, l), sgm_eps) &
-                                     *((pres_relax + pres_inf(i))/(pres_K_init(i) + &
-                                                                   pres_inf(i)))**(1._wp/gamma_min(i))
+                                     *((pres_relax + ps_inf(i))/(pres_K_init(i) + &
+                                                                 ps_inf(i)))**(1._wp/gs_min(i))
                         f_pres = f_pres + q_cons_vf(i + contxb - 1)%sf(j, k, l)/rho_K_s(i)
                         df_pres = df_pres - q_cons_vf(i + contxb - 1)%sf(j, k, l) &
-                                  /(gamma_min(i)*rho_K_s(i)*(pres_relax + pres_inf(i)))
+                                  /(gs_min(i)*rho_K_s(i)*(pres_relax + ps_inf(i)))
                     end if
                 end do
             end if

src/simulation/m_start_up.fpp

@@ -1052,9 +1052,8 @@ contains
                                             dyn_pres, pi_inf, gamma, rho, qv, rhoYks, pres, T, pres_mag=pres_mag)
 
                     do i = 1, num_fluids
-                        v_vf(i + internalEnergies_idx%beg - 1)%sf(j, k, l) = v_vf(i + adv_idx%beg - 1)%sf(j, k, l)* &
-                                                                             (fluid_pp(i)%gamma*pres + fluid_pp(i)%pi_inf) &
-                                                                             + v_vf(i + cont_idx%beg - 1)%sf(j, k, l)*fluid_pp(i)%qv
+                        v_vf(i + intxb - 1)%sf(j, k, l) = v_vf(i + advxb - 1)%sf(j, k, l)*(gammas(i)*pres + pi_infs(i)) &
+                                                          + v_vf(i + contxb - 1)%sf(j, k, l)*qvs(i)
                     end do
 
                 end do