diff --git a/docs/documentation/case.md b/docs/documentation/case.md index 23463f15b..249e6c15a 100644 --- a/docs/documentation/case.md +++ b/docs/documentation/case.md @@ -414,7 +414,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 48bba55fb..c5af1e123 100644 --- a/src/simulation/m_global_parameters.fpp +++ b/src/simulation/m_global_parameters.fpp @@ -486,7 +486,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 35b97ea31..0a845c8fe 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 @@ -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 = 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) + 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 = 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 + + 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 de5b62206..7a55635eb 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -898,7 +898,9 @@ contains type(ib_patch_parameters), allocatable :: particle_cloud_ibs(:) integer :: num_particle_cloud_ibs - call get_neighbor_bounds() + call s_instantiate_STL_models() + call s_initialize_ib_airfoils() + call s_get_neighbor_bounds() if (cfl_dt .and. n_start > 0) then call s_read_ib_restart_data(n_start) @@ -911,8 +913,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) end block @@ -1492,10 +1492,10 @@ contains 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, min_rank_width + integer :: k, send_neighbor, recv_neighbor, ierr, temporary_radius ! Default: unbounded in all directions (covers single-rank and no-MPI cases) @@ -1507,6 +1507,32 @@ 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 + + ! determine the maximum length of space that needs to be contained by the neighborhood + max_ib_bound = -1._wp + 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 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 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_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, 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 + ! 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')] @@ -1541,6 +1567,6 @@ contains #:endfor #endif - end subroutine get_neighbor_bounds + end subroutine s_get_neighbor_bounds end module m_start_up diff --git a/toolchain/mfc/params/definitions.py b/toolchain/mfc/params/definitions.py index 971685f09..d04b533d7 100644 --- a/toolchain/mfc/params/definitions.py +++ b/toolchain/mfc/params/definitions.py @@ -376,7 +376,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}, "nb": {"min": 1},