diff --git a/CMakeLists.txt b/CMakeLists.txt index f61394f79..c1f228d77 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -168,11 +168,16 @@ macro(HANDLE_FYPP target) set(f90_filepath "${f90_dirpath}/autogen/${f90_filename}") + set(depends_on "${fpp_filepath}") + if (EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/src/${target}/case.fpp") + list(APPEND depends_on "${CMAKE_CURRENT_SOURCE_DIR}/src/${target}/case.fpp") + endif() + add_custom_command( OUTPUT "${f90_filepath}" COMMAND "${FYPP_EXE}" "${fpp_filepath}" "${f90_filepath}" - DEPENDS "${fpp_filepath}" "${CMAKE_CURRENT_SOURCE_DIR}/src/${target}/case.fpp" - COMMENT "Preprocessing ${fpp_filename}" + DEPENDS "${depends_on}" + COMMENT "Preprocessing (Fypp) ${fpp_filename}" VERBATIM ) diff --git a/src/simulation/m_variables_conversion.f90 b/src/common/m_variables_conversion.fpp similarity index 51% rename from src/simulation/m_variables_conversion.f90 rename to src/common/m_variables_conversion.fpp index 3cef9d119..60ad56977 100644 --- a/src/simulation/m_variables_conversion.f90 +++ b/src/common/m_variables_conversion.fpp @@ -2,15 +2,10 @@ !! @file m_variables_conversion.f90 !! @brief Contains module m_variables_conversion -!> @brief This module features a database of subroutines that allow for the -!! conversion of state variables from one type into another. At this -!! time, the state variables type conversions below are available: -!! 1) Mixture => Mixture -!! 2) Species => Mixture -!! 3) Conservative => Primitive -!! 5) Conservative => Flux -!! 6) Primitive => Conservative -!! 8) Primitive => Flux +!> @brief This module consists of subroutines used in the conversion of the +!! conservative variables into the primitive ones and vice versa. In +!! addition, the module also contains the subroutines used to obtain +!! the mixture variables and the subroutines used to compute pressure. module m_variables_conversion ! Dependencies ============================================================= @@ -19,8 +14,6 @@ module m_variables_conversion use m_global_parameters !< Definitions of the global parameters use m_mpi_proxy !< Message passing interface (MPI) module proxy - - use nvtx ! ========================================================================== implicit none @@ -35,41 +28,45 @@ module m_variables_conversion s_convert_conservative_to_primitive_variables, & s_convert_primitive_to_conservative_variables, & s_convert_primitive_to_flux_variables, & + s_compute_pressure, & s_finalize_variables_conversion_module - abstract interface ! ======================================================= + !> Abstract interface to two subroutines designed for the transfer/conversion + !! of the mixture/species variables to the mixture variables - !> The abstract interface to the procedures that are utilized to convert - !! the mixture and the species variables into the mixture variables. For - !! more information, refer to: - !! 1) s_convert_mixture_to_mixture_variables - !! 2) s_convert_species_to_mixture_variables - !! @param qK_vf Conservative/primitive variables - !! @param rho_K Mixture density - !! @param gamma_K Mixture sp. heat ratio - !! @param pi_inf_K Mixture stiffness function - !! @param Re_K Reynolds number - !! @param i Cell location first index - !! @param j Cell location second index - !! @param k Cell location third index - subroutine s_convert_abstract_to_mixture_variables(qK_vf, rho_K, & - gamma_K, pi_inf_K, & - Re_K, i, j, k, G_K, G) + abstract interface ! ======================================================= + !> Structure of the s_convert_mixture_to_mixture_variables + !! and s_convert_species_to_mixture_variables subroutines + !! @param q_vf Conservative or primitive variables + !! @param i First-coordinate cell index + !! @param j First-coordinate cell index + !! @param k First-coordinate cell index + !! @param rho Density + !! @param gamma Specific heat ratio function + !! @param pi_inf Liquid stiffness function + subroutine s_convert_xxxxx_to_mixture_variables(q_vf, i, j, k, & + rho, gamma, pi_inf, Re_K, G_K, G) + + ! Importing the derived type scalar_field from m_derived_types.f90 + ! and global variable sys_size, from m_global_variables.f90, as + ! the abstract interface does not inherently have access to them import :: scalar_field, sys_size, num_fluids - type(scalar_field), dimension(sys_size), intent(IN) :: qK_vf + type(scalar_field), dimension(sys_size), intent(IN) :: q_vf - real(kind(0d0)), intent(OUT) :: rho_K, gamma_K, pi_inf_K + integer, intent(IN) :: i, j, k - real(kind(0d0)), dimension(2), intent(OUT) :: Re_K + real(kind(0d0)), intent(OUT), target :: rho + real(kind(0d0)), intent(OUT), target :: gamma + real(kind(0d0)), intent(OUT), target :: pi_inf - integer, intent(IN) :: i, j, k + real(kind(0d0)), optional, dimension(2), intent(OUT) :: Re_K real(kind(0d0)), optional, intent(OUT) :: G_K real(kind(0d0)), optional, dimension(num_fluids), intent(IN) :: G - end subroutine s_convert_abstract_to_mixture_variables + end subroutine s_convert_xxxxx_to_mixture_variables !> The abstract interface to the procedures that are used to compute the !! Roe and the arithmetic average states. For additional information see: @@ -103,23 +100,27 @@ end subroutine s_compute_abstract_average_state !> @} integer :: ixb, ixe, iyb, iye, izb, ize - real(kind(0d0)), allocatable, dimension(:) :: Gs + !! In simulation, gammas and pi_infs is already declared in m_global_variables +#ifndef MFC_SIMULATION + real(kind(0d0)), allocatable, dimension(:) :: gammas, pi_infs +#endif + real(kind(0d0)), allocatable, dimension(:) :: Gs integer, allocatable, dimension(:) :: bubrs real(kind(0d0)), allocatable, dimension(:, :) :: Res -!$acc declare create(ixb, ixe, iyb, iye, izb, ize, bubrs, Gs, Res) +!$acc declare create(ixb, ixe, iyb, iye, izb, ize, momxb, momxe, bubxb, bubxe, contxb, contxe, advxb, advxe, strxb, strxe, gammas, pi_infs, bubrs, Gs, Res) integer :: is1b, is2b, is3b, is1e, is2e, is3e !$acc declare create(is1b, is2b, is3b, is1e, is2e, is3e) - ! Density, dynamic pressure, surface energy, specific heat ratio - ! function, liquid stiffness function, shear and volume Reynolds - ! numbers and the Weber numbers + real(kind(0d0)), allocatable, dimension(:, :, :), target, public :: rho_sf !< Scalar density function + real(kind(0d0)), allocatable, dimension(:, :, :), target, public :: gamma_sf !< Scalar sp. heat ratio function + real(kind(0d0)), allocatable, dimension(:, :, :), target, public :: pi_inf_sf !< Scalar liquid stiffness function - procedure(s_convert_abstract_to_mixture_variables), & + procedure(s_convert_xxxxx_to_mixture_variables), & pointer :: s_convert_to_mixture_variables => null() !< - !! Pointer to the procedure utilized to convert either the mixture or the - !! species variables into the mixture variables, based on model equations + !! Pointer referencing the subroutine s_convert_mixture_to_mixture_variables + !! or s_convert_species_to_mixture_variables, based on model equations choice procedure(s_compute_abstract_average_state), & pointer :: s_compute_average_state => null() !< @@ -128,39 +129,154 @@ end subroutine s_compute_abstract_average_state contains - !> This procedure is used alongside with the gamma/pi_inf - !! model to transfer the density, the specific heat ratio - !! function and liquid stiffness function from the vector - !! of conservative or primitive variables to their scalar - !! counterparts. - !! @param qK_vf conservative or primitive variables + !> This procedure calculates the pressure based on energy when + !! there are no bubbles present and model_eqns != 4 + !! @param energy Energy + !! @param alf Void Fraction + !! @param dyn_p Dynamic Pressure + !! @param pi_inf Liquid Stiffness + !! @param gamma Specific Heat Ratio + !! @param pres Pressure to calculate + subroutine s_compute_pressure_from_energy(energy, alf, dyn_p, pi_inf, gamma, pres) +!$acc routine seq + + real(kind(0d0)) :: energy, alf + + real(kind(0d0)), intent(IN) :: dyn_p + real(kind(0d0)), intent(OUT) :: pres + + real(kind(0d0)) :: pi_inf, gamma + + pres = (energy - dyn_p - pi_inf)/gamma + + end subroutine s_compute_pressure_from_energy + + !> This procedure calculates the pressure when there + !! are bubbles present and model_eqns != 4 + !! @param energy Energy + !! @param alf Void Fraction + !! @param dyn_p Dynamic Pressure + !! @param pi_inf Liquid Stiffness + !! @param gamma Specific Heat Ratio + !! @param pres Pressure to calculate + subroutine s_compute_pressure_from_bubbles(energy, alf, dyn_p, pi_inf, gamma, pres) +!$acc routine seq + + real(kind(0d0)) :: energy, alf + + real(kind(0d0)), intent(IN) :: dyn_p + real(kind(0d0)), intent(OUT) :: pres + + real(kind(0d0)) :: pi_inf, gamma + + pres = ((energy - dyn_p)/(1.d0 - alf) - pi_inf)/gamma + + end subroutine s_compute_pressure_from_bubbles + + !> This procedure calculates the pressure when model_eqns = 4 + !! @param energy Energy + !! @param alf Void Fraction + !! @param dyn_p Dynamic Pressure + !! @param pi_inf Liquid Stiffness + !! @param gamma Specific Heat Ratio + !! @param pres Pressure to calculate + subroutine s_compute_pressure_4eqns(energy, alf, dyn_p, pi_inf, gamma, pres) +!$acc routine seq + + real(kind(0d0)) :: energy, alf + + real(kind(0d0)), intent(IN) :: dyn_p + real(kind(0d0)), intent(OUT) :: pres + + real(kind(0d0)) :: pi_inf, gamma + + pres = (pref + pi_inf)* & + (energy/ & + (rhoref*(1 - alf)) & + )**(1/gamma + 1) - pi_inf + + end subroutine s_compute_pressure_4eqns + + !> This procedure conditionally calls the appropriate pressure-computing + !! subroutine. + !! @param energy Energy + !! @param alf Void Fraction + !! @param dyn_p Dynamic Pressure + !! @param pi_inf Liquid Stiffness + !! @param gamma Specific Heat Ratio + !! @param pres Pressure to calculate + subroutine s_compute_pressure(energy, alf, dyn_p, pi_inf, gamma, pres) +!$acc routine seq + + real(kind(0d0)) :: energy, alf + + real(kind(0d0)), intent(IN) :: dyn_p + real(kind(0d0)), intent(OUT) :: pres + + real(kind(0d0)) :: pi_inf, gamma + + ! Depending on model_eqns and bubbles, the appropriate procedure + ! for computing pressure is targeted by the procedure pointer + + if ((model_eqns /= 4) .and. (bubbles .neqv. .true.)) then + call s_compute_pressure_from_energy(energy, alf, dyn_p, pi_inf, gamma, pres) + else if ((model_eqns /= 4) .and. bubbles) then + call s_compute_pressure_from_bubbles(energy, alf, dyn_p, pi_inf, gamma, pres) + else + call s_compute_pressure_4eqns(energy, alf, dyn_p, pi_inf, gamma, pres) + end if + + end subroutine s_compute_pressure + + !> This subroutine is designed for the gamma/pi_inf model + !! and provided a set of either conservative or primitive + !! variables, transfers the density, specific heat ratio + !! function and the liquid stiffness function from q_vf to + !! rho, gamma and pi_inf. + !! @param q_vf conservative or primitive variables !! @param i cell index to transfer mixture variables !! @param j cell index to transfer mixture variables !! @param k cell index to transfer mixture variables - !! @param rho_K density - !! @param gamma_K specific heat ratio function - !! @param pi_inf_K liquid stiffness - !! @param Re_k Reynolds number - subroutine s_convert_mixture_to_mixture_variables(qK_vf, rho_K, & - gamma_K, pi_inf_K, & - Re_K, i, j, k, G_K, G) + !! @param rho density + !! @param gamma specific heat ratio function + !! @param pi_inf liquid stiffness + subroutine s_convert_mixture_to_mixture_variables(q_vf, i, j, k, & + rho, gamma, pi_inf, Re_K, G_K, G) - type(scalar_field), dimension(sys_size), intent(IN) :: qK_vf + type(scalar_field), dimension(sys_size), intent(IN) :: q_vf - real(kind(0d0)), intent(OUT) :: rho_K, gamma_K, pi_inf_K + integer, intent(IN) :: i, j, k - real(kind(0d0)), dimension(2), intent(OUT) :: Re_K + real(kind(0d0)), intent(OUT), target :: rho + real(kind(0d0)), intent(OUT), target :: gamma + real(kind(0d0)), intent(OUT), target :: pi_inf - integer, intent(IN) :: i, j, k + real(kind(0d0)), optional, dimension(2), intent(OUT) :: Re_K real(kind(0d0)), optional, intent(OUT) :: G_K real(kind(0d0)), optional, dimension(num_fluids), intent(IN) :: G - ! Performing the transfer of the density, the specific heat ratio - ! function as well as the liquid stiffness function, respectively - rho_K = qK_vf(1)%sf(i, j, k) - gamma_K = qK_vf(gamma_idx)%sf(i, j, k) - pi_inf_K = qK_vf(pi_inf_idx)%sf(i, j, k) + real(kind(0d0)), pointer :: rho_K, gamma_K, pi_inf_K + + !> Post process requires rho_sf/gamma_sf/pi_inf_sf to be + !! updated instead of rho/gamma/pi_inf. Therefore, the + !! versions of these variables appended with '_K' represent + !! pointers that target the correct variable. +#ifdef MFC_POST_PROCESS + rho_K => rho_sf(i, j, k) + gamma_K => gamma_sf(i, j, k) + pi_inf_K => pi_inf_sf(i, j, k) +#else + rho_K => rho + gamma_K => gamma + pi_inf_K => pi_inf +#endif + + ! Transfering the density, the specific heat ratio function and the + ! liquid stiffness function, respectively + rho_K = q_vf(1)%sf(i, j, k) + gamma_K = q_vf(gamma_idx)%sf(i, j, k) + pi_inf_K = q_vf(pi_inf_idx)%sf(i, j, k) end subroutine s_convert_mixture_to_mixture_variables ! ---------------- @@ -174,29 +290,43 @@ end subroutine s_convert_mixture_to_mixture_variables ! ---------------- !! @param rho_K density !! @param gamma_K specific heat ratio !! @param pi_inf_K liquid stiffness - !! @param Re_K mixture Reynolds number - !! @param i Cell index !! @param j Cell index !! @param k Cell index - subroutine s_convert_species_to_mixture_variables_bubbles(qK_vf, rho_K, & - gamma_K, pi_inf_K, & - Re_K, i, j, k, G_K, G) - - type(scalar_field), dimension(sys_size), intent(IN) :: qK_vf + !! @param l Cell index + subroutine s_convert_species_to_mixture_variables_bubbles(q_vf, j, k, l, & + rho, gamma, pi_inf, Re_K, G_K, G) - real(kind(0d0)), intent(OUT) :: rho_K, gamma_K, pi_inf_K + type(scalar_field), dimension(sys_size), intent(IN) :: q_vf - real(kind(0d0)), dimension(2), intent(OUT) :: Re_K + integer, intent(IN) :: j, k, l - real(kind(0d0)), dimension(num_fluids) :: alpha_rho_K, alpha_K !< - !! Partial densities and volume fractions + real(kind(0d0)), intent(OUT), target :: rho + real(kind(0d0)), intent(OUT), target :: gamma + real(kind(0d0)), intent(OUT), target :: pi_inf - integer, intent(IN) :: i, j, k - integer :: l + real(kind(0d0)), optional, dimension(2), intent(OUT) :: Re_K real(kind(0d0)), optional, intent(OUT) :: G_K real(kind(0d0)), optional, dimension(num_fluids), intent(IN) :: G + integer :: i + + real(kind(0d0)), pointer :: rho_K, gamma_K, pi_inf_K + + !> Post process requires rho_sf/gamma_sf/pi_inf_sf to be + !! updated instead of rho/gamma/pi_inf. Therefore, the + !! versions of these variables appended with '_K' represent + !! pointers that target the correct variable. +#ifdef MFC_POST_PROCESS + rho_K => rho_sf(j, k, l) + gamma_K => gamma_sf(j, k, l) + pi_inf_K => pi_inf_sf(j, k, l) +#else + rho_K => rho + gamma_K => gamma + pi_inf_K => pi_inf +#endif + ! Constraining the partial densities and the volume fractions within ! their physical bounds to make sure that any mixture variables that ! are derived from them result within the limits that are set by the @@ -208,43 +338,39 @@ subroutine s_convert_species_to_mixture_variables_bubbles(qK_vf, rho_K, & ! function as well as the liquid stiffness function, respectively if (model_eqns == 4) then - rho_K = qK_vf(1)%sf(i, j, k) - gamma_K = gammas(1)!qK_vf(gamma_idx)%sf(i,j,k) - pi_inf_K = pi_infs(1) !qK_vf(pi_inf_idx)%sf(i,j,k) + rho_K = q_vf(1)%sf(j, k, l) + gamma_K = fluid_pp(1)%gamma !qK_vf(gamma_idx)%sf(i,j,k) + pi_inf_K = fluid_pp(1)%pi_inf !qK_vf(pi_inf_idx)%sf(i,j,k) else if ((model_eqns == 2) .and. bubbles) then - rho_k = 0d0; gamma_k = 0d0; pi_inf_k = 0d0 + rho_K = 0d0; gamma_K = 0d0; pi_inf_K = 0d0 if (mpp_lim .and. (num_fluids > 2)) then - do l = 1, num_fluids - rho_k = rho_k + qK_vf(l)%sf(i, j, k) - gamma_k = gamma_k + qK_vf(l + E_idx)%sf(i, j, k)*gammas(l) - pi_inf_k = pi_inf_k + qK_vf(l + E_idx)%sf(i, j, k)*pi_infs(l) + do i = 1, num_fluids + rho_K = rho_K + q_vf(i)%sf(j, k, l) + gamma_K = gamma_K + q_vf(i + E_idx)%sf(j, k, l)*fluid_pp(i)%gamma + pi_inf_K = pi_inf_K + q_vf(i + E_idx)%sf(j, k, l)*fluid_pp(i)%pi_inf end do else if (num_fluids == 2) then - rho_K = qK_vf(1)%sf(i, j, k) - gamma_K = gammas(1) - pi_inf_K = pi_infs(1) + rho_K = q_vf(1)%sf(j, k, l) + gamma_K = fluid_pp(1)%gamma + pi_inf_K = fluid_pp(1)%pi_inf else if (num_fluids > 2) then - do l = 1, num_fluids - 1 !leave out bubble part of mixture - rho_k = rho_k + qK_vf(l)%sf(i, j, k) - gamma_k = gamma_k + qK_vf(l + E_idx)%sf(i, j, k)*gammas(l) - pi_inf_k = pi_inf_k + qK_vf(l + E_idx)%sf(i, j, k)*pi_infs(l) + !TODO: This may need fixing for hypo + bubbles + do i = 1, num_fluids - 1 !leave out bubble part of mixture + rho_K = rho_K + q_vf(i)%sf(j, k, l) + gamma_K = gamma_K + q_vf(i + E_idx)%sf(j, k, l)*fluid_pp(i)%gamma + pi_inf_K = pi_inf_K + q_vf(i + E_idx)%sf(j, k, l)*fluid_pp(i)%pi_inf end do + !rho_K = qK_vf(1)%sf(j,k,l) + !gamma_K = fluid_pp(1)%gamma + !pi_inf_K = fluid_pp(1)%pi_inf else - rho_K = qK_vf(1)%sf(i, j, k) - gamma_K = gammas(1) - pi_inf_K = pi_infs(1) + rho_K = q_vf(1)%sf(j, k, l) + gamma_K = fluid_pp(1)%gamma + pi_inf_K = fluid_pp(1)%pi_inf end if end if - if (present(G_K)) then - G_K = 0d0 - do l = 1, num_fluids - G_K = G_K + alpha_K(l)*G(l) - end do - G_K = max(0d0, G_K) - end if - end subroutine s_convert_species_to_mixture_variables_bubbles ! ---------------- !> This subroutine is designed for the volume fraction model @@ -252,23 +378,25 @@ end subroutine s_convert_species_to_mixture_variables_bubbles ! ---------------- !! variables, computes the density, the specific heat ratio !! function and the liquid stiffness function from q_vf and !! stores the results into rho, gamma and pi_inf. - !! @param qK_vf primitive variables - !! @param rho_K density - !! @param gamma_K specific heat ratio - !! @param pi_inf_K liquid stiffness - !! @param Re_K mixture Reynolds number + !! @param q_vf primitive variables + !! @param rho density + !! @param gamma specific heat ratio + !! @param pi_inf liquid stiffness + !! @param j Cell index !! @param k Cell index !! @param l Cell index - !! @param r Cell index - subroutine s_convert_species_to_mixture_variables(qK_vf, rho_K, & - gamma_K, pi_inf_K, & - Re_K, k, l, r, G_K, G) + subroutine s_convert_species_to_mixture_variables(q_vf, k, l, r, & + rho, gamma, pi_inf, Re_K, G_K, G) - type(scalar_field), dimension(sys_size), intent(IN) :: qK_vf + type(scalar_field), dimension(sys_size), intent(IN) :: q_vf - real(kind(0d0)), intent(OUT) :: rho_K, gamma_K, pi_inf_K + integer, intent(IN) :: k, l, r - real(kind(0d0)), dimension(2), intent(OUT) :: Re_K + real(kind(0d0)), intent(OUT), target :: rho + real(kind(0d0)), intent(OUT), target :: gamma + real(kind(0d0)), intent(OUT), target :: pi_inf + + real(kind(0d0)), optional, dimension(2), intent(OUT) :: Re_K real(kind(0d0)), dimension(num_fluids) :: alpha_rho_K, alpha_K !< !! Partial densities and volume fractions @@ -276,18 +404,30 @@ subroutine s_convert_species_to_mixture_variables(qK_vf, rho_K, & real(kind(0d0)), optional, intent(OUT) :: G_K real(kind(0d0)), optional, dimension(num_fluids), intent(IN) :: G - integer, intent(IN) :: k, l, r + integer :: i, j !< Generic loop iterator - integer :: i, j !< Generic loop iterators + real(kind(0d0)), pointer :: rho_K, gamma_K, pi_inf_K - ! Constraining the partial densities and the volume fractions within - ! their physical bounds to make sure that any mixture variables that - ! are derived from them result within the limits that are set by the - ! fluids physical parameters that make up the mixture + !> Post process requires rho_sf/gamma_sf/pi_inf_sf to be + !! updated instead of rho/gamma/pi_inf. Therefore, the + !! versions of these variables appended with '_K' represent + !! pointers that target the correct variable. +#ifdef MFC_POST_PROCESS + rho_K => rho_sf(k, l, r) + gamma_K => gamma_sf(k, l, r) + pi_inf_K => pi_inf_sf(k, l, r) +#else + rho_K => rho + gamma_K => gamma + pi_inf_K => pi_inf +#endif + + ! Computing the density, the specific heat ratio function and the + ! liquid stiffness function, respectively do i = 1, num_fluids - alpha_rho_K(i) = qK_vf(i)%sf(k, l, r) - alpha_K(i) = qK_vf(advxb + i - 1)%sf(k, l, r) + alpha_rho_K(i) = q_vf(i)%sf(k, l, r) + alpha_K(i) = q_vf(advxb + i - 1)%sf(k, l, r) end do if (mpp_lim) then @@ -311,6 +451,7 @@ subroutine s_convert_species_to_mixture_variables(qK_vf, rho_K, & pi_inf_K = pi_inf_K + alpha_K(i)*pi_infs(i) end do +#ifdef MFC_SIMULATION ! Computing the shear and bulk Reynolds numbers from species analogs do i = 1, 2 @@ -324,6 +465,7 @@ subroutine s_convert_species_to_mixture_variables(qK_vf, rho_K, & Re_K(i) = 1d0/max(Re_K(i), sgm_eps) end do +#endif if (present(G_K)) then G_K = 0d0 @@ -355,6 +497,7 @@ subroutine s_convert_species_to_mixture_variables_acc(rho_K, & integer :: i, j !< Generic loop iterators real(kind(0d0)) :: alpha_K_sum +#ifdef MFC_SIMULATION ! Constraining the partial densities and the volume fractions within ! their physical bounds to make sure that any mixture variables that ! are derived from them result within the limits that are set by the @@ -407,6 +550,7 @@ subroutine s_convert_species_to_mixture_variables_acc(rho_K, & end do end if +#endif end subroutine s_convert_species_to_mixture_variables_acc ! ---------------- @@ -422,6 +566,7 @@ subroutine s_convert_species_to_mixture_variables_bubbles_acc(rho_K, & integer, intent(IN) :: k, l, r integer :: i, j !< Generic loop iterators +#ifdef MFC_SIMULATION rho_K = 0d0 gamma_K = 0d0 pi_inf_K = 0d0 @@ -443,6 +588,7 @@ subroutine s_convert_species_to_mixture_variables_bubbles_acc(rho_K, & gamma_K = gammas(1) pi_inf_K = pi_infs(1) end if +#endif end subroutine s_convert_species_to_mixture_variables_bubbles_acc @@ -452,7 +598,12 @@ end subroutine s_convert_species_to_mixture_variables_bubbles_acc subroutine s_initialize_variables_conversion_module() ! ---------------- integer :: i, j +!$acc update device(momxb, momxe, bubxb, bubxe, advxb, advxe, contxb, contxe, strxb, strxe) +#ifdef MFC_PRE_PROCESS + ixb = 0; iyb = 0; izb = 0; + ixe = m; iye = n; ize = p; +#else ixb = -buff_size ixe = m - ixb @@ -466,25 +617,32 @@ subroutine s_initialize_variables_conversion_module() ! ---------------- izb = -buff_size ize = p - izb end if +#endif !$acc update device(ixb, ixe, iyb, iye, izb, ize) + allocate (gammas(1:num_fluids)) + allocate (pi_infs(1:num_fluids)) allocate (Gs(1:num_fluids)) - if (any(Re_size > 0)) then - allocate (Res(1:2, 1:maxval(Re_size))) - end if - allocate (bubrs(1:nb)) do i = 1, num_fluids + gammas(i) = fluid_pp(i)%gamma + pi_infs(i) = fluid_pp(i)%pi_inf Gs(i) = fluid_pp(i)%G end do -!$acc update device(Gs) +!$acc update device(gammas, pi_infs, Gs) !TODO: this update was in previous version: (no longer needed?) !!!$acc update device(small_alf, dflt_real, dflt_int, pi, nnode, sgm_eps) +#ifdef MFC_SIMULATION + + if (any(Re_size > 0)) then + allocate (Res(1:2, 1:maxval(Re_size))) + end if + if (any(Re_size > 0)) then do i = 1, 2 do j = 1, Re_size(i) @@ -493,6 +651,7 @@ subroutine s_initialize_variables_conversion_module() ! ---------------- end do !$acc update device(Res, Re_idx, Re_size) end if +#endif if (bubbles) then @@ -504,7 +663,7 @@ subroutine s_initialize_variables_conversion_module() ! ---------------- !$acc update device(bubrs) !$acc update device(dt, sys_size, pref, rhoref, gamma_idx, pi_inf_idx, E_idx, alf_idx, stress_idx, mpp_lim, bubbles, hypoelasticity, alt_soundspeed, avg_state, num_fluids, model_eqns, num_dims, mixture_err, nb, weight, grid_geometry, cyl_coord, mapped_weno, mp_weno, weno_eps) -!$acc update device(nb, R0ref, Ca, Web, Re_inv, weight, R0, V0, bubbles, polytropic, polydisperse, qbmm, nmom, nmomsp, nmomtot, R0_type, ptil, bubble_model, thermal, poly_sigma) +!$acc update device(nb, R0ref, Ca, Web, Re_inv, weight, R0, V0, bubbles, polytropic, polydisperse, qbmm, R0_type, ptil, bubble_model, thermal, poly_sigma) !$acc update device(R_n, R_v, phi_vn, phi_nv, Pe_c, Tw, pv, M_n, M_v, k_n, k_v, pb0, mass_n0, mass_v0, Pe_T, Re_trans_T, Re_trans_c, Im_trans_T, Im_trans_c, omegaN , mul0, ss, gamma_v, mu_v, gamma_m, gamma_n, mu_n, gam) @@ -522,15 +681,66 @@ subroutine s_initialize_variables_conversion_module() ! ---------------- !$acc update device(monopole, num_mono) - ! Associating the procedural pointer to the appropriate subroutine - ! that will be utilized in the conversion to the mixture variables - if (model_eqns == 1) then ! gamma/pi_inf model +#ifdef MFC_POST_PROCESS + ! Allocating the density, the specific heat ratio function and the + ! liquid stiffness function, respectively + + ! Simulation is at least 2D + if (n > 0) then + + ! Simulation is 3D + if (p > 0) then + + allocate (rho_sf(-buff_size:m + buff_size, & + -buff_size:n + buff_size, & + -buff_size:p + buff_size)) + allocate (gamma_sf(-buff_size:m + buff_size, & + -buff_size:n + buff_size, & + -buff_size:p + buff_size)) + allocate (pi_inf_sf(-buff_size:m + buff_size, & + -buff_size:n + buff_size, & + -buff_size:p + buff_size)) + + ! Simulation is 2D + else + + allocate (rho_sf(-buff_size:m + buff_size, & + -buff_size:n + buff_size, & + 0:0)) + allocate (gamma_sf(-buff_size:m + buff_size, & + -buff_size:n + buff_size, & + 0:0)) + allocate (pi_inf_sf(-buff_size:m + buff_size, & + -buff_size:n + buff_size, & + 0:0)) + + end if + + ! Simulation is 1D + else + + allocate (rho_sf(-buff_size:m + buff_size, & + 0:0, & + 0:0)) + allocate (gamma_sf(-buff_size:m + buff_size, & + 0:0, & + 0:0)) + allocate (pi_inf_sf(-buff_size:m + buff_size, & + 0:0, & + 0:0)) + + end if +#endif + + if (model_eqns == 1) then ! Gamma/pi_inf model s_convert_to_mixture_variables => & s_convert_mixture_to_mixture_variables - elseif (bubbles) then ! Volume fraction for bubbles + + else if (bubbles) then s_convert_to_mixture_variables => & s_convert_species_to_mixture_variables_bubbles - else ! Volume fraction model + else + ! Volume fraction model s_convert_to_mixture_variables => & s_convert_species_to_mixture_variables end if @@ -559,105 +769,88 @@ subroutine s_convert_conservative_to_primitive_variables(qK_cons_vf, & intent(INOUT) :: qK_prim_vf type(scalar_field), & - allocatable, dimension(:), & + allocatable, optional, dimension(:), & intent(IN) :: gm_alphaK_vf - type(int_bounds_info), intent(IN) :: ix, iy, iz + type(int_bounds_info), optional, intent(IN) :: ix, iy, iz real(kind(0d0)), dimension(num_fluids) :: alpha_K, alpha_rho_K real(kind(0d0)), dimension(2) :: Re_K real(kind(0d0)) :: rho_K, gamma_K, pi_inf_K, dyn_pres_K - real(kind(0d0)), dimension(nb) :: nRtmp + + real(kind(0d0)), dimension(:), allocatable :: nRtmp real(kind(0d0)) :: vftmp, nR3, nbub_sc real(kind(0d0)) :: G_K - integer :: i, j, k, l !< Generic loop iterators - - - if ((model_eqns /= 4) .and. (bubbles .neqv. .true.)) then -!$acc parallel loop collapse(3) gang vector default(present) private( alpha_K, alpha_rho_K, Re_K) - do l = izb, ize - do k = iyb, iye - do j = ixb, ixe - dyn_pres_K = 0d0 + real(kind(0d0)) :: pres + integer :: i, j, k, l !< Generic loop iterators + + if (bubbles) then + allocate(nRtmp(nb)) + else + allocate(nRtmp(0)) + endif + +!$acc parallel loop collapse(3) gang vector default(present) private(alpha_K, alpha_rho_K, Re_K, nRtmp, rho_K, gamma_K, pi_inf_K, dyn_pres_K) + do l = izb, ize + do k = iyb, iye + do j = ixb, ixe + dyn_pres_K = 0d0 !$acc loop seq - do i = 1, num_fluids - alpha_rho_K(i) = qK_cons_vf(i)%sf(j, k, l) - alpha_K(i) = qK_cons_vf(advxb + i - 1)%sf(j, k, l) - end do + do i = 1, num_fluids + alpha_rho_K(i) = qK_cons_vf(i)%sf(j, k, l) + alpha_K(i) = qK_cons_vf(advxb + i - 1)%sf(j, k, l) + end do + if (model_eqns /= 4) then +#ifdef MFC_SIMULATION + ! If in simulation, use acc mixture subroutines if (hypoelasticity) then call s_convert_species_to_mixture_variables_acc(rho_K, gamma_K, pi_inf_K, alpha_K, & alpha_rho_K, Re_K, j, k, l, G_K, Gs) + else if (bubbles) then + call s_convert_species_to_mixture_variables_bubbles_acc(rho_K, gamma_K, pi_inf_K, & + alpha_K, alpha_rho_K, j, k, l) + else + call s_convert_species_to_mixture_variables_acc(rho_K, gamma_K, pi_inf_K, & + alpha_K, alpha_rho_K, Re_K, j, k, l) + end if +#else + ! If pre-processing, use non acc mixture subroutines + if (hypoelasticity) then + call s_convert_to_mixture_variables(qK_cons_vf, j, k, l, & + rho_K, gamma_K, pi_inf_K, Re_K, G_K) else - call s_convert_species_to_mixture_variables_acc(rho_K, gamma_K, pi_inf_K, alpha_K, alpha_rho_K, & - Re_K, j, k, l) + call s_convert_to_mixture_variables(qK_cons_vf, j, k, l, & + rho_K, gamma_K, pi_inf_K) end if +#endif + end if +#ifdef MFC_SIMULATION + rho_K = max(rho_K, sgm_eps) +#endif !$acc loop seq - do i = momxb, momxe + do i = momxb, momxe + if (model_eqns /= 4) then qK_prim_vf(i)%sf(j, k, l) = qK_cons_vf(i)%sf(j, k, l) & - /max(rho_K, sgm_eps) + /rho_K dyn_pres_K = dyn_pres_K + 5d-1*qK_cons_vf(i)%sf(j, k, l) & *qK_prim_vf(i)%sf(j, k, l) - end do - - qK_prim_vf(E_idx)%sf(j, k, l) = (qK_cons_vf(E_idx)%sf(j, k, l) & - - dyn_pres_K - pi_inf_K)/gamma_K - - if (hypoelasticity) then -!$acc loop seq - do i = strxb, strxe - qK_prim_vf(i)%sf(j, k, l) = qK_cons_vf(i)%sf(j, k, l) & - /max(rho_K, sgm_eps) - ! subtracting elastic contribution for pressure calculation - if (G_K > 1000) then !TODO: check if stable for >0 - qK_prim_vf(E_idx)%sf(j, k, l) = qK_prim_vf(E_idx)%sf(j, k, l) - & - ((qK_prim_vf(i)%sf(j, k, l)**2d0)/(4d0*G_K))/gamma_K - ! extra terms in 2 and 3D - if ((i == strxb + 1) .or. & - (i == strxb + 3) .or. & - (i == strxb + 4)) then - qK_prim_vf(E_idx)%sf(j, k, l) = qK_prim_vf(E_idx)%sf(j, k, l) - & - ((qK_prim_vf(i)%sf(j, k, l)**2d0)/(4d0*G_K))/gamma_K - end if - end if - end do + else + 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 - end do - end do - end do -!$acc end parallel loop - - elseif ((model_eqns /= 4)) then !TODO: add routine below for bubble + hypo -!$acc parallel loop collapse(3) gang vector default(present) private(alpha_rho_K, alpha_K, nRtmp) - do l = izb, ize - do k = iyb, iye - do j = ixb, ixe - dyn_pres_K = 0d0 -!$acc loop seq - do i = 1, num_fluids - alpha_rho_K(i) = qK_cons_vf(i)%sf(j, k, l) - alpha_K(i) = qK_cons_vf(advxb + i - 1)%sf(j, k, l) - end do - - call s_convert_species_to_mixture_variables_bubbles_acc(rho_K, gamma_K, & - pi_inf_K, alpha_K, alpha_rho_K, j, k, l) - -!$acc loop seq - do i = momxb, momxe - qK_prim_vf(i)%sf(j, k, l) = qK_cons_vf(i)%sf(j, k, l) & - /max(rho_K, sgm_eps) - dyn_pres_K = dyn_pres_K + 5d-1*qK_cons_vf(i)%sf(j, k, l) & - *qK_prim_vf(i)%sf(j, k, l) - end do - qK_prim_vf(E_idx)%sf(j, k, l) = (((qK_cons_vf(E_idx)%sf(j, k, l) - dyn_pres_K)/ & - (1.d0 - qK_cons_vf(alf_idx)%sf(j, k, l))) - pi_inf_K)/gamma_K + call s_compute_pressure(qK_cons_vf(E_idx)%sf(j, k, l), & + qK_cons_vf(alf_idx)%sf(j, k, l), & + dyn_pres_K, pi_inf_K, gamma_K, pres) + qK_prim_vf(E_idx)%sf(j, k, l) = pres + if (bubbles) then !$acc loop seq do i = 1, nb nRtmp(i) = qK_cons_vf(bubrs(i))%sf(j, k, l) @@ -666,36 +859,35 @@ subroutine s_convert_conservative_to_primitive_variables(qK_cons_vf, & vftmp = qK_cons_vf(alf_idx)%sf(j, k, l) call s_comp_n_from_cons(vftmp, nRtmp, nbub_sc) - !$acc loop seq do i = bubxb, bubxe qK_prim_vf(i)%sf(j, k, l) = qK_cons_vf(i)%sf(j, k, l)/nbub_sc end do + end if - end do - end do - end do -!$acc end parallel loop - else -!$acc parallel loop collapse(3) gang vector default(present) private(rho_K, gamma_K, pi_inf_K, dyn_pres_K) - do l = izb, ize - do k = iyb, iye - do j = ixb, ixe - dyn_pres_K = 0d0 -!$acc loop - do i = momxb, momxe + if (hypoelasticity) then +!$acc loop seq + do i = strxb, strxe qK_prim_vf(i)%sf(j, k, l) = qK_cons_vf(i)%sf(j, k, l) & - /qK_cons_vf(1)%sf(j, k, l) + /rho_K + ! subtracting elastic contribution for pressure calculation + if (G_K > 1000) then !TODO: check if stable for >0 + qK_prim_vf(E_idx)%sf(j, k, l) = qK_prim_vf(E_idx)%sf(j, k, l) - & + ((qK_prim_vf(i)%sf(j, k, l)**2d0)/(4d0*G_K))/gamma_K + ! extra terms in 2 and 3D + if ((i == strxb + 1) .or. & + (i == strxb + 3) .or. & + (i == strxb + 4)) then + qK_prim_vf(E_idx)%sf(j, k, l) = qK_prim_vf(E_idx)%sf(j, k, l) - & + ((qK_prim_vf(i)%sf(j, k, l)**2d0)/(4d0*G_K))/gamma_K + end if + end if end do -!$acc end loop - qK_prim_vf(E_idx)%sf(j, k, l) = (pref + pi_infs(1))* & - ((qK_cons_vf(1)%sf(j, k, l)/(rhoref*(1.d0 - qK_cons_vf(E_idx + 1)%sf(j, k, l))) & - )**(1.d0/gammas(1) + 1.d0)) - pi_infs(1) - end do + end if end do end do + end do !$acc end parallel loop - end if end subroutine s_convert_conservative_to_primitive_variables ! --------- @@ -707,43 +899,139 @@ end subroutine s_convert_conservative_to_primitive_variables ! --------- !! @param ix Index bounds in the first coordinate direction !! @param iy Index bounds in the second coordinate direction !! @param iz Index bounds in the third coordinate direction - subroutine s_convert_primitive_to_conservative_variables(qK_prim_vf, & - qK_cons_vf, & - gm_alphaK_vf, & - ix, iy, iz) + subroutine s_convert_primitive_to_conservative_variables(q_prim_vf, & + q_cons_vf) type(scalar_field), & dimension(sys_size), & - intent(IN) :: qK_prim_vf + intent(IN) :: q_prim_vf type(scalar_field), & dimension(sys_size), & - intent(INOUT) :: qK_cons_vf + intent(INOUT) :: q_cons_vf + + ! Density, specific heat ratio function, liquid stiffness function + ! and dynamic pressure, as defined in the incompressible flow sense, + ! respectively + real(kind(0d0)) :: rho + real(kind(0d0)) :: gamma + real(kind(0d0)) :: pi_inf + real(kind(0d0)) :: dyn_pres + real(kind(0d0)) :: nbub, R3, vftmp + real(kind(0d0)), dimension(nb) :: Rtmp + + real(kind(0d0)) :: G + + integer :: i, j, k, l, q !< Generic loop iterators + +#ifndef MFC_SIMULATION + ! Converting the primitive variables to the conservative variables + do l = 0, p + do k = 0, n + do j = 0, m + + ! Obtaining the density, specific heat ratio function + ! and the liquid stiffness function, respectively + call s_convert_to_mixture_variables(q_prim_vf, j, k, l, & + rho, gamma, pi_inf) + + ! Transferring the continuity equation(s) variable(s) + do i = 1, cont_idx%end + q_cons_vf(i)%sf(j, k, l) = q_prim_vf(i)%sf(j, k, l) + end do - type(scalar_field), & - allocatable, dimension(:), & - intent(IN) :: gm_alphaK_vf + ! Zeroing out the dynamic pressure since it is computed + ! iteratively by cycling through the velocity equations + dyn_pres = 0d0 + + ! Computing momenta and dynamic pressure from velocity + do i = mom_idx%beg, mom_idx%end + q_cons_vf(i)%sf(j, k, l) = rho*q_prim_vf(i)%sf(j, k, l) + dyn_pres = dyn_pres + q_cons_vf(i)%sf(j, k, l)* & + q_prim_vf(i)%sf(j, k, l)/2d0 + end do - type(int_bounds_info), intent(IN) :: ix, iy, iz + ! Computing the energy from the pressure + if ((model_eqns /= 4) .and. (bubbles .neqv. .true.)) then + ! E = Gamma*P + \rho u u /2 + \pi_inf + q_cons_vf(E_idx)%sf(j, k, l) = & + gamma*q_prim_vf(E_idx)%sf(j, k, l) + dyn_pres + pi_inf + else if ((model_eqns /= 4) .and. (bubbles)) then + ! \tilde{E} = dyn_pres + (1-\alf)(\Gamma p_l + \Pi_inf) + q_cons_vf(E_idx)%sf(j, k, l) = dyn_pres + & + (1.d0 - q_prim_vf(alf_idx)%sf(j, k, l))* & + (gamma*q_prim_vf(E_idx)%sf(j, k, l) + pi_inf) + else + !Tait EOS, no conserved energy variable + q_cons_vf(E_idx)%sf(j, k, l) = 0. + end if - integer :: j, k, l !< Generic loop iterators + ! Computing the internal energies from the pressure and continuities + if (model_eqns == 3) then + do i = internalEnergies_idx%beg, internalEnergies_idx%end + q_cons_vf(i)%sf(j, k, l) = q_cons_vf(i - adv_idx%end)%sf(j, k, l)* & + fluid_pp(i - adv_idx%end)%gamma* & + q_prim_vf(E_idx)%sf(j, k, l) + & + fluid_pp(i - adv_idx%end)%pi_inf + end do + end if - ! Calculating the momentum and energy from the velocity and pressure - do l = iz%beg, iz%end - do k = iy%beg, iy%end - do j = ix%beg, ix%end + ! Transferring the advection equation(s) variable(s) + do i = adv_idx%beg, adv_idx%end + q_cons_vf(i)%sf(j, k, l) = q_prim_vf(i)%sf(j, k, l) + end do - if (proc_rank == 0) then - print '(A)', 'Conversion from primitive to '// & - 'conservative variables not '// & - 'implemented. Exiting ...' - call s_mpi_abort() + if (bubbles) then + ! From prim: Compute nbub = (3/4pi) * \alpha / \bar{R^3} + do i = 1, nb + Rtmp(i) = q_prim_vf(bub_idx%rs(i))%sf(j, k, l) + end do + !call s_comp_n_from_prim_cpu(q_prim_vf(alf_idx)%sf(j, k, l), Rtmp, nbub) + vftmp = q_prim_vf(alf_idx)%sf(j, k, l) + R3 = 0d0 + do q = 1, nb + R3 = R3 + weight(q)*(Rtmp(q)**3d0) + end do + nbub = (3.d0/(4.d0*pi))*vftmp/R3 + if (j == 0 .and. k == 0 .and. l == 0) print *, 'In convert, nbub:', nbub + do i = bub_idx%beg, bub_idx%end + q_cons_vf(i)%sf(j, k, l) = q_prim_vf(i)%sf(j, k, l)*nbub + ! IF( j==0 .and. k==0 .and. l==0) THEN + ! PRINT*, 'nmom', i, q_cons_vf(i)%sf(j,k,l) + ! END IF + end do end if + if (hypoelasticity) then + do i = stress_idx%beg, stress_idx%end + q_cons_vf(i)%sf(j, k, l) = rho*q_prim_vf(i)%sf(j, k, l) + ! adding elastic contribution + if (G > 1000) then + q_cons_vf(E_idx)%sf(j, k, l) = q_cons_vf(E_idx)%sf(j, k, l) + & + (q_prim_vf(i)%sf(j, k, l)**2d0)/(4d0*G) + ! extra terms in 2 and 3D + if ((i == stress_idx%beg + 1) .or. & + (i == stress_idx%beg + 3) .or. & + (i == stress_idx%beg + 4)) then + q_cons_vf(E_idx)%sf(j, k, l) = q_cons_vf(E_idx)%sf(j, k, l) + & + (q_prim_vf(i)%sf(j, k, l)**2d0)/(4d0*G) + end if + end if + end do + end if end do end do end do +#else + if (proc_rank == 0) then + print '(A)', 'Conversion from primitive to '// & + 'conservative variables not '// & + 'implemented. Exiting ...' + call s_mpi_abort() + end if +#endif + end subroutine s_convert_primitive_to_conservative_variables ! --------- !> The following subroutine handles the conversion between @@ -791,7 +1079,7 @@ subroutine s_convert_primitive_to_flux_variables(qK_prim_vf, & ! ------ ! Computing the flux variables from the primitive variables, without ! accounting for the contribution of either viscosity or capillarity - +#ifdef MFC_SIMULATION !$acc parallel loop collapse(3) gang vector default(present) private(alpha_rho_K, vel_K, alpha_K, Re_K) do l = is3b, is3e do k = is2b, is2e @@ -875,16 +1163,28 @@ subroutine s_convert_primitive_to_flux_variables(qK_prim_vf, & ! ------ end do end do end do +#endif end subroutine s_convert_primitive_to_flux_variables ! ----------------- subroutine s_finalize_variables_conversion_module() ! ------------------ - deallocate ( Gs) + ! Deallocating the density, the specific heat ratio function and the + ! liquid stiffness function +#ifdef MFC_POST_PROCESS + deallocate (rho_sf) + deallocate (gamma_sf) + deallocate (pi_inf_sf) +#endif + + deallocate (gammas, pi_infs, Gs) deallocate (bubrs) + ! Nullifying the procedure pointer to the subroutine transfering/ + ! computing the mixture/species variables to the mixture variables s_convert_to_mixture_variables => null() end subroutine s_finalize_variables_conversion_module ! ---------------- -end module m_variables_conversion + +end module m_variables_conversion \ No newline at end of file diff --git a/src/post_process/m_global_parameters.f90 b/src/post_process/m_global_parameters.f90 index 4051535ad..aea7c5aad 100644 --- a/src/post_process/m_global_parameters.f90 +++ b/src/post_process/m_global_parameters.f90 @@ -235,6 +235,16 @@ module m_global_parameters real(kind(0d0)) :: poly_sigma !> @} + !> @name Index variables used for m_variables_conversion + !> @{ + integer :: momxb, momxe + integer :: advxb, advxe + integer :: contxb, contxe + integer :: intxb, intxe + integer :: bubxb, bubxe + integer :: strxb, strxe + !> @} + ! Mathematical and Physical Constants ====================================== real(kind(0d0)), parameter :: pi = 3.141592653589793d0 ! ========================================================================== @@ -516,6 +526,19 @@ subroutine s_initialize_global_parameters_module() ! ---------------------- end if end if end if + + momxb = mom_idx%beg + momxe = mom_idx%end + advxb = adv_idx%beg + advxe = adv_idx%end + contxb = cont_idx%beg + contxe = cont_idx%end + bubxb = bub_idx%beg + bubxe = bub_idx%end + strxb = stress_idx%beg + strxe = stress_idx%end + intxb = internalEnergies_idx%beg + intxe = internalEnergies_idx%end ! ================================================================== #ifdef MFC_MPI diff --git a/src/post_process/m_variables_conversion.f90 b/src/post_process/m_variables_conversion.f90 deleted file mode 100644 index 918661499..000000000 --- a/src/post_process/m_variables_conversion.f90 +++ /dev/null @@ -1,453 +0,0 @@ -!> -!! @file m_variables_conversion.f90 -!! @brief Contains module m_variables_conversion - -!> @brief This module consists of subroutines used in the conversion of the -!! conservative variables into the primitive variables. In addition, -!! the module also contains subroutines used to compute the mixture -!! variables. For a specific time-step undergoing the post-process, -!! the mixture variables are stored everywhere on the grid so that -!! they may possibly be outputted or used to compute other derived -!! flow quantities. -module m_variables_conversion - - ! Dependencies ============================================================= - use m_derived_types !< Definitions of the derived types - - use m_global_parameters !< Global parameters for the code - ! ========================================================================== - - implicit none - - private; public :: s_initialize_variables_conversion_module, & - s_convert_to_mixture_variables, & - s_convert_mixture_to_mixture_variables, & - s_convert_species_to_mixture_variables_bubbles, & - s_convert_species_to_mixture_variables, & - s_convert_conservative_to_primitive_variables, & - s_finalize_variables_conversion_module - - !> Abstract interface to two subroutines designed for the transfer/conversion - !! of the mixture/species variables to the mixture variables - abstract interface - - !> Structure of the s_convert_mixture_to_mixture_variables - !! and s_convert_species_to_mixture_variables subroutines - !! @param q_cons_vf Conservative variables - !! @param i cell index to transfer mixture variables - !! @param j cell index to transfer mixture variables - !! @param k cell index to transfer mixture variables - subroutine s_convert_xxxxx_to_mixture_variables(q_cons_vf, i, j, k, G_K, G) - - ! Importing the derived type scalar_field from m_derived_types.f90 - ! and global variable sys_size, from m_global_variables.f90, as - ! the abstract interface does not inherently have access to them - import :: scalar_field, sys_size, num_fluids - - type(scalar_field), & - dimension(sys_size), & - intent(IN) :: q_cons_vf - - integer, intent(IN) :: i, j, k - - real(kind(0d0)), optional, intent(out) :: G_K - real(kind(0d0)), optional, dimension(num_fluids), intent(in) :: G - - end subroutine s_convert_xxxxx_to_mixture_variables - - end interface - - ! NOTE: These abstract interfaces allow for the declaration of a pointer to - ! a procedure such that the choice of the model equations does not have to - ! be queried every time the mixture quantities are needed. - - real(kind(0d0)), allocatable, dimension(:, :, :), public :: rho_sf !< Scalar density function - real(kind(0d0)), allocatable, dimension(:, :, :), public :: gamma_sf !< Scalar sp. heat ratio function - real(kind(0d0)), allocatable, dimension(:, :, :), public :: pi_inf_sf !< Scalar liquid stiffness function - - procedure(s_convert_xxxxx_to_mixture_variables), & - pointer :: s_convert_to_mixture_variables => null() !< - !! Pointer referencing the subroutine s_convert_mixture_to_mixture_variables - !! or s_convert_species_to_mixture_variables, based on model equations choice - - integer, private :: flg !< - !! Flagging (flg) variable used to annotate the dimensionality of the dataset - !! that is undergoing the post-process. A flag value of 1 indicates that the - !! dataset is 3D, while a flag value of 0 indicates that it is not. This flg - !! variable is necessary to avoid cycling through the third dimension of the - !! flow variable(s) when the simulation is not 3D and the size of the buffer - !! is non-zero. Note that a similar procedure does not have to be applied to - !! the second dimension since in 1D, the buffer size is always zero. - - integer, private :: flg_2d - -contains - - !> This subroutine is constructed for the gamma/pi_inf model - !! and provided a set of conservative variables, transfers - !! the density, specific heat ratio function and the liquid - !! stiffness function from q_cons_vf to rho_sf, gamma_sf and - !! pi_inf_sf. - !! @param q_cons_vf Conservative variables - !! @param i cell index to transfer mixture variables - !! @param j cell index to transfer mixture variables - !! @param k cell index to transfer mixture variables - subroutine s_convert_mixture_to_mixture_variables(q_cons_vf, i, j, k, G_K, G) ! -- - - type(scalar_field), & - dimension(sys_size), & - intent(IN) :: q_cons_vf - - integer, intent(IN) :: i, j, k - - real(kind(0d0)), optional, intent(out) :: G_K - real(kind(0d0)), optional, dimension(num_fluids), intent(in) :: G - - ! Transfering the density, the specific heat ratio function and the - ! liquid stiffness function, respectively - rho_sf(i, j, k) = q_cons_vf(1)%sf(i, j, k) - gamma_sf(i, j, k) = q_cons_vf(gamma_idx)%sf(i, j, k) - pi_inf_sf(i, j, k) = q_cons_vf(pi_inf_idx)%sf(i, j, k) - - end subroutine s_convert_mixture_to_mixture_variables ! ---------------- - - !> This subroutine is designed for the volume fraction model - !! with sub-grid ensemble averaged bubbles (Bryngelson 2019) - !! and provided a set of conservative variables, calculates - !! the density, the specific heat ratio function and liquid - !! stiffness function from q_cons_vf and stores the results - !! into rho_sf, gamma_sf and pi_inf_sf. - !! @param qK_vf Conservative variables - !! @param j cell index to transfer mixture variables - !! @param k cell index to transfer mixture variables - !! @param l cell index to transfer mixture variables - subroutine s_convert_species_to_mixture_variables_bubbles(q_cons_vf, j, k, l, G_K, G) - - type(scalar_field), dimension(sys_size), intent(IN) :: q_cons_vf - integer, intent(IN) :: j, k, l - - real(kind(0d0)), optional, intent(out) :: G_K - real(kind(0d0)), optional, dimension(num_fluids), intent(in) :: G - - integer :: i !< Generic loop iterator - - if (model_eqns == 4) then - rho_sf(j, k, l) = q_cons_vf(1)%sf(j, k, l) - gamma_sf(j, k, l) = fluid_pp(1)%gamma - pi_inf_sf(j, k, l) = fluid_pp(1)%pi_inf - else if (model_eqns == 2 .and. bubbles) then - rho_sf(j, k, l) = 0d0; gamma_sf(j, k, l) = 0d0; pi_inf_sf(j, k, l) = 0d0 - - if (mpp_lim .and. num_fluids > 1) then - do i = 1, num_fluids - rho_sf(j, k, l) = rho_sf(j, k, l) + q_cons_vf(i + E_idx)%sf(j, k, l)*q_cons_vf(i)%sf(j, k, l) - gamma_sf(j, k, l) = gamma_sf(j, k, l) + q_cons_vf(i + E_idx)%sf(j, k, l)*fluid_pp(i)%gamma - pi_inf_sf(j, k, l) = pi_inf_sf(j, k, l) + q_cons_vf(i + E_idx)%sf(j, k, l)*fluid_pp(i)%pi_inf - end do - else if (num_fluids > 1) then - do i = 1, num_fluids - 1 !leave out bubble part of mixture - rho_sf(j, k, l) = rho_sf(j, k, l) + q_cons_vf(i + E_idx)%sf(j, k, l)*q_cons_vf(i)%sf(j, k, l) - gamma_sf(j, k, l) = gamma_sf(j, k, l) + q_cons_vf(i + E_idx)%sf(j, k, l)*fluid_pp(i)%gamma - pi_inf_sf(j, k, l) = pi_inf_sf(j, k, l) + q_cons_vf(i + E_idx)%sf(j, k, l)*fluid_pp(i)%pi_inf - end do - else - rho_sf(j, k, l) = q_cons_vf(1)%sf(j, k, l) - gamma_sf(j, k, l) = fluid_pp(1)%gamma - pi_inf_sf(j, k, l) = fluid_pp(1)%pi_inf - end if - end if - - end subroutine s_convert_species_to_mixture_variables_bubbles ! ---------------- - - !> This subroutine is designed for the volume fraction model - !! and provided a set of conservative variables, calculates - !! the density, the specific heat ratio function and liquid - !! stiffness function from q_cons_vf and stores the results - !! into rho_sf, gamma_sf and pi_inf_sf. - !! @param q_cons_vf Conservative variables - !! @param j cell index to transfer mixture variables - !! @param k cell index to transfer mixture variables - !! @param l cell index to transfer mixture variables - subroutine s_convert_species_to_mixture_variables(q_cons_vf, j, k, l, G_K, G) ! -- - - type(scalar_field), & - dimension(sys_size), & - intent(IN) :: q_cons_vf - - real(kind(0d0)), optional, intent(out) :: G_K - real(kind(0d0)), optional, dimension(num_fluids), intent(in) :: G - - integer, intent(IN) :: j, k, l - - integer :: i !< Generic loop iterator - - ! Computing the density, the specific heat ratio function and the - ! liquid stiffness function, respectively - if (bubbles .neqv. .true.) then - rho_sf(j, k, l) = 0d0 - gamma_sf(j, k, l) = 0d0 - pi_inf_sf(j, k, l) = 0d0 - - do i = 1, num_fluids - rho_sf(j, k, l) = rho_sf(j, k, l) & - + q_cons_vf(i)%sf(j, k, l) - gamma_sf(j, k, l) = gamma_sf(j, k, l) & - + q_cons_vf(i + E_idx)%sf(j, k, l) & - *fluid_pp(i)%gamma - pi_inf_sf(j, k, l) = pi_inf_sf(j, k, l) & - + q_cons_vf(i + E_idx)%sf(j, k, l) & - *fluid_pp(i)%pi_inf - end do - else - rho_sf(j, k, l) = q_cons_vf(1)%sf(j, k, l) - gamma_sf(j, k, l) = fluid_pp(1)%gamma - pi_inf_sf(j, k, l) = fluid_pp(1)%pi_inf - end if - - if (present(G_K)) then - G_K = 0d0 - do i = 1, num_fluids - G_K = G_K + q_cons_vf(i + E_idx)%sf(j, k, l)*G(i) - end do - G_K = max(0d0, G_K) - end if - - end subroutine s_convert_species_to_mixture_variables ! ---------------- - - !> Computation of parameters, allocation procedures, and/or - !! any other tasks needed to properly setup the module - subroutine s_initialize_variables_conversion_module() ! ------------------- - - ! Allocating the density, the specific heat ratio function and the - ! liquid stiffness function, respectively - - ! Simulation is at least 2D - if (n > 0) then - - ! Simulation is 3D - if (p > 0) then - - allocate (rho_sf(-buff_size:m + buff_size, & - -buff_size:n + buff_size, & - -buff_size:p + buff_size)) - allocate (gamma_sf(-buff_size:m + buff_size, & - -buff_size:n + buff_size, & - -buff_size:p + buff_size)) - allocate (pi_inf_sf(-buff_size:m + buff_size, & - -buff_size:n + buff_size, & - -buff_size:p + buff_size)) - - ! Simulation is 2D - else - - allocate (rho_sf(-buff_size:m + buff_size, & - -buff_size:n + buff_size, & - 0:0)) - allocate (gamma_sf(-buff_size:m + buff_size, & - -buff_size:n + buff_size, & - 0:0)) - allocate (pi_inf_sf(-buff_size:m + buff_size, & - -buff_size:n + buff_size, & - 0:0)) - - end if - - ! Simulation is 1D - else - - allocate (rho_sf(-buff_size:m + buff_size, & - 0:0, & - 0:0)) - allocate (gamma_sf(-buff_size:m + buff_size, & - 0:0, & - 0:0)) - allocate (pi_inf_sf(-buff_size:m + buff_size, & - 0:0, & - 0:0)) - - end if - - ! Depending on the model selection for the equations of motion, the - ! appropriate procedure for the conversion to the mixture variables - ! is targeted by the procedure pointer - - if (model_eqns == 1) then ! Gamma/pi_inf model - s_convert_to_mixture_variables => & - s_convert_mixture_to_mixture_variables - elseif (bubbles) then ! Volume fraction model - s_convert_to_mixture_variables => & - s_convert_species_to_mixture_variables_bubbles - else ! Volume fraction model - s_convert_to_mixture_variables => & - s_convert_species_to_mixture_variables - end if - - ! Annotating the dimensionality of the dataset undergoing the post- - ! process. A flag value of 1 indicates that the dataset is 3D, while - ! a flag value of 0 indicates that it is not. - if (p > 0) then - flg = 1 - else - flg = 0 - end if - - end subroutine s_initialize_variables_conversion_module ! ----------------- - - !> Converts the conservative variables to the primitive ones - !! @param q_cons_vf Conservative variabels - !! @param q_prim_vf Primitive variables - subroutine s_convert_conservative_to_primitive_variables(q_cons_vf, & - q_prim_vf) - - type(scalar_field), & - dimension(sys_size), & - intent(IN) :: q_cons_vf - - type(scalar_field), & - dimension(sys_size), & - intent(INOUT) :: q_prim_vf - - ! Dynamic pressure, as defined in the incompressible flow sense - real(kind(0d0)) :: dyn_pres - - ! Bubble parameters - real(kind(0d0)) :: nbub - real(kind(0d0)), dimension(:), allocatable :: nRtmp - - real(kind(0d0)) :: G_K - - integer :: i, j, k, l !< Generic loop iterators - - if (bubbles) then - allocate (nRtmp(nb)) - end if - - ! Converting the conservative variables to the primitive variables - do l = -buff_size*flg, (p + buff_size)*flg - do k = -buff_size, n + buff_size - do j = -buff_size, m + buff_size - - ! Obtaining the density, specific heat ratio function - ! and the liquid stiffness function, respectively - if (hypoelasticity) then - call s_convert_to_mixture_variables(q_cons_vf, j, k, l, & - G_K, fluid_pp(:)%G) - else - call s_convert_to_mixture_variables(q_cons_vf, j, k, l) - end if - - ! Transferring the continuity equation(s) variable(s) - do i = 1, cont_idx%end - q_prim_vf(i)%sf(j, k, l) = q_cons_vf(i)%sf(j, k, l) - end do - - ! Zeroing out the dynamic pressure since it is computed - ! iteratively by cycling through the momentum equations - dyn_pres = 0d0 - - ! Computing velocity and dynamic pressure from momenta - do i = mom_idx%beg, mom_idx%end - if (model_eqns == 4) then - !u = \rho u/ \rho - q_prim_vf(i)%sf(j, k, l) = q_cons_vf(i)%sf(j, k, l)/ & - q_cons_vf(1)%sf(j, k, l) - else - q_prim_vf(i)%sf(j, k, l) = q_cons_vf(i)%sf(j, k, l)/ & - rho_sf(j, k, l) - dyn_pres = dyn_pres + q_cons_vf(i)%sf(j, k, l)* & - q_prim_vf(i)%sf(j, k, l)/2d0 - end if - end do - - if (model_eqns == 4) then - ! Computing the pressure from the energy - ! Tait EOS - ! p_l = (pl0 + pi_infty)(rho/(rho_l0(1-alf)))^gamma - ! - pi_infty - q_prim_vf(E_idx)%sf(j, k, l) = & - (pref + fluid_pp(1)%pi_inf)* & - (( & - q_cons_vf(1)%sf(j, k, l)/ & - (rhoref*(1.d0 - q_cons_vf(alf_idx)%sf(j, k, l))) & - )**(1.d0/fluid_pp(1)%gamma + 1.d0)) - fluid_pp(1)%pi_inf - else if ((model_eqns /= 4) .and. bubbles .neqv. .true.) then - ! Computing the pressure from the energy - q_prim_vf(E_idx)%sf(j, k, l) = & - (q_cons_vf(E_idx)%sf(j, k, l) & - - dyn_pres - pi_inf_sf(j, k, l))/gamma_sf(j, k, l) - else - ! p = ( E/(1-alf) - 0.5 rho u u/(1-alf) - pi_inf_k )/gamma_k - q_prim_vf(E_idx)%sf(j, k, l) = & - ((q_cons_vf(E_idx)%sf(j, k, l) & - - dyn_pres)/(1.d0 - q_cons_vf(alf_idx)%sf(j, k, l)) & - - pi_inf_sf(j, k, l) & - )/gamma_sf(j, k, l) - end if - - ! Set partial pressures to mixture pressure - if (model_eqns == 3) then - do i = internalEnergies_idx%beg, internalEnergies_idx%end - q_prim_vf(i)%sf(j, k, l) = q_prim_vf(E_idx)%sf(j, k, l) - end do - end if - - ! Transferring the advection equation(s) variable(s) - do i = adv_idx%beg, adv_idx%end - q_prim_vf(i)%sf(j, k, l) = q_cons_vf(i)%sf(j, k, l) - end do - - ! \phi = (n\phi)/n (n = nbub) - if (bubbles) then - ! n = sqrt( 4pi/(3 alpha) * (nR)**3 ) - do i = 1, nb - nRtmp(i) = q_cons_vf(bub_idx%rs(i))%sf(j, k, l) - end do - call s_comp_n_from_cons(q_cons_vf(alf_idx)%sf(j, k, l), nRtmp, nbub) - do i = bub_idx%beg, sys_size - q_prim_vf(i)%sf(j, k, l) = q_cons_vf(i)%sf(j, k, l)/nbub - end do - end if - - if (hypoelasticity) then - do i = stress_idx%beg, stress_idx%end - q_prim_vf(i)%sf(j, k, l) = q_cons_vf(i)%sf(j, k, l)/ & - rho_sf(j, k, l) - if (G_K > 1000) then - q_prim_vf(E_idx)%sf(j, k, l) = q_prim_vf(E_idx)%sf(j, k, l) - & - ((q_prim_vf(i)%sf(j, k, l)**2d0)/(4d0*G_K))/gamma_sf(j, k, l) - ! 2D and 3D terms - if ((i == stress_idx%beg + 1) .or. & - (i == stress_idx%beg + 3) .or. & - (i == stress_idx%beg + 4)) then - q_prim_vf(E_idx)%sf(j, k, l) = q_prim_vf(E_idx)%sf(j, k, l) - & - ((q_prim_vf(i)%sf(j, k, l)**2d0)/(4d0*G_K))/gamma_sf(j, k, l) - end if - end if - end do - end if - end do - end do - end do - - if (bubbles) then - deallocate (nRtmp) - end if - - end subroutine s_convert_conservative_to_primitive_variables ! --------- - - !> Deallocation procedures for the module - subroutine s_finalize_variables_conversion_module() ! ---------------- - - ! Deallocating the density, the specific heat ratio function and the - ! liquid stiffness function - deallocate (rho_sf) - deallocate (gamma_sf) - deallocate (pi_inf_sf) - - ! Nullifying the procedure pointer to the subroutine transfering/ - ! computing the mixture/species variables to the mixture variables - s_convert_to_mixture_variables => null() - - end subroutine s_finalize_variables_conversion_module ! -------------- - -end module m_variables_conversion diff --git a/src/pre_process/m_data_output.f90 b/src/pre_process/m_data_output.f90 index 42aed533a..706f81d9f 100644 --- a/src/pre_process/m_data_output.f90 +++ b/src/pre_process/m_data_output.f90 @@ -86,6 +86,7 @@ subroutine s_write_serial_data_files(q_cons_vf) ! ----------- real(kind(0d0)) :: nbub !< Temporary bubble number density real(kind(0d0)) :: gamma, lit_gamma, pi_inf !< Temporary EOS params real(kind(0d0)) :: rho !< Temporary density + real(kind(0d0)) :: pres !< Temporary pressure t_step = 0 @@ -168,33 +169,9 @@ subroutine s_write_serial_data_files(q_cons_vf) ! ----------- else if (i == stress_idx%beg) then !tau_e write (2, FMT) x_cb(j), q_cons_vf(stress_idx%beg)%sf(j, 0, 0)/rho else if (i == E_idx) then !p - if (model_eqns == 4) then - !Tait pressure from density - write (2, FMT) x_cb(j), & - (pref + pi_inf)*( & - (q_cons_vf(1)%sf(j, 0, 0)/ & - (rhoref*(1.d0 - q_cons_vf(4)%sf(j, 0, 0))) & - )**lit_gamma) & - - pi_inf - !TODO: add elastic contribution here for cases with initial tau_e nonzero - else if (model_eqns == 2 .and. (bubbles .neqv. .true.)) then - !Stiffened gas pressure from energy - write (2, FMT) x_cb(j), & - ( & - q_cons_vf(E_idx)%sf(j, 0, 0) - & - 0.5d0*(q_cons_vf(mom_idx%beg)%sf(j, 0, 0)**2.d0)/rho - & - pi_inf & - )/gamma - else - !Stiffened gas pressure from energy with bubbles - write (2, FMT) x_cb(j), & - ( & - (q_cons_vf(E_idx)%sf(j, 0, 0) - & - 0.5d0*(q_cons_vf(mom_idx%beg)%sf(j, 0, 0)**2.d0)/rho)/ & - (1.d0 - q_cons_vf(alf_idx)%sf(j, 0, 0)) - & - pi_inf & - )/gamma - end if + call s_compute_pressure(q_cons_vf(E_idx)%sf(j, 0, 0), q_cons_vf(alf_idx)%sf(j, 0, 0), & + 0.5d0*(q_cons_vf(mom_idx%beg)%sf(j, 0, 0)**2.d0)/rho, pi_inf, gamma, pres) + write (2, FMT) x_cb(j), pres else if ((i >= bub_idx%beg) .and. (i <= bub_idx%end) .and. bubbles) then do k = 1, nb nRtmp(k) = q_cons_vf(bub_idx%rs(k))%sf(j, 0, 0) diff --git a/src/pre_process/m_global_parameters.fpp b/src/pre_process/m_global_parameters.fpp index c80f1bb7f..2b7fc03a2 100644 --- a/src/pre_process/m_global_parameters.fpp +++ b/src/pre_process/m_global_parameters.fpp @@ -174,6 +174,16 @@ module m_global_parameters integer :: R0_type !1 = simpson !> @} + !> @name Index variables used for m_variables_conversion + !> @{ + integer :: momxb, momxe + integer :: advxb, advxe + integer :: contxb, contxe + integer :: intxb, intxe + integer :: bubxb, bubxe + integer :: strxb, strxe + !> @} + integer, allocatable, dimension(:, :, :) :: logic_grid ! Mathematical and Physical Constants ====================================== @@ -560,6 +570,20 @@ contains end if end if end if + + momxb = mom_idx%beg + momxe = mom_idx%end + advxb = adv_idx%beg + advxe = adv_idx%end + contxb = cont_idx%beg + contxe = cont_idx%end + bubxb = bub_idx%beg + bubxe = bub_idx%end + strxb = stress_idx%beg + strxe = stress_idx%end + intxb = internalEnergies_idx%beg + intxe = internalEnergies_idx%end + ! ================================================================== #ifdef MFC_MPI @@ -843,7 +867,7 @@ contains !! @param mom is the computed moment subroutine s_quad(func, mom) - real(kind(0.d0)), dimension(nb), intent(IN) :: func + real(kind(0.d0)), dimension(nb), intent(IN) :: func real(kind(0.d0)), intent(OUT) :: mom mom = dot_product(weight, func) diff --git a/src/pre_process/m_initial_condition.fpp b/src/pre_process/m_initial_condition.fpp index 0913b92e7..809107930 100644 --- a/src/pre_process/m_initial_condition.fpp +++ b/src/pre_process/m_initial_condition.fpp @@ -100,7 +100,7 @@ module m_initial_condition !! Depending on the multicomponent flow model, this variable is a pointer to !! either the subroutine s_assign_patch_mixture_primitive_variables, or the !! subroutine s_assign_patch_species_primitive_variables - + contains !> This subroutine assigns the mixture primitive variables @@ -757,6 +757,7 @@ contains s_assign_patch_species_primitive_variables end if + end subroutine s_initialize_initial_condition_module ! ----------------- !> This subroutine peruses the patches and depending on the diff --git a/src/pre_process/m_variables_conversion.f90 b/src/pre_process/m_variables_conversion.f90 deleted file mode 100644 index 4e117f63a..000000000 --- a/src/pre_process/m_variables_conversion.f90 +++ /dev/null @@ -1,523 +0,0 @@ -!> -!! @file m_variables_conversion.f90 -!! @brief Contains module m_variables_conversion - -!> @brief This module consists of subroutines used in the conversion of the -!! conservative variables into the primitive ones and vice versa. In -!! addition, the module also contains the subroutines used to obtain -!! the mixture variables. -module m_variables_conversion - - ! Dependencies ============================================================= - use m_derived_types !< Definitions of the derived types - - use m_global_parameters !< Global parameters for the code - ! ========================================================================== - - implicit none - - !> Abstract interface to two subroutines designed for the transfer/conversion - !! of the mixture/species variables to the mixture variables - abstract interface - - !> Structure of the s_convert_mixture_to_mixture_variables - !! and s_convert_species_to_mixture_variables subroutines - !! @param q_vf Conservative or primitive variables - !! @param i First-coordinate cell index - !! @param j First-coordinate cell index - !! @param k First-coordinate cell index - !! @param rho Density - !! @param gamma Specific heat ratio function - !! @param pi_inf Liquid stiffness function - subroutine s_convert_xxxxx_to_mixture_variables(q_vf, i, j, k, & - rho, gamma, pi_inf, G) - - ! Importing the derived type scalar_field from m_derived_types.f90 - ! and global variable sys_size, from m_global_variables.f90, as - ! the abstract interface does not inherently have access to them - import :: scalar_field, sys_size - - type(scalar_field), dimension(sys_size), intent(IN) :: q_vf - - integer, intent(IN) :: i, j, k - - real(kind(0d0)), intent(OUT) :: rho - real(kind(0d0)), intent(OUT) :: gamma - real(kind(0d0)), intent(OUT) :: pi_inf - - real(kind(0d0)), optional, intent(OUT) :: G - - end subroutine s_convert_xxxxx_to_mixture_variables - - end interface - - ! NOTE: These abstract interfaces allow for the declaration of a pointer to - ! a procedure such that the choice of the model equations does not have to - ! be queried every time the mixture quantities are needed. - - procedure(s_convert_xxxxx_to_mixture_variables), & - pointer :: s_convert_to_mixture_variables => null() !< - !! Pointer referencing the subroutine s_convert_mixture_to_mixture_variables - !! or s_convert_species_to_mixture_variables, based on model equations choice - -contains - - !> This subroutine is designed for the gamma/pi_inf model - !! and provided a set of either conservative or primitive - !! variables, transfers the density, specific heat ratio - !! function and the liquid stiffness function from q_vf to - !! rho, gamma and pi_inf. - !! @param q_vf conservative or primitive variables - !! @param i cell index to transfer mixture variables - !! @param j cell index to transfer mixture variables - !! @param k cell index to transfer mixture variables - !! @param rho density - !! @param gamma specific heat ratio function - !! @param pi_inf liquid stiffness - subroutine s_convert_mixture_to_mixture_variables(q_vf, i, j, k, & - rho, gamma, pi_inf, G) - - type(scalar_field), dimension(sys_size), intent(IN) :: q_vf - integer, intent(IN) :: i, j, k - - real(kind(0d0)), intent(OUT) :: rho - real(kind(0d0)), intent(OUT) :: gamma - real(kind(0d0)), intent(OUT) :: pi_inf - - real(kind(0d0)), optional, intent(OUT) :: G - - ! Transfering the density, the specific heat ratio function and the - ! liquid stiffness function, respectively - rho = q_vf(1)%sf(i, j, k) - gamma = q_vf(gamma_idx)%sf(i, j, k) - pi_inf = q_vf(pi_inf_idx)%sf(i, j, k) - - end subroutine s_convert_mixture_to_mixture_variables ! ---------------- - - !> This procedure is used alongside with the gamma/pi_inf - !! model to transfer the density, the specific heat ratio - !! function and liquid stiffness function from the vector - !! of conservative or primitive variables to their scalar - !! counterparts. Specifially designed for when subgrid bubbles - !! must be included. - !! @param qK_vf primitive variables - !! @param rho_K density - !! @param gamma_K specific heat ratio - !! @param pi_inf_K liquid stiffness - !! @param j Cell index - !! @param k Cell index - !! @param l Cell index - subroutine s_convert_species_to_mixture_variables_bubbles(qK_vf, & - j, k, l, & - rho_K, gamma_K, pi_inf_K, & - G) - - type(scalar_field), dimension(sys_size), intent(IN) :: qK_vf - - real(kind(0d0)), intent(OUT) :: rho_K, gamma_K, pi_inf_K - - real(kind(0d0)), dimension(num_fluids) :: alpha_rho_K, alpha_K !< - !! Partial densities and volume fractions - - real(kind(0d0)), optional, intent(OUT) :: G - - integer, intent(IN) :: j, k, l - integer :: i - - ! Constraining the partial densities and the volume fractions within - ! their physical bounds to make sure that any mixture variables that - ! are derived from them result within the limits that are set by the - ! fluids physical parameters that make up the mixture - ! alpha_rho_K(1) = qK_vf(i)%sf(i,j,k) - ! alpha_K(1) = qK_vf(E_idx+i)%sf(i,j,k) - - ! Performing the transfer of the density, the specific heat ratio - ! function as well as the liquid stiffness function, respectively - - if (model_eqns == 4) then - rho_K = qK_vf(1)%sf(j, k, l) - gamma_K = fluid_pp(1)%gamma !qK_vf(gamma_idx)%sf(i,j,k) - pi_inf_K = fluid_pp(1)%pi_inf !qK_vf(pi_inf_idx)%sf(i,j,k) - else if ((model_eqns == 2) .and. bubbles) then - rho_k = 0d0; gamma_k = 0d0; pi_inf_k = 0d0 - - if (mpp_lim .and. (num_fluids > 2)) then - do i = 1, num_fluids - rho_k = rho_k + qK_vf(i)%sf(j, k, l) - gamma_k = gamma_k + qK_vf(i + E_idx)%sf(j, k, l)*fluid_pp(i)%gamma - pi_inf_k = pi_inf_k + qK_vf(i + E_idx)%sf(j, k, l)*fluid_pp(i)%pi_inf - end do - else if (num_fluids == 2) then - rho_K = qK_vf(1)%sf(j, k, l) - gamma_K = fluid_pp(1)%gamma - pi_inf_K = fluid_pp(1)%pi_inf - else if (num_fluids > 2) then - !TODO: This may need fixing for hypo + bubbles - do i = 1, num_fluids - 1 !leave out bubble part of mixture - rho_k = rho_k + qK_vf(i)%sf(j, k, l) - gamma_k = gamma_k + qK_vf(i + E_idx)%sf(j, k, l)*fluid_pp(i)%gamma - pi_inf_k = pi_inf_k + qK_vf(i + E_idx)%sf(j, k, l)*fluid_pp(i)%pi_inf - end do - !rho_K = qK_vf(1)%sf(j,k,l) - !gamma_K = fluid_pp(1)%gamma - !pi_inf_K = fluid_pp(1)%pi_inf - else - rho_K = qK_vf(1)%sf(j, k, l) - gamma_K = fluid_pp(1)%gamma - pi_inf_K = fluid_pp(1)%pi_inf - end if - end if - - end subroutine s_convert_species_to_mixture_variables_bubbles ! ---------------- - - !> This subroutine is designed for the volume fraction model - !! and provided a set of either conservative or primitive - !! variables, computes the density, the specific heat ratio - !! function and the liquid stiffness function from q_vf and - !! stores the results into rho, gamma and pi_inf. - !! @param q_vf primitive variables - !! @param rho density - !! @param gamma specific heat ratio - !! @param pi_inf liquid stiffness - !! @param j Cell index - !! @param k Cell index - !! @param l Cell index - subroutine s_convert_species_to_mixture_variables(q_vf, j, k, l, & - rho, gamma, pi_inf, G) - - type(scalar_field), dimension(sys_size), intent(IN) :: q_vf - - integer, intent(IN) :: j, k, l !< - ! Indices of the cell for which to compute the mixture variables - - real(kind(0d0)), intent(OUT) :: rho - real(kind(0d0)), intent(OUT) :: gamma - real(kind(0d0)), intent(OUT) :: pi_inf - - real(kind(0d0)), optional, intent(OUT) :: G - - integer :: i !< Generic loop iterator - - ! Computing the density, the specific heat ratio function and the - ! liquid stiffness function, respectively - - rho = 0d0; gamma = 0d0; pi_inf = 0d0 - - do i = 1, num_fluids - rho = rho + q_vf(i)%sf(j, k, l) - gamma = gamma + q_vf(i + E_idx)%sf(j, k, l)*fluid_pp(i)%gamma - pi_inf = pi_inf + q_vf(i + E_idx)%sf(j, k, l)*fluid_pp(i)%pi_inf - end do - - if (present(G)) then - G = 0d0 - do i = 1, num_fluids - G = G + q_vf(i + E_idx)%sf(j, k, l)*fluid_pp(i)%G - end do - G = max(0d0, G) - end if - - end subroutine s_convert_species_to_mixture_variables ! ---------------- - - !> Computation of parameters, allocation procedures, and/or - !! any other tasks needed to properly setup the module - subroutine s_initialize_variables_conversion_module() ! ------------------- - - ! Depending on the model selection for the equations of motion, the - ! appropriate procedure for the conversion to the mixture variables - ! is targeted by the procedure pointer - - if (model_eqns == 1) then ! Gamma/pi_inf model - s_convert_to_mixture_variables => & - s_convert_mixture_to_mixture_variables - - else if (bubbles) then - s_convert_to_mixture_variables => & - s_convert_species_to_mixture_variables_bubbles - else - ! Volume fraction model - s_convert_to_mixture_variables => & - s_convert_species_to_mixture_variables - end if - - end subroutine s_initialize_variables_conversion_module ! ----------------- - - !> Converts the conservative variables to the primitive ones - !! @param q_cons_vf Conservative variables - !! @param q_prim_vf Primitive variables - subroutine s_convert_conservative_to_primitive_variables(q_cons_vf, & - q_prim_vf) - - type(scalar_field), & - dimension(sys_size), & - intent(IN) :: q_cons_vf - - type(scalar_field), & - dimension(sys_size), & - intent(INOUT) :: q_prim_vf - - ! Density, specific heat ratio function, liquid stiffness function - ! and dynamic pressure, as defined in the incompressible flow sense, - ! respectively - real(kind(0d0)) :: rho - real(kind(0d0)) :: gamma - real(kind(0d0)) :: pi_inf - real(kind(0d0)) :: dyn_pres - real(kind(0d0)) :: nbub, nR3, vftmp - real(kind(0d0)), dimension(nb) :: nRtmp - - real(kind(0d0)) :: G - - ! Generic loop iterators - integer :: i, j, k, l, q - - ! Converting the conservative variables to the primitive variables - do l = 0, p - do k = 0, n - do j = 0, m - - ! Obtaining the density, specific heat ratio function - ! and the liquid stiffness function, respectively - if (model_eqns /= 4) then - if (hypoelasticity) then - call s_convert_to_mixture_variables(q_cons_vf, j, k, l, & - rho, gamma, pi_inf, G) - else - call s_convert_to_mixture_variables(q_cons_vf, j, k, l, & - rho, gamma, pi_inf) - end if - end if - - ! Transferring the continuity equation(s) variable(s) - do i = 1, cont_idx%end - q_prim_vf(i)%sf(j, k, l) = q_cons_vf(i)%sf(j, k, l) - end do - - ! Zeroing out the dynamic pressure since it is computed - ! iteratively by cycling through the momentum equations - dyn_pres = 0d0 - - ! Computing velocity and dynamic pressure from momenta - do i = mom_idx%beg, mom_idx%end - if (model_eqns /= 4) then - q_prim_vf(i)%sf(j, k, l) = q_cons_vf(i)%sf(j, k, l)/rho - dyn_pres = dyn_pres + q_cons_vf(i)%sf(j, k, l)* & - q_prim_vf(i)%sf(j, k, l)/2d0 - else - q_prim_vf(i)%sf(j, k, l) = q_cons_vf(i)%sf(j, k, l) & - /q_cons_vf(1)%sf(j, k, l) - end if - end do - - if ((model_eqns /= 4) .and. (bubbles .neqv. .true.)) then - ! Computing the pressure from the energy - q_prim_vf(E_idx)%sf(j, k, l) = & - (q_cons_vf(E_idx)%sf(j, k, l) - dyn_pres - pi_inf)/gamma - else if ((model_eqns /= 4) .and. bubbles) then - print *, 'getting model_eqns 2 with bubbles. cons to prim' - ! p = ( E/(1-alf) - 0.5 rho u u/(1-alf) - pi_inf_k )/gamma_k - q_prim_vf(E_idx)%sf(j, k, l) = & - ((q_cons_vf(E_idx)%sf(j, k, l) - dyn_pres)/(1.d0 - q_cons_vf(alf_idx)%sf(j, k, l)) & - - pi_inf)/gamma - else - ! Tait EOS - ! p = (pl0 + pi_infty)(rho/(rho_l0(1-alf)))^gamma - pi_infty - q_prim_vf(E_idx)%sf(j, k, l) = & - (pref + fluid_pp(1)%pi_inf)* & - ( & - q_prim_vf(1)%sf(j, k, l)/ & - (rhoref*(1 - q_prim_vf(alf_idx)%sf(j, k, l))) & - )**(1/fluid_pp(1)%gamma + 1) - fluid_pp(1)%pi_inf - end if - - ! Set partial pressures to mixture pressure - if (model_eqns == 3) then - do i = internalEnergies_idx%beg, internalEnergies_idx%end - q_prim_vf(i)%sf(j, k, l) = q_prim_vf(E_idx)%sf(j, k, l) - end do - end if - - ! Transfer the advection equation(s) variable(s) - do i = adv_idx%beg, adv_idx%end - q_prim_vf(i)%sf(j, k, l) = q_cons_vf(i)%sf(j, k, l) - end do - - if (bubbles) then - ! From cons: ntmp = DSQRT( (4.d0*pi/3.d0)*nR3/vftmp ) - do i = 1, nb - nRtmp(i) = q_cons_vf(bub_idx%rs(i))%sf(j, k, l) - end do - vftmp = q_cons_vf(alf_idx)%sf(j, k, l) - nR3 = 0d0 - do q = 1, nb - nR3 = nR3 + weight(q)*(nRtmp(q)**3d0) - end do - nbub = DSQRT((4.d0*pi/3.d0)*nR3/vftmp) - do i = bub_idx%beg, bub_idx%end - q_prim_vf(i)%sf(j, k, l) = q_cons_vf(i)%sf(j, k, l)/nbub - end do - end if - - if (hypoelasticity) then - do i = stress_idx%beg, stress_idx%end - q_prim_vf(i)%sf(j, k, l) = q_cons_vf(i)%sf(j, k, l)/rho - ! subtracting elastic contribution for pressure calculation - if (G > 1000) then !TODO: Change to >0 if stable - q_prim_vf(E_idx)%sf(j, k, l) = q_prim_vf(E_idx)%sf(j, k, l) - & - ((q_prim_vf(i)%sf(j, k, l)**2d0)/(4d0*G))/gamma - ! extra terms in 2 and 3D - if ((i == stress_idx%beg + 1) .or. & - (i == stress_idx%beg + 3) .or. & - (i == stress_idx%beg + 4)) then - q_prim_vf(E_idx)%sf(j, k, l) = q_prim_vf(E_idx)%sf(j, k, l) - & - ((q_prim_vf(i)%sf(j, k, l)**2d0)/(4d0*G))/gamma - end if - end if - end do - end if - - end do - end do - end do - - end subroutine s_convert_conservative_to_primitive_variables ! --------- - - !> Converts the primitive variables to the conservative ones. - !! Used when initializing patches. - !! @param q_cons_vf Conservative variables - !! @param q_prim_vf Primitive variables - subroutine s_convert_primitive_to_conservative_variables(q_prim_vf, & - q_cons_vf) - - type(scalar_field), & - dimension(sys_size), & - intent(IN) :: q_prim_vf - - type(scalar_field), & - dimension(sys_size), & - intent(INOUT) :: q_cons_vf - - ! Density, specific heat ratio function, liquid stiffness function - ! and dynamic pressure, as defined in the incompressible flow sense, - ! respectively - real(kind(0d0)) :: rho - real(kind(0d0)) :: gamma - real(kind(0d0)) :: pi_inf - real(kind(0d0)) :: dyn_pres - real(kind(0d0)) :: nbub, R3, vftmp - real(kind(0d0)), dimension(nb) :: Rtmp - - real(kind(0d0)) :: G - - integer :: i, j, k, l, q !< Generic loop iterators - - ! Converting the primitive variables to the conservative variables - do l = 0, p - do k = 0, n - do j = 0, m - - ! Obtaining the density, specific heat ratio function - ! and the liquid stiffness function, respectively - call s_convert_to_mixture_variables(q_prim_vf, j, k, l, & - rho, gamma, pi_inf) - - ! Transferring the continuity equation(s) variable(s) - do i = 1, cont_idx%end - q_cons_vf(i)%sf(j, k, l) = q_prim_vf(i)%sf(j, k, l) - end do - - ! Zeroing out the dynamic pressure since it is computed - ! iteratively by cycling through the velocity equations - dyn_pres = 0d0 - - ! Computing momenta and dynamic pressure from velocity - do i = mom_idx%beg, mom_idx%end - q_cons_vf(i)%sf(j, k, l) = rho*q_prim_vf(i)%sf(j, k, l) - dyn_pres = dyn_pres + q_cons_vf(i)%sf(j, k, l)* & - q_prim_vf(i)%sf(j, k, l)/2d0 - end do - - ! Computing the energy from the pressure - if ((model_eqns /= 4) .and. (bubbles .neqv. .true.)) then - ! E = Gamma*P + \rho u u /2 + \pi_inf - q_cons_vf(E_idx)%sf(j, k, l) = & - gamma*q_prim_vf(E_idx)%sf(j, k, l) + dyn_pres + pi_inf - else if ((model_eqns /= 4) .and. (bubbles)) then - ! \tilde{E} = dyn_pres + (1-\alf)(\Gamma p_l + \Pi_inf) - q_cons_vf(E_idx)%sf(j, k, l) = dyn_pres + & - (1.d0 - q_prim_vf(alf_idx)%sf(j, k, l))* & - (gamma*q_prim_vf(E_idx)%sf(j, k, l) + pi_inf) - else - !Tait EOS, no conserved energy variable - q_cons_vf(E_idx)%sf(j, k, l) = 0. - end if - - ! Computing the internal energies from the pressure and continuities - if (model_eqns == 3) then - do i = internalEnergies_idx%beg, internalEnergies_idx%end - q_cons_vf(i)%sf(j, k, l) = q_cons_vf(i - adv_idx%end)%sf(j, k, l)* & - fluid_pp(i - adv_idx%end)%gamma* & - q_prim_vf(E_idx)%sf(j, k, l) + & - fluid_pp(i - adv_idx%end)%pi_inf - end do - end if - - ! Transferring the advection equation(s) variable(s) - do i = adv_idx%beg, adv_idx%end - q_cons_vf(i)%sf(j, k, l) = q_prim_vf(i)%sf(j, k, l) - end do - - if (bubbles) then - ! From prim: Compute nbub = (3/4pi) * \alpha / \bar{R^3} - do i = 1, nb - Rtmp(i) = q_prim_vf(bub_idx%rs(i))%sf(j, k, l) - end do - !call s_comp_n_from_prim_cpu(q_prim_vf(alf_idx)%sf(j, k, l), Rtmp, nbub) - vftmp = q_prim_vf(alf_idx)%sf(j, k, l) - R3 = 0d0 - do q = 1, nb - R3 = R3 + weight(q)*(Rtmp(q)**3d0) - end do - nbub = (3.d0/(4.d0*pi))*vftmp/R3 - if (j == 0 .and. k == 0 .and. l == 0) print *, 'In convert, nbub:', nbub - do i = bub_idx%beg, bub_idx%end - q_cons_vf(i)%sf(j, k, l) = q_prim_vf(i)%sf(j, k, l)*nbub - ! IF( j==0 .and. k==0 .and. l==0) THEN - ! PRINT*, 'nmom', i, q_cons_vf(i)%sf(j,k,l) - ! END IF - end do - end if - - if (hypoelasticity) then - do i = stress_idx%beg, stress_idx%end - q_cons_vf(i)%sf(j, k, l) = rho*q_prim_vf(i)%sf(j, k, l) - ! adding elastic contribution - if (G > 1000) then - q_cons_vf(E_idx)%sf(j, k, l) = q_cons_vf(E_idx)%sf(j, k, l) + & - (q_prim_vf(i)%sf(j, k, l)**2d0)/(4d0*G) - ! extra terms in 2 and 3D - if ((i == stress_idx%beg + 1) .or. & - (i == stress_idx%beg + 3) .or. & - (i == stress_idx%beg + 4)) then - q_cons_vf(E_idx)%sf(j, k, l) = q_cons_vf(E_idx)%sf(j, k, l) + & - (q_prim_vf(i)%sf(j, k, l)**2d0)/(4d0*G) - end if - end if - end do - end if - end do - end do - end do - - end subroutine s_convert_primitive_to_conservative_variables ! --------- - - !> Deallocation procedures for the module - subroutine s_finalize_variables_conversion_module() ! ---------------- - - ! Nullifying the procedure pointer to the subroutine transfering/ - ! computing the mixture/species variables to the mixture variables - s_convert_to_mixture_variables => null() - - end subroutine s_finalize_variables_conversion_module ! -------------- - -end module m_variables_conversion diff --git a/src/simulation/m_data_output.f90 b/src/simulation/m_data_output.f90 index d174b9e11..0c5be03c0 100644 --- a/src/simulation/m_data_output.f90 +++ b/src/simulation/m_data_output.f90 @@ -691,7 +691,7 @@ subroutine s_write_serial_data_files(q_cons_vf, t_step) ! --------------------- open (2, FILE=trim(file_path)) do j = 0, m - call s_convert_to_mixture_variables(q_cons_vf, rho, gamma, pi_inf, Re, j, 0, 0, & + call s_convert_to_mixture_variables(q_cons_vf, j, 0, 0, rho, gamma, pi_inf, Re, & G, fluid_pp(:)%G) lit_gamma = 1d0/gamma + 1d0 @@ -1743,10 +1743,9 @@ subroutine s_write_probe_files(t_step, q_cons_vf, accel_mag) ! ----------- l = 0 ! Computing/Sharing necessary state variables - call s_convert_to_mixture_variables(q_cons_vf, rho, & - gamma, pi_inf, & - Re, j - 2, k, l, & - G, fluid_pp(:)%G) + call s_convert_to_mixture_variables(q_cons_vf, j - 2, k, l, & + rho, gamma, pi_inf, & + Re, G, fluid_pp(:)%G) do s = 1, num_dims vel(s) = q_cons_vf(cont_idx%end + s)%sf(j - 2, k, l)/rho end do @@ -1885,10 +1884,9 @@ subroutine s_write_probe_files(t_step, q_cons_vf, accel_mag) ! ----------- l = 0 ! Computing/Sharing necessary state variables - call s_convert_to_mixture_variables(q_cons_vf, rho, & - gamma, pi_inf, & - Re, j - 2, k - 2, l, & - G, fluid_pp(:)%G) + call s_convert_to_mixture_variables(q_cons_vf, j - 2, k - 2, l, & + rho, gamma, pi_inf, & + Re, G, fluid_pp(:)%G) do s = 1, num_dims vel(s) = q_cons_vf(cont_idx%end + s)%sf(j - 2, k - 2, l)/rho end do @@ -2014,10 +2012,9 @@ subroutine s_write_probe_files(t_step, q_cons_vf, accel_mag) ! ----------- if (l == 1) l = 2 ! Pick first point if probe is at edge ! Computing/Sharing necessary state variables - call s_convert_to_mixture_variables(q_cons_vf, rho, & - gamma, pi_inf, & - Re, j - 2, k - 2, l - 2, & - G, fluid_pp(:)%G) + call s_convert_to_mixture_variables(q_cons_vf, j - 2, k - 2, l - 2, & + rho, gamma, pi_inf, & + Re, G, fluid_pp(:)%G) do s = 1, num_dims vel(s) = q_cons_vf(cont_idx%end + s)%sf(j - 2, k - 2, l - 2)/rho end do @@ -2257,9 +2254,8 @@ subroutine s_write_probe_files(t_step, q_cons_vf, accel_mag) ! ----------- if ((integral(i)%xmin <= x_cb(j)) .and. (integral(i)%xmax >= x_cb(j))) then npts = npts + 1 - call s_convert_to_mixture_variables(q_cons_vf, rho, & - gamma, pi_inf, & - Re, j, k, l) + call s_convert_to_mixture_variables(q_cons_vf, j, k, l, & + rho, gamma, pi_inf, Re) do s = 1, num_dims vel(s) = q_cons_vf(cont_idx%end + s)%sf(j, k, l)/rho end do @@ -2330,9 +2326,8 @@ subroutine s_write_probe_files(t_step, q_cons_vf, accel_mag) ! ----------- if (trigger) then npts = npts + 1 - call s_convert_to_mixture_variables(q_cons_vf, rho, & - gamma, pi_inf, & - Re, j, k, l) + call s_convert_to_mixture_variables(q_cons_vf, j, k, l, & + rho, gamma, pi_inf, Re) do s = 1, num_dims vel(s) = q_cons_vf(cont_idx%end + s)%sf(j, k, l)/rho end do diff --git a/src/simulation/m_fftw.fpp b/src/simulation/m_fftw.fpp index 1e15dba44..cb1bd253e 100644 --- a/src/simulation/m_fftw.fpp +++ b/src/simulation/m_fftw.fpp @@ -17,7 +17,6 @@ module m_fftw #if defined(_OPENACC) && defined(__PGI) use cufft #endif - ! ========================================================================== implicit none diff --git a/src/simulation/m_global_parameters.fpp b/src/simulation/m_global_parameters.fpp index a4e8589c9..534400ddd 100644 --- a/src/simulation/m_global_parameters.fpp +++ b/src/simulation/m_global_parameters.fpp @@ -818,14 +818,6 @@ contains !$acc update device(momxb, momxe, advxb, advxe, contxb, contxe, bubxb, bubxe, intxb, intxe, sys_size, buff_size, E_idx, alf_idx, strxb, strxe) - allocate (gammas(1:num_fluids), pi_infs(1:num_fluids)) - - do i = 1, num_fluids - gammas(i) = fluid_pp(i)%gamma - pi_infs(i) = fluid_pp(i)%pi_inf - end do -!$acc update device(gammas, pi_infs) - ! Allocating grid variables for the x-, y- and z-directions allocate (x_cb(-1 - buff_size:m + buff_size)) allocate (x_cc(-buff_size:m + buff_size)) diff --git a/src/simulation/m_riemann_solvers.fpp b/src/simulation/m_riemann_solvers.fpp index 1492f3738..c8c264366 100644 --- a/src/simulation/m_riemann_solvers.fpp +++ b/src/simulation/m_riemann_solvers.fpp @@ -4660,7 +4660,6 @@ contains deallocate (V0_L, V0_R) end if - deallocate (gammas, pi_infs, Gs) ! Disassociating procedural pointer to the subroutine which was ! utilized to calculate the solution of a given Riemann problem s_riemann_solver => null() diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index d5c835d81..15f9c5977 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -1118,7 +1118,7 @@ contains do k = 0, n do l = 0, p - call s_convert_to_mixture_variables(v_vf, rho, gamma, pi_inf, Re, j, k, l) + call s_convert_to_mixture_variables(v_vf, j, k, l, rho, gamma, pi_inf, Re) dyn_pres = 0d0 do i = mom_idx%beg, mom_idx%end