183 lines
9.7 KiB
Diff
183 lines
9.7 KiB
Diff
--- /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()
|