Give the finite-difference coefficients the buffer range the IB force asks for

examples/2D_ibm_force_decomposition/README.md

@@ -0,0 +1,43 @@
+# Immersed-boundary force under a changing MPI decomposition
+
+A static cylinder in uniform flow at Re 40, Ma 0.1, on a uniform 200 x 200 grid. The body sits at the
+origin, which is also the center of the domain, so any decomposition with an even number of ranks in a
+direction puts a subdomain edge straight through it.
+
+The force on a rigid body is a property of the flow. Running the same case on a different number of ranks
+must not change it. This case measures whether it does.
+
+## Running it
+
+```
+./mfc.sh run examples/2D_ibm_force_decomposition/case.py -n 1
+./mfc.sh run examples/2D_ibm_force_decomposition/case.py -n 4
+```
+
+Each run writes `restart_data/ib_state_200.dat`: 20 doubles, `[time, Fx, Fy, Fz, Tx, Ty, Tz, ...]`.
+
+```
+python3 -c "import numpy as np; a = np.fromfile('restart_data/ib_state_200.dat'); print(a[1], a[2])"
+```
+
+## What it shows
+
+`fd_order = 4` here, so the force integral's stencil reaches two cells, and a body cell two cells from a
+subdomain edge asks for finite-difference coefficients outside the interior. Before the coefficients were
+defined over that range they were read off the end of the array:
+
+| ranks | Fx before | Fx after |
+| --- | --- | --- |
+| 1 | 1.02352019 | 1.02352019 |
+| 4 | 1.02933438 | 1.02351629 |
+
+0.57 percent of drag appearing out of adjacent memory, against 0.0004 percent after. The single-rank
+answer is unchanged, because a single rank never reaches outside its own interior.
+
+The size of the discrepancy is set by whatever is adjacent in memory and is not bounded by anything: on a
+3D sphere at Re 100 on 64 ranks the same read produced a transverse force of 1.08 times the drag on a body
+that has none.
+
+Lift is a second, independent check: the cylinder is symmetric about y = 0, so Fy must be zero. It is
+about 8e-6 here in every configuration, which is the level at which the staircase representation of a
+circle on this grid is symmetric, and it does not move.

examples/2D_ibm_force_decomposition/case.py

