558 lines
30 KiB
Python
558 lines
30 KiB
Python
"""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 check_clock(stepper, release_clock, completed):
|
|
if stepper.step_count != release_clock+completed:
|
|
raise RuntimeError('observed solver clock differs from recorded response clock')
|
|
|
|
|
|
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)
|
|
first = sim.bodies.count
|
|
for point in CENTRES:
|
|
sim.add_body('circle', point, 20.) # runtime additions return staged ID -1
|
|
sim.sync_bodies()
|
|
bodies = list(range(first, first+len(CENTRES)))
|
|
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
|
|
check_clock(sim.stepper, plan['release_clock'], completed)
|
|
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()
|