"""Low-Re startup flow in an x-periodic confined domain; nozzle remains coarse. From this checkout, after checking the selected physical GPU is idle: PYTHONPATH="$PWD/src" CUDA_VISIBLE_DEVICES=0 python examples/fixed_patch_jet.py --output /tmp/patch-jet The fixed-amplitude start includes acoustic transients, not a developed jet. Only native leaf cells are plotted. No interpolation, restart or flow oracle. """ import argparse from dataclasses import replace import json import os from pathlib import Path import tempfile import time import numpy as np import CelerisLab from CelerisLab.config import load_lbm_config from CelerisLab.cuda.context import CudaContext from CelerisLab.lbm.equilibrium import compute_feq_d2q9, compute_macro_from_ddf from CelerisLab.lbm.refinement import ACTIVE from CelerisLab.lbm.refinement_stepper import FixedPatchStepper BOUNDS, AMPLITUDE, WIDTH = (16, 4, 32, 28), .02, 12 TIMES, MAX_MACH, MAX_SECONDS = (0, 300, 1200), .1, 120 def measure(arrays, masks, centers, initial_mass): fields = [] for array, mask in zip(arrays, masks): if not np.isfinite(array[:, mask]).all(): raise RuntimeError("Nonfinite leaf populations") field = np.full((*mask.shape, 3), np.nan) field[mask] = [compute_macro_from_ddf(f) for f in array[:, mask].T] fields.append(field) leaves = np.concatenate([f[m] for f, m in zip(fields, masks)]) speed = np.hypot(leaves[:, 1], leaves[:, 2]) mass = sum(area*a[:, m].sum(dtype=np.float64) for area, a, m in zip((1, .25), arrays, masks)) record = dict(finite=bool(np.isfinite(leaves).all()), rho_min=float(leaves[:, 0].min()), rho_max=float(leaves[:, 0].max()), speed_max=float(speed.max()), measured_Ma=float(speed.max()*np.sqrt(3)), density_deviation=float(np.abs(leaves[:, 0]-1).max()), discrete_leaf_mass=float(mass), mass_change=float(mass-initial_mass)) if not record['finite'] or record['rho_min'] <= 0 or record['measured_Ma'] > MAX_MACH: raise RuntimeError(f"Finite/positive-density/Ma<={MAX_MACH} guard failed: {record}") walls = fields[0][[1, -2]] record['macro_sum_rho_v_bottom_top'] = (walls[:, :, 0]*walls[:, :, 2]).sum(axis=1).tolist() record['centerline'] = [] for level, y in ((0, 3.5), (1, 8.25), (1, 16.25), (0, 32.5)): xs, ys = centers[level] columns = np.argsort(np.abs(xs-24))[:2] row = int(np.flatnonzero(ys == y)[0]) values = fields[level][row, columns].mean(axis=0) record['centerline'].append(dict(level=level, x=24, y=y, sampled_x=sorted(xs[columns].tolist()), rule='arithmetic mean of two neighboring native columns', rho_ux_uy=values.tolist())) return fields, record def draw(out, snapshots, masks, centers): import matplotlib matplotlib.use('Agg') import matplotlib.pyplot as plt from matplotlib.collections import PolyCollection def decorate(ax): ax.set(xlim=(0, 48), ylim=(1.5, 46.5), aspect='equal', xlabel='X', ylabel='Y') ax.plot([16, 32, 32, 16, 16], [4, 4, 28, 28, 4], color='tab:red', lw=1) ax.axhline(1.5, color='black', lw=1) ax.axhline(46.5, color='black', lw=1, ls='--') ax.plot([18, 30], [1.5, 1.5], color='tab:orange', lw=4, clip_on=False) fig, ax = plt.subplots(figsize=(6, 6), constrained_layout=True) for mask, (xs, ys), dx in zip(masks, centers, (1, .5)): yy, xx = np.meshgrid(ys, xs, indexing='ij') cells = np.stack((xx[mask], yy[mask]), axis=-1) vertices = cells[:, None, :] + dx*np.array([[-.5, -.5], [.5, -.5], [.5, .5], [-.5, .5]]) ax.add_collection(PolyCollection(vertices, facecolors='none', edgecolors='0.55', linewidths=.3)) decorate(ax) ax.set_title('Native leaf grid: coarse dx=1, fine dx=0.5\nRed: internal patch; orange: coarse slot; dashed: reservoir') fig.savefig(out/'layout.png', dpi=150) plt.close(fig) vmax = max(1., *(np.nanmax(np.hypot(f[..., 1], f[..., 2]))/AMPLITUDE for fields in snapshots for f in fields)) fig, axes = plt.subplots(1, 3, figsize=(12, 4.8), constrained_layout=True) for ax, step, fields in zip(axes, TIMES, snapshots): for field, mask, (xs, ys), dx in zip(fields, masks, centers, (1, .5)): speed = np.ma.masked_where(~mask, np.hypot(field[..., 1], field[..., 2])/AMPLITUDE) image = ax.pcolormesh(np.r_[xs-dx/2, xs[-1]+dx/2], np.r_[ys-dx/2, ys[-1]+dx/2], speed, shading='flat', cmap='viridis', vmin=0, vmax=vmax, rasterized=True) stride = int(4/dx) yy, xx = np.meshgrid(ys, xs, indexing='ij') selection = np.zeros_like(mask) selection[::stride, ::stride] = mask[::stride, ::stride] ax.quiver(xx[selection], yy[selection], field[..., 1][selection], field[..., 2][selection], angles='xy', scale_units='xy', scale=AMPLITUDE/3, width=.005, color='white') decorate(ax) ax.set_title(f'coarse t={step}, fine steps={2*step}') fig.colorbar(image, ax=axes, label='speed / imposed amplitude') fig.suptitle('Low-Re startup flow in an x-periodic confined domain\nFixed-amplitude start includes acoustic transients; not steady/developed') fig.savefig(out/'velocity.png', dpi=150) plt.close(fig) def main(): parser = argparse.ArgumentParser(description=__doc__) parser.add_argument('--output', type=Path, required=True) out = parser.parse_args().output.resolve() root = Path(__file__).resolve().parents[1] if Path(CelerisLab.__file__).resolve() != root/'src/CelerisLab/__init__.py': raise RuntimeError(f'Set PYTHONPATH to {root / "src"}') out.mkdir(parents=True, exist_ok=False) cfg = replace(load_lbm_config(str(root/'src/CelerisLab/configs/config_lbm.json')), nx=48, ny=48, viscosity=.1, velocity=0., rho=1., wall_slot_width=WIDTH, wall_normal_velocity=AMPLITUDE, y_wall_bc='on_site_velocity', top_wall_bc='reservoir', lattice_model='D2Q9', collision='SRT', streaming='double_buffer', store_precision='FP32', ddf_shifting=False, les_enabled=False) report = dict(config=dict(nx=48, ny=48, bounds=BOUNDS, W=WIDTH, A=AMPLITUDE, nu=.1, rho=1, tau_coarse=.8, tau_fine=1.1, coarse_boundary='slot_reservoir'), guards=dict(check_every=50, positive_density=True, finite=True, max_Mach=MAX_MACH, max_seconds=MAX_SECONDS), notes='Node mass uses areas 1 and 1/4, not an exact physical volume integral. ' 'Boundary rho*v is a macro measurement, not a closed population budget.', samples=[], status='running', last_synchronized_coarse_step=0) start, run = time.monotonic(), None try: with CudaContext(0) as context, FixedPatchStepper(cfg, BOUNDS, context=context, coarse_boundary='slot_reservoir') as run: p = run.patch masks = [np.ones(p.coarse_shape, bool), p.fine_kind == ACTIVE] masks[0][[0, -1]] = False masks[0][4:28, 16:32] = False centers = [(origin[0]+dx*(np.arange(shape[1])+.5), origin[1]+dx*(np.arange(shape[0])+.5)) for shape, origin, dx in ((p.coarse_shape, (0, 0), 1), (p.fine_shape, p.fine_origin, .5))] initial = [np.broadcast_to(compute_feq_d2q9(1., 0., 0.)[:, None, None], (9, *m.shape)).copy() for m in masks] for x in range(cfg.nx): v = AMPLITUDE*max(0., 1-((2*x-(cfg.nx-1))/WIDTH)**2)**2 initial[0][:, 1, x] = compute_feq_d2q9(1., 0., v) run.initialize(*initial) mass0 = sum(area*a[:, m].sum(dtype=np.float64) for area, a, m in zip((1, .25), initial, masks)) snapshots = [] for step in range(0, TIMES[-1]+1, 50): if time.monotonic()-start > MAX_SECONDS: raise TimeoutError('120 second demonstration limit reached') if step: run.step(50) report['last_synchronized_coarse_step'] = run.step_count fields, sample = measure(run.download(), masks, centers, mass0) print(f'step={step} rho=[{sample["rho_min"]:.7g},{sample["rho_max"]:.7g}] ' f'Ma={sample["measured_Ma"]:.7g}', flush=True) if step in TIMES: sample.update(coarse_steps=step, fine_steps=2*step) report['samples'].append(sample) snapshots.append(fields) with tempfile.TemporaryDirectory(prefix='.plot-', dir=out) as cache: os.environ['MPLCONFIGDIR'] = cache draw(out, snapshots, masks, centers) report['status'] = 'complete' except BaseException as exc: report.update(status='failed', error=f'{type(exc).__name__}: {exc}', last_synchronized_coarse_step=run.step_count if run is not None else 0) raise finally: report['elapsed_seconds'] = time.monotonic()-start (out/'summary.json').write_text(json.dumps(report, indent=2)+'\n') print(f'{report["status"]}: last synchronized coarse step ' f'{report["last_synchronized_coarse_step"]}; output={out}', flush=True) if __name__ == '__main__': main()