From 40af4f1dbeaee7366dbd6227456a991c8d75102e Mon Sep 17 00:00:00 2001 From: "Daniel J. Vickers" Date: Thu, 10 Sep 2026 10:16:20 -0400 Subject: [PATCH 1/5] Initial move of code from simulation to pre_process --- src/pre_process/m_data_output.fpp | 46 ++++- src/pre_process/m_global_parameters.fpp | 21 ++ .../m_particle_cloud.fpp | 88 +++----- src/pre_process/m_start_up.fpp | 17 +- src/simulation/m_start_up.fpp | 193 +++++++++++++----- toolchain/mfc/params/definitions.py | 3 + toolchain/mfc/test/cases.py | 32 +++ 7 files changed, 286 insertions(+), 114 deletions(-) rename src/{simulation => pre_process}/m_particle_cloud.fpp (83%) diff --git a/src/pre_process/m_data_output.fpp b/src/pre_process/m_data_output.fpp index e6b4847a06..f8f8c59e41 100644 --- a/src/pre_process/m_data_output.fpp +++ b/src/pre_process/m_data_output.fpp @@ -28,7 +28,7 @@ module m_data_output private public :: s_write_serial_data_files, s_write_parallel_data_files, s_write_data_files, s_initialize_data_output_module, & - & s_finalize_data_output_module + & s_finalize_data_output_module, s_write_ib_state_0 type(scalar_field), allocatable, dimension(:) :: q_cons_temp @@ -733,6 +733,50 @@ contains end subroutine s_initialize_data_output_module + !> @brief Writes restart_data/ib_state_0.dat: the initial IB layout (namelist patch_ib entries, then any generated + !! particle-cloud beds). Rank-0-only. Read back by simulation at startup via s_read_ib_restart_data(0, ...) + !! (src/simulation/m_start_up.fpp). Uses the same 20-field record layout s_write_serial_ib_state + !! (src/simulation/m_data_output.fpp) writes, so simulation's own writers can freely overwrite this file later at t_step_start + !! == 0 without changing its layout - only position (fields 17:19) and radius (field 20) are populated here; everything else + !! (time, force, torque, vel, angular_vel, angles) is zero for a freshly generated IB. + impure subroutine s_write_ib_state_0(particle_cloud_ibs, num_particle_cloud_ibs) + + type(ib_patch_parameters), dimension(:), intent(in) :: particle_cloud_ibs + integer, intent(in) :: num_particle_cloud_ibs + character(LEN=len_trim(case_dir) + 2*name_len) :: file_loc + integer :: i, ios, file_unit + integer, parameter :: NFIELDS_PER_IB = 20 + real(wp) :: ib_buf(NFIELDS_PER_IB) + + call s_create_directory(trim(case_dir) // '/restart_data') + + file_loc = trim(case_dir) // '/restart_data/ib_state_0.dat' + + open (newunit=file_unit, file=trim(file_loc), form='unformatted', access='stream', status='replace', iostat=ios) + if (ios /= 0) call s_mpi_abort('Cannot open IB state output file: ' // trim(file_loc)) + + ib_buf = 0._wp + + do i = 1, num_ibs + ib_buf(17) = patch_ib(i)%x_centroid + ib_buf(18) = patch_ib(i)%y_centroid + ib_buf(19) = patch_ib(i)%z_centroid + ib_buf(20) = patch_ib(i)%radius + write (file_unit) ib_buf + end do + + do i = 1, num_particle_cloud_ibs + ib_buf(17) = particle_cloud_ibs(i)%x_centroid + ib_buf(18) = particle_cloud_ibs(i)%y_centroid + ib_buf(19) = particle_cloud_ibs(i)%z_centroid + ib_buf(20) = particle_cloud_ibs(i)%radius + write (file_unit) ib_buf + end do + + close (file_unit) + + end subroutine s_write_ib_state_0 + !> Resets s_write_data_files pointer impure subroutine s_finalize_data_output_module diff --git a/src/pre_process/m_global_parameters.fpp b/src/pre_process/m_global_parameters.fpp index 7d85c2ac95..a3d7b835ee 100644 --- a/src/pre_process/m_global_parameters.fpp +++ b/src/pre_process/m_global_parameters.fpp @@ -361,6 +361,27 @@ contains patch_ib(i)%rotation_matrix_inverse = patch_ib(i)%rotation_matrix end do + num_particle_clouds = 0 + do i = 1, num_particle_clouds_max + particle_cloud(i)%x_centroid = 0._wp + particle_cloud(i)%y_centroid = 0._wp + particle_cloud(i)%z_centroid = 0._wp + particle_cloud(i)%length_x = dflt_real + particle_cloud(i)%length_y = dflt_real + particle_cloud(i)%length_z = dflt_real + particle_cloud(i)%num_particles = 0 + particle_cloud(i)%radius = dflt_real + particle_cloud(i)%mass = dflt_real + particle_cloud(i)%min_spacing = 0._wp + particle_cloud(i)%shell_inner_radius = dflt_real + particle_cloud(i)%shell_outer_radius = dflt_real + particle_cloud(i)%moving_ibm = 0 + particle_cloud(i)%seed = 0 + particle_cloud(i)%cloud_geometry = 1 + particle_cloud(i)%packing_method = dflt_int + particle_cloud(i)%periodic = 0 + end do + do i = 1, num_ib_airfoils_max ib_airfoil(i)%c = dflt_real ib_airfoil(i)%p = dflt_real diff --git a/src/simulation/m_particle_cloud.fpp b/src/pre_process/m_particle_cloud.fpp similarity index 83% rename from src/simulation/m_particle_cloud.fpp rename to src/pre_process/m_particle_cloud.fpp index 6e08433523..be26fde59b 100644 --- a/src/simulation/m_particle_cloud.fpp +++ b/src/pre_process/m_particle_cloud.fpp @@ -1,18 +1,19 @@ !> !! @file m_particle_cloud.fpp !! @brief Generates particle beds: converts particle_cloud specifications into -!! individual sphere/circle particle_cloud_ibs entries before reduction. +!! individual sphere/circle particle_cloud_ibs entries, written to +!! restart_data/ib_state_0.dat for simulation to read at startup. #:include 'macros.fpp' -!> @brief Generates particle beds by converting particle_cloud patch specifications into individual immersed boundary patches before -!! domain reduction. Each rank runs the same deterministic placement so no MPI broadcast of particle positions is needed. +!> @brief Generates particle beds by converting particle_cloud patch specifications into individual immersed boundary patches. Runs +!! once, on rank 0 only, producing the full global (unfiltered) particle list - simulation applies its own neighborhood filtering +!! when it reads the result back from restart_data/ib_state_0.dat. module m_particle_cloud use m_global_parameters use m_constants use m_mpi_common - use m_collisions implicit none @@ -22,11 +23,10 @@ module m_particle_cloud contains - !> Generate all particle beds and fill particle_cloud_ibs. Called on all ranks before s_reduce_ib_patch_array. Each packing - !! method owns and allocates its own per-cloud working array (see s_particle_cloud_lattice / s_particle_cloud_rejection_pack) - !! and hands back only the entries that fall within this rank's IB neighborhood. Only the first num_particle_cloud_ibs of them - !! are actually written - callers must use that count, not size(particle_cloud_ibs), since the remainder of the array is left - !! uninitialized. + !> Generate all particle beds and fill particle_cloud_ibs with the full, global (unfiltered) particle list. Called once, on rank + !! 0 only. Each packing method owns and allocates its own per-cloud working array (see s_particle_cloud_lattice / + !! s_particle_cloud_rejection_pack). Only the first num_particle_cloud_ibs of them are actually written - callers must use that + !! count, not size(particle_cloud_ibs), since the remainder of the array is left uninitialized. impure subroutine s_generate_particle_clouds(particle_cloud_ibs, num_particle_cloud_ibs) type(ib_patch_parameters), allocatable, intent(out), dimension(:) :: particle_cloud_ibs @@ -78,12 +78,10 @@ contains end subroutine s_generate_particle_clouds - !> Rejection-samples particle centres into a box or hemisphere-shell region with a minimum centre-to-centre spacing. Rejection - !! sampling needs every placed particle tracked (regardless of which rank's neighborhood it falls in) to detect overlaps - !! deterministically, so cloud_ibs is allocated here to the cloud's full requested particle count and only pared down to this - !! rank's neighborhood afterwards, via s_reduce_particle_cloud_ibs. Only the per-candidate geometry sampling differs between box + !> Rejection-samples particle centres into a box or hemisphere-shell region with a minimum centre-to-centre spacing. cloud_ibs + !! is allocated here to the cloud's full requested particle count. Only the per-candidate geometry sampling differs between box !! and hemisphere shell; it is delegated to s_sample_cloud_candidate, and every other step (overlap rejection via the spatial - !! hash, acceptance, reduction) is geometry-independent. + !! hash, acceptance) is geometry-independent. subroutine s_particle_cloud_rejection_pack(cloud_idx, glbl_idx, cloud_ibs, num_cloud_ibs) integer, intent(in) :: cloud_idx @@ -174,7 +172,6 @@ contains deallocate (placed, hash_head, chain_next) - call s_reduce_particle_cloud_ibs(cloud_ibs, ib_idx) num_cloud_ibs = ib_idx end subroutine s_particle_cloud_rejection_pack @@ -243,9 +240,7 @@ contains !> Places particles on the optimally dense lattice for the cloud region: a triangular lattice in 2D, a face-centered cubic !! lattice in 3D. The lattice spacing is set by the particle density (num_particles over the region area/volume); if that !! spacing falls below the required centre-to-centre distance (2*radius + min_spacing), the region is too dense and the run is - !! aborted. No two lattice sites can overlap, so unlike rejection packing each site's IB neighborhood membership - !! (get_neighbor_bounds() must already have run) is checked as it is generated and only in-neighborhood sites are stored; - !! cloud_ibs is therefore allocated to the neighborhood-sized cap rather than the cloud's full particle count. + !! aborted. No two lattice sites can overlap, so cloud_ibs is allocated to the cloud's full particle count. subroutine s_particle_cloud_lattice(cloud_idx, glbl_idx, cloud_ibs, num_cloud_ibs) integer, intent(in) :: cloud_idx @@ -259,7 +254,7 @@ contains real(wp), dimension(4) :: bx_off, by_off, bz_off real(wp), dimension(3) :: centroid - allocate (cloud_ibs(min(num_ib_patches_max_namelist, particle_cloud(cloud_idx)%num_particles))) + allocate (cloud_ibs(particle_cloud(cloud_idx)%num_particles)) ib_idx = 0 xmin = particle_cloud(cloud_idx)%x_centroid - 0.5_wp*particle_cloud(cloud_idx)%length_x @@ -301,10 +296,7 @@ contains do while (px <= xmax .and. n_placed < n_target) glbl_idx = glbl_idx + 1 centroid = [px, py, particle_cloud(cloud_idx)%z_centroid] - if (f_neighborhood_ranks_own_location(centroid)) then - call s_add_cloud_particle(cloud_idx, ib_idx, glbl_idx, geom, centroid(1), centroid(2), centroid(3), & - & cloud_ibs) - end if + call s_add_cloud_particle(cloud_idx, ib_idx, glbl_idx, geom, centroid(1), centroid(2), centroid(3), cloud_ibs) n_placed = n_placed + 1 col = col + 1 px = x0 + real(col, wp)*spacing @@ -328,10 +320,8 @@ contains centroid = [xmin + real(ix, wp)*cell + bx_off(b), ymin + real(jy, wp)*cell + by_off(b), & & zmin + real(kz, wp)*cell + bz_off(b)] glbl_idx = glbl_idx + 1 - if (f_neighborhood_ranks_own_location(centroid)) then - call s_add_cloud_particle(cloud_idx, ib_idx, glbl_idx, geom, centroid(1), centroid(2), & - & centroid(3), cloud_ibs) - end if + call s_add_cloud_particle(cloud_idx, ib_idx, glbl_idx, geom, centroid(1), centroid(2), centroid(3), & + & cloud_ibs) n_placed = n_placed + 1 end do end do @@ -344,11 +334,11 @@ contains end subroutine s_particle_cloud_lattice - !> Writes a single placed particle into particle_cloud_ibs at the next free slot, advancing ib_idx. The caller decides whether - !! this particle belongs in the array (neighborhood membership, for lattice packing, or unconditionally for rejection packing - - !! see s_particle_cloud_lattice / s_particle_cloud_rejection_pack) and supplies its already-assigned, absolute global patch id - !! via glbl_idx - s_reduce_ib_patch_array copies gbl_patch_id as-is. Shared by all packing methods so the per-particle - !! ib_patch_parameters setup stays in one place. + !> Writes a single placed particle into particle_cloud_ibs at the next free slot, advancing ib_idx, tagged with its + !! already-assigned, absolute global patch id via glbl_idx. Shared by all packing methods so the per-particle + !! ib_patch_parameters setup stays in one place. Only x/y/z_centroid and radius are ever read back out of this struct (by + !! s_write_ib_state_0); the rest is filled in for completeness and because leaving fields uninitialized here would otherwise + !! surface as garbage if this struct is ever read further downstream. subroutine s_add_cloud_particle(cloud_idx, ib_idx, glbl_idx, geom, px, py, pz, particle_cloud_ibs) integer, intent(in) :: cloud_idx, glbl_idx, geom @@ -358,7 +348,7 @@ contains ib_idx = ib_idx + 1 @:PROHIBIT(ib_idx > size(particle_cloud_ibs), & - & "Too many particle-cloud IBs in one rank's neighborhood. Modify case file or increase num_ib_patches_max_namelist.") + & "Too many particle-cloud IBs. Modify case file or increase num_ib_patches_max_namelist.") particle_cloud_ibs(ib_idx)%gbl_patch_id = glbl_idx particle_cloud_ibs(ib_idx)%geometry = geom @@ -388,12 +378,9 @@ contains particle_cloud_ibs(ib_idx)%moving_ibm = particle_cloud(cloud_idx)%moving_ibm particle_cloud_ibs(ib_idx)%slip = .false. - ! Particles are inert surfaces. These must be set explicitly: particle_cloud_ibs is - ! allocated (not default-initialized) and s_reduce_ib_patch_array copies the whole - ! struct into patch_ib, overwriting the defaults from - ! s_assign_default_values_to_user_inputs -- so anything left unset here reaches the - ! solver as uninitialized memory (a nonzero v_blow injects a garbage wall-normal - ! velocity and NaNs the field). + ! Particles are inert surfaces. particle_cloud_ibs is allocated (not default-initialized), + ! so these must be set explicitly even though only x/y/z_centroid and radius are read back + ! out of this struct downstream. particle_cloud_ibs(ib_idx)%v_blow = 0._wp particle_cloud_ibs(ib_idx)%inj_species = 0 particle_cloud_ibs(ib_idx)%burn_rate_exp = 0._wp @@ -500,29 +487,6 @@ contains end subroutine s_check_cloud_particle_overlap - !> Compacts cloud_ibs(1:num_ibs) in place, discarding entries outside this rank's IB neighborhood (get_neighbor_bounds() must - !! already have run) and updating num_ibs to the retained count. Used by rejection packing, which cannot filter as it places - !! particles (see s_particle_cloud_rejection_pack), to pare its full, unfiltered placement down to this rank's neighborhood. - subroutine s_reduce_particle_cloud_ibs(cloud_ibs, num_cloud_ibs) - - type(ib_patch_parameters), intent(inout), dimension(:) :: cloud_ibs - integer, intent(inout) :: num_cloud_ibs - integer :: i, write_idx - real(wp), dimension(3) :: centroid - - write_idx = 0 - do i = 1, num_cloud_ibs - centroid = [cloud_ibs(i)%x_centroid, cloud_ibs(i)%y_centroid, 0._wp] - if (num_dims == 3) centroid(3) = cloud_ibs(i)%z_centroid - if (f_neighborhood_ranks_own_location(centroid)) then - write_idx = write_idx + 1 - if (write_idx /= i) cloud_ibs(write_idx) = cloud_ibs(i) - end if - end do - num_cloud_ibs = write_idx - - end subroutine s_reduce_particle_cloud_ibs - !> Xorshift PRNG. Advances seed in-place and returns a value in [0, 1). function f_xorshift(seed) result(rval) diff --git a/src/pre_process/m_start_up.fpp b/src/pre_process/m_start_up.fpp index 04e3997378..bbf0a2ddbf 100644 --- a/src/pre_process/m_start_up.fpp +++ b/src/pre_process/m_start_up.fpp @@ -28,6 +28,7 @@ module m_start_up use m_check_patches use m_check_ib_patches + use m_particle_cloud use m_helper use m_checker_common use m_checker @@ -133,7 +134,21 @@ contains call s_check_patches() - if (ib) call s_check_ib_patches() + if (ib) then + call s_check_ib_patches() + + if (proc_rank == 0) then + block + type(ib_patch_parameters), allocatable :: particle_cloud_ibs(:) + integer :: num_particle_cloud_ibs + + call s_generate_particle_clouds(particle_cloud_ibs, num_particle_cloud_ibs) + call s_write_ib_state_0(particle_cloud_ibs, num_particle_cloud_ibs) + deallocate (particle_cloud_ibs) + end block + end if + call s_mpi_barrier() + end if end subroutine s_check_input_file diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index 485643aff5..26199cd414 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -43,7 +43,6 @@ module m_start_up use m_ibm use m_ib_patches use m_model - use m_particle_cloud use m_collisions use m_compile_specific use m_checker_common @@ -909,16 +908,12 @@ contains call s_initialize_ib_airfoils() call s_get_neighbor_bounds() + call s_assign_particle_cloud_ib_defaults(particle_cloud_ibs, num_particle_cloud_ibs) + call s_read_ib_restart_data(0, particle_cloud_ibs, num_particle_cloud_ibs) if (cfl_dt .and. n_start > 0) then - call s_read_ib_restart_data(n_start) - allocate (particle_cloud_ibs(0)) - num_particle_cloud_ibs = 0 + call s_read_ib_restart_data(n_start, particle_cloud_ibs, num_particle_cloud_ibs) else if (t_step_start > 0) then - call s_read_ib_restart_data(t_step_start) - allocate (particle_cloud_ibs(0)) - num_particle_cloud_ibs = 0 - else - call s_generate_particle_clouds(particle_cloud_ibs, num_particle_cloud_ibs) + call s_read_ib_restart_data(t_step_start, particle_cloud_ibs, num_particle_cloud_ibs) end if call s_reduce_ib_patch_array(particle_cloud_ibs, num_particle_cloud_ibs) deallocate (particle_cloud_ibs) @@ -1148,20 +1143,90 @@ contains end subroutine s_finalize_modules - !> @brief Reads IB kinematic state from restart_data/ib_state.dat on restart. Rank 0 reads the last num_ibs records and - !! broadcasts to all ranks. Overwrites patch_ib vel, angular_vel, angles, and centroid. - impure subroutine s_read_ib_restart_data(t_step) + !> @brief Allocates particle_cloud_ibs to sum(particle_cloud(:)%num_particles) and fills the properties that are constant per + !! cloud (geometry, mass, moving_ibm, slip, v_blow, inj_species, burn rate, identity rotation matrix, zeroed kinematics) - + !! mirroring the defaults s_add_cloud_particle (src/pre_process/m_particle_cloud.fpp) sets at generation time, in the same + !! cloud-then-particle-index order, so entry ib_idx here corresponds to global id num_ibs + ib_idx exactly as pre_process + !! assigned it. Position and radius are left at zero; s_read_ib_restart_data fills them in from restart_data/ib_state_0.dat. + impure subroutine s_assign_particle_cloud_ib_defaults(particle_cloud_ibs, num_particle_cloud_ibs) - integer, intent(in) :: t_step - character(len=path_len + 2*name_len) :: file_loc - integer :: i, ios, file_unit, ierr - integer :: r, nlocal, gbl_id - integer, parameter :: NFIELDS_PER_IB = 20 - real(wp) :: ib_buf(NFIELDS_PER_IB) - logical :: file_exist - character(len=10) :: t_step_string + type(ib_patch_parameters), allocatable, intent(out), dimension(:) :: particle_cloud_ibs + integer, intent(out) :: num_particle_cloud_ibs + integer :: cloud_idx, i, ib_idx, geom - if (file_per_process) then + num_particle_cloud_ibs = 0 + do cloud_idx = 1, num_particle_clouds + num_particle_cloud_ibs = num_particle_cloud_ibs + particle_cloud(cloud_idx)%num_particles + end do + allocate (particle_cloud_ibs(num_particle_cloud_ibs)) + + geom = merge(2, 8, num_dims < 3) + + ib_idx = 0 + do cloud_idx = 1, num_particle_clouds + do i = 1, particle_cloud(cloud_idx)%num_particles + ib_idx = ib_idx + 1 + + particle_cloud_ibs(ib_idx)%gbl_patch_id = num_ibs + ib_idx + particle_cloud_ibs(ib_idx)%geometry = geom + particle_cloud_ibs(ib_idx)%x_centroid = 0._wp + particle_cloud_ibs(ib_idx)%y_centroid = 0._wp + particle_cloud_ibs(ib_idx)%z_centroid = 0._wp + particle_cloud_ibs(ib_idx)%step_x_centroid = 0._wp + particle_cloud_ibs(ib_idx)%step_y_centroid = 0._wp + particle_cloud_ibs(ib_idx)%step_z_centroid = 0._wp + particle_cloud_ibs(ib_idx)%angles(:) = 0._wp + particle_cloud_ibs(ib_idx)%step_angles(:) = 0._wp + particle_cloud_ibs(ib_idx)%vel(:) = 0._wp + particle_cloud_ibs(ib_idx)%step_vel(:) = 0._wp + particle_cloud_ibs(ib_idx)%angular_vel(:) = 0._wp + particle_cloud_ibs(ib_idx)%step_angular_vel(:) = 0._wp + particle_cloud_ibs(ib_idx)%force(:) = 0._wp + particle_cloud_ibs(ib_idx)%torque(:) = 0._wp + particle_cloud_ibs(ib_idx)%centroid_offset(:) = 0._wp + particle_cloud_ibs(ib_idx)%rotation_matrix = 0._wp + particle_cloud_ibs(ib_idx)%rotation_matrix(1, 1) = 1._wp + particle_cloud_ibs(ib_idx)%rotation_matrix(2, 2) = 1._wp + particle_cloud_ibs(ib_idx)%rotation_matrix(3, 3) = 1._wp + particle_cloud_ibs(ib_idx)%rotation_matrix_inverse = particle_cloud_ibs(ib_idx)%rotation_matrix + particle_cloud_ibs(ib_idx)%radius = 0._wp + particle_cloud_ibs(ib_idx)%mass = particle_cloud(cloud_idx)%mass + particle_cloud_ibs(ib_idx)%moment = dflt_real + particle_cloud_ibs(ib_idx)%moving_ibm = particle_cloud(cloud_idx)%moving_ibm + particle_cloud_ibs(ib_idx)%slip = .false. + particle_cloud_ibs(ib_idx)%v_blow = 0._wp + particle_cloud_ibs(ib_idx)%inj_species = 0 + particle_cloud_ibs(ib_idx)%burn_rate_exp = 0._wp + particle_cloud_ibs(ib_idx)%burn_rate_pref = 0._wp + end do + end do + + end subroutine s_assign_particle_cloud_ib_defaults + + !> @brief Reads IB kinematic state (plus, for particle-cloud beds, radius) from restart_data/ib_state_.dat. Records with + !! gbl_id <= num_ibs (namelist patches) are applied directly to patch_ib, exactly as before. Records with gbl_id > num_ibs + !! (particle-cloud beds) are applied to particle_cloud_ibs(gbl_id - num_ibs) instead: patch_ib is GPU-declared at a fixed size + !! and, elsewhere in this module, only ever written in its final, already-neighborhood-filtered form (by + !! s_reduce_ib_patch_array) - bed entries stay in this plain staging array until that merge. Called with t_step = 0 to load the + !! initial layout pre_process generated (src/pre_process/m_data_output.fpp:s_write_ib_state_0) - since pre_process has no + !! file_per_process concept, t_step == 0 always uses the single-global-file layout below, regardless of this run's own + !! file_per_process setting - and again with the true step on an actual restart, overlaying live kinematics on top. + impure subroutine s_read_ib_restart_data(t_step, particle_cloud_ibs, num_particle_cloud_ibs) + + integer, intent(in) :: t_step + type(ib_patch_parameters), dimension(:), intent(inout) :: particle_cloud_ibs + integer, intent(in) :: num_particle_cloud_ibs + character(len=path_len + 2*name_len) :: file_loc + integer :: i, ios, file_unit, ierr + integer :: r, nlocal, gbl_id, n_records + integer, parameter :: NFIELDS_PER_IB = 20 + real(wp) :: ib_buf(NFIELDS_PER_IB) + logical :: file_exist + character(len=10) :: t_step_string + + n_records = num_ibs + num_particle_cloud_ibs + + if (file_per_process .and. t_step /= 0) then call s_int_to_str(t_step, t_step_string) do r = 0, num_procs - 1 @@ -1183,12 +1248,11 @@ contains read (file_unit, iostat=ios) ib_buf if (ios /= 0) call s_mpi_abort('Error reading IB state data: ' // trim(file_loc)) - patch_ib(gbl_id)%vel = ib_buf(8:10) - patch_ib(gbl_id)%angular_vel = ib_buf(11:13) - patch_ib(gbl_id)%angles = ib_buf(14:16) - patch_ib(gbl_id)%x_centroid = ib_buf(17) - patch_ib(gbl_id)%y_centroid = ib_buf(18) - patch_ib(gbl_id)%z_centroid = ib_buf(19) + if (gbl_id <= num_ibs) then + call s_apply_ib_buf(patch_ib(gbl_id), ib_buf) + else + call s_apply_ib_buf(particle_cloud_ibs(gbl_id - num_ibs), ib_buf) + end if end do close (file_unit) @@ -1206,16 +1270,15 @@ contains open (newunit=file_unit, file=trim(file_loc), form='unformatted', access='stream', status='old', iostat=ios) if (ios /= 0) call s_mpi_abort('Error opening IB state restart file: ' // trim(file_loc)) - do i = 1, num_ibs + do i = 1, n_records read (file_unit, iostat=ios) ib_buf if (ios /= 0) call s_mpi_abort('Error reading IB state restart file') - patch_ib(i)%vel = ib_buf(8:10) - patch_ib(i)%angular_vel = ib_buf(11:13) - patch_ib(i)%angles = ib_buf(14:16) - patch_ib(i)%x_centroid = ib_buf(17) - patch_ib(i)%y_centroid = ib_buf(18) - patch_ib(i)%z_centroid = ib_buf(19) + if (i <= num_ibs) then + call s_apply_ib_buf(patch_ib(i), ib_buf) + else + call s_apply_ib_buf(particle_cloud_ibs(i - num_ibs), ib_buf) + end if end do close (file_unit) @@ -1229,18 +1292,45 @@ contains call MPI_BCAST(patch_ib(i)%x_centroid, 1, mpi_p, 0, MPI_COMM_WORLD, ierr) call MPI_BCAST(patch_ib(i)%y_centroid, 1, mpi_p, 0, MPI_COMM_WORLD, ierr) call MPI_BCAST(patch_ib(i)%z_centroid, 1, mpi_p, 0, MPI_COMM_WORLD, ierr) + call MPI_BCAST(patch_ib(i)%radius, 1, mpi_p, 0, MPI_COMM_WORLD, ierr) + end do + do i = 1, num_particle_cloud_ibs + call MPI_BCAST(particle_cloud_ibs(i)%vel, 3, mpi_p, 0, MPI_COMM_WORLD, ierr) + call MPI_BCAST(particle_cloud_ibs(i)%angular_vel, 3, mpi_p, 0, MPI_COMM_WORLD, ierr) + call MPI_BCAST(particle_cloud_ibs(i)%angles, 3, mpi_p, 0, MPI_COMM_WORLD, ierr) + call MPI_BCAST(particle_cloud_ibs(i)%x_centroid, 1, mpi_p, 0, MPI_COMM_WORLD, ierr) + call MPI_BCAST(particle_cloud_ibs(i)%y_centroid, 1, mpi_p, 0, MPI_COMM_WORLD, ierr) + call MPI_BCAST(particle_cloud_ibs(i)%z_centroid, 1, mpi_p, 0, MPI_COMM_WORLD, ierr) + call MPI_BCAST(particle_cloud_ibs(i)%radius, 1, mpi_p, 0, MPI_COMM_WORLD, ierr) end do #endif end if + contains + + !> Applies one 20-field IB state record to a single patch_ib/particle_cloud_ibs entry. + subroutine s_apply_ib_buf(entry, buf) + + type(ib_patch_parameters), intent(inout) :: entry + real(wp), dimension(NFIELDS_PER_IB), intent(in) :: buf + + entry%vel = buf(8:10) + entry%angular_vel = buf(11:13) + entry%angles = buf(14:16) + entry%x_centroid = buf(17) + entry%y_centroid = buf(18) + entry%z_centroid = buf(19) + entry%radius = buf(20) + + end subroutine s_apply_ib_buf + end subroutine s_read_ib_restart_data - !> @brief Merges patch_ib (namelist patches, fixed at num_ib_patches_max_namelist) with particle_cloud_ibs (already filtered by - !! s_generate_particle_clouds to this rank's IB neighborhood, each entry already tagged with its final, absolute gbl_patch_id) - !! and reduces to only the patches in or near the local computational domain. patch_ib is never reallocated; the local subset is - !! written in-place from the front. particle_cloud_ibs is owned by the caller and freed there after this returns. - !! num_particle_cloud_ibs is the number of entries s_generate_particle_clouds actually wrote into particle_cloud_ibs - it may be - !! allocated to a larger worst-case capacity, so size() of it must never be used as the valid-entry count. + !> @brief Merges patch_ib (namelist patches, fixed at num_ib_patches_max_namelist) with particle_cloud_ibs (the full, global, + !! unfiltered bed particle set - see s_assign_particle_cloud_ib_defaults/s_read_ib_restart_data - each entry already tagged with + !! its final, absolute gbl_patch_id) and reduces to only the patches in or near the local computational domain. patch_ib is + !! never reallocated; the local subset is written in-place from the front. particle_cloud_ibs is owned by the caller and freed + !! there after this returns. subroutine s_reduce_ib_patch_array(particle_cloud_ibs, num_particle_cloud_ibs) type(ib_patch_parameters), intent(in), dimension(:) :: particle_cloud_ibs @@ -1308,19 +1398,22 @@ contains end if end if end do - ! particle_cloud_ibs entries already passed the neighborhood check at generation time so no need to recheck it here. + ! particle_cloud_ibs holds the full global bed set; filter to this rank's neighborhood exactly as above for namelist + ! patches. do i = 1, num_particle_cloud_ibs centroid = [particle_cloud_ibs(i)%x_centroid, particle_cloud_ibs(i)%y_centroid, 0._wp] if (num_dims == 3) centroid(3) = particle_cloud_ibs(i)%z_centroid - num_ibs = num_ibs + 1 - @:PROHIBIT(num_ibs > num_ib_patches_max_namelist, & - & "Local IB count exceeds patch_ib capacity. Increase num_ib_patches_max_namelist.") - patch_ib(num_ibs) = particle_cloud_ibs(i) - if (f_local_rank_owns_location(centroid)) then - num_local_ibs = num_local_ibs + 1 - @:PROHIBIT(num_local_ibs > num_local_ibs_max, & - & "Too many IBs on a single processor rank. Modify case file or increase limit of num_local_ibs_max to resolve.") - local_ib_patch_ids(num_local_ibs) = num_ibs + if (f_neighborhood_ranks_own_location(centroid)) then + num_ibs = num_ibs + 1 + @:PROHIBIT(num_ibs > num_ib_patches_max_namelist, & + & "Local IB count exceeds patch_ib capacity. Increase num_ib_patches_max_namelist.") + patch_ib(num_ibs) = particle_cloud_ibs(i) + if (f_local_rank_owns_location(centroid)) then + num_local_ibs = num_local_ibs + 1 + @:PROHIBIT(num_local_ibs > num_local_ibs_max, & + & "Too many IBs on a single processor rank. Modify case file or increase limit of num_local_ibs_max to resolve.") + local_ib_patch_ids(num_local_ibs) = num_ibs + end if end if end do end if diff --git a/toolchain/mfc/params/definitions.py b/toolchain/mfc/params/definitions.py index ca888ab52f..d91656456f 100644 --- a/toolchain/mfc/params/definitions.py +++ b/toolchain/mfc/params/definitions.py @@ -1346,6 +1346,9 @@ def _decl(targets: set, *names: str) -> None: "avg_state", "alt_soundspeed", "mixture_err", +) +_nv( + _ALL, "num_particle_clouds", "particle_cloud", ) diff --git a/toolchain/mfc/test/cases.py b/toolchain/mfc/test/cases.py index ff81664ff4..fbce5f63a5 100644 --- a/toolchain/mfc/test/cases.py +++ b/toolchain/mfc/test/cases.py @@ -1309,6 +1309,38 @@ def alter_ib(dimInfo, six_eqn_model=False, viscous=False): ) ) + # Restart roundtrip regression: particle-cloud beds must survive a restart, not just namelist patch_ib patches - + # pre_process now generates them once and simulation reads that layout back on every start (fresh or restart). + cases.append( + define_case_d( + stack, + "IBM -> Particle Cloud -> Box -> Restart", + { + "ib": "T", + "num_ibs": 0, + "num_particle_clouds": 1, + "fd_order": 2, + "n": 49, + "particle_cloud(1)%cloud_geometry": 1, + "particle_cloud(1)%packing_method": 1, + "particle_cloud(1)%x_centroid": 0.5, + "particle_cloud(1)%y_centroid": 0.5, + "particle_cloud(1)%length_x": 0.6, + "particle_cloud(1)%length_y": 0.6, + "particle_cloud(1)%num_particles": 4, + "particle_cloud(1)%radius": 0.02, + "particle_cloud(1)%mass": 1.0, + "particle_cloud(1)%min_spacing": 0.005, + "particle_cloud(1)%moving_ibm": 0, + "particle_cloud(1)%seed": 12345, + "patch_icpp(1)%vel(1)": 0.001, + "patch_icpp(2)%vel(1)": 0.001, + "patch_icpp(3)%vel(1)": 0.001, + }, + restart_check=True, + ) + ) + if len(dimInfo[0]) == 3 and not viscous: cases.append( define_case_d( From d7e05742d251df95a32d6caf59b0bca34cd822d8 Mon Sep 17 00:00:00 2001 From: "Daniel J. Vickers" Date: Thu, 10 Sep 2026 11:44:37 -0400 Subject: [PATCH 2/5] Intermitent commit, --- src/common/m_helper.fpp | 45 ++++++- src/pre_process/m_data_output.fpp | 181 ++++++++++++++++++++++++--- src/pre_process/m_particle_cloud.fpp | 85 +++++++++---- src/pre_process/m_start_up.fpp | 41 +++--- src/pre_process/p_main.f90 | 2 + src/simulation/m_collisions.fpp | 41 +----- src/simulation/m_ibm.fpp | 2 +- src/simulation/m_start_up.fpp | 12 +- 8 files changed, 310 insertions(+), 99 deletions(-) diff --git a/src/common/m_helper.fpp b/src/common/m_helper.fpp index 4a456d6615..5a046c769f 100644 --- a/src/common/m_helper.fpp +++ b/src/common/m_helper.fpp @@ -10,6 +10,7 @@ module m_helper use m_derived_types use m_global_parameters + use m_constants use ieee_arithmetic !< For checking NaN implicit none @@ -18,7 +19,8 @@ module m_helper public :: s_comp_n_from_prim, s_comp_n_from_cons, s_initialize_bubbles_model, s_initialize_nonpoly, s_simpson, s_transcoeff, & & s_int_to_str, s_transform_vec, s_transform_triangle, s_transform_model, s_swap, f_cross, f_create_transform_matrix, & & f_create_bbox, s_print_2D_array, f_xor, f_logical_to_int, associated_legendre, real_ylm, double_factorial, factorial, & - & f_cut_on, f_cut_off, s_downsample_data, s_upsample_data, s_cross_product, f_unit_vector, s_prng, modmul + & f_cut_on, f_cut_off, s_downsample_data, s_upsample_data, s_cross_product, f_unit_vector, s_prng, modmul, & + & f_local_rank_owns_location contains @@ -699,4 +701,45 @@ contains end subroutine s_upsample_data + !> @brief True if `location` falls within this rank's own subdomain (a strict partition - each location is owned by exactly one + !! rank, unlike an overlapping multi-rank ghost-stencil neighborhood). Used by both pre_process (to decide which + !! generated/namelist IBs a rank writes to its own restart_data/ib_state_0.dat chunk) and simulation (to decide which IBs a rank + !! owns for its own periodic IB-state writes) so the two stay consistent. glb_bounds_in is the global domain extent, used only + !! to project a location that falls just outside the domain (floating-point edge case) onto the domain so some rank still claims + !! it. + function f_local_rank_owns_location(location, glb_bounds_in) result(owns_location) + + $:GPU_ROUTINE(parallelism='[seq]') + + real(wp), dimension(3), intent(in) :: location + type(bounds_info), dimension(3), intent(in) :: glb_bounds_in + logical :: owns_location + real(wp), dimension(3) :: projected_location + + owns_location = .true. + +#ifdef MFC_MPI + if (num_procs > 1) then + projected_location(:) = location(:) + + ! catch the edge case where the location lies just outside the computational domain + #:for X, ID, DIM in [('x', 1, 'm'), ('y', 2, 'n'), ('z', 3, 'p')] + if (num_dims >= ${ID}$) then + if (bc_${X}$%beg /= BC_PERIODIC) then + ! if it is outside the domain in one direction, project it somewhere inside so at least one rank owns it + if (location(${ID}$) < glb_bounds_in(${ID}$)%beg) then + projected_location(${ID}$) = glb_bounds_in(${ID}$)%beg + else if (glb_bounds_in(${ID}$)%end < location(${ID}$)) then + projected_location(${ID}$) = glb_bounds_in(${ID}$)%end - 1.0e-10_wp + end if + end if + owns_location = owns_location .and. ${X}$_cb(-1) <= projected_location(${ID}$) & + & .and. projected_location(${ID}$) < ${X}$_cb(${DIM}$) + end if + #:endfor + end if +#endif + + end function f_local_rank_owns_location + end module m_helper diff --git a/src/pre_process/m_data_output.fpp b/src/pre_process/m_data_output.fpp index f8f8c59e41..1519d286a3 100644 --- a/src/pre_process/m_data_output.fpp +++ b/src/pre_process/m_data_output.fpp @@ -28,7 +28,7 @@ module m_data_output private public :: s_write_serial_data_files, s_write_parallel_data_files, s_write_data_files, s_initialize_data_output_module, & - & s_finalize_data_output_module, s_write_ib_state_0 + & s_finalize_data_output_module, s_write_ib_state_0_file type(scalar_field), allocatable, dimension(:) :: q_cons_temp @@ -733,14 +733,59 @@ contains end subroutine s_initialize_data_output_module - !> @brief Writes restart_data/ib_state_0.dat: the initial IB layout (namelist patch_ib entries, then any generated - !! particle-cloud beds). Rank-0-only. Read back by simulation at startup via s_read_ib_restart_data(0, ...) - !! (src/simulation/m_start_up.fpp). Uses the same 20-field record layout s_write_serial_ib_state - !! (src/simulation/m_data_output.fpp) writes, so simulation's own writers can freely overwrite this file later at t_step_start - !! == 0 without changing its layout - only position (fields 17:19) and radius (field 20) are populated here; everything else - !! (time, force, torque, vel, angular_vel, angles) is zero for a freshly generated IB. - impure subroutine s_write_ib_state_0(particle_cloud_ibs, num_particle_cloud_ibs) + !> @brief Writes restart_data/ib_state_0.dat (or, under file_per_process, one restart_data/lustre_0/ib_state_0_.dat chunk + !! per rank): the initial IB layout - namelist patch_ib entries, then any generated particle-cloud beds - with each rank writing + !! only the entries f_local_rank_owns_location says are its own. Read back by simulation at startup via + !! s_read_ib_restart_data(0, ...) (src/simulation/m_start_up.fpp), which dispatches on the same file_per_process flag, so the + !! two writers must stay format-compatible: this mirrors s_write_parallel_ib_state/s_write_serial_ib_state + !! (src/simulation/m_data_output.fpp) exactly, just sourcing entries from this rank's local subset instead of patch_ib/ + !! local_ib_patch_ids post-reduce. Only position (fields 17:19) and radius (field 20) are populated here; everything else (time, + !! force, torque, vel, angular_vel, angles) is zero for a freshly generated IB. + impure subroutine s_write_ib_state_0_file(glb_bounds, particle_cloud_ibs, num_particle_cloud_ibs) + + type(bounds_info), dimension(3), intent(in) :: glb_bounds + type(ib_patch_parameters), dimension(:), intent(in) :: particle_cloud_ibs + integer, intent(in) :: num_particle_cloud_ibs + integer, allocatable :: local_namelist_ids(:) + integer :: num_local_namelist + integer :: i + real(wp), dimension(3) :: centroid + + allocate (local_namelist_ids(max(1, num_ibs))) + num_local_namelist = 0 + do i = 1, num_ibs + centroid = [patch_ib(i)%x_centroid, patch_ib(i)%y_centroid, 0._wp] + if (num_dims == 3) centroid(3) = patch_ib(i)%z_centroid + if (f_local_rank_owns_location(centroid, glb_bounds)) then + num_local_namelist = num_local_namelist + 1 + local_namelist_ids(num_local_namelist) = i + end if + end do + + if (.not. parallel_io) then + ! Mirrors s_write_serial_ib_state (src/simulation/m_data_output.fpp): "must be called only on rank 0" - parallel_io + ! = F is not MPI-aware for ib_state, so this combined with num_procs > 1 is a pre-existing limitation, not one + ! introduced here. + if (proc_rank == 0) call s_write_ib_state_0_serial(local_namelist_ids, num_local_namelist, particle_cloud_ibs, & + & num_particle_cloud_ibs) + else if (file_per_process) then + call s_write_ib_state_0_file_per_process(local_namelist_ids, num_local_namelist, particle_cloud_ibs, & + & num_particle_cloud_ibs) + else + call s_write_ib_state_0_shared(local_namelist_ids, num_local_namelist, particle_cloud_ibs, num_particle_cloud_ibs) + end if + deallocate (local_namelist_ids) + + end subroutine s_write_ib_state_0_file + + !> Writes this rank's local entries into its own restart_data/lustre_0/ib_state_0_.dat chunk file - mirrors + !! s_write_parallel_ib_state's file_per_process branch (src/simulation/m_data_output.fpp). + subroutine s_write_ib_state_0_file_per_process(local_namelist_ids, num_local_namelist, particle_cloud_ibs, & + & num_particle_cloud_ibs) + + integer, dimension(:), intent(in) :: local_namelist_ids + integer, intent(in) :: num_local_namelist type(ib_patch_parameters), dimension(:), intent(in) :: particle_cloud_ibs integer, intent(in) :: num_particle_cloud_ibs character(LEN=len_trim(case_dir) + 2*name_len) :: file_loc @@ -748,8 +793,114 @@ contains integer, parameter :: NFIELDS_PER_IB = 20 real(wp) :: ib_buf(NFIELDS_PER_IB) - call s_create_directory(trim(case_dir) // '/restart_data') + if (proc_rank == 0) call s_create_directory(trim(case_dir) // '/restart_data/lustre_0') + call s_mpi_barrier() + call s_delay_file_access(proc_rank) + + write (file_loc, '(A,i7.7,A)') 'ib_state_0_', proc_rank, '.dat' + file_loc = trim(case_dir) // '/restart_data/lustre_0/' // trim(file_loc) + + open (newunit=file_unit, file=trim(file_loc), form='unformatted', access='stream', status='replace', iostat=ios) + if (ios /= 0) call s_mpi_abort('Cannot open IB state output file: ' // trim(file_loc)) + + ib_buf = 0._wp + + write (file_unit) num_local_namelist + num_particle_cloud_ibs + + do i = 1, num_local_namelist + ib_buf(17) = patch_ib(local_namelist_ids(i))%x_centroid + ib_buf(18) = patch_ib(local_namelist_ids(i))%y_centroid + ib_buf(19) = patch_ib(local_namelist_ids(i))%z_centroid + ib_buf(20) = patch_ib(local_namelist_ids(i))%radius + write (file_unit) local_namelist_ids(i) + write (file_unit) ib_buf + end do + + do i = 1, num_particle_cloud_ibs + ib_buf(17) = particle_cloud_ibs(i)%x_centroid + ib_buf(18) = particle_cloud_ibs(i)%y_centroid + ib_buf(19) = particle_cloud_ibs(i)%z_centroid + ib_buf(20) = particle_cloud_ibs(i)%radius + write (file_unit) particle_cloud_ibs(i)%gbl_patch_id + write (file_unit) ib_buf + end do + + close (file_unit) + + end subroutine s_write_ib_state_0_file_per_process + !> Writes this rank's local entries into the single shared restart_data/ib_state_0.dat, each rank placing its own records at + !! their gbl_patch_id-based offset via MPI-IO - mirrors s_write_parallel_ib_state's non-file_per_process branch + !! (src/simulation/m_data_output.fpp). Only reached when parallel_io = T (see s_write_ib_state_0_file), matching that + !! routine's own precondition of running under MFC_MPI. + subroutine s_write_ib_state_0_shared(local_namelist_ids, num_local_namelist, particle_cloud_ibs, num_particle_cloud_ibs) + + integer, dimension(:), intent(in) :: local_namelist_ids + integer, intent(in) :: num_local_namelist + type(ib_patch_parameters), dimension(:), intent(in) :: particle_cloud_ibs + integer, intent(in) :: num_particle_cloud_ibs + integer, parameter :: NFIELDS_PER_IB = 20 + real(wp) :: ib_buf(NFIELDS_PER_IB) + integer :: i +#ifdef MFC_MPI + character(LEN=len_trim(case_dir) + 2*name_len) :: file_loc + integer(kind=MPI_OFFSET_KIND) :: disp, WP_MOK + integer :: ifile, ierr + integer, dimension(MPI_STATUS_SIZE) :: status + logical :: file_exist + + if (proc_rank == 0) call s_create_directory(trim(case_dir) // '/restart_data') + call s_mpi_barrier() + + WP_MOK = int(storage_size(0._wp)/8, MPI_OFFSET_KIND) + file_loc = trim(case_dir) // '/restart_data/ib_state_0.dat' + + inquire (FILE=trim(file_loc), EXIST=file_exist) + if (file_exist .and. proc_rank == 0) call MPI_FILE_DELETE(file_loc, mpi_info_int, ierr) + call s_mpi_barrier() + + call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, ior(MPI_MODE_WRONLY, MPI_MODE_CREATE), mpi_info_int, ifile, ierr) + + ib_buf = 0._wp + + do i = 1, num_local_namelist + ib_buf(17) = patch_ib(local_namelist_ids(i))%x_centroid + ib_buf(18) = patch_ib(local_namelist_ids(i))%y_centroid + ib_buf(19) = patch_ib(local_namelist_ids(i))%z_centroid + ib_buf(20) = patch_ib(local_namelist_ids(i))%radius + disp = int(local_namelist_ids(i) - 1, MPI_OFFSET_KIND)*int(NFIELDS_PER_IB, MPI_OFFSET_KIND)*WP_MOK + call MPI_FILE_WRITE_AT(ifile, disp, ib_buf, NFIELDS_PER_IB, mpi_p, status, ierr) + end do + + do i = 1, num_particle_cloud_ibs + ib_buf(17) = particle_cloud_ibs(i)%x_centroid + ib_buf(18) = particle_cloud_ibs(i)%y_centroid + ib_buf(19) = particle_cloud_ibs(i)%z_centroid + ib_buf(20) = particle_cloud_ibs(i)%radius + disp = int(particle_cloud_ibs(i)%gbl_patch_id - 1, MPI_OFFSET_KIND)*int(NFIELDS_PER_IB, MPI_OFFSET_KIND)*WP_MOK + call MPI_FILE_WRITE_AT(ifile, disp, ib_buf, NFIELDS_PER_IB, mpi_p, status, ierr) + end do + + call MPI_FILE_CLOSE(ifile, ierr) +#endif + + end subroutine s_write_ib_state_0_shared + + !> Writes ALL of this rank's local entries directly into restart_data/ib_state_0.dat with plain Fortran I/O - mirrors + !! s_write_serial_ib_state (src/simulation/m_data_output.fpp). Used only when parallel_io = F; the caller (rank 0 only) mirrors + !! that routine's "must be called only on rank 0" precondition, so this is a plain sequential write, no MPI needed. + subroutine s_write_ib_state_0_serial(local_namelist_ids, num_local_namelist, particle_cloud_ibs, num_particle_cloud_ibs) + + integer, dimension(:), intent(in) :: local_namelist_ids + integer, intent(in) :: num_local_namelist + type(ib_patch_parameters), dimension(:), intent(in) :: particle_cloud_ibs + integer, intent(in) :: num_particle_cloud_ibs + character(LEN=len_trim(case_dir) + 2*name_len) :: file_loc + integer :: i, ios, file_unit + integer, parameter :: NFIELDS_PER_IB = 20 + real(wp) :: ib_buf(NFIELDS_PER_IB) + + call s_create_directory(trim(case_dir) // '/restart_data') file_loc = trim(case_dir) // '/restart_data/ib_state_0.dat' open (newunit=file_unit, file=trim(file_loc), form='unformatted', access='stream', status='replace', iostat=ios) @@ -757,11 +908,11 @@ contains ib_buf = 0._wp - do i = 1, num_ibs - ib_buf(17) = patch_ib(i)%x_centroid - ib_buf(18) = patch_ib(i)%y_centroid - ib_buf(19) = patch_ib(i)%z_centroid - ib_buf(20) = patch_ib(i)%radius + do i = 1, num_local_namelist + ib_buf(17) = patch_ib(local_namelist_ids(i))%x_centroid + ib_buf(18) = patch_ib(local_namelist_ids(i))%y_centroid + ib_buf(19) = patch_ib(local_namelist_ids(i))%z_centroid + ib_buf(20) = patch_ib(local_namelist_ids(i))%radius write (file_unit) ib_buf end do @@ -775,7 +926,7 @@ contains close (file_unit) - end subroutine s_write_ib_state_0 + end subroutine s_write_ib_state_0_serial !> Resets s_write_data_files pointer impure subroutine s_finalize_data_output_module diff --git a/src/pre_process/m_particle_cloud.fpp b/src/pre_process/m_particle_cloud.fpp index be26fde59b..37bc7c1896 100644 --- a/src/pre_process/m_particle_cloud.fpp +++ b/src/pre_process/m_particle_cloud.fpp @@ -6,14 +6,17 @@ #:include 'macros.fpp' -!> @brief Generates particle beds by converting particle_cloud patch specifications into individual immersed boundary patches. Runs -!! once, on rank 0 only, producing the full global (unfiltered) particle list - simulation applies its own neighborhood filtering -!! when it reads the result back from restart_data/ib_state_0.dat. +!> @brief Generates particle beds by converting particle_cloud patch specifications into individual immersed boundary patches before +!! writing them to restart_data/ib_state_0.dat. Runs on every rank: each rank computes the same deterministic placement (so no MPI +!! broadcast of particle positions is needed) and keeps only the particles f_local_rank_owns_location says are its own, so every +!! generated particle is written by exactly one rank - consistent with how simulation partitions patch_ib for its own restart writes +!! (src/simulation/m_start_up.fpp:s_reduce_ib_patch_array). module m_particle_cloud use m_global_parameters use m_constants use m_mpi_common + use m_helper implicit none @@ -23,12 +26,14 @@ module m_particle_cloud contains - !> Generate all particle beds and fill particle_cloud_ibs with the full, global (unfiltered) particle list. Called once, on rank - !! 0 only. Each packing method owns and allocates its own per-cloud working array (see s_particle_cloud_lattice / - !! s_particle_cloud_rejection_pack). Only the first num_particle_cloud_ibs of them are actually written - callers must use that - !! count, not size(particle_cloud_ibs), since the remainder of the array is left uninitialized. - impure subroutine s_generate_particle_clouds(particle_cloud_ibs, num_particle_cloud_ibs) + !> Generate all particle beds and fill particle_cloud_ibs, keeping only the particles this rank owns (see module docs). Called + !! on every rank. Each packing method owns and allocates its own per-cloud working array (see s_particle_cloud_lattice / + !! s_particle_cloud_rejection_pack) and hands back only the entries that fall within this rank's own subdomain. Only the first + !! num_particle_cloud_ibs of them are actually written - callers must use that count, not size(particle_cloud_ibs), since the + !! remainder of the array is left uninitialized. + impure subroutine s_generate_particle_clouds(glb_bounds, particle_cloud_ibs, num_particle_cloud_ibs) + type(bounds_info), dimension(3), intent(in) :: glb_bounds type(ib_patch_parameters), allocatable, intent(out), dimension(:) :: particle_cloud_ibs integer, intent(out) :: num_particle_cloud_ibs type(ib_patch_parameters), allocatable :: cloud_ibs(:) @@ -58,9 +63,9 @@ contains ! packing is box-only - the hemisphere-shell + lattice combination is rejected in case_validator.py. select case (particle_cloud(cloud_idx)%packing_method) case (1) ! rejection (random) packing method - call s_particle_cloud_rejection_pack(cloud_idx, glbl_idx, cloud_ibs, num_cloud_ibs) + call s_particle_cloud_rejection_pack(cloud_idx, glbl_idx, glb_bounds, cloud_ibs, num_cloud_ibs) case (2) ! lattice packing method - call s_particle_cloud_lattice(cloud_idx, glbl_idx, cloud_ibs, num_cloud_ibs) + call s_particle_cloud_lattice(cloud_idx, glbl_idx, glb_bounds, cloud_ibs, num_cloud_ibs) case default call s_mpi_abort("Particle cloud packing method is not a known packing method of MFC. Exiting.") end select @@ -78,14 +83,17 @@ contains end subroutine s_generate_particle_clouds - !> Rejection-samples particle centres into a box or hemisphere-shell region with a minimum centre-to-centre spacing. cloud_ibs - !! is allocated here to the cloud's full requested particle count. Only the per-candidate geometry sampling differs between box - !! and hemisphere shell; it is delegated to s_sample_cloud_candidate, and every other step (overlap rejection via the spatial - !! hash, acceptance) is geometry-independent. - subroutine s_particle_cloud_rejection_pack(cloud_idx, glbl_idx, cloud_ibs, num_cloud_ibs) + !> Rejection-samples particle centres into a box or hemisphere-shell region with a minimum centre-to-centre spacing. Rejection + !! sampling needs every placed particle tracked (regardless of which rank owns it) to detect overlaps deterministically, so + !! cloud_ibs is allocated here to the cloud's full requested particle count and only pared down to this rank's own subdomain + !! afterwards, via s_reduce_particle_cloud_ibs. Only the per-candidate geometry sampling differs between box and hemisphere + !! shell; it is delegated to s_sample_cloud_candidate, and every other step (overlap rejection via the spatial hash, acceptance, + !! reduction) is geometry-independent. + subroutine s_particle_cloud_rejection_pack(cloud_idx, glbl_idx, glb_bounds, cloud_ibs, num_cloud_ibs) integer, intent(in) :: cloud_idx integer, intent(inout) :: glbl_idx + type(bounds_info), dimension(3), intent(in) :: glb_bounds type(ib_patch_parameters), allocatable, intent(out), dimension(:) :: cloud_ibs integer, intent(out) :: num_cloud_ibs integer :: ib_idx, n_placed, geom, seed, alloc_stat @@ -172,6 +180,7 @@ contains deallocate (placed, hash_head, chain_next) + call s_reduce_particle_cloud_ibs(cloud_ibs, glb_bounds, ib_idx) num_cloud_ibs = ib_idx end subroutine s_particle_cloud_rejection_pack @@ -240,11 +249,14 @@ contains !> Places particles on the optimally dense lattice for the cloud region: a triangular lattice in 2D, a face-centered cubic !! lattice in 3D. The lattice spacing is set by the particle density (num_particles over the region area/volume); if that !! spacing falls below the required centre-to-centre distance (2*radius + min_spacing), the region is too dense and the run is - !! aborted. No two lattice sites can overlap, so cloud_ibs is allocated to the cloud's full particle count. - subroutine s_particle_cloud_lattice(cloud_idx, glbl_idx, cloud_ibs, num_cloud_ibs) + !! aborted. No two lattice sites can overlap, so unlike rejection packing each site's local ownership + !! (f_local_rank_owns_location) is checked as it is generated and only this rank's sites are stored; cloud_ibs is therefore + !! allocated to a worst-case cap rather than the cloud's full particle count. + subroutine s_particle_cloud_lattice(cloud_idx, glbl_idx, glb_bounds, cloud_ibs, num_cloud_ibs) integer, intent(in) :: cloud_idx integer, intent(inout) :: glbl_idx + type(bounds_info), dimension(3), intent(in) :: glb_bounds type(ib_patch_parameters), allocatable, intent(out), dimension(:) :: cloud_ibs integer, intent(out) :: num_cloud_ibs integer :: ib_idx, n_placed, n_target, geom @@ -254,7 +266,7 @@ contains real(wp), dimension(4) :: bx_off, by_off, bz_off real(wp), dimension(3) :: centroid - allocate (cloud_ibs(particle_cloud(cloud_idx)%num_particles)) + allocate (cloud_ibs(min(num_ib_patches_max_namelist, particle_cloud(cloud_idx)%num_particles))) ib_idx = 0 xmin = particle_cloud(cloud_idx)%x_centroid - 0.5_wp*particle_cloud(cloud_idx)%length_x @@ -296,7 +308,10 @@ contains do while (px <= xmax .and. n_placed < n_target) glbl_idx = glbl_idx + 1 centroid = [px, py, particle_cloud(cloud_idx)%z_centroid] - call s_add_cloud_particle(cloud_idx, ib_idx, glbl_idx, geom, centroid(1), centroid(2), centroid(3), cloud_ibs) + if (f_local_rank_owns_location(centroid, glb_bounds)) then + call s_add_cloud_particle(cloud_idx, ib_idx, glbl_idx, geom, centroid(1), centroid(2), centroid(3), & + & cloud_ibs) + end if n_placed = n_placed + 1 col = col + 1 px = x0 + real(col, wp)*spacing @@ -320,8 +335,10 @@ contains centroid = [xmin + real(ix, wp)*cell + bx_off(b), ymin + real(jy, wp)*cell + by_off(b), & & zmin + real(kz, wp)*cell + bz_off(b)] glbl_idx = glbl_idx + 1 - call s_add_cloud_particle(cloud_idx, ib_idx, glbl_idx, geom, centroid(1), centroid(2), centroid(3), & - & cloud_ibs) + if (f_local_rank_owns_location(centroid, glb_bounds)) then + call s_add_cloud_particle(cloud_idx, ib_idx, glbl_idx, geom, centroid(1), centroid(2), & + & centroid(3), cloud_ibs) + end if n_placed = n_placed + 1 end do end do @@ -334,6 +351,30 @@ contains end subroutine s_particle_cloud_lattice + !> Compacts cloud_ibs(1:num_cloud_ibs) in place, discarding entries this rank does not own (f_local_rank_owns_location) and + !! updating num_cloud_ibs to the retained count. Used by rejection packing, which cannot filter as it places particles (see + !! s_particle_cloud_rejection_pack), to pare its full, unfiltered placement down to this rank's own subdomain. + subroutine s_reduce_particle_cloud_ibs(cloud_ibs, glb_bounds, num_cloud_ibs) + + type(ib_patch_parameters), intent(inout), dimension(:) :: cloud_ibs + type(bounds_info), dimension(3), intent(in) :: glb_bounds + integer, intent(inout) :: num_cloud_ibs + integer :: i, write_idx + real(wp), dimension(3) :: centroid + + write_idx = 0 + do i = 1, num_cloud_ibs + centroid = [cloud_ibs(i)%x_centroid, cloud_ibs(i)%y_centroid, 0._wp] + if (num_dims == 3) centroid(3) = cloud_ibs(i)%z_centroid + if (f_local_rank_owns_location(centroid, glb_bounds)) then + write_idx = write_idx + 1 + if (write_idx /= i) cloud_ibs(write_idx) = cloud_ibs(i) + end if + end do + num_cloud_ibs = write_idx + + end subroutine s_reduce_particle_cloud_ibs + !> Writes a single placed particle into particle_cloud_ibs at the next free slot, advancing ib_idx, tagged with its !! already-assigned, absolute global patch id via glbl_idx. Shared by all packing methods so the per-particle !! ib_patch_parameters setup stays in one place. Only x/y/z_centroid and radius are ever read back out of this struct (by @@ -348,7 +389,7 @@ contains ib_idx = ib_idx + 1 @:PROHIBIT(ib_idx > size(particle_cloud_ibs), & - & "Too many particle-cloud IBs. Modify case file or increase num_ib_patches_max_namelist.") + & "Too many particle-cloud IBs on one rank. Modify case file or increase num_ib_patches_max_namelist.") particle_cloud_ibs(ib_idx)%gbl_patch_id = glbl_idx particle_cloud_ibs(ib_idx)%geometry = geom diff --git a/src/pre_process/m_start_up.fpp b/src/pre_process/m_start_up.fpp index bbf0a2ddbf..bf9a14f1ba 100644 --- a/src/pre_process/m_start_up.fpp +++ b/src/pre_process/m_start_up.fpp @@ -40,7 +40,8 @@ module m_start_up private public :: s_read_input_file, s_check_input_file, s_read_grid_data_files, s_read_ic_data_files, s_read_serial_grid_data_files, & & s_read_serial_ic_data_files, s_read_parallel_grid_data_files, s_read_parallel_ic_data_files, s_check_grid_data_files, & - & s_initialize_modules, s_initialize_mpi_domain, s_finalize_modules, s_apply_initial_condition, s_save_data, s_read_grid + & s_initialize_modules, s_initialize_mpi_domain, s_finalize_modules, s_apply_initial_condition, s_save_data, s_read_grid, & + & s_write_ib_state_0 abstract interface @@ -134,24 +135,32 @@ contains call s_check_patches() - if (ib) then - call s_check_ib_patches() - - if (proc_rank == 0) then - block - type(ib_patch_parameters), allocatable :: particle_cloud_ibs(:) - integer :: num_particle_cloud_ibs - - call s_generate_particle_clouds(particle_cloud_ibs, num_particle_cloud_ibs) - call s_write_ib_state_0(particle_cloud_ibs, num_particle_cloud_ibs) - deallocate (particle_cloud_ibs) - end block - end if - call s_mpi_barrier() - end if + if (ib) call s_check_ib_patches() end subroutine s_check_input_file + !> @brief Generates the particle-cloud beds (if any) and writes restart_data/ib_state_0.dat: the initial IB layout that + !! simulation reads back at startup (src/simulation/m_start_up.fpp:s_read_ib_restart_data, called there with t_step = 0). Must + !! run after the computational domain is decomposed (s_initialize_mpi_domain) and the grid is populated (s_read_grid): every + !! rank computes the same deterministic candidate stream and keeps only the particles (and namelist patch_ib entries) + !! f_local_rank_owns_location says are its own, using this rank's x_cb/y_cb/z_cb - exactly mirroring how simulation's own + !! restart writers partition patch_ib, so file_per_process output stays consistent between the two executables. + impure subroutine s_write_ib_state_0() + + type(ib_patch_parameters), allocatable :: particle_cloud_ibs(:) + integer :: num_particle_cloud_ibs + type(bounds_info), dimension(3) :: glb_bounds + + if (.not. ib) return + + glb_bounds = (/x_domain_glb, y_domain_glb, z_domain_glb/) + + call s_generate_particle_clouds(glb_bounds, particle_cloud_ibs, num_particle_cloud_ibs) + call s_write_ib_state_0_file(glb_bounds, particle_cloud_ibs, num_particle_cloud_ibs) + deallocate (particle_cloud_ibs) + + end subroutine s_write_ib_state_0 + !> The goal of this subroutine is to read in any preexisting grid data as well as based on the imported grid, complete the !! necessary global computational domain parameters. impure subroutine s_read_serial_grid_data_files diff --git a/src/pre_process/p_main.f90 b/src/pre_process/p_main.f90 index e16ff68571..08cdc75766 100644 --- a/src/pre_process/p_main.f90 +++ b/src/pre_process/p_main.f90 @@ -24,6 +24,8 @@ program p_main call s_read_grid() + call s_write_ib_state_0() + allocate (proc_time(0:num_procs - 1)) call s_apply_initial_condition(start, finish) diff --git a/src/simulation/m_collisions.fpp b/src/simulation/m_collisions.fpp index 188b03977c..3cf7d311e8 100644 --- a/src/simulation/m_collisions.fpp +++ b/src/simulation/m_collisions.fpp @@ -21,7 +21,7 @@ module m_collisions implicit none private; public :: s_apply_collision_forces, s_initialize_collisions_module, s_finalize_collisions_module, & - & f_local_rank_owns_location, f_neighborhood_ranks_own_location, ib_gbl_idx_lookup + & f_neighborhood_ranks_own_location, ib_gbl_idx_lookup ! overlap distances for computing collisions integer, allocatable, dimension(:,:) :: collision_lookup real(wp), allocatable, dimension(:,:) :: wall_overlap_distances @@ -120,7 +120,7 @@ contains overlap_distance = patch_ib(pid1)%radius + patch_ib(pid2)%radius - norm2(normal_vector) if (overlap_distance > 0._wp) then ! if the two patches are close enough to collide normal_vector = normal_vector/norm2(normal_vector) - if (f_local_rank_owns_location(centroid_1)) then + if (f_local_rank_owns_location(centroid_1, glb_bounds)) then ! compute constants of the collision effective_mass = 1.0_wp/((1.0_wp/patch_ib(pid1)%mass) + (1._wp/(patch_ib(pid2)%mass))) k = spring_stiffness*effective_mass @@ -199,7 +199,7 @@ contains ! ensure the local rank owns that collision before proceeding collision_location = [patch_ib(patch_id)%x_centroid, patch_ib(patch_id)%y_centroid, 0._wp] if (num_dims == 3) collision_location(3) = patch_ib(patch_id)%z_centroid - if (f_local_rank_owns_location(collision_location)) then + if (f_local_rank_owns_location(collision_location, glb_bounds)) then k = spring_stiffness*patch_ib(patch_id)%mass eta = damping_parameter*patch_ib(patch_id)%mass @@ -372,41 +372,6 @@ contains end subroutine s_detect_wall_collisions - !> @brief function checks if this local MPI processor owns this specific collision - function f_local_rank_owns_location(location) result(owns_collision) - - $:GPU_ROUTINE(parallelism='[seq]') - - real(wp), dimension(3), intent(in) :: location - logical :: owns_collision - real(wp), dimension(3) :: projected_location - - owns_collision = .true. - -#ifdef MFC_MPI - if (num_procs > 1) then - projected_location(:) = location(:) - - ! catch the edge case where th collision lies just outside the computational domain - #:for X, ID, DIM in [('x', 1, 'm'), ('y', 2, 'n'), ('z', 3, 'p')] - if (num_dims >= ${ID}$) then - if (ib_bc_${X}$%beg /= BC_PERIODIC) then - ! if it is outside the domain in one direction, project it somewhere inside so at least one rank owns it - if (location(${ID}$) < glb_bounds(${ID}$)%beg) then - projected_location(${ID}$) = glb_bounds(${ID}$)%beg - else if (glb_bounds(${ID}$)%end < location(${ID}$)) then - projected_location(${ID}$) = glb_bounds(${ID}$)%end - 1.0e-10_wp - end if - end if - owns_collision = owns_collision .and. ${X}$_cb(-1) <= projected_location(${ID}$) & - & .and. projected_location(${ID}$) < ${X}$_cb(${DIM}$) - end if - #:endfor - end if -#endif - - end function f_local_rank_owns_location - !> @brief function checks if this local MPI processor owns this specific collision function f_neighborhood_ranks_own_location(location) result(owns_collision) diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp index 99a2122605..16e76e616d 100644 --- a/src/simulation/m_ibm.fpp +++ b/src/simulation/m_ibm.fpp @@ -1445,7 +1445,7 @@ contains end if ! check if in local domain - if (f_local_rank_owns_location(centroid)) then + if (f_local_rank_owns_location(centroid, glb_bounds)) then local_output_idx = local_output_idx + 1 local_ib_patch_ids(local_output_idx) = output_idx end if diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index 26199cd414..d9fd387411 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -1208,9 +1208,9 @@ contains !! (particle-cloud beds) are applied to particle_cloud_ibs(gbl_id - num_ibs) instead: patch_ib is GPU-declared at a fixed size !! and, elsewhere in this module, only ever written in its final, already-neighborhood-filtered form (by !! s_reduce_ib_patch_array) - bed entries stay in this plain staging array until that merge. Called with t_step = 0 to load the - !! initial layout pre_process generated (src/pre_process/m_data_output.fpp:s_write_ib_state_0) - since pre_process has no - !! file_per_process concept, t_step == 0 always uses the single-global-file layout below, regardless of this run's own - !! file_per_process setting - and again with the true step on an actual restart, overlaying live kinematics on top. + !! initial layout pre_process generated (src/pre_process/m_data_output.fpp:s_write_ib_state_0_file, which dispatches on the same + !! file_per_process flag and is kept format-compatible with this reader) - and again with the true step on an actual restart, + !! overlaying live kinematics on top. impure subroutine s_read_ib_restart_data(t_step, particle_cloud_ibs, num_particle_cloud_ibs) integer, intent(in) :: t_step @@ -1226,7 +1226,7 @@ contains n_records = num_ibs + num_particle_cloud_ibs - if (file_per_process .and. t_step /= 0) then + if (file_per_process) then call s_int_to_str(t_step, t_step_string) do r = 0, num_procs - 1 @@ -1390,7 +1390,7 @@ contains num_ibs = num_ibs + 1 patch_ib(num_ibs) = patch_ib(i) patch_ib(num_ibs)%gbl_patch_id = i - if (f_local_rank_owns_location(centroid)) then + if (f_local_rank_owns_location(centroid, glb_bounds)) then num_local_ibs = num_local_ibs + 1 @:PROHIBIT(num_local_ibs > num_local_ibs_max, & & "Too many IBs on a single processor rank. Modify case file or increase limit of num_local_ibs_max to resolve.") @@ -1408,7 +1408,7 @@ contains @:PROHIBIT(num_ibs > num_ib_patches_max_namelist, & & "Local IB count exceeds patch_ib capacity. Increase num_ib_patches_max_namelist.") patch_ib(num_ibs) = particle_cloud_ibs(i) - if (f_local_rank_owns_location(centroid)) then + if (f_local_rank_owns_location(centroid, glb_bounds)) then num_local_ibs = num_local_ibs + 1 @:PROHIBIT(num_local_ibs > num_local_ibs_max, & & "Too many IBs on a single processor rank. Modify case file or increase limit of num_local_ibs_max to resolve.") From 2675a363f8dbf9c87f3ff5ecfe8b72e6c06557c6 Mon Sep 17 00:00:00 2001 From: danieljvickers Date: Sun, 13 Sep 2026 01:42:40 -0400 Subject: [PATCH 3/5] Updated the files to have a file-per-process path that lets them read only the local ranks and then share. This also reduces the memroy reading required in exascale cases, which optimizes read times --- .claude/rules/common-pitfalls.md | 11 +- src/post_process/m_global_parameters.fpp | 2 +- src/pre_process/m_data_output.fpp | 209 +++-------- src/pre_process/m_mpi_proxy.fpp | 14 + src/pre_process/m_particle_cloud.fpp | 19 +- src/pre_process/m_start_up.fpp | 11 +- src/simulation/m_start_up.fpp | 427 ++++++++--------------- 7 files changed, 230 insertions(+), 463 deletions(-) diff --git a/.claude/rules/common-pitfalls.md b/.claude/rules/common-pitfalls.md index b3d7b0e47b..fdd5ead8b2 100644 --- a/.claude/rules/common-pitfalls.md +++ b/.claude/rules/common-pitfalls.md @@ -110,12 +110,11 @@ covered in `docs/documentation/contributing.md`. type silently drops their broadcast, so every non-root rank keeps the `dflt_real` sentinel; single-rank goldens cannot see this, so pair such a change with a `ppn=2` test and confirm it fails without the emitter. -- A `patch_ib` member that any `m_ibm` ghost-point code reads must ALSO be set in - `s_add_cloud_particle` (`src/simulation/m_particle_cloud.fpp`): `particle_cloud_ibs` is - `allocate`d without default initialization, and `s_reduce_ib_patch_array` copies the whole - struct into `patch_ib`, overwriting the defaults from - `s_assign_default_values_to_user_inputs`. Anything left unset reaches the solver as - uninitialized memory, and only where the allocation is not already zero-filled — a +- A `patch_ib` member that any `m_ibm` ghost-point code reads must ALSO be set for particle-cloud + IBs in `s_assign_particle_cloud_ib_defaults` (`src/simulation/m_start_up.fpp`): only position, + kinematics and radius come from the IB state file, and the routine writes into a reused + `patch_ib` slot, so anything it leaves unset keeps whatever that slot held (a namelist patch's + members, or memory that was never initialized). Such bugs show only where that memory is not already zero-filled — a 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 diff --git a/src/post_process/m_global_parameters.fpp b/src/post_process/m_global_parameters.fpp index 4c3a6712d9..a6b8052693 100644 --- a/src/post_process/m_global_parameters.fpp +++ b/src/post_process/m_global_parameters.fpp @@ -364,7 +364,7 @@ contains ! Particle clouds expand into individual IB patches at simulation startup, so num_ibs as read ! from the case file counts only the namelist patches. Match the global count the simulation - ! arrives at (s_reduce_ib_patch_array) so the IB state records can be read back. + ! arrives at (s_read_ib_restart_data) so the IB state records can be read back. do i = 1, num_particle_clouds num_ibs = num_ibs + particle_cloud(i)%num_particles end do diff --git a/src/pre_process/m_data_output.fpp b/src/pre_process/m_data_output.fpp index e3b4eb35fc..59717c5068 100644 --- a/src/pre_process/m_data_output.fpp +++ b/src/pre_process/m_data_output.fpp @@ -733,201 +733,74 @@ contains end subroutine s_initialize_data_output_module - !> @brief Writes restart_data/ib_state_0.dat (or, under file_per_process, one restart_data/lustre_0/ib_state_0_.dat chunk - !! per rank): the initial IB layout - namelist patch_ib entries, then any generated particle-cloud beds - with each rank writing - !! only the entries f_local_rank_owns_location says are its own. Read back by simulation at startup via s_read_ib_restart_data - !! at step 0 (src/simulation/m_start_up.fpp), which dispatches on the same file_per_process flag, so the two writers must stay - !! format-compatible: this mirrors s_write_parallel_ib_state/s_write_serial_ib_state (src/simulation/m_data_output.fpp) exactly, - !! just sourcing entries from this rank's local subset instead of patch_ib/ local_ib_patch_ids post-reduce. Only position - !! (fields 17:19) and radius (field 20) are populated here; everything else (time, force, torque, vel, angular_vel, angles) is - !! zero for a freshly generated IB. + !> @brief Writes the initial IB layout (namelist patch_ib entries, then generated particle-cloud beds) that simulation reads back + !! at startup (s_read_ib_restart_data, src/simulation/m_start_up.fpp), in the layouts simulation's own IB state writers use. + !! Under file_per_process each rank writes only the IBs it owns to restart_data/lustre_0/ib_state_0_.dat; otherwise only + !! rank 0 calls this and writes every IB, in global-id order, to restart_data/ib_state_0.dat. impure subroutine s_write_ib_state_0_file(glb_bounds, particle_cloud_ibs, num_particle_cloud_ibs) type(bounds_info), dimension(3), intent(in) :: glb_bounds type(ib_patch_parameters), dimension(:), intent(in) :: particle_cloud_ibs integer, intent(in) :: num_particle_cloud_ibs - integer, allocatable :: local_namelist_ids(:) - integer :: num_local_namelist - integer :: i + character(LEN=len_trim(case_dir) + 2*name_len) :: file_loc + integer :: i, ios, file_unit + logical, dimension(num_ibs) :: owned real(wp), dimension(3) :: centroid - allocate (local_namelist_ids(max(1, num_ibs))) - num_local_namelist = 0 - do i = 1, num_ibs - centroid = [patch_ib(i)%x_centroid, patch_ib(i)%y_centroid, 0._wp] - if (num_dims == 3) centroid(3) = patch_ib(i)%z_centroid - if (f_local_rank_owns_location(centroid, glb_bounds)) then - num_local_namelist = num_local_namelist + 1 - local_namelist_ids(num_local_namelist) = i - end if - end do + if (file_per_process) then + do i = 1, num_ibs + centroid = [patch_ib(i)%x_centroid, patch_ib(i)%y_centroid, patch_ib(i)%z_centroid] + owned(i) = f_local_rank_owns_location(centroid, glb_bounds) + end do - if (.not. parallel_io) then - ! Mirrors s_write_serial_ib_state (src/simulation/m_data_output.fpp): "must be called only on rank 0" - parallel_io - ! = F is not MPI-aware for ib_state, so this combined with num_procs > 1 is a pre-existing limitation, not one - ! introduced here. - if (proc_rank == 0) call s_write_ib_state_0_serial(local_namelist_ids, num_local_namelist, particle_cloud_ibs, & - & num_particle_cloud_ibs) - else if (file_per_process) then - call s_write_ib_state_0_file_per_process(local_namelist_ids, num_local_namelist, particle_cloud_ibs, & - & num_particle_cloud_ibs) + if (proc_rank == 0) call s_create_directory(trim(case_dir) // '/restart_data/lustre_0') + call s_mpi_barrier() + call s_delay_file_access(proc_rank) + write (file_loc, '(A,i7.7,A)') '/restart_data/lustre_0/ib_state_0_', proc_rank, '.dat' else - call s_write_ib_state_0_shared(local_namelist_ids, num_local_namelist, particle_cloud_ibs, num_particle_cloud_ibs) + owned = .true. + call s_create_directory(trim(case_dir) // '/restart_data') + file_loc = '/restart_data/ib_state_0.dat' end if - - deallocate (local_namelist_ids) - - end subroutine s_write_ib_state_0_file - - !> Writes this rank's local entries into its own restart_data/lustre_0/ib_state_0_.dat chunk file - mirrors - !! s_write_parallel_ib_state's file_per_process branch (src/simulation/m_data_output.fpp). - subroutine s_write_ib_state_0_file_per_process(local_namelist_ids, num_local_namelist, particle_cloud_ibs, & - & num_particle_cloud_ibs) - - integer, dimension(:), intent(in) :: local_namelist_ids - integer, intent(in) :: num_local_namelist - type(ib_patch_parameters), dimension(:), intent(in) :: particle_cloud_ibs - integer, intent(in) :: num_particle_cloud_ibs - character(LEN=len_trim(case_dir) + 2*name_len) :: file_loc - integer :: i, ios, file_unit - integer, parameter :: NFIELDS_PER_IB = 20 - real(wp) :: ib_buf(NFIELDS_PER_IB) - - if (proc_rank == 0) call s_create_directory(trim(case_dir) // '/restart_data/lustre_0') - call s_mpi_barrier() - call s_delay_file_access(proc_rank) - - write (file_loc, '(A,i7.7,A)') 'ib_state_0_', proc_rank, '.dat' - file_loc = trim(case_dir) // '/restart_data/lustre_0/' // trim(file_loc) + file_loc = trim(case_dir) // trim(file_loc) open (newunit=file_unit, file=trim(file_loc), form='unformatted', access='stream', status='replace', iostat=ios) if (ios /= 0) call s_mpi_abort('Cannot open IB state output file: ' // trim(file_loc)) - ib_buf = 0._wp - - write (file_unit) num_local_namelist + num_particle_cloud_ibs - - do i = 1, num_local_namelist - ib_buf(17) = patch_ib(local_namelist_ids(i))%x_centroid - ib_buf(18) = patch_ib(local_namelist_ids(i))%y_centroid - ib_buf(19) = patch_ib(local_namelist_ids(i))%z_centroid - ib_buf(20) = patch_ib(local_namelist_ids(i))%radius - write (file_unit) local_namelist_ids(i) - write (file_unit) ib_buf + if (file_per_process) write (file_unit) count(owned) + num_particle_cloud_ibs + do i = 1, num_ibs + if (owned(i)) call s_write_ib_state_record(patch_ib(i), i) end do - do i = 1, num_particle_cloud_ibs - ib_buf(17) = particle_cloud_ibs(i)%x_centroid - ib_buf(18) = particle_cloud_ibs(i)%y_centroid - ib_buf(19) = particle_cloud_ibs(i)%z_centroid - ib_buf(20) = particle_cloud_ibs(i)%radius - write (file_unit) particle_cloud_ibs(i)%gbl_patch_id - write (file_unit) ib_buf + call s_write_ib_state_record(particle_cloud_ibs(i), particle_cloud_ibs(i)%gbl_patch_id) end do close (file_unit) - end subroutine s_write_ib_state_0_file_per_process - - !> Writes this rank's local entries into the single shared restart_data/ib_state_0.dat, each rank placing its own records at - !! their gbl_patch_id-based offset via MPI-IO - mirrors s_write_parallel_ib_state's non-file_per_process branch - !! (src/simulation/m_data_output.fpp). Only reached when parallel_io = T (see s_write_ib_state_0_file), matching that routine's - !! own precondition of running under MFC_MPI. - subroutine s_write_ib_state_0_shared(local_namelist_ids, num_local_namelist, particle_cloud_ibs, num_particle_cloud_ibs) - - integer, dimension(:), intent(in) :: local_namelist_ids - integer, intent(in) :: num_local_namelist - type(ib_patch_parameters), dimension(:), intent(in) :: particle_cloud_ibs - integer, intent(in) :: num_particle_cloud_ibs - integer, parameter :: NFIELDS_PER_IB = 20 - real(wp) :: ib_buf(NFIELDS_PER_IB) - integer :: i - -#ifdef MFC_MPI - character(LEN=len_trim(case_dir) + 2*name_len) :: file_loc - integer(kind=MPI_OFFSET_KIND) :: disp, WP_MOK - integer :: ifile, ierr - integer, dimension(MPI_STATUS_SIZE) :: status - logical :: file_exist - - if (proc_rank == 0) call s_create_directory(trim(case_dir) // '/restart_data') - call s_mpi_barrier() - - WP_MOK = int(storage_size(0._wp)/8, MPI_OFFSET_KIND) - file_loc = trim(case_dir) // '/restart_data/ib_state_0.dat' - - inquire (FILE=trim(file_loc), EXIST=file_exist) - if (file_exist .and. proc_rank == 0) call MPI_FILE_DELETE(file_loc, mpi_info_int, ierr) - call s_mpi_barrier() - - call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, ior(MPI_MODE_WRONLY, MPI_MODE_CREATE), mpi_info_int, ifile, ierr) - - ib_buf = 0._wp - - do i = 1, num_local_namelist - ib_buf(17) = patch_ib(local_namelist_ids(i))%x_centroid - ib_buf(18) = patch_ib(local_namelist_ids(i))%y_centroid - ib_buf(19) = patch_ib(local_namelist_ids(i))%z_centroid - ib_buf(20) = patch_ib(local_namelist_ids(i))%radius - disp = int(local_namelist_ids(i) - 1, MPI_OFFSET_KIND)*int(NFIELDS_PER_IB, MPI_OFFSET_KIND)*WP_MOK - call MPI_FILE_WRITE_AT(ifile, disp, ib_buf, NFIELDS_PER_IB, mpi_p, status, ierr) - end do - - do i = 1, num_particle_cloud_ibs - ib_buf(17) = particle_cloud_ibs(i)%x_centroid - ib_buf(18) = particle_cloud_ibs(i)%y_centroid - ib_buf(19) = particle_cloud_ibs(i)%z_centroid - ib_buf(20) = particle_cloud_ibs(i)%radius - disp = int(particle_cloud_ibs(i)%gbl_patch_id - 1, MPI_OFFSET_KIND)*int(NFIELDS_PER_IB, MPI_OFFSET_KIND)*WP_MOK - call MPI_FILE_WRITE_AT(ifile, disp, ib_buf, NFIELDS_PER_IB, mpi_p, status, ierr) - end do - - call MPI_FILE_CLOSE(ifile, ierr) -#endif - - end subroutine s_write_ib_state_0_shared - - !> Writes ALL of this rank's local entries directly into restart_data/ib_state_0.dat with plain Fortran I/O - mirrors - !! s_write_serial_ib_state (src/simulation/m_data_output.fpp). Used only when parallel_io = F; the caller (rank 0 only) mirrors - !! that routine's "must be called only on rank 0" precondition, so this is a plain sequential write, no MPI needed. - subroutine s_write_ib_state_0_serial(local_namelist_ids, num_local_namelist, particle_cloud_ibs, num_particle_cloud_ibs) - - integer, dimension(:), intent(in) :: local_namelist_ids - integer, intent(in) :: num_local_namelist - type(ib_patch_parameters), dimension(:), intent(in) :: particle_cloud_ibs - integer, intent(in) :: num_particle_cloud_ibs - character(LEN=len_trim(case_dir) + 2*name_len) :: file_loc - integer :: i, ios, file_unit - integer, parameter :: NFIELDS_PER_IB = 20 - real(wp) :: ib_buf(NFIELDS_PER_IB) + contains - call s_create_directory(trim(case_dir) // '/restart_data') - file_loc = trim(case_dir) // '/restart_data/ib_state_0.dat' + !> Writes one 20-field IB state record, prefixed by its global id under file_per_process. + subroutine s_write_ib_state_record(ib_patch, gbl_id) - open (newunit=file_unit, file=trim(file_loc), form='unformatted', access='stream', status='replace', iostat=ios) - if (ios /= 0) call s_mpi_abort('Cannot open IB state output file: ' // trim(file_loc)) + type(ib_patch_parameters), intent(in) :: ib_patch + integer, intent(in) :: gbl_id + real(wp), dimension(20) :: ib_buf - ib_buf = 0._wp + ib_buf = 0._wp + ib_buf(8:10) = ib_patch%vel + ib_buf(11:13) = ib_patch%angular_vel + ib_buf(14:16) = ib_patch%angles + ib_buf(17) = ib_patch%x_centroid + ib_buf(18) = ib_patch%y_centroid + ib_buf(19) = ib_patch%z_centroid + ib_buf(20) = ib_patch%radius - do i = 1, num_local_namelist - ib_buf(17) = patch_ib(local_namelist_ids(i))%x_centroid - ib_buf(18) = patch_ib(local_namelist_ids(i))%y_centroid - ib_buf(19) = patch_ib(local_namelist_ids(i))%z_centroid - ib_buf(20) = patch_ib(local_namelist_ids(i))%radius + if (file_per_process) write (file_unit) gbl_id write (file_unit) ib_buf - end do - do i = 1, num_particle_cloud_ibs - ib_buf(17) = particle_cloud_ibs(i)%x_centroid - ib_buf(18) = particle_cloud_ibs(i)%y_centroid - ib_buf(19) = particle_cloud_ibs(i)%z_centroid - ib_buf(20) = particle_cloud_ibs(i)%radius - write (file_unit) ib_buf - end do + end subroutine s_write_ib_state_record - close (file_unit) - - end subroutine s_write_ib_state_0_serial + end subroutine s_write_ib_state_0_file !> Resets s_write_data_files pointer impure subroutine s_finalize_data_output_module diff --git a/src/pre_process/m_mpi_proxy.fpp b/src/pre_process/m_mpi_proxy.fpp index 45c1a3b505..826af94eb4 100644 --- a/src/pre_process/m_mpi_proxy.fpp +++ b/src/pre_process/m_mpi_proxy.fpp @@ -135,6 +135,20 @@ contains call MPI_BCAST(patch_ib(i)%slip, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr) end do + ! manual: particle_cloud (runtime loop to num_particle_clouds; irregular member subset) + do i = 1, num_particle_clouds + #:for VAR in ['x_centroid', 'y_centroid', 'z_centroid', 'length_x', 'length_y', 'length_z', & + & 'radius', 'mass', 'min_spacing', 'shell_inner_radius', 'shell_outer_radius'] + call MPI_BCAST(particle_cloud(i)%${VAR}$, 1, mpi_p, 0, MPI_COMM_WORLD, ierr) + #:endfor + call MPI_BCAST(particle_cloud(i)%num_particles, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) + call MPI_BCAST(particle_cloud(i)%moving_ibm, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) + call MPI_BCAST(particle_cloud(i)%seed, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) + call MPI_BCAST(particle_cloud(i)%cloud_geometry, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) + call MPI_BCAST(particle_cloud(i)%packing_method, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) + call MPI_BCAST(particle_cloud(i)%periodic, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) + end do + ! manual: ib_airfoil (kept manual alongside patch_ib) do i = 1, num_ib_airfoils_max #:for VAR in ['c', 'p', 't', 'm'] diff --git a/src/pre_process/m_particle_cloud.fpp b/src/pre_process/m_particle_cloud.fpp index 37bc7c1896..acd1cb3038 100644 --- a/src/pre_process/m_particle_cloud.fpp +++ b/src/pre_process/m_particle_cloud.fpp @@ -7,10 +7,9 @@ #:include 'macros.fpp' !> @brief Generates particle beds by converting particle_cloud patch specifications into individual immersed boundary patches before -!! writing them to restart_data/ib_state_0.dat. Runs on every rank: each rank computes the same deterministic placement (so no MPI -!! broadcast of particle positions is needed) and keeps only the particles f_local_rank_owns_location says are its own, so every -!! generated particle is written by exactly one rank - consistent with how simulation partitions patch_ib for its own restart writes -!! (src/simulation/m_start_up.fpp:s_reduce_ib_patch_array). +!! writing them to the initial IB state file. Under file_per_process it runs on every rank: each rank computes the same deterministic +!! placement (so no MPI broadcast of particle positions is needed) and keeps only the particles f_local_rank_owns_location says are +!! its own, so every generated particle is written by exactly one rank. Otherwise only rank 0 runs it and keeps every particle. module m_particle_cloud use m_global_parameters @@ -26,9 +25,9 @@ module m_particle_cloud contains - !> Generate all particle beds and fill particle_cloud_ibs, keeping only the particles this rank owns (see module docs). Called - !! on every rank. Each packing method owns and allocates its own per-cloud working array (see s_particle_cloud_lattice / - !! s_particle_cloud_rejection_pack) and hands back only the entries that fall within this rank's own subdomain. Only the first + !> Generate all particle beds and fill particle_cloud_ibs, keeping only the particles this rank owns under file_per_process (see + !! module docs). Each packing method owns and allocates its own per-cloud working array (see s_particle_cloud_lattice / + !! s_particle_cloud_rejection_pack) and hands back only the entries this rank keeps. Only the first !! num_particle_cloud_ibs of them are actually written - callers must use that count, not size(particle_cloud_ibs), since the !! remainder of the array is left uninitialized. impure subroutine s_generate_particle_clouds(glb_bounds, particle_cloud_ibs, num_particle_cloud_ibs) @@ -180,7 +179,7 @@ contains deallocate (placed, hash_head, chain_next) - call s_reduce_particle_cloud_ibs(cloud_ibs, glb_bounds, ib_idx) + if (file_per_process) call s_reduce_particle_cloud_ibs(cloud_ibs, glb_bounds, ib_idx) num_cloud_ibs = ib_idx end subroutine s_particle_cloud_rejection_pack @@ -308,7 +307,7 @@ contains do while (px <= xmax .and. n_placed < n_target) glbl_idx = glbl_idx + 1 centroid = [px, py, particle_cloud(cloud_idx)%z_centroid] - if (f_local_rank_owns_location(centroid, glb_bounds)) then + if (.not. file_per_process .or. f_local_rank_owns_location(centroid, glb_bounds)) then call s_add_cloud_particle(cloud_idx, ib_idx, glbl_idx, geom, centroid(1), centroid(2), centroid(3), & & cloud_ibs) end if @@ -335,7 +334,7 @@ contains centroid = [xmin + real(ix, wp)*cell + bx_off(b), ymin + real(jy, wp)*cell + by_off(b), & & zmin + real(kz, wp)*cell + bz_off(b)] glbl_idx = glbl_idx + 1 - if (f_local_rank_owns_location(centroid, glb_bounds)) then + if (.not. file_per_process .or. f_local_rank_owns_location(centroid, glb_bounds)) then call s_add_cloud_particle(cloud_idx, ib_idx, glbl_idx, geom, centroid(1), centroid(2), & & centroid(3), cloud_ibs) end if diff --git a/src/pre_process/m_start_up.fpp b/src/pre_process/m_start_up.fpp index bf9a14f1ba..5c96b5ef37 100644 --- a/src/pre_process/m_start_up.fpp +++ b/src/pre_process/m_start_up.fpp @@ -139,12 +139,10 @@ contains end subroutine s_check_input_file - !> @brief Generates the particle-cloud beds (if any) and writes restart_data/ib_state_0.dat: the initial IB layout that - !! simulation reads back at startup (src/simulation/m_start_up.fpp:s_read_ib_restart_data, called there with t_step = 0). Must - !! run after the computational domain is decomposed (s_initialize_mpi_domain) and the grid is populated (s_read_grid): every - !! rank computes the same deterministic candidate stream and keeps only the particles (and namelist patch_ib entries) - !! f_local_rank_owns_location says are its own, using this rank's x_cb/y_cb/z_cb - exactly mirroring how simulation's own - !! restart writers partition patch_ib, so file_per_process output stays consistent between the two executables. + !> @brief Generates the particle-cloud beds (if any) and writes the initial IB state file that simulation reads back at startup + !! (src/simulation/m_start_up.fpp:s_read_ib_restart_data). Must run after the domain is decomposed (s_initialize_mpi_domain) and + !! the grid is populated (s_read_grid). Under file_per_process every rank computes the same deterministic placement and keeps + !! only the IBs f_local_rank_owns_location says are its own; otherwise rank 0 alone generates and writes every IB. impure subroutine s_write_ib_state_0() type(ib_patch_parameters), allocatable :: particle_cloud_ibs(:) @@ -152,6 +150,7 @@ contains type(bounds_info), dimension(3) :: glb_bounds if (.not. ib) return + if (.not. file_per_process .and. proc_rank /= 0) return glb_bounds = (/x_domain_glb, y_domain_glb, z_domain_glb/) diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index 795e0db561..185f53adc7 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -912,24 +912,17 @@ contains if (model_eqns == model_eqns_6eq) call s_initialize_internal_energy_equations(q_cons_ts(1)%vf) if (ib) then - block - type(ib_patch_parameters), allocatable :: particle_cloud_ibs(:) - integer :: num_particle_cloud_ibs - - call s_instantiate_STL_models() - call s_initialize_ib_airfoils() - call s_get_neighbor_bounds() - - call s_assign_particle_cloud_ib_defaults(particle_cloud_ibs, num_particle_cloud_ibs) - call s_read_ib_restart_data(0, particle_cloud_ibs, num_particle_cloud_ibs) - if (cfl_dt .and. n_start > 0) then - call s_read_ib_restart_data(n_start, particle_cloud_ibs, num_particle_cloud_ibs) - else if (t_step_start > 0) then - call s_read_ib_restart_data(t_step_start, particle_cloud_ibs, num_particle_cloud_ibs) - end if - call s_reduce_ib_patch_array(particle_cloud_ibs, num_particle_cloud_ibs) - deallocate (particle_cloud_ibs) - end block + call s_instantiate_STL_models() + call s_initialize_ib_airfoils() + call s_get_neighbor_bounds() + if (cfl_dt .and. n_start > 0) then + call s_read_ib_restart_data(n_start) + else if (t_step_start > 0) then + call s_read_ib_restart_data(t_step_start) + else + call s_read_ib_restart_data(0) + end if + call s_build_ib_neighborhood() call s_ibm_setup() if (t_step_start == 0 .or. (cfl_dt .and. n_start == 0)) then call s_write_ib_data_file(0) @@ -1155,298 +1148,188 @@ contains end subroutine s_finalize_modules - !> @brief Allocates particle_cloud_ibs to sum(particle_cloud(:)%num_particles) and fills the properties that are constant per - !! cloud (geometry, mass, moving_ibm, slip, v_blow, inj_species, burn rate, identity rotation matrix, zeroed kinematics) - - !! mirroring the defaults s_add_cloud_particle (src/pre_process/m_particle_cloud.fpp) sets at generation time, in the same - !! cloud-then-particle-index order, so entry ib_idx here corresponds to global id num_ibs + ib_idx exactly as pre_process - !! assigned it. Position and radius are left at zero; s_read_ib_restart_data fills them in from restart_data/ib_state_0.dat. - impure subroutine s_assign_particle_cloud_ib_defaults(particle_cloud_ibs, num_particle_cloud_ibs) + !> @brief Fills the properties of a generated particle-cloud IB that are constant per cloud (geometry, mass, moving_ibm, inert + !! surface, identity rotation matrix, zeroed kinematics), mirroring s_add_cloud_particle (src/pre_process/m_particle_cloud.fpp). + !! cloud_ib_idx is the global patch id minus the number of namelist patches; pre_process numbers particles cloud by cloud. + subroutine s_assign_particle_cloud_ib_defaults(cloud_ib_idx, ib_patch) - type(ib_patch_parameters), allocatable, intent(out), dimension(:) :: particle_cloud_ibs - integer, intent(out) :: num_particle_cloud_ibs - integer :: cloud_idx, i, ib_idx, geom + integer, intent(in) :: cloud_ib_idx + type(ib_patch_parameters), intent(inout) :: ib_patch + integer :: cloud_idx, idx_in_cloud - num_particle_cloud_ibs = 0 - do cloud_idx = 1, num_particle_clouds - num_particle_cloud_ibs = num_particle_cloud_ibs + particle_cloud(cloud_idx)%num_particles - end do - allocate (particle_cloud_ibs(num_particle_cloud_ibs)) - - geom = merge(2, 8, num_dims < 3) - - ib_idx = 0 - do cloud_idx = 1, num_particle_clouds - do i = 1, particle_cloud(cloud_idx)%num_particles - ib_idx = ib_idx + 1 - - particle_cloud_ibs(ib_idx)%gbl_patch_id = num_ibs + ib_idx - particle_cloud_ibs(ib_idx)%geometry = geom - particle_cloud_ibs(ib_idx)%x_centroid = 0._wp - particle_cloud_ibs(ib_idx)%y_centroid = 0._wp - particle_cloud_ibs(ib_idx)%z_centroid = 0._wp - particle_cloud_ibs(ib_idx)%step_x_centroid = 0._wp - particle_cloud_ibs(ib_idx)%step_y_centroid = 0._wp - particle_cloud_ibs(ib_idx)%step_z_centroid = 0._wp - particle_cloud_ibs(ib_idx)%angles(:) = 0._wp - particle_cloud_ibs(ib_idx)%step_angles(:) = 0._wp - particle_cloud_ibs(ib_idx)%vel(:) = 0._wp - particle_cloud_ibs(ib_idx)%step_vel(:) = 0._wp - particle_cloud_ibs(ib_idx)%angular_vel(:) = 0._wp - particle_cloud_ibs(ib_idx)%step_angular_vel(:) = 0._wp - particle_cloud_ibs(ib_idx)%force(:) = 0._wp - particle_cloud_ibs(ib_idx)%torque(:) = 0._wp - particle_cloud_ibs(ib_idx)%centroid_offset(:) = 0._wp - particle_cloud_ibs(ib_idx)%rotation_matrix = 0._wp - particle_cloud_ibs(ib_idx)%rotation_matrix(1, 1) = 1._wp - particle_cloud_ibs(ib_idx)%rotation_matrix(2, 2) = 1._wp - particle_cloud_ibs(ib_idx)%rotation_matrix(3, 3) = 1._wp - particle_cloud_ibs(ib_idx)%rotation_matrix_inverse = particle_cloud_ibs(ib_idx)%rotation_matrix - particle_cloud_ibs(ib_idx)%radius = 0._wp - particle_cloud_ibs(ib_idx)%mass = particle_cloud(cloud_idx)%mass - particle_cloud_ibs(ib_idx)%moment = dflt_real - particle_cloud_ibs(ib_idx)%moving_ibm = particle_cloud(cloud_idx)%moving_ibm - particle_cloud_ibs(ib_idx)%slip = .false. - particle_cloud_ibs(ib_idx)%v_blow = 0._wp - particle_cloud_ibs(ib_idx)%inj_species = 0 - particle_cloud_ibs(ib_idx)%burn_rate_exp = 0._wp - particle_cloud_ibs(ib_idx)%burn_rate_pref = 0._wp - end do + idx_in_cloud = cloud_ib_idx + do cloud_idx = 1, num_particle_clouds - 1 + if (idx_in_cloud <= particle_cloud(cloud_idx)%num_particles) exit + idx_in_cloud = idx_in_cloud - particle_cloud(cloud_idx)%num_particles end do + ib_patch%geometry = merge(2, 8, num_dims < 3) + ib_patch%step_x_centroid = 0._wp + ib_patch%step_y_centroid = 0._wp + ib_patch%step_z_centroid = 0._wp + ib_patch%step_angles(:) = 0._wp + ib_patch%step_vel(:) = 0._wp + ib_patch%step_angular_vel(:) = 0._wp + ib_patch%force(:) = 0._wp + ib_patch%torque(:) = 0._wp + ib_patch%centroid_offset(:) = 0._wp + ib_patch%rotation_matrix = 0._wp + ib_patch%rotation_matrix(1, 1) = 1._wp + ib_patch%rotation_matrix(2, 2) = 1._wp + ib_patch%rotation_matrix(3, 3) = 1._wp + ib_patch%rotation_matrix_inverse = ib_patch%rotation_matrix + ib_patch%mass = particle_cloud(cloud_idx)%mass + ib_patch%moment = dflt_real + ib_patch%moving_ibm = particle_cloud(cloud_idx)%moving_ibm + ib_patch%slip = .false. + ib_patch%v_blow = 0._wp + ib_patch%inj_species = 0 + ib_patch%burn_rate_exp = 0._wp + ib_patch%burn_rate_pref = 0._wp + end subroutine s_assign_particle_cloud_ib_defaults - !> @brief Reads IB kinematic state (plus, for particle-cloud beds, radius) from restart_data/ib_state_.dat. Records with - !! gbl_id <= num_ibs (namelist patches) are applied directly to patch_ib, exactly as before. Records with gbl_id > num_ibs - !! (particle-cloud beds) are applied to particle_cloud_ibs(gbl_id - num_ibs) instead: patch_ib is GPU-declared at a fixed size - !! and, elsewhere in this module, only ever written in its final, already-neighborhood-filtered form (by - !! s_reduce_ib_patch_array) - bed entries stay in this plain staging array until that merge. Called with t_step = 0 to load the - !! initial layout pre_process generated (src/pre_process/m_data_output.fpp:s_write_ib_state_0_file, which dispatches on the same - !! file_per_process flag and is kept format-compatible with this reader) - and again with the true step on an actual restart, - !! overlaying live kinematics on top. - impure subroutine s_read_ib_restart_data(t_step, particle_cloud_ibs, num_particle_cloud_ibs) - - integer, intent(in) :: t_step - type(ib_patch_parameters), dimension(:), intent(inout) :: particle_cloud_ibs - integer, intent(in) :: num_particle_cloud_ibs - character(len=path_len + 2*name_len) :: file_loc - integer :: i, ios, file_unit, ierr - integer :: r, nlocal, gbl_id, n_records - integer, parameter :: NFIELDS_PER_IB = 20 - real(wp) :: ib_buf(NFIELDS_PER_IB) - logical :: file_exist - character(len=10) :: t_step_string - - n_records = num_ibs + num_particle_cloud_ibs + !> @brief Loads the IBs this rank owns from the IB state file for t_step into patch_ib(1:num_ibs), all of which are local. + !! Under file_per_process the rank reads only its own restart_data/lustre_/ib_state__.dat, which holds + !! exactly its IBs; otherwise every rank reads every record of restart_data/ib_state_.dat and keeps the ones + !! f_local_rank_owns_location assigns it. Records carry kinematics, position and radius; every other property comes from the + !! namelist patch (global id <= num_ibs) or the particle cloud the IB was generated from. Written by pre_process at t_step = 0 + !! (src/pre_process/m_data_output.fpp:s_write_ib_state_0_file) and by s_write_ib_state_file on later steps. + impure subroutine s_read_ib_restart_data(t_step) + + integer, intent(in) :: t_step + type(ib_patch_parameters), allocatable, dimension(:) :: namelist_ibs + character(len=path_len + 2*name_len) :: file_loc + integer :: i, ios, file_unit, gbl_id, n_records + integer, parameter :: NFIELDS_PER_IB = 20 + real(wp) :: ib_buf(NFIELDS_PER_IB) + character(len=10) :: t_step_string + + moving_immersed_boundary_flag = any(patch_ib(1:num_ibs)%moving_ibm /= 0) & + & .or. any(particle_cloud(1:num_particle_clouds)%moving_ibm /= 0) + + allocate (namelist_ibs(num_ibs)) + namelist_ibs(:) = patch_ib(1:num_ibs) + num_gbl_ibs = num_ibs + sum(particle_cloud(1:num_particle_clouds)%num_particles) if (file_per_process) then call s_int_to_str(t_step, t_step_string) - - do r = 0, num_procs - 1 - write (file_loc, '(A,I0,A,i7.7,A)') 'ib_state_', t_step, '_', r, '.dat' - file_loc = trim(case_dir) // '/restart_data/lustre_' // trim(t_step_string) // '/' // trim(file_loc) - - inquire (FILE=trim(file_loc), EXIST=file_exist) - if (.not. file_exist) call s_mpi_abort('Cannot open IB state file for restart: ' // trim(file_loc)) - - open (newunit=file_unit, file=trim(file_loc), form='unformatted', access='stream', status='old', iostat=ios) - if (ios /= 0) call s_mpi_abort('Error opening IB state restart file: ' // trim(file_loc)) - - read (file_unit, iostat=ios) nlocal - if (ios /= 0) call s_mpi_abort('Error reading IB state file header: ' // trim(file_loc)) - - do i = 1, nlocal - read (file_unit, iostat=ios) gbl_id - if (ios /= 0) call s_mpi_abort('Error reading IB patch ID: ' // trim(file_loc)) - read (file_unit, iostat=ios) ib_buf - if (ios /= 0) call s_mpi_abort('Error reading IB state data: ' // trim(file_loc)) - - if (gbl_id <= num_ibs) then - call s_apply_ib_buf(patch_ib(gbl_id), ib_buf) - else - call s_apply_ib_buf(particle_cloud_ibs(gbl_id - num_ibs), ib_buf) - end if - end do - - close (file_unit) - end do + write (file_loc, '(A,I0,A,i7.7,A)') 'ib_state_', t_step, '_', proc_rank, '.dat' + file_loc = trim(case_dir) // '/restart_data/lustre_' // trim(t_step_string) // '/' // trim(file_loc) else write (file_loc, '(A,I0,A)') '/restart_data/ib_state_', t_step, '.dat' file_loc = trim(case_dir) // trim(file_loc) + end if - if (proc_rank == 0) then - inquire (FILE=trim(file_loc), EXIST=file_exist) - if (.not. file_exist) then - call s_mpi_abort('Cannot open IB state file for restart: ' // trim(file_loc)) - end if - - open (newunit=file_unit, file=trim(file_loc), form='unformatted', access='stream', status='old', iostat=ios) - if (ios /= 0) call s_mpi_abort('Error opening IB state restart file: ' // trim(file_loc)) - - do i = 1, n_records - read (file_unit, iostat=ios) ib_buf - if (ios /= 0) call s_mpi_abort('Error reading IB state restart file') - - if (i <= num_ibs) then - call s_apply_ib_buf(patch_ib(i), ib_buf) - else - call s_apply_ib_buf(particle_cloud_ibs(i - num_ibs), ib_buf) - end if - end do - - close (file_unit) - end if + open (newunit=file_unit, file=trim(file_loc), form='unformatted', access='stream', status='old', action='read', iostat=ios) + if (ios /= 0) call s_mpi_abort('Error opening IB state file: ' // trim(file_loc)) -#ifdef MFC_MPI - do i = 1, num_ibs - call MPI_BCAST(patch_ib(i)%vel, 3, mpi_p, 0, MPI_COMM_WORLD, ierr) - call MPI_BCAST(patch_ib(i)%angular_vel, 3, mpi_p, 0, MPI_COMM_WORLD, ierr) - call MPI_BCAST(patch_ib(i)%angles, 3, mpi_p, 0, MPI_COMM_WORLD, ierr) - call MPI_BCAST(patch_ib(i)%x_centroid, 1, mpi_p, 0, MPI_COMM_WORLD, ierr) - call MPI_BCAST(patch_ib(i)%y_centroid, 1, mpi_p, 0, MPI_COMM_WORLD, ierr) - call MPI_BCAST(patch_ib(i)%z_centroid, 1, mpi_p, 0, MPI_COMM_WORLD, ierr) - call MPI_BCAST(patch_ib(i)%radius, 1, mpi_p, 0, MPI_COMM_WORLD, ierr) - end do - do i = 1, num_particle_cloud_ibs - call MPI_BCAST(particle_cloud_ibs(i)%vel, 3, mpi_p, 0, MPI_COMM_WORLD, ierr) - call MPI_BCAST(particle_cloud_ibs(i)%angular_vel, 3, mpi_p, 0, MPI_COMM_WORLD, ierr) - call MPI_BCAST(particle_cloud_ibs(i)%angles, 3, mpi_p, 0, MPI_COMM_WORLD, ierr) - call MPI_BCAST(particle_cloud_ibs(i)%x_centroid, 1, mpi_p, 0, MPI_COMM_WORLD, ierr) - call MPI_BCAST(particle_cloud_ibs(i)%y_centroid, 1, mpi_p, 0, MPI_COMM_WORLD, ierr) - call MPI_BCAST(particle_cloud_ibs(i)%z_centroid, 1, mpi_p, 0, MPI_COMM_WORLD, ierr) - call MPI_BCAST(particle_cloud_ibs(i)%radius, 1, mpi_p, 0, MPI_COMM_WORLD, ierr) - end do -#endif + n_records = num_gbl_ibs + if (file_per_process) then + read (file_unit, iostat=ios) n_records + if (ios /= 0) call s_mpi_abort('Error reading IB state file header: ' // trim(file_loc)) end if - contains + num_ibs = 0 + do i = 1, n_records + gbl_id = i + if (file_per_process) read (file_unit, iostat=ios) gbl_id + if (ios == 0) read (file_unit, iostat=ios) ib_buf + if (ios /= 0) call s_mpi_abort('Error reading IB state file: ' // trim(file_loc)) - !> Applies one 20-field IB state record to a single patch_ib/particle_cloud_ibs entry. - subroutine s_apply_ib_buf(entry, buf) + if (.not. file_per_process) then + if (.not. f_local_rank_owns_location(ib_buf(17:19), glb_bounds)) cycle + end if - type(ib_patch_parameters), intent(inout) :: entry - real(wp), dimension(NFIELDS_PER_IB), intent(in) :: buf + num_ibs = num_ibs + 1 + @:PROHIBIT(num_ibs > num_local_ibs_max, & + & "Too many IBs on a single processor rank. Modify case file or increase limit of num_local_ibs_max to resolve.") + if (gbl_id <= size(namelist_ibs)) then + patch_ib(num_ibs) = namelist_ibs(gbl_id) + else + call s_assign_particle_cloud_ib_defaults(gbl_id - size(namelist_ibs), patch_ib(num_ibs)) + end if + patch_ib(num_ibs)%gbl_patch_id = gbl_id + patch_ib(num_ibs)%vel = ib_buf(8:10) + patch_ib(num_ibs)%angular_vel = ib_buf(11:13) + patch_ib(num_ibs)%angles = ib_buf(14:16) + patch_ib(num_ibs)%x_centroid = ib_buf(17) + patch_ib(num_ibs)%y_centroid = ib_buf(18) + patch_ib(num_ibs)%z_centroid = ib_buf(19) + patch_ib(num_ibs)%radius = ib_buf(20) + local_ib_patch_ids(num_ibs) = num_ibs + end do - entry%vel = buf(8:10) - entry%angular_vel = buf(11:13) - entry%angles = buf(14:16) - entry%x_centroid = buf(17) - entry%y_centroid = buf(18) - entry%z_centroid = buf(19) - entry%radius = buf(20) + close (file_unit) + deallocate (namelist_ibs) - end subroutine s_apply_ib_buf + num_local_ibs = num_ibs end subroutine s_read_ib_restart_data - !> @brief Merges patch_ib (namelist patches, fixed at num_ib_patches_max_namelist) with particle_cloud_ibs (the full, global, - !! unfiltered bed particle set - see s_assign_particle_cloud_ib_defaults/s_read_ib_restart_data - each entry already tagged with - !! its final, absolute gbl_patch_id) and reduces to only the patches in or near the local computational domain. patch_ib is - !! never reallocated; the local subset is written in-place from the front. particle_cloud_ibs is owned by the caller and freed - !! there after this returns. - subroutine s_reduce_ib_patch_array(particle_cloud_ibs, num_particle_cloud_ibs) - - type(ib_patch_parameters), intent(in), dimension(:) :: particle_cloud_ibs - integer, intent(in) :: num_particle_cloud_ibs - real(wp), dimension(3) :: centroid - integer :: i - integer :: num_namelist_ibs, num_bed_ibs - - num_namelist_ibs = num_ibs - num_bed_ibs = 0 - do i = 1, num_particle_clouds - num_bed_ibs = num_bed_ibs + particle_cloud(i)%num_particles - end do + !> @brief Completes this rank's IB neighborhood once s_read_ib_restart_data has loaded only the IBs it owns: every rank sends its + !! own IBs to, and receives the owned IBs of, each distinct rank in ib_neighbor_ranks, appending them to patch_ib after its own. + subroutine s_build_ib_neighborhood() - ! Check for moving IBs across both namelist and particle cloud patches. - moving_immersed_boundary_flag = .false. - do i = 1, num_namelist_ibs - if (patch_ib(i)%moving_ibm /= 0) then - moving_immersed_boundary_flag = .true. - exit - end if - end do - if (.not. moving_immersed_boundary_flag) then - do i = 1, num_particle_clouds - if (particle_cloud(i)%moving_ibm /= 0) then - moving_immersed_boundary_flag = .true. - exit - end if - end do - end if +#ifdef MFC_MPI + integer, allocatable, dimension(:) :: nbr_ranks, recv_counts, requests + type(ib_patch_parameters), allocatable, dimension(:,:) :: recv_ibs + integer :: i, n_nbrs, nreqs, patch_bytes, ierr +#endif call s_compute_ib_neighbor_ranks() #ifdef MFC_MPI - if (num_procs == 1) then - ! single-rank: all patches are local; append particle bed entries directly into patch_ib. - do i = 1, num_particle_cloud_ibs - patch_ib(num_namelist_ibs + i) = particle_cloud_ibs(i) + if (num_procs > 1) then + ! A rank can fill several table slots (periodicity, few ranks) or be its own neighbor; exchange once per distinct rank + allocate (nbr_ranks(size(ib_neighbor_ranks))) + nbr_ranks = reshape(ib_neighbor_ranks, [size(ib_neighbor_ranks)]) + n_nbrs = 0 + do i = 1, size(nbr_ranks) + if (nbr_ranks(i) < 0 .or. nbr_ranks(i) == proc_rank) cycle + if (any(nbr_ranks(1:n_nbrs) == nbr_ranks(i))) cycle + n_nbrs = n_nbrs + 1 + nbr_ranks(n_nbrs) = nbr_ranks(i) end do - num_gbl_ibs = num_namelist_ibs + num_particle_cloud_ibs - @:PROHIBIT(num_gbl_ibs > num_ib_patches_max_namelist, & - & "Total IB count exceeds patch_ib capacity. Increase num_ib_patches_max_namelist.") - num_ibs = num_gbl_ibs - num_local_ibs = num_gbl_ibs - do i = 1, num_gbl_ibs - local_ib_patch_ids(i) = i + + allocate (recv_counts(n_nbrs), requests(2*n_nbrs)) + do i = 1, n_nbrs + call MPI_IRECV(recv_counts(i), 1, MPI_INTEGER, nbr_ranks(i), 500, MPI_COMM_WORLD, requests(2*i - 1), ierr) + call MPI_ISEND(num_local_ibs, 1, MPI_INTEGER, nbr_ranks(i), 500, MPI_COMM_WORLD, requests(2*i), ierr) end do - else - ! multi-rank: compact namelist patches in-place (write_idx <= read_idx, no aliasing), then append local particle beds. - num_ibs = 0 - num_local_ibs = 0 - num_gbl_ibs = num_namelist_ibs + num_bed_ibs - do i = 1, num_namelist_ibs - centroid = [patch_ib(i)%x_centroid, patch_ib(i)%y_centroid, 0._wp] - if (num_dims == 3) centroid(3) = patch_ib(i)%z_centroid - if (f_neighborhood_ranks_own_location(centroid)) then - num_ibs = num_ibs + 1 - patch_ib(num_ibs) = patch_ib(i) - patch_ib(num_ibs)%gbl_patch_id = i - if (f_local_rank_owns_location(centroid, glb_bounds)) then - num_local_ibs = num_local_ibs + 1 - @:PROHIBIT(num_local_ibs > num_local_ibs_max, & - & "Too many IBs on a single processor rank. Modify case file or increase limit of num_local_ibs_max to resolve.") - local_ib_patch_ids(num_local_ibs) = num_ibs - end if + call MPI_WAITALL(2*n_nbrs, requests, MPI_STATUSES_IGNORE, ierr) + + patch_bytes = storage_size(patch_ib(1))/8 + allocate (recv_ibs(max(1, maxval(recv_counts)), n_nbrs)) + nreqs = 0 + do i = 1, n_nbrs + if (recv_counts(i) > 0) then + nreqs = nreqs + 1 + call MPI_IRECV(recv_ibs(:,i), recv_counts(i)*patch_bytes, MPI_BYTE, nbr_ranks(i), 501, MPI_COMM_WORLD, & + & requests(nreqs), ierr) end if - end do - ! particle_cloud_ibs holds the full global bed set; filter to this rank's neighborhood exactly as above for namelist - ! patches. - do i = 1, num_particle_cloud_ibs - centroid = [particle_cloud_ibs(i)%x_centroid, particle_cloud_ibs(i)%y_centroid, 0._wp] - if (num_dims == 3) centroid(3) = particle_cloud_ibs(i)%z_centroid - if (f_neighborhood_ranks_own_location(centroid)) then - num_ibs = num_ibs + 1 - @:PROHIBIT(num_ibs > num_ib_patches_max_namelist, & - & "Local IB count exceeds patch_ib capacity. Increase num_ib_patches_max_namelist.") - patch_ib(num_ibs) = particle_cloud_ibs(i) - if (f_local_rank_owns_location(centroid, glb_bounds)) then - num_local_ibs = num_local_ibs + 1 - @:PROHIBIT(num_local_ibs > num_local_ibs_max, & - & "Too many IBs on a single processor rank. Modify case file or increase limit of num_local_ibs_max to resolve.") - local_ib_patch_ids(num_local_ibs) = num_ibs - end if + if (num_local_ibs > 0) then + nreqs = nreqs + 1 + call MPI_ISEND(patch_ib, num_local_ibs*patch_bytes, MPI_BYTE, nbr_ranks(i), 501, MPI_COMM_WORLD, & + & requests(nreqs), ierr) end if end do + call MPI_WAITALL(nreqs, requests, MPI_STATUSES_IGNORE, ierr) + + do i = 1, n_nbrs + @:PROHIBIT(num_ibs + recv_counts(i) > num_ib_patches_max_namelist, & + & "IB neighborhood exceeds patch_ib capacity. Increase num_ib_patches_max_namelist.") + patch_ib(num_ibs + 1:num_ibs + recv_counts(i)) = recv_ibs(1:recv_counts(i), i) + num_ibs = num_ibs + recv_counts(i) + end do + + deallocate (nbr_ranks, recv_counts, requests, recv_ibs) end if -#else - ! no-MPI: all patches are local; append particle bed entries directly into patch_ib. - do i = 1, num_particle_cloud_ibs - patch_ib(num_namelist_ibs + i) = particle_cloud_ibs(i) - end do - num_gbl_ibs = num_namelist_ibs + num_particle_cloud_ibs - @:PROHIBIT(num_gbl_ibs > num_ib_patches_max_namelist, & - & "Total IB count exceeds patch_ib capacity. Increase num_ib_patches_max_namelist.") - num_ibs = num_gbl_ibs - num_local_ibs = num_gbl_ibs - do i = 1, num_gbl_ibs - local_ib_patch_ids(i) = i - end do #endif @:ALLOCATE(ib_gbl_idx_lookup(1:num_gbl_ibs)) - end subroutine s_reduce_ib_patch_array + end subroutine s_build_ib_neighborhood !> Build ib_neighbor_ranks(-1:1,-1:1,-1:1): MPI ranks of all neighbor domains. Uses two rounds of MPI_SENDRECV cascades - face !! neighbors are known from bc_*, edge neighbors are obtained in round 1, and (3D) corner neighbors in round 2. From 87ce523ce79d03fb5f39b2392dcbf99e0a390763 Mon Sep 17 00:00:00 2001 From: danieljvickers Date: Sun, 13 Sep 2026 01:52:57 -0400 Subject: [PATCH 4/5] Addressed review commetns --- src/common/m_helper.fpp | 4 +- src/post_process/m_global_parameters.fpp | 1 + src/post_process/m_start_up.fpp | 6 ++ src/pre_process/m_data_output.fpp | 8 +-- src/pre_process/m_global_parameters.fpp | 3 +- src/pre_process/m_particle_cloud.fpp | 73 ++++++------------------ src/pre_process/m_start_up.fpp | 6 ++ src/simulation/m_start_up.fpp | 31 +++++----- toolchain/mfc/case_validator.py | 2 + 9 files changed, 56 insertions(+), 78 deletions(-) diff --git a/src/common/m_helper.fpp b/src/common/m_helper.fpp index 5a046c769f..17a986d6a9 100644 --- a/src/common/m_helper.fpp +++ b/src/common/m_helper.fpp @@ -10,7 +10,7 @@ module m_helper use m_derived_types use m_global_parameters - use m_constants + use m_constants, only: BC_PERIODIC use ieee_arithmetic !< For checking NaN implicit none @@ -725,7 +725,7 @@ contains ! catch the edge case where the location lies just outside the computational domain #:for X, ID, DIM in [('x', 1, 'm'), ('y', 2, 'n'), ('z', 3, 'p')] if (num_dims >= ${ID}$) then - if (bc_${X}$%beg /= BC_PERIODIC) then + if (ib_bc_${X}$%beg /= BC_PERIODIC) then ! if it is outside the domain in one direction, project it somewhere inside so at least one rank owns it if (location(${ID}$) < glb_bounds_in(${ID}$)%beg) then projected_location(${ID}$) = glb_bounds_in(${ID}$)%beg diff --git a/src/post_process/m_global_parameters.fpp b/src/post_process/m_global_parameters.fpp index a6b8052693..05bed7f795 100644 --- a/src/post_process/m_global_parameters.fpp +++ b/src/post_process/m_global_parameters.fpp @@ -93,6 +93,7 @@ module m_global_parameters !> @name Boundary conditions in the x-, y- and z-coordinate directions !> @{ type(int_bounds_info) :: bc_x, bc_y, bc_z + type(int_bounds_info) :: ib_bc_x, ib_bc_y, ib_bc_z !< bc_x/y/z before decomposition overwrites them with MPI neighbor ranks type(bc_xyz_info) :: bc !> @} diff --git a/src/post_process/m_start_up.fpp b/src/post_process/m_start_up.fpp index 86a00c0018..c97133fff3 100644 --- a/src/post_process/m_start_up.fpp +++ b/src/post_process/m_start_up.fpp @@ -952,6 +952,12 @@ contains end if call s_mpi_bcast_user_inputs() + + ! Save original BCs before decomposition overwrites them with MPI neighbor ranks + ib_bc_x = bc_x + ib_bc_y = bc_y + ib_bc_z = bc_z + call s_initialize_parallel_io() 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., & diff --git a/src/pre_process/m_data_output.fpp b/src/pre_process/m_data_output.fpp index 59717c5068..85e9300ed4 100644 --- a/src/pre_process/m_data_output.fpp +++ b/src/pre_process/m_data_output.fpp @@ -733,10 +733,10 @@ contains end subroutine s_initialize_data_output_module - !> @brief Writes the initial IB layout (namelist patch_ib entries, then generated particle-cloud beds) that simulation reads back - !! at startup (s_read_ib_restart_data, src/simulation/m_start_up.fpp), in the layouts simulation's own IB state writers use. - !! Under file_per_process each rank writes only the IBs it owns to restart_data/lustre_0/ib_state_0_.dat; otherwise only - !! rank 0 calls this and writes every IB, in global-id order, to restart_data/ib_state_0.dat. + !> @brief Writes the initial IB layout (namelist patch_ib entries, then generated particle-cloud beds) that simulation reads + !! back at startup (s_read_ib_restart_data, src/simulation/m_start_up.fpp), in the layouts simulation's own IB state writers + !! use. Under file_per_process each rank writes only the IBs it owns to restart_data/lustre_0/ib_state_0_.dat; otherwise + !! only rank 0 calls this and writes every IB, in global-id order, to restart_data/ib_state_0.dat. impure subroutine s_write_ib_state_0_file(glb_bounds, particle_cloud_ibs, num_particle_cloud_ibs) type(bounds_info), dimension(3), intent(in) :: glb_bounds diff --git a/src/pre_process/m_global_parameters.fpp b/src/pre_process/m_global_parameters.fpp index a3d7b835ee..551f1432e7 100644 --- a/src/pre_process/m_global_parameters.fpp +++ b/src/pre_process/m_global_parameters.fpp @@ -52,7 +52,8 @@ module m_global_parameters ! Cell indices (InDices With BUFFer): includes buffer except in pre_process type(int_bounds_info) :: idwbuff(1:3) type(int_bounds_info) :: bc_x, bc_y, bc_z !< Boundary conditions in the x-, y- and z-coordinate directions - type(bc_xyz_info) :: bc !< Combined BC storage (used by the shared beta-buffer routines; pre-process-local) + type(int_bounds_info) :: ib_bc_x, ib_bc_y, ib_bc_z !< bc_x/y/z before decomposition overwrites them with MPI neighbor ranks + type(bc_xyz_info) :: bc !< Combined BC storage (used by the shared beta-buffer routines; pre-process-local) ! simplex_params: auto-generated in generated_decls.fpp ! shear_num/shear_indices/shear_BC_flip_*, bc: in m_global_parameters_common integer :: fd_order !< Finite-difference order for CoM/probe derivative approximations diff --git a/src/pre_process/m_particle_cloud.fpp b/src/pre_process/m_particle_cloud.fpp index acd1cb3038..fba6452ca9 100644 --- a/src/pre_process/m_particle_cloud.fpp +++ b/src/pre_process/m_particle_cloud.fpp @@ -7,9 +7,10 @@ #:include 'macros.fpp' !> @brief Generates particle beds by converting particle_cloud patch specifications into individual immersed boundary patches before -!! writing them to the initial IB state file. Under file_per_process it runs on every rank: each rank computes the same deterministic -!! placement (so no MPI broadcast of particle positions is needed) and keeps only the particles f_local_rank_owns_location says are -!! its own, so every generated particle is written by exactly one rank. Otherwise only rank 0 runs it and keeps every particle. +!! writing them to the initial IB state file. Under file_per_process it runs on every rank: each rank computes the same +!! deterministic placement (so no MPI broadcast of particle positions is needed) and keeps only the particles +!! f_local_rank_owns_location says are its own, so every generated particle is written by exactly one rank. Otherwise only rank 0 +!! runs it and keeps every particle. module m_particle_cloud use m_global_parameters @@ -27,9 +28,9 @@ contains !> Generate all particle beds and fill particle_cloud_ibs, keeping only the particles this rank owns under file_per_process (see !! module docs). Each packing method owns and allocates its own per-cloud working array (see s_particle_cloud_lattice / - !! s_particle_cloud_rejection_pack) and hands back only the entries this rank keeps. Only the first - !! num_particle_cloud_ibs of them are actually written - callers must use that count, not size(particle_cloud_ibs), since the - !! remainder of the array is left uninitialized. + !! s_particle_cloud_rejection_pack) and hands back only the entries this rank keeps. Only the first num_particle_cloud_ibs of + !! them are actually written - callers must use that count, not size(particle_cloud_ibs), since the remainder of the array is + !! left uninitialized. impure subroutine s_generate_particle_clouds(glb_bounds, particle_cloud_ibs, num_particle_cloud_ibs) type(bounds_info), dimension(3), intent(in) :: glb_bounds @@ -95,7 +96,7 @@ contains type(bounds_info), dimension(3), intent(in) :: glb_bounds type(ib_patch_parameters), allocatable, intent(out), dimension(:) :: cloud_ibs integer, intent(out) :: num_cloud_ibs - integer :: ib_idx, n_placed, geom, seed, alloc_stat + integer :: ib_idx, n_placed, seed, alloc_stat integer(8) :: n_attempts, max_attempts real(wp) :: min_dist, rx, ry, rz logical :: overlaps, reject, periodic_pack @@ -125,12 +126,6 @@ contains nz_bins = max(1, ceiling(length_z/min_dist)) if (num_dims < 3) nz_bins = 1 - if (num_dims < 3) then - geom = 2 ! circle for 2D - else - geom = 8 ! sphere for 3D - end if - max_attempts = int(particle_cloud(cloud_idx)%num_particles, 8)*1000_8 n_placed = 0 n_attempts = 0 @@ -169,7 +164,7 @@ contains hash_head(slot) = n_placed glbl_idx = glbl_idx + 1 - call s_add_cloud_particle(cloud_idx, ib_idx, glbl_idx, geom, rx, ry, rz, cloud_ibs) + call s_add_cloud_particle(cloud_idx, ib_idx, glbl_idx, rx, ry, rz, cloud_ibs) end if end do @@ -258,7 +253,7 @@ contains type(bounds_info), dimension(3), intent(in) :: glb_bounds type(ib_patch_parameters), allocatable, intent(out), dimension(:) :: cloud_ibs integer, intent(out) :: num_cloud_ibs - integer :: ib_idx, n_placed, n_target, geom + integer :: ib_idx, n_placed, n_target integer :: row, col, ncx, ncy, ix, jy, kz, b real(wp) :: xmin, xmax, ymin, ymax, zmin, zmax, min_dist real(wp) :: spacing, row_dy, cell, x0, px, py @@ -280,11 +275,9 @@ contains n_placed = 0 if (num_dims < 3) then - geom = 2 ! circle for 2D ! Triangular lattice: area per particle = (sqrt(3)/2)*spacing**2. spacing = sqrt(2._wp*(xmax - xmin)*(ymax - ymin)/(sqrt(3._wp)*real(n_target, wp))) else - geom = 8 ! sphere for 3D ! Face-centered cubic lattice: volume per particle = spacing**3/sqrt(2). spacing = (sqrt(2._wp)*(xmax - xmin)*(ymax - ymin)*(zmax - zmin)/real(n_target, wp))**(1._wp/3._wp) end if @@ -308,8 +301,7 @@ contains glbl_idx = glbl_idx + 1 centroid = [px, py, particle_cloud(cloud_idx)%z_centroid] if (.not. file_per_process .or. f_local_rank_owns_location(centroid, glb_bounds)) then - call s_add_cloud_particle(cloud_idx, ib_idx, glbl_idx, geom, centroid(1), centroid(2), centroid(3), & - & cloud_ibs) + call s_add_cloud_particle(cloud_idx, ib_idx, glbl_idx, centroid(1), centroid(2), centroid(3), cloud_ibs) end if n_placed = n_placed + 1 col = col + 1 @@ -335,8 +327,8 @@ contains & zmin + real(kz, wp)*cell + bz_off(b)] glbl_idx = glbl_idx + 1 if (.not. file_per_process .or. f_local_rank_owns_location(centroid, glb_bounds)) then - call s_add_cloud_particle(cloud_idx, ib_idx, glbl_idx, geom, centroid(1), centroid(2), & - & centroid(3), cloud_ibs) + call s_add_cloud_particle(cloud_idx, ib_idx, glbl_idx, centroid(1), centroid(2), centroid(3), & + & cloud_ibs) end if n_placed = n_placed + 1 end do @@ -375,13 +367,11 @@ contains end subroutine s_reduce_particle_cloud_ibs !> Writes a single placed particle into particle_cloud_ibs at the next free slot, advancing ib_idx, tagged with its - !! already-assigned, absolute global patch id via glbl_idx. Shared by all packing methods so the per-particle - !! ib_patch_parameters setup stays in one place. Only x/y/z_centroid and radius are ever read back out of this struct (by - !! s_write_ib_state_0); the rest is filled in for completeness and because leaving fields uninitialized here would otherwise - !! surface as garbage if this struct is ever read further downstream. - subroutine s_add_cloud_particle(cloud_idx, ib_idx, glbl_idx, geom, px, py, pz, particle_cloud_ibs) + !! already-assigned, absolute global patch id via glbl_idx. Only the fields s_write_ib_state_0_file writes are set; simulation + !! fills every other property in s_assign_particle_cloud_ib_defaults (src/simulation/m_start_up.fpp). + subroutine s_add_cloud_particle(cloud_idx, ib_idx, glbl_idx, px, py, pz, particle_cloud_ibs) - integer, intent(in) :: cloud_idx, glbl_idx, geom + integer, intent(in) :: cloud_idx, glbl_idx integer, intent(inout) :: ib_idx real(wp), intent(in) :: px, py, pz type(ib_patch_parameters), intent(inout), dimension(:) :: particle_cloud_ibs @@ -391,40 +381,13 @@ contains & "Too many particle-cloud IBs on one rank. Modify case file or increase num_ib_patches_max_namelist.") particle_cloud_ibs(ib_idx)%gbl_patch_id = glbl_idx - particle_cloud_ibs(ib_idx)%geometry = geom particle_cloud_ibs(ib_idx)%x_centroid = px particle_cloud_ibs(ib_idx)%y_centroid = py particle_cloud_ibs(ib_idx)%z_centroid = pz - particle_cloud_ibs(ib_idx)%step_x_centroid = px - particle_cloud_ibs(ib_idx)%step_y_centroid = py - particle_cloud_ibs(ib_idx)%step_z_centroid = pz - particle_cloud_ibs(ib_idx)%angles(:) = 0._wp - particle_cloud_ibs(ib_idx)%step_angles(:) = 0._wp particle_cloud_ibs(ib_idx)%vel(:) = 0._wp - particle_cloud_ibs(ib_idx)%step_vel(:) = 0._wp particle_cloud_ibs(ib_idx)%angular_vel(:) = 0._wp - particle_cloud_ibs(ib_idx)%step_angular_vel(:) = 0._wp - particle_cloud_ibs(ib_idx)%force(:) = 0._wp - particle_cloud_ibs(ib_idx)%torque(:) = 0._wp - particle_cloud_ibs(ib_idx)%centroid_offset(:) = 0._wp - particle_cloud_ibs(ib_idx)%rotation_matrix = 0._wp - particle_cloud_ibs(ib_idx)%rotation_matrix(1, 1) = 1._wp - particle_cloud_ibs(ib_idx)%rotation_matrix(2, 2) = 1._wp - particle_cloud_ibs(ib_idx)%rotation_matrix(3, 3) = 1._wp - particle_cloud_ibs(ib_idx)%rotation_matrix_inverse = particle_cloud_ibs(ib_idx)%rotation_matrix + particle_cloud_ibs(ib_idx)%angles(:) = 0._wp particle_cloud_ibs(ib_idx)%radius = particle_cloud(cloud_idx)%radius - particle_cloud_ibs(ib_idx)%mass = particle_cloud(cloud_idx)%mass - particle_cloud_ibs(ib_idx)%moment = dflt_real - particle_cloud_ibs(ib_idx)%moving_ibm = particle_cloud(cloud_idx)%moving_ibm - particle_cloud_ibs(ib_idx)%slip = .false. - - ! Particles are inert surfaces. particle_cloud_ibs is allocated (not default-initialized), - ! so these must be set explicitly even though only x/y/z_centroid and radius are read back - ! out of this struct downstream. - particle_cloud_ibs(ib_idx)%v_blow = 0._wp - particle_cloud_ibs(ib_idx)%inj_species = 0 - particle_cloud_ibs(ib_idx)%burn_rate_exp = 0._wp - particle_cloud_ibs(ib_idx)%burn_rate_pref = 0._wp end subroutine s_add_cloud_particle diff --git a/src/pre_process/m_start_up.fpp b/src/pre_process/m_start_up.fpp index 5c96b5ef37..fd21f95d97 100644 --- a/src/pre_process/m_start_up.fpp +++ b/src/pre_process/m_start_up.fpp @@ -646,6 +646,12 @@ contains ! Broadcasting the user inputs to all of the processors and performing the parallel computational domain decomposition. ! Neither procedure has to be carried out if pre-process is in fact not truly executed in parallel. call s_mpi_bcast_user_inputs() + + ! Save original BCs before decomposition overwrites them with MPI neighbor ranks + ib_bc_x = bc_x + ib_bc_y = bc_y + ib_bc_z = bc_z + call s_initialize_parallel_io() ! Save the global domain bounds before decomposition overwrites x/y/z_domain with each processor's local sub-domain bounds diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index 185f53adc7..6c9ef4a72b 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -924,10 +924,7 @@ contains end if call s_build_ib_neighborhood() call s_ibm_setup() - if (t_step_start == 0 .or. (cfl_dt .and. n_start == 0)) then - call s_write_ib_data_file(0) - call s_write_ib_state_file(0) - end if + if (t_step_start == 0 .or. (cfl_dt .and. n_start == 0)) call s_write_ib_data_file(0) end if if (bodyForces .or. synthetic_turbulence) call s_initialize_body_forces_module() if (acoustic_source) call s_precalculate_acoustic_spatial_sources() @@ -1148,9 +1145,10 @@ contains end subroutine s_finalize_modules - !> @brief Fills the properties of a generated particle-cloud IB that are constant per cloud (geometry, mass, moving_ibm, inert - !! surface, identity rotation matrix, zeroed kinematics), mirroring s_add_cloud_particle (src/pre_process/m_particle_cloud.fpp). - !! cloud_ib_idx is the global patch id minus the number of namelist patches; pre_process numbers particles cloud by cloud. + !> @brief Fills the properties of a generated particle-cloud IB that the IB state file does not carry (geometry, mass, + !! moving_ibm, inert surface, identity rotation matrix, zeroed step state). This is the only place they are set - pre_process + !! (s_add_cloud_particle) generates only position, kinematics and radius. cloud_ib_idx is the global patch id minus the number + !! of namelist patches; pre_process numbers particles cloud by cloud. subroutine s_assign_particle_cloud_ib_defaults(cloud_ib_idx, ib_patch) integer, intent(in) :: cloud_ib_idx @@ -1189,9 +1187,9 @@ contains end subroutine s_assign_particle_cloud_ib_defaults - !> @brief Loads the IBs this rank owns from the IB state file for t_step into patch_ib(1:num_ibs), all of which are local. - !! Under file_per_process the rank reads only its own restart_data/lustre_/ib_state__.dat, which holds - !! exactly its IBs; otherwise every rank reads every record of restart_data/ib_state_.dat and keeps the ones + !> @brief Loads the IBs this rank owns from the IB state file for t_step into patch_ib(1:num_ibs), all of which are local. Under + !! file_per_process the rank reads only its own restart_data/lustre_/ib_state__.dat, which holds exactly + !! its IBs; otherwise every rank reads every record of restart_data/ib_state_.dat and keeps the ones !! f_local_rank_owns_location assigns it. Records carry kinematics, position and radius; every other property comes from the !! namelist patch (global id <= num_ibs) or the particle cloud the IB was generated from. Written by pre_process at t_step = 0 !! (src/pre_process/m_data_output.fpp:s_write_ib_state_0_file) and by s_write_ib_state_file on later steps. @@ -1206,7 +1204,7 @@ contains character(len=10) :: t_step_string moving_immersed_boundary_flag = any(patch_ib(1:num_ibs)%moving_ibm /= 0) & - & .or. any(particle_cloud(1:num_particle_clouds)%moving_ibm /= 0) + & .or. any(particle_cloud(1:num_particle_clouds)%moving_ibm /= 0) allocate (namelist_ibs(num_ibs)) namelist_ibs(:) = patch_ib(1:num_ibs) @@ -1267,14 +1265,15 @@ contains end subroutine s_read_ib_restart_data - !> @brief Completes this rank's IB neighborhood once s_read_ib_restart_data has loaded only the IBs it owns: every rank sends its - !! own IBs to, and receives the owned IBs of, each distinct rank in ib_neighbor_ranks, appending them to patch_ib after its own. + !> @brief Completes this rank's IB neighborhood once s_read_ib_restart_data has loaded only the IBs it owns: every rank sends + !! its own IBs to, and receives the owned IBs of, each distinct rank in ib_neighbor_ranks, appending them to patch_ib after its + !! own. subroutine s_build_ib_neighborhood() #ifdef MFC_MPI - integer, allocatable, dimension(:) :: nbr_ranks, recv_counts, requests + integer, allocatable, dimension(:) :: nbr_ranks, recv_counts, requests type(ib_patch_parameters), allocatable, dimension(:,:) :: recv_ibs - integer :: i, n_nbrs, nreqs, patch_bytes, ierr + integer :: i, n_nbrs, nreqs, patch_bytes, ierr #endif call s_compute_ib_neighbor_ranks() @@ -1319,7 +1318,7 @@ contains do i = 1, n_nbrs @:PROHIBIT(num_ibs + recv_counts(i) > num_ib_patches_max_namelist, & & "IB neighborhood exceeds patch_ib capacity. Increase num_ib_patches_max_namelist.") - patch_ib(num_ibs + 1:num_ibs + recv_counts(i)) = recv_ibs(1:recv_counts(i), i) + patch_ib(num_ibs + 1:num_ibs + recv_counts(i)) = recv_ibs(1:recv_counts(i),i) num_ibs = num_ibs + recv_counts(i) end do diff --git a/toolchain/mfc/case_validator.py b/toolchain/mfc/case_validator.py index 1da790a90e..c048f8029a 100644 --- a/toolchain/mfc/case_validator.py +++ b/toolchain/mfc/case_validator.py @@ -1872,6 +1872,8 @@ def check_parallel_io_pre_process(self): m = self.get("m", 0) n = self.get("n", 0) + self.prohibit(file_per_process and not parallel_io, "file_per_process requires parallel_io = T") + if down_sample: self.prohibit(not parallel_io, "down sample requires parallel_io = T") self.prohibit(not igr, "down sample requires igr = T") From 120b94e44a63f042ab646f9180876097ada1f9c5 Mon Sep 17 00:00:00 2001 From: danieljvickers Date: Sun, 13 Sep 2026 10:05:33 -0400 Subject: [PATCH 5/5] Update docs --- docs/documentation/contributing.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/documentation/contributing.md b/docs/documentation/contributing.md index dd3c0f33c6..366ea1d844 100644 --- a/docs/documentation/contributing.md +++ b/docs/documentation/contributing.md @@ -218,7 +218,7 @@ Both human reviewers and AI code reviewers reference this section. ### Parameter Plumbing - **Derived-type parameters are not auto-broadcast.** `generated_bcast.fpp` covers namelist *scalars* only. Each derived type (`chem_params`, `lag_params`, `rburn`) needs a hand-written `_emit_` in `toolchain/mfc/params/generators/fortran_gen.py` plus its call site in that generator's simulation branch, and, if it is read on device, an explicit ``$:GPU_UPDATE(device='[name]')`` in both the target's `m_global_parameters.fpp` and `src/simulation/m_start_up.fpp` — `GPU_DECLARE` alone does not make it device-resident. Regrouping existing scalars into a derived type silently drops their broadcast, leaving every non-root rank holding the `dflt_real` sentinel. Single-rank golden files cannot catch this, so pair such a change with a `ppn=2` test and confirm it fails without the emitter. -- **A `patch_ib` member that immersed-boundary ghost-point code reads must also be set in `s_add_cloud_particle`** (`src/simulation/m_particle_cloud.fpp`). `particle_cloud_ibs` is allocated without default initialization, and `s_reduce_ib_patch_array` copies the whole struct into `patch_ib`, overwriting the defaults assigned in `s_assign_default_values_to_user_inputs`. Anything left unset reaches the solver as uninitialized memory, and only where the allocation is not already zero-filled. A platform-only NaN is the signature of this class: a garbage `v_blow` once failed an AMD lane with `ICFL is NaN` while every NVIDIA lane and all local runs passed. +- **A `patch_ib` member that immersed-boundary ghost-point code reads must also be set for particle-cloud IBs in `s_assign_particle_cloud_ib_defaults`** (`src/simulation/m_start_up.fpp`). Pre-process writes only position, kinematics and radius to the IB state file; simulation builds every other property there, writing into a reused `patch_ib` slot. Anything it leaves unset keeps whatever that slot held, which may be a namelist patch's value or uninitialized memory, and shows up only where that memory is not already zero-filled. A platform-only NaN is the signature of this class: a garbage `v_blow` once failed an AMD lane with `ICFL is NaN` while every NVIDIA lane and all local runs passed. - **Runtime checks go where they run.** Shared constraints belong in `src/common/m_checker_common.fpp`, simulation-only ones in `src/simulation/m_checker.fpp`, and pre- and post-process ones in their own `m_checker.fpp`. Those two `s_check_inputs` are currently empty; that is still the correct home for their checks, not `m_checker_common`. - **Analytic initial conditions are compiled into the binary** and their expressions are AST-validated at case load, so syntax errors and unknown variables surface immediately and by name. Each IC variable maps to an `eqn_idx` expression in `QPVF_IDX_VARS` (`toolchain/mfc/case.py`); adding a patch-settable conserved variable means updating that map and the Fortran `eqn_idx` builder together, because a mismatch is a silent wrong index. - **Under `--case-optimization` the baked-in constants are dropped from the namelist**, so changing one requires a rebuild rather than a case-file edit.