From 04cb847effbdfc29800c258d2046f5418721dc15 Mon Sep 17 00:00:00 2001 From: Daniel Vickers Date: Wed, 8 Jul 2026 19:45:00 -0400 Subject: [PATCH 1/8] Swapped the particle bed generation to be neighborhood aware during instantiation, saving memroy. --- src/common/m_constants.fpp | 6 +-- src/simulation/m_collisions.fpp | 4 +- src/simulation/m_particle_cloud.fpp | 58 ++++++++++++++++++----------- src/simulation/m_start_up.fpp | 55 +++++++++++++++------------ 4 files changed, 74 insertions(+), 49 deletions(-) diff --git a/src/common/m_constants.fpp b/src/common/m_constants.fpp index 4857d93249..ed67683917 100644 --- a/src/common/m_constants.fpp +++ b/src/common/m_constants.fpp @@ -27,10 +27,10 @@ module m_constants integer, parameter :: num_ib_airfoils_max = 5 !< Maximum number of ib_airfoil instances integer, parameter :: num_stl_models_max = 10 !> Maximum number of immersed boundary patches (legacy, not used for patch_ib sizing) - integer, parameter :: num_ib_patches_max = 2050000 + integer, parameter :: num_ib_patches_max = 1500000000 !> Fixed capacity of patch_ib (namelist patches + local particle bed subset after reduction) - integer, parameter :: num_ib_patches_max_namelist = 54000 - integer, parameter :: num_local_ibs_max = 2000 !< Maximum number of immersed boundary patches (patch_ib) + integer, parameter :: num_ib_patches_max_namelist = 810000 + integer, parameter :: num_local_ibs_max = 30000 !< Maximum number of immersed boundary patches (patch_ib) integer, parameter :: num_particle_clouds_max = 10 !< Maximum number of particle bed patch specifications integer, parameter :: num_bc_patches_max = 10 !< Maximum number of boundary condition patches integer, parameter :: max_2d_fourier_modes = 10 !< Max Fourier mode index for 2D modal patch (geometry 13) diff --git a/src/simulation/m_collisions.fpp b/src/simulation/m_collisions.fpp index 5498da29c2..2a03e9a616 100644 --- a/src/simulation/m_collisions.fpp +++ b/src/simulation/m_collisions.fpp @@ -65,8 +65,8 @@ contains ! get is distance used in the force calculation with each IB and each wall call s_detect_wall_collisions() - ! call s_detect_ib_collisions(ghost_points, ib_markers, num_gps, num_considered_collisions) - call s_detect_ib_collisions_n2(num_considered_collisions) + call s_detect_ib_collisions(ghost_points, ib_markers, num_gps, num_considered_collisions) + ! call s_detect_ib_collisions_n2(num_considered_collisions) select case (collision_model) case (1) ! soft sphere model diff --git a/src/simulation/m_particle_cloud.fpp b/src/simulation/m_particle_cloud.fpp index 081c5b3463..4fd8970894 100644 --- a/src/simulation/m_particle_cloud.fpp +++ b/src/simulation/m_particle_cloud.fpp @@ -3,6 +3,8 @@ !! @brief Generates particle beds: converts particle_cloud specifications into !! individual sphere/circle particle_cloud_ibs entries before reduction. +#: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. module m_particle_cloud @@ -10,6 +12,7 @@ module m_particle_cloud use m_global_parameters use m_constants use m_mpi_common + use m_collisions implicit none @@ -19,11 +22,13 @@ module m_particle_cloud contains - !> Generate all particle beds and fill particle_cloud_ibs. Called on all ranks before s_reduce_ib_patch_array. + !> Generate all particle beds and fill particle_cloud_ibs. Called on all ranks before s_reduce_ib_patch_array. Particles + !! outside this rank's IB neighborhood (get_neighbor_bounds() must already have been called) are discarded as they are + !! generated, so particle_cloud_ibs never holds more than num_local_ibs_max entries regardless of the global particle count. impure subroutine s_generate_particle_clouds(particle_cloud_ibs) type(ib_patch_parameters), allocatable, intent(out), dimension(:) :: particle_cloud_ibs - integer :: cloud_idx, ib_idx, n_total_particles + integer :: cloud_idx, ib_idx, glbl_idx, n_total_particles real(wp) :: t_start, t_end if (num_particle_clouds == 0) then @@ -33,34 +38,34 @@ contains call cpu_time(t_start) - ! Pre-count total particles across all beds so particle_cloud_ibs can be allocated exactly once. n_total_particles = 0 do cloud_idx = 1, num_particle_clouds n_total_particles = n_total_particles + particle_cloud(cloud_idx)%num_particles end do - allocate (particle_cloud_ibs(n_total_particles)) + allocate (particle_cloud_ibs(min(num_local_ibs_max, n_total_particles))) - ib_idx = 0 ! index into particle_cloud_ibs + ib_idx = 0 ! index into particle_cloud_ibs (this rank's in-neighborhood particles only) + glbl_idx = 0 ! running index across all generated particles, kept regardless of locality do cloud_idx = 1, num_particle_clouds select case (particle_cloud(cloud_idx)%packing_method) case (1) ! random box packing method - call s_particle_cloud_random_box(cloud_idx, ib_idx, particle_cloud_ibs) + call s_particle_cloud_random_box(cloud_idx, ib_idx, glbl_idx, particle_cloud_ibs) case (2) ! lattice packing method - call s_particle_cloud_lattice(cloud_idx, ib_idx, particle_cloud_ibs) + call s_particle_cloud_lattice(cloud_idx, ib_idx, glbl_idx, particle_cloud_ibs) end select end do call cpu_time(t_end) - if (proc_rank == 0) print '(a,i0,a,f0.3,a)', 'Particle beds placed ', ib_idx, ' particles in ', t_end - t_start, ' seconds.' + if (proc_rank == 0) print '(a,i0,a,f0.3,a)', 'Particle beds placed ', glbl_idx, ' particles in ', t_end - t_start, ' seconds.' end subroutine s_generate_particle_clouds !> Generates a random distributions of particles in a box with a minimum spacing - subroutine s_particle_cloud_random_box(cloud_idx, ib_idx, particle_cloud_ibs) + subroutine s_particle_cloud_random_box(cloud_idx, ib_idx, glbl_idx, particle_cloud_ibs) integer, intent(in) :: cloud_idx - integer, intent(inout) :: ib_idx + integer, intent(inout) :: ib_idx, glbl_idx type(ib_patch_parameters), intent(inout), dimension(:) :: particle_cloud_ibs integer :: n_placed, geom, seed integer(8) :: n_attempts, max_attempts @@ -161,7 +166,7 @@ contains chain_next(n_placed) = hash_head(slot) hash_head(slot) = n_placed - call s_add_cloud_particle(cloud_idx, ib_idx, geom, rx, ry, rz, particle_cloud_ibs) + call s_add_cloud_particle(cloud_idx, ib_idx, glbl_idx, geom, rx, ry, rz, particle_cloud_ibs) end if end do @@ -177,10 +182,10 @@ contains !! 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. - subroutine s_particle_cloud_lattice(cloud_idx, ib_idx, particle_cloud_ibs) + subroutine s_particle_cloud_lattice(cloud_idx, ib_idx, glbl_idx, particle_cloud_ibs) integer, intent(in) :: cloud_idx - integer, intent(inout) :: ib_idx + integer, intent(inout) :: ib_idx, glbl_idx type(ib_patch_parameters), intent(inout), dimension(:) :: particle_cloud_ibs integer :: n_placed, n_target, geom integer :: row, col, ncx, ncy, ix, jy, kz, b @@ -225,7 +230,7 @@ contains col = 0 px = x0 do while (px <= xmax .and. n_placed < n_target) - call s_add_cloud_particle(cloud_idx, ib_idx, geom, px, py, particle_cloud(cloud_idx)%z_centroid, & + call s_add_cloud_particle(cloud_idx, ib_idx, glbl_idx, geom, px, py, particle_cloud(cloud_idx)%z_centroid, & & particle_cloud_ibs) n_placed = n_placed + 1 col = col + 1 @@ -247,7 +252,7 @@ contains do ix = 0, ncx - 1 do b = 1, 4 if (n_placed >= n_target) exit - call s_add_cloud_particle(cloud_idx, ib_idx, geom, xmin + real(ix, wp)*cell + bx_off(b), & + call s_add_cloud_particle(cloud_idx, ib_idx, glbl_idx, geom, xmin + real(ix, wp)*cell + bx_off(b), & & ymin + real(jy, wp)*cell + by_off(b), zmin + real(kz, & & wp)*cell + bz_off(b), particle_cloud_ibs) n_placed = n_placed + 1 @@ -260,19 +265,30 @@ contains end subroutine s_particle_cloud_lattice - !> Writes a single placed particle into particle_cloud_ibs at the next free slot, advancing ib_idx. Shared by all packing - !! methods so the per-particle ib_patch_parameters setup stays in one place. - subroutine s_add_cloud_particle(cloud_idx, ib_idx, geom, px, py, pz, particle_cloud_ibs) + !> Writes a single placed particle into particle_cloud_ibs at the next free slot if it falls within this rank's IB + !! neighborhood (get_neighbor_bounds() must already have run), advancing ib_idx; particles outside the neighborhood are + !! discarded here rather than stored. glbl_idx always advances so gbl_patch_id records the particle's position in the full, + !! unfiltered enumeration - s_reduce_ib_patch_array offsets it by the namelist IB count for final global indexing. Shared by + !! all packing methods so the per-particle ib_patch_parameters setup stays in one place. + subroutine s_add_cloud_particle(cloud_idx, ib_idx, glbl_idx, geom, px, py, pz, particle_cloud_ibs) integer, intent(in) :: cloud_idx, geom - integer, intent(inout) :: ib_idx + integer, intent(inout) :: ib_idx, glbl_idx real(wp), intent(in) :: px, py, pz type(ib_patch_parameters), intent(inout), dimension(:) :: particle_cloud_ibs + real(wp), dimension(3) :: centroid + + glbl_idx = glbl_idx + 1 + + centroid = [px, py, 0._wp] + if (num_dims == 3) centroid(3) = pz + if (.not. f_neighborhood_ranks_own_location(centroid)) return ib_idx = ib_idx + 1 + @:PROHIBIT(ib_idx > num_local_ibs_max, & + & "Too many particle-cloud IBs in one rank's neighborhood. Modify case file or increase num_local_ibs_max.") - ! gbl_patch_id is relative within particle_cloud_ibs here; s_reduce_ib_patch_array adjusts to global indexing. - particle_cloud_ibs(ib_idx)%gbl_patch_id = ib_idx + 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 diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index 488c858535..7fe7c1d637 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -876,6 +876,10 @@ contains block type(ib_patch_parameters), allocatable :: particle_cloud_ibs(:) + ! Neighborhood bounds must exist before s_generate_particle_clouds so it can discard out-of-neighborhood + ! particles as they are generated, instead of materializing the full global particle-cloud on every rank. + call get_neighbor_bounds() + if (cfl_dt .and. n_start > 0) then call s_read_ib_restart_data(n_start) allocate (particle_cloud_ibs(0)) @@ -1200,23 +1204,28 @@ contains end subroutine s_read_ib_restart_data - !> @brief Merges patch_ib (namelist patches, fixed at num_ib_patches_max_namelist) with particle_cloud_ibs (CPU-only, exact - !! size) 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. + !> @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, using the gbl_patch_id each entry was tagged with at + !! generation time to recover its position in the full, unfiltered particle ordering) 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) type(ib_patch_parameters), intent(in), dimension(:) :: particle_cloud_ibs real(wp), dimension(3) :: centroid integer :: i - integer :: num_namelist_ibs, num_bed_ibs + integer :: num_namelist_ibs, num_bed_ibs, num_local_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 + num_local_bed_ibs = size(particle_cloud_ibs) - ! Check for moving IBs across both namelist and particle bed patches. + ! Check for moving IBs across both namelist and particle bed patches. moving_ibm is a per-cloud property (every particle + ! in a cloud shares particle_cloud(:)%moving_ibm), so checking the cloud list gives the same, globally-consistent answer + ! on every rank without needing this rank's full (neighborhood-filtered) particle_cloud_ibs. moving_immersed_boundary_flag = .false. do i = 1, num_namelist_ibs if (patch_ib(i)%moving_ibm /= 0) then @@ -1225,15 +1234,14 @@ contains end if end do if (.not. moving_immersed_boundary_flag) then - do i = 1, num_bed_ibs - if (particle_cloud_ibs(i)%moving_ibm /= 0) 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 - call get_neighbor_bounds() call s_compute_ib_neighbor_ranks() num_gbl_ibs = num_namelist_ibs + num_bed_ibs @@ -1243,9 +1251,9 @@ contains ! single-rank: all patches are local; append particle bed entries directly into patch_ib. @:PROHIBIT(num_gbl_ibs > num_ib_patches_max_namelist, & & "Total IB count exceeds patch_ib capacity. Increase num_ib_patches_max_namelist.") - do i = 1, num_bed_ibs + do i = 1, num_local_bed_ibs patch_ib(num_namelist_ibs + i) = particle_cloud_ibs(i) - patch_ib(num_namelist_ibs + i)%gbl_patch_id = num_namelist_ibs + i + patch_ib(num_namelist_ibs + i)%gbl_patch_id = num_namelist_ibs + particle_cloud_ibs(i)%gbl_patch_id end do num_ibs = num_gbl_ibs num_local_ibs = num_gbl_ibs @@ -1269,19 +1277,20 @@ contains end if end if end do - do i = 1, num_bed_ibs + ! particle_cloud_ibs entries already passed the neighborhood check at generation time (against the same + ! neighbor_domain_x/y/z computed by get_neighbor_bounds() before s_generate_particle_clouds ran), so no need to + ! recheck it here. + do i = 1, num_local_bed_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) - patch_ib(num_ibs)%gbl_patch_id = num_namelist_ibs + i - if (f_local_rank_owns_location(centroid)) then - num_local_ibs = num_local_ibs + 1 - local_ib_patch_ids(num_local_ibs) = num_ibs - end if + 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) + patch_ib(num_ibs)%gbl_patch_id = num_namelist_ibs + particle_cloud_ibs(i)%gbl_patch_id + if (f_local_rank_owns_location(centroid)) then + num_local_ibs = num_local_ibs + 1 + local_ib_patch_ids(num_local_ibs) = num_ibs end if end do @:PROHIBIT(num_local_ibs > num_local_ibs_max, & @@ -1291,9 +1300,9 @@ contains ! no-MPI: all patches are local; append particle bed entries directly into patch_ib. @:PROHIBIT(num_gbl_ibs > num_ib_patches_max_namelist, & & "Total IB count exceeds patch_ib capacity. Increase num_ib_patches_max_namelist.") - do i = 1, num_bed_ibs + do i = 1, num_local_bed_ibs patch_ib(num_namelist_ibs + i) = particle_cloud_ibs(i) - patch_ib(num_namelist_ibs + i)%gbl_patch_id = num_namelist_ibs + i + patch_ib(num_namelist_ibs + i)%gbl_patch_id = num_namelist_ibs + particle_cloud_ibs(i)%gbl_patch_id end do num_ibs = num_gbl_ibs num_local_ibs = num_gbl_ibs From a22bea913d378ffe1e4f74ed3340cc3710943ae6 Mon Sep 17 00:00:00 2001 From: Daniel Vickers Date: Wed, 8 Jul 2026 23:14:54 -0400 Subject: [PATCH 2/8] Fixed an out-of-bounds error caused by reading a bad length for the number of IBs. --- src/common/m_constants.fpp | 2 +- src/simulation/m_particle_cloud.fpp | 39 ++++++++++++++++++----------- src/simulation/m_start_up.fpp | 37 ++++++++++++++------------- 3 files changed, 45 insertions(+), 33 deletions(-) diff --git a/src/common/m_constants.fpp b/src/common/m_constants.fpp index ed67683917..662a8d62fc 100644 --- a/src/common/m_constants.fpp +++ b/src/common/m_constants.fpp @@ -30,7 +30,7 @@ module m_constants integer, parameter :: num_ib_patches_max = 1500000000 !> Fixed capacity of patch_ib (namelist patches + local particle bed subset after reduction) integer, parameter :: num_ib_patches_max_namelist = 810000 - integer, parameter :: num_local_ibs_max = 30000 !< Maximum number of immersed boundary patches (patch_ib) + integer, parameter :: num_local_ibs_max = 30000 !< Maximum number of immersed boundary patches (patch_ib) integer, parameter :: num_particle_clouds_max = 10 !< Maximum number of particle bed patch specifications integer, parameter :: num_bc_patches_max = 10 !< Maximum number of boundary condition patches integer, parameter :: max_2d_fourier_modes = 10 !< Max Fourier mode index for 2D modal patch (geometry 13) diff --git a/src/simulation/m_particle_cloud.fpp b/src/simulation/m_particle_cloud.fpp index 4fd8970894..6a23bb56ad 100644 --- a/src/simulation/m_particle_cloud.fpp +++ b/src/simulation/m_particle_cloud.fpp @@ -22,17 +22,23 @@ module m_particle_cloud contains - !> Generate all particle beds and fill particle_cloud_ibs. Called on all ranks before s_reduce_ib_patch_array. Particles - !! outside this rank's IB neighborhood (get_neighbor_bounds() must already have been called) are discarded as they are - !! generated, so particle_cloud_ibs never holds more than num_local_ibs_max entries regardless of the global particle count. - impure subroutine s_generate_particle_clouds(particle_cloud_ibs) + !> Generate all particle beds and fill particle_cloud_ibs. Called on all ranks before s_reduce_ib_patch_array. Particles outside + !! this rank's IB neighborhood (get_neighbor_bounds() must already have been called) are discarded as they are generated, so + !! particle_cloud_ibs is allocated to a worst-case capacity of num_ib_patches_max_namelist entries (the same neighborhood-sized + !! cap patch_ib is bound by in s_reduce_ib_patch_array) regardless of the global particle count, but 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. Each entry's gbl_patch_id is already the final, absolute global patch id + !! (namelist patches occupy 1..num_ibs, so the running index here starts at num_ibs) - s_reduce_ib_patch_array copies it as-is. + impure subroutine s_generate_particle_clouds(particle_cloud_ibs, num_particle_cloud_ibs) type(ib_patch_parameters), allocatable, intent(out), dimension(:) :: particle_cloud_ibs + integer, intent(out) :: num_particle_cloud_ibs integer :: cloud_idx, ib_idx, glbl_idx, n_total_particles real(wp) :: t_start, t_end if (num_particle_clouds == 0) then allocate (particle_cloud_ibs(0)) + num_particle_cloud_ibs = 0 return end if @@ -42,10 +48,11 @@ contains do cloud_idx = 1, num_particle_clouds n_total_particles = n_total_particles + particle_cloud(cloud_idx)%num_particles end do - allocate (particle_cloud_ibs(min(num_local_ibs_max, n_total_particles))) + allocate (particle_cloud_ibs(min(num_ib_patches_max_namelist, n_total_particles))) - ib_idx = 0 ! index into particle_cloud_ibs (this rank's in-neighborhood particles only) - glbl_idx = 0 ! running index across all generated particles, kept regardless of locality + ib_idx = 0 ! index into particle_cloud_ibs (this rank's in-neighborhood particles only) + glbl_idx = num_ibs ! running global patch id across all generated particles, kept regardless of locality; starts after + ! the namelist patches (1..num_ibs) so each particle's gbl_patch_id is already its final global id do cloud_idx = 1, num_particle_clouds select case (particle_cloud(cloud_idx)%packing_method) @@ -55,9 +62,11 @@ contains call s_particle_cloud_lattice(cloud_idx, ib_idx, glbl_idx, particle_cloud_ibs) end select end do + num_particle_cloud_ibs = ib_idx call cpu_time(t_end) - if (proc_rank == 0) print '(a,i0,a,f0.3,a)', 'Particle beds placed ', glbl_idx, ' particles in ', t_end - t_start, ' seconds.' + if (proc_rank == 0) print '(a,i0,a,f0.3,a)', 'Particle beds placed ', glbl_idx - num_ibs, ' particles in ', & + & t_end - t_start, ' seconds.' end subroutine s_generate_particle_clouds @@ -265,11 +274,11 @@ contains end subroutine s_particle_cloud_lattice - !> Writes a single placed particle into particle_cloud_ibs at the next free slot if it falls within this rank's IB - !! neighborhood (get_neighbor_bounds() must already have run), advancing ib_idx; particles outside the neighborhood are - !! discarded here rather than stored. glbl_idx always advances so gbl_patch_id records the particle's position in the full, - !! unfiltered enumeration - s_reduce_ib_patch_array offsets it by the namelist IB count for final global indexing. 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 if it falls within this rank's IB neighborhood + !! (get_neighbor_bounds() must already have run), advancing ib_idx; particles outside the neighborhood are discarded here rather + !! than stored. glbl_idx always advances (starting from num_ibs, see s_generate_particle_clouds) so gbl_patch_id is already the + !! particle's final, absolute global patch id - s_reduce_ib_patch_array copies it as-is. Shared by all packing methods so the + !! per-particle ib_patch_parameters setup stays in one place. subroutine s_add_cloud_particle(cloud_idx, ib_idx, glbl_idx, geom, px, py, pz, particle_cloud_ibs) integer, intent(in) :: cloud_idx, geom @@ -285,8 +294,8 @@ contains if (.not. f_neighborhood_ranks_own_location(centroid)) return ib_idx = ib_idx + 1 - @:PROHIBIT(ib_idx > num_local_ibs_max, & - & "Too many particle-cloud IBs in one rank's neighborhood. Modify case file or increase num_local_ibs_max.") + @:PROHIBIT(ib_idx > num_ib_patches_max_namelist, & + & "Too many particle-cloud IBs in one rank's neighborhood. 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/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index 7fe7c1d637..0853605e5f 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -875,6 +875,7 @@ contains if (ib) then block type(ib_patch_parameters), allocatable :: particle_cloud_ibs(:) + integer :: num_particle_cloud_ibs ! Neighborhood bounds must exist before s_generate_particle_clouds so it can discard out-of-neighborhood ! particles as they are generated, instead of materializing the full global particle-cloud on every rank. @@ -883,15 +884,17 @@ contains 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 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) + call s_generate_particle_clouds(particle_cloud_ibs, num_particle_cloud_ibs) end if call s_instantiate_STL_models() call s_initialize_ib_airfoils() - call s_reduce_ib_patch_array(particle_cloud_ibs) + call s_reduce_ib_patch_array(particle_cloud_ibs, num_particle_cloud_ibs) deallocate (particle_cloud_ibs) end block call s_ibm_setup() @@ -1205,23 +1208,24 @@ contains 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, using the gbl_patch_id each entry was tagged with at - !! generation time to recover its position in the full, unfiltered particle ordering) 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) + !! 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. + 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_local_bed_ibs + 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 - num_local_bed_ibs = size(particle_cloud_ibs) ! Check for moving IBs across both namelist and particle bed patches. moving_ibm is a per-cloud property (every particle ! in a cloud shares particle_cloud(:)%moving_ibm), so checking the cloud list gives the same, globally-consistent answer @@ -1251,9 +1255,8 @@ contains ! single-rank: all patches are local; append particle bed entries directly into patch_ib. @:PROHIBIT(num_gbl_ibs > num_ib_patches_max_namelist, & & "Total IB count exceeds patch_ib capacity. Increase num_ib_patches_max_namelist.") - do i = 1, num_local_bed_ibs + do i = 1, num_particle_cloud_ibs patch_ib(num_namelist_ibs + i) = particle_cloud_ibs(i) - patch_ib(num_namelist_ibs + i)%gbl_patch_id = num_namelist_ibs + particle_cloud_ibs(i)%gbl_patch_id end do num_ibs = num_gbl_ibs num_local_ibs = num_gbl_ibs @@ -1273,6 +1276,8 @@ contains patch_ib(num_ibs)%gbl_patch_id = 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 @@ -1280,29 +1285,27 @@ contains ! particle_cloud_ibs entries already passed the neighborhood check at generation time (against the same ! neighbor_domain_x/y/z computed by get_neighbor_bounds() before s_generate_particle_clouds ran), so no need to ! recheck it here. - do i = 1, num_local_bed_ibs + 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) - patch_ib(num_ibs)%gbl_patch_id = num_namelist_ibs + particle_cloud_ibs(i)%gbl_patch_id 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 do - @: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.") end if #else ! no-MPI: all patches are local; append particle bed entries directly into patch_ib. @:PROHIBIT(num_gbl_ibs > num_ib_patches_max_namelist, & & "Total IB count exceeds patch_ib capacity. Increase num_ib_patches_max_namelist.") - do i = 1, num_local_bed_ibs + do i = 1, num_particle_cloud_ibs patch_ib(num_namelist_ibs + i) = particle_cloud_ibs(i) - patch_ib(num_namelist_ibs + i)%gbl_patch_id = num_namelist_ibs + particle_cloud_ibs(i)%gbl_patch_id end do num_ibs = num_gbl_ibs num_local_ibs = num_gbl_ibs From a1bec2f9414781813ac19bc9732a3fdb1f4a08d3 Mon Sep 17 00:00:00 2001 From: Daniel Vickers Date: Wed, 8 Jul 2026 23:28:29 -0400 Subject: [PATCH 3/8] Cleaned up comments --- src/simulation/m_particle_cloud.fpp | 9 ++++----- src/simulation/m_start_up.fpp | 10 ++-------- 2 files changed, 6 insertions(+), 13 deletions(-) diff --git a/src/simulation/m_particle_cloud.fpp b/src/simulation/m_particle_cloud.fpp index 6a23bb56ad..698c9c62d8 100644 --- a/src/simulation/m_particle_cloud.fpp +++ b/src/simulation/m_particle_cloud.fpp @@ -50,9 +50,8 @@ contains end do allocate (particle_cloud_ibs(min(num_ib_patches_max_namelist, n_total_particles))) - ib_idx = 0 ! index into particle_cloud_ibs (this rank's in-neighborhood particles only) - glbl_idx = num_ibs ! running global patch id across all generated particles, kept regardless of locality; starts after - ! the namelist patches (1..num_ibs) so each particle's gbl_patch_id is already its final global id + ib_idx = 0 ! index into particle_cloud_ibs + glbl_idx = num_ibs do cloud_idx = 1, num_particle_clouds select case (particle_cloud(cloud_idx)%packing_method) @@ -213,7 +212,7 @@ contains n_target = particle_cloud(cloud_idx)%num_particles n_placed = 0 - if (p == 0) then + 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))) @@ -228,7 +227,7 @@ contains & // "reduce num_particles or min_spacing, or enlarge the cloud region") end if - if (p == 0) then + if (num_dims < 3) then ! Triangular lattice: rows pitched by spacing*sqrt(3)/2, odd rows shifted by half a spacing. row_dy = spacing*sqrt(3._wp)/2._wp row = 0 diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index 0853605e5f..9b4aab8127 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -877,8 +877,6 @@ contains type(ib_patch_parameters), allocatable :: particle_cloud_ibs(:) integer :: num_particle_cloud_ibs - ! Neighborhood bounds must exist before s_generate_particle_clouds so it can discard out-of-neighborhood - ! particles as they are generated, instead of materializing the full global particle-cloud on every rank. call get_neighbor_bounds() if (cfl_dt .and. n_start > 0) then @@ -1227,9 +1225,7 @@ contains num_bed_ibs = num_bed_ibs + particle_cloud(i)%num_particles end do - ! Check for moving IBs across both namelist and particle bed patches. moving_ibm is a per-cloud property (every particle - ! in a cloud shares particle_cloud(:)%moving_ibm), so checking the cloud list gives the same, globally-consistent answer - ! on every rank without needing this rank's full (neighborhood-filtered) particle_cloud_ibs. + ! 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 @@ -1282,9 +1278,7 @@ contains end if end if end do - ! particle_cloud_ibs entries already passed the neighborhood check at generation time (against the same - ! neighbor_domain_x/y/z computed by get_neighbor_bounds() before s_generate_particle_clouds ran), so no need to - ! recheck it here. + ! particle_cloud_ibs entries already passed the neighborhood check at generation time so no need to recheck it here. 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 From ba5a7a5247b5211f02dbd25899105eadc5b05047 Mon Sep 17 00:00:00 2001 From: Daniel Vickers Date: Mon, 27 Jul 2026 10:46:37 -0400 Subject: [PATCH 4/8] file access stagger and cleanup from timings in full-scale run --- src/common/m_boundary_io.fpp | 4 +-- src/common/m_delay_file_access.f90 | 6 ++--- src/pre_process/m_data_output.fpp | 2 +- src/simulation/m_data_output.fpp | 8 ++++-- src/simulation/m_ibm.fpp | 1 + src/simulation/m_start_up.fpp | 40 +++++++++++++++++++++++++++--- 6 files changed, 50 insertions(+), 11 deletions(-) diff --git a/src/common/m_boundary_io.fpp b/src/common/m_boundary_io.fpp index e926c39570..cf17f74b22 100644 --- a/src/common/m_boundary_io.fpp +++ b/src/common/m_boundary_io.fpp @@ -130,7 +130,7 @@ contains call s_mpi_barrier() - call DelayFileAccess(proc_rank) + call s_delay_file_access(proc_rank) write (proc_rank_str, '(I7.7)') proc_rank file_path = trim(file_loc) // '/bc_' // trim(proc_rank_str) // '.dat' @@ -233,7 +233,7 @@ contains call s_mpi_barrier() - call DelayFileAccess(proc_rank) + call s_delay_file_access(proc_rank) write (proc_rank_str, '(I7.7)') proc_rank file_path = trim(file_loc) // '/bc_' // trim(proc_rank_str) // '.dat' diff --git a/src/common/m_delay_file_access.f90 b/src/common/m_delay_file_access.f90 index 0f154af46c..a5d110bb08 100644 --- a/src/common/m_delay_file_access.f90 +++ b/src/common/m_delay_file_access.f90 @@ -11,14 +11,14 @@ module m_delay_file_access private - public :: DelayFileAccess + public :: s_delay_file_access integer, private, parameter :: N_PROCESSES_FILE_ACCESS = 128, FILE_ACCESS_DELAY_UNIT = 10000 contains !> Introduce a rank-dependent busy-wait delay to stagger parallel file access and reduce I/O contention. - impure subroutine DelayFileAccess(ProcessRank) + impure subroutine s_delay_file_access(ProcessRank) integer, intent(in) :: ProcessRank integer :: iDelay, nFileAccessDelayIterations @@ -31,6 +31,6 @@ impure subroutine DelayFileAccess(ProcessRank) Dummy = Number*Number end do - end subroutine DelayFileAccess + end subroutine s_delay_file_access end module m_delay_file_access diff --git a/src/pre_process/m_data_output.fpp b/src/pre_process/m_data_output.fpp index 54c6e84b1a..3f148cdfdd 100644 --- a/src/pre_process/m_data_output.fpp +++ b/src/pre_process/m_data_output.fpp @@ -454,7 +454,7 @@ contains call s_create_directory(trim(file_loc)) end if call s_mpi_barrier() - call DelayFileAccess(proc_rank) + call s_delay_file_access(proc_rank) if (down_sample) then call s_initialize_mpi_data_ds(q_cons_temp) diff --git a/src/simulation/m_data_output.fpp b/src/simulation/m_data_output.fpp index f6b58f4e39..a492c4c638 100644 --- a/src/simulation/m_data_output.fpp +++ b/src/simulation/m_data_output.fpp @@ -712,7 +712,7 @@ contains call s_create_directory(trim(file_loc)) end if call s_mpi_barrier() - call DelayFileAccess(proc_rank) + call s_delay_file_access(proc_rank) call s_initialize_mpi_data(q_cons_vf) @@ -896,6 +896,10 @@ contains write (file_loc, '(A)') 'ib.dat' file_loc = trim(case_dir) // '/restart_data' // trim(mpiiofs) // trim(file_loc) + + call s_mpi_barrier() + call s_delay_file_access(proc_rank) + call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, ior(MPI_MODE_WRONLY, MPI_MODE_CREATE), mpi_info_int, ifile, ierr) var_MOK = int(sys_size + 1, MPI_OFFSET_KIND) @@ -955,7 +959,7 @@ contains call s_create_directory(trim(file_loc)) end if call s_mpi_barrier() - call DelayFileAccess(proc_rank) + call s_delay_file_access(proc_rank) 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) diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp index d2e598ed22..25ebddf718 100644 --- a/src/simulation/m_ibm.fpp +++ b/src/simulation/m_ibm.fpp @@ -983,6 +983,7 @@ contains real(wp), dimension(1:3,1:3) :: viscous_stress real(wp), dimension(1:3) :: local_force_contribution, radial_vector, local_torque_contribution real(wp) :: cell_volume, dynamic_viscosity + real(wp) :: t_start, t_end #:if not MFC_CASE_OPTIMIZATION and USING_AMD real(wp), dimension(3) :: dynamic_viscosities diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index 9b4aab8127..e30bf3a0ed 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -796,8 +796,9 @@ contains !> Initialize all simulation sub-modules in the required dependency order impure subroutine s_initialize_modules - integer :: m_ds, n_ds, p_ds - integer :: i + integer :: m_ds, n_ds, p_ds + integer :: i + real(wp) :: t_start, t_end call s_initialize_global_parameters_module() #:if USING_AMD @@ -813,9 +814,14 @@ contains if (bubbles_euler .or. bubbles_lagrange) then call s_initialize_bubbles_model() end if + if (proc_rank == 0) print *, "Starting s_initialize_mpi_common_module" + call cpu_time(t_start) call s_initialize_mpi_common_module() call s_initialize_mpi_proxy_module() call s_initialize_variables_conversion_module() + call cpu_time(t_end) + if (proc_rank == 0) print *, "s_initialize_mpi_common_module finished in ", t_end - t_start + if (grid_geometry == 3) call s_initialize_fftw_module() if (bubbles_euler) call s_initialize_bubbles_EE_module() @@ -838,9 +844,13 @@ contains if (relax) call s_initialize_phasechange_module() + if (proc_rank == 0) print *, "Starting s_initialize_data_output_module" + call cpu_time(t_start) call s_initialize_data_output_module() call s_initialize_derived_variables_module() call s_initialize_time_steppers_module() + call cpu_time(t_end) + if (proc_rank == 0) print *, "s_initialize_data_output_module finished in ", t_end - t_start call s_initialize_boundary_common_module() @@ -866,7 +876,11 @@ contains end do deallocate (q_cons_temp) else + if (proc_rank == 0) print *, "Starting s_read_data_files" + call cpu_time(t_start) call s_read_data_files(q_cons_ts(1)%vf) + call cpu_time(t_end) + if (proc_rank == 0) print *, "s_read_data_files finished in ", t_end - t_start end if call s_populate_grid_variables_buffers() @@ -877,7 +891,11 @@ contains type(ib_patch_parameters), allocatable :: particle_cloud_ibs(:) integer :: num_particle_cloud_ibs + if (proc_rank == 0) print *, "Starting get_neighbor_bounds" + call cpu_time(t_start) call get_neighbor_bounds() + call cpu_time(t_end) + if (proc_rank == 0) print *, "get_neighbor_bounds finished in ", t_end - t_start if (cfl_dt .and. n_start > 0) then call s_read_ib_restart_data(n_start) @@ -888,14 +906,23 @@ contains allocate (particle_cloud_ibs(0)) num_particle_cloud_ibs = 0 else + if (proc_rank == 0) print *, "Starting s_generate_particle_clouds" call s_generate_particle_clouds(particle_cloud_ibs, num_particle_cloud_ibs) end if call s_instantiate_STL_models() call s_initialize_ib_airfoils() + if (proc_rank == 0) print *, "Starting s_reduce_ib_patch_array" + call cpu_time(t_start) call s_reduce_ib_patch_array(particle_cloud_ibs, num_particle_cloud_ibs) + call cpu_time(t_end) + if (proc_rank == 0) print *, "s_reduce_ib_patch_array finished in ", t_end - t_start deallocate (particle_cloud_ibs) end block + if (proc_rank == 0) print *, "Starting s_ibm_setup" + call cpu_time(t_start) call s_ibm_setup() + call cpu_time(t_end) + if (proc_rank == 0) print *, "s_ibm_setup finished in ", t_end - t_start 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) @@ -927,12 +954,15 @@ contains if (hypoelasticity) call s_initialize_hypoelastic_module() + if (proc_rank == 0) print *, "exiting s_initialize_modules" + end subroutine s_initialize_modules !> Set up the MPI execution environment, bind GPUs, and decompose the computational domain impure subroutine s_initialize_mpi_domain - integer :: ierr + integer :: ierr + real(wp) :: t_start, t_end #ifdef MFC_GPU real(wp) :: starttime, endtime @@ -991,7 +1021,11 @@ contains #endif end if + if (proc_rank == 0) print *, "Starting s_mpi_bcast_user_inputs" + call cpu_time(t_start) call s_mpi_bcast_user_inputs() + call cpu_time(t_end) + if (proc_rank == 0) print *, "s_mpi_bcast_user_inputs finished in ", t_end - t_start ! Save original BCs before decomposition overwrites them with MPI neighbor ranks ib_bc_x = bc_x From efec23a27ad4ecdb41a27ceddf9a03004d13ea22 Mon Sep 17 00:00:00 2001 From: Daniel Vickers Date: Mon, 27 Jul 2026 10:51:46 -0400 Subject: [PATCH 5/8] Missed a few prints --- src/simulation/m_ibm.fpp | 1 - src/simulation/m_start_up.fpp | 39 +++-------------------------------- 2 files changed, 3 insertions(+), 37 deletions(-) diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp index 25ebddf718..d2e598ed22 100644 --- a/src/simulation/m_ibm.fpp +++ b/src/simulation/m_ibm.fpp @@ -983,7 +983,6 @@ contains real(wp), dimension(1:3,1:3) :: viscous_stress real(wp), dimension(1:3) :: local_force_contribution, radial_vector, local_torque_contribution real(wp) :: cell_volume, dynamic_viscosity - real(wp) :: t_start, t_end #:if not MFC_CASE_OPTIMIZATION and USING_AMD real(wp), dimension(3) :: dynamic_viscosities diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index e30bf3a0ed..9e14f8a0e4 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -796,9 +796,8 @@ contains !> Initialize all simulation sub-modules in the required dependency order impure subroutine s_initialize_modules - integer :: m_ds, n_ds, p_ds - integer :: i - real(wp) :: t_start, t_end + integer :: m_ds, n_ds, p_ds + integer :: i call s_initialize_global_parameters_module() #:if USING_AMD @@ -814,13 +813,9 @@ contains if (bubbles_euler .or. bubbles_lagrange) then call s_initialize_bubbles_model() end if - if (proc_rank == 0) print *, "Starting s_initialize_mpi_common_module" - call cpu_time(t_start) call s_initialize_mpi_common_module() call s_initialize_mpi_proxy_module() call s_initialize_variables_conversion_module() - call cpu_time(t_end) - if (proc_rank == 0) print *, "s_initialize_mpi_common_module finished in ", t_end - t_start if (grid_geometry == 3) call s_initialize_fftw_module() @@ -844,13 +839,9 @@ contains if (relax) call s_initialize_phasechange_module() - if (proc_rank == 0) print *, "Starting s_initialize_data_output_module" - call cpu_time(t_start) call s_initialize_data_output_module() call s_initialize_derived_variables_module() call s_initialize_time_steppers_module() - call cpu_time(t_end) - if (proc_rank == 0) print *, "s_initialize_data_output_module finished in ", t_end - t_start call s_initialize_boundary_common_module() @@ -876,11 +867,7 @@ contains end do deallocate (q_cons_temp) else - if (proc_rank == 0) print *, "Starting s_read_data_files" - call cpu_time(t_start) call s_read_data_files(q_cons_ts(1)%vf) - call cpu_time(t_end) - if (proc_rank == 0) print *, "s_read_data_files finished in ", t_end - t_start end if call s_populate_grid_variables_buffers() @@ -891,11 +878,7 @@ contains type(ib_patch_parameters), allocatable :: particle_cloud_ibs(:) integer :: num_particle_cloud_ibs - if (proc_rank == 0) print *, "Starting get_neighbor_bounds" - call cpu_time(t_start) call get_neighbor_bounds() - call cpu_time(t_end) - if (proc_rank == 0) print *, "get_neighbor_bounds finished in ", t_end - t_start if (cfl_dt .and. n_start > 0) then call s_read_ib_restart_data(n_start) @@ -906,23 +889,14 @@ contains allocate (particle_cloud_ibs(0)) num_particle_cloud_ibs = 0 else - if (proc_rank == 0) print *, "Starting s_generate_particle_clouds" call s_generate_particle_clouds(particle_cloud_ibs, num_particle_cloud_ibs) end if call s_instantiate_STL_models() call s_initialize_ib_airfoils() - if (proc_rank == 0) print *, "Starting s_reduce_ib_patch_array" - call cpu_time(t_start) call s_reduce_ib_patch_array(particle_cloud_ibs, num_particle_cloud_ibs) - call cpu_time(t_end) - if (proc_rank == 0) print *, "s_reduce_ib_patch_array finished in ", t_end - t_start deallocate (particle_cloud_ibs) end block - if (proc_rank == 0) print *, "Starting s_ibm_setup" - call cpu_time(t_start) call s_ibm_setup() - call cpu_time(t_end) - if (proc_rank == 0) print *, "s_ibm_setup finished in ", t_end - t_start 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) @@ -954,15 +928,12 @@ contains if (hypoelasticity) call s_initialize_hypoelastic_module() - if (proc_rank == 0) print *, "exiting s_initialize_modules" - end subroutine s_initialize_modules !> Set up the MPI execution environment, bind GPUs, and decompose the computational domain impure subroutine s_initialize_mpi_domain - integer :: ierr - real(wp) :: t_start, t_end + integer :: ierr #ifdef MFC_GPU real(wp) :: starttime, endtime @@ -1021,11 +992,7 @@ contains #endif end if - if (proc_rank == 0) print *, "Starting s_mpi_bcast_user_inputs" - call cpu_time(t_start) call s_mpi_bcast_user_inputs() - call cpu_time(t_end) - if (proc_rank == 0) print *, "s_mpi_bcast_user_inputs finished in ", t_end - t_start ! Save original BCs before decomposition overwrites them with MPI neighbor ranks ib_bc_x = bc_x From 3a0d013c150245014719a40e91fc751e142875d2 Mon Sep 17 00:00:00 2001 From: Daniel Vickers Date: Wed, 29 Jul 2026 22:55:14 -0400 Subject: [PATCH 6/8] Fixing issue with rejection packing algorithm not having enough allocated memory to pair down collision list --- src/simulation/m_particle_cloud.fpp | 175 ++++++++++++++++++---------- 1 file changed, 115 insertions(+), 60 deletions(-) diff --git a/src/simulation/m_particle_cloud.fpp b/src/simulation/m_particle_cloud.fpp index 698c9c62d8..b7eb3b9515 100644 --- a/src/simulation/m_particle_cloud.fpp +++ b/src/simulation/m_particle_cloud.fpp @@ -22,18 +22,20 @@ module m_particle_cloud contains - !> Generate all particle beds and fill particle_cloud_ibs. Called on all ranks before s_reduce_ib_patch_array. Particles outside - !! this rank's IB neighborhood (get_neighbor_bounds() must already have been called) are discarded as they are generated, so - !! particle_cloud_ibs is allocated to a worst-case capacity of num_ib_patches_max_namelist entries (the same neighborhood-sized - !! cap patch_ib is bound by in s_reduce_ib_patch_array) regardless of the global particle count, but 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. Each entry's gbl_patch_id is already the final, absolute global patch id - !! (namelist patches occupy 1..num_ibs, so the running index here starts at num_ibs) - s_reduce_ib_patch_array copies it as-is. + !> 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_random_box) and + !! hands back only the entries that fall within this rank's IB neighborhood (get_neighbor_bounds() must already have been + !! called); those are copied here into particle_cloud_ibs, allocated once to a worst-case capacity of + !! min(num_ib_patches_max_namelist, total particles requested across all clouds). 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. Each entry's gbl_patch_id is already the final, absolute global patch id (namelist patches occupy 1..num_ibs, + !! so the running index here starts at num_ibs) - s_reduce_ib_patch_array copies it as-is. impure subroutine s_generate_particle_clouds(particle_cloud_ibs, num_particle_cloud_ibs) type(ib_patch_parameters), allocatable, intent(out), dimension(:) :: particle_cloud_ibs integer, intent(out) :: num_particle_cloud_ibs - integer :: cloud_idx, ib_idx, glbl_idx, n_total_particles + type(ib_patch_parameters), allocatable :: cloud_ibs(:) + integer :: cloud_idx, glbl_idx, num_cloud_ibs, n_total_particles real(wp) :: t_start, t_end if (num_particle_clouds == 0) then @@ -50,18 +52,23 @@ contains end do allocate (particle_cloud_ibs(min(num_ib_patches_max_namelist, n_total_particles))) - ib_idx = 0 ! index into particle_cloud_ibs + num_particle_cloud_ibs = 0 glbl_idx = num_ibs do cloud_idx = 1, num_particle_clouds select case (particle_cloud(cloud_idx)%packing_method) case (1) ! random box packing method - call s_particle_cloud_random_box(cloud_idx, ib_idx, glbl_idx, particle_cloud_ibs) + call s_particle_cloud_random_box(cloud_idx, glbl_idx, cloud_ibs, num_cloud_ibs) case (2) ! lattice packing method - call s_particle_cloud_lattice(cloud_idx, ib_idx, glbl_idx, particle_cloud_ibs) + call s_particle_cloud_lattice(cloud_idx, glbl_idx, cloud_ibs, num_cloud_ibs) end select + + @:PROHIBIT(num_particle_cloud_ibs + num_cloud_ibs > num_ib_patches_max_namelist, & + & "Too many particle-cloud IBs in one rank's neighborhood. Modify case file or increase num_ib_patches_max_namelist.") + particle_cloud_ibs(num_particle_cloud_ibs + 1:num_particle_cloud_ibs + num_cloud_ibs) = cloud_ibs(1:num_cloud_ibs) + num_particle_cloud_ibs = num_particle_cloud_ibs + num_cloud_ibs + deallocate (cloud_ibs) end do - num_particle_cloud_ibs = ib_idx call cpu_time(t_end) if (proc_rank == 0) print '(a,i0,a,f0.3,a)', 'Particle beds placed ', glbl_idx - num_ibs, ' particles in ', & @@ -69,22 +76,33 @@ contains end subroutine s_generate_particle_clouds - !> Generates a random distributions of particles in a box with a minimum spacing - subroutine s_particle_cloud_random_box(cloud_idx, ib_idx, glbl_idx, particle_cloud_ibs) - - integer, intent(in) :: cloud_idx - integer, intent(inout) :: ib_idx, glbl_idx - type(ib_patch_parameters), intent(inout), dimension(:) :: particle_cloud_ibs - integer :: n_placed, geom, seed - integer(8) :: n_attempts, max_attempts - real(wp) :: xmin, xmax, ymin, ymax, zmin, zmax, min_dist - real(wp) :: rx, ry, rz, dist - logical :: overlaps - real(wp), allocatable :: placed(:,:) - integer :: hash_size, slot - integer :: bx, by, bz, nbx, nby, nbz - integer :: dx_b, dy_b, dz_b, dz_lo, dz_hi, j - integer, allocatable :: hash_head(:), chain_next(:) + !> Generates a random distributions of particles in a box with a minimum 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. + subroutine s_particle_cloud_random_box(cloud_idx, glbl_idx, cloud_ibs, num_cloud_ibs) + + integer, intent(in) :: cloud_idx + integer, intent(inout) :: glbl_idx + 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(8) :: n_attempts, max_attempts + real(wp) :: xmin, xmax, ymin, ymax, zmin, zmax, min_dist + real(wp) :: rx, ry, rz, dist + logical :: overlaps + real(wp), allocatable :: placed(:,:) + integer :: hash_size, slot + integer :: bx, by, bz, nbx, nby, nbz + integer :: dx_b, dy_b, dz_b, dz_lo, dz_hi, j + integer, allocatable :: hash_head(:), chain_next(:) + + allocate (cloud_ibs(particle_cloud(cloud_idx)%num_particles), stat=alloc_stat) + if (alloc_stat /= 0) then + call s_mpi_abort("Error :: Ran out of CPU memory trying to allocate particle cloud IB array. " & + & // "Current system resources cannot perform rejection packing with the specified number of particles.") + end if + ib_idx = 0 xmin = particle_cloud(cloud_idx)%x_centroid - 0.5_wp*particle_cloud(cloud_idx)%length_x xmax = particle_cloud(cloud_idx)%x_centroid + 0.5_wp*particle_cloud(cloud_idx)%length_x @@ -174,7 +192,8 @@ contains chain_next(n_placed) = hash_head(slot) hash_head(slot) = n_placed - call s_add_cloud_particle(cloud_idx, ib_idx, glbl_idx, geom, rx, ry, rz, particle_cloud_ibs) + glbl_idx = glbl_idx + 1 + call s_add_cloud_particle(cloud_idx, ib_idx, glbl_idx, geom, rx, ry, rz, cloud_ibs) end if end do @@ -184,22 +203,32 @@ 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_random_box !> 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. - subroutine s_particle_cloud_lattice(cloud_idx, ib_idx, glbl_idx, particle_cloud_ibs) - - integer, intent(in) :: cloud_idx - integer, intent(inout) :: ib_idx, glbl_idx - type(ib_patch_parameters), intent(inout), dimension(:) :: particle_cloud_ibs - integer :: n_placed, n_target, geom - 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 - real(wp), dimension(4) :: bx_off, by_off, bz_off + !! 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. + subroutine s_particle_cloud_lattice(cloud_idx, glbl_idx, cloud_ibs, num_cloud_ibs) + + integer, intent(in) :: cloud_idx + integer, intent(inout) :: glbl_idx + 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 :: 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 + 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))) + ib_idx = 0 xmin = particle_cloud(cloud_idx)%x_centroid - 0.5_wp*particle_cloud(cloud_idx)%length_x xmax = particle_cloud(cloud_idx)%x_centroid + 0.5_wp*particle_cloud(cloud_idx)%length_x @@ -238,8 +267,12 @@ contains col = 0 px = x0 do while (px <= xmax .and. n_placed < n_target) - call s_add_cloud_particle(cloud_idx, ib_idx, glbl_idx, geom, px, py, particle_cloud(cloud_idx)%z_centroid, & - & particle_cloud_ibs) + 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 n_placed = n_placed + 1 col = col + 1 px = x0 + real(col, wp)*spacing @@ -260,9 +293,13 @@ contains do ix = 0, ncx - 1 do b = 1, 4 if (n_placed >= n_target) exit - call s_add_cloud_particle(cloud_idx, ib_idx, glbl_idx, geom, xmin + real(ix, wp)*cell + bx_off(b), & - & ymin + real(jy, wp)*cell + by_off(b), zmin + real(kz, & - & wp)*cell + bz_off(b), particle_cloud_ibs) + 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 n_placed = n_placed + 1 end do end do @@ -271,29 +308,24 @@ contains end do end if + num_cloud_ibs = ib_idx + end subroutine s_particle_cloud_lattice - !> Writes a single placed particle into particle_cloud_ibs at the next free slot if it falls within this rank's IB neighborhood - !! (get_neighbor_bounds() must already have run), advancing ib_idx; particles outside the neighborhood are discarded here rather - !! than stored. glbl_idx always advances (starting from num_ibs, see s_generate_particle_clouds) so gbl_patch_id is already the - !! particle's final, absolute global patch id - s_reduce_ib_patch_array copies it 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. 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_random_box) 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. subroutine s_add_cloud_particle(cloud_idx, ib_idx, glbl_idx, geom, px, py, pz, particle_cloud_ibs) - integer, intent(in) :: cloud_idx, geom - integer, intent(inout) :: ib_idx, glbl_idx + integer, intent(in) :: cloud_idx, glbl_idx, geom + integer, intent(inout) :: ib_idx real(wp), intent(in) :: px, py, pz type(ib_patch_parameters), intent(inout), dimension(:) :: particle_cloud_ibs - real(wp), dimension(3) :: centroid - - glbl_idx = glbl_idx + 1 - - centroid = [px, py, 0._wp] - if (num_dims == 3) centroid(3) = pz - if (.not. f_neighborhood_ranks_own_location(centroid)) return ib_idx = ib_idx + 1 - @:PROHIBIT(ib_idx > num_ib_patches_max_namelist, & + @: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.") particle_cloud_ibs(ib_idx)%gbl_patch_id = glbl_idx @@ -337,6 +369,29 @@ contains end subroutine s_add_cloud_particle + !> 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_random_box), to pare its full, unfiltered placement down to this rank's neighborhood. + subroutine s_reduce_particle_cloud_ibs(cloud_ibs, num_ibs) + + type(ib_patch_parameters), intent(inout), dimension(:) :: cloud_ibs + integer, intent(inout) :: num_ibs + integer :: i, write_idx + real(wp), dimension(3) :: centroid + + write_idx = 0 + do i = 1, num_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_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) From 2349818f5f440b058eec388164b7b4a4a2c755b5 Mon Sep 17 00:00:00 2001 From: Daniel Vickers Date: Thu, 30 Jul 2026 20:12:29 -0700 Subject: [PATCH 7/8] Update to the tuo mako file for tme limits value in the mako --- toolchain/templates/tuo.mako | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/toolchain/templates/tuo.mako b/toolchain/templates/tuo.mako index 1b3247675e..74c6dc4fc5 100644 --- a/toolchain/templates/tuo.mako +++ b/toolchain/templates/tuo.mako @@ -8,7 +8,7 @@ # flux: --job-name="${name}" # flux: --output="${name}.out" # flux: --error="${name}.err" -# flux: --time=${walltime} +# flux: --time-limit=${walltime} # flux: --exclusive # flux: --setattr=thp=always % if account: From 75c7c61e4ec8f61fe3460ea4b9d02af91d763b66 Mon Sep 17 00:00:00 2001 From: "Daniel J. Vickers" Date: Thu, 30 Jul 2026 23:42:30 -0400 Subject: [PATCH 8/8] Addressing some review comments --- src/common/m_constants.fpp | 5 ++- src/simulation/m_collisions.fpp | 50 ----------------------------- src/simulation/m_particle_cloud.fpp | 23 ++++++------- 3 files changed, 12 insertions(+), 66 deletions(-) diff --git a/src/common/m_constants.fpp b/src/common/m_constants.fpp index 662a8d62fc..e01531eb8d 100644 --- a/src/common/m_constants.fpp +++ b/src/common/m_constants.fpp @@ -27,10 +27,9 @@ module m_constants integer, parameter :: num_ib_airfoils_max = 5 !< Maximum number of ib_airfoil instances integer, parameter :: num_stl_models_max = 10 !> Maximum number of immersed boundary patches (legacy, not used for patch_ib sizing) - integer, parameter :: num_ib_patches_max = 1500000000 !> Fixed capacity of patch_ib (namelist patches + local particle bed subset after reduction) - integer, parameter :: num_ib_patches_max_namelist = 810000 - integer, parameter :: num_local_ibs_max = 30000 !< Maximum number of immersed boundary patches (patch_ib) + integer, parameter :: num_ib_patches_max_namelist = 54000 + integer, parameter :: num_local_ibs_max = 2000 !< Maximum number of immersed boundary patches (patch_ib) integer, parameter :: num_particle_clouds_max = 10 !< Maximum number of particle bed patch specifications integer, parameter :: num_bc_patches_max = 10 !< Maximum number of boundary condition patches integer, parameter :: max_2d_fourier_modes = 10 !< Max Fourier mode index for 2D modal patch (geometry 13) diff --git a/src/simulation/m_collisions.fpp b/src/simulation/m_collisions.fpp index 2a03e9a616..f9ed5b732e 100644 --- a/src/simulation/m_collisions.fpp +++ b/src/simulation/m_collisions.fpp @@ -66,7 +66,6 @@ contains ! get is distance used in the force calculation with each IB and each wall call s_detect_wall_collisions() call s_detect_ib_collisions(ghost_points, ib_markers, num_gps, num_considered_collisions) - ! call s_detect_ib_collisions_n2(num_considered_collisions) select case (collision_model) case (1) ! soft sphere model @@ -334,55 +333,6 @@ contains end subroutine s_detect_ib_collisions - subroutine s_detect_ib_collisions_n2(num_considered_collisions) - - integer, intent(out) :: num_considered_collisions - integer :: pid1, pid2, encoded_pid2, current_collisions - integer :: xp_lower, xp_upper, yp_lower, yp_upper, zp_lower, zp_upper, xp, yp, zp - real(wp), dimension(3) :: centroid_1, centroid_2, distance_vec - - num_considered_collisions = 0 - - call s_get_periodicities(xp_lower, xp_upper, yp_lower, yp_upper, zp_lower, zp_upper) - - $:GPU_PARALLEL_LOOP(private='[pid1, pid2, encoded_pid2, centroid_1, centroid_2, xp, yp, zp, distance_vec, & - & current_collisions]', copyin='[xp_lower, xp_upper, yp_lower, yp_upper, zp_lower, zp_upper]', copy='[num_considered_collisions]') - do pid1 = 1, num_ibs - 1 - centroid_1 = [patch_ib(pid1)%x_centroid, patch_ib(pid1)%y_centroid, 0._wp] - if (num_dims == 3) centroid_1(3) = patch_ib(pid1)%z_centroid - do pid2 = pid1 + 1, num_ibs - periodic_search: do xp = xp_lower, xp_upper - do yp = yp_lower, yp_upper - do zp = zp_lower, zp_upper - centroid_2(1) = patch_ib(pid2)%x_centroid + real(xp, wp)*(glb_bounds(1)%end - glb_bounds(1)%beg) - centroid_2(2) = patch_ib(pid2)%y_centroid + real(yp, wp)*(glb_bounds(2)%end - glb_bounds(2)%beg) - if (num_dims == 3) centroid_2(3) = patch_ib(pid2)%z_centroid + real(zp, & - & wp)*(glb_bounds(3)%end - glb_bounds(3)%beg) - distance_vec = centroid_2 - centroid_1 - - if (norm2(distance_vec) < patch_ib(pid1)%radius + patch_ib(pid2)%radius) then - $:GPU_ATOMIC(atomic='capture') - num_considered_collisions = num_considered_collisions + 1 - current_collisions = num_considered_collisions - $:END_GPU_ATOMIC_CAPTURE() - - call s_encode_patch_periodicity(patch_ib(pid2)%gbl_patch_id, xp, yp, zp, encoded_pid2) - - collision_lookup(current_collisions, 1) = pid1 - collision_lookup(current_collisions, 2) = pid2 - collision_lookup(current_collisions, 3) = patch_ib(pid1)%gbl_patch_id - collision_lookup(current_collisions, 4) = encoded_pid2 - exit periodic_search - end if - end do - end do - end do periodic_search - end do - end do - $:END_GPU_PARALLEL_LOOP() - - end subroutine s_detect_ib_collisions_n2 - !> @brief uses boundary conditions and particle locations to check for wall conditions subroutine s_detect_wall_collisions() diff --git a/src/simulation/m_particle_cloud.fpp b/src/simulation/m_particle_cloud.fpp index b7eb3b9515..51976316f7 100644 --- a/src/simulation/m_particle_cloud.fpp +++ b/src/simulation/m_particle_cloud.fpp @@ -24,12 +24,7 @@ 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_random_box) and - !! hands back only the entries that fall within this rank's IB neighborhood (get_neighbor_bounds() must already have been - !! called); those are copied here into particle_cloud_ibs, allocated once to a worst-case capacity of - !! min(num_ib_patches_max_namelist, total particles requested across all clouds). 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. Each entry's gbl_patch_id is already the final, absolute global patch id (namelist patches occupy 1..num_ibs, - !! so the running index here starts at num_ibs) - s_reduce_ib_patch_array copies it as-is. + !! hands back only the entries that fall within this rank's IB neighborhood impure subroutine s_generate_particle_clouds(particle_cloud_ibs, num_particle_cloud_ibs) type(ib_patch_parameters), allocatable, intent(out), dimension(:) :: particle_cloud_ibs @@ -61,6 +56,8 @@ contains call s_particle_cloud_random_box(cloud_idx, glbl_idx, cloud_ibs, num_cloud_ibs) case (2) ! lattice packing method call s_particle_cloud_lattice(cloud_idx, glbl_idx, 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 @:PROHIBIT(num_particle_cloud_ibs + num_cloud_ibs > num_ib_patches_max_namelist, & @@ -113,7 +110,7 @@ contains min_dist = 2._wp*particle_cloud(cloud_idx)%radius + particle_cloud(cloud_idx)%min_spacing - if (p == 0) then + if (num_dims < 3) then geom = 2 ! circle for 2D dz_lo = 0 dz_hi = 0 @@ -144,7 +141,7 @@ contains rx = xmin + f_xorshift(seed)*(xmax - xmin) ry = ymin + f_xorshift(seed)*(ymax - ymin) - if (p == 0) then + if (num_dims < 3) then rz = particle_cloud(cloud_idx)%z_centroid else rz = zmin + f_xorshift(seed)*(zmax - zmin) @@ -166,7 +163,7 @@ contains slot = f_bin_hash(nbx, nby, nbz, hash_size) j = hash_head(slot) do while (j > 0) - if (p == 0) then + if (num_dims < 3) then dist = sqrt((rx - placed(1, j))**2 + (ry - placed(2, j))**2) else dist = sqrt((rx - placed(1, j))**2 + (ry - placed(2, j))**2 + (rz - placed(3, j))**2) @@ -372,15 +369,15 @@ contains !> 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_random_box), to pare its full, unfiltered placement down to this rank's neighborhood. - subroutine s_reduce_particle_cloud_ibs(cloud_ibs, num_ibs) + subroutine s_reduce_particle_cloud_ibs(cloud_ibs, num_cloud_ibs) type(ib_patch_parameters), intent(inout), dimension(:) :: cloud_ibs - integer, intent(inout) :: num_ibs + integer, intent(inout) :: num_cloud_ibs integer :: i, write_idx real(wp), dimension(3) :: centroid write_idx = 0 - do i = 1, num_ibs + 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 @@ -388,7 +385,7 @@ contains if (write_idx /= i) cloud_ibs(write_idx) = cloud_ibs(i) end if end do - num_ibs = write_idx + num_cloud_ibs = write_idx end subroutine s_reduce_particle_cloud_ibs