diff --git a/docs/steady_pinball_fiv.md b/docs/steady_pinball_fiv.md new file mode 100644 index 0000000..714c9c3 --- /dev/null +++ b/docs/steady_pinball_fiv.md @@ -0,0 +1,69 @@ +# Uncontrolled steady pinball: development protocol + +Frozen before results: 1280 x 512 D2Q9/FP32/MRT, unshifted double buffer; +rho=1, mean U=.01, nu=.004, D=20, Re_D=50. Centres are (600,255.5), +(626,270.5), (626,240.5), radius10. No upstream cylinder. Passive probes stay +at (800,295.5), (800,255.5), (800,215.5), radius5. Every omega is zero. +Parabolic inlet, regularized NEQ damping1; y bounce-back planes .5/510.5; +NEQ-extrapolated outlet with existing backflow clamp/blending. Initial fluid is +parabolic equilibrium. Fixed warmup is 4 NX/U=512000 steps. + +Geometry and fluid parameters follow DynamisLab 609116ec26f19766f5e9ce2b0df427d4e24faae8, +`legacy_test/core/legacy_env_builder.py:build_steady_cloak` and the retained +`config_lbm_pinball_legacy_re50d_u001.json`. This is a modern Celeris reconstruction, +not a bitwise Legacy reproduction: Legacy uses shifted incompressible MRT, a +different curved interpolation and prescribed parabolic outlet velocity. Current +Celeris uses standard MRT, two-branch BFL and extrapolated outlet velocity. +The fixed-hold control uses the SAME modern settings as the free case. + +## Shared translation and ownership + +`Translation(..., shared=True)` is explicitly transverse only. Reusing one +envelope INSTANCE declares common motion; equal but separate envelopes are not +linked. Before initialization and every positive Simulation.run, group members +must have identical applied offsets, effective target offsets and interval wall +velocities. Stage each member with set_body, then call run(1) once. An incomplete +or inconsistent target group is rejected before device uploads and can be fixed. +Existing pending-target direct-stepper restrictions remain. + +Only members of that group may have overlapping whole swept bands. Their fixed +surface clearance must exceed4 lattice cells. This conservative bound exceeds +2 sqrt(2)+MAX_STEP = 3.329: M1 changed cells lie in old/new disks and read one-link +neighbors; M2 reads at most two links outside its new disk. Old/new relative +displacement is at most .5. Thus reads cannot reach another member's old/new +solid or M1 write region. M1 writes are disjoint. M2 writes opposite population +slots inside its own new disk, also disjoint. Separate first_slot ranges still +allocate each body's logical band/debug storage. Shared band cells alone do not +write anything: M1 returns unless membership changes; M2 returns unless a cut +link exists. Unrelated moving/static body and probe overlaps still reject. + +This is a sufficient restricted-group condition, not support for near-contact +or independent motion. The r4, distance8.3, dy=.5 counterexample remains rejected. +CPU tests enumerate candidate donors/write footprints over subcell phases; +GPU tests compare fixed vs hold and member-order permutations through crossings. +CUDA kernels and arithmetic are unchanged. + +## New structural parameters (not inherited from the paper) + +One transverse state y,v; total assembly M=10 rho (3 pi D^2/4)=9424.77796, +fn=U/(5D)=.0001, Tn=10000, zeta=.01, K=M(2 pi fn)^2, +C=2 zeta M(2 pi fn). Initial y=v=0, no seed or later kick. Force is the sum +of the three instantaneous CME_y+Wen_y values from the previous fluid step; +node impulses are logged separately and never added. The damped oscillator is +integrated exactly for that piecewise-constant force over dt=1. Wall velocity +is y[n+1]-y[n], not the endpoint oscillator velocity. No filtering or retuning. + +Observe600000 steps=60Tn=300 D/U after warmup, in each of free and fixed-hold +cases. Final20Tn is a descriptive comparison window; ten-Tn blocks expose drift. +No settled/periodic claim follows merely from reaching the cap. Stop on |y|>D, +nonfinite values, no_donor, rejected displacement or speed>.05. Report any quiet +response honestly. No automatic extension. --smoke uses32 warmup+32 steps only. + +Native-step rows retain y,v, each body's CME/Wen/node terms and refill counters; +each600000x15 FP64 file is72MB. No full-field output or policy training. Compile +and run on one device serially, keep original outputs and dirty trees untouched. +This stage makes no claim of physical grid convergence or vibration suppression. + +Run from checkout with PYTHONPATH=src and the existing Python environment: +`python research/steady_pinball/run.py --smoke --out ` +then separately `--out ` and `--fixed --out `. diff --git a/research/steady_pinball/config.json b/research/steady_pinball/config.json new file mode 100644 index 0000000..2baf1b9 --- /dev/null +++ b/research/steady_pinball/config.json @@ -0,0 +1,49 @@ +{ + "_doc": "Legacy geometry Re_D=50, U0=0.01, SI=800 baseline (kar_legacy_re50d_u001)", + "grid": { + "lattice_model": "D2Q9", + "nx": 1280, + "ny": 512, + "nz": 1 + }, + "physics": { + "data_type": "FP32", + "viscosity": 0.004, + "velocity": 0.01, + "rho": 1.0 + }, + "method": { + "collision": "MRT", + "streaming": "double_buffer", + "store_precision": "FP32", + "ddf_shifting": false, + "les": { + "enabled": false, + "cs": 0.16, + "closed_form": true + }, + "trt": { + "magic_param": 0.1875 + }, + "inlet": { + "profile": "parabolic", + "scheme": "regularized", + "regularized_neq_damp": 1.0 + }, + "outlet": { + "mode": "neq_extrap", + "backflow_clamp": true, + "blend_alpha": 0.7, + "srt_neq_damp": 0.5 + }, + "y_wall_bc": "bounce_back", + "omega_guard": { + "min": 0.01, + "max": 1.99 + } + }, + "cuda": { + "threads_per_block": 256, + "compute_capability": "auto" + } +} diff --git a/research/steady_pinball/run.py b/research/steady_pinball/run.py new file mode 100644 index 0000000..2af149a --- /dev/null +++ b/research/steady_pinball/run.py @@ -0,0 +1,157 @@ +"""Uncontrolled paper-geometry steady Re_D=50, one transverse assembly DOF. + +One fixed protocol; --smoke performs 32 warmup + 32 observation steps only. +See docs/steady_pinball_fiv.md for historical solver differences and new structure. +""" +from __future__ import annotations + +import argparse +import json +import math +from pathlib import Path +import subprocess +import time + +import numpy as np + +D, U, NU = 20., .01, .004 +CENTRES = ((600., 255.5), (626., 270.5), (626., 240.5)) +PROBES = ((800., 295.5), (800., 255.5), (800., 215.5)) +WARMUP, STEPS, PERIOD = 512000, 600000, 10000 +COLUMNS = ('step', 'y', 'v', 'cme_F', 'cme_U', 'cme_L', 'wen_F', 'wen_U', 'wen_L', + 'node_F', 'node_U', 'node_L', 'eq_only', 'no_donor', 'crossings') +ROOT = Path(__file__).resolve().parents[2] + + +def config(): + return json.loads(Path(__file__).with_name('config.json').read_text()) + + +class Structure: + """Exact damped oscillator update for a force held over dt=1; no force filter.""" + def __init__(self): + self.mass = 10. * 3. * math.pi * D * D / 4. + self.omega = 2. * math.pi / PERIOD + self.zeta = .01 + self.stiffness = self.mass * self.omega**2 + self.damping = 2. * self.zeta * self.mass * self.omega + self.y = self.v = 0. + a = self.zeta * self.omega + wd = self.omega * math.sqrt(1. - self.zeta**2) + self.ec = math.exp(-a) * math.cos(wd) + self.es = math.exp(-a) * math.sin(wd) / wd + self.a = a + + def candidate(self, force): + equilibrium = force / self.stiffness + z = self.y - equilibrium + return (equilibrium + (self.ec + self.a*self.es)*z + self.es*self.v, + -self.omega**2*self.es*z + (self.ec - self.a*self.es)*self.v) + + +def summarize(rows): + windows = [] + for start in range(0, len(rows), 10*PERIOD): + block = rows[start:start + 10*PERIOD] + fy = block[:, 3:9].sum(axis=1) + windows.append(dict(start_step=int(block[0, 0]), end_step=int(block[-1, 0]), + y_mean_D=float(block[:, 1].mean()/D), + y_rms_D=float(np.sqrt(np.mean(block[:, 1]**2))/D), + y_peak_D=float(np.abs(block[:, 1]).max()/D), + fy_mean=float(fy.mean()), fy_rms=float(np.sqrt(np.mean(fy**2))), + mean_power=float(np.mean(fy*block[:, 2])))) + return windows + + +def main(): + ap = argparse.ArgumentParser(description=__doc__) + ap.add_argument('--out', required=True, type=Path) + ap.add_argument('--fixed', action='store_true', help='same modern moving-path geometry held fixed') + ap.add_argument('--smoke', action='store_true', help='32+32 steps, no physical result') + args = ap.parse_args() + out = args.out.resolve() + if out == ROOT or ROOT in out.parents: + ap.error('output must be outside the checkout') + out.mkdir(parents=True, exist_ok=False) + cfg = config() + (out/'config.json').write_text(json.dumps(cfg, indent=2)+'\n') + source = subprocess.check_output(['git', '-C', str(ROOT), 'rev-parse', 'HEAD'], text=True).strip() + s = Structure() + warmup, count = (32, 32) if args.smoke else (WARMUP, STEPS) + rows = np.lib.format.open_memmap(out/'steps.npy', mode='w+', dtype=np.float64, + shape=(count, len(COLUMNS))) + rows[:] = np.nan + from CelerisLab import Simulation + from CelerisLab.body.motion import MovingTables, Translation + names = MovingTables.DIAG + col = {k: names.index(k) for k in names} + status, completed, f_last = 'running', 0, 0. + started = time.monotonic() + report = dict(source=source, fixed=args.fixed, smoke=args.smoke, columns=COLUMNS, + structure=dict(mass=s.mass, stiffness=s.stiffness, damping=s.damping, + zeta=s.zeta, fn=1./PERIOD, y0=0., v0=0.), + warmup_steps=warmup, requested_steps=count, dt=1., D=D, U=U, nu=NU, + duration_Tn=count/PERIOD, duration_tU_D=count*U/D, + rhs='sum of three instantaneous CME_y + Wen_y, one-step lag, no filtering') + try: + with Simulation(str(out/'config.json')) as sim: + # Match Legacy object order; probes remain passive and fixed in the lab. + for c in PROBES: + sim.add_body('sensor', c, 5.) + envelope = Translation((0., 0., -D, D), .05, shared=True) + bodies = [sim.add_body('circle', c, D/2, translation=envelope) for c in CENTRES] + sim.initialize() + # At release obtain the last instantaneous force, not a warmup average. + remaining = warmup - 1 + while remaining: + chunk = min(10000, remaining) + sim.run(chunk) + remaining -= chunk + sim.run(1) + d = sim.moving_diagnostics() + f_last = float(d[:, col['cme_y']].sum(dtype=np.float64) + d[:, col['wen_y']].sum(dtype=np.float64)) + report['warmup_elapsed_s'] = time.monotonic()-started + print(json.dumps(dict(event='release', **report)), flush=True) + for n in range(count): + y1, v1 = (0., 0.) if args.fixed else s.candidate(f_last) + if not all(map(math.isfinite, (y1, v1))) or abs(y1) > D: + status = 'stopped_nonfinite_or_travel_bound' + break + if not args.fixed: + for bid in bodies: + sim.set_body(bid, offset=(0., y1), vy=y1-s.y) + sim.run(1) + d = sim.moving_diagnostics() + cme, wen = d[:, col['cme_y']], d[:, col['wen_y']] + node = d[:, col['cover_jy']] - d[:, col['uncover_jy']] + row = (n+1, y1, v1, *cme, *wen, *node, + d[:, col['eq_only']].sum(), d[:, col['no_donor']].sum(), d[:, col['crossings']].sum()) + if not np.isfinite(row).all() or row[-2] != 0: + status = 'stopped_nonfinite_or_no_donor' + break + rows[n] = row + completed = n+1 + s.y, s.v = y1, v1 + f_last = float(cme.sum(dtype=np.float64)+wen.sum(dtype=np.float64)) + if completed % PERIOD == 0: + rows.flush() + print(json.dumps(dict(event='period', step=completed, t_Tn=completed/PERIOD, + y_D=s.y/D, elapsed_s=time.monotonic()-started)), flush=True) + else: + status = 'smoke_completed' if args.smoke else 'duration_completed_not_qualified' + except BaseException as exc: + status = f'error: {type(exc).__name__}: {exc}' + raise + finally: + rows.flush() + report.update(status=status, completed_steps=completed, elapsed_s=time.monotonic()-started, + windows=summarize(rows[:completed]), + totals={key: float(rows[:completed, i].sum()) for key, i in + [('eq_only', 12), ('no_donor', 13), ('crossings', 14)]}, + limitation='Uncontrolled development response only. No automatic settling, grid convergence, control or historical-solver equivalence claim.') + (out/'report.json').write_text(json.dumps(report, indent=2)+'\n') + print(json.dumps(report), flush=True) + + +if __name__ == '__main__': + main() diff --git a/src/CelerisLab/body/motion.py b/src/CelerisLab/body/motion.py index c82e25b..b55f0f1 100644 --- a/src/CelerisLab/body/motion.py +++ b/src/CelerisLab/body/motion.py @@ -14,6 +14,7 @@ from typing import Callable, Optional, Tuple MAX_STEP = 0.5 # Euclidean centre displacement per step (lattice cells) BAND_MARGIN = 2 # cells added around the swept circle AABB +SHARED_CLEARANCE = 4.0 # > two diagonal links + MAX_STEP; docs/steady_pinball_fiv.md R_MIN = 4.0 # smallest radius covered by the donor enumeration Vec = Tuple[float, float] @@ -31,6 +32,7 @@ class Translation: """Declared motion envelope: offset box (x_min, x_max, y_min, y_max) and speed bound.""" offset_box: Tuple[float, float, float, float] max_speed: float + shared: bool = False # same envelope instance: synchronous transverse motion def __post_init__(self): box = tuple(float(v) for v in self.offset_box) @@ -41,10 +43,29 @@ class Translation: raise ValueError("offset_box must contain the zero offset") if not 0.0 < speed <= MAX_STEP: raise ValueError(f"max_speed must satisfy 0 < max_speed <= {MAX_STEP}") + if self.shared and box[:2] != (0.0, 0.0): + raise ValueError("shared translation is transverse only (x offsets must be zero)") object.__setattr__(self, "offset_box", box) object.__setattr__(self, "max_speed", speed) +def validate_shared_translations(circles): + """Check common old/target offsets and wall velocities before device work. + + Per-body setters may stage members in any order. Envelope identity (not + value equality) defines a group; a missing target means hold at applied. + """ + groups = {} + for motion in circles: + if not motion.translation.shared: + continue + target = motion.pending if motion.pending is not None else (motion.applied, (0.0, 0.0)) + state = (motion.applied, target) + key = id(motion.translation) + if groups.setdefault(key, state) != state: + raise ValueError("shared translation needs identical applied/target offsets and velocities") + + class MovingCircle: """Applied offset, pending target and the fail-closed commit of one circle.""" @@ -72,6 +93,8 @@ class MovingCircle: x0, x1, y0, y1 = self.translation.offset_box if not (x0 <= offset[0] <= x1 and y0 <= offset[1] <= y1): raise ValueError(f"offset {offset} outside declared box {self.translation.offset_box}") + if self.translation.shared and velocity[0] != 0.0: + raise ValueError("shared translation needs zero x velocity") step = math.hypot(offset[0] - self.applied[0], offset[1] - self.applied[1]) if step > MAX_STEP: raise ValueError(f"per-step displacement {step} exceeds {MAX_STEP}") diff --git a/src/CelerisLab/simulation.py b/src/CelerisLab/simulation.py index a7e1195..600e871 100644 --- a/src/CelerisLab/simulation.py +++ b/src/CelerisLab/simulation.py @@ -259,6 +259,8 @@ class Simulation: moving = self._moving_bodies() if not moving: return + from .body.motion import SHARED_CLEARANCE, validate_shared_translations + validate_shared_translations([obj.motion for obj in moving]) cfg = self.lbm_cfg if (not cfg.is_d2q9 or cfg.streaming != "double_buffer" or cfg.ddf_shifting or cfg.store_precision != "FP32"): @@ -277,6 +279,13 @@ class Simulation: if other is obj: continue (cx, cy), r = other.center[:2], float(other.radius) + if (other.motion is not None and obj.motion.translation.shared + and obj.motion.translation is other.motion.translation): + # Whole envelopes may overlap; simultaneous local stencils may not. + gap = np.hypot(obj.center[0] - cx, obj.center[1] - cy) - obj.radius - r + if gap <= SHARED_CLEARANCE: + raise ValueError("shared translation needs surface clearance > 4 cells") + continue if other.motion is not None: ox0, ox1, oy0, oy1 = other.motion.band() else: @@ -719,6 +728,9 @@ class Simulation: if steps < 0 or checkpoint_interval < 0: raise ValueError("steps and checkpoint_interval must be nonnegative") moving = self._moving_bodies() + if steps > 0 and moving: + from .body.motion import validate_shared_translations + validate_shared_translations([obj.motion for obj in moving]) if steps > 0 and moving and self.stepper.moving is None: raise RuntimeError(MOVING_BLOCKED) if moving and checkpoint_interval > 0: diff --git a/tests/integration/test_shared_translation_gpu.py b/tests/integration/test_shared_translation_gpu.py new file mode 100644 index 0000000..090d9e2 --- /dev/null +++ b/tests/integration/test_shared_translation_gpu.py @@ -0,0 +1,60 @@ +"""Small serial GPU gates for shared transverse circles; no physical benchmark.""" +import json +from pathlib import Path + +import numpy as np + +from CelerisLab import Simulation +from CelerisLab.body.motion import MovingTables, Translation + +CENTRES = ((70., 79.5), (96., 94.5), (96., 64.5)) + + +def simulate(tmp_path, *, moving, offsets, reverse=False): + config = json.loads((Path(__file__).resolve().parents[2]/'research/steady_pinball/config.json').read_text()) + config['grid'].update(nx=192, ny=160) + p = tmp_path/'config.json' + p.write_text(json.dumps(config)) + with Simulation(str(p)) as sim: + envelope = Translation((0., 0., -20., 20.), .5, shared=True) + centres = list(reversed(CENTRES)) if reverse else list(CENTRES) + ids = [sim.add_body('circle', c, 10., translation=envelope if moving else None) for c in centres] + sim.initialize() + previous = 0. + crossings = 0. + for y in offsets: + if moving: + for bid in ids: + sim.set_body(bid, offset=(0., float(y)), vy=float(y-previous)) + sim.run(1) + if moving: + diag = sim.moving_diagnostics() + assert np.isfinite(diag).all() + assert np.all(diag[:, MovingTables.DIAG.index('no_donor')] == 0) + crossings += diag[:, MovingTables.DIAG.index('crossings')].sum() + previous = y + forces = np.array([sim.read_force(bid) for bid in ids]) + if reverse: + forces = forces[::-1] + return sim.get_ddf(), sim.get_flags(), forces, crossings + + +def test_fixed_vs_shared_hold(tmp_path): + fixed = simulate(tmp_path, moving=False, offsets=[0.]*32) + held = simulate(tmp_path, moving=True, offsets=[0.]*32) + assert np.array_equal(fixed[0], held[0], equal_nan=True) + assert np.array_equal(fixed[1], held[1]) + # FP32 atomic force reduction order only; populations remain exactly equal. + np.testing.assert_allclose(fixed[2], held[2], rtol=0., atol=1e-4) + assert held[3] == 0 + + +def test_common_motion_member_order_has_no_population_effect(tmp_path): + offsets = np.concatenate((np.arange(.5, 20.5, .5), np.arange(19.5, -20.5, -.5), + np.arange(-19.5, .5, .5))) + forward = simulate(tmp_path, moving=True, offsets=offsets) + reverse = simulate(tmp_path, moving=True, offsets=offsets, reverse=True) + assert forward[3] > 0 and forward[3] == reverse[3] + assert np.array_equal(forward[0], reverse[0], equal_nan=True) + assert np.array_equal(forward[1], reverse[1]) + np.testing.assert_allclose(forward[2], reverse[2], rtol=0., atol=1e-4) diff --git a/tests/unit/test_shared_translation.py b/tests/unit/test_shared_translation.py new file mode 100644 index 0000000..2155326 --- /dev/null +++ b/tests/unit/test_shared_translation.py @@ -0,0 +1,122 @@ +"""Host admission and independent lattice-stencil checks; no CUDA context.""" +from types import SimpleNamespace + +import numpy as np +import pytest + +from CelerisLab.body.motion import MovingCircle, MovingTables, Translation, validate_shared_translations +from CelerisLab.lbm.descriptors import FLUID +from tests.unit.test_moving_bodies import STREAM, body, moving_simulation, routed_simulation + + +def scene(monkeypatch, envelope=None, gap=10.): + sim = routed_simulation(monkeypatch) + envelope = envelope or Translation((0., 0., -6., 6.), .5, shared=True) + for y in (40. - (20. + gap) / 2., 40. + (20. + gap) / 2.): + sim.bodies.add(body(60., y, 10., envelope)) + sim.field = SimpleNamespace(nx=120, ny=80, build_channel_flags=lambda: np.full((80, 120), FLUID)) + return sim + + +def test_shared_bands_are_admitted_but_other_groups_are_not(monkeypatch): + sim = scene(monkeypatch) + sim._validate_moving_bodies() + a, b = [o.motion for o in sim.bodies.objects] + assert a.band()[3] >= b.band()[2] # whole envelopes genuinely overlap + band, geom = MovingTables.layout([(a, 0), (b, 1)]) + assert band[10] == band[2] * band[3] * 8 + assert len(geom) == 14 + b.translation = Translation((0., 0., -6., 6.), .5, shared=True) + with pytest.raises(ValueError, match='overlaps body'): + sim._validate_moving_bodies() # equal values do not link independent motion + b.translation = a.translation + sim.bodies.add(body(60., 40., 2.)) + with pytest.raises(ValueError, match='overlaps body'): + sim._validate_moving_bodies() + + +@pytest.mark.parametrize('gap', (.3, 4.)) +def test_small_surface_gaps_rejected(monkeypatch, gap): + with pytest.raises(ValueError, match='clearance'): + scene(monkeypatch, gap=gap)._validate_moving_bodies() + + +def test_shared_envelope_is_transverse_only(): + with pytest.raises(ValueError, match='transverse'): + Translation((-.1, .1, -1., 1.), .5, shared=True) + + envelope = Translation((0., 0., -1., 1.), .5, shared=True) + with pytest.raises(ValueError, match="zero x velocity"): + MovingCircle((0., 0.), 10., envelope).check_target((0., .1), (.1, .1)) + +def test_group_rejection_is_before_writes_and_can_be_corrected(monkeypatch): + sim, writes = moving_simulation(monkeypatch) + env = Translation((0., 0., -2., 2.), .5, shared=True) + sim.bodies.objects = [body(40., y, 10., env) for y in (20., 50.)] + a, b = [o.motion for o in sim.bodies.objects] + sim.stepper.moving = object() + a.set_target((0., .2), (0., .2)) + with pytest.raises(ValueError, match='identical'): + sim.run(1, stream=STREAM) + assert writes == [] and a.applied == b.applied == (0., 0.) + b.set_target((0., .2), (0., .1)) + with pytest.raises(ValueError, match='identical'): + sim.run(1, stream=STREAM) + assert writes == [] + b.set_target((0., .2), (0., .2)) + validate_shared_translations([a, b]) + b.applied = (0., -.1) + with pytest.raises(ValueError, match='identical'): + validate_shared_translations([a, b]) + + +C = np.array([(x, y) for x in (-1, 0, 1) for y in (-1, 0, 1) if (x, y) != (0, 0)]) + + +def footprint_conflicts(centres, radius, dy): + """Enumerate geometric write/read supports, not refill/BFL formulas. + + Tests every candidate donor, conservatively including M2's second link. + Body-indexed slot numbers are separate from physical (cell,direction) writes. + """ + centres = np.asarray(centres, float) + lo = np.floor(centres.min(0) - radius - 4).astype(int) + hi = np.ceil(centres.max(0) + radius + 4).astype(int) + yy, xx = np.mgrid[lo[1]:hi[1]+1, lo[0]:hi[0]+1] + points = np.column_stack((xx.ravel(), yy.ravel())) + def inside(p, c): + return ((p - c)**2).sum(axis=-1) < radius**2 + old = np.array([inside(points, c) for c in centres]) + new_centres = centres + (0., dy) + new = np.array([inside(points, c) for c in new_centres]) + changed = old ^ new + conflicts = int(np.count_nonzero(changed.sum(0) > 1)) + writes = set() + for i, (oc, nc) in enumerate(zip(centres, new_centres)): + others = [j for j in range(len(centres)) if j != i] + fresh = points[old[i] & ~new[i]] + for e in C: + donor = fresh + e + valid = ~inside(donor, oc) & ~inside(donor, nc) + for j in others: + conflicts += int(np.count_nonzero(valid & (inside(donor, centres[j]) | inside(donor, new_centres[j])))) + link = ~new[i] & inside(points + e, nc) + fluid = points[link] + solid = fluid + e + for p in solid: + slot = (int(p[0]), int(p[1]), int(e[0]), int(e[1])) + conflicts += int(slot in writes) + writes.add(slot) + for read in (fluid, fluid - e): + for j in others: + conflicts += int(np.count_nonzero(inside(read, centres[j]) | inside(read, new_centres[j]))) + return conflicts + + +def test_counterexample_and_pinball_local_ownership(): + # Close but disjoint circles are insufficient: parent/Grok independent example. + assert footprint_conflicts(((.25, .25), (.25, 8.55)), 4., .5) > 0 + for phase in (0., .25, .5, .75): + for dy in (-.5, -.1, .1, .5): + centres = np.array(((0., 0.), (26., 15.), (26., -15.))) + phase + assert footprint_conflicts(centres, 10., dy) == 0 diff --git a/tests/unit/test_steady_pinball.py b/tests/unit/test_steady_pinball.py new file mode 100644 index 0000000..dffd806 --- /dev/null +++ b/tests/unit/test_steady_pinball.py @@ -0,0 +1,54 @@ +"""Real config loading and pure structural update for the frozen development case.""" +import importlib.util +import math +from pathlib import Path + +import numpy as np +import pytest + +path = Path(__file__).resolve().parents[2] / 'research/steady_pinball/run.py' +spec = importlib.util.spec_from_file_location('steady_pinball_run', path) +runner = importlib.util.module_from_spec(spec) +spec.loader.exec_module(runner) + + +def test_actual_loader_and_paper_geometry(tmp_path): + import json + from CelerisLab.config import load_lbm_config + p = tmp_path/'config.json' + p.write_text(json.dumps(runner.config())) + cfg = load_lbm_config(str(p)) + cfg.validate() + assert (cfg.nx, cfg.ny, cfg.viscosity, cfg.velocity) == (1280, 512, .004, .01) + assert cfg.collision == 'MRT' and cfg.y_wall_bc == 'bounce_back' + assert runner.U * runner.D / runner.NU == 50. + assert runner.CENTRES == ((600., 255.5), (626., 270.5), (626., 240.5)) + + +def test_structure_rest_equilibrium_and_energy_decay(): + s = runner.Structure() + assert s.candidate(0.) == (0., 0.) + assert s.mass == pytest.approx(10*3*math.pi*100) + s.y, s.v = .2, 0. + # One damped period: independent continuous solution at total elapsed time. + for _ in range(runner.PERIOD): + s.y, s.v = s.candidate(0.) + t = runner.PERIOD + a = s.zeta*s.omega + w = s.omega*math.sqrt(1-s.zeta**2) + expected = .2*math.exp(-a*t)*(math.cos(w*t)+a/w*math.sin(w*t)) + assert s.y == pytest.approx(expected, abs=2e-12) + assert .5*s.mass*s.v*s.v+.5*s.stiffness*s.y*s.y < .5*s.stiffness*.2**2 + s.y, s.v = .3, 0. + assert s.candidate(s.stiffness*.3) == pytest.approx((.3, 0.), abs=1e-15) + + +def test_report_uses_completed_rows_and_net_force(): + rows = np.zeros((4, 15)) + rows[:, 0] = np.arange(1, 5) + rows[:, 1] = 2. + rows[:, 3:9] = 1. + got = runner.summarize(rows) + assert got[0]['end_step'] == 4 + assert got[0]['y_rms_D'] == .1 + assert got[0]['fy_rms'] == 6.