Add the thin-plate force case that measures issue 1849

examples/2D_ibm_thin_plate_force/README.md

@@ -0,0 +1,60 @@
+# Immersed-boundary force on a thin plate (issue #1849)
+
+A 2D flat plate pitching about its leading edge, 0 → 45° on an Eldredge smoothed ramp (a = 21), K = π/8
+(case C1), Re_c = 300, Ma 0.2, plate thickness 2.5 % of chord. There is a published measurement to compare
+against:
+
+> Jantzen, Taira, Granlund & Ol, *Phys. Fluids* **26**, 053606 (2014), Fig. 10, 2D panel, curve C1.
+
+`C_L = F_y / (½ ρ U² c) = 2 F_y` here, with ρ = U = c = 1 and MFC's 2D force being per unit depth.
+
+This case exists to **measure** the defect in #1849, not to fix it. Whoever does fix it needs a number to move.
+
+## Running the sweep
+
+`NCELL` sets how many cells lie across the plate thickness. The physical problem is identical at every level —
+same chord, same thickness, same domain, same times — so the sweep isolates the resolution requirement from
+any change of geometry, which a thickness sweep would confound.
+
+```
+NCELL=2  ./mfc.sh run examples/2D_ibm_thin_plate_force/case.py     # dx = 0.0125 c,   0.22 M cells
+NCELL=4  ./mfc.sh run examples/2D_ibm_thin_plate_force/case.py     # dx = 0.00625 c,  0.90 M cells
+NCELL=8  ./mfc.sh run examples/2D_ibm_thin_plate_force/case.py     # dx = 0.003125 c, 3.58 M cells
+NCELL=16 ./mfc.sh run examples/2D_ibm_thin_plate_force/case.py     # dx = 0.0015625 c, 14.3 M cells
+```
+
+`SUMMARY=1 python3 case.py` prints the grid and step count without running anything.
+
+## What the sweep shows
+
+| cells across the thickness | peak C_L | C_L at t = 4 | rms difference from the reference |
+| --- | --- | --- | --- |
+| 2 | 6.413 | 1.209 | 20.9 % |
+| 4 | 4.682 | 1.271 | 33.5 % |
+| 8 | 4.458 | 1.111 | 27.9 % |
+| 16 | 4.833 | 1.005 | 26.9 % |
+| reference | 6.998 | 1.521 | — |
+
+Two separate readings:
+
+**The peak converges from four cells up.** 4.682, 4.458 and 4.833 across a fourfold refinement — a spread of
+8 % with no trend. Two cells is genuinely under-resolved and 35 % out; four is enough. So a thin body does not
+need ten or more cells before the immersed boundary resolves it, which is worth knowing on its own for cost
+estimates.
+
+**The disagreement with the reference does not shrink.** The last column sits at 27–34 % at every resolution
+and shows no sign of falling as the grid refines. **That is the measurement that matters for #1849**: it rules
+out under-resolution as the explanation and points at the force computation itself.
+
+## Why a zero-reference check is not available here
+
+The obvious cheap test — a symmetric body whose true force is exactly zero — does not apply to a pitching
+plate, whose lift is large and unknown. For that style of check see
+`examples/2D_ibm_force_decomposition` (a cylinder, whose lift must be zero) and
+`examples/3D_ibm_neighborhood_radius`. This case trades the exact reference for a published one.
+
+## Related
+
+- #1849 — the force integral over body cells on a thin plate, which this measures
+- #1859 — an out-of-bounds coefficient read in the same force path, fixed
+- #1863 — a spurious transverse force on multi-rank runs, open

examples/2D_ibm_thin_plate_force/case.py