@@ -0,0 +1,93 @@
+#!/usr/bin/env python3
+"""The smallest case that shows MFC's immersed-boundary force depending on the MPI decomposition.
+
+A static cylinder in uniform flow, Re 40, on a grid coarse enough to run on a laptop. The body sits at the
+origin, which is also the domain center, so a decomposition with an even number of ranks in a direction puts
+a subdomain edge straight through the body. Run it at one rank and at four and compare `restart_data/ib_state_200.dat`:
+the cylinder is symmetric about y = 0 and the lift must be zero, and the drag must not care how the domain
+was cut up.
+"""
+
+import json
+import os
+
+Re, Ma, d, gamma, U, rho = 40.0, 0.1, 1.0, 1.4, 1.0, 1.0
+P = rho * U**2 / (gamma * Ma**2)
+L = float(os.environ.get("L", 5.0))
+dx = float(os.environ.get("DX", 0.1))
+N = int(2 * L / dx)
+NSTEP = int(os.environ.get("NSTEP", 200))
+
+case = {
+    "run_time_info": "T",
+    "parallel_io": "T",
+    "prim_vars_wrt": "T",
+    "ib_state_wrt": "T",
+    "format": "silo",
+    "precision": "double",
+    "x_domain%beg": -L,
+    "x_domain%end": L,
+    "y_domain%beg": -L,
+    "y_domain%end": L,
+    "m": N - 1,
+    "n": N - 1,
+    "p": 0,
+    "cyl_coord": "F",
+    "dt": 0.3 * dx / (U + U / Ma),
+    "t_step_start": 0,
+    "t_step_stop": NSTEP,
+    "t_step_save": NSTEP,
+    "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": int(os.environ.get("FDORDER", 4)),
+    "patch_icpp(1)%geometry": 3,
+    "patch_icpp(1)%x_centroid": 0.0,
+    "patch_icpp(1)%y_centroid": 0.0,
+    "patch_icpp(1)%length_x": 2 * L,
+    "patch_icpp(1)%length_y": 2 * L,
+    "patch_icpp(1)%vel(1)": U,
+    "patch_icpp(1)%vel(2)": 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 / (d * 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%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,
+    "ib": "T",
+    "num_ibs": 1,
+    "ib_neighborhood_radius": 3,
+    "patch_ib(1)%geometry": 2,
+    "patch_ib(1)%x_centroid": 0.0,
+    "patch_ib(1)%y_centroid": 0.0,
+    "patch_ib(1)%radius": d / 2,
+    "patch_ib(1)%slip": "F",
+    "patch_ib(1)%moving_ibm": 0,
+}
+print(json.dumps(case, indent=4))

src/simulation/m_derived_variables.fpp

@@ -38,14 +38,16 @@ contains
         ! higher than fourth-order accuracy coefficients are wanted, the formulae required to compute these coefficients will have
         ! to be implemented in the subroutine s_compute_finite_difference_coefficients.
 
-        ! Allocating centered finite-difference coefficients
+        ! Allocating centered finite-difference coefficients. The IB force integral evaluates the viscous stress
+        ! at every cell its own stencil reaches, so it asks for coefficients up to fd_number outside the interior
+        ! whenever a body touches a subdomain edge; that range is included here rather than read off the end.
         if (probe_wrt .or. ib) then
-            @:ALLOCATE(fd_coeff_x(-fd_number:fd_number, 0:m))
+            @:ALLOCATE(fd_coeff_x(-fd_number:fd_number, -fd_number:m + fd_number))
             if (n > 0) then
-                @:ALLOCATE(fd_coeff_y(-fd_number:fd_number, 0:n))
+                @:ALLOCATE(fd_coeff_y(-fd_number:fd_number, -fd_number:n + fd_number))
             end if
             if (p > 0) then
-                @:ALLOCATE(fd_coeff_z(-fd_number:fd_number, 0:p))
+                @:ALLOCATE(fd_coeff_z(-fd_number:fd_number, -fd_number:p + fd_number))
             end if
 
             @:ALLOCATE(accel_mag(0:m, 0:n, 0:p))
@@ -63,22 +65,25 @@ contains
     !> Allocate and open derived variables. Computing FD coefficients.
     impure subroutine s_initialize_derived_variables
 
+        type(int_bounds_info) :: fd_offset
+
         if (probe_wrt .or. ib) then
             ! Opening and writing header of flow probe files
             if (proc_rank == 0 .and. probe_wrt) then
                 call s_open_probe_files()
                 call s_open_com_files()
             end if
             ! Computing centered finite difference coefficients
-            call s_compute_finite_difference_coefficients(m, x_cc, fd_coeff_x, buff_size, fd_number, fd_order)
+            fd_offset%beg = fd_number; fd_offset%end = fd_number
+            call s_compute_finite_difference_coefficients(m, x_cc, fd_coeff_x, buff_size, fd_number, fd_order, fd_offset)
             $:GPU_UPDATE(device='[fd_coeff_x]')
 
             if (n > 0) then
-                call s_compute_finite_difference_coefficients(n, y_cc, fd_coeff_y, buff_size, fd_number, fd_order)
+                call s_compute_finite_difference_coefficients(n, y_cc, fd_coeff_y, buff_size, fd_number, fd_order, fd_offset)
                 $:GPU_UPDATE(device='[fd_coeff_y]')
             end if
             if (p > 0) then
-                call s_compute_finite_difference_coefficients(p, z_cc, fd_coeff_z, buff_size, fd_number, fd_order)
+                call s_compute_finite_difference_coefficients(p, z_cc, fd_coeff_z, buff_size, fd_number, fd_order, fd_offset)
                 $:GPU_UPDATE(device='[fd_coeff_z]')
             end if
         end if