From 76873d911c4e2caba1bacc25aaa1495ea768453f Mon Sep 17 00:00:00 2001 From: Frank14f <1515444314@qq.com> Date: Mon, 2 Nov 2026 21:23:38 +0800 Subject: [PATCH] Add continuous OFF ON OFF VIV video research driver --- docs/continuous_viv_video.md | 86 +++++++++ research/continuous_viv/render.py | 220 +++++++++++++++++++++++ research/continuous_viv/run.py | 222 ++++++++++++++++++++++++ tests/unit/test_continuous_viv_video.py | 148 ++++++++++++++++ 4 files changed, 676 insertions(+) create mode 100644 docs/continuous_viv_video.md create mode 100644 research/continuous_viv/render.py create mode 100644 research/continuous_viv/run.py create mode 100644 tests/unit/test_continuous_viv_video.py diff --git a/docs/continuous_viv_video.md b/docs/continuous_viv_video.md new file mode 100644 index 0000000..e92f50b --- /dev/null +++ b/docs/continuous_viv_video.md @@ -0,0 +1,86 @@ +# Continuous VIV videos + +Research-only, two separate settings, each one uninterrupted OFF -> ON -> OFF +fluid/oscillator trajectory. No production changes, topology switching, moving +checkpoint feature, retraining, force filtering, or parameter search. + +The user selected separate videos because runtime topology changes are unsupported +with translating circles. This is deliberately different from the old Legacy +multi-setting video, which rebuilt solvers and transferred full state. + +## Frozen physical and control contract + +- Source baseline `73d3f6214e0d28ced4744665aa3bb05522aa2fa8`; D40, U=.005, + nu=.004, MRT/FP32/unshifted double-buffer. Same parabolic regularized inlet, + NEQ outlet and bounce-back walls as the retained m*=4 pair. +- Common transverse assembly: m*=4 (all three disks), fn D/U=.2, zeta=0. + Structure uses previous native-step sum CME_y+Wen_y. No policy-force averaging + is applied to this RHS. Same +/-D travel, speed, finite-state and no-donor stops. +- Each phase: Karman400 complete SI of3200, steady320 SI of4000 native steps + =160 D/U. Steady's SI follows the actual old video1000-step interval under + diffusive D40 scaling, rather than borrowing PPO's800-step interval. Total480 D/U. + No extension, retuning or replacement if steady motion/control is weak. +- Steady: canonical new-reproduce constant surface velocities + (0,-5.1,+5.1) U. This is NOT a PPO transfer. No model/normalizer is loaded. + Existing 2,048,000-step at-rest warmup precedes the free response. +- Karman: actual pinball Re_D50, historically named `karman_re100`, frozen + model `148240c1dcb8b11d8e5a0a9e991f6874c4e799c9744707bf0fb32a8efeb0468d` + and norm `534fda818c030a56960e95d0f3ac1e350d887d03db7f80f1d6453195a502b7f0`. + Use the verified static release once, at video start (clock4576000). +- Both switches occur on SI boundaries. The oscillator y/v, previous RHS, + both DDF buffers, moving geometry and solver clock are never reset. OFF still + gathers sensor/CME observations. First PPO action uses the immediately previous + complete OFF interval, unlike an independent arm's initial all-zero observation. + Deterministic feed-forward PPO has no recurrent state. EMA resets to zero each + SI, as in training/reproduction; original video's SplitAdvance also reset the + constant-steady actuation every interval. Physical OFF is exactly zero wall speed. + There is no extra transition blend. The first1 D/U after either switch is + explicitly labeled a switch transient. ON from an already moving state differs + from the previous independent pair; +/-D stops are retained without retuning. + A visual ON label means enabled control, + not a proven suppression/cloaking result. + +## Video provenance and presentation + +Read-only source: DynamisLab `src/drl_pinball/reproduce/continuous_video.py`, SHA256 +`3bd9c777171c9e8932c78298accbb6f61a8323b0eed5938e0f56e167b3df0743`. +Approved original MP4 SHA256 +`578d2cf5113429006b6919bbb8ef8ff6ce65ebcd0c056ac478f8ce2515e4cbab`. +The actual new-reproduce conversation confirms fixed-position plain-language +labels, white bars, OFF/ON toggle, rotation sectors, and 1/2/5/10/5/2/1 gears. + +Only the display idiom/gear schedule are reused. Legacy population decoding, +physics, topology handoff and policy runner are not imported. Modern macroscopic +fields feed centered vorticity with all stencil nodes fluid. Fixed color limits +[-3,+3], 21 RdBu_r bins; per-frame clipping fraction retained. Box downsampling +is display-only. Both main field and enlarged inset have fixed lab coordinates; +the inset shows the equilibrium line. Actual moving centers and applied surface +speeds drive the body overlays. y/D trace has fixed +/-1 limits. No motion gain. + +1280x800,30fps; first80 D/U of each phase uses the original display gear ramps, +last80 D/U uses normal cadence (800 native steps/frame =.1D/U). Every native +step and every SI is still executed. Streaming ffmpeg avoids stored PNG sequences. +Keep raw scalar histories, SI observations/actions, per-frame time/motion records, +selected keyframes, source/input hashes and terminal status. These are new +continuous runs, not an exact replay of the old independent paired histories. + +## Gates and limits + +CPU: switch-state/clock continuity, first ON observation, exactly zero OFF action, +SI timing despite sampling, physical scaling, schedule monotonicity, vorticity +sign/mask, fixed-coordinate mapping. Real smoke:32 warmup +3x32 response steps +on D40 (smoke SI32), all switches/encoder; not a physical comparison. Grok reviews +the exact research diff before full runs. Capture no-donor, finite state, native +RHS and all stop outcomes. Report per-phase late-half motion statistics and +raw-force limits; appearance and successful duration are not physics qualification. + +Example (outputs must be new, outside checkout; select a safe GPU): + +```sh +CUDA_VISIBLE_DEVICES=0 PYTHONPATH=src python research/continuous_viv/run.py --scene steady --out /tmp/celeris-results-continuous-viv/steady +CUDA_VISIBLE_DEVICES=0 PYTHONPATH=src python research/continuous_viv/run.py --scene karman --model /tmp/celeris-karman-viv-assets/model.zip --norm /tmp/celeris-karman-viv-assets/norm.json --release-dir /tmp/celeris-results-karman-viv-516975e/pair --out /tmp/celeris-results-continuous-viv/karman +``` + +Append `--smoke` and use separate output directories for the bounded gates. +The Karman smoke does not read/install the physical release. Source and review +SHA, interpreter/version and GPU resource readings belong in the execution receipt. diff --git a/research/continuous_viv/render.py b/research/continuous_viv/render.py new file mode 100644 index 0000000..6bc16a6 --- /dev/null +++ b/research/continuous_viv/render.py @@ -0,0 +1,220 @@ +"""Sparse scientific frames and the approved new-reproduce presentation idiom. + +Frame gears/UI colors/rotation sectors derive from DynamisLab continuous_video.py +(read-only SHA in docs). Modern velocity decoding and moving coordinates are used. +All camera limits, vorticity limits, and motion scale are fixed, never following y. +""" +from functools import lru_cache +import json +from pathlib import Path +import subprocess +import time + +import numpy as np + +D,U,CY=40.,.005,510.5 +SIZE=(1280,800) +CW,CCW='#8854B8','#138F78' +MAIN=(0,120,1280,440) +INSET=(960,492,1240,772) +MAIN_LIMITS=(6.4,57.6,-6.4,6.4) +INSET_LIMITS=(28.75,33.75,-2.5,2.5) + + +def frame_schedule(steps,base_interval,accelerated=False): + """Original 15-frame gear holds; only readback cadence changes.""" + if accelerated: + ramp=15*(1+2+5)*base_interval + if steps<=2*ramp: + raise ValueError('too short for both display ramps') + blocks=[(s,15*s*base_interval) for s in (1,2,5)] + blocks += [(10,steps-2*ramp)] + blocks += [(s,15*s*base_interval) for s in (5,2,1)] + else: + blocks=[(1,steps)] + result,done=[],0 + for speed,duration in blocks: + end=done+duration + while done1.00001: + raise RuntimeError('applied rotation exceeds fixed display reference') + r2=radius*.72 + inner=(cx-r2,cy-r2,cx+r2,cy+r2) + draw.ellipse(inner,fill='#E4E9ED') + angle=min(fraction,1)*360 + start,end=(-90,-90+angle) if wall_speed<0 else (-90-angle,-90) + if angle>0: + draw.pieslice(inner,start,end,fill=CW if wall_speed<0 else CCW) + + +class Video: + def __init__(self,out,scene,bodies,phase_steps,smoke=False,source='unknown'): + from PIL import Image + from matplotlib import colormaps + self.out,self.scene,self.bodies=out,scene,bodies + self.schedule=video_schedule(phase_steps,smoke) + self.phase_steps=phase_steps + self.source=source + self.records=[] + self.palette=np.rint(colormaps['RdBu_r'](np.linspace(0,1,21))[:,:3]*255).astype('uint8') + self.log=open(out/'ffmpeg.log','wb') + self.proc=subprocess.Popen(['ffmpeg','-hide_banner','-loglevel','warning','-n', + '-f','rawvideo','-pix_fmt','rgb24','-s','1280x800','-r','30','-i','pipe:0', + '-an','-c:v','libx264','-threads','4','-preset','fast','-crf','18', + '-pix_fmt','yuv420p','-movflags','+faststart',str(out/'continuous_viv.mp4')], + stdin=subprocess.PIPE,stderr=self.log) + self.closed=False + self.seconds=0. + + def capture(self,sim,step,phase,eff,rows): + speed=self.schedule.get(step) + if speed is None: + return + from PIL import Image,ImageDraw + start=time.monotonic() + macro=sim.get_macroscopic() + flags=sim.get_flags().reshape(sim.field.ny,sim.field.nx) + fluid=(flags&1)!=0 + if any(not np.isfinite(macro[k][fluid]).all() for k in ('rho','ux','uy')): + raise RuntimeError('nonfinite fluid field at video sample') + omega=vorticity(macro['ux'],macro['uy'],flags) + finite=np.isfinite(omega) + index=np.floor((np.clip(np.nan_to_num(omega),-3,3)+3)/6*21).astype(int).clip(0,20) + rgb=self.palette[index] + rgb[~finite]=[232,235,237] + field=Image.fromarray(rgb[::-1]) + canvas=Image.new('RGB',SIZE,'white') + draw=ImageDraw.Draw(canvas) + for box,limits in ((MAIN,MAIN_LIMITS),(INSET,INSET_LIMITS)): + x0,x1,y0,y1=limits + # Lattice cell-centre coordinates -> image pixel-edge coordinates. + crop=(x0*D+.5,sim.field.ny-(CY+y1*D+.5), + x1*D+.5,sim.field.ny-(CY+y0*D+.5)) + picture=field.crop(crop).resize((box[2]-box[0],box[3]-box[1]),Image.Resampling.BOX) + canvas.paste(picture,box[:2]) + if box==INSET: + _,eq=pixel(30,0,box,limits) + for x in range(box[0],box[2],10): + draw.line((x,eq,min(x+5,box[2]),eq),fill='#657887',width=1) + if self.scene=='karman': + draw_circle(draw,10,0,1,None,box,limits) if box==MAIN else None + for bid,w in zip(self.bodies,eff): + obj=sim.bodies.objects[bid] + x,y=obj.motion.centre(obj.motion.applied) + draw_circle(draw,x/D,(y-CY)/D,obj.radius/D,float(w),box,limits) + x0,y0=pixel(INSET_LIMITS[0],INSET_LIMITS[2],MAIN,MAIN_LIMITS) + x1,y1=pixel(INSET_LIMITS[1],INSET_LIMITS[3],MAIN,MAIN_LIMITS) + draw.rectangle((x0,y1,x1,y0),outline='#647C8C',width=1) + draw.rectangle(INSET,outline='#647C8C',width=1) + def text(x,y,value,size=16,color='#263D4B',bold=False): + draw.text((x,y),value,font=font(size,bold),fill=color) + on=phase=='on' + text(28,15,'INCOMING FLOW',12,'#667887',True) + text(28,39,'Steady inflow' if self.scene=='steady' else 'Periodic cylinder wake',29,bold=True) + text(28,83,'One continuous flow and oscillator trajectory | Re_D = 50 | zero structural damping',15) + text(850,18,'CYLINDER ROTATION',12,bold=True) + text(850,43,'Clockwise',14,CW);text(1020,43,'Counterclockwise',14,CCW) + text(850,70,'Sector angle = actual surface speed',13,'#667887') + text(850,90,'Full circle = 12 x reference inflow speed',13,'#667887') + draw.line((24,453,1256,453),fill='#E5EBEF') + text(28,465,'CLOAKING CONTROL',12,'#667887',True) + draw.rounded_rectangle((28,490,92,524),radius=17,fill='#137F73' if on else '#A8B4BC') + cx=75 if on else 45 + draw.ellipse((cx-12,495,cx+12,519),fill='white') + text(108,491,'ON' if on else 'OFF',22,'#137F73' if on else '#667887',True) + text(195,493,('Fixed rotations' if self.scene=='steady' else 'Frozen PPO') if on else 'Zero rotation',18) + text(570,467,f't U/D = {step*U/D:.2f}',18,bold=True) + local=(step-1)%self.phase_steps+1 + transition=(phase!='off_before' and local*U/D<=1.) + text(570,495,f'Playback {speed}x | '+('switch transient' if transition else phase.replace('_',' ')),16) + text(960,467,f'Fixed lab-frame inset | y/D={rows[-1,1]/D:+.3f}',13,bold=True) + # Motion trace: fixed +/-D limits, a fixed 60 D/U trailing time span. + box=(55,560,895,732) + draw.rectangle(box,outline='#CCD6DD') + middle=(box[1]+box[3])/2 + draw.line((box[0],middle,box[2],middle),fill='#A8B4BC') + t=step*U/D + left=max(0.,t-60) + first=max(0,int(left*D/U)) + sample=rows[first::max(1,(len(rows)-first)//1200)] + xy=[(box[0]+(r[0]*U/D-left)/60*(box[2]-box[0]),middle-r[1]/D*(box[3]-box[1])/2) for r in sample] + if len(xy)>1:draw.line(xy,fill='#263D4B',width=2) + text(28,535,'Transverse displacement y/D (fixed scale; dashed inset line = equilibrium)',14) + text(24,554,'+1',12);text(25,642,'0',12);text(24,722,'-1',12) + text(55,739,f'{left:.1f}',12);text(835,739,f'{left+60:.1f}',12) + text(350,739,'response time t U/D',12) + text(28,774,f'Vorticity D/U: fixed [-3, +3] | exploratory, not convergence proof | source {self.source[:8]}',12) + self.proc.stdin.write(canvas.tobytes()) + local=step%self.phase_steps + if not self.records or local in (0,self.phase_steps//2) or self.records[-1]['phase']!=phase: + canvas.save(self.out/f'keyframe-{step:09d}.png') + elapsed=time.monotonic()-start + self.seconds+=elapsed + self.records.append(dict(frame=len(self.records),step=step,absolute_clock=sim.stepper.step_count, + tU_D=t,phase=phase,playback_speed=speed,y_D=float(rows[-1,1]/D), + surface_speeds=eff.tolist(),render_s=elapsed, + fraction_finite_outside_color_limits=float(np.mean(np.abs(omega[finite])>3)))) + + def close(self): + if self.closed:return + self.closed=True + self.proc.stdin.close() + code=self.proc.wait(timeout=60) + self.log.close() + receipt=dict(frames=self.records,frame_count=len(self.records),fps=30,size=SIZE, + duration_s=len(self.records)/30,render_total_s=self.seconds,ffmpeg_exit=code, + normalization='centered dx=1 vorticity D/U, all 5 stencil cells fluid', + color_limits=[-3,3],color_bins=21,main_limits_D=MAIN_LIMITS,inset_limits_D=INSET_LIMITS, + y_origin_lattice=CY,trace_limits_D=[-1,1],inset_camera='fixed lab frame', + acceleration='frame cadence only; native solver and SI unchanged') + (self.out/'video.json').write_text(json.dumps(receipt,indent=2)+'\n') + if code:raise RuntimeError(f'ffmpeg exit {code}') diff --git a/research/continuous_viv/run.py b/research/continuous_viv/run.py new file mode 100644 index 0000000..84b4f81 --- /dev/null +++ b/research/continuous_viv/run.py @@ -0,0 +1,222 @@ +"""Continuous OFF -> ON -> OFF VIV illustration; docs/continuous_viv_video.md. + +Each invocation owns ONE Simulation and ONE oscillator. Rendering never advances +the solver or changes the policy clock. No topology changes or moving restart. +""" +from __future__ import annotations + +import argparse +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.karman_viv import run as base +from research.continuous_viv.render import Video, frame_schedule + +D, U, SCALE, SI = 40., .005, 2., 3200 +PHASES = ('off_before', 'on', 'off_after') +STEADY_WALL = np.array([0., -5.1, 5.1]) * U + + +def timing(scene, smoke): + if smoke: + return 1,32,32 + # Prior video's steady interval was 1000 at D20; diffusive D40 uses 4000. + return (320,4000,2048000) if scene=='steady' else (400,SI,2048000) + + +def command(scene, phase, obs, model): + """Physical surface speeds; OFF is not a zero normalized PPO action.""" + if phase != 'on': + return np.zeros(3), np.array([0., .5, -.5]) + if scene == 'steady': + return STEADY_WALL.copy(), np.full(3, np.nan) + action = 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 PPO action') + return base.wall_target(action, SCALE), action + + +def add_scene(sim, scene): + if scene == 'karman': + return base.add_scene(sim, SCALE, moving=True) + from CelerisLab.body.motion import Translation + probes = [sim.add_body('sensor', c, 10.) for c in base.scaled_points(base.PROBES, SCALE)] + envelope = Translation((0., 0., -D, D), .025, shared=True) + bodies = [sim.add_body('circle', c, D/2, translation=envelope) + for c in base.scaled_points(base.CENTRES, SCALE)] + return probes, bodies + + +def evolve(sim, probes, bodies, structure, force, *, scene, model, norm, + intervals, si, out, render, start, wall_cap): + """One native-step loop across all phases; completed SI feeds the next SI.""" + from CelerisLab.body.motion import MovingTables + cols = {k:i for i,k in enumerate(MovingTables.DIAG)} + count = len(PHASES)*intervals*si + rows = np.lib.format.open_memmap(out/'steps.npy', mode='w+', dtype=np.float64, + shape=(count, len(base.STEP_COLUMNS))) + # The completed count is authoritative after any interrupted run. + completed = 0 + obs = np.zeros(12, np.float32) + signals, switches = [], [] + origin_clock = sim.stepper.step_count + status = 'running' + try: + for phase in PHASES: + switches.append(dict(phase=phase, step=completed, absolute_clock=sim.stepper.step_count, + y=structure.y, v=structure.v, previous_rhs=force, + input_observation=obs.tolist())) + for i in range(intervals): + target, action = command(scene, phase, obs, model) + # Match the frozen inference contract: reset the EMA every SI. + eff = np.zeros(3) + sensor_sum, force_sum = np.zeros(6), np.zeros(6) + for k in range(si): + if k % 1000 == 0: + base.check_time(start, wall_cap) + y1, v1 = structure.candidate(force) + if not np.isfinite([y1,v1]).all() or abs(y1)>D: + raise RuntimeError('nonfinite structure or |y| > D; no retuning') + eff = base.ema(eff, target, SCALE) + omega = eff/(D/2) + for bid,w in zip(bodies,omega): + sim.set_body(bid, 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') + sensor_sum += s + force_sum += f + force=float(cme.sum(dtype=np.float64)+wen.sum(dtype=np.float64)) + rows[completed]=row + completed += 1 + structure.y,structure.v=y1,v1 + if sim.stepper.step_count != origin_clock+completed: + raise RuntimeError('noncontinuous solver clock') + render(sim, completed, phase, eff, rows[:completed]) + if scene == 'karman': + raw, obs, clipped = base.observation(sensor_sum/si,force_sum/si,SCALE,norm) + else: + # Steady uses prescribed rotations, not a model or borrowed norm. + raw=np.r_[sensor_sum/si,force_sum/si] + clipped=0 + signals.append(dict(phase=phase, step=completed, raw=raw.tolist(), + observation=obs.tolist(), clipped=clipped, + action=action.tolist() if np.isfinite(action).all() else None, + target_surface_speed=target.tolist())) + if (i+1)%10==0 or i+1==intervals: + rows.flush() + progress=dict(phase=phase, interval=i+1, completed_steps=completed, + tU_D=completed*U/D, y_D=structure.y/D, + elapsed_s=time.monotonic()-start) + (out/'progress.json').write_text(json.dumps(progress)+'\n') + print(json.dumps(progress),flush=True) + status='duration_completed_not_physics_qualified' + except BaseException as exc: + status=f'error: {type(exc).__name__}: {exc}' + raise + finally: + rows.flush() + (out/'signals.json').write_text(json.dumps(signals,allow_nan=False)+'\n') + report=dict(status=status,completed_steps=completed,origin_clock=origin_clock, + final_clock=sim.stepper.step_count,switches=switches, + elapsed_s=time.monotonic()-start,columns=base.STEP_COLUMNS, + M=structure.mass,K=structure.stiffness,fn=structure.frequency,damping=0., + max_abs_y_D=float(np.abs(rows[:completed,1]).max()/D) if completed else None) + (out/'report.json').write_text(json.dumps(report,indent=2,allow_nan=False)+'\n') + + +def main(): + ap=argparse.ArgumentParser(description=__doc__) + ap.add_argument('--scene',choices=('steady','karman'),required=True) + ap.add_argument('--out',type=Path,required=True) + ap.add_argument('--model',type=Path) + ap.add_argument('--norm',type=Path) + ap.add_argument('--release-dir',type=Path) + ap.add_argument('--smoke',action='store_true',help='3 x 32 steps, all switches, no physical result') + ap.add_argument('--wall-seconds',type=float,default=10800.) + args=ap.parse_args() + if not np.isfinite(args.wall_seconds) or args.wall_seconds<=0: + ap.error('finite positive wall cap required') + if args.scene=='karman': + if args.model is None or args.norm is None: + ap.error('Karman requires the frozen model and norm') + if base.sha(args.model)!=base.MODEL_SHA or base.sha(args.norm)!=base.NORM_SHA: + ap.error('frozen Karman model/normalizer hash mismatch') + if not args.smoke and args.release_dir is None: + ap.error('physical Karman requires its verified static release') + elif any(p is not None for p in (args.model,args.norm,args.release_dir)): + ap.error('steady uses constant control and its own warmup, no borrowed policy or release') + 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) + intervals,si,warm=timing(args.scene,args.smoke) + cfgpath=out/'config.json' + cfgpath.write_text(json.dumps(base.config(40),indent=2)+'\n') + source=subprocess.check_output(['git','-C',str(ROOT),'rev-parse','HEAD'],text=True).strip() + case=dict(source=source,scene=args.scene,smoke=args.smoke,D=D,U=U,nu=.004, + phases=PHASES,SI=si,intervals_per_phase=intervals, + phase_duration_tU_D=intervals*si*U/D, mass_ratio=4.,fn_D_U=.2,zeta=0., + first_ppo_observation='previous complete OFF interval; no reset at switches', + rhs='previous native step CME_y+Wen_y; no filtering; policy uses SI-mean CME', + control='fixed (0,-5.1,+5.1)U' if args.scene=='steady' else 'frozen karman_re100 PPO; actual Re_D50', + model_sha256=base.sha(args.model) if args.model else None, + norm_sha256=base.sha(args.norm) if args.norm else None, + config_sha256=base.sha(cfgpath),render='sparse post-step snapshots only; no step skipping') + model=base.load_policy(args.model) if args.model else None + norm=json.loads(args.norm.read_text()) if args.norm else None + seed=None + if args.scene=='karman' and not args.smoke: + seed,force=base.saved_release(args.release_dir) + case['release_sha256']=base.sha(args.release_dir/'release.npz') + (out/'case.json').write_text(json.dumps(case,indent=2)+'\n') + from CelerisLab import Simulation + from CelerisLab.body.motion import MovingTables + started=time.monotonic() + structure=base.Structure(4.,2.5e-5,0.,D) + video=None + try: + with Simulation(str(cfgpath)) as sim: + probes,bodies=add_scene(sim,args.scene) + sim.initialize() + if seed is not None: + case['release_transfer']=base.install_seed(sim,seed) + del seed + else: + base.warmup(sim,warm-1,started,args.wall_seconds) + sim.run(1) + d=sim.moving_diagnostics() + force=float(d[:,MovingTables.DIAG.index('cme_y')].sum(dtype=np.float64) + +d[:,MovingTables.DIAG.index('wen_y')].sum(dtype=np.float64)) + case.update(warmup_executed=0 if 'release_transfer' in case else warm, + origin_clock=sim.stepper.step_count) + (out/'case.json').write_text(json.dumps(case,indent=2)+'\n') + phase_steps=intervals*si + video=Video(out,args.scene,bodies,phase_steps,smoke=args.smoke,source=source) + evolve(sim,probes,bodies,structure,force,scene=args.scene,model=model,norm=norm, + intervals=intervals,si=si,out=out,render=video.capture, + start=started,wall_cap=args.wall_seconds) + finally: + if video is not None: + video.close() + + +if __name__=='__main__': + main() diff --git a/tests/unit/test_continuous_viv_video.py b/tests/unit/test_continuous_viv_video.py new file mode 100644 index 0000000..f6c7733 --- /dev/null +++ b/tests/unit/test_continuous_viv_video.py @@ -0,0 +1,148 @@ +"""Research video contracts: continuous state and display-only sampling.""" +import importlib.util +import json +from pathlib import Path +from types import SimpleNamespace + +import numpy as np +import pytest + +ROOT=Path(__file__).resolve().parents[2] +spec=importlib.util.spec_from_file_location('continuous_viv_video_driver',ROOT/'research/continuous_viv/run.py') +m=importlib.util.module_from_spec(spec) +spec.loader.exec_module(m) +from research.continuous_viv import render as r + + +def test_physical_zero_and_independent_steady_control(): + model=SimpleNamespace(predict=lambda *a,**k:pytest.fail('steady/OFF must not invoke PPO')) + for scene in ('steady','karman'): + for phase in ('off_before','off_after'): + target,action=m.command(scene,phase,None,model) + np.testing.assert_array_equal(target,[0,0,0]) + np.testing.assert_array_equal(m.base.wall_target(action,2),target) + target,action=m.command('steady','on',None,model) + np.testing.assert_array_equal(target,np.array([0,-5.1,5.1])*.005) + assert np.isnan(action).all() + + +@pytest.mark.parametrize('bad',[[2,0,0],[np.nan,0,0],[0,0]]) +def test_invalid_policy_output(bad): + model=SimpleNamespace(predict=lambda *a,**k:(bad,None)) + with pytest.raises(RuntimeError,match='invalid PPO'): + m.command('karman','on',np.zeros(12),model) + + +def test_continuous_switches_previous_observation_and_native_rhs(tmp_path): + from CelerisLab.body.motion import MovingTables + calls,forces,seen,frames=[],[],[],[] + class Structure: + y=v=0.;mass=2.;stiffness=3.;frequency=.1 + def candidate(self,f): + forces.append(f) + return self.y+.001,self.v+.002 + class Sim: + stepper=SimpleNamespace(step_count=100) + n=0 + def set_body(self,p,**kw):calls.append((p,kw)) + def run(self,n,**kw): + assert n==1 and kw==dict(zero_obs=True,sync_obs=True) + self.n+=1;self.stepper.step_count+=1 + def moving_diagnostics(self): + d=np.zeros((3,len(MovingTables.DIAG)),np.float32) + d[:,MovingTables.DIAG.index('cme_y')]=self.n + d[:,MovingTables.DIAG.index('wen_y')]=10*self.n + d[:,MovingTables.DIAG.index('cover_jy')]=1000 + return d + def read_force(self,p):return [0.,float(self.n),0.] + def read_sensor(self,p):return [.001,0.,0.] + model=SimpleNamespace(predict=lambda obs,**kw:(seen.append(obs.copy()) or np.ones(3)*.1,None)) + norm=dict(force_norm_fact=100.,sens_deviation=np.zeros(6),sens_norm_fact=np.ones(6)) + sim,structure=Sim(),Structure() + def render(sim,step,phase,eff,rows): + if step%3==0:frames.append((step,sim.stepper.step_count,eff.copy())) + m.evolve(sim,[1,2,3],[4,5,6],structure,7.,scene='karman',model=model,norm=norm, + intervals=2,si=2,out=tmp_path,render=render,start=m.time.monotonic(),wall_cap=60) + assert sim.n==12 and sim.stepper.step_count==112 + assert forces==[7.,*list(33.*np.arange(1,12))] # no reset, SI averaging, or node impulse in RHS + rows=np.load(tmp_path/'steps.npy') + np.testing.assert_allclose(rows[:,1],.001*np.arange(1,13)) + np.testing.assert_allclose(rows[:,2],.002*np.arange(1,13)) + np.testing.assert_array_equal(rows[:4,15:],0) + np.testing.assert_array_equal(rows[8:,15:],0) + np.testing.assert_array_equal(rows[4,15:],rows[6,15:]) # per-SI EMA reset, not per-frame + signals=json.loads((tmp_path/'signals.json').read_text()) + np.testing.assert_array_equal(seen[0],np.asarray(signals[1]['observation'],np.float32)) + np.testing.assert_array_equal(seen[1],np.asarray(signals[2]['observation'],np.float32)) + assert len(seen)==2 and seen[0].any() + report=json.loads((tmp_path/'report.json').read_text()) + assert [x['absolute_clock'] for x in report['switches']]==[100,104,108] + assert [x['previous_rhs'] for x in report['switches']]==[7.,132.,264.] + np.testing.assert_allclose([x['y'] for x in report['switches']],[0,.004,.008]) + assert [x[0] for x in frames]==[3,6,9,12] + assert all(x[1]==100+x[0] for x in frames) + assert len(calls)==36 and all(c[1]['vy']==pytest.approx(.001) for c in calls) + + +def test_gear_schedule_and_complete_phase_endpoints(): + n=1280000 + schedule=r.video_schedule(n) + times=list(schedule) + assert times==sorted(set(times)) and times[-1]==3*n + assert all(p*n in schedule for p in (1,2,3)) + assert set(schedule.values())=={1,2,5,10} + assert all(speed==1 for step,speed in schedule.items() if step%n>n//2) + part=r.frame_schedule(n//2,800,True) + assert [s for _,s in part[:15]]==[1]*15 + assert [s for _,s in part[15:30]]==[2]*15 + assert [s for _,s in part[-15:]]==[1]*15 + assert part[-1][0]==n//2 + + +def test_smoke_schedule_covers_all_switches(): + schedule=r.video_schedule(32,True) + assert list(schedule)==list(range(8,97,8)) + assert set(schedule.values())=={1} + + +def test_vorticity_sign_scale_and_solid_stencil_mask(): + y,x=np.mgrid[:9,:11] + ux=-.0001*y;uy=.0001*x + flags=np.ones(x.shape,np.uint16) + omega=r.vorticity(ux,uy,flags) + np.testing.assert_allclose(omega[1:-1,1:-1],.0002*40/.005) + assert np.isnan(omega[0]).all() + flags[4,5]=2 + omega=r.vorticity(ux,uy,flags) + for yy,xx in ((4,5),(3,5),(5,5),(4,4),(4,6)): + assert np.isnan(omega[yy,xx]) + assert np.isfinite(omega[3,4]) + + +def test_fixed_inset_coordinate_and_no_motion_gain(): + a=np.array(r.pixel(30,0,r.INSET,r.INSET_LIMITS)) + b=np.array(r.pixel(30,1,r.INSET,r.INSET_LIMITS)) + c=np.array(r.pixel(31,0,r.INSET,r.INSET_LIMITS)) + assert b[0]==a[0] and a[1]-b[1]==56 + assert c[0]-a[0]==56 + assert r.pixel(r.INSET_LIMITS[0],r.INSET_LIMITS[3],r.INSET,r.INSET_LIMITS)==r.INSET[:2] + + +def test_config_real_loader_and_structure_unchanged(tmp_path): + from CelerisLab.config import load_lbm_config + p=tmp_path/'config.json';p.write_text(json.dumps(m.base.config(40))) + cfg=load_lbm_config(str(p));cfg.validate() + assert (cfg.nx,cfg.ny,cfg.velocity,cfg.viscosity)==(2559,1022,.005,.004) + s=m.base.Structure(4.,2.5e-5,0.,40.) + assert s.mass==pytest.approx(15079.644737231007) + assert s.stiffness==pytest.approx(.00037207532016359776) + assert s.damping==0 and s.frequency*40/.005==.2 + + +def test_scene_timing_keeps_actual_legacy_intervals_and_same_physical_duration(): + assert m.timing('steady',False)==(320,4000,2048000) + assert m.timing('karman',False)==(400,3200,2048000) + for scene in ('steady','karman'): + n,si,warm=m.timing(scene,False) + assert n*si*.005/40==160 + assert m.timing(scene,True)==(1,32,32)