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_constants.fpp b/src/common/m_constants.fpp index 4857d93249..e01531eb8d 100644 --- a/src/common/m_constants.fpp +++ b/src/common/m_constants.fpp @@ -27,7 +27,6 @@ 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 !> 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) 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 2c4f93b0bc..2585017969 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(m_ds, n_ds, p_ds) diff --git a/src/simulation/m_collisions.fpp b/src/simulation/m_collisions.fpp index 5498da29c2..f9ed5b732e 100644 --- a/src/simulation/m_collisions.fpp +++ b/src/simulation/m_collisions.fpp @@ -65,8 +65,7 @@ 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) 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_data_output.fpp b/src/simulation/m_data_output.fpp index 998ab07b43..547a336c15 100644 --- a/src/simulation/m_data_output.fpp +++ b/src/simulation/m_data_output.fpp @@ -713,7 +713,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, qbmm_pb=pb_ts(1), qbmm_mv=mv_ts(1)) @@ -898,6 +898,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) @@ -957,7 +961,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_particle_cloud.fpp b/src/simulation/m_particle_cloud.fpp index 081c5b3463..51976316f7 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,59 +22,84 @@ module m_particle_cloud contains - !> Generate all particle beds and fill particle_cloud_ibs. Called on all ranks before s_reduce_ib_patch_array. - 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. 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 + 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 :: cloud_idx, ib_idx, n_total_particles + integer, intent(out) :: num_particle_cloud_ibs + 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 allocate (particle_cloud_ibs(0)) + num_particle_cloud_ibs = 0 return end if 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_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, 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, particle_cloud_ibs) + 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, & + & "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 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 - num_ibs, ' 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) - - integer, intent(in) :: cloud_idx - integer, intent(inout) :: ib_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 @@ -82,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 @@ -113,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) @@ -135,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) @@ -161,7 +189,8 @@ 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) + 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 @@ -171,22 +200,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, particle_cloud_ibs) - - integer, intent(in) :: cloud_idx - integer, intent(inout) :: ib_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 @@ -199,7 +238,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))) @@ -214,7 +253,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 @@ -225,8 +264,12 @@ 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, & - & 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 @@ -247,9 +290,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, 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 @@ -258,21 +305,27 @@ 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, 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, 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(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 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.") - ! 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 @@ -313,6 +366,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_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/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index 7b7657ecee..135001a5db 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -896,19 +896,24 @@ contains if (ib) then block type(ib_patch_parameters), allocatable :: particle_cloud_ibs(:) + integer :: num_particle_cloud_ibs + + 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)) + 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() @@ -1221,12 +1226,16 @@ 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. - subroutine s_reduce_ib_patch_array(particle_cloud_ibs) + !> @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. + 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 @@ -1237,7 +1246,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. + ! 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 @@ -1246,15 +1255,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 @@ -1264,9 +1272,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_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 + i end do num_ibs = num_gbl_ibs num_local_ibs = num_gbl_ibs @@ -1286,35 +1293,34 @@ 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 end do - do i = 1, num_bed_ibs + ! 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 - 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) + 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_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 + i end do num_ibs = num_gbl_ibs num_local_ibs = num_gbl_ibs 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: