Changes to through the code to make fluid variables uniform

The following PR comprises changes to the code in order to call fluid variables uniformly through the code. The main changes for this PR is to change all the fluid_pp(i)% to their respective, which are defined on src/common/m_variables_conversion.fpp. For instance, substituting fluid_pp(i)%pi_inf by pi_infs(i), and so on.

In certain cases, variable indices are also replaced to be uniform

The importance of this PR is that, if these variables change through the code for whatever reason, they will be followed through the rest of the code

Author: JRChreim

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,8 +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 = gs_min(1)
+        B_tait = pi_infs(1)
 
         n_tait = 1._wp/n_tait + 1._wp
         B_tait = B_tait*(n_tait - 1._wp)/n_tait
@@ -563,9 +563,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
@@ -797,7 +799,7 @@ contains
                             !zero density, reassign according to Tait EOS
                             q_prim_vf(1)%sf(i, j, 0) = &
                                 (((q_prim_vf(E_idx)%sf(i, j, 0) + pi_inf)/(pref + pi_inf))**(1._wp/lit_gamma))* &
-                                rhoref*(1._wp - q_prim_vf(alf_idx)%sf(i, j, 0))
+                                (1._wp - q_prim_vf(alf_idx)%sf(i, j, 0))
                         end if
 
                         ! Updating the patch identities bookkeeping variable
@@ -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._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