diff --git a/README.md b/README.md index 191cb5a..7a59a88 100644 --- a/README.md +++ b/README.md @@ -665,6 +665,7 @@ streak.render("streakline.png") ## Validation campaigns +- [Periodic Taylor–Green AMR](docs/validation/taylor_green.md): analytic periodic-decay pilot and frozen baseline-relative assessment. - [Matched refined DFG](docs/validation/refined_dfg.md): frozen SRT uniform/refined comparison, requiring channel GPU gates before execution. - [DFG 2D-2 Re100](docs/validation/dfg_re100.md): frozen periodicity and official-reference comparisons, with a declared noncanonical outlet. - [Uniform LES shear decay](docs/validation/uniform_les_flow.md): independent continuum-oracle flow diagnostic, baseline and isolated corrected source variants. diff --git a/docs/validation/taylor_green.md b/docs/validation/taylor_green.md new file mode 100644 index 0000000..d710f4d --- /dev/null +++ b/docs/validation/taylor_green.md @@ -0,0 +1,161 @@ +# Periodic 2D Taylor–Green AMR campaign (prepared, not GPU-qualified) + +This validation harness complements DFG Re100 and Jet/body tests. It measures +smooth 2D periodic decay and errors crossing fixed refinement interfaces at +Re_L=100 and 1000. It does **not** certify 3D turbulence, LES, Jet performance, +wall conditions, dynamic refinement or general high-Re flows. + +## Reference and conventions + +Primary implementation-author documentation, checked 2026-10-18: +- TensorMesh, exact velocity and pressure with a consistent phase: + https://docs.tensor-mesh.com/example_gallery/fluid/taylor_green.html +- AeroSim Nassu, 2D decay, energy and lattice-length Re definition: + https://docs.aerosim.io/nassu/validation/analytical/04_taylor_green_vortex_2d/README.html + +We use the TensorMesh phase, rescaled to length L=N root cells: + + k = 2*pi/L; A(t) = U*exp(-2*nu*k*k*t) + ux = -A*cos(k*x)*sin(k*y) + uy = A*sin(k*x)*cos(k*y) + p/rho0 = -A*A/4 * (cos(2*k*x)+cos(2*k*y)) + E = <|u|^2>/2 = U*U/4 * exp(-4*nu*k*k*t) + +Pressure is kinematic, with rho0=1. A global velocity sign reversal leaves +pressure unchanged; exchanging sine/cosine without the corresponding pressure +phase can reverse its sign. CPU tests independently check momentum residuals. +The nonlinear acceleration is a pressure gradient, **not zero**; the AeroSim +page's prose calling advection zero is not adopted. Do not mix phase variants +from different references. The formulas above are verified by substitution, +not by trusting a copied pressure expression. + +Re_L=U*L/nu, not Re_k=U/(nu*k)=Re_L/(2*pi). U=.03 in lattice units, +Ma=U/cs=.0519615, cs^2=1/3, nu_root=.03*N/Re_L. Acoustic refinement uses +h_l=dt_l=2^-l: lattice velocity stays U, nu_l=2^l*nu_root and + tau_l=.5+3*nu_l. The worst root tau is .50288 (N32/Re1000), omega=1.98855, +inside the existing config guard. No viscosity clamping or guard override. + +All initial arrays are physical FP64 equilibrium populations, supplied to the +public shifted-FP32 postcollision-G initializer. rho=1+3p at each true cell +center. These are valid equilibrium postcollision states, but do not include +Chapman–Enskog nonequilibrium stress; startup acoustic/kinetic error remains in +the measurements and is shared with uniform baselines. The incompressible +reference is not an exact finite-Mach LBM solution. This is why resolution alone +need not remove the full error floor. + +## Geometry and workload + +All patches use level-wide union form, so origins are zero, including fine +arrays with holes. Coordinates are (index+.5)*2^-level. Metrics integrate only +ACTIVE leaves, each weighted by 4^-level. Covered cells and halo cells never +enter norms or mass. Geometry tests prove exact single coverage at root-cell +resolution. All blocks are 8x8, and both root dimensions are periodic. + +At each root N=32,64: +- 2 levels: L1 rectangle (N/8,N/8,7N/8,3N/4), in L0 indices. +- 3 levels: same L1 plus L2 rectangle (3N/4,N/2,3N/2,5N/4), in L1 indices. + +These off-axis boundaries cross nonzero flow and gradients rather than sitting +only on symmetry or stagnation lines. Whole geometry scales with N. + +16 cases: at each Re, uniform N=32,64,128,256 plus AMR (root32/root64)x(depth2/depth3). +Each AMR case therefore has both a same-root uniform baseline and a uniform +baseline at its finest spacing. Uniform N is its own lattice length, so nu and +step count are scaled to represent the same nondimensional solution. +Samples t*=tU/L = 0,.12,.24,.36,.48. End time is 16*N base steps, a meaningful +.48 domain turnover or 3.016 vortex-scale turnover times (scale 1/k). +At Re1000 energy decays about 7.3%; at Re100 about 53.1%. No plots or full fields; +each case writes five diagnostic rows and provenance, incrementally. + +## Running (GPU owner must schedule launches) + +CPU checks, from the repository root: + + PYTHONPATH=src python -m unittest discover -s tests/validation/taylor_green -p 'test_*.py' -v + +One GPU case, on the scheduled worker after checking occupancy: + + CUDA_VISIBLE_DEVICES= python tests/validation/taylor_green/run_case.py --repo /path/to/reviewed/checkout --re 100 --n 32 --depth 1 --output output/taylor_green/re100_n32_d1.json + +Depth 1 is uniform. Run uniform cases before AMR for each Re; inspect baseline +screening failures before proceeding. For a complete matrix, repeat the uniform +call with N=32,64,128,256 and AMR calls with N=32,64 / depth=2,3, at each Re. +Use serial process startup on the parent's assigned Horsea GPUs. This harness +does not schedule or acquire GPUs. No GPU was used during preparation. + + python tests/validation/taylor_green/assess.py output/taylor_green/re*_n*_d*.json --output output/taylor_green/assessment.json + +The local API/geometry review used clean staging commit +cd3a657f9d3833e812c3cca173525f165b48710a. This is **not** a claim of results on d79 +or any historical commit. Each runtime output records actual Git HEAD/dirty +state, an aggregate source-tree hash, plus SHA256 of all harness Python files. Assessment refuses dirty checkouts and mixed source +or harness identities. Keep reports outside the checkout or .gitignore them so +new outputs do not change its dirty flag between cases. The source-tree hash excludes bytecode caches; production qualification requires +a clean checkout. +The runner uses the existing PeriodicBlockStepper / FixedPatchStepper periodic +paths only. It does not depend on the new channel-root implementation. + +## Preregistered assessment, before any GPU results + +All sampled times must satisfy finite screening ceilings: +- Relative composite velocity L2 <= 5% +- Relative globally gauge-adjusted pressure L2 <= 15% +- Relative kinetic-energy error <= 5% +- Composite mass drift from downloaded t0 <= 1e-5 + +These are screening budgets, not literature-certified tolerances. Acceptance +also requires all 16 cases, common source/harness, correct nondimensional times +and viscosity, and baseline-relative checks on worst post-startup-sample errors: +- AMR error <= max(2 * same-root uniform error, declared numerical floor). +- Root64 AMR error <= max(.9 * root32 AMR error, floor), for each depth/Re. +- Uniform256 error <= max(.9 * uniform32 error, floor). + +Floors: velocity .002, pressure .02, energy .004. They are declared low-Mach / +equilibrium-startup budgets, not fitted to observed results. Ma^2=.0027 motivates +their scale; it is an asymptotic expectation, not a proof of those coefficients. +The coarse baseline sets the AMR budget because much of the domain remains +coarse. The finest uniform error and AMR/fine ratio are always reported, with no +claim that partial refinement should match a fully fine domain. A pass below a +floor is bounded accuracy, not a measured convergence order. Report errors and +ratios even for a pass. Failure requires diagnosis or a separately preregistered +follow-up, never retroactive tolerance tuning. Higher Re is deliberately absent; +Re100/1000 must first be stable, resolved and reviewed. + +Diagnostics use one global pressure gauge, direct cell-center analytic values, +kinetic energy per unit reference mass, and actual leaf cell areas. Composite +midpoint quadrature introduces a small initial energy integration error, retained +rather than subtracted. Mass drift is relative to downloaded t0, so density +initialization and mixed-grid quadrature do not masquerade as time drift. + +## Pilot scheduling and transfer-consistency scope + +Start with one uniform case and one matched AMR case; the 16-case matrix is a +qualification plan, not authorization to launch every case automatically. The +GPU owner reviews the pilot before expanding it. Each process prints a flushed +READY only after initialization, a real base step exercising every level's +kernel/halo/transfer path, and a finite positive-density download check. That +step belongs to the measured trajectory; it is not discarded warmup. PROGRESS +messages contain each of the five diagnostic samples. READY alone is not an +accuracy pass. Reports include per-level active cells, hops and collision updates. + +The complete matrix costs 768,540,672 active collision updates and 21,504 root +intervals, excluding ghost streaming, halos, transfers and diagnostic work. +Uniform N256 is the largest case (268,435,456 updates / 4096 intervals per Re). +AMR N64/depth3 needs 50,987,008 updates / 1024 intervals. These are work counts, +not GPU wall-time estimates. Baseline comparisons remain available on failed +screens and partial pilots; absent/nonfinite operands are null, never passing. + +Acoustic viscosity scaling and unchanged lattice velocity do not by themselves +prove nonequilibrium-stress consistency across refinement interfaces. A smooth +collocated Chapman–Enskog postcollision nonequilibrium has a (tau-1)*h factor; +the existing implementation instead uses volumetric piecewise-constant explosion +and stream-only ghosts. Whether that actual transfer construction achieves the +needed stress accuracy is an unresolved question for this pilot to measure. +It is not a proven defect, and a collocated-state scaling factor must not be +inserted into this different transfer operator without a separate derivation +and validation. This campaign measures the current operator and makes no +theoretical stress-consistency certification or production-kernel change. + +## Repository integration + +The four Python harness files remain together under tests/validation/taylor_green, so runtime provenance hashes precisely this harness directory. Formatting changes produce a new harness identity; preparation and integrated receipts must not be mixed. All analytic definitions, numerical floors, assessment gates, sample times, and source-import assertions are unchanged. Production source is unchanged from the cd3 snapshot; these periodic tests do not waive the separate channel GPU gate. diff --git a/tests/validation/taylor_green/assess.py b/tests/validation/taylor_green/assess.py new file mode 100644 index 0000000..cffc9fe --- /dev/null +++ b/tests/validation/taylor_green/assess.py @@ -0,0 +1,124 @@ +"""Frozen assessment of the complete 16-case matrix, independent of GPU launch.""" + +import argparse +import json +import math +from pathlib import Path + +from oracle import TIMES, LIMITS, U + +# Numerical floors reflect finite-Mach/FP32 and equilibrium startup, not observed runs. +FLOORS = dict(velocity_l2=0.002, pressure_l2=0.02, energy_relative=0.004) +BASELINE_FACTOR = 2.0 +REFINEMENT_FACTOR = 0.9 + + +def assess(records): + indexed = {} + for r in records: + key = (r["re_L"], r["n"], r["depth"]) + if key in indexed: + raise ValueError(f"duplicate case {key}") + indexed[key] = r + expected = {(re, n, 1) for re in (100, 1000) for n in (32, 64, 128, 256)} + expected |= {(re, n, d) for re in (100, 1000) for n in (32, 64) for d in (2, 3)} + missing = sorted(expected - indexed.keys()) + checks = { + "complete_matrix": not missing, + "no_extra_cases": set(indexed) <= expected, + } + identities = {json.dumps(r["source"], sort_keys=True) for r in records} + checks["same_source_and_harness"] = len(identities) == 1 + checks["clean_source"] = bool(records) and all( + not r["source"]["dirty"] for r in records + ) + for key, r in indexed.items(): + label = str(key) + rows = r["rows"] + checks[label + " complete"] = r["completed"] and len(rows) == len(TIMES) + checks[label + " scales"] = ( + r["U"] == U and abs(r["viscosity"] - U * r["n"] / r["re_L"]) < 1e-15 + ) + checks[label + " times"] = len(rows) == len(TIMES) and all( + abs(v["tstar"] - t) < 1e-12 for v, t in zip(rows, TIMES) + ) + checks[label + " absolute"] = bool(rows) and all( + 0 <= row[k] <= limit for row in rows for k, limit in LIMITS.items() + ) + comparisons = [] + + def worst(key, metric): + values = [r[metric] for r in indexed.get(key, {}).get("rows", [])[1:]] + return max(values) if values and all(math.isfinite(v) for v in values) else None + + # Diagnostics survive failed screens and partial pilot matrices. Missing or + # nonfinite operands are explicit nulls; they cannot acquire a passing gate. + for re in (100, 1000): + for depth in (2, 3): + for n in (32, 64): + if (re, n, depth) not in indexed: + continue + for metric, floor in FLOORS.items(): + amr = worst((re, n, depth), metric) + coarse = worst((re, n, 1), metric) + fine = worst((re, n * 2 ** (depth - 1), 1), metric) + ok = ( + amr is not None + and coarse is not None + and amr <= max(BASELINE_FACTOR * coarse, floor) + ) + comparisons.append( + dict( + re_L=re, + n=n, + depth=depth, + metric=metric, + amr=amr, + uniform_coarse=coarse, + uniform_finest=fine, + amr_to_finest=( + None + if amr is None or fine is None + else amr / max(fine, 1e-30) + ), + passed=ok, + ) + ) + checks[f"baseline {re}/{n}/{depth}/{metric}"] = ok + for metric, floor in FLOORS.items(): + low = worst((re, 32, depth), metric) + high = worst((re, 64, depth), metric) + checks[f"refinement {re}/{depth}/{metric}"] = ( + low is not None + and high is not None + and high <= max(REFINEMENT_FACTOR * low, floor) + ) + for metric, floor in FLOORS.items(): + low = worst((re, 32, 1), metric) + high = worst((re, 256, 1), metric) + checks[f"uniform refinement {re}/{metric}"] = ( + low is not None + and high is not None + and high <= max(REFINEMENT_FACTOR * low, floor) + ) + return dict( + passed=all(checks.values()), + missing=missing, + checks=checks, + comparisons=comparisons, + floors=FLOORS, + baseline_factor=BASELINE_FACTOR, + refinement_factor=REFINEMENT_FACTOR, + scope="2D periodic smooth-decay AMR verification only; no turbulence, LES, wall or Jet certification", + ) + + +if __name__ == "__main__": + p = argparse.ArgumentParser(description=__doc__) + p.add_argument("files", type=Path, nargs="+") + p.add_argument("--output", type=Path, required=True) + a = p.parse_args() + result = assess([json.loads(path.read_text()) for path in a.files]) + a.output.write_text(json.dumps(result, indent=2, allow_nan=False) + "\n") + print(json.dumps(result, indent=2)) + raise SystemExit(0 if result["passed"] else 1) diff --git a/tests/validation/taylor_green/oracle.py b/tests/validation/taylor_green/oracle.py new file mode 100644 index 0000000..8e82ab6 --- /dev/null +++ b/tests/validation/taylor_green/oracle.py @@ -0,0 +1,103 @@ +"""CPU-only 2D Taylor–Green reference and composite-leaf diagnostics.""" + +import numpy as np + +CX = np.array([0, 1, -1, 0, 0, 1, -1, 1, -1]) +CY = np.array([0, 0, 0, 1, -1, 1, -1, -1, 1]) +W = np.array([4 / 9] + [1 / 9] * 4 + [1 / 36] * 4) +U = 0.03 +HORIZON = 0.48 # t* = t U/L; exactly 16*N base steps +TIMES = [0.0, 0.12, 0.24, 0.36, 0.48] +# Frozen screening ceilings; baseline-relative qualification is in assess.py. +LIMITS = dict( + velocity_l2=0.05, pressure_l2=0.15, energy_relative=0.05, mass_relative=1e-5 +) + + +def exact(x, y, t, length, viscosity, amplitude=U): + k = 2 * np.pi / length + a = amplitude * np.exp(-2 * viscosity * k * k * t) + # TensorMesh phase: negative cos(x)sin(y), positive sin(x)cos(y). + ux = -a * np.cos(k * x) * np.sin(k * y) + uy = a * np.sin(k * x) * np.cos(k * y) + p = -0.25 * a * a * (np.cos(2 * k * x) + np.cos(2 * k * y)) + return ux, uy, p # kinematic pressure, mean zero + + +def centers(shape, level=0): + # Campaign deliberately uses level-wide union layouts, all origins zero. + y, x = np.indices(shape, dtype=float) + h = 2.0 ** (-level) + return h * (x + 0.5), h * (y + 0.5) + + +def initial(shape, level, length, viscosity): + x, y = centers(shape, level) + ux, uy, p = exact(x, y, 0.0, length, viscosity) + rho = 1 + 3 * p # cs^2=1/3, rho0=1 + cu = CX[:, None, None] * ux + CY[:, None, None] * uy + # Exact FP64 weights; public shifted initializer encodes before FP32 rounding. + return ( + W[:, None, None] + * rho + * (1 + 3 * cu + 4.5 * cu * cu - 1.5 * (ux * ux + uy * uy)) + ) + + +def rectangles(n, depth): + if n < 32 or n % 32 or depth not in (2, 3): + raise ValueError("AMR requires N >= 32 divisible by 32 and depth 2 or 3") + # Broad off-axis rectangle cuts nonzero velocity and both gradients. + out = [[(n // 8, n // 8, 7 * n // 8, 3 * n // 4)]] + if depth == 3: + # L1 index coordinates; global physical bounds (3N/8,N/4,3N/4,5N/8). + out.append([(3 * n // 4, n // 2, 3 * n // 2, 5 * n // 4)]) + return out + + +def metrics(arrays, masks, t, length, viscosity, initial_mass=None): + accum = [] + for level, (g, mask) in enumerate(zip(arrays, masks)): + if g.shape != (9, *mask.shape): + raise ValueError("population/mask shape mismatch") + values = g[:, mask].astype(float) + if not np.isfinite(values).all(): + raise FloatingPointError("nonfinite leaf populations") + rho = values.sum(axis=0) + if (rho <= 0).any(): + raise FloatingPointError("nonpositive leaf density") + ux, uy = CX @ values / rho, CY @ values / rho + x, y = centers(mask.shape, level) + ex, ey, ep = exact(x[mask], y[mask], t, length, viscosity) + area = np.full(rho.size, 0.25**level) + accum.append((rho, ux, uy, (rho - 1) / 3, ex, ey, ep, area)) + rho, ux, uy, p, ex, ey, ep, area = [np.concatenate(v) for v in zip(*accum)] + volume = area.sum() + if not np.isclose(volume, length**2, rtol=0, atol=1e-10): + raise ValueError("leaf partition does not cover domain once") + mass = float(area @ rho) + # Remove one global pressure gauge, never one gauge per level. + p = p - area @ p / volume + ep = ep - area @ ep / volume + energy = float(area @ (ux * ux + uy * uy) / 2 / volume) + exact_energy = U * U / 4 * np.exp(-4 * viscosity * (2 * np.pi / length) ** 2 * t) + return dict( + tstar=t * U / length, + mass=mass, + mass_relative=( + 0.0 if initial_mass is None else abs(mass - initial_mass) / initial_mass + ), + velocity_l2=float( + np.sqrt( + area @ ((ux - ex) ** 2 + (uy - ey) ** 2) / (area @ (ex * ex + ey * ey)) + ) + ), + pressure_l2=float(np.sqrt(area @ ((p - ep) ** 2) / (area @ (ep * ep)))), + energy=energy, + exact_energy=exact_energy, + energy_relative=abs(energy - exact_energy) / exact_energy, + rho_min=float(rho.min()), + rho_max=float(rho.max()), + max_speed=float(np.max(np.hypot(ux, uy))), + leaf_cells=int(len(rho)), + ) diff --git a/tests/validation/taylor_green/run_case.py b/tests/validation/taylor_green/run_case.py new file mode 100644 index 0000000..9f0acb2 --- /dev/null +++ b/tests/validation/taylor_green/run_case.py @@ -0,0 +1,205 @@ +"""One bounded GPU case. No production edits, plots, checkpoints or field dumps.""" + +import argparse +import hashlib +import json +import os +from pathlib import Path +import subprocess +import sys +import time + +import numpy as np +from oracle import U, TIMES, LIMITS, initial, metrics, rectangles + + +def provenance(repo): + git = lambda *a: subprocess.check_output( + ["git", "-C", str(repo), *a], text=True + ).strip() + source_hash = hashlib.sha256() + for path in sorted((repo / "src").rglob("*")): + if path.is_file() and "__pycache__" not in path.parts and path.suffix != ".pyc": + source_hash.update( + str(path.relative_to(repo)).encode() + b"\0" + path.read_bytes() + b"\0" + ) + return dict( + revision=git("rev-parse", "HEAD"), + dirty=bool(git("status", "--porcelain")), + source_tree_sha256=source_hash.hexdigest(), + harness_sha256={ + p.name: hashlib.sha256(p.read_bytes()).hexdigest() + for p in Path(__file__).parent.glob("*.py") + }, + ) + + +def config(n, re): + from CelerisLab.config import LBMConfig + + cfg = LBMConfig( + nx=n, + ny=n, + viscosity=U * n / re, + velocity=0.0, + rho=1.0, + lattice_model="D2Q9", + collision="SRT", + streaming="double_buffer", + store_precision="FP32", + ddf_shifting=True, + les_enabled=False, + ) + cfg.validate() + return cfg + + +def workload(masks, base_steps): + cells = [int(np.count_nonzero(m)) for m in masks] + hops = [base_steps * 2**level for level in range(len(masks))] + updates = [count * hop for count, hop in zip(cells, hops)] + return dict( + active_cells_by_level=cells, + level_hops=hops, + active_collision_updates_by_level=updates, + active_collision_updates=sum(updates), + root_intervals=base_steps, + excludes="ghost streaming, halo copies, transfers, initialization and diagnostics", + ) + + +def main(): + p = argparse.ArgumentParser(description=__doc__) + p.add_argument("--repo", type=Path, required=True) + p.add_argument("--re", type=int, choices=[100, 1000], required=True) + p.add_argument("--n", type=int, choices=[32, 64, 128, 256], required=True) + p.add_argument("--depth", type=int, choices=[1, 2, 3], required=True) + p.add_argument("--output", type=Path, required=True) + args = p.parse_args() + if not os.environ.get("CUDA_VISIBLE_DEVICES"): + p.error("Explicit physical CUDA_VISIBLE_DEVICES selection required") + repo = args.repo.resolve() + sys.path.insert(0, str(repo / "src")) + import CelerisLab + + if Path(CelerisLab.__file__).resolve() != repo / "src/CelerisLab/__init__.py": + raise RuntimeError("Wrong imported source checkout") + from CelerisLab.cuda.context import CudaContext + from CelerisLab.lbm.blocks import BlockLayout2D + from CelerisLab.lbm.block_stepper import PeriodicBlockStepper + from CelerisLab.lbm.refinement_stepper import FixedPatchStepper + from CelerisLab.lbm.refinement import ACTIVE + from CelerisLab.lbm.hierarchy import build_hierarchy, leaf_partition + + cfg = config(args.n, args.re) + refinements = [] if args.depth == 1 else rectangles(args.n, args.depth) + h = build_hierarchy(cfg, refinements) + np.testing.assert_array_equal(leaf_partition(h), np.ones((args.n, args.n))) + masks = [k == ACTIVE for k in h.kinds] + report = dict( + schema=1, + source=provenance(repo), + n=args.n, + re_L=args.re, + depth=args.depth, + refinements=refinements, + U=U, + mach=U * np.sqrt(3), + viscosity=cfg.viscosity, + tau_by_level=[0.5 + 3 * c.viscosity for c in h.configs], + initialization="physical FP64 equilibrium G at actual cell centers; rho=1+3p", + cuda_visible_devices=os.environ["CUDA_VISIBLE_DEVICES"], + device_id=0, + limits=LIMITS, + workload=workload(masks, round(TIMES[-1] * args.n / U)), + rows=[], + completed=False, + ) + args.output.parent.mkdir(parents=True, exist_ok=True) + save = lambda: args.output.write_text( + json.dumps(report, indent=2, allow_nan=False) + "\n" + ) + save() + started = time.monotonic() + with CudaContext(0) as ctx: + if args.depth == 1: + layout = BlockLayout2D.build(args.n, args.n, 8, 8, periodic=(True, True)) + run = PeriodicBlockStepper(cfg, layout, context=ctx) + else: + run = FixedPatchStepper( + cfg, + context=ctx, + coarse_boundary="periodic", + refinements=refinements, + level_blocks=[(8, 8)] * args.depth, + ) + with run: + run.initialize( + *[ + initial(m.shape, level, args.n, cfg.viscosity) + for level, m in enumerate(masks) + ] + ) + initial_mass = None + for tstar in TIMES: + step = round(tstar * args.n / U) + run.step(step - run.step_count) + arrays = (run.download(),) if args.depth == 1 else run.download() + row = metrics(arrays, masks, step, args.n, cfg.viscosity, initial_mass) + row["step"] = step + if initial_mass is None: + initial_mass = row["mass"] + report["rows"].append(row) + save() + print( + json.dumps( + dict( + event="PROGRESS", + n=args.n, + re_L=args.re, + depth=args.depth, + **row, + ) + ), + flush=True, + ) + if step == 0: + # Exercise real kernels/halos/transfers on every level before READY. + # This is the first measured trajectory step, not discarded warmup. + run.step(1) + probe = (run.download(),) if args.depth == 1 else run.download() + metrics(probe, masks, 1, args.n, cfg.viscosity, initial_mass) + report["ready_after_step"] = 1 + save() + print( + json.dumps( + dict( + event="READY", + step=1, + n=args.n, + re_L=args.re, + depth=args.depth, + workload=report["workload"], + ) + ), + flush=True, + ) + report.update( + completed=True, + seconds=time.monotonic() - started, + screening_pass=all( + row[key] <= limit for row in report["rows"] for key, limit in LIMITS.items() + ), + ) + save() + print( + json.dumps( + {k: report[k] for k in ("n", "re_L", "depth", "screening_pass", "seconds")} + ), + flush=True, + ) + return 0 if report["screening_pass"] else 1 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/tests/validation/taylor_green/test_oracle.py b/tests/validation/taylor_green/test_oracle.py new file mode 100644 index 0000000..cd0231a --- /dev/null +++ b/tests/validation/taylor_green/test_oracle.py @@ -0,0 +1,170 @@ +"""Standalone CPU oracle/geometry checks. Set PYTHONPATH=/src.""" + +import unittest +import numpy as np +from oracle import CX, CY, U, exact, centers, initial, rectangles, metrics +from run_case import config + + +class OracleTests(unittest.TestCase): + def test_phase_pressure_and_decay(self): + u, v, p = exact(0.0, 8.0, 0.0, 32.0, 0.01) + self.assertAlmostEqual(u, -U) + self.assertAlmostEqual(v, 0.0) + self.assertAlmostEqual(p, 0.0) + self.assertAlmostEqual(exact(0.0, 0.0, 0.0, 32.0, 0.01)[2], -U * U / 2) + a = np.array(exact(3.0, 5.0, 0.0, 32.0, 0.01)) + b = np.array(exact(3.0, 5.0, 20.0, 32.0, 0.01)) + f = np.exp(-2 * 0.01 * (2 * np.pi / 32) ** 2 * 20) + np.testing.assert_allclose(b, a * np.array([f, f, f * f]), atol=1e-17) + + def test_navier_stokes_residual(self): + # Independent finite-difference derivatives catch sign/phase mistakes. + x = np.array([2.7, 8.1, 13.3]) + y = np.array([3.2, 7.4, 21.6]) + t = 10.0 + h = 0.002 + n = 32.0 + nu = 0.01 + f = lambda x, y, t: np.array(exact(x, y, t, n, nu)) + z = f(x, y, t) + dx = (f(x + h, y, t) - f(x - h, y, t)) / (2 * h) + dy = (f(x, y + h, t) - f(x, y - h, t)) / (2 * h) + dt = (f(x, y, t + h) - f(x, y, t - h)) / (2 * h) + lap = ( + f(x + h, y, t) + f(x - h, y, t) + f(x, y + h, t) + f(x, y - h, t) - 4 * z + ) / h**2 + self.assertLess(np.max(np.abs(dx[0] + dy[1])), 1e-11) + residual = ( + dt[:2] + + z[0] * dx[:2] + + z[1] * dy[:2] + + np.array([dx[2], dy[2]]) + - nu * lap[:2] + ) + self.assertLess(np.max(np.abs(residual)), 1e-10) + + def test_equilibrium_moments_and_metrics(self): + n = 32 + nu = 0.01 + g = initial((n, n), 0, n, nu) + rho = g.sum(axis=0) + x, y = centers((n, n)) + u, v, p = exact(x, y, 0, n, nu) + np.testing.assert_allclose(rho, 1 + 3 * p, atol=1e-15) + np.testing.assert_allclose(np.einsum("q,qyx->yx", CX, g) / rho, u, atol=1e-15) + np.testing.assert_allclose(np.einsum("q,qyx->yx", CY, g) / rho, v, atol=1e-15) + m = metrics((g,), [np.ones((n, n), bool)], 0, n, nu) + for k in ("velocity_l2", "pressure_l2", "energy_relative"): + self.assertLess(m[k], 1e-11) + with self.assertRaises(FloatingPointError): + g[:, 2, 2] = np.nan + metrics((g,), [np.ones((n, n), bool)], 0, n, nu) + + def test_geometry_cell_centers_and_scaling(self): + from CelerisLab.lbm.hierarchy import build_hierarchy, leaf_partition + from CelerisLab.lbm.refinement import ACTIVE + + for n in (32, 64): + for re in (100, 1000): + for depth in (2, 3): + cfg = config(n, re) + h = build_hierarchy(cfg, rectangles(n, depth)) + np.testing.assert_array_equal(leaf_partition(h), np.ones((n, n))) + for level, cfg_l in enumerate(h.configs): + self.assertAlmostEqual( + cfg_l.viscosity, cfg.viscosity * 2**level + ) + self.assertLessEqual(cfg_l.omega, cfg_l.omega_max) + x, y = centers(h.kinds[level].shape, level) + self.assertEqual(x[0, 0], 0.5 / 2**level) + self.assertEqual(x[0, -1], n - 0.5 / 2**level) + masks = [k == ACTIVE for k in h.kinds] + g = [ + initial(k.shape, l, n, cfg.viscosity) + for l, k in enumerate(h.kinds) + ] + result = metrics(g, masks, 0, n, cfg.viscosity) + self.assertLess(result["velocity_l2"], 1e-12) + self.assertLess(result["pressure_l2"], 1e-10) + self.assertLess(result["energy_relative"], 0.002) + + def test_complete_assessment_and_corruptions(self): + from copy import deepcopy + from assess import assess + from oracle import TIMES + + records = [] + for re in (100, 1000): + cases = [(n, 1) for n in (32, 64, 128, 256)] + cases += [(n, d) for n in (32, 64) for d in (2, 3)] + for n, d in cases: + rows = [ + dict( + tstar=t, + velocity_l2=0.0, + pressure_l2=0.0, + energy_relative=0.0, + mass_relative=0.0, + ) + for t in TIMES + ] + records.append( + dict( + re_L=re, + n=n, + depth=d, + U=U, + viscosity=U * n / re, + source=dict(revision="synthetic", dirty=False), + completed=True, + rows=rows, + ) + ) + self.assertTrue(assess(records)["passed"]) + bad = deepcopy(records) + bad[-1]["rows"][-1]["pressure_l2"] = 0.5 + failed = assess(bad) + self.assertFalse(failed["passed"]) + self.assertEqual(len(failed["comparisons"]), 24) + partial = assess(records[0:1] + records[4:5]) + self.assertFalse(partial["passed"]) + self.assertEqual(len(partial["comparisons"]), 3) + self.assertIsNone(partial["comparisons"][0]["uniform_finest"]) + + bad = deepcopy(records) + bad[-1]["rows"][-1]["pressure_l2"] = float("nan") + self.assertFalse(assess(bad)["passed"]) + bad = deepcopy(records) + bad[-1]["source"]["dirty"] = True + self.assertFalse(assess(bad)["passed"]) + bad = deepcopy(records) + bad[-1]["rows"][-1]["tstar"] = 0.49 + self.assertFalse(assess(bad)["passed"]) + with self.assertRaises(ValueError): + assess(records + [records[0]]) + + def test_workload(self): + from run_case import workload + from CelerisLab.lbm.hierarchy import build_hierarchy + from CelerisLab.lbm.refinement import ACTIVE + + cases = [(n, 1) for n in (32, 64, 128, 256)] + cases += [(n, d) for n in (32, 64) for d in (2, 3)] + costs = [] + for n, d in cases: + h = build_hierarchy(config(n, 100), [] if d == 1 else rectangles(n, d)) + costs.append(workload([k == ACTIVE for k in h.kinds], 16 * n)) + self.assertEqual(sum(c["active_collision_updates"] for c in costs), 384270336) + self.assertEqual(sum(c["root_intervals"] for c in costs), 10752) + self.assertEqual(costs[-1]["active_collision_updates"], 50987008) + + def test_incomplete_assessment_fails(self): + from assess import assess + + self.assertFalse(assess([])["passed"]) + self.assertEqual(len(assess([])["missing"]), 16) + + +if __name__ == "__main__": + unittest.main()