From 1f9af078a03e3cc81fd3f7d5481cd6f7e9fdbe1f Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Sat, 12 Sep 2026 14:58:07 -0400 Subject: [PATCH 1/2] Add the thin-plate force case that measures issue 1849 A 2D plate pitching about its leading edge against Jantzen et al. (2014) Fig. 10 curve C1, with the number of cells across the 2.5 percent thickness set by NCELL so the physical problem is fixed and only the resolution moves. This measures the defect rather than fixing it. The reading that matters is that the disagreement with the reference does not shrink with refinement: 20.9, 33.5, 27.9 and 26.9 percent rms at 2, 4, 8 and 16 cells across. That rules out under-resolution and points at the force computation. A second reading falls out of the same sweep: the peak converges from four cells up (4.682, 4.458, 4.833 across a fourfold refinement), while two cells is 35 percent out. A thin body does not need ten or more cells for the immersed boundary to resolve it, which matters for cost estimates. --- examples/2D_ibm_thin_plate_force/README.md | 60 +++++++++ examples/2D_ibm_thin_plate_force/case.py | 138 +++++++++++++++++++++ 2 files changed, 198 insertions(+) create mode 100644 examples/2D_ibm_thin_plate_force/README.md create mode 100644 examples/2D_ibm_thin_plate_force/case.py diff --git a/examples/2D_ibm_thin_plate_force/README.md b/examples/2D_ibm_thin_plate_force/README.md new file mode 100644 index 000000000..9d6b0e0c7 --- /dev/null +++ b/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 diff --git a/examples/2D_ibm_thin_plate_force/case.py b/examples/2D_ibm_thin_plate_force/case.py new file mode 100644 index 000000000..b984ba5de --- /dev/null +++ b/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)) From 617b7f192fa7be4f6e373d19dff2a81801a5af11 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Sat, 12 Sep 2026 23:04:08 -0500 Subject: [PATCH 2/2] Skip the thin-plate force example in the test suite; it is a resolution study --- 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..5d48275aa 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", + # 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