From b6bf034c0daafd191e6399647009fca2a7f1d6bb Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Sat, 25 Jul 2026 20:52:08 -0400 Subject: [PATCH 01/15] Refactor variable conversion into common code --- src/common/m_global_parameters_common.fpp | 31 ++ src/common/m_variables_conversion.fpp | 268 +++++++----------- src/post_process/m_global_parameters.fpp | 6 - src/post_process/m_start_up.fpp | 2 +- src/pre_process/m_global_parameters.fpp | 2 - src/simulation/m_bubbles_EL.fpp | 2 +- src/simulation/m_cbc.fpp | 2 +- src/simulation/m_global_parameters.fpp | 29 -- src/simulation/m_hyperelastic.fpp | 2 +- src/simulation/m_ibm.fpp | 6 +- src/simulation/m_sim_helpers.fpp | 4 +- src/simulation/m_start_up.fpp | 2 +- toolchain/mfc/params/definitions.py | 24 +- .../mfc/params/generators/fortran_gen.py | 10 +- .../mfc/params_tests/test_fortran_gen.py | 25 ++ 15 files changed, 198 insertions(+), 217 deletions(-) diff --git a/src/common/m_global_parameters_common.fpp b/src/common/m_global_parameters_common.fpp index 4e65b887e0..dfcb0e1b3b 100644 --- a/src/common/m_global_parameters_common.fpp +++ b/src/common/m_global_parameters_common.fpp @@ -64,6 +64,29 @@ module m_global_parameters_common integer, dimension(3, 2) :: shear_BC_flip_indices !< Shear stress BC reflection indices (1:3, 1:shear_BC_flip_num) !> @} + !> @name Material properties derived from fluid_pp + !> @{ One declaration is shared by all executables and initialized by m_variables_conversion after the case parameters have been + !! read. + real(wp), allocatable, dimension(:) :: gammas, gs_min, pi_infs, ps_inf, cvs, qvs, qvps + $:GPU_DECLARE(create='[gammas, gs_min, pi_infs, ps_inf, cvs, qvs, qvps]') + !> @} + + !> @name Fluids participating in shear and bulk viscosity + !> @{ + integer, dimension(2) :: Re_size = 0 + integer :: Re_size_max = 0 + integer, allocatable, dimension(:,:) :: Re_idx + $:GPU_DECLARE(create='[Re_size, Re_size_max, Re_idx]') + !> @} + + !> @name Working-direction index mapping + !> @{ + integer, dimension(3) :: dir_idx + real(wp), dimension(3) :: dir_flg + integer, dimension(3) :: dir_idx_tau + $:GPU_DECLARE(create='[dir_idx, dir_flg, dir_idx_tau]') + !> @} + #ifdef MFC_SIMULATION $:GPU_DECLARE(create='[sys_size, eqn_idx, b_size, tensor_size]') $:GPU_DECLARE(create='[shear_num, shear_indices, shear_BC_flip_num, shear_BC_flip_indices]') @@ -434,6 +457,8 @@ contains igr_order = dflt_int mhd = .false. relativity = .false. + viscous = .false. + riemann_solver = dflt_int #:endif ! Tait EOS @@ -464,6 +489,12 @@ contains down_sample = .false. fft_wrt = .false. + ! Mixture conversion and sound-speed behavior + avg_state = dflt_int + alt_soundspeed = .false. + mixture_err = .false. + sigR = dflt_real + end subroutine s_assign_common_defaults end module m_global_parameters_common diff --git a/src/common/m_variables_conversion.fpp b/src/common/m_variables_conversion.fpp index fe17b89c27..e764e5bf46 100644 --- a/src/common/m_variables_conversion.fpp +++ b/src/common/m_variables_conversion.fpp @@ -21,29 +21,12 @@ module m_variables_conversion implicit none private - public :: s_initialize_variables_conversion_module, & - s_initialize_pb, & - s_initialize_mv, & - s_convert_to_mixture_variables, & - s_convert_mixture_to_mixture_variables, & - s_convert_species_to_mixture_variables, & - s_convert_species_to_mixture_variables_acc, & - s_convert_conservative_to_primitive_variables, & - s_convert_primitive_to_conservative_variables, & - s_convert_primitive_to_flux_variables, & - s_compute_pressure, & - s_compute_species_fraction, & -#ifndef MFC_PRE_PROCESS - s_compute_speed_of_sound, & - s_compute_fast_magnetosonic_speed, & -#endif - s_finalize_variables_conversion_module - - ! In simulation, gammas, pi_infs, and qvs are already declared in m_global_variables -#ifndef MFC_SIMULATION - real(wp), allocatable, public, dimension(:) :: gammas, gs_min, pi_infs, ps_inf, cvs, qvs, qvps - $:GPU_DECLARE(create='[gammas, gs_min, pi_infs, ps_inf, cvs, qvs, qvps]') -#endif + public :: s_initialize_variables_conversion_module, s_initialize_pb, s_initialize_mv, s_convert_to_mixture_variables, & + & s_convert_mixture_to_mixture_variables, s_convert_species_to_mixture_variables, & + & s_convert_species_to_mixture_variables_kernel, s_convert_conservative_to_primitive_variables, & + & s_convert_primitive_to_conservative_variables, s_convert_primitive_to_flux_variables, s_compute_pressure, & + & s_compute_species_fraction, s_compute_speed_of_sound, s_compute_fast_magnetosonic_speed, & + & s_finalize_variables_conversion_module, gammas, gs_min, pi_infs, ps_inf, cvs, qvs, qvps real(wp), allocatable, dimension(:) :: Gs_vc integer, allocatable, dimension(:) :: bubrs_vc @@ -53,6 +36,11 @@ module m_variables_conversion integer :: is1b, is2b, is3b, is1e, is2e, is3e $:GPU_DECLARE(create='[is1b, is2b, is3b, is1e, is2e, is3e]') + logical :: enforce_density_floor_vc = .false. + logical :: preserve_qbmm_number_vc = .false. + integer :: lagrange_beta_index_vc = 0 + $:GPU_DECLARE(create='[enforce_density_floor_vc, preserve_qbmm_number_vc, lagrange_beta_index_vc]') + real(wp), allocatable, dimension(:,:,:), public :: rho_sf !< Scalar density function real(wp), allocatable, dimension(:,:,:), public :: gamma_sf !< Scalar sp. heat ratio function real(wp), allocatable, dimension(:,:,:), public :: pi_inf_sf !< Scalar liquid stiffness function @@ -165,13 +153,13 @@ contains pi_inf = q_vf(eqn_idx%pi_inf)%sf(i, j, k) qv = 0._wp ! keep this value nil for now. For future adjustment - ! Post process requires rho_sf/gamma_sf/pi_inf_sf/qv_sf to also be updated -#ifdef MFC_POST_PROCESS - rho_sf(i, j, k) = rho - gamma_sf(i, j, k) = gamma - pi_inf_sf(i, j, k) = pi_inf - qv_sf(i, j, k) = qv -#endif + ! Store derived mixture fields when requested during module initialization. + if (allocated(rho_sf)) then + rho_sf(i, j, k) = rho + gamma_sf(i, j, k) = gamma + pi_inf_sf(i, j, k) = pi_inf + qv_sf(i, j, k) = qv + end if end subroutine s_convert_mixture_to_mixture_variables @@ -195,60 +183,33 @@ contains call s_compute_species_fraction(q_vf, k, l, r, alpha_rho_K, alpha_K) - ! Calculating the density, the specific heat ratio function, the liquid stiffness function, and the energy reference - ! function, respectively, from the species analogs - if (num_fluids == 1 .and. bubbles_euler) then - rho = alpha_rho_K(1) - gamma = gammas(1) - pi_inf = pi_infs(1) - qv = qvs(1) + ! Use the same scalar kernel on host and device so mixture semantics do not depend on the executable or accelerator backend. + if (present(G_K)) then + if (present(Re_K)) then + call s_convert_species_to_mixture_variables_kernel(rho, gamma, pi_inf, qv, alpha_K, alpha_rho_K, Re_K, G_K, G) + else + call s_convert_species_to_mixture_variables_kernel(rho, gamma, pi_inf, qv, alpha_K, alpha_rho_K, G_K=G_K, G=G) + end if + else if (present(Re_K)) then + call s_convert_species_to_mixture_variables_kernel(rho, gamma, pi_inf, qv, alpha_K, alpha_rho_K, Re_K) else - rho = 0._wp; gamma = 0._wp; pi_inf = 0._wp; qv = 0._wp - do i = 1, num_fluids - rho = rho + alpha_rho_K(i) - gamma = gamma + alpha_K(i)*gammas(i) - pi_inf = pi_inf + alpha_K(i)*pi_infs(i) - qv = qv + alpha_rho_K(i)*qvs(i) - end do + call s_convert_species_to_mixture_variables_kernel(rho, gamma, pi_inf, qv, alpha_K, alpha_rho_K) end if -#ifdef MFC_SIMULATION - ! Computing the shear and bulk Reynolds numbers from species analogs - if (viscous) then - do i = 1, 2 - Re_K(i) = dflt_real; if (Re_size(i) > 0) Re_K(i) = 0._wp - - do j = 1, Re_size(i) - Re_K(i) = alpha_K(Re_idx(i, j))/fluid_pp(Re_idx(i, j))%Re(i) + Re_K(i) - end do - - Re_K(i) = 1._wp/max(Re_K(i), sgm_eps) - end do + ! Store derived mixture fields when requested during module initialization. + if (allocated(rho_sf)) then + rho_sf(k, l, r) = rho + gamma_sf(k, l, r) = gamma + pi_inf_sf(k, l, r) = pi_inf + qv_sf(k, l, r) = qv end if -#endif - - if (present(G_K)) then - G_K = 0._wp - do i = 1, num_fluids - G_K = G_K + alpha_K(i)*G(i) - end do - G_K = max(0._wp, G_K) - end if - - ! Post process requires rho_sf/gamma_sf/pi_inf_sf/qv_sf to also be updated -#ifdef MFC_POST_PROCESS - rho_sf(k, l, r) = rho - gamma_sf(k, l, r) = gamma - pi_inf_sf(k, l, r) = pi_inf - qv_sf(k, l, r) = qv -#endif end subroutine s_convert_species_to_mixture_variables - !> GPU-accelerated conversion of species volume fractions and partial densities to mixture density, gamma, pi_inf, and qv. - subroutine s_convert_species_to_mixture_variables_acc(rho_K, gamma_K, pi_inf_K, qv_K, alpha_K, alpha_rho_K, Re_K, G_K, G) + !> Host- and device-callable conversion kernel for species and mixture variables. + subroutine s_convert_species_to_mixture_variables_kernel(rho_K, gamma_K, pi_inf_K, qv_K, alpha_K, alpha_rho_K, Re_K, G_K, G) - $:GPU_ROUTINE(function_name='s_convert_species_to_mixture_variables_acc', parallelism='[seq]', cray_noinline=True) + $: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 @@ -258,19 +219,18 @@ contains real(wp), dimension(num_fluids), intent(inout) :: alpha_rho_K, alpha_K real(wp), optional, dimension(num_fluids), intent(in) :: G #:endif - real(wp), 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), 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 pi_inf_K = 0._wp qv_K = 0._wp - Re_K = dflt_real + if (present(Re_K)) Re_K = dflt_real if (present(G_K)) G_K = 0._wp -#ifdef MFC_SIMULATION ! Constrain partial densities and volume fractions within physical bounds if (num_fluids == 1 .and. bubbles_euler) then rho_K = alpha_rho_K(1) @@ -305,7 +265,7 @@ contains G_K = max(0._wp, G_K) end if - if (viscous) then + if (viscous .and. present(Re_K)) then do i = 1, 2 Re_K(i) = dflt_real @@ -318,16 +278,30 @@ contains Re_K(i) = 1._wp/max(Re_K(i), sgm_eps) end do end if -#endif - end subroutine s_convert_species_to_mixture_variables_acc + end subroutine s_convert_species_to_mixture_variables_kernel !> Initialize the variables conversion module. - impure subroutine s_initialize_variables_conversion_module - - integer :: i, j + impure subroutine s_initialize_variables_conversion_module(store_mixture_fields, enforce_density_floor, preserve_qbmm_number, & + & lagrange_beta_index) + + integer :: i, j + logical, optional, intent(in) :: store_mixture_fields + logical, optional, intent(in) :: enforce_density_floor, preserve_qbmm_number + integer, optional, intent(in) :: lagrange_beta_index + logical :: allocate_mixture_fields + + allocate_mixture_fields = .false. + if (present(store_mixture_fields)) allocate_mixture_fields = store_mixture_fields + enforce_density_floor_vc = .false. + if (present(enforce_density_floor)) enforce_density_floor_vc = enforce_density_floor + preserve_qbmm_number_vc = .false. + if (present(preserve_qbmm_number)) preserve_qbmm_number_vc = preserve_qbmm_number + lagrange_beta_index_vc = 0 + if (present(lagrange_beta_index)) lagrange_beta_index_vc = lagrange_beta_index $:GPU_ENTER_DATA(copyin='[is1b, is1e, is2b, is2e, is3b, is3e]') + $:GPU_UPDATE(device='[enforce_density_floor_vc, preserve_qbmm_number_vc, lagrange_beta_index_vc]') @:ALLOCATE(gammas (1:num_fluids)) @:ALLOCATE(gs_min (1:num_fluids)) @@ -350,18 +324,17 @@ contains end do $:GPU_UPDATE(device='[gammas, gs_min, pi_infs, ps_inf, cvs, qvs, qvps, Gs_vc]') -#ifdef MFC_SIMULATION - if (viscous) then - @:ALLOCATE(Res_vc(1:2, 1:Re_size_max)) + @:ALLOCATE(Res_vc(1:2, 1:max(1, Re_size_max))) + Res_vc = dflt_real + if (allocated(Re_idx)) then do i = 1, 2 do j = 1, Re_size(i) Res_vc(i, j) = fluid_pp(Re_idx(i, j))%Re(i) end do end do - - $:GPU_UPDATE(device='[Res_vc, Re_idx, Re_size]') + $:GPU_UPDATE(device='[Re_idx]') end if -#endif + $:GPU_UPDATE(device='[Res_vc, Re_size]') if (bubbles_euler) then @:ALLOCATE(bubrs_vc(1:nb)) @@ -371,34 +344,27 @@ contains $:GPU_UPDATE(device='[bubrs_vc]') end if -#ifdef MFC_POST_PROCESS - ! Allocating the density, the specific heat ratio function and the liquid stiffness function, respectively - - ! Simulation is at least 2D - if (n > 0) then - ! Simulation is 3D - if (p > 0) then - allocate (rho_sf(-buff_size:m + buff_size,-buff_size:n + buff_size,-buff_size:p + buff_size)) - allocate (gamma_sf(-buff_size:m + buff_size,-buff_size:n + buff_size,-buff_size:p + buff_size)) - allocate (pi_inf_sf(-buff_size:m + buff_size,-buff_size:n + buff_size,-buff_size:p + buff_size)) - allocate (qv_sf(-buff_size:m + buff_size,-buff_size:n + buff_size,-buff_size:p + buff_size)) - - ! Simulation is 2D + if (allocate_mixture_fields) then + ! Allocate derived mixture fields over the available grid storage. + if (n > 0) then + if (p > 0) then + allocate (rho_sf(-buff_size:m + buff_size,-buff_size:n + buff_size,-buff_size:p + buff_size)) + allocate (gamma_sf(-buff_size:m + buff_size,-buff_size:n + buff_size,-buff_size:p + buff_size)) + allocate (pi_inf_sf(-buff_size:m + buff_size,-buff_size:n + buff_size,-buff_size:p + buff_size)) + allocate (qv_sf(-buff_size:m + buff_size,-buff_size:n + buff_size,-buff_size:p + buff_size)) + else + allocate (rho_sf(-buff_size:m + buff_size,-buff_size:n + buff_size,0:0)) + allocate (gamma_sf(-buff_size:m + buff_size,-buff_size:n + buff_size,0:0)) + allocate (pi_inf_sf(-buff_size:m + buff_size,-buff_size:n + buff_size,0:0)) + allocate (qv_sf(-buff_size:m + buff_size,-buff_size:n + buff_size,0:0)) + end if else - allocate (rho_sf(-buff_size:m + buff_size,-buff_size:n + buff_size,0:0)) - allocate (gamma_sf(-buff_size:m + buff_size,-buff_size:n + buff_size,0:0)) - allocate (pi_inf_sf(-buff_size:m + buff_size,-buff_size:n + buff_size,0:0)) - allocate (qv_sf(-buff_size:m + buff_size,-buff_size:n + buff_size,0:0)) + allocate (rho_sf(-buff_size:m + buff_size,0:0,0:0)) + allocate (gamma_sf(-buff_size:m + buff_size,0:0,0:0)) + allocate (pi_inf_sf(-buff_size:m + buff_size,0:0,0:0)) + allocate (qv_sf(-buff_size:m + buff_size,0:0,0:0)) end if - - ! Simulation is 1D - else - allocate (rho_sf(-buff_size:m + buff_size,0:0,0:0)) - allocate (gamma_sf(-buff_size:m + buff_size,0:0,0:0)) - allocate (pi_inf_sf(-buff_size:m + buff_size,0:0,0:0)) - allocate (qv_sf(-buff_size:m + buff_size,0:0,0:0)) end if -#endif end subroutine s_initialize_variables_conversion_module @@ -513,17 +479,17 @@ contains call s_compute_species_fraction(qK_cons_vf, j, k, l, alpha_rho_K, alpha_K) if (model_eqns /= model_eqns_4eq) then -#ifdef MFC_SIMULATION - ! If in simulation, use acc mixture subroutines +#ifdef MFC_GPU + ! Device regions call the device-compiled scalar kernel directly. if (elasticity) then - call s_convert_species_to_mixture_variables_acc(rho_K, gamma_K, pi_inf_K, qv_K, alpha_K, alpha_rho_K, & - & Re_K, G_K, Gs_vc) + call s_convert_species_to_mixture_variables_kernel(rho_K, gamma_K, pi_inf_K, qv_K, alpha_K, & + & alpha_rho_K, Re_K, G_K, Gs_vc) else - call s_convert_species_to_mixture_variables_acc(rho_K, gamma_K, pi_inf_K, qv_K, alpha_K, alpha_rho_K, & - & Re_K) + call s_convert_species_to_mixture_variables_kernel(rho_K, gamma_K, pi_inf_K, qv_K, alpha_K, & + & alpha_rho_K, Re_K) end if #else - ! If pre-processing, use non acc mixture subroutines + ! Host execution uses the wrapper, which also stores requested diagnostics. if (elasticity) then call s_convert_to_mixture_variables(qK_cons_vf, j, k, l, rho_K, gamma_K, pi_inf_K, qv_K, Re_K, G_K, & & fluid_pp(:)%G) @@ -635,9 +601,7 @@ contains end do end if -#ifdef MFC_SIMULATION - rho_K = max(rho_K, sgm_eps) -#endif + if (enforce_density_floor_vc) rho_K = max(rho_K, sgm_eps) ! Recover velocity from momentum: u = rho*u / rho, and accumulate dynamic pressure 0.5*rho*|u|^2 $:GPU_LOOP(parallelism='[seq]') @@ -700,9 +664,9 @@ contains qK_prim_vf(i)%sf(j, k, l) = qK_cons_vf(i)%sf(j, k, l)/nbub_sc end do ! Need to keep track of nb in the primitive variable list (converted back to true value before output) -#ifdef MFC_SIMULATION - qK_prim_vf(eqn_idx%bub%beg)%sf(j, k, l) = qK_cons_vf(eqn_idx%bub%beg)%sf(j, k, l) -#endif + if (preserve_qbmm_number_vc) then + qK_prim_vf(eqn_idx%bub%beg)%sf(j, k, l) = qK_cons_vf(eqn_idx%bub%beg)%sf(j, k, l) + end if else if (adv_n) then qK_prim_vf(eqn_idx%n)%sf(j, k, l) = qK_cons_vf(eqn_idx%n)%sf(j, k, l) @@ -770,9 +734,9 @@ contains if (cont_damage) qK_prim_vf(eqn_idx%damage)%sf(j, k, l) = qK_cons_vf(eqn_idx%damage)%sf(j, k, l) if (hyper_cleaning) qK_prim_vf(eqn_idx%psi)%sf(j, k, l) = qK_cons_vf(eqn_idx%psi)%sf(j, k, l) -#ifdef MFC_POST_PROCESS - if (bubbles_lagrange) qK_prim_vf(beta_idx)%sf(j, k, l) = qK_cons_vf(beta_idx)%sf(j, k, l) -#endif + if (bubbles_lagrange .and. lagrange_beta_index_vc > 0) then + qK_prim_vf(lagrange_beta_index_vc)%sf(j, k, l) = qK_cons_vf(lagrange_beta_index_vc)%sf(j, k, l) + end if end do end do end do @@ -812,7 +776,6 @@ contains G = 0._wp -#ifndef MFC_SIMULATION ! Converting the primitive variables to the conservative variables do l = 0, p do k = 0, n @@ -1019,11 +982,6 @@ contains end do end do end do -#else - if (proc_rank == 0) then - call s_mpi_abort('Conversion from primitive to ' // 'conservative variables not ' // 'implemented. Exiting.') - end if -#endif end subroutine s_convert_primitive_to_conservative_variables @@ -1070,7 +1028,6 @@ contains ! Computing the flux variables from the primitive variables, without accounting for the contribution of either viscosity or ! capillarity -#ifdef MFC_SIMULATION $:GPU_PARALLEL_LOOP(collapse=3, private='[alpha_rho_K, vel_K, alpha_K, Re_K, Y_K, rho_K, vel_K_sum, pres_K, E_K, gamma_K, & & pi_inf_K, qv_K, G_K, T_K, mix_mol_weight, R_gas]') do l = is3b, is3e @@ -1099,10 +1056,11 @@ contains pres_K = qK_prim_vf(j, k, l, eqn_idx%E) if (elasticity) then - call s_convert_species_to_mixture_variables_acc(rho_K, gamma_K, pi_inf_K, qv_K, alpha_K, alpha_rho_K, & + call s_convert_species_to_mixture_variables_kernel(rho_K, gamma_K, pi_inf_K, qv_K, alpha_K, alpha_rho_K, & & Re_K, G_K, Gs_vc) else - call s_convert_species_to_mixture_variables_acc(rho_K, gamma_K, pi_inf_K, qv_K, alpha_K, alpha_rho_K, Re_K) + call s_convert_species_to_mixture_variables_kernel(rho_K, gamma_K, pi_inf_K, qv_K, alpha_K, alpha_rho_K, & + & Re_K) end if ! Computing the energy from the pressure @@ -1169,7 +1127,6 @@ contains end do end do $:END_GPU_PARALLEL_LOOP() -#endif end subroutine s_convert_primitive_to_flux_variables @@ -1227,29 +1184,18 @@ contains !> Deallocate fluid property arrays and post-processing fields allocated during module initialization. impure subroutine s_finalize_variables_conversion_module() - ! Deallocating the density, the specific heat ratio function and the liquid stiffness function -#ifdef MFC_POST_PROCESS - deallocate (rho_sf, gamma_sf, pi_inf_sf, qv_sf) -#endif + if (allocated(rho_sf)) deallocate (rho_sf, gamma_sf, pi_inf_sf, qv_sf) -#ifdef MFC_SIMULATION @:DEALLOCATE(gammas, gs_min, pi_infs, ps_inf, cvs, qvs, qvps, Gs_vc) - if (bubbles_euler) then + if (allocated(bubrs_vc)) then @:DEALLOCATE(bubrs_vc) end if - if (viscous) then + if (allocated(Res_vc)) then @:DEALLOCATE(Res_vc) end if -#else - @:DEALLOCATE(gammas, gs_min, pi_infs, ps_inf, cvs, qvs, qvps, Gs_vc) - if (bubbles_euler) then - @:DEALLOCATE(bubrs_vc) - end if -#endif end subroutine s_finalize_variables_conversion_module -#ifndef MFC_PRE_PROCESS !> Compute the speed of sound from thermodynamic state variables, supporting multiple equation-of-state models. subroutine s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, adv, vel_sum, c_c, c, qv) @@ -1309,9 +1255,7 @@ contains end if end subroutine s_compute_speed_of_sound -#endif -#ifndef MFC_PRE_PROCESS !> Compute the fast magnetosonic wave speed from the sound speed, density, and magnetic field components. subroutine s_compute_fast_magnetosonic_speed(rho, c, B, norm, c_fast, h) @@ -1348,5 +1292,5 @@ contains c_fast = sqrt(0.5_wp*(term + sqrt(disc))) end subroutine s_compute_fast_magnetosonic_speed -#endif + end module m_variables_conversion diff --git a/src/post_process/m_global_parameters.fpp b/src/post_process/m_global_parameters.fpp index 7e1c804700..f90bbcb0d1 100644 --- a/src/post_process/m_global_parameters.fpp +++ b/src/post_process/m_global_parameters.fpp @@ -190,10 +190,6 @@ contains t_save = dflt_real t_stop = dflt_real - ! Simulation algorithm parameters (post-specific) - mixture_err = .false. - alt_soundspeed = .false. - bc_io = .false. num_bc_patches = dflt_int @@ -298,11 +294,9 @@ contains schlieren_alpha = dflt_real fd_order = dflt_int - avg_state = dflt_int ! Bubble modeling (post-specific) nb = dflt_int - sigR = dflt_real ! Output partial domain (post-specific) output_partial_domain = .false. diff --git a/src/post_process/m_start_up.fpp b/src/post_process/m_start_up.fpp index 3b294ff075..ab4181a84d 100644 --- a/src/post_process/m_start_up.fpp +++ b/src/post_process/m_start_up.fpp @@ -781,7 +781,7 @@ contains call s_initialize_mpi_common_module() end if call s_initialize_boundary_common_module() - call s_initialize_variables_conversion_module() + call s_initialize_variables_conversion_module(store_mixture_fields=.true., lagrange_beta_index=beta_idx) call s_initialize_data_input_module() call s_initialize_derived_variables_module() call s_initialize_data_output_module() diff --git a/src/pre_process/m_global_parameters.fpp b/src/pre_process/m_global_parameters.fpp index bbf2444803..748e4e2a41 100644 --- a/src/pre_process/m_global_parameters.fpp +++ b/src/pre_process/m_global_parameters.fpp @@ -184,7 +184,6 @@ contains pre_stress = .false. precision = 2 - viscous = .false. mixlayer_vel_profile = .false. mixlayer_vel_coef = 1._wp mixlayer_perturb = .false. @@ -318,7 +317,6 @@ contains Web = dflt_real nmom = 1 - sigR = dflt_real sigV = dflt_real rhoRV = 0._wp dist_type = dflt_int diff --git a/src/simulation/m_bubbles_EL.fpp b/src/simulation/m_bubbles_EL.fpp index a7057867bf..fb4ed20c39 100644 --- a/src/simulation/m_bubbles_EL.fpp +++ b/src/simulation/m_bubbles_EL.fpp @@ -659,7 +659,7 @@ contains ! Obtain liquid density and computing speed of sound from pinf call s_compute_species_fraction(q_prim_vf, cell(1), cell(2), cell(3), myalpha_rho, myalpha) - call s_convert_species_to_mixture_variables_acc(myRho, gamma, pi_inf, qv, myalpha, myalpha_rho, Re) + call s_convert_species_to_mixture_variables_kernel(myRho, gamma, pi_inf, qv, myalpha, myalpha_rho, Re) call s_compute_cson_from_pinf(q_prim_vf, myPinf, cell, myRho, gamma, pi_inf, myCson) ! Adaptive time stepping diff --git a/src/simulation/m_cbc.fpp b/src/simulation/m_cbc.fpp index 839485464c..20dac4e200 100644 --- a/src/simulation/m_cbc.fpp +++ b/src/simulation/m_cbc.fpp @@ -624,7 +624,7 @@ contains adv_local(i) = q_prim_rs${XYZ}$_vf(0, k, r, eqn_idx%E + i) end do - call s_convert_species_to_mixture_variables_acc(rho, gamma, pi_inf, qv, adv_local, alpha_rho, Re_cbc) + call s_convert_species_to_mixture_variables_kernel(rho, gamma, pi_inf, qv, adv_local, alpha_rho, Re_cbc) $:GPU_LOOP(parallelism='[seq]') do i = 1, eqn_idx%cont%end diff --git a/src/simulation/m_global_parameters.fpp b/src/simulation/m_global_parameters.fpp index 5db4b0f72e..1ef48eaf3d 100644 --- a/src/simulation/m_global_parameters.fpp +++ b/src/simulation/m_global_parameters.fpp @@ -150,16 +150,6 @@ module m_global_parameters type(int_bounds_info) :: idwbuff(1:3) $:GPU_DECLARE(create='[idwbuff]') - !> @name The number of fluids, along with their identifying indexes, respectively, for which viscous effects, e.g. the shear - !! and/or the volume Reynolds (Re) numbers, will be non-negligible. - !> @{ - integer, dimension(2) :: Re_size - integer :: Re_size_max - integer, allocatable, dimension(:,:) :: Re_idx - !> @} - - $:GPU_DECLARE(create='[Re_size, Re_size_max, Re_idx]') - !> @name Herschel-Bulkley non-Newtonian viscosity: per-fluid flags and parameter arrays. !> @{ logical :: any_non_newtonian !< .true. if any fluid is non-Newtonian @@ -178,17 +168,6 @@ module m_global_parameters $:GPU_DECLARE(create='[wa_flg]') - !> @name The coordinate direction indexes and flags (flg), respectively, for which the configurations will be determined with - !! respect to a working direction and that will be used to isolate the contributions, in that direction, in the dimensionally - !! split system of equations. - !> @{ - integer, dimension(3) :: dir_idx - real(wp), dimension(3) :: dir_flg - integer, dimension(3) :: dir_idx_tau !< used for hypoelasticity=true - !> @} - - $:GPU_DECLARE(create='[dir_idx, dir_flg, dir_idx_tau]') - integer :: buff_size !< Number of ghost cells for boundary condition storage $:GPU_DECLARE(create='[buff_size]') @@ -280,9 +259,6 @@ module m_global_parameters !> @{ !> @} - real(wp), allocatable, dimension(:) :: gammas, gs_min, pi_infs, ps_inf, cvs, qvs, qvps - $:GPU_DECLARE(create='[gammas, gs_min, pi_infs, ps_inf, cvs, qvs, qvps]') - real(wp) :: mytime !< Current simulation time real(wp) :: finaltime !< Final simulation time type(pres_field), allocatable, dimension(:) :: pb_ts @@ -378,13 +354,9 @@ contains mp_weno = .false. weno_avg = .false. weno_Re_flux = .false. - riemann_solver = dflt_int low_Mach = 0 wave_speeds = dflt_int - avg_state = dflt_int - alt_soundspeed = .false. null_weights = .false. - mixture_err = .false. precision = 2 palpha_eps = dflt_real ptgalpha_eps = dflt_real @@ -407,7 +379,6 @@ contains wenoz_q = dflt_real igr_order = dflt_int igr_pres_lim = .false. - viscous = .false. igr_iter_solver = 1 #:endif diff --git a/src/simulation/m_hyperelastic.fpp b/src/simulation/m_hyperelastic.fpp index e0017a3685..1a4d2bc518 100644 --- a/src/simulation/m_hyperelastic.fpp +++ b/src/simulation/m_hyperelastic.fpp @@ -101,7 +101,7 @@ contains call s_compute_species_fraction(q_cons_vf, j, k, l, alpha_rho_k, alpha_k) ! If in simulation, use acc mixture subroutines - call s_convert_species_to_mixture_variables_acc(rho, gamma, pi_inf, qv, alpha_k, alpha_rho_k, Re, G_local, & + call s_convert_species_to_mixture_variables_kernel(rho, gamma, pi_inf, qv, alpha_k, alpha_rho_k, Re, G_local, & & Gs_hyper) rho = max(rho, sgm_eps) G_local = max(G_local, sgm_eps) diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp index 624784ee38..dd6477c8f3 100644 --- a/src/simulation/m_ibm.fpp +++ b/src/simulation/m_ibm.fpp @@ -296,10 +296,10 @@ contains if (model_eqns /= model_eqns_4eq) then ! If in simulation, use acc mixture subroutines if (elasticity) then - call s_convert_species_to_mixture_variables_acc(rho, gamma, pi_inf, qv_K, alpha_IP, alpha_rho_IP, Re_K, & - & G_K, Gs) + call s_convert_species_to_mixture_variables_kernel(rho, gamma, pi_inf, qv_K, alpha_IP, alpha_rho_IP, & + & Re_K, G_K, Gs) else - call s_convert_species_to_mixture_variables_acc(rho, gamma, pi_inf, qv_K, alpha_IP, alpha_rho_IP, Re_K) + call s_convert_species_to_mixture_variables_kernel(rho, gamma, pi_inf, qv_K, alpha_IP, alpha_rho_IP, Re_K) end if end if diff --git a/src/simulation/m_sim_helpers.fpp b/src/simulation/m_sim_helpers.fpp index 3e380a8081..29f07c7ce6 100644 --- a/src/simulation/m_sim_helpers.fpp +++ b/src/simulation/m_sim_helpers.fpp @@ -97,9 +97,9 @@ contains call s_compute_species_fraction(q_prim_vf, j, k, l, alpha_rho, alpha) if (elasticity) then - call s_convert_species_to_mixture_variables_acc(rho, gamma, pi_inf, qv, alpha, alpha_rho, Re, G_local, Gs) + call s_convert_species_to_mixture_variables_kernel(rho, gamma, pi_inf, qv, alpha, alpha_rho, Re, G_local, Gs) else - call s_convert_species_to_mixture_variables_acc(rho, gamma, pi_inf, qv, alpha, alpha_rho, Re) + call s_convert_species_to_mixture_variables_kernel(rho, gamma, pi_inf, qv, alpha, alpha_rho, Re) end if if (igr) then diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index 97a1c76afd..88c13c1005 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -816,7 +816,7 @@ contains end if call s_initialize_mpi_common_module() call s_initialize_mpi_proxy_module() - call s_initialize_variables_conversion_module() + call s_initialize_variables_conversion_module(enforce_density_floor=.true., preserve_qbmm_number=.true.) if (grid_geometry == 3) call s_initialize_fftw_module() if (bubbles_euler) call s_initialize_bubbles_EE_module() diff --git a/toolchain/mfc/params/definitions.py b/toolchain/mfc/params/definitions.py index 8eb3002f44..93fd27263a 100644 --- a/toolchain/mfc/params/definitions.py +++ b/toolchain/mfc/params/definitions.py @@ -1160,6 +1160,11 @@ def _init_registry(): NAMELIST_VARS: dict[str, set[str]] = {} +# Some common modules are compiled into every executable even though the +# corresponding user input remains meaningful in only a subset of namelists. +# Keep declaration visibility separate from input acceptance. +DECLARATION_TARGETS: dict[str, set[str]] = {} + # Maps indexed-family base names to their Fortran dimension expression. # The generator emits `{type}, dimension({dim}) :: {name}` for each entry. # Add here whenever a new array param needs no manual Fortran declaration. @@ -1207,12 +1212,19 @@ def _nv(targets: set, *names: str) -> None: NAMELIST_VARS[n] = set(targets) +def _decl(targets: set, *names: str) -> None: + for n in names: + DECLARATION_TARGETS[n] = set(targets) + + _ALL, _PRE_SIM, _SIM_POST = {"pre", "sim", "post"}, {"pre", "sim"}, {"sim", "post"} _PRE_POST = {"pre", "post"} _SIM = {"sim"} _PRE = {"pre"} _POST = {"post"} +_decl(_ALL, "avg_state", "alt_soundspeed", "mixture_err", "sigR", "viscous", "riemann_solver") + _nv( _ALL, "m", @@ -1271,12 +1283,12 @@ def _nv(targets: set, *names: str) -> None: "t_stop", "t_save", "cfl_target", - "avg_state", "prim_vars_wrt", - "alt_soundspeed", - "mixture_err", "fd_order", "ib_state_wrt", + "avg_state", + "alt_soundspeed", + "mixture_err", ) _nv( _PRE_SIM, @@ -1286,7 +1298,7 @@ def _nv(targets: set, *names: str) -> None: "patch_ib", "pi_fac", ) -_nv(_PRE_POST, "num_fluids", "weno_order", "recon_type", "muscl_order", "mhd", "nb", "sigR", "igr", "igr_order") +_nv(_PRE_POST, "num_fluids", "weno_order", "recon_type", "muscl_order", "mhd", "nb", "igr", "igr_order", "sigR") _nv(_ALL, "reactive_burn", "rburn") _nv(_PRE_SIM, "ib_airfoil") _nv(_PRE_SIM, "stl_models", "num_stl_models") @@ -1309,7 +1321,6 @@ def _nv(targets: set, *names: str) -> None: "int_comp", "ic_eps", "ic_beta", - "riemann_solver", "wave_speeds", "low_Mach", "hyper_cleaning_speed", @@ -1374,6 +1385,7 @@ def _nv(targets: set, *names: str) -> None: "nv_uvm_out_of_core", "nv_uvm_igr_temps_on_gpu", "nv_uvm_pref_gpu", + "riemann_solver", ) _nv( _PRE, @@ -1402,7 +1414,6 @@ def _nv(targets: set, *names: str) -> None: "sigV", "dist_type", "rhoRV", - "viscous", "old_grid", "old_ic", "perturb_flow", @@ -1423,6 +1434,7 @@ def _nv(targets: set, *names: str) -> None: "simplex_params", "files_dir", "file_extension", + "viscous", ) _nv( _POST, diff --git a/toolchain/mfc/params/generators/fortran_gen.py b/toolchain/mfc/params/generators/fortran_gen.py index f302db8bb0..d5d6bc02ea 100644 --- a/toolchain/mfc/params/generators/fortran_gen.py +++ b/toolchain/mfc/params/generators/fortran_gen.py @@ -3,7 +3,7 @@ from pathlib import Path from typing import List, Tuple -from ..definitions import CASE_OPT_PARAMS, FORTRAN_ARRAY_DIMS, NAMELIST_VARS, TYPED_DECLS # noqa: F401 - triggers registry population +from ..definitions import CASE_OPT_PARAMS, DECLARATION_TARGETS, FORTRAN_ARRAY_DIMS, NAMELIST_VARS, TYPED_DECLS # noqa: F401 - triggers registry population from ..registry import REGISTRY from ..schema import ParamDef, ParamType @@ -59,6 +59,12 @@ def _vars_for_target(target: str) -> List[str]: return sorted(v for v, ts in NAMELIST_VARS.items() if target in ts) +def _decl_vars_for_target(target: str) -> List[str]: + namelist_vars = {v for v, targets in NAMELIST_VARS.items() if target in targets} + common_module_vars = {v for v, targets in DECLARATION_TARGETS.items() if target in targets} + return sorted(namelist_vars | common_module_vars) + + def _pack_namelist(vars_list: List[str], first_prefix: str, cont_prefix: str, max_line: int) -> List[str]: """Pack variable names into Fortran continuation lines; all but last end with ', &'.""" if not vars_list: @@ -111,7 +117,7 @@ def generate_decls_fpp(target: str) -> str: """Return Fortran declarations (scalars + known arrays) for a target.""" _check_target(target) lines = [_HEADER.rstrip()] - for name in _vars_for_target(target): + for name in _decl_vars_for_target(target): if not _is_simple_scalar(name): continue if target == "sim" and name in CASE_OPT_PARAMS: diff --git a/toolchain/mfc/params_tests/test_fortran_gen.py b/toolchain/mfc/params_tests/test_fortran_gen.py index 58c1f8525a..0771f0f3c7 100644 --- a/toolchain/mfc/params_tests/test_fortran_gen.py +++ b/toolchain/mfc/params_tests/test_fortran_gen.py @@ -117,6 +117,31 @@ def test_decls_dt_for_sim(): assert "real(wp) :: dt" in generate_decls_fpp("sim") +def test_common_module_decls_do_not_expand_namelists(): + from mfc.params.definitions import NAMELIST_VARS + from mfc.params.generators.fortran_gen import generate_decls_fpp, generate_namelist_fpp + + common_decls = ("avg_state", "alt_soundspeed", "mixture_err", "sigR", "riemann_solver") + for target in ("pre", "sim", "post"): + declarations = generate_decls_fpp(target) + for name in common_decls: + assert f":: {name}" in declarations + + # Viscous is a case-optimization parameter, so simulation declares it in + # generated_case_opt_decls.fpp instead of generated_decls.fpp. + assert ":: viscous" in generate_decls_fpp("pre") + assert ":: viscous" in generate_decls_fpp("post") + + assert "riemann_solver" not in generate_namelist_fpp("pre") + assert "riemann_solver" not in generate_namelist_fpp("post") + assert NAMELIST_VARS["avg_state"] == {"sim", "post"} + assert NAMELIST_VARS["alt_soundspeed"] == {"sim", "post"} + assert NAMELIST_VARS["mixture_err"] == {"sim", "post"} + assert NAMELIST_VARS["sigR"] == {"pre", "post"} + assert NAMELIST_VARS["viscous"] == {"pre", "sim"} + assert NAMELIST_VARS["riemann_solver"] == {"sim"} + + def test_decls_no_percent_vars(): from mfc.params.generators.fortran_gen import generate_decls_fpp From 9cbdd68d76c2627e524461403662bde09dfee444 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Sat, 25 Jul 2026 21:11:43 -0400 Subject: [PATCH 02/15] Make grid buffer handling stage-independent --- src/common/m_boundary_common.fpp | 142 +++++++++++++++++++------------ src/common/m_mpi_common.fpp | 49 +++++++++++ src/post_process/m_start_up.fpp | 8 +- src/simulation/m_start_up.fpp | 18 +++- 4 files changed, 159 insertions(+), 58 deletions(-) diff --git a/src/common/m_boundary_common.fpp b/src/common/m_boundary_common.fpp index f407bdbcf1..07776d1c10 100644 --- a/src/common/m_boundary_common.fpp +++ b/src/common/m_boundary_common.fpp @@ -342,70 +342,68 @@ contains !> Populate the buffers of the grid variables, which are constituted of the cell-boundary locations and cell-width !! distributions, based on the boundary conditions. - subroutine s_populate_grid_variables_buffers - -#ifdef MFC_SIMULATION - ! Simulation allocates the cell-boundary arrays with buff_size ghost layers, so the - ! ghost extrapolation extends over the full buffer. In post-process the module-level - ! offset_x/y/z (which the x_cb/y_cb/z_cb allocations are sized to) are used instead. - type(int_bounds_info) :: offset_x, offset_y, offset_z - - offset_x%beg = buff_size; offset_x%end = buff_size - offset_y%beg = buff_size; offset_y%end = buff_size - offset_z%beg = buff_size; offset_z%end = buff_size - - ! Global domain bounds + subroutine s_populate_grid_variables_buffers(x_cb_in, x_cc_in, dx_in, x_offset, y_offset, z_offset, y_cb_in, y_cc_in, dy_in, & + & z_cb_in, z_cc_in, dz_in, global_bounds) + + type(int_bounds_info), intent(in) :: x_offset, y_offset, z_offset + real(wp), intent(inout) :: x_cb_in(-1 - x_offset%beg:) + real(wp), intent(inout) :: x_cc_in(-buff_size:), dx_in(-buff_size:) + real(wp), optional, intent(inout) :: y_cb_in(-1 - y_offset%beg:), z_cb_in(-1 - z_offset%beg:) + real(wp), optional, intent(inout) :: y_cc_in(-buff_size:), dy_in(-buff_size:) + real(wp), optional, intent(inout) :: z_cc_in(-buff_size:), dz_in(-buff_size:) + type(bounds_info), optional, dimension(3), intent(out) :: global_bounds + + if (present(global_bounds)) then #ifdef MFC_MPI - call s_mpi_allreduce_min(x_cb(-1), glb_bounds(1)%beg) - call s_mpi_allreduce_max(x_cb(m), glb_bounds(1)%end) - if (n > 0) then - call s_mpi_allreduce_min(y_cb(-1), glb_bounds(2)%beg) - call s_mpi_allreduce_max(y_cb(n), glb_bounds(2)%end) - if (p > 0) then - call s_mpi_allreduce_min(z_cb(-1), glb_bounds(3)%beg) - call s_mpi_allreduce_max(z_cb(p), glb_bounds(3)%end) + call s_mpi_allreduce_min(x_cb_in(-1), global_bounds(1)%beg) + call s_mpi_allreduce_max(x_cb_in(m), global_bounds(1)%end) + if (n > 0) then + call s_mpi_allreduce_min(y_cb_in(-1), global_bounds(2)%beg) + call s_mpi_allreduce_max(y_cb_in(n), global_bounds(2)%end) + if (p > 0) then + call s_mpi_allreduce_min(z_cb_in(-1), global_bounds(3)%beg) + call s_mpi_allreduce_max(z_cb_in(p), global_bounds(3)%end) + end if end if - end if #else - glb_bounds(1)%beg = x_cb(-1); glb_bounds(1)%end = x_cb(m) - if (n > 0) then - glb_bounds(2)%beg = y_cb(-1); glb_bounds(2)%end = y_cb(n) - if (p > 0) then - glb_bounds(3)%beg = z_cb(-1); glb_bounds(3)%end = z_cb(p) + global_bounds(1)%beg = x_cb_in(-1); global_bounds(1)%end = x_cb_in(m) + if (n > 0) then + global_bounds(2)%beg = y_cb_in(-1); global_bounds(2)%end = y_cb_in(n) + if (p > 0) then + global_bounds(3)%beg = z_cb_in(-1); global_bounds(3)%end = z_cb_in(p) + end if end if - end if -#endif - $:GPU_UPDATE(device='[glb_bounds]') #endif + end if -#ifndef MFC_PRE_PROCESS - call s_populate_grid_bc_direction(1, -1, bc_x, offset_x) - call s_populate_grid_bc_direction(1, 1, bc_x, offset_x) + call s_populate_grid_bc_direction(x_cb_in, x_cc_in, dx_in, m, 1, -1, bc_x, x_offset) + call s_populate_grid_bc_direction(x_cb_in, x_cc_in, dx_in, m, 1, 1, bc_x, x_offset) if (n == 0) return #:if not MFC_CASE_OPTIMIZATION or num_dims > 1 - call s_populate_grid_bc_direction(2, -1, bc_y, offset_y) - call s_populate_grid_bc_direction(2, 1, bc_y, offset_y) + call s_populate_grid_bc_direction(y_cb_in, y_cc_in, dy_in, n, 2, -1, bc_y, y_offset) + call s_populate_grid_bc_direction(y_cb_in, y_cc_in, dy_in, n, 2, 1, bc_y, y_offset) #:endif if (p == 0) return #:if not MFC_CASE_OPTIMIZATION or num_dims > 2 - call s_populate_grid_bc_direction(3, -1, bc_z, offset_z) - call s_populate_grid_bc_direction(3, 1, bc_z, offset_z) + call s_populate_grid_bc_direction(z_cb_in, z_cc_in, dz_in, p, 3, -1, bc_z, z_offset) + call s_populate_grid_bc_direction(z_cb_in, z_cc_in, dz_in, p, 3, 1, bc_z, z_offset) #:endif -#endif end subroutine s_populate_grid_variables_buffers -#ifndef MFC_PRE_PROCESS - !> Populate grid variable buffers (cell widths and centers) for one direction and location. - subroutine s_populate_grid_bc_direction(bc_dir, bc_loc, bc_bounds, offset_dir) + !> Populate cell-boundary, cell-center, and cell-width buffers for one coordinate direction. + subroutine s_populate_grid_bc_direction(cell_boundaries, cell_centers, cell_widths, num_cells, bc_dir, bc_loc, bc_bounds, & + & offset) - integer, intent(in) :: bc_dir, bc_loc - type(int_bounds_info), intent(in) :: bc_bounds, offset_dir - integer :: bc_edge + integer, intent(in) :: num_cells, bc_dir, bc_loc + type(int_bounds_info), intent(in) :: bc_bounds, offset + real(wp), intent(inout) :: cell_boundaries(-1 - offset%beg:) + real(wp), intent(inout) :: cell_centers(-buff_size:), cell_widths(-buff_size:) + integer :: bc_edge, i, source_index if (bc_loc == -1) then bc_edge = bc_bounds%beg @@ -414,23 +412,55 @@ contains end if if (bc_edge >= 0) then - call s_mpi_sendrecv_grid_variables_buffers(bc_dir, bc_loc, offset_dir) + call s_mpi_sendrecv_grid_variable_buffer(cell_boundaries, cell_centers, cell_widths, num_cells, bc_bounds, bc_loc, & + & offset) return end if - select case (bc_edge) - case (BC_PERIODIC) - call s_grid_periodic_bc(bc_dir, bc_loc, offset_dir) - case (BC_REFLECTIVE) - call s_grid_reflective_bc(bc_dir, bc_loc, offset_dir) - case (BC_AXIS) - if (bc_dir == 2) call s_grid_axis_bc(bc_loc, offset_dir) - case default - call s_grid_ghost_cell_extrapolation_bc(bc_dir, bc_loc, offset_dir) - end select + if (bc_edge == BC_AXIS .and. (bc_dir /= 2 .or. bc_loc == 1)) return + + do i = 1, buff_size + if (bc_loc == -1) then + select case (bc_edge) + case (BC_PERIODIC) + source_index = num_cells - i + 1 + case (BC_REFLECTIVE, BC_AXIS) + source_index = i - 1 + case default + source_index = 0 + end select + cell_widths(-i) = cell_widths(source_index) + else + select case (bc_edge) + case (BC_PERIODIC) + source_index = i - 1 + case (BC_REFLECTIVE) + source_index = num_cells - i + 1 + case default + source_index = num_cells + end select + cell_widths(num_cells + i) = cell_widths(source_index) + end if + end do + + if (bc_loc == -1) then + do i = 1, offset%beg + cell_boundaries(-1 - i) = cell_boundaries(-i) - cell_widths(-i) + end do + do i = 1, buff_size + cell_centers(-i) = cell_centers(1 - i) - (cell_widths(1 - i) + cell_widths(-i))/2._wp + end do + else + do i = 1, offset%end + cell_boundaries(num_cells + i) = cell_boundaries(num_cells + i - 1) + cell_widths(num_cells + i) + end do + do i = 1, buff_size + cell_centers(num_cells + i) = cell_centers(num_cells + i - 1) + (cell_widths(num_cells + i - 1) & + & + cell_widths(num_cells + i))/2._wp + end do + end if end subroutine s_populate_grid_bc_direction -#endif !> Deallocate boundary condition buffer arrays allocated during module initialization. subroutine s_finalize_boundary_common_module() diff --git a/src/common/m_mpi_common.fpp b/src/common/m_mpi_common.fpp index 3f6e34f6c5..efe44c4fe7 100644 --- a/src/common/m_mpi_common.fpp +++ b/src/common/m_mpi_common.fpp @@ -1769,6 +1769,55 @@ contains !> The goal of this procedure is to populate the buffers of the grid variables by communicating with the neighboring processors. !! Note that only the buffers of the cell-width distributions are handled in such a way. This is because the buffers of !! cell-boundary locations may be calculated directly from those of the cell-width distributions. + subroutine s_mpi_sendrecv_grid_variable_buffer(cell_boundaries, cell_centers, cell_widths, num_cells, bc_bounds, pbc_loc, & + & offset) + + integer, intent(in) :: num_cells, pbc_loc + type(int_bounds_info), intent(in) :: bc_bounds, offset + real(wp), intent(inout) :: cell_boundaries(-1 - offset%beg:) + real(wp), intent(inout) :: cell_centers(-buff_size:) + real(wp), intent(inout) :: cell_widths(-buff_size:) + +#ifdef MFC_MPI + integer :: ierr + integer :: i + + if (pbc_loc == -1) then + if (bc_bounds%end >= 0) then + call MPI_SENDRECV(cell_widths(num_cells - buff_size + 1), buff_size, mpi_p, bc_bounds%end, 0, & + & cell_widths(-buff_size), buff_size, mpi_p, bc_bounds%beg, 0, MPI_COMM_WORLD, & + & MPI_STATUS_IGNORE, ierr) + else + call MPI_SENDRECV(cell_widths(0), buff_size, mpi_p, bc_bounds%beg, 1, cell_widths(-buff_size), buff_size, mpi_p, & + & bc_bounds%beg, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr) + end if + do i = 1, offset%beg + cell_boundaries(-1 - i) = cell_boundaries(-i) - cell_widths(-i) + end do + do i = 1, buff_size + cell_centers(-i) = cell_centers(1 - i) - (cell_widths(1 - i) + cell_widths(-i))/2._wp + end do + else + if (bc_bounds%beg >= 0) then + call MPI_SENDRECV(cell_widths(0), buff_size, mpi_p, bc_bounds%beg, 1, cell_widths(num_cells + 1), buff_size, & + & mpi_p, bc_bounds%end, 1, MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr) + else + call MPI_SENDRECV(cell_widths(num_cells - buff_size + 1), buff_size, mpi_p, bc_bounds%end, 0, & + & cell_widths(num_cells + 1), buff_size, mpi_p, bc_bounds%end, 1, MPI_COMM_WORLD, & + & MPI_STATUS_IGNORE, ierr) + end if + do i = 1, offset%end + cell_boundaries(num_cells + i) = cell_boundaries(num_cells + i - 1) + cell_widths(num_cells + i) + end do + do i = 1, buff_size + cell_centers(num_cells + i) = cell_centers(num_cells + i - 1) + (cell_widths(num_cells + i - 1) & + & + cell_widths(num_cells + i))/2._wp + end do + end if +#endif + + end subroutine s_mpi_sendrecv_grid_variable_buffer + #ifndef MFC_PRE_PROCESS subroutine s_mpi_sendrecv_grid_variables_buffers(mpi_dir, pbc_loc, offset) diff --git a/src/post_process/m_start_up.fpp b/src/post_process/m_start_up.fpp index ab4181a84d..021ea86705 100644 --- a/src/post_process/m_start_up.fpp +++ b/src/post_process/m_start_up.fpp @@ -162,7 +162,13 @@ contains if (chemistry) call s_compute_q_T_sf(q_T_sf, q_cons_vf, idwint) if (buff_size > 0) then - call s_populate_grid_variables_buffers() + if (n == 0) then + call s_populate_grid_variables_buffers(x_cb, x_cc, dx, offset_x, offset_y, offset_z) + else if (p == 0) then + call s_populate_grid_variables_buffers(x_cb, x_cc, dx, offset_x, offset_y, offset_z, y_cb, y_cc, dy) + else + call s_populate_grid_variables_buffers(x_cb, x_cc, dx, offset_x, offset_y, offset_z, y_cb, y_cc, dy, z_cb, z_cc, dz) + end if call s_populate_variables_buffers(bc_type, q_cons_vf, q_T_sf=q_T_sf) end if diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index 88c13c1005..082a471aae 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -870,7 +870,23 @@ contains call s_read_data_files(q_cons_ts(1)%vf) end if - call s_populate_grid_variables_buffers() + block + type(int_bounds_info), dimension(3) :: grid_offsets + + grid_offsets(:)%beg = buff_size + grid_offsets(:)%end = buff_size + if (n == 0) then + call s_populate_grid_variables_buffers(x_cb, x_cc, dx, grid_offsets(1), grid_offsets(2), grid_offsets(3), & + & global_bounds=glb_bounds) + else if (p == 0) then + call s_populate_grid_variables_buffers(x_cb, x_cc, dx, grid_offsets(1), grid_offsets(2), grid_offsets(3), y_cb, & + & y_cc, dy, global_bounds=glb_bounds) + else + call s_populate_grid_variables_buffers(x_cb, x_cc, dx, grid_offsets(1), grid_offsets(2), grid_offsets(3), y_cb, & + & y_cc, dy, z_cb, z_cc, dz, glb_bounds) + end if + end block + $:GPU_UPDATE(device='[glb_bounds]') if (model_eqns == model_eqns_6eq) call s_initialize_internal_energy_equations(q_cons_ts(1)%vf) if (ib) then From 83a05fdb415273cead5271d3e464bc1f06f56474 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Sat, 25 Jul 2026 21:22:30 -0400 Subject: [PATCH 03/15] Separate minimum spacing from grid width arrays --- src/common/include/2dHardcodedIC.fpp | 2 +- src/common/include/3dHardcodedIC.fpp | 4 +- src/common/m_boundary_primitives.fpp | 253 ---------------------- src/common/m_global_parameters_common.fpp | 6 + src/common/m_model.fpp | 14 +- src/common/m_mpi_common.fpp | 143 ++---------- src/pre_process/m_global_parameters.fpp | 1 - src/pre_process/m_grid.f90 | 58 ++--- src/pre_process/m_icpp_patches.fpp | 40 ++-- src/pre_process/m_start_up.fpp | 24 +- src/simulation/m_start_up.fpp | 3 + 11 files changed, 94 insertions(+), 454 deletions(-) diff --git a/src/common/include/2dHardcodedIC.fpp b/src/common/include/2dHardcodedIC.fpp index c982636779..360c49b718 100644 --- a/src/common/include/2dHardcodedIC.fpp +++ b/src/common/include/2dHardcodedIC.fpp @@ -161,7 +161,7 @@ ! Smoothening function to smooth out sharp discontinuity in the interface if (x_cc(i) <= 0.7_wp*lam) then d = x_cc(i) - lam*(0.4_wp - 0.1_wp*sin(2.0_wp*pi*(y_cc(j)/lam + 0.25_wp))) - fsm = 0.5_wp*(1.0_wp + erf(d/(ei*sqrt(dx*dy)))) + fsm = 0.5_wp*(1.0_wp + erf(d/(ei*sqrt(dx_min*dy_min)))) alpha_air = eps + (1.0_wp - 2.0_wp*eps)*fsm alpha_sf6 = 1.0_wp - alpha_air q_prim_vf(eqn_idx%cont%beg)%sf(i, j, 0) = alpha_sf6*5.04_wp diff --git a/src/common/include/3dHardcodedIC.fpp b/src/common/include/3dHardcodedIC.fpp index 3542b8eddf..76835a1aa4 100644 --- a/src/common/include/3dHardcodedIC.fpp +++ b/src/common/include/3dHardcodedIC.fpp @@ -233,7 +233,7 @@ q_prim_vf(eqn_idx%E)%sf(i, j, k) = p_th*rcut*xcut + p_am case (304) ! 3D Interface from file cartesian - alph = 0.5_wp*(1 + (1._wp - 2._wp*eps)*tanh((ih(start_idx(2) + j, start_idx(3) + k) - x_cc(i))*(0.5_wp/dx))) + alph = 0.5_wp*(1 + (1._wp - 2._wp*eps)*tanh((ih(start_idx(2) + j, start_idx(3) + k) - x_cc(i))*(0.5_wp/dx_min))) q_prim_vf(eqn_idx%adv%beg)%sf(i, j, k) = alph q_prim_vf(eqn_idx%adv%end)%sf(i, j, k) = 1._wp - alph @@ -246,7 +246,7 @@ if (surface_tension) q_prim_vf(eqn_idx%c)%sf(i, j, k) = alph case (305) ! 3D Interface from file axisymmetric - alph = 0.5_wp*(1 + (1._wp - 2._wp*eps)*tanh((ih(start_idx(2) + j, 0) - x_cc(i))*(0.01_wp/dx))) + alph = 0.5_wp*(1 + (1._wp - 2._wp*eps)*tanh((ih(start_idx(2) + j, 0) - x_cc(i))*(0.01_wp/dx_min))) q_prim_vf(eqn_idx%adv%beg)%sf(i, j, k) = alph q_prim_vf(eqn_idx%adv%end)%sf(i, j, k) = 1._wp - alph diff --git a/src/common/m_boundary_primitives.fpp b/src/common/m_boundary_primitives.fpp index 2fce625617..00a61bd4a5 100644 --- a/src/common/m_boundary_primitives.fpp +++ b/src/common/m_boundary_primitives.fpp @@ -1619,257 +1619,4 @@ contains end subroutine s_beta_reflective -#ifndef MFC_PRE_PROCESS - !> Apply periodic boundary conditions to grid variables by copying cell widths from opposite domain boundary. - subroutine s_grid_periodic_bc(bc_dir, bc_loc, offset_dir) - - integer, intent(in) :: bc_dir, bc_loc - type(int_bounds_info), intent(in) :: offset_dir - integer :: i - - if (bc_dir == 1) then - if (bc_loc == -1) then - do i = 1, buff_size - dx(-i) = dx(m - (i - 1)) - end do - do i = 1, offset_dir%beg - x_cb(-1 - i) = x_cb(-i) - dx(-i) - end do - do i = 1, buff_size - x_cc(-i) = x_cc(1 - i) - (dx(1 - i) + dx(-i))/2._wp - end do - else - do i = 1, buff_size - dx(m + i) = dx(i - 1) - end do - do i = 1, offset_dir%end - x_cb(m + i) = x_cb(m + (i - 1)) + dx(m + i) - end do - do i = 1, buff_size - x_cc(m + i) = x_cc(m + (i - 1)) + (dx(m + (i - 1)) + dx(m + i))/2._wp - end do - end if - else if (bc_dir == 2) then - if (bc_loc == -1) then - do i = 1, buff_size - dy(-i) = dy(n - (i - 1)) - end do - do i = 1, offset_dir%beg - y_cb(-1 - i) = y_cb(-i) - dy(-i) - end do - do i = 1, buff_size - y_cc(-i) = y_cc(1 - i) - (dy(1 - i) + dy(-i))/2._wp - end do - else - do i = 1, buff_size - dy(n + i) = dy(i - 1) - end do - do i = 1, offset_dir%end - y_cb(n + i) = y_cb(n + (i - 1)) + dy(n + i) - end do - do i = 1, buff_size - y_cc(n + i) = y_cc(n + (i - 1)) + (dy(n + (i - 1)) + dy(n + i))/2._wp - end do - end if - else - if (bc_loc == -1) then - do i = 1, buff_size - dz(-i) = dz(p - (i - 1)) - end do - do i = 1, offset_dir%beg - z_cb(-1 - i) = z_cb(-i) - dz(-i) - end do - do i = 1, buff_size - z_cc(-i) = z_cc(1 - i) - (dz(1 - i) + dz(-i))/2._wp - end do - else - do i = 1, buff_size - dz(p + i) = dz(i - 1) - end do - do i = 1, offset_dir%end - z_cb(p + i) = z_cb(p + (i - 1)) + dz(p + i) - end do - do i = 1, buff_size - z_cc(p + i) = z_cc(p + (i - 1)) + (dz(p + (i - 1)) + dz(p + i))/2._wp - end do - end if - end if - - end subroutine s_grid_periodic_bc - - !> Apply reflective boundary conditions to grid variables by mirroring cell widths across the boundary. - subroutine s_grid_reflective_bc(bc_dir, bc_loc, offset_dir) - - integer, intent(in) :: bc_dir, bc_loc - type(int_bounds_info), intent(in) :: offset_dir - integer :: i - - if (bc_dir == 1) then - if (bc_loc == -1) then - do i = 1, buff_size - dx(-i) = dx(i - 1) - end do - do i = 1, offset_dir%beg - x_cb(-1 - i) = x_cb(-i) - dx(-i) - end do - do i = 1, buff_size - x_cc(-i) = x_cc(1 - i) - (dx(1 - i) + dx(-i))/2._wp - end do - else - do i = 1, buff_size - dx(m + i) = dx(m - (i - 1)) - end do - do i = 1, offset_dir%end - x_cb(m + i) = x_cb(m + (i - 1)) + dx(m + i) - end do - do i = 1, buff_size - x_cc(m + i) = x_cc(m + (i - 1)) + (dx(m + (i - 1)) + dx(m + i))/2._wp - end do - end if - else if (bc_dir == 2) then - if (bc_loc == -1) then - do i = 1, buff_size - dy(-i) = dy(i - 1) - end do - do i = 1, offset_dir%beg - y_cb(-1 - i) = y_cb(-i) - dy(-i) - end do - do i = 1, buff_size - y_cc(-i) = y_cc(1 - i) - (dy(1 - i) + dy(-i))/2._wp - end do - else - do i = 1, buff_size - dy(n + i) = dy(n - (i - 1)) - end do - do i = 1, offset_dir%end - y_cb(n + i) = y_cb(n + (i - 1)) + dy(n + i) - end do - do i = 1, buff_size - y_cc(n + i) = y_cc(n + (i - 1)) + (dy(n + (i - 1)) + dy(n + i))/2._wp - end do - end if - else - if (bc_loc == -1) then - do i = 1, buff_size - dz(-i) = dz(i - 1) - end do - do i = 1, offset_dir%beg - z_cb(-1 - i) = z_cb(-i) - dz(-i) - end do - do i = 1, buff_size - z_cc(-i) = z_cc(1 - i) - (dz(1 - i) + dz(-i))/2._wp - end do - else - do i = 1, buff_size - dz(p + i) = dz(p - (i - 1)) - end do - do i = 1, offset_dir%end - z_cb(p + i) = z_cb(p + (i - 1)) + dz(p + i) - end do - do i = 1, buff_size - z_cc(p + i) = z_cc(p + (i - 1)) + (dz(p + (i - 1)) + dz(p + i))/2._wp - end do - end if - end if - - end subroutine s_grid_reflective_bc - - !> Extrapolate grid variables by copying boundary cell width into ghost cells. - subroutine s_grid_ghost_cell_extrapolation_bc(bc_dir, bc_loc, offset_dir) - - integer, intent(in) :: bc_dir, bc_loc - type(int_bounds_info), intent(in) :: offset_dir - integer :: i - - if (bc_dir == 1) then - if (bc_loc == -1) then - do i = 1, buff_size - dx(-i) = dx(0) - end do - do i = 1, offset_dir%beg - x_cb(-1 - i) = x_cb(-i) - dx(-i) - end do - do i = 1, buff_size - x_cc(-i) = x_cc(1 - i) - (dx(1 - i) + dx(-i))/2._wp - end do - else - do i = 1, buff_size - dx(m + i) = dx(m) - end do - do i = 1, offset_dir%end - x_cb(m + i) = x_cb(m + (i - 1)) + dx(m + i) - end do - do i = 1, buff_size - x_cc(m + i) = x_cc(m + (i - 1)) + (dx(m + (i - 1)) + dx(m + i))/2._wp - end do - end if - else if (bc_dir == 2) then - if (bc_loc == -1) then - do i = 1, buff_size - dy(-i) = dy(0) - end do - do i = 1, offset_dir%beg - y_cb(-1 - i) = y_cb(-i) - dy(-i) - end do - do i = 1, buff_size - y_cc(-i) = y_cc(1 - i) - (dy(1 - i) + dy(-i))/2._wp - end do - else - do i = 1, buff_size - dy(n + i) = dy(n) - end do - do i = 1, offset_dir%end - y_cb(n + i) = y_cb(n + (i - 1)) + dy(n + i) - end do - do i = 1, buff_size - y_cc(n + i) = y_cc(n + (i - 1)) + (dy(n + (i - 1)) + dy(n + i))/2._wp - end do - end if - else - if (bc_loc == -1) then - do i = 1, buff_size - dz(-i) = dz(0) - end do - do i = 1, offset_dir%beg - z_cb(-1 - i) = z_cb(-i) - dz(-i) - end do - do i = 1, buff_size - z_cc(-i) = z_cc(1 - i) - (dz(1 - i) + dz(-i))/2._wp - end do - else - do i = 1, buff_size - dz(p + i) = dz(p) - end do - do i = 1, offset_dir%end - z_cb(p + i) = z_cb(p + (i - 1)) + dz(p + i) - end do - do i = 1, buff_size - z_cc(p + i) = z_cc(p + (i - 1)) + (dz(p + (i - 1)) + dz(p + i))/2._wp - end do - end if - end if - - end subroutine s_grid_ghost_cell_extrapolation_bc - - !> Apply axis boundary conditions to grid variables for cylindrical coordinates. - subroutine s_grid_axis_bc(bc_loc, offset_dir) - - integer, intent(in) :: bc_loc - type(int_bounds_info), intent(in) :: offset_dir - integer :: i - - if (bc_loc == -1) then - do i = 1, buff_size - dy(-i) = dy(i - 1) - end do - do i = 1, offset_dir%beg - y_cb(-1 - i) = y_cb(-i) - dy(-i) - end do - do i = 1, buff_size - y_cc(-i) = y_cc(1 - i) - (dy(1 - i) + dy(-i))/2._wp - end do - end if - - end subroutine s_grid_axis_bc -#endif end module m_boundary_primitives diff --git a/src/common/m_global_parameters_common.fpp b/src/common/m_global_parameters_common.fpp index dfcb0e1b3b..fc396b3b3a 100644 --- a/src/common/m_global_parameters_common.fpp +++ b/src/common/m_global_parameters_common.fpp @@ -87,6 +87,9 @@ module m_global_parameters_common $:GPU_DECLARE(create='[dir_idx, dir_flg, dir_idx_tau]') !> @} + !> Minimum cell widths. These are distinct from the per-cell width arrays named dx, dy, and dz in simulation and post-process. + real(wp) :: dx_min, dy_min, dz_min + #ifdef MFC_SIMULATION $:GPU_DECLARE(create='[sys_size, eqn_idx, b_size, tensor_size]') $:GPU_DECLARE(create='[shear_num, shear_indices, shear_BC_flip_num, shear_BC_flip_indices]') @@ -494,6 +497,9 @@ contains alt_soundspeed = .false. mixture_err = .false. sigR = dflt_real + dx_min = dflt_real + dy_min = dflt_real + dz_min = dflt_real end subroutine s_assign_common_defaults diff --git a/src/common/m_model.fpp b/src/common/m_model.fpp index 76770a4a62..0183aab55f 100644 --- a/src/common/m_model.fpp +++ b/src/common/m_model.fpp @@ -858,7 +858,6 @@ contains ! Variables for IBM+STL integer :: boundary_vertex_count, boundary_edge_count !< Boundary vertex real(wp), allocatable, dimension(:,:,:) :: boundary_v !< Boundary vertex buffer - real(wp) :: dx_local, dy_local, dz_local !< Levelset distance buffer integer :: i, j, k !< Generic loop iterators integer :: stl_id type(t_bbox) :: bbox, bbox_old @@ -868,13 +867,6 @@ contains real(wp) :: grid_mm(1:3,1:2) real(wp), dimension(1:4,1:4) :: transform, transform_n -#ifdef MFC_SIMULATION - dx_local = minval(dx); dy_local = minval(dy) - if (p /= 0) dz_local = minval(dz) -#else - dx_local = dx; dy_local = dy - if (p /= 0) dz_local = dz -#endif if (num_stl_models == 0) return @:ALLOCATE(stl_bounding_boxes(num_stl_models,1:3,1:3)) @@ -924,11 +916,11 @@ contains write (*, "(A, 3(2X, F20.10))") " > Cen:", (bbox%min(1:3) + bbox%max(1:3))/2._wp write (*, "(A, 3(2X, F20.10))") " > Max:", bbox%max(1:3) - grid_mm(1,:) = (/minval(x_cc(0:m)) - 0.5_wp*dx_local, maxval(x_cc(0:m)) + 0.5_wp*dx_local/) - grid_mm(2,:) = (/minval(y_cc(0:n)) - 0.5_wp*dy_local, maxval(y_cc(0:n)) + 0.5_wp*dy_local/) + grid_mm(1,:) = (/minval(x_cc(0:m)) - 0.5_wp*dx_min, maxval(x_cc(0:m)) + 0.5_wp*dx_min/) + grid_mm(2,:) = (/minval(y_cc(0:n)) - 0.5_wp*dy_min, maxval(y_cc(0:n)) + 0.5_wp*dy_min/) if (p > 0) then - grid_mm(3,:) = (/minval(z_cc(0:p)) - 0.5_wp*dz_local, maxval(z_cc(0:p)) + 0.5_wp*dz_local/) + grid_mm(3,:) = (/minval(z_cc(0:p)) - 0.5_wp*dz_min, maxval(z_cc(0:p)) + 0.5_wp*dz_min/) else grid_mm(3,:) = (/0._wp, 0._wp/) end if diff --git a/src/common/m_mpi_common.fpp b/src/common/m_mpi_common.fpp index efe44c4fe7..b239183e99 100644 --- a/src/common/m_mpi_common.fpp +++ b/src/common/m_mpi_common.fpp @@ -1537,15 +1537,15 @@ contains else #ifdef MFC_PRE_PROCESS if (old_grid .neqv. .true.) then - dz = (z_domain%end - z_domain%beg)/real(p_glb + 1, wp) + dz_min = (z_domain%end - z_domain%beg)/real(p_glb + 1, wp) if (proc_coords(3) < rem_cells) then - z_domain%beg = z_domain%beg + dz*real((p + 1)*proc_coords(3)) - z_domain%end = z_domain%end - dz*real((p + 1)*(num_procs_z - proc_coords(3) - 1) - (num_procs_z & - & - rem_cells)) + z_domain%beg = z_domain%beg + dz_min*real((p + 1)*proc_coords(3)) + z_domain%end = z_domain%end - dz_min*real((p + 1)*(num_procs_z - proc_coords(3) - 1) - (num_procs_z & + & - rem_cells)) else - z_domain%beg = z_domain%beg + dz*real((p + 1)*proc_coords(3) + rem_cells) - z_domain%end = z_domain%end - dz*real((p + 1)*(num_procs_z - proc_coords(3) - 1)) + z_domain%beg = z_domain%beg + dz_min*real((p + 1)*proc_coords(3) + rem_cells) + z_domain%end = z_domain%end - dz_min*real((p + 1)*(num_procs_z - proc_coords(3) - 1)) end if end if #endif @@ -1651,15 +1651,15 @@ contains else #ifdef MFC_PRE_PROCESS if (old_grid .neqv. .true.) then - dy = (y_domain%end - y_domain%beg)/real(n_glb + 1, wp) + dy_min = (y_domain%end - y_domain%beg)/real(n_glb + 1, wp) if (proc_coords(2) < rem_cells) then - y_domain%beg = y_domain%beg + dy*real((n + 1)*proc_coords(2)) - y_domain%end = y_domain%end - dy*real((n + 1)*(num_procs_y - proc_coords(2) - 1) - (num_procs_y & - & - rem_cells)) + y_domain%beg = y_domain%beg + dy_min*real((n + 1)*proc_coords(2)) + y_domain%end = y_domain%end - dy_min*real((n + 1)*(num_procs_y - proc_coords(2) - 1) - (num_procs_y & + & - rem_cells)) else - y_domain%beg = y_domain%beg + dy*real((n + 1)*proc_coords(2) + rem_cells) - y_domain%end = y_domain%end - dy*real((n + 1)*(num_procs_y - proc_coords(2) - 1)) + y_domain%beg = y_domain%beg + dy_min*real((n + 1)*proc_coords(2) + rem_cells) + y_domain%end = y_domain%end - dy_min*real((n + 1)*(num_procs_y - proc_coords(2) - 1)) end if end if #endif @@ -1736,14 +1736,15 @@ contains else #ifdef MFC_PRE_PROCESS if (old_grid .neqv. .true.) then - dx = (x_domain%end - x_domain%beg)/real(m_glb + 1, wp) + dx_min = (x_domain%end - x_domain%beg)/real(m_glb + 1, wp) if (proc_coords(1) < rem_cells) then - x_domain%beg = x_domain%beg + dx*real((m + 1)*proc_coords(1)) - x_domain%end = x_domain%end - dx*real((m + 1)*(num_procs_x - proc_coords(1) - 1) - (num_procs_x - rem_cells)) + x_domain%beg = x_domain%beg + dx_min*real((m + 1)*proc_coords(1)) + x_domain%end = x_domain%end - dx_min*real((m + 1)*(num_procs_x - proc_coords(1) - 1) - (num_procs_x & + & - rem_cells)) else - x_domain%beg = x_domain%beg + dx*real((m + 1)*proc_coords(1) + rem_cells) - x_domain%end = x_domain%end - dx*real((m + 1)*(num_procs_x - proc_coords(1) - 1)) + x_domain%beg = x_domain%beg + dx_min*real((m + 1)*proc_coords(1) + rem_cells) + x_domain%end = x_domain%end - dx_min*real((m + 1)*(num_procs_x - proc_coords(1) - 1)) end if end if #endif @@ -1818,113 +1819,6 @@ contains end subroutine s_mpi_sendrecv_grid_variable_buffer -#ifndef MFC_PRE_PROCESS - subroutine s_mpi_sendrecv_grid_variables_buffers(mpi_dir, pbc_loc, offset) - - integer, intent(in) :: mpi_dir - integer, intent(in) :: pbc_loc - !> Ghost layers for cell-boundary arrays (buff_size in simulation, module offset_* in post-process) - type(int_bounds_info), intent(in) :: offset - -#ifdef MFC_MPI - integer :: ierr !< Generic flag used to identify and report MPI errors - integer :: i - - if (mpi_dir == 1) then - if (pbc_loc == -1) then ! PBC at the beginning - if (bc_x%end >= 0) then ! PBC at the beginning and end - call MPI_SENDRECV(dx(m - buff_size + 1), buff_size, mpi_p, bc_x%end, 0, dx(-buff_size), buff_size, mpi_p, & - & bc_x%beg, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr) - else ! PBC at the beginning only - call MPI_SENDRECV(dx(0), buff_size, mpi_p, bc_x%beg, 1, dx(-buff_size), buff_size, mpi_p, bc_x%beg, 0, & - & MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr) - end if - do i = 1, offset%beg - x_cb(-1 - i) = x_cb(-i) - dx(-i) - end do - do i = 1, buff_size - x_cc(-i) = x_cc(1 - i) - (dx(1 - i) + dx(-i))/2._wp - end do - else ! PBC at the end - if (bc_x%beg >= 0) then ! PBC at the end and beginning - call MPI_SENDRECV(dx(0), buff_size, mpi_p, bc_x%beg, 1, dx(m + 1), buff_size, mpi_p, bc_x%end, 1, & - & MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr) - else ! PBC at the end only - call MPI_SENDRECV(dx(m - buff_size + 1), buff_size, mpi_p, bc_x%end, 0, dx(m + 1), buff_size, mpi_p, & - & bc_x%end, 1, MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr) - end if - do i = 1, offset%end - x_cb(m + i) = x_cb(m + (i - 1)) + dx(m + i) - end do - do i = 1, buff_size - x_cc(m + i) = x_cc(m + (i - 1)) + (dx(m + (i - 1)) + dx(m + i))/2._wp - end do - end if - else if (mpi_dir == 2) then - if (pbc_loc == -1) then ! PBC at the beginning - if (bc_y%end >= 0) then ! PBC at the beginning and end - call MPI_SENDRECV(dy(n - buff_size + 1), buff_size, mpi_p, bc_y%end, 0, dy(-buff_size), buff_size, mpi_p, & - & bc_y%beg, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr) - else ! PBC at the beginning only - call MPI_SENDRECV(dy(0), buff_size, mpi_p, bc_y%beg, 1, dy(-buff_size), buff_size, mpi_p, bc_y%beg, 0, & - & MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr) - end if - do i = 1, offset%beg - y_cb(-1 - i) = y_cb(-i) - dy(-i) - end do - do i = 1, buff_size - y_cc(-i) = y_cc(1 - i) - (dy(1 - i) + dy(-i))/2._wp - end do - else ! PBC at the end - if (bc_y%beg >= 0) then ! PBC at the end and beginning - call MPI_SENDRECV(dy(0), buff_size, mpi_p, bc_y%beg, 1, dy(n + 1), buff_size, mpi_p, bc_y%end, 1, & - & MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr) - else ! PBC at the end only - call MPI_SENDRECV(dy(n - buff_size + 1), buff_size, mpi_p, bc_y%end, 0, dy(n + 1), buff_size, mpi_p, & - & bc_y%end, 1, MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr) - end if - do i = 1, offset%end - y_cb(n + i) = y_cb(n + (i - 1)) + dy(n + i) - end do - do i = 1, buff_size - y_cc(n + i) = y_cc(n + (i - 1)) + (dy(n + (i - 1)) + dy(n + i))/2._wp - end do - end if - else - if (pbc_loc == -1) then ! PBC at the beginning - if (bc_z%end >= 0) then ! PBC at the beginning and end - call MPI_SENDRECV(dz(p - buff_size + 1), buff_size, mpi_p, bc_z%end, 0, dz(-buff_size), buff_size, mpi_p, & - & bc_z%beg, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr) - else ! PBC at the beginning only - call MPI_SENDRECV(dz(0), buff_size, mpi_p, bc_z%beg, 1, dz(-buff_size), buff_size, mpi_p, bc_z%beg, 0, & - & MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr) - end if - do i = 1, offset%beg - z_cb(-1 - i) = z_cb(-i) - dz(-i) - end do - do i = 1, buff_size - z_cc(-i) = z_cc(1 - i) - (dz(1 - i) + dz(-i))/2._wp - end do - else ! PBC at the end - if (bc_z%beg >= 0) then ! PBC at the end and beginning - call MPI_SENDRECV(dz(0), buff_size, mpi_p, bc_z%beg, 1, dz(p + 1), buff_size, mpi_p, bc_z%end, 1, & - & MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr) - else ! PBC at the end only - call MPI_SENDRECV(dz(p - buff_size + 1), buff_size, mpi_p, bc_z%end, 0, dz(p + 1), buff_size, mpi_p, & - & bc_z%end, 1, MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr) - end if - do i = 1, offset%end - z_cb(p + i) = z_cb(p + (i - 1)) + dz(p + i) - end do - do i = 1, buff_size - z_cc(p + i) = z_cc(p + (i - 1)) + (dz(p + (i - 1)) + dz(p + i))/2._wp - end do - end if - end if -#endif - - end subroutine s_mpi_sendrecv_grid_variables_buffers - !> Populate the local cell-boundary, cell-center, and cell-width arrays in one direction directly from the global cell-boundary !! array. This guarantees that every rank sees bitwise-identical values at any shared physical cell or boundary subroutine s_apply_grid_from_global_dim(x_cb_glb, m_dim_glb, m_dim, sidx, bc_beg, bc_end, cb_lo, cb_hi, cw_lo, cw_hi, & @@ -2003,7 +1897,6 @@ contains end do end subroutine s_apply_grid_from_global_dim -#endif !> Module deallocation and/or disassociation procedures impure subroutine s_finalize_mpi_common_module diff --git a/src/pre_process/m_global_parameters.fpp b/src/pre_process/m_global_parameters.fpp index 748e4e2a41..44dd397f90 100644 --- a/src/pre_process/m_global_parameters.fpp +++ b/src/pre_process/m_global_parameters.fpp @@ -36,7 +36,6 @@ module m_global_parameters real(wp), allocatable, dimension(:) :: x_cc, y_cc, z_cc !> Locations of cell-boundaries (cb) in x-, y- and z-directions, respectively real(wp), allocatable, dimension(:) :: x_cb, y_cb, z_cb - real(wp) :: dx, dy, dz !< Minimum cell-widths in the x-, y- and z-coordinate directions type(bounds_info) :: x_domain, y_domain, z_domain !< Locations of the domain bounds in the x-, y- and z-coordinate directions !> Global (pre-decomposition) domain bounds, needed by s_generate_serial_grid to stretch the grid using the full domain length !! rather than a local processor's sub-domain length diff --git a/src/pre_process/m_grid.f90 b/src/pre_process/m_grid.f90 index ea81c1f1ec..aeb1d2d93e 100644 --- a/src/pre_process/m_grid.f90 +++ b/src/pre_process/m_grid.f90 @@ -36,13 +36,13 @@ impure subroutine s_generate_serial_grid ! Generic loop iterator integer :: i, j !< generic loop operators real(wp) :: length !< domain lengths - ! Uniform grid: dx = (x_end - x_beg) / (m + 1) + ! Uniform grid: dx_min = (x_end - x_beg) / (m + 1) - dx = (x_domain%end - x_domain%beg)/real(m + 1, wp) + dx_min = (x_domain%end - x_domain%beg)/real(m + 1, wp) do i = 0, m - x_cc(i) = x_domain%beg + 5.e-1_wp*dx*real(2*i + 1, wp) - x_cb(i - 1) = x_domain%beg + dx*real(i, wp) + x_cc(i) = x_domain%beg + 5.e-1_wp*dx_min*real(2*i + 1, wp) + x_cb(i - 1) = x_domain%beg + dx_min*real(i, wp) end do x_cb(m) = x_domain%end @@ -64,9 +64,9 @@ impure subroutine s_generate_serial_grid x_cc(0:m) = (x_cb(0:m) + x_cb(-1:m - 1))/2._wp - dx = minval(x_cb(0:m) - x_cb(-1:m - 1)) + dx_min = minval(x_cb(0:m) - x_cb(-1:m - 1)) print *, 'Stretched grid: min/max x grid: ', minval(x_cc(:)), maxval(x_cc(:)) - if (num_procs > 1) call s_mpi_reduce_min(dx) + if (num_procs > 1) call s_mpi_reduce_min(dx_min) end if ! Grid Generation in the y-direction @@ -74,21 +74,21 @@ impure subroutine s_generate_serial_grid ! Axisymmetric cylindrical grid (r, z): half-cell offset at r=0 axis if (grid_geometry == 2 .and. f_approx_equal(y_domain%beg, 0.0_wp)) then - dy = (y_domain%end - y_domain%beg)/real(2*n + 1, wp) + dy_min = (y_domain%end - y_domain%beg)/real(2*n + 1, wp) - y_cc(0) = y_domain%beg + 5.e-1_wp*dy + y_cc(0) = y_domain%beg + 5.e-1_wp*dy_min y_cb(-1) = y_domain%beg do i = 1, n - y_cc(i) = y_domain%beg + 2._wp*dy*real(i, wp) - y_cb(i - 1) = y_domain%beg + dy*real(2*i - 1, wp) + y_cc(i) = y_domain%beg + 2._wp*dy_min*real(i, wp) + y_cb(i - 1) = y_domain%beg + dy_min*real(2*i - 1, wp) end do else - dy = (y_domain%end - y_domain%beg)/real(n + 1, wp) + dy_min = (y_domain%end - y_domain%beg)/real(n + 1, wp) do i = 0, n - y_cc(i) = y_domain%beg + 5.e-1_wp*dy*real(2*i + 1, wp) - y_cb(i - 1) = y_domain%beg + dy*real(i, wp) + y_cc(i) = y_domain%beg + 5.e-1_wp*dy_min*real(2*i + 1, wp) + y_cb(i - 1) = y_domain%beg + dy_min*real(i, wp) end do end if @@ -111,19 +111,19 @@ impure subroutine s_generate_serial_grid y_cb = y_cb*length y_cc(0:n) = (y_cb(0:n) + y_cb(-1:n - 1))/2._wp - dy = minval(y_cb(0:n) - y_cb(-1:n - 1)) + dy_min = minval(y_cb(0:n) - y_cb(-1:n - 1)) - if (num_procs > 1) call s_mpi_reduce_min(dy) + if (num_procs > 1) call s_mpi_reduce_min(dy_min) end if ! Grid Generation in the z-direction if (p == 0) return - dz = (z_domain%end - z_domain%beg)/real(p + 1, wp) + dz_min = (z_domain%end - z_domain%beg)/real(p + 1, wp) do i = 0, p - z_cc(i) = z_domain%beg + 5.e-1_wp*dz*real(2*i + 1, wp) - z_cb(i - 1) = z_domain%beg + dz*real(i, wp) + z_cc(i) = z_domain%beg + 5.e-1_wp*dz_min*real(2*i + 1, wp) + z_cb(i - 1) = z_domain%beg + dz_min*real(i, wp) end do z_cb(p) = z_domain%end @@ -145,9 +145,9 @@ impure subroutine s_generate_serial_grid z_cb = z_cb*length z_cc(0:p) = (z_cb(0:p) + z_cb(-1:p - 1))/2._wp - dz = minval(z_cb(0:p) - z_cb(-1:p - 1)) + dz_min = minval(z_cb(0:p) - z_cb(-1:p - 1)) - if (num_procs > 1) call s_mpi_reduce_min(dz) + if (num_procs > 1) call s_mpi_reduce_min(dz_min) end if end subroutine s_generate_serial_grid @@ -168,10 +168,10 @@ impure subroutine s_generate_parallel_grid allocate (y_cb_glb(-1:n_glb)) allocate (z_cb_glb(-1:p_glb)) - ! Uniform grid: dx = (x_end - x_beg) / (m_glb + 1) - dx = (x_domain%end - x_domain%beg)/real(m_glb + 1, wp) + ! Uniform grid: dx_min = (x_end - x_beg) / (m_glb + 1) + dx_min = (x_domain%end - x_domain%beg)/real(m_glb + 1, wp) do i = 0, m_glb - x_cb_glb(i - 1) = x_domain%beg + dx*real(i, wp) + x_cb_glb(i - 1) = x_domain%beg + dx_min*real(i, wp) end do x_cb_glb(m_glb) = x_domain%end ! Hyperbolic tangent grid stretching in x-direction (parallel version) @@ -197,15 +197,15 @@ impure subroutine s_generate_parallel_grid if (n_glb > 0) then ! Axisymmetric cylindrical grid (r, z): half-cell offset at r=0 axis if (grid_geometry == 2 .and. f_approx_equal(y_domain%beg, 0.0_wp)) then - dy = (y_domain%end - y_domain%beg)/real(2*n_glb + 1, wp) + dy_min = (y_domain%end - y_domain%beg)/real(2*n_glb + 1, wp) y_cb_glb(-1) = y_domain%beg do i = 1, n_glb - y_cb_glb(i - 1) = y_domain%beg + dy*real(2*i - 1, wp) + y_cb_glb(i - 1) = y_domain%beg + dy_min*real(2*i - 1, wp) end do else - dy = (y_domain%end - y_domain%beg)/real(n_glb + 1, wp) + dy_min = (y_domain%end - y_domain%beg)/real(n_glb + 1, wp) do i = 0, n_glb - y_cb_glb(i - 1) = y_domain%beg + dy*real(i, wp) + y_cb_glb(i - 1) = y_domain%beg + dy_min*real(i, wp) end do end if y_cb_glb(n_glb) = y_domain%end @@ -229,9 +229,9 @@ impure subroutine s_generate_parallel_grid ! Grid generation in the z-direction if (p_glb > 0) then - dz = (z_domain%end - z_domain%beg)/real(p_glb + 1, wp) + dz_min = (z_domain%end - z_domain%beg)/real(p_glb + 1, wp) do i = 0, p_glb - z_cb_glb(i - 1) = z_domain%beg + dz*real(i, wp) + z_cb_glb(i - 1) = z_domain%beg + dz_min*real(i, wp) end do z_cb_glb(p_glb) = z_domain%end if (stretch_z) then diff --git a/src/pre_process/m_icpp_patches.fpp b/src/pre_process/m_icpp_patches.fpp index d80949037a..286addf567 100644 --- a/src/pre_process/m_icpp_patches.fpp +++ b/src/pre_process/m_icpp_patches.fpp @@ -320,8 +320,8 @@ contains do i = 0, m if (patch_icpp(patch_id)%smoothen) then ! Smooth Heaviside via hyperbolic tangent; smooth_coeff controls interface sharpness - eta = tanh(smooth_coeff/min(dx, & - & dy)*(sqrt((x_cc(i) - x_centroid)**2 + (y_cc(j) - y_centroid)**2) - radius))*(-0.5_wp) + 0.5_wp + eta = tanh(smooth_coeff/min(dx_min, & + & dy_min)*(sqrt((x_cc(i) - x_centroid)**2 + (y_cc(j) - y_centroid)**2) - radius))*(-0.5_wp) + 0.5_wp end if if ((f_is_inside_cylinder(x_cc(i) - x_centroid, y_cc(j) - y_centroid, 0._wp, radius, & @@ -493,8 +493,8 @@ contains do j = 0, n do i = 0, m if (patch_icpp(patch_id)%smoothen) then - eta = tanh(smooth_coeff/min(dx, & - & dy)*(sqrt(((x_cc(i) - x_centroid)/a)**2 + ((y_cc(j) - y_centroid)/b)**2) - 1._wp))*(-0.5_wp) & + eta = tanh(smooth_coeff/min(dx_min, & + & dy_min)*(sqrt(((x_cc(i) - x_centroid)/a)**2 + ((y_cc(j) - y_centroid)/b)**2) - 1._wp))*(-0.5_wp) & & + 0.5_wp end if @@ -563,8 +563,8 @@ contains end if if (patch_icpp(patch_id)%smoothen) then - eta = tanh(smooth_coeff/min(dx, dy, & - & dz)*(sqrt(((x_cc(i) - x_centroid)/a)**2 + ((cart_y - y_centroid)/b)**2 + ((cart_z & + eta = tanh(smooth_coeff/min(dx_min, dy_min, & + & dz_min)*(sqrt(((x_cc(i) - x_centroid)/a)**2 + ((cart_y - y_centroid)/b)**2 + ((cart_z & & - z_centroid)/c)**2) - 1._wp))*(-0.5_wp) + 0.5_wp end if @@ -695,7 +695,7 @@ contains do j = 0, n do i = 0, m if (patch_icpp(patch_id)%smoothen) then - eta = 5.e-1_wp + 5.e-1_wp*tanh(smooth_coeff/min(dx, dy)*(a*x_cc(i) + b*y_cc(j) + c)/sqrt(a**2 + b**2)) + eta = 5.e-1_wp + 5.e-1_wp*tanh(smooth_coeff/min(dx_min, dy_min)*(a*x_cc(i) + b*y_cc(j) + c)/sqrt(a**2 + b**2)) end if if ((a*x_cc(i) + b*y_cc(j) + c >= 0._wp .and. patch_icpp(patch_id)%alter_patch(patch_id_fp(i, j, & @@ -885,7 +885,7 @@ contains end if end if if (patch_icpp(patch_id)%smoothen) then - eta = 0.5_wp + 0.5_wp*tanh(smooth_coeff/min(dx, dy)*(R_boundary - r)) + eta = 0.5_wp + 0.5_wp*tanh(smooth_coeff/min(dx_min, dy_min)*(R_boundary - r)) end if if ((r <= R_boundary .and. patch_icpp(patch_id)%alter_patch(patch_id_fp(i, j, 0))) .or. patch_id_fp(i, j, & & 0) == smooth_patch_id) then @@ -947,7 +947,7 @@ contains end do end do if (patch_icpp(patch_id)%smoothen) then - eta_local = 0.5_wp + 0.5_wp*tanh(smooth_coeff/min(dx, dy, dz)*(R_surface - r)) + eta_local = 0.5_wp + 0.5_wp*tanh(smooth_coeff/min(dx_min, dy_min, dz_min)*(R_surface - r)) end if if ((r <= R_surface .and. patch_icpp(patch_id)%alter_patch(patch_id_fp(i, j, k))) .or. patch_id_fp(i, j, & & k) == smooth_patch_id) then @@ -1005,9 +1005,9 @@ contains end if if (patch_icpp(patch_id)%smoothen) then - eta = tanh(smooth_coeff/min(dx, dy, & - & dz)*(sqrt((x_cc(i) - x_centroid)**2 + (cart_y - y_centroid)**2 + (cart_z - z_centroid)**2) & - & - radius))*(-0.5_wp) + 0.5_wp + eta = tanh(smooth_coeff/min(dx_min, dy_min, & + & dz_min)*(sqrt((x_cc(i) - x_centroid)**2 + (cart_y - y_centroid)**2 + (cart_z - z_centroid) & + & **2) - radius))*(-0.5_wp) + 0.5_wp end if if ((f_is_inside_sphere(x_cc(i) - x_centroid, cart_y - y_centroid, cart_z - z_centroid, & @@ -1153,16 +1153,16 @@ contains if (patch_icpp(patch_id)%smoothen) then if (.not. f_is_default(length_x)) then - eta = tanh(smooth_coeff/min(dy, & - & dz)*(sqrt((cart_y - y_centroid)**2 + (cart_z - z_centroid)**2) - radius))*(-0.5_wp) & + eta = tanh(smooth_coeff/min(dy_min, & + & dz_min)*(sqrt((cart_y - y_centroid)**2 + (cart_z - z_centroid)**2) - radius))*(-0.5_wp) & & + 0.5_wp else if (.not. f_is_default(length_y)) then - eta = tanh(smooth_coeff/min(dx, & - & dz)*(sqrt((x_cc(i) - x_centroid)**2 + (cart_z - z_centroid)**2) - radius))*(-0.5_wp) & + eta = tanh(smooth_coeff/min(dx_min, & + & dz_min)*(sqrt((x_cc(i) - x_centroid)**2 + (cart_z - z_centroid)**2) - radius))*(-0.5_wp) & & + 0.5_wp else - eta = tanh(smooth_coeff/min(dx, & - & dy)*(sqrt((x_cc(i) - x_centroid)**2 + (cart_y - y_centroid)**2) - radius))*(-0.5_wp) & + eta = tanh(smooth_coeff/min(dx_min, & + & dy_min)*(sqrt((x_cc(i) - x_centroid)**2 + (cart_y - y_centroid)**2) - radius))*(-0.5_wp) & & + 0.5_wp end if end if @@ -1240,8 +1240,8 @@ contains end if if (patch_icpp(patch_id)%smoothen) then - eta = 5.e-1_wp + 5.e-1_wp*tanh(smooth_coeff/min(dx, dy, & - & dz)*(a*x_cc(i) + b*cart_y + c*cart_z + d)/sqrt(a**2 + b**2 + c**2)) + eta = 5.e-1_wp + 5.e-1_wp*tanh(smooth_coeff/min(dx_min, dy_min, & + & dz_min)*(a*x_cc(i) + b*cart_y + c*cart_z + d)/sqrt(a**2 + b**2 + c**2)) end if if ((a*x_cc(i) + b*cart_y + c*cart_z + d >= 0._wp .and. patch_icpp(patch_id)%alter_patch(patch_id_fp(i, j, & diff --git a/src/pre_process/m_start_up.fpp b/src/pre_process/m_start_up.fpp index a4fef6d5b6..cbf8dd8ace 100644 --- a/src/pre_process/m_start_up.fpp +++ b/src/pre_process/m_start_up.fpp @@ -171,8 +171,8 @@ contains x_cc(0:m) = (x_cb(0:m) + x_cb(-1:(m - 1)))/2._wp - dx = minval(x_cb(0:m) - x_cb(-1:m - 1)) - if (num_procs > 1) call s_mpi_reduce_min(dx) + dx_min = minval(x_cb(0:m) - x_cb(-1:m - 1)) + if (num_procs > 1) call s_mpi_reduce_min(dx_min) x_domain%beg = x_cb(-1) x_domain%end = x_cb(m) @@ -191,8 +191,8 @@ contains y_cc(0:n) = (y_cb(0:n) + y_cb(-1:(n - 1)))/2._wp - dy = minval(y_cb(0:n) - y_cb(-1:n - 1)) - if (num_procs > 1) call s_mpi_reduce_min(dy) + dy_min = minval(y_cb(0:n) - y_cb(-1:n - 1)) + if (num_procs > 1) call s_mpi_reduce_min(dy_min) y_domain%beg = y_cb(-1) y_domain%end = y_cb(n) @@ -211,8 +211,8 @@ contains z_cc(0:p) = (z_cb(0:p) + z_cb(-1:(p - 1)))/2._wp - dz = minval(z_cb(0:p) - z_cb(-1:p - 1)) - if (num_procs > 1) call s_mpi_reduce_min(dz) + dz_min = minval(z_cb(0:p) - z_cb(-1:p - 1)) + if (num_procs > 1) call s_mpi_reduce_min(dz_min) z_domain%beg = z_cb(-1) z_domain%end = z_cb(p) @@ -344,8 +344,8 @@ contains x_cb(-1:m) = x_cb_glb((start_idx(1) - 1):(start_idx(1) + m)) x_cc(0:m) = (x_cb(0:m) + x_cb(-1:(m - 1)))/2._wp - dx = minval(x_cb(0:m) - x_cb(-1:(m - 1))) - if (num_procs > 1) call s_mpi_reduce_min(dx) + dx_min = minval(x_cb(0:m) - x_cb(-1:(m - 1))) + if (num_procs > 1) call s_mpi_reduce_min(dx_min) x_domain%beg = x_cb(-1) x_domain%end = x_cb(m) @@ -364,8 +364,8 @@ contains y_cb(-1:n) = y_cb_glb((start_idx(2) - 1):(start_idx(2) + n)) y_cc(0:n) = (y_cb(0:n) + y_cb(-1:(n - 1)))/2._wp - dy = minval(y_cb(0:n) - y_cb(-1:(n - 1))) - if (num_procs > 1) call s_mpi_reduce_min(dy) + dy_min = minval(y_cb(0:n) - y_cb(-1:(n - 1))) + if (num_procs > 1) call s_mpi_reduce_min(dy_min) y_domain%beg = y_cb(-1) y_domain%end = y_cb(n) @@ -384,8 +384,8 @@ contains z_cb(-1:p) = z_cb_glb((start_idx(3) - 1):(start_idx(3) + p)) z_cc(0:p) = (z_cb(0:p) + z_cb(-1:(p - 1)))/2._wp - dz = minval(z_cb(0:p) - z_cb(-1:(p - 1))) - if (num_procs > 1) call s_mpi_reduce_min(dz) + dz_min = minval(z_cb(0:p) - z_cb(-1:(p - 1))) + if (num_procs > 1) call s_mpi_reduce_min(dz_min) z_domain%beg = z_cb(-1) z_domain%end = z_cb(p) end if diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index 082a471aae..6d54a914a1 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -887,6 +887,9 @@ contains end if end block $:GPU_UPDATE(device='[glb_bounds]') + dx_min = minval(dx) + if (n > 0) dy_min = minval(dy) + if (p > 0) dz_min = minval(dz) if (model_eqns == model_eqns_6eq) call s_initialize_internal_energy_equations(q_cons_ts(1)%vf) if (ib) then From a80e777c91eac6e140637ad8cac0ba36ae73d5f7 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Sat, 25 Jul 2026 21:35:27 -0400 Subject: [PATCH 04/15] Make MPI I/O field binding stage-independent --- src/common/m_mpi_common.fpp | 69 ++++++++++++++----------------- src/post_process/m_data_input.f90 | 6 +-- src/pre_process/m_data_output.fpp | 6 +-- src/pre_process/m_start_up.fpp | 2 +- src/simulation/m_data_output.fpp | 16 +++---- src/simulation/m_start_up.fpp | 12 +++--- 6 files changed, 54 insertions(+), 57 deletions(-) diff --git a/src/common/m_mpi_common.fpp b/src/common/m_mpi_common.fpp index b239183e99..c8a208e299 100644 --- a/src/common/m_mpi_common.fpp +++ b/src/common/m_mpi_common.fpp @@ -109,17 +109,25 @@ contains end subroutine s_mpi_initialize !> Set up MPI I/O data views and variable pointers for parallel file output. - impure subroutine s_initialize_mpi_data(q_cons_vf, ib_markers, beta) + impure subroutine s_initialize_mpi_data(q_cons_vf, ib_markers, ib_mpi_data, beta, qbmm_pb, qbmm_mv) type(scalar_field), dimension(sys_size), intent(in) :: q_cons_vf type(integer_field), optional, intent(in) :: ib_markers + type(mpi_io_ib_var), optional, intent(inout) :: ib_mpi_data type(scalar_field), intent(in), optional :: beta + type(pres_field), intent(in), optional :: qbmm_pb, qbmm_mv integer, dimension(num_dims) :: sizes_glb, sizes_loc #ifdef MFC_MPI integer :: i, j integer :: ierr !< Generic flag used to identify and report MPI errors integer :: alt_sys + logical :: bind_qbmm_fields + + if (present(qbmm_pb) .neqv. present(qbmm_mv)) then + call s_mpi_abort('QBMM MPI I/O requires both pressure and moment fields.') + end if + bind_qbmm_fields = qbmm .and. .not. polytropic .and. present(qbmm_pb) .and. present(qbmm_mv) if (present(beta)) then alt_sys = sys_size + 1 @@ -136,16 +144,11 @@ contains end if ! Additional variables pb and mv for non-polytropic qbmm - if (qbmm .and. .not. polytropic) then + if (bind_qbmm_fields) then do i = 1, nb do j = 1, nnode -#ifdef MFC_PRE_PROCESS - MPI_IO_DATA%var(sys_size + (i - 1)*nnode + j)%sf => pb%sf(0:m,0:n,0:p,j, i) - MPI_IO_DATA%var(sys_size + (i - 1)*nnode + j + nb*nnode)%sf => mv%sf(0:m,0:n,0:p,j, i) -#elif defined (MFC_SIMULATION) - MPI_IO_DATA%var(sys_size + (i - 1)*nnode + j)%sf => pb_ts(1)%sf(0:m,0:n,0:p,j, i) - MPI_IO_DATA%var(sys_size + (i - 1)*nnode + j + nb*nnode)%sf => mv_ts(1)%sf(0:m,0:n,0:p,j, i) -#endif + MPI_IO_DATA%var(sys_size + (i - 1)*nnode + j)%sf => qbmm_pb%sf(0:m,0:n,0:p,j, i) + MPI_IO_DATA%var(sys_size + (i - 1)*nnode + j + nb*nnode)%sf => qbmm_mv%sf(0:m,0:n,0:p,j, i) end do end do end if @@ -166,56 +169,46 @@ contains call MPI_TYPE_COMMIT(MPI_IO_DATA%view(i), ierr) end do -#ifndef MFC_POST_PROCESS - if (qbmm .and. .not. polytropic) then + if (bind_qbmm_fields) then do i = sys_size + 1, sys_size + 2*nb*nnode call MPI_TYPE_CREATE_SUBARRAY(num_dims, sizes_glb, sizes_loc, start_idx, MPI_ORDER_FORTRAN, mpi_p, & & MPI_IO_DATA%view(i), ierr) call MPI_TYPE_COMMIT(MPI_IO_DATA%view(i), ierr) end do end if -#endif -#ifndef MFC_PRE_PROCESS - if (present(ib_markers)) then - MPI_IO_IB_DATA%var%sf => ib_markers%sf(0:m,0:n,0:p) + if (present(ib_markers) .neqv. present(ib_mpi_data)) then + call s_mpi_abort('Immersed-boundary MPI I/O requires both marker and descriptor fields.') + end if + if (present(ib_markers)) then + ib_mpi_data%var%sf => ib_markers%sf(0:m,0:n,0:p) call MPI_TYPE_CREATE_SUBARRAY(num_dims, sizes_glb, sizes_loc, start_idx, MPI_ORDER_FORTRAN, MPI_INTEGER, & - & MPI_IO_IB_DATA%view, ierr) - call MPI_TYPE_COMMIT(MPI_IO_IB_DATA%view, ierr) + & ib_mpi_data%view, ierr) + call MPI_TYPE_COMMIT(ib_mpi_data%view, ierr) end if -#endif #endif end subroutine s_initialize_mpi_data !> Set up MPI I/O data views for downsampled (coarsened) parallel file output. - subroutine s_initialize_mpi_data_ds(q_cons_vf) + subroutine s_initialize_mpi_data_ds(m_ds, n_ds, p_ds, q_cons_vf) - type(scalar_field), dimension(sys_size), intent(in) :: q_cons_vf - integer, dimension(num_dims) :: sizes_loc - integer, dimension(3) :: sf_start_idx + integer, intent(in) :: m_ds, n_ds, p_ds + type(scalar_field), dimension(sys_size), intent(in), optional :: q_cons_vf + integer, dimension(num_dims) :: sizes_loc + integer, dimension(3) :: sf_start_idx #ifdef MFC_MPI - integer :: i, m_ds, n_ds, p_ds, ierr + integer :: i, ierr sf_start_idx = (/0, 0, 0/) -#ifndef MFC_POST_PROCESS - m_ds = int((m + 1)/3) - 1 - n_ds = int((n + 1)/3) - 1 - p_ds = int((p + 1)/3) - 1 -#else - m_ds = m - n_ds = n - p_ds = p -#endif - -#ifdef MFC_POST_PROCESS - do i = 1, sys_size - MPI_IO_DATA%var(i)%sf => q_cons_vf(i)%sf(-1:m_ds + 1,-1:n_ds + 1,-1:p_ds + 1) - end do -#endif + if (present(q_cons_vf)) then + do i = 1, sys_size + MPI_IO_DATA%var(i)%sf => q_cons_vf(i)%sf(-1:m_ds + 1,-1:n_ds + 1,-1:p_ds + 1) + end do + end if ! Define global(g) and local(l) sizes for flow variables sizes_loc(1) = m_ds + 3 if (n > 0) then diff --git a/src/post_process/m_data_input.f90 b/src/post_process/m_data_input.f90 index b24941c46c..da037910ac 100644 --- a/src/post_process/m_data_input.f90 +++ b/src/post_process/m_data_input.f90 @@ -84,7 +84,7 @@ impure subroutine s_setup_mpi_io_params(data_size, m_MOK, n_MOK, p_MOK, WP_MOK, integer(KIND=MPI_OFFSET_KIND), intent(out) :: WP_MOK, MOK, str_MOK, NVARS_MOK if (ib) then - call s_initialize_mpi_data(q_cons_vf, ib_markers) + call s_initialize_mpi_data(q_cons_vf, ib_markers=ib_markers, ib_mpi_data=MPI_IO_IB_DATA) else call s_initialize_mpi_data(q_cons_vf) end if @@ -384,10 +384,10 @@ impure subroutine s_read_parallel_conservative_data(t_step, m_MOK, n_MOK, p_MOK, call MPI_FILE_OPEN(MPI_COMM_SELF, file_loc, MPI_MODE_RDONLY, mpi_info_int, ifile, ierr) if (down_sample) then - call s_initialize_mpi_data_ds(q_cons_temp) + call s_initialize_mpi_data_ds(m, n, p, q_cons_temp) else if (ib) then - call s_initialize_mpi_data(q_cons_vf, ib_markers) + call s_initialize_mpi_data(q_cons_vf, ib_markers=ib_markers, ib_mpi_data=MPI_IO_IB_DATA) else call s_initialize_mpi_data(q_cons_vf) end if diff --git a/src/pre_process/m_data_output.fpp b/src/pre_process/m_data_output.fpp index 0376529957..0fb852b56b 100644 --- a/src/pre_process/m_data_output.fpp +++ b/src/pre_process/m_data_output.fpp @@ -457,9 +457,9 @@ contains call DelayFileAccess(proc_rank) if (down_sample) then - call s_initialize_mpi_data_ds(q_cons_temp) + call s_initialize_mpi_data_ds(m_ds, n_ds, p_ds) else - call s_initialize_mpi_data(q_cons_vf) + call s_initialize_mpi_data(q_cons_vf, qbmm_pb=pb, qbmm_mv=mv) end if if (cfl_dt) then @@ -544,7 +544,7 @@ contains call MPI_FILE_CLOSE(ifile, ierr) else - call s_initialize_mpi_data(q_cons_vf) + call s_initialize_mpi_data(q_cons_vf, qbmm_pb=pb, qbmm_mv=mv) if (cfl_dt) then write (file_loc, '(I0,A)') n_start, '.dat' diff --git a/src/pre_process/m_start_up.fpp b/src/pre_process/m_start_up.fpp index cbf8dd8ace..e36cd571ab 100644 --- a/src/pre_process/m_start_up.fpp +++ b/src/pre_process/m_start_up.fpp @@ -424,7 +424,7 @@ contains if (file_exist) then call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, MPI_MODE_RDONLY, mpi_info_int, ifile, ierr) - call s_initialize_mpi_data(q_cons_vf_in) + call s_initialize_mpi_data(q_cons_vf_in, qbmm_pb=pb, qbmm_mv=mv) data_size = (m + 1)*(n + 1)*(p + 1) diff --git a/src/simulation/m_data_output.fpp b/src/simulation/m_data_output.fpp index ad621c796f..347726c397 100644 --- a/src/simulation/m_data_output.fpp +++ b/src/simulation/m_data_output.fpp @@ -694,12 +694,13 @@ contains call s_int_to_str(t_step, t_step_string) if (down_sample) then - call s_initialize_mpi_data_ds(q_cons_temp_ds) + call s_initialize_mpi_data_ds(m_ds, n_ds, p_ds) else if (ib) then - call s_initialize_mpi_data(q_cons_vf, ib_markers) + call s_initialize_mpi_data(q_cons_vf, ib_markers=ib_markers, ib_mpi_data=MPI_IO_IB_DATA, qbmm_pb=pb_ts(1), & + & qbmm_mv=mv_ts(1)) else - call s_initialize_mpi_data(q_cons_vf) + call s_initialize_mpi_data(q_cons_vf, qbmm_pb=pb_ts(1), qbmm_mv=mv_ts(1)) end if end if @@ -714,7 +715,7 @@ contains call s_mpi_barrier() call DelayFileAccess(proc_rank) - call s_initialize_mpi_data(q_cons_vf) + call s_initialize_mpi_data(q_cons_vf, qbmm_pb=pb_ts(1), qbmm_mv=mv_ts(1)) write (file_loc, '(I0,A,i7.7,A)') t_step, '_', proc_rank, '.dat' file_loc = trim(case_dir) // '/restart_data/lustre_' // trim(t_step_string) // trim(mpiiofs) // trim(file_loc) @@ -780,11 +781,12 @@ contains end if else if (ib) then - call s_initialize_mpi_data(q_cons_vf, ib_markers) + call s_initialize_mpi_data(q_cons_vf, ib_markers=ib_markers, ib_mpi_data=MPI_IO_IB_DATA, qbmm_pb=pb_ts(1), & + & qbmm_mv=mv_ts(1)) else if (present(beta)) then - call s_initialize_mpi_data(q_cons_vf, beta=beta) + call s_initialize_mpi_data(q_cons_vf, beta=beta, qbmm_pb=pb_ts(1), qbmm_mv=mv_ts(1)) else - call s_initialize_mpi_data(q_cons_vf) + call s_initialize_mpi_data(q_cons_vf, qbmm_pb=pb_ts(1), qbmm_mv=mv_ts(1)) end if write (file_loc, '(I0,A)') t_step, '.dat' diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index 6d54a914a1..e9577cc43b 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -367,12 +367,13 @@ contains call MPI_FILE_OPEN(MPI_COMM_SELF, file_loc, MPI_MODE_RDONLY, mpi_info_int, ifile, ierr) if (down_sample) then - call s_initialize_mpi_data_ds(q_cons_vf) + call s_initialize_mpi_data_ds(m_ds, n_ds, p_ds) else if (ib) then - call s_initialize_mpi_data(q_cons_vf, ib_markers) + call s_initialize_mpi_data(q_cons_vf, ib_markers=ib_markers, ib_mpi_data=MPI_IO_IB_DATA, & + & qbmm_pb=pb_ts(1), qbmm_mv=mv_ts(1)) else - call s_initialize_mpi_data(q_cons_vf) + call s_initialize_mpi_data(q_cons_vf, qbmm_pb=pb_ts(1), qbmm_mv=mv_ts(1)) end if end if @@ -443,9 +444,10 @@ contains call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, MPI_MODE_RDONLY, mpi_info_int, ifile, ierr) if (ib) then - call s_initialize_mpi_data(q_cons_vf, ib_markers) + call s_initialize_mpi_data(q_cons_vf, ib_markers=ib_markers, ib_mpi_data=MPI_IO_IB_DATA, qbmm_pb=pb_ts(1), & + & qbmm_mv=mv_ts(1)) else - call s_initialize_mpi_data(q_cons_vf) + call s_initialize_mpi_data(q_cons_vf, qbmm_pb=pb_ts(1), qbmm_mv=mv_ts(1)) end if data_size = (m + 1)*(n + 1)*(p + 1) From b2dfcadf8dce917f1df97090249431756b214a54 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Sat, 25 Jul 2026 21:46:20 -0400 Subject: [PATCH 05/15] Make Dirichlet boundary policy explicit --- src/common/m_boundary_common.fpp | 9 +++++++-- src/common/m_boundary_primitives.fpp | 12 ++++++++---- src/simulation/m_start_up.fpp | 2 +- 3 files changed, 16 insertions(+), 7 deletions(-) diff --git a/src/common/m_boundary_common.fpp b/src/common/m_boundary_common.fpp index 07776d1c10..3dfa30528b 100644 --- a/src/common/m_boundary_common.fpp +++ b/src/common/m_boundary_common.fpp @@ -34,9 +34,14 @@ module m_boundary_common contains !> Allocate and set up boundary condition buffer arrays for all coordinate directions. - impure subroutine s_initialize_boundary_common_module() + impure subroutine s_initialize_boundary_common_module(use_dirichlet_buffers) - integer :: i, j, sys_size_alloc + integer :: i, j, sys_size_alloc + logical, intent(in), optional :: use_dirichlet_buffers + + dirichlet_from_buffers = .false. + if (present(use_dirichlet_buffers)) dirichlet_from_buffers = use_dirichlet_buffers + $:GPU_UPDATE(device='[dirichlet_from_buffers]') @:ALLOCATE(bc_buffers(1:3, 1:2)) diff --git a/src/common/m_boundary_primitives.fpp b/src/common/m_boundary_primitives.fpp index 00a61bd4a5..f284283b2d 100644 --- a/src/common/m_boundary_primitives.fpp +++ b/src/common/m_boundary_primitives.fpp @@ -18,6 +18,9 @@ module m_boundary_primitives type(scalar_field), dimension(:,:), allocatable :: bc_buffers $:GPU_DECLARE(create='[bc_buffers]') + logical :: dirichlet_from_buffers + $:GPU_DECLARE(create='[dirichlet_from_buffers]') + contains !> Fill ghost cells by copying the nearest boundary cell value along the specified direction. @@ -904,7 +907,11 @@ contains integer :: j, i type(scalar_field), optional, intent(inout) :: q_T_sf -#ifdef MFC_SIMULATION + if (.not. dirichlet_from_buffers) then + call s_ghost_cell_extrapolation(q_prim_vf, bc_dir, bc_loc, k, l, q_T_sf) + return + end if + if (bc_dir == 1) then !< x-direction if (bc_loc == -1) then ! bc_x%beg do i = 1, sys_size @@ -982,9 +989,6 @@ contains end if #:endif end if -#else - call s_ghost_cell_extrapolation(q_prim_vf, bc_dir, bc_loc, k, l, q_T_sf) -#endif end subroutine s_dirichlet diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index e9577cc43b..66a7c57a99 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -845,7 +845,7 @@ contains call s_initialize_derived_variables_module() call s_initialize_time_steppers_module() - call s_initialize_boundary_common_module() + call s_initialize_boundary_common_module(use_dirichlet_buffers=.true.) if (down_sample) then m_ds = int((m + 1)/3) - 1 From 6f7bf444ceabde78d5df1d3f09e7bd154d5dbcc6 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Sat, 25 Jul 2026 21:55:13 -0400 Subject: [PATCH 06/15] Make common validation policies explicit --- src/common/m_checker_common.fpp | 22 ++++++++++------------ src/common/m_finite_differences.fpp | 5 ----- src/post_process/m_start_up.fpp | 2 +- src/pre_process/m_start_up.fpp | 2 +- src/simulation/m_start_up.fpp | 2 +- 5 files changed, 13 insertions(+), 20 deletions(-) diff --git a/src/common/m_checker_common.fpp b/src/common/m_checker_common.fpp index dfdcff7678..f7b7b42d74 100644 --- a/src/common/m_checker_common.fpp +++ b/src/common/m_checker_common.fpp @@ -20,33 +20,33 @@ module m_checker_common contains !> Checks compatibility of parameters in the input file. Used by all three stages - impure subroutine s_check_inputs_common + impure subroutine s_check_inputs_common(check_total_cells, n_global) -#ifndef MFC_SIMULATION - call s_check_total_cells -#endif + logical, intent(in) :: check_total_cells + 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 #:endif end subroutine s_check_inputs_common -#ifndef MFC_SIMULATION !> Verify that the total number of grid cells meets the minimum required by the number of dimensions and MPI ranks. - impure subroutine s_check_total_cells + impure subroutine s_check_total_cells(n_global) - character(len=18) :: numStr !< for int to string conversion - integer(kind=8) :: min_cells + character(len=18) :: numStr !< for int to string conversion + integer(kind=8) :: min_cells + integer(kind=8), intent(in) :: n_global min_cells = int(2, kind=8)**int(min(1, m) + min(1, n) + min(1, p), kind=8)*int(num_procs, kind=8) call s_int_to_str(2**(min(1, m) + min(1, n) + min(1, p))*num_procs, numStr) - @:PROHIBIT(nGlobal < min_cells, & + @:PROHIBIT(n_global < min_cells, & & "Total number of cells must be at least (2^[number of dimensions])*num_procs, " // "which is currently " & & // trim(numStr)) end subroutine s_check_total_cells -#endif !> Check that simulation parameters stay within AMD GPU compiler limits when case optimization is disabled. impure subroutine s_check_amd @@ -59,6 +59,4 @@ contains end subroutine s_check_amd -#ifndef MFC_POST_PROCESS -#endif end module m_checker_common diff --git a/src/common/m_finite_differences.fpp b/src/common/m_finite_differences.fpp index 00cf6a0c99..5639890daa 100644 --- a/src/common/m_finite_differences.fpp +++ b/src/common/m_finite_differences.fpp @@ -38,11 +38,6 @@ contains lE = q end if -#ifdef MFC_POST_PROCESS - if (allocated(fd_coeff_s)) deallocate (fd_coeff_s) - allocate (fd_coeff_s(-fd_number_in:fd_number_in,lB:lE)) -#endif - ! Computing the 1st order finite-difference coefficients if (fd_order_in == 1) then do i = lB, lE diff --git a/src/post_process/m_start_up.fpp b/src/post_process/m_start_up.fpp index 021ea86705..97d67a293a 100644 --- a/src/post_process/m_start_up.fpp +++ b/src/post_process/m_start_up.fpp @@ -120,7 +120,7 @@ contains call s_mpi_abort('Unsupported choice for the value of ' // 'case_dir. Exiting.') end if - call s_check_inputs_common() + call s_check_inputs_common(check_total_cells=.true., n_global=nGlobal) call s_check_inputs() end subroutine s_check_input_file diff --git a/src/pre_process/m_start_up.fpp b/src/pre_process/m_start_up.fpp index e36cd571ab..a7676b3104 100644 --- a/src/pre_process/m_start_up.fpp +++ b/src/pre_process/m_start_up.fpp @@ -128,7 +128,7 @@ contains call s_mpi_abort('Unsupported choice for the value of case_dir.' // 'Exiting.') end if - call s_check_inputs_common() + call s_check_inputs_common(check_total_cells=.true., n_global=nGlobal) call s_check_inputs() call s_check_patches() diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index 66a7c57a99..764bb7e253 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -144,7 +144,7 @@ contains call s_mpi_abort(trim(file_path) // ' is missing. Exiting.') end if - call s_check_inputs_common() + call s_check_inputs_common(check_total_cells=.false., n_global=0_8) call s_check_inputs() end subroutine s_check_input_file From d086e7a1dde08a7b4a86359bd05a38f2907b1f13 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Sat, 25 Jul 2026 22:06:45 -0400 Subject: [PATCH 07/15] Make common global ownership stage-independent --- src/common/m_global_parameters_common.fpp | 73 +---------- src/post_process/m_global_parameters.fpp | 2 +- src/pre_process/m_global_parameters.fpp | 2 +- src/simulation/m_global_parameters.fpp | 2 +- .../mfc/params/generators/fortran_gen.py | 119 +++++++++++++++++- .../mfc/params_tests/test_fortran_gen.py | 19 +++ 6 files changed, 141 insertions(+), 76 deletions(-) diff --git a/src/common/m_global_parameters_common.fpp b/src/common/m_global_parameters_common.fpp index fc396b3b3a..a5585202a0 100644 --- a/src/common/m_global_parameters_common.fpp +++ b/src/common/m_global_parameters_common.fpp @@ -27,20 +27,6 @@ module m_global_parameters_common ! num_dims, num_vels, weno_polyn, muscl_polyn, weno_num_stencils, wenojs, igr, etc. #:include 'generated_case_opt_decls.fpp' - ! For pre_process and post_process: num_dims and num_vels are declared manually here - ! (sim gets them from generated_case_opt_decls.fpp above) -#ifndef MFC_SIMULATION - integer :: num_dims !< Number of spatial dimensions - integer :: num_vels !< Number of velocity components (different from num_dims for mhd) -#endif - - ! For pre_process: weno_polyn and muscl_polyn are declared manually here - ! (sim gets them from generated_case_opt_decls.fpp; post does not use them) -#ifdef MFC_PRE_PROCESS - integer :: weno_polyn !< Degree of the WENO polynomials - integer :: muscl_polyn !< Degree of the MUSCL polynomials -#endif - !> @name Annotations of the structure of the state and flux vectors in terms of the size and configuration of the system of !! equations !> @{ @@ -90,49 +76,9 @@ module m_global_parameters_common !> Minimum cell widths. These are distinct from the per-cell width arrays named dx, dy, and dz in simulation and post-process. real(wp) :: dx_min, dy_min, dz_min -#ifdef MFC_SIMULATION $:GPU_DECLARE(create='[sys_size, eqn_idx, b_size, tensor_size]') $:GPU_DECLARE(create='[shear_num, shear_indices, shear_BC_flip_num, shear_BC_flip_indices]') - ! Device residency for namelist/case-opt state declared above via the generated - ! includes: declare directives must live in the declaring module (Cray ftn rejects - ! declare-target on use-associated names), so these moved here from simulation. - $:GPU_DECLARE(create='[cyl_coord]') - $:GPU_DECLARE(create='[dt, m, n, p]') - $:GPU_DECLARE(create='[cfl_target]') - $:GPU_DECLARE(create='[int_comp, ic_eps, ic_beta]') - $:GPU_DECLARE(create='[muscl_eps]') - $:GPU_DECLARE(create='[mpp_lim, model_eqns, mixture_err, alt_soundspeed]') - $:GPU_DECLARE(create='[avg_state, mp_weno, weno_eps, teno_CT, hypoelasticity]') - $:GPU_DECLARE(create='[hyperelasticity, elasticity, low_Mach]') - $:GPU_DECLARE(create='[cont_damage, hyper_cleaning]') - $:GPU_DECLARE(create='[relax, relax_model, palpha_eps, ptgalpha_eps]') - $:GPU_DECLARE(create='[down_sample]') - $:GPU_DECLARE(create='[fd_order]') - $:GPU_DECLARE(create='[rhoref, pref]') - $:GPU_DECLARE(create='[ib, num_ibs]') - $:GPU_DECLARE(create='[ib_coefficient_of_friction]') - $:GPU_DECLARE(create='[Ca, Web, Re_inv]') - $:GPU_DECLARE(create='[bubbles_euler, polytropic, polydisperse]') - $:GPU_DECLARE(create='[adv_n, adap_dt, adap_dt_tol, adap_dt_max_iters]') - $:GPU_DECLARE(create='[bubble_model, thermal]') - $:GPU_DECLARE(create='[poly_sigma]') - $:GPU_DECLARE(create='[qbmm, pi_fac]') - $:GPU_DECLARE(create='[R0ref]') - $:GPU_DECLARE(create='[acoustic_source, num_source]') - $:GPU_DECLARE(create='[sigma, surface_tension]') - $:GPU_DECLARE(create='[bubbles_lagrange]') - $:GPU_DECLARE(create='[Bx0]') - $:GPU_DECLARE(create='[tau_star, cont_damage_s, alpha_bar]') - $:GPU_DECLARE(create='[hyper_cleaning_speed, hyper_cleaning_tau]') - $:GPU_DECLARE(create='[synthetic_turbulence, num_turbulent_sources, synth_U_inf]') - #:if not MFC_CASE_OPTIMIZATION - $:GPU_DECLARE(create='[num_dims, num_vels, weno_polyn, weno_order]') - $:GPU_DECLARE(create='[weno_num_stencils, num_fluids, wenojs]') - $:GPU_DECLARE(create='[mapped_weno, wenoz, teno, wenoz_q, mhd, relativity]') - $:GPU_DECLARE(create='[igr_iter_solver, igr_order, viscous, igr_pres_lim, igr]') - $:GPU_DECLARE(create='[recon_type, muscl_order, muscl_polyn, muscl_lim]') - #:endif -#endif + $:GPU_DECLARE(create='[elasticity]') !> @name Processor coordinates and parallel-IO addressing (identical declaration across all three targets) !> @{ @@ -163,10 +109,11 @@ contains !! - qbmm_idx allocations and fills (diverge between pre vs sim/post) !! - sim-only: gam = bub_pp%gam_g, nmomsp/nmomtot, Re_idx allocation, GPU_UPDATE calls !! - post-only: beta_idx increment (bubbles_lagrange path), offset/grid allocations - impure subroutine s_initialize_eqn_idx(nmom_in, nb_in) + impure subroutine s_initialize_eqn_idx(nmom_in, nb_in, six_eqn_alf_is_advected) integer, intent(in) :: nmom_in integer, intent(in) :: nb_in + logical, intent(in) :: six_eqn_alf_is_advected ! Gamma/Pi_inf Model @@ -252,9 +199,7 @@ contains 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 -#ifdef MFC_SIMULATION - eqn_idx%alf = eqn_idx%adv%end -#endif + 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 @@ -352,7 +297,6 @@ contains integer :: ierr !< Generic flag used to identify and report MPI errors #endif -#ifdef MFC_SIMULATION ! Under case-optimization, num_dims and num_vels are compile-time parameters; skip assignment. #:if not MFC_CASE_OPTIMIZATION num_dims = 1 + min(1, n) + min(1, p) @@ -363,15 +307,6 @@ contains num_vels = num_dims end if #:endif -#else - num_dims = 1 + min(1, n) + min(1, p) - - if (mhd) then - num_vels = 3 - else - num_vels = num_dims - end if -#endif allocate (proc_coords(1:num_dims)) diff --git a/src/post_process/m_global_parameters.fpp b/src/post_process/m_global_parameters.fpp index f90bbcb0d1..ac67023f74 100644 --- a/src/post_process/m_global_parameters.fpp +++ b/src/post_process/m_global_parameters.fpp @@ -326,7 +326,7 @@ contains if (model_eqns == model_eqns_5eq .and. qbmm) nmom = 6 ! Populate eqn_idx, sys_size, b_size, tensor_size, elasticity, shear_* (shared logic) - call s_initialize_eqn_idx(nmom, nb) + call s_initialize_eqn_idx(nmom, nb, six_eqn_alf_is_advected=.false.) ! post-only: 6eq alf is a dummy (no void fraction in 6eq) if (model_eqns == model_eqns_6eq) eqn_idx%alf = 1 diff --git a/src/pre_process/m_global_parameters.fpp b/src/pre_process/m_global_parameters.fpp index 44dd397f90..b13f17b4db 100644 --- a/src/pre_process/m_global_parameters.fpp +++ b/src/pre_process/m_global_parameters.fpp @@ -447,7 +447,7 @@ contains if (model_eqns == model_eqns_5eq .and. bubbles_euler .and. qbmm .and. nnode == 4) nmom = 6 ! Populate eqn_idx, sys_size, b_size, tensor_size, elasticity, shear_* (shared logic) - call s_initialize_eqn_idx(nmom, nb) + call s_initialize_eqn_idx(nmom, nb, six_eqn_alf_is_advected=.false.) ! Per-target (pre_process): qbmm_idx allocations and fills if (model_eqns == model_eqns_5eq .and. bubbles_euler) then diff --git a/src/simulation/m_global_parameters.fpp b/src/simulation/m_global_parameters.fpp index 1ef48eaf3d..6f88f99c13 100644 --- a/src/simulation/m_global_parameters.fpp +++ b/src/simulation/m_global_parameters.fpp @@ -708,7 +708,7 @@ contains Re_size_max = 0 ! Populate eqn_idx, sys_size, b_size, tensor_size, elasticity, shear_* (shared logic) - call s_initialize_eqn_idx(nmom, nb) + call s_initialize_eqn_idx(nmom, nb, six_eqn_alf_is_advected=.true.) ! sim-only: GPU update for shear state after s_initialize_eqn_idx populated it if (model_eqns == model_eqns_5eq .or. model_eqns == model_eqns_6eq) then diff --git a/toolchain/mfc/params/generators/fortran_gen.py b/toolchain/mfc/params/generators/fortran_gen.py index d5d6bc02ea..a77d0d42d8 100644 --- a/toolchain/mfc/params/generators/fortran_gen.py +++ b/toolchain/mfc/params/generators/fortran_gen.py @@ -28,6 +28,97 @@ _VALID_TARGETS = ("pre", "sim", "post") +# Scalar simulation parameters whose generated declarations own device storage. +# Case-optimization parameters are emitted by generated_case_opt_decls.fpp; all +# other names are emitted by generated_decls.fpp. +SIM_GPU_DECL_VARS = { + "Bx0", + "Ca", + "R0ref", + "Re_inv", + "Web", + "acoustic_source", + "adap_dt", + "adap_dt_max_iters", + "adap_dt_tol", + "adv_n", + "alpha_bar", + "alt_soundspeed", + "avg_state", + "bubble_model", + "bubbles_euler", + "bubbles_lagrange", + "cfl_target", + "cont_damage", + "cont_damage_s", + "cyl_coord", + "down_sample", + "dt", + "fd_order", + "hyper_cleaning", + "hyper_cleaning_speed", + "hyper_cleaning_tau", + "hyperelasticity", + "hypoelasticity", + "ib", + "ib_coefficient_of_friction", + "ic_beta", + "ic_eps", + "igr", + "igr_iter_solver", + "igr_order", + "igr_pres_lim", + "int_comp", + "low_Mach", + "m", + "mapped_weno", + "mixture_err", + "model_eqns", + "mp_weno", + "mpp_lim", + "muscl_eps", + "muscl_lim", + "muscl_order", + "muscl_polyn", + "n", + "num_dims", + "num_fluids", + "num_ibs", + "num_source", + "num_turbulent_sources", + "num_vels", + "p", + "palpha_eps", + "pi_fac", + "poly_sigma", + "polydisperse", + "polytropic", + "pref", + "ptgalpha_eps", + "qbmm", + "recon_type", + "relativity", + "relax", + "relax_model", + "rhoref", + "sigma", + "surface_tension", + "synth_U_inf", + "synthetic_turbulence", + "tau_star", + "teno", + "teno_CT", + "thermal", + "viscous", + "weno_eps", + "weno_num_stencils", + "weno_order", + "weno_polyn", + "wenojs", + "wenoz", + "wenoz_q", +} + def _check_target(target: str) -> None: if target not in _VALID_TARGETS: @@ -117,6 +208,7 @@ def generate_decls_fpp(target: str) -> str: """Return Fortran declarations (scalars + known arrays) for a target.""" _check_target(target) lines = [_HEADER.rstrip()] + declared_names = set() for name in _decl_vars_for_target(target): if not _is_simple_scalar(name): continue @@ -132,6 +224,7 @@ def generate_decls_fpp(target: str) -> str: ftype = fortran_type_decl(member) dim = FORTRAN_ARRAY_DIMS[name] lines.append(f"{(ftype + ', dimension(' + dim + ')').ljust(_ARRAY_DECL_COL)}:: {name}") + declared_names.add(name) continue param = REGISTRY.all_params.get(name) if param is None: @@ -139,6 +232,7 @@ def generate_decls_fpp(target: str) -> str: if any(k.startswith(f"{name}(") for k in REGISTRY.all_params): raise ValueError(f"{name!r} has indexed variants (e.g. {name}(1)) but is missing from " "FORTRAN_ARRAY_DIMS. Add it there with its Fortran dimension expression.") lines.append(f"{fortran_type_decl(param).ljust(_DECL_COL)}:: {name}") + declared_names.add(name) for name, (ftype, dim, gpu, desc) in TYPED_DECLS.items(): if name not in NAMELIST_VARS or target not in NAMELIST_VARS[name]: continue @@ -148,8 +242,12 @@ def generate_decls_fpp(target: str) -> str: padded += " " doc = f" !< {desc}" if desc else "" lines.append(f"{padded}:: {name}{doc}") + declared_names.add(name) if gpu and target == "sim": lines.append(f"$:GPU_DECLARE(create='[{name}]')") + if target == "sim": + for name in sorted(SIM_GPU_DECL_VARS & declared_names): + lines.append(f"$:GPU_DECLARE(create='[{name}]')") return "\n".join(lines) + "\n" @@ -179,6 +277,7 @@ def generate_constants_fpp() -> str: ("weno_num_stencils", "integer", "Number of stencils for WENO reconstruction"), ("wenojs", "logical", "WENO-JS (default)"), ] +COMMON_CASE_OPT_EXTRA_NAMES = {"num_dims", "num_vels", "weno_polyn", "muscl_polyn"} _CASE_OPT_DECL_COL = 24 # '::' alignment for case-opt declarations @@ -227,6 +326,10 @@ def generate_case_opt_decls_fpp() -> str: if_lines.append(f" {ftype}, parameter :: {name} = ${{{name}}}$ !< {desc}") else_lines.append(f" {ftype.ljust(_CASE_OPT_DECL_COL)}:: {name}") + declared_names = {name for name, _, _ in CASE_OPT_EXTRA_LINES} | set(params_to_emit) + for name in sorted(SIM_GPU_DECL_VARS & declared_names): + else_lines.append(f" $:GPU_DECLARE(create='[{name}]')") + parts = [_HEADER.rstrip(), "#:if MFC_CASE_OPTIMIZATION"] parts.extend(if_lines) parts.append("#:else") @@ -235,6 +338,15 @@ def generate_case_opt_decls_fpp() -> str: return "\n".join(parts) + "\n" +def generate_common_extra_decls_fpp() -> str: + """Return computed-scalar declarations needed by common pre/post code.""" + lines = [_HEADER.rstrip()] + for name, ftype, _ in CASE_OPT_EXTRA_LINES: + if name in COMMON_CASE_OPT_EXTRA_NAMES: + lines.append(f"{ftype.ljust(_CASE_OPT_DECL_COL)}:: {name}") + return "\n".join(lines) + "\n" + + # Struct roots in NAMELIST_VARS whose member-level broadcasts are irregular # (per-target member subsets, grouped array members, nested loops, etc.). # These are kept in the manual residue of m_mpi_proxy.fpp. @@ -572,9 +684,8 @@ def get_generated_files(build_dir: Path) -> List[Tuple[Path, str]]: Paths match the cmake include directory structure: build_dir/include/{full_target}/generated_{namelist,decls,constants,case_opt_decls,bcast}.fpp - Every target gets generated_case_opt_decls.fpp: real content for simulation, - a header-only stub for the others (Fypp resolves #:include at parse time, so - the file must exist for every target even inside a dead conditional). + 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. """ result = [] @@ -585,7 +696,7 @@ def get_generated_files(build_dir: Path) -> List[Tuple[Path, str]]: result.append((inc / "generated_constants.fpp", generate_constants_fpp())) for short, full in TARGETS: inc = build_dir / "include" / full - content = generate_case_opt_decls_fpp() if short == "sim" else _HEADER + "! (no case-optimization declarations for this target)\n" + content = generate_case_opt_decls_fpp() if short == "sim" else generate_common_extra_decls_fpp() result.append((inc / "generated_case_opt_decls.fpp", content)) for short, full in TARGETS: inc = build_dir / "include" / full diff --git a/toolchain/mfc/params_tests/test_fortran_gen.py b/toolchain/mfc/params_tests/test_fortran_gen.py index 0771f0f3c7..9752e6bdd9 100644 --- a/toolchain/mfc/params_tests/test_fortran_gen.py +++ b/toolchain/mfc/params_tests/test_fortran_gen.py @@ -323,6 +323,25 @@ def test_generate_case_opt_decls_fpp(): assert "AUTO-GENERATED" in out +def test_common_computed_scalars_are_declared_for_every_target(): + from pathlib import Path + + from mfc.params.generators.fortran_gen import get_generated_files + + files = {path.parent.name: content for path, content in get_generated_files(Path("/build")) if path.name == "generated_case_opt_decls.fpp"} + for target in ("pre_process", "simulation", "post_process"): + for name in ("num_dims", "num_vels", "weno_polyn", "muscl_polyn"): + assert f":: {name}" in files[target] + + +def test_simulation_generated_scalars_own_gpu_declarations(): + from mfc.params.generators.fortran_gen import SIM_GPU_DECL_VARS, generate_case_opt_decls_fpp, generate_decls_fpp + + generated = generate_decls_fpp("sim") + generate_case_opt_decls_fpp() + for name in SIM_GPU_DECL_VARS: + assert generated.count(f"$:GPU_DECLARE(create='[{name}]')") == 1 + + # ── generate_bcast_fpp tests ────────────────────────────────────────────────── From b50464544cc59c90a587dc021a62e36b6473fb70 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Sat, 25 Jul 2026 22:17:18 -0400 Subject: [PATCH 08/15] Make common MPI transport stage-independent --- src/common/include/macros.fpp | 2 -- src/common/m_mpi_common.fpp | 39 ++++++++++++++------------------- src/post_process/m_start_up.fpp | 2 +- src/pre_process/m_start_up.fpp | 2 +- src/simulation/m_start_up.fpp | 2 +- 5 files changed, 19 insertions(+), 28 deletions(-) diff --git a/src/common/include/macros.fpp b/src/common/include/macros.fpp index da55c016c3..a221fcf8c5 100644 --- a/src/common/include/macros.fpp +++ b/src/common/include/macros.fpp @@ -16,7 +16,6 @@ ! rank. That's because for both cudaMemAdvise (preferred location) and cudaMemPrefetchAsync we use location = device_id = 0. For an ! example see misc/nvidia_uvm/bind.sh. #:def PREFER_GPU(*args) -#ifdef MFC_SIMULATION #ifdef __NVCOMPILER_GPU_UNIFIED_MEM block ! NVIDIA CUDA Fortran 25.3+: uses submodules (cuda_runtime_api, gpu_reductions, sort) See @@ -52,7 +51,6 @@ end if end block #endif -#endif #:enddef ! Allocate and create GPU device memory diff --git a/src/common/m_mpi_common.fpp b/src/common/m_mpi_common.fpp index c8a208e299..a75951c885 100644 --- a/src/common/m_mpi_common.fpp +++ b/src/common/m_mpi_common.fpp @@ -40,10 +40,19 @@ module m_mpi_common integer(kind=8) :: halo_size $:GPU_DECLARE(create='[halo_size]') + logical, private :: exchange_all_chemistry_temperatures + logical, private :: use_rdma_transport + contains !> Initialize the module. - impure subroutine s_initialize_mpi_common_module + impure subroutine s_initialize_mpi_common_module(exchange_all_chemistry_temperatures_in, use_rdma_transport_in) + + logical, intent(in) :: exchange_all_chemistry_temperatures_in + logical, intent(in) :: use_rdma_transport_in + + exchange_all_chemistry_temperatures = exchange_all_chemistry_temperatures_in + use_rdma_transport = use_rdma_transport_in #ifdef MFC_MPI ! Allocating buff_send/recv and. Please note that for the sake of simplicity, both variables are provided sufficient storage @@ -51,7 +60,7 @@ contains if (qbmm .and. .not. polytropic) then v_size = sys_size + 2*nb*nnode - else if (chemistry .and. chem_params%diffusion) then + else if (chemistry .and. (chem_params%diffusion .or. exchange_all_chemistry_temperatures)) then v_size = sys_size + 1 else v_size = sys_size @@ -311,7 +320,6 @@ contains Rc_min_glb = Rc_min_loc ccfl_max_glb = ccfl_max_loc -#ifdef MFC_SIMULATION #ifdef MFC_MPI block integer :: ierr @@ -346,7 +354,6 @@ contains end if if (bubbles_lagrange) bubs_glb = bubs_loc -#endif #endif end subroutine s_mpi_reduce_stability_criteria_extrema @@ -564,14 +571,10 @@ contains v_size = nVar + 2*nb*nnode buffer_counts = (/buff_size*v_size*(n + 1)*(p + 1), buff_size*v_size*(m + 2*buff_size + 1)*(p + 1), & & buff_size*v_size*(m + 2*buff_size + 1)*(n + 2*buff_size + 1)/) -#ifdef MFC_SIMULATION - else if (present(q_T_sf) .and. chemistry .and. chem_params%diffusion) then -#else - else if (present(q_T_sf) .and. chemistry) then - ! post_process converts cons->prim over the ghost-inclusive bounds, so the temperature - ! Newton guess must be valid at rank seams for EVERY chemistry run (not only diffusion): + else if (present(q_T_sf) .and. chemistry .and. (chem_params%diffusion .or. exchange_all_chemistry_temperatures)) then + ! Consumers that convert over ghost-inclusive bounds request temperature exchange for every chemistry run. + ! The temperature Newton guess must be valid at rank seams even when diffusion is disabled: ! an unexchanged seam ghost is an uninitialized guess -> NaN T/pres/c in the output -#endif chem_diff_comm = .true. v_size = nVar + 1 buffer_counts = (/buff_size*v_size*(n + 1)*(p + 1), buff_size*v_size*(m + 2*buff_size + 1)*(p + 1), & @@ -800,9 +803,8 @@ contains call nvtxEndRange ! Packbuf ! Send/Recv -#ifdef MFC_SIMULATION #:for rdma_mpi in [False, True] - if (rdma_mpi .eqv. ${'.true.' if rdma_mpi else '.false.'}$) then + if (use_rdma_transport .eqv. ${'.true.' if rdma_mpi else '.false.'}$) then #:if rdma_mpi #:call GPU_HOST_DATA(use_device_addr='[buff_send, buff_recv]') call nvtxStartRange("RHS-COMM-SENDRECV-RDMA") @@ -830,10 +832,6 @@ contains #:endif end if #:endfor -#else - call MPI_SENDRECV(buff_send, buffer_count, mpi_p, dst_proc, send_tag, buff_recv, buffer_count, mpi_p, src_proc, recv_tag, & - & MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr) -#endif ! Unpack Received Buffer call nvtxStartRange("RHS-COMM-UNPACKBUF") @@ -1190,9 +1188,8 @@ contains call nvtxEndRange ! Packbuf ! Send/Recv -#ifdef MFC_SIMULATION #:for rdma_mpi in [False, True] - if (rdma_mpi .eqv. ${'.true.' if rdma_mpi else '.false.'}$) then + if (use_rdma_transport .eqv. ${'.true.' if rdma_mpi else '.false.'}$) then #:if rdma_mpi #:call GPU_HOST_DATA(use_device_addr='[buff_send, buff_recv]') call nvtxStartRange("BETA-COMM-SENDRECV-RDMA") @@ -1220,10 +1217,6 @@ contains #:endif end if #:endfor -#else - call MPI_SENDRECV(buff_send, buffer_count, mpi_p, dst_proc, send_tag, buff_recv, buffer_count, mpi_p, src_proc, recv_tag, & - & MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr) -#endif ! Unpack Received Buffer (skip if no source rank) call nvtxStartRange("BETA-COMM-UNPACKBUF") diff --git a/src/post_process/m_start_up.fpp b/src/post_process/m_start_up.fpp index 97d67a293a..707d5b117a 100644 --- a/src/post_process/m_start_up.fpp +++ b/src/post_process/m_start_up.fpp @@ -784,7 +784,7 @@ contains end if if (num_procs > 1) then call s_initialize_mpi_proxy_module() - call s_initialize_mpi_common_module() + call s_initialize_mpi_common_module(exchange_all_chemistry_temperatures_in=.true., use_rdma_transport_in=.false.) end if call s_initialize_boundary_common_module() call s_initialize_variables_conversion_module(store_mixture_fields=.true., lagrange_beta_index=beta_idx) diff --git a/src/pre_process/m_start_up.fpp b/src/pre_process/m_start_up.fpp index a7676b3104..cc60ea5fe4 100644 --- a/src/pre_process/m_start_up.fpp +++ b/src/pre_process/m_start_up.fpp @@ -476,7 +476,7 @@ contains if (bubbles_euler .or. bubbles_lagrange) then call s_initialize_bubbles_model() end if - call s_initialize_mpi_common_module() + call s_initialize_mpi_common_module(exchange_all_chemistry_temperatures_in=.false., use_rdma_transport_in=.false.) call s_initialize_data_output_module() call s_initialize_variables_conversion_module() call s_initialize_grid_module() diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index 764bb7e253..c39f52bc80 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -816,7 +816,7 @@ contains if (bubbles_euler .or. bubbles_lagrange) then call s_initialize_bubbles_model() end if - call s_initialize_mpi_common_module() + call s_initialize_mpi_common_module(exchange_all_chemistry_temperatures_in=.false., use_rdma_transport_in=rdma_mpi) call s_initialize_mpi_proxy_module() call s_initialize_variables_conversion_module(enforce_density_floor=.true., preserve_qbmm_number=.true.) if (grid_geometry == 3) call s_initialize_fftw_module() From 913573c5005a09220a7816983cc0bd4e3fb007e3 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Sat, 25 Jul 2026 22:47:28 -0400 Subject: [PATCH 09/15] Make MPI decomposition stage-independent --- src/common/m_mpi_common.fpp | 176 +++++++++++++++----------------- src/post_process/m_start_up.fpp | 11 +- src/pre_process/m_start_up.fpp | 9 +- src/simulation/m_start_up.fpp | 2 +- 4 files changed, 99 insertions(+), 99 deletions(-) diff --git a/src/common/m_mpi_common.fpp b/src/common/m_mpi_common.fpp index a75951c885..7220823906 100644 --- a/src/common/m_mpi_common.fpp +++ b/src/common/m_mpi_common.fpp @@ -17,10 +17,12 @@ module m_mpi_common use m_helper use ieee_arithmetic use m_nvtx - use m_constants, only: recon_type_weno, format_silo + use m_constants, only: recon_type_weno implicit none + private :: s_apply_decomposition_policies + integer, private :: v_size $:GPU_DECLARE(create='[v_size]') @@ -1314,7 +1316,12 @@ contains !> The purpose of this procedure is to optimally decompose the computational domain among the available processors. This is !! performed by attempting to award each processor, in each of the coordinate directions, approximately the same number of !! cells, and then recomputing the affected global parameters. - subroutine s_mpi_decompose_computational_domain + subroutine s_mpi_decompose_computational_domain(write_silo_ghost_offsets, adjust_local_domains, output_offsets, local_domains) + + logical, intent(in) :: write_silo_ghost_offsets + logical, intent(in) :: adjust_local_domains + type(int_bounds_info), dimension(3), intent(inout), optional :: output_offsets + type(bounds_info), dimension(3), intent(inout), optional :: local_domains #ifdef MFC_MPI !> Non-optimal number of processors in the x-, y- and z-directions @@ -1322,6 +1329,7 @@ contains real(wp) :: fct_min !< Processor factorization (fct) minimization parameter integer :: MPI_COMM_CART !< Cartesian processor topology communicator integer :: rem_cells !< Remaining cells after distribution among processors + integer :: rem_cells_by_dim(3) integer :: recon_order !< WENO or MUSCL reconstruction order integer :: i, j, k !< Generic loop iterators integer :: ierr !< Generic flag used to identify and report MPI errors @@ -1329,6 +1337,8 @@ contains ! temp array to store neighbor rank coordinates integer, dimension(1:num_dims) :: neighbor_coords + rem_cells_by_dim = 0 + ! Zeroing out communication needs for moving EL bubbles/particles nidx(1)%beg = 0; nidx(1)%end = 0 nidx(2)%beg = 0; nidx(2)%end = 0 @@ -1470,6 +1480,7 @@ contains ! Number of remaining cells rem_cells = mod(p + 1, num_procs_z) + rem_cells_by_dim(3) = rem_cells ! Optimal number of cells per processor p = (p + 1)/num_procs_z - 1 @@ -1497,22 +1508,6 @@ contains nidx(3)%end = 1 end if -#ifdef MFC_POST_PROCESS - ! Ghost zone at the beginning - if (proc_coords(3) > 0 .and. format == format_silo) then - offset_z%beg = 2 - else - offset_z%beg = 0 - end if - - ! Ghost zone at the end - if (proc_coords(3) < num_procs_z - 1 .and. format == format_silo) then - offset_z%end = 2 - else - offset_z%end = 0 - end if -#endif - ! Beginning and end sub-domain boundary locations if (parallel_io) then if (proc_coords(3) < rem_cells) then @@ -1520,21 +1515,6 @@ contains else start_idx(3) = (p + 1)*proc_coords(3) + rem_cells end if - else -#ifdef MFC_PRE_PROCESS - if (old_grid .neqv. .true.) then - dz_min = (z_domain%end - z_domain%beg)/real(p_glb + 1, wp) - - if (proc_coords(3) < rem_cells) then - z_domain%beg = z_domain%beg + dz_min*real((p + 1)*proc_coords(3)) - z_domain%end = z_domain%end - dz_min*real((p + 1)*(num_procs_z - proc_coords(3) - 1) - (num_procs_z & - & - rem_cells)) - else - z_domain%beg = z_domain%beg + dz_min*real((p + 1)*proc_coords(3) + rem_cells) - z_domain%end = z_domain%end - dz_min*real((p + 1)*(num_procs_z - proc_coords(3) - 1)) - end if - end if -#endif end if ! 2D Cartesian Processor Topology @@ -1584,6 +1564,7 @@ contains ! Number of remaining cells rem_cells = mod(n + 1, num_procs_y) + rem_cells_by_dim(2) = rem_cells ! Optimal number of cells per processor n = (n + 1)/num_procs_y - 1 @@ -1611,22 +1592,6 @@ contains nidx(2)%end = 1 end if -#ifdef MFC_POST_PROCESS - ! Ghost zone at the beginning - if (proc_coords(2) > 0 .and. format == format_silo) then - offset_y%beg = 2 - else - offset_y%beg = 0 - end if - - ! Ghost zone at the end - if (proc_coords(2) < num_procs_y - 1 .and. format == format_silo) then - offset_y%end = 2 - else - offset_y%end = 0 - end if -#endif - ! Beginning and end sub-domain boundary locations if (parallel_io) then if (proc_coords(2) < rem_cells) then @@ -1634,21 +1599,6 @@ contains else start_idx(2) = (n + 1)*proc_coords(2) + rem_cells end if - else -#ifdef MFC_PRE_PROCESS - if (old_grid .neqv. .true.) then - dy_min = (y_domain%end - y_domain%beg)/real(n_glb + 1, wp) - - if (proc_coords(2) < rem_cells) then - y_domain%beg = y_domain%beg + dy_min*real((n + 1)*proc_coords(2)) - y_domain%end = y_domain%end - dy_min*real((n + 1)*(num_procs_y - proc_coords(2) - 1) - (num_procs_y & - & - rem_cells)) - else - y_domain%beg = y_domain%beg + dy_min*real((n + 1)*proc_coords(2) + rem_cells) - y_domain%end = y_domain%end - dy_min*real((n + 1)*(num_procs_y - proc_coords(2) - 1)) - end if - end if -#endif end if ! 1D Cartesian Processor Topology @@ -1667,6 +1617,7 @@ contains ! Number of remaining cells rem_cells = mod(m + 1, num_procs_x) + rem_cells_by_dim(1) = rem_cells ! Optimal number of cells per processor m = (m + 1)/num_procs_x - 1 @@ -1696,22 +1647,6 @@ contains nidx(1)%end = 1 end if -#ifdef MFC_POST_PROCESS - ! Ghost zone at the beginning - if (proc_coords(1) > 0 .and. format == format_silo) then - offset_x%beg = 2 - else - offset_x%beg = 0 - end if - - ! Ghost zone at the end - if (proc_coords(1) < num_procs_x - 1 .and. format == format_silo) then - offset_x%end = 2 - else - offset_x%end = 0 - end if -#endif - ! Beginning and end sub-domain boundary locations if (parallel_io) then if (proc_coords(1) < rem_cells) then @@ -1719,23 +1654,12 @@ contains else start_idx(1) = (m + 1)*proc_coords(1) + rem_cells end if - else -#ifdef MFC_PRE_PROCESS - if (old_grid .neqv. .true.) then - dx_min = (x_domain%end - x_domain%beg)/real(m_glb + 1, wp) - - if (proc_coords(1) < rem_cells) then - x_domain%beg = x_domain%beg + dx_min*real((m + 1)*proc_coords(1)) - x_domain%end = x_domain%end - dx_min*real((m + 1)*(num_procs_x - proc_coords(1) - 1) - (num_procs_x & - & - rem_cells)) - else - x_domain%beg = x_domain%beg + dx_min*real((m + 1)*proc_coords(1) + rem_cells) - x_domain%end = x_domain%end - dx_min*real((m + 1)*(num_procs_x - proc_coords(1) - 1)) - end if - end if -#endif end if + call s_apply_decomposition_policies((/num_procs_x, num_procs_y, num_procs_z/), rem_cells_by_dim, (/m, n, p/), (/m_glb, & + & n_glb, p_glb/), write_silo_ghost_offsets, & + & adjust_local_domains .and. (.not. parallel_io), output_offsets, local_domains) + @:ALLOCATE(neighbor_ranks(nidx(1)%beg:nidx(1)%end, nidx(2)%beg:nidx(2)%end, nidx(3)%beg:nidx(3)%end)) do k = nidx(3)%beg, nidx(3)%end do j = nidx(2)%beg, nidx(2)%end @@ -1753,6 +1677,68 @@ contains end subroutine s_mpi_decompose_computational_domain + !> Apply executable-configured output and local-domain policies after the shared Cartesian decomposition. + subroutine s_apply_decomposition_policies(proc_counts, remainders, local_cells, global_cells, write_silo_ghost_offsets, & + & adjust_local_domains, output_offsets, local_domains) + + integer, dimension(3), intent(in) :: proc_counts, remainders, local_cells, global_cells + logical, intent(in) :: write_silo_ghost_offsets, adjust_local_domains + type(int_bounds_info), dimension(3), intent(inout), optional :: output_offsets + type(bounds_info), dimension(3), intent(inout), optional :: local_domains + integer :: dim + real(wp) :: domain_beg, domain_end, spacing + + if (write_silo_ghost_offsets .and. .not. present(output_offsets)) then + call s_mpi_abort('Silo ghost-offset policy requires output offset storage.') + end if + if (adjust_local_domains .and. .not. present(local_domains)) then + call s_mpi_abort('Local-domain adjustment policy requires domain storage.') + end if + + if (present(output_offsets)) then + do dim = 1, 3 + output_offsets(dim)%beg = 0 + output_offsets(dim)%end = 0 + end do + + if (write_silo_ghost_offsets) then + do dim = 1, num_dims + if (proc_coords(dim) > 0) output_offsets(dim)%beg = 2 + if (proc_coords(dim) < proc_counts(dim) - 1) output_offsets(dim)%end = 2 + end do + end if + end if + + if (adjust_local_domains) then + do dim = 1, num_dims + domain_beg = local_domains(dim)%beg + domain_end = local_domains(dim)%end + spacing = (domain_end - domain_beg)/real(global_cells(dim) + 1, wp) + + select case (dim) + case (1) + dx_min = spacing + case (2) + dy_min = spacing + case (3) + dz_min = spacing + end select + + if (proc_coords(dim) < remainders(dim)) then + local_domains(dim)%beg = domain_beg + spacing*real((local_cells(dim) + 1)*proc_coords(dim), wp) + local_domains(dim)%end = domain_end - spacing*real((local_cells(dim) + 1)*(proc_counts(dim) - proc_coords(dim) & + & - 1) - (proc_counts(dim) - remainders(dim)), wp) + else + local_domains(dim)%beg = domain_beg + spacing*real((local_cells(dim) + 1)*proc_coords(dim) + remainders(dim), & + & wp) + local_domains(dim)%end = domain_end - spacing*real((local_cells(dim) + 1)*(proc_counts(dim) - proc_coords(dim) & + & - 1), wp) + end if + end do + end if + + end subroutine s_apply_decomposition_policies + !> The goal of this procedure is to populate the buffers of the grid variables by communicating with the neighboring processors. !! Note that only the buffers of the cell-width distributions are handled in such a way. This is because the buffers of !! cell-boundary locations may be calculated directly from those of the cell-width distributions. diff --git a/src/post_process/m_start_up.fpp b/src/post_process/m_start_up.fpp index 707d5b117a..a18c113d2b 100644 --- a/src/post_process/m_start_up.fpp +++ b/src/post_process/m_start_up.fpp @@ -26,7 +26,7 @@ module m_start_up use m_checker use m_thermochem, only: num_species, species_names use m_finite_differences - use m_constants, only: model_eqns_gamma_law, model_eqns_5eq, model_eqns_6eq, model_eqns_4eq + use m_constants, only: model_eqns_gamma_law, model_eqns_5eq, model_eqns_6eq, model_eqns_4eq, format_silo use m_chemistry #ifdef MFC_MPI @@ -932,6 +932,8 @@ contains !> Set up the MPI environment, read and broadcast user inputs, and decompose the computational domain. impure subroutine s_initialize_mpi_domain + type(int_bounds_info), dimension(3) :: output_offsets + num_dims = 1 + min(1, n) + min(1, p) call s_mpi_initialize() @@ -946,7 +948,12 @@ contains call s_mpi_bcast_user_inputs() call s_initialize_parallel_io() - call s_mpi_decompose_computational_domain() + output_offsets = (/offset_x, offset_y, offset_z/) + call s_mpi_decompose_computational_domain(write_silo_ghost_offsets=format == format_silo, adjust_local_domains=.false., & + & output_offsets=output_offsets) + offset_x = output_offsets(1) + offset_y = output_offsets(2) + offset_z = output_offsets(3) call s_check_inputs_fft() bc = bc_xyz_info(bc_x, bc_y, bc_z) diff --git a/src/pre_process/m_start_up.fpp b/src/pre_process/m_start_up.fpp index cc60ea5fe4..04e3997378 100644 --- a/src/pre_process/m_start_up.fpp +++ b/src/pre_process/m_start_up.fpp @@ -608,6 +608,8 @@ contains !> Initialize MPI, read and validate user inputs on rank 0, and decompose the computational domain. impure subroutine s_initialize_mpi_domain + type(bounds_info), dimension(3) :: local_domains + call s_mpi_initialize() if (proc_rank == 0) then @@ -628,7 +630,12 @@ contains y_domain_glb = y_domain z_domain_glb = z_domain - call s_mpi_decompose_computational_domain() + local_domains = (/x_domain, y_domain, z_domain/) + call s_mpi_decompose_computational_domain(write_silo_ghost_offsets=.false., adjust_local_domains= .not. old_grid, & + & local_domains=local_domains) + x_domain = local_domains(1) + y_domain = local_domains(2) + z_domain = local_domains(3) bc = bc_xyz_info(bc_x, bc_y, bc_z) diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index c39f52bc80..7f20004800 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -1018,7 +1018,7 @@ contains call s_initialize_parallel_io() - call s_mpi_decompose_computational_domain() + call s_mpi_decompose_computational_domain(write_silo_ghost_offsets=.false., adjust_local_domains=.false.) bc = bc_xyz_info(bc_x, bc_y, bc_z) From bf43fd2c64be216d8e117bc45bafe675eafad5f2 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Sat, 25 Jul 2026 22:52:06 -0400 Subject: [PATCH 10/15] Make STL model setup stage-independent --- src/common/m_model.fpp | 4 ---- toolchain/mfc/params/definitions.py | 2 +- toolchain/mfc/params/generators/fortran_gen.py | 3 ++- toolchain/mfc/params_tests/test_fortran_gen.py | 12 ++++++++++++ 4 files changed, 15 insertions(+), 6 deletions(-) diff --git a/src/common/m_model.fpp b/src/common/m_model.fpp index 0183aab55f..7fa14cdfa4 100644 --- a/src/common/m_model.fpp +++ b/src/common/m_model.fpp @@ -24,9 +24,7 @@ module m_model public :: s_check_boundary, s_register_edge, f_model_is_inside, s_distance_normals_3D, s_distance_normals_2D, & & s_pack_model_for_gpu -#ifndef MFC_POST_PROCESS public :: s_instantiate_STL_models -#endif type(t_model_array), allocatable, target :: models(:) !< STL/OBJ models for IB markers and levelset integer, allocatable :: gpu_ntrs(:) !< GPU-friendly flat arrays for STL model data @@ -851,7 +849,6 @@ contains end subroutine s_distance_normals_2D -#ifndef MFC_POST_PROCESS !> Load, transform, and register STL/OBJ immersed-boundary models onto the simulation grid. subroutine s_instantiate_STL_models() @@ -1002,7 +999,6 @@ contains end block end subroutine s_instantiate_STL_models -#endif !> Pack triangle vertices and normals from a model into flat arrays for GPU transfer. subroutine s_pack_model_for_gpu(ma) diff --git a/toolchain/mfc/params/definitions.py b/toolchain/mfc/params/definitions.py index 93fd27263a..54cdef6ca2 100644 --- a/toolchain/mfc/params/definitions.py +++ b/toolchain/mfc/params/definitions.py @@ -1223,7 +1223,7 @@ def _decl(targets: set, *names: str) -> None: _PRE = {"pre"} _POST = {"post"} -_decl(_ALL, "avg_state", "alt_soundspeed", "mixture_err", "sigR", "viscous", "riemann_solver") +_decl(_ALL, "avg_state", "alt_soundspeed", "mixture_err", "sigR", "viscous", "riemann_solver", "stl_models", "num_stl_models") _nv( _ALL, diff --git a/toolchain/mfc/params/generators/fortran_gen.py b/toolchain/mfc/params/generators/fortran_gen.py index a77d0d42d8..f4f36d6a89 100644 --- a/toolchain/mfc/params/generators/fortran_gen.py +++ b/toolchain/mfc/params/generators/fortran_gen.py @@ -234,7 +234,8 @@ def generate_decls_fpp(target: str) -> str: lines.append(f"{fortran_type_decl(param).ljust(_DECL_COL)}:: {name}") declared_names.add(name) for name, (ftype, dim, gpu, desc) in TYPED_DECLS.items(): - if name not in NAMELIST_VARS or target not in NAMELIST_VARS[name]: + declaration_targets = NAMELIST_VARS.get(name, set()) | DECLARATION_TARGETS.get(name, set()) + if target not in declaration_targets: continue decl = f"{ftype}, dimension({dim})" if dim else ftype padded = decl.ljust(_ARRAY_DECL_COL) diff --git a/toolchain/mfc/params_tests/test_fortran_gen.py b/toolchain/mfc/params_tests/test_fortran_gen.py index 9752e6bdd9..7ff5d8b493 100644 --- a/toolchain/mfc/params_tests/test_fortran_gen.py +++ b/toolchain/mfc/params_tests/test_fortran_gen.py @@ -142,6 +142,18 @@ def test_common_module_decls_do_not_expand_namelists(): assert NAMELIST_VARS["riemann_solver"] == {"sim"} +def test_stl_model_declarations_do_not_expand_post_namelist(): + from mfc.params.generators.fortran_gen import generate_decls_fpp, generate_namelist_fpp + + post_declarations = generate_decls_fpp("post") + post_namelist = generate_namelist_fpp("post") + + assert ":: num_stl_models" in post_declarations + assert ":: stl_models" in post_declarations + assert "num_stl_models" not in post_namelist + assert "stl_models" not in post_namelist + + def test_decls_no_percent_vars(): from mfc.params.generators.fortran_gen import generate_decls_fpp From 395b0224d1bf723d17a58cd7b94e3b58a1e312ce Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Sat, 25 Jul 2026 23:00:02 -0400 Subject: [PATCH 11/15] Make phase change stage-independent --- src/common/m_phase_change.fpp | 3 +-- toolchain/mfc/params/definitions.py | 14 +++++++++++++- toolchain/mfc/params_tests/test_fortran_gen.py | 11 +++++++++++ 3 files changed, 25 insertions(+), 3 deletions(-) diff --git a/src/common/m_phase_change.fpp b/src/common/m_phase_change.fpp index e5b601331a..b56c9a6ecf 100644 --- a/src/common/m_phase_change.fpp +++ b/src/common/m_phase_change.fpp @@ -8,7 +8,6 @@ !> @brief Phase transition relaxation solvers for liquid-vapor flows with cavitation and boiling module m_phase_change -#ifndef MFC_POST_PROCESS use m_derived_types use m_global_parameters use m_mpi_proxy @@ -470,5 +469,5 @@ contains impure subroutine s_finalize_relaxation_solver_module end subroutine s_finalize_relaxation_solver_module -#endif + end module m_phase_change diff --git a/toolchain/mfc/params/definitions.py b/toolchain/mfc/params/definitions.py index 54cdef6ca2..649c81fc0f 100644 --- a/toolchain/mfc/params/definitions.py +++ b/toolchain/mfc/params/definitions.py @@ -1223,7 +1223,19 @@ def _decl(targets: set, *names: str) -> None: _PRE = {"pre"} _POST = {"post"} -_decl(_ALL, "avg_state", "alt_soundspeed", "mixture_err", "sigR", "viscous", "riemann_solver", "stl_models", "num_stl_models") +_decl( + _ALL, + "avg_state", + "alt_soundspeed", + "mixture_err", + "sigR", + "viscous", + "riemann_solver", + "stl_models", + "num_stl_models", + "palpha_eps", + "ptgalpha_eps", +) _nv( _ALL, diff --git a/toolchain/mfc/params_tests/test_fortran_gen.py b/toolchain/mfc/params_tests/test_fortran_gen.py index 7ff5d8b493..a57175514f 100644 --- a/toolchain/mfc/params_tests/test_fortran_gen.py +++ b/toolchain/mfc/params_tests/test_fortran_gen.py @@ -154,6 +154,17 @@ def test_stl_model_declarations_do_not_expand_post_namelist(): assert "stl_models" not in post_namelist +def test_phase_change_declarations_do_not_expand_post_namelist(): + from mfc.params.generators.fortran_gen import generate_decls_fpp, generate_namelist_fpp + + post_declarations = generate_decls_fpp("post") + post_namelist = generate_namelist_fpp("post") + + for name in ("palpha_eps", "ptgalpha_eps"): + assert f":: {name}" in post_declarations + assert name not in post_namelist + + def test_decls_no_percent_vars(): from mfc.params.generators.fortran_gen import generate_decls_fpp From d875530f5681876b3cdad980a16f00f8dd627995 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Sun, 26 Jul 2026 09:52:43 -0400 Subject: [PATCH 12/15] Restore mhd GPU declaration --- toolchain/mfc/params/generators/fortran_gen.py | 1 + 1 file changed, 1 insertion(+) diff --git a/toolchain/mfc/params/generators/fortran_gen.py b/toolchain/mfc/params/generators/fortran_gen.py index f4f36d6a89..6c2a1ad97c 100644 --- a/toolchain/mfc/params/generators/fortran_gen.py +++ b/toolchain/mfc/params/generators/fortran_gen.py @@ -75,6 +75,7 @@ "mixture_err", "model_eqns", "mp_weno", + "mhd", "mpp_lim", "muscl_eps", "muscl_lim", From 5ae5c6804c1e6ead244c30a11e65998e252a4ab5 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Sun, 26 Jul 2026 13:45:11 -0400 Subject: [PATCH 13/15] Simplify conversion dispatch and guard generated GPU declarations Collapse the optional-argument dispatch in s_convert_species_to_mixture_variables: absent optional dummies forward as absent, so the four branches were unnecessary. Default dirichlet_from_buffers so it is defined before module initialization. Fail generation when SIM_GPU_DECL_VARS names a variable simulation does not declare, instead of silently dropping its device residency. Drop a stray build log and refresh the stage-independence notes. --- .claude/rules/common-pitfalls.md | 14 +++-- build_merge.log | 62 ------------------- src/common/m_boundary_primitives.fpp | 2 +- src/common/m_variables_conversion.fpp | 13 +--- .../mfc/params/generators/fortran_gen.py | 9 +++ 5 files changed, 20 insertions(+), 80 deletions(-) delete mode 100644 build_merge.log diff --git a/.claude/rules/common-pitfalls.md b/.claude/rules/common-pitfalls.md index c4838f4c62..71bf6e452f 100644 --- a/.claude/rules/common-pitfalls.md +++ b/.claude/rules/common-pitfalls.md @@ -77,15 +77,17 @@ covered in `docs/documentation/contributing.md`. garbage `v_blow` failed Frontier AMD with `ICFL is NaN` while every NVIDIA lane and all local CPU/GPU runs passed. A platform-only NaN is the signature of this class. - Shared-state pattern: namelist declarations (`#:include 'generated_decls.fpp'`), the - `eqn_idx`/`sys_size`/`b_size`/`tensor_size` state variables, and the common defaults + `eqn_idx`/`sys_size` state variables, and the common defaults core all live in `src/common/m_global_parameters_common.fpp`. Each per-target `m_global_parameters.fpp` does `use m_global_parameters_common` (default-public), so `use m_global_parameters` continues to work for all downstream modules without change. - Sim-only declarations (GPU_DECLARE, Re_idx allocation) stay in - `m_global_parameters_common` behind `#ifdef MFC_SIMULATION`. Generated includes - (`generated_decls.fpp`, `generated_bcast.fpp`, `generated_case_opt_decls.fpp`) must exist for every target — the build - emits stubs where the content is sim-only, so a common file that includes one will - compile for pre/post too. + `src/common/` carries no `MFC_PRE_PROCESS`/`MFC_SIMULATION`/`MFC_POST_PROCESS` guards: + stage-varying behavior is passed in as an explicit argument or initialization policy, and + device residency for generated simulation scalars is emitted from `SIM_GPU_DECL_VARS` + (`toolchain/mfc/params/generators/fortran_gen.py`). Generated includes + (`generated_decls.fpp`, `generated_bcast.fpp`, `generated_case_opt_decls.fpp`) must exist for every target — pre/post + get the common computed scalars (`num_dims`, `num_vels`, `weno_polyn`, `muscl_polyn`), so + a common file that includes one will compile for pre/post too. - Runtime checks (`@:PROHIBIT`) go where they run: shared → `src/common/m_checker_common.fpp`; simulation-only → `src/simulation/m_checker.fpp`; pre/post-only → `src/{pre,post}_process/m_checker.fpp` (their `s_check_inputs` are diff --git a/build_merge.log b/build_merge.log deleted file mode 100644 index 16b4b538b2..0000000000 --- a/build_merge.log +++ /dev/null @@ -1,62 +0,0 @@ -mfc: Starting... -mfc: OK > Created a Python virtual environment (venv). -mfc: OK > (venv) Entered the Python 3.12.3 virtual environment (>= 3.9). -mfc: (venv) Installing Python packages... -mfc: (venv) Installing uv package manager for fast installation... -mfc: OK > (venv) Installed uv. -mfc: (venv) Using uv for fast installation... -Using Python 3.12.3 environment at: build/venv -Resolved 96 packages in 2.90s - Building mfc @ file:///fastscratch/sbryngelson3/mfc-pr1679-merge/toolchain - Built mfc @ file:///fastscratch/sbryngelson3/mfc-pr1679-merge/toolchain -Prepared 1 package in 464ms -Installed 96 packages in 1.69s -mfc: OK > (venv) Installation succeeded. -mfc: OK > (venv) Applied fypp linemarker-resync patch. - - Auto-regenerating completion scripts... - Tab completions updated. Run: source /nethome/sbryngelson3/.local/share/mfc/completions/mfc.bash - .=++*: -+*+=. | sbryngelson3@wingtip-gpu3.cc.gatech.edu [Linux] - :+ -*- == =* . | ----------------------------------------------- - :*+ == ++ .+- | - :*##-.....:*+ .#%+++=--+=:::. | --jobs 24 - -=-++-======#=--**+++==+*++=::-:. | --mpi --no-gpu --no-debug --no-reldebug --no-gcov --no-unified --no-single --no-mixed --no-fastmath - .:++=----------====+*= ==..:%..... | --targets pre_process, simulation, and post_process - .:-=++++===--==+=-+= +. := | - +#=::::::::=%=. -+: =+ *: | ------------------------------------------------------------------ - .*=-=*=.. :=+*+: -...-- | $ ./mfc.sh (build, run, test, clean, new, validate, params) --help - - Build | syscheck, pre_process, simulation, and post_process - - Generating case.fpp. - Writing a (new) custom case.fpp file. - Configuring syscheck... - ✓ Configured syscheck - - Generating case.fpp. - Writing a (new) custom case.fpp file. - Building syscheck... - (build took 8.1s) - ✓ Built syscheck - - Installing syscheck... - ✓ Installed syscheck - - Generating case.fpp. - Writing a (new) custom case.fpp file. - Configuring pre_process... - ✓ Configured pre_process - - Generating case.fpp. - Writing a (new) custom case.fpp file. - Building pre_process... - (build took 55.7s) - ✓ Built pre_process - - Installing pre_process... - ✓ Installed pre_process - - Configuring fftw... - ✓ Configured fftw - - Building fftw... diff --git a/src/common/m_boundary_primitives.fpp b/src/common/m_boundary_primitives.fpp index b947b3b274..90deaf09a7 100644 --- a/src/common/m_boundary_primitives.fpp +++ b/src/common/m_boundary_primitives.fpp @@ -18,7 +18,7 @@ module m_boundary_primitives type(scalar_field), dimension(:,:), allocatable :: bc_buffers $:GPU_DECLARE(create='[bc_buffers]') - logical :: dirichlet_from_buffers + logical :: dirichlet_from_buffers = .false. $:GPU_DECLARE(create='[dirichlet_from_buffers]') contains diff --git a/src/common/m_variables_conversion.fpp b/src/common/m_variables_conversion.fpp index 61c7328486..695bc89a2c 100644 --- a/src/common/m_variables_conversion.fpp +++ b/src/common/m_variables_conversion.fpp @@ -184,17 +184,8 @@ contains call s_compute_species_fraction(q_vf, k, l, r, alpha_rho_K, alpha_K) ! Use the same scalar kernel on host and device so mixture semantics do not depend on the executable or accelerator backend. - if (present(G_K)) then - if (present(Re_K)) then - call s_convert_species_to_mixture_variables_kernel(rho, gamma, pi_inf, qv, alpha_K, alpha_rho_K, Re_K, G_K, G) - else - call s_convert_species_to_mixture_variables_kernel(rho, gamma, pi_inf, qv, alpha_K, alpha_rho_K, G_K=G_K, G=G) - end if - else if (present(Re_K)) then - call s_convert_species_to_mixture_variables_kernel(rho, gamma, pi_inf, qv, alpha_K, alpha_rho_K, Re_K) - else - call s_convert_species_to_mixture_variables_kernel(rho, gamma, pi_inf, qv, alpha_K, alpha_rho_K) - end if + ! Absent optional dummies forward as absent, so the optional arguments need no dispatch here. + call s_convert_species_to_mixture_variables_kernel(rho, gamma, pi_inf, qv, alpha_K, alpha_rho_K, Re_K, G_K, G) ! Store derived mixture fields when requested during module initialization. if (allocated(rho_sf)) then diff --git a/toolchain/mfc/params/generators/fortran_gen.py b/toolchain/mfc/params/generators/fortran_gen.py index 8debc46881..0c8aba9687 100644 --- a/toolchain/mfc/params/generators/fortran_gen.py +++ b/toolchain/mfc/params/generators/fortran_gen.py @@ -695,10 +695,19 @@ def get_generated_files(build_dir: Path) -> List[Tuple[Path, str]]: result.append((inc / "generated_namelist.fpp", generate_namelist_fpp(short))) result.append((inc / "generated_decls.fpp", generate_decls_fpp(short))) result.append((inc / "generated_constants.fpp", generate_constants_fpp())) + sim_gpu_decls = "" for short, full in TARGETS: inc = build_dir / "include" / full content = generate_case_opt_decls_fpp() if short == "sim" else generate_common_extra_decls_fpp() + if short == "sim": + sim_gpu_decls = content result.append((inc / "generated_case_opt_decls.fpp", content)) + # A name here that simulation never declares emits no GPU_DECLARE, so the variable + # silently loses device residency. Fail on the stale entry instead. + sim_gpu_decls += generate_decls_fpp("sim") + stale = sorted(n for n in SIM_GPU_DECL_VARS if f"$:GPU_DECLARE(create='[{n}]')" not in sim_gpu_decls) + if stale: + raise ValueError(f"SIM_GPU_DECL_VARS names that simulation does not declare: {', '.join(stale)}. Remove them or restore the parameter.") for short, full in TARGETS: inc = build_dir / "include" / full result.append((inc / "generated_bcast.fpp", generate_bcast_fpp(short))) From bad33421e700478aae7e596c6bb9eeff535d2798 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Sun, 26 Jul 2026 14:07:36 -0400 Subject: [PATCH 14/15] Fail generation on an orphaned common case-optimization declaration COMMON_CASE_OPT_EXTRA_NAMES silently emitted nothing for a name that CASE_OPT_EXTRA_LINES does not define, leaving pre_process and post_process without a declaration that common code references. Raise instead, matching the SIM_GPU_DECL_VARS guard. --- toolchain/mfc/params/generators/fortran_gen.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/toolchain/mfc/params/generators/fortran_gen.py b/toolchain/mfc/params/generators/fortran_gen.py index 0c8aba9687..8209334e2b 100644 --- a/toolchain/mfc/params/generators/fortran_gen.py +++ b/toolchain/mfc/params/generators/fortran_gen.py @@ -341,6 +341,11 @@ def generate_case_opt_decls_fpp() -> str: def generate_common_extra_decls_fpp() -> str: """Return computed-scalar declarations needed by common pre/post code.""" + # A name here that CASE_OPT_EXTRA_LINES does not define emits nothing, leaving pre/post + # without a declaration that common code references. Fail on the orphan instead. + orphans = sorted(COMMON_CASE_OPT_EXTRA_NAMES - {name for name, _, _ in CASE_OPT_EXTRA_LINES}) + if orphans: + raise ValueError(f"COMMON_CASE_OPT_EXTRA_NAMES entries missing from CASE_OPT_EXTRA_LINES: {', '.join(orphans)}.") lines = [_HEADER.rstrip()] for name, ftype, _ in CASE_OPT_EXTRA_LINES: if name in COMMON_CASE_OPT_EXTRA_NAMES: From 9599c5bc7d64cc8fb140762952f92c8ad757fae1 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Sun, 26 Jul 2026 17:01:49 -0400 Subject: [PATCH 15/15] Fix case-optimization riemann_solver default and undefined reads riemann_solver is not a case-optimization parameter, so its default belongs outside the case-optimization guard; case-optimized simulation builds otherwise lost the dflt_int sentinel. post_process declares viscous for the common conversion kernel but neither accepts it as input nor broadcasts it, and assigns defaults on rank 0 only, so give it a value on every rank. Initialize the MPI transport policies, keep global_bounds intent(inout) so inactive dimensions retain their defaults, and mark the grid-width buffer contiguous where an element is passed to MPI. --- src/common/m_boundary_common.fpp | 14 +++++++------- src/common/m_global_parameters_common.fpp | 4 +++- src/common/m_mpi_common.fpp | 8 +++++--- src/post_process/m_global_parameters.fpp | 4 ++++ 4 files changed, 19 insertions(+), 11 deletions(-) diff --git a/src/common/m_boundary_common.fpp b/src/common/m_boundary_common.fpp index 3dfa30528b..06681e0d50 100644 --- a/src/common/m_boundary_common.fpp +++ b/src/common/m_boundary_common.fpp @@ -350,13 +350,13 @@ contains subroutine s_populate_grid_variables_buffers(x_cb_in, x_cc_in, dx_in, x_offset, y_offset, z_offset, y_cb_in, y_cc_in, dy_in, & & z_cb_in, z_cc_in, dz_in, global_bounds) - type(int_bounds_info), intent(in) :: x_offset, y_offset, z_offset - real(wp), intent(inout) :: x_cb_in(-1 - x_offset%beg:) - real(wp), intent(inout) :: x_cc_in(-buff_size:), dx_in(-buff_size:) - real(wp), optional, intent(inout) :: y_cb_in(-1 - y_offset%beg:), z_cb_in(-1 - z_offset%beg:) - real(wp), optional, intent(inout) :: y_cc_in(-buff_size:), dy_in(-buff_size:) - real(wp), optional, intent(inout) :: z_cc_in(-buff_size:), dz_in(-buff_size:) - type(bounds_info), optional, dimension(3), intent(out) :: global_bounds + type(int_bounds_info), intent(in) :: x_offset, y_offset, z_offset + real(wp), intent(inout) :: x_cb_in(-1 - x_offset%beg:) + real(wp), intent(inout) :: x_cc_in(-buff_size:), dx_in(-buff_size:) + real(wp), optional, intent(inout) :: y_cb_in(-1 - y_offset%beg:), z_cb_in(-1 - z_offset%beg:) + real(wp), optional, intent(inout) :: y_cc_in(-buff_size:), dy_in(-buff_size:) + real(wp), optional, intent(inout) :: z_cc_in(-buff_size:), dz_in(-buff_size:) + type(bounds_info), optional, dimension(3), intent(inout) :: global_bounds if (present(global_bounds)) then #ifdef MFC_MPI diff --git a/src/common/m_global_parameters_common.fpp b/src/common/m_global_parameters_common.fpp index 88d0dcd9f6..b43bbddecc 100644 --- a/src/common/m_global_parameters_common.fpp +++ b/src/common/m_global_parameters_common.fpp @@ -376,9 +376,11 @@ contains mhd = .false. relativity = .false. viscous = .false. - riemann_solver = dflt_int #:endif + ! Not a case-optimization parameter, so it stays a runtime variable in every build. + riemann_solver = dflt_int + ! Tait EOS rhoref = dflt_real pref = dflt_real diff --git a/src/common/m_mpi_common.fpp b/src/common/m_mpi_common.fpp index 7220823906..1154b68f2e 100644 --- a/src/common/m_mpi_common.fpp +++ b/src/common/m_mpi_common.fpp @@ -42,8 +42,8 @@ module m_mpi_common integer(kind=8) :: halo_size $:GPU_DECLARE(create='[halo_size]') - logical, private :: exchange_all_chemistry_temperatures - logical, private :: use_rdma_transport + logical, private :: exchange_all_chemistry_temperatures = .false. + logical, private :: use_rdma_transport = .false. contains @@ -1749,7 +1749,9 @@ contains type(int_bounds_info), intent(in) :: bc_bounds, offset real(wp), intent(inout) :: cell_boundaries(-1 - offset%beg:) real(wp), intent(inout) :: cell_centers(-buff_size:) - real(wp), intent(inout) :: cell_widths(-buff_size:) + ! Contiguous so that passing an element to MPI is a plain address, with no descriptor + ! or copy-in/copy-out. Every actual argument is a whole module array. + real(wp), contiguous, intent(inout) :: cell_widths(-buff_size:) #ifdef MFC_MPI integer :: ierr diff --git a/src/post_process/m_global_parameters.fpp b/src/post_process/m_global_parameters.fpp index 3ec6eb98ca..1423c62feb 100644 --- a/src/post_process/m_global_parameters.fpp +++ b/src/post_process/m_global_parameters.fpp @@ -318,6 +318,10 @@ contains if (n == 0) m_root = m_glb + ! Declared for the common conversion kernel but not a post_process input, so it is neither + ! defaulted on non-root ranks nor broadcast. Post-process never carries viscous stresses. + viscous = .false. + ! Gamma/Pi_inf: force num_fluids=1 (post_process-specific side effect of the gamma-law model) if (model_eqns == model_eqns_gamma_law) num_fluids = 1