--- /dev/null +++ b/examples/fixed_patch_jet.py @@ -0,0 +1,179 @@ +"""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()