From 04a33246a1464f8f403259cbd0eae65cff5b2e72 Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 30 Oct 2026 10:59:55 +0000 Subject: [PATCH] =?UTF-8?q?research:=20D=C3=BCtsch=20driver=20=E2=80=94=20?= =?UTF-8?q?exact-point=20body=20mask,=20periodic=20geometry=20from=20step?= =?UTF-8?q?=20mod=20T?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-Authored-By: Claude Opus 5.5 --- research/duetsch_re100_kc5/protocol.py | 21 ++++++++++++++++++--- research/duetsch_re100_kc5/run.py | 5 +++-- tests/unit/test_duetsch_protocol.py | 22 ++++++++++++++++++++++ 3 files changed, 43 insertions(+), 5 deletions(-) diff --git a/research/duetsch_re100_kc5/protocol.py b/research/duetsch_re100_kc5/protocol.py index df65e35..c070535 100644 --- a/research/duetsch_re100_kc5/protocol.py +++ b/research/duetsch_re100_kc5/protocol.py @@ -25,12 +25,15 @@ def amplitude(D): def position(s, T, A): - """Centre position relative to the mean: x = A cos(2 pi s / T).""" - return A * math.cos(2.0 * math.pi * s / T) + """Centre position relative to the mean: x = A cos(2 pi s / T), evaluated from the + native step modulo T so every cycle repeats the same geometry bit for bit.""" + return A * math.cos(2.0 * math.pi * (int(s) % T) / T) def phase(s, T): - return (360.0 * s / T - 90.0) % 360.0 + """Paper phase in degrees from the native step modulo T (the unwrapped cycle + index is kept separately by ``cycle``).""" + return (360.0 * (int(s) % T) / T - 90.0) % 360.0 def cycle(s, T): @@ -63,6 +66,18 @@ def morison_fit(f_norm, b1, b2): return float(c_d), float(c_i) +def inside_circle(x, y, centre, radius): + return (x - centre[0]) ** 2 + (y - centre[1]) ** 2 < radius * radius + + +def sample(field, fluid, x, y, centre, radius): + """Profile sample: MISSING if the exact point is inside the body (checked before + any corner interpolation), else ``bilinear_fluid``.""" + if inside_circle(x, y, centre, radius): + return None, (), 0.0 + return bilinear_fluid(field, fluid, x, y) + + def bilinear_fluid(field, fluid, x, y): """Bilinear value at (x, y) in cell coordinates from fluid cells only. diff --git a/research/duetsch_re100_kc5/run.py b/research/duetsch_re100_kc5/run.py index f50cce9..72017d3 100644 --- a/research/duetsch_re100_kc5/run.py +++ b/research/duetsch_re100_kc5/run.py @@ -70,11 +70,12 @@ def main(): fields = sim.get_macroscopic() fluid = (sim.get_flags().reshape(n, n) & 1) != 0 ux, uy = (np.asarray(fields[k]).reshape(n, n) for k in ("ux", "uy")) + centre = (xm + p.position(s, T, A), ym) # analytic centre after the step for xs in p.STATIONS: for x2 in p.X2: x, y = xm + xs * D, ym + x2 * D - u, wu, wsum = p.bilinear_fluid(ux, fluid, x, y) - v, _, _ = p.bilinear_fluid(uy, fluid, x, y) + u, wu, wsum = p.sample(ux, fluid, x, y, centre, R) + v, _, _ = p.sample(uy, fluid, x, y, centre, R) samples.append((s, p.cycle(s, T), ph, xs, x2, math.nan if u is None else u, math.nan if v is None else v, wsum, json.dumps(wu))) diff --git a/tests/unit/test_duetsch_protocol.py b/tests/unit/test_duetsch_protocol.py index 2de5e34..0751be7 100644 --- a/tests/unit/test_duetsch_protocol.py +++ b/tests/unit/test_duetsch_protocol.py @@ -62,3 +62,25 @@ def test_interval_bases_and_morison_fit_recover_known_coefficients(): assert np.allclose(b[:, 1], np.sin(phi), atol=1e-5) # midpoint ~ interval mean f = 2.09 * b[:, 0] - (math.pi ** 2 * 1.45 / p.KC) * b[:, 1] assert p.morison_fit(f, b[:, 0], b[:, 1]) == pytest.approx((2.09, 1.45), abs=1e-9) + + +def test_geometry_repeats_exactly_across_cycles(): + A = p.amplitude(20) + for m in (0, 1, T // 4, T // 2 - 1, 7 * T // 12): + for c in (1, 5, 11): + assert p.position(m + c * T, T, A) == p.position(m, T, A) + assert p.phase(m + c * T, T) == p.phase(m, T) + v0 = p.position(m + 1, T, A) - p.position(m, T, A) + assert p.position(m + 1 + c * T, T, A) - p.position(m + c * T, T, A) == v0 + assert p.cycle(T // 4 + 5 * T, T) == 6 # unwrapped index stays separate + + +def test_exact_point_inside_body_is_missing_even_with_fluid_corners(): + field = np.zeros((6, 6)) + fluid = np.ones((6, 6), bool) + fluid[3, 3] = False # only one corner solid + centre, radius = (3.0, 3.0), 0.3 + value, used, total = p.sample(field, fluid, 3.2, 3.1, centre, radius) # inside the circle + assert value is None and used == () and total == 0.0 + value, used, total = p.sample(field, fluid, 3.35, 3.0, centre, radius) # just outside + assert value is not None and total > 0.0