From 841f863e32e6cfcfee6fc8c2e3ae7ec32977ef8a Mon Sep 17 00:00:00 2001 From: JRChreim Date: Tue, 25 Nov 2025 19:31:18 -0800 Subject: [PATCH 1/4] 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 --- src/common/m_variables_conversion.fpp | 13 +++++------ src/post_process/m_data_output.fpp | 10 +++++---- src/post_process/m_derived_variables.fpp | 8 +++---- src/pre_process/m_assign_variables.fpp | 10 ++++----- src/pre_process/m_data_output.fpp | 8 +++---- src/pre_process/m_icpp_patches.fpp | 28 +++++++++++++----------- src/simulation/m_data_output.fpp | 12 +++++----- src/simulation/m_global_parameters.fpp | 2 -- src/simulation/m_pressure_relaxation.fpp | 26 ++++++---------------- src/simulation/m_start_up.fpp | 5 ++--- 10 files changed, 54 insertions(+), 68 deletions(-) diff --git a/src/common/m_variables_conversion.fpp b/src/common/m_variables_conversion.fpp index f965f8ef3c..307406e6bd 100644 --- a/src/common/m_variables_conversion.fpp +++ b/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 diff --git a/src/post_process/m_data_output.fpp b/src/post_process/m_data_output.fpp index 41363136ca..410ab9ebcf 100644 --- a/src/post_process/m_data_output.fpp +++ b/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 diff --git a/src/post_process/m_derived_variables.fpp b/src/post_process/m_derived_variables.fpp index bc1242a97d..497178ccf3 100644 --- a/src/post_process/m_derived_variables.fpp +++ b/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 diff --git a/src/pre_process/m_assign_variables.fpp b/src/pre_process/m_assign_variables.fpp index 3eb3be17e6..8a9f5c2124 100644 --- a/src/pre_process/m_assign_variables.fpp +++ b/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) = & diff --git a/src/pre_process/m_data_output.fpp b/src/pre_process/m_data_output.fpp index fea0502496..ca5fb3c56a 100644 --- a/src/pre_process/m_data_output.fpp +++ b/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)" diff --git a/src/pre_process/m_icpp_patches.fpp b/src/pre_process/m_icpp_patches.fpp index 79fe1ac6e3..285f9ffd9e 100644 --- a/src/pre_process/m_icpp_patches.fpp +++ b/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 diff --git a/src/simulation/m_data_output.fpp b/src/simulation/m_data_output.fpp index 1c9f2e2933..2ddae9bc30 100644 --- a/src/simulation/m_data_output.fpp +++ b/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 diff --git a/src/simulation/m_global_parameters.fpp b/src/simulation/m_global_parameters.fpp index 31230ad40b..ab0d9a5c46 100644 --- a/src/simulation/m_global_parameters.fpp +++ b/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 diff --git a/src/simulation/m_pressure_relaxation.fpp b/src/simulation/m_pressure_relaxation.fpp index 4c6618dd95..a8452876f4 100644 --- a/src/simulation/m_pressure_relaxation.fpp +++ b/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 diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index 5e56dd1562..eb5fd3179e 100644 --- a/src/simulation/m_start_up.fpp +++ b/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 From 9e05105a3ebee74dad95af913b0572f33a6e7216 Mon Sep 17 00:00:00 2001 From: JRChreim Date: Tue, 25 Nov 2025 20:16:43 -0800 Subject: [PATCH 2/4] Update src/simulation/m_pressure_relaxation.fpp Co-authored-by: qodo-merge-pro[bot] <151058649+qodo-merge-pro[bot]@users.noreply.github.com> --- src/simulation/m_pressure_relaxation.fpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/simulation/m_pressure_relaxation.fpp b/src/simulation/m_pressure_relaxation.fpp index a8452876f4..c78fbe90ae 100644 --- a/src/simulation/m_pressure_relaxation.fpp +++ b/src/simulation/m_pressure_relaxation.fpp @@ -177,7 +177,7 @@ contains ! Enforce pressure bounds do i = 1, num_fluids 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 + pres_relax = -(1._wp - 1.e-8_wp)*ps_inf(i) + 1.e-8_wp end do ! Newton-Raphson step From 151504ea504e7b544c8e3816caf97e8a6d04dc82 Mon Sep 17 00:00:00 2001 From: JRChreim Date: Tue, 25 Nov 2025 20:21:52 -0800 Subject: [PATCH 3/4] Update m_assign_variables.fpp --- src/pre_process/m_assign_variables.fpp | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/src/pre_process/m_assign_variables.fpp b/src/pre_process/m_assign_variables.fpp index 8a9f5c2124..0dd0deab77 100644 --- a/src/pre_process/m_assign_variables.fpp +++ b/src/pre_process/m_assign_variables.fpp @@ -214,10 +214,7 @@ contains pres_mag = 1.e-1_wp loc = x_cc(177) 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 + 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) From 41db1d9e68d27a8c4ec4c868a90dccb23f1ceb2c Mon Sep 17 00:00:00 2001 From: JRChreim Date: Tue, 25 Nov 2025 20:25:21 -0800 Subject: [PATCH 4/4] Update src/pre_process/m_icpp_patches.fpp Co-authored-by: cubic-dev-ai[bot] <191113872+cubic-dev-ai[bot]@users.noreply.github.com> --- src/pre_process/m_icpp_patches.fpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pre_process/m_icpp_patches.fpp b/src/pre_process/m_icpp_patches.fpp index 285f9ffd9e..fc0afdeec6 100644 --- a/src/pre_process/m_icpp_patches.fpp +++ b/src/pre_process/m_icpp_patches.fpp @@ -799,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))* & - (1._wp - q_prim_vf(alf_idx)%sf(i, j, 0)) + rhoref*(1._wp - q_prim_vf(alf_idx)%sf(i, j, 0)) end if ! Updating the patch identities bookkeeping variable