Add streakline post-processing and exp control-matrix batch runners (vorticity + streakline)

This commit is contained in:
Frank14f 2026-05-27 19:19:51 +08:00
parent 50b2b6a7ca
commit 4758eb3215
30 changed files with 4128 additions and 4600 deletions

View File

@ -217,6 +217,31 @@ Tested on Tesla V100-SXM2-16GB (CUDA 12.4):
| Re100 EsoPull SRT | 384×192 | ~3900 |
| Re3000 MRT+LES | 384×192 | ~4360 |
### Performance methodology
For a "kernel-dominant" baseline (closest to FluidX3D-style throughput testing),
use the dedicated script:
```bash
conda run -n pycuda_3_10 python tests/run_perf_baseline.py \
--lattice-model D2Q9 --nx 384 --ny 192 --nz 1 \
--steps 4000 --warmup-steps 400 --batch-steps 100
```
This path times GPU stepping (`stepper.step`) and reports MLUPS and batch latency
percentiles. By default it avoids host readbacks inside the timed loop.
APIs that trigger device-to-host transfers (DTOH) and can reduce MLUPS:
- `Simulation.get_macroscopic()` / `LBMField.get_macroscopic()` (downloads full DDF)
- `Simulation.get_ddf()` / `LBMField.download_ddf()`
- `Simulation.save_checkpoint()` (downloads field/state buffers)
- Body observation downloads (e.g. `ObjectManager.download_obs_full_async(...)`)
Use `tests/run_perf_baseline.py` switches (`--macro-every`, `--ddf-every`,
`--checkpoint-every`, `--obs-every`) to quantify each overhead path against the
pure-step baseline.
## Citation
If you use CelerisLab in your research, please cite:

View File

@ -4,5 +4,6 @@ Common utilities and preprocessing functions.
"""
from . import preprocess
from . import streakline
__all__ = ['preprocess']
__all__ = ["preprocess", "streakline"]

View File

@ -0,0 +1,778 @@
# CelerisLab/common/streakline.py
"""Streakline utilities for online/offline particle-based flow visualization.
Responsibilities:
- Build dense release seeds (point/line/strip).
- Integrate particles with RK4 and space-time interpolated velocity fields.
- Support both online (in-memory frames from Simulation) and offline playback.
- Render density-like streakline images or accumulated path trails.
"""
from __future__ import annotations
import math
import os
from concurrent.futures import ThreadPoolExecutor
from dataclasses import dataclass, field
from typing import Any, List, Optional, Sequence, Tuple
import numpy as np
CylinderSpec = Tuple[Tuple[float, float], float]
@dataclass(frozen=True)
class FlowFrame:
"""One sampled velocity frame."""
step: int
ux: np.ndarray
uy: np.ndarray
@dataclass(frozen=True)
class ReleaseConfig:
"""Release strategy configuration for dense streakline seeding."""
mode: str = "strip" # point | line | strip
line_span: float = 12.0
line_count: int = 5
downstream_count: int = 4
downstream_spacing: float = 3.0
inject_per_seed: int = 2
@dataclass(frozen=True)
class IntegratorConfig:
"""Particle integration and filtering settings."""
alpha_t: float = 0.2
alpha_x: float = 0.4
age_decay_steps: float = 20000.0
max_particle_age: Optional[float] = None
diffusion_coeff: float = 0.0
random_seed: int = 1234
num_threads: int = 0
@dataclass
class ParticleTrailSet:
"""Growing polylines per particle (pathline history — not physical streakline)."""
trails: List[np.ndarray] = field(default_factory=list)
active: List[bool] = field(default_factory=list)
def inject(self, points: np.ndarray) -> None:
pts = np.asarray(points, dtype=np.float64)
if pts.ndim != 2 or pts.shape[1] != 2:
raise ValueError("points must be [N,2].")
for p in pts:
self.trails.append(np.asarray([[p[0], p[1]]], dtype=np.float64))
self.active.append(True)
def keep_indices(self, survived_mask: np.ndarray) -> None:
"""Deactivate trails whose particle did not survive the last advance."""
mask = np.asarray(survived_mask, dtype=bool).ravel()
j = 0
for i, is_active in enumerate(self.active):
if not is_active:
continue
if j >= mask.size:
break
if not bool(mask[j]):
self.active[i] = False
j += 1
def keep_mask(self, survived_mask: np.ndarray) -> None:
"""Alias for keep_indices (boolean mask over active particles)."""
self.keep_indices(survived_mask)
def append_positions(self, particles: np.ndarray) -> None:
pts = np.asarray(particles, dtype=np.float64)
j = 0
for i, is_active in enumerate(self.active):
if not is_active:
continue
if j >= pts.shape[0]:
break
self.trails[i] = np.vstack([self.trails[i], pts[j]])
j += 1
def configure_compute_threads(num_threads: int = 0) -> int:
"""Set BLAS/OpenMP thread env vars for host-side numpy work."""
if num_threads <= 0:
num_threads = min(32, os.cpu_count() or 4)
for var in (
"OMP_NUM_THREADS",
"OPENBLAS_NUM_THREADS",
"MKL_NUM_THREADS",
"NUMEXPR_NUM_THREADS",
"VECLIB_MAXIMUM_THREADS",
):
os.environ[var] = str(int(num_threads))
return int(num_threads)
def cylinders_from_triangle_layout(layout: dict) -> List[CylinderSpec]:
"""Build cylinder list from triangle-layout dict used in experiment notebook."""
radius = float(layout["radius_lb"])
return [
((float(layout["x_apex"]), float(layout["y_center"])), radius),
((float(layout["x_rear"]), float(layout["y_lower"])), radius),
((float(layout["x_rear"]), float(layout["y_upper"])), radius),
]
def estimate_sampling_plan(
*,
st_ref: float,
diameter: float,
u_ref: float,
snapshots_per_period: float = 24.0,
periods: int = 5,
) -> dict:
"""Estimate save interval and snapshot count from nominal shedding scale."""
period_steps = float(diameter) / (float(st_ref) * float(u_ref))
save_every = int(max(20, round(period_steps / snapshots_per_period / 10.0) * 10))
n_snapshots = int(max(20, round(float(periods) * float(snapshots_per_period))))
return {
"st_ref": float(st_ref),
"period_steps_est": float(period_steps),
"snapshots_per_period_target": float(snapshots_per_period),
"save_every_recommended": int(save_every),
"snapshot_count_recommended": int(n_snapshots),
}
def build_release_points(
base_points: np.ndarray,
release_cfg: ReleaseConfig,
) -> np.ndarray:
"""Expand base release points to dense points by line/strip policy."""
pts = np.asarray(base_points, dtype=np.float64)
if pts.ndim != 2 or pts.shape[1] != 2:
raise ValueError("base_points must be [N,2].")
mode = str(release_cfg.mode).lower()
if mode not in ("point", "line", "strip"):
raise ValueError("release mode must be point|line|strip.")
line_count = max(1, int(release_cfg.line_count))
ds_count = max(1, int(release_cfg.downstream_count))
line_offsets = np.linspace(
-0.5 * float(release_cfg.line_span),
0.5 * float(release_cfg.line_span),
line_count,
)
ds_offsets = np.arange(ds_count, dtype=np.float64) * float(release_cfg.downstream_spacing)
out: list[list[float]] = []
for x0, y0 in pts:
line_iter = [0.0] if mode == "point" else line_offsets
ds_iter = [0.0] if mode in ("point", "line") else ds_offsets
for dy in line_iter:
for dx in ds_iter:
out.append([float(x0 + dx), float(y0 + dy)])
return np.asarray(out, dtype=np.float64)
def _bilinear_sample(field: np.ndarray, x: np.ndarray, y: np.ndarray) -> np.ndarray:
ny, nx = field.shape
x0 = np.floor(x).astype(np.int64)
y0 = np.floor(y).astype(np.int64)
x1 = x0 + 1
y1 = y0 + 1
inside = (x0 >= 0) & (x1 < nx) & (y0 >= 0) & (y1 < ny)
out = np.full_like(x, np.nan, dtype=np.float64)
if not np.any(inside):
return out
xi = x[inside] - x0[inside]
yi = y[inside] - y0[inside]
f00 = field[y0[inside], x0[inside]]
f10 = field[y0[inside], x1[inside]]
f01 = field[y1[inside], x0[inside]]
f11 = field[y1[inside], x1[inside]]
out[inside] = (
(1.0 - xi) * (1.0 - yi) * f00
+ xi * (1.0 - yi) * f10
+ (1.0 - xi) * yi * f01
+ xi * yi * f11
)
return out
def _velocity_at(pos: np.ndarray, theta: float, f0: FlowFrame, f1: FlowFrame) -> np.ndarray:
x = pos[:, 0]
y = pos[:, 1]
vx0 = _bilinear_sample(f0.ux, x, y)
vy0 = _bilinear_sample(f0.uy, x, y)
vx1 = _bilinear_sample(f1.ux, x, y)
vy1 = _bilinear_sample(f1.uy, x, y)
vx = (1.0 - theta) * vx0 + theta * vx1
vy = (1.0 - theta) * vy0 + theta * vy1
return np.stack([vx, vy], axis=1)
def _particles_valid_mask(
p_new: np.ndarray,
*,
nx: int,
ny: int,
solid_center: Optional[Tuple[float, float]] = None,
solid_radius: Optional[float] = None,
cylinders: Optional[Sequence[CylinderSpec]] = None,
) -> np.ndarray:
valid = np.all(np.isfinite(p_new), axis=1)
valid &= (p_new[:, 0] >= 0.0) & (p_new[:, 0] <= nx - 1)
valid &= (p_new[:, 1] >= 0.0) & (p_new[:, 1] <= ny - 1)
specs: List[CylinderSpec] = []
if cylinders:
specs.extend(cylinders)
elif solid_center is not None and solid_radius is not None:
specs.append((solid_center, float(solid_radius)))
for (cx, cy), radius in specs:
dx = p_new[:, 0] - float(cx)
dy = p_new[:, 1] - float(cy)
valid &= (dx * dx + dy * dy) > (float(radius) * float(radius))
return valid
def _inject_particles(
particles: np.ndarray,
ages: np.ndarray,
release_points: np.ndarray,
*,
inject_per_seed: int,
) -> Tuple[np.ndarray, np.ndarray]:
new_pts = np.repeat(release_points, max(1, int(inject_per_seed)), axis=0)
if particles.size == 0:
return new_pts.copy(), np.zeros((new_pts.shape[0],), dtype=np.float64)
return (
np.vstack([particles, new_pts]),
np.concatenate([ages, np.zeros((new_pts.shape[0],), dtype=np.float64)]),
)
def advance_particles_between_frames(
particles: np.ndarray,
ages: np.ndarray,
frame0: FlowFrame,
frame1: FlowFrame,
*,
nx: int,
ny: int,
cfg: IntegratorConfig,
solid_center: Optional[Tuple[float, float]] = None,
solid_radius: Optional[float] = None,
cylinders: Optional[Sequence[CylinderSpec]] = None,
rng: Optional[np.random.Generator] = None,
) -> Tuple[np.ndarray, np.ndarray, dict, np.ndarray]:
"""Advance particles from frame0.step to frame1.step with RK4.
Returns:
particles, ages, diagnostics, survived_mask (bool, length = input count).
"""
if particles.size == 0:
return particles, ages, {"n_substeps": 0, "max_speed": 0.0}, np.zeros((0,), dtype=bool)
dt_save = float(frame1.step - frame0.step)
max_speed = float(
max(
np.nanmax(np.sqrt(frame0.ux * frame0.ux + frame0.uy * frame0.uy)),
np.nanmax(np.sqrt(frame1.ux * frame1.ux + frame1.uy * frame1.uy)),
)
)
umax = max(max_speed, 1e-8)
dt_trace = min(float(cfg.alpha_t) * dt_save, float(cfg.alpha_x) / umax)
n_sub = max(1, int(math.ceil(dt_save / dt_trace)))
dt_sub = dt_save / float(n_sub)
p = np.asarray(particles, dtype=np.float64)
a = np.asarray(ages, dtype=np.float64)
track = np.ones((p.shape[0],), dtype=bool)
local_rng = rng or np.random.default_rng(int(cfg.random_seed))
diffusion = float(cfg.diffusion_coeff)
for sub_idx in range(n_sub):
if p.size == 0:
break
theta = float(sub_idx) / float(n_sub)
k1 = _velocity_at(p, theta, frame0, frame1)
k2 = _velocity_at(p + 0.5 * dt_sub * k1, min(1.0, theta + 0.5 / n_sub), frame0, frame1)
k3 = _velocity_at(p + 0.5 * dt_sub * k2, min(1.0, theta + 0.5 / n_sub), frame0, frame1)
k4 = _velocity_at(p + dt_sub * k3, min(1.0, theta + 1.0 / n_sub), frame0, frame1)
p_new = p + (k1 + 2.0 * k2 + 2.0 * k3 + k4) * (dt_sub / 6.0)
if diffusion > 0.0:
sigma = math.sqrt(max(0.0, 2.0 * diffusion * dt_sub))
p_new += local_rng.normal(0.0, sigma, size=p_new.shape)
valid = _particles_valid_mask(
p_new,
nx=nx,
ny=ny,
solid_center=solid_center,
solid_radius=solid_radius,
cylinders=cylinders,
)
track[track] = valid
p = p_new[valid]
a = a[valid] + dt_sub
if cfg.max_particle_age is not None:
age_keep = a <= float(cfg.max_particle_age)
track[track] = age_keep
p = p[age_keep]
a = a[age_keep]
survived_mask = track
return p, a, {"dt_trace": float(dt_sub), "n_substeps": int(n_sub), "max_speed": float(max_speed)}, survived_mask
def run_streakline_offline(
frames: Sequence[FlowFrame],
*,
nx: int,
ny: int,
release_points: np.ndarray,
release_cfg: ReleaseConfig,
integrator_cfg: IntegratorConfig,
solid_center: Optional[Tuple[float, float]] = None,
solid_radius: Optional[float] = None,
cylinders: Optional[Sequence[CylinderSpec]] = None,
) -> Tuple[np.ndarray, np.ndarray, dict]:
"""Reconstruct streakline from pre-collected velocity frames."""
if len(frames) < 2:
raise ValueError("At least two frames are required.")
pts = build_release_points(release_points, release_cfg)
particles = np.zeros((0, 2), dtype=np.float64)
ages = np.zeros((0,), dtype=np.float64)
rng = np.random.default_rng(int(integrator_cfg.random_seed))
max_speed = 0.0
last_step = int(frames[0].step)
total_substeps = 0
for i in range(len(frames) - 1):
particles, ages = _inject_particles(
particles,
ages,
pts,
inject_per_seed=release_cfg.inject_per_seed,
)
particles, ages, diag, _ = advance_particles_between_frames(
particles,
ages,
frames[i],
frames[i + 1],
nx=nx,
ny=ny,
cfg=integrator_cfg,
solid_center=solid_center,
solid_radius=solid_radius,
cylinders=cylinders,
rng=rng,
)
max_speed = max(max_speed, float(diag["max_speed"]))
total_substeps += int(diag["n_substeps"])
last_step = int(frames[i + 1].step)
return particles, ages, {
"frames_used": int(len(frames)),
"last_step": int(last_step),
"n_particles_alive": int(particles.shape[0]),
"max_speed_seen": float(max_speed),
"total_substeps": int(total_substeps),
}
def run_streakline_online(
sim: Any,
*,
start_step: int,
sample_every: int,
n_samples: int,
release_points: np.ndarray,
release_cfg: ReleaseConfig,
integrator_cfg: IntegratorConfig,
solid_center: Optional[Tuple[float, float]] = None,
solid_radius: Optional[float] = None,
cylinders: Optional[Sequence[CylinderSpec]] = None,
) -> Tuple[np.ndarray, np.ndarray, dict]:
"""Compute streakline online from a running Simulation without disk snapshots."""
if sample_every < 1 or n_samples < 2:
raise ValueError("sample_every must be >=1 and n_samples must be >=2.")
if not getattr(sim, "_initialized", False):
raise RuntimeError("Simulation must be initialized before online streakline run.")
nx = int(sim.lbm_cfg.nx)
ny = int(sim.lbm_cfg.ny)
frames: list[FlowFrame] = []
target_last = int(start_step + sample_every * (n_samples - 1))
while int(sim.stepper.step_count) < target_last:
sim.step(1)
step = int(sim.stepper.step_count)
if step < int(start_step):
continue
if (step - int(start_step)) % int(sample_every) != 0:
continue
macro = sim.get_macroscopic()
frames.append(
FlowFrame(
step=step,
ux=np.asarray(macro["ux"], dtype=np.float64),
uy=np.asarray(macro["uy"], dtype=np.float64),
)
)
if len(frames) >= int(n_samples):
break
if len(frames) < 2:
raise RuntimeError("Online run did not collect enough frames.")
particles, ages, diag = run_streakline_offline(
frames,
nx=nx,
ny=ny,
release_points=release_points,
release_cfg=release_cfg,
integrator_cfg=integrator_cfg,
solid_center=solid_center,
solid_radius=solid_radius,
cylinders=cylinders,
)
diag["mode"] = "online"
diag["sample_every"] = int(sample_every)
diag["start_step"] = int(start_step)
return particles, ages, diag
def compute_vorticity(ux: np.ndarray, uy: np.ndarray) -> np.ndarray:
"""Central-difference vorticity ω = ∂uy/∂x ∂ux/∂y on interior cells."""
ux_arr = np.asarray(ux, dtype=np.float64)
uy_arr = np.asarray(uy, dtype=np.float64)
vort = np.zeros_like(ux_arr, dtype=np.float64)
dudy = 0.5 * (ux_arr[2:, 1:-1] - ux_arr[:-2, 1:-1])
dvdx = 0.5 * (uy_arr[1:-1, 2:] - uy_arr[1:-1, :-2])
vort[1:-1, 1:-1] = dvdx - dudy
return vort
def render_vorticity_field(
vort: np.ndarray,
*,
nx: int,
ny: int,
out_path: str,
cylinders: Optional[Sequence[CylinderSpec]] = None,
vmin: Optional[float] = None,
vmax: Optional[float] = None,
percentile: float = 99.0,
minimal_axes: bool = True,
) -> dict:
"""Render vorticity with white at zero, green negative, purple positive."""
try:
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from matplotlib.colors import LinearSegmentedColormap, TwoSlopeNorm
except ImportError as exc:
raise RuntimeError("render_vorticity_field requires matplotlib.") from exc
field = np.asarray(vort, dtype=np.float64)
vmax_field = float(np.nanmax(np.abs(field)))
abs_pct = float(np.nanpercentile(np.abs(field), percentile))
# Freestream dominates |ω|; do not let a near-zero percentile wash out wake structure.
abs_max = max(abs_pct, 0.35 * vmax_field, 1e-6)
if vmin is None:
vmin = -abs_max
if vmax is None:
vmax = abs_max
if vmin >= 0.0 or vmax <= 0.0:
raise ValueError(f"vmin/vmax must bracket zero: vmin={vmin}, vmax={vmax}")
cmap = LinearSegmentedColormap.from_list(
"vort_green_white_purple",
["#1a9850", "#f7f7f7", "#762a83"],
N=256,
)
norm = TwoSlopeNorm(vmin=float(vmin), vcenter=0.0, vmax=float(vmax))
rgba = cmap(norm(field))
if cylinders:
yy, xx = np.mgrid[0:ny, 0:nx]
for (cx, cy), radius in cylinders:
mask = (xx - float(cx)) ** 2 + (yy - float(cy)) ** 2 <= float(radius) ** 2
rgba[mask] = (0.0, 0.0, 0.0, 1.0)
os.makedirs(os.path.dirname(os.path.abspath(out_path)), exist_ok=True)
fig = plt.figure(figsize=(nx / 100.0, ny / 100.0), dpi=100)
ax = fig.add_axes([0.0, 0.0, 1.0, 1.0])
ax.imshow(rgba, origin="lower", interpolation="nearest", aspect="auto")
if minimal_axes:
ax.axis("off")
fig.savefig(out_path, dpi=100, bbox_inches="tight", pad_inches=0)
plt.close(fig)
return {
"image_path": out_path,
"vmin": float(vmin),
"vmax": float(vmax),
"vort_min": float(np.nanmin(field)),
"vort_max": float(np.nanmax(field)),
}
def gaussian_blur2d(img: np.ndarray, sigma: float = 1.0) -> np.ndarray:
"""Separable Gaussian blur helper."""
if float(sigma) <= 0.0:
return np.asarray(img, dtype=np.float64)
radius = max(1, int(math.ceil(3.0 * float(sigma))))
x = np.arange(-radius, radius + 1, dtype=np.float64)
k = np.exp(-(x**2) / (2.0 * sigma * sigma))
k /= np.sum(k)
arr = np.asarray(img, dtype=np.float64)
arr = np.apply_along_axis(lambda m: np.convolve(m, k, mode="same"), 1, arr)
arr = np.apply_along_axis(lambda m: np.convolve(m, k, mode="same"), 0, arr)
return arr
def _rasterize_segment(
intensity: np.ndarray,
x0: float,
y0: float,
x1: float,
y1: float,
weight: float,
) -> None:
"""Stamp one segment into a scalar intensity map (for white-bg red compose)."""
ny, nx = intensity.shape
n = int(max(abs(x1 - x0), abs(y1 - y0), 1.0)) + 1
xs = np.linspace(x0, x1, n)
ys = np.linspace(y0, y1, n)
xi = np.clip(np.rint(xs).astype(np.int64), 0, nx - 1)
yi = np.clip(np.rint(ys).astype(np.int64), 0, ny - 1)
w = float(max(0.0, min(1.0, weight)))
np.add.at(intensity, (yi, xi), w)
def _draw_trail(
intensity: np.ndarray,
trail: np.ndarray,
*,
fade_along_trail: bool,
) -> None:
if trail.shape[0] < 2:
return
nseg = trail.shape[0] - 1
for i in range(nseg):
if fade_along_trail:
t = float(i + 1) / float(nseg)
weight = 0.15 + 0.85 * t
else:
weight = 1.0
_rasterize_segment(
intensity,
float(trail[i, 0]),
float(trail[i, 1]),
float(trail[i + 1, 0]),
float(trail[i + 1, 1]),
weight,
)
def _draw_cylinders(rgb: np.ndarray, cylinders: Sequence[CylinderSpec]) -> None:
ny, nx = rgb.shape[0], rgb.shape[1]
yy, xx = np.mgrid[0:ny, 0:nx]
for (cx, cy), radius in cylinders:
mask = (xx - float(cx)) ** 2 + (yy - float(cy)) ** 2 <= float(radius) ** 2
rgb[mask] = 0.0
def render_snapshot_trails(
trail_set: ParticleTrailSet,
*,
nx: int,
ny: int,
out_path: str,
cylinders: Optional[Sequence[CylinderSpec]] = None,
blur_sigma: float = 0.0,
fade_along_trail: bool = True,
num_threads: int = 0,
background_color: Tuple[float, float, float] = (1.0, 1.0, 1.0),
streak_color: Tuple[float, float, float] = (1.0, 0.0, 0.0),
) -> dict:
"""Render accumulated path trails (pathline collage — debug only, not streakline)."""
n_threads = configure_compute_threads(num_threads)
intensity = np.zeros((ny, nx), dtype=np.float64)
trails = trail_set.trails
if not trails:
density_max = 0.0
else:
chunk = max(1, (len(trails) + n_threads - 1) // n_threads)
def _worker(batch: List[np.ndarray]) -> np.ndarray:
local = np.zeros((ny, nx), dtype=np.float64)
for tr in batch:
_draw_trail(local, tr, fade_along_trail=fade_along_trail)
return local
batches = [trails[i : i + chunk] for i in range(0, len(trails), chunk)]
if len(batches) == 1:
intensity = _worker(batches[0])
else:
with ThreadPoolExecutor(max_workers=n_threads) as pool:
parts = list(pool.map(_worker, batches))
for part in parts:
intensity = np.maximum(intensity, part)
if np.any(intensity > 0):
vmax = float(np.percentile(intensity[intensity > 0], 99.0))
intensity = np.clip(intensity / max(vmax, 1e-12), 0.0, 1.0)
density_max = float(np.max(intensity))
rgb = np.ones((ny, nx, 3), dtype=np.float64)
rgb[:, :, 0] = background_color[0]
rgb[:, :, 1] = background_color[1] * (1.0 - intensity) + streak_color[1] * intensity
rgb[:, :, 2] = background_color[2] * (1.0 - intensity) + streak_color[2] * intensity
if streak_color[0] < 0.999:
rgb[:, :, 0] = background_color[0] * (1.0 - intensity) + streak_color[0] * intensity
if cylinders:
_draw_cylinders(rgb, cylinders)
if float(blur_sigma) > 0.0:
for c in range(3):
rgb[:, :, c] = gaussian_blur2d(rgb[:, :, c], sigma=float(blur_sigma))
rgb = np.clip(rgb, 0.0, 1.0)
os.makedirs(os.path.dirname(os.path.abspath(out_path)), exist_ok=True)
try:
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
fig = plt.figure(figsize=(nx / 100.0, ny / 100.0), dpi=100)
ax = fig.add_axes([0.0, 0.0, 1.0, 1.0])
ax.axis("off")
ax.imshow(rgb, origin="lower", interpolation="nearest", aspect="auto")
fig.savefig(out_path, dpi=100, bbox_inches="tight", pad_inches=0)
plt.close(fig)
except ImportError:
from PIL import Image
Image.fromarray((np.clip(rgb, 0.0, 1.0) * 255.0).astype(np.uint8)).save(out_path)
return {
"image_path": out_path,
"trail_count": int(len(trails)),
"density_max": density_max,
"num_threads": int(n_threads),
}
def render_streakline_density(
positions: np.ndarray,
ages: np.ndarray,
*,
nx: int,
ny: int,
out_path: str,
release_points: Optional[np.ndarray] = None,
solid_center: Optional[Tuple[float, float]] = None,
solid_radius: Optional[float] = None,
cylinders: Optional[Sequence[CylinderSpec]] = None,
age_decay_steps: float = 20000.0,
blur_sigma: float = 1.2,
title: str = "Streakline",
minimal_axes: bool = False,
background_color: Tuple[float, float, float] = (0.0, 0.0, 0.0),
streak_color: Tuple[float, float, float] = (1.0, 0.5, 0.0),
show_release_points: bool = True,
) -> dict:
"""Render streakline as current particle positions (density cloud at observation time)."""
try:
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from matplotlib.colors import LinearSegmentedColormap
except ImportError as exc:
raise RuntimeError("render_streakline_density requires matplotlib.") from exc
if positions.size == 0:
density = np.zeros((ny, nx), dtype=np.float64)
else:
w = np.exp(-np.asarray(ages, dtype=np.float64) / max(1.0, float(age_decay_steps)))
y_idx = np.clip(np.floor(positions[:, 1]).astype(np.int64), 0, ny - 1)
x_idx = np.clip(np.floor(positions[:, 0]).astype(np.int64), 0, nx - 1)
density = np.zeros((ny, nx), dtype=np.float64)
np.add.at(density, (y_idx, x_idx), w)
density = gaussian_blur2d(density, sigma=float(blur_sigma))
if minimal_axes:
rgb = np.ones((ny, nx, 3), dtype=np.float64)
rgb[:, :, 0] = background_color[0]
rgb[:, :, 1] = background_color[1]
rgb[:, :, 2] = background_color[2]
if np.any(density > 0):
vmax = float(np.percentile(density[density > 0], 99.0))
alpha = np.clip(density / max(vmax, 1e-12), 0.0, 1.0)
for c, val in enumerate(streak_color):
rgb[:, :, c] = np.clip(rgb[:, :, c] * (1.0 - alpha) + val * alpha, 0.0, 1.0)
specs = list(cylinders or [])
if solid_center is not None and solid_radius is not None:
specs.append((solid_center, float(solid_radius)))
if specs:
_draw_cylinders(rgb, specs)
os.makedirs(os.path.dirname(os.path.abspath(out_path)), exist_ok=True)
fig = plt.figure(figsize=(nx / 100.0, ny / 100.0), dpi=100)
ax = fig.add_axes([0.0, 0.0, 1.0, 1.0])
ax.axis("off")
ax.imshow(rgb, origin="lower", interpolation="nearest", aspect="auto")
fig.savefig(out_path, dpi=100, bbox_inches="tight", pad_inches=0)
plt.close(fig)
return {
"image_path": out_path,
"n_particles": int(positions.shape[0]) if positions.ndim == 2 else 0,
"density_max": float(np.max(density)) if density.size else 0.0,
}
os.makedirs(os.path.dirname(os.path.abspath(out_path)), exist_ok=True)
fig, ax = plt.subplots(figsize=(12, 5))
vmax = np.percentile(density[density > 0], 99.0) if np.any(density > 0) else 1.0
cmap = LinearSegmentedColormap.from_list("streak", [(0.0, 0.0, 0.0), streak_color])
ax.imshow(
density,
origin="lower",
cmap=cmap,
extent=(0, nx - 1, 0, ny - 1),
vmin=0.0,
vmax=max(vmax, 1e-12),
aspect="equal",
)
if solid_center is not None and solid_radius is not None:
circ = plt.Circle(solid_center, radius=solid_radius, fill=False, color="cyan", linewidth=1.4)
ax.add_patch(circ)
if show_release_points and release_points is not None and np.asarray(release_points).size:
rp = np.asarray(release_points)
ax.scatter(rp[:, 0], rp[:, 1], c="white", s=10, marker="x")
ax.set_title(title)
ax.set_xlabel("x (lattice)")
ax.set_ylabel("y (lattice)")
ax.set_xlim(0, nx - 1)
ax.set_ylim(0, ny - 1)
fig.tight_layout()
fig.savefig(out_path, dpi=170, bbox_inches="tight")
plt.close(fig)
return {
"image_path": out_path,
"n_particles": int(positions.shape[0]) if positions.ndim == 2 else 0,
"density_max": float(np.max(density)) if density.size else 0.0,
}

View File

@ -2,8 +2,8 @@
"_doc": "Three rotating cylinders in confined channel (triangle layout, free-slip walls).",
"grid": {
"lattice_model": "D2Q9",
"nx": 6000,
"ny": 1200,
"nx": 3000,
"ny": 300,
"nz": 1
},
"physics": {
@ -27,7 +27,9 @@
},
"inlet": {
"profile": "uniform",
"trt_neq_damp": 0.5
"scheme": "regularized",
"trt_neq_damp": 0.5,
"regularized_neq_damp": 0.5
},
"outlet": {
"mode": "neq_extrap",

View File

@ -69,6 +69,12 @@ class LBMField:
# Snapshot
self._ddf_snap: np.ndarray | None = None
self._host_ddf_step: int | None = None
# Cached lattice constants mirrored from CUDA __constant__ memory.
self._cached_w: np.ndarray | None = None
self._cached_vectors_2d: tuple[np.ndarray, np.ndarray] | None = None
self._cached_vectors_3d: tuple[np.ndarray, np.ndarray, np.ndarray] | None = None
# Host mirror of runtime parameters so partial updates preserve state.
self._runtime_params = {
@ -82,6 +88,16 @@ class LBMField:
# Upload d_params immediately
self._upload_params()
def invalidate_host_ddf_cache(self) -> None:
"""Mark host DDF mirror as stale after GPU-side updates."""
self._host_ddf_step = None
def invalidate_module_constant_cache(self) -> None:
"""Drop cached __constant__ mirrors after module replacement."""
self._cached_w = None
self._cached_vectors_2d = None
self._cached_vectors_3d = None
def build_channel_flags(self) -> np.ndarray:
"""Return a fresh host flag array for the base channel domain.
@ -123,9 +139,15 @@ class LBMField:
return self.sensors.count
# -- DDF download with precision decode ----------------------------------
def download_ddf(self):
def download_ddf(self, *, step_id: int | None = None, force: bool = False):
"""Copy DDF from GPU to host self.ddf (float32), decoding FP16S if needed.
If DDF shifting is enabled, un-shifts by adding back lattice weights."""
# Fast path for repeated host queries within one solver step.
# Callers pass step_id from LBMStepper; force=True bypasses cache
# for snapshot-style semantics.
if not force and step_id is not None and self._host_ddf_step == int(step_id):
return
if self.store_bytes == 2:
buf = np.empty(self.n * self.nq, dtype=self.store_dtype)
cuda.memcpy_dtoh(buf, self.ddf_gpu)
@ -143,13 +165,17 @@ class LBMField:
ddf_2d = self.ddf.reshape(self.nq, -1)
for i in range(self.nq):
ddf_2d[i] += w[i]
self._host_ddf_step = int(step_id) if step_id is not None else None
def _read_lattice_weights(self):
"""Read d_w[] from CUDA __constant__ memory (single source of truth in descriptors.cuh)."""
if self._cached_w is not None:
return self._cached_w
ptr, size = self.module.get_global("d_w")
w = np.empty(self.nq, dtype=np.float32)
cuda.memcpy_dtoh(w, ptr)
return w
self._cached_w = w
return self._cached_w
# -- d_params upload -----------------------------------------------------
def _upload_params(self):
@ -209,6 +235,7 @@ class LBMField:
else:
cuda.memcpy_htod(self.ddf_gpu, self.ddf)
cuda.memcpy_htod(self.temp_gpu, self.ddf)
self.invalidate_host_ddf_cache()
# The canonical download_ddf is defined above in the class body.
@ -223,6 +250,11 @@ class LBMField:
# -- Read lattice descriptors from CUDA module ---------------------------
def _read_lattice_vectors(self):
"""Read d_cx, d_cy [, d_cz] from CUDA __constant__ memory."""
if self.cfg.is_d2q9 and self._cached_vectors_2d is not None:
return self._cached_vectors_2d
if self.cfg.is_d3q19 and self._cached_vectors_3d is not None:
return self._cached_vectors_3d
nq = self.nq
cx = np.empty(nq, dtype=np.int32)
cy = np.empty(nq, dtype=np.int32)
@ -234,13 +266,17 @@ class LBMField:
cz = np.empty(nq, dtype=np.int32)
ptr_cz, _ = self.module.get_global("d_cz")
cuda.memcpy_dtoh(cz, ptr_cz)
return cx, cy, cz
return cx, cy
self._cached_vectors_3d = (cx, cy, cz)
return self._cached_vectors_3d
self._cached_vectors_2d = (cx, cy)
return self._cached_vectors_2d
# -- Macroscopic field extraction ----------------------------------------
def get_macroscopic(self):
def get_macroscopic(self, *, step_id: int | None = None):
"""Download DDF and compute rho, ux, uy [, uz] on host."""
self.download_ddf()
# Reuse the DDF host mirror when step_id is unchanged to avoid
# duplicate DTOH in repeated diagnostics of the same timestep.
self.download_ddf(step_id=step_id)
nq = self.nq
if self.cfg.is_d2q9:
@ -266,7 +302,7 @@ class LBMField:
# -- Snapshots -----------------------------------------------------------
def snapshot(self):
self.download_ddf()
self.download_ddf(force=True)
self._ddf_snap = self.ddf.copy()
def restore(self):

