Read IB viscous FD coefficients at the nearest interior cell

examples/2D_ibm_force_decomposition/README.md

@@ -0,0 +1,45 @@
+# 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 "after" column was
+measured with the coefficients defined over the range that is read; clamping to the nearest interior cell
+gives the same answer on this uniform grid, where the two stencils coincide.)
+
+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_viscous.fpp

@@ -1127,24 +1127,33 @@ contains
         real(wp), dimension(1:3,1:3)                          :: velocity_gradient_tensor
         real(wp)                                              :: divergence
         real(wp)                                              :: mu_eff, gamma_dot_c
-        integer                                               :: l, q  !< iterators
+        integer                                               :: l, q              !< iterators
         integer                                               :: fl
         integer                                               :: r
+        integer                                               :: i_fd, j_fd, k_fd  !< sample clamped to the interior
 
         ! zero the viscous stress and collection of velocity derivatives
         viscous_stress_tensor = 0._wp
         velocity_gradient_tensor = 0._wp
 
+        ! fd_coeff_x/y/z are computed for interior cells only (0:m, 0:n, 0:p), but s_compute_ib_forces samples this routine
+        ! fd_number cells out from an interior cell, so a body near a domain boundary asks for a coefficient that was never
+        ! computed. Read the nearest interior cell's coefficients there: on a uniform grid they are the same, and on a
+        ! stretched one this is the stencil the boundary cell itself uses.
+        i_fd = min(max(i, 0), m)
+        j_fd = min(max(j, 0), n)
+        k_fd = min(max(k, 0), p)
+
         ! compute the velocity gradient tensor with the same fd_order-respecting stencil as the stress-divergence outer derivative
         do l = 1, num_dims
             do r = -fd_number, fd_number
                 velocity_gradient_tensor(l, 1) = velocity_gradient_tensor(l, 1) + fd_coeff_x(r, &
-                                         & i)*q_prim_vf(eqn_idx%mom%beg + l - 1)%sf(i + r, j, k)
+                                         & i_fd)*q_prim_vf(eqn_idx%mom%beg + l - 1)%sf(i + r, j, k)
                 velocity_gradient_tensor(l, 2) = velocity_gradient_tensor(l, 2) + fd_coeff_y(r, &
-                                         & j)*q_prim_vf(eqn_idx%mom%beg + l - 1)%sf(i, j + r, k)
+                                         & j_fd)*q_prim_vf(eqn_idx%mom%beg + l - 1)%sf(i, j + r, k)
                 if (num_dims == 3) then
                     velocity_gradient_tensor(l, 3) = velocity_gradient_tensor(l, 3) + fd_coeff_z(r, &
-                                             & k)*q_prim_vf(eqn_idx%mom%beg + l - 1)%sf(i, j, k + r)
+                                             & k_fd)*q_prim_vf(eqn_idx%mom%beg + l - 1)%sf(i, j, k + r)
                 end if
             end do
         end do