From 733f4aff3576906936ec26c9ec04ae7d87130577 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Sat, 12 Sep 2026 14:20:33 -0400 Subject: [PATCH 1/2] Choose ib_neighborhood_radius from the thinnest rank, not the widest ib_neighborhood_radius counts rank hops: a rank keeps an immersed-boundary patch only while the centroid lies inside its subdomain grown by that many hops, and drops it otherwise. The radius therefore has to span the body in every direction, and the hop that needs the most is the one crossing the thinnest rank. The automatic choice assembled its rank width the wrong way round -- each rank reported its widest extent and the minimum was taken across ranks. Where every rank is long in one direction and thin in another, the width reported is the long one and the radius comes out too small. examples/3D_ibm_neighborhood_radius is a thin plate in a long narrow channel, sized so the topology search settles on 16 x 2 x 2 at 64 ranks: x ranks 1.250 wide against y and z at 4.000 and 3.000. The plate's half-extent is 1.3010, so crossing it takes two hops of 1.250 and one of 4.000. MFC prints radius 1 before this change and 2 after; both runs complete. The two production grids that motivated it agree: a gust case at 16 x 2 x 4 moves 1 -> 2, and a flapping wing at 8 x 4 x 4 moves 1 -> 2. The new width is never larger than the old, so the radius never decreases -- the change can only make the neighbourhood more conservative. Cases that set the parameter explicitly are untouched, as is any near-cubic decomposition, where the widest and narrowest extents coincide. --- examples/3D_ibm_neighborhood_radius/README.md | 77 ++++++++++ examples/3D_ibm_neighborhood_radius/case.py | 132 ++++++++++++++++++ src/simulation/m_start_up.fpp | 10 +- 3 files changed, 216 insertions(+), 3 deletions(-) create mode 100644 examples/3D_ibm_neighborhood_radius/README.md create mode 100644 examples/3D_ibm_neighborhood_radius/case.py diff --git a/examples/3D_ibm_neighborhood_radius/README.md b/examples/3D_ibm_neighborhood_radius/README.md new file mode 100644 index 000000000..f635700ef --- /dev/null +++ b/examples/3D_ibm_neighborhood_radius/README.md @@ -0,0 +1,77 @@ +# Automatic `ib_neighborhood_radius` when ranks are not cubes + +`ib_neighborhood_radius` is a count of **rank hops**. A rank keeps an immersed-boundary patch only while the +patch's centroid lies inside its own subdomain grown outward by that many hops +(`s_get_neighbor_bounds`, `f_neighborhood_ranks_own_location`), and drops it otherwise. So the radius has to be +large enough that the body is reachable within that many hops **in every direction** — and the hop that needs +the most is the one stepping across the *thinnest* rank. + +When the radius is not set in the case file, MFC chooses it from the body's half-extent divided by a rank +width. The width it used was assembled the wrong way round: + +```fortran +local_rank_width = -1._wp +do each direction: + local_rank_width = max(local_rank_width, ) +call s_mpi_allreduce_min(local_rank_width, min_rank_width) +``` + +Each rank reports its **widest** extent, and the minimum is taken over ranks. On a decomposition where every +rank is long in one direction and thin in another, the reported width is the long one, the radius comes out too +small, and ranks that should have kept the patch drop it. + +## The case + +A thin plate in a long, narrow channel: 20 chords by 8 by 6, with 400 x 50 x 50 cells chosen so MFC's topology +search settles on **16 x 2 x 2** at 64 ranks. That is an ordinary shape for a wake, a jet or a channel — the +flow direction resolved far more finely than the cross-stream ones — and it makes the ranks strongly +anisotropic: + +| direction | ranks | extent per rank | +| --- | --- | --- | +| x | 16 | **1.250** | +| y | 2 | 4.000 | +| z | 2 | 3.000 | + +The plate is 1.0 x 2.4 x 0.1, so `s_get_ib_bound` (geometry 9, the cuboid's half-diagonal) returns **1.3010**. +Crossing that at 1.250 per hop needs `ceil(1.1 * 1.3010 / 1.250) = 2` hops. The old width of 4.000 gives +`ceil(1.1 * 1.3010 / 4.000) = 1`. + +## Running it + +``` +./mfc.sh run examples/3D_ibm_neighborhood_radius/case.py -n 64 +``` + +and read the line MFC prints at start-up: + +``` +Automatic choice of ib_neighborhood_radius selected: N +``` + +| | printed radius | +| --- | --- | +| before | **1** | +| after | **2** | + +Both runs complete; the case is 1 M cells and takes a few minutes on two CPU nodes. `SUMMARY=1 python3 +case.py` prints the half-extent, the rank extents and the arithmetic above without running anything. + +## It is not only this case + +The same two production grids that motivated the fix, measured from their own `lustre_*_cb.dat`: + +| case | topology | old width | old radius | new width | new radius | +| --- | --- | --- | --- | --- | --- | +| gust encounter, 128 ranks | 16 x 2 x 4 | 2.051 | 1 | 1.052 | **2** | +| flapping wing, 128 ranks | 8 x 4 x 4 | 1.745 | 1 | 1.027 | **2** | + +Both pick 1 where 2 is required. + +## Scope + +The new width is never larger than the old one, so the chosen radius never decreases: the change can only make +the neighbourhood more conservative, at the cost of more hops in the force reduction. Cases that set +`ib_neighborhood_radius` explicitly are untouched, and so is any decomposition whose ranks are close to cubic, +where the widest and narrowest extents coincide — which is why a uniform grid with a balanced topology shows no +difference. diff --git a/examples/3D_ibm_neighborhood_radius/case.py b/examples/3D_ibm_neighborhood_radius/case.py new file mode 100644 index 000000000..1669458c6 --- /dev/null +++ b/examples/3D_ibm_neighborhood_radius/case.py @@ -0,0 +1,132 @@ +#!/usr/bin/env python3 +"""Automatic ib_neighborhood_radius on a decomposition whose ranks are not cubes. + +A thin plate held in a long, narrow channel. The domain is 20 chords long and 8 by 6 across, and the cell +counts (400 x 50 x 50) are chosen so MFC's topology search settles on 16 x 2 x 2 at 64 ranks. That makes the +x ranks 1.25 chords wide while the y and z ranks are 4.0 and 3.0 -- an ordinary situation for a wake, a jet or +a channel, where the flow direction is resolved far more finely than the cross-stream ones. + +`ib_neighborhood_radius` is deliberately left unset, so MFC chooses it at start-up and prints + + Automatic choice of ib_neighborhood_radius selected: N + +The radius counts rank hops, and the body must be reachable within that many hops in *every* direction, so +the hop that matters is the one crossing the thinnest rank. The plate's half-extent is 1.301 chords and the +thinnest rank is 1.250 wide, so it needs 2 hops, not 1. + +See README.md for the measured before/after and the argument. +""" + +import json +import math +import os + +Re, Ma, gamma, U, rho = 1000.0, 0.1, 1.4, 1.0, 1.0 +cs = U / Ma +P = rho * cs**2 / gamma +c = 1.0 # chord +SPAN = float(os.environ.get("SPAN", 2.4)) # spanwise length of the plate +THICK = 0.1 * c + +M = int(os.environ.get("M", 400)) # cells in x; 400/16 ranks = 25, the stencil minimum +N = int(os.environ.get("N", 50)) # cells in y; 50/2 ranks = 25 +PZ = int(os.environ.get("PZ", 50)) # cells in z; 50/2 ranks = 25 +TSTOP = float(os.environ.get("TSTOP", 0.5)) + +x0, x1 = -10.0 * c, 10.0 * c +y0, y1 = -4.0 * c, 4.0 * c +z0, z1 = -3.0 * c, 3.0 * c + +case = { + "run_time_info": "T", + "parallel_io": "T", + "prim_vars_wrt": "T", + "ib_state_wrt": "T", + "format": "silo", + "precision": "double", + "x_domain%beg": x0, + "x_domain%end": x1, + "y_domain%beg": y0, + "y_domain%end": y1, + "z_domain%beg": z0, + "z_domain%end": z1, + "m": M - 1, + "n": N - 1, + "p": PZ - 1, + "cyl_coord": "F", + "cfl_adap_dt": "T", + "cfl_target": 0.4, + "n_start": 0, + "t_save": TSTOP, + "t_stop": TSTOP, + "num_patches": 1, + "num_fluids": 1, + "model_eqns": "5eq", + "alt_soundspeed": "F", + "mpp_lim": "F", + "mixture_err": "T", + "time_stepper": "rk3", + "weno_order": 5, + "weno_eps": 1.0e-10, + "weno_Re_flux": "T", + "weno_avg": "T", + "avg_state": "arithmetic", + "mapped_weno": "T", + "null_weights": "F", + "mp_weno": "F", + "riemann_solver": "hllc", + "low_Mach": 2, + "wave_speeds": "direct", + "viscous": "T", + "fd_order": 4, + "patch_icpp(1)%geometry": 9, + "patch_icpp(1)%x_centroid": 0.0, + "patch_icpp(1)%y_centroid": 0.0, + "patch_icpp(1)%z_centroid": 0.0, + "patch_icpp(1)%length_x": x1 - x0, + "patch_icpp(1)%length_y": y1 - y0, + "patch_icpp(1)%length_z": z1 - z0, + "patch_icpp(1)%vel(1)": U, + "patch_icpp(1)%vel(2)": 0.0, + "patch_icpp(1)%vel(3)": 0.0, + "patch_icpp(1)%pres": P, + "patch_icpp(1)%alpha_rho(1)": rho, + "patch_icpp(1)%alpha(1)": 1.0, + "fluid_pp(1)%gamma": 1.0 / (gamma - 1.0), + "fluid_pp(1)%eos": "ideal_gas", + "fluid_pp(1)%Re(1)": Re / (c * U), + "bc_x%beg": -7, + "bc_x%grcbc_in": "T", + "bc_x%vel_in(1)": U, + "bc_x%vel_in(2)": 0.0, + "bc_x%vel_in(3)": 0.0, + "bc_x%pres_in": P, + "bc_x%alpha_rho_in(1)": rho, + "bc_x%alpha_in(1)": 1.0, + "bc_x%end": -8, + "bc_y%beg": -8, + "bc_y%end": -8, + "bc_z%beg": -8, + "bc_z%end": -8, + # ib_neighborhood_radius deliberately NOT set: this case exists to exercise the automatic choice + "ib": "T", + "num_ibs": 1, + "patch_ib(1)%geometry": 9, # cuboid + "patch_ib(1)%x_centroid": 0.0, + "patch_ib(1)%y_centroid": 0.0, + "patch_ib(1)%z_centroid": 0.0, + "patch_ib(1)%length_x": c, + "patch_ib(1)%length_y": SPAN, + "patch_ib(1)%length_z": THICK, + "patch_ib(1)%slip": "F", + "patch_ib(1)%moving_ibm": 0, +} + +if __name__ == "__main__": + if os.environ.get("SUMMARY"): + bound = 0.5 * math.sqrt(c**2 + SPAN**2 + THICK**2) + print(f"plate half-extent (s_get_ib_bound, geometry 9): {bound:.4f}") + print(f"rank extents at 64 ranks (16 x 2 x 2): x {(x1 - x0) / 16:.3f}, " f"y {(y1 - y0) / 2:.3f}, z {(z1 - z0) / 2:.3f}") + print(f"hops needed across the thinnest rank: ceil(1.1 * {bound:.4f} / {(x1 - x0) / 16:.3f}) = " f"{max(1, math.ceil(1.1 * bound / ((x1 - x0) / 16)))}") + else: + print(json.dumps(case, indent=4)) diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index 6a5434f5b..f1faba9af 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -1542,10 +1542,14 @@ contains 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 + ! Narrowest rank extent, over every direction as well as every rank. The radius is a count of rank + ! hops, so the distance one hop covers is the extent of the rank it steps over, and the direction + ! needing the most hops to span the body is the one whose ranks are thinnest. Reducing over each + ! rank's widest extent first reports the wrong number whenever ranks are anisotropic, which is the + ! norm on a stretched grid or an elongated domain. + local_rank_width = huge(0._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))) + if (num_dims >= ${ID}$) local_rank_width = min(local_rank_width, abs(${X}$_cb(${DIM}$) - ${X}$_cb(-1))) #:endfor call s_mpi_allreduce_min(local_rank_width, min_rank_width) From 15460d3ef60277a2200d7b2b3e14cf7cc4d2a4d6 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Sat, 12 Sep 2026 23:00:09 -0500 Subject: [PATCH 2/2] Skip the neighborhood-radius example in the test suite; it needs its rank topology --- toolchain/mfc/test/cases.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/toolchain/mfc/test/cases.py b/toolchain/mfc/test/cases.py index fff3b4552..1639a8d84 100644 --- a/toolchain/mfc/test/cases.py +++ b/toolchain/mfc/test/cases.py @@ -3167,6 +3167,8 @@ def foreach_example(): # the transverse momentum drifts past the 1e-3 Example tolerance across compilers # (nvhpc passes; Intel and CCE disagree by ~2e-3 absolute). No single golden is portable. "2D_hybrid_slab", + # Needs its 16 x 2 x 2 rank topology; the Example suite runs it on one rank and a shrunken grid. + "3D_ibm_neighborhood_radius", ] if path in casesToSkip: continue