From 7b81584458904809117466f46d06b66d53937b59 Mon Sep 17 00:00:00 2001 From: Frank14f <1515444314@qq.com> Date: Sat, 31 Oct 2026 15:51:39 +0800 Subject: [PATCH] Predeclare zero-damping pinball development cases --- docs/steady_pinball_fiv.md | 33 ++++++++++++++++++++++++++++++- research/steady_pinball/run.py | 30 ++++++++++++++++++---------- tests/unit/test_steady_pinball.py | 32 +++++++++++++++++++++++++++++- 3 files changed, 83 insertions(+), 12 deletions(-) diff --git a/docs/steady_pinball_fiv.md b/docs/steady_pinball_fiv.md index a2eb3a3..977fbda 100644 --- a/docs/steady_pinball_fiv.md +++ b/docs/steady_pinball_fiv.md @@ -5,7 +5,8 @@ rho=1, mean U=.01, nu=.004, D=20, Re_D=50. Centres are (600,255.5), (626,270.5), (626,240.5), radius10. No upstream cylinder. Passive probes stay at (800,295.5), (800,255.5), (800,215.5), radius5. Every omega is zero. Parabolic inlet, regularized NEQ damping1; y bounce-back planes .5/510.5; -NEQ-extrapolated outlet with existing backflow clamp/blending. Initial fluid is +NEQ-extrapolated outlet with backflow clamp (the configured blend is inactive in +this MRT route). Initial fluid is parabolic equilibrium. Fixed warmup is 4 NX/U=512000 steps. Geometry and fluid parameters follow DynamisLab 609116ec26f19766f5e9ce2b0df427d4e24faae8, @@ -71,3 +72,33 @@ This stage makes no claim of physical grid convergence or vibration suppression. Run from checkout with PYTHONPATH=src and the existing Python environment: `python research/steady_pinball/run.py --smoke --out ` then separately `--out ` and `--fixed --out `. + +## Zero-damping development follow-up (predeclared 2026-10-31) + +The original d0bcefa run and its zeta=.01 results remain historical evidence. +The current runner exposes three fixed `--case` choices, all with C=0 as requested: + +| Case | M / (rho 3 pi D^2/4) | fn | M | K | Purpose | +|---|---:|---:|---:|---:|---| +| Z0 | 10 | .0001 | 9424.77796 | .00372075320 | Remove damping only | +| M2 | 2 | .0001 | 1884.95559 | .000744150640 | Reduce mass and stiffness together at fixed fn | +| R2 | 2 | .00006 | 1884.95559 | .000267894231 | Change stiffness only relative to M2 | + +R2 deliberately places fn at the retained fixed-body lift peak, f=.00006 +(St=.12 on mean U). This is an informed development choice, not independent +physical validation or a promise of lock-in. All geometries, fluid settings, +512000-step fixed warmup, y0=v0=0, force time level and unfiltered RHS are unchanged. +No run result may select a replacement parameter or an automatic extension. + +Each case has a 600000-step cap, the same 300 D/U physical duration. This equals +60 natural periods for Z0/M2 and 36 for R2. Compare the common final200000 native +steps, and six100000-step blocks, retaining actual natural/response frequencies +and ongoing drift. Do not call those blocks ten R2 periods. No settled-VIV claim +follows from a larger displacement. Report all three outcomes, including quiet +responses or stopped cases, alongside the existing fixed and damped-free records. + +Unchanged stops: |y|>D, invalid move/speed, nonfinite values or no_donor. A 1800s +per-case wall cap and GPU headroom check bound resource use. Preserve partial rows; +do not add damping, clip trajectories, filter forces, or relax any stop. Each +case adds72MB of scalar history, no full fields. The separate Legacy/modern +baseline diagnosis does not change this FSI protocol. diff --git a/research/steady_pinball/run.py b/research/steady_pinball/run.py index 2af149a..67e66b2 100644 --- a/research/steady_pinball/run.py +++ b/research/steady_pinball/run.py @@ -1,6 +1,6 @@ """Uncontrolled paper-geometry steady Re_D=50, one transverse assembly DOF. -One fixed protocol; --smoke performs 32 warmup + 32 observation steps only. +Three predeclared zero-damping structures; --smoke uses only 32 + 32 steps. See docs/steady_pinball_fiv.md for historical solver differences and new structure. """ from __future__ import annotations @@ -18,6 +18,7 @@ D, U, NU = 20., .01, .004 CENTRES = ((600., 255.5), (626., 270.5), (626., 240.5)) PROBES = ((800., 295.5), (800., 255.5), (800., 215.5)) WARMUP, STEPS, PERIOD = 512000, 600000, 10000 +STRUCTURES = {'Z0': (10., 1.e-4), 'M2': (2., 1.e-4), 'R2': (2., 6.e-5)} COLUMNS = ('step', 'y', 'v', 'cme_F', 'cme_U', 'cme_L', 'wen_F', 'wen_U', 'wen_L', 'node_F', 'node_U', 'node_L', 'eq_only', 'no_donor', 'crossings') ROOT = Path(__file__).resolve().parents[2] @@ -29,10 +30,15 @@ def config(): class Structure: """Exact damped oscillator update for a force held over dt=1; no force filter.""" - def __init__(self): - self.mass = 10. * 3. * math.pi * D * D / 4. - self.omega = 2. * math.pi / PERIOD - self.zeta = .01 + def __init__(self, mass_ratio=10., frequency=1.e-4, zeta=0.): + if not (math.isfinite(mass_ratio) and mass_ratio > 0. and + math.isfinite(frequency) and frequency > 0. and + math.isfinite(zeta) and 0. <= zeta < 1.): + raise ValueError('positive finite mass ratio/frequency and 0 <= zeta < 1 required') + self.mass_ratio, self.frequency = mass_ratio, frequency + self.mass = mass_ratio * 3. * math.pi * D * D / 4. + self.omega = 2. * math.pi * frequency + self.zeta = zeta self.stiffness = self.mass * self.omega**2 self.damping = 2. * self.zeta * self.mass * self.omega self.y = self.v = 0. @@ -68,6 +74,8 @@ def main(): ap.add_argument('--out', required=True, type=Path) ap.add_argument('--fixed', action='store_true', help='same modern moving-path geometry held fixed') ap.add_argument('--smoke', action='store_true', help='32+32 steps, no physical result') + ap.add_argument('--case', choices=STRUCTURES, default='Z0', + help='predeclared structure; all cases have zero structural damping') args = ap.parse_args() out = args.out.resolve() if out == ROOT or ROOT in out.parents: @@ -76,7 +84,7 @@ def main(): cfg = config() (out/'config.json').write_text(json.dumps(cfg, indent=2)+'\n') source = subprocess.check_output(['git', '-C', str(ROOT), 'rev-parse', 'HEAD'], text=True).strip() - s = Structure() + s = Structure(*STRUCTURES[args.case], zeta=0.) warmup, count = (32, 32) if args.smoke else (WARMUP, STEPS) rows = np.lib.format.open_memmap(out/'steps.npy', mode='w+', dtype=np.float64, shape=(count, len(COLUMNS))) @@ -87,11 +95,13 @@ def main(): col = {k: names.index(k) for k in names} status, completed, f_last = 'running', 0, 0. started = time.monotonic() - report = dict(source=source, fixed=args.fixed, smoke=args.smoke, columns=COLUMNS, + report = dict(source=source, case=args.case, fixed=args.fixed, smoke=args.smoke, columns=COLUMNS, structure=dict(mass=s.mass, stiffness=s.stiffness, damping=s.damping, - zeta=s.zeta, fn=1./PERIOD, y0=0., v0=0.), + mass_ratio_three_disks=s.mass_ratio, + zeta=s.zeta, fn=s.frequency, Tn=1./s.frequency, y0=0., v0=0.), warmup_steps=warmup, requested_steps=count, dt=1., D=D, U=U, nu=NU, - duration_Tn=count/PERIOD, duration_tU_D=count*U/D, + duration_Tn=count*s.frequency, duration_tU_D=count*U/D, + summary_block_steps=10*PERIOD, rhs='sum of three instantaneous CME_y + Wen_y, one-step lag, no filtering') try: with Simulation(str(out/'config.json')) as sim: @@ -135,7 +145,7 @@ def main(): f_last = float(cme.sum(dtype=np.float64)+wen.sum(dtype=np.float64)) if completed % PERIOD == 0: rows.flush() - print(json.dumps(dict(event='period', step=completed, t_Tn=completed/PERIOD, + print(json.dumps(dict(event='progress', step=completed, t_Tn=completed*s.frequency, y_D=s.y/D, elapsed_s=time.monotonic()-started)), flush=True) else: status = 'smoke_completed' if args.smoke else 'duration_completed_not_qualified' diff --git a/tests/unit/test_steady_pinball.py b/tests/unit/test_steady_pinball.py index dffd806..d76394f 100644 --- a/tests/unit/test_steady_pinball.py +++ b/tests/unit/test_steady_pinball.py @@ -26,7 +26,7 @@ def test_actual_loader_and_paper_geometry(tmp_path): def test_structure_rest_equilibrium_and_energy_decay(): - s = runner.Structure() + s = runner.Structure(zeta=.01) assert s.candidate(0.) == (0., 0.) assert s.mass == pytest.approx(10*3*math.pi*100) s.y, s.v = .2, 0. @@ -52,3 +52,33 @@ def test_report_uses_completed_rows_and_net_force(): assert got[0]['end_step'] == 4 assert got[0]['y_rms_D'] == .1 assert got[0]['fy_rms'] == 6. + + +@pytest.mark.parametrize('case', ['Z0', 'M2', 'R2']) +def test_zero_damping_cases_conserve_unforced_energy_and_match_forced_solution(case): + s = runner.Structure(*runner.STRUCTURES[case]) + assert s.zeta == s.damping == 0. + s.y, s.v = .2, 1.e-5 + y0, v0 = s.y, s.v + energy = .5*s.mass*v0**2 + .5*s.stiffness*y0**2 + for _ in range(12345): + s.y, s.v = s.candidate(0.) + assert .5*s.mass*s.v**2 + .5*s.stiffness*s.y**2 == pytest.approx(energy, rel=2e-11) + force, elapsed = 1.e-4, 3210 + s.y, s.v = y0, v0 + for _ in range(elapsed): + s.y, s.v = s.candidate(force) + phase = s.omega*elapsed + equilibrium = force/s.stiffness + expected_y = equilibrium+(y0-equilibrium)*math.cos(phase)+v0*math.sin(phase)/s.omega + expected_v = -(y0-equilibrium)*s.omega*math.sin(phase)+v0*math.cos(phase) + assert (s.y, s.v) == pytest.approx((expected_y, expected_v), abs=2e-11) + + +def test_predeclared_comparisons_change_structure_only(): + z, m, r = [runner.Structure(*runner.STRUCTURES[k]) for k in ['Z0', 'M2', 'R2']] + assert m.mass/z.mass == pytest.approx(.2) + assert m.stiffness/z.stiffness == pytest.approx(.2) + assert r.mass == m.mass + assert r.stiffness/m.stiffness == pytest.approx(.36) + assert [runner.STEPS*s.frequency for s in [z, m, r]] == pytest.approx([60., 60., 36.])