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..0dd0deab77 100644 --- a/src/pre_process/m_assign_variables.fpp +++ b/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) = & 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..fc0afdeec6 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 @@ -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..c78fbe90ae 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.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 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