From 63e5c3828b9bf45ec375b750ab24e26eb093975c Mon Sep 17 00:00:00 2001 From: Frank14f <1515444314@qq.com> Date: Tue, 3 Nov 2026 01:17:45 +0800 Subject: [PATCH] Collect paired moving-pinball Fig2 research material --- docs/fig2_collection.md | 109 ++++++ research/fig2_collection/analyze.py | 130 +++++++ research/fig2_collection/plot.py | 93 +++++ research/fig2_collection/run.py | 548 ++++++++++++++++++++++++++++ research/fig2_materials/metrics.py | 94 +++++ tests/unit/test_fig2_analysis.py | 36 ++ tests/unit/test_fig2_collection.py | 225 ++++++++++++ tests/unit/test_fig2_metrics.py | 49 +++ 8 files changed, 1284 insertions(+) create mode 100644 docs/fig2_collection.md create mode 100644 research/fig2_collection/analyze.py create mode 100644 research/fig2_collection/plot.py create mode 100644 research/fig2_collection/run.py create mode 100644 research/fig2_materials/metrics.py create mode 100644 tests/unit/test_fig2_analysis.py create mode 100644 tests/unit/test_fig2_collection.py create mode 100644 tests/unit/test_fig2_metrics.py diff --git a/docs/fig2_collection.md b/docs/fig2_collection.md new file mode 100644 index 0000000..2cf4cd3 --- /dev/null +++ b/docs/fig2_collection.md @@ -0,0 +1,109 @@ +# Fig. 2 material collection: modern moving plant + +Research-only child of `76873d911c4e2caba1bacc25aaa1495ea768453f`. +No production/kernel/manuscript or DynamisLab edits. These are fresh common-release +paired experiments, separate from the continuous OFF/ON/OFF demonstration videos. + +## Fixed plant and cases + +All cases use D40, Uref=.005, nu=.004, MRT/FP32/unshifted double buffering, +the existing mean-parabolic/regularized inlet, NEQ outlet and bounce-back walls. +One transverse assembly DOF: m*=4 based on all three disks, fnD/U=.2, zeta=0. +No force filter or retuning. Existing +/-D, speed.025, nonfinite and no-donor stops. +Structure RHS is previous native-step sum of pinball CME_y+Wen_y. Policy uses +SI-mean packed CME and fixed lab probe signals; these are different force paths. +SI=3200; EMA alpha=.9**.25 resets every SI as in the acquisition contract. +Steady video used SI4000; its descriptive numbers are not the new collection. + +| case | prelude/release clock | controller | collection | +|---|---|---|---| +| steady | new static pinball, 2,048,000 steps | constant wall speeds (0,-5.1,5.1)U, no PPO | 400 SI,160 D/U | +| Karman | verified saved static release at4,576,000 | historical `karman_re100` PPO; actual pinball Re_D50, scale8/bias(0,-4,4) | 400 SI,160 D/U | +| Taylor y0 | clean512,000 then static pinball512,000; clock1,024,000 | `vortex_taylor` PPO, scale4/samebias | 150 SI,60 D/U | + +Zero/control share exact current/spare populations, flags, clock and zero initial +structural state. Target is a separate no-pinball reference at the SAME absolute +clocks, not a copy of a solid-containing release. Its different geometry is the +definition of the reference, not identical microscopic initial conditions. +Karman target must run to4,576,000; old target.npy is from earlier clocks and is +not reused for same-clock field errors. Model/norm files are verified by SHA256 +in the driver and written with all source hashes into case.json. + +## Taylor injection + +Exactly one event at release in every role, inside the owning synchronized CUDA +context: centre(600,510.5), R80, strength **.000075**. Reuse the existing initializer +unchanged. It rewrites ALL interior cells to equilibrium and uploads BOTH buffers; +this removes existing nonequilibrium globally. The impulse and first two SI of +native force/displacement remain in the record. This is a declared limitation, +not an independently qualified physical vortex generator. + +Strength has inverse-time units: u_theta=(s*pi/2)r exp(-r²/R²). D20 R40/U.01/s.0003 +maps to D40 R80/U.005/s.000075 under diffusive scaling. Using .03*U at D40 would +double u/U and quadruple p/U²; actual-function CPU positive/negative controls lock it. + +Current Horsea `acquisition.py` SHA +`ea674080501407d00ded6d23467423c516cd37e04d166d5adda16465100204d3` +defines EIGHT offsets(-20,-15,-10,-5,0,5,10,15). `legacy_test/acquire.py` SHA +`91cd29a368dfc53a0804ec37397d2e5f25f0f8721f01b560a0f4328b8fd1dd9c` +selects the controlled centre-uy event: median first20, peak over38..87 inclusive. +All roles use that common event index; own target/role peak delays are separately +reported. The earlier five-offset remote summary is superseded by actual source. + +## Estimands and storage + +Primary: unaligned same-clock/event velocity-vector errors, Eq1 +`sqrt(mean_ROI((dux²+duy²)/Uref²))`. No target-energy denominator. Lab ROI is +front-equilibrium-centre-relative x/D=[4,24], y/D=[-5,5], equivalent to sensor-plane +x/D=[-6,14]; D40 inclusive columns1360..2160 and rows311..710. All ROI cells must +be fluid and finite in every compared field; record support coverage. + +Motion uses native last80 D/U for steady/Karman, whole60 D/U for Taylor; mean, +centered RMS and half-range remain distinct. Karman complete-cycle means and +eight independently phased nearest snapshots reuse actual acquisition helpers, +centre-uy column3 and minimum crossing gap. Missing crossings yield unavailable +phase diagnostics, never fabricated phase slots. Steady retains its last snapshot +and separately labeled fixed-window mean. Field collection for steady/Karman is +the last160 SI boundaries (96.4..160 D/U); motion window remains80..160 D/U. +Taylor retains its eight common event snapshots. Fixed-window mean error differs +from mean instantaneous error; both are labeled. Native six-channel DTW is a +supplementary training-space metric, never field recovery or physical delay; +its actual legacy lag uses flat sensor column1, unlike the phase helper's column3. + +Use one FP32 full-grid RAM cache per role, plus compact ROI RAM histories for +same-clock comparisons. Admit only with >=16GiB MemAvailable; no3.35GB disk field +cache. Save only selected fields, means, explicit masks, scalars and source receipts. +Nonfluid placeholders are zero-filled for helpers but always masked in plots and +scores; nonfinite FLUID is a hard failure. Serialize compilation/GPU work. Verify +compressed native histories before removing only redundant task scratch copies. + +Plot standards: vorticity curl(u)D/U, five-point fluid stencil, common RdBu_r21bins +with fixed[-3,3] limits and nearest interpolation; lab-fixed coordinates and ROI. +Mean plots exclude cells that were not fluid at every collected boundary. Selected +snapshots use actual step geometry. Raw arrays and scalar errors are not clipped. + +## V5 boundary + +Standard V5 `kar_re100` seeds41..45 DO exist (actual pinball Re_D50), but use +2000x600 uniform/free-slip and different absolute geometry/normalization/action +contracts. These are not silently transferred onto the present Legacy-contract +plant. No verified retained trained V5 steady/Taylor bundle was found; a vortex +environment prototype is not a trained policy. Existing V5 within-case results +may be indexed separately; cross-generation scores cannot imply an identical plant. + +## Execution + +First run focused CPU contracts, independent source review, then each scene with +`--smoke` (three32-step SI per role; no physical scores). Only the reviewed exact +SHA may run physical cases, one at a time, new outputs outside the checkout: + +``` +CUDA_VISIBLE_DEVICES=3 PYTHONPATH=src python research/fig2_collection/run.py --scene karman --out /tmp/celeris-results-fig2-collection/karman +python research/fig2_collection/analyze.py /tmp/celeris-results-fig2-collection/karman +python research/fig2_collection/plot.py /tmp/celeris-results-fig2-collection/karman +``` + +Use the same invocation for steady/taylor with new output directories. No automatic +extension or parameter sweep. All results are developmental material, not asserted +settled/physically qualified. Preserve partial failures and distinguish missing +evidence from passing gates. diff --git a/research/fig2_collection/analyze.py b/research/fig2_collection/analyze.py new file mode 100644 index 0000000..1c2d359 --- /dev/null +++ b/research/fig2_collection/analyze.py @@ -0,0 +1,130 @@ +"""Offline Fig2 estimands; preserve raw same-clock errors beside registration. + +No thresholds are fitted here. Probe DTW is the actual Legacy helper, not a +field metric; independent wake phase is a diagnostic, not a time correction. +""" +import argparse +import gc +import importlib.util +import json +from pathlib import Path +import sys +import types + +import numpy as np + +ROOT=Path(__file__).resolve().parents[2] +sys.path.insert(0,str(ROOT)) +from research.fig2_collection import run +from research.fig2_materials.metrics import field_error, motion_summary + + +def frequency(signal,dt): + """Largest non-DC rectangular-window FFT bin, descriptive resolution only.""" + a=np.asarray(signal,float);a=a-a.mean() + if len(a)<3 or np.max(np.abs(a))<1e-12:return None + power=np.abs(np.fft.rfft(a))**2 + i=1+int(np.argmax(power[1:])) + return dict(f_D_U=float(np.fft.rfftfreq(len(a),dt)[i]),bin_width_D_U=1/(len(a)*dt), + method='demeaned rectangular window; dominant non-DC bin; not a confidence interval') + + +def native_summary(rows,scene): + start=max(0,len(rows)-640000) if scene!='taylor' else 0 + a=rows[start:] + y=a[:,1]/40.;cy=a[:,3:9].sum(axis=1)/(.5*.005**2*40) + return dict(window_D_U=[start*.005/40,len(rows)*.005/40], + motion=motion_summary(y),motion_frequency=frequency(y,.005/40), + structural_RHS_Cy_RMS=float(np.sqrt(np.mean(cy**2))), + structural_RHS_Cy_mean=float(cy.mean()), + max_native_delta_Cy=float(np.abs(np.diff(cy)).max()) if len(cy)>1 else None, + force_frequency=frequency(cy,.005/40), + no_donor=float(rows[:,13].sum()),eq_only=float(rows[:,12].sum()), + crossings=float(rows[:,14].sum()), + force_convention='body force, sum pinball CME_y+Wen_y; native previous-step RHS, no filtering') + + +def compare_selected(target,current,u0=.005): + """Fixed lab ROI; all actual compared roles must be fluid and finite here.""" + ys,xs=slice(311,711),slice(1360,2161) + t,c=target['ux'][:,ys,xs],current['ux'][:,ys,xs] + tv,cv=target['uy'][:,ys,xs],current['uy'][:,ys,xs] + if t.shape!=c.shape or tv.shape!=cv.shape:raise ValueError('selected shape mismatch') + for fields in (target,current): + if not np.all((fields['flags'][:,ys,xs]&1)!=0):raise ValueError('nonfluid selected ROI') + if not all(np.isfinite(a).all() for a in (t,c,tv,cv)):raise ValueError('nonfinite selected ROI') + e2=((c.astype(float)-t)/u0)**2+((cv.astype(float)-tv)/u0)**2 + return dict(per_slot_E=np.sqrt(e2.mean(axis=(1,2))),aggregate_E=float(np.sqrt(e2.mean()))) + + +def main(): + ap=argparse.ArgumentParser(description=__doc__);ap.add_argument('directory',type=Path) + args=ap.parse_args();p=args.directory + case=json.loads((p/'case.json').read_text());term=json.loads((p/'terminal.json').read_text()) + if case['plan']['smoke'] or term['status']!='three_roles_completed_not_qualified': + raise ValueError('physical analysis requires all three complete roles; smoke/partial is unscored') + helpers=run.load_helpers() + # Reuse the canonical pure module under the same isolated read-only package. + native=run.load_file('_fig2_readonly_drl.legacy_test.metrics', + run.ACQUISITION.parent/'legacy_test/metrics.py') + reports={k:json.loads((p/k/'report.json').read_text()) for k in ('target','zero','controlled')} + sig={} + result=dict(scene=case['scene'],source=case['source'],Re_D=50.,roles={}, + primary='unaligned same-clock/event U-normalized vector errors; common lab ROI', + periodic_phase='independently registered centre-uy phase diagnostic, not a shared-clock comparison', + notes=['Modern moving-plant data; no Legacy numerical equivalence or physical qualification', + 'Motion, native structural force, policy SI means and probe DTW are different estimands']) + for role in ('target','zero','controlled'): + with np.load(p/role/'signals.npz') as z:sig[role]={k:z[k] for k in z.files} + if role!='target': + with np.load(p/role/'native.npz') as z: + result['roles'][role]=native_summary(z['rows'],case['scene']) + t=sig['target']['sensors'] + for role in ('zero','controlled'): + s=sig[role]['sensors'] + if not np.array_equal(sig[role]['absolute_steps'],sig['target']['absolute_steps']): + raise ValueError('probe absolute clocks mismatch') + result['roles'][role]['unaligned_probe_area_velocity_RMS_over_U']=float(np.sqrt(np.mean(((s/156-t/156)/.005)**2))) + result['roles'][role]['probe_channel_correlations']=native.channel_correlations(t,s) + if case['scene']=='taylor': + spec=types.SimpleNamespace(scene='vortex',conv_len=30) + dtw=native.native_similarity(spec,t,s,len(s)-1) + result['roles'][role]['native_probe_DTW']=dict(similarity=dtw,definition='source vortex training-window roll; supplementary, not event field recovery',conv_len=30) + own=helpers.select_vortex_event(s,'vortex_taylor') + target_event=helpers.select_vortex_event(t,'vortex_taylor') + result['roles'][role]['own_event_peak']=int(own['peak_index']) + result['roles'][role]['event_peak_delay_D_U']=(int(own['peak_index'])-int(target_event['peak_index']))*.4 + else: + result['roles'][role]['native_probe_DTW']=native.cycle_similarity(t,s,30) + result['roles'][role]['native_probe_DTW'].update(definition='canonical six-channel legacy cycle DTW, training D20 units', + lag_column=1,target_template_indices=[30,60],state_window_indices=[len(s)-30,len(s)], + registration_only=True,lag_identifiable=bool(case['scene']=='karman' and np.std(t[:,1])>1e-10)) + data={} + for role in ('target','zero','controlled'): + with np.load(p/role/'fields.npz') as z:data[role]={k:z[k] for k in z.files} + for role in ('zero','controlled'): + name={'karman':'independent_phase8_diagnostic','steady':'late_snapshot','taylor':'common_event8'}[case['scene']] + if case['scene']!='karman' or all(len(data[k]['ux'])==8 for k in (role,'target')): + result['roles'][role][name]=compare_selected(data['target'],data[role]) + else: + result['roles'][role][name]=dict(status='unavailable_insufficient_periodic_crossings') + result['roles'][role]['selection_metadata']=reports[role]['selection'] + result['roles'][role]['field_window_tU_D']=reports[role]['field_window_tU_D'] + if case['scene']=='karman' and all('cycle_mean_ux' in data[k] for k in (role,'target')): + support=np.zeros(data[role]['cycle_mean_ux'].shape,bool);support[311:711,1360:2161]=True + result['roles'][role]['E_complete_cycle_mean']=field_error( + data[role]['cycle_mean_ux'],data[role]['cycle_mean_uy'], + data['target']['cycle_mean_ux'],data['target']['cycle_mean_uy'],.005,support) + if case['scene']=='taylor': + clocks=[data[k]['absolute_steps'] for k in data] + if not all(np.array_equal(clocks[0],v) for v in clocks[1:]):raise ValueError('event snapshot clocks differ') + result['Taylor_injection']=dict(operator='all interior equilibrium rewrite; both buffers, identical event clock', + strength=.000075,center=[600.,510.5],radius=80.,impulse_window='first two SI; retained native history', + source_offsets=list(helpers.VORTEX_EVENT_OFFSETS)) + result['same_clock_field_errors']=json.loads((p/'same_clock_errors.json').read_text()) + run.write_json(p/'statistics.json',result) + del data;gc.collect() + print(json.dumps(dict(scene=case['scene'],statistics=str(p/'statistics.json')))) + + +if __name__=='__main__':main() diff --git a/research/fig2_collection/plot.py b/research/fig2_collection/plot.py new file mode 100644 index 0000000..c79fd5b --- /dev/null +++ b/research/fig2_collection/plot.py @@ -0,0 +1,93 @@ +"""Shared scientific plots of retained fields; no solver or synthetic fields.""" +import argparse +import json +from pathlib import Path +import sys + +import numpy as np +import matplotlib +matplotlib.use('Agg') +import matplotlib.pyplot as plt +from matplotlib.colors import BoundaryNorm +from matplotlib.patches import Circle,Rectangle + +ROOT=Path(__file__).resolve().parents[2] +sys.path.insert(0,str(ROOT)) +from research.continuous_viv.render import vorticity + + +def decorate(ax,scene,role,y,*,mean=False): + if not mean: + if scene=='karman':ax.add_patch(Circle((-20,0),1,facecolor='white',edgecolor='#253544',lw=.7)) + if role!='target': + for x,dy in ((0,0),(1.3,.75),(1.3,-.75)): + ax.add_patch(Circle((x,dy+y),.5,facecolor='white',edgecolor='#253544',lw=.7)) + ax.add_patch(Rectangle((4,-5),20,10,fill=False,ls='--',lw=.6,edgecolor='#555')) + ax.set_xlim(-23.6,27.6);ax.set_ylim(-6.4,6.4) + ax.set_xlabel('(x - pinball-front equilibrium centre) / D',fontsize=8) + ax.set_ylabel('(y - channel centre) / D',fontsize=8) + ax.tick_params(labelsize=8) + + +def main(): + ap=argparse.ArgumentParser(description=__doc__);ap.add_argument('directory',type=Path) + args=ap.parse_args();p=args.directory + case=json.loads((p/'case.json').read_text());scene=case['scene'] + cmap=plt.get_cmap('RdBu_r',21).copy();cmap.set_bad('#DDE2E5') + norm=BoundaryNorm(np.linspace(-3,3,22),21,clip=True) + prefix='SMOKE — plumbing only | ' if case['plan']['smoke'] else '' + roles=('target','zero','controlled') + data={};motion={} + for role in roles: + with np.load(p/role/'fields.npz') as z:data[role]={k:z[k] for k in z.files} + with np.load(p/role/'native.npz') as z: + motion[role]=z['rows'][:,1].copy()/40. + shape=data['target']['ux'].shape[-2:] + extent=((-0.5-1200)/40,(shape[1]-.5-1200)/40, + (-.5-510.5)/40,(shape[0]-.5-510.5)/40) + fig,axs=plt.subplots(3,1,figsize=(13,10),constrained_layout=True) + for ax,role in zip(axs,roles): + d=data[role];flags=d['common_fluid_all_boundaries'].astype(np.uint16) + w=vorticity(d['fixed_window_mean_ux'],d['fixed_window_mean_uy'],flags) + image=ax.imshow(w,origin='lower',extent=extent,cmap=cmap,norm=norm,interpolation='nearest') + decorate(ax,scene,role,0,mean=True) + label='no pinball reference' if role=='target' else ('zero rotation' if role=='zero' else ('constant rotations' if scene=='steady' else 'frozen PPO')) + ax.set_title(f'{role}: {label} | arithmetic fixed-window mean',fontsize=10) + fig.suptitle(prefix+f'{scene}: modern Re_D=50, m*=4, zero structural damping\n' + 'Gray: cells outside the common fluid support; mean fields are not instantaneous geometry',fontsize=12) + fig.colorbar(image,ax=list(axs),shrink=.65,label='vorticity D/U, fixed limits [-3,3]') + fig.savefig(p/'mean-flow.png',dpi=160);plt.close(fig) + slots=max(len(data[r]['ux']) for r in roles) + fig,axs=plt.subplots(3,slots,figsize=(max(8,3*slots),7),squeeze=False,constrained_layout=True) + for row,role in enumerate(roles): + d=data[role] + for col,ax in enumerate(axs[row]): + if col>=len(d['ux']):ax.axis('off');continue + w=vorticity(d['ux'][col],d['uy'][col],d['flags'][col]) + image=ax.imshow(w,origin='lower',extent=extent,cmap=cmap,norm=norm,interpolation='nearest') + step=int(d['selected_response_steps'][col]);y=float(motion[role][step-1]) + decorate(ax,scene,role,y) + ax.set_xlim(-3,14);ax.set_ylim(-3,3) + ax.set_xlabel('x/D relative to equilibrium front',fontsize=7) + ax.set_title(f'{role} | slot {col+1}\ntU/D={step*.005/40:.2f}',fontsize=8) + family={'karman':'independently phased diagnostic snapshots (different real times)', + 'steady':'late snapshots at matching real times', + 'taylor':'common event-relative snapshots at matching real times'}[scene] + fig.suptitle(prefix+f'{scene}: {family}\nFixed lab coordinates; same vorticity normalization/limits; source '+case['source'][:8],fontsize=12) + fig.colorbar(image,ax=list(axs.ravel()),shrink=.7,label='vorticity D/U [-3,3]') + fig.savefig(p/'selected-flow.png',dpi=160);plt.close(fig) + fig,ax=plt.subplots(figsize=(11,3.5),constrained_layout=True) + for role,color in (('zero','#5A6570'),('controlled','#147D73')): + y=motion[role];stride=max(1,len(y)//12000) + time=(np.arange(len(y))+1)*.005/40 + ax.plot(time[::stride],y[::stride],lw=.8,label=role,color=color) + ax.set(xlabel='response time tU/D',ylabel='transverse y/D',ylim=(-1,1)) + ax.legend();ax.grid(alpha=.2) + ax.set_title(prefix+scene+': independent paired releases; full retained duration') + fig.savefig(p/'motion.png',dpi=160);plt.close(fig) + print(json.dumps(dict(scene=scene,plots=['mean-flow.png','selected-flow.png','motion.png'], + vorticity_limits=[-3,3],normalization='curl(u) D/U; five-node fluid stencil', + mean_mask='intersection over collected boundary fields',display_interpolation='nearest'))) + + +if __name__=='__main__':main() diff --git a/research/fig2_collection/run.py b/research/fig2_collection/run.py new file mode 100644 index 0000000..f432bca --- /dev/null +++ b/research/fig2_collection/run.py @@ -0,0 +1,548 @@ +"""Research-only Fig. 2 acquisition on the modern moving plant. + +These paired releases are new experiments, not Legacy numerical reproductions. +Nothing here qualifies a steady state or automatically extends/retunes a run. +GPU dependencies are lazy so the collection contracts can be checked on CPU. +""" +from __future__ import annotations + +import argparse +import ast +import gc +import hashlib +import importlib.util +import json +import os +from pathlib import Path +import platform +import subprocess +import sys +import time +import types + +import numpy as np + +sys.dont_write_bytecode = True +ROOT = Path(__file__).resolve().parents[2] +D, U, NU, SI = 40., .005, .004, 3200 +BASELINE = '76873d911c4e2caba1bacc25aaa1495ea768453f' +CENTRES = ((1200., 510.5), (1252., 540.5), (1252., 480.5)) +PROBES = ((1600., 590.5), (1600., 510.5), (1600., 430.5)) +BIAS = np.array([0., -4., 4.]) +STEADY_SURFACE = np.array([0., -5.1, 5.1])*U +MIN_AVAILABLE = 16*1024**3 +ACQUISITION = Path('/home/frank14f/DynamisLab/src/drl_pinball/acquisition.py') +K_RELEASE = Path('/tmp/celeris-results-karman-viv-516975e/pair') +ASSETS = { + 'karman': (Path('/tmp/celeris-karman-viv-assets/model.zip'), + '148240c1dcb8b11d8e5a0a9e991f6874c4e799c9744707bf0fb32a8efeb0468d', + Path('/tmp/celeris-karman-viv-assets/norm.json'), + '534fda818c030a56960e95d0f3ac1e350d887d03db7f80f1d6453195a502b7f0'), + 'taylor': (Path('/home/frank14f/DynamisLab/models/old/vortex_taylor.zip'), + '58882721900f1eb0779f64aea0a2e233d307a94b640b9aaf33efa036db3c79e3', + Path('/home/frank14f/DynamisLab/src/SR_analysis/data/vortex/vortex_taylor/norm.json'), + 'bd38b4cf2b6da56b5cf378e646d673b97a606531c68163460313b330c9957b86'), +} +NATIVE_COLUMNS = ('response_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', 'omega_F', 'omega_U', 'omega_L') + + +def sha(path): + digest = hashlib.sha256() + with Path(path).open('rb') as stream: + for block in iter(lambda: stream.read(1024*1024), b''): + digest.update(block) + return digest.hexdigest() + + +def json_value(value): + if isinstance(value, dict): + return {str(k): json_value(v) for k, v in value.items()} + if isinstance(value, (tuple, list)): + return [json_value(v) for v in value] + if isinstance(value, np.ndarray): + return json_value(value.tolist()) + if isinstance(value, np.generic): + return value.item() + if isinstance(value, Path): + return str(value.resolve()) + return value + + +def write_json(path, value): + Path(path).write_text(json.dumps(json_value(value), indent=2, allow_nan=False)+'\n') + + +def load_file(name, path): + spec = importlib.util.spec_from_file_location(name, path) + module = importlib.util.module_from_spec(spec) + sys.modules[name] = module + spec.loader.exec_module(module) + return module + + +def load_helpers(path=ACQUISITION): + """Load read-only CPU helpers under an isolated package; no original imports. + + Empty namespace packages avoid executing application __init__ modules. The + acquisition module and its pure DTW dependency execute directly from source. + Event functions are AST-selected verbatim to avoid acquire.py's heavy imports; + offsets come from that acquisition module, never a copied/fabricated selector. + """ + path = Path(path) + package = '_fig2_readonly_drl' + for suffix, directory in (('', path.parent), ('.legacy_test', path.parent/'legacy_test'), + ('.legacy_test.core', path.parent/'legacy_test/core')): + module = types.ModuleType(package+suffix) + module.__path__ = [str(directory)] + sys.modules[package+suffix] = module + acquisition = load_file(package+'.acquisition', path) + event_path = path.parent/'legacy_test/acquire.py' + source = ast.parse(event_path.read_text(), filename=str(event_path)) + names = {'_vortex_event_channel', '_select_vortex_event'} + nodes = [node for node in source.body if isinstance(node, ast.FunctionDef) and node.name in names] + if {node.name for node in nodes} != names: + raise ValueError('authoritative event selector functions missing') + namespace = dict(np=np, VORTEX_EVENT_OFFSETS=acquisition.VORTEX_EVENT_OFFSETS) + exec(compile(ast.Module(body=nodes, type_ignores=[]), str(event_path), 'exec'), namespace) + acquisition.select_vortex_event = namespace['_select_vortex_event'] + return acquisition + + +def schedule(scene, smoke=False): + if scene not in ('karman', 'steady', 'taylor'): + raise ValueError('unknown scene') + intervals = 150 if scene == 'taylor' else 400 + release = {'karman': 4576000, 'steady': 2048000, 'taylor': 1024000}[scene] + return dict(si=32 if smoke else SI, intervals=3 if smoke else intervals, + field_start=0 if smoke or scene == 'taylor' else 240, + release_clock=(64 if scene == 'taylor' else 32) if smoke else release, + smoke=bool(smoke)) + + +def available_memory(path=Path('/proc/meminfo')): + for line in Path(path).read_text().splitlines(): + if line.startswith('MemAvailable:'): + return int(line.split()[1])*1024 + raise RuntimeError('MemAvailable unavailable; RAM admission cannot be verified') + + +def admit_memory(smoke=False): + if smoke: + return None + available = available_memory() + if available < MIN_AVAILABLE: + raise MemoryError(f'MemAvailable {available} < required {MIN_AVAILABLE}') + return available + + +def wall_surface(scene, role, action=None): + if role == 'target' or role == 'zero': + return np.zeros(3, np.float64) # exactly zero physical rotation + if scene == 'steady': + return STEADY_SURFACE.copy() + action = np.asarray(action) + if action.shape != (3,) or not np.isfinite(action).all() or (np.abs(action) > 1).any(): + raise ValueError('invalid deterministic policy action') + return ((8. if scene == 'karman' else 4.)*action+BIAS)*U + + +def check_structure(y, v, previous_y): + if not np.isfinite([y, v]).all() or abs(y) > D or abs(v) > .025 or abs(y-previous_y) > .025: + raise RuntimeError('nonfinite structure, travel > D or speed > .025; no retuning') + + +def clean_field(fields, flags): + """Canonical modern bit-0 fluid mask; hard-fail before nonfluid zero fill.""" + ux, uy = np.asarray(fields['ux']), np.asarray(fields['uy']) + rho = np.asarray(fields['rho']) + if ux.ndim != 2 or uy.shape != ux.shape or rho.shape != ux.shape: + raise ValueError('canonical shared (ny,nx) macroscopic arrays required') + flags = np.asarray(flags) + if flags.size != ux.size or flags.dtype != np.uint16: + raise ValueError('modern uint16 flags matching field required') + flags = flags.reshape(ux.shape) + fluid = (flags & 1) != 0 + if any(not np.isfinite(value[fluid]).all() for value in (rho, ux, uy)): + raise RuntimeError('nonfinite FLUID macroscopic field') + result = np.zeros((2, *ux.shape), np.float32) + result[0, fluid], result[1, fluid] = ux[fluid], uy[fluid] + return result, flags.copy(), fluid + + +def select_fields(scene, plan, helpers, times, sensors, cache, event=None): + """Fixed-window means remain separate from independently phased diagnostics.""" + mean = helpers.accumulate_mean_fields(cache[:, 0], cache[:, 1], np.arange(len(cache))) + result = dict(fixed_window_mean_ux=mean['mean_ux'], fixed_window_mean_uy=mean['mean_uy']) + metadata = dict(fixed_window_count=len(cache), periodic_phase_status='not_requested') + indices = np.array([len(cache)-1], np.int64) + if plan['smoke']: + metadata['selection'] = 'smoke_last_snapshot_not_physical' + elif scene == 'taylor': + if event is None: + event = helpers.select_vortex_event(sensors, 'vortex_taylor') + indices = np.asarray(event['field_indices'], np.int64) + if (indices < 0).any() or (indices >= len(cache)).any(): + raise ValueError('authoritative event offsets exceed collected fields') + metadata.update(selection='source_vortex_event', event=event) + elif scene == 'karman': + try: + phase = helpers.assign_periodic_phase(times, sensors, minimum_crossing_interval=.5*30*SI) + complete = helpers.complete_cycle_field_indices(times, phase['crossing_times']) + chosen = helpers.select_phase_fields(times, phase['crossing_times']) + cycle_mean = helpers.accumulate_mean_fields(cache[:, 0], cache[:, 1], complete['field_indices']) + result.update(cycle_mean_ux=cycle_mean['mean_ux'], cycle_mean_uy=cycle_mean['mean_uy']) + indices = chosen['field_indices'] + metadata.update(periodic_phase_status='available_diagnostic_only', phase_selection=chosen, + complete_cycles=complete, crossing_times=phase['crossing_times'], + rejected_crossing_times=phase['rejected_crossing_times']) + except ValueError as exc: + metadata.update(periodic_phase_status='insufficient_crossings', phase_reason=str(exc)) + metadata['selection'] = 'independent_center_probe_uy_diagnostic' + else: + metadata.update(selection='last_snapshot', periodic_phase_status='not_assumed_steady') + result.update(ux=cache[indices, 0], uy=cache[indices, 1], selected_index=indices, + selected_response_steps=times[indices]) + return result, metadata, event + + +def add_scene(sim, scene, *, moving=False, pinball=True): + from CelerisLab.body.motion import Translation + if scene == 'karman': + sim.add_body('circle', (400., 510.5), 40.) + probes = [sim.add_body('sensor', point, 10.) for point in PROBES] + envelope = Translation((0., 0., -D, D), .025, shared=True) if moving else None + bodies = [sim.add_body('circle', point, 20., translation=envelope) for point in CENTRES] if pinball else [] + return probes, bodies + + +def inject_taylor(sim): + from CelerisLab.lbm.initializers import add_vortex + clock = int(sim.stepper.step_count) + with sim.ctx.activate(): + sim.stream.synchronize() + add_vortex(sim.field, (600., 510.5), 80., .000075, 'taylor') + sim.stream.synchronize() + if sim.stepper.step_count != clock: + raise RuntimeError('vortex initializer advanced response clock') + return dict(absolute_clock=clock, response_clock=0, center=[600., 510.5], radius=80., + strength=.000075, type='taylor', count=1, + operator='unmodified add_vortex: all interior reset to equilibrium; both buffers uploaded', + impulse='full equilibrium reset discards pre-existing non-equilibrium on all interior cells') + + +def make_release(scene, plan, cfgpath, base, out, started, cap, release_dir): + from CelerisLab import Simulation + if scene == 'karman' and not plan['smoke']: + seed, force = base.saved_release(release_dir) + return seed, force, dict(kind='verified_saved_static_release', path=release_dir, + release_npz_sha256=sha(release_dir/'release.npz')) + with Simulation(str(cfgpath)) as sim: + _, bodies = add_scene(sim, scene, pinball=scene != 'taylor') + sim.initialize() + warm = 32 if plan['smoke'] else (512000 if scene == 'taylor' else 2048000) + if scene == 'taylor': + base.warmup(sim, warm, started, cap) + bodies = [sim.add_body('circle', point, 20.) for point in CENTRES] + sim.sync_bodies() + base.warmup(sim, warm-1, started, cap) + sim.run(1, zero_obs=True, sync_obs=True) + force = sum(float(sim.read_force(body)[1]) for body in bodies) + seed = base.capture_seed(sim) + if seed[3] != plan['release_clock'] or not np.isfinite(force): + raise RuntimeError('unexpected static release clock/force') + np.savez_compressed(out/'release.npz', current=seed[0], spare=seed[1], flags=seed[2], + clock=seed[3], force=force) + return seed, force, dict(kind='fresh_static_full_two_buffer_release', + release_npz_sha256=sha(out/'release.npz')) + + +def compress_native(path, rows, completed): + """Bounded disk memmap; preserve raw file until parent elects removal.""" + rows.flush() + compressed = path.with_suffix('.npz') + np.savez_compressed(compressed, rows=rows[:completed]) + with np.load(compressed) as data: + retained = data['rows'] + if retained.shape != (completed, len(NATIVE_COLUMNS)): + raise RuntimeError('compressed native shape mismatch') + for start in range(0, completed, 10000): + end = min(start+10000, completed) + if not np.array_equal(retained[start:end], rows[start:end]): + raise RuntimeError('compressed native verification failed') + return dict(raw_path=path.name, compressed_path=compressed.name, verified=True, + compressed_sha256=sha(compressed), raw_preserved_for_parent_cleanup=True) + + +def run_role(scene, role, plan, cfgpath, base, helpers, metrics, seed, initial_force, + model, norm, out, started, cap, event): + from CelerisLab import Simulation + from CelerisLab.body.motion import MovingTables + directory = out/role + directory.mkdir() + steps = plan['si']*plan['intervals'] + rows = np.lib.format.open_memmap(directory/'native.npy', mode='w+', dtype=np.float64, + shape=(steps, len(NATIVE_COLUMNS))) + signals = {key: np.zeros((plan['intervals'], width), np.float32) for key, width in + [('sensor_area_mean', 6), ('sensors', 6), ('forces', 6), ('actions', 3), + ('observations', 12), ('policy_input', 12), ('raw_observation', 12)]} + structure = base.Structure(4., 2.5e-5, 0., D) + columns = {name: i for i, name in enumerate(MovingTables.DIAG)} + completed = completed_si = 0 + duration = steps*U/D + report = dict(status='running', role=role, target_is_free_body=False if role == 'target' else None, + columns=NATIVE_COLUMNS, + motion_window=dict(start_tU_D=max(0., duration-80.), end_tU_D=duration, + interpretation='smoke_short_not_physical' if plan['smoke'] else + ('entire_60DU_event_not_settled' if scene == 'taylor' else 'declared_last_80DU')), + structure=None if role == 'target' else dict(mass=structure.mass, fn=structure.frequency, + mass_ratio=4., zeta=0., y0=0., v0=0.)) + cache, common_fluid, selected_flags, roi_cache = None, None, {}, None + obs, force = np.zeros(12, np.float32), initial_force + try: + report['MemAvailable_admission_bytes'] = admit_memory(plan['smoke']) + with Simulation(str(cfgpath)) as sim: + probes, bodies = add_scene(sim, scene, moving=role != 'target', pinball=role != 'target') + sim.initialize() + if role == 'target': + base.warmup(sim, plan['release_clock'], started, cap) + report['initialization'] = 'own no-pinball initialization at matching absolute clock' + else: + report['release_transfer'] = base.install_seed(sim, seed) + if sim.stepper.step_count != plan['release_clock']: + raise RuntimeError('role release clocks differ') + if scene == 'taylor': + report['injection'] = inject_taylor(sim) + count = plan['intervals']-plan['field_start'] + shape = (sim.field.ny, sim.field.nx) + cache = np.empty((count, 2, *shape), np.float32) + common_fluid = np.ones(shape, bool) + roi = metrics.roi_slices(shape, D=D, xs=1600., y0=510.5, + downstream_solid_extent=1272.) + # Only compact ROI histories survive a role; full-grid histories never reach disk. + roi_shape = cache[(slice(None), slice(None), *roi)].shape + roi_cache = np.empty(roi_shape, np.float32) + retain = set(event['field_indices'].tolist()) if event is not None else set() + if scene == 'taylor' and event is None and not plan['smoke']: + # Authoritative selector searches 38..87; any imported offset can be selected. + retain = {i+int(offset) for i in range(38, 88) for offset in helpers.VORTEX_EVENT_OFFSETS} + # K phase unknown until the final signal: retain flags for every boundary in RAM. + flag_cache = np.empty((count, *shape), np.uint16) if scene == 'karman' else None + for interval in range(plan['intervals']): + base.check_time(started, cap) + action = np.zeros(3, np.float32) + signals['policy_input'][interval] = obs + if role == 'controlled' and scene != 'steady': + action = np.asarray(model.predict(obs, deterministic=True)[0], np.float32) + target, effective = wall_surface(scene, role, action), np.zeros(3) + sensor_sum, force_sum = np.zeros(6), np.zeros(6) + for native in range(plan['si']): + if native % 1000 == 0: + base.check_time(started, cap) + y, v = (0., 0.) if role == 'target' else structure.candidate(force) + check_structure(y, v, structure.y) + effective = .9**.25*effective + (1.-.9**.25)*target + omega = effective/20. + for body, w in zip(bodies, omega): + sim.set_body(body, omega=float(w), offset=(0., y), vy=y-structure.y) + sim.run(1, zero_obs=True, sync_obs=True) + if bodies: + diag = sim.moving_diagnostics() + cme, wen = diag[:, columns['cme_y']], diag[:, columns['wen_y']] + node = diag[:, columns['cover_jy']]-diag[:, columns['uncover_jy']] + counters = [diag[:, columns[key]].sum() for key in ('eq_only', 'no_donor', 'crossings')] + packed = np.array([value for body in bodies for value in sim.read_force(body)[:2]]) + else: + cme, wen, node, counters, packed = np.zeros(3), np.zeros(3), np.zeros(3), [0., 0., 0.], np.zeros(6) + sensors = np.array([value for probe in probes for value in sim.read_sensor(probe)[:2]]) + row = np.array([completed+1, y, v, *cme, *wen, *node, *counters, *omega]) + if not np.isfinite(row).all() or counters[1] != 0 or not np.isfinite([packed, sensors]).all(): + raise RuntimeError('nonfinite native diagnostics/observations or no_donor') + rows[completed] = row + completed += 1 + structure.y, structure.v = y, v + force = float(cme.sum(dtype=np.float64)+wen.sum(dtype=np.float64)) + sensor_sum += sensors + force_sum += packed + sensors, packed = sensor_sum/plan['si'], force_sum/plan['si'] + raw = np.r_[sensors*156., packed*2.].astype(np.float32) + if norm is not None: + raw, obs, clipped = base.observation(sensors, packed, 2., norm) + report['clipped_channels'] = report.get('clipped_channels', 0)+clipped + for key, value in [('sensor_area_mean', sensors), ('sensors', sensors*156.), + ('forces', packed*2.), ('actions', action), ('observations', obs), + ('raw_observation', raw)]: + signals[key][interval] = value + if interval >= plan['field_start']: + index = interval-plan['field_start'] + clean, flags, fluid = clean_field(sim.get_macroscopic(), sim.get_flags()) + if not fluid[roi].all(): + raise RuntimeError('prescribed downstream ROI contains nonfluid cells') + cache[index] = clean + common_fluid &= fluid + roi_cache[index] = clean[(slice(None), *roi)] + if flag_cache is not None: + flag_cache[index] = flags + elif index in retain or index == count-1: + selected_flags[index] = flags + completed_si += 1 + if completed_si % 10 == 0 or completed_si == plan['intervals']: + rows.flush() + print(json.dumps(dict(event='progress', scene=scene, role=role, interval=completed_si, + absolute_clock=sim.stepper.step_count, y_D=y/D)), flush=True) + times = np.arange(plan['field_start']+1, plan['intervals']+1, dtype=np.int64)*plan['si'] + selected, selection, event = select_fields(scene, plan, helpers, times, + signals['sensors'][plan['field_start']:], cache, event) + indices = selected['selected_index'] + selected['flags'] = flag_cache[indices] if flag_cache is not None else np.stack([selected_flags[int(i)] for i in indices]) + selected['common_fluid_all_boundaries'] = common_fluid + selected['absolute_steps'] = selected['selected_response_steps']+plan['release_clock'] + np.savez_compressed(directory/'fields.npz', **selected) + report.update(selection=selection, field_response_steps=times, + field_absolute_steps=times+plan['release_clock'], + field_window_tU_D=[int(times[0])*U/D, int(times[-1])*U/D], + status='smoke_completed_not_physical' if plan['smoke'] else 'duration_completed_not_qualified') + roi_mask = common_fluid[roi].copy() + report['motion_summary'] = metrics.motion_summary(rows[max(0, completed-640000):completed, 1]/D) + return dict(roi=roi_cache, mask=roi_mask, times=times, event=event) + except BaseException as exc: + report['status'] = f'error: {type(exc).__name__}: {exc}' + raise + finally: + report.update(completed_native_steps=completed, completed_intervals=completed_si) + np.savez_compressed(directory/'signals.npz', response_steps=np.arange(1, completed_si+1)*plan['si'], + absolute_steps=plan['release_clock']+np.arange(1, completed_si+1)*plan['si'], + **{key: value[:completed_si] for key, value in signals.items()}) + if completed: + report['native_archive'] = compress_native(directory/'native.npy', rows, completed) + write_json(directory/'report.json', report) + del cache, rows + gc.collect() + + +def score_roles(roles, metrics): + target = roles['target'] + common = np.logical_and.reduce([role['mask'] for role in roles.values()]) + if not common.all(): + raise ValueError('prescribed downstream ROI must be entirely common fluid') + result = dict(definition='sqrt(mean common ROI (dux^2+duy^2)/U^2)', + alignment='same absolute native clock; no lag/phase alignment', + roi_xy_inclusive=[1360, 2160, 311, 710], support_cells=int(common.sum()), + support_coverage_fraction=1.) + for name in ('zero', 'controlled'): + role = roles[name] + if not np.array_equal(role['times'], target['times']): + raise ValueError('unmatched response clocks') + errors = np.array([metrics.field_error(*a, *b, U, common) for a, b in zip(role['roi'], target['roi'])]) + # FP64 reduction without duplicating entire history. + mean = [a.mean(axis=0, dtype=np.float64) for a in (role['roi'], target['roi'])] + result[name] = dict(same_clock_instantaneous_E=errors, + mean_instantaneous_E=float(errors.mean()), + E_of_fixed_window_mean=metrics.field_error(*mean[0], *mean[1], U, common)) + return result, common + + +def source_manifest(base, helpers): + files = [Path(__file__), ROOT/'research/fig2_materials/metrics.py', Path(base.__file__), + ROOT/'research/steady_pinball/run.py', ROOT/'research/steady_pinball/config.json', + Path(helpers.__file__), Path(helpers.__file__).parent/'legacy_test/acquire.py', + Path(helpers.__file__).parent/'legacy_test/core/dtw_metrics.py'] + files += sorted((ROOT/'src/CelerisLab').rglob('*.py')) + for suffix in ('*.cu', '*.cuh', '*.h'): + files += sorted((ROOT/'src/CelerisLab').rglob(suffix)) + return {str(path.resolve()): sha(path) for path in files} + + +def main(argv=None): + ap = argparse.ArgumentParser(description=__doc__) + ap.add_argument('--scene', choices=('karman', 'steady', 'taylor'), required=True) + ap.add_argument('--out', type=Path, required=True) + ap.add_argument('--smoke', action='store_true', help='32-step warmups; 3 SI of 32; no physical score') + ap.add_argument('--D', type=int, choices=(40,), default=40) + ap.add_argument('--model', type=Path) + ap.add_argument('--norm', type=Path) + ap.add_argument('--release-dir', type=Path, default=K_RELEASE) + ap.add_argument('--acquisition', type=Path, default=ACQUISITION) + ap.add_argument('--wall-seconds', type=float, default=43200.) + args = ap.parse_args(argv) + if not np.isfinite(args.wall_seconds) or args.wall_seconds <= 0: + ap.error('positive finite wall cap required') + if args.scene == 'steady' and (args.model or args.norm): + ap.error('steady constant control has no PPO or normalizer') + out = args.out.resolve() + if out == ROOT or ROOT in out.parents: + ap.error('fresh output must be outside checkout') + sys.path.insert(0, str(ROOT)) + from research.karman_viv import run as base + metrics = load_file('_fig2_metrics', ROOT/'research/fig2_materials/metrics.py') + helpers = load_helpers(args.acquisition) + model = norm = None + inputs = {} + if args.scene != 'steady': + model_path, model_sha, norm_path, norm_sha = ASSETS[args.scene] + model_path, norm_path = args.model or model_path, args.norm or norm_path + if sha(model_path) != model_sha or sha(norm_path) != norm_sha: + ap.error('model/normalizer hashes differ from frozen inputs') + norm = json.loads(norm_path.read_text()) + inputs = dict(model_path=model_path, model_sha256=model_sha, + norm_path=norm_path, norm_sha256=norm_sha) + model = base.load_policy(model_path) + sources = source_manifest(base, helpers) + source = subprocess.check_output(['git', '-C', str(ROOT), 'rev-parse', 'HEAD'], text=True).strip() + dirty = subprocess.check_output(['git', '-C', str(ROOT), 'status', '--porcelain'], text=True) + plan = schedule(args.scene, args.smoke) + admit_memory(args.smoke) + out.mkdir(parents=True, exist_ok=False) + cfgpath = out/'config.json' + write_json(cfgpath, base.config(40)) + receipt = dict(scene=args.scene, plan=plan, D=D, U=U, nu=NU, actual_Re_D=50., + historical_case='karman_re100' if args.scene == 'karman' else args.scene, + expected_baseline=BASELINE, source=source, git_dirty=bool(dirty), git_status=dirty, + source_sha256=sources, inputs=inputs, centres=CENTRES, probes=PROBES, + environment=dict(python=sys.version, platform=platform.platform(), numpy=np.__version__, + conda_env=os.environ.get('CONDA_DEFAULT_ENV'), + cuda_visible_devices=os.environ.get('CUDA_VISIBLE_DEVICES')), + acquisition_event_offsets=list(helpers.VORTEX_EVENT_OFFSETS), + event_selector='source _select_vortex_event AST verbatim; actual acquisition offsets', + structure=dict(mass_ratio=4., fn_native=2.5e-5, fn_D_U=.2, zeta=0., travel_D=1., speed=.025), + rhs='previous native-step sum CME_y+Wen_y; initial force last static native step', + policy='SI mean packed CME *2; area-mean sensors *156; frozen norm; first input zeros', + control='steady constant (0,-5.1,5.1)*U; K scale8 or Taylor scale4, bias(0,-4,4)*U', + ema=dict(alpha=.9**.25, reset='every SI'), + storage='one full-grid float32 RAM cache at a time; compact ROI histories in RAM; no full history field disk', + limitation='Modern moving-plant experiment; no Legacy numerical equivalence or steady/convergence claim') + write_json(out/'case.json', receipt) + started, status = time.monotonic(), 'running' + try: + seed, force, release = make_release(args.scene, plan, cfgpath, base, out, started, + args.wall_seconds, args.release_dir) + receipt['release'] = release + receipt['release'].update(clock=int(seed[3]), initial_force=force, + current_sha256=base.array_sha(seed[0]), + spare_sha256=base.array_sha(seed[1]), flags_sha256=base.array_sha(seed[2])) + write_json(out/'case.json', receipt) + roles, event = {}, None + order = ('controlled', 'zero', 'target') if args.scene == 'taylor' else ('target', 'zero', 'controlled') + for role in order: + roles[role] = run_role(args.scene, role, plan, cfgpath, base, helpers, metrics, + seed, force, model, norm, out, started, args.wall_seconds, event) + if args.scene == 'taylor': + event = roles[role]['event'] + if not args.smoke: + scores, support = score_roles(roles, metrics) + write_json(out/'same_clock_errors.json', scores) + np.savez_compressed(out/'common_roi_support.npz', support=support) + if source_manifest(base, helpers) != sources: + raise RuntimeError('source files changed during collection') + status = 'smoke_completed_not_physical' if args.smoke else 'three_roles_completed_not_qualified' + except BaseException as exc: + status = f'error: {type(exc).__name__}: {exc}' + raise + finally: + write_json(out/'terminal.json', dict(status=status, elapsed_s=time.monotonic()-started, + physical_score_claim=False if args.smoke else None)) + + +if __name__ == '__main__': + main() diff --git a/research/fig2_materials/metrics.py b/research/fig2_materials/metrics.py new file mode 100644 index 0000000..84d0f57 --- /dev/null +++ b/research/fig2_materials/metrics.py @@ -0,0 +1,94 @@ +"""Pure NumPy Fig2 primitives; no solver, phase selection, or probe DTW. + +Adapted formula and inclusive ROI bounds from DynamisLab +src/drl_pinball/eval/wake_l2.py::mean_field_metric, roi_slices, inspected +2026-11-02 at HEAD 609116ec26f19766f5e9ce2b0df427d4e24faae8 +(working-tree source, not asserted clean). The explicit common-fluid/finite +mask and motion summaries are additions. legacy_test/metrics.py::cycle_similarity +was inspected only: its six probe channels are not a field recovery metric. +""" +import math +import numpy as np + + +def roi_slices(shape, *, D, xs, y0=None, downstream_solid_extent): + """Lab-fixed ROI: sensor-relative x/D=-6..14, y/D=-5..5, inclusive. + + Equivalent to pinball-front-relative x/D=4..24 when xs=front+10D. + Supply actual geometry, e.g. D=40, xs=1600, y0=510.5. + """ + ny, nx = shape + y0 = (ny - 1) / 2 if y0 is None else float(y0) + if D <= 0 or not all(np.isfinite(v) for v in (D, xs, y0, downstream_solid_extent)): + raise ValueError("invalid geometry") + x0, x1 = math.ceil(xs - 6 * D), math.floor(xs + 14 * D) + lo, hi = math.ceil(y0 - 5 * D), math.floor(y0 + 5 * D) + if x0 < 0 or lo < 0 or x1 >= nx or hi >= ny: + raise ValueError("ROI outside grid") + if x0 <= downstream_solid_extent: + raise ValueError("ROI not downstream of solid extent") + return slice(lo, hi + 1), slice(x0, x1 + 1) + + +def common_fluid_roi_mask(fluid_masks, fields, roi): + """Intersect explicit boolean fluid masks, ROI, and finite 2-D fields. + + Pass every compared role/component/sample to freeze a common support over + time. Zero velocity never implies solid. Masks must share the lab grid. + """ + fields = [np.asarray(a, dtype=float) for a in fields] + masks = [np.asarray(a) for a in fluid_masks] + if not fields or not masks: + raise ValueError("fields and explicit fluid masks required") + shape = fields[0].shape + if len(shape) != 2 or any(a.shape != shape for a in fields + masks): + raise ValueError("equal-shape 2-D arrays required") + if any(a.dtype != np.dtype(bool) for a in masks): + raise ValueError("fluid masks must be boolean") + support = np.zeros(shape, dtype=bool) + ys, xs = roi + for s, n in zip((ys, xs), shape): + if not isinstance(s, slice) or s.step not in (None, 1) or not 0 <= (s.start or 0) < (n if s.stop is None else s.stop) <= n: + raise ValueError("invalid ROI slices") + support[ys, xs] = True + for a in masks: + support &= a + for a in fields: + support &= np.isfinite(a) + if not support.any(): + raise ValueError("empty common finite fluid ROI") + return support + + +def field_error(rux, ruy, tux, tuy, u0, support): + """Eq1: sqrt(mean_support((dux**2 + duy**2)/U0**2)). + + For E(mean field), pass arithmetic mean fields from matched windows. + Mean instantaneous E requires separate calls per snapshot and then a mean; + these operations are not interchangeable. No target-energy denominator. + """ + a = [np.asarray(v, dtype=float) for v in (rux, ruy, tux, tuy)] + support = np.asarray(support) + if a[0].ndim != 2 or any(v.shape != a[0].shape for v in a) or support.shape != a[0].shape or support.dtype != np.dtype(bool): + raise ValueError("equal-shape 2-D fields and boolean support required") + if not np.isfinite(u0) or u0 <= 0 or not support.any(): + raise ValueError("positive finite U0 and nonempty support required") + if any(not np.isfinite(v[support]).all() for v in a): + raise ValueError("nonfinite field on fixed support") + du = (a[0][support] - a[2][support]) / u0 + dv = (a[1][support] - a[3][support]) / u0 + return float(np.sqrt(np.mean(du * du + dv * dv))) + + +def motion_summary(displacement): + """Arithmetic mean, centered RMS, half peak-to-peak along time axis 0. + + Input is displacement in caller-declared units; no implicit D scaling. + Use the same retained window as field evaluation. No sinusoid assumption. + """ + a = np.asarray(displacement, dtype=float) + if a.ndim == 0 or a.shape[0] == 0 or not np.isfinite(a).all(): + raise ValueError("nonempty finite time samples required") + mean = a.mean(axis=0) + return {"mean": mean, "centered_rms": np.sqrt(np.mean((a - mean) ** 2, axis=0)), + "half_range": (a.max(axis=0) - a.min(axis=0)) / 2} diff --git a/tests/unit/test_fig2_analysis.py b/tests/unit/test_fig2_analysis.py new file mode 100644 index 0000000..2a1a619 --- /dev/null +++ b/tests/unit/test_fig2_analysis.py @@ -0,0 +1,36 @@ +"""Independent analytic data checks for offline estimands, no CFD context.""" +import unittest +import numpy as np +from research.fig2_collection.analyze import frequency, native_summary, compare_selected + + +class AnalysisTests(unittest.TestCase): + def test_frequency_units_and_no_dc(self): + dt=.005/40;n=640000 + t=np.arange(n)*dt + result=frequency(4+np.sin(2*np.pi*.2*t),dt) + self.assertAlmostEqual(result['f_D_U'],.2) + self.assertAlmostEqual(result['bin_width_D_U'],1/80) + self.assertIsNone(frequency(np.ones(20),dt)) + + def test_motion_and_native_rhs_are_distinct_from_node_impulse(self): + a=np.zeros((10,18));a[:,0]=np.arange(1,11);a[:,1]=40*np.tile([.5,-.5],5) + a[:,3]=.001;a[:,6]=.002;a[:,9]=999 + result=native_summary(a,'taylor') + self.assertEqual(result['motion']['centered_rms'],.5) + self.assertEqual(result['motion']['half_range'],.5) + self.assertAlmostEqual(result['structural_RHS_Cy_RMS'],6.) + self.assertEqual(result['window_D_U'],[0.,10*.005/40]) + + def test_selected_vector_error_uses_common_roi_and_actual_flags(self): + shape=(1,712,2162) + target=dict(ux=np.zeros(shape,np.float32),uy=np.zeros(shape,np.float32),flags=np.ones(shape,np.uint16)) + current={k:v.copy() for k,v in target.items()} + current['ux'][:,311:711,1360:2161]=.003 + current['uy'][:,311:711,1360:2161]=.004 + self.assertAlmostEqual(compare_selected(target,current)['aggregate_E'],1.,places=6) + current['flags'][0,400,1500]=2 + with self.assertRaises(ValueError):compare_selected(target,current) + + +if __name__=='__main__':unittest.main() diff --git a/tests/unit/test_fig2_collection.py b/tests/unit/test_fig2_collection.py new file mode 100644 index 0000000..9f1036f --- /dev/null +++ b/tests/unit/test_fig2_collection.py @@ -0,0 +1,225 @@ +"""CPU contracts for Fig. 2 acquisition; no CUDA or simulation import required.""" +import importlib.util +import json +from pathlib import Path +import sys +import tempfile +import types +import unittest +from unittest.mock import patch + +import numpy as np + +ROOT = Path(__file__).resolve().parents[2] +SPEC = importlib.util.spec_from_file_location('fig2_collection_tested', ROOT/'research/fig2_collection/run.py') +driver = importlib.util.module_from_spec(SPEC) +SPEC.loader.exec_module(driver) +METRIC_SPEC = importlib.util.spec_from_file_location('fig2_collection_metrics', ROOT/'research/fig2_materials/metrics.py') +metrics = importlib.util.module_from_spec(METRIC_SPEC) +METRIC_SPEC.loader.exec_module(metrics) + + +class CollectionContracts(unittest.TestCase): + def test_frozen_clocks_and_windows(self): + for scene, release, count, first in [('karman', 4576000, 400, 240), + ('steady', 2048000, 400, 240), + ('taylor', 1024000, 150, 0)]: + plan = driver.schedule(scene) + self.assertEqual((plan['release_clock'], plan['intervals'], plan['field_start']), + (release, count, first)) + self.assertEqual(plan['si'], 3200) + self.assertEqual((400-240)*3200*driver.U/driver.D, 64.) + self.assertEqual(400*3200*driver.U/driver.D, 160.) + self.assertEqual(150*3200*driver.U/driver.D, 60.) + + def test_smoke_has_no_physical_timing(self): + for scene in ('karman', 'steady', 'taylor'): + plan = driver.schedule(scene, True) + self.assertEqual((plan['si'], plan['intervals'], plan['field_start']), (32, 3, 0)) + self.assertEqual(plan['release_clock'], 64 if scene == 'taylor' else 32) + + def test_zero_is_exact_physical_zero_despite_bias(self): + for scene in ('karman', 'steady', 'taylor'): + np.testing.assert_array_equal(driver.wall_surface(scene, 'zero'), np.zeros(3)) + np.testing.assert_array_equal(driver.wall_surface(scene, 'target'), np.zeros(3)) + np.testing.assert_allclose(driver.wall_surface('karman', 'controlled', [.5, -.5, 1]), + [.02, -.04, .06]) + np.testing.assert_allclose(driver.wall_surface('taylor', 'controlled', [.5, -.5, 1]), + [.01, -.03, .04]) + np.testing.assert_allclose(driver.wall_surface('steady', 'controlled'), [0, -.0255, .0255]) + with self.assertRaises(ValueError): + driver.wall_surface('karman', 'controlled', [0, 0, np.nan]) + + def test_structure_hard_stops(self): + driver.check_structure(40., 0., 40.) + for y, v, old in [(40.01, 0, 40), (0, .026, 0), (.026, 0, 0), (np.nan, 0, 0)]: + with self.subTest(y=y, v=v), self.assertRaises(RuntimeError): + driver.check_structure(y, v, old) + + def test_fluid_nonfinite_fails_nonfluid_only_zero_filled(self): + flags = np.array([[1, 257, 2]], np.uint16) + fields = dict(rho=np.array([[1., 1., np.nan]]), ux=np.array([[0., .1, np.nan]]), + uy=np.array([[0., -.1, np.inf]])) + clean, retained, fluid = driver.clean_field(fields, flags) + np.testing.assert_array_equal(fluid, [[True, True, False]]) + np.testing.assert_array_equal(retained, flags) + self.assertTrue(np.isfinite(clean).all()) + self.assertEqual(clean[0, 0, 0], 0.) # zero flow remains fluid + fields['ux'][0, 1] = np.nan + with self.assertRaisesRegex(RuntimeError, 'FLUID'): + driver.clean_field(fields, flags) + + def test_memory_admission_uses_available_not_free(self): + with tempfile.TemporaryDirectory() as directory: + path = Path(directory)/'meminfo' + path.write_text('MemFree: 1 kB\nMemAvailable: 16777216 kB\n') + self.assertEqual(driver.available_memory(path), driver.MIN_AVAILABLE) + with patch.object(driver, 'available_memory', return_value=driver.MIN_AVAILABLE-1): + with self.assertRaises(MemoryError): + driver.admit_memory() + self.assertIsNone(driver.admit_memory(True)) + + def test_roi_bounds_match_full_downstream_support(self): + roi = metrics.roi_slices((1022, 2559), D=40, xs=1600, y0=510.5, + downstream_solid_extent=1272.) + self.assertEqual(roi, (slice(311, 711), slice(1360, 2161))) + + def test_same_clock_error_and_mean_are_distinct(self): + shape = (2, 2, 2, 3) + target = np.zeros(shape, np.float32) + changed = target.copy() + changed[0, 0], changed[1, 0] = driver.U, -driver.U + role = lambda fields: dict(roi=fields, mask=np.ones((2, 3), bool), times=np.array([3200, 6400])) + roles = dict(target=role(target), zero=role(changed), controlled=role(target)) + result, mask = driver.score_roles(roles, metrics) + self.assertAlmostEqual(result['zero']['mean_instantaneous_E'], 1., places=7) + self.assertEqual(result['zero']['E_of_fixed_window_mean'], 0.) + roles['zero']['times'] = np.array([6400, 9600]) + with self.assertRaisesRegex(ValueError, 'clocks'): + driver.score_roles(roles, metrics) + roles['zero']['times'] = roles['target']['times'] + roles['zero']['mask'][0, 0] = False + with self.assertRaisesRegex(ValueError, 'entirely common fluid'): + driver.score_roles(roles, metrics) + + def test_native_archive_verified_raw_retained(self): + with tempfile.TemporaryDirectory() as directory: + path = Path(directory)/'native.npy' + rows = np.lib.format.open_memmap(path, mode='w+', dtype=np.float64, + shape=(9, len(driver.NATIVE_COLUMNS))) + rows[:] = np.arange(rows.size).reshape(rows.shape) + report = driver.compress_native(path, rows, 7) + self.assertTrue(report['verified']) + self.assertTrue(path.exists()) + with np.load(path.with_suffix('.npz')) as data: + np.testing.assert_array_equal(data['rows'], rows[:7]) + + def test_readonly_import_reuses_actual_offsets_and_functions(self): + with tempfile.TemporaryDirectory() as directory: + root = Path(directory) + (root/'legacy_test/core').mkdir(parents=True) + # Package __init__ must never be executed by the isolated source loader. + for folder in (root, root/'legacy_test', root/'legacy_test/core'): + (folder/'__init__.py').write_text("raise RuntimeError('application import forbidden')\n") + (root/'legacy_test/core/dtw_metrics.py').write_text('calc_dtw_sim = None\ncalc_lag = None\n') + (root/'acquisition.py').write_text( + 'from .legacy_test.core.dtw_metrics import calc_dtw_sim, calc_lag\n' + 'VORTEX_EVENT_OFFSETS = (-20,-15,-10,-5,0,5,10,15)\n') + (root/'legacy_test/acquire.py').write_text( + "raise RuntimeError('heavy imports forbidden')\n" + 'def _vortex_event_channel(case_name):\n return 3, "center_sensor_uy"\n' + 'def _select_vortex_event(sensors, case_name):\n' + ' return dict(offsets=VORTEX_EVENT_OFFSETS, channel=_vortex_event_channel(case_name))\n') + helper = driver.load_helpers(root/'acquisition.py') + result = helper.select_vortex_event(None, 'vortex_taylor') + self.assertEqual(result['offsets'], (-20,-15,-10,-5,0,5,10,15)) + self.assertEqual(result['channel'], (3, 'center_sensor_uy')) + self.assertFalse(list(root.rglob('*.pyc'))) + + def test_geometry_and_shared_translation_envelope(self): + calls, translations = [], [] + def translation(*args, **kwargs): + item = types.SimpleNamespace(args=args, kwargs=kwargs) + translations.append(item) + return item + def add_body(*args, **kwargs): + calls.append((args, kwargs)) + return len(calls)-1 + motion = types.ModuleType('CelerisLab.body.motion') + motion.Translation = translation + with patch.dict(sys.modules, {'CelerisLab.body.motion': motion}): + probes, bodies = driver.add_scene(types.SimpleNamespace(add_body=add_body), 'karman', moving=True) + self.assertEqual(calls[0][0], ('circle', (400.,510.5),40.)) + self.assertEqual(probes, [1,2,3]) + self.assertEqual(bodies, [4,5,6]) + self.assertEqual(translations[0].args, ((0.,0.,-40.,40.), .025)) + self.assertTrue(translations[0].kwargs['shared']) + self.assertTrue(all(calls[i][1]['translation'] is translations[0] for i in bodies)) + calls.clear() + with patch.dict(sys.modules, {'CelerisLab.body.motion': motion}): + driver.add_scene(types.SimpleNamespace(add_body=add_body), 'steady', pinball=False) + self.assertEqual(len(calls), 3) + self.assertTrue(all(call[0][0] == 'sensor' for call in calls)) + + def test_taylor_operator_called_once_inside_context_at_response_zero(self): + events = [] + class Context: + def __enter__(self): + events.append('enter') + def __exit__(self, *args): + events.append('exit') + initializer = types.ModuleType('CelerisLab.lbm.initializers') + initializer.add_vortex = lambda *args: events.append(args) + field = object() + sim = types.SimpleNamespace(field=field, stepper=types.SimpleNamespace(step_count=1024000), + ctx=types.SimpleNamespace(activate=Context), + stream=types.SimpleNamespace(synchronize=lambda: events.append('sync'))) + with patch.dict(sys.modules, {'CelerisLab.lbm.initializers': initializer}): + receipt = driver.inject_taylor(sim) + self.assertEqual(events, ['enter', 'sync', (field, (600.,510.5),80.,.000075,'taylor'), 'sync', 'exit']) + self.assertEqual((receipt['response_clock'], receipt['absolute_clock'], receipt['count']), + (0,1024000,1)) + + +class SelectionContracts(unittest.TestCase): + @staticmethod + def mean(ux, uy, indices): + return dict(mean_ux=ux[indices].mean(axis=0), mean_uy=uy[indices].mean(axis=0)) + + def test_karman_missing_cycles_preserves_fixed_mean_and_last(self): + helper = types.SimpleNamespace(accumulate_mean_fields=self.mean) + def no_crossing(*args, **kwargs): + raise ValueError('at least two accepted rising crossings are required') + helper.assign_periodic_phase = no_crossing + cache = np.arange(24, dtype=np.float32).reshape(3, 2, 2, 2) + arrays, metadata, event = driver.select_fields('karman', dict(smoke=False), helper, + np.array([3200, 6400, 9600]), np.zeros((3, 6)), cache) + self.assertEqual(metadata['periodic_phase_status'], 'insufficient_crossings') + self.assertNotIn('cycle_mean_ux', arrays) + np.testing.assert_array_equal(arrays['ux'][0], cache[-1, 0]) + np.testing.assert_allclose(arrays['fixed_window_mean_ux'], cache[:, 0].mean(axis=0)) + + def test_taylor_common_controlled_event_not_reselected(self): + def forbidden(*args): + raise AssertionError('other roles must not select their own event') + helper = types.SimpleNamespace(accumulate_mean_fields=self.mean, select_vortex_event=forbidden) + cache = np.zeros((150, 2, 2, 2), np.float32) + event = dict(field_indices=np.array([20,25,30,35,40,45,50,55])) + arrays, meta, returned = driver.select_fields('taylor', dict(smoke=False), helper, + np.arange(1,151)*3200, np.zeros((150,6)), cache, event) + self.assertIs(returned, event) + np.testing.assert_array_equal(arrays['selected_index'], event['field_indices']) + self.assertEqual(len(arrays['ux']), 8) + + def test_smoke_selection_has_no_event_or_phase_claim(self): + helper = types.SimpleNamespace(accumulate_mean_fields=self.mean) + for scene in ('karman', 'steady', 'taylor'): + arrays, metadata, event = driver.select_fields(scene, dict(smoke=True), helper, + np.array([32,64,96]), np.zeros((3,6)), + np.zeros((3,2,2,2), np.float32)) + self.assertEqual(metadata['selection'], 'smoke_last_snapshot_not_physical') + self.assertIsNone(event) + + +if __name__ == '__main__': + unittest.main() diff --git a/tests/unit/test_fig2_metrics.py b/tests/unit/test_fig2_metrics.py new file mode 100644 index 0000000..f5555af --- /dev/null +++ b/tests/unit/test_fig2_metrics.py @@ -0,0 +1,49 @@ +"""CPU-only analytic checks, runnable with stdlib unittest (or pytest).""" +import unittest +import numpy as np +from research.fig2_materials.metrics import ( + common_fluid_roi_mask, field_error, motion_summary, roi_slices, +) + + +class MetricsTests(unittest.TestCase): + def test_u0_normalization_not_target_energy(self): + zero = np.zeros((2, 2)); support = np.ones((2, 2), bool) + self.assertEqual(field_error(zero + 3, zero + 4, zero, zero, 2, support), 2.5) + self.assertEqual(field_error(zero + 13, zero + 4, zero + 10, zero, 2, support), 2.5) + for u0 in (0, -1, np.nan): + with self.assertRaises(ValueError): + field_error(zero, zero, zero, zero, u0, support) + + def test_common_mask_roi_finite_and_fluid(self): + a = np.zeros((3, 4)); b = a.copy(); b[1, 2] = np.nan + m = np.ones(a.shape, bool); n = m.copy(); n[1, 1] = False + support = common_fluid_roi_mask([m, n], [a, b], (slice(1, 3), slice(1, 4))) + self.assertEqual(int(support.sum()), 4) + self.assertFalse(support[1, 1]); self.assertFalse(support[1, 2]) + self.assertTrue(support[2, 1]) # zero velocity is fluid + self.assertEqual(field_error(a, a, b, a, 1, support), 0) + with self.assertRaises(ValueError): + common_fluid_roi_mask([m & False], [a], (slice(None), slice(None))) + with self.assertRaises(ValueError): + field_error(a, a, b, a, 1, m) # fixed support must be finite + + def test_mean_field_differs_from_mean_instantaneous_error(self): + snapshots = np.array([[[2.]], [[-2.]]]); z = np.zeros((1, 1)); mask = np.ones((1, 1), bool) + self.assertEqual(field_error(snapshots.mean(0), z, z, z, 1, mask), 0) + self.assertEqual(np.mean([field_error(s, z, z, z, 1, mask) for s in snapshots]), 2) + + def test_actual_geometry_and_motion_bias(self): + roi = roi_slices((1022, 2400), D=40, xs=1600, y0=510.5, downstream_solid_extent=1350) + self.assertEqual(roi, (slice(311, 711), slice(1360, 2161))) + m = motion_summary([8, 10, 12]) + self.assertEqual(m['mean'], 10) + self.assertAlmostEqual(m['centered_rms'], np.sqrt(8 / 3)) + self.assertEqual(m['half_range'], 2) + for a in ([], [np.nan]): + with self.assertRaises(ValueError): + motion_summary(a) + + +if __name__ == '__main__': + unittest.main()