From 8737ba96376f0f2c78cb3588052cc3c320960a2c Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Fri, 2 Oct 2026 23:38:37 -0500 Subject: [PATCH 1/5] Replace the USING_AMD array guards with BOUND() and a fixed-bounds switch Per-thread arrays in GPU kernels get compile-time extents because runtime-sized private arrays spill to scratch (4-5x on amdflang). The 83 copy-pasted `#:if [not MFC_CASE_OPTIMIZATION and] USING_AMD` blocks become single declarations using ${BOUND('name')}$, which returns the fixed maximum from one table (num_fluids 3, nb 3, sys_size 70 or 10 + species, ...) or the runtime extent. A converter checked every literal against that table. MFC_FIXED_BOUNDS (CMake, default ON) is decided per target: only the offloaded simulation uses fixed bounds, so pre/post_process callers always match the routines they call. This commit keeps it amdflang-only, and the generated code for amdflang simulation and for every non-fixed build matches master. s_check_amd becomes s_check_fixed_bounds and absorbs the CBC limits. --- CMakeLists.txt | 1 + cmake/Fypp.cmake | 8 + src/common/include/shared_parallel_macros.fpp | 18 +- src/common/m_checker_common.fpp | 19 +- src/common/m_eos.fpp | 78 ++---- src/common/m_phase_change.fpp | 13 +- src/common/m_variables_conversion.fpp | 135 ++++------ src/simulation/m_acoustic_src.fpp | 43 ++-- src/simulation/m_bubbles_EE.fpp | 18 +- src/simulation/m_bubbles_EL.fpp | 26 +- src/simulation/m_cbc.fpp | 60 ++--- src/simulation/m_compute_cbc.fpp | 240 +++++------------- src/simulation/m_data_output.fpp | 52 ++-- src/simulation/m_ibm.fpp | 59 ++--- src/simulation/m_pressure_relaxation.fpp | 30 +-- src/simulation/m_qbmm.fpp | 86 ++----- src/simulation/m_riemann_solver_hll.fpp | 19 +- src/simulation/m_riemann_solver_hllc.fpp | 150 +++++------ src/simulation/m_riemann_solver_hlld.fpp | 18 +- src/simulation/m_riemann_solver_hypo_hlld.fpp | 98 ++++--- src/simulation/m_riemann_solver_lf.fpp | 18 +- src/simulation/m_riemann_state.fpp | 112 +++----- src/simulation/m_sim_helpers.fpp | 77 +++--- src/simulation/m_start_up.fpp | 12 +- src/simulation/m_surface_tension.fpp | 29 +-- src/simulation/m_time_steppers.fpp | 36 ++- src/simulation/m_viscous.fpp | 29 +-- src/simulation/m_weno.fpp | 34 +-- toolchain/mfc/lint_source.py | 2 +- 29 files changed, 572 insertions(+), 948 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 5495f62088..5d3443d550 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -20,6 +20,7 @@ project(MFC LANGUAGES C CXX Fortran) option(MFC_MPI "Build with MPI" ON) option(MFC_OpenACC "Build with OpenACC" OFF) option(MFC_OpenMP "Build with OpenMP" OFF) +option(MFC_FIXED_BOUNDS "Compile-time bounds for per-thread GPU arrays" ON) option(MFC_GCov "Build with GCov" OFF) option(MFC_Unified "Build with unified CPU & GPU memory (GH-200 only)" OFF) option(MFC_Fastmath "Build with -gpu=fastmath on NV GPUs" OFF) diff --git a/cmake/Fypp.cmake b/cmake/Fypp.cmake index 03cc237af7..4a4b3a24b8 100644 --- a/cmake/Fypp.cmake +++ b/cmake/Fypp.cmake @@ -101,6 +101,13 @@ macro(HANDLE_SOURCES target useCommon) list(APPEND ${target}_incs ${common_incs}) endif() + # Fixed per-thread array bounds (see shared_parallel_macros.fpp): only simulation is offloaded. + set(_fixed_bounds False) + if (MFC_FIXED_BOUNDS AND (MFC_OpenACC OR MFC_OpenMP) AND "${target}" STREQUAL "simulation" + AND CMAKE_Fortran_COMPILER_ID STREQUAL "LLVMFlang") + set(_fixed_bounds True) + endif() + # /path/to/*.fpp (used by ) -> /fypp//*.f90 file(MAKE_DIRECTORY "${CMAKE_BINARY_DIR}/fypp/${target}") foreach(fpp ${${target}_FPPs}) @@ -118,6 +125,7 @@ macro(HANDLE_SOURCES target useCommon) -D MFC_${CMAKE_Fortran_COMPILER_ID} -D MFC_${${target}_UPPER} -D MFC_COMPILER="${CMAKE_Fortran_COMPILER_ID}" + -D MFC_FIXED_BOUNDS=${_fixed_bounds} -D MFC_CASE_OPTIMIZATION=False -D chemistry=False --line-numbering diff --git a/src/common/include/shared_parallel_macros.fpp b/src/common/include/shared_parallel_macros.fpp index 8c05cd1dec..1bd13b650e 100644 --- a/src/common/include/shared_parallel_macros.fpp +++ b/src/common/include/shared_parallel_macros.fpp @@ -11,9 +11,23 @@ #! NUM_SPECIES (m_thermochem's species count) and CHEMISTRY, written per build by the toolchain. #:include 'thermochem.fpp' -#! Fallback extent the USING_AMD guards substitute for sys_size arrays when case optimization is off. +#! Fixed bounds: in GPU simulation builds (MFC_FIXED_BOUNDS, set per target by CMake) per-thread +#! arrays get compile-time extents, since runtime-sized private arrays spill to scratch. Case +#! optimization makes most of these compile-time anyway; sys_size never is, so it is always fixed. #! Chemistry pins num_fluids to 1, leaving at most 10 flow variables beside the species. -#:set AMD_SYS_SIZE_MAX = 10 + NUM_SPECIES if CHEMISTRY else 70 +#:set MFC_FIXED_BOUNDS = defined('MFC_FIXED_BOUNDS') and MFC_FIXED_BOUNDS +#:set FIXED_BOUNDS = MFC_FIXED_BOUNDS and not MFC_CASE_OPTIMIZATION +#:set NUM_FLUIDS_MAX = 3 +#:set NB_MAX = 3 +#:set SYS_SIZE_MAX = 10 + NUM_SPECIES if CHEMISTRY else 70 +#:set BOUND_MAX = {'num_fluids': NUM_FLUIDS_MAX, 'nb': NB_MAX, 'num_dims': 3, 'num_vels': 3, & + & 'weno_polyn': 3, 'weno_num_stencils': 4, 'nterms': 32, 'n_stress': 6, 'sys_size': SYS_SIZE_MAX} +#:set BOUND_RUNTIME = {'n_stress': 'eqn_idx%stress%end - eqn_idx%stress%beg + 1'} + +#! Extent for a per-thread array sized by `name`: its fixed maximum or the runtime expression. +#:def BOUND(name) + $:BOUND_MAX[name] if (MFC_FIXED_BOUNDS if name == 'sys_size' else FIXED_BOUNDS) else BOUND_RUNTIME.get(name, name) +#:enddef #:def ASSERT_LIST(data, datatype) #:assert data is not None diff --git a/src/common/m_checker_common.fpp b/src/common/m_checker_common.fpp index 83e1adf4aa..40b2f49186 100644 --- a/src/common/m_checker_common.fpp +++ b/src/common/m_checker_common.fpp @@ -26,8 +26,8 @@ contains integer(kind=8), intent(in) :: n_global if (check_total_cells) call s_check_total_cells(n_global) - #:if USING_AMD - call s_check_amd + #:if MFC_FIXED_BOUNDS + call s_check_fixed_bounds #:endif end subroutine s_check_inputs_common @@ -48,17 +48,16 @@ contains end subroutine s_check_total_cells - !> Check that simulation parameters stay within AMD GPU compiler limits when case optimization is disabled. - impure subroutine s_check_amd + !> Check that the case fits the fixed per-thread array bounds of a GPU simulation build. + impure subroutine s_check_fixed_bounds #:if not MFC_CASE_OPTIMIZATION - @:PROHIBIT(num_fluids > 3, "num_fluids <= 3 for AMDFLang when Case optimization is off") - @:PROHIBIT((bubbles_euler .or. bubbles_lagrange) .and. nb > 3, "nb <= 3 for AMDFLang when Case optimization is off") - ! HLLC's star states and the CBC L vectors have no other bound. - @:PROHIBIT(sys_size > ${AMD_SYS_SIZE_MAX}$, & - & "sys_size <= ${AMD_SYS_SIZE_MAX}$ for AMDFLang when Case optimization is off") + @:PROHIBIT(num_fluids > ${NUM_FLUIDS_MAX}$, & + & "num_fluids <= ${NUM_FLUIDS_MAX}$ in GPU builds; rebuild with --case-optimization") + @:PROHIBIT((bubbles_euler .or. bubbles_lagrange) .and. nb > ${NB_MAX}$, & + & "nb <= ${NB_MAX}$ in GPU builds; rebuild with --case-optimization") #:endif - end subroutine s_check_amd + end subroutine s_check_fixed_bounds end module m_checker_common diff --git a/src/common/m_eos.fpp b/src/common/m_eos.fpp index 16e126297b..ea8721f390 100644 --- a/src/common/m_eos.fpp +++ b/src/common/m_eos.fpp @@ -305,15 +305,11 @@ contains $:GPU_ROUTINE(function_name='f_mixture_temperature', parallelism='[seq]', cray_inline=True) - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3), intent(in) :: alpha_rho_K - #:else - real(wp), dimension(num_fluids), intent(in) :: alpha_rho_K - #:endif - real(wp), intent(in) :: pres, gamma_K, pi_inf_K - real(wp) :: T - real(wp) :: mCP !< sum of alpha_rho_i*cp_i; cp_i = n_i*cv_i - integer :: i + real(wp), dimension(${BOUND('num_fluids')}$), intent(in) :: alpha_rho_K + real(wp), intent(in) :: pres, gamma_K, pi_inf_K + real(wp) :: T + real(wp) :: mCP !< sum of alpha_rho_i*cp_i; cp_i = n_i*cv_i + integer :: i mCP = 0._wp $:GPU_LOOP(parallelism='[seq]') @@ -569,15 +565,11 @@ contains $:GPU_ROUTINE(function_name='s_compute_mixture_coefficients', parallelism='[seq]', cray_inline=True) - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3), intent(in) :: alpha_rho_K, alpha_K - #:else - real(wp), dimension(num_fluids), intent(in) :: alpha_rho_K, alpha_K - #:endif - real(wp), intent(out) :: rho_K, gamma_K, pi_inf_K, qv_K - real(wp) :: gamma_i, pi_inf_i, dpi_i, dgamma_i - real(wp) :: rho_i, alpha_i, alpha_rho_i - integer :: i !< Loop iterator over fluids + real(wp), dimension(${BOUND('num_fluids')}$), intent(in) :: alpha_rho_K, alpha_K + real(wp), intent(out) :: rho_K, gamma_K, pi_inf_K, qv_K + real(wp) :: gamma_i, pi_inf_i, dpi_i, dgamma_i + real(wp) :: rho_i, alpha_i, alpha_rho_i + integer :: i !< Loop iterator over fluids ! The bubbly closure is written for one carrier liquid, which keeps its own coefficients ! undiluted: Gamma_l*p_l = (E - rho|u|^2/2)/(1 - alf) - Pi_inf_l, the void entering only through @@ -614,14 +606,10 @@ contains $:GPU_ROUTINE(function_name='s_compute_mixture_coefficients_dt', parallelism='[seq]', cray_inline=True) - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3), intent(in) :: dalpha_rho_dt, dadv_dt, alpha_rho, adv - #:else - real(wp), dimension(num_fluids), intent(in) :: dalpha_rho_dt, dadv_dt, alpha_rho, adv - #:endif - real(wp), intent(out) :: drho_dt, dgamma_dt, dpi_inf_dt, dqv_dt - real(wp) :: rho_i, gamma_i, pi_inf_i, dpi_i, dgamma_i, alpha_i, alpha_rho_i - integer :: i !< Loop iterator over fluids + real(wp), dimension(${BOUND('num_fluids')}$), intent(in) :: dalpha_rho_dt, dadv_dt, alpha_rho, adv + real(wp), intent(out) :: drho_dt, dgamma_dt, dpi_inf_dt, dqv_dt + real(wp) :: rho_i, gamma_i, pi_inf_i, dpi_i, dgamma_i, alpha_i, alpha_rho_i + integer :: i !< Loop iterator over fluids dgamma_dt = 0._wp dpi_inf_dt = 0._wp @@ -654,21 +642,13 @@ contains $:GPU_ROUTINE(parallelism='[seq]') - real(wp), intent(in) :: pres, rho, gamma, pi_inf - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3), intent(in) :: adv - #:else - real(wp), dimension(num_fluids), intent(in) :: adv - #:endif - real(wp), intent(out) :: c - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3), intent(in), optional :: alpha_rho - #:else - real(wp), dimension(num_fluids), intent(in), optional :: alpha_rho - #:endif - real(wp) :: alf !< Subgrid void fraction; dilute by construction - real(wp) :: blkmod_q, alpha_q, alpha_rho_q, gamma_q, pi_inf_q - integer :: q + real(wp), intent(in) :: pres, rho, gamma, pi_inf + real(wp), dimension(${BOUND('num_fluids')}$), intent(in) :: adv + real(wp), intent(out) :: c + real(wp), dimension(${BOUND('num_fluids')}$), intent(in), optional :: alpha_rho + real(wp) :: alf !< Subgrid void fraction; dilute by construction + real(wp) :: blkmod_q, alpha_q, alpha_rho_q, gamma_q, pi_inf_q + integer :: q if (chemistry) then ! Reacting mixture sound speed c = sqrt((1.0_wp + 1.0_wp/gamma)*pres/rho) @@ -741,18 +721,10 @@ contains $:GPU_ROUTINE(parallelism='[seq]') - real(wp), intent(in) :: pres, rho, gamma, pi_inf, qv, vel_sum, H, c_c - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3), intent(in) :: adv - #:else - real(wp), dimension(num_fluids), intent(in) :: adv - #:endif - real(wp), intent(out) :: c - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3), intent(in), optional :: alpha_rho - #:else - real(wp), dimension(num_fluids), intent(in), optional :: alpha_rho - #:endif + real(wp), intent(in) :: pres, rho, gamma, pi_inf, qv, vel_sum, H, c_c + real(wp), dimension(${BOUND('num_fluids')}$), intent(in) :: adv + real(wp), intent(out) :: c + real(wp), dimension(${BOUND('num_fluids')}$), intent(in), optional :: alpha_rho if (chemistry) then ! Reacting mixture sound speed if (avg_state == avg_state_roe .and. abs(c_c) > verysmall) then diff --git a/src/common/m_phase_change.fpp b/src/common/m_phase_change.fpp index c0c6c81f1e..582da96325 100644 --- a/src/common/m_phase_change.fpp +++ b/src/common/m_phase_change.fpp @@ -59,12 +59,7 @@ contains real(wp) :: rhoe, dynE, rhos !< total internal energy, kinetic energy, and total entropy real(wp) :: rho, rM, m1, m2, MCT !< total density, total reacting mass, individual reacting masses real(wp) :: TvF !< total volume fraction - - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3) :: p_infpT, sk, hk, gk, ek, rhok - #:else - real(wp), dimension(num_fluids) :: p_infpT, sk, hk, gk, ek, rhok - #:endif + real(wp), dimension(${BOUND('num_fluids')}$) :: p_infpT, sk, hk, gk, ek, rhok !> Generic loop iterators integer :: i, j, k, l @@ -214,10 +209,10 @@ contains mQ = mQ + q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, l)*qvs(i) end do - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - if (num_fluids < 3) then + #:if FIXED_BOUNDS + if (num_fluids < ${NUM_FLUIDS_MAX}$) then $:GPU_LOOP(parallelism='[seq]') - do i = num_fluids + 1, 3 + do i = num_fluids + 1, ${NUM_FLUIDS_MAX}$ p_infpT(i) = p_infpT_sum end do end if diff --git a/src/common/m_variables_conversion.fpp b/src/common/m_variables_conversion.fpp index 87c3d2654c..0fabe8af0a 100644 --- a/src/common/m_variables_conversion.fpp +++ b/src/common/m_variables_conversion.fpp @@ -186,18 +186,13 @@ contains $:GPU_ROUTINE(function_name='s_convert_species_to_mixture_variables_kernel', parallelism='[seq]', cray_noinline=True) - real(wp), intent(out) :: rho_K, gamma_K, pi_inf_K, qv_K - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3), intent(inout) :: alpha_rho_K, alpha_K - real(wp), optional, dimension(3), intent(in) :: G - #:else - real(wp), dimension(num_fluids), intent(inout) :: alpha_rho_K, alpha_K - real(wp), optional, dimension(num_fluids), intent(in) :: G - #:endif - real(wp), optional, dimension(2), intent(out) :: Re_K - real(wp), optional, intent(out) :: G_K - real(wp) :: alpha_K_sum - integer :: i, j !< Generic loop iterators + real(wp), intent(out) :: rho_K, gamma_K, pi_inf_K, qv_K + real(wp), dimension(${BOUND('num_fluids')}$), intent(inout) :: alpha_rho_K, alpha_K + real(wp), optional, dimension(${BOUND('num_fluids')}$), intent(in) :: G + real(wp), optional, dimension(2), intent(out) :: Re_K + real(wp), optional, intent(out) :: G_K + real(wp) :: alpha_K_sum + integer :: i, j !< Generic loop iterators rho_K = 0._wp gamma_K = 0._wp @@ -385,37 +380,31 @@ contains use m_global_parameters_common, only: shear_indices ! Performance fix with AMDFlang - type(scalar_field), dimension(sys_size), intent(in) :: qK_cons_vf - type(scalar_field), intent(inout) :: q_T_sf + type(scalar_field), dimension(sys_size), intent(in) :: qK_cons_vf + type(scalar_field), intent(inout) :: q_T_sf type(scalar_field), dimension(sys_size), intent(inout) :: qK_prim_vf - type(int_bounds_info), dimension(1:3), intent(in) :: ibounds - - #:if USING_AMD and not MFC_CASE_OPTIMIZATION - real(wp), dimension(3) :: alpha_K, alpha_rho_K - real(wp), dimension(3) :: nRtmp - #:else - real(wp), dimension(num_fluids) :: alpha_K, alpha_rho_K - real(wp), dimension(nb) :: nRtmp - #:endif - real(wp) :: rhoYks(1:${NUM_SPECIES}$) + type(int_bounds_info), dimension(1:3), intent(in) :: ibounds + real(wp), dimension(${BOUND('num_fluids')}$) :: alpha_K, alpha_rho_K + real(wp), dimension(${BOUND('nb')}$) :: nRtmp + real(wp) :: rhoYks(1:${NUM_SPECIES}$) real(wp), dimension(2) :: Re_K - real(wp) :: rho_K, gamma_K, pi_inf_K, qv_K, dyn_pres_K - real(wp) :: vftmp, nbub_sc - real(wp) :: G_K - real(wp) :: solid_partial_density - real(wp) :: pres - integer :: i, j, k, l !< Generic loop iterators - real(wp) :: T - real(wp) :: pres_mag - real(wp) :: Ga !< Lorentz factor (gamma in relativity) - real(wp) :: B2 !< Magnetic field magnitude squared - real(wp) :: B(3) !< Magnetic field components - real(wp) :: m2 !< Relativistic momentum magnitude squared - real(wp) :: S !< Dot product of the magnetic field and the relativistic momentum - real(wp) :: W, dW !< W := rho*v*Ga**2; f = f(W) in Newton-Raphson - real(wp) :: E, D !< Prim/Cons variables within Newton-Raphson iteration - real(wp) :: f, dGa_dW, dp_dW, df_dW !< Functions within Newton-Raphson iteration - integer :: iter !< Newton-Raphson iteration counter + real(wp) :: rho_K, gamma_K, pi_inf_K, qv_K, dyn_pres_K + real(wp) :: vftmp, nbub_sc + real(wp) :: G_K + real(wp) :: solid_partial_density + real(wp) :: pres + integer :: i, j, k, l !< Generic loop iterators + real(wp) :: T + real(wp) :: pres_mag + real(wp) :: Ga !< Lorentz factor (gamma in relativity) + real(wp) :: B2 !< Magnetic field magnitude squared + real(wp) :: B(3) !< Magnetic field components + real(wp) :: m2 !< Relativistic momentum magnitude squared + real(wp) :: S !< Dot product of the magnetic field and the relativistic momentum + real(wp) :: W, dW !< W := rho*v*Ga**2; f = f(W) in Newton-Raphson + real(wp) :: E, D !< Prim/Cons variables within Newton-Raphson iteration + real(wp) :: f, dGa_dW, dp_dW, df_dW !< Functions within Newton-Raphson iteration + integer :: iter !< Newton-Raphson iteration counter $:GPU_PARALLEL_LOOP(collapse=3, private='[alpha_K, alpha_rho_K, Re_K, nRtmp, rho_K, gamma_K, pi_inf_K, qv_K, dyn_pres_K, & & rhoYks, B, pres, vftmp, nbub_sc, G_K, solid_partial_density, T, pres_mag, Ga, B2, m2, S, W, dW, E, & @@ -946,28 +935,22 @@ contains ! Partial densities, density, velocity, pressure, energy, advection variables, the specific heat ratio and liquid stiffness ! functions, the shear and volume Reynolds numbers and the Weber numbers - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3) :: alpha_rho_K - real(wp), dimension(3) :: alpha_K - real(wp), dimension(3) :: vel_K - #:else - real(wp), dimension(num_fluids) :: alpha_rho_K - real(wp), dimension(num_fluids) :: alpha_K - real(wp), dimension(num_vels) :: vel_K - #:endif - real(wp), dimension(${NUM_SPECIES}$) :: Y_K - real(wp) :: rho_K - real(wp) :: vel_K_sum - real(wp) :: pres_K - real(wp) :: E_K - real(wp) :: gamma_K - real(wp) :: pi_inf_K - real(wp) :: qv_K - real(wp), dimension(2) :: Re_K - real(wp) :: G_K - real(wp) :: blkmod1_K, blkmod2_K, K_K - real(wp) :: T_K, mix_mol_weight, R_gas - integer :: i, j, k, l !< Generic loop iterators + real(wp), dimension(${BOUND('num_fluids')}$) :: alpha_rho_K + real(wp), dimension(${BOUND('num_fluids')}$) :: alpha_K + real(wp), dimension(${BOUND('num_vels')}$) :: vel_K + real(wp), dimension(${NUM_SPECIES}$) :: Y_K + real(wp) :: rho_K + real(wp) :: vel_K_sum + real(wp) :: pres_K + real(wp) :: E_K + real(wp) :: gamma_K + real(wp) :: pi_inf_K + real(wp) :: qv_K + real(wp), dimension(2) :: Re_K + real(wp) :: G_K + real(wp) :: blkmod1_K, blkmod2_K, K_K + real(wp) :: T_K, mix_mol_weight, R_gas + integer :: i, j, k, l !< Generic loop iterators is1b = is1%beg; is1e = is1%end is2b = is2%beg; is2e = is2%end @@ -1119,15 +1102,11 @@ contains subroutine s_compute_species_fraction(q_vf, k, l, r, alpha_rho_K, alpha_K) $:GPU_ROUTINE(function_name='s_compute_species_fraction', parallelism='[seq]', cray_noinline=True) - type(scalar_field), dimension(sys_size), intent(in) :: q_vf - integer, intent(in) :: k, l, r - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3), intent(out) :: alpha_rho_K, alpha_K - #:else - real(wp), dimension(num_fluids), intent(out) :: alpha_rho_K, alpha_K - #:endif - integer :: i - real(wp) :: alpha_K_sum + type(scalar_field), dimension(sys_size), intent(in) :: q_vf + integer, intent(in) :: k, l, r + real(wp), dimension(${BOUND('num_fluids')}$), intent(out) :: alpha_rho_K, alpha_K + integer :: i + real(wp) :: alpha_K_sum if (num_fluids == 1) then alpha_rho_K(1) = q_vf(eqn_idx%cont%beg)%sf(k, l, r) @@ -1187,14 +1166,10 @@ contains $:GPU_ROUTINE(function_name='s_compute_energy', parallelism='[seq]', cray_inline=True) - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3), intent(in) :: alpha_rho_K, alpha_K - #:else - real(wp), dimension(num_fluids), intent(in) :: alpha_rho_K, alpha_K - #:endif - real(wp), intent(in) :: pres, vel_sum - real(wp), intent(out) :: E - real(wp) :: rho, gamma, pi_inf, qv + real(wp), dimension(${BOUND('num_fluids')}$), intent(in) :: alpha_rho_K, alpha_K + real(wp), intent(in) :: pres, vel_sum + real(wp), intent(out) :: E + real(wp) :: rho, gamma, pi_inf, qv call s_compute_mixture_coefficients(alpha_rho_K, alpha_K, rho, gamma, pi_inf, qv) diff --git a/src/simulation/m_acoustic_src.fpp b/src/simulation/m_acoustic_src.fpp index 5ca556ede1..71d453b80c 100644 --- a/src/simulation/m_acoustic_src.fpp +++ b/src/simulation/m_acoustic_src.fpp @@ -126,31 +126,26 @@ contains !> Compute mass, momentum, and energy acoustic source terms and add to the RHS impure subroutine s_acoustic_src_calculations(q_cons_vf, q_prim_vf, rhs_vf) - type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf !< Conservative variables - type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf !< Primitive variables + type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf !< Conservative variables + type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf !< Primitive variables type(scalar_field), dimension(sys_size), intent(inout) :: rhs_vf - - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3) :: myalpha, myalpha_rho - #:else - real(wp), dimension(num_fluids) :: myalpha, myalpha_rho - #:endif - real(wp) :: myRho, pi_inf_mix, qv_dummy - real(wp) :: sim_time, c, gamma_mix - real(wp) :: blkmod_q, pres_q, alpha_q, alpha_rho_q - real(wp) :: frequency_local, gauss_sigma_time_local - real(wp) :: mass_src_diff, mom_src_diff - real(wp) :: source_temporal - real(wp) :: period_BB !< period of each sine wave in broadband source - real(wp) :: sl_BB !< spectral level at each frequency - real(wp) :: ffre_BB !< source term corresponding to each frequency - real(wp) :: sum_BB !< total source term for the broadband wave - real(wp), allocatable, dimension(:) :: phi_rn !< random phase shift for each frequency - integer :: i, j, k, l, q !< generic loop variables - integer :: ai !< acoustic source index - integer :: num_points - logical :: freq_conv_flag, gauss_conv_flag - integer, parameter :: mass_label = 1, mom_label = 2 + real(wp), dimension(${BOUND('num_fluids')}$) :: myalpha, myalpha_rho + real(wp) :: myRho, pi_inf_mix, qv_dummy + real(wp) :: sim_time, c, gamma_mix + real(wp) :: blkmod_q, pres_q, alpha_q, alpha_rho_q + real(wp) :: frequency_local, gauss_sigma_time_local + real(wp) :: mass_src_diff, mom_src_diff + real(wp) :: source_temporal + real(wp) :: period_BB !< period of each sine wave in broadband source + real(wp) :: sl_BB !< spectral level at each frequency + real(wp) :: ffre_BB !< source term corresponding to each frequency + real(wp) :: sum_BB !< total source term for the broadband wave + real(wp), allocatable, dimension(:) :: phi_rn !< random phase shift for each frequency + integer :: i, j, k, l, q !< generic loop variables + integer :: ai !< acoustic source index + integer :: num_points + logical :: freq_conv_flag, gauss_conv_flag + integer, parameter :: mass_label = 1, mom_label = 2 sim_time = mytime ! Accumulated time, correct under adaptive dt diff --git a/src/simulation/m_bubbles_EE.fpp b/src/simulation/m_bubbles_EE.fpp index c1f9547422..bc0f9e15a0 100644 --- a/src/simulation/m_bubbles_EE.fpp +++ b/src/simulation/m_bubbles_EE.fpp @@ -149,19 +149,13 @@ contains real(wp) :: pb_local, mv_local, vflux, pbdot real(wp) :: n_tait, B_tait, qv_bub real(wp) :: chi_vw_l, k_mw_l, rho_mw_l !< Per-thread bubble-wall scratch (avoid module-scalar race) - - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3) :: Rtmp, Vtmp - real(wp), dimension(3) :: myalpha, myalpha_rho - #:else - real(wp), dimension(nb) :: Rtmp, Vtmp - real(wp), dimension(num_fluids) :: myalpha, myalpha_rho - #:endif + real(wp), dimension(${BOUND('nb')}$) :: Rtmp, Vtmp + real(wp), dimension(${BOUND('num_fluids')}$) :: myalpha, myalpha_rho real(wp) :: myR, myV, alf, myP, myRho, R2Vav, R3 - real(wp) :: nbub !< Bubble number density - integer :: i, j, k, l, q, ii !< Loop variables - integer :: adap_dt_stop_sum, adap_dt_stop !< Fail-safe exit if max iteration count reached - integer :: dmBub_id !< Dummy variables for unified subgrid bubble subroutines + real(wp) :: nbub !< Bubble number density + integer :: i, j, k, l, q, ii !< Loop variables + integer :: adap_dt_stop_sum, adap_dt_stop !< Fail-safe exit if max iteration count reached + integer :: dmBub_id !< Dummy variables for unified subgrid bubble subroutines real(wp) :: dmMass_v, dmMass_n, dmBeta_c, dmBeta_t, dmCson $:GPU_PARALLEL_LOOP(private='[j, k, l, q]', collapse=3) diff --git a/src/simulation/m_bubbles_EL.fpp b/src/simulation/m_bubbles_EL.fpp index 0f40e6fe6f..7de920f933 100644 --- a/src/simulation/m_bubbles_EL.fpp +++ b/src/simulation/m_bubbles_EL.fpp @@ -576,22 +576,17 @@ contains !> Contains the bubble dynamics subroutines. subroutine s_compute_bubble_EL_dynamics(q_prim_vf, bc_type, stage) - type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf + type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf type(integer_field), dimension(1:num_dims,1:2), intent(in) :: bc_type - integer, intent(in) :: stage - real(wp) :: myVapFlux - real(wp) :: preterm1, term2, paux, pint, Romega, term1_fac - real(wp) :: myR_m, mygamma_m, myPb, myMass_n, myMass_v - real(wp) :: myR, myV, myBeta_c, myBeta_t, myR0, myPbdot, myMvdot - real(wp) :: myPinf, aux1, aux2, myCson, myRho - real(wp), dimension(3) :: myPos, myVel - real(wp) :: gamma, pi_inf, qv, f_b, myRe - - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3) :: myalpha_rho, myalpha - #:else - real(wp), dimension(num_fluids) :: myalpha_rho, myalpha - #:endif + integer, intent(in) :: stage + real(wp) :: myVapFlux + real(wp) :: preterm1, term2, paux, pint, Romega, term1_fac + real(wp) :: myR_m, mygamma_m, myPb, myMass_n, myMass_v + real(wp) :: myR, myV, myBeta_c, myBeta_t, myR0, myPbdot, myMvdot + real(wp) :: myPinf, aux1, aux2, myCson, myRho + real(wp), dimension(3) :: myPos, myVel + real(wp) :: gamma, pi_inf, qv, f_b, myRe + real(wp), dimension(${BOUND('num_fluids')}$) :: myalpha_rho, myalpha real(wp), dimension(2) :: Re integer, dimension(3) :: cell integer :: adap_dt_stop_sum, adap_dt_stop !< Fail-safe exit if max iteration count reached @@ -599,6 +594,7 @@ contains integer :: k, l ! Subgrid p_inf model based on Maeda and Colonius (2018). + if (lag_params%pressure_corrector) then call nvtxStartRange("LAGRANGE-BUBBLE-PINF-CORRECTION") ! Calculate velocity potentials (valid for one bubble per cell) diff --git a/src/simulation/m_cbc.fpp b/src/simulation/m_cbc.fpp index 49cc12c707..4884ce3f91 100644 --- a/src/simulation/m_cbc.fpp +++ b/src/simulation/m_cbc.fpp @@ -475,43 +475,29 @@ contains real(wp) :: dpi_inf_dt real(wp) :: dqv_dt real(wp) :: dpres_ds - real(wp) :: ramp !< inflow ramp factor; unity unless a ramp is set - - #:if USING_AMD - real(wp), dimension(${AMD_SYS_SIZE_MAX}$) :: L - #:else - real(wp), dimension(sys_size) :: L - #:endif - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3) :: alpha_rho, dalpha_rho_ds, mf - real(wp), dimension(3) :: vel, dvel_ds - real(wp), dimension(3) :: adv_local, dadv_ds - real(wp), dimension(3) :: dadv_dt - real(wp), dimension(3) :: dvel_dt - real(wp), dimension(3) :: dalpha_rho_dt - #:else - real(wp), dimension(num_fluids) :: alpha_rho, dalpha_rho_ds, mf - real(wp), dimension(num_vels) :: vel, dvel_ds - real(wp), dimension(num_fluids) :: adv_local, dadv_ds - real(wp), dimension(num_fluids) :: dadv_dt - real(wp), dimension(num_dims) :: dvel_dt - real(wp), dimension(num_fluids) :: dalpha_rho_dt - #:endif - real(wp), dimension(${NUM_SPECIES}$) :: Ys, h_k, dYs_dt, dYs_ds, Xs, Gamma_i, Cp_i - real(wp), dimension(2) :: Re_cbc - real(wp), dimension(3) :: lambda - real(wp) :: rho !< Cell averaged density - real(wp) :: pres !< Cell averaged pressure - real(wp) :: E !< Cell averaged energy - real(wp) :: gamma !< Cell averaged specific heat ratio - real(wp) :: pi_inf !< Cell averaged liquid stiffness - real(wp) :: qv !< Cell averaged fluid reference energy - real(wp) :: c - real(wp) :: Ma - real(wp) :: T, sum_Enthalpies - real(wp) :: Cv, Cp, e_mix, Mw, R_gas - real(wp) :: vel_K_sum, vel_dv_dt_sum - integer :: i, j, k, r !< Generic loop iterators + real(wp) :: ramp !< inflow ramp factor; unity unless a ramp is set + real(wp), dimension(${BOUND('sys_size')}$) :: L + real(wp), dimension(${BOUND('num_fluids')}$) :: alpha_rho, dalpha_rho_ds, mf + real(wp), dimension(${BOUND('num_vels')}$) :: vel, dvel_ds + real(wp), dimension(${BOUND('num_fluids')}$) :: adv_local, dadv_ds + real(wp), dimension(${BOUND('num_fluids')}$) :: dadv_dt + real(wp), dimension(${BOUND('num_dims')}$) :: dvel_dt + real(wp), dimension(${BOUND('num_fluids')}$) :: dalpha_rho_dt + real(wp), dimension(${NUM_SPECIES}$) :: Ys, h_k, dYs_dt, dYs_ds, Xs, Gamma_i, Cp_i + real(wp), dimension(2) :: Re_cbc + real(wp), dimension(3) :: lambda + real(wp) :: rho !< Cell averaged density + real(wp) :: pres !< Cell averaged pressure + real(wp) :: E !< Cell averaged energy + real(wp) :: gamma !< Cell averaged specific heat ratio + real(wp) :: pi_inf !< Cell averaged liquid stiffness + real(wp) :: qv !< Cell averaged fluid reference energy + real(wp) :: c + real(wp) :: Ma + real(wp) :: T, sum_Enthalpies + real(wp) :: Cv, Cp, e_mix, Mw, R_gas + real(wp) :: vel_K_sum, vel_dv_dt_sum + integer :: i, j, k, r !< Generic loop iterators ! Reshaping of inputted data and association of the FD and PI coefficients, or CBC coefficients, respectively, hinging on ! selected CBC coordinate direction diff --git a/src/simulation/m_compute_cbc.fpp b/src/simulation/m_compute_cbc.fpp index 0a24c524ed..0c539f0783 100644 --- a/src/simulation/m_compute_cbc.fpp +++ b/src/simulation/m_compute_cbc.fpp @@ -21,14 +21,10 @@ contains function f_base_L1(lambda, rho, c, dpres_ds, dvel_ds) result(L1) $:GPU_ROUTINE(parallelism='[seq]') - real(wp), dimension(3), intent(in) :: lambda - real(wp), intent(in) :: rho, c, dpres_ds - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3), intent(in) :: dvel_ds - #:else - real(wp), dimension(num_dims), intent(in) :: dvel_ds - #:endif - real(wp) :: L1 + real(wp), dimension(3), intent(in) :: lambda + real(wp), intent(in) :: rho, c, dpres_ds + real(wp), dimension(${BOUND('num_dims')}$), intent(in) :: dvel_ds + real(wp) :: L1 L1 = lambda(1)*(dpres_ds - rho*c*dvel_ds(dir_idx(1))) end function f_base_L1 @@ -37,19 +33,11 @@ contains subroutine s_fill_density_L(L, lambda_factor, lambda2, c, mf, dalpha_rho_ds, dpres_ds) $:GPU_ROUTINE(parallelism='[seq]') - #:if USING_AMD - real(wp), dimension(${AMD_SYS_SIZE_MAX}$), intent(inout) :: L - #:else - real(wp), dimension(sys_size), intent(inout) :: L - #:endif - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3), intent(in) :: mf, dalpha_rho_ds - #:else - real(wp), dimension(num_fluids), intent(in) :: mf, dalpha_rho_ds - #:endif - real(wp), intent(in) :: lambda_factor, lambda2, c - real(wp), intent(in) :: dpres_ds - integer :: i + real(wp), dimension(${BOUND('sys_size')}$), intent(inout) :: L + real(wp), dimension(${BOUND('num_fluids')}$), intent(in) :: mf, dalpha_rho_ds + real(wp), intent(in) :: lambda_factor, lambda2, c + real(wp), intent(in) :: dpres_ds + integer :: i do i = 2, eqn_idx%mom%beg L(i) = lambda_factor*lambda2*(c*c*dalpha_rho_ds(i - 1) - mf(i - 1)*dpres_ds) @@ -61,18 +49,10 @@ contains subroutine s_fill_velocity_L(L, lambda_factor, lambda2, dvel_ds) $:GPU_ROUTINE(parallelism='[seq]') - #:if USING_AMD - real(wp), dimension(${AMD_SYS_SIZE_MAX}$), intent(inout) :: L - #:else - real(wp), dimension(sys_size), intent(inout) :: L - #:endif - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3), intent(in) :: dvel_ds - #:else - real(wp), dimension(num_dims), intent(in) :: dvel_ds - #:endif - real(wp), intent(in) :: lambda_factor, lambda2 - integer :: i + real(wp), dimension(${BOUND('sys_size')}$), intent(inout) :: L + real(wp), dimension(${BOUND('num_dims')}$), intent(in) :: dvel_ds + real(wp), intent(in) :: lambda_factor, lambda2 + integer :: i do i = eqn_idx%mom%beg + 1, eqn_idx%mom%end L(i) = lambda_factor*lambda2*dvel_ds(dir_idx(i - eqn_idx%cont%end)) @@ -84,18 +64,10 @@ contains subroutine s_fill_advection_L(L, lambda_factor, lambda2, dadv_ds) $:GPU_ROUTINE(parallelism='[seq]') - #:if USING_AMD - real(wp), dimension(${AMD_SYS_SIZE_MAX}$), intent(inout) :: L - #:else - real(wp), dimension(sys_size), intent(inout) :: L - #:endif - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3), intent(in) :: dadv_ds - #:else - real(wp), dimension(num_fluids), intent(in) :: dadv_ds - #:endif - real(wp), intent(in) :: lambda_factor, lambda2 - integer :: i + real(wp), dimension(${BOUND('sys_size')}$), intent(inout) :: L + real(wp), dimension(${BOUND('num_fluids')}$), intent(in) :: dadv_ds + real(wp), intent(in) :: lambda_factor, lambda2 + integer :: i do i = eqn_idx%E, eqn_idx%adv%end - 1 L(i) = lambda_factor*lambda2*dadv_ds(i - eqn_idx%mom%end) @@ -107,14 +79,10 @@ contains subroutine s_fill_chemistry_L(L, lambda_factor, lambda2, dYs_ds) $:GPU_ROUTINE(parallelism='[seq]') - #:if USING_AMD - real(wp), dimension(${AMD_SYS_SIZE_MAX}$), intent(inout) :: L - #:else - real(wp), dimension(sys_size), intent(inout) :: L - #:endif - real(wp), dimension(${NUM_SPECIES}$), intent(in) :: dYs_ds - real(wp), intent(in) :: lambda_factor, lambda2 - integer :: i + real(wp), dimension(${BOUND('sys_size')}$), intent(inout) :: L + real(wp), dimension(${NUM_SPECIES}$), intent(in) :: dYs_ds + real(wp), intent(in) :: lambda_factor, lambda2 + integer :: i if (.not. chemistry) return @@ -129,18 +97,10 @@ contains $:GPU_ROUTINE(function_name='s_compute_slip_wall_L',parallelism='[seq]', cray_inline=True) - real(wp), dimension(3), intent(in) :: lambda - #:if USING_AMD - real(wp), dimension(${AMD_SYS_SIZE_MAX}$), intent(inout) :: L - #:else - real(wp), dimension(sys_size), intent(inout) :: L - #:endif - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3), intent(in) :: dvel_ds - #:else - real(wp), dimension(num_dims), intent(in) :: dvel_ds - #:endif - real(wp), intent(in) :: rho, c, dpres_ds + real(wp), dimension(3), intent(in) :: lambda + real(wp), dimension(${BOUND('sys_size')}$), intent(inout) :: L + real(wp), dimension(${BOUND('num_dims')}$), intent(in) :: dvel_ds + real(wp), intent(in) :: rho, c, dpres_ds L(1) = f_base_L1(lambda, rho, c, dpres_ds, dvel_ds) L(2:eqn_idx%adv%end - 1) = 0._wp @@ -153,25 +113,15 @@ contains $:GPU_ROUTINE(function_name='s_compute_nonreflecting_subsonic_buffer_L', parallelism='[seq]', cray_inline=True) - real(wp), dimension(3), intent(in) :: lambda - #:if USING_AMD - real(wp), dimension(${AMD_SYS_SIZE_MAX}$), intent(inout) :: L - #:else - real(wp), dimension(sys_size), intent(inout) :: L - #:endif - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3), intent(in) :: mf, dalpha_rho_ds - real(wp), dimension(3), intent(in) :: dvel_ds - real(wp), dimension(3), intent(in) :: dadv_ds - #:else - real(wp), dimension(num_fluids), intent(in) :: mf, dalpha_rho_ds - real(wp), dimension(num_dims), intent(in) :: dvel_ds - real(wp), dimension(num_fluids), intent(in) :: dadv_ds - #:endif - real(wp), dimension(${NUM_SPECIES}$), intent(in) :: dYs_ds - real(wp), intent(in) :: rho, c - real(wp), intent(in) :: dpres_ds - real(wp) :: lambda_factor + real(wp), dimension(3), intent(in) :: lambda + real(wp), dimension(${BOUND('sys_size')}$), intent(inout) :: L + real(wp), dimension(${BOUND('num_fluids')}$), intent(in) :: mf, dalpha_rho_ds + real(wp), dimension(${BOUND('num_dims')}$), intent(in) :: dvel_ds + real(wp), dimension(${BOUND('num_fluids')}$), intent(in) :: dadv_ds + real(wp), dimension(${NUM_SPECIES}$), intent(in) :: dYs_ds + real(wp), intent(in) :: rho, c + real(wp), intent(in) :: dpres_ds + real(wp) :: lambda_factor lambda_factor = (5.e-1_wp - 5.e-1_wp*sign(1._wp, lambda(1))) L(1) = lambda_factor*lambda(1)*(dpres_ds - rho*c*dvel_ds(dir_idx(1))) @@ -192,18 +142,10 @@ contains $:GPU_ROUTINE(function_name='s_compute_nonreflecting_subsonic_inflow_L', parallelism='[seq]', cray_inline=True) - real(wp), dimension(3), intent(in) :: lambda - #:if USING_AMD - real(wp), dimension(${AMD_SYS_SIZE_MAX}$), intent(inout) :: L - #:else - real(wp), dimension(sys_size), intent(inout) :: L - #:endif - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3), intent(in) :: dvel_ds - #:else - real(wp), dimension(num_dims), intent(in) :: dvel_ds - #:endif - real(wp), intent(in) :: rho, c, dpres_ds + real(wp), dimension(3), intent(in) :: lambda + real(wp), dimension(${BOUND('sys_size')}$), intent(inout) :: L + real(wp), dimension(${BOUND('num_dims')}$), intent(in) :: dvel_ds + real(wp), intent(in) :: rho, c, dpres_ds L(1) = f_base_L1(lambda, rho, c, dpres_ds, dvel_ds) L(2:eqn_idx%adv%end) = 0._wp @@ -216,24 +158,14 @@ contains $:GPU_ROUTINE(function_name='s_compute_nonreflecting_subsonic_outflow_L', parallelism='[seq]', cray_inline=True) - real(wp), dimension(3), intent(in) :: lambda - #:if USING_AMD - real(wp), dimension(${AMD_SYS_SIZE_MAX}$), intent(inout) :: L - #:else - real(wp), dimension(sys_size), intent(inout) :: L - #:endif - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3), intent(in) :: mf, dalpha_rho_ds - real(wp), dimension(3), intent(in) :: dvel_ds - real(wp), dimension(3), intent(in) :: dadv_ds - #:else - real(wp), dimension(num_fluids), intent(in) :: mf, dalpha_rho_ds - real(wp), dimension(num_dims), intent(in) :: dvel_ds - real(wp), dimension(num_fluids), intent(in) :: dadv_ds - #:endif - real(wp), dimension(${NUM_SPECIES}$), intent(in) :: dYs_ds - real(wp), intent(in) :: rho, c - real(wp), intent(in) :: dpres_ds + real(wp), dimension(3), intent(in) :: lambda + real(wp), dimension(${BOUND('sys_size')}$), intent(inout) :: L + real(wp), dimension(${BOUND('num_fluids')}$), intent(in) :: mf, dalpha_rho_ds + real(wp), dimension(${BOUND('num_dims')}$), intent(in) :: dvel_ds + real(wp), dimension(${BOUND('num_fluids')}$), intent(in) :: dadv_ds + real(wp), dimension(${NUM_SPECIES}$), intent(in) :: dYs_ds + real(wp), intent(in) :: rho, c + real(wp), intent(in) :: dpres_ds L(1) = f_base_L1(lambda, rho, c, dpres_ds, dvel_ds) call s_fill_density_L(L, 1._wp, lambda(2), c, mf, dalpha_rho_ds, dpres_ds) @@ -249,23 +181,13 @@ contains $:GPU_ROUTINE(function_name='s_compute_force_free_subsonic_outflow_L', parallelism='[seq]', cray_inline=True) - real(wp), dimension(3), intent(in) :: lambda - #:if USING_AMD - real(wp), dimension(${AMD_SYS_SIZE_MAX}$), intent(inout) :: L - #:else - real(wp), dimension(sys_size), intent(inout) :: L - #:endif - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3), intent(in) :: mf, dalpha_rho_ds - real(wp), dimension(3), intent(in) :: dvel_ds - real(wp), dimension(3), intent(in) :: dadv_ds - #:else - real(wp), dimension(num_fluids), intent(in) :: mf, dalpha_rho_ds - real(wp), dimension(num_dims), intent(in) :: dvel_ds - real(wp), dimension(num_fluids), intent(in) :: dadv_ds - #:endif - real(wp), intent(in) :: rho, c - real(wp), intent(in) :: dpres_ds + real(wp), dimension(3), intent(in) :: lambda + real(wp), dimension(${BOUND('sys_size')}$), intent(inout) :: L + real(wp), dimension(${BOUND('num_fluids')}$), intent(in) :: mf, dalpha_rho_ds + real(wp), dimension(${BOUND('num_dims')}$), intent(in) :: dvel_ds + real(wp), dimension(${BOUND('num_fluids')}$), intent(in) :: dadv_ds + real(wp), intent(in) :: rho, c + real(wp), intent(in) :: dpres_ds L(1) = f_base_L1(lambda, rho, c, dpres_ds, dvel_ds) call s_fill_density_L(L, 1._wp, lambda(2), c, mf, dalpha_rho_ds, dpres_ds) @@ -280,23 +202,13 @@ contains $:GPU_ROUTINE(function_name='s_compute_constant_pressure_subsonic_outflow_L', parallelism='[seq]', cray_inline=True) - real(wp), dimension(3), intent(in) :: lambda - #:if USING_AMD - real(wp), dimension(${AMD_SYS_SIZE_MAX}$), intent(inout) :: L - #:else - real(wp), dimension(sys_size), intent(inout) :: L - #:endif - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3), intent(in) :: mf, dalpha_rho_ds - real(wp), dimension(3), intent(in) :: dvel_ds - real(wp), dimension(3), intent(in) :: dadv_ds - #:else - real(wp), dimension(num_fluids), intent(in) :: mf, dalpha_rho_ds - real(wp), dimension(num_dims), intent(in) :: dvel_ds - real(wp), dimension(num_fluids), intent(in) :: dadv_ds - #:endif - real(wp), intent(in) :: rho, c - real(wp), intent(in) :: dpres_ds + real(wp), dimension(3), intent(in) :: lambda + real(wp), dimension(${BOUND('sys_size')}$), intent(inout) :: L + real(wp), dimension(${BOUND('num_fluids')}$), intent(in) :: mf, dalpha_rho_ds + real(wp), dimension(${BOUND('num_dims')}$), intent(in) :: dvel_ds + real(wp), dimension(${BOUND('num_fluids')}$), intent(in) :: dadv_ds + real(wp), intent(in) :: rho, c + real(wp), intent(in) :: dpres_ds L(1) = f_base_L1(lambda, rho, c, dpres_ds, dvel_ds) call s_fill_density_L(L, 1._wp, lambda(2), c, mf, dalpha_rho_ds, dpres_ds) @@ -310,11 +222,7 @@ contains subroutine s_compute_supersonic_inflow_L(L) $:GPU_ROUTINE(function_name='s_compute_supersonic_inflow_L', parallelism='[seq]', cray_inline=True) - #:if USING_AMD - real(wp), dimension(${AMD_SYS_SIZE_MAX}$), intent(inout) :: L - #:else - real(wp), dimension(sys_size), intent(inout) :: L - #:endif + real(wp), dimension(${BOUND('sys_size')}$), intent(inout) :: L L(1:eqn_idx%adv%end) = 0._wp if (chemistry) L(eqn_idx%species%beg:eqn_idx%species%end) = 0._wp @@ -325,24 +233,14 @@ contains $:GPU_ROUTINE(function_name='s_compute_supersonic_outflow_L', parallelism='[seq]', cray_inline=True) - real(wp), dimension(3), intent(in) :: lambda - #:if USING_AMD - real(wp), dimension(${AMD_SYS_SIZE_MAX}$), intent(inout) :: L - #:else - real(wp), dimension(sys_size), intent(inout) :: L - #:endif - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3), intent(in) :: mf, dalpha_rho_ds - real(wp), dimension(3), intent(in) :: dvel_ds - real(wp), dimension(3), intent(in) :: dadv_ds - #:else - real(wp), dimension(num_fluids), intent(in) :: mf, dalpha_rho_ds - real(wp), dimension(num_dims), intent(in) :: dvel_ds - real(wp), dimension(num_fluids), intent(in) :: dadv_ds - #:endif - real(wp), dimension(${NUM_SPECIES}$), intent(in) :: dYs_ds - real(wp), intent(in) :: rho, c - real(wp), intent(in) :: dpres_ds + real(wp), dimension(3), intent(in) :: lambda + real(wp), dimension(${BOUND('sys_size')}$), intent(inout) :: L + real(wp), dimension(${BOUND('num_fluids')}$), intent(in) :: mf, dalpha_rho_ds + real(wp), dimension(${BOUND('num_dims')}$), intent(in) :: dvel_ds + real(wp), dimension(${BOUND('num_fluids')}$), intent(in) :: dadv_ds + real(wp), dimension(${NUM_SPECIES}$), intent(in) :: dYs_ds + real(wp), intent(in) :: rho, c + real(wp), intent(in) :: dpres_ds L(1) = f_base_L1(lambda, rho, c, dpres_ds, dvel_ds) call s_fill_density_L(L, 1._wp, lambda(2), c, mf, dalpha_rho_ds, dpres_ds) diff --git a/src/simulation/m_data_output.fpp b/src/simulation/m_data_output.fpp index 7902ff5abc..ceb06ffe49 100644 --- a/src/simulation/m_data_output.fpp +++ b/src/simulation/m_data_output.fpp @@ -159,35 +159,29 @@ contains impure subroutine s_write_run_time_information(q_prim_vf, t_step) type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf - integer, intent(in) :: t_step - real(wp) :: rho !< Cell-avg. density - - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3) :: alpha, alpha_rho !< Cell-avg. volume fraction, partial density - real(wp), dimension(3) :: vel !< Cell-avg. velocity - #:else - real(wp), dimension(num_fluids) :: alpha, alpha_rho !< Cell-avg. volume fraction, partial density - real(wp), dimension(num_vels) :: vel !< Cell-avg. velocity - #:endif - real(wp) :: vel_sum !< Cell-avg. velocity sum - real(wp) :: pres !< Cell-avg. pressure - real(wp) :: gamma !< Cell-avg. sp. heat ratio - real(wp) :: pi_inf !< Cell-avg. liquid stiffness function - real(wp) :: qv !< Cell-avg. internal energy reference value - real(wp) :: c !< Cell-avg. sound speed - real(wp), dimension(2) :: Re !< Cell-avg. Reynolds numbers - integer :: j, k, l - real(wp) :: icfl_max_loc, icfl_max_glb !< ICFL stability extrema on local and global grids - real(wp) :: vcfl_max_loc, vcfl_max_glb !< VCFL stability extrema on local and global grids - real(wp) :: ccfl_max_loc, ccfl_max_glb !< CCFL stability extrema on local and global grids - real(wp) :: tcfl_max_loc, tcfl_max_glb !< TCFL stability extrema on local and global grids - real(wp) :: Rc_min_loc, Rc_min_glb !< Rc stability extrema on local and global grids - real(wp) :: icfl, vcfl, ccfl, tcfl, Rc - real(wp) :: mu_frac, mu_frac_max_loc, mu_frac_max_glb !< Compression as a fraction of the EOS limit - integer :: fl !< Fluid loop iterator - logical :: include_cell !< Cell is fluid, not ghost/inside an IB - real(wp), dimension(4) :: stab_max_loc, stab_max_glb !< Max-reduced criteria (ICFL, VCFL, CCFL, TCFL), packed - real(wp), dimension(1) :: stab_min_loc, stab_min_glb !< Min-reduced criteria (Rc), packed + integer, intent(in) :: t_step + real(wp) :: rho !< Cell-avg. density + real(wp), dimension(${BOUND('num_fluids')}$) :: alpha, alpha_rho !< Cell-avg. volume fraction, partial density + real(wp), dimension(${BOUND('num_vels')}$) :: vel !< Cell-avg. velocity + real(wp) :: vel_sum !< Cell-avg. velocity sum + real(wp) :: pres !< Cell-avg. pressure + real(wp) :: gamma !< Cell-avg. sp. heat ratio + real(wp) :: pi_inf !< Cell-avg. liquid stiffness function + real(wp) :: qv !< Cell-avg. internal energy reference value + real(wp) :: c !< Cell-avg. sound speed + real(wp), dimension(2) :: Re !< Cell-avg. Reynolds numbers + integer :: j, k, l + real(wp) :: icfl_max_loc, icfl_max_glb !< ICFL stability extrema on local and global grids + real(wp) :: vcfl_max_loc, vcfl_max_glb !< VCFL stability extrema on local and global grids + real(wp) :: ccfl_max_loc, ccfl_max_glb !< CCFL stability extrema on local and global grids + real(wp) :: tcfl_max_loc, tcfl_max_glb !< TCFL stability extrema on local and global grids + real(wp) :: Rc_min_loc, Rc_min_glb !< Rc stability extrema on local and global grids + real(wp) :: icfl, vcfl, ccfl, tcfl, Rc + real(wp) :: mu_frac, mu_frac_max_loc, mu_frac_max_glb !< Compression as a fraction of the EOS limit + integer :: fl !< Fluid loop iterator + logical :: include_cell !< Cell is fluid, not ghost/inside an IB + real(wp), dimension(4) :: stab_max_loc, stab_max_glb !< Max-reduced criteria (ICFL, VCFL, CCFL, TCFL), packed + real(wp), dimension(1) :: stab_min_loc, stab_min_glb !< Min-reduced criteria (Rc), packed icfl_max_loc = 0._wp vcfl_max_loc = 0._wp diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp index d31c4addb7..7b3d0beeea 100644 --- a/src/simulation/m_ibm.fpp +++ b/src/simulation/m_ibm.fpp @@ -244,32 +244,23 @@ contains type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf !< Primitive Variables type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf !< Primitive Variables real(stp), dimension(idwbuff(1)%beg:,idwbuff(2)%beg:,idwbuff(3)%beg:,1:,1:), intent(inout) :: pb_in, mv_in - integer :: i, j, k, l, q, r !< Iterator variables - integer :: jj, kk, ll !< Neighbor iterators - integer :: rad, rad_z, num_nbrs !< Neighbor stencil radius and population - integer :: patch_id, patch_id_temp !< Patch ID of ghost point - real(wp) :: rho, gamma, pi_inf, dyn_pres !< Mixture variables - real(wp) :: vel_sum_g, E_ghost !< Ghost-point velocity magnitude and energy + integer :: i, j, k, l, q, r !< Iterator variables + integer :: jj, kk, ll !< Neighbor iterators + integer :: rad, rad_z, num_nbrs !< Neighbor stencil radius and population + integer :: patch_id, patch_id_temp !< Patch ID of ghost point + real(wp) :: rho, gamma, pi_inf, dyn_pres !< Mixture variables + real(wp) :: vel_sum_g, E_ghost !< Ghost-point velocity magnitude and energy real(wp), dimension(2) :: Re_K real(wp) :: G_K real(wp) :: qv_K real(wp) :: pres_IP, pres_GP real(wp), dimension(3) :: vel_IP real(wp) :: c_IP - - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3) :: Gs - real(wp), dimension(3) :: alpha_rho_IP, alpha_IP, alpha_rho_GP - real(wp), dimension(3) :: r_IP, v_IP, pb_IP, mv_IP - real(wp), dimension(18) :: nmom_IP - real(wp), dimension(12) :: presb_IP, massv_IP - #:else - real(wp), dimension(num_fluids) :: Gs - real(wp), dimension(num_fluids) :: alpha_rho_IP, alpha_IP, alpha_rho_GP - real(wp), dimension(nb) :: r_IP, v_IP, pb_IP, mv_IP - real(wp), dimension(nb*nmom) :: nmom_IP - real(wp), dimension(nb*nnode) :: presb_IP, massv_IP - #:endif + real(wp), dimension(${BOUND('num_fluids')}$) :: Gs + real(wp), dimension(${BOUND('num_fluids')}$) :: alpha_rho_IP, alpha_IP, alpha_rho_GP + real(wp), dimension(${BOUND('nb')}$) :: r_IP, v_IP, pb_IP, mv_IP + real(wp), dimension(${BOUND('nb')}$*nmom) :: nmom_IP + real(wp), dimension(${BOUND('nb')}$*nnode) :: presb_IP, massv_IP real(wp), dimension(${NUM_SPECIES}$) :: Ys_IP real(wp) :: alpha_q, alpha_rho_q, e_q real(wp) :: T_IP, mw_IP, e_IP !< Image-point temperature, mixture MW, and mass-specific internal energy (chemistry) @@ -283,6 +274,7 @@ contains type(ghost_point) :: gp ! set the Moving IBM interior conservative variables + $:GPU_PARALLEL_LOOP(private='[i, j, k, patch_id, rho, patch_id_temp]', collapse=3) do l = 0, p do k = 0, n @@ -906,18 +898,14 @@ contains real(wp), intent(inout) :: pres_IP real(wp), dimension(3), intent(inout) :: vel_IP real(wp), intent(inout) :: c_IP - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3), intent(inout) :: alpha_IP, alpha_rho_IP - #:else - real(wp), dimension(num_fluids), intent(inout) :: alpha_IP, alpha_rho_IP - #:endif + real(wp), dimension(${BOUND('num_fluids')}$), intent(inout) :: alpha_IP, alpha_rho_IP real(wp), optional, dimension(:), intent(inout) :: r_IP, v_IP, pb_IP, mv_IP real(wp), optional, dimension(:), intent(inout) :: nmom_IP real(wp), optional, dimension(:), intent(inout) :: presb_IP, massv_IP - real(wp), optional, dimension(:), intent(inout) :: Ys_IP !< Interpolated species mass fractions (chemistry) - integer :: i, j, k, l, q !< Iterator variables - integer :: i1, i2, j1, j2, k1, k2 !< Iterator variables - real(wp) :: coeff + real(wp), optional, dimension(:), intent(inout) :: Ys_IP !< Interpolated species mass fractions (chemistry) + integer :: i, j, k, l, q !< Iterator variables + integer :: i1, i2, j1, j2, k1, k2 !< Iterator variables + real(wp) :: coeff i1 = gp%ip_grid(1); i2 = i1 + 1 j1 = gp%ip_grid(2); j2 = j1 + 1 @@ -1172,15 +1160,10 @@ contains integer :: i, j, k, l, encoded_ib_idx, xp, yp, zp, ib_idx, ib_idx_temp, fluid_idx real(wp), dimension(num_ibs, 3) :: forces, torques ! viscous stress tensor with temp vectors to hold divergence calculations - real(wp), dimension(1:3,1:3) :: viscous_stress - real(wp), dimension(1:3) :: local_force_contribution, radial_vector, local_torque_contribution - real(wp) :: cell_volume, dynamic_viscosity - - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3) :: dynamic_viscosities - #:else - real(wp), dimension(num_fluids) :: dynamic_viscosities - #:endif + real(wp), dimension(1:3,1:3) :: viscous_stress + real(wp), dimension(1:3) :: local_force_contribution, radial_vector, local_torque_contribution + real(wp) :: cell_volume, dynamic_viscosity + real(wp), dimension(${BOUND('num_fluids')}$) :: dynamic_viscosities call nvtxStartRange("COMPUTE-IB-FORCES") diff --git a/src/simulation/m_pressure_relaxation.fpp b/src/simulation/m_pressure_relaxation.fpp index 6c6dfe6ddb..d286ed5f3d 100644 --- a/src/simulation/m_pressure_relaxation.fpp +++ b/src/simulation/m_pressure_relaxation.fpp @@ -31,18 +31,14 @@ contains type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf integer :: i, j, k, l - - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3) :: alpha_rho, alpha - #:else - real(wp), dimension(num_fluids) :: alpha_rho, alpha - #:endif - real(wp) :: rho, gamma, pi_inf, qv_mix - integer :: hit_cap, hit_cap_sum, unusable, unusable_sum - real(wp) :: resid, resid_max + real(wp), dimension(${BOUND('num_fluids')}$) :: alpha_rho, alpha + real(wp) :: rho, gamma, pi_inf, qv_mix + integer :: hit_cap, hit_cap_sum, unusable, unusable_sum + real(wp) :: resid, resid_max ! Formed here, not one call deeper: CCE OpenACC accepts a num_fluids-sized array passed to a device routine from a ! parallel-loop body, and rejects the same call from inside another acc routine seq. + hit_cap_sum = 0 unusable_sum = 0 resid_max = 0._wp @@ -153,16 +149,12 @@ contains $:GPU_ROUTINE(parallelism='[seq]') type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf - integer, intent(in) :: j, k, l - integer, intent(out) :: hit_cap, unusable - real(wp), intent(out) :: resid - real(wp) :: pres_relax, f_pres, df_pres - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3) :: pres_K_init, rho_K_init, rho_K_s - #:else - real(wp), dimension(num_fluids) :: pres_K_init, rho_K_init, rho_K_s - #:endif - real(wp) :: gamma_K, pi_inf_K, dpi_K, dgamma_K, c2_K, alpha_i, alpha_rho_i, rho_i, p_i, rho_s_i + integer, intent(in) :: j, k, l + integer, intent(out) :: hit_cap, unusable + real(wp), intent(out) :: resid + real(wp) :: pres_relax, f_pres, df_pres + real(wp), dimension(${BOUND('num_fluids')}$) :: pres_K_init, rho_K_init, rho_K_s + real(wp) :: gamma_K, pi_inf_K, dpi_K, dgamma_K, c2_K, alpha_i, alpha_rho_i, rho_i, p_i, rho_s_i integer, parameter :: MAX_ITER = 50 ! Pressure relaxation convergence tolerance real(wp), parameter :: TOLERANCE = 1.e-10_wp diff --git a/src/simulation/m_qbmm.fpp b/src/simulation/m_qbmm.fpp index cde6830a7b..1c801838bf 100644 --- a/src/simulation/m_qbmm.fpp +++ b/src/simulation/m_qbmm.fpp @@ -588,14 +588,9 @@ contains $:GPU_ROUTINE(function_name='s_coeff_nonpoly',parallelism='[seq]', cray_inline=True) - real(wp), intent(in) :: pres, rho, c - #:if USING_AMD - real(wp), dimension(32,0:2,0:2), intent(out) :: coeffs - #:else - real(wp), dimension(nterms,0:2,0:2), intent(out) :: coeffs - #:endif - - integer :: i1, i2 + real(wp), intent(in) :: pres, rho, c + real(wp), dimension(${BOUND('nterms')}$,0:2,0:2), intent(out) :: coeffs + integer :: i1, i2 coeffs(:,:,:) = 0._wp @@ -667,14 +662,9 @@ contains $:GPU_ROUTINE(function_name='s_coeff',parallelism='[seq]', cray_inline=True) - real(wp), intent(in) :: pres, rho, c - #:if USING_AMD - real(wp), dimension(32,0:2,0:2), intent(out) :: coeffs - #:else - real(wp), dimension(nterms,0:2,0:2), intent(out) :: coeffs - #:endif - - integer :: i1, i2 + real(wp), intent(in) :: pres, rho, c + real(wp), dimension(${BOUND('nterms')}$,0:2,0:2), intent(out) :: coeffs + integer :: i1, i2 coeffs(:,:,:) = 0._wp @@ -734,29 +724,19 @@ contains !> Perform moment inversion to recover quadrature weights and abscissas and evaluate bubble source terms subroutine s_mom_inv(q_cons_vf, q_prim_vf, momsp, moms3d, pb, rhs_pb, mv, rhs_mv, ix, iy, iz) - type(scalar_field), dimension(:), intent(inout) :: q_cons_vf, q_prim_vf - type(scalar_field), dimension(:), intent(inout) :: momsp - type(scalar_field), dimension(0:,0:,:), intent(inout) :: moms3d + type(scalar_field), dimension(:), intent(inout) :: q_cons_vf, q_prim_vf + type(scalar_field), dimension(:), intent(inout) :: momsp + type(scalar_field), dimension(0:,0:,:), intent(inout) :: moms3d real(stp), dimension(idwbuff(1)%beg:,idwbuff(2)%beg:,idwbuff(3)%beg:,1:,1:), intent(inout) :: pb - real(wp), dimension(idwbuff(1)%beg:,idwbuff(2)%beg:,idwbuff(3)%beg:,1:,1:), intent(inout) :: rhs_pb + real(wp), dimension(idwbuff(1)%beg:,idwbuff(2)%beg:,idwbuff(3)%beg:,1:,1:), intent(inout) :: rhs_pb real(stp), dimension(idwbuff(1)%beg:,idwbuff(2)%beg:,idwbuff(3)%beg:,1:,1:), intent(inout) :: mv - real(wp), dimension(idwbuff(1)%beg:,idwbuff(2)%beg:,idwbuff(3)%beg:,1:,1:), intent(inout) :: rhs_mv - type(int_bounds_info), intent(in) :: ix, iy, iz - - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(6) :: moms, msum - real(wp), dimension(4, 3) :: wght, abscX, abscY, wght_pb, wght_mv, wght_ht, ht - #:else - real(wp), dimension(nmom) :: moms, msum - real(wp), dimension(nnode, nb) :: wght, abscX, abscY, wght_pb, wght_mv, wght_ht, ht - #:endif - #:if USING_AMD - real(wp), dimension(32,0:2,0:2) :: coeff - #:else - real(wp), dimension(nterms,0:2,0:2) :: coeff - #:endif + real(wp), dimension(idwbuff(1)%beg:,idwbuff(2)%beg:,idwbuff(3)%beg:,1:,1:), intent(inout) :: rhs_mv + type(int_bounds_info), intent(in) :: ix, iy, iz + real(wp), dimension(nmom) :: moms, msum + real(wp), dimension(nnode, ${BOUND('nb')}$) :: wght, abscX, abscY, wght_pb, wght_mv, wght_ht, ht + real(wp), dimension(${BOUND('nterms')}$,0:2,0:2) :: coeff real(wp) :: pres, rho, nbub, c, alf, momsum, drdt, drdt2, chi_vw, x_vw, rho_mw, k_mw, grad_T - integer :: id1, id2, id3, i1, i2, j, q, r + integer :: id1, id2, id3, i1, i2, j, q, r is1_qbmm = ix; is2_qbmm = iy; is3_qbmm = iz $:GPU_UPDATE(device='[is1_qbmm, is2_qbmm, is3_qbmm]') @@ -929,13 +909,9 @@ contains subroutine s_coeff_selector(pres, rho, c, coeff, polytropic) $:GPU_ROUTINE(function_name='s_coeff_selector',parallelism='[seq]', cray_inline=True) - real(wp), intent(in) :: pres, rho, c - #:if USING_AMD - real(wp), dimension(32,0:2,0:2), intent(out) :: coeff - #:else - real(wp), dimension(nterms,0:2,0:2), intent(out) :: coeff - #:endif - logical, intent(in) :: polytropic + real(wp), intent(in) :: pres, rho, c + real(wp), dimension(${BOUND('nterms')}$,0:2,0:2), intent(out) :: coeff + logical, intent(in) :: polytropic if (polytropic) then call s_coeff(pres, rho, c, coeff) else @@ -1027,14 +1003,10 @@ contains function f_quad(abscX, abscY, wght_in, q, r, s) $:GPU_ROUTINE(parallelism='[seq]') - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(4, 3), intent(in) :: abscX, abscY, wght_in - #:else - real(wp), dimension(nnode, nb), intent(in) :: abscX, abscY, wght_in - #:endif - real(wp), intent(in) :: q, r, s - real(wp) :: f_quad_RV, f_quad - integer :: i, i1 + real(wp), dimension(nnode, ${BOUND('nb')}$), intent(in) :: abscX, abscY, wght_in + real(wp), intent(in) :: q, r, s + real(wp) :: f_quad_RV, f_quad + integer :: i, i1 f_quad = 0._wp $:GPU_LOOP(parallelism='[seq]') @@ -1053,14 +1025,10 @@ contains function f_quad2D(abscX, abscY, wght_in, pow) $:GPU_ROUTINE(parallelism='[seq]') - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(4), intent(in) :: abscX, abscY, wght_in - #:else - real(wp), dimension(nnode), intent(in) :: abscX, abscY, wght_in - #:endif - real(wp), dimension(3), intent(in) :: pow - real(wp) :: f_quad2D - integer :: i + real(wp), dimension(nnode), intent(in) :: abscX, abscY, wght_in + real(wp), dimension(3), intent(in) :: pow + real(wp) :: f_quad2D + integer :: i f_quad2D = 0._wp $:GPU_LOOP(parallelism='[seq]') diff --git a/src/simulation/m_riemann_solver_hll.fpp b/src/simulation/m_riemann_solver_hll.fpp index 5cddef70c3..5416bb90aa 100644 --- a/src/simulation/m_riemann_solver_hll.fpp +++ b/src/simulation/m_riemann_solver_hll.fpp @@ -36,19 +36,12 @@ contains ! Intercell fluxes type(scalar_field), dimension(sys_size), intent(inout) :: flux_vf, flux_src_vf, flux_gsrc_vf - real(wp) :: flux_tau_L, flux_tau_R - integer, intent(in) :: norm_dir - type(int_bounds_info), intent(in) :: ix, iy, iz - - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3) :: alpha_rho_L, alpha_rho_R - real(wp), dimension(3) :: vel_L, vel_R - real(wp), dimension(3) :: alpha_L, alpha_R - #:else - real(wp), dimension(num_fluids) :: alpha_rho_L, alpha_rho_R - real(wp), dimension(num_vels) :: vel_L, vel_R - real(wp), dimension(num_fluids) :: alpha_L, alpha_R - #:endif + real(wp) :: flux_tau_L, flux_tau_R + integer, intent(in) :: norm_dir + type(int_bounds_info), intent(in) :: ix, iy, iz + real(wp), dimension(${BOUND('num_fluids')}$) :: alpha_rho_L, alpha_rho_R + real(wp), dimension(${BOUND('num_vels')}$) :: vel_L, vel_R + real(wp), dimension(${BOUND('num_fluids')}$) :: alpha_L, alpha_R real(wp), dimension(${NUM_SPECIES}$) :: Ys_L, Ys_R, R_species, h_iL, h_iR real(wp), dimension(${NUM_SPECIES}$) :: Cp_iL, Cp_iR, Xs_L, Xs_R, Gamma_iL, Gamma_iR real(wp) :: rho_L, rho_R diff --git a/src/simulation/m_riemann_solver_hllc.fpp b/src/simulation/m_riemann_solver_hllc.fpp index a94badf27a..5b80e838da 100644 --- a/src/simulation/m_riemann_solver_hllc.fpp +++ b/src/simulation/m_riemann_solver_hllc.fpp @@ -41,60 +41,43 @@ contains type(scalar_field), dimension(sys_size), intent(inout) :: flux_vf, flux_src_vf, flux_gsrc_vf integer, intent(in) :: norm_dir type(int_bounds_info), intent(in) :: ix, iy, iz - - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3) :: alpha_rho_L, alpha_rho_R - real(wp), dimension(3) :: alpha_L, alpha_R - real(wp), dimension(3) :: alpha_lim_L, alpha_lim_R - real(wp), dimension(3) :: vel_L, vel_R - #:else - real(wp), dimension(num_fluids) :: alpha_rho_L, alpha_rho_R - real(wp), dimension(num_fluids) :: alpha_L, alpha_R - !> Post-limiter volume fractions (alpha_L/R retain the pre-limiter loads used downstream) - real(wp), dimension(num_fluids) :: alpha_lim_L, alpha_lim_R - real(wp), dimension(num_dims) :: vel_L, vel_R - #:endif - - real(wp) :: rho_L, rho_R - real(wp) :: pres_L, pres_R - real(wp) :: E_L, E_R - real(wp) :: H_L, H_R + real(wp), dimension(${BOUND('num_fluids')}$) :: alpha_rho_L, alpha_rho_R + real(wp), dimension(${BOUND('num_fluids')}$) :: alpha_L, alpha_R + !> Post-limiter volume fractions (alpha_L/R retain the pre-limiter loads used downstream) + real(wp), dimension(${BOUND('num_fluids')}$) :: alpha_lim_L, alpha_lim_R + real(wp), dimension(${BOUND('num_dims')}$) :: vel_L, vel_R + real(wp) :: rho_L, rho_R + real(wp) :: pres_L, pres_R + real(wp) :: E_L, E_R + real(wp) :: H_L, H_R real(wp), dimension(${NUM_SPECIES}$) :: Ys_L, Ys_R, Xs_L, Xs_R, Gamma_iL, Gamma_iR, Cp_iL, Cp_iR, R_species, h_iL, h_iR - real(wp) :: c_sum_Yi_Phi - real(wp) :: T_L, T_R - real(wp) :: MW_L, MW_R - real(wp) :: R_gas_L, R_gas_R - real(wp) :: Cp_L, Cp_R - real(wp) :: Cv_L, Cv_R - real(wp) :: Gamm_L, Gamm_R - real(wp) :: Y_L, Y_R - real(wp) :: gamma_L, gamma_R - real(wp) :: pi_inf_L, pi_inf_R - real(wp) :: qv_L, qv_R - real(wp) :: c_L, c_R - real(wp), dimension(2) :: Re_L, Re_R - real(wp) :: rho_avg - real(wp) :: H_avg - real(wp) :: gamma_avg - real(wp) :: qv_avg - real(wp) :: c_avg - real(wp) :: s_L, s_R, s_M, s_P, s_S - real(wp) :: xi_L, xi_R !< Left and right wave speeds functions - real(wp) :: xi_L_m1, xi_R_m1 !< xi_L/R - 1, computed without cancellation - real(wp) :: xi_M, xi_P - real(wp) :: xi_MP, xi_PP - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3) :: R0_L, R0_R - real(wp), dimension(3) :: V0_L, V0_R - real(wp), dimension(3) :: P0_L, P0_R - real(wp), dimension(3) :: pbw_L, pbw_R - #:else - real(wp), dimension(nb) :: R0_L, R0_R - real(wp), dimension(nb) :: V0_L, V0_R - real(wp), dimension(nb) :: P0_L, P0_R - real(wp), dimension(nb) :: pbw_L, pbw_R - #:endif - + real(wp) :: c_sum_Yi_Phi + real(wp) :: T_L, T_R + real(wp) :: MW_L, MW_R + real(wp) :: R_gas_L, R_gas_R + real(wp) :: Cp_L, Cp_R + real(wp) :: Cv_L, Cv_R + real(wp) :: Gamm_L, Gamm_R + real(wp) :: Y_L, Y_R + real(wp) :: gamma_L, gamma_R + real(wp) :: pi_inf_L, pi_inf_R + real(wp) :: qv_L, qv_R + real(wp) :: c_L, c_R + real(wp), dimension(2) :: Re_L, Re_R + real(wp) :: rho_avg + real(wp) :: H_avg + real(wp) :: gamma_avg + real(wp) :: qv_avg + real(wp) :: c_avg + real(wp) :: s_L, s_R, s_M, s_P, s_S + real(wp) :: xi_L, xi_R !< Left and right wave speeds functions + real(wp) :: xi_L_m1, xi_R_m1 !< xi_L/R - 1, computed without cancellation + real(wp) :: xi_M, xi_P + real(wp) :: xi_MP, xi_PP + real(wp), dimension(${BOUND('nb')}$) :: R0_L, R0_R + real(wp), dimension(${BOUND('nb')}$) :: V0_L, V0_R + real(wp), dimension(${BOUND('nb')}$) :: P0_L, P0_R + real(wp), dimension(${BOUND('nb')}$) :: pbw_L, pbw_R real(wp) :: alpha_L_sum, alpha_R_sum, nbub_L, nbub_R real(wp) :: ptilde_L, ptilde_R real(wp) :: PbwR3Lbar, PbwR3Rbar @@ -113,43 +96,34 @@ contains integer :: Re_size_loc1, Re_size_loc2 !< host copy of Re_size; amdflang reads the declare-target original stale cross-TU ! HLLC star-state helpers - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(${AMD_SYS_SIZE_MAX}$) :: U_L, U_R - real(wp), dimension(${AMD_SYS_SIZE_MAX}$) :: F_L, F_R, F_star_L, F_star_R, F_HLLC - #:else - real(wp), dimension(sys_size) :: U_L, U_R - real(wp), dimension(sys_size) :: F_L, F_R, F_star_L, F_star_R, F_HLLC - #:endif - real(wp) :: u_n_HLLC, u_t_HLLC, u_t2_HLLC - real(wp) :: pres_tot_L, pres_tot_R - real(wp) :: u_n_L, u_n_R, u_t_L, u_t_R - real(wp) :: u_t2_L, u_t2_R - real(wp) :: tau_nn_L, tau_nn_R, tau_nt_L, tau_nt_R, tau_tt_L, tau_tt_R - real(wp) :: tau_nt2_L, tau_nt2_R, tau_t2t2_L, tau_t2t2_R, tau_t1t2_L, tau_t1t2_R - real(wp) :: tau_qq_L, tau_qq_R - real(wp) :: p_face, tau_qq_face - real(wp) :: A_L, A_R, denom_A - real(wp) :: u_t_star, tau_nt_star - real(wp) :: u_t2_star, tau_nt2_star - real(wp) :: pres_tot_star - integer :: idx_phys + real(wp), dimension(${BOUND('sys_size')}$) :: U_L, U_R + real(wp), dimension(${BOUND('sys_size')}$) :: F_L, F_R, F_star_L, F_star_R, F_HLLC + real(wp) :: u_n_HLLC, u_t_HLLC, u_t2_HLLC + real(wp) :: pres_tot_L, pres_tot_R + real(wp) :: u_n_L, u_n_R, u_t_L, u_t_R + real(wp) :: u_t2_L, u_t2_R + real(wp) :: tau_nn_L, tau_nn_R, tau_nt_L, tau_nt_R, tau_tt_L, tau_tt_R + real(wp) :: tau_nt2_L, tau_nt2_R, tau_t2t2_L, tau_t2t2_R, tau_t1t2_L, tau_t1t2_R + real(wp) :: tau_qq_L, tau_qq_R + real(wp) :: p_face, tau_qq_face + real(wp) :: A_L, A_R, denom_A + real(wp) :: u_t_star, tau_nt_star + real(wp) :: u_t2_star, tau_nt2_star + real(wp) :: pres_tot_star + integer :: idx_phys ! ADC (HLL -> HLLC) - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(${AMD_SYS_SIZE_MAX}$) :: F_HLL - #:else - real(wp), dimension(sys_size) :: F_HLL - #:endif - real(wp) :: u_n_HLL_trace, u_t_HLL_trace - real(wp) :: u_t2_HLL_trace - real(wp) :: p_face_HLL, tau_qq_face_HLL, tau_nn_HLL - real(wp) :: phi - real(wp) :: Sigma_L, Sigma_R, dSigma, Sigma_ref - real(wp) :: a_L_ref, a_R_ref, a_ref - real(wp) :: du_t, dtau_nt - real(wp) :: du_t2, dtau_nt2 - real(wp) :: sensor_ptot, sensor_vt, sensor_tnt, sensor_combined - real(wp), parameter :: ADC_power = 1.0_wp + real(wp), dimension(${BOUND('sys_size')}$) :: F_HLL + real(wp) :: u_n_HLL_trace, u_t_HLL_trace + real(wp) :: u_t2_HLL_trace + real(wp) :: p_face_HLL, tau_qq_face_HLL, tau_nn_HLL + real(wp) :: phi + real(wp) :: Sigma_L, Sigma_R, dSigma, Sigma_ref + real(wp) :: a_L_ref, a_R_ref, a_ref + real(wp) :: du_t, dtau_nt + real(wp) :: du_t2, dtau_nt2 + real(wp) :: sensor_ptot, sensor_vt, sensor_tnt, sensor_combined + real(wp), parameter :: ADC_power = 1.0_wp ! Populating the buffers of the left and right Riemann problem states variables, based on the choice of boundary conditions diff --git a/src/simulation/m_riemann_solver_hlld.fpp b/src/simulation/m_riemann_solver_hlld.fpp index bce95254d8..de284bee54 100644 --- a/src/simulation/m_riemann_solver_hlld.fpp +++ b/src/simulation/m_riemann_solver_hlld.fpp @@ -34,17 +34,13 @@ contains ! Local variables: - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3) :: alpha_L, alpha_R, alpha_rho_L, alpha_rho_R - #:else - real(wp), dimension(num_fluids) :: alpha_L, alpha_R, alpha_rho_L, alpha_rho_R - #:endif - type(riemann_states_vec3) :: vel - type(riemann_states) :: rho, pres, E, H_no_mag - type(riemann_states) :: gamma, pi_inf, qv - type(riemann_states) :: vel_rms - type(riemann_states_vec3) :: B - type(riemann_states) :: c, c_fast, pres_mag + real(wp), dimension(${BOUND('num_fluids')}$) :: alpha_L, alpha_R, alpha_rho_L, alpha_rho_R + type(riemann_states_vec3) :: vel + type(riemann_states) :: rho, pres, E, H_no_mag + type(riemann_states) :: gamma, pi_inf, qv + type(riemann_states) :: vel_rms + type(riemann_states_vec3) :: B + type(riemann_states) :: c, c_fast, pres_mag ! HLLD speeds and intermediate state variables: real(wp) :: s_L, s_R, s_M, s_starL, s_starR diff --git a/src/simulation/m_riemann_solver_hypo_hlld.fpp b/src/simulation/m_riemann_solver_hypo_hlld.fpp index 2fb798e6a8..dbc600fe28 100644 --- a/src/simulation/m_riemann_solver_hypo_hlld.fpp +++ b/src/simulation/m_riemann_solver_hypo_hlld.fpp @@ -83,10 +83,8 @@ contains #:if MFC_CASE_OPTIMIZATION real(wp), dimension(${max(num_fluids, 2)}$) :: alpha_L, alpha_R, alpha_rho_L, alpha_rho_R - #:elif USING_AMD - real(wp), dimension(3) :: alpha_L, alpha_R, alpha_rho_L, alpha_rho_R #:else - real(wp), dimension(num_fluids) :: alpha_L, alpha_R, alpha_rho_L, alpha_rho_R + real(wp), dimension(${BOUND('num_fluids')}$) :: alpha_L, alpha_R, alpha_rho_L, alpha_rho_R #:endif type(riemann_states_vec3) :: vel type(riemann_states) :: rho, pres, E @@ -109,62 +107,52 @@ contains ! HLLD Hypo variables - real(wp) :: G_eff, G_eff_tol, C_NC, sqrtC_NC - real(wp) :: A_L, A_R, denomA, fac_L, fac_R - real(wp) :: u_n_L, u_t_L, u_n_R, u_t_R - real(wp) :: u_t2_L, u_t2_R - real(wp) :: tau_nn_L, tau_nt_L, tau_tt_L, tau_nn_R, tau_nt_R, tau_tt_R - real(wp) :: tau_nt2_L, tau_nt2_R, tau_t2t2_L, tau_t2t2_R, tau_t1t2_L, tau_t1t2_R - real(wp) :: tau_qq_L, tau_qq_R - real(wp) :: G_L, G_R - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(6) :: tau_e_L, tau_e_R - #:else - real(wp), dimension(eqn_idx%stress%end - eqn_idx%stress%beg + 1) :: tau_e_L, tau_e_R - #:endif - - real(wp) :: alpha1_L_star, alpha1_R_star, alpha2_L_star, alpha2_R_star - real(wp) :: u_t_star, tau_nt_star - real(wp) :: u_t2_star, tau_nt2_star - real(wp) :: tau_nn_L_star, tau_nn_R_star, tau_tt_L_star, tau_tt_R_star - real(wp) :: tau_tt_L_starstar, tau_tt_R_starstar - real(wp) :: tau_t2t2_L_star, tau_t2t2_R_star - real(wp) :: tau_t2t2_L_starstar, tau_t2t2_R_starstar - real(wp) :: tau_t1t2_L_star, tau_t1t2_R_star - real(wp) :: tau_t1t2_L_starstar, tau_t1t2_R_starstar - real(wp) :: tau_qq_L_star, tau_qq_R_star - real(wp) :: pTot_star - real(wp) :: E_L_star, E_R_star - real(wp) :: E_L_starstar, E_R_starstar - real(wp) :: p_face, tau_qq_face - real(wp) :: u_n_face, u_t_face - real(wp) :: G_hat - real(wp) :: rho_hat - real(wp) :: tau_nn_hat, tau_nt_hat, tau_tt_hat, tau_qq_hat - real(wp) :: tau_nt2_hat, tau_t2t2_hat, tau_t1t2_hat + real(wp) :: G_eff, G_eff_tol, C_NC, sqrtC_NC + real(wp) :: A_L, A_R, denomA, fac_L, fac_R + real(wp) :: u_n_L, u_t_L, u_n_R, u_t_R + real(wp) :: u_t2_L, u_t2_R + real(wp) :: tau_nn_L, tau_nt_L, tau_tt_L, tau_nn_R, tau_nt_R, tau_tt_R + real(wp) :: tau_nt2_L, tau_nt2_R, tau_t2t2_L, tau_t2t2_R, tau_t1t2_L, tau_t1t2_R + real(wp) :: tau_qq_L, tau_qq_R + real(wp) :: G_L, G_R + real(wp), dimension(${BOUND('n_stress')}$) :: tau_e_L, tau_e_R + real(wp) :: alpha1_L_star, alpha1_R_star, alpha2_L_star, alpha2_R_star + real(wp) :: u_t_star, tau_nt_star + real(wp) :: u_t2_star, tau_nt2_star + real(wp) :: tau_nn_L_star, tau_nn_R_star, tau_tt_L_star, tau_tt_R_star + real(wp) :: tau_tt_L_starstar, tau_tt_R_starstar + real(wp) :: tau_t2t2_L_star, tau_t2t2_R_star + real(wp) :: tau_t2t2_L_starstar, tau_t2t2_R_starstar + real(wp) :: tau_t1t2_L_star, tau_t1t2_R_star + real(wp) :: tau_t1t2_L_starstar, tau_t1t2_R_starstar + real(wp) :: tau_qq_L_star, tau_qq_R_star + real(wp) :: pTot_star + real(wp) :: E_L_star, E_R_star + real(wp) :: E_L_starstar, E_R_starstar + real(wp) :: p_face, tau_qq_face + real(wp) :: u_n_face, u_t_face + real(wp) :: G_hat + real(wp) :: rho_hat + real(wp) :: tau_nn_hat, tau_nt_hat, tau_tt_hat, tau_qq_hat + real(wp) :: tau_nt2_hat, tau_t2t2_hat, tau_t1t2_hat ! alpha_hat/alpha_rho_hat: same max(num_fluids, 2) reason as alpha_* above #:if MFC_CASE_OPTIMIZATION - real(wp), dimension(${max(num_fluids, 2)}$) :: alpha_hat, alpha_rho_hat - real(wp), dimension(eqn_idx%stress%end - eqn_idx%stress%beg + 1) :: tau_e_hat - #:elif USING_AMD - real(wp), dimension(3) :: alpha_hat, alpha_rho_hat - real(wp), dimension(6) :: tau_e_hat + real(wp), dimension(${max(num_fluids, 2)}$) :: alpha_hat, alpha_rho_hat #:else - real(wp), dimension(num_fluids) :: alpha_hat, alpha_rho_hat - real(wp), dimension(eqn_idx%stress%end - eqn_idx%stress%beg + 1) :: tau_e_hat + real(wp), dimension(${BOUND('num_fluids')}$) :: alpha_hat, alpha_rho_hat #:endif - - real(wp) :: pres_hat, blkmod1_hat, blkmod2_hat, K_hat, alpha_hat_q, alpha_rho_hat_q - real(wp) :: C_hat_1, C_hat_2 - real(wp) :: Sigma_L, Sigma_R, dSigma, Sigma_ref - real(wp) :: a_L_ref, a_R_ref, a_ref - real(wp) :: du_t, dtau_nt, du_t2, dtau_nt2 - real(wp) :: sensor_ptot, sensor_vt, sensor_tnt, sensor_combined - real(wp) :: phi - real(wp), parameter :: ADC_power = 1.0_wp - real(wp) :: alpha_L_sum, alpha_R_sum - logical :: degenerate, shear_degenerate, fan_fallback, shear_cond - integer :: i, j, k, l, ipass, zone + real(wp), dimension(${BOUND('n_stress')}$) :: tau_e_hat + real(wp) :: pres_hat, blkmod1_hat, blkmod2_hat, K_hat, alpha_hat_q, alpha_rho_hat_q + real(wp) :: C_hat_1, C_hat_2 + real(wp) :: Sigma_L, Sigma_R, dSigma, Sigma_ref + real(wp) :: a_L_ref, a_R_ref, a_ref + real(wp) :: du_t, dtau_nt, du_t2, dtau_nt2 + real(wp) :: sensor_ptot, sensor_vt, sensor_tnt, sensor_combined + real(wp) :: phi + real(wp), parameter :: ADC_power = 1.0_wp + real(wp) :: alpha_L_sum, alpha_R_sum + logical :: degenerate, shear_degenerate, fan_fallback, shear_cond + integer :: i, j, k, l, ipass, zone call s_populate_riemann_states_variables_buffers(qL_prim_rsx_vf, dqL_prim_dx_vf, dqL_prim_dy_vf, dqL_prim_dz_vf, & & qR_prim_rsx_vf, dqR_prim_dx_vf, dqR_prim_dy_vf, dqR_prim_dz_vf, norm_dir, ix, iy, iz) diff --git a/src/simulation/m_riemann_solver_lf.fpp b/src/simulation/m_riemann_solver_lf.fpp index 6893a31683..80531fe60c 100644 --- a/src/simulation/m_riemann_solver_lf.fpp +++ b/src/simulation/m_riemann_solver_lf.fpp @@ -35,19 +35,11 @@ contains type(scalar_field), dimension(sys_size), intent(inout) :: flux_vf, flux_src_vf, flux_gsrc_vf integer, intent(in) :: norm_dir type(int_bounds_info), intent(in) :: ix, iy, iz - - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3) :: alpha_rho_L, alpha_rho_R - real(wp), dimension(3) :: vel_L, vel_R - real(wp), dimension(3) :: alpha_L, alpha_R - real(wp), dimension(3, 3) :: vel_grad_L, vel_grad_R !< Averaged velocity gradient tensor `d(vel_i)/d(coord_j)`. - #:else - real(wp), dimension(num_fluids) :: alpha_rho_L, alpha_rho_R - real(wp), dimension(num_vels) :: vel_L, vel_R - real(wp), dimension(num_fluids) :: alpha_L, alpha_R - !> Averaged velocity gradient tensor `d(vel_i)/d(coord_j)`. - real(wp), dimension(num_dims, num_dims) :: vel_grad_L, vel_grad_R - #:endif + real(wp), dimension(${BOUND('num_fluids')}$) :: alpha_rho_L, alpha_rho_R + real(wp), dimension(${BOUND('num_vels')}$) :: vel_L, vel_R + real(wp), dimension(${BOUND('num_fluids')}$) :: alpha_L, alpha_R + !> Averaged velocity gradient tensor `d(vel_i)/d(coord_j)`. + real(wp), dimension(${BOUND('num_dims')}$, ${BOUND('num_dims')}$) :: vel_grad_L, vel_grad_R real(wp), dimension(${NUM_SPECIES}$) :: Ys_L, Ys_R real(wp), dimension(${NUM_SPECIES}$) :: Cp_iL, Cp_iR, Xs_L, Xs_R, Gamma_iL, Gamma_iR real(wp) :: rho_L, rho_R diff --git a/src/simulation/m_riemann_state.fpp b/src/simulation/m_riemann_state.fpp index aeb34c8a9f..98c3f3bfcf 100644 --- a/src/simulation/m_riemann_state.fpp +++ b/src/simulation/m_riemann_state.fpp @@ -175,14 +175,10 @@ contains real(wp), intent(in) :: H_L, H_R !< Left and right total enthalpies real(wp), intent(in) :: gamma_L, gamma_R !< Left and right specific heat ratio functions real(wp), intent(in) :: qv_L, qv_R !< Left and right reference energies - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3), intent(in) :: vel_L, vel_R - #:else - real(wp), dimension(num_vels), intent(in) :: vel_L, vel_R - #:endif + real(wp), dimension(${BOUND('num_vels')}$), intent(in) :: vel_L, vel_R real(wp), intent(out) :: rho_avg, H_avg, gamma_avg, qv_avg - real(wp), intent(out) :: vel_avg_rms !< Squared magnitude of the averaged velocity, summed over all components - integer :: i + real(wp), intent(out) :: vel_avg_rms !< Squared magnitude of the averaged velocity, summed over all components + integer :: i vel_avg_rms = 0._wp @@ -779,25 +775,18 @@ contains ! Local variables - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3) :: avg_v_int !< Averaged interface velocity (\f$v_x, v_y, v_z\f$) (grid directions). - real(wp), dimension(3) :: avg_dvdx_int !< Averaged interface \f$\partial v_i/\partial x\f$ (grid dir 1). - real(wp), dimension(3) :: avg_dvdy_int !< Averaged interface \f$\partial v_i/\partial y\f$ (grid dir 2). - real(wp), dimension(3) :: avg_dvdz_int !< Averaged interface \f$\partial v_i/\partial z\f$ (grid dir 3). - real(wp), dimension(3) :: vel_src_int !< Interface velocity (\f$v_1,v_2,v_3\f$) (grid directions) for viscous work. - - !> Shear stress vector (\f$\sigma_{N1}, \sigma_{N2}, \sigma_{N3}\f$) on N-face (grid directions). - real(wp), dimension(3) :: stress_vector_shear - #:else - real(wp), dimension(num_dims) :: avg_v_int !< Averaged interface velocity (\f$v_x, v_y, v_z\f$) (grid directions). - real(wp), dimension(num_dims) :: avg_dvdx_int !< Averaged interface \f$\partial v_i/\partial x\f$ (grid dir 1). - real(wp), dimension(num_dims) :: avg_dvdy_int !< Averaged interface \f$\partial v_i/\partial y\f$ (grid dir 2). - real(wp), dimension(num_dims) :: avg_dvdz_int !< Averaged interface \f$\partial v_i/\partial z\f$ (grid dir 3). - !> Interface velocity (\f$v_1,v_2,v_3\f$) (grid directions) for viscous work. - real(wp), dimension(num_dims) :: vel_src_int - !> Shear stress vector (\f$\sigma_{N1}, \sigma_{N2}, \sigma_{N3}\f$) on N-face (grid directions). - real(wp), dimension(num_dims) :: stress_vector_shear - #:endif + !> Averaged interface velocity (\f$v_x, v_y, v_z\f$) (grid directions). + real(wp), dimension(${BOUND('num_dims')}$) :: avg_v_int + !> Averaged interface \f$\partial v_i/\partial x\f$ (grid dir 1). + real(wp), dimension(${BOUND('num_dims')}$) :: avg_dvdx_int + !> Averaged interface \f$\partial v_i/\partial y\f$ (grid dir 2). + real(wp), dimension(${BOUND('num_dims')}$) :: avg_dvdy_int + !> Averaged interface \f$\partial v_i/\partial z\f$ (grid dir 3). + real(wp), dimension(${BOUND('num_dims')}$) :: avg_dvdz_int + !> Interface velocity (\f$v_1,v_2,v_3\f$) (grid directions) for viscous work. + real(wp), dimension(${BOUND('num_dims')}$) :: vel_src_int + !> Shear stress vector (\f$\sigma_{N1}, \sigma_{N2}, \sigma_{N3}\f$) on N-face (grid directions). + real(wp), dimension(${BOUND('num_dims')}$) :: stress_vector_shear real(wp) :: stress_normal_bulk !< Normal bulk stress component \f$\sigma_{NN}\f$ on N-face. real(wp) :: Re_s, Re_b !< Effective interface shear and bulk Reynolds numbers. real(wp) :: r_eff !< Effective radius at interface for cylindrical terms. @@ -1000,30 +989,24 @@ contains ! Local variables - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3, 3) :: vel_grad_avg !< Averaged velocity gradient tensor `d(vel_i)/d(coord_j)`. - real(wp), dimension(3, 3) :: current_tau_shear !< Current shear stress tensor. - real(wp), dimension(3, 3) :: current_tau_bulk !< Current bulk stress tensor. - real(wp), dimension(3) :: vel_src_at_interface !< Interface velocities (u,v,w) for viscous work. - #:else - real(wp), dimension(num_dims, num_dims) :: vel_grad_avg !< Averaged velocity gradient tensor `d(vel_i)/d(coord_j)`. - real(wp), dimension(num_dims, num_dims) :: current_tau_shear !< Current shear stress tensor. - real(wp), dimension(num_dims, num_dims) :: current_tau_bulk !< Current bulk stress tensor. - real(wp), dimension(num_dims) :: vel_src_at_interface !< Interface velocities (u,v,w) for viscous work. - #:endif - integer, dimension(3) :: idx_right_phys !< Physical (j,k,l) indices for right state. - real(wp) :: Re_shear !< Interface shear Reynolds number. - real(wp) :: Re_bulk !< Interface bulk Reynolds number. - integer :: j_loop !< Physical x-index loop iterator. - integer :: k_loop !< Physical y-index loop iterator. - integer :: l_loop !< Physical z-index loop iterator. - integer :: i_dim !< Generic dimension/component iterator. - integer :: vel_comp_idx !< Velocity component iterator (1=u, 2=v, 3=w). - real(wp) :: divergence_v !< Velocity divergence at interface. - real(wp) :: gamma_dot, D_xx, D_yy, D_zz, D_xy, D_xz, D_yz - real(wp), dimension(2) :: Re_nn + !> Averaged velocity gradient tensor `d(vel_i)/d(coord_j)`. + real(wp), dimension(${BOUND('num_dims')}$, ${BOUND('num_dims')}$) :: vel_grad_avg + real(wp), dimension(${BOUND('num_dims')}$, ${BOUND('num_dims')}$) :: current_tau_shear !< Current shear stress tensor. + real(wp), dimension(${BOUND('num_dims')}$, ${BOUND('num_dims')}$) :: current_tau_bulk !< Current bulk stress tensor. + real(wp), dimension(${BOUND('num_dims')}$) :: vel_src_at_interface !< Interface velocities (u,v,w) for viscous work. + integer, dimension(3) :: idx_right_phys !< Physical (j,k,l) indices for right state. + real(wp) :: Re_shear !< Interface shear Reynolds number. + real(wp) :: Re_bulk !< Interface bulk Reynolds number. + integer :: j_loop !< Physical x-index loop iterator. + integer :: k_loop !< Physical y-index loop iterator. + integer :: l_loop !< Physical z-index loop iterator. + integer :: i_dim !< Generic dimension/component iterator. + integer :: vel_comp_idx !< Velocity component iterator (1=u, 2=v, 3=w). + real(wp) :: divergence_v !< Velocity divergence at interface. + real(wp) :: gamma_dot, D_xx, D_yy, D_zz, D_xy, D_xz, D_yz + real(wp), dimension(2) :: Re_nn real(wp), dimension(num_fluids) :: alpha_avg - integer :: fl + integer :: fl $:GPU_PARALLEL_LOOP(collapse=3, private='[idx_right_phys, vel_grad_avg, current_tau_shear, current_tau_bulk, & & vel_src_at_interface, Re_shear, Re_bulk, divergence_v, i_dim, vel_comp_idx, gamma_dot, D_xx, D_yy, & @@ -1163,15 +1146,10 @@ contains $:GPU_ROUTINE(parallelism='[seq]') ! Arguments - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3, 3), intent(in) :: vel_grad_avg - real(wp), dimension(3, 3), intent(out) :: tau_shear_out - #:else - real(wp), dimension(num_dims, num_dims), intent(in) :: vel_grad_avg - real(wp), dimension(num_dims, num_dims), intent(out) :: tau_shear_out - #:endif - real(wp), intent(in) :: Re_shear - real(wp), intent(in) :: divergence_v + real(wp), dimension(${BOUND('num_dims')}$, ${BOUND('num_dims')}$), intent(in) :: vel_grad_avg + real(wp), dimension(${BOUND('num_dims')}$, ${BOUND('num_dims')}$), intent(out) :: tau_shear_out + real(wp), intent(in) :: Re_shear + real(wp), intent(in) :: divergence_v ! Local variables integer :: i_dim !< Loop iterator for face normal. @@ -1195,13 +1173,9 @@ contains $:GPU_ROUTINE(parallelism='[seq]') ! Arguments - real(wp), intent(in) :: Re_bulk - real(wp), intent(in) :: divergence_v - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3, 3), intent(out) :: tau_bulk_out - #:else - real(wp), dimension(num_dims, num_dims), intent(out) :: tau_bulk_out - #:endif + real(wp), intent(in) :: Re_bulk + real(wp), intent(in) :: divergence_v + real(wp), dimension(${BOUND('num_dims')}$, ${BOUND('num_dims')}$), intent(out) :: tau_bulk_out ! Local variables integer :: i_dim !< Loop iterator for diagonal components. @@ -1219,12 +1193,8 @@ contains $:GPU_ROUTINE(function_name='s_compute_interface_reynolds', parallelism='[seq]', cray_inline=True) - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3), intent(in) :: alpha_K - #:else - real(wp), dimension(num_fluids), intent(in) :: alpha_K - #:endif - real(wp), dimension(2), intent(out) :: Re_K + real(wp), dimension(${BOUND('num_fluids')}$), intent(in) :: alpha_K + real(wp), dimension(2), intent(out) :: Re_K !> host copies of Re_size; amdflang reads the declare-target original stale cross-TU integer, intent(in) :: Re_size_loc1, Re_size_loc2 integer :: i, q !< Loop iterators diff --git a/src/simulation/m_sim_helpers.fpp b/src/simulation/m_sim_helpers.fpp index 82e56e943a..c2d84f75a8 100644 --- a/src/simulation/m_sim_helpers.fpp +++ b/src/simulation/m_sim_helpers.fpp @@ -50,25 +50,16 @@ contains $:GPU_ROUTINE(function_name='s_compute_cell_state',parallelism='[seq]', cray_inline=True) - type(scalar_field), intent(in), dimension(sys_size) :: q_prim_vf - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), intent(inout), dimension(3) :: alpha, alpha_rho - real(wp), intent(inout), dimension(3) :: vel - #:else - real(wp), intent(inout), dimension(num_fluids) :: alpha, alpha_rho - real(wp), intent(inout), dimension(num_vels) :: vel - #:endif - real(wp), intent(inout) :: rho, gamma, pi_inf, vel_sum, pres - real(wp), intent(out) :: qv - integer, intent(in) :: j, k, l - real(wp), dimension(2), intent(inout) :: Re - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3) :: Gs - #:else - real(wp), dimension(num_fluids) :: Gs - #:endif - real(wp) :: G_local - integer :: i + type(scalar_field), intent(in), dimension(sys_size) :: q_prim_vf + real(wp), intent(inout), dimension(${BOUND('num_fluids')}$) :: alpha, alpha_rho + real(wp), intent(inout), dimension(${BOUND('num_vels')}$) :: vel + real(wp), intent(inout) :: rho, gamma, pi_inf, vel_sum, pres + real(wp), intent(out) :: qv + integer, intent(in) :: j, k, l + real(wp), dimension(2), intent(inout) :: Re + real(wp), dimension(${BOUND('num_fluids')}$) :: Gs + real(wp) :: G_local + integer :: i call s_compute_species_fraction(q_prim_vf, j, k, l, alpha_rho, alpha) @@ -108,20 +99,16 @@ contains subroutine s_compute_stability_from_dt(vel, c, rho, Re_l, alpha, alpha_rho, j, k, l, icfl, vcfl, Rc, ccfl, tcfl) $:GPU_ROUTINE(parallelism='[seq]') - real(wp), intent(in), dimension(num_vels) :: vel - real(wp), intent(in) :: c, rho - real(wp), intent(inout) :: icfl - real(wp), intent(inout) :: vcfl, Rc, ccfl, tcfl - real(wp), dimension(2), intent(in) :: Re_l - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3), intent(in) :: alpha, alpha_rho - #:else - real(wp), dimension(num_fluids), intent(in) :: alpha, alpha_rho - #:endif - integer, intent(in) :: j, k, l - real(wp) :: fltr_dtheta - real(wp) :: k_mix, rho_cv - integer :: i + real(wp), intent(in), dimension(num_vels) :: vel + real(wp), intent(in) :: c, rho + real(wp), intent(inout) :: icfl + real(wp), intent(inout) :: vcfl, Rc, ccfl, tcfl + real(wp), dimension(2), intent(in) :: Re_l + real(wp), dimension(${BOUND('num_fluids')}$), intent(in) :: alpha, alpha_rho + integer, intent(in) :: j, k, l + real(wp) :: fltr_dtheta + real(wp) :: k_mix, rho_cv + integer :: i ! Inviscid CFL calculation ! The multi-dimensional CFL terms are written out here rather than @@ -214,20 +201,16 @@ contains subroutine s_compute_dt_from_cfl(vel, c, max_dt, rho, Re_l, alpha, alpha_rho, j, k, l) $:GPU_ROUTINE(parallelism='[seq]') - real(wp), dimension(num_vels), intent(in) :: vel - real(wp), intent(in) :: c, rho - real(wp), dimension(4), intent(out) :: max_dt - real(wp), dimension(2), intent(in) :: Re_l - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3), intent(in) :: alpha, alpha_rho - #:else - real(wp), dimension(num_fluids), intent(in) :: alpha, alpha_rho - #:endif - integer, intent(in) :: j, k, l - real(wp) :: vcfl_dt, ccfl_dt, tcfl_dt - real(wp) :: fltr_dtheta - real(wp) :: k_mix, rho_cv - integer :: i + real(wp), dimension(num_vels), intent(in) :: vel + real(wp), intent(in) :: c, rho + real(wp), dimension(4), intent(out) :: max_dt + real(wp), dimension(2), intent(in) :: Re_l + real(wp), dimension(${BOUND('num_fluids')}$), intent(in) :: alpha, alpha_rho + integer, intent(in) :: j, k, l + real(wp) :: vcfl_dt, ccfl_dt, tcfl_dt + real(wp) :: fltr_dtheta + real(wp) :: k_mix, rho_cv + integer :: i max_dt(2) = huge(1._wp) max_dt(3) = huge(1._wp) diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index df311bca2b..ee620a415d 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -826,15 +826,9 @@ contains integer :: i call s_initialize_global_parameters_module() - #:if USING_AMD - #:for BC in [-5, -6, -7, -8, -9, -10, -11, -12, -13] - @:PROHIBIT(any((/bc_x%beg, bc_x%end, bc_y%beg, bc_y%end, bc_z%beg, & - & bc_z%end/) == ${BC}$) .and. eqn_idx%adv%end > 70 .and. (.not. chemistry), & - & "CBC module with AMD compiler requires eqn_idx%adv%end <= 70 when case optimization is turned off") - @:PROHIBIT(any((/bc_x%beg, bc_x%end, bc_y%beg, bc_y%end, bc_z%beg, & - & bc_z%end/) == ${BC}$) .and. sys_size > 70 .and. (chemistry), & - & "CBC module with AMD compiler and chemistry requires sys_size <= 70 when case optimization is turned off") - #:endfor + #:if MFC_FIXED_BOUNDS + ! sys_size exists only from here on, so this cannot live in s_check_fixed_bounds. + @:PROHIBIT(sys_size > ${SYS_SIZE_MAX}$, "sys_size <= ${SYS_SIZE_MAX}$ in GPU builds") #:endif if (bubbles_euler .or. bubbles_lagrange) then call s_initialize_bubbles_model() diff --git a/src/simulation/m_surface_tension.fpp b/src/simulation/m_surface_tension.fpp index 2fbe62261b..573336fa3b 100644 --- a/src/simulation/m_surface_tension.fpp +++ b/src/simulation/m_surface_tension.fpp @@ -63,12 +63,8 @@ contains $:GPU_ROUTINE(function_name='s_compute_capillary_stress_tensor', parallelism='[seq]', cray_inline=True) - real(wp), intent(in) :: sigma_c, w1, w2, w3, normW - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3, 3), intent(inout) :: Omega - #:else - real(wp), dimension(num_dims, num_dims), intent(inout) :: Omega - #:endif + real(wp), intent(in) :: sigma_c, w1, w2, w3, normW + real(wp), dimension(${BOUND('num_dims')}$, ${BOUND('num_dims')}$), intent(inout) :: Omega Omega(1, 1) = -sigma_c*(w2*w2 + w3*w3)/normW #:if not MFC_CASE_OPTIMIZATION or num_dims > 1 @@ -94,19 +90,14 @@ contains subroutine s_compute_capillary_source_flux(vSrc_rsx_vf, flux_src_vf, id, isx, isy, isz) - real(wp), dimension(-1:,-1:,-1:,1:), intent(in) :: vSrc_rsx_vf - type(scalar_field), dimension(sys_size), intent(inout) :: flux_src_vf - integer, intent(in) :: id - type(int_bounds_info), intent(in) :: isx, isy, isz - - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3, 3) :: Omega - #:else - real(wp), dimension(num_dims, num_dims) :: Omega - #:endif - real(wp) :: w1L, w1R, w2L, w2R, w3L, w3R, w1, w2, w3 - real(wp) :: normWL, normWR, normW - integer :: j, k, l, i + real(wp), dimension(-1:,-1:,-1:,1:), intent(in) :: vSrc_rsx_vf + type(scalar_field), dimension(sys_size), intent(inout) :: flux_src_vf + integer, intent(in) :: id + type(int_bounds_info), intent(in) :: isx, isy, isz + real(wp), dimension(${BOUND('num_dims')}$, ${BOUND('num_dims')}$) :: Omega + real(wp) :: w1L, w1R, w2L, w2R, w3L, w3R, w1, w2, w3 + real(wp) :: normWL, normWR, normW + integer :: j, k, l, i if (id == 1) then $:GPU_PARALLEL_LOOP(collapse=3, private='[Omega, w1L, w2L, w3L, w1R, w2R, w3R, w1, w2, w3, normWL, normWR, normW]') diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp index 5fddc3f132..6f2c6a5d09 100644 --- a/src/simulation/m_time_steppers.fpp +++ b/src/simulation/m_time_steppers.fpp @@ -663,29 +663,23 @@ contains impure subroutine s_compute_dt() real(wp) :: rho !< Cell-avg. density - - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3) :: vel !< Cell-avg. velocity - real(wp), dimension(3) :: alpha, alpha_rho !< Cell-avg. volume fraction, partial density - #:else - real(wp), dimension(num_vels) :: vel !< Cell-avg. velocity - real(wp), dimension(num_fluids) :: alpha, alpha_rho !< Cell-avg. volume fraction, partial density - #:endif - real(wp) :: vel_sum !< Cell-avg. velocity sum - real(wp) :: pres !< Cell-avg. pressure - real(wp) :: gamma !< Cell-avg. sp. heat ratio - real(wp) :: pi_inf !< Cell-avg. liquid stiffness function - real(wp) :: qv !< Cell-avg. fluid reference energy - real(wp) :: c !< Cell-avg. sound speed - real(wp), dimension(2) :: Re !< Cell-avg. Reynolds numbers - real(wp), dimension(4) :: max_dt !< Cell dt candidates (inviscid, viscous, capillary, thermal) - real(wp) :: icfl_dt_local, vcfl_dt_local, ccfl_dt_local, tcfl_dt_local, coll_dt_local + real(wp), dimension(${BOUND('num_vels')}$) :: vel !< Cell-avg. velocity + real(wp), dimension(${BOUND('num_fluids')}$) :: alpha, alpha_rho !< Cell-avg. volume fraction, partial density + real(wp) :: vel_sum !< Cell-avg. velocity sum + real(wp) :: pres !< Cell-avg. pressure + real(wp) :: gamma !< Cell-avg. sp. heat ratio + real(wp) :: pi_inf !< Cell-avg. liquid stiffness function + real(wp) :: qv !< Cell-avg. fluid reference energy + real(wp) :: c !< Cell-avg. sound speed + real(wp), dimension(2) :: Re !< Cell-avg. Reynolds numbers + real(wp), dimension(4) :: max_dt !< Cell dt candidates (inviscid, viscous, capillary, thermal) + real(wp) :: icfl_dt_local, vcfl_dt_local, ccfl_dt_local, tcfl_dt_local, coll_dt_local real(wp), dimension(5) :: dt_candidates_loc !< Rank-local dt candidates (ICFL, VCFL, CCFL, TCFL, collision cap) real(wp), dimension(5) :: dt_candidates_glb !< Global dt candidates (ICFL, VCFL, CCFL, TCFL, collision cap) - real(wp) :: dt_prev - logical :: is_fluid_cell !< Cell lies outside every immersed boundary - integer :: j, k, l !< Generic loop iterators - integer :: fl !< Fluid loop iterator + real(wp) :: dt_prev + logical :: is_fluid_cell !< Cell lies outside every immersed boundary + integer :: j, k, l !< Generic loop iterators + integer :: fl !< Fluid loop iterator if (.not. igr) then call s_convert_conservative_to_primitive_variables(q_cons_ts(1)%vf, q_T_sf, q_prim_vf, idwint) diff --git a/src/simulation/m_viscous.fpp b/src/simulation/m_viscous.fpp index 37afd8988e..b4c51ac030 100644 --- a/src/simulation/m_viscous.fpp +++ b/src/simulation/m_viscous.fpp @@ -58,16 +58,12 @@ contains $:GPU_ROUTINE(function_name='s_compute_axis_inv_re', parallelism='[seq]', cray_inline=True) - type(scalar_field), dimension(num_dims), intent(in) :: grad_x_vf, grad_y_vf, grad_z_vf - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3), intent(in) :: alpha_visc - #:else - real(wp), dimension(num_fluids), intent(in) :: alpha_visc - #:endif - integer, intent(in) :: j, k, l - real(wp), dimension(2), intent(out) :: Re_visc - real(wp) :: gamma_dot_c - integer :: i, q + type(scalar_field), dimension(num_dims), intent(in) :: grad_x_vf, grad_y_vf, grad_z_vf + real(wp), dimension(${BOUND('num_fluids')}$), intent(in) :: alpha_visc + integer, intent(in) :: j, k, l + real(wp), dimension(2), intent(out) :: Re_visc + real(wp) :: gamma_dot_c + integer :: i, q if (any_non_newtonian) then gamma_dot_c = f_compute_shear_rate_from_components(grad_x_vf(1)%sf(j, k, l), grad_y_vf(2)%sf(j, k, l), 0._wp, & @@ -106,16 +102,9 @@ contains type(int_bounds_info), intent(in) :: ix, iy, iz real(wp) :: rho_visc, gamma_visc, pi_inf_visc, qv_visc, alpha_visc_sum !< Mixture variables real(wp), dimension(2) :: Re_visc - - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(3) :: alpha_visc, alpha_rho_visc - real(wp), dimension(3, 3) :: tau_Re - #:else - real(wp), dimension(num_fluids) :: alpha_visc, alpha_rho_visc - real(wp), dimension(num_dims, num_dims) :: tau_Re - #:endif - - integer :: i, j, k, l, q !< Generic loop iterator + real(wp), dimension(${BOUND('num_fluids')}$) :: alpha_visc, alpha_rho_visc + real(wp), dimension(${BOUND('num_dims')}$, ${BOUND('num_dims')}$) :: tau_Re + integer :: i, j, k, l, q !< Generic loop iterator is1_viscous = ix; is2_viscous = iy; is3_viscous = iz diff --git a/src/simulation/m_weno.fpp b/src/simulation/m_weno.fpp index ef18917abf..69d20c8be3 100644 --- a/src/simulation/m_weno.fpp +++ b/src/simulation/m_weno.fpp @@ -912,31 +912,21 @@ contains !> Perform WENO reconstruction of left and right cell-boundary values from cell-averaged variables subroutine s_weno(v_vf, vL_rs_vf_x, vR_rs_vf_x, weno_dir, is1_weno_d, is2_weno_d, is3_weno_d) - type(scalar_field), dimension(1:), intent(in) :: v_vf + type(scalar_field), dimension(1:), intent(in) :: v_vf real(wp), dimension(idwbuff(1)%beg:,idwbuff(2)%beg:,idwbuff(3)%beg:,1:), intent(inout) :: vL_rs_vf_x real(wp), dimension(idwbuff(1)%beg:,idwbuff(2)%beg:,idwbuff(3)%beg:,1:), intent(inout) :: vR_rs_vf_x - integer, intent(in) :: weno_dir - type(int_bounds_info), intent(in) :: is1_weno_d, is2_weno_d, is3_weno_d - - #:if not MFC_CASE_OPTIMIZATION and USING_AMD - real(wp), dimension(-3:2) :: dvd - real(wp), dimension(0:4) :: poly - real(wp), dimension(0:4) :: alpha - real(wp), dimension(0:4) :: omega - real(wp), dimension(0:4) :: beta - real(wp), dimension(0:4) :: delta - #:else - real(wp), dimension(-weno_polyn:weno_polyn - 1) :: dvd - real(wp), dimension(0:weno_num_stencils) :: poly - real(wp), dimension(0:weno_num_stencils) :: alpha - real(wp), dimension(0:weno_num_stencils) :: omega - real(wp), dimension(0:weno_num_stencils) :: beta - real(wp), dimension(0:weno_num_stencils) :: delta - #:endif + integer, intent(in) :: weno_dir + type(int_bounds_info), intent(in) :: is1_weno_d, is2_weno_d, is3_weno_d + real(wp), dimension(-${BOUND('weno_polyn')}$:${BOUND('weno_polyn')}$ - 1) :: dvd + real(wp), dimension(0:${BOUND('weno_num_stencils')}$) :: poly + real(wp), dimension(0:${BOUND('weno_num_stencils')}$) :: alpha + real(wp), dimension(0:${BOUND('weno_num_stencils')}$) :: omega + real(wp), dimension(0:${BOUND('weno_num_stencils')}$) :: beta + real(wp), dimension(0:${BOUND('weno_num_stencils')}$) :: delta real(wp), dimension(-3:3) :: v !< temporary field value array for clarity (WENO7 only) - real(wp) :: tau - integer :: i, j, k, l, q - real(wp) :: vp0, vp1, vp2, vp3, vm1, vm2, vm3 + real(wp) :: tau + integer :: i, j, k, l, q + real(wp) :: vp0, vp1, vp2, vp3, vm1, vm2, vm3 is1_weno = is1_weno_d is2_weno = is2_weno_d diff --git a/toolchain/mfc/lint_source.py b/toolchain/mfc/lint_source.py index 4093fcf7d4..4112b69d9a 100644 --- a/toolchain/mfc/lint_source.py +++ b/toolchain/mfc/lint_source.py @@ -46,7 +46,7 @@ # RUNTIME_CHECK_MARKER instead, so anything new added there is still flagged. RUNTIME_CHECKER_SUBROUTINES = { # Compiler conditionals (#ifdef / #if guarded). - "s_check_amd", + "s_check_fixed_bounds", "s_check_inputs_compilers", "s_check_inputs_nvidia_uvm", # MPI decomposition: n_global, num_procs_y/z. From 77bf251f9d2b8e5859b10fec00a3a6fe0c441dba Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Fri, 2 Oct 2026 23:51:49 -0500 Subject: [PATCH 2/5] Generate the state-vector layout (eqn_idx, sys_size) from one Python table toolchain/mfc/params/eqn_layout.py lists each model's fields in order with their size and condition as small Python expressions. The params codegen translates them into generated_eqn_idx.fpp, which s_initialize_eqn_idx now includes in place of 140 hand-written lines (the hypoelastic shear tables stay hand-written). evaluate_layout() computes the same layout for a case, which case optimization will use for a compile-time sys_size. Checked against master's routine over all 7,962,624 combinations of model, the 13 layout flags, 1D/multi-D, and the fluid/velocity/dimension/bubble/ species counts: every eqn_idx field and sys_size match, and evaluate_layout agrees on sys_size for a 300,000-case sample. --- cmake/ParamsCodegen.cmake | 5 +- docs/documentation/contributing.md | 4 +- src/common/m_global_parameters_common.fpp | 165 +++-------------- toolchain/mfc/params/eqn_layout.py | 168 ++++++++++++++++++ .../mfc/params/generators/fortran_gen.py | 6 +- toolchain/mfc/params_tests/test_eqn_layout.py | 54 ++++++ .../mfc/params_tests/test_fortran_gen.py | 8 +- 7 files changed, 261 insertions(+), 149 deletions(-) create mode 100644 toolchain/mfc/params/eqn_layout.py create mode 100644 toolchain/mfc/params_tests/test_eqn_layout.py diff --git a/cmake/ParamsCodegen.cmake b/cmake/ParamsCodegen.cmake index a27b081882..a5f83c2dc3 100644 --- a/cmake/ParamsCodegen.cmake +++ b/cmake/ParamsCodegen.cmake @@ -9,7 +9,7 @@ file(GLOB_RECURSE _mfc_gen_inputs "${CMAKE_CURRENT_SOURCE_DIR}/toolchain/mfc/params/*.py" ) -# Enumerate the 18 generated .fpp files explicitly so ninja can track them as +# Enumerate the 21 generated .fpp files explicitly so ninja can track them as # build-time outputs and so HANDLE_SOURCES does not need a configure-time GLOB # of ${CMAKE_BINARY_DIR}/include// (which fails when the dir is empty). set(_mfc_gen_inc "${CMAKE_BINARY_DIR}/include") @@ -18,6 +18,7 @@ set(_mfc_gen_files_pre_process "${_mfc_gen_inc}/pre_process/generated_decls.fpp" "${_mfc_gen_inc}/pre_process/generated_constants.fpp" "${_mfc_gen_inc}/pre_process/generated_eos.fpp" + "${_mfc_gen_inc}/pre_process/generated_eqn_idx.fpp" "${_mfc_gen_inc}/pre_process/generated_bcast.fpp" "${_mfc_gen_inc}/pre_process/generated_case_opt_decls.fpp" ) @@ -26,6 +27,7 @@ set(_mfc_gen_files_simulation "${_mfc_gen_inc}/simulation/generated_decls.fpp" "${_mfc_gen_inc}/simulation/generated_constants.fpp" "${_mfc_gen_inc}/simulation/generated_eos.fpp" + "${_mfc_gen_inc}/simulation/generated_eqn_idx.fpp" "${_mfc_gen_inc}/simulation/generated_bcast.fpp" "${_mfc_gen_inc}/simulation/generated_case_opt_decls.fpp" ) @@ -34,6 +36,7 @@ set(_mfc_gen_files_post_process "${_mfc_gen_inc}/post_process/generated_decls.fpp" "${_mfc_gen_inc}/post_process/generated_constants.fpp" "${_mfc_gen_inc}/post_process/generated_eos.fpp" + "${_mfc_gen_inc}/post_process/generated_eqn_idx.fpp" "${_mfc_gen_inc}/post_process/generated_bcast.fpp" "${_mfc_gen_inc}/post_process/generated_case_opt_decls.fpp" ) diff --git a/docs/documentation/contributing.md b/docs/documentation/contributing.md index 9767705a0e..b190b701a2 100644 --- a/docs/documentation/contributing.md +++ b/docs/documentation/contributing.md @@ -87,7 +87,7 @@ between them is narrow. mfc.sh (env bootstrap, venv, module loading, lock) └─ toolchain/mfc/build.py (config slugs, cmake invocation) └─ CMakeLists.txt + cmake/{GPU,Fypp,ParamsCodegen,MFCTargets}.cmake - └─ toolchain/mfc/params/generators/cmake_gen.py (writes 18 generated .fpp includes) + └─ toolchain/mfc/params/generators/cmake_gen.py (writes 21 generated .fpp includes) ``` **`mfc.sh` → `build.py`.** `mfc.sh` is a thin shell wrapper that activates the Python @@ -100,7 +100,7 @@ mode, debug, chemistry, MPI. Staging and install trees are namespaced by slug u **CMake layer.** `cmake/Fypp.cmake` defines `HANDLE_SOURCES`, which sets up one `add_custom_command` per `.fpp` file to run Fypp at build time. `cmake/ParamsCodegen.cmake` registers a single ninja-tracked `add_custom_command` (DEPENDS all `params/*.py`) that -invokes `cmake_gen.py` and writes the 18 generated includes under +invokes `cmake_gen.py` and writes the 21 generated includes under `build/include//`. There is no configure-time generation: all 18 files are build outputs, so changing any `params/*.py` triggers only a targeted rebuild, not a full reconfigure. diff --git a/src/common/m_global_parameters_common.fpp b/src/common/m_global_parameters_common.fpp index 80e3c2923a..46e4ce573f 100644 --- a/src/common/m_global_parameters_common.fpp +++ b/src/common/m_global_parameters_common.fpp @@ -100,8 +100,8 @@ module m_global_parameters_common contains - !> Initialize equation-index state (eqn_idx and sys_size) from the namelist parameters. This is the shared skeleton: it covers - !! the model_eqns dispatch, all eqn_idx field assignments, and the hypoelastic/surface-tension/chemistry extensions. + !> Initialize equation-index state (eqn_idx and sys_size) from the namelist parameters. The layout itself is generated from + !! toolchain/mfc/params/eqn_layout.py, which case optimization also evaluates; edit fields there, not here. !! !! @param nmom_in Number of carried moments per R0 location (per-target: pre/post pass an !! integer variable; sim passes its integer parameter nmom = 6). Used only in the 5eq @@ -117,146 +117,29 @@ contains integer, intent(in) :: nb_in logical, intent(in) :: six_eqn_alf_is_advected - ! Gamma/Pi_inf Model - - if (model_eqns == model_eqns_gamma_law) then - ! Annotating structure of the state and flux vectors belonging to the system of - ! equations defined by the selected number of spatial dimensions and the gamma/pi_inf model - eqn_idx%cont%beg = 1 - eqn_idx%cont%end = eqn_idx%cont%beg - eqn_idx%mom%beg = eqn_idx%cont%end + 1 - eqn_idx%mom%end = eqn_idx%cont%end + num_vels - eqn_idx%E = eqn_idx%mom%end + 1 - eqn_idx%adv%beg = eqn_idx%E + 1 - eqn_idx%adv%end = eqn_idx%adv%beg + 1 - eqn_idx%gamma = eqn_idx%adv%beg - eqn_idx%pi_inf = eqn_idx%adv%end - sys_size = eqn_idx%adv%end - - ! Volume Fraction Model (5-equation model) - else if (model_eqns == model_eqns_5eq) then - ! Annotating structure of the state and flux vectors belonging to the system of - ! equations defined by the selected number of spatial dimensions and the volume fraction model - eqn_idx%cont%beg = 1 - eqn_idx%cont%end = num_fluids - eqn_idx%mom%beg = eqn_idx%cont%end + 1 - eqn_idx%mom%end = eqn_idx%cont%end + num_vels - eqn_idx%E = eqn_idx%mom%end + 1 - - if (igr) then - ! IGR: volume fractions after energy (N-1 for N fluids; skipped when num_fluids=1) - eqn_idx%adv%beg = eqn_idx%E + 1 - eqn_idx%adv%end = eqn_idx%E + num_fluids - 1 - else - ! WENO/MUSCL + Riemann tracks a total of (N) volume fractions for N fluids - eqn_idx%adv%beg = eqn_idx%E + 1 - eqn_idx%adv%end = eqn_idx%E + num_fluids - end if - - sys_size = eqn_idx%adv%end - - if (bubbles_euler) then - eqn_idx%alf = eqn_idx%adv%end - else - eqn_idx%alf = 1 - end if - - if (bubbles_euler) then - eqn_idx%bub%beg = sys_size + 1 - if (qbmm) then - eqn_idx%bub%end = eqn_idx%adv%end + nb_in*nmom_in - else - if (.not. polytropic) then - eqn_idx%bub%end = sys_size + 4*nb_in - else - eqn_idx%bub%end = sys_size + 2*nb_in - end if - end if - sys_size = eqn_idx%bub%end - - if (adv_n) then - eqn_idx%n = eqn_idx%bub%end + 1 - sys_size = eqn_idx%n - end if - end if - - if (mhd) then - eqn_idx%B%beg = sys_size + 1 - if (n == 0) then - eqn_idx%B%end = sys_size + 2 ! 1D: By, Bz - else - eqn_idx%B%end = sys_size + 3 ! 2D/3D: Bx, By, Bz - end if - sys_size = eqn_idx%B%end - end if - - ! Volume Fraction Model (6-equation model) - else if (model_eqns == model_eqns_6eq) then - ! Annotating structure of the state and flux vectors belonging to the system of - ! equations defined by the selected number of spatial dimensions and the volume fraction model - eqn_idx%cont%beg = 1 - eqn_idx%cont%end = num_fluids - eqn_idx%mom%beg = eqn_idx%cont%end + 1 - eqn_idx%mom%end = eqn_idx%cont%end + num_vels - eqn_idx%E = eqn_idx%mom%end + 1 - eqn_idx%adv%beg = eqn_idx%E + 1 - eqn_idx%adv%end = eqn_idx%E + num_fluids - if (six_eqn_alf_is_advected) eqn_idx%alf = eqn_idx%adv%end - eqn_idx%int_en%beg = eqn_idx%adv%end + 1 - eqn_idx%int_en%end = eqn_idx%adv%end + num_fluids - sys_size = eqn_idx%int_en%end - end if - - if (model_eqns == model_eqns_5eq .or. model_eqns == model_eqns_6eq) then - if (hypoelasticity) then - eqn_idx%stress%beg = sys_size + 1 - eqn_idx%stress%end = sys_size + (num_dims*(num_dims + 1))/2 - if (cyl_coord) eqn_idx%stress%end = eqn_idx%stress%end + 1 - ! number of stresses is 1 in 1D, 3 in 2D, 4 in 2D-Axisym, 6 in 3D - sys_size = eqn_idx%stress%end - - ! shear stress index is 2 for 2D and 2,4,5 for 3D. Readers test the whole array - ! rather than the first shear_num entries, so unused slots must not be garbage. - shear_indices = 0 - if (num_dims == 1) then - shear_num = 0 - else if (num_dims == 2) then - shear_num = 1 - shear_indices(1) = eqn_idx%stress%beg - 1 + 2 - shear_BC_flip_num = 1 - shear_BC_flip_indices(1:2,1) = shear_indices(1) - ! Both x-dir and y-dir: flip tau_xy only - else if (num_dims == 3) then - shear_num = 3 - shear_indices(1:3) = eqn_idx%stress%beg - 1 + (/2, 4, 5/) - shear_BC_flip_num = 2 - shear_BC_flip_indices(1,1:2) = shear_indices((/1, 2/)) - shear_BC_flip_indices(2,1:2) = shear_indices((/1, 3/)) - shear_BC_flip_indices(3,1:2) = shear_indices((/2, 3/)) - ! x-dir: flip tau_xy and tau_xz; y-dir: flip tau_xy and tau_yz; z-dir: flip tau_xz and tau_yz - end if + #:include 'generated_eqn_idx.fpp' + + ! shear stress index is 2 for 2D and 2,4,5 for 3D. Readers test the whole array + ! rather than the first shear_num entries, so unused slots must not be garbage. + if ((model_eqns == model_eqns_5eq .or. model_eqns == model_eqns_6eq) .and. hypoelasticity) then + shear_indices = 0 + if (num_dims == 1) then + shear_num = 0 + else if (num_dims == 2) then + shear_num = 1 + shear_indices(1) = eqn_idx%stress%beg - 1 + 2 + shear_BC_flip_num = 1 + shear_BC_flip_indices(1:2,1) = shear_indices(1) + ! Both x-dir and y-dir: flip tau_xy only + else if (num_dims == 3) then + shear_num = 3 + shear_indices(1:3) = eqn_idx%stress%beg - 1 + (/2, 4, 5/) + shear_BC_flip_num = 2 + shear_BC_flip_indices(1,1:2) = shear_indices((/1, 2/)) + shear_BC_flip_indices(2,1:2) = shear_indices((/1, 3/)) + shear_BC_flip_indices(3,1:2) = shear_indices((/2, 3/)) + ! x-dir: flip tau_xy and tau_xz; y-dir: flip tau_xy and tau_yz; z-dir: flip tau_xz and tau_yz end if - - if (surface_tension) then - eqn_idx%c = sys_size + 1 - sys_size = eqn_idx%c - end if - - if (cont_damage) then - eqn_idx%damage = sys_size + 1 - sys_size = eqn_idx%damage - end if - - if (hyper_cleaning) then - eqn_idx%psi = sys_size + 1 - sys_size = eqn_idx%psi - end if - end if - - if (chemistry) then - eqn_idx%species%beg = sys_size + 1 - eqn_idx%species%end = sys_size + num_species - sys_size = eqn_idx%species%end end if ! Resolved here, not with the other fluid properties, because the MPI halo buffers are sized before diff --git a/toolchain/mfc/params/eqn_layout.py b/toolchain/mfc/params/eqn_layout.py new file mode 100644 index 0000000000..bb748f3d1b --- /dev/null +++ b/toolchain/mfc/params/eqn_layout.py @@ -0,0 +1,168 @@ +"""Layout of the state vector (eqn_idx and sys_size), defined once. + +Each model appends fields in order; a field occupies `size` slots (one if `size` is None) when its +`when` condition holds. Sizes, conditions, and alias values are Python expressions over the +s_initialize_eqn_idx inputs and over earlier fields (`adv.end`). They are both translated to the +Fortran of s_initialize_eqn_idx (fortran_layout) and evaluated for a given case (evaluate_layout), +so the runtime layout and a case-optimized sys_size cannot disagree. +""" + +import ast +from dataclasses import dataclass +from typing import Dict, List, Optional + + +@dataclass(frozen=True) +class Field: + name: str + size: Optional[str] = None # None: a single index (eqn_idx%name); else a range (eqn_idx%name%beg:end) + when: Optional[str] = None + + +@dataclass(frozen=True) +class Alias: + name: str + value: str + when: Optional[str] = None + + +HEAD = [Field("mom", "num_vels"), Field("E")] + +MODELS: Dict[str, List] = { + "gamma_law": [Field("cont", "1"), *HEAD, Field("adv", "2"), Alias("gamma", "adv.beg"), Alias("pi_inf", "adv.end")], + "5eq": [ + Field("cont", "num_fluids"), + *HEAD, + Field("adv", "num_fluids - 1 if igr else num_fluids"), + Alias("alf", "adv.end if bubbles_euler else 1"), + Field("bub", "nb_in*nmom_in if qbmm else (2 if polytropic else 4)*nb_in", "bubbles_euler"), + Field("n", None, "bubbles_euler and adv_n"), + Field("B", "2 if n == 0 else 3", "mhd"), + ], + "6eq": [ + Field("cont", "num_fluids"), + *HEAD, + Field("adv", "num_fluids"), + Alias("alf", "adv.end", "six_eqn_alf_is_advected"), + Field("int_en", "num_fluids"), + ], +} + +# Appended after the model fields: for the five- and six-equation models, then for every model. +MULTIPHASE = [ + Field("stress", "num_dims*(num_dims + 1)//2 + (1 if cyl_coord else 0)", "hypoelasticity"), + Field("c", None, "surface_tension"), + Field("damage", None, "cont_damage"), + Field("psi", None, "hyper_cleaning"), +] +ALL_MODELS = [Field("species", "num_species", "chemistry")] + + +class _ToFortran(ast.NodeVisitor): + OPS = {ast.Add: "+", ast.Sub: "-", ast.Mult: "*", ast.FloorDiv: "/", ast.Eq: "==", ast.NotEq: "/="} + + def visit_Expression(self, node): + return self.visit(node.body) + + def visit_Name(self, node): + return {"True": ".true.", "False": ".false."}.get(node.id, node.id) + + def visit_Constant(self, node): + return {True: ".true.", False: ".false."}[node.value] if isinstance(node.value, bool) else str(node.value) + + def visit_Attribute(self, node): + return f"eqn_idx%{node.value.id}%{node.attr}" + + def visit_BinOp(self, node): + return f"({self.visit(node.left)} {self.OPS[type(node.op)]} {self.visit(node.right)})" + + def visit_Compare(self, node): + return f"({self.visit(node.left)} {self.OPS[type(node.ops[0])]} {self.visit(node.comparators[0])})" + + def visit_BoolOp(self, node): + op = " .and. " if isinstance(node.op, ast.And) else " .or. " + return "(" + op.join(self.visit(v) for v in node.values) + ")" + + def visit_UnaryOp(self, node): + return f"(.not. {self.visit(node.operand)})" + + def visit_IfExp(self, node): + return f"merge({self.visit(node.body)}, {self.visit(node.orelse)}, {self.visit(node.test)})" + + def generic_visit(self, node): + raise ValueError(f"eqn_layout: unsupported expression node {type(node).__name__}") + + +def _f(expr: str) -> str: + """Translate a layout expression to Fortran.""" + return _ToFortran().visit(ast.parse(expr, mode="eval")) + + +def _fields_fortran(fields, ind: str) -> List[str]: + out = [] + for f in fields: + body = [] + if isinstance(f, Alias): + body.append(f"eqn_idx%{f.name} = {_f(f.value)}") + elif f.size is None: + body += [f"eqn_idx%{f.name} = sys_size + 1", f"sys_size = eqn_idx%{f.name}"] + else: + body += [f"eqn_idx%{f.name}%beg = sys_size + 1", f"eqn_idx%{f.name}%end = sys_size + {_f(f.size)}", f"sys_size = eqn_idx%{f.name}%end"] + if f.when: + cond = _f(f.when) + out += [f"{ind}if ({cond[1:-1] if cond.startswith('(') else cond}) then", *(f"{ind} {b}" for b in body), f"{ind}end if"] + else: + out += [f"{ind}{b}" for b in body] + return out + + +def fortran_layout() -> str: + """Fortran statements setting eqn_idx and sys_size (included by s_initialize_eqn_idx).""" + lines = ["! Generated by toolchain/mfc/params/eqn_layout.py. Do not edit.", "sys_size = 0"] + for i, (model, fields) in enumerate(MODELS.items()): + lines.append(f"{'if' if i == 0 else 'else if'} (model_eqns == model_eqns_{model}) then") + lines += _fields_fortran(fields, " ") + lines.append("end if") + lines.append("if (model_eqns == model_eqns_5eq .or. model_eqns == model_eqns_6eq) then") + lines += _fields_fortran(MULTIPHASE, " ") + lines.append("end if") + lines += _fields_fortran(ALL_MODELS, "") + return "\n".join(lines) + "\n" + + +class _Ranges(dict): + """Lets layout expressions read earlier ranges as `adv.beg` / `adv.end`.""" + + def __getattr__(self, key): + return self[key] + + +def evaluate_layout(inputs: Dict) -> Dict: + """Evaluate the layout for given inputs (model_eqns as 'gamma_law'/'5eq'/'6eq', flags, counts). + + Returns {"sys_size": int, field: index or (beg, end), alias: index}. + """ + env = dict(inputs) + out: Dict = {} + size = 0 + + def run(fields): + nonlocal size + for f in fields: + if f.when and not eval(f.when, {}, env): + continue + if isinstance(f, Alias): + out[f.name] = eval(f.value, {}, env) + continue + beg, end = size + 1, size + (1 if f.size is None else int(eval(f.size, {}, env))) + out[f.name] = beg if f.size is None else (beg, end) + if f.size is not None: # only ranges are referenced (adv.end); `n` is also the grid size + env[f.name] = _Ranges(beg=beg, end=end) + size = end + + run(MODELS[inputs["model_eqns"]]) + if inputs["model_eqns"] in ("5eq", "6eq"): + run(MULTIPHASE) + run(ALL_MODELS) + out["sys_size"] = size + return out diff --git a/toolchain/mfc/params/generators/fortran_gen.py b/toolchain/mfc/params/generators/fortran_gen.py index 605c97a9fa..c29b441b9a 100644 --- a/toolchain/mfc/params/generators/fortran_gen.py +++ b/toolchain/mfc/params/generators/fortran_gen.py @@ -5,6 +5,7 @@ from typing import List, Tuple from ..definitions import CASE_OPT_PARAMS, DECLARATION_TARGETS, FORTRAN_ARRAY_DIMS, NAMELIST_VARS, TYPED_DECLS # noqa: F401 - triggers registry population +from ..eqn_layout import fortran_layout from ..registry import REGISTRY from ..schema import ParamDef, ParamType @@ -852,10 +853,10 @@ def resolve_namelist_content(fpp_path: Path) -> str: def get_generated_files(build_dir: Path) -> List[Tuple[Path, str]]: - """Return (path, content) for all 18 generated .fpp files under build_dir. + """Return (path, content) for all 21 generated .fpp files under build_dir. Paths match the cmake include directory structure: - build_dir/include/{full_target}/generated_{namelist,decls,constants,eos,case_opt_decls,bcast}.fpp + build_dir/include/{full_target}/generated_{namelist,decls,constants,eos,eqn_idx,case_opt_decls,bcast}.fpp Every target gets generated_case_opt_decls.fpp: the full case-optimization block for simulation and common computed-scalar declarations for pre/post. Every target gets generated_bcast.fpp with its MPI broadcast statements. @@ -867,6 +868,7 @@ def get_generated_files(build_dir: Path) -> List[Tuple[Path, str]]: result.append((inc / "generated_decls.fpp", generate_decls_fpp(short))) result.append((inc / "generated_constants.fpp", generate_constants_fpp())) result.append((inc / "generated_eos.fpp", generate_eos_fpp())) + result.append((inc / "generated_eqn_idx.fpp", fortran_layout())) sim_gpu_decls = "" for short, full in TARGETS: inc = build_dir / "include" / full diff --git a/toolchain/mfc/params_tests/test_eqn_layout.py b/toolchain/mfc/params_tests/test_eqn_layout.py new file mode 100644 index 0000000000..ec80b276d6 --- /dev/null +++ b/toolchain/mfc/params_tests/test_eqn_layout.py @@ -0,0 +1,54 @@ +BASE = dict( + num_fluids=2, + num_vels=2, + num_dims=2, + nb_in=1, + nmom_in=6, + num_species=10, + n=5, + igr=False, + bubbles_euler=False, + qbmm=False, + polytropic=True, + adv_n=False, + mhd=False, + hypoelasticity=False, + cyl_coord=False, + surface_tension=False, + cont_damage=False, + hyper_cleaning=False, + chemistry=False, + six_eqn_alf_is_advected=False, +) + + +def layout(**kw): + from mfc.params.eqn_layout import evaluate_layout + + return evaluate_layout({**BASE, "model_eqns": "5eq", **kw}) + + +def test_five_equation_layout(): + out = layout() + assert (out["cont"], out["mom"], out["E"], out["adv"]) == ((1, 2), (3, 4), 5, (6, 7)) + assert out["sys_size"] == 7 and out["alf"] == 1 + + +def test_extensions_append_in_order(): + out = layout(bubbles_euler=True, polytropic=False, nb_in=3, hypoelasticity=True, chemistry=True) + assert out["bub"] == (8, 19) and out["alf"] == 7 + assert out["stress"] == (20, 22) and out["species"] == (23, 32) and out["sys_size"] == 32 + + +def test_grid_n_is_not_the_bubble_index(): + # eqn_idx%n (bubble number density) must not shadow the grid size n used by the MHD field count + out = layout(bubbles_euler=True, adv_n=True, mhd=True, n=0) + assert out["B"][1] - out["B"][0] + 1 == 2 + + +def test_fortran_translation(): + from mfc.params.eqn_layout import _f, fortran_layout + + assert _f("adv.end if bubbles_euler else 1") == "merge(eqn_idx%adv%end, 1, bubbles_euler)" + assert _f("num_dims*(num_dims + 1)//2") == "((num_dims * (num_dims + 1)) / 2)" + assert "eqn_idx%species%end = sys_size + num_species" in fortran_layout() diff --git a/toolchain/mfc/params_tests/test_fortran_gen.py b/toolchain/mfc/params_tests/test_fortran_gen.py index 0c82351333..f17783ceaf 100644 --- a/toolchain/mfc/params_tests/test_fortran_gen.py +++ b/toolchain/mfc/params_tests/test_fortran_gen.py @@ -258,13 +258,13 @@ def test_check_target_raises_on_bad_target(): generate_decls_fpp("bad") -def test_get_generated_files_returns_eighteen(): +def test_get_generated_files_returns_twenty_one(): from pathlib import Path from mfc.params.generators.fortran_gen import get_generated_files files = get_generated_files(Path("/build")) - assert len(files) == 18 + assert len(files) == 21 paths = [str(p) for p, _ in files] assert any("pre_process/generated_namelist.fpp" in p for p in paths) assert any("simulation/generated_decls.fpp" in p for p in paths) @@ -273,6 +273,7 @@ def test_get_generated_files_returns_eighteen(): assert any("pre_process/generated_bcast.fpp" in p for p in paths) assert any("simulation/generated_bcast.fpp" in p for p in paths) assert any("post_process/generated_bcast.fpp" in p for p in paths) + assert any("simulation/generated_eqn_idx.fpp" in p for p in paths) def test_generate_constants_fpp_content(): @@ -299,10 +300,11 @@ def test_get_generated_files_includes_bcast(): "generated_decls.fpp", "generated_constants.fpp", "generated_eos.fpp", + "generated_eqn_idx.fpp", "generated_case_opt_decls.fpp", "generated_bcast.fpp", } - assert len(files) == 18 + assert len(files) == 21 def test_generate_case_opt_decls_fpp(): From d93b45b644b4d4ba7886989eb832eca848f4398e Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Sat, 3 Oct 2026 00:00:00 -0500 Subject: [PATCH 3/5] Bake the case-optimized sys_size into per-thread array extents Case optimization now evaluates the layout table for the case and writes CASE_OPT_SIZES (sys_size and the hypoelastic stress count) into case.fpp. BOUND() returns those exact values, or the case-optimized parameters, under case optimization on any compiler; before, sys_size arrays fell back to the runtime value or, on amdflang, to the fixed maximum of 70. The sys_size variable itself stays a runtime variable filled by the generated layout, so nothing that declares or updates it on the device changes. Simulation checks at startup that the runtime layout matches the baked-in values. CASE_OPT_SIZES is part of case.fpp, so it is already in the build fingerprint: a case with a different layout gets its own build. --- src/common/include/shared_parallel_macros.fpp | 12 ++++--- src/simulation/m_start_up.fpp | 8 +++-- toolchain/mfc/case.py | 31 +++++++++++++++++++ 3 files changed, 44 insertions(+), 7 deletions(-) diff --git a/src/common/include/shared_parallel_macros.fpp b/src/common/include/shared_parallel_macros.fpp index 1bd13b650e..4c5a1ff7c8 100644 --- a/src/common/include/shared_parallel_macros.fpp +++ b/src/common/include/shared_parallel_macros.fpp @@ -11,9 +11,9 @@ #! NUM_SPECIES (m_thermochem's species count) and CHEMISTRY, written per build by the toolchain. #:include 'thermochem.fpp' -#! Fixed bounds: in GPU simulation builds (MFC_FIXED_BOUNDS, set per target by CMake) per-thread -#! arrays get compile-time extents, since runtime-sized private arrays spill to scratch. Case -#! optimization makes most of these compile-time anyway; sys_size never is, so it is always fixed. +#! Compile-time extents for per-thread arrays, since runtime-sized private arrays spill to scratch +#! on GPUs. Case optimization gives exact ones (parameters, plus CASE_OPT_SIZES from the layout +#! table); otherwise GPU simulation builds (MFC_FIXED_BOUNDS, set per target by CMake) use maxima. #! Chemistry pins num_fluids to 1, leaving at most 10 flow variables beside the species. #:set MFC_FIXED_BOUNDS = defined('MFC_FIXED_BOUNDS') and MFC_FIXED_BOUNDS #:set FIXED_BOUNDS = MFC_FIXED_BOUNDS and not MFC_CASE_OPTIMIZATION @@ -24,9 +24,11 @@ & 'weno_polyn': 3, 'weno_num_stencils': 4, 'nterms': 32, 'n_stress': 6, 'sys_size': SYS_SIZE_MAX} #:set BOUND_RUNTIME = {'n_stress': 'eqn_idx%stress%end - eqn_idx%stress%beg + 1'} -#! Extent for a per-thread array sized by `name`: its fixed maximum or the runtime expression. +#! Extent for a per-thread array sized by `name`: exact under case optimization, else the fixed +#! maximum in GPU simulation builds, else the runtime expression. #:def BOUND(name) - $:BOUND_MAX[name] if (MFC_FIXED_BOUNDS if name == 'sys_size' else FIXED_BOUNDS) else BOUND_RUNTIME.get(name, name) + $:CASE_OPT_SIZES.get(name, & + & name) if MFC_CASE_OPTIMIZATION else BOUND_MAX[name] if MFC_FIXED_BOUNDS else BOUND_RUNTIME.get(name, name) #:enddef #:def ASSERT_LIST(data, datatype) diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index ee620a415d..ccf1a43f1c 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -826,8 +826,12 @@ contains integer :: i call s_initialize_global_parameters_module() - #:if MFC_FIXED_BOUNDS - ! sys_size exists only from here on, so this cannot live in s_check_fixed_bounds. + ! sys_size exists only from here on, so these cannot live in the input checks. + #:if MFC_CASE_OPTIMIZATION + @:PROHIBIT(sys_size /= ${CASE_OPT_SIZES['sys_size']}$, "sys_size differs from the case-optimized build; rebuild it") + @:PROHIBIT(hypoelasticity .and. eqn_idx%stress%end - eqn_idx%stress%beg + 1 /= ${CASE_OPT_SIZES['n_stress']}$, & + & "stress count differs from the case-optimized build; rebuild it") + #:elif MFC_FIXED_BOUNDS @:PROHIBIT(sys_size > ${SYS_SIZE_MAX}$, "sys_size <= ${SYS_SIZE_MAX}$ in GPU builds") #:endif if (bubbles_euler .or. bubbles_lagrange) then diff --git a/toolchain/mfc/case.py b/toolchain/mfc/case.py index 0359e45230..1947446c60 100644 --- a/toolchain/mfc/case.py +++ b/toolchain/mfc/case.py @@ -334,6 +334,34 @@ def __get_analytic_mib_fpp(self, print: bool) -> str: """ return content + def __case_opt_sizes(self, num_dims: int, num_vels: int) -> dict: + """Compile-time extents case optimization bakes in beside the parameters: sys_size and the + hypoelastic stress count, from the same layout table s_initialize_eqn_idx is generated from.""" + from .params.eqn_layout import evaluate_layout + + def flag(name: str, default: str = "F") -> bool: + return self.params.get(name, default) == "T" + + chemistry = flag("chemistry") + layout = evaluate_layout( + { + "model_eqns": {1: "gamma_law", 2: "5eq", 3: "6eq"}[int(self.params.get("model_eqns", 2))], + **{k: flag(k) for k in ("igr", "bubbles_euler", "qbmm", "adv_n", "mhd", "hypoelasticity", "cyl_coord", "surface_tension", "cont_damage", "hyper_cleaning")}, + "polytropic": flag("polytropic", "T"), + "chemistry": chemistry, + "six_eqn_alf_is_advected": True, + "num_fluids": int(self.params.get("num_fluids", 1)), + "num_dims": num_dims, + "num_vels": num_vels, + "nb_in": int(self.params.get("nb", 1)), + "nmom_in": 6, + "n": int(self.params.get("n", 0)), + "num_species": self.get_cantera_solution().n_species if chemistry else 0, + } + ) + stress = layout.get("stress", (1, 0)) + return {"sys_size": layout["sys_size"], "n_stress": max(stress[1] - stress[0] + 1, 1)} + def __get_sim_fpp(self, print: bool) -> str: if ARG("case_optimization"): if print: @@ -388,9 +416,12 @@ def __get_sim_fpp(self, print: bool) -> str: igr = 1 if self.params.get("igr", "F") == "T" else 0 igr_pres_lim = 1 if self.params.get("igr_pres_lim", "F") == "T" else 0 + sizes = self.__case_opt_sizes(num_dims, num_vels) + # Throw error if wenoz_q is required but not set out = f"""\ #:set MFC_CASE_OPTIMIZATION = {ARG("case_optimization")} +#:set CASE_OPT_SIZES = {sizes} #:set recon_type = {recon_type} #:set weno_order = {weno_order} #:set weno_polyn = {weno_polyn} From b507fe520f8b156e16c35f1f1452d61672bc92e2 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Sat, 3 Oct 2026 00:07:31 -0500 Subject: [PATCH 4/5] Use fixed per-thread array bounds in every GPU simulation build MFC_FIXED_BOUNDS now applies to NVHPC and CCE (OpenACC and OpenMP), not just amdflang. CPU builds keep runtime extents: with GNU the fixed bounds cost up to 39% on the 5-equation cases. In #1938's CI benchmarks the fixed bounds sped up CCE by 12-36% and NVHPC OpenMP by 9-62%, with AMD unchanged. GPU simulation builds without case optimization now take at most 3 fluids, nb <= 3, and sys_size <= 70 (or 10 + species); larger cases rebuild with --case-optimization, and s_check_fixed_bounds says so. Docs describe BOUND() and why runtime-sized per-thread arrays are slow on GPUs. --- cmake/Fypp.cmake | 3 +-- docs/documentation/gpuParallelization.md | 27 ++++++++++-------------- 2 files changed, 12 insertions(+), 18 deletions(-) diff --git a/cmake/Fypp.cmake b/cmake/Fypp.cmake index 4a4b3a24b8..19f3a18173 100644 --- a/cmake/Fypp.cmake +++ b/cmake/Fypp.cmake @@ -103,8 +103,7 @@ macro(HANDLE_SOURCES target useCommon) # Fixed per-thread array bounds (see shared_parallel_macros.fpp): only simulation is offloaded. set(_fixed_bounds False) - if (MFC_FIXED_BOUNDS AND (MFC_OpenACC OR MFC_OpenMP) AND "${target}" STREQUAL "simulation" - AND CMAKE_Fortran_COMPILER_ID STREQUAL "LLVMFlang") + if (MFC_FIXED_BOUNDS AND (MFC_OpenACC OR MFC_OpenMP) AND "${target}" STREQUAL "simulation") set(_fixed_bounds True) endif() diff --git a/docs/documentation/gpuParallelization.md b/docs/documentation/gpuParallelization.md index 4a03b2f23c..b661a11f9a 100644 --- a/docs/documentation/gpuParallelization.md +++ b/docs/documentation/gpuParallelization.md @@ -650,14 +650,12 @@ helper receives only scalars or small arrays with explicit-shape dimensioning. T is required because the `SF` indexing lambda used in solver loops is defined locally inside each solver's `#:for NORM_DIR` block and cannot be referenced from a helper. -**AMD case-opt compatibility.** Under `--case-optimization` with the AMD backend, -arrays that are sized by runtime parameters at compile time must be declared with an -explicit constant bound. Use an explicit `n` argument (e.g., `integer, intent(in) :: -nf`) and dimension helpers as `dimension(nf)` rather than `dimension(num_fluids)`. -See `s_compute_interface_reynolds` in `src/simulation/m_riemann_state.fpp` for the -`#:if not MFC_CASE_OPTIMIZATION and USING_AMD` guard pattern: the guard sits on the -dummy-argument declaration in the helper's definition, with matching guards on the -callers' own local declarations so the actual and dummy bounds agree. +**Per-thread array extents.** Size per-thread arrays with ``${BOUND('name')}$`` +(`src/common/include/shared_parallel_macros.fpp`), e.g. ``dimension(${BOUND('num_fluids')}$)``, +not `dimension(num_fluids)`. Use it on both a helper's dummy arguments and its callers' locals +so the bounds agree. `BOUND` gives the exact value under `--case-optimization`, a fixed maximum +in GPU simulation builds (`MFC_FIXED_BOUNDS`), and the runtime extent otherwise. A new extent +needs an entry in `BOUND_MAX` and, if it can exceed it, a check in `s_check_fixed_bounds`. **Declare scoping.** The `GPU_ROUTINE` directive must appear in the source file that defines the routine. Helpers added to `m_riemann_state.fpp` are automatically @@ -940,14 +938,11 @@ answer is wrong, or one backend diverges from all the others. passes a `parameter` array from `m_thermochem`, such as `molecular_weights`, into a declare-target routine. Read such arrays directly in the kernel, or pass a plain local computed from them. -- **The `USING_AMD` fypp guards are load-bearing, not a stale workaround.** They swap a - device-global array bound for a literal in `src/common/include/shared_parallel_macros.fpp` - and its 86 use sites. Setting `USING_AMD = False` and rebuilding amdflang `--gpu mp` - without case optimization compiles completely clean, then produces NaNs in CBC, the - `wave_speeds=2` Riemann path, immersed boundaries, surface tension, QBMM and viscous - cases, and MHD HLLD, while both Lagrange bubble cases complete with out-of-tolerance - answers. A compile-only check returns green, so any attempt to remove these must run the - tests rather than just build. +- **Runtime-sized per-thread arrays are slow on GPUs.** A private array sized by a runtime value + (`num_fluids`, `sys_size`) cannot live in registers and spills to scratch: on amdflang (AFAR + 24.3) the WENO kernel ran 11x slower and HLLC 4.9x. This is why `BOUND` gives fixed maxima in + GPU builds; turning them off (`-DMFC_FIXED_BOUNDS=OFF`) needs a benchmark, not just the tests. + On AFAR 23.2.x the runtime-sized form also gave NaNs, so passing tests alone is not enough. - `@:ACC_SETUP_VFs` and `@:ACC_SETUP_SFs` compile only under Cray. Around MPI, use `GPU_UPDATE(host=...)` before a send and `GPU_UPDATE(device=...)` after a receive. From 41a73aed8cfb4d319da71cddc0ffb89694109e52 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Sun, 4 Oct 2026 23:56:45 -0400 Subject: [PATCH 5/5] Size host callers of fixed-bound routines with BOUND In GPU simulation builds BOUND('num_fluids') is 3, but three host paths still passed num_fluids-sized arrays to routines with BOUND dummies: s_convert_species_to_mixture_variables (its alpha_K/alpha_rho_K locals and the forwarded G), s_report_icfl_violation, and s_write_probe_files. With one or two fluids and mpp_lim, the kernel's whole-array normalization of alpha_K wrote past the end of the host local. Declare those with BOUND (every G caller passes fluid_pp(:)%G, sized num_fluids_max), and normalize only alpha_K(1:num_fluids). Also: fix the generated-file count in contributing.md (21, not 18), and tell users to rebuild with --case-optimization when sys_size exceeds the GPU maximum. --- docs/documentation/contributing.md | 2 +- src/common/m_variables_conversion.fpp | 36 ++++++++-------- src/simulation/m_data_output.fpp | 62 +++++++++++++-------------- src/simulation/m_start_up.fpp | 2 +- 4 files changed, 51 insertions(+), 51 deletions(-) diff --git a/docs/documentation/contributing.md b/docs/documentation/contributing.md index b190b701a2..1c9a265444 100644 --- a/docs/documentation/contributing.md +++ b/docs/documentation/contributing.md @@ -101,7 +101,7 @@ mode, debug, chemistry, MPI. Staging and install trees are namespaced by slug u `add_custom_command` per `.fpp` file to run Fypp at build time. `cmake/ParamsCodegen.cmake` registers a single ninja-tracked `add_custom_command` (DEPENDS all `params/*.py`) that invokes `cmake_gen.py` and writes the 21 generated includes under -`build/include//`. There is no configure-time generation: all 18 files are build +`build/include//`. There is no configure-time generation: all 21 files are build outputs, so changing any `params/*.py` triggers only a targeted rebuild, not a full reconfigure. diff --git a/src/common/m_variables_conversion.fpp b/src/common/m_variables_conversion.fpp index 0fabe8af0a..2452565b5f 100644 --- a/src/common/m_variables_conversion.fpp +++ b/src/common/m_variables_conversion.fpp @@ -53,12 +53,12 @@ contains !! procedure pointer. subroutine s_convert_to_mixture_variables(q_vf, i, j, k, rho, gamma, pi_inf, qv, Re_K, G_K, G) - type(scalar_field), dimension(sys_size), intent(in) :: q_vf - integer, intent(in) :: i, j, k - real(wp), intent(out), target :: rho, gamma, pi_inf, qv - real(wp), optional, dimension(2), intent(out) :: Re_K - real(wp), optional, intent(out) :: G_K - real(wp), optional, dimension(num_fluids), intent(in) :: G + type(scalar_field), dimension(sys_size), intent(in) :: q_vf + integer, intent(in) :: i, j, k + real(wp), intent(out), target :: rho, gamma, pi_inf, qv + real(wp), optional, dimension(2), intent(out) :: Re_K + real(wp), optional, intent(out) :: G_K + real(wp), optional, dimension(${BOUND('num_fluids')}$), intent(in) :: G if (model_eqns == model_eqns_gamma_law) then ! Gamma/pi_inf model call s_convert_mixture_to_mixture_variables(q_vf, i, j, k, rho, gamma, pi_inf, qv) @@ -153,17 +153,17 @@ contains !! stores the results into rho, gamma and pi_inf. subroutine s_convert_species_to_mixture_variables(q_vf, k, l, r, rho, gamma, pi_inf, qv, Re_K, G_K, G) - type(scalar_field), dimension(sys_size), intent(in) :: q_vf - integer, intent(in) :: k, l, r - real(wp), intent(out), target :: rho - real(wp), intent(out), target :: gamma - real(wp), intent(out), target :: pi_inf - real(wp), intent(out), target :: qv - real(wp), optional, dimension(2), intent(out) :: Re_K - real(wp), optional, intent(out) :: G_K - real(wp), dimension(num_fluids) :: alpha_rho_K, alpha_K - real(wp), optional, dimension(num_fluids), intent(in) :: G - integer :: i, j !< Generic loop iterator + type(scalar_field), dimension(sys_size), intent(in) :: q_vf + integer, intent(in) :: k, l, r + real(wp), intent(out), target :: rho + real(wp), intent(out), target :: gamma + real(wp), intent(out), target :: pi_inf + real(wp), intent(out), target :: qv + real(wp), optional, dimension(2), intent(out) :: Re_K + real(wp), optional, intent(out) :: G_K + real(wp), dimension(${BOUND('num_fluids')}$) :: alpha_rho_K, alpha_K + real(wp), optional, dimension(${BOUND('num_fluids')}$), intent(in) :: G + integer :: i, j !< Generic loop iterator ! Computing the density, the specific heat ratio function and the liquid stiffness function, respectively call s_compute_species_fraction(q_vf, k, l, r, alpha_rho_K, alpha_K) @@ -209,7 +209,7 @@ contains alpha_K(i) = min(max(0._wp, alpha_K(i)), 1._wp) alpha_K_sum = alpha_K_sum + alpha_K(i) end do - alpha_K = alpha_K/max(alpha_K_sum, sgm_eps) + alpha_K(1:num_fluids) = alpha_K(1:num_fluids)/max(alpha_K_sum, sgm_eps) end if call s_compute_mixture_coefficients(alpha_rho_K, alpha_K, rho_K, gamma_K, pi_inf_K, qv_K) diff --git a/src/simulation/m_data_output.fpp b/src/simulation/m_data_output.fpp index ceb06ffe49..4e17abd400 100644 --- a/src/simulation/m_data_output.fpp +++ b/src/simulation/m_data_output.fpp @@ -347,7 +347,7 @@ contains impure subroutine s_report_icfl_violation(q_prim_vf) type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf - real(wp), dimension(num_fluids) :: alpha, alpha_rho + real(wp), dimension(${BOUND('num_fluids')}$) :: alpha, alpha_rho real(wp), dimension(num_vels) :: vel, vel_hit real(wp), dimension(2) :: Re real(wp) :: rho, vel_sum, pres, gamma, pi_inf, qv, c @@ -1435,36 +1435,36 @@ contains ! The cell-averaged partial densities, density, velocity, pressure, volume fractions, specific heat ratio function, liquid ! stiffness function, and sound speed. - real(wp) :: lit_gamma, nbub - real(wp) :: rho - real(wp), dimension(num_vels) :: vel - real(wp) :: pres - real(wp) :: ptilde - real(wp) :: ptot - real(wp) :: alf - real(wp) :: alfgr - real(wp), dimension(num_fluids) :: alpha, alpha_rho - real(wp) :: gamma - real(wp) :: pi_inf - real(wp) :: qv - real(wp) :: c - real(wp) :: M00, M10, M01, M20, M02 - real(wp) :: varR, varV - real(wp), dimension(Nb) :: nR, R, nRdot, Rdot - real(wp) :: nR3 - real(wp) :: accel - real(wp) :: int_pres - real(wp) :: max_pres - real(wp), dimension(2) :: Re - real(wp), dimension(6) :: tau_e - real(wp) :: G_undamaged, G_damaged - real(wp) :: dyn_p, T - real(wp) :: damage_state - real(wp) :: solid_partial_density !< damageable-solid partial density at the probe cell - integer :: i, j, k, l, s, d !< Generic loop iterator - real(wp) :: nondim_time !< Non-dimensional time - real(wp) :: tmp !< Temporary variable to store quantity for mpi_allreduce - real(wp) :: rhoYks(1:num_species) + real(wp) :: lit_gamma, nbub + real(wp) :: rho + real(wp), dimension(num_vels) :: vel + real(wp) :: pres + real(wp) :: ptilde + real(wp) :: ptot + real(wp) :: alf + real(wp) :: alfgr + real(wp), dimension(${BOUND('num_fluids')}$) :: alpha, alpha_rho + real(wp) :: gamma + real(wp) :: pi_inf + real(wp) :: qv + real(wp) :: c + real(wp) :: M00, M10, M01, M20, M02 + real(wp) :: varR, varV + real(wp), dimension(Nb) :: nR, R, nRdot, Rdot + real(wp) :: nR3 + real(wp) :: accel + real(wp) :: int_pres + real(wp) :: max_pres + real(wp), dimension(2) :: Re + real(wp), dimension(6) :: tau_e + real(wp) :: G_undamaged, G_damaged + real(wp) :: dyn_p, T + real(wp) :: damage_state + real(wp) :: solid_partial_density !< damageable-solid partial density at the probe cell + integer :: i, j, k, l, s, d !< Generic loop iterator + real(wp) :: nondim_time !< Non-dimensional time + real(wp) :: tmp !< Temporary variable to store quantity for mpi_allreduce + real(wp) :: rhoYks(1:num_species) T = dflt_T_guess diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index ccf1a43f1c..536400bd7f 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -832,7 +832,7 @@ contains @:PROHIBIT(hypoelasticity .and. eqn_idx%stress%end - eqn_idx%stress%beg + 1 /= ${CASE_OPT_SIZES['n_stress']}$, & & "stress count differs from the case-optimized build; rebuild it") #:elif MFC_FIXED_BOUNDS - @:PROHIBIT(sys_size > ${SYS_SIZE_MAX}$, "sys_size <= ${SYS_SIZE_MAX}$ in GPU builds") + @:PROHIBIT(sys_size > ${SYS_SIZE_MAX}$, "sys_size <= ${SYS_SIZE_MAX}$ in GPU builds; rebuild with --case-optimization") #:endif if (bubbles_euler .or. bubbles_lagrange) then call s_initialize_bubbles_model()