From 516975ef6482f983ace347cb10de7450d2fbd680 Mon Sep 17 00:00:00 2001 From: Frank14f <1515444314@qq.com> Date: Sat, 31 Oct 2026 23:21:27 +0800 Subject: [PATCH] Add bounded paired Karman pinball VIV research driver --- docs/karman_viv.md | 83 +++++ docs/research/moving_boundary_force_spikes.md | 99 ++++++ research/karman_viv/run.py | 292 ++++++++++++++++++ tests/unit/test_karman_viv.py | 115 +++++++ 4 files changed, 589 insertions(+) create mode 100644 docs/karman_viv.md create mode 100644 docs/research/moving_boundary_force_spikes.md create mode 100644 research/karman_viv/run.py create mode 100644 tests/unit/test_karman_viv.py diff --git a/docs/karman_viv.md b/docs/karman_viv.md new file mode 100644 index 0000000..add7e77 --- /dev/null +++ b/docs/karman_viv.md @@ -0,0 +1,83 @@ +# Re_D = 50 Kármán inflow: paired pinball VIV / frozen PPO demo + +Development experiment, frozen before results, 2026-10-31. No new production +method, structural damping, force filter, policy training, or main merge. +[Force-spike research index](research/moving_boundary_force_spikes.md). + +## Plant and structure + +Base `e9213dd2cc9cfed0e8487a3cf4ebdd94d23dc64d`; fixed-policy contract from +`580397eaf83c85b233495791ecc51662a374a247`. Reuse `steady_pinball` config, +coordinate map and exact held-force oscillator. D40, U=.005, nu=.004, MRT, +2559 x 1022, FP32 unshifted double buffers. Regularized parabolic inlet, +NEQ extrapolation outlet, bounce-back y walls (.5,1020.5); no BC changes. +D20 coordinates map x'=2x, y'=.5+2(y-.5). +Upstream cylinder (400,510.5), radius40, is static. Pinballs radius20 at +(1200,510.5), (1252,540.5), (1252,480.5) share only transverse translation. +Lab probes radius10 at (1600,590.5), (1600,510.5), (1600,430.5) stay fixed. + +One shared `Translation`, travel ±D, max wall speed .025, existing clearance +and per-step move guards unchanged. M2: m*=2 based on all three disks, +M=7539.8223686, fn=.000025, K=.000186037660, C=zeta=0; release y=v=0. +Each native step uses previous-step sum(CME_y+Wen_y), then exact oscillator +update with held force. Wall speed is interval displacement y[n+1]-y[n]. +Cover/uncover impulses are logged, never added to the selected RHS. This +estimator remains provisional; existing spike/convergence limitations remain. + +## Same release state, two clocks + +A static precursor retains the accepted Kármán sequence: upstream+probes, +2,048,000-step warmup, 150 x 3,200-step no-pinball target samples, add static +pinballs/sync, second 2,048,000-step warmup. Save both DDF buffers, flags, +clock and last instantaneous force. Construct each moving arm fresh at the +same zero-offset geometry, check flags, copy both buffers and clock, verify +exact copy hashes. This one-time research initialization is not a moving +checkpoint or FSI rollback. Initial policy observation and spin are zero. + +Frozen model SHA256 `148240c1dcb8b11d8e5a0a9e991f6874c4e799c9744707bf0fb32a8efeb0468d`; +normalizer `534fda818c030a56960e95d0f3ac1e350d887d03db7f80f1d6453195a502b7f0`. +Only these explicit supplied assets are accepted. Policy is CPU deterministic. +Native policy signals: packed CME (same `read_force` path), lab area-mean +velocity, each averaged over SI3200. Convert to original D20 training units: +forces x2, area-mean velocity x2 x78, then frozen force-first normalization +and clipping. Fine-grid probe pixel footprint is not exactly identical; +using the original78 cells avoids a hidden area-count normalization change. +The moving plant, motion-induced force, lab-frame wake displacement and D40 +resolution remain out of training distribution. No transfer guarantee. + +A=(8a+(0,-4,4))*.005; EMA alpha=.9^(1/4), reset0 each SI; omega=+A_eff/20. +This preserves historical physical actuator timescale under dt/4. +The zero arm substitutes a=(0,.5,-.5), so all spins are identically zero. +Policy interval averaging does **not** average the native structural RHS. + +## Bound, output, and interpretation + +One serial pair, zero then controlled: 400SI=1,280,000 native steps=160D/U +=32Tn per arm. No retuning/automatic extension. Stop on nonfinite state, +no_donor, travel/speed/move rejection, or four-hour whole-pair wall cap. +A failure stops the remaining run; retain valid rows and completed intervals. +NaN trailing rows are unexecuted, not zero-valued observations. + +Report first/last80D/U windows, mean offset separately from centered RMS, +half-range, peak, frequency and drift; native versus SI-averaged force, +CME/Wen/node/crossing terms and action/observation clipping. Last80D/U is +the principal descriptive window, not an asserted settled window. If zero +arm half-range <.1D there, classify weak response and do not claim vibration +suppression. Report continuous ratios regardless. Cloaking uses retained +canonical last30 lagged sensor DTW against this precursor's no-pinball target, +with fixed-body D20 results contextual only. These are development criteria, +not paper uncertainty, grid convergence, or proof of validated FSI. + +Save native18-column rows (~184MB/arm), small SI signals, target, config, +release two buffers (~188MiB), hashes and terminal reports outside checkout. +No per-step field output. Expected new data <650MB, check actual resources. + +``` +CUDA_VISIBLE_DEVICES=0 PYTHONPATH=src python research/karman_viv/run.py \ + --model /tmp/celeris-karman-viv-assets/model.zip \ + --norm /tmp/celeris-karman-viv-assets/norm.json --out /tmp/NEW_OUTPUT +``` + +First `--smoke`: same D40 geometry, warmup32/target2x8/warmup32, +4x8 response steps each arm; strictly plumbing only. CPU source checks and +independent exact-tip review precede GPU. Preserve original DynamisLab files. diff --git a/docs/research/moving_boundary_force_spikes.md b/docs/research/moving_boundary_force_spikes.md new file mode 100644 index 0000000..2d2f3ff --- /dev/null +++ b/docs/research/moving_boundary_force_spikes.md @@ -0,0 +1,99 @@ +# Moving-boundary force spikes: evidence and source index + +Updated 2026-10-31. + +This note indexes supplied audited observations and primary-source extractions. +It does not report a newly measured cure or implement an alternative method. +Paper details below were supplied by the parent investigation; no new source +retrieval or independent paper verification was performed for this note. + +## Current implementation and interpretation + +- Celeris M1 reconstructs a fresh node as `feq(rho_bar, u_wall) + donor NEQ`. +- In exact arithmetic, the donor non-equilibrium component has zero mass and + momentum when formed against the donor's own equilibrium moments. +- Re-projecting that same reconstruction onto the same wall velocity therefore + does not, by itself, change its mass or momentum or supply a cure for spikes. +- Floating-point residuals should be measured separately from this algebraic fact. +- Peng's alternative differs upstream: extrapolate PDFs from three donors, then + replace selected moments. It is not simply another projection of current M1. +- PPO interval-force observations and native per-step FSI right-hand-side forces + are distinct signals; comparisons must identify which signal is being used. + +## Audited run observations supplied to this note + +- D40 revision: `e9213dd2cc9cfed0e8487a3cf4ebdd94d23dc64d`. +- Window: first `97.303 D/U` of physical time. +- Maximum displacement: `max |y|/D = 0.304893`. +- Maximum native force-coefficient increment: `max |ΔCy| = 0.66688`. +- D20 over the same physical-time window: `max |ΔCy| = 1.80735`. +- The largest D40 event involved 30 nodes. +- These are supplied audited observations, not measurements made for this note. +- Time steps and trajectories differ between runs. This comparison establishes + neither spatial convergence nor a unique cause for the force spikes. + +Evidence locations supplied by the investigation: + +- Horsea: `/tmp/celeris-results-pinball-resolution`. +- Mac: `validation-evidence/pinball-force-diagnosis`. + +These locations index run evidence, not local copies of the cited papers. + +## Tao, Hu, and Guo (2016): time-step decoupling + +Computers & Fluids **133**, 1–14. DOI: +[10.1016/j.compfluid.2016.04.009](https://doi.org/10.1016/j.compfluid.2016.04.009). +[Publisher abstract](https://www.sciencedirect.com/science/article/pii/S0045793016301116). + +- Abstract-level support only: non-equilibrium refill reduces oscillations. +- With local iteration, the compared refill methods give similar oscillation magnitudes. +- Time-step decoupling reduces oscillations further in the reported sedimentation, + Couette, and Poiseuille examples. +- The exact decoupling schedule and time-step ratios have not been verified. +- This abstract does not establish equivalence to an arbitrary force filter. +- Do not implement a claimed reproduction until the actual methods are available. + +## Marson et al. (2021): local iterations + +Physical Review E **103**, 053308. DOI: +[10.1103/PhysRevE.103.053308](https://doi.org/10.1103/PhysRevE.103.053308). +[Open repository record](https://archive-ouverte.unige.ch/unige:151868). + +- §VI, p. 18 discusses pressure waves associated with fresh nodes. +- Appendix E.3, p. 25 and Fig. 18, p. 24 describe the local procedure. +- Copy a local sublattice, perform local collision and streaming, write the fresh + nodes, and repeat before the global collision and streaming operation. +- One to three iterations performed best in the reported tests; this is not a + universal iteration-count recommendation. +- The reported smoothed torque-noise assessment is qualitative here. +- Its thin-shell setting differs from the current moving-body configuration. + +## Peng (2016): extrapolation and moment replacement + +Computers & Mathematics with Applications **72**, 349–374. DOI: +[10.1016/j.camwa.2015.08.027](https://doi.org/10.1016/j.camwa.2015.08.027). + +- §5.1.2, Eq. (50), p. 368: extrapolate PDFs along an outward normal direction: + `f_new(x) = 3 f(x + e) - 3 f(x + 2e) + f(x + 3e)`. +- §5.1.3, Eq. (51): wall equilibrium plus donor non-equilibrium is the family + containing the current M1 reconstruction. +- §5.1.4, Eqs. (52)–(54): in MRT space, replace momentum with `j = rho0 u_wall`, + preserve the other moments, then apply the inverse transform. +- A compressible Celeris adaptation must choose the appropriate density; + copying the paper's `rho0` blindly is not justified. +- Three valid donors require new eligibility checks and fallback behavior beyond + the current one-link shared-body clearance rule. +- §§5.2–5.3 retain nonzero fluctuations in specific configurations; these results + do not promise elimination of force spikes in the present problem. + +## Bounded next investigation + +1. Retain current M1 as the baseline and compare matched prescribed trajectories + over matched physical-time windows. +2. Record raw force RHS, node crossings, donor moments, mass and momentum changes, + impulse, spectra, and trajectory response together. +3. Evaluate local iterations, Peng-style reconstruction, and force filtering as + separate interventions, with explicit definitions and diagnostics. +4. Obtain the actual time-step-decoupling methods before implementing that method. + +The current task indexes evidence only; it makes no alternative-method change. diff --git a/research/karman_viv/run.py b/research/karman_viv/run.py new file mode 100644 index 0000000..d899369 --- /dev/null +++ b/research/karman_viv/run.py @@ -0,0 +1,292 @@ +"""One paired zero/PPO Karman pinball VIV experiment; see docs/karman_viv.md. + +Research-only adapter: static precursor -> two fresh, identical moving plants. +No moving checkpoint/rollback, policy retraining, or force filtering. +""" +from __future__ import annotations + +import argparse +import hashlib +import json +from pathlib import Path +import subprocess +import sys +import time + +import numpy as np + +ROOT = Path(__file__).resolve().parents[2] +sys.path.insert(0, str(ROOT)) +from research.steady_pinball.run import Structure, config, scaled_points, CENTRES, PROBES, COLUMNS + +BIAS = np.array([0., -4., 4.]) +MODEL_SHA = '148240c1dcb8b11d8e5a0a9e991f6874c4e799c9744707bf0fb32a8efeb0468d' +NORM_SHA = '534fda818c030a56960e95d0f3ac1e350d887d03db7f80f1d6453195a502b7f0' +STEP_COLUMNS = (*COLUMNS, 'omega_F', 'omega_U', 'omega_L') + + +def sha(path): + h = hashlib.sha256() + with open(path, 'rb') as f: + for block in iter(lambda: f.read(1024*1024), b''): + h.update(block) + return h.hexdigest() + + +def array_sha(a): + return hashlib.sha256(memoryview(np.ascontiguousarray(a)).cast('B')).hexdigest() + + +def observation(sensor_mean, force_mean, scale, norm): + """Fine-grid area means/CME -> trained D20 sensor sums/forces -> frozen norm.""" + sensors = np.asarray(sensor_mean, np.float64)*scale*78. + forces = np.asarray(force_mean, np.float64)*scale + raw = np.asarray([*sensors, *forces], np.float32) + preclip = np.hstack((raw[6:]/norm['force_norm_fact'], + (raw[:6]-np.asarray(norm['sens_deviation'], np.float32)) / + np.asarray(norm['sens_norm_fact'], np.float32))) + if not np.isfinite(preclip).all(): + raise RuntimeError('nonfinite policy observation') + return raw, np.clip(preclip, -1, 1).astype(np.float32), int((np.abs(preclip)>1).sum()) + + +def wall_target(action, scale): + return (8.*np.asarray(action)+BIAS)*(.01/scale) + + +def ema(eff, target, scale): + alpha = .9**(1./scale**2) + return alpha*eff + (1.-alpha)*target + + +def load_policy(model): + # Same inference spaces as legacy_test.core.model_loader, without importing it. + import gymnasium as gym + from stable_baselines3 import PPO + class InferenceEnv(gym.Env): + observation_space = gym.spaces.Box(-1., 1., (12,), np.float32) + action_space = gym.spaces.Box(-1., 1., (3,), np.float32) + def reset(self, *, seed=None, options=None): + super().reset(seed=seed) + return np.zeros(12, np.float32), {} + def step(self, action): + return np.zeros(12, np.float32), 0., False, False, {} + return PPO.load(str(model), env=InferenceEnv(), device='cpu') + + +def add_scene(sim, scale, moving=False, pinball=True): + from CelerisLab.body.motion import Translation + sim.add_body('circle', scaled_points(((200., 255.5),), scale)[0], 20.*scale) + probes = [sim.add_body('sensor', c, 5.*scale) for c in scaled_points(PROBES, scale)] + bodies = [] + if pinball: + envelope = Translation((0., 0., -20.*scale, 20.*scale), .05/scale, shared=True) if moving else None + bodies = [sim.add_body('circle', c, 10.*scale, translation=envelope) + for c in scaled_points(CENTRES, scale)] + return probes, bodies + + +def capture_seed(sim): + """Capture both FP32 buffers at fixed geometry; never a moving-state checkpoint.""" + import pycuda.driver as cuda + with sim.ctx.activate(): + sim.stream.synchronize() + buffers = [np.empty_like(sim.field.ddf) for _ in range(2)] + for dst, src in zip(buffers, (sim.field.ddf_gpu, sim.field.temp_gpu)): + cuda.memcpy_dtoh(dst, src) + return (*buffers, sim.get_flags(), sim.stepper.step_count) + + +def install_seed(sim, seed): + """Only a new, zero-offset instance may receive the static release state.""" + import pycuda.driver as cuda + current, spare, flags, clock = seed + if sim.stepper.step_count != 0 or any(a.shape != sim.field.ddf.shape or a.dtype != sim.field.ddf.dtype + for a in (current, spare)): + raise ValueError('release transfer requires a fresh matching field') + if not np.array_equal(sim.get_flags(), flags): + raise ValueError('static/moving release flags differ') + for obj in sim.bodies.objects: + if obj.motion is not None and (obj.motion.applied != (0., 0.) or obj.motion.pending is not None): + raise ValueError('release transfer requires zero geometry and no target') + with sim.ctx.activate(): + sim.stream.synchronize() + cuda.memcpy_htod(sim.field.ddf_gpu, current) + cuda.memcpy_htod(sim.field.temp_gpu, spare) + sim.stepper._step_count = clock # one-time research initialization, not restore() + sim.field.invalidate_host_ddf_cache() + copied = capture_seed(sim) + if copied[3] != clock or any(not np.array_equal(a, b, equal_nan=True) for a, b in zip(copied[:3], seed[:3])): + raise RuntimeError('release transfer did not preserve both buffers/flags') + return dict(current=array_sha(copied[0]), spare=array_sha(copied[1]), flags=array_sha(copied[2]), clock=copied[3]) + + +def check_time(start, cap): + if time.monotonic()-start > cap: + raise TimeoutError('paired experiment wall cap') + + +def warmup(sim, count, start, cap): + while count: + check_time(start, cap) + n = min(count, 10000) + sim.run(n) + count -= n + + +def precursor(cfgpath, scale, warm, si, targets, out, start, cap): + from CelerisLab import Simulation + with Simulation(str(cfgpath)) as sim: + probes, _ = add_scene(sim, scale, pinball=False) + sim.initialize() + warmup(sim, warm, start, cap) + target = [] + for _ in range(targets): + check_time(start, cap) + sim.run(si) + target.append([v*scale*78. for p in probes for v in sim.read_sensor(p)[:2]]) + np.save(out/'target.npy', np.asarray(target, np.float32)) + first = sim.bodies.count + for c in scaled_points(CENTRES, scale): + sim.add_body('circle', c, 10.*scale) + sim.sync_bodies() + bodies = list(range(first, first+3)) + warmup(sim, warm-1, start, cap) + sim.run(1) + force = sum(float(sim.read_force(p)[1]) for p in bodies) + seed = capture_seed(sim) + np.savez(out/'release.npz', current=seed[0], spare=seed[1], flags=seed[2], clock=seed[3], force=force) + print(json.dumps(dict(event='precursor_complete', clock=seed[3], force=force, + elapsed_s=time.monotonic()-start)), flush=True) + return seed, force + + +def run_arm(name, cfgpath, scale, si, intervals, seed, initial_force, model, norm, out, start, cap): + from CelerisLab import Simulation + from CelerisLab.body.motion import MovingTables + cols = {k:i for i, k in enumerate(MovingTables.DIAG)} + diameter = 20.*scale + structure = Structure(2., 1.e-4/scale**2, 0., diameter) + arm = out/name + arm.mkdir() + rows = np.lib.format.open_memmap(arm/'steps.npy', mode='w+', dtype=np.float64, + shape=(si*intervals, len(STEP_COLUMNS))) + rows[:] = np.nan + signals = {k:np.full((intervals, width), np.nan, np.float32) for k, width in + [('sensors',6), ('forces',6), ('actions',3), ('observations',12)]} + completed = completed_intervals = clipped = 0 + status, transfer = 'running', None + force = initial_force + obs = np.zeros(12, np.float32) + arm_start = time.monotonic() + try: + with Simulation(str(cfgpath)) as sim: + probes, bodies = add_scene(sim, scale, moving=True) + sim.initialize() + transfer = install_seed(sim, seed) + for i in range(intervals): + action = -BIAS/8 if name == 'zero' else np.asarray(model.predict(obs, deterministic=True)[0], np.float32) + if action.shape != (3,) or not np.isfinite(action).all() or (np.abs(action)>1).any(): + raise RuntimeError('invalid policy action') + target, eff = wall_target(action, scale), np.zeros(3) + sensor_sum, force_sum = np.zeros(6), np.zeros(6) + for k in range(si): + if k % 1000 == 0: + check_time(start, cap) + y1, v1 = structure.candidate(force) + if not np.isfinite([y1,v1]).all() or abs(y1)>diameter: + raise RuntimeError('nonfinite structure or |y| > D; no retuning') + eff = ema(eff, target, scale) + omega = eff/(diameter/2) + for p, w in zip(bodies, omega): + sim.set_body(p, omega=float(w), offset=(0., y1), vy=y1-structure.y) + sim.run(1, zero_obs=True, sync_obs=True) + d = sim.moving_diagnostics() + cme, wen = d[:,cols['cme_y']], d[:,cols['wen_y']] + node = d[:,cols['cover_jy']]-d[:,cols['uncover_jy']] + row = (completed+1,y1,v1,*cme,*wen,*node, + d[:,cols['eq_only']].sum(),d[:,cols['no_donor']].sum(),d[:,cols['crossings']].sum(),*omega) + if not np.isfinite(row).all() or row[13] != 0: + raise RuntimeError('nonfinite moving diagnostic or no_donor') + f = np.array([v for p in bodies for v in sim.read_force(p)[:2]]) + s = np.array([v for p in probes for v in sim.read_sensor(p)[:2]]) + if not np.isfinite([f,s]).all(): + raise RuntimeError('nonfinite force/sensor observation') + # Policy sees packed CME as in fixed-body training; structure sees CME+Wen. + force_sum += f + sensor_sum += s + force = float(cme.sum(dtype=np.float64)+wen.sum(dtype=np.float64)) + rows[completed] = row + completed += 1 + structure.y, structure.v = y1,v1 + raw, obs, nclip = observation(sensor_sum/si, force_sum/si, scale, norm) + clipped += nclip + for key, val in [('sensors',raw[:6]),('forces',raw[6:]),('actions',action),('observations',obs)]: + signals[key][i] = val + completed_intervals = i+1 + if (i+1)%10 == 0 or i+1 == intervals: + rows.flush() + np.savez(arm/'signals.npz', **signals) + print(json.dumps(dict(event='progress', arm=name, interval=i+1, step=completed, + tU_D=completed*.01/scale/diameter, y_D=y1/diameter, + elapsed_s=time.monotonic()-arm_start)), flush=True) + status = 'duration_completed_not_qualified' + except BaseException as exc: + status = f'error: {type(exc).__name__}: {exc}' + raise + finally: + rows.flush() + np.savez(arm/'signals.npz', **signals) + report = dict(status=status, completed_steps=completed, completed_intervals=completed_intervals, + elapsed_s=time.monotonic()-arm_start, release_transfer=transfer, clipped_observation_channels=clipped, + columns=STEP_COLUMNS, M=structure.mass,K=structure.stiffness,fn=structure.frequency,damping=0., + max_abs_y_D=float(np.max(np.abs(rows[:completed,1]))/diameter) if completed else None) + (arm/'report.json').write_text(json.dumps(report,indent=2)+'\n') + + +def main(): + ap = argparse.ArgumentParser(description=__doc__) + ap.add_argument('--model', type=Path, required=True) + ap.add_argument('--norm', type=Path, required=True) + ap.add_argument('--out', type=Path, required=True) + ap.add_argument('--smoke', action='store_true', help='short plumbing only; no physical conclusion') + ap.add_argument('--wall-seconds', type=float, default=14400.) + args = ap.parse_args() + if not np.isfinite(args.wall_seconds) or args.wall_seconds<=0: + ap.error('positive finite wall cap required') + if sha(args.model)!=MODEL_SHA or sha(args.norm)!=NORM_SHA: + ap.error('model/normalizer must match frozen provenance') + out = args.out.resolve() + if out == ROOT or ROOT in out.parents: + ap.error('output must be outside checkout') + out.mkdir(parents=True, exist_ok=False) + norm = json.loads(args.norm.read_text()) + model = load_policy(args.model) + scale = 2. + warm, si, targets, intervals = (32,8,2,4) if args.smoke else (2048000,3200,150,400) + cfgpath = out/'config.json' + cfgpath.write_text(json.dumps(config(40),indent=2)+'\n') + receipt = dict(source=subprocess.check_output(['git','-C',str(ROOT),'rev-parse','HEAD'],text=True).strip(), + model_sha256=sha(args.model),norm_sha256=sha(args.norm),smoke=args.smoke,D=40,U=.005,nu=.004, + warmup_each=warm,target_intervals=targets,SI=si,response_intervals=intervals, + response_steps=si*intervals,response_tU_D=si*intervals*.005/40, + rhs='previous native step sum CME_y+Wen_y; policy uses interval packed CME', + probe_frame='lab-fixed; fine area mean * 2 * trained 78 cells', + ema_native_alpha=.9**.25, structural_damping=0., columns=STEP_COLUMNS) + (out/'case.json').write_text(json.dumps(receipt,indent=2)+'\n') + start = time.monotonic() + status = 'running' + try: + seed, force = precursor(cfgpath,scale,warm,si,targets,out,start,args.wall_seconds) + for name in ('zero','controlled'): + run_arm(name,cfgpath,scale,si,intervals,seed,force,model,norm,out,start,args.wall_seconds) + status = 'smoke_completed' if args.smoke else 'pair_completed_not_qualified' + except BaseException as exc: + status = f'error: {type(exc).__name__}: {exc}' + raise + finally: + (out/'terminal.json').write_text(json.dumps(dict(status=status,elapsed_s=time.monotonic()-start),indent=2)+'\n') + + +if __name__ == '__main__': + main() diff --git a/tests/unit/test_karman_viv.py b/tests/unit/test_karman_viv.py new file mode 100644 index 0000000..efab574 --- /dev/null +++ b/tests/unit/test_karman_viv.py @@ -0,0 +1,115 @@ +"""Research interface contracts; no CUDA context or physical validation.""" +import importlib.util +from pathlib import Path +from types import SimpleNamespace +from unittest.mock import patch + +import numpy as np +import pytest + +ROOT = Path(__file__).resolve().parents[2] +spec = importlib.util.spec_from_file_location('karman_viv_research_driver', ROOT/'research/karman_viv/run.py') +m = importlib.util.module_from_spec(spec) +spec.loader.exec_module(m) + + +def test_observation_unit_mapping_and_order(): + norm = dict(force_norm_fact=2.,sens_deviation=np.arange(6)/10.,sens_norm_fact=np.ones(6)*3) + s=np.arange(6)*.001; f=np.arange(6)*.01 + raw,obs,nclip=m.observation(s,f,2,norm) + np.testing.assert_array_equal(raw,np.asarray([*(s*156),*(f*2)],np.float32)) + expected=np.hstack((raw[6:]/2,(raw[:6]-np.asarray(norm['sens_deviation'],np.float32))/3)) + np.testing.assert_array_equal(obs,expected.astype(np.float32)) + assert nclip==0 + + +def test_same_physical_signal_d20_d40(): + norm=dict(force_norm_fact=1.,sens_deviation=np.zeros(6),sens_norm_fact=np.ones(6)) + a=m.observation(np.ones(6)*.001,np.arange(6)*.01,1,norm) + b=m.observation(np.ones(6)*.0005,np.arange(6)*.005,2,norm) + np.testing.assert_array_equal(a[0],b[0]); np.testing.assert_array_equal(a[1],b[1]) + + +def test_zero_action_and_ccw_scale(): + np.testing.assert_array_equal(m.wall_target(-m.BIAS/8,2),np.zeros(3)) + np.testing.assert_allclose(m.wall_target(np.ones(3),2),[.04,.02,.06]) + assert m.wall_target(np.ones(3),2)[0]/20 == .002 + + +def test_ema_equal_physical_time_and_interval_reset(): + target=np.array([.01,-.02,.03]); coarse=np.zeros(3); fine=np.zeros(3) + for _ in range(13): coarse=m.ema(coarse,target,1) + for _ in range(52): fine=m.ema(fine,target,2) + np.testing.assert_allclose(fine,coarse,rtol=2e-15,atol=1e-17) + np.testing.assert_allclose(m.ema(np.zeros(3),target,2),(1-.9**.25)*target) + + +def test_config_real_parser(tmp_path): + import json + from CelerisLab.config import load_lbm_config + path=tmp_path/'config.json';path.write_text(json.dumps(m.config(40))) + cfg=load_lbm_config(str(path));cfg.validate() + assert (cfg.nx,cfg.ny,cfg.velocity,cfg.viscosity)==(2559,1022,.005,.004) + + +def test_structure_physical_scaling_and_no_damping(): + coarse=m.Structure(2,1e-4,0,20);fine=m.Structure(2,2.5e-5,0,40) + assert fine.mass==4*coarse.mass + assert fine.stiffness==coarse.stiffness/4 + assert fine.damping==0 + assert fine.candidate(0)==(0.,0.) + + +def test_observation_clipping_and_nonfinite(): + norm=dict(force_norm_fact=1.,sens_deviation=np.zeros(6),sens_norm_fact=np.ones(6)) + _,obs,n=m.observation(np.ones(6),np.ones(6),2,norm) + assert n==12 and np.array_equal(obs,np.ones(12,np.float32)) + with pytest.raises(RuntimeError,match='nonfinite'): + m.observation(np.full(6,np.nan),np.zeros(6),2,norm) + + +def test_arm_step_order_and_native_rhs(tmp_path): + import CelerisLab + from CelerisLab.body.motion import MovingTables + from contextlib import nullcontext + log=[];candidates=[] + class Structure: + y=v=0.;mass=2.;stiffness=3.;frequency=.1 + def __init__(self,*args): pass + def candidate(self,f): + candidates.append(f) + return self.y+.001,self.v+.002 + class Sim: + def __init__(self,*args): self.n=0 + def __enter__(self): return self + def __exit__(self,*args): pass + def initialize(self): pass + def set_body(self,p,**kw): log.append(('set',p,kw)) + def run(self,n,**kw): self.n+=1;log.append(('run',n,kw)) + def moving_diagnostics(self): + a=np.zeros((3,len(MovingTables.DIAG)),np.float32) + a[:,MovingTables.DIAG.index('cme_y')]=self.n + a[:,MovingTables.DIAG.index('wen_y')]=10*self.n + a[:,MovingTables.DIAG.index('cover_jy')]=1000 + return a + def read_force(self,p): return [0.,float(self.n),0.] + def read_sensor(self,p): return [.001,0.,0.] + seen=[] + model=SimpleNamespace(predict=lambda obs,**kw:(seen.append(obs.copy()) or np.ones(3)*.1,None)) + norm=dict(force_norm_fact=1.,sens_deviation=np.zeros(6),sens_norm_fact=np.ones(6)) + with patch.object(CelerisLab,'Simulation',Sim),patch.object(m,'Structure',Structure), \ + patch.object(m,'add_scene',return_value=([1,2,3],[4,5,6])),patch.object(m,'install_seed',return_value={}), \ + patch.object(m,'check_time'): + m.run_arm('controlled','cfg',2,2,2,None,7.,model,norm,tmp_path,0,1) + assert candidates==[7.,33.,66.,99.] # previous native force, NOT SI mean or node impulse + assert not seen[0].any() and len(seen)==2 + assert [e[0] for e in log]==['set','set','set','run']*4 + for e in log: + if e[0]=='set': assert e[2]['vy']==pytest.approx(.001) and e[2]['omega']!=0 + else: assert e[2]==dict(zero_obs=True,sync_obs=True) + rows=np.load(tmp_path/'controlled/steps.npy') + np.testing.assert_allclose(rows[:,1],[.001,.002,.003,.004]) + assert rows.shape==(4,18) + np.testing.assert_array_equal(rows[0,15:],rows[2,15:]) # EMA resets at SI boundary + signals=np.load(tmp_path/'controlled/signals.npz') + np.testing.assert_array_equal(signals['forces'][:,1],[3.,7.]) # interval CME *2, not CME+Wen