1.8 MiB
1.8 MiB
In [1]:
# ==================== CLEAN 1. Imports ====================
from pathlib import Path
import math
import json
import tempfile
import numpy as np
import matplotlib.pyplot as plt
import pycuda.driver as cuda
from CelerisLab import Simulation
from CelerisLab.common.streakline import Streakline, ReleaseConfig, IntegratorConfig
from CelerisLab.common.render import compute_vorticity, render_vorticity_field
from CelerisLab.common.preprocess import cylinders_from_triangle_layout
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 # experiment positive(clockwise) -> kernel omega(CCW positive)
In [2]:
# ==================== CLEAN 2. Parameters (edit here) ====================
cwd = Path.cwd()
WORKDIR = (cwd / "tests") if (cwd / "src").exists() else cwd
WORKDIR.mkdir(parents=True, exist_ok=True)
CFD_PARAMS = {
"config_path": str((cwd / "src/CelerisLab/configs/config_lbm_three_cylinder_triangle.json") if (cwd / "src").exists() else Path("../src/CelerisLab/configs/config_lbm_three_cylinder_triangle.json")),
"inlet_scheme": "regularized",
"checkpoint_path": str(WORKDIR / "triangle_state_clean.h5"),
"output_streak_png": str(WORKDIR / "streakline_triangle_clean.png"),
"output_streak_meta": str(WORKDIR / "streakline_triangle_clean.json"),
"output_vort_png": str(WORKDIR / "vorticity_triangle_clean.png"),
"report_every": 10000,
}
RUN_PARAMS = {
"steady": False,
"steps_this_chunk": 0, # 0 => 2 * nx / (3 * U_lb)
"load_checkpoint": False,
}
STREAK_PARAMS = {
"sample_every": 200,
"inject_per_seed": 1,
# Per seed: 5 particles along +x with small spacing (strip mode).
"release_mode": "strip",
"downstream_count": 5,
"downstream_spacing": 1.0,
"max_particle_age": None,
"age_decay_steps": 80000.0,
"blur_sigma": 0.8,
"num_threads": 0,
}
ACTION_FEATURES = {
"action1": {
"mean": 0.0,
"components": [
# (0.2708, 0.0044, -2.6190),
# (0.8124, 0.0012, 2.8560),
(0.2708, 0.000, -2.6190),
],
},
"action2": {
# "mean": -0.0140,
"mean": -0.0,
"components": [
# (0.2708, 0.0075, -0.9330),
# (0.8124, 0.0015, 0.8010),
# (0.5418, 0.0005, 1.8540),
(0.2708, 0.000, -0.9330),
],
},
"action3": {
# "mean": 0.0140,
"mean": 0.0,
"components": [
# (0.2708, 0.0075, -1.3980),
# (0.8124, 0.0015, 2.2080),
# (0.5418, 0.0005, 1.8100),
(0.2708, 0.000, -1.3980),
],
},
}
print("WORKDIR:", WORKDIR)
print("CONFIG:", CFD_PARAMS["config_path"])
WORKDIR: /home/frank14f/CelerisLab/tests CONFIG: ../src/CelerisLab/configs/config_lbm_three_cylinder_triangle.json
In [3]:
# ==================== CLEAN 3. Helpers ====================
def generate_signal(t_phys, feature):
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):
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_m_s_at_time(t_phys, steady=False):
if steady:
return (
float(ACTION_FEATURES["action1"]["mean"]),
float(ACTION_FEATURES["action2"]["mean"]),
float(ACTION_FEATURES["action3"]["mean"]),
)
s1 = generate_signal(t_phys, ACTION_FEATURES["action1"])
s2 = generate_signal(t_phys, ACTION_FEATURES["action2"])
s3 = generate_signal(t_phys, ACTION_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, u_lb):
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, omega0, omega1, omega2):
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 ensure_compat_config(config_path, preferred_scheme="regularized"):
config_file = Path(config_path)
with config_file.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_file)
tmp = tempfile.NamedTemporaryFile(
mode="w",
suffix="_compat_config_lbm.json",
delete=False,
encoding="utf-8",
)
with tmp:
json.dump(cfg, tmp, indent=4)
return tmp.name
def triangle_layout(cfg):
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
y_upper = y_center + 0.5 * spacing_lb
y_lower = y_center - 0.5 * spacing_lb
return {
"x_apex": x_apex,
"x_rear": x_rear,
"y_center": y_center,
"y_upper": y_upper,
"y_lower": y_lower,
"radius_lb": radius_lb,
}
def add_triangle_cylinders(sim):
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 build_release_points_for_triangle(layout):
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.float32)
return np.column_stack([np.full(4, release_x, dtype=np.float32), ys]).astype(np.float32)
def filter_particles_outside_cylinders(points, layout):
centers = np.array(
[
[layout["x_apex"], layout["y_center"]],
[layout["x_rear"], layout["y_lower"]],
[layout["x_rear"], layout["y_upper"]],
],
dtype=np.float32,
)
dist2_min = np.full(points.shape[0], np.inf, dtype=np.float32)
for center in centers:
d2 = (points[:, 0] - center[0]) ** 2 + (points[:, 1] - center[1]) ** 2
dist2_min = np.minimum(dist2_min, d2)
return dist2_min > (layout["radius_lb"] ** 2)
In [4]:
# ==================== CLEAN 4. Init/restore session ====================
if "CLEAN_STATE" not in globals():
CLEAN_STATE = {"sim": None, "steps_done": 0, "u_lb": None, "dt_phys": None}
if CLEAN_STATE["sim"] is None:
compat_config_path = ensure_compat_config(
CFD_PARAMS["config_path"], preferred_scheme=CFD_PARAMS["inlet_scheme"]
)
sim = Simulation(compat_config_path)
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)
CLEAN_STATE["sim"] = sim
CLEAN_STATE["layout"] = layout
CLEAN_STATE["u_lb"] = u_lb
CLEAN_STATE["dt_phys"] = dt_phys
CLEAN_STATE["steps_done"] = 0
ckpt = Path(CFD_PARAMS["checkpoint_path"])
if RUN_PARAMS["load_checkpoint"] and ckpt.exists():
sim.load_checkpoint(str(ckpt))
CLEAN_STATE["steps_done"] = int(sim.stepper.step_count)
print(f"Restored checkpoint: {ckpt}, step={CLEAN_STATE['steps_done']}")
else:
print("Initialized from scratch.")
print(f"Using config: {compat_config_path}")
print(f"Inlet scheme: {sim.lbm_cfg.inlet_scheme}")
else:
print("Session already initialized.")
print(f"u_lb={CLEAN_STATE['u_lb']}, dt_phys={CLEAN_STATE['dt_phys']}")Initialized from scratch. u_lb=0.04, dt_phys=0.0022153300841825435
In [6]:
# ==================== CLEAN 5. Compute chunk only ====================
sim = CLEAN_STATE["sim"]
if sim is None:
raise RuntimeError("Run CLEAN 4 first.")
layout = CLEAN_STATE["layout"]
steps = int(RUN_PARAMS["steps_this_chunk"])
if steps <= 0:
steps = int(round(2.0 * sim.lbm_cfg.nx / (3.0 * CLEAN_STATE["u_lb"])))
steady = bool(RUN_PARAMS["steady"])
report_every = max(1, int(CFD_PARAMS["report_every"]))
sample_every = max(1, int(STREAK_PARAMS["sample_every"]))
n_threads = configure_compute_threads(int(STREAK_PARAMS.get("num_threads", 0)))
print(f"Host compute threads (BLAS/blur): {n_threads}")
release_cfg = ReleaseConfig(
mode=str(STREAK_PARAMS.get("release_mode", "strip")),
line_count=1,
line_span=0.0,
downstream_count=max(1, int(STREAK_PARAMS.get("downstream_count", 5))),
downstream_spacing=float(STREAK_PARAMS.get("downstream_spacing", 1.0)),
inject_per_seed=max(1, int(STREAK_PARAMS["inject_per_seed"])),
)
max_age = STREAK_PARAMS.get("max_particle_age")
integrator_cfg = IntegratorConfig(
alpha_t=0.25,
alpha_x=0.40,
max_particle_age=None if max_age is None else float(max_age),
num_threads=n_threads,
)
base_release = build_release_points_for_triangle(layout)
release_points = build_release_points(base_release, release_cfg).astype(np.float64)
cylinders = cylinders_from_triangle_layout(layout)
print(
f"Release seeds: {base_release.shape[0]} base -> {release_points.shape[0]} "
f"(+x x{release_cfg.downstream_count}, spacing={release_cfg.downstream_spacing})"
)
# Streakline = continuous release at fixed seeds + current particle cloud at end (not path history).
particles = np.empty((0, 2), dtype=np.float64)
ages = np.empty((0,), dtype=np.float64)
prev_frame = None
chunk_start_step = int(CLEAN_STATE["steps_done"])
for i in range(steps):
t_phys = CLEAN_STATE["steps_done"] * CLEAN_STATE["dt_phys"]
a1, a2, a3 = actions_m_s_at_time(t_phys, steady=steady)
w0 = action_to_omega_lb(a1, CLEAN_STATE["u_lb"])
w2 = action_to_omega_lb(a2, CLEAN_STATE["u_lb"])
w3 = action_to_omega_lb(a3, CLEAN_STATE["u_lb"])
# action2 -> rear-upper, action3 -> rear-lower (swap vs body add order)
set_body_omegas(sim, w0, w3, w2)
sim.run(1)
CLEAN_STATE["steps_done"] += 1
if CLEAN_STATE["steps_done"] % sample_every == 0:
macro = sim.get_macroscopic()
curr_frame = FlowFrame(
step=int(CLEAN_STATE["steps_done"]),
ux=macro["ux"].astype(np.float64, copy=False),
uy=macro["uy"].astype(np.float64, copy=False),
)
# Continuous dye injection: new particles always born at upstream seeds.
new_particles = np.repeat(release_points, release_cfg.inject_per_seed, axis=0).astype(np.float64)
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=sim.lbm_cfg.nx,
ny=sim.lbm_cfg.ny,
cfg=integrator_cfg,
cylinders=cylinders,
)
prev_frame = curr_frame
if (i + 1) % report_every == 0:
print(
f"chunk_step={i+1}/{steps} global={CLEAN_STATE['steps_done']} "
f"a=({a1:+.5f},{a2:+.5f},{a3:+.5f}) particles={particles.shape[0]}"
)
sim.save_checkpoint(CFD_PARAMS["checkpoint_path"])
if particles.size == 0:
raise RuntimeError("No alive particles at end. Increase steps or lower sample_every.")
streak_png = Path(CFD_PARAMS["output_streak_png"])
render_info = render_streakline_density(
particles,
ages,
nx=sim.lbm_cfg.nx,
ny=sim.lbm_cfg.ny,
out_path=str(streak_png),
cylinders=cylinders,
age_decay_steps=float(STREAK_PARAMS["age_decay_steps"]),
blur_sigma=float(STREAK_PARAMS["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,
)
streak_meta = {
"chunk_start_step": chunk_start_step,
"chunk_end_step": int(CLEAN_STATE["steps_done"]),
"sample_every": sample_every,
"release_points_base": base_release.tolist(),
"release_points_dense": release_points.tolist(),
"release_x_4D_upstream": float(base_release[0, 0]),
"particle_count": int(particles.shape[0]),
"max_particle_age": max_age,
"age_decay_steps": float(STREAK_PARAMS["age_decay_steps"]),
"render_mode": "particle_cloud_streakline",
"render": render_info,
}
meta_path = Path(CFD_PARAMS["output_streak_meta"])
meta_path.parent.mkdir(parents=True, exist_ok=True)
macro_end = sim.get_macroscopic()
ux_end = macro_end["ux"].astype(np.float32, copy=False)
uy_end = macro_end["uy"].astype(np.float32, copy=False)
vort_end = compute_vorticity(ux_end, uy_end)
vort_png = Path(CFD_PARAMS["output_vort_png"])
vort_info = render_vorticity_field(
vort_end,
nx=sim.lbm_cfg.nx,
ny=sim.lbm_cfg.ny,
out_path=str(vort_png),
cylinders=cylinders,
minimal_axes=True,
)
streak_meta["vorticity"] = vort_info
CLEAN_STATE["streak_particles"] = particles
CLEAN_STATE["streak_age"] = ages
CLEAN_STATE["streak_image"] = plt.imread(streak_png)
CLEAN_STATE["vort"] = vort_end
CLEAN_STATE["vort_image"] = plt.imread(vort_png)
with meta_path.open("w", encoding="utf-8") as f:
json.dump(streak_meta, f, indent=2)
print(f"Compute done. step={CLEAN_STATE['steps_done']}")
print(f"Saved streakline image: {streak_png}")
print(f"Saved vorticity image: {vort_png}")
print(
f"Vorticity stats: min={vort_end.min():+.6e}, max={vort_end.max():+.6e} "
f"(plot range [{vort_info['vmin']:+.3e}, {vort_info['vmax']:+.3e}])"
)
print(f"Saved streakline meta: {meta_path}")
print(f"Saved checkpoint: {CFD_PARAMS['checkpoint_path']}")chunk_step=10000/100000 global=10000 a=(+0.00000,+0.00000,+0.00000) chunk_step=20000/100000 global=20000 a=(+0.00000,+0.00000,+0.00000) chunk_step=30000/100000 global=30000 a=(+0.00000,+0.00000,+0.00000) chunk_step=40000/100000 global=40000 a=(+0.00000,+0.00000,+0.00000) chunk_step=50000/100000 global=50000 a=(+0.00000,+0.00000,+0.00000) chunk_step=60000/100000 global=60000 a=(+0.00000,+0.00000,+0.00000) chunk_step=70000/100000 global=70000 a=(+0.00000,+0.00000,+0.00000) chunk_step=80000/100000 global=80000 a=(+0.00000,+0.00000,+0.00000) chunk_step=90000/100000 global=90000 a=(+0.00000,+0.00000,+0.00000) chunk_step=100000/100000 global=100000 a=(+0.00000,+0.00000,+0.00000) Compute done. step=100000 Vorticity stats: min=-2.888157e-02, max=+2.866180e-02 Saved cache: /home/frank14f/CelerisLab/tests/triangle_latest_fields_clean.npz Saved checkpoint: /home/frank14f/CelerisLab/tests/triangle_state_clean.h5
In [7]:
# ==================== CLEAN 6. Plot only (repeat) ====================
# Streakline
img_path = Path(CFD_PARAMS["output_streak_png"])
if "streak_image" in CLEAN_STATE:
img = CLEAN_STATE["streak_image"]
else:
if not img_path.exists():
raise RuntimeError("No streakline image. Run CLEAN 5 first.")
img = plt.imread(img_path)
step = int(CLEAN_STATE.get("steps_done", -1))
sim = CLEAN_STATE.get("sim")
if sim is None:
raise RuntimeError("Simulation state not found. Run CLEAN 4 first.")
print(f"Plotting streakline at step={step}, image={img_path}")
fig = plt.figure(figsize=(sim.lbm_cfg.nx / 100.0, sim.lbm_cfg.ny / 100.0), dpi=100)
ax = fig.add_axes([0.0, 0.0, 1.0, 1.0])
ax.axis("off")
ax.imshow(img, origin="lower", interpolation="nearest", aspect="auto")
plt.show()
vort_path = Path(CFD_PARAMS["output_vort_png"])
if "vort_image" in CLEAN_STATE:
vort_img = CLEAN_STATE["vort_image"]
elif vort_path.exists():
vort_img = plt.imread(vort_path)
else:
raise RuntimeError("No vorticity image. Run CLEAN 5 first.")
print(f"Vorticity image: {vort_path}")
fig2 = plt.figure(figsize=(sim.lbm_cfg.nx / 100.0, sim.lbm_cfg.ny / 100.0), dpi=100)
ax2 = fig2.add_axes([0.0, 0.0, 1.0, 1.0])
ax2.axis("off")
ax2.imshow(vort_img, origin="lower", interpolation="nearest", aspect="auto")
plt.show()Plotting step=100000, vort min=-2.888157e-02, max=+2.866180e-02 Saved figure: /home/frank14f/CelerisLab/tests/vorticity_triangle_clean.png
In [ ]:
# ==================== CLEAN 7. Optional close ====================
# if CLEAN_STATE.get("sim") is not None:
# CLEAN_STATE["sim"].close()
# CLEAN_STATE["sim"] = None
# print("Simulation closed.")