View File

@ -6,8 +6,8 @@
#define NT 256
#define MULT_GPU 0
#define NX 401
#define NY 201
#define NX 3000
#define NY 300
#define NZ 1
// ---- Lattice model (single source of truth) ----

View File

@ -17,7 +17,7 @@
#define OUTLET_MODE 0
#define OUTLET_BLEND_ALPHA 0.700f
#define OUTLET_BACKFLOW_CLAMP 1
#define Y_WALL_BC 0
#define Y_WALL_BC 1
#define OMEGA_COLLISION_MIN 0.01f
#define OMEGA_COLLISION_MAX 1.960f

View File

@ -3,6 +3,6 @@
#ifndef CELERIS_CONFIG_OBJECTS_H
#define CELERIS_CONFIG_OBJECTS_H
#define N_OBJS 0
#define N_OBJS 3
#endif

View File

@ -4,9 +4,9 @@
#define CELERIS_CONFIG_PHYSICS_H
#define LBtype float
#define VIS 0.0090000000
#define VIS 0.0040000000
#define RHO 1.0
#define U0 0.03
#define U0 0.04
#define PI 3.141592653589793238

View File

@ -51,7 +51,7 @@ class LBMStepper:
_gpu_bytes = f.n * f.nq * f.store_bytes
cuda.memcpy_dtod(f.temp_gpu, f.ddf_gpu, _gpu_bytes)
cuda.memcpy_dtoh(f.flag, f.flag_gpu)
f.download_ddf()
f.download_ddf(step_id=self._step_count, force=True)
def step(
self,
@ -88,6 +88,9 @@ class LBMStepper:
self._launch_sensor(obs_gpu, **launch_kw)
self._step_count += 1
if n > 0:
# Any completed GPU step invalidates host DDF mirror.
f.invalidate_host_ddf_cache()
def _launch_curved(self, action_gpu, obs_gpu, **launch_kw):
f = self.field

View File

@ -125,6 +125,7 @@ class Simulation:
self._module = compiler.load_module(self._ptx_path)
# Reconnect field and stepper to new module
self.field.module = self._module
self.field.invalidate_module_constant_cache()
self.field._upload_params()
_prev_step_count = self.stepper._step_count
self.stepper = LBMStepper(
@ -194,10 +195,10 @@ class Simulation:
# -- Data access ---------------------------------------------------------
def get_macroscopic(self) -> Dict[str, np.ndarray]:
"""Download DDF and return rho, ux, uy [, uz]."""
return self.field.get_macroscopic()
return self.field.get_macroscopic(step_id=self.stepper.step_count)
def get_ddf(self) -> np.ndarray:
self.field.download_ddf()
self.field.download_ddf(step_id=self.stepper.step_count)
return self.field.ddf.copy()
def get_flags(self) -> np.ndarray:

View File

@ -1,64 +1,52 @@
# Rotating cylinder validation plan
# Rotating cylinder validation against [Kan99b]
##### [**Undermind**](https://undermind.ai)
## Goal
---
This validation should stay small, direct, and defensible.
## Rotating cylinder validation against \[Kan99b\]
The main design rules are:
This plan defines a practical validation campaign for 2D flow past a rotating circular cylinder using the current CelerisLab solver. The reference is the rotating cylinder study of Kang, Choi, and Lee \[Kan99b\]. The goal is not to reproduce their O type far field mesh exactly. The goal is to show that the present rectangular domain with uniform inflow, slip top and bottom boundaries, curved moving wall treatment, and open outlet can recover the same force, shedding, and suppression trends with controlled boundary sensitivity.
- use the paper's direct numeric anchor at `Re = 100, alpha = 1.0` as the main hard benchmark
- use a low-rotation case to test the lift trend
- use suppression cases to test flow classification, not exact threshold fitting
- do not treat values read from figures near `alpha_L` as tight amplitude targets
The validation focuses on three questions:
This keeps the matrix representative without overfitting to sensitive threshold points.
- whether the moving curved wall treatment gives the correct mean force trend under rotation
- whether the wake frequency and fluctuation amplitudes are credible before shedding is suppressed
- whether the predicted critical rotation rate for shedding suppression is close to the literature value
## Strong numeric anchors from [Kan99b]
## Reference targets from \[Kan99b\]
The paper defines the Reynolds number and spin ratio as
``` math
Re = \frac{U_\infty D}{\nu}, \qquad \alpha = \frac{\omega_{body} D}{2 U_\infty}
```
and uses the standard force coefficients
``` math
C_D = \frac{2 F_x}{\rho U_\infty^2 D}, \qquad C_L = \frac{2 F_y}{\rho U_\infty^2 D}
```
The fluctuation amplitudes are half peak to peak values over one shedding period, not RMS values \[Kan99b\].
The most useful exact anchor point is the tabulated convergence case at $`Re=100`$ and $`\alpha=1.0`$ \[Kan99b\].
The strongest exact benchmark in the paper is the convergence case at `Re = 100, alpha = 1.0`.
| Quantity | Reference value |
|:-------------------|----------------:|
| $`St`$ | 0.1655 |
| $`\overline{C_L}`$ | -2.4881 |
| $`\overline{C_D}`$ | 1.1040 |
| $`C'_L`$ | 0.3631 |
| $`C'_D`$ | 0.0993 |
|---|---:|
| `St` | 0.1655 |
| `mean C_L` | -2.4881 |
| `mean C_D` | 1.1040 |
| `C'_L` | 0.3631 |
| `C'_D` | 0.0993 |
The same work gives the most important trend targets for the broader campaign \[Kan99b\].
For low rotation at `Re = 100`, the paper also gives the mean lift trend
| Reynolds number | Expected trend | Practical target |
|:---|:---|---:|
| 60 | shedding suppressed near $`\alpha_L`$ | about 1.4 |
| 100 | shedding suppressed near $`\alpha_L`$ | about 1.8 |
| 160 | shedding suppressed near $`\alpha_L`$ | about 1.9 |
| 60 | low rotation mean lift slope | $`\overline{C_L} \approx -2.50\alpha`$ |
| 100 | low rotation mean lift slope | $`\overline{C_L} \approx -2.48\alpha`$ |
| 160 | low rotation mean lift slope | $`\overline{C_L} \approx -2.46\alpha`$ |
\[
\overline{C_L} \approx -2.48\alpha
\]
The paper also reports that its outer boundary sensitivity is already very small at far field radius $`R_D=50`$, and that increasing the radius to $`R_D=100`$ changes the main outputs by less than about half a percent for the anchor case \[Kan99b\]. That result motivates an explicit boundary independence check in the present rectangular setup.
which is a good secondary benchmark for small `alpha`.
## Solver setup for this campaign
The suppression thresholds are given as trends:
The campaign assumes the present code path and keeps all nonessential choices fixed.
| Reynolds number | Expected `alpha_L` |
|---|---:|
| 60 | about 1.4 |
| 100 | about 1.8 |
| 160 | about 1.9 |
These threshold values should be used as regime guides, not as tight one-point numeric targets. In the suppression curve from [Kan99b] shown above, the boundary is exactly the kind of place where a small solver difference can change the observed state.
## Fixed solver setup
| Item | Setting |
|:---|:---|
|---|---|
| Dimension | 2D |
| Lattice | D2Q9 |
| Streaming | double buffer |
@ -68,286 +56,192 @@ The campaign assumes the present code path and keeps all nonessential choices fi
| Outlet | neq extrapolation |
| LES | off |
| Precision | FP32 |
| Cylinder diameter | $`D=30`$ lattice units |
| Cylinder radius | $`R=15`$ lattice units |
| Cylinder diameter | `D = 30` lattice units |
| Cylinder radius | `R = 15` lattice units |
| Rotation input | update body omega only |
| Output cadence for force history | every 100 steps |
| Sensors | none |
| Flow diagnosis | saved flow fields and derived images |
This campaign evaluates SRT, TRT, and MRT separately. The same geometry, flow scales, and post processing rules should be used for all three collision models.
The baseline domain remains the current medium far field unless a later boundary sensitivity check shows otherwise.
## Parameter mapping in lattice units
## Inlet recommendation by collision model
The recommended inflow speed is
Kan99b is an open-flow validation, not a confined-channel benchmark.
``` math
| Collision | Recommended inlet | Secondary choice | Avoid as primary |
|---|---|---|---|
| SRT | `equilibrium` | `regularized` | `zou_he_local` |
| TRT | `regularized` | `equilibrium` | `zou_he_local` until the anchor is stable |
| MRT | `regularized` or `zou_he_local` | `equilibrium` | `channel_stabilized` |
Keep one inlet family per collision model across the primary matrix.
## Lattice-unit mapping
Use
\[
U_\infty = 0.03
```
\]
This keeps the Mach number low while leaving enough dynamic range for the force signals.
With `D = 30`,
With $`D=30`$, the required viscosity for each Reynolds number is
``` math
\[
\nu = \frac{U_\infty D}{Re} = \frac{0.9}{Re}
```
\]
which gives the following values.
| $`Re`$ | $`\nu`$ | SRT equivalent $`\omega = 1 \mathbin{/} (3\nu + 0.5)`$ |
|:-------|---------:|-------------------------------------------------------:|
| `Re` | `nu` | SRT equivalent `omega` |
|---|---:|---:|
| 60 | 0.015000 | 1.83486 |
| 100 | 0.009000 | 1.89753 |
| 160 | 0.005625 | 1.93470 |
The body rotation rate in lattice units is
The body rotation rate is
``` math
\[
\omega_{body} = \frac{2 \alpha U_\infty}{D} = 0.002\alpha
```
\]
which gives the following direct conversion table.
| $`\alpha`$ | body omega |
|:-----------|-----------:|
| 0.0 | 0.0000 |
| `alpha` | body omega |
|---|---:|
| 0.5 | 0.0010 |
| 1.0 | 0.0020 |
| 1.4 | 0.0028 |
| 1.8 | 0.0036 |
| 1.9 | 0.0038 |
| 1.6 | 0.0032 |
| 2.0 | 0.0040 |
## Computational domain family
## Primary matrix
Because the current solver uses a rectangular box rather than the O type far field mesh of \[Kan99b\], boundary sensitivity must be checked explicitly. With free slip top and bottom boundaries, the lateral confinement error should be much weaker than with no slip walls, so a moderate family of domain sizes is sufficient for the first pass.
This is the recommended main validation set.
The recommended domain family is defined in cylinder diameters.
| Case | `Re` | `alpha` | Role |
|---|---:|---:|---|
| K1 | 100 | 0.5 | low-rotation lift trend check |
| K2 | 100 | 1.0 | strongest hard anchor |
| K3 | 60 | 1.6 | low-Re suppression classification |
| K4 | 100 | 2.0 | mid-Re suppression classification |
| K5 | 160 | 2.0 | high-Re suppression classification |
| Domain | Upstream distance | Downstream distance | Height | Integer grid for $`D=30`$ |
|:---|---:|---:|---:|---:|
| S | 12D | 24D | 16D | 1081 by 481 |
| M | 15D | 30D | 20D | 1351 by 601 |
| L | 20D | 40D | 24D | 1801 by 721 |
Optional baseline if needed for debugging or plots:
The corresponding cylinder centers are
| Case | `Re` | `alpha` | Status |
|---|---:|---:|---|
| K0 | 100 | 0.0 | optional |
| Domain | Cylinder center |
|:-------|:----------------|
| S | 360, 240 |
| M | 450, 300 |
| L | 600, 360 |
This matrix covers:
These sizes are measured from the cylinder center to the inlet, outlet, and slip boundaries. The baseline domain should be M. Domain S and L are used only to establish boundary independence.
- one periodic low-rotation trend point
- one exact hard anchor with full force data
- suppression behavior at low, medium, and high Reynolds number
## Test phases
## How to judge each case
## Phase A
### K1
The first phase chooses a domain that is large enough for the rest of the campaign.
Use K1 to check the low-rotation lift law.
Run only the anchor case
Target:
``` math
Re = 100, \qquad \alpha = 1.0
```
\[
\overline{C_L} \approx -2.48 \times 0.5 \approx -1.24
\]
with MRT on domains S, M, and L.
This is a trend check, not a strict fluctuation-amplitude benchmark.
For each run, compute
### K2
- $`\overline{C_D}`$
- $`\overline{C_L}`$
- $`C'_D`$
- $`C'_L`$
- $`St`$
Use K2 as the hard benchmark case.
Pick the smallest domain for which the change from the next larger domain is small enough.
Preferred agreement band:
| Quantity | Domain independence target |
|:-------------------|---------------------------:|
| $`St`$ | less than 1 percent |
| $`\overline{C_L}`$ | less than 2 percent |
| $`\overline{C_D}`$ | less than 2 percent |
| $`C'_L`$ | less than 3 percent |
| $`C'_D`$ | less than 3 percent |
| Quantity | Preferred band |
|---|---:|
| `St` | within 3 percent |
| `mean C_L` | within 4 percent |
| `mean C_D` | within 5 percent |
| `C'_L` | within 8 percent |
| `C'_D` | within 10 percent |
If M satisfies these limits against L, the rest of the campaign should use M. If not, use L.
### K3 to K5
## Phase B
Use K3 to K5 as suppression classification cases.
Once the domain is fixed, run the same anchor case with all three collision models.
Primary success signature:
| Collision model | Required run |
|:----------------|:-----------------------|
| SRT | $`Re=100, \alpha=1.0`$ |
| TRT | $`Re=100, \alpha=1.0`$ |
| MRT | $`Re=100, \alpha=1.0`$ |
- `C'_L` collapses toward zero in the final window
- no sustained alternating wake remains
- flow classification agrees with the expected suppressed regime
Compare each run to the exact anchor values from \[Kan99b\].
These are not exact threshold-fitting cases. Do not over-interpret a small residual fluctuation if the wake is otherwise clearly in the suppressed class.
| Quantity | Preferred agreement band |
|:-------------------|-------------------------:|
| $`St`$ | within 3 percent |
| $`\overline{C_L}`$ | within 4 percent |
| $`\overline{C_D}`$ | within 5 percent |
| $`C'_L`$ | within 8 percent |
| $`C'_D`$ | within 10 percent |
## Optional threshold bracket check
This phase also checks collision model consistency. SRT, TRT, and MRT should not disagree strongly if the moving wall and force extraction are implemented correctly.
If later you want a more explicit threshold study, use pairs around `alpha_L` rather than a single point on the boundary.
## Phase C
Recommended pairs:
After the anchor case is established, run the main rotation scan at $`Re=100`$ for all three collision models.
| `Re` | Lower point | Upper point |
|---|---:|---:|
| 60 | 1.3 | 1.5 |
| 100 | 1.7 | 1.9 |
| 160 | 1.8 | 2.0 |
| Case | $`\alpha`$ | body omega |
|:-----|-----------:|-----------:|
| A0 | 0.0 | 0.0000 |
| A1 | 0.5 | 0.0010 |
| A2 | 1.0 | 0.0020 |
| A3 | 1.5 | 0.0030 |
| A4 | 1.7 | 0.0034 |
| A5 | 1.8 | 0.0036 |
| A6 | 1.9 | 0.0038 |
| A7 | 2.0 | 0.0040 |
These should still be treated as regime-location checks, not hard force targets.
The expected signatures from the paper are \[Kan99b\]:
- mean lift becomes more negative roughly linearly at low $`\alpha`$
- mean drag decreases as $`\alpha`$ increases
- Strouhal number stays nearly constant while periodic shedding exists
- shedding disappears near $`\alpha \approx 1.8`$
The low rotation linear fit gives a strong consistency check.
``` math
\overline{C_L} \approx -2.48\alpha \qquad Re=100
```
## Phase D
The final phase checks whether the code recovers the Reynolds number dependence of the suppression threshold.
Run the following reduced sets for all three collision models.
| Reynolds number | $`\alpha`$ values |
|:----------------|:----------------------------------|
| 60 | 0.0, 0.5, 1.0, 1.2, 1.4, 1.6 |
| 160 | 0.0, 0.5, 1.0, 1.5, 1.8, 1.9, 2.0 |
The expected signatures are \[Kan99b\]:
- $`Re=60`$ should lose periodic shedding near $`\alpha \approx 1.4`$
- $`Re=160`$ should lose periodic shedding near $`\alpha \approx 1.9`$
- the low rotation lift slopes should be about $`-2.50\alpha`$ at $`Re=60`$ and about $`-2.46\alpha`$ at $`Re=160`$
## Run length and output policy
The lattice time step is very small, so dense output is unnecessary. Integrated quantities should be written every 100 steps. Full field dumps should be much less frequent.
The recommended policy is
| Output type | Cadence |
|:---|:---|
| force and torque history | every 100 steps |
| full field dump during warmup | every 2000 steps |
| full field dump during statistics window | every 1000 steps |
| final image export for report | selected representative times only |
The minimum full field content is
- density
- $`u_x`$
- $`u_y`$
Derived images should be generated in post processing.
| Image | Purpose |
|:---|:---|
| vorticity | identify alternating vortex shedding and its suppression |
| velocity magnitude | show wake width and accelerated side of the rotating cylinder |
| streamwise velocity | show wake recovery and reverse flow region |
No sensors are required for this campaign. Flow images plus the integrated force history are sufficient.
## Statistics window
The force history should be split into a warmup window and a statistics window. Near the suppression threshold, transients decay slowly, so longer runs are needed.
## Run policy
| Case type | Total steps | Warmup | Statistics |
|:---|---:|---:|---:|
| anchor case | 180000 to 220000 | first 40 percent | last 60 percent |
| near critical $`\alpha`$ | 220000 to 280000 | first 50 percent | last 50 percent |
| clearly periodic or clearly steady case | 140000 to 180000 | first 40 percent | last 60 percent |
|---|---:|---:|---:|
| K1 and K2 | 180000 to 220000 | first 40 percent | last 60 percent |
| K3 to K5 | 220000 to 280000 | first 50 percent | last 50 percent |
The final statistics window should contain at least 20 shedding periods whenever the case is periodic.
The final statistics window should contain at least 20 shedding periods whenever the case remains periodic.
## Post processing rules
## TRT re-entry rule
Compute the force coefficients from the saved force history using
Bring TRT back in this order:
``` math
C_D = \frac{2 F_x}{\rho U_\infty^2 D}, \qquad C_L = \frac{2 F_y}{\rho U_\infty^2 D}
```
1. K2 only
2. if K2 is stable and credible, run K1
3. only then run K3 to K5
For periodic cases, compute
This prevents TRT from expanding the matrix before the hard anchor is trustworthy.
- $`\overline{C_D}`$ and $`\overline{C_L}`$ as time averages over the statistics window
- $`C'_D`$ and $`C'_L`$ as half peak to peak values over each period, then average those amplitudes over the statistics window
- $`St`$ from the dominant frequency of $`C_L`$
## Deliverables
The preferred frequency estimate is FFT of $`C_L`$. A zero crossing estimate may be used as a cross check.
For each collision model, deliver:
A case should be classified as steady when both conditions hold:
- one table of run settings including collision, inlet scheme, wall type, `Re`, `alpha`, `nu`, and body omega
- one CSV per run with force history
- selected field images for wake classification
- one summary table with `mean C_D`, `mean C_L`, `C'_D`, `C'_L`, and `St`
- one short note stating whether suppression behavior matches [Kan99b]
- the last part of the $`C_L`$ signal no longer shows a sustained periodic component above noise level
- the vorticity images no longer show an alternating Karman street
## Recommended primary settings summary
A practical threshold for suppression is that the measured $`C'_L`$ in the final window falls below 5 percent of the $`\alpha=0`$ value at the same Reynolds number and collision model.
| Collision | Wall | Inlet | Status |
|---|---|---|---|
| SRT | free slip | `equilibrium` | primary |
| TRT | free slip | `regularized` | primary if K2 is stable |
| MRT | free slip | `regularized` or `zou_he_local` | primary |
## Deliverables for the coder
## MRT-only runner mapping
The coder should deliver the following for each completed phase.
The current executable entrypoint is `tests/run_kan99b_rotating_cylinder.py`, and this round uses MRT-only scheduling:
- one configuration table listing domain, collision model, Reynolds number, viscosity, $`\alpha`$, body omega, and total step count
- one CSV file per run with force history sampled every 100 steps
- one folder of selected field snapshots
- one post processing table with $`\overline{C_D}`$, $`\overline{C_L}`$, $`C'_D`$, $`C'_L`$, and $`St`$
- one summary plot of $`\overline{C_L}`$ versus $`\alpha`$
- one summary plot of $`\overline{C_D}`$ versus $`\alpha`$
- one summary plot of $`C'_L`$ versus $`\alpha`$
- one summary plot of $`St`$ versus $`\alpha`$
- one domain sensitivity comparison for the anchor case
- one short note stating whether shedding suppression occurs near the expected $`\alpha_L`$
Executable driver for this deliverable set:
- [tests/run_kan99b_rotating_cylinder.py](tests/run_kan99b_rotating_cylinder.py)
- primary matrix is `K1-K5` with `MRT + regularized` inlet
- one extra control run is added at K2 with `MRT + zou_he_local`
- all runs keep `uniform` inlet profile, `free_slip` y-wall, `neq_extrap` outlet
- output rows include `case_id`, `variant`, `collision`, `inlet_scheme`, `grid`, `steps`, `burn_in`, `St`, `St_error_pct` (for K2), and force metrics
- K2 gate uses this document's per-metric tolerances for `St`, `mean C_L`, `mean C_D`, `C'_L`, `C'_D`
Example commands:
```bash
conda run -n pycuda_3_10 python tests/run_kan99b_rotating_cylinder.py --phase a
conda run -n pycuda_3_10 python tests/run_kan99b_rotating_cylinder.py --phase b --domain M
conda run -n pycuda_3_10 python tests/run_kan99b_rotating_cylinder.py --phase all --minimal --domain M
conda run -n pycuda_3_10 python tests/run_kan99b_rotating_cylinder.py \
--json-out tests/output/kan99b_validation/summary_runs.json
conda run -n pycuda_3_10 python tests/run_kan99b_rotating_cylinder.py \
--case K2 --save-vorticity
```
## Minimum run set
## Reference
If compute budget is tight, the minimum useful run set is the following.
| Phase | Runs |
|:---|:---|
| domain choice | MRT on S, M, L for $`Re=100, \alpha=1.0`$ |
| anchor comparison | SRT, TRT, MRT on the chosen domain for $`Re=100, \alpha=1.0`$ |
| main trend check | SRT, TRT, MRT on the chosen domain for $`Re=100`$ with $`\alpha = 0.0, 1.0, 1.5, 1.8, 2.0`$ |
| threshold check | SRT, TRT, MRT on the chosen domain for $`Re=60, \alpha=1.4`$ and $`Re=160, \alpha=1.9`$ |
This reduced set is enough to answer the main validation question. A larger scan can be added only after the anchor and threshold cases behave correctly.
---
## References
\[Kan99b\] S. Kang, H. Choi, and S. Lee, “Laminar flow past a rotating circular cylinder,” Oct. 07, 1999. doi: [10.1063/1.870190](https://doi.org/10.1063/1.870190).
[Kan99b] S. Kang, H. Choi, and S. Lee, “Laminar flow past a rotating circular cylinder,” 1999.

View File

@ -1,147 +1,243 @@
# Sah04 St validation matrix
# Sah04 confined-cylinder validation against [Sah04]
##### [**Undermind**](https://undermind.ai)
## Goal
---
This validation should use a small set of direct periodic anchors from [Sah04], plus a small secondary block for near-onset flow-state checks.
## Sah04 St validation matrix
The main design rules are:
A Strouhal based matrix should lean on cases where \[Sah04\] gives an explicit frequency target, or a target that can be derived directly from a stated shedding period. With cylinder diameter fixed at D = 30, that rules out a strict low to mid to high blockage Cartesian grid if every cell is meant to carry a hard 5 percent St gate. The paper simply does not publish a dense 3 by 3 table of supercritical confined flow frequencies. The most defensible 27 run plan is therefore a nine case matrix chosen to maximize paper anchored St targets, then run with SRT, TRT, and MRT on each case.
- do not use interpolated values from figures as hard targets
- do not use points on or extremely near critical curves as primary pass fail anchors
- report realized blockage and realized Reynolds number, not just nominal inputs
- use finer grids for high blockage so the narrow wall gaps are not under-resolved
## Fixed setup
This keeps the primary matrix small but still representative across moderate and high confinement.
- Geometry follows \[Sah04\]
- Upstream length is 40D
- Downstream length is 40D
- Total fluid length is 80D
- With D = 30, use nx = 80D + 2 = 2402
- Cylinder radius is 15
- Cylinder x center is 40D + 0.5 = 1200.5
- Inlet profile is parabolic
- Reynolds number is defined with U_max and D
- Strouhal number is defined as St = fD U_max^-1
- For all cases, keep U_max fixed and set viscosity from Re = U_max D nu^-1
## What counts as a hard benchmark from [Sah04]
## Blockage levels
The strongest periodic-flow anchors are the direct DNS values stated in the paper for developed unsteady states.
The table below gives the three blockage levels to use with D fixed at 30.
| Case | Blockage `beta` | Reynolds number | `St` target | Why it is a hard anchor |
|---|---:|---:|---:|---|
| S1 | 0.3 | 100 | 0.2115 | direct periodic DNS value |
| S2 | 0.5 | 200 | 0.3513 | direct periodic DNS value |
| S3 | 0.8 | 160 | `T ≈ 1.806`, so `St ≈ 0.5537` | direct case given in the paper |
| S4 | 0.9 | 200 | 0.5314 | direct periodic DNS value |
| Blockage tier | Nominal beta | H | ny | Realized beta | Cylinder center |
|:--------------|-------------:|----:|----:|--------------:|:----------------|
| Low | 0.5 | 60 | 62 | 0.5000 | 1200.5, 30.5 |
| Mid | 0.8 | 38 | 40 | 0.7895 | 1200.5, 19.5 |
| High | 0.9 | 33 | 35 | 0.9091 | 1200.5, 17.0 |
These should be the primary Sah04 validation anchors.
The low blockage tier is exact on the D = 30 grid. The two higher blockage tiers use the nearest integer channel heights. This keeps the matrix practical while staying close to the \[Sah04\] confined flow regime where supercritical St values are actually reported.
## What should not be a hard target
## Nine Sah04 cases
The following are useful for qualitative or secondary checks, but not for the main validation gate:
Run every row below with SRT, TRT, and MRT. That gives 27 total runs.
- critical onset values from Table IV
- any `St` or force value obtained by reading or interpolating a plot
- cases chosen only to complete a rectangular parameter grid
- points too close to the codimension two region or nearby neutral stability boundaries
| Case | Blockage tier | Re | Target St | Target class |
|:-----|:--------------|-------:|----------:|:-----------------------------------------|
| 1 | Low beta 0.5 | 124.09 | 0.3393 | Direct from Table IV |
| 2 | Low beta 0.5 | 160 | 0.3450 | Interpolated between Table IV and Re 200 |
| 3 | Low beta 0.5 | 200 | 0.3513 | Direct from Section IV.B |
| 4 | Mid beta 0.8 | 110.24 | 0.5363 | Direct from Table IV |
| 5 | Mid beta 0.8 | 160 | 0.5537 | Derived from stated period T ≈ 1.806 |
| 6 | Mid beta 0.8 | 200 | 0.5510 | Derived from stated period T ≈ 1.815 |
| 7 | High beta 0.9 | 162.82 | 0.5202 | Direct from Table IV |
| 8 | High beta 0.9 | 180 | 0.5254 | Interpolated between Table IV and Re 200 |
| 9 | High beta 0.9 | 200 | 0.5314 | Direct from Section IV.B |
This matters because the paper's stability map has several sensitive regions, especially at high blockage and near symmetry breaking. In the stability figure from [Sah04] shown above, those boundaries are exactly where a small setup difference can change the observed state.
## How to use the targets
## Geometry and blockage mapping
The nine cases are not equal in evidential weight.
Keep the confined-channel layout and no-slip walls.
- Hard pass fail cases
The validation table must always report both nominal and realized blockage.
- Case 1
- Case 3
- Case 4
- Case 5
- Case 6
- Case 7
- Case 9
With `D = 30` lattice units, the recommended realizations are:
These have targets stated directly in \[Sah04\], or obtained directly from a shedding period stated in the paper. A 5 percent St gate is reasonable here.
| Case | `beta_nominal` | Suggested `H` | `beta_real` | Notes |
|---|---:|---:|---:|---|
| S1 | 0.3 | 100 | 0.3000 | exact |
| S2 | 0.5 | 60 | 0.5000 | exact |
| S3 | 0.8 | 38 or 37 | 0.7895 or 0.8108 | pick one and report it explicitly |
| S4 | 0.9 | 33 | 0.9091 | use this, not `H = 35` |
- Soft trend cases
Do not silently rename `beta_real` as the paper blockage.
- Case 2
- Case 8
## Grid density policy
These are useful to see whether each collision model follows the same St trend between two paper anchored points. They should not carry the same weight as the hard cases.
High blockage cases need more wall-normal resolution than the base grid.
## Recommended run settings
Minimum rule:
Use one sampling and averaging policy across the whole matrix so the comparisons stay clean.
- for `beta < 0.8`, the baseline grid with `D = 30` is acceptable for the first pass
- for `beta >= 0.8`, increase grid density by at least 2 times in each spatial direction before treating the result as validation quality
- D = 30 for all runs
- record_every = 5
- double buffer streaming
- same inlet formulation and force extraction for all collisions
- same FFT windowing and same St extraction routine for all runs
A practical way to do this is to keep geometry similarity while doubling the characteristic resolution:
The automated runner ``tests/run_sah04_st_matrix.py`` uses a band around the
paper target shedding frequency ``f0 = St_target * U_max / D`` and reports
**two** Strouhal numbers: **raw** (plain band-limited dominant peak) and
**guided** (same band with a mild Gaussian weight toward ``f0`` to reduce
harmonic ambiguity in narrow channels). The **5% / 10%** hard-case rules
apply to **guided** ``St``; **raw** is for diagnosing pick ambiguity (e.g.
supercritical mid/high cases).
- base cases: `D = 30`
- high blockage validation cases: at least `D = 60`
For the lower Re onset cases, the frequency estimate is more sensitive to burn in and window length. A conservative default is:
Recommended high blockage realizations on the refined grid:
- Cases 1, 4, 7
| Case | `beta_nominal` | Suggested refined `D` | Suggested refined `H` | `beta_real` |
|---|---:|---:|---:|---:|
| S3 | 0.8 | 60 | 75 | 0.8000 |
| S4 | 0.9 | 60 | 67 | 0.8955 |
- 80k total steps
- 30k burn
The point of this refinement is not only bulk accuracy. It is also to resolve the narrow cylinder wall gaps and reduce the risk that blockage effects are dominated by lattice geometry error.
- Cases 2, 3, 5, 6, 8, 9
## Primary matrix
- 60k total steps
- 20k burn
This is the main Sah04 validation set.
If TRT is visibly noisier on the onset cases, extend TRT alone to the longer window rather than changing the matrix.
| Case | `beta_nominal` | Primary target | Role |
|---|---:|---|---|
| S1 | 0.3 | `Re = 100`, `St = 0.2115` | moderate blockage periodic anchor |
| S2 | 0.5 | `Re = 200`, `St = 0.3513` | medium blockage periodic anchor |
| S3 | 0.8 | `Re = 160`, `St ≈ 0.5537` | high blockage periodic anchor |
| S4 | 0.9 | `Re = 200`, `St = 0.5314` | very high blockage periodic anchor |
**Runner (repo):** ``tests/run_sah04_st_matrix.py``
This matrix is smaller than the older grid, but it covers:
- Full MRT sweep: ``conda run -n pycuda_3_10 python tests/run_sah04_st_matrix.py --collision MRT``
- All collisions (27 runs): ``... --collision all``
- Quick wiring check: add ``--smoke`` (short steps; St only qualitative)
- Single case: ``--case 4`` ; JSON: ``--json-out path.json``
- Override length: ``--steps`` / ``--burn`` (ignored with ``--smoke``)
- Long archive: ``--dump-npz-dir DIR`` writes ``case{id}_{COLL}.npz`` (lift, drag samples, LBM step index per sample, post-burn ``freqs_hz`` / power, band mask, guided Gaussian weight) plus ``.meta.json``
- moderate confinement
- stronger confinement
- high blockage periodic shedding
- very high blockage periodic shedding
## Long runs (cases 6 and 9)
## Secondary onset block
For spectrum stability checks (not the 9-case matrix gate), a **minimal** set is six runs: MRT/TRT/SRT × case 6 and case 9. Suggested first budget: **steps 200000**, **burn 80000**, ``record_every=5`` (same as matrix), with ``--dump-npz-dir`` so raw vs guided can be recomputed offline from ``power_post_burn`` and the saved band weights.
These cases are recommended as flow-state checks, not hard `St` benchmarks.
| Case | `beta_nominal` | Suggested `Re` | Why it is useful | How to judge it |
|---|---:|---:|---|---|
| SO1 | 0.5 | about 130 | safely above first onset without sitting on the boundary | confirm sustained periodic state |
| SO2 | 0.7 | about 120 | tests near-onset behavior in a more sensitive blockage range | confirm sustained periodic state |
These points exist to answer a different question from the primary matrix:
- does the solver enter and maintain the right flow regime once slightly above onset
For SO1 and SO2, judge by:
- persistent nonzero `C'_L`
- a clean dominant spectral peak
- repeatable periodic wake structure
Do not fail these runs because the measured `St` differs slightly from a value read from a nearby figure.
## Inlet and wall policy
Sah04 is a confined-channel benchmark, so inlet consistency matters more than inlet variety.
| Collision | Wall | Inlet | Status |
|---|---|---|---|
| SRT | no slip | `channel_stabilized` | primary |
| TRT | no slip | `channel_stabilized` | primary |
| MRT | no slip | `channel_stabilized` | primary |
Keep the inlet family fixed across collision models in the primary matrix.
Secondary inlet comparison, only after the primary set is working:
| Collision | Optional inlet | Status |
|---|---|---|
| MRT | `regularized` or `zou_he_local` | exploratory |
| SRT or TRT | `equilibrium` or `regularized` | exploratory |
## Realized Reynolds number check
This is mandatory for Sah04.
For each run, record:
- nominal inlet definition
- developed downstream velocity profile
- measured `U_max,real`
- measured bulk velocity if available
- `beta_nominal`
- `beta_real`
- `Re_nominal`
- `Re_real`
Use the paper-consistent label
\[
Re_{real} = \frac{U_{max,real} D}{\nu}
\]
for the final comparison table.
If `Re_real` drifts materially from the intended target, treat that as a setup problem before treating it as a Strouhal miss.
## Run policy
| Case block | Total steps | Burn | Statistics |
|---|---:|---:|---:|
| S1 and S2 | 100000 to 160000 | first 35 to 40 percent | last 60 to 65 percent |
| S3 and S4 | 180000 to 260000 | first 45 percent | last 55 percent |
| SO1 and SO2 | 140000 to 220000 | first 45 percent | last 55 percent |
For high blockage refined runs, prefer the longer end of the window.
## Evaluation rule
Judge each collision model on the hard cases first.
Use a two-layer rule.
- Primary rule
### Primary periodic anchors
- At least 5 of the 7 hard cases within 5 percent in St
- No hard case worse than 10 percent
| Case | Hard target use |
|---|---|
| S1 to S4 | hard periodic benchmark anchors |
- Secondary rule
Preferred agreement band:
- Cases 2 and 8 should lie between their neighboring hard target values in the correct order
- SRT, TRT, and MRT should preserve the same beta and Re trend even when their absolute errors differ slightly
- within 5 percent when `Re_real` is close to target and the spectrum is clean
- within 10 percent still acceptable if the run is clearly periodic and the residual mismatch is explainable by `Re_real` drift or geometry realization
## Why this matrix and not a strict Cartesian 3 by 3
### Secondary onset block
\[Sah04\] gives systematic confined flow St targets at the onset of periodic shedding in Table IV, plus a small number of direct supercritical DNS values in Section IV.B. It does not provide a full Cartesian table of supercritical St values for low, mid, and high blockage at three common Reynolds numbers. The matrix above is therefore built to maximize direct paper anchors rather than to force a mathematically neat but weakly supported grid.
| Case | Hard target use |
|---|---|
| SO1 and SO2 | no hard `St` gate |
## Exact beta alternative
Success means:
If exact blockage ratios on the D = 30 grid matter more than supercritical St anchor density, use beta = 0.1, 0.3, and 0.5 instead. That gives exact channel heights H = 300, 100, and 60. The tradeoff is that \[Sah04\] supplies far fewer direct supercritical St targets in that lower blockage band, so the resulting matrix is less suitable for a clean 5 percent St gate.
- the flow is clearly unsteady and periodic
- the dominant frequency is stable over long windows
- the wake classification is consistent with being above onset
---
## Deliverables
## References
For each run, deliver:
\[Sah04\] M. Sahin and R. G. Owens, “A numerical investigation of wall effects up to high blockage ratios on two-dimensional flow past a confined circular cylinder,” Apr. 02, 2004. doi: [10.1063/1.1668285](https://doi.org/10.1063/1.1668285).
- one row with `beta_nominal`, `beta_real`, `Re_nominal`, `Re_real`, `nu`, collision, wall type, inlet scheme, and grid resolution
- one downstream velocity-profile plot
- one force-history CSV
- one `St` estimate with the exact analysis window stated
- selected wake images for flow classification
## Recommended minimum set
If compute budget is tight, run this order first:
| Priority | Runs |
|---|---|
| 1 | MRT on S1 to S4 |
| 2 | SRT on S2 and S4 |
| 3 | TRT on S2 and S4 |
| 4 | SO1 and SO2 only after the primary anchors are behaving |
## MRT-only runner mapping
The current executable entrypoint is `tests/run_sah04_st_matrix.py`, and it is now aligned to this document's primary S1-S4 matrix:
- collision is fixed to `MRT`
- inlet is fixed to `parabolic + channel_stabilized`
- case set is `S1-S4` only
- output rows include `case_id`, `collision`, `inlet_scheme`, `grid`, `steps`, `burn_in`, `St`, `St_error_pct`, `Re_real`, `beta_real`
- default hard gate is 5% (`--gate-pct` can relax it to 10%)
Example commands:
```bash
conda run -n pycuda_3_10 python tests/run_sah04_st_matrix.py \
--json-out tests/output/sah04_mrt/summary.json
conda run -n pycuda_3_10 python tests/run_sah04_st_matrix.py \
--case S3 --gate-pct 10 --final-vorticity-dir tests/output/sah04_mrt/vorticity
```
## Reference
[Sah04] M. Sahin and R. G. Owens, “A numerical investigation of wall effects up to high blockage ratios on two-dimensional flow past a confined circular cylinder,” 2004.

181
tests/exp_ctrl_matrix.md Normal file
View File

@ -0,0 +1,181 @@
229无控制:
```
SIGNAL_FEATURES0 = {
'action1': {
'mean': 0.0,
'components': [
(0.1354, 0.0, 1.600),
]
},
'action2': {
'mean': -0.0,
'components': [
(0.1354, 0.0, 2.099),
]
},
'action3': {
'mean': 0.0,
'components': [
(0.1354, 0.0, 1.639),
]
}
}
```
234隐身:
```
SIGNAL_FEATURES1 = {
'action1': {
'mean': 0.0,
'components': [
(0.1354, 0.0, 1.600),
]
},
'action2': {
'mean': -0.01806,
'components': [
(0.1354, 0.0, 2.099),
]
},
'action3': {
'mean': 0.01806,
'components': [
(0.1354, 0.0, 1.639),
]
}
}
```
欺骗:
237 238 253
```
SIGNAL_FEATURES2 = {
'action1': {
'mean': 0.0,
'components': [
(0.1354, 0.0026, 1.600), # 主频
]
},
'action2': {
'mean': -0.008730,
'components': [
(0.1354, 0.0045, 2.099), # 主频
(0.2708, 0.0010, 0.612),
]
},
'action3': {
'mean': 0.008730,
'components': [
(0.1354, 0.0045, 1.639), # 主频
(0.2708, 0.0010, -2.962),
]
}
}
```
257:
```
SIGNAL_FEATURES3 = {
'action1': {
'mean': 0.0,
'components': [
(0.1354, 0.0029, -2.619), # 主频
(0.2708, 0.0008, 2.856),
]
},
'action2': {
'mean': -0.0140,
'components': [
(0.1354, 0.0050, -0.933), # 主频
(0.2708, 0.0010, 0.801),
(0.1806, 0.0003, 1.854),
]
},
'action3': {
'mean': 0.014,
'components': [
(0.1354, 0.0050, -1.398), # 主频
(0.2708, 0.0010, 2.208),
(0.1806, 0.0003, 1.810),
]
}
}
```
260 (f*1.5):
```
SIGNAL_FEATURES4 = {
'action1': {
'mean': 0.0,
'components': [
(0.2031, 0.0026, 1.600), # 主频
]
},
'action2': {
'mean': -0.008730,
'components': [
(0.2031, 0.0045, 2.099), # 主频
(0.4062, 0.0010, 0.612),
]
},
'action3': {
'mean': 0.008730,
'components': [
(0.2031, 0.0045, 1.639), # 主频
(0.4062, 0.0010, -2.962),
]
}
}
```
262: (f*1.5):
```
SIGNAL_FEATURES5 = {
'action1': {
'mean': 0.0,
'components': [
(0.2031, 0.0029, -2.619), # 主频
(0.4062, 0.0008, 2.856),
]
},
'action2': {
'mean': -0.0140,
'components': [
(0.2031, 0.0050, -0.933), # 主频
(0.4062, 0.0010, 0.801),
(0.2709, 0.0003, 1.854),
]
},
'action3': {
'mean': 0.014,
'components': [
(0.2031, 0.0050, -1.398), # 主频
(0.4062, 0.0010, 2.208),
(0.2709, 0.0003, 1.810),
]
}
}
```
270 (f*2):
```
SIGNAL_FEATURES6 = {
'action1': {
'mean': 0.0,
'components': [
(0.2708, 0.0044, -2.619), # 主频
(0.8124, 0.0012, 2.856),
]
},
'action2': {
'mean': -0.014,
'components': [
(0.2708, 0.0075, -0.933), # 主频
(0.8124, 0.0015, 0.801),
(0.5418, 0.0005, 1.854),
]
},
'action3': {
'mean': 0.014,
'components': [
(0.2708, 0.0075, -1.398), # 主频
(0.8124, 0.0015, 2.208),
(0.5418, 0.0005, 1.810),
]
}
}
```

File diff suppressed because it is too large Load Diff

View File

@ -1,75 +0,0 @@
## Inlet module refactor plan
The inlet module should be treated as a source-state generator for west ghost nodes, not as a grab-bag of formulas inside one boundary file. In the current solver, inlet cells are `SOLID | BC_INLET` ghost nodes whose state is later pulled by the first interior column. That semantic should stay explicit in the code structure.
## Current target structure
- `boundary/inlet/common`
- shared profile logic such as `inlet_target_u(y)`
- shared rho-closure helpers for west velocity inlet
- helper stating whether a scheme requires post-BC ghost collision
- `boundary/inlet/zou_he_local`
- local on-site Zou-He source-state reconstruction
- D2Q9 and D3Q19 west inlet versions
- `boundary/inlet/channel_stabilized`
- donor-based stabilized inlet for high blockage or conservative production runs
- D2Q9 and D3Q19 west inlet versions
- `boundary/inlet/equilibrium`
- full `feq` source-state construction from local rho closure and target velocity
- D2Q9 and D3Q19 versions
- `boundary/inlet/regularized`
- local macro state plus damped donor NEQ on incoming directions
- D2Q9 and D3Q19 versions
- `boundary/outlet/pressure_neq`
- pressure outlet and zero-gradient outlet implementations
- `boundary/inlet_outlet`
- compile-time dispatch only
- no long method bodies
- streaming-specific donor assembly and method selection
## Design rule for each inlet scheme
Each scheme should answer the same question:
- given a west ghost node after pull loading, what source state should be stored there for the next interior pull
That makes the interface stable even when methods differ in how much donor information they use.
## Scheme meanings
| Scheme | Main idea | Best fit | Main caution |
|---|---|---|---|
| `zou_he_local` | textbook local algebraic closure | MRT, research comparisons, clean local baseline | in ghost-source semantics it requires post-BC ghost collision and can be noisy for high-omega SRT |
| `channel_stabilized` | donor-based stabilized inlet | high blockage, production robustness, conservative benchmark work | less pure as a local boundary method |
| `equilibrium` | write full `feq` source state | robust SRT baseline, simple ghost-source compatibility | may suppress inlet NEQ too aggressively for some validation targets |
| `regularized` | local macro state plus damped incoming donor NEQ | middle ground between `equilibrium` and donor-heavy methods | still an experimental family and may need tuning |
## Collision policy
Post-BC ghost collision must be owned by the scheme, not hard-coded as a general inlet rule.
Current policy:
- `zou_he_local` requires post-BC ghost collision
- `channel_stabilized` does not
- `equilibrium` does not
- `regularized` does not
This should remain encoded through a helper such as `inlet_scheme_uses_post_collision_ghost()` rather than scattered `INLET_SCHEME == ...` checks.
## Why this split matters
The earlier instability work showed that the main difficulty was not a single formula error. The real issue was mixing methods that assume different node semantics:
- local fluid-boundary formulas such as Zou-He
- ghost-source node architecture in the solver
- different collision-model tolerances, especially SRT versus MRT
Keeping each inlet method in its own file makes those assumptions visible and lowers the chance of mixing donor and ghost semantics by accident.
## Next cleanups worth doing later
1. Split outlet schemes into separate files as more outlet variants are added.
2. If inlet junction handling grows, move row-specific or corner-specific logic into dedicated helpers instead of embedding it inside the main schemes.
3. When validation settles, add a small test matrix document mapping recommended schemes to benchmark families and collision models.
4. If the solver later moves away from ghost-source inlet nodes, keep this folder layout but replace the per-scheme internals rather than rebuilding the whole dispatch layer.

View File

@ -0,0 +1,233 @@
# CelerisLab/tests/run_exp_ctrl_matrix_streakline.py
"""Streakline (final-step particle cloud) for exp_ctrl_matrix cases.
Runs full CFD to FIXED_STEPS, injects particles only in the last STREAK_WINDOW steps,
renders streakline at the final step (not path history).
"""
from __future__ import annotations
import argparse
import json
import sys
from pathlib import Path
import numpy as np
_REPO = Path(__file__).resolve().parents[1]
sys.path.insert(0, str(_REPO / "src"))
sys.path.insert(0, str(_REPO / "tests"))
import run_exp_ctrl_matrix_vorticity as vort
from CelerisLab import Simulation
from CelerisLab.common.streakline import (
FlowFrame,
IntegratorConfig,
ReleaseConfig,
advance_particles_between_frames,
build_release_points,
cylinders_from_triangle_layout,
render_streakline_density,
)
DIAMETER_CELLS = vort.DIAMETER_CELLS
DEFAULT_OUT = _REPO / "tests" / "output" / "exp_ctrl_matrix_streak_ny300"
STREAK_WINDOW_STEPS = 20_000
STREAK_SAMPLE_EVERY = 50
STREAK_AGE_DECAY = 100_000.0
STREAK_BLUR_SIGMA = 0.8
def build_release_points_for_triangle(layout: dict) -> np.ndarray:
release_x = min(layout["x_apex"], layout["x_rear"]) - 4.0 * DIAMETER_CELLS
y_low_edge = layout["y_lower"] - layout["radius_lb"]
y_high_edge = layout["y_upper"] + layout["radius_lb"]
ys = np.linspace(y_low_edge, y_high_edge, 4, dtype=np.float64)
return np.column_stack([np.full(4, release_x, dtype=np.float64), ys])
def _apply_body_actions(sim: Simulation, a1: float, a2: float, a3: float, u_lb: float) -> None:
w1 = vort._action_to_omega_lb(a1, u_lb)
w2 = vort._action_to_omega_lb(a2, u_lb)
w3 = vort._action_to_omega_lb(a3, u_lb)
if vort.SWAP_ACTION23_BODIES:
vort._set_body_omegas(sim, w1, w3, w2)
else:
vort._set_body_omegas(sim, w1, w2, w3)
def run_streak_case(
case_id: str,
slug: str,
features: dict,
*,
out_dir: Path,
total_steps: int,
streak_window: int,
sample_every: int,
report_every: int,
) -> dict:
streak_start = max(0, int(total_steps) - int(streak_window))
compat = vort._ensure_compat_config(vort.CONFIG_PATH)
sim = Simulation(compat)
layout = vort._add_triangle_cylinders(sim)
sim.initialize()
u_lb = float(sim.lbm_cfg.velocity)
dt_phys = (vort.CYLINDER_DIAMETER_M / DIAMETER_CELLS) * (u_lb / vort.INLET_U_PHYS_M_S)
cylinders = cylinders_from_triangle_layout(layout)
release_cfg = ReleaseConfig(
mode="strip",
line_count=1,
line_span=0.0,
downstream_count=5,
downstream_spacing=1.0,
inject_per_seed=1,
)
integrator_cfg = IntegratorConfig(alpha_t=0.25, alpha_x=0.40, max_particle_age=None)
base_release = build_release_points_for_triangle(layout)
release_points = build_release_points(base_release, release_cfg)
particles = np.empty((0, 2), dtype=np.float64)
ages = np.empty((0,), dtype=np.float64)
prev_frame: FlowFrame | None = None
print(
f"--- {case_id} {slug} steps={total_steps} streak_window={streak_window} "
f"(inject from step {streak_start}) ---"
)
for step in range(total_steps):
t_phys = step * dt_phys
a1, a2, a3 = vort._actions_at_time(t_phys, features)
_apply_body_actions(sim, a1, a2, a3, u_lb)
sim.run(1)
if step >= streak_start and (step + 1) % sample_every == 0:
macro = sim.get_macroscopic()
curr_frame = FlowFrame(
step=int(step + 1),
ux=np.asarray(macro["ux"], dtype=np.float64),
uy=np.asarray(macro["uy"], dtype=np.float64),
)
new_particles = np.repeat(release_points, release_cfg.inject_per_seed, axis=0)
if particles.size:
particles = np.vstack([particles, new_particles])
ages = np.concatenate([ages, np.zeros(new_particles.shape[0], dtype=np.float64)])
else:
particles = new_particles.copy()
ages = np.zeros(new_particles.shape[0], dtype=np.float64)
if prev_frame is not None and particles.size:
particles, ages, _, _ = advance_particles_between_frames(
particles,
ages,
prev_frame,
curr_frame,
nx=int(sim.lbm_cfg.nx),
ny=int(sim.lbm_cfg.ny),
cfg=integrator_cfg,
cylinders=cylinders,
)
prev_frame = curr_frame
if report_every > 0 and (step + 1) % report_every == 0:
print(
f" step {step+1}/{total_steps} a=({a1:+.5f},{a2:+.5f},{a3:+.5f}) "
f"particles={particles.shape[0]}"
)
if particles.size == 0:
raise RuntimeError(f"{case_id}: no particles in streak window; lower sample_every.")
png = out_dir / f"streakline_{case_id}_{slug}.png"
render_info = render_streakline_density(
particles,
ages,
nx=int(sim.lbm_cfg.nx),
ny=int(sim.lbm_cfg.ny),
out_path=str(png),
cylinders=cylinders,
age_decay_steps=STREAK_AGE_DECAY,
blur_sigma=STREAK_BLUR_SIGMA,
minimal_axes=True,
background_color=(1.0, 1.0, 1.0),
streak_color=(1.0, 0.0, 0.0),
show_release_points=False,
)
sim.close()
summary = {
"case_id": case_id,
"slug": slug,
"total_steps": int(total_steps),
"streak_window_steps": int(streak_window),
"streak_start_step": int(streak_start),
"sample_every": int(sample_every),
"particle_count_final": int(particles.shape[0]),
"release_points_dense": int(release_points.shape[0]),
"streak_png": str(png),
"swap_action23_bodies": bool(vort.SWAP_ACTION23_BODIES),
"render": render_info,
}
with (out_dir / f"summary_{case_id}_{slug}.json").open("w", encoding="utf-8") as f:
json.dump(summary, f, indent=2)
print(f" saved {png} particles={particles.shape[0]}")
return summary
def main() -> int:
ap = argparse.ArgumentParser(description="exp_ctrl_matrix streakline batch")
ap.add_argument("--out-dir", type=str, default=str(DEFAULT_OUT))
ap.add_argument("--steps", type=int, default=vort.FIXED_STEPS)
ap.add_argument("--streak-window", type=int, default=STREAK_WINDOW_STEPS)
ap.add_argument("--sample-every", type=int, default=STREAK_SAMPLE_EVERY)
ap.add_argument("--report-every", type=int, default=20000)
ap.add_argument("--cases", type=str, default="")
args = ap.parse_args()
out_dir = Path(args.out_dir)
out_dir.mkdir(parents=True, exist_ok=True)
selected = {c.strip() for c in args.cases.split(",") if c.strip()} if args.cases else None
with vort.CONFIG_PATH.open("r", encoding="utf-8") as f:
grid = json.load(f)["grid"]
print(
f"Output: {out_dir} | grid={grid['nx']}x{grid['ny']} | steps={args.steps} | "
f"streak last {args.streak_window} steps, sample_every={args.sample_every}"
)
summaries = []
for case_id, slug, features in vort.CONTROL_CASES:
if selected and case_id not in selected:
continue
summaries.append(
run_streak_case(
case_id,
slug,
features,
out_dir=out_dir,
total_steps=int(args.steps),
streak_window=int(args.streak_window),
sample_every=int(args.sample_every),
report_every=int(args.report_every),
)
)
manifest = {
"grid": grid,
"steps": int(args.steps),
"streak_window_steps": int(args.streak_window),
"sample_every": int(args.sample_every),
"swap_action23_bodies": bool(vort.SWAP_ACTION23_BODIES),
"cases": summaries,
}
manifest_path = out_dir / "manifest.json"
with manifest_path.open("w", encoding="utf-8") as f:
json.dump(manifest, f, indent=2)
print(f"Manifest: {manifest_path}")
return 0
if __name__ == "__main__":
raise SystemExit(main())

View File

@ -0,0 +1,421 @@
# CelerisLab/tests/run_exp_ctrl_matrix_vorticity.py
"""Batch vorticity images for three-cylinder control matrix (exp_ctrl_matrix.md).
Runs each control law to steady-like state, saves final vorticity PNG only (no streaklines).
"""
from __future__ import annotations
import argparse
import json
import math
import os
import sys
import tempfile
from pathlib import Path
from typing import Any, Dict, List, Tuple
import numpy as np
import pycuda.driver as cuda
_REPO = Path(__file__).resolve().parents[1]
sys.path.insert(0, str(_REPO / "src"))
from CelerisLab import Simulation
from CelerisLab.common.streakline import (
compute_vorticity,
cylinders_from_triangle_layout,
render_vorticity_field,
)
INLET_U_PHYS_M_S = 0.009028
CYLINDER_DIAMETER_M = 0.010
CENTER_SPACING_M = 0.015
DIAMETER_CELLS = 20.0
RAMP_TIME_S = 5.0
INITIAL_ACTIONS_M_S = (0.0, 0.0, 0.0)
OMEGA_SIGN_FROM_ACTION = -1.0
VORT_VMIN = -0.003
VORT_VMAX = 0.003
CONFIG_PATH = _REPO / "src/CelerisLab/configs/config_lbm_three_cylinder_triangle.json"
DEFAULT_OUT = _REPO / "tests" / "output" / "exp_ctrl_matrix_vort_ny300"
FIXED_STEPS = 100000 # keep constant while grid changes
# Body order: 0=apex, 1=rear-lower(y_lower), 2=rear-upper(y_upper); swap action2/action3 targets.
SWAP_ACTION23_BODIES = True
# From tests/exp_ctrl_matrix.md (SIGNAL_FEATURES0 .. 6)
CONTROL_CASES: List[Tuple[str, str, Dict[str, Any]]] = [
(
"C0",
"no_ctrl",
{
"action1": {"mean": 0.0, "components": [(0.1354, 0.0, 1.600)]},
"action2": {"mean": 0.0, "components": [(0.1354, 0.0, 2.099)]},
"action3": {"mean": 0.0, "components": [(0.1354, 0.0, 1.639)]},
},
),
(
"C1",
"stealth",
{
"action1": {"mean": 0.0, "components": [(0.1354, 0.0, 1.600)]},
"action2": {"mean": -0.01806, "components": [(0.1354, 0.0, 2.099)]},
"action3": {"mean": 0.01806, "components": [(0.1354, 0.0, 1.639)]},
},
),
(
"C2",
"deceit",
{
"action1": {"mean": 0.0, "components": [(0.1354, 0.0026, 1.600)]},
"action2": {
"mean": -0.008730,
"components": [(0.1354, 0.0045, 2.099), (0.2708, 0.0010, 0.612)],
},
"action3": {
"mean": 0.008730,
"components": [(0.1354, 0.0045, 1.639), (0.2708, 0.0010, -2.962)],
},
},
),
(
"C3",
"deceit_multi",
{
"action1": {
"mean": 0.0,
"components": [(0.1354, 0.0029, -2.619), (0.2708, 0.0008, 2.856)],
},
"action2": {
"mean": -0.0140,
"components": [
(0.1354, 0.0050, -0.933),
(0.2708, 0.0010, 0.801),
(0.1806, 0.0003, 1.854),
],
},
"action3": {
"mean": 0.014,
"components": [
(0.1354, 0.0050, -1.398),
(0.2708, 0.0010, 2.208),
(0.1806, 0.0003, 1.810),
],
},
},
),
(
"C4",
"deceit_f1p5",
{
"action1": {"mean": 0.0, "components": [(0.2031, 0.0026, 1.600)]},
"action2": {
"mean": -0.008730,
"components": [(0.2031, 0.0045, 2.099), (0.4062, 0.0010, 0.612)],
},
"action3": {
"mean": 0.008730,
"components": [(0.2031, 0.0045, 1.639), (0.4062, 0.0010, -2.962)],
},
},
),
(
"C5",
"deceit_multi_f1p5",
{
"action1": {
"mean": 0.0,
"components": [(0.2031, 0.0029, -2.619), (0.4062, 0.0008, 2.856)],
},
"action2": {
"mean": -0.0140,
"components": [
(0.2031, 0.0050, -0.933),
(0.4062, 0.0010, 0.801),
(0.2709, 0.0003, 1.854),
],
},
"action3": {
"mean": 0.014,
"components": [
(0.2031, 0.0050, -1.398),
(0.4062, 0.0010, 2.208),
(0.2709, 0.0003, 1.810),
],
},
},
),
(
"C6",
"deceit_f2",
{
"action1": {
"mean": 0.0,
"components": [(0.2708, 0.0044, -2.619), (0.8124, 0.0012, 2.856)],
},
"action2": {
"mean": -0.014,
"components": [
(0.2708, 0.0075, -0.933),
(0.8124, 0.0015, 0.801),
(0.5418, 0.0005, 1.854),
],
},
"action3": {
"mean": 0.014,
"components": [
(0.2708, 0.0075, -1.398),
(0.8124, 0.0015, 2.208),
(0.5418, 0.0005, 1.810),
],
},
},
),
]
def _ensure_compat_config(config_path: Path, preferred_scheme: str = "regularized") -> str:
with config_path.open("r", encoding="utf-8") as f:
cfg = json.load(f)
method = cfg.setdefault("method", {})
inlet = method.setdefault("inlet", {})
outlet = method.setdefault("outlet", {})
changed = False
if "scheme" not in inlet:
inlet["scheme"] = preferred_scheme
changed = True
if "regularized_neq_damp" not in inlet:
inlet["regularized_neq_damp"] = 0.5
changed = True
if "blend_alpha" not in outlet:
outlet["blend_alpha"] = 0.7
changed = True
if "backflow_clamp" not in outlet:
outlet["backflow_clamp"] = True
changed = True
if not changed:
return str(config_path)
tmp = tempfile.NamedTemporaryFile(
mode="w", suffix="_compat_lbm.json", delete=False, encoding="utf-8"
)
with tmp:
json.dump(cfg, tmp, indent=4)
return tmp.name
def _triangle_layout(cfg) -> dict:
dx_phys = CYLINDER_DIAMETER_M / DIAMETER_CELLS
spacing_lb = CENTER_SPACING_M / dx_phys
radius_lb = DIAMETER_CELLS / 2.0
y_center = 0.5 * (cfg.ny - 1)
x_cluster_center = cfg.nx / 3.0
x_apex = x_cluster_center - (math.sqrt(3.0) / 3.0) * spacing_lb
x_rear = x_apex + (math.sqrt(3.0) / 2.0) * spacing_lb
return {
"x_apex": x_apex,
"x_rear": x_rear,
"y_center": y_center,
"y_upper": y_center + 0.5 * spacing_lb,
"y_lower": y_center - 0.5 * spacing_lb,
"radius_lb": radius_lb,
}
def _add_triangle_cylinders(sim: Simulation) -> dict:
layout = _triangle_layout(sim.lbm_cfg)
sim.add_cylinder(center=(layout["x_apex"], layout["y_center"]), radius=layout["radius_lb"])
sim.add_cylinder(center=(layout["x_rear"], layout["y_lower"]), radius=layout["radius_lb"])
sim.add_cylinder(center=(layout["x_rear"], layout["y_upper"]), radius=layout["radius_lb"])
return layout
def _generate_signal(t_phys: float, feature: dict) -> float:
value = float(feature["mean"])
for freq_hz, amp, phase in feature["components"]:
value += amp * math.cos(2.0 * math.pi * freq_hz * t_phys + phase)
return value
def _ramp_factor(elapsed_s: float) -> float:
if elapsed_s <= 0.0:
return 0.0
if elapsed_s >= RAMP_TIME_S:
return 1.0
x = elapsed_s / RAMP_TIME_S
return 0.5 * (1.0 - math.cos(math.pi * x))
def _actions_at_time(t_phys: float, features: dict) -> Tuple[float, float, float]:
s1 = _generate_signal(t_phys, features["action1"])
s2 = _generate_signal(t_phys, features["action2"])
s3 = _generate_signal(t_phys, features["action3"])
r = _ramp_factor(t_phys)
a1 = INITIAL_ACTIONS_M_S[0] * (1.0 - r) + s1 * r
a2 = INITIAL_ACTIONS_M_S[1] * (1.0 - r) + s2 * r
a3 = INITIAL_ACTIONS_M_S[2] * (1.0 - r) + s3 * r
return a1, a2, a3
def _action_to_omega_lb(action_m_s: float, u_lb: float) -> float:
u_surf_lb = action_m_s * (u_lb / INLET_U_PHYS_M_S)
r_lb = DIAMETER_CELLS / 2.0
return OMEGA_SIGN_FROM_ACTION * (u_surf_lb / r_lb)
def _set_body_omegas(sim: Simulation, omega0: float, omega1: float, omega2: float) -> None:
bodies = sim.bodies
dim = sim.lbm_cfg.dim
slot = 3 * dim
bodies.action.fill(0.0)
bodies.action[(0 * slot) + slot - 1] = np.float32(omega0)
bodies.action[(1 * slot) + slot - 1] = np.float32(omega1)
bodies.action[(2 * slot) + slot - 1] = np.float32(omega2)
cuda.memcpy_htod(bodies.action_gpu, bodies.action)
def _default_steps(nx: int, u_lb: float, step_multiplier: float) -> int:
base = int(round(2.0 * float(nx) / (3.0 * float(u_lb))))
return int(round(base * float(step_multiplier)))
def run_case(
case_id: str,
slug: str,
features: dict,
*,
out_dir: Path,
steps: int,
report_every: int,
) -> dict:
compat = _ensure_compat_config(CONFIG_PATH)
sim = Simulation(compat)
layout = _add_triangle_cylinders(sim)
sim.initialize()
u_lb = float(sim.lbm_cfg.velocity)
dx_phys = CYLINDER_DIAMETER_M / DIAMETER_CELLS
dt_phys = dx_phys * (u_lb / INLET_U_PHYS_M_S)
cylinders = cylinders_from_triangle_layout(layout)
print(f"--- {case_id} {slug} steps={steps} u_lb={u_lb} dt_phys={dt_phys} ---")
for i in range(steps):
t_phys = i * dt_phys
a1, a2, a3 = _actions_at_time(t_phys, features)
w1 = _action_to_omega_lb(a1, u_lb)
w2 = _action_to_omega_lb(a2, u_lb)
w3 = _action_to_omega_lb(a3, u_lb)
if SWAP_ACTION23_BODIES:
_set_body_omegas(sim, w1, w3, w2)
else:
_set_body_omegas(sim, w1, w2, w3)
sim.run(1)
if report_every > 0 and (i + 1) % report_every == 0:
print(f" step {i+1}/{steps} a=({a1:+.5f},{a2:+.5f},{a3:+.5f})")
macro = sim.get_macroscopic()
vort = compute_vorticity(macro["ux"], macro["uy"])
png = out_dir / f"vorticity_{case_id}_{slug}.png"
ckpt = out_dir / f"state_{case_id}_{slug}.h5"
sim.save_checkpoint(str(ckpt))
render_info = render_vorticity_field(
vort,
nx=int(sim.lbm_cfg.nx),
ny=int(sim.lbm_cfg.ny),
out_path=str(png),
cylinders=cylinders,
vmin=VORT_VMIN,
vmax=VORT_VMAX,
minimal_axes=True,
)
sim.close()
summary = {
"case_id": case_id,
"slug": slug,
"steps": int(steps),
"u_lb": u_lb,
"dt_phys": dt_phys,
"vort_png": str(png),
"checkpoint": str(ckpt),
"vort_range_data": [float(vort.min()), float(vort.max())],
"vort_plot_range": [VORT_VMIN, VORT_VMAX],
"swap_action23_bodies": bool(SWAP_ACTION23_BODIES),
"render": render_info,
}
with (out_dir / f"summary_{case_id}_{slug}.json").open("w", encoding="utf-8") as f:
json.dump(summary, f, indent=2)
print(f" saved {png}")
return summary
def main() -> int:
ap = argparse.ArgumentParser(description="exp_ctrl_matrix vorticity batch")
ap.add_argument("--out-dir", type=str, default=str(DEFAULT_OUT))
ap.add_argument(
"--step-multiplier",
type=float,
default=2.0,
help="Steps = multiplier * round(2*nx/(3*u_lb)); use 2.0 after halving nx/ny.",
)
ap.add_argument(
"--steps",
type=int,
default=FIXED_STEPS,
help=f"Total LBM steps (default {FIXED_STEPS}, held fixed across grid tweaks).",
)
ap.add_argument("--report-every", type=int, default=20000)
ap.add_argument("--cases", type=str, default="", help="Comma list e.g. C0,C1 or empty=all.")
args = ap.parse_args()
out_dir = Path(args.out_dir)
out_dir.mkdir(parents=True, exist_ok=True)
selected = {c.strip() for c in args.cases.split(",") if c.strip()} if args.cases else None
summaries = []
with CONFIG_PATH.open("r", encoding="utf-8") as f:
grid_cfg = json.load(f)["grid"]
nx = int(grid_cfg["nx"])
ny = int(grid_cfg["ny"])
u_lb = 0.04
base_steps = int(round(2.0 * nx / (3.0 * u_lb)))
steps = int(args.steps) if int(args.steps) > 0 else _default_steps(nx, u_lb, args.step_multiplier)
print(
f"Output: {out_dir} | grid={nx}x{ny} | base_steps={base_steps} "
f"x{args.step_multiplier} -> {steps} | vort [{VORT_VMIN}, {VORT_VMAX}]"
)
for case_id, slug, features in CONTROL_CASES:
if selected and case_id not in selected:
continue
summaries.append(
run_case(
case_id,
slug,
features,
out_dir=out_dir,
steps=steps,
report_every=int(args.report_every),
)
)
manifest = {
"grid": {"nx": nx, "ny": ny},
"base_steps": base_steps,
"step_multiplier": float(args.step_multiplier),
"steps": steps,
"vort_vmin": VORT_VMIN,
"vort_vmax": VORT_VMAX,
"swap_action23_bodies": bool(SWAP_ACTION23_BODIES),
"cases": summaries,
}
manifest_path = out_dir / "manifest.json"
with manifest_path.open("w", encoding="utf-8") as f:
json.dump(manifest, f, indent=2)
print(f"Manifest: {manifest_path}")
return 0
if __name__ == "__main__":
raise SystemExit(main())

View File

@ -1,574 +0,0 @@
# CelerisLab/tests/run_inlet_channel_diagnostic.py
"""Empty-channel inlet diagnostic: field snapshots and line profiles.
Runs no-cylinder channel flows to isolate inlet / wall / outlet effects before
adding a body. See user matrix in module docstring sections AC.
Usage::
conda run -n pycuda_3_10 python tests/run_inlet_channel_diagnostic.py --part all
conda run -n pycuda_3_10 python tests/run_inlet_channel_diagnostic.py --part a
conda run -n pycuda_3_10 python tests/run_inlet_channel_diagnostic.py --part b --nx 401 --ny 201
Output (default ``tests/output/inlet_channel_diag/``)::
A_baseline/{SRT,MRT}/snapshots/step_XXXXXX.{npz,png}
B_matrix/caseNN_.../snapshots/...
B_matrix/caseNN_.../lines/ux_lines_stepXXXXXX.png
summary.csv
"""
from __future__ import annotations
import argparse
import csv
import json
import os
import sys
import tempfile
from dataclasses import dataclass
from typing import Any, Dict, Iterable, List, Optional, Sequence, Tuple
import numpy as np
_REPO = os.path.abspath(os.path.join(os.path.dirname(__file__), ".."))
_DEFAULT_LBM = os.path.join(_REPO, "src", "CelerisLab", "configs", "config_lbm.json")
# Default snapshot steps for parts A and B.
DEFAULT_SNAPSHOT_STEPS: Tuple[int, ...] = (100, 500, 1000, 1500, 2000)
@dataclass(frozen=True)
class CaseSpec:
"""One empty-channel configuration."""
case_id: str
label: str
inlet_scheme: str
y_wall_bc: str
outlet_mode: str
collision: str = "MRT"
def _load_json(path: str) -> dict:
with open(path, "r", encoding="utf-8") as f:
return json.load(f)
def _write_json(path: str, payload: dict) -> None:
os.makedirs(os.path.dirname(path) or ".", exist_ok=True)
with open(path, "w", encoding="utf-8") as f:
json.dump(payload, f, indent=2)
def vorticity_z(ux: np.ndarray, uy: np.ndarray) -> np.ndarray:
"""ωz = ∂uy/∂x ∂ux/∂y on (ny, nx) arrays."""
ux = np.asarray(ux, dtype=np.float64)
uy = np.asarray(uy, dtype=np.float64)
return np.gradient(uy, axis=1) - np.gradient(ux, axis=0)
def _line_y_indices(ny: int) -> List[Tuple[int, str]]:
return [
(1, "y1"),
(ny // 2, f"y{ny // 2}"),
(ny - 2, f"y{ny - 2}"),
]
def _build_cfg(
base: dict,
*,
nx: int,
ny: int,
collision: str,
inlet_scheme: str,
inlet_profile: str,
y_wall_bc: str,
outlet_mode: str,
velocity: float,
viscosity: float,
) -> dict:
cfg = json.loads(json.dumps(base))
cfg["grid"]["nx"] = int(nx)
cfg["grid"]["ny"] = int(ny)
cfg["grid"]["nz"] = 1
cfg["physics"]["velocity"] = float(velocity)
cfg["physics"]["viscosity"] = float(viscosity)
cfg["physics"]["rho"] = 1.0
cfg["method"]["collision"] = str(collision).upper()
cfg["method"]["streaming"] = "double_buffer"
cfg["method"]["store_precision"] = "FP32"
cfg["method"]["les"]["enabled"] = False
cfg["method"]["inlet"]["profile"] = str(inlet_profile)
cfg["method"]["inlet"]["scheme"] = str(inlet_scheme)
cfg["method"]["y_wall_bc"] = str(y_wall_bc)
cfg["method"]["outlet"]["mode"] = str(outlet_mode)
return cfg
def _field_stats(rho: np.ndarray, ux: np.ndarray, vort: np.ndarray) -> Dict[str, float]:
def _f(a: np.ndarray) -> float:
b = a[np.isfinite(a)]
return float("nan") if b.size == 0 else float(np.max(np.abs(b)))
return {
"rho_min": float(np.nanmin(rho)) if np.isfinite(rho).any() else float("nan"),
"rho_max": float(np.nanmax(rho)) if np.isfinite(rho).any() else float("nan"),
"ux_max": _f(ux),
"vort_max": _f(vort),
"finite": bool(np.isfinite(rho).all() and np.isfinite(ux).all()),
}
def _save_snapshot_npz(
path: str,
*,
step: int,
rho: np.ndarray,
ux: np.ndarray,
uy: np.ndarray,
vort: np.ndarray,
meta: dict,
) -> None:
os.makedirs(os.path.dirname(path), exist_ok=True)
np.savez_compressed(
path,
rho=np.asarray(rho, dtype=np.float32),
ux=np.asarray(ux, dtype=np.float32),
uy=np.asarray(uy, dtype=np.float32),
vort=np.asarray(vort, dtype=np.float32),
step=np.int32(step),
meta_json=np.asarray(json.dumps(meta)),
)
def _save_field_pngs(
out_dir: str,
prefix: str,
*,
rho: np.ndarray,
ux: np.ndarray,
vort: np.ndarray,
title: str,
) -> List[str]:
try:
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
except ImportError:
return []
os.makedirs(out_dir, exist_ok=True)
ny, nx = rho.shape
extent = (0, nx - 1, 0, ny - 1)
paths: List[str] = []
def _one(arr: np.ndarray, name: str, cmap: str, sym: bool) -> None:
a = np.asarray(arr, dtype=np.float64)
fin = a[np.isfinite(a)]
if fin.size == 0:
vmin, vmax = -1.0, 1.0
elif sym:
v = float(np.percentile(np.abs(fin), 99.5)) or 1.0
vmin, vmax = -v, v
else:
vmin = float(np.percentile(fin, 0.5))
vmax = float(np.percentile(fin, 99.5))
if vmax <= vmin:
vmax = vmin + 1.0
fig, ax = plt.subplots(figsize=(min(16.0, max(8.0, nx / 80.0)), min(8.0, max(3.0, ny / 50.0))))
im = ax.imshow(a, origin="lower", aspect="auto", cmap=cmap, vmin=vmin, vmax=vmax, extent=extent)
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_title(f"{title}{name}")
fig.colorbar(im, ax=ax, fraction=0.046, pad=0.04)
fig.tight_layout()
p = os.path.join(out_dir, f"{prefix}_{name}.png")
fig.savefig(p, dpi=140, bbox_inches="tight")
plt.close(fig)
paths.append(p)
_one(rho, "rho", "viridis", sym=False)
_one(ux, "ux", "RdBu_r", sym=True)
_one(vort, "vort", "RdBu_r", sym=True)
return paths
def _save_line_plots(
path: str,
*,
rho: np.ndarray,
ux: np.ndarray,
step: int,
case_label: str,
y_rows: Sequence[Tuple[int, str]],
) -> None:
try:
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
except ImportError:
return
ny, nx = rho.shape
x = np.arange(nx, dtype=np.float64)
fig, axes = plt.subplots(2, 1, figsize=(min(14.0, max(8.0, nx / 60.0)), 7.0), sharex=True)
for y_idx, y_lab in y_rows:
y_idx = int(np.clip(y_idx, 0, ny - 1))
axes[0].plot(x, ux[y_idx, :], label=y_lab, linewidth=1.2)
axes[1].plot(x, rho[y_idx, :], label=y_lab, linewidth=1.2)
axes[0].set_ylabel("u_x")
axes[0].legend(loc="best", fontsize=8)
axes[0].grid(True, alpha=0.3)
axes[1].set_ylabel("rho")
axes[1].set_xlabel("x (lattice)")
axes[1].legend(loc="best", fontsize=8)
axes[1].grid(True, alpha=0.3)
fig.suptitle(f"{case_label} — line profiles at step {step}")
fig.tight_layout()
os.makedirs(os.path.dirname(path) or ".", exist_ok=True)
fig.savefig(path, dpi=150, bbox_inches="tight")
plt.close(fig)
def _run_channel(
case: CaseSpec,
*,
base_cfg: dict,
nx: int,
ny: int,
velocity: float,
viscosity: float,
out_root: str,
snapshot_steps: Sequence[int],
max_step: int,
save_png: bool,
save_lines: bool,
stop_on_nan: bool,
) -> List[Dict[str, Any]]:
"""Run one case; write snapshots and optional line plots. Return summary rows."""
sys.path.insert(0, os.path.join(_REPO, "src"))
from CelerisLab import Simulation # noqa: WPS433
cfg = _build_cfg(
base_cfg,
nx=nx,
ny=ny,
collision=case.collision,
inlet_scheme=case.inlet_scheme,
inlet_profile="uniform",
y_wall_bc=case.y_wall_bc,
outlet_mode=case.outlet_mode,
velocity=velocity,
viscosity=viscosity,
)
bdoc = {"objects": []}
run_dir = os.path.join(out_root, case.case_id)
snap_dir = os.path.join(run_dir, "snapshots")
line_dir = os.path.join(run_dir, "lines")
os.makedirs(snap_dir, exist_ok=True)
if save_lines:
os.makedirs(line_dir, exist_ok=True)
tmpd = tempfile.mkdtemp(prefix="inlet_diag_")
lbm_tmp = os.path.join(tmpd, "config_lbm.json")
body_tmp = os.path.join(tmpd, "config_body.json")
_write_json(lbm_tmp, cfg)
_write_json(body_tmp, bdoc)
meta_base = {
"case_id": case.case_id,
"label": case.label,
"nx": nx,
"ny": ny,
"inlet_scheme": case.inlet_scheme,
"inlet_profile": "uniform",
"y_wall_bc": case.y_wall_bc,
"outlet_mode": case.outlet_mode,
"collision": case.collision,
"velocity": velocity,
"viscosity": viscosity,
}
_write_json(os.path.join(run_dir, "case_meta.json"), meta_base)
sim = Simulation(lbm_config_path=lbm_tmp, body_config_path=body_tmp)
sim.initialize()
y_rows = _line_y_indices(ny)
want_steps = sorted({int(s) for s in snapshot_steps if 0 < int(s) <= max_step})
next_snap = 0
rows: List[Dict[str, Any]] = []
for step in range(1, max_step + 1):
sim.step(1)
if next_snap < len(want_steps) and step == want_steps[next_snap]:
macro = sim.get_macroscopic()
rho = np.asarray(macro["rho"], dtype=np.float64)
ux = np.asarray(macro["ux"], dtype=np.float64)
uy = np.asarray(macro["uy"], dtype=np.float64)
vort = vorticity_z(ux, uy)
stats = _field_stats(rho, ux, vort)
stem = f"step_{step:06d}"
meta = {**meta_base, "step": step, **stats}
npz_path = os.path.join(snap_dir, f"{stem}.npz")
_save_snapshot_npz(
npz_path,
step=step,
rho=rho,
ux=ux,
uy=uy,
vort=vort,
meta=meta,
)
if save_png:
_save_field_pngs(
snap_dir,
stem,
rho=rho,
ux=ux,
vort=vort,
title=f"{case.label} step {step}",
)
if save_lines:
line_png = os.path.join(line_dir, f"lines_{stem}.png")
_save_line_plots(
line_png,
rho=rho,
ux=ux,
step=step,
case_label=case.label,
y_rows=y_rows,
)
# Also save raw 1D data for replotting.
line_npz = os.path.join(line_dir, f"lines_{stem}.npz")
payload = {"x": np.arange(nx, dtype=np.float32)}
for y_idx, y_lab in y_rows:
payload[f"ux_{y_lab}"] = ux[y_idx, :].astype(np.float32)
payload[f"rho_{y_lab}"] = rho[y_idx, :].astype(np.float32)
payload["step"] = np.int32(step)
np.savez_compressed(line_npz, **payload)
rows.append(
{
"case_id": case.case_id,
"label": case.label,
"collision": case.collision,
"inlet_scheme": case.inlet_scheme,
"y_wall_bc": case.y_wall_bc,
"outlet_mode": case.outlet_mode,
"step": step,
**stats,
"npz": npz_path,
}
)
if stop_on_nan and not stats["finite"]:
sim.close()
rows[-1]["stopped_early"] = True
return rows
next_snap += 1
sim.close()
return rows
def _part_a_cases() -> List[CaseSpec]:
# Kan99b-style baseline: zou_he + free_slip + neq_extrap; SRT and MRT.
base = CaseSpec(
case_id="",
label="",
inlet_scheme="zou_he_local",
y_wall_bc="free_slip",
outlet_mode="neq_extrap",
)
out: List[CaseSpec] = []
for coll in ("SRT", "MRT"):
cid = f"A_{coll.lower()}_zouhe_fs_neq"
out.append(
CaseSpec(
case_id=cid,
label=f"A baseline {coll} zou_he / free_slip / neq_extrap",
inlet_scheme=base.inlet_scheme,
y_wall_bc=base.y_wall_bc,
outlet_mode=base.outlet_mode,
collision=coll,
)
)
return out
def _part_b_cases() -> List[CaseSpec]:
return [
CaseSpec(
case_id="B_case01_zouhe_fs_neq",
label="1 zou_he / free_slip / neq_extrap",
inlet_scheme="zou_he_local",
y_wall_bc="free_slip",
outlet_mode="neq_extrap",
collision="MRT",
),
CaseSpec(
case_id="B_case02_zouhe_bb_neq",
label="2 zou_he / bounce_back / neq_extrap",
inlet_scheme="zou_he_local",
y_wall_bc="bounce_back",
outlet_mode="neq_extrap",
collision="MRT",
),
CaseSpec(
case_id="B_case03_zouhe_fs_zgrad",
label="3 zou_he / free_slip / zero_gradient",
inlet_scheme="zou_he_local",
y_wall_bc="free_slip",
outlet_mode="zero_gradient",
collision="MRT",
),
CaseSpec(
case_id="B_case04_stab_fs_neq",
label="4 channel_stabilized / free_slip / neq_extrap",
inlet_scheme="channel_stabilized",
y_wall_bc="free_slip",
outlet_mode="neq_extrap",
collision="MRT",
),
]
def _write_summary_csv(path: str, rows: Sequence[Dict[str, Any]]) -> None:
if not rows:
return
keys = [
"case_id",
"label",
"collision",
"inlet_scheme",
"y_wall_bc",
"outlet_mode",
"step",
"finite",
"rho_min",
"rho_max",
"ux_max",
"vort_max",
"stopped_early",
]
os.makedirs(os.path.dirname(path) or ".", exist_ok=True)
with open(path, "w", encoding="utf-8", newline="") as f:
w = csv.DictWriter(f, fieldnames=keys, extrasaction="ignore")
w.writeheader()
for r in rows:
w.writerow(r)
def main() -> int:
ap = argparse.ArgumentParser(description="Empty-channel inlet diagnostic (fields + line profiles)")
ap.add_argument(
"--part",
choices=("a", "b", "all"),
default="all",
help="A=baseline SRT/MRT; B=4-case matrix; all=both",
)
ap.add_argument("--nx", type=int, default=401, help="Channel length (lattice)")
ap.add_argument("--ny", type=int, default=201, help="Channel height (lattice)")
ap.add_argument("--velocity", type=float, default=0.03, help="Uniform inlet U0")
ap.add_argument("--viscosity", type=float, default=0.009, help="Kinematic viscosity")
ap.add_argument(
"--out-dir",
type=str,
default=os.path.join(_REPO, "tests", "output", "inlet_channel_diag"),
)
ap.add_argument(
"--steps",
type=str,
default="",
help="Comma-separated snapshot steps (default: 100,500,1000,1500,2000)",
)
ap.add_argument("--no-png", action="store_true", help="Skip rho/ux/vort PNG maps")
ap.add_argument("--no-lines", action="store_true", help="Skip ux/rho line-profile plots")
ap.add_argument(
"--continue-on-nan",
action="store_true",
help="Keep stepping after non-finite fields (default: stop case early)",
)
args = ap.parse_args()
if not os.path.isfile(_DEFAULT_LBM):
print(f"Missing config: {_DEFAULT_LBM}", file=sys.stderr)
return 2
if args.steps.strip():
snap_steps = tuple(int(s.strip()) for s in args.steps.split(",") if s.strip())
else:
snap_steps = DEFAULT_SNAPSHOT_STEPS
max_step = max(snap_steps)
base_cfg = _load_json(_DEFAULT_LBM)
out_dir = os.path.abspath(args.out_dir)
os.makedirs(out_dir, exist_ok=True)
cases: List[CaseSpec] = []
if args.part in ("a", "all"):
cases.extend(_part_a_cases())
if args.part in ("b", "all"):
cases.extend(_part_b_cases())
save_png = not args.no_png
# Part A: field maps only; Part B: fields + line plots.
all_rows: List[Dict[str, Any]] = []
for case in cases:
part = "A_baseline" if case.case_id.startswith("A_") else "B_matrix"
root = os.path.join(out_dir, part)
save_lines = not args.no_lines and part == "B_matrix"
print(f"--- {case.case_id}: {case.label} ({case.collision}) ---", flush=True)
rows = _run_channel(
case,
base_cfg=base_cfg,
nx=args.nx,
ny=args.ny,
velocity=args.velocity,
viscosity=args.viscosity,
out_root=root,
snapshot_steps=snap_steps,
max_step=max_step,
save_png=save_png,
save_lines=save_lines,
stop_on_nan=not args.continue_on_nan,
)
all_rows.extend(rows)
for r in rows:
fin = "OK" if r.get("finite") else "NONFINITE"
print(
f" step {r['step']:5d} {fin} rho=[{r['rho_min']:.4f},{r['rho_max']:.4f}] "
f"ux_max={r['ux_max']:.4f} vort_max={r['vort_max']:.4f}",
flush=True,
)
if rows and rows[-1].get("stopped_early"):
print(" (stopped early: non-finite fields)", flush=True)
summary_path = os.path.join(out_dir, "summary.csv")
_write_summary_csv(summary_path, all_rows)
manifest = {
"snapshot_steps": list(snap_steps),
"nx": args.nx,
"ny": args.ny,
"velocity": args.velocity,
"viscosity": args.viscosity,
"line_y_indices": [{"y": y, "label": lab} for y, lab in _line_y_indices(args.ny)],
"cases": [c.case_id for c in cases],
}
_write_json(os.path.join(out_dir, "manifest.json"), manifest)
print(f"Wrote: {summary_path}", flush=True)
print(f"Wrote: {os.path.join(out_dir, 'manifest.json')}", flush=True)
print(f"Output root: {out_dir}", flush=True)
return 0
if __name__ == "__main__":
raise SystemExit(main())

View File

@ -1,460 +0,0 @@
# CelerisLab/tests/run_inlet_ghost_timing_experiment.py
"""Minimal ghost-inlet timing experiments (post field-diagnostic).
Experiment 1 DDF time series at inlet center vs first interior column:
(x=0, y=NY/2) and (x=1, y=NY/2), steps 1..N.
Populations f1,f2,f5,f6,f7,f8 plus rho, ux (macro and sum f).
Experiment 2 Same channel with ``inlet.collide=false`` vs ``true``:
When collide is on, inlet ghost nodes undergo collision after Zou-He BC.
Compare rho_max / ux_max vs step to test ghost-source timing hypothesis.
Usage::
conda run -n pycuda_3_10 python tests/run_inlet_ghost_timing_experiment.py
conda run -n pycuda_3_10 python tests/run_inlet_ghost_timing_experiment.py --exp 1 --steps 50
conda run -n pycuda_3_10 python tests/run_inlet_ghost_timing_experiment.py --exp 2 --steps 500
"""
from __future__ import annotations
import argparse
import csv
import json
import os
import sys
import tempfile
from typing import Any, Dict, List, Optional, Sequence, Tuple
import numpy as np
_REPO = os.path.abspath(os.path.join(os.path.dirname(__file__), ".."))
_DEFAULT_LBM = os.path.join(_REPO, "src", "CelerisLab", "configs", "config_lbm.json")
# D2Q9 indices logged (see zou_he_velocity.cuh).
POP_IDX: Tuple[int, ...] = (1, 2, 5, 6, 7, 8)
POP_NAMES: Tuple[str, ...] = tuple(f"f{i}" for i in POP_IDX)
# cx, cy for macroscopic ux, uy from local f (matches descriptors.cuh D2Q9).
_CX = np.array([0, 1, -1, 0, 0, 1, -1, 1, -1], dtype=np.float64)
_CY = np.array([0, 0, 0, 1, -1, 1, -1, -1, 1], dtype=np.float64)
def _load_json(path: str) -> dict:
with open(path, "r", encoding="utf-8") as f:
return json.load(f)
def _write_json(path: str, payload: dict) -> None:
os.makedirs(os.path.dirname(path) or ".", exist_ok=True)
with open(path, "w", encoding="utf-8") as f:
json.dump(payload, f, indent=2)
def _build_cfg(
base: dict,
*,
nx: int,
ny: int,
collision: str,
inlet_collide: bool,
velocity: float,
viscosity: float,
) -> dict:
cfg = json.loads(json.dumps(base))
cfg["grid"]["nx"] = int(nx)
cfg["grid"]["ny"] = int(ny)
cfg["grid"]["nz"] = 1
cfg["physics"]["velocity"] = float(velocity)
cfg["physics"]["viscosity"] = float(viscosity)
cfg["physics"]["rho"] = 1.0
cfg["method"]["collision"] = str(collision).upper()
cfg["method"]["streaming"] = "double_buffer"
cfg["method"]["store_precision"] = "FP32"
cfg["method"]["les"]["enabled"] = False
cfg["method"]["inlet"]["profile"] = "uniform"
cfg["method"]["inlet"]["scheme"] = "zou_he_local"
cfg["method"]["inlet"]["collide"] = bool(inlet_collide)
cfg["method"]["y_wall_bc"] = "free_slip"
cfg["method"]["outlet"]["mode"] = "neq_extrap"
return cfg
def _macro_from_f(f: np.ndarray) -> Tuple[float, float, float]:
f = np.asarray(f, dtype=np.float64)
rho = float(np.sum(f))
if abs(rho) < 1e-14:
return rho, 0.0, 0.0
ux = float(np.dot(f, _CX) / rho)
uy = float(np.dot(f, _CY) / rho)
return rho, ux, uy
def _sample_node(ddf_qnyx: np.ndarray, x: int, y: int) -> Dict[str, float]:
f = ddf_qnyx[:, y, x].astype(np.float64)
rho_m, ux_m, uy_m = _macro_from_f(f)
out: Dict[str, float] = {
"rho_sum": rho_m,
"ux_macro": ux_m,
"uy_macro": uy_m,
}
for i, name in zip(POP_IDX, POP_NAMES):
out[name] = float(f[i])
return out
def _run_steps(
cfg: dict,
*,
steps: int,
y_mid: int,
probe_x: Sequence[int] = (0, 1),
) -> Tuple[List[Dict[str, Any]], Dict[str, Any]]:
sys.path.insert(0, os.path.join(_REPO, "src"))
from CelerisLab import Simulation # noqa: WPS433
bdoc = {"objects": []}
tmpd = tempfile.mkdtemp(prefix="ghost_timing_")
lbm_tmp = os.path.join(tmpd, "config_lbm.json")
body_tmp = os.path.join(tmpd, "config_body.json")
with open(lbm_tmp, "w", encoding="utf-8") as f:
json.dump(cfg, f)
with open(body_tmp, "w", encoding="utf-8") as f:
json.dump(bdoc, f)
sim = Simulation(lbm_config_path=lbm_tmp, body_config_path=body_tmp)
sim.initialize()
nx, ny = sim.lbm_cfg.nx, sim.lbm_cfg.ny
y_mid = int(np.clip(y_mid, 1, ny - 2))
rows: List[Dict[str, Any]] = []
for step in range(1, int(steps) + 1):
sim.step(1)
sim.field.download_ddf()
farr = sim.field.ddf.reshape(sim.lbm_cfg.nq, ny, nx)
rec: Dict[str, Any] = {"step": step}
for x in probe_x:
tag = "inlet" if x == 0 else "interior"
s = _sample_node(farr, int(x), y_mid)
for k, v in s.items():
rec[f"{tag}_{k}"] = v
# Pull semantics at interior: f[2] is read from stored f[2] at x=0 (same link index).
rec["cross_f2_match"] = abs(rec["interior_f2"] - rec["inlet_f2"]) < 1e-5
rec["cross_f1_inlet_to_int_pull"] = float(
farr[1, y_mid, 1]
) # after step, what x=1 holds in f1
rec["delta_inlet_f1"] = (
float("nan") if step == 1 else rec["inlet_f1"] - rows[-1]["inlet_f1"]
)
rec["delta_inlet_ux"] = (
float("nan")
if step == 1
else rec["inlet_ux_macro"] - rows[-1]["inlet_ux_macro"]
)
macro = sim.get_macroscopic()
rho_f = np.asarray(macro["rho"], dtype=np.float64)
ux_f = np.asarray(macro["ux"], dtype=np.float64)
rec["domain_rho_max"] = float(np.nanmax(rho_f))
rec["domain_rho_min"] = float(np.nanmin(rho_f))
rec["domain_ux_max"] = float(np.nanmax(np.abs(ux_f)))
rec["finite"] = bool(np.isfinite(rho_f).all() and np.isfinite(ux_f).all())
rows.append(rec)
meta = {
"nx": nx,
"ny": ny,
"y_mid": y_mid,
"probe_x": list(probe_x),
"inlet_collide": bool(cfg["method"]["inlet"].get("collide", False)),
"collision": cfg["method"]["collision"],
"inlet_scheme": cfg["method"]["inlet"]["scheme"],
"y_wall_bc": cfg["method"]["y_wall_bc"],
"outlet_mode": cfg["method"]["outlet"]["mode"],
"velocity": cfg["physics"]["velocity"],
"viscosity": cfg["physics"]["viscosity"],
}
sim.close()
return rows, meta
def _write_csv(path: str, rows: Sequence[Dict[str, Any]]) -> None:
if not rows:
return
keys: List[str] = []
for r in rows:
for k in r:
if k not in keys:
keys.append(k)
os.makedirs(os.path.dirname(path) or ".", exist_ok=True)
with open(path, "w", encoding="utf-8", newline="") as f:
w = csv.DictWriter(f, fieldnames=keys)
w.writeheader()
w.writerows(rows)
def _plot_exp1(out_dir: str, rows: Sequence[Dict[str, Any]], y_mid: int) -> List[str]:
try:
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
except ImportError:
return []
steps = [int(r["step"]) for r in rows]
paths: List[str] = []
def _ts(key: str, label: str, ax, **kw):
ax.plot(steps, [r[key] for r in rows], label=label, **kw)
# Populations
fig, axes = plt.subplots(2, 1, figsize=(11, 7), sharex=True)
for name in POP_NAMES:
_ts(f"inlet_{name}", f"inlet {name}", axes[0], linewidth=1.0)
axes[0].set_ylabel("f at x=0")
axes[0].legend(ncol=3, fontsize=7, loc="best")
axes[0].grid(True, alpha=0.3)
for name in POP_NAMES:
_ts(f"interior_{name}", f"int {name}", axes[1], linewidth=1.0)
axes[1].set_ylabel("f at x=1")
axes[1].set_xlabel("step")
axes[1].legend(ncol=3, fontsize=7, loc="best")
axes[1].grid(True, alpha=0.3)
fig.suptitle(f"Exp1 populations (y={y_mid})")
fig.tight_layout()
p1 = os.path.join(out_dir, "exp1_populations.png")
fig.savefig(p1, dpi=150, bbox_inches="tight")
plt.close(fig)
paths.append(p1)
# rho / ux + step-to-step deltas
fig, axes = plt.subplots(3, 1, figsize=(11, 8), sharex=True)
_ts("inlet_ux_macro", "inlet ux", axes[0])
_ts("interior_ux_macro", "interior ux", axes[0])
axes[0].axhline(0.03, color="k", ls="--", lw=0.8, label="U0")
axes[0].set_ylabel("u_x")
axes[0].legend(fontsize=8)
axes[0].grid(True, alpha=0.3)
_ts("inlet_rho_sum", "inlet rho", axes[1])
_ts("interior_rho_sum", "interior rho", axes[1])
axes[1].set_ylabel("rho")
axes[1].grid(True, alpha=0.3)
_ts("delta_inlet_f1", "|Δf1| inlet", axes[2])
axes[2].set_ylabel("Δf1")
axes[2].set_xlabel("step")
axes[2].grid(True, alpha=0.3)
fig.suptitle(f"Exp1 macro / inlet f1 increment (y={y_mid})")
fig.tight_layout()
p2 = os.path.join(out_dir, "exp1_macro_delta.png")
fig.savefig(p2, dpi=150, bbox_inches="tight")
plt.close(fig)
paths.append(p2)
return paths
def _plot_exp2(out_dir: str, rows_a: Sequence[Dict[str, Any]], rows_b: Sequence[Dict[str, Any]]) -> List[str]:
try:
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
except ImportError:
return []
fig, axes = plt.subplots(2, 1, figsize=(10, 6), sharex=True)
for rows, lab, c in (
(rows_a, "ghost (no collide)", "C0"),
(rows_b, "inlet collide", "C1"),
):
steps = [int(r["step"]) for r in rows]
axes[0].plot(steps, [r["domain_rho_max"] for r in rows], label=lab, color=c)
axes[1].plot(steps, [r["domain_ux_max"] for r in rows], label=lab, color=c)
axes[0].set_ylabel("rho_max")
axes[0].legend()
axes[0].grid(True, alpha=0.3)
axes[1].set_ylabel("|ux|_max")
axes[1].set_xlabel("step")
axes[1].legend()
axes[1].grid(True, alpha=0.3)
fig.suptitle("Exp2: ghost vs inlet collide")
fig.tight_layout()
p = os.path.join(out_dir, "exp2_stability_compare.png")
fig.savefig(p, dpi=150, bbox_inches="tight")
plt.close(fig)
return [p]
def _oscillation_summary(rows: Sequence[Dict[str, Any]], *, last_n: int = 20) -> Dict[str, float]:
"""High-frequency proxy: std of inlet f1 and ux over the last *last_n* steps."""
if len(rows) < 2:
return {}
tail = rows[-min(last_n, len(rows)) :]
f1 = np.array([r["inlet_f1"] for r in tail], dtype=np.float64)
ux = np.array([r["inlet_ux_macro"] for r in tail], dtype=np.float64)
d1 = np.array([r["delta_inlet_f1"] for r in tail if np.isfinite(r["delta_inlet_f1"])], dtype=np.float64)
return {
"std_inlet_f1_last": float(np.std(f1)),
"std_inlet_ux_last": float(np.std(ux)),
"mean_abs_delta_f1_last": float(np.mean(np.abs(d1))) if d1.size else float("nan"),
}
def run_exp1(
base: dict,
*,
out_dir: str,
nx: int,
ny: int,
steps: int,
collision: str,
velocity: float,
viscosity: float,
) -> None:
y_mid = ny // 2
cfg = _build_cfg(
base,
nx=nx,
ny=ny,
collision=collision,
inlet_collide=False,
velocity=velocity,
viscosity=viscosity,
)
print(f"Exp1: zou_he ghost inlet, y_mid={y_mid}, steps={steps}", flush=True)
rows, meta = _run_steps(cfg, steps=steps, y_mid=y_mid)
meta["oscillation"] = _oscillation_summary(rows)
exp_dir = os.path.join(out_dir, "exp1_ddf_timeseries")
os.makedirs(exp_dir, exist_ok=True)
_write_csv(os.path.join(exp_dir, "timeseries.csv"), rows)
_write_json(os.path.join(exp_dir, "meta.json"), meta)
plots = _plot_exp1(exp_dir, rows, y_mid)
for p in plots:
print(f" plot: {p}", flush=True)
# Console summary for quick read
print(" last 5 steps (inlet center):", flush=True)
for r in rows[-5:]:
print(
f" step {r['step']:3d} f1={r['inlet_f1']:.6f} ux={r['inlet_ux_macro']:.6f} "
f"Δf1={r['delta_inlet_f1']:.2e} rho_max={r['domain_rho_max']:.4f} finite={r['finite']}",
flush=True,
)
print(f" oscillation: {meta['oscillation']}", flush=True)
print(f"Wrote: {exp_dir}/timeseries.csv", flush=True)
def run_exp2(
base: dict,
*,
out_dir: str,
nx: int,
ny: int,
steps: int,
collision: str,
velocity: float,
viscosity: float,
) -> None:
y_mid = ny // 2
exp_dir = os.path.join(out_dir, "exp2_inlet_collide")
os.makedirs(exp_dir, exist_ok=True)
summaries: Dict[str, Any] = {}
all_rows: Dict[str, List[Dict[str, Any]]] = {}
for collide, tag in ((False, "ghost_no_collide"), (True, "ghost_with_collide")):
cfg = _build_cfg(
base,
nx=nx,
ny=ny,
collision=collision,
inlet_collide=collide,
velocity=velocity,
viscosity=viscosity,
)
print(f"Exp2 [{tag}]: inlet.collide={collide}, steps={steps}", flush=True)
rows, meta = _run_steps(cfg, steps=steps, y_mid=y_mid)
all_rows[tag] = rows
_write_csv(os.path.join(exp_dir, f"{tag}.csv"), rows)
last_finite = next(
(int(r["step"]) for r in rows if not r.get("finite", True)),
None,
)
summaries[tag] = {
**meta,
"first_nonfinite_step": last_finite,
"final_rho_max": rows[-1]["domain_rho_max"] if rows else None,
"final_finite": rows[-1].get("finite") if rows else None,
"oscillation": _oscillation_summary(rows),
}
print(
f" final rho_max={summaries[tag]['final_rho_max']:.4f} "
f"finite={summaries[tag]['final_finite']} "
f"first_nonfinite={summaries[tag]['first_nonfinite_step']}",
flush=True,
)
_write_json(os.path.join(exp_dir, "summary.json"), summaries)
plots = _plot_exp2(exp_dir, all_rows["ghost_no_collide"], all_rows["ghost_with_collide"])
for p in plots:
print(f" plot: {p}", flush=True)
print(f"Wrote: {exp_dir}/summary.json", flush=True)
def main() -> int:
ap = argparse.ArgumentParser(description="Ghost inlet timing experiments")
ap.add_argument("--exp", choices=("1", "2", "all"), default="all")
ap.add_argument("--steps", type=int, default=50, help="Steps for exp1 (default 50)")
ap.add_argument("--steps2", type=int, default=500, help="Steps for exp2 (default 500)")
ap.add_argument("--nx", type=int, default=401)
ap.add_argument("--ny", type=int, default=201)
ap.add_argument("--collision", default="MRT", choices=("SRT", "TRT", "MRT"))
ap.add_argument("--velocity", type=float, default=0.03)
ap.add_argument("--viscosity", type=float, default=0.009)
ap.add_argument(
"--out-dir",
default=os.path.join(_REPO, "tests", "output", "inlet_ghost_timing"),
)
args = ap.parse_args()
if not os.path.isfile(_DEFAULT_LBM):
print(f"Missing {_DEFAULT_LBM}", file=sys.stderr)
return 2
base = _load_json(_DEFAULT_LBM)
out_dir = os.path.abspath(args.out_dir)
os.makedirs(out_dir, exist_ok=True)
if args.exp in ("1", "all"):
run_exp1(
base,
out_dir=out_dir,
nx=args.nx,
ny=args.ny,
steps=args.steps,
collision=args.collision,
velocity=args.velocity,
viscosity=args.viscosity,
)
if args.exp in ("2", "all"):
run_exp2(
base,
out_dir=out_dir,
nx=args.nx,
ny=args.ny,
steps=args.steps2,
collision=args.collision,
velocity=args.velocity,
viscosity=args.viscosity,
)
return 0
if __name__ == "__main__":
raise SystemExit(main())

View File

@ -1,587 +0,0 @@
# CelerisLab/tests/run_inlet_scenario_fields.py
"""Three-scenario inlet field export: final or last-stable step at 5000 LBM steps.
Scenarios:
- empty_channel: zou_he_local × {SRT,MRT} × {free_slip,bounce_back}
- empty_channel_inlet_matrix: all inlet schemes × {SRT,MRT}, bounce_back only
- kan99b: zou_he_local × {SRT,MRT} × {free_slip,bounce_back}, Re=100, domain M
- sah04_case9: channel_stabilized × {SRT,MRT}, high-blockage case 9 geometry
Outputs per run (under ``--out-dir``):
fields/final_{rho,ux,vort}.png, fields/final.npz
lines/lines_ux_rho.png, lines/lines.npz
run_meta.json
Usage::
conda run -n pycuda_3_10 python tests/run_inlet_scenario_fields.py
conda run -n pycuda_3_10 python tests/run_inlet_scenario_fields.py --scenario empty_channel
conda run -n pycuda_3_10 python tests/run_inlet_scenario_fields.py --scenario empty_channel_inlet_matrix
"""
from __future__ import annotations
import argparse
import csv
import json
import os
import sys
import tempfile
from dataclasses import dataclass, replace
from typing import Any, Dict, List, Optional, Sequence, Tuple
import numpy as np
_REPO = os.path.abspath(os.path.join(os.path.dirname(__file__), ".."))
_DEFAULT_LBM = os.path.join(_REPO, "src", "CelerisLab", "configs", "config_lbm.json")
# Kan99b lattice contract
_KAN_U_INF = 0.03
_KAN_D = 30.0
_KAN_R = 15.0
_KAN_RE = 100.0
_KAN_ALPHA = 1.0
# Sah04 case 9 (high tier)
_SAH_D = 30
_SAH_NX = 80 * _SAH_D + 2
_SAH_NY = 35
_SAH_CX = 40.0 * _SAH_D + 0.5
_SAH_CY = 17.0
_SAH_RE = 200.0
_SAH_U_MAX = 0.1
@dataclass(frozen=True)
class RunSpec:
"""One simulation run specification."""
scenario: str
run_id: str
label: str
nx: int
ny: int
collision: str
inlet_scheme: str
inlet_profile: str
y_wall_bc: str
outlet_mode: str
velocity: float
viscosity: float
steps: int
has_cylinder: bool
cylinder_center: Tuple[float, float] = (0.0, 0.0)
cylinder_radius: float = 0.0
cylinder_omega: float = 0.0
def _load_json(path: str) -> dict:
with open(path, "r", encoding="utf-8") as f:
return json.load(f)
def _write_json(path: str, payload: dict) -> None:
os.makedirs(os.path.dirname(path) or ".", exist_ok=True)
with open(path, "w", encoding="utf-8") as f:
json.dump(payload, f, indent=2)
def vorticity_z(ux: np.ndarray, uy: np.ndarray) -> np.ndarray:
ux = np.asarray(ux, dtype=np.float64)
uy = np.asarray(uy, dtype=np.float64)
return np.gradient(uy, axis=1) - np.gradient(ux, axis=0)
def _line_y_indices(ny: int) -> List[Tuple[int, str]]:
return [(1, "y1"), (ny // 2, f"y{ny // 2}"), (ny - 2, f"y{ny - 2}")]
_INLET_SCHEMES = (
"zou_he_local",
"channel_stabilized",
"equilibrium",
"regularized",
)
def _empty_channel_inlet_matrix_specs() -> List[RunSpec]:
"""All inlet schemes on empty channel, bounce_back, SRT/MRT (5000-step field export)."""
specs: List[RunSpec] = []
for coll in ("SRT", "MRT"):
for scheme in _INLET_SCHEMES:
run_id = f"{coll.lower()}_{scheme}"
specs.append(
RunSpec(
scenario="empty_channel_inlet_matrix",
run_id=run_id,
label=f"empty {scheme} {coll} bounce_back",
nx=401,
ny=201,
collision=coll,
inlet_scheme=scheme,
inlet_profile="uniform",
y_wall_bc="bounce_back",
outlet_mode="neq_extrap",
velocity=0.03,
viscosity=0.009,
steps=5000,
has_cylinder=False,
)
)
return specs
def _all_specs() -> List[RunSpec]:
specs: List[RunSpec] = []
for coll in ("SRT", "MRT"):
for wall in ("free_slip", "bounce_back"):
wid = f"{coll.lower()}_{wall}"
specs.append(
RunSpec(
scenario="empty_channel",
run_id=wid,
label=f"empty zou_he {coll} {wall}",
nx=401,
ny=201,
collision=coll,
inlet_scheme="zou_he_local",
inlet_profile="uniform",
y_wall_bc=wall,
outlet_mode="neq_extrap",
velocity=0.03,
viscosity=0.009,
steps=5000,
has_cylinder=False,
)
)
dom_m = (1351, 601, (450.0, 300.0))
nu_k = _KAN_U_INF * _KAN_D / _KAN_RE
omega = 2.0 * _KAN_ALPHA * _KAN_U_INF / _KAN_D
for coll in ("SRT", "MRT"):
for wall in ("free_slip", "bounce_back"):
wid = f"{coll.lower()}_{wall}"
specs.append(
RunSpec(
scenario="kan99b",
run_id=wid,
label=f"kan99b zou_he {coll} {wall}",
nx=dom_m[0],
ny=dom_m[1],
collision=coll,
inlet_scheme="zou_he_local",
inlet_profile="uniform",
y_wall_bc=wall,
outlet_mode="neq_extrap",
velocity=_KAN_U_INF,
viscosity=nu_k,
steps=5000,
has_cylinder=True,
cylinder_center=dom_m[2],
cylinder_radius=_KAN_R,
cylinder_omega=omega,
)
)
u0_mean = _SAH_U_MAX / 1.5
nu_s = _SAH_U_MAX * _SAH_D / _SAH_RE
for coll in ("SRT", "MRT"):
wid = f"{coll.lower()}_channel_stab"
specs.append(
RunSpec(
scenario="sah04_case9",
run_id=wid,
label=f"sah04 case9 channel_stab {coll}",
nx=_SAH_NX,
ny=_SAH_NY,
collision=coll,
inlet_scheme="channel_stabilized",
inlet_profile="parabolic",
y_wall_bc="bounce_back",
outlet_mode="neq_extrap",
velocity=u0_mean,
viscosity=nu_s,
steps=5000,
has_cylinder=True,
cylinder_center=(_SAH_CX, _SAH_CY),
cylinder_radius=0.5 * _SAH_D,
cylinder_omega=0.0,
)
)
return specs
def _build_cfg(base: dict, spec: RunSpec) -> dict:
cfg = json.loads(json.dumps(base))
cfg["grid"]["nx"] = spec.nx
cfg["grid"]["ny"] = spec.ny
cfg["grid"]["nz"] = 1
cfg["physics"]["velocity"] = float(spec.velocity)
cfg["physics"]["viscosity"] = float(spec.viscosity)
cfg["physics"]["rho"] = 1.0
cfg["method"]["collision"] = spec.collision.upper()
cfg["method"]["streaming"] = "double_buffer"
cfg["method"]["store_precision"] = "FP32"
cfg["method"]["les"]["enabled"] = False
cfg["method"]["inlet"]["profile"] = spec.inlet_profile
cfg["method"]["inlet"]["scheme"] = spec.inlet_scheme
cfg["method"]["y_wall_bc"] = spec.y_wall_bc
cfg["method"]["outlet"]["mode"] = spec.outlet_mode
return cfg
def _body_doc(spec: RunSpec) -> dict:
if not spec.has_cylinder:
return {"objects": []}
return {
"objects": [
{
"type": "cylinder",
"center": [float(spec.cylinder_center[0]), float(spec.cylinder_center[1])],
"radius": float(spec.cylinder_radius),
"omega": float(spec.cylinder_omega),
}
]
}
def _save_field_pngs(
out_dir: str,
prefix: str,
*,
rho: np.ndarray,
ux: np.ndarray,
vort: np.ndarray,
title: str,
) -> List[str]:
try:
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
except ImportError:
return []
os.makedirs(out_dir, exist_ok=True)
ny, nx = rho.shape
extent = (0, nx - 1, 0, ny - 1)
paths: List[str] = []
def _one(arr: np.ndarray, name: str, cmap: str, sym: bool) -> None:
a = np.asarray(arr, dtype=np.float64)
fin = a[np.isfinite(a)]
if fin.size == 0:
vmin, vmax = -1.0, 1.0
elif sym:
v = float(np.percentile(np.abs(fin), 99.5)) or 1.0
vmin, vmax = -v, v
else:
vmin = float(np.percentile(fin, 0.5))
vmax = float(np.percentile(fin, 99.5))
if vmax <= vmin:
vmax = vmin + 1.0
fw = min(18.0, max(8.0, nx / 70.0))
fh = min(10.0, max(3.0, ny / 45.0))
fig, ax = plt.subplots(figsize=(fw, fh))
im = ax.imshow(a, origin="lower", aspect="auto", cmap=cmap, vmin=vmin, vmax=vmax, extent=extent)
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_title(f"{title}{name}")
fig.colorbar(im, ax=ax, fraction=0.046, pad=0.04)
fig.tight_layout()
p = os.path.join(out_dir, f"{prefix}_{name}.png")
fig.savefig(p, dpi=150, bbox_inches="tight")
plt.close(fig)
paths.append(p)
_one(rho, "rho", "viridis", sym=False)
_one(ux, "ux", "RdBu_r", sym=True)
_one(vort, "vort", "RdBu_r", sym=True)
return paths
def _save_line_plots(
path: str,
*,
rho: np.ndarray,
ux: np.ndarray,
step: int,
label: str,
y_rows: Sequence[Tuple[int, str]],
) -> None:
try:
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
except ImportError:
return
ny, nx = rho.shape
x = np.arange(nx, dtype=np.float64)
fig, axes = plt.subplots(2, 1, figsize=(min(14.0, max(8.0, nx / 55.0)), 7.0), sharex=True)
for y_idx, y_lab in y_rows:
yi = int(np.clip(y_idx, 0, ny - 1))
axes[0].plot(x, ux[yi, :], label=y_lab, linewidth=1.0)
axes[1].plot(x, rho[yi, :], label=y_lab, linewidth=1.0)
axes[0].set_ylabel("u_x")
axes[0].legend(loc="best", fontsize=8)
axes[0].grid(True, alpha=0.3)
axes[1].set_ylabel("rho")
axes[1].set_xlabel("x (lattice)")
axes[1].legend(loc="best", fontsize=8)
axes[1].grid(True, alpha=0.3)
fig.suptitle(f"{label} — ux/rho lines at step {step}")
fig.tight_layout()
os.makedirs(os.path.dirname(path) or ".", exist_ok=True)
fig.savefig(path, dpi=150, bbox_inches="tight")
plt.close(fig)
def _snapshot_from_sim(sim) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
macro = sim.get_macroscopic()
rho = np.asarray(macro["rho"], dtype=np.float64)
ux = np.asarray(macro["ux"], dtype=np.float64)
uy = np.asarray(macro["uy"], dtype=np.float64)
vort = vorticity_z(ux, uy)
return rho, ux, uy, vort
def _is_stable_fields(
rho: np.ndarray,
ux: np.ndarray,
*,
rho_lo: float = 0.85,
rho_hi: float = 1.25,
ux_cap: float = 0.15,
) -> bool:
"""Finite fields within a physically plausible band (reject pre-blow-up states)."""
if not (np.isfinite(rho).all() and np.isfinite(ux).all()):
return False
r0 = float(np.min(rho))
r1 = float(np.max(rho))
umax = float(np.max(np.abs(ux)))
return (rho_lo <= r0) and (r1 <= rho_hi) and (umax <= ux_cap)
def run_one(spec: RunSpec, base_cfg: dict, out_root: str) -> Dict[str, Any]:
sys.path.insert(0, os.path.join(_REPO, "src"))
import pycuda.driver as cuda
from CelerisLab import Simulation # noqa: WPS433
run_dir = os.path.join(out_root, spec.scenario, spec.run_id)
field_dir = os.path.join(run_dir, "fields")
line_dir = os.path.join(run_dir, "lines")
os.makedirs(field_dir, exist_ok=True)
os.makedirs(line_dir, exist_ok=True)
cfg = _build_cfg(base_cfg, spec)
tmpd = tempfile.mkdtemp(prefix="inlet_scenario_")
lbm_tmp = os.path.join(tmpd, "config_lbm.json")
body_tmp = os.path.join(tmpd, "config_body.json")
_write_json(lbm_tmp, cfg)
_write_json(body_tmp, _body_doc(spec))
sim = Simulation(lbm_config_path=lbm_tmp, body_config_path=body_tmp)
if spec.has_cylinder and spec.cylinder_omega != 0.0:
sim.bodies.get(0).state.omega = np.float32(spec.cylinder_omega)
sim.initialize()
stream = cuda.Stream()
y_rows = _line_y_indices(spec.ny)
last_good: Optional[Dict[str, Any]] = None
first_bad_step: Optional[int] = None
force_bad_step: Optional[int] = None
print(f" [{spec.scenario}/{spec.run_id}] {spec.label} steps={spec.steps}", flush=True)
for step in range(1, spec.steps + 1):
if spec.has_cylinder:
sim.bodies.zero_force_segment_async(stream)
sim.stepper.step(
1,
action_gpu=sim.bodies.action_gpu,
obs_gpu=sim.bodies.obs_gpu,
stream=stream,
)
if step % 100 == 0 or step == spec.steps:
stream.synchronize()
sim.bodies.download_obs_full_async(stream)
stream.synchronize()
fvec = sim.bodies.read_force(0)
if not (np.isfinite(fvec[0]) and np.isfinite(fvec[1])):
if force_bad_step is None:
force_bad_step = step
else:
sim.step(1)
rho, ux, uy, vort = _snapshot_from_sim(sim)
if _is_stable_fields(rho, ux):
last_good = {
"step": step,
"rho": rho.copy(),
"ux": ux.copy(),
"uy": uy.copy(),
"vort": vort.copy(),
}
elif first_bad_step is None:
first_bad_step = step
sim.close()
if last_good is None:
raise RuntimeError(f"No finite snapshot for {spec.run_id}")
out_step = int(last_good["step"])
rho = last_good["rho"]
ux = last_good["ux"]
uy = last_good["uy"]
vort = last_good["vort"]
requested_final = spec.steps
used_last_stable = out_step < requested_final
meta = {
"scenario": spec.scenario,
"run_id": spec.run_id,
"label": spec.label,
"nx": spec.nx,
"ny": spec.ny,
"collision": spec.collision,
"inlet_scheme": spec.inlet_scheme,
"inlet_profile": spec.inlet_profile,
"y_wall_bc": spec.y_wall_bc,
"outlet_mode": spec.outlet_mode,
"velocity": spec.velocity,
"viscosity": spec.viscosity,
"requested_steps": requested_final,
"output_step": out_step,
"used_last_stable": used_last_stable,
"first_nonfinite_step": first_bad_step,
"first_force_nonfinite_step": force_bad_step,
"rho_min": float(np.min(rho)),
"rho_max": float(np.max(rho)),
"ux_max": float(np.max(np.abs(ux))),
"vort_max": float(np.max(np.abs(vort[np.isfinite(vort)]))) if np.isfinite(vort).any() else float("nan"),
}
_write_json(os.path.join(run_dir, "run_meta.json"), meta)
stem = f"step_{out_step:06d}"
np.savez_compressed(
os.path.join(field_dir, "final.npz"),
rho=rho.astype(np.float32),
ux=ux.astype(np.float32),
uy=uy.astype(np.float32),
vort=vort.astype(np.float32),
step=np.int32(out_step),
)
title = f"{spec.label} (step {out_step}" + (", last stable" if used_last_stable else ", final") + ")"
pngs = _save_field_pngs(field_dir, "final", rho=rho, ux=ux, vort=vort, title=title)
_save_line_plots(
os.path.join(line_dir, "lines_ux_rho.png"),
rho=rho,
ux=ux,
step=out_step,
label=spec.label,
y_rows=y_rows,
)
line_payload: Dict[str, Any] = {"x": np.arange(spec.nx, dtype=np.float32), "step": np.int32(out_step)}
for y_idx, y_lab in y_rows:
yi = int(np.clip(y_idx, 0, spec.ny - 1))
line_payload[f"ux_{y_lab}"] = ux[yi, :].astype(np.float32)
line_payload[f"rho_{y_lab}"] = rho[yi, :].astype(np.float32)
np.savez_compressed(os.path.join(line_dir, "lines.npz"), **line_payload)
status = "last_stable" if used_last_stable else "final"
print(
f" -> {status} step {out_step} rho=[{meta['rho_min']:.4f},{meta['rho_max']:.4f}] "
f"ux_max={meta['ux_max']:.4f} force_bad={force_bad_step}",
flush=True,
)
return {**meta, "field_pngs": pngs, "run_dir": run_dir}
def main() -> int:
ap = argparse.ArgumentParser(description="Three-scenario inlet field export (5000 steps)")
ap.add_argument(
"--scenario",
choices=(
"empty_channel",
"empty_channel_inlet_matrix",
"kan99b",
"sah04_case9",
"all",
),
default="all",
)
ap.add_argument(
"--collision",
default="",
help="Optional filter: SRT or MRT only",
)
ap.add_argument("--steps", type=int, default=5000)
ap.add_argument(
"--out-dir",
default=os.path.join(_REPO, "tests", "output", "inlet_scenario_fields"),
)
args = ap.parse_args()
if not os.path.isfile(_DEFAULT_LBM):
print(f"Missing {_DEFAULT_LBM}", file=sys.stderr)
return 2
base = _load_json(_DEFAULT_LBM)
if args.scenario == "empty_channel_inlet_matrix":
specs = _empty_channel_inlet_matrix_specs()
elif args.scenario == "all":
specs = _all_specs() + _empty_channel_inlet_matrix_specs()
else:
specs = _all_specs()
if args.scenario != "all":
specs = [s for s in specs if s.scenario == args.scenario]
if args.collision.strip():
coll = args.collision.strip().upper()
specs = [s for s in specs if s.collision.upper() == coll]
specs = [replace(s, steps=int(args.steps)) for s in specs]
out_dir = os.path.abspath(args.out_dir)
os.makedirs(out_dir, exist_ok=True)
rows: List[Dict[str, Any]] = []
for spec in specs:
try:
row = run_one(spec, base, out_dir)
rows.append(row)
except Exception as e: # noqa: BLE001
print(f" FAILED {spec.scenario}/{spec.run_id}: {e}", flush=True)
rows.append(
{
"scenario": spec.scenario,
"run_id": spec.run_id,
"label": spec.label,
"error": str(e),
}
)
summary_path = os.path.join(out_dir, "summary.csv")
if rows:
keys: List[str] = []
for r in rows:
for k in r:
if k not in keys and k != "field_pngs":
keys.append(k)
with open(summary_path, "w", encoding="utf-8", newline="") as f:
w = csv.DictWriter(f, fieldnames=keys, extrasaction="ignore")
w.writeheader()
w.writerows(rows)
_write_json(
os.path.join(out_dir, "manifest.json"),
{"steps": args.steps, "scenario_filter": args.scenario, "runs": [s.run_id for s in specs]},
)
print(f"Wrote: {summary_path}", flush=True)
print(f"Output: {out_dir}", flush=True)
return 0
if __name__ == "__main__":
raise SystemExit(main())

View File

@ -1,30 +1,11 @@
# CelerisLab/tests/run_kan99b_rotating_cylinder.py
"""Kan99b rotating-cylinder validation driver.
"""Kan99b MRT-only rotating-cylinder validation runner.
This script executes the rotating-cylinder campaign in
``tests/Rotating_cylinder_validation_plan.md`` against Kan99b anchors.
This script follows ``tests/Kan99b_validation.md`` for the current round:
Core lattice mapping (fixed by campaign contract):
- D = 30, R = 15
- U_inf = 0.03
- nu = U_inf * D / Re = 0.9 / Re
- omega_body = 2 * alpha * U_inf / D = 0.002 * alpha
- inlet.profile = uniform
- y_wall_bc = free_slip
- outlet.mode = neq_extrap
- streaming = double_buffer
Phases:
- A: domain independence at Re=100, alpha=1.0 (MRT, domains S/M/L)
- B: anchor collision sweep at Re=100, alpha=1.0 (SRT/TRT/MRT)
- C: Re=100 alpha scan
- D: Re=60 and Re=160 threshold scan
Usage examples::
conda run -n pycuda_3_10 python tests/run_kan99b_rotating_cylinder.py --phase a
conda run -n pycuda_3_10 python tests/run_kan99b_rotating_cylinder.py --phase b --domain M
conda run -n pycuda_3_10 python tests/run_kan99b_rotating_cylinder.py --phase c --minimal
conda run -n pycuda_3_10 python tests/run_kan99b_rotating_cylinder.py --phase all --minimal
- Primary matrix: K1-K5 with collision fixed to MRT.
- Primary inlet: regularized (uniform profile).
- Extra control: K2 with ``zou_he_local`` inlet for sensitivity only.
"""
from __future__ import annotations
@ -36,7 +17,7 @@ import os
import sys
import tempfile
from dataclasses import dataclass
from typing import Any, Dict, Iterable, List, Optional, Sequence, Tuple
from typing import Any, Dict, List, Optional, Sequence, Tuple
import numpy as np
import pycuda.driver as cuda
@ -48,7 +29,7 @@ U_INF = 0.03
D_LATTICE = 30.0
R_LATTICE = 15.0
# Kan99b Table I anchor (Re=100, alpha=1.0).
KAN99B_ANCHOR = {
"St": 0.1655,
"mean_cl": -2.4881,
@ -57,7 +38,6 @@ KAN99B_ANCHOR = {
"amp_cd": 0.0993,
}
# Preferred agreement bands from the validation plan (fractional errors).
ANCHOR_BANDS = {
"St": 0.03,
"mean_cl": 0.04,
@ -66,20 +46,9 @@ ANCHOR_BANDS = {
"amp_cd": 0.10,
}
# Domain sensitivity thresholds vs domain L (fractional errors).
DOMAIN_THRESH = {
"St": 0.01,
"mean_cl": 0.02,
"mean_cd": 0.02,
"amp_cl": 0.03,
"amp_cd": 0.03,
}
@dataclass(frozen=True)
class DomainSpec:
"""Rectangular domain defined in lattice units."""
key: str
nx: int
ny: int
@ -87,18 +56,45 @@ class DomainSpec:
@dataclass(frozen=True)
class RunSpec:
"""One executable run specification."""
phase: str
collision: str
domain: str
class KanCase:
case_id: str
re: float
alpha: float
steps: int
burn: int
@dataclass(frozen=True)
class RunSpec:
case_id: str
variant: str
domain: str
re: float
alpha: float
inlet_scheme: str
steps: int
burn: int
CASES: Tuple[KanCase, ...] = (
KanCase("K1", 100.0, 0.5, 200_000, 80_000),
KanCase("K2", 100.0, 1.0, 200_000, 80_000),
KanCase("K3", 60.0, 1.6, 240_000, 120_000),
KanCase("K4", 100.0, 2.0, 240_000, 120_000),
KanCase("K5", 160.0, 2.0, 240_000, 120_000),
)
CASE_K0 = KanCase("K0", 100.0, 0.0, 180_000, 72_000)
def _domain_specs() -> Dict[str, DomainSpec]:
return {
"S": DomainSpec("S", 1081, 481, (360.0, 240.0)),
"M": DomainSpec("M", 1351, 601, (450.0, 300.0)),
"L": DomainSpec("L", 1801, 721, (600.0, 360.0)),
}
def _load_json(path: str) -> dict:
with open(path, "r", encoding="utf-8") as f:
return json.load(f)
@ -109,14 +105,6 @@ def _write_json(path: str, payload: dict) -> None:
json.dump(payload, f, indent=2)
def _domain_specs() -> Dict[str, DomainSpec]:
return {
"S": DomainSpec("S", 1081, 481, (360.0, 240.0)),
"M": DomainSpec("M", 1351, 601, (450.0, 300.0)),
"L": DomainSpec("L", 1801, 721, (600.0, 360.0)),
}
def _nu_from_re(re: float) -> float:
return U_INF * D_LATTICE / float(re)
@ -127,10 +115,20 @@ def _omega_body(alpha: float) -> float:
def _run_id(spec: RunSpec) -> str:
a = f"{spec.alpha:.3f}".replace(".", "p")
return f"phase{spec.phase}_dom{spec.domain}_re{int(spec.re)}_a{a}_{spec.collision.lower()}"
return (
f"{spec.case_id.lower()}_{spec.variant}_dom{spec.domain}_re{int(spec.re)}_a{a}_"
f"{spec.inlet_scheme.lower()}_mrt"
)
def _build_cfg(base_cfg: dict, *, nx: int, ny: int, collision: str, re: float) -> dict:
def _build_cfg(
base_cfg: dict,
*,
nx: int,
ny: int,
re: float,
inlet_scheme: str,
) -> dict:
cfg = json.loads(json.dumps(base_cfg))
cfg["grid"]["nx"] = int(nx)
cfg["grid"]["ny"] = int(ny)
@ -138,18 +136,19 @@ def _build_cfg(base_cfg: dict, *, nx: int, ny: int, collision: str, re: float) -
cfg["physics"]["velocity"] = float(U_INF)
cfg["physics"]["viscosity"] = float(_nu_from_re(re))
cfg["physics"]["rho"] = 1.0
cfg["method"]["collision"] = str(collision).upper()
cfg["method"]["collision"] = "MRT"
cfg["method"]["streaming"] = "double_buffer"
cfg["method"]["store_precision"] = "FP32"
cfg["method"]["ddf_shifting"] = False
cfg["method"]["les"]["enabled"] = False
cfg["method"]["inlet"]["profile"] = "uniform"
cfg["method"]["inlet"]["scheme"] = str(inlet_scheme)
cfg["method"]["outlet"]["mode"] = "neq_extrap"
cfg["method"]["y_wall_bc"] = "free_slip"
return cfg
def _body_doc(center: Tuple[float, float], *, alpha: float) -> dict:
def _body_doc(center: Tuple[float, float], alpha: float) -> dict:
return {
"objects": [
{
@ -163,13 +162,12 @@ def _body_doc(center: Tuple[float, float], *, alpha: float) -> dict:
def _rfft_spectrum(x: np.ndarray, sample_dt: float) -> Tuple[np.ndarray, np.ndarray]:
v = np.asarray(x, dtype=np.float64)
if v.size < 64:
arr = np.asarray(x, dtype=np.float64)
if arr.size < 64:
return np.zeros(0, dtype=np.float64), np.zeros(0, dtype=np.float64)
v = v - np.mean(v)
win = np.hanning(v.size)
spec = np.abs(np.fft.rfft(v * win)) ** 2
freqs = np.fft.rfftfreq(v.size, d=float(sample_dt))
arr = arr - np.mean(arr)
spec = np.abs(np.fft.rfft(arr * np.hanning(arr.size))) ** 2
freqs = np.fft.rfftfreq(arr.size, d=float(sample_dt))
return freqs.astype(np.float64), spec.astype(np.float64)
@ -183,39 +181,34 @@ def _peak_freq_parabolic(freqs: np.ndarray, spec: np.ndarray, idx: int) -> float
den = y0 - 2.0 * y1 + y2
if abs(den) < 1e-20:
return float(freqs[i])
delta = 0.5 * (y0 - y2) / den
delta = float(np.clip(delta, -1.0, 1.0))
df = float(freqs[i + 1] - freqs[i])
return float(freqs[i]) + delta * df
delta = float(np.clip(0.5 * (y0 - y2) / den, -1.0, 1.0))
return float(freqs[i]) + delta * float(freqs[i + 1] - freqs[i])
def _st_from_lift(lift: np.ndarray, sample_dt: float) -> float:
freqs, spec = _rfft_spectrum(lift, sample_dt=sample_dt)
if freqs.size <= 1:
return float("nan")
# Ignore DC bin.
idx = int(np.argmax(spec[1:])) + 1
f_peak = _peak_freq_parabolic(freqs, spec, idx)
return float(f_peak * D_LATTICE / U_INF)
def _cycle_half_p2p(y: np.ndarray) -> float:
"""Mean half peak-to-peak amplitude over cycles of demeaned signal."""
s = np.asarray(y, dtype=np.float64)
if s.size < 8:
arr = np.asarray(y, dtype=np.float64)
if arr.size < 8:
return float("nan")
d = s - np.mean(s)
crossing = np.where((d[:-1] <= 0.0) & (d[1:] > 0.0))[0]
centered = arr - np.mean(arr)
crossing = np.where((centered[:-1] <= 0.0) & (centered[1:] > 0.0))[0]
if crossing.size >= 2:
amps: List[float] = []
for i in range(crossing.size - 1):
seg = s[crossing[i] + 1 : crossing[i + 1] + 1]
if seg.size < 3:
continue
seg = arr[crossing[i] + 1 : crossing[i + 1] + 1]
if seg.size >= 3:
amps.append(0.5 * (float(np.max(seg)) - float(np.min(seg))))
if amps:
return float(np.mean(amps))
return 0.5 * (float(np.max(s)) - float(np.min(s)))
return 0.5 * (float(np.max(arr)) - float(np.min(arr)))
def _vorticity_z(ux: np.ndarray, uy: np.ndarray) -> np.ndarray:
@ -232,7 +225,6 @@ def _save_vorticity_png(path: str, ux: np.ndarray, uy: np.ndarray, title: str) -
import matplotlib.pyplot as plt
except ImportError:
return
omega = _vorticity_z(ux, uy)
abs_o = np.abs(omega[np.isfinite(omega)])
vmax = float(np.percentile(abs_o, 99.5)) if abs_o.size else 1.0
@ -268,19 +260,24 @@ def _run_one(
field_every: int,
save_vorticity: bool,
) -> Dict[str, Any]:
cfg = _build_cfg(base_cfg, nx=domain.nx, ny=domain.ny, collision=spec.collision, re=spec.re)
bdoc = _body_doc(domain.center, alpha=spec.alpha)
cfg = _build_cfg(
base_cfg,
nx=domain.nx,
ny=domain.ny,
re=spec.re,
inlet_scheme=spec.inlet_scheme,
)
body = _body_doc(domain.center, alpha=spec.alpha)
tmpd = tempfile.mkdtemp(prefix="celeris_kan99b_")
tmpd = tempfile.mkdtemp(prefix="celeris_kan99b_mrt_")
lbm_tmp = os.path.join(tmpd, "config_lbm.json")
body_tmp = os.path.join(tmpd, "config_body.json")
_write_json(lbm_tmp, cfg)
_write_json(body_tmp, bdoc)
_write_json(body_tmp, body)
from CelerisLab import Simulation # noqa: WPS433
sim = Simulation(lbm_config_path=lbm_tmp, body_config_path=body_tmp)
# Contract: body omega is host-side runtime state from alpha conversion.
if sim.bodies.count < 1:
sim.close()
raise RuntimeError("Expected one cylinder in body config.")
@ -294,7 +291,7 @@ def _run_one(
sim.close()
raise ValueError("burn + steps must be >= 1")
steps: List[int] = []
step_hist: List[int] = []
fx_hist: List[float] = []
fy_hist: List[float] = []
field_snapshots: List[str] = []
@ -315,26 +312,26 @@ def _run_one(
stream.synchronize()
sim.bodies.download_obs_full_async(stream)
stream.synchronize()
fvec = sim.bodies.read_force(0)
fx = float(fvec[0])
fy = float(fvec[1])
steps.append(step)
fx_hist.append(fx)
fy_hist.append(fy)
force = sim.bodies.read_force(0)
fx = float(force[0])
fy = float(force[1])
if not np.isfinite(fx) or not np.isfinite(fy):
sim.close()
raise RuntimeError(f"NaN/Inf force at step {step}")
step_hist.append(step)
fx_hist.append(fx)
fy_hist.append(fy)
if field_every > 0 and (step % int(field_every) == 0 or step == total):
stream.synchronize()
macro = sim.get_macroscopic()
save_p = os.path.join(snap_dir, f"macro_step{step:08d}.npz")
snap_path = os.path.join(snap_dir, f"macro_step{step:08d}.npz")
np.savez_compressed(
save_p,
snap_path,
rho=np.asarray(macro["rho"], dtype=np.float32),
ux=np.asarray(macro["ux"], dtype=np.float32),
uy=np.asarray(macro["uy"], dtype=np.float32),
)
field_snapshots.append(save_p)
field_snapshots.append(snap_path)
stream.synchronize()
macro_last = sim.get_macroscopic()
@ -343,15 +340,15 @@ def _run_one(
rho_last = np.asarray(macro_last["rho"], dtype=np.float64).reshape(domain.ny, domain.nx)
sim.close()
step_arr = np.asarray(steps, dtype=np.int64)
step_arr = np.asarray(step_hist, dtype=np.int64)
fx_arr = np.asarray(fx_hist, dtype=np.float64)
fy_arr = np.asarray(fy_hist, dtype=np.float64)
burn_mask = step_arr >= int(spec.burn)
if not np.any(burn_mask):
burn_mask = np.ones_like(step_arr, dtype=bool)
cl = 2.0 * fy_arr / (1.0 * (U_INF ** 2) * D_LATTICE)
cd = 2.0 * fx_arr / (1.0 * (U_INF ** 2) * D_LATTICE)
cl = 2.0 * fy_arr / (U_INF**2 * D_LATTICE)
cd = 2.0 * fx_arr / (U_INF**2 * D_LATTICE)
cl_tail = cl[burn_mask]
cd_tail = cd[burn_mask]
st = _st_from_lift(cl_tail, sample_dt=float(rec))
@ -375,249 +372,126 @@ def _run_one(
ux_last,
uy_last,
title=(
f"Kan99b {spec.phase.upper()} {spec.collision} dom={spec.domain} "
f"Re={spec.re:.0f} alpha={spec.alpha:.3f}"
f"Kan99b {spec.case_id} {spec.variant} MRT dom={spec.domain} "
f"Re={spec.re:.0f} alpha={spec.alpha:.3f} inlet={spec.inlet_scheme}"
),
)
return {
"run_id": run_id,
"phase": spec.phase,
"collision": spec.collision,
"case_id": spec.case_id,
"variant": spec.variant,
"collision": "MRT",
"inlet_scheme": spec.inlet_scheme,
"inlet_profile": "uniform",
"domain": spec.domain,
"re": float(spec.re),
"alpha": float(spec.alpha),
"omega_body": float(_omega_body(spec.alpha)),
"nu": float(_nu_from_re(spec.re)),
"steps": int(spec.steps),
"burn": int(spec.burn),
"burn_in": int(spec.burn),
"total_steps": int(total),
"record_every": int(rec),
"n_samples": int(step_arr.size),
"mean_cd": float(np.mean(cd_tail)),
"mean_cl": float(np.mean(cl_tail)),
"amp_cd": float(amp_cd),
"amp_cl": float(amp_cl),
"St": float(st),
"st": float(st),
"mean_cl": float(np.mean(cl_tail)),
"mean_cd": float(np.mean(cd_tail)),
"amp_cl": float(amp_cl),
"amp_cd": float(amp_cd),
"rho_min_final": float(np.min(rho_last)),
"rho_max_final": float(np.max(rho_last)),
"grid": {"nx": int(domain.nx), "ny": int(domain.ny), "diameter": int(D_LATTICE)},
"beta_real": None,
"Re_real": None,
"re_real": None,
"force_csv": csv_path,
"field_snapshots": field_snapshots,
}
def _alpha_list_from_str(text: str) -> List[float]:
vals: List[float] = []
for t in text.split(","):
t = t.strip()
if t:
vals.append(float(t))
return vals
def _phase_runs(
phase: str,
*,
minimal: bool,
domain_key: str,
collisions: Sequence[str],
alpha_override: Optional[List[float]],
steps: int,
burn: int,
) -> List[RunSpec]:
runs: List[RunSpec] = []
def add_many(
p: str,
ds: Iterable[str],
cs: Iterable[str],
res: Iterable[float],
alphas: Iterable[float],
*,
phase_steps: Optional[int] = None,
phase_burn: Optional[int] = None,
) -> None:
for d in ds:
for c in cs:
for re in res:
for a in alphas:
runs.append(
RunSpec(
phase=p,
collision=str(c).upper(),
domain=d,
re=float(re),
alpha=float(a),
steps=int(phase_steps if phase_steps is not None else steps),
burn=int(phase_burn if phase_burn is not None else burn),
)
)
# Plan-driven defaults.
alpha_c = [0.0, 0.5, 1.0, 1.5, 1.7, 1.8, 1.9, 2.0]
alpha_c_min = [0.0, 1.0, 1.5, 1.8, 2.0]
alpha_d_60 = [0.0, 0.5, 1.0, 1.2, 1.4, 1.6]
alpha_d_160 = [0.0, 0.5, 1.0, 1.5, 1.8, 1.9, 2.0]
alpha_d_min = {60.0: [1.4], 160.0: [1.9]}
anchor_steps = 200_000
anchor_burn = 80_000
near_steps = 240_000
near_burn = 120_000
periodic_steps = 160_000
periodic_burn = 64_000
if phase in ("a", "all"):
add_many("a", ["S", "M", "L"], ["MRT"], [100.0], [1.0], phase_steps=anchor_steps, phase_burn=anchor_burn)
if phase in ("anchor", "b", "all"):
add_many("b", [domain_key], collisions, [100.0], [1.0], phase_steps=anchor_steps, phase_burn=anchor_burn)
if phase in ("c", "all"):
alphas = alpha_override if alpha_override is not None else (alpha_c_min if minimal else alpha_c)
# Near-critical values need longer windows.
for a in alphas:
ps = near_steps if abs(a - 1.8) < 0.11 else periodic_steps
pb = near_burn if abs(a - 1.8) < 0.11 else periodic_burn
add_many("c", [domain_key], collisions, [100.0], [a], phase_steps=ps, phase_burn=pb)
if phase in ("d", "all"):
if minimal:
for re, alphas in alpha_d_min.items():
add_many("d", [domain_key], collisions, [re], alphas, phase_steps=near_steps, phase_burn=near_burn)
else:
add_many("d", [domain_key], collisions, [60.0], alpha_d_60, phase_steps=periodic_steps, phase_burn=periodic_burn)
add_many("d", [domain_key], collisions, [160.0], alpha_d_160, phase_steps=periodic_steps, phase_burn=periodic_burn)
# CLI override for quick tests.
if steps > 0:
for i in range(len(runs)):
runs[i] = RunSpec(
phase=runs[i].phase,
collision=runs[i].collision,
domain=runs[i].domain,
re=runs[i].re,
alpha=runs[i].alpha,
steps=steps,
burn=burn,
)
return runs
def _rel_err(meas: float, ref: float) -> Optional[float]:
if not np.isfinite(meas) or ref == 0:
def _rel_err(measured: float, ref: float) -> Optional[float]:
if not np.isfinite(measured) or ref == 0.0:
return None
return abs(float(meas) - float(ref)) / abs(float(ref))
return abs(float(measured) - float(ref)) / abs(float(ref))
def _phase_a_gate(rows: List[Dict[str, Any]]) -> Dict[str, Any]:
dom = {r["domain"]: r for r in rows}
out: Dict[str, Any] = {"phase": "a", "pass": False}
if not all(k in dom for k in ("S", "M", "L")):
out["error"] = "Phase A needs S, M, L rows."
return out
l = dom["L"]
compare: Dict[str, Any] = {}
for k in ("S", "M"):
r = dom[k]
compare[k] = {
"St": _rel_err(r["st"], l["st"]),
"mean_cl": _rel_err(r["mean_cl"], l["mean_cl"]),
"mean_cd": _rel_err(r["mean_cd"], l["mean_cd"]),
"amp_cl": _rel_err(r["amp_cl"], l["amp_cl"]),
"amp_cd": _rel_err(r["amp_cd"], l["amp_cd"]),
}
choose = "L"
m_ok = all(
(compare["M"][metric] is not None and compare["M"][metric] <= DOMAIN_THRESH[metric])
for metric in DOMAIN_THRESH
)
if m_ok:
choose = "M"
else:
s_ok = all(
(compare["S"][metric] is not None and compare["S"][metric] <= DOMAIN_THRESH[metric])
for metric in DOMAIN_THRESH
)
if s_ok:
choose = "S"
out["compare_vs_L"] = compare
out["recommended_domain"] = choose
out["pass"] = True
return out
def _phase_b_anchor_eval(rows: List[Dict[str, Any]]) -> Dict[str, Any]:
by_coll = {r["collision"]: r for r in rows}
out: Dict[str, Any] = {"phase": "b", "rows": {}}
for coll in ("SRT", "TRT", "MRT"):
row = by_coll.get(coll)
if row is None:
def _k2_anchor_gate(rows: Sequence[Dict[str, Any]]) -> List[Dict[str, Any]]:
"""Evaluate K2 rows against Kan99b anchor tolerances."""
out: List[Dict[str, Any]] = []
for row in rows:
if row.get("case_id") != "K2" or "error" in row:
continue
metrics = {
"St": _rel_err(row["st"], KAN99B_ANCHOR["St"]),
rel = {
"St": _rel_err(row["St"], KAN99B_ANCHOR["St"]),
"mean_cl": _rel_err(row["mean_cl"], KAN99B_ANCHOR["mean_cl"]),
"mean_cd": _rel_err(row["mean_cd"], KAN99B_ANCHOR["mean_cd"]),
"amp_cl": _rel_err(row["amp_cl"], KAN99B_ANCHOR["amp_cl"]),
"amp_cd": _rel_err(row["amp_cd"], KAN99B_ANCHOR["amp_cd"]),
}
out["rows"][coll] = {
"rel_err": metrics,
"pass_bands": {
m: (metrics[m] is not None and metrics[m] <= ANCHOR_BANDS[m]) for m in ANCHOR_BANDS
},
pass_bands = {
key: (rel[key] is not None and rel[key] <= ANCHOR_BANDS[key]) for key in ANCHOR_BANDS
}
out.append(
{
"run_id": row["run_id"],
"variant": row["variant"],
"inlet_scheme": row["inlet_scheme"],
"rel_err": rel,
"pass_bands": pass_bands,
"pass_all": bool(all(pass_bands.values())),
}
)
return out
def _save_summary_plots(rows: List[Dict[str, Any]], out_dir: str) -> None:
try:
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
except ImportError:
return
summary_dir = os.path.join(out_dir, "summary_plots")
os.makedirs(summary_dir, exist_ok=True)
def plot_metric(metric: str, ylabel: str, filename: str) -> None:
fig, ax = plt.subplots(figsize=(8, 5))
for coll in ("SRT", "TRT", "MRT"):
coll_rows = [r for r in rows if r["collision"] == coll]
if not coll_rows:
continue
# Use Re=100 sweep first if present, else all points sorted by alpha.
target = [r for r in coll_rows if abs(r["re"] - 100.0) < 1e-9]
data = target if target else coll_rows
data = sorted(data, key=lambda r: (r["re"], r["alpha"]))
ax.plot(
[r["alpha"] for r in data],
[r[metric] for r in data],
marker="o",
linewidth=1.4,
label=coll,
def _build_runs(
cases: Sequence[KanCase],
*,
domain: str,
include_k2_control: bool,
steps_override: int,
burn_override: int,
) -> List[RunSpec]:
runs: List[RunSpec] = []
for case in cases:
steps = int(steps_override) if steps_override > 0 else int(case.steps)
burn = int(burn_override) if burn_override > 0 else int(case.burn)
runs.append(
RunSpec(
case_id=case.case_id,
variant="baseline",
domain=domain,
re=case.re,
alpha=case.alpha,
inlet_scheme="regularized",
steps=steps,
burn=burn,
)
ax.set_xlabel("alpha")
ax.set_ylabel(ylabel)
ax.set_title(f"{ylabel} vs alpha")
ax.grid(True, alpha=0.3)
ax.legend(loc="best")
fig.tight_layout()
fig.savefig(os.path.join(summary_dir, filename), dpi=150, bbox_inches="tight")
plt.close(fig)
plot_metric("mean_cl", "mean C_L", "mean_cl_vs_alpha.png")
plot_metric("mean_cd", "mean C_D", "mean_cd_vs_alpha.png")
plot_metric("amp_cl", "C'_L (half peak-to-peak)", "amp_cl_vs_alpha.png")
plot_metric("st", "St", "st_vs_alpha.png")
)
if include_k2_control and case.case_id == "K2":
runs.append(
RunSpec(
case_id=case.case_id,
variant="k2_inlet_control",
domain=domain,
re=case.re,
alpha=case.alpha,
inlet_scheme="zou_he_local",
steps=steps,
burn=burn,
)
)
return runs
def main() -> int:
ap = argparse.ArgumentParser(description="Kan99b rotating-cylinder validation driver")
ap.add_argument("--phase", default="all", choices=("anchor", "a", "b", "c", "d", "all"))
ap.add_argument("--minimal", action="store_true", help="Run reduced minimum set from the plan.")
ap.add_argument("--domain", default="M", choices=("S", "M", "L"), help="Default domain for phases B/C/D.")
ap.add_argument("--collision", default="all", choices=("SRT", "TRT", "MRT", "all"))
ap.add_argument("--alpha", type=float, default=None, help="Single alpha override (for c/d phases).")
ap.add_argument("--alpha-list", type=str, default="", help="Comma-separated alpha list override.")
ap = argparse.ArgumentParser(description="Kan99b MRT-only primary matrix runner")
ap.add_argument("--case", default="all", help='Case id K1-K5/K0 or "all"')
ap.add_argument("--include-k0", action="store_true", help="Include optional K0 baseline.")
ap.add_argument("--no-k2-control", action="store_true", help="Disable K2 zou_he_local control run.")
ap.add_argument("--domain", default="M", choices=("S", "M", "L"))
ap.add_argument("--steps", type=int, default=0, help="Override run steps for all selected runs.")
ap.add_argument("--burn", type=int, default=0, help="Override burn steps for all selected runs.")
ap.add_argument("--record-every", type=int, default=100)
@ -632,112 +506,112 @@ def main() -> int:
print(f"Missing base config: {_DEFAULT_LBM}", file=sys.stderr)
return 2
base_cfg = _load_json(_DEFAULT_LBM)
domains = _domain_specs()
out_dir = os.path.abspath(args.out_dir)
os.makedirs(out_dir, exist_ok=True)
sel = str(args.case).upper()
allowed = {case.case_id for case in CASES} | {"K0", "ALL"}
if sel not in allowed:
print("--case must be K0,K1,K2,K3,K4,K5,all", file=sys.stderr)
return 2
collisions = ["SRT", "TRT", "MRT"] if args.collision == "all" else [str(args.collision).upper()]
alpha_override: Optional[List[float]] = None
if args.alpha is not None:
alpha_override = [float(args.alpha)]
elif args.alpha_list.strip():
alpha_override = _alpha_list_from_str(args.alpha_list)
if args.smoke:
o_steps = 2000
o_burn = 800
selected: List[KanCase] = []
if sel == "ALL":
selected.extend(CASES)
if args.include_k0:
selected.insert(0, CASE_K0)
elif sel == "K0":
selected.append(CASE_K0)
else:
o_steps = max(0, int(args.steps))
o_burn = max(0, int(args.burn))
selected.extend(case for case in CASES if case.case_id == sel)
runs = _phase_runs(
args.phase,
minimal=bool(args.minimal),
domain_key=args.domain,
collisions=collisions,
alpha_override=alpha_override,
steps=o_steps,
burn=o_burn,
)
if not runs:
if not selected:
print("No runs selected.", file=sys.stderr)
return 2
domains = _domain_specs()
steps_override = 2000 if args.smoke else max(0, int(args.steps))
burn_override = 800 if args.smoke else max(0, int(args.burn))
runs = _build_runs(
selected,
domain=args.domain,
include_k2_control=not bool(args.no_k2_control),
steps_override=steps_override,
burn_override=burn_override,
)
out_dir = os.path.abspath(args.out_dir)
os.makedirs(out_dir, exist_ok=True)
rows: List[Dict[str, Any]] = []
contract = {
"U_inf": U_INF,
"D_lattice": D_LATTICE,
"R_lattice": R_LATTICE,
"nu_formula": "nu = U_inf * D / Re = 0.9 / Re",
"omega_formula": "omega_body = 2 * alpha * U_inf / D = 0.002 * alpha",
"method_contract": {
"inlet_profile": "uniform",
"y_wall_bc": "free_slip",
"outlet_mode": "neq_extrap",
"streaming": "double_buffer",
"store_precision": "FP32",
"les_enabled": False,
},
}
for spec in runs:
dspec = domains[spec.domain]
print(
f"--- {spec.phase.upper()} {spec.collision} dom={spec.domain} Re={spec.re:.0f} "
f"alpha={spec.alpha:.3f} burn={spec.burn} steps={spec.steps} ---",
f"--- {spec.case_id} {spec.variant} MRT dom={spec.domain} Re={spec.re:.0f} "
f"alpha={spec.alpha:.3f} inlet={spec.inlet_scheme} burn={spec.burn} steps={spec.steps} ---",
flush=True,
)
try:
row = _run_one(
spec,
domain=dspec,
domain=domains[spec.domain],
base_cfg=base_cfg,
out_dir=out_dir,
record_every=max(1, int(args.record_every)),
field_every=max(0, int(args.field_every)),
save_vorticity=bool(args.save_vorticity),
)
except Exception as e: # noqa: BLE001
rows.append(
{
"run_id": _run_id(spec),
"phase": spec.phase,
"collision": spec.collision,
"domain": spec.domain,
"re": float(spec.re),
"alpha": float(spec.alpha),
"error": str(e),
}
)
print(f"FAILED: {e}", flush=True)
continue
if spec.case_id == "K2":
rel = _rel_err(row["St"], KAN99B_ANCHOR["St"])
row["St_error_pct"] = 100.0 * rel if rel is not None else None
else:
row["St_error_pct"] = None
rows.append(row)
print(
" "
f"St={row['st']:.5f} mean_CL={row['mean_cl']:.4f} mean_CD={row['mean_cd']:.4f} "
f" St={row['St']:.5f} mean_CL={row['mean_cl']:.4f} mean_CD={row['mean_cd']:.4f} "
f"C'L={row['amp_cl']:.4f} C'D={row['amp_cd']:.4f}",
flush=True,
)
except Exception as exc: # noqa: BLE001
rows.append(
{
"run_id": _run_id(spec),
"case_id": spec.case_id,
"variant": spec.variant,
"collision": "MRT",
"inlet_scheme": spec.inlet_scheme,
"inlet_profile": "uniform",
"domain": spec.domain,
"re": float(spec.re),
"alpha": float(spec.alpha),
"steps": int(spec.steps),
"burn_in": int(spec.burn),
"error": str(exc),
}
)
print(f"FAILED: {exc}", flush=True)
k2_gate = _k2_anchor_gate(rows)
print("\n=== Kan99b K2 gate summary ===", flush=True)
print(json.dumps({"k2_runs": k2_gate}, indent=2), flush=True)
# Summary table outputs
summary_csv = os.path.join(out_dir, "summary_runs.csv")
csv_keys = [
"run_id",
"phase",
"case_id",
"variant",
"collision",
"inlet_scheme",
"inlet_profile",
"domain",
"re",
"alpha",
"omega_body",
"nu",
"burn",
"burn_in",
"steps",
"total_steps",
"record_every",
"n_samples",
"St",
"st",
"St_error_pct",
"mean_cl",
"mean_cd",
"amp_cl",
@ -748,30 +622,31 @@ def main() -> int:
"error",
]
with open(summary_csv, "w", newline="", encoding="utf-8") as f:
w = csv.DictWriter(f, fieldnames=csv_keys)
w.writeheader()
for r in rows:
w.writerow({k: r.get(k, "") for k in csv_keys})
phase_reports: Dict[str, Any] = {}
phase_a_rows = [r for r in rows if r.get("phase") == "a" and "error" not in r]
if phase_a_rows:
phase_reports["a"] = _phase_a_gate(phase_a_rows)
phase_b_rows = [r for r in rows if r.get("phase") == "b" and "error" not in r]
if phase_b_rows:
phase_reports["b"] = _phase_b_anchor_eval(phase_b_rows)
_save_summary_plots([r for r in rows if "error" not in r], out_dir)
writer = csv.DictWriter(f, fieldnames=csv_keys)
writer.writeheader()
for row in rows:
writer.writerow({k: row.get(k, "") for k in csv_keys})
summary = {
"contract": contract,
"contract": {
"collision": "MRT",
"primary_inlet_scheme": "regularized",
"k2_control_inlet_scheme": "zou_he_local",
"inlet_profile": "uniform",
"y_wall_bc": "free_slip",
"outlet_mode": "neq_extrap",
"streaming": "double_buffer",
"store_precision": "FP32",
"les_enabled": False,
},
"requested": {
"phase": args.phase,
"minimal": bool(args.minimal),
"case": args.case,
"include_k0": bool(args.include_k0),
"include_k2_control": not bool(args.no_k2_control),
"domain": args.domain,
"collision": args.collision,
"steps_override": int(o_steps),
"burn_override": int(o_burn),
"smoke": bool(args.smoke),
"steps_override": int(steps_override),
"burn_override": int(burn_override),
"record_every": int(args.record_every),
"field_every": int(args.field_every),
"save_vorticity": bool(args.save_vorticity),
@ -781,7 +656,7 @@ def main() -> int:
"completed_runs": sum(1 for r in rows if "error" not in r),
"failed_runs": sum(1 for r in rows if "error" in r),
},
"phase_reports": phase_reports,
"k2_gate": k2_gate,
"rows": rows,
}
json_out = (
@ -789,11 +664,13 @@ def main() -> int:
if args.json_out.strip()
else os.path.join(out_dir, "summary_runs.json")
)
json_out_dir = os.path.dirname(json_out)
if json_out_dir:
os.makedirs(json_out_dir, exist_ok=True)
_write_json(json_out, summary)
print(f"Wrote: {summary_csv}", flush=True)
print(f"Wrote: {json_out}", flush=True)
print(f"Wrote: {os.path.join(out_dir, 'summary_plots')}", flush=True)
return 0

View File

@ -0,0 +1,309 @@
# CelerisLab/tests/run_kan99b_streakline.py
"""Kan99b streakline demo using CelerisLab common streakline engine.
Modes:
- online: run CFD and compute streakline in memory (no velocity dump).
- offline: load velocity snapshots and reconstruct streakline.
"""
from __future__ import annotations
import argparse
import json
import os
import tempfile
from dataclasses import dataclass
from typing import List, Tuple
import numpy as np
from CelerisLab import Simulation
from CelerisLab.common.streakline import (
FlowFrame,
IntegratorConfig,
ReleaseConfig,
build_release_points,
estimate_sampling_plan,
render_streakline_density,
run_streakline_offline,
run_streakline_online,
)
_REPO = os.path.abspath(os.path.join(os.path.dirname(__file__), ".."))
_DEFAULT_LBM = os.path.join(_REPO, "src", "CelerisLab", "configs", "config_lbm.json")
U_INF = 0.03
D_LATTICE = 30.0
R_LATTICE = 15.0
KAN99B_ST_REF = 0.1655
@dataclass(frozen=True)
class DomainSpec:
key: str
nx: int
ny: int
center: Tuple[float, float]
def _domain_specs() -> dict:
return {
"S": DomainSpec("S", 1081, 481, (360.0, 240.0)),
"M": DomainSpec("M", 1351, 601, (450.0, 300.0)),
"L": DomainSpec("L", 1801, 721, (600.0, 360.0)),
}
def _load_json(path: str) -> dict:
with open(path, "r", encoding="utf-8") as f:
return json.load(f)
def _write_json(path: str, payload: dict) -> None:
with open(path, "w", encoding="utf-8") as f:
json.dump(payload, f, indent=2)
def _nu_from_re(reynolds: float) -> float:
return U_INF * D_LATTICE / float(reynolds)
def _omega_body(alpha: float) -> float:
return 2.0 * float(alpha) * U_INF / D_LATTICE
def _build_cfg(base_cfg: dict, *, nx: int, ny: int, re: float, inlet_scheme: str) -> dict:
cfg = json.loads(json.dumps(base_cfg))
cfg["grid"]["nx"] = int(nx)
cfg["grid"]["ny"] = int(ny)
cfg["grid"]["nz"] = 1
cfg["physics"]["velocity"] = float(U_INF)
cfg["physics"]["viscosity"] = float(_nu_from_re(re))
cfg["physics"]["rho"] = 1.0
cfg["method"]["collision"] = "MRT"
cfg["method"]["streaming"] = "double_buffer"
cfg["method"]["store_precision"] = "FP32"
cfg["method"]["ddf_shifting"] = False
cfg["method"]["les"]["enabled"] = False
cfg["method"]["inlet"]["profile"] = "uniform"
cfg["method"]["inlet"]["scheme"] = str(inlet_scheme)
cfg["method"]["outlet"]["mode"] = "neq_extrap"
cfg["method"]["y_wall_bc"] = "free_slip"
return cfg
def _build_simulation(*, domain: DomainSpec, re: float, alpha: float, inlet_scheme: str) -> Simulation:
base_cfg = _load_json(_DEFAULT_LBM)
cfg = _build_cfg(base_cfg, nx=domain.nx, ny=domain.ny, re=re, inlet_scheme=inlet_scheme)
body_doc = {
"objects": [
{
"type": "cylinder",
"center": [float(domain.center[0]), float(domain.center[1])],
"radius": float(R_LATTICE),
"omega": float(_omega_body(alpha)),
}
]
}
tmpd = tempfile.mkdtemp(prefix="celeris_streakline_online_")
lbm_tmp = os.path.join(tmpd, "config_lbm.json")
body_tmp = os.path.join(tmpd, "config_body.json")
_write_json(lbm_tmp, cfg)
_write_json(body_tmp, body_doc)
sim = Simulation(lbm_config_path=lbm_tmp, body_config_path=body_tmp)
sim.bodies.get(0).state.omega = np.float32(_omega_body(alpha))
sim.initialize()
return sim
def _default_base_release_points(center: Tuple[float, float]) -> np.ndarray:
x_rel = float(center[0] - 6.0 * D_LATTICE)
y0 = float(center[1])
return np.array(
[
[x_rel, y0 - 18.0],
[x_rel, y0 - 6.0],
[x_rel, y0 + 6.0],
[x_rel, y0 + 18.0],
],
dtype=np.float64,
)
def _load_offline_frames(snapshot_dir: str) -> List[FlowFrame]:
files = sorted(
[
os.path.join(snapshot_dir, name)
for name in os.listdir(snapshot_dir)
if name.endswith(".npz") and name.startswith("vel_step")
]
)
frames: List[FlowFrame] = []
for path in files:
data = np.load(path)
step = int(np.asarray(data["step"]).reshape(-1)[0])
frames.append(
FlowFrame(
step=step,
ux=np.asarray(data["ux"], dtype=np.float64),
uy=np.asarray(data["uy"], dtype=np.float64),
)
)
return frames
def main() -> int:
ap = argparse.ArgumentParser(description="Kan99b streakline demo (online/offline)")
ap.add_argument("--mode", default="online", choices=("online", "offline"))
ap.add_argument("--domain", default="M", choices=("S", "M", "L"))
ap.add_argument("--re", type=float, default=100.0)
ap.add_argument("--alpha", type=float, default=1.0)
ap.add_argument("--inlet-scheme", default="regularized", choices=("regularized", "zou_he_local"))
ap.add_argument("--start-step", type=int, default=60_000)
ap.add_argument("--sample-every", type=int, default=0, help="0 uses recommended value.")
ap.add_argument("--n-snapshots", type=int, default=0, help="0 uses recommended value.")
ap.add_argument("--snapshot-dir", type=str, default="", help="Required in offline mode.")
ap.add_argument("--release-mode", default="strip", choices=("point", "line", "strip"))
ap.add_argument("--line-span", type=float, default=0.0)
ap.add_argument("--line-count", type=int, default=1)
ap.add_argument(
"--downstream-count",
type=int,
default=5,
help="Particles per seed along +x (denser streak).",
)
ap.add_argument(
"--downstream-spacing",
type=float,
default=1.0,
help="Lattice spacing between x-staggered release points.",
)
ap.add_argument("--inject-per-seed", type=int, default=2)
ap.add_argument("--alpha-t", type=float, default=0.2)
ap.add_argument("--alpha-x", type=float, default=0.4)
ap.add_argument("--diffusion-coeff", type=float, default=0.0)
ap.add_argument(
"--out-dir",
type=str,
default=os.path.join(_REPO, "tests", "output", "streakline", "kan99b_k2_online"),
)
args = ap.parse_args()
domain = _domain_specs()[args.domain]
out_dir = os.path.abspath(args.out_dir)
os.makedirs(out_dir, exist_ok=True)
plan = estimate_sampling_plan(st_ref=KAN99B_ST_REF, diameter=D_LATTICE, u_ref=U_INF)
sample_every = int(args.sample_every) if int(args.sample_every) > 0 else int(plan["save_every_recommended"])
n_snapshots = int(args.n_snapshots) if int(args.n_snapshots) > 0 else int(plan["snapshot_count_recommended"])
release_cfg = ReleaseConfig(
mode=args.release_mode,
line_span=float(args.line_span),
line_count=max(1, int(args.line_count)),
downstream_count=max(1, int(args.downstream_count)),
downstream_spacing=float(args.downstream_spacing),
inject_per_seed=max(1, int(args.inject_per_seed)),
)
integrator_cfg = IntegratorConfig(
alpha_t=float(args.alpha_t),
alpha_x=float(args.alpha_x),
diffusion_coeff=float(args.diffusion_coeff),
)
base_release_points = _default_base_release_points(domain.center)
dense_release_points = build_release_points(base_release_points, release_cfg)
if args.mode == "online":
sim = _build_simulation(
domain=domain,
re=float(args.re),
alpha=float(args.alpha),
inlet_scheme=args.inlet_scheme,
)
try:
particles, ages, diag = run_streakline_online(
sim,
start_step=int(args.start_step),
sample_every=int(sample_every),
n_samples=int(n_snapshots),
release_points=base_release_points,
release_cfg=release_cfg,
integrator_cfg=integrator_cfg,
solid_center=domain.center,
solid_radius=R_LATTICE,
)
finally:
sim.close()
frame_count = int(n_snapshots)
frame_source = "in-memory online sampling"
else:
if not args.snapshot_dir:
raise ValueError("--snapshot-dir is required in offline mode.")
frames = _load_offline_frames(os.path.abspath(args.snapshot_dir))
if len(frames) < 2:
raise RuntimeError("Offline mode needs at least two velocity snapshots.")
particles, ages, diag = run_streakline_offline(
frames,
nx=domain.nx,
ny=domain.ny,
release_points=base_release_points,
release_cfg=release_cfg,
integrator_cfg=integrator_cfg,
solid_center=domain.center,
solid_radius=R_LATTICE,
)
frame_count = len(frames)
frame_source = os.path.abspath(args.snapshot_dir)
render_info = render_streakline_density(
particles,
ages,
nx=domain.nx,
ny=domain.ny,
out_path=os.path.join(out_dir, "streakline.png"),
release_points=dense_release_points,
solid_center=domain.center,
solid_radius=R_LATTICE,
age_decay_steps=integrator_cfg.age_decay_steps,
blur_sigma=1.2,
title=f"Kan99b streakline ({args.mode}, {args.release_mode})",
)
meta = {
"case": {
"domain": args.domain,
"re": float(args.re),
"alpha": float(args.alpha),
"inlet_scheme": args.inlet_scheme,
"collision": "MRT",
},
"mode": args.mode,
"sampling_estimate": plan,
"sampling_used": {
"start_step": int(args.start_step),
"sample_every": int(sample_every),
"n_snapshots": int(n_snapshots),
"frame_source": frame_source,
"frames_used": int(frame_count),
},
"release": {
"config": vars(args),
"base_points": base_release_points.tolist(),
"dense_point_count": int(dense_release_points.shape[0]),
"dense_points_preview": dense_release_points[: min(12, dense_release_points.shape[0])].tolist(),
},
"diagnostics": diag,
"render": render_info,
}
_write_json(os.path.join(out_dir, "streakline_meta.json"), meta)
print(f"Recommended sample_every: {plan['save_every_recommended']} steps")
print(f"Recommended snapshots: {plan['snapshot_count_recommended']}")
print(f"Mode: {args.mode} | frames used: {frame_count}")
print(f"Dense release points: {dense_release_points.shape[0]}")
print(f"Output image: {render_info['image_path']}")
return 0
if __name__ == "__main__":
raise SystemExit(main())

269
tests/run_perf_baseline.py Normal file
View File

@ -0,0 +1,269 @@
"""Performance baseline for minimal host-interference LBM stepping.
This script builds a temporary config, runs warmup + measured batches, and
reports MLUPS under a FluidX3D-like benchmark mindset:
- keep the main loop on GPU (`stepper.step`)
- avoid host downloads by default
- selectively enable host-touch paths to quantify overhead
- sweep `inlet.scheme` to check stability-sensitive combinations
Usage::
python tests/run_perf_baseline.py
python tests/run_perf_baseline.py --lattice-model D2Q9 --nx 384 --ny 192 --steps 30000
python tests/run_perf_baseline.py --macro-every 500 --ddf-every 1000
python tests/run_perf_baseline.py --with-cylinder --obs-every 50
"""
from __future__ import annotations
import argparse
import json
import os
import tempfile
import time
from typing import Any, Dict, List
import pycuda.driver as cuda
_REPO = os.path.abspath(os.path.join(os.path.dirname(__file__), ".."))
_DEFAULT_LBM = os.path.join(_REPO, "src", "CelerisLab", "configs", "config_lbm.json")
def _load_json(path: str) -> dict:
with open(path, "r", encoding="utf-8") as f:
return json.load(f)
def _write_json(path: str, payload: dict) -> None:
with open(path, "w", encoding="utf-8") as f:
json.dump(payload, f, indent=2)
def _build_lbm_cfg(base: dict, args: argparse.Namespace) -> dict:
cfg = json.loads(json.dumps(base))
cfg["grid"]["lattice_model"] = args.lattice_model
cfg["grid"]["nx"] = int(args.nx)
cfg["grid"]["ny"] = int(args.ny)
cfg["grid"]["nz"] = int(args.nz)
cfg["physics"]["viscosity"] = float(args.viscosity)
cfg["physics"]["velocity"] = float(args.velocity)
cfg["physics"]["rho"] = float(args.rho)
cfg["method"]["collision"] = str(args.collision).upper()
cfg["method"]["streaming"] = str(args.streaming)
cfg["method"]["store_precision"] = str(args.store_precision).upper()
cfg["method"]["ddf_shifting"] = bool(args.ddf_shifting)
cfg["method"]["les"]["enabled"] = bool(args.enable_les)
cfg["method"]["outlet"]["mode"] = str(args.outlet_mode)
cfg["method"]["inlet"]["profile"] = str(args.inlet_profile)
# Expose inlet scheme as a benchmark axis; useful when stability depends
# on collision/inlet coupling.
cfg["method"]["inlet"]["scheme"] = str(args.inlet_scheme)
cfg["method"]["y_wall_bc"] = str(args.y_wall_bc)
cfg["cuda"]["threads_per_block"] = int(args.threads_per_block)
cfg["cuda"]["compute_capability"] = str(args.compute_capability)
return cfg
def _build_body_cfg(args: argparse.Namespace) -> dict:
if not args.with_cylinder:
return {"objects": []}
cx = 0.5 * float(args.nx)
cy = 0.5 * float(args.ny)
radius = max(2.0, min(float(args.nx), float(args.ny)) * 0.08)
return {
"objects": [
{
"type": "cylinder",
"center": [cx, cy],
"radius": radius,
}
]
}
def _maybe_probe_macroscopic(sim: Any, step: int, every: int, repeat: int) -> None:
if every > 0 and step % every == 0:
for _ in range(max(1, int(repeat))):
_ = sim.get_macroscopic()
def _maybe_probe_ddf(sim: Any, step: int, every: int, repeat: int) -> None:
if every > 0 and step % every == 0:
for _ in range(max(1, int(repeat))):
_ = sim.get_ddf()
def _maybe_checkpoint(sim: Any, step: int, every: int, out_dir: str) -> None:
if every > 0 and step % every == 0:
path = os.path.join(out_dir, f"checkpoint_step_{step:09d}.h5")
sim.save_checkpoint(path)
def _maybe_probe_obs(sim: Any, stream: cuda.Stream, step: int, every: int) -> None:
if every <= 0 or step % every != 0:
return
if sim.bodies.count <= 0:
return
sim.bodies.download_obs_full_async(stream)
stream.synchronize()
_ = sim.bodies.read_force(0)
def run(args: argparse.Namespace) -> Dict[str, Any]:
if not os.path.isfile(_DEFAULT_LBM):
raise FileNotFoundError(f"Base config missing: {_DEFAULT_LBM}")
base = _load_json(_DEFAULT_LBM)
lbm_cfg = _build_lbm_cfg(base, args)
body_cfg = _build_body_cfg(args)
tmpd = tempfile.mkdtemp(prefix="celeris_perf_baseline_")
lbm_path = os.path.join(tmpd, "config_lbm.json")
body_path = os.path.join(tmpd, "config_body.json")
ckpt_dir = os.path.join(tmpd, "checkpoints")
os.makedirs(ckpt_dir, exist_ok=True)
_write_json(lbm_path, lbm_cfg)
_write_json(body_path, body_cfg)
from CelerisLab import Simulation # noqa: WPS433
sim = Simulation(lbm_config_path=lbm_path, body_config_path=body_path, device_id=args.device_id)
sim.initialize()
stream = cuda.Stream()
total_cells = int(args.nx) * int(args.ny) * int(args.nz)
# Warmup before measurement window.
warmup_done = 0
while warmup_done < int(args.warmup_steps):
chunk = min(int(args.batch_steps), int(args.warmup_steps) - warmup_done)
sim.stepper.step(chunk, action_gpu=sim.bodies.action_gpu, obs_gpu=sim.bodies.obs_gpu, stream=stream)
warmup_done += chunk
stream.synchronize()
measured_batch_s: List[float] = []
measured_steps = int(args.steps)
done = 0
t0 = time.perf_counter()
while done < measured_steps:
chunk = min(int(args.batch_steps), measured_steps - done)
step_start = time.perf_counter()
sim.stepper.step(chunk, action_gpu=sim.bodies.action_gpu, obs_gpu=sim.bodies.obs_gpu, stream=stream)
stream.synchronize()
step_end = time.perf_counter()
measured_batch_s.append(step_end - step_start)
done += chunk
global_step = sim.stepper.step_count
_maybe_probe_macroscopic(sim, global_step, int(args.macro_every), int(args.macro_repeat))
_maybe_probe_ddf(sim, global_step, int(args.ddf_every), int(args.ddf_repeat))
_maybe_checkpoint(sim, global_step, int(args.checkpoint_every), ckpt_dir)
_maybe_probe_obs(sim, stream, global_step, int(args.obs_every))
stream.synchronize()
elapsed_s = time.perf_counter() - t0
# Optional final sanity readback (outside core timing path by default).
if args.final_macro_snapshot:
_ = sim.get_macroscopic()
sim.close()
mlups = (total_cells * measured_steps) / max(elapsed_s, 1e-12) / 1.0e6
batch_ms = [1000.0 * x for x in measured_batch_s]
batch_ms_sorted = sorted(batch_ms)
p50 = batch_ms_sorted[len(batch_ms_sorted) // 2] if batch_ms_sorted else 0.0
p90 = batch_ms_sorted[min(len(batch_ms_sorted) - 1, int(0.9 * (len(batch_ms_sorted) - 1)))] if batch_ms_sorted else 0.0
return {
"benchmark": "celerislab_stepper_baseline",
"device_id": int(args.device_id),
"lattice_model": args.lattice_model,
"grid": {"nx": int(args.nx), "ny": int(args.ny), "nz": int(args.nz)},
"collision": str(args.collision).upper(),
"inlet_scheme": str(args.inlet_scheme),
"streaming": str(args.streaming),
"store_precision": str(args.store_precision).upper(),
"steps": measured_steps,
"warmup_steps": int(args.warmup_steps),
"batch_steps": int(args.batch_steps),
"mlups": float(mlups),
"elapsed_s": float(elapsed_s),
"batch_ms_p50": float(p50),
"batch_ms_p90": float(p90),
"overhead_switches": {
"macro_every": int(args.macro_every),
"macro_repeat": int(args.macro_repeat),
"ddf_every": int(args.ddf_every),
"ddf_repeat": int(args.ddf_repeat),
"checkpoint_every": int(args.checkpoint_every),
"obs_every": int(args.obs_every),
"with_cylinder": bool(args.with_cylinder),
"final_macro_snapshot": bool(args.final_macro_snapshot),
},
}
def parse_args() -> argparse.Namespace:
ap = argparse.ArgumentParser(description="CelerisLab pure-step performance baseline")
ap.add_argument("--device-id", type=int, default=0)
ap.add_argument("--lattice-model", choices=("D2Q9", "D3Q19"), default="D3Q19")
ap.add_argument("--nx", type=int, default=256)
ap.add_argument("--ny", type=int, default=256)
ap.add_argument("--nz", type=int, default=256)
ap.add_argument("--steps", type=int, default=3000)
ap.add_argument("--warmup-steps", type=int, default=400)
ap.add_argument("--batch-steps", type=int, default=100)
ap.add_argument("--collision", choices=("SRT", "TRT", "MRT"), default="SRT")
ap.add_argument("--streaming", choices=("double_buffer", "esopull"), default="double_buffer")
ap.add_argument("--store-precision", choices=("FP32", "FP16S"), default="FP32")
ap.add_argument("--ddf-shifting", action="store_true")
ap.add_argument("--enable-les", action="store_true")
ap.add_argument("--outlet-mode", choices=("neq_extrap", "zero_gradient", "blended"), default="neq_extrap")
ap.add_argument("--inlet-profile", choices=("uniform", "parabolic"), default="uniform")
ap.add_argument(
"--inlet-scheme",
choices=("zou_he_local", "channel_stabilized", "equilibrium", "regularized"),
default="zou_he_local",
)
ap.add_argument("--y-wall-bc", choices=("bounce_back", "free_slip"), default="bounce_back")
ap.add_argument("--viscosity", type=float, default=0.0035)
ap.add_argument("--velocity", type=float, default=0.03)
ap.add_argument("--rho", type=float, default=1.0)
ap.add_argument("--threads-per-block", type=int, default=256)
ap.add_argument("--compute-capability", type=str, default="auto")
# Overhead attribution toggles.
ap.add_argument("--macro-every", type=int, default=0, help="Call get_macroscopic() every N steps (0=off)")
ap.add_argument("--macro-repeat", type=int, default=1, help="Repeat get_macroscopic() calls per probe step")
ap.add_argument("--ddf-every", type=int, default=0, help="Call get_ddf() every N steps (0=off)")
ap.add_argument("--ddf-repeat", type=int, default=1, help="Repeat get_ddf() calls per probe step")
ap.add_argument("--checkpoint-every", type=int, default=0, help="Save checkpoint every N steps (0=off)")
ap.add_argument("--obs-every", type=int, default=0, help="Download object obs every N steps (0=off)")
ap.add_argument("--with-cylinder", action="store_true", help="Inject one cylinder object (needed for obs probes)")
ap.add_argument("--final-macro-snapshot", action="store_true")
ap.add_argument("--json-out", type=str, default="", help="Optional path to save metrics JSON")
return ap.parse_args()
def main() -> int:
args = parse_args()
result = run(args)
print(json.dumps(result, indent=2))
if args.json_out.strip():
out = os.path.abspath(args.json_out)
os.makedirs(os.path.dirname(out), exist_ok=True)
_write_json(out, result)
print(f"Wrote: {out}")
return 0
if __name__ == "__main__":
raise SystemExit(main())

View File

@ -1,238 +0,0 @@
# tests/run_sah04_case9_flow_diagnostics.py
"""Sah04 matrix case 9: velocity profiles at several x and full-field vorticity.
Loads geometry and setup from ``tests/tests/run_sah04_st_matrix.py``, runs one
simulation to the last LBM step, then writes PNG + JSON under
``tests/output/sah04_case9_flow_diag/`` (default).
Usage::
conda run -n pycuda_3_10 python tests/run_sah04_case9_flow_diagnostics.py
conda run -n pycuda_3_10 python tests/run_sah04_case9_flow_diagnostics.py --steps 60000 --burn 20000
Design::
Intended for inlet / channel diagnostics after boundary fixes; not part of
the St matrix gate.
"""
from __future__ import annotations
import argparse
import importlib.util
import json
import os
import shutil
import sys
import tempfile
from typing import Any, Dict, List, Sequence
import numpy as np
import pycuda.driver as cuda
_REPO = os.path.abspath(os.path.join(os.path.dirname(__file__), ".."))
_mtx_path = os.path.join(_REPO, "tests", "tests", "run_sah04_st_matrix.py")
_spec = importlib.util.spec_from_file_location("sah04_mtx", _mtx_path)
if _spec is None or _spec.loader is None:
raise RuntimeError(f"Cannot load {_mtx_path}")
sah = importlib.util.module_from_spec(_spec)
sys.modules["sah04_mtx"] = sah
_spec.loader.exec_module(sah)
def _inlet_parabolic_ref_u(y: np.ndarray, ny: int, u0_mean: float) -> np.ndarray:
"""Match ``inlet_target_u`` in ``inlet_outlet.cuh`` (parabolic, mean ``u0_mean``)."""
h = max(float(ny - 2), 1.0)
yc = np.clip(y.astype(np.float64), 1.0, float(ny - 2))
eta = (yc - 0.5) / h
shape = np.maximum(0.0, 4.0 * eta * (1.0 - eta))
return u0_mean * 1.5 * shape
def _default_x_indices(nx: int) -> List[int]:
raw = [1, 2, 3, 5, 10, 20, 40, 80, 160, 320, 640, 960, 1200, 1600, 2000, 2200]
out: List[int] = []
for x in raw:
if 1 <= x < nx - 1:
out.append(int(x))
mid = max(1, (nx - 2) // 2)
if mid not in out:
out.append(mid)
out = sorted(set(out))
return out
def _plot_profiles(
out_png: str,
*,
ux: np.ndarray,
uy: np.ndarray,
x_indices: Sequence[int],
ny: int,
u0_mean: float,
title: str,
) -> None:
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
y = np.arange(ny, dtype=np.float64)
u_ref = _inlet_parabolic_ref_u(y, ny, u0_mean)
fig, axes = plt.subplots(1, 2, figsize=(12, 6), sharey=True)
for xi in x_indices:
axes[0].plot(ux[:, xi], y, label=f"x={xi}", linewidth=1.2)
axes[1].plot(uy[:, xi], y, label=f"x={xi}", linewidth=1.2)
axes[0].plot(u_ref, y, "k--", linewidth=1.5, label="inlet parabolic ref")
axes[0].set_xlabel("u_x")
axes[0].set_ylabel("y (lattice)")
axes[0].set_title("u_x (y) at selected x")
axes[0].legend(loc="best", fontsize=7)
axes[0].grid(True, alpha=0.3)
axes[1].set_xlabel("u_y")
axes[1].set_title("u_y (y) at selected x")
axes[1].legend(loc="best", fontsize=7)
axes[1].grid(True, alpha=0.3)
fig.suptitle(title)
fig.tight_layout()
os.makedirs(os.path.dirname(os.path.abspath(out_png)) or ".", exist_ok=True)
fig.savefig(out_png, dpi=150, bbox_inches="tight")
plt.close(fig)
def main() -> int:
ap = argparse.ArgumentParser(description="Sah04 case 9 flow: profiles + vorticity")
ap.add_argument("--collision", default="MRT", help="SRT, TRT, or MRT")
ap.add_argument("--steps", type=int, default=60_000)
ap.add_argument("--burn", type=int, default=20_000)
ap.add_argument("--record-every", type=int, default=5)
ap.add_argument("--outlet", default="neq_extrap")
ap.add_argument(
"--out-dir",
default=os.path.join(_REPO, "tests", "output", "sah04_case9_flow_diag"),
help="Output directory for PNG + JSON",
)
ap.add_argument(
"--x-list",
type=str,
default="",
help='Comma-separated x indices (default: built-in list); e.g. "1,2,5,40,200"',
)
args = ap.parse_args()
mc = next(c for c in sah.MATRIX if c.case_id == 9)
tier = sah._TIERS[mc.tier]
ny = int(tier["ny"])
nx = int(sah._NX)
center = (sah._CX, float(tier["center_y"]))
u_max = 0.1
u0_mean = u_max / 1.5
nu = u_max * float(sah._D) / float(mc.re)
if args.x_list.strip():
x_indices = [int(s.strip()) for s in args.x_list.split(",") if s.strip()]
x_indices = [x for x in x_indices if 1 <= x < nx - 1]
else:
x_indices = _default_x_indices(nx)
lbm_path = sah._DEFAULT_LBM
cfg = sah._load_json(lbm_path)
cfg["grid"]["nx"] = nx
cfg["grid"]["ny"] = ny
cfg["grid"]["nz"] = 1
cfg["physics"]["viscosity"] = float(nu)
cfg["physics"]["velocity"] = float(u0_mean)
cfg["physics"]["rho"] = 1.0
cfg["method"]["collision"] = str(args.collision).upper()
cfg["method"]["streaming"] = "double_buffer"
cfg["method"]["les"]["enabled"] = False
cfg["method"]["outlet"]["mode"] = args.outlet
body_doc = {
"objects": [
{
"type": "cylinder",
"center": list(center),
"radius": float(sah._R_CYL),
}
]
}
tmpd = tempfile.mkdtemp(prefix="celeris_sah09_diag_")
lbm_tmp = os.path.join(tmpd, "config_lbm.json")
body_tmp = os.path.join(tmpd, "config_body.json")
sah._write_json(lbm_tmp, cfg)
sah._write_json(body_tmp, body_doc)
os.makedirs(args.out_dir, exist_ok=True)
png_profiles = os.path.join(args.out_dir, "case9_profiles_ux_uy.png")
png_vort = os.path.join(args.out_dir, "case9_omega_z_full.png")
json_path = os.path.join(args.out_dir, "case9_profiles.json")
try:
from CelerisLab import Simulation # noqa: WPS433
sim = Simulation(lbm_config_path=lbm_tmp, body_config_path=body_tmp)
sim.initialize()
stream = cuda.Stream()
rec_every = max(1, int(args.record_every))
steps = int(args.steps)
for step in range(1, steps + 1):
sim.bodies.zero_force_segment_async(stream)
sim.stepper.step(
1,
action_gpu=sim.bodies.action_gpu,
obs_gpu=sim.bodies.obs_gpu,
stream=stream,
)
stream.synchronize()
macro = sim.get_macroscopic()
ux = np.asarray(macro["ux"], dtype=np.float64)
uy = np.asarray(macro["uy"], dtype=np.float64)
sim.close()
y = np.arange(ny, dtype=np.float64)
u_ref = _inlet_parabolic_ref_u(y, ny, u0_mean)
title = (
f"Sah04 case 9 {args.collision} Re={mc.re} ny={ny} "
f"steps={steps} (last step)"
)
_plot_profiles(png_profiles, ux=ux, uy=uy, x_indices=x_indices, ny=ny, u0_mean=u0_mean, title=title)
sah.save_final_vorticity_png(png_vort, macro["ux"], macro["uy"], title=title + " omega_z")
profiles: Dict[str, Any] = {
"case_id": 9,
"collision": str(args.collision).upper(),
"Re": mc.re,
"nx": nx,
"ny": ny,
"steps": steps,
"burn_requested": int(args.burn),
"u_max": u_max,
"u0_mean_parabolic": u0_mean,
"x_indices": list(x_indices),
"y_lattice": y.tolist(),
"inlet_parabolic_ref_ux": u_ref.tolist(),
"profiles": {},
}
for xi in x_indices:
profiles["profiles"][str(xi)] = {
"ux": ux[:, xi].tolist(),
"uy": uy[:, xi].tolist(),
}
with open(json_path, "w", encoding="utf-8") as f:
json.dump(profiles, f, indent=2)
print("Wrote:", png_profiles, flush=True)
print("Wrote:", png_vort, flush=True)
print("Wrote:", json_path, flush=True)
finally:
shutil.rmtree(tmpd, ignore_errors=True)
return 0
if __name__ == "__main__":
raise SystemExit(main())

View File

@ -1,452 +0,0 @@
# tests/run_sah04_case9_grid_blockage_compare.py
"""Sah04 matrix case 9: compare inlet-near flow vs grid refinement vs blockage.
Runs three D2Q9 setups anchored to case 9 (high tier, Re=200, ``u_max``=0.1):
1. **baseline** matrix geometry: ``D=30``, ``nx=80*D+2``, ``ny=35``,
cylinder ``(40*D+0.5, 17)``, ``r=D/2``.
2. **grid_2x** double lattice resolution with the same **blockage** ``D/H``
and the same ``Lx/D=80`` convention: ``D=60``, ``nx=80*D+2``, ``ny=68``,
``center_y`` doubled with the channel, ``r=D/2``, ``Re`` still based on ``D``.
3. **radius_half** same channel height as baseline but **cylinder diameter
halved** in lattice units: ``D=15``, ``nx=80*D+2``, ``ny=35``, same
``center_y``, ``r=D/2``, ``Re`` based on the smaller ``D`` (lower blockage).
Outputs under ``tests/output/sah04_case9_compare/`` by default:
- ``compare_meta.json`` per-variant grid, ``nu``, estimated ``beta=D/H`` (no
full arrays; see ``compare_fields.npz``).
- ``compare_fields.npz`` full ``ux`` / ``uy`` per variant.
- ``ux_vs_y_inlet.png`` ``u_x`` vs normalized wall-normal ``eta`` at several
``x`` indices (all variants overlaid per panel).
- ``omega_z_<variant>.png`` final-step vorticity ``omega_z = d u_y/dx - d u_x/dy``
for each variant.
Forward schedule: ``burn`` warm-up steps (may be ``0``), then ``steps`` more
steps; the reported field is the state **after** ``burn + steps`` LBM steps.
Usage::
conda run -n pycuda_3_10 python tests/run_sah04_case9_grid_blockage_compare.py --smoke
conda run -n pycuda_3_10 python tests/run_sah04_case9_grid_blockage_compare.py \\
--steps 20000 --burn 0 --collision MRT
Design::
Isolates whether an inlet-channel artifact scales with **mesh** (grid_2x)
or with **blockage** (radius_half) while keeping the Sah04-style confined
channel layout. Requires **matplotlib** for PNG output.
"""
from __future__ import annotations
import argparse
import json
import os
import sys
import tempfile
from dataclasses import dataclass
from typing import Any, Dict, List, Tuple
import numpy as np
import pycuda.driver as cuda
_REPO = os.path.abspath(os.path.join(os.path.dirname(__file__), ".."))
_DEFAULT_LBM = os.path.join(
_REPO, "src", "CelerisLab", "configs", "config_lbm.json",
)
def _load_json(path: str) -> dict:
with open(path, "r", encoding="utf-8") as f:
return json.load(f)
def _write_json(path: str, d: dict) -> None:
with open(path, "w", encoding="utf-8") as f:
json.dump(d, f, indent=2)
@dataclass(frozen=True)
class Variant:
"""One lattice-resolved Sah04-like case (case 9 Re / u_max convention)."""
key: str
nx: int
ny: int
center: Tuple[float, float]
d_lattice: float
r_cyl: float
re: float
note: str
def _variants() -> Tuple[Variant, Variant, Variant]:
d0 = 30
ny0 = 35
cy0 = 17.0
v0 = Variant(
key="baseline",
nx=80 * d0 + 2,
ny=ny0,
center=(40.0 * d0 + 0.5, cy0),
d_lattice=float(d0),
r_cyl=0.5 * d0,
re=200.0,
note="Matrix case 9 (high tier): D=30, ny=35.",
)
d2 = 60
ny2 = 2 * ny0 - 2
v1 = Variant(
key="grid_2x",
nx=80 * d2 + 2,
ny=ny2,
center=(40.0 * d2 + 0.5, 2.0 * cy0),
d_lattice=float(d2),
r_cyl=0.5 * d2,
re=200.0,
note="Double D and ny-2; Lx/D=80 unchanged; blockage D/H ~ baseline.",
)
dh = 15
v2 = Variant(
key="radius_half",
nx=80 * dh + 2,
ny=ny0,
center=(40.0 * dh + 0.5, cy0),
d_lattice=float(dh),
r_cyl=0.5 * dh,
re=200.0,
note="Half cylinder diameter in lattice units; same ny as baseline.",
)
return v0, v1, v2
def _vorticity_z_from_velocity(ux: np.ndarray, uy: np.ndarray) -> np.ndarray:
"""``omega_z = d u_y/dx - d u_x/dy`` on a 2D ``(ny, nx)`` slice (unit lattice spacing)."""
ux = np.asarray(ux, dtype=np.float64)
uy = np.asarray(uy, dtype=np.float64)
duy_dx = np.gradient(uy, axis=1)
dux_dy = np.gradient(ux, axis=0)
return duy_dx - dux_dy
def _save_vorticity_png(
path: str,
ux: np.ndarray,
uy: np.ndarray,
*,
title: str,
) -> None:
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
omega = _vorticity_z_from_velocity(ux, uy)
abs_o = np.abs(omega[np.isfinite(omega)])
if abs_o.size:
vmax = float(np.percentile(abs_o, 99.5))
if vmax <= 0.0:
vmax = float(np.max(abs_o)) or 1.0
else:
vmax = 1.0
ny, nx = omega.shape
fw = min(18.0, max(8.0, nx / 100.0))
fh = min(12.0, max(3.0, ny / 40.0))
fig, ax = plt.subplots(figsize=(fw, fh))
im = ax.imshow(
omega,
origin="lower",
aspect="equal",
cmap="RdBu_r",
vmin=-vmax,
vmax=vmax,
extent=(0, nx - 1, 0, ny - 1),
)
ax.set_xlabel("x (lattice)")
ax.set_ylabel("y (lattice)")
ax.set_title(title)
fig.colorbar(im, ax=ax, fraction=0.046, pad=0.04, label="omega_z")
fig.tight_layout()
os.makedirs(os.path.dirname(os.path.abspath(path)) or ".", exist_ok=True)
fig.savefig(path, dpi=150, bbox_inches="tight")
plt.close(fig)
def _run_variant(
v: Variant,
*,
collision: str,
outlet: str,
u_max: float,
burn: int,
steps: int,
) -> Dict[str, Any]:
u0_mean = u_max / 1.5
nu = u_max * float(v.d_lattice) / float(v.re)
if not os.path.isfile(_DEFAULT_LBM):
raise FileNotFoundError(_DEFAULT_LBM)
cfg = _load_json(_DEFAULT_LBM)
cfg["grid"]["nx"] = int(v.nx)
cfg["grid"]["ny"] = int(v.ny)
cfg["grid"]["nz"] = 1
cfg["physics"]["viscosity"] = float(nu)
cfg["physics"]["velocity"] = float(u0_mean)
cfg["physics"]["rho"] = 1.0
cfg["method"]["collision"] = collision
cfg["method"]["streaming"] = "double_buffer"
cfg["method"]["les"]["enabled"] = False
cfg["method"]["outlet"]["mode"] = outlet
body_doc = {
"objects": [
{
"type": "cylinder",
"center": list(v.center),
"radius": float(v.r_cyl),
}
]
}
tmpd = tempfile.mkdtemp(prefix="celeris_sah09cmp_")
lbm_tmp = os.path.join(tmpd, "config_lbm.json")
body_tmp = os.path.join(tmpd, "config_body.json")
_write_json(lbm_tmp, cfg)
_write_json(body_tmp, body_doc)
from CelerisLab import Simulation # noqa: WPS433
sim = Simulation(lbm_config_path=lbm_tmp, body_config_path=body_tmp)
sim.initialize()
stream = cuda.Stream()
total = int(burn) + int(steps)
if total < 1:
sim.close()
raise ValueError("burn + steps must be >= 1")
for _step in range(1, total + 1):
sim.bodies.zero_force_segment_async(stream)
sim.stepper.step(
1,
action_gpu=sim.bodies.action_gpu,
obs_gpu=sim.bodies.obs_gpu,
stream=stream,
)
stream.synchronize()
macro = sim.get_macroscopic()
ux = np.asarray(macro["ux"], dtype=np.float64).reshape(v.ny, v.nx)
uy = np.asarray(macro["uy"], dtype=np.float64).reshape(v.ny, v.nx)
rho = np.asarray(macro["rho"], dtype=np.float64).reshape(v.ny, v.nx)
sim.close()
h_fluid = max(int(v.ny) - 2, 1)
beta_est = float(v.d_lattice) / float(h_fluid)
return {
"key": v.key,
"nx": int(v.nx),
"ny": int(v.ny),
"center": [float(v.center[0]), float(v.center[1])],
"d_lattice": float(v.d_lattice),
"r_cyl": float(v.r_cyl),
"re": float(v.re),
"nu": float(nu),
"u_max": float(u_max),
"u0_mean": float(u0_mean),
"beta_est": beta_est,
"note": v.note,
"burn": int(burn),
"steps": int(steps),
"total_lbm_steps": int(total),
"ux_shape": list(ux.shape),
"rho_min": float(np.min(rho)),
"rho_max": float(np.max(rho)),
"ux": ux.tolist(),
"uy": uy.tolist(),
}
def _eta_coords(ny: int) -> np.ndarray:
"""Wall-normal fractional coordinate in (0,1) for fluid rows y=1..ny-2."""
y = np.arange(ny, dtype=np.float64)
h = max(float(ny - 2), 1.0)
out = np.zeros_like(y)
for yi in range(1, max(ny - 1, 1)):
out[yi] = (float(yi) - 0.5) / h
return out
def _plot_ux_vs_y_inlet(
path: str,
results: List[Dict[str, Any]],
x_indices: List[int],
) -> None:
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
n_x = len(x_indices)
fig, axes = plt.subplots(1, n_x, figsize=(4.2 * n_x, 5.0), sharey=True)
if n_x == 1:
axes = [axes]
colors = {"baseline": "C0", "grid_2x": "C1", "radius_half": "C2"}
for ax, xi in zip(axes, x_indices):
for r in results:
nx = int(r["nx"])
if xi < 1 or xi >= nx - 1:
continue
ux = np.asarray(r["ux"], dtype=np.float64).reshape(r["ny"], r["nx"])
eta = _eta_coords(int(r["ny"]))
y = np.arange(int(r["ny"]))
ax.plot(ux[:, xi], eta, label=r["key"], color=colors.get(r["key"], "k"))
ax.set_title(f"x = {xi}")
ax.set_xlabel("u_x")
ax.grid(True, alpha=0.3)
axes[0].set_ylabel("eta = (y-0.5)/H (fluid band)")
fig.suptitle("u_x vs wall-normal coordinate near inlet (several x)")
axes[0].set_ylim(0.0, 1.0)
axes[0].legend(loc="best", fontsize=8)
fig.tight_layout()
os.makedirs(os.path.dirname(os.path.abspath(path)) or ".", exist_ok=True)
fig.savefig(path, dpi=150, bbox_inches="tight")
plt.close(fig)
def main() -> int:
ap = argparse.ArgumentParser(
description="Case 9 baseline vs finer grid vs half-radius: inlet-near u_x",
)
ap.add_argument("--collision", default="MRT", help="SRT, TRT, or MRT")
ap.add_argument("--outlet", default="neq_extrap")
ap.add_argument("--u-max", type=float, default=0.1, dest="u_max")
ap.add_argument(
"--burn",
type=int,
default=0,
help="Warm-up LBM steps before the --steps production segment (0 is valid).",
)
ap.add_argument(
"--steps",
type=int,
default=20_000,
help="LBM steps after burn-in; final field is after burn+steps.",
)
ap.add_argument(
"--smoke",
action="store_true",
help="Short run: burn=0, steps=2500 (overrides --burn and --steps)",
)
ap.add_argument(
"--out-dir",
default=os.path.join(_REPO, "tests", "output", "sah04_case9_compare"),
)
ap.add_argument(
"--inlet-x",
type=str,
default="1,2,3,5,8",
help="Comma-separated x indices for u_x(y) panels (must exist for all nx)",
)
args = ap.parse_args()
if args.smoke:
burn = 0
steps = 2500
else:
burn = max(0, int(args.burn))
steps = max(1, int(args.steps))
out_dir = os.path.abspath(args.out_dir)
os.makedirs(out_dir, exist_ok=True)
v0, v1, v2 = _variants()
variants = (v0, v1, v2)
x_list = [int(s.strip()) for s in args.inlet_x.split(",") if s.strip()]
min_nx = min(v.nx for v in variants)
x_list = [x for x in x_list if 1 <= x < min_nx - 1]
if not x_list:
print("No valid --inlet-x indices common to all variants.", file=sys.stderr)
return 2
results: List[Dict[str, Any]] = []
for v in variants:
print(f"--- {v.key}: nx={v.nx} ny={v.ny} D={v.d_lattice} r={v.r_cyl} ---", flush=True)
row = _run_variant(
v,
collision=str(args.collision).upper(),
outlet=str(args.outlet),
u_max=float(args.u_max),
burn=burn,
steps=steps,
)
results.append(row)
meta_path = os.path.join(out_dir, "compare_meta.json")
slim = []
for r in results:
slim.append({k: v for k, v in r.items() if k not in ("ux", "uy")})
_write_json(
meta_path,
{
"variants_summary": slim,
"inlet_x_used": x_list,
"burn": burn,
"steps": steps,
"total_lbm_steps": burn + steps,
"collision": str(args.collision).upper(),
},
)
try:
_plot_ux_vs_y_inlet(
os.path.join(out_dir, "ux_vs_y_inlet.png"),
results,
x_list,
)
for r in results:
key = str(r["key"])
ux = np.asarray(r["ux"], dtype=np.float64).reshape(r["ny"], r["nx"])
uy = np.asarray(r["uy"], dtype=np.float64).reshape(r["ny"], r["nx"])
_save_vorticity_png(
os.path.join(out_dir, f"omega_z_{key}.png"),
ux,
uy,
title=f"omega_z final ({key}) nx={r['nx']} ny={r['ny']} "
f"D={r['d_lattice']} burn={r['burn']} steps={r['steps']}",
)
except ImportError:
print(
"matplotlib not installed; skipped PNG. Install matplotlib for figures.",
file=sys.stderr,
)
full_npz = os.path.join(out_dir, "compare_fields.npz")
np.savez_compressed(
full_npz,
keys=np.array([r["key"] for r in results]),
**{f"ux_{r['key']}": np.asarray(r["ux"], dtype=np.float64) for r in results},
**{f"uy_{r['key']}": np.asarray(r["uy"], dtype=np.float64) for r in results},
**{
f"omega_z_{r['key']}": _vorticity_z_from_velocity(
np.asarray(r["ux"], dtype=np.float64).reshape(r["ny"], r["nx"]),
np.asarray(r["uy"], dtype=np.float64).reshape(r["ny"], r["nx"]),
)
for r in results
},
)
print(f"Wrote: {meta_path}", flush=True)
print(f"Wrote: {full_npz}", flush=True)
print(f"Wrote: {os.path.join(out_dir, 'ux_vs_y_inlet.png')}", flush=True)
for r in results:
oz_path = os.path.join(out_dir, f"omega_z_{r['key']}.png")
print(f"Wrote: {oz_path}", flush=True)
return 0
if __name__ == "__main__":
raise SystemExit(main())

File diff suppressed because it is too large Load Diff

493
tests/streakline_notes.md Normal file
View File

@ -0,0 +1,493 @@
## Streakline postprocessing design
这份说明面向一个最小可实现的离线后处理程序。输入是按时间保存的二维速度场序列,输出是接近实验染色图的 streakline 图像。核心路线是连续释粒,因为 streakline 正是“固定位置持续释放染料后,在某个时刻所有已释放粒子的位置集合” [Lan96]。对非定常尾迹,这比单帧 streamline 更贴近实验图像 [Lan96]。
搜索和全文阅读给出的信息已经足够支撑第一版实现。最直接的方案是按保存时刻连续注入粒子,用时空插值后的速度场推进粒子,再将粒子云渲染为图像 [Lan96, Ken96]。如果纯粒子结果过于尖锐,再在推进后加入一个很小的随机扩散项,作为对真实染色液扩散的近似。更完整的被动标量法当然更物理,但实现成本更高,不适合作为第一版 [Kim04]。
## Core recommendation
| 方法 | 与实验染色图的对应 | 实现复杂度 | 适合作为第一版 |
|---|---|---:|---|
| 单帧 streamline | 弱 | 很低 | 不推荐 |
| 连续释粒 streakline | 强 | 低 | 最推荐 |
| streakline 加少量扩散 | 很强 | 低到中 | 推荐作为第二步 |
| 被动标量对流扩散 | 最强 | 中到高 | 暂不作为第一版 |
最小可实现路线如下:
- 从 CelerisLab 导出一串二维速度场快照 `u_x(x,y,t_k), u_y(x,y,t_k)`
- 选定一个固定释放点或一小段释放线
- 在每个保存时刻都注入新粒子
- 在两个相邻快照之间,用时间插值和空间插值计算粒子速度
- 用二阶或四阶时间积分推进粒子
- 删去出域粒子与进入固体的粒子
- 将当前时刻存活粒子按位置和年龄渲染成图像
这个程序逻辑直接对应 [Lan96] 的 streakline 算法,插值与推进细节由 [Ken96] 和 [Dar96] 支撑。
## Data contract
默认输入是规则 Cartesian 网格上的二维时间序列速度场。对当前项目,最自然的输入约定是:
| 名称 | 含义 |
|---|---|
| `t_k` | 第 `k` 个保存时刻 |
| `ux[k, j, i]` | 时刻 `t_k` 的 x 方向速度 |
| `uy[k, j, i]` | 时刻 `t_k` 的 y 方向速度 |
| `mask[j, i]` | 可选,固体与流体标记 |
| `x_i, y_j` | 网格坐标 |
第一版最好满足下面三条:
- 时间快照间隔固定或至少已知
- 网格坐标固定不动
- 固体几何位置已知,至少能判断粒子是否进入圆柱内部
如果输出来自 CelerisLab本质上只需要把每个保存时刻的 `ux, uy` 和对应时间写出来即可。算法本身是通用的,并不依赖 LBM 本身。
## Mathematical model for the base streakline
### Particle motion
粒子位置满足拉格朗日运动方程 [Lan96]
\[
\frac{d \boldsymbol{x}}{dt} = \boldsymbol{u}(\boldsymbol{x}, t)
\]
其中
\[
\boldsymbol{x} = (x, y), \qquad \boldsymbol{u} = (u_x, u_y)
\]
积分形式为 [Lan96]
\[
\boldsymbol{x}(t + \Delta t) = \boldsymbol{x}(t) + \int_t^{t+\Delta t} \boldsymbol{u}(\boldsymbol{x}(\tau), \tau) \, d\tau
\]
streakline 在观察时刻 \(t_n\) 的定义是:从固定释放位置 \(\boldsymbol{x}_s\) 在过去各时刻持续注入的粒子,在 \(t_n\) 时刻的全部位置集合 [Lan96]。
### Time interpolation
速度场只在离散时刻存储,所以必须做时间插值。对位于 \(t_k \le t \le t_{k+1}\) 的任意子步,最简单的做法是线性时间插值 [Ken96]
\[
\delta = \frac{t - t_k}{t_{k+1} - t_k}
\]
\[
\boldsymbol{u}(\boldsymbol{x}, t) = (1 - \delta) \, \boldsymbol{u}_k(\boldsymbol{x}) + \delta \, \boldsymbol{u}_{k+1}(\boldsymbol{x})
\]
这里 \(\boldsymbol{u}_k(\boldsymbol{x})\) 和 \(\boldsymbol{u}_{k+1}(\boldsymbol{x})\) 仍需通过空间插值得到。
### Spatial interpolation
在规则网格上,第一版直接用双线性插值即可。虽然 [Ken96] 讨论的是非结构网格上的点定位与线性插值,但其一般原则完全适用于规则网格:
- 先定位粒子所在单元
- 再用单元顶点速度插值得到粒子点速度 [Ken96]
若粒子位于单元局部坐标 \((\xi, \eta) \in [0,1]^2\),四个角点速度为 \(\boldsymbol{u}_{00}, \boldsymbol{u}_{10}, \boldsymbol{u}_{01}, \boldsymbol{u}_{11}\),则
\[
\boldsymbol{u}(\xi, \eta) = (1-\xi)(1-\eta) \boldsymbol{u}_{00}
+ \xi (1-\eta) \boldsymbol{u}_{10}
+ (1-\xi)\eta \boldsymbol{u}_{01}
+ \xi\eta \boldsymbol{u}_{11}
\]
## Time integration choice
[Lan96] 和 [Ken96] 都直接使用四阶 Runge Kutta。对你的第一版程序推荐两个选项
| 积分器 | 优点 | 缺点 | 建议 |
|---|---|---|---|
| RK2 | 简洁,容易调试 | 精度一般 | 最小原型可用 |
| RK4 | 文献最一致,精度更稳 | 每步插值次数更多 | 默认推荐 |
RK4 的更新式为 [Lan96, Ken96]
\[
\boldsymbol{x}_{n+1} = \boldsymbol{x}_n + \frac{1}{6}(\boldsymbol{a} + 2\boldsymbol{b} + 2\boldsymbol{c} + \boldsymbol{d})
\]
其中
\[
\boldsymbol{a} = \Delta t \, \boldsymbol{u}(\boldsymbol{x}_n, t_n)
\]
\[
\boldsymbol{b} = \Delta t \, \boldsymbol{u}(\boldsymbol{x}_n + \tfrac{1}{2}\boldsymbol{a}, t_n + \tfrac{1}{2}\Delta t)
\]
\[
\boldsymbol{c} = \Delta t \, \boldsymbol{u}(\boldsymbol{x}_n + \tfrac{1}{2}\boldsymbol{b}, t_n + \tfrac{1}{2}\Delta t)
\]
\[
\boldsymbol{d} = \Delta t \, \boldsymbol{u}(\boldsymbol{x}_n + \boldsymbol{c}, t_n + \Delta t)
\]
[Ken96] 特别强调RK4 的每一个子步都要重新做点定位与时空插值。这在规则网格上不难实现,只是意味着一次完整步进需要四次速度查询 [Ken96]。
## Practical timestep guidance
第一版不必上自适应步长,但步长不能随意设大。搜索里最重要的精度提醒来自 [Dar96]
- 非定常粒子积分的误差常常由时间离散控制,而不是概念本身
- 合理步长必须不大于流动主要非定常时间尺度的同量级 [Dar96]
- 若时间步过大,粒子轨迹会在拓扑上都出错,而不仅是位置略偏 [Dar96]
对当前项目,最简单的工程准则是:
\[
\Delta t_{trace} \le \min \bigl(\alpha_t \, \Delta t_{save},\; \alpha_x \, \frac{\Delta x}{\max |u|} \bigr)
\]
其中 \(\Delta t_{save}\) 是两个保存快照之间的时间间隔,\(\Delta x\) 是网格尺度。第一版可取:
- \(\alpha_t = 0.1 \sim 0.25\)
- \(\alpha_x = 0.25 \sim 0.5\)
这不是文献中的严格上界,而是结合 [Dar96] 的结论给出的实现准则:粒子推进子步既要显著小于保存间隔,也不要一子步跨过太多网格。
[Ken96] 还给了一个很实用的自适应思路:根据相邻速度方向夹角调节步长。如果速度方向变化太快就减半,变化很小就加倍 [Ken96]。这很适合后续增强版,但不属于第一版必需项。
## Base streakline algorithm
### Release strategy
最符合染色实验的是连续释粒 [Lan96]。第一版可用两种释放方式:
| 方式 | 适用场景 | 建议 |
|---|---|---|
| 单点释放 | 针头式染料入口 | 最简单 |
| 短线段释放 | 更像细缝或细带注入 | 更稳健 |
如果实验是在圆柱上游某一点持续释放染色液,第一版就直接使用单点释放。
### Particle state
每个粒子只需保存:
| 字段 | 含义 |
|---|---|
| `x, y` | 当前坐标 |
| `t_birth` | 释放时刻 |
| `age` | 当前年龄 |
| `alive` | 是否仍在域内 |
第一版不需要保存整条历史轨迹。因为你要的是当前时刻的 streakline 图,而不是每个粒子的 pathline 曲线。
**常见误区(务必区分)**
| 概念 | 画什么 | 是否等于水洞染色迹线 |
|---|---|---|
| **Streakline迹线** | 固定源点**持续释粒**,在观察时刻 \(t_n\) 画出**所有仍存活粒子的当前位置**(可带年龄衰减权重) | 是 |
| **Pathline迹线/轨道)** | 单个粒子从释放到当前的**整条历史轨迹** | 否 |
| **错误实现** | 每隔一段时间放一个粒子,再把**所有粒子走过的路径段**全部叠加到图上 | 否(这是 pathline 叠加,不是 streakline |
水洞实验:针头在固定点连续注 dye → 某一时刻拍照 → 看到的是“此刻染料粒子在流场里的分布”,一条下游色带由**不同释放时刻、当前仍在本流场中的粒子**共同构成,而不是把每个粒子从出生到现在的轨迹都画出来。
### Main loop pseudocode
```text
given velocity snapshots U[k] = {ux[k], uy[k]} at times t[k]
given seeding point or seeding segment S
initialize empty particle list P
for k = 0 to N-2:
inject new particles at source S at time t[k]
set t_local = t[k]
while t_local < t[k+1]:
dt = min(dt_trace, t[k+1] - t_local)
for each particle p in P with p.alive:
v = interpolate_velocity(p.position, t_local)
advance p by one integration step using RK2 or RK4
update p.age
if p leaves domain:
p.alive = false
if p enters solid body:
p.alive = false
t_local = t_local + dt
optionally remove very old particles
optionally render current particle cloud
```
这个结构与 [Lan96] 的离散 streakline 算法一致,只是把“从 `t_k``t_{k+1}` 的一次推进”细化成多个更小的粒子子步,以满足 [Dar96] 对精度的要求。
### Velocity query pseudocode
```text
function interpolate_velocity(position x, time t):
find k such that t[k] <= t <= t[k+1]
compute delta = (t - t[k]) / (t[k+1] - t[k])
u_k = bilinear_interpolation(U[k], x)
u_k1 = bilinear_interpolation(U[k+1], x)
return (1 - delta) * u_k + delta * u_k1
```
这是第一版最核心的数值部件。只要这部分实现正确streakline 程序主体就很直接。
## Small diffusion particle model
当纯 streakline 太细、太锐利、不像染色液图像时,可以给每个粒子加一个很小的随机扩散项。这不是严格的被动标量求解,但能以很低代价增加染色带宽度。
### Guiding idea
[Kim04] 表明,染料图像更接近一个被动标量浓度场,其基本控制方程是对流扩散方程:
\[
\frac{\partial c}{\partial t} + \boldsymbol{u} \cdot \nabla c = D \nabla^2 c
\]
其中 \(c\) 是染料浓度,\(D\) 是扩散系数。对第一版粒子法,一个常见近似是将每个粒子的位置更新写成“对流加随机扩散”:
\[
\boldsymbol{x}_{n+1} = \boldsymbol{x}_n + \Delta \boldsymbol{x}_{adv} + \Delta \boldsymbol{x}_{diff}
\]
其中 \(\Delta \boldsymbol{x}_{adv}\) 由 RK2 或 RK4 给出,扩散项取二维各向同性随机增量:
\[
\Delta \boldsymbol{x}_{diff} = \sqrt{2 D \Delta t}
\begin{bmatrix}
\eta_x \\
\eta_y
\end{bmatrix}
\]
这里 \(\eta_x, \eta_y \sim \mathcal{N}(0,1)\)。
这个形式本身是对扩散过程的标准随机游走近似。搜索结果提醒,随机游走模型如果处理不当会产生假漂移与错误浓度偏置 [Mac92]。因此第一版的使用原则应当很克制:
- 只用于加入少量模糊和厚度
- 不把粒子密度当成严格浓度
- 不在强近壁统计上过度解读结果
### Diffusive particle update pseudocode
```text
for each particle p in P with p.alive:
x_adv = RK4_step(p.position, t_local, dt)
sigma = sqrt(2 * D * dt)
dx_rand = sigma * normal(0, 1)
dy_rand = sigma * normal(0, 1)
x_new = x_adv + [dx_rand, dy_rand]
if x_new leaves domain:
p.alive = false
else if x_new enters solid:
p.alive = false
else:
p.position = x_new
p.age += dt
```
### Choosing the diffusion level
第一版不必试图从真实染料物性严格标定 \(D\)。更实用的做法是把 \(D\) 当作视觉匹配参数,并保持它足够小,使图像结构仍主要由对流控制。
[Kim04] 用 Schmidt 数控制扩散强弱。若已有参考速度 \(U\) 和长度尺度 \(L\),则
\[
Sc = \frac{\nu}{D}
\]
也可写成
\[
D = \frac{\nu}{Sc}
\]
对第一版,可以用下面的思路选扩散强度:
| 目标效果 | 建议 |
|---|---|
| 只想让线条略微变厚 | 取较大 `Sc`,即很小的 `D` |
| 想模拟明显洗开和模糊 | 取较小 `Sc`,即较大的 `D` |
| 不确定 | 先从几乎看不出的弱扩散开始 |
因为当前目标是“少量扩散”,所以推荐先把扩散当成弱修饰,而不是主导机制。
## Rendering logic
最终图像不需要把每个粒子的完整轨迹都画出来。更像实验染色图的做法是把当前存活粒子投影到图像网格上,生成粒子密度图或带年龄权重的强度图。
### Minimal rendering choices
| 方法 | 图像风格 | 实现难度 |
|---|---|---:|
| 直接散点 | 最简陋 | 很低 |
| 网格计数直方图 | 像浓度图 | 低 |
| 高斯核累积 | 更平滑 | 低到中 |
第一版推荐:
- 把每个粒子投到像素网格
- 对像素做计数或加权累积
- 最后做一次轻微 Gaussian blur
这与先做完整被动标量相比便宜很多,但视觉上已经会很接近实验染色图。
### Optional particle weighting
可以给粒子一个简单权重:
\[
I = \sum_p w_p K(\boldsymbol{x} - \boldsymbol{x}_p)
\]
其中 \(K\) 是像素核或 Gaussian 核。第一版里,权重 \(w_p\) 可直接取 1也可对年龄做衰减例如
\[
w_p = \exp\left(- \frac{\mathrm{age}_p}{\tau_f} \right)
\]
这样旧粒子会逐渐淡出,图像不会无限堆积。
## Boundary handling
这部分先按最小原则处理即可。
| 情况 | 第一版处理 |
|---|---|
| 粒子出计算域 | 直接删除 |
| 粒子进入圆柱内部 | 直接删除 |
| 粒子贴近边界滑动 | 暂不专门处理 |
| 扩散后跨入固体 | 直接删除 |
这样做的优点是简单稳妥。后续若发现近壁 streakline 误差明显,再考虑反射、投影或更物理的壁面处理。
## Minimal implementation plan
### Version 1
目标是尽快得到可信的 streakline 图:
- 读取 `ux, uy, t`
- 单点连续释粒
- 双线性空间插值
- 线性时间插值
- RK4 粒子推进
- 出域删除与入固体删除
- 粒子计数成图
### Version 2
在不改变主体架构的前提下增强视觉效果:
- 增加短线段释放
- 增加年龄衰减
- 增加弱随机扩散
- 用 Gaussian 核代替简单计数
### Version 3
若后续发现实验图像明显受扩散与混合主导,再转向更重的模型:
- 被动标量对流扩散
- 用 \(Sc\) 或 \(D\) 做参数标定 [Kim04]
## Design decisions that are already justified by the literature
下面这些设计已经有足够文献支撑,可以直接采用:
- 用 streakline 而不是 streamline 对应持续染色实验 [Lan96]
- 用连续释放粒子离线重构 streakline [Lan96]
- 用时空插值从离散速度快照查询粒子速度 [Ken96]
- 用 RK4 作为默认推进器 [Lan96, Ken96]
- 把步长取得明显小于保存时间间隔,并避免一次跨过太多网格 [Dar96]
- 若需要更像染料图,可以加入弱扩散,或后续上升到被动标量 [Kim04]
## Recommended default choices
| 项目 | 默认选择 |
|---|---|
| 释放方式 | 单点连续释粒 |
| 空间插值 | 双线性 |
| 时间插值 | 线性 |
| 积分器 | RK4 |
| 粒子子步 | 保存步长的 0.1 到 0.25 |
| 固体处理 | 入固体即删除 |
| 输出图像 | 粒子密度图加轻微模糊 |
| 扩散 | 默认关闭,作为第二步 |
## Bottom line
对当前项目,最简洁且足够准确的程序逻辑不是去解新的染料场,而是先做一个离线连续释粒 streakline 后处理器。它只依赖时间序列速度场,程序结构清楚,数值风险也集中在可控的几个环节:时空插值、积分步长和边界删除 [Lan96, Dar96, Ken96]。在此基础上,再增加一个弱随机扩散项,就能以很小代价把图像从“几何上正确的粒子线”推进到“更像实验染色照片”的粒子云 [Kim04]。
## CelerisLab implementation mapping
Current implementation lives in:
- library: `src/CelerisLab/common/streakline.py`
- Kan99b demo CLI: `tests/run_kan99b_streakline.py`
- experiment notebook: `tests/experiment.ipynb` (CLEAN 5/6)
Available runtime modes:
- `online`: sample `Simulation.get_macroscopic()` directly in memory (no velocity snapshot files).
- `offline`: replay from an existing snapshot directory and run the same streakline integrator.
Dense release support:
- base upstream points are expanded by `release_mode` (`point|line|strip`).
- `strip` mode densifies both cross-stream (`line_count`, `line_span`) and downstream (`downstream_count`, `downstream_spacing`) directions.
### Performance notes (why it can feel slow)
Streakline post-processing is **CPU-side** in the current stack:
1. **GPU→host velocity copies** dominate for large grids (`nx=6000, ny=1200`). Each `get_macroscopic()` pulls full `ux/uy` arrays.
2. Particle RK4 + bilinear interpolation runs in NumPy on CPU (vectorized over particles, not GPU).
3. Rendering accumulates polylines into an image (optionally multi-threaded via `ThreadPoolExecutor`).
Practical tuning:
| knob | effect |
|---|---|
| `sample_every` | larger → fewer host copies, faster, coarser streaks |
| `max_particle_age` | `None` → no hard cutoff; trails survive until boundary/solid |
| `num_threads` | `0` = auto; sets OpenBLAS/OMP threads for blur/host ops |
| `blur_sigma` | `0` disables Gaussian blur pass |
For `experiment.ipynb` triangle case, use **`render_streakline_density(..., minimal_axes=True)`** on the **final alive particle cloud** (`positions`, `ages`). Do **not** accumulate and draw full path history (`render_snapshot_trails` is pathline-style debug only).
### Clean render style (experiment)
`render_snapshot_trails` defaults:
- white background
- red streaks, brighter toward the downstream end (`fade_along_trail=True`)
- black filled cylinders
- no axes, labels, colorbar, or release-point markers
Example commands:
```bash
# Online in-memory streakline (no snapshot write)
conda run -n pycuda_3_10 python tests/run_kan99b_streakline.py \
--mode online --domain M --re 100 --alpha 1.0 \
--sample-every 300 --n-snapshots 20 \
--release-mode strip --line-count 7 --downstream-count 6
# Offline replay from existing velocity snapshots
conda run -n pycuda_3_10 python tests/run_kan99b_streakline.py \
--mode offline \
--snapshot-dir tests/output/final_validation_round/streakline_kan99b_k2/velocity_snapshots \
--release-mode strip --line-count 7 --downstream-count 6
```

Binary file not shown.

After

Width:  |  Height:  |  Size: 316 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 134 KiB