From fc6a308dae7a9ef143bde801b1c7a46dff146c20 Mon Sep 17 00:00:00 2001 From: danieljvickers Date: Tue, 18 Aug 2026 10:50:51 -0400 Subject: [PATCH 1/6] Neighborhood radius now starts at 0. Documentation updated. There is a bounds check to get the temporary values --- docs/documentation/case.md | 2 +- src/simulation/m_global_parameters.fpp | 2 +- src/simulation/m_ib_patches.fpp | 91 ++++++++++------ src/simulation/m_start_up.fpp | 139 +++++++++++++++---------- toolchain/mfc/params/definitions.py | 2 +- 5 files changed, 148 insertions(+), 88 deletions(-) diff --git a/docs/documentation/case.md b/docs/documentation/case.md index ec82b5bbd8..c7a8b71e4c 100644 --- a/docs/documentation/case.md +++ b/docs/documentation/case.md @@ -416,7 +416,7 @@ Additional details on this specification can be found in [NACA airfoil](https:// - `ib_coefficient_of_friction` is the coefficient of friction used in IB collisions. -- `ib_neighborhood_radius` controls the size of the neighborhood size. This value defaults to 1, which indicates that any given rank is aware of IBs up to 1 ranks away. This parameter is required to strong-scale a case when IBs eventually grow to be larger than one full processor domain wide. +- `ib_neighborhood_radius` controls the size of the neighborhood size. A value of $r$ indicates that any given rank is aware of IBs up to $r$ ranks away. This value defaults to 0, which leaves the radius unset so that it is selected automatically. This parameter is required to strong-scale a case when IBs eventually grow to be larger than one full processor domain wide. #### Particle Clouds diff --git a/src/simulation/m_global_parameters.fpp b/src/simulation/m_global_parameters.fpp index 6ec9c3107b..eb03727812 100644 --- a/src/simulation/m_global_parameters.fpp +++ b/src/simulation/m_global_parameters.fpp @@ -479,7 +479,7 @@ contains bub_pp%R_g = dflt_real; R_g = dflt_real ! Immersed Boundaries (sim-specific extras) - ib_neighborhood_radius = 1 + ib_neighborhood_radius = 0 collision_model = 0 coefficient_of_restitution = dflt_real collision_time = dflt_real diff --git a/src/simulation/m_ib_patches.fpp b/src/simulation/m_ib_patches.fpp index 35b97ea31d..e271284867 100644 --- a/src/simulation/m_ib_patches.fpp +++ b/src/simulation/m_ib_patches.fpp @@ -65,7 +65,7 @@ contains call s_encode_patch_periodicity(patch_ib(patch_id)%gbl_patch_id, xp, yp, zp, encoded_patch_id) ! find the indices to the left and right of the IB in i, j, k - call get_bounding_indices(patch_ib(patch_id), center, il, ir, jl, jr, kl, kr) + call s_get_bounding_indices(patch_ib(patch_id), center, il, ir, jl, jr, kl, kr) ! skip patches whose bounding box does not overlap this rank's domain if (ir < il .or. jr < jl .or. kr < kl) cycle @@ -137,7 +137,7 @@ contains call s_encode_patch_periodicity(patch_ib(patch_id)%gbl_patch_id, xp, yp, 0, encoded_patch_id) ! find the indices to the left and right of the IB in i, j, k - call get_bounding_indices(patch_ib(patch_id), center, il, ir, jl, jr, kl, kr) + call s_get_bounding_indices(patch_ib(patch_id), center, il, ir, jl, jr, kl, kr) ! skip patches whose bounding box does not overlap this rank's domain if (ir < il .or. jr < jl) cycle @@ -218,7 +218,7 @@ contains call s_encode_patch_periodicity(patch_ib(patch_id)%gbl_patch_id, xp, yp, zp, encoded_patch_id) ! find the indices to the left and right of the IB in i, j, k - call get_bounding_indices(patch_ib(patch_id), center, il, ir, jl, jr, kl, kr) + call s_get_bounding_indices(patch_ib(patch_id), center, il, ir, jl, jr, kl, kr) do k = kl, kr do j = jl, jr @@ -287,7 +287,7 @@ contains call s_encode_patch_periodicity(patch_ib(patch_id)%gbl_patch_id, xp, yp, 0, encoded_patch_id) ! find the indices to the left and right of the IB in i, j, k - call get_bounding_indices(patch_ib(patch_id), center, il, ir, jl, jr, kl, kr) + call s_get_bounding_indices(patch_ib(patch_id), center, il, ir, jl, jr, kl, kr) do j = jl, jr do i = il, ir @@ -458,31 +458,67 @@ contains end subroutine s_update_ib_rotation_matrix - subroutine get_bounding_indices(patch, center, il, ir, jl, jr, kl, kr) + subroutine s_get_ib_bound(patch, bound) $:GPU_ROUTINE(parallelism='[seq]') type(ib_patch_parameters), intent(in) :: patch - real(wp), dimension(3), intent(in) :: center - integer, intent(out) :: il, ir, jl, jr, kl, kr - real(wp), dimension(3) :: bbox_min, bbox_max, local_corner, world_corner + real(wp), intent(out) :: bound real(wp), dimension(2) :: lx, ly, lz - integer :: cx, cy, cz - logical :: outside_domain if (patch%geometry == 2 .or. patch%geometry == 8) then ! circle and sphere geometries - bbox_min = center - patch%radius - bbox_max = center + patch%radius + bound = patch%radius else if (patch%geometry == 3) then - ! rectangular geometries - bbox_min = center - 0.5_wp*sqrt(patch%length_x**2 + patch%length_y**2) - bbox_max = center + 0.5_wp*sqrt(patch%length_x**2 + patch%length_y**2) + bound = 0.5_wp*sqrt(patch%length_x**2 + patch%length_y**2) else if (patch%geometry == 4 .or. patch%geometry == 11) then - ! airfoil geometries TODO :: This can be better optimized, since airfoils are typically very long in one dimension - bbox_min = center - ib_airfoil(patch%airfoil_id)%c - bbox_max = center + ib_airfoil(patch%airfoil_id)%c + ! rectangular geometries + bound = ib_airfoil(patch%airfoil_id)%c else if (patch%geometry == 5) then + ! STL model geometry + lx(1) = stl_bounding_boxes(patch%model_id, 1, 1) + lx(2) = stl_bounding_boxes(patch%model_id, 1, 3) + ly(1) = stl_bounding_boxes(patch%model_id, 2, 1) + ly(2) = stl_bounding_boxes(patch%model_id, 2, 3) + + bound = sqrt((lx(2) - lx(1))**2 + (ly(2) - ly(1))**2) + else if (patch%geometry == 6) then + ! ellipse geometry + bound = 0.5_wp*max(patch%length_x, patch%length_y) + else if (patch%geometry == 9) then + ! cuboid geometries + bound = 0.5_wp*sqrt(patch%length_x**2 + patch%length_y**2 + patch%length_z**2) + else if (patch%geometry == 10) then + ! cylinder geometry + bound = sqrt(patch%radius**2 + patch%length_x**2) + else if (patch%geometry == 12) then + ! Local-space bounding box extents (min=1, max=2 in the third index) + lx(1) = stl_bounding_boxes(patch%model_id, 1, 1) + patch%centroid_offset(1) + lx(2) = stl_bounding_boxes(patch%model_id, 1, 3) + patch%centroid_offset(1) + ly(1) = stl_bounding_boxes(patch%model_id, 2, 1) + patch%centroid_offset(2) + ly(2) = stl_bounding_boxes(patch%model_id, 2, 3) + patch%centroid_offset(2) + lz(1) = stl_bounding_boxes(patch%model_id, 3, 1) + patch%centroid_offset(3) + lz(2) = stl_bounding_boxes(patch%model_id, 3, 3) + patch%centroid_offset(3) + + bound = sqrt((lx(2) - lx(1))**2 + (ly(2) - ly(1))**2 + (lz(2) - lz(1))**2) + end if + + end subroutine s_get_ib_bound + + subroutine s_get_bounding_indices(patch, center, il, ir, jl, jr, kl, kr) + + $:GPU_ROUTINE(parallelism='[seq]') + + type(ib_patch_parameters), intent(in) :: patch + real(wp), dimension(3), intent(in) :: center + integer, intent(out) :: il, ir, jl, jr, kl, kr + real(wp), dimension(3) :: bbox_min, bbox_max, local_corner, world_corner + real(wp), dimension(2) :: lx, ly, lz + real(wp) :: bound + integer :: cx, cy, cz + logical :: outside_domain + + if (patch%geometry == 5) then ! STL model geometry lx(1) = stl_bounding_boxes(patch%model_id, 1, 1) + patch%centroid_offset(1) lx(2) = stl_bounding_boxes(patch%model_id, 1, 3) + patch%centroid_offset(1) @@ -502,18 +538,6 @@ contains bbox_max(2) = max(bbox_max(2), world_corner(2)) end do end do - else if (patch%geometry == 6) then - ! ellipse geometry - bbox_min = center - 0.5_wp*max(patch%length_x, patch%length_y) - bbox_max = center + 0.5_wp*max(patch%length_x, patch%length_y) - else if (patch%geometry == 9) then - ! cuboid geometries - bbox_min = center - 0.5_wp*sqrt(patch%length_x**2 + patch%length_y**2 + patch%length_z**2) - bbox_max = center + 0.5_wp*sqrt(patch%length_x**2 + patch%length_y**2 + patch%length_z**2) - else if (patch%geometry == 10) then - ! cylinder geometry - bbox_min = center - sqrt(patch%radius**2 + patch%length_x**2) - bbox_max = center + sqrt(patch%radius**2 + patch%length_x**2) else if (patch%geometry == 12) then ! Local-space bounding box extents (min=1, max=2 in the third index) lx(1) = stl_bounding_boxes(patch%model_id, 1, 1) + patch%centroid_offset(1) @@ -540,6 +564,11 @@ contains end do end do end do + else + ! All other IBs + call s_get_ib_bound(patch, bound) + bbox_min = center - bound + bbox_max = center + bound end if ! completely skip patches whose bounding box does not overlap this rank's domain @@ -566,7 +595,7 @@ contains call get_indices_from_bounds(bbox_min(2), bbox_max(2), y_cc, jl, jr) if (num_dims == 3) call get_indices_from_bounds(bbox_min(3), bbox_max(3), z_cc, kl, kr) - end subroutine get_bounding_indices + end subroutine s_get_bounding_indices subroutine get_indices_from_bounds(left_bound, right_bound, cell_centers, left_index, right_index) diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index 488c858535..0533fb8e01 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -1310,32 +1310,54 @@ contains !! neighbors are known from bc_*, edge neighbors are obtained in round 1, and (3D) corner neighbors in round 2. subroutine s_compute_ib_neighbor_ranks() - integer :: ax, k, nbr_idx, nreqs, sx, sy, sz, dx, dy, dz - integer, allocatable :: send_table(:,:,:), recv_tables(:,:,:,:) + integer :: i, j, k, nbr_idx, nreqs, sx, sy, sz, dx, dy, dz + integer, allocatable :: send_table(:,:,:), recv_tables(:,:,:,:), temporary_neighbors integer, dimension(52) :: requests #ifdef MFC_MPI integer :: ierr integer, dimension(4) :: buf4, rbuf4 integer, dimension(2) :: buf2, rbuf2 + logical :: automatic_neighborhood_radius + real(wp) :: maximum_ib_bound + + ! perform setup if we are doing automatic radius checking + if (ib_neighborhood_radius < 1) then + ! set up a temporary radius array to hold the indices + automatic_neighborhood_radius = .true. + temporary_radius = max(num_procs_x - 1, 1) ! ensure that we get at least 1 neighbor rank + if (num_dims > 1) temporary_radius = max(temporary_radius, num_procs_y - 1) + if (num_dims > 2) temporary_radius = max(temporary_radius, num_procs_z - 1) + + ! determine the maximum length of space that needs to be contained by the neighborhood + maximum_ib_bound = -1._wp + do i = 1, num_ibs + call s_get_ib_bound(patch_ib(i), bound) + maximum_ib_bound = max(maximum_ib_bound, bound) + end do + do i = 1, num_particle_clouds + maximum_ib_bound = max(maximum_ib_bound, particle_cloud(i)%radius) + end do + else + automatic_neighborhood_radius = .false. + temporary_radius = ib_neighborhood_radius + end if - ax = ib_neighborhood_radius - - if (allocated(ib_neighbor_ranks)) deallocate (ib_neighbor_ranks) - allocate (ib_neighbor_ranks(-ax:ax,-ax:ax,-ax:ax)) - ib_neighbor_ranks = MPI_PROC_NULL - ib_neighbor_ranks(0, 0, 0) = proc_rank + allocate (temporary_neighbors(-ib_neighborhood_radius:ib_neighborhood_radius, & + & -ib_neighborhood_radius:ib_neighborhood_radius,-ib_neighborhood_radius:ib_neighborhood_radius)) + temporary_neighbors = MPI_PROC_NULL + temporary_neighbors(0, 0, 0) = proc_rank ! Fill radius-1 entries: face neighbors are known from domain decomposition - ib_neighbor_ranks(-1, 0, 0) = bc_x%beg - ib_neighbor_ranks(+1, 0, 0) = bc_x%end + temporary_neighbors(-1, 0, 0) = bc_x%beg + temporary_neighbors(+1, 0, 0) = bc_x%end if (num_dims >= 2) then - ib_neighbor_ranks(0, -1, 0) = bc_y%beg - ib_neighbor_ranks(0, +1, 0) = bc_y%end + temporary_neighbors(0, -1, 0) = bc_y%beg + temporary_neighbors(0, +1, 0) = bc_y%end end if if (num_dims == 3) then - ib_neighbor_ranks(0, 0, -1) = bc_z%beg - ib_neighbor_ranks(0, 0, +1) = bc_z%end + temporary_neighbors(0, 0, -1) = bc_z%beg + temporary_neighbors(0, 0, +1) = bc_z%end end if if (num_dims >= 2) then @@ -1346,19 +1368,19 @@ contains call MPI_SENDRECV(buf4, 4, MPI_INTEGER, merge(bc_x%beg, MPI_PROC_NULL, bc_x%beg >= 0), 310, rbuf4, 4, MPI_INTEGER, & & merge(bc_x%end, MPI_PROC_NULL, bc_x%end >= 0), 310, MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr) if (bc_x%end >= 0) then - ib_neighbor_ranks(+1, -1, 0) = rbuf4(1) - ib_neighbor_ranks(+1, +1, 0) = rbuf4(2) - ib_neighbor_ranks(+1, 0, -1) = rbuf4(3) - ib_neighbor_ranks(+1, 0, +1) = rbuf4(4) + temporary_neighbors(+1, -1, 0) = rbuf4(1) + temporary_neighbors(+1, +1, 0) = rbuf4(2) + temporary_neighbors(+1, 0, -1) = rbuf4(3) + temporary_neighbors(+1, 0, +1) = rbuf4(4) end if call MPI_SENDRECV(buf4, 4, MPI_INTEGER, merge(bc_x%end, MPI_PROC_NULL, bc_x%end >= 0), 311, rbuf4, 4, MPI_INTEGER, & & merge(bc_x%beg, MPI_PROC_NULL, bc_x%beg >= 0), 311, MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr) if (bc_x%beg >= 0) then - ib_neighbor_ranks(-1, -1, 0) = rbuf4(1) - ib_neighbor_ranks(-1, +1, 0) = rbuf4(2) - ib_neighbor_ranks(-1, 0, -1) = rbuf4(3) - ib_neighbor_ranks(-1, 0, +1) = rbuf4(4) + temporary_neighbors(-1, -1, 0) = rbuf4(1) + temporary_neighbors(-1, +1, 0) = rbuf4(2) + temporary_neighbors(-1, 0, -1) = rbuf4(3) + temporary_neighbors(-1, 0, +1) = rbuf4(4) end if end if @@ -1369,41 +1391,43 @@ contains call MPI_SENDRECV(buf2, 2, MPI_INTEGER, merge(bc_y%beg, MPI_PROC_NULL, bc_y%beg >= 0), 312, rbuf2, 2, MPI_INTEGER, & & merge(bc_y%end, MPI_PROC_NULL, bc_y%end >= 0), 312, MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr) if (bc_y%end >= 0) then - ib_neighbor_ranks(0, +1, -1) = rbuf2(1) - ib_neighbor_ranks(0, +1, +1) = rbuf2(2) + temporary_neighbors(0, +1, -1) = rbuf2(1) + temporary_neighbors(0, +1, +1) = rbuf2(2) end if call MPI_SENDRECV(buf2, 2, MPI_INTEGER, merge(bc_y%end, MPI_PROC_NULL, bc_y%end >= 0), 313, rbuf2, 2, MPI_INTEGER, & & merge(bc_y%beg, MPI_PROC_NULL, bc_y%beg >= 0), 313, MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr) if (bc_y%beg >= 0) then - ib_neighbor_ranks(0, -1, -1) = rbuf2(1) - ib_neighbor_ranks(0, -1, +1) = rbuf2(2) + temporary_neighbors(0, -1, -1) = rbuf2(1) + temporary_neighbors(0, -1, +1) = rbuf2(2) end if ! Round 2: exchange z face ranks with xy-diagonal edge neighbors -> corner ranks. Each of the 4 xy diagonals gives 2 ! corners (the +/-z variants). Pattern: send buf2 to mirror diagonal, receive from this diagonal -> that edge's z face ! ranks. #:for DX, DY, MDX, MDY, TAG in [(1,1,-1,-1,320), (1,-1,-1,1,321), (-1,1,1,-1,322), (-1,-1,1,1,323)] - call MPI_SENDRECV(buf2, 2, MPI_INTEGER, merge(ib_neighbor_ranks(${MDX}$, ${MDY}$, 0), MPI_PROC_NULL, & - & ib_neighbor_ranks(${MDX}$, ${MDY}$, 0) >= 0), ${TAG}$, rbuf2, 2, MPI_INTEGER, & - & merge(ib_neighbor_ranks(${DX}$, ${DY}$, 0), MPI_PROC_NULL, ib_neighbor_ranks(${DX}$, ${DY}$, & - & 0) >= 0), ${TAG}$, MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr) - if (ib_neighbor_ranks(${DX}$, ${DY}$, 0) >= 0) then - ib_neighbor_ranks(${DX}$, ${DY}$, -1) = rbuf2(1) - ib_neighbor_ranks(${DX}$, ${DY}$, +1) = rbuf2(2) + call MPI_SENDRECV(buf2, 2, MPI_INTEGER, merge(temporary_neighbors(${MDX}$, ${MDY}$, 0), MPI_PROC_NULL, & + & temporary_neighbors(${MDX}$, ${MDY}$, 0) >= 0), ${TAG}$, rbuf2, 2, MPI_INTEGER, & + & merge(temporary_neighbors(${DX}$, ${DY}$, 0), MPI_PROC_NULL, temporary_neighbors(${DX}$, & + & ${DY}$, 0) >= 0), ${TAG}$, MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr) + if (temporary_neighbors(${DX}$, ${DY}$, 0) >= 0) then + temporary_neighbors(${DX}$, ${DY}$, -1) = rbuf2(1) + temporary_neighbors(${DX}$, ${DY}$, +1) = rbuf2(2) end if #:endfor end if ! For radius > 1: extend the table by iterative 26-neighbor full-table exchanges. In each round, every rank broadcasts its ! current table to all 26 immediate neighbors. Their entry at offset (dx,dy,dz) from them = our entry at - ! (dx+sx,dy+sy,dz+sz). One extension round fills the entire next shell, so ax-1 rounds suffice. - if (ax > 1) then - allocate (send_table(-ax:ax,-ax:ax,-ax:ax)) - allocate (recv_tables(-ax:ax,-ax:ax,-ax:ax,1:26)) + ! (dx+sx,dy+sy,dz+sz). One extension round fills the entire next shell, so ib_neighborhood_radius-1 rounds suffice. + if (ib_neighborhood_radius > 1) then + allocate (send_table(-ib_neighborhood_radius:ib_neighborhood_radius,-ib_neighborhood_radius:ib_neighborhood_radius, & + & -ib_neighborhood_radius:ib_neighborhood_radius)) + allocate (recv_tables(-ib_neighborhood_radius:ib_neighborhood_radius,-ib_neighborhood_radius:ib_neighborhood_radius, & + & -ib_neighborhood_radius:ib_neighborhood_radius,1:26)) - do k = 2, ax - send_table = ib_neighbor_ranks + do k = 2, ib_neighborhood_radius + send_table = temporary_neighbors nreqs = 0 nbr_idx = 0 @@ -1412,10 +1436,10 @@ contains do sx = -1, 1 if (sx == 0 .and. sy == 0 .and. sz == 0) cycle nbr_idx = nbr_idx + 1 - if (ib_neighbor_ranks(sx, sy, sz) < 0) cycle + if (temporary_neighbors(sx, sy, sz) < 0) cycle nreqs = nreqs + 1 - call MPI_IRECV(recv_tables(:,:,:,nbr_idx), (2*ax + 1)**3, MPI_INTEGER, ib_neighbor_ranks(sx, sy, sz), & - & 400, MPI_COMM_WORLD, requests(nreqs), ierr) + call MPI_IRECV(recv_tables(:,:,:,nbr_idx), (2*ib_neighborhood_radius + 1)**3, MPI_INTEGER, & + & temporary_neighbors(sx, sy, sz), 400, MPI_COMM_WORLD, requests(nreqs), ierr) end do end do end do @@ -1424,10 +1448,10 @@ contains do sy = -1, 1 do sx = -1, 1 if (sx == 0 .and. sy == 0 .and. sz == 0) cycle - if (ib_neighbor_ranks(sx, sy, sz) < 0) cycle + if (temporary_neighbors(sx, sy, sz) < 0) cycle nreqs = nreqs + 1 - call MPI_ISEND(send_table, (2*ax + 1)**3, MPI_INTEGER, ib_neighbor_ranks(sx, sy, sz), 400, & - & MPI_COMM_WORLD, requests(nreqs), ierr) + call MPI_ISEND(send_table, (2*ib_neighborhood_radius + 1)**3, MPI_INTEGER, temporary_neighbors(sx, & + & sy, sz), 400, MPI_COMM_WORLD, requests(nreqs), ierr) end do end do end do @@ -1440,16 +1464,16 @@ contains do sx = -1, 1 if (sx == 0 .and. sy == 0 .and. sz == 0) cycle nbr_idx = nbr_idx + 1 - if (ib_neighbor_ranks(sx, sy, sz) < 0) cycle - do dz = -ax, ax - do dy = -ax, ax - do dx = -ax, ax + if (temporary_neighbors(sx, sy, sz) < 0) cycle + do dz = -ib_neighborhood_radius, ib_neighborhood_radius + do dy = -ib_neighborhood_radius, ib_neighborhood_radius + do dx = -ib_neighborhood_radius, ib_neighborhood_radius if (recv_tables(dx, dy, dz, nbr_idx) == MPI_PROC_NULL) cycle - if (dx + sx < -ax .or. dx + sx > ax) cycle - if (dy + sy < -ax .or. dy + sy > ax) cycle - if (dz + sz < -ax .or. dz + sz > ax) cycle - if (ib_neighbor_ranks(dx + sx, dy + sy, dz + sz) /= MPI_PROC_NULL) cycle - ib_neighbor_ranks(dx + sx, dy + sy, dz + sz) = recv_tables(dx, dy, dz, nbr_idx) + if (dx + sx < -ib_neighborhood_radius .or. dx + sx > ib_neighborhood_radius) cycle + if (dy + sy < -ib_neighborhood_radius .or. dy + sy > ib_neighborhood_radius) cycle + if (dz + sz < -ib_neighborhood_radius .or. dz + sz > ib_neighborhood_radius) cycle + if (temporary_neighbors(dx + sx, dy + sy, dz + sz) /= MPI_PROC_NULL) cycle + temporary_neighbors(dx + sx, dy + sy, dz + sz) = recv_tables(dx, dy, dz, nbr_idx) end do end do end do @@ -1460,6 +1484,13 @@ contains deallocate (send_table, recv_tables) end if + + if (allocated(ib_neighbor_ranks)) deallocate (ib_neighbor_ranks) + allocate (ib_neighbor_ranks(-ib_neighborhood_radius:ib_neighborhood_radius,-ib_neighborhood_radius:ib_neighborhood_radius, & + & -ib_neighborhood_radius:ib_neighborhood_radius)) + ib_neighbor_ranks = temporary_neighbors(-ib_neighborhood_radius, ib_neighborhood_radius, -ib_neighborhood_radius, & + & ib_neighborhood_radius, -ib_neighborhood_radius, ib_neighborhood_radius) + deallocate (temporary_neighbors) #endif end subroutine s_compute_ib_neighbor_ranks diff --git a/toolchain/mfc/params/definitions.py b/toolchain/mfc/params/definitions.py index 886d4252db..5675e4f7cd 100644 --- a/toolchain/mfc/params/definitions.py +++ b/toolchain/mfc/params/definitions.py @@ -377,7 +377,7 @@ def get_value_label(param_name: str, value: int) -> str: "num_fluids": {"min": 1, "max": NF}, "num_patches": {"min": 0, "max": NUM_PATCHES_MAX}, "num_ibs": {"min": 0}, - "ib_neighborhood_radius": {"min": 1}, + "ib_neighborhood_radius": {"min": 0}, "num_source": {"min": 1}, "num_probes": {"min": 1}, "num_integrals": {"min": 1}, From 36507d63be57b0d5cdf84662ba419e133e8df1cd Mon Sep 17 00:00:00 2001 From: danieljvickers Date: Tue, 18 Aug 2026 14:36:08 -0400 Subject: [PATCH 2/6] Initial Implementation done --- src/simulation/m_start_up.fpp | 178 +++++++++++++++++----------------- 1 file changed, 89 insertions(+), 89 deletions(-) diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index 0533fb8e01..2495cd5411 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -1233,7 +1233,7 @@ contains end do end if - call get_neighbor_bounds() + call s_get_neighbor_bounds() call s_compute_ib_neighbor_ranks() num_gbl_ibs = num_namelist_ibs + num_bed_ibs @@ -1310,54 +1310,32 @@ contains !! neighbors are known from bc_*, edge neighbors are obtained in round 1, and (3D) corner neighbors in round 2. subroutine s_compute_ib_neighbor_ranks() - integer :: i, j, k, nbr_idx, nreqs, sx, sy, sz, dx, dy, dz - integer, allocatable :: send_table(:,:,:), recv_tables(:,:,:,:), temporary_neighbors + integer :: ax, k, nbr_idx, nreqs, sx, sy, sz, dx, dy, dz + integer, allocatable :: send_table(:,:,:), recv_tables(:,:,:,:) integer, dimension(52) :: requests #ifdef MFC_MPI integer :: ierr integer, dimension(4) :: buf4, rbuf4 integer, dimension(2) :: buf2, rbuf2 - logical :: automatic_neighborhood_radius - real(wp) :: maximum_ib_bound - ! perform setup if we are doing automatic radius checking - if (ib_neighborhood_radius < 1) then - ! set up a temporary radius array to hold the indices - automatic_neighborhood_radius = .true. - temporary_radius = max(num_procs_x - 1, 1) ! ensure that we get at least 1 neighbor rank - if (num_dims > 1) temporary_radius = max(temporary_radius, num_procs_y - 1) - if (num_dims > 2) temporary_radius = max(temporary_radius, num_procs_z - 1) - - ! determine the maximum length of space that needs to be contained by the neighborhood - maximum_ib_bound = -1._wp - do i = 1, num_ibs - call s_get_ib_bound(patch_ib(i), bound) - maximum_ib_bound = max(maximum_ib_bound, bound) - end do - do i = 1, num_particle_clouds - maximum_ib_bound = max(maximum_ib_bound, particle_cloud(i)%radius) - end do - else - automatic_neighborhood_radius = .false. - temporary_radius = ib_neighborhood_radius - end if + ax = ib_neighborhood_radius - allocate (temporary_neighbors(-ib_neighborhood_radius:ib_neighborhood_radius, & - & -ib_neighborhood_radius:ib_neighborhood_radius,-ib_neighborhood_radius:ib_neighborhood_radius)) - temporary_neighbors = MPI_PROC_NULL - temporary_neighbors(0, 0, 0) = proc_rank + if (allocated(ib_neighbor_ranks)) deallocate (ib_neighbor_ranks) + allocate (ib_neighbor_ranks(-ax:ax,-ax:ax,-ax:ax)) + ib_neighbor_ranks = MPI_PROC_NULL + ib_neighbor_ranks(0, 0, 0) = proc_rank ! Fill radius-1 entries: face neighbors are known from domain decomposition - temporary_neighbors(-1, 0, 0) = bc_x%beg - temporary_neighbors(+1, 0, 0) = bc_x%end + ib_neighbor_ranks(-1, 0, 0) = bc_x%beg + ib_neighbor_ranks(+1, 0, 0) = bc_x%end if (num_dims >= 2) then - temporary_neighbors(0, -1, 0) = bc_y%beg - temporary_neighbors(0, +1, 0) = bc_y%end + ib_neighbor_ranks(0, -1, 0) = bc_y%beg + ib_neighbor_ranks(0, +1, 0) = bc_y%end end if if (num_dims == 3) then - temporary_neighbors(0, 0, -1) = bc_z%beg - temporary_neighbors(0, 0, +1) = bc_z%end + ib_neighbor_ranks(0, 0, -1) = bc_z%beg + ib_neighbor_ranks(0, 0, +1) = bc_z%end end if if (num_dims >= 2) then @@ -1368,19 +1346,19 @@ contains call MPI_SENDRECV(buf4, 4, MPI_INTEGER, merge(bc_x%beg, MPI_PROC_NULL, bc_x%beg >= 0), 310, rbuf4, 4, MPI_INTEGER, & & merge(bc_x%end, MPI_PROC_NULL, bc_x%end >= 0), 310, MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr) if (bc_x%end >= 0) then - temporary_neighbors(+1, -1, 0) = rbuf4(1) - temporary_neighbors(+1, +1, 0) = rbuf4(2) - temporary_neighbors(+1, 0, -1) = rbuf4(3) - temporary_neighbors(+1, 0, +1) = rbuf4(4) + ib_neighbor_ranks(+1, -1, 0) = rbuf4(1) + ib_neighbor_ranks(+1, +1, 0) = rbuf4(2) + ib_neighbor_ranks(+1, 0, -1) = rbuf4(3) + ib_neighbor_ranks(+1, 0, +1) = rbuf4(4) end if call MPI_SENDRECV(buf4, 4, MPI_INTEGER, merge(bc_x%end, MPI_PROC_NULL, bc_x%end >= 0), 311, rbuf4, 4, MPI_INTEGER, & & merge(bc_x%beg, MPI_PROC_NULL, bc_x%beg >= 0), 311, MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr) if (bc_x%beg >= 0) then - temporary_neighbors(-1, -1, 0) = rbuf4(1) - temporary_neighbors(-1, +1, 0) = rbuf4(2) - temporary_neighbors(-1, 0, -1) = rbuf4(3) - temporary_neighbors(-1, 0, +1) = rbuf4(4) + ib_neighbor_ranks(-1, -1, 0) = rbuf4(1) + ib_neighbor_ranks(-1, +1, 0) = rbuf4(2) + ib_neighbor_ranks(-1, 0, -1) = rbuf4(3) + ib_neighbor_ranks(-1, 0, +1) = rbuf4(4) end if end if @@ -1391,43 +1369,41 @@ contains call MPI_SENDRECV(buf2, 2, MPI_INTEGER, merge(bc_y%beg, MPI_PROC_NULL, bc_y%beg >= 0), 312, rbuf2, 2, MPI_INTEGER, & & merge(bc_y%end, MPI_PROC_NULL, bc_y%end >= 0), 312, MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr) if (bc_y%end >= 0) then - temporary_neighbors(0, +1, -1) = rbuf2(1) - temporary_neighbors(0, +1, +1) = rbuf2(2) + ib_neighbor_ranks(0, +1, -1) = rbuf2(1) + ib_neighbor_ranks(0, +1, +1) = rbuf2(2) end if call MPI_SENDRECV(buf2, 2, MPI_INTEGER, merge(bc_y%end, MPI_PROC_NULL, bc_y%end >= 0), 313, rbuf2, 2, MPI_INTEGER, & & merge(bc_y%beg, MPI_PROC_NULL, bc_y%beg >= 0), 313, MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr) if (bc_y%beg >= 0) then - temporary_neighbors(0, -1, -1) = rbuf2(1) - temporary_neighbors(0, -1, +1) = rbuf2(2) + ib_neighbor_ranks(0, -1, -1) = rbuf2(1) + ib_neighbor_ranks(0, -1, +1) = rbuf2(2) end if ! Round 2: exchange z face ranks with xy-diagonal edge neighbors -> corner ranks. Each of the 4 xy diagonals gives 2 ! corners (the +/-z variants). Pattern: send buf2 to mirror diagonal, receive from this diagonal -> that edge's z face ! ranks. #:for DX, DY, MDX, MDY, TAG in [(1,1,-1,-1,320), (1,-1,-1,1,321), (-1,1,1,-1,322), (-1,-1,1,1,323)] - call MPI_SENDRECV(buf2, 2, MPI_INTEGER, merge(temporary_neighbors(${MDX}$, ${MDY}$, 0), MPI_PROC_NULL, & - & temporary_neighbors(${MDX}$, ${MDY}$, 0) >= 0), ${TAG}$, rbuf2, 2, MPI_INTEGER, & - & merge(temporary_neighbors(${DX}$, ${DY}$, 0), MPI_PROC_NULL, temporary_neighbors(${DX}$, & - & ${DY}$, 0) >= 0), ${TAG}$, MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr) - if (temporary_neighbors(${DX}$, ${DY}$, 0) >= 0) then - temporary_neighbors(${DX}$, ${DY}$, -1) = rbuf2(1) - temporary_neighbors(${DX}$, ${DY}$, +1) = rbuf2(2) + call MPI_SENDRECV(buf2, 2, MPI_INTEGER, merge(ib_neighbor_ranks(${MDX}$, ${MDY}$, 0), MPI_PROC_NULL, & + & ib_neighbor_ranks(${MDX}$, ${MDY}$, 0) >= 0), ${TAG}$, rbuf2, 2, MPI_INTEGER, & + & merge(ib_neighbor_ranks(${DX}$, ${DY}$, 0), MPI_PROC_NULL, ib_neighbor_ranks(${DX}$, ${DY}$, & + & 0) >= 0), ${TAG}$, MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr) + if (ib_neighbor_ranks(${DX}$, ${DY}$, 0) >= 0) then + ib_neighbor_ranks(${DX}$, ${DY}$, -1) = rbuf2(1) + ib_neighbor_ranks(${DX}$, ${DY}$, +1) = rbuf2(2) end if #:endfor end if ! For radius > 1: extend the table by iterative 26-neighbor full-table exchanges. In each round, every rank broadcasts its ! current table to all 26 immediate neighbors. Their entry at offset (dx,dy,dz) from them = our entry at - ! (dx+sx,dy+sy,dz+sz). One extension round fills the entire next shell, so ib_neighborhood_radius-1 rounds suffice. - if (ib_neighborhood_radius > 1) then - allocate (send_table(-ib_neighborhood_radius:ib_neighborhood_radius,-ib_neighborhood_radius:ib_neighborhood_radius, & - & -ib_neighborhood_radius:ib_neighborhood_radius)) - allocate (recv_tables(-ib_neighborhood_radius:ib_neighborhood_radius,-ib_neighborhood_radius:ib_neighborhood_radius, & - & -ib_neighborhood_radius:ib_neighborhood_radius,1:26)) + ! (dx+sx,dy+sy,dz+sz). One extension round fills the entire next shell, so ax-1 rounds suffice. + if (ax > 1) then + allocate (send_table(-ax:ax,-ax:ax,-ax:ax)) + allocate (recv_tables(-ax:ax,-ax:ax,-ax:ax,1:26)) - do k = 2, ib_neighborhood_radius - send_table = temporary_neighbors + do k = 2, ax + send_table = ib_neighbor_ranks nreqs = 0 nbr_idx = 0 @@ -1436,10 +1412,10 @@ contains do sx = -1, 1 if (sx == 0 .and. sy == 0 .and. sz == 0) cycle nbr_idx = nbr_idx + 1 - if (temporary_neighbors(sx, sy, sz) < 0) cycle + if (ib_neighbor_ranks(sx, sy, sz) < 0) cycle nreqs = nreqs + 1 - call MPI_IRECV(recv_tables(:,:,:,nbr_idx), (2*ib_neighborhood_radius + 1)**3, MPI_INTEGER, & - & temporary_neighbors(sx, sy, sz), 400, MPI_COMM_WORLD, requests(nreqs), ierr) + call MPI_IRECV(recv_tables(:,:,:,nbr_idx), (2*ax + 1)**3, MPI_INTEGER, ib_neighbor_ranks(sx, sy, sz), & + & 400, MPI_COMM_WORLD, requests(nreqs), ierr) end do end do end do @@ -1448,10 +1424,10 @@ contains do sy = -1, 1 do sx = -1, 1 if (sx == 0 .and. sy == 0 .and. sz == 0) cycle - if (temporary_neighbors(sx, sy, sz) < 0) cycle + if (ib_neighbor_ranks(sx, sy, sz) < 0) cycle nreqs = nreqs + 1 - call MPI_ISEND(send_table, (2*ib_neighborhood_radius + 1)**3, MPI_INTEGER, temporary_neighbors(sx, & - & sy, sz), 400, MPI_COMM_WORLD, requests(nreqs), ierr) + call MPI_ISEND(send_table, (2*ax + 1)**3, MPI_INTEGER, ib_neighbor_ranks(sx, sy, sz), 400, & + & MPI_COMM_WORLD, requests(nreqs), ierr) end do end do end do @@ -1464,16 +1440,16 @@ contains do sx = -1, 1 if (sx == 0 .and. sy == 0 .and. sz == 0) cycle nbr_idx = nbr_idx + 1 - if (temporary_neighbors(sx, sy, sz) < 0) cycle - do dz = -ib_neighborhood_radius, ib_neighborhood_radius - do dy = -ib_neighborhood_radius, ib_neighborhood_radius - do dx = -ib_neighborhood_radius, ib_neighborhood_radius + if (ib_neighbor_ranks(sx, sy, sz) < 0) cycle + do dz = -ax, ax + do dy = -ax, ax + do dx = -ax, ax if (recv_tables(dx, dy, dz, nbr_idx) == MPI_PROC_NULL) cycle - if (dx + sx < -ib_neighborhood_radius .or. dx + sx > ib_neighborhood_radius) cycle - if (dy + sy < -ib_neighborhood_radius .or. dy + sy > ib_neighborhood_radius) cycle - if (dz + sz < -ib_neighborhood_radius .or. dz + sz > ib_neighborhood_radius) cycle - if (temporary_neighbors(dx + sx, dy + sy, dz + sz) /= MPI_PROC_NULL) cycle - temporary_neighbors(dx + sx, dy + sy, dz + sz) = recv_tables(dx, dy, dz, nbr_idx) + if (dx + sx < -ax .or. dx + sx > ax) cycle + if (dy + sy < -ax .or. dy + sy > ax) cycle + if (dz + sz < -ax .or. dz + sz > ax) cycle + if (ib_neighbor_ranks(dx + sx, dy + sy, dz + sz) /= MPI_PROC_NULL) cycle + ib_neighbor_ranks(dx + sx, dy + sy, dz + sz) = recv_tables(dx, dy, dz, nbr_idx) end do end do end do @@ -1484,21 +1460,15 @@ contains deallocate (send_table, recv_tables) end if - - if (allocated(ib_neighbor_ranks)) deallocate (ib_neighbor_ranks) - allocate (ib_neighbor_ranks(-ib_neighborhood_radius:ib_neighborhood_radius,-ib_neighborhood_radius:ib_neighborhood_radius, & - & -ib_neighborhood_radius:ib_neighborhood_radius)) - ib_neighbor_ranks = temporary_neighbors(-ib_neighborhood_radius, ib_neighborhood_radius, -ib_neighborhood_radius, & - & ib_neighborhood_radius, -ib_neighborhood_radius, ib_neighborhood_radius) - deallocate (temporary_neighbors) #endif end subroutine s_compute_ib_neighbor_ranks - subroutine get_neighbor_bounds() + subroutine s_get_neighbor_bounds() - real(wp) :: beg_val, end_val, recv_val - integer :: k, send_neighbor, recv_neighbor, ierr + real(wp) :: beg_val, end_val, recv_val, bound, max_ib_bound, local_rank_width, max_rank_width + integer :: k, send_neighbor, recv_neighbor, ierr, temporary_radius + logical :: automatic_neighborhood_radius ! Default: unbounded in all directions (covers single-rank and no-MPI cases) @@ -1510,6 +1480,36 @@ contains neighbor_domain_z%end = huge(0._wp) #ifdef MFC_MPI + ! perform setup if we are doing automatic radius checking + if (ib_neighborhood_radius < 1) then + ib_neighborhood_radius = 0 ! ensure we are starting with 0 neighborhood radius + + ! set up a temporary radius array to hold the indices + automatic_neighborhood_radius = .true. + + ! determine the maximum length of space that needs to be contained by the neighborhood + max_ib_bound = -1._wp + do i = 1, num_ibs + call s_get_ib_bound(patch_ib(i), bound) + max_ib_bound = max(max_ib_bound, bound) + end do + do i = 1, num_particle_clouds + max_ib_bound = max(max_ib_bound, particle_cloud(i)%radius) + end do + + ! determine the upper bound on the size + local_rank_width = -1._wp + #:for X, ID, DIM [('x', 1, 'm'), ('y', 2, 'n'), ('z', 3, 'p')] + local_rank_width = max(local_rank_width, abs(${X}$_cb(${DIM}$) - ${X}$_cb(-1))) + #:endfor + call s_mpi_allreduce_max(local_rank_width, max_rank_width) + + ! approximate the size of the neighborhood with a local 1.1x fudge factor for safety + ib_neighborhood_radius = floor(0.5_wp*max_rank_width/(1.1_wp*max_ib_bound)) + else + automatic_neighborhood_radius = .false. + end if + ! For each direction, propagate the left/right boundary edges outward ib_neighborhood_radius hops. After k rounds: beg_val = ! left edge of the rank k hops to the left; end_val = right edge of the rank k hops to the right. #:for X, ID, TAG, DIM in [('x', 1, 100, 'm'), ('y', 2, 102, 'n'), ('z', 3, 104, 'p')] @@ -1544,6 +1544,6 @@ contains #:endfor #endif - end subroutine get_neighbor_bounds + end subroutine s_get_neighbor_bounds end module m_start_up From e448caab7cb5e7dce9df6578d2d57267c187350b Mon Sep 17 00:00:00 2001 From: danieljvickers Date: Tue, 18 Aug 2026 14:48:14 -0400 Subject: [PATCH 3/6] Compiled NVHPC GPU --- src/simulation/m_ib_patches.fpp | 2 +- src/simulation/m_start_up.fpp | 22 ++++++++-------------- 2 files changed, 9 insertions(+), 15 deletions(-) diff --git a/src/simulation/m_ib_patches.fpp b/src/simulation/m_ib_patches.fpp index e271284867..917b19ed28 100644 --- a/src/simulation/m_ib_patches.fpp +++ b/src/simulation/m_ib_patches.fpp @@ -23,7 +23,7 @@ module m_ib_patches implicit none private; public :: s_apply_ib_patches, s_update_ib_rotation_matrix, s_instantiate_STL_models, s_decode_patch_periodicity, & - & s_encode_patch_periodicity, s_initialize_ib_airfoils, s_get_periodicities + & s_encode_patch_periodicity, s_initialize_ib_airfoils, s_get_periodicities, s_get_ib_bound contains diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index 2495cd5411..2c86e1d59e 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -1468,7 +1468,6 @@ contains real(wp) :: beg_val, end_val, recv_val, bound, max_ib_bound, local_rank_width, max_rank_width integer :: k, send_neighbor, recv_neighbor, ierr, temporary_radius - logical :: automatic_neighborhood_radius ! Default: unbounded in all directions (covers single-rank and no-MPI cases) @@ -1484,30 +1483,25 @@ contains if (ib_neighborhood_radius < 1) then ib_neighborhood_radius = 0 ! ensure we are starting with 0 neighborhood radius - ! set up a temporary radius array to hold the indices - automatic_neighborhood_radius = .true. - ! determine the maximum length of space that needs to be contained by the neighborhood max_ib_bound = -1._wp - do i = 1, num_ibs - call s_get_ib_bound(patch_ib(i), bound) + do k = 1, num_ibs + call s_get_ib_bound(patch_ib(k), bound) max_ib_bound = max(max_ib_bound, bound) end do - do i = 1, num_particle_clouds - max_ib_bound = max(max_ib_bound, particle_cloud(i)%radius) + do k = 1, num_particle_clouds + max_ib_bound = max(max_ib_bound, particle_cloud(k)%radius) end do ! determine the upper bound on the size local_rank_width = -1._wp - #:for X, ID, DIM [('x', 1, 'm'), ('y', 2, 'n'), ('z', 3, 'p')] - local_rank_width = max(local_rank_width, abs(${X}$_cb(${DIM}$) - ${X}$_cb(-1))) + #:for X, ID, DIM in [('x', 1, 'm'), ('y', 2, 'n'), ('z', 3, 'p')] + if (num_dims >= ${ID}$) local_rank_width = max(local_rank_width, abs(${X}$_cb(${DIM}$) - ${X}$_cb(-1))) #:endfor call s_mpi_allreduce_max(local_rank_width, max_rank_width) - ! approximate the size of the neighborhood with a local 1.1x fudge factor for safety - ib_neighborhood_radius = floor(0.5_wp*max_rank_width/(1.1_wp*max_ib_bound)) - else - automatic_neighborhood_radius = .false. + ! approximate the size of the neighborhood with a local 1.1x fudge factor for safety, lower bound of 1 + ib_neighborhood_radius = max(1, floor(1.1_wp*max_ib_bound/(0.5_wp*max_rank_width))) end if ! For each direction, propagate the left/right boundary edges outward ib_neighborhood_radius hops. After k rounds: beg_val = From bc3e9ac52ffe7f9f65458fcb67ab5779d5e927dd Mon Sep 17 00:00:00 2001 From: danieljvickers Date: Wed, 19 Aug 2026 12:04:26 -0400 Subject: [PATCH 4/6] Fixed STL issues --- src/simulation/m_ib_patches.fpp | 4 ++-- src/simulation/m_start_up.fpp | 3 ++- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/src/simulation/m_ib_patches.fpp b/src/simulation/m_ib_patches.fpp index 917b19ed28..0a845c8fe1 100644 --- a/src/simulation/m_ib_patches.fpp +++ b/src/simulation/m_ib_patches.fpp @@ -481,7 +481,7 @@ contains ly(1) = stl_bounding_boxes(patch%model_id, 2, 1) ly(2) = stl_bounding_boxes(patch%model_id, 2, 3) - bound = sqrt((lx(2) - lx(1))**2 + (ly(2) - ly(1))**2) + bound = 0.5_wp*sqrt((lx(2) - lx(1))**2 + (ly(2) - ly(1))**2) else if (patch%geometry == 6) then ! ellipse geometry bound = 0.5_wp*max(patch%length_x, patch%length_y) @@ -500,7 +500,7 @@ contains lz(1) = stl_bounding_boxes(patch%model_id, 3, 1) + patch%centroid_offset(3) lz(2) = stl_bounding_boxes(patch%model_id, 3, 3) + patch%centroid_offset(3) - bound = sqrt((lx(2) - lx(1))**2 + (ly(2) - ly(1))**2 + (lz(2) - lz(1))**2) + bound = 0.5_wp*sqrt((lx(2) - lx(1))**2 + (ly(2) - ly(1))**2 + (lz(2) - lz(1))**2) end if end subroutine s_get_ib_bound diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index 3b6a194ba3..7e33d2690e 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -898,6 +898,7 @@ contains type(ib_patch_parameters), allocatable :: particle_cloud_ibs(:) integer :: num_particle_cloud_ibs + call s_instantiate_STL_models() call s_get_neighbor_bounds() if (cfl_dt .and. n_start > 0) then @@ -911,7 +912,6 @@ contains else 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, num_particle_cloud_ibs) deallocate (particle_cloud_ibs) @@ -1530,6 +1530,7 @@ contains ! approximate the size of the neighborhood with a local 1.1x fudge factor for safety, lower bound of 1 ib_neighborhood_radius = max(1, floor(1.1_wp*max_ib_bound/(0.5_wp*max_rank_width))) + if (proc_rank == 0) print *, "Automatic choice of ib_neighborhood_radius selected: ", ib_neighborhood_radius, end if ! For each direction, propagate the left/right boundary edges outward ib_neighborhood_radius hops. After k rounds: beg_val = From 255096685230df7c7fd00ecf10bbab609130e5ce Mon Sep 17 00:00:00 2001 From: danieljvickers Date: Wed, 19 Aug 2026 14:10:23 -0400 Subject: [PATCH 5/6] Fixed trailing comma --- src/simulation/m_start_up.fpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index 7e33d2690e..63cebda913 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -1530,7 +1530,7 @@ contains ! approximate the size of the neighborhood with a local 1.1x fudge factor for safety, lower bound of 1 ib_neighborhood_radius = max(1, floor(1.1_wp*max_ib_bound/(0.5_wp*max_rank_width))) - if (proc_rank == 0) print *, "Automatic choice of ib_neighborhood_radius selected: ", ib_neighborhood_radius, + if (proc_rank == 0) print *, "Automatic choice of ib_neighborhood_radius selected: ", ib_neighborhood_radius end if ! For each direction, propagate the left/right boundary edges outward ib_neighborhood_radius hops. After k rounds: beg_val = From dd2fdf7290c29c170bcc8098f30d67d7bc7aaae8 Mon Sep 17 00:00:00 2001 From: danieljvickers Date: Wed, 19 Aug 2026 19:46:24 -0400 Subject: [PATCH 6/6] Addressed AI comments --- src/simulation/m_start_up.fpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index 63cebda913..7a55635eb2 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -899,6 +899,7 @@ contains integer :: num_particle_cloud_ibs call s_instantiate_STL_models() + call s_initialize_ib_airfoils() call s_get_neighbor_bounds() if (cfl_dt .and. n_start > 0) then @@ -912,7 +913,6 @@ contains else call s_generate_particle_clouds(particle_cloud_ibs, num_particle_cloud_ibs) end if - call s_initialize_ib_airfoils() call s_reduce_ib_patch_array(particle_cloud_ibs, num_particle_cloud_ibs) deallocate (particle_cloud_ibs) end block @@ -1494,7 +1494,7 @@ contains subroutine s_get_neighbor_bounds() - real(wp) :: beg_val, end_val, recv_val, bound, max_ib_bound, local_rank_width, max_rank_width + real(wp) :: beg_val, end_val, recv_val, bound, max_ib_bound, local_rank_width, min_rank_width integer :: k, send_neighbor, recv_neighbor, ierr, temporary_radius ! Default: unbounded in all directions (covers single-rank and no-MPI cases) @@ -1526,10 +1526,10 @@ contains #:for X, ID, DIM in [('x', 1, 'm'), ('y', 2, 'n'), ('z', 3, 'p')] if (num_dims >= ${ID}$) local_rank_width = max(local_rank_width, abs(${X}$_cb(${DIM}$) - ${X}$_cb(-1))) #:endfor - call s_mpi_allreduce_max(local_rank_width, max_rank_width) + call s_mpi_allreduce_min(local_rank_width, min_rank_width) ! approximate the size of the neighborhood with a local 1.1x fudge factor for safety, lower bound of 1 - ib_neighborhood_radius = max(1, floor(1.1_wp*max_ib_bound/(0.5_wp*max_rank_width))) + ib_neighborhood_radius = max(1, ceiling(1.1_wp*max_ib_bound/(min_rank_width))) if (proc_rank == 0) print *, "Automatic choice of ib_neighborhood_radius selected: ", ib_neighborhood_radius end if