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.

Author: sbryngelson

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, <this rank's extent in that direction>)
+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.

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))

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)