@@ -0,0 +1,138 @@
+#!/usr/bin/env python3
+"""Immersed-boundary force on a thin plate, against a published measurement (issue 1849).
+
+A 2D flat plate pitching about its leading edge, 0 -> 45 degrees on an Eldredge smoothed ramp (a = 21),
+K = pi/8 (case C1), Re_c = 300, Ma 0.2. The plate is 2.5 percent of the chord thick, matching the experiment.
+
+    Jantzen, Taira, Granlund & Ol, Phys. Fluids 26, 053606 (2014), Fig. 10, 2D panel, curve C1.
+
+`NCELL` sets how many cells lie across the plate thickness; the physical problem does not change with it, so
+the sweep isolates the immersed-boundary resolution requirement from any change of geometry. See README.md
+for what the sweep shows and why it matters.
+
+    NCELL=2 python3 case.py      dx = 0.0125 c
+    NCELL=4 python3 case.py      dx = 0.00625 c   (default)
+    NCELL=8 python3 case.py      dx = 0.003125 c
+    NCELL=16 python3 case.py     dx = 0.0015625 c
+"""
+
+import json
+import math
+import os
+
+U, rho, Ma, gamma, Re = 1.0, 1.0, 0.2, 1.4, 300.0
+P = rho * U**2 / (gamma * Ma**2)
+cs = math.sqrt(gamma * P / rho)
+K = math.pi / 8
+Omega = 2 * K * U  # rad per c/U; pitch time = 45 deg / Omega = 1 c/U
+th_max = math.radians(45.0)
+t_p = th_max / Omega
+a_smooth = 21.0
+THICK = 0.025
+t0 = 2.0  # settle at 0 deg before the ramp
+t_end = t0 + t_p + 4.0
+x0, x1, y0, y1 = -2.0, 5.0, -2.5, 2.5
+
+# L4 added after the first three failed to converge: refining 2 -> 4 -> 8 cells across the thickness moved
+# the peak lift 6.41 -> 4.68 -> 4.46 against a reference of 7.00, i.e. away from it and then stalling.
+# Two under-resolved answers landing near each other is not convergence. If a few cells across a thin
+# body is simply too few for the immersed boundary, 16 should move back toward the reference; if the
+# finite thickness is genuinely the difference, it should stay near 4.5.
+LEVELS = {"L1": 0.0125, "L2": 0.00625, "L3": 0.003125, "L4": 0.0015625}
+
+
+NCELL = int(os.environ.get("NCELL", 4))
+dx = THICK / NCELL
+m, n = int((x1 - x0) / dx) - 1, int((y1 - y0) / dx) - 1
+dt = 0.4 * dx / (U + cs)
+nt = int(t_end / dt)
+
+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,
+    "m": m,
+    "n": n,
+    "p": 0,
+    "cyl_coord": "F",
+    "dt": dt,
+    "t_step_start": 0,
+    "t_step_stop": nt,
+    "t_step_save": max(1, nt // 40),
+    "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": 3,
+    "patch_icpp(1)%x_centroid": 0.5 * (x0 + x1),
+    "patch_icpp(1)%y_centroid": 0.5 * (y0 + y1),
+    "patch_icpp(1)%length_x": x1 - x0,
+    "patch_icpp(1)%length_y": y1 - y0,
+    "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,
+    "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_x%grcbc_out": "T",
+    "bc_x%pres_out": P,
+    "bc_y%beg": -9,
+    "bc_y%end": -9,
+    "ib": "T",
+    "num_ibs": 1,
+    "patch_ib(1)%geometry": 3,
+    "patch_ib(1)%x_centroid": 0.5,
+    "patch_ib(1)%y_centroid": 0.0,
+    "patch_ib(1)%length_x": 1.0,
+    "patch_ib(1)%length_y": THICK,
+    "patch_ib(1)%slip": "F",
+    "patch_ib(1)%moving_ibm": 1,
+    "patch_ib(1)%angles(3)": 0.0,
+    "patch_ib(1)%angular_vel(3)": 0.0,
+    "omega_wrt(3)": "T",
+}
+tau = f"(t - {t0})"
+th = f"(0.5*{th_max}*(1.0 + (log(cosh({a_smooth}*{tau})) - log(cosh({a_smooth}*({tau} - {t_p}))))/{a_smooth * t_p}))"
+thd = f"(0.5*{Omega}*(tanh({a_smooth}*{tau}) - tanh({a_smooth}*({tau} - {t_p}))))"
+case["patch_ib(1)%angular_vel(3)"] = thd
+case["patch_ib(1)%vel(1)"] = f"-0.5*{thd}*sin({th})"
+case["patch_ib(1)%vel(2)"] = f"0.5*{thd}*cos({th})"
+
+if __name__ == "__main__":
+    if os.environ.get("SUMMARY"):
+        print(f"{NCELL} cells across the {THICK:g} c thickness: dx = {dx:g} c, " f"{m + 1} x {n + 1} = {(m + 1) * (n + 1) / 1e6:.2f} M cells, {nt} steps of dt = {dt:.2e}")
+    else:
+        print(json.dumps(case, indent=4))

toolchain/mfc/test/cases.py

@@ -3224,6 +3224,8 @@ def foreach_example():
                 # that no code change can perturb. The same kinematics are covered meaningfully by the
                 # "IBM -> Prescribed Kinematics -> Pitch Ramp" case below, whose plate is four cells thick.
                 "3D_ibm_pitchup_plate",
+                # A grid-resolution study; the Example suite's 25-cell cap removes the resolution it measures.
+                "2D_ibm_thin_plate_force",
             ]
             if path in casesToSkip:
                 continue