""" TorsoCAE Validation -- cfd/laminar_flow case 2: Karman vortex street (2D cylinder in crossflow, Re=100). Classic vortex-shedding benchmark: above the critical Reynolds number (~Re=47 for an unconfined circular cylinder), the wake behind a bluff body sheds alternating vortices at a frequency f. The dimensionless Strouhal number St = f * D / U is one of the most widely reproduced numbers in CFD validation; at Re=100 the literature value sits in a narrow, well-established band, St ~ 0.16-0.17 (e.g. Williamson 1996 review; Roshko 1954 experiments extrapolate to St~0.16-0.17 at Re=100 on the lower-Re branch of the well-known St-Re curve). New backend capability needed and added for this case ------------------------------------------------------- Every prior CFD case in this suite reads a STEADY final state -- Strouhal number needs a per-step scalar TIME SERIES, which nothing in the transient laminar_flow driver produced (only the final step's field was ever persisted). Added an opt-in `force_tags` solver option to `_solve_cfd_laminar_flow_fields` (dolfinx_backend.py): when set, each converged time step calls the SAME `_compute_cfd_surface_forces` integral already used by the stokes driver (cfd/stokes/case_02) -- no duplicate force-computation code, per DRY -- on the requested surface tag(s), and accumulates (t, Fx, Fy) into `force_time`/`force_tags`/`force_values` arrays saved into the result npz. Threads through with zero new plumbing: solver options are already a generic passthrough dict (`session.set_solver_options(**opts)` only blocklists physics/contact/ coupling keys), so `force_tags=[cylinder_tag]` just works. Face tags verified fresh via `gmsh.model.occ.getCenterOfMass` on the actual exported STEP file for this exact geometry -- NOT assumed from either prior case's convention, both of which turned out to disagree with each other (case_01's 1x1x0.1 box: sequential 1=x-min,2=x-max,3=y-min,4=y-max,5=z-min, 6=z-max; cfd/stokes case_02's ~L^3 box-minus-sphere: interleaved 1=x-min,6=x-max,2=y-min,4=y-max,3=z-min,5=z-max). This box (20Dx10Dx0.2D, box-minus-cylinder, cylinder axis spanning the full z-thickness so its end caps coincide with the box's z faces rather than floating inside like the sphere) gave YET a third ordering: 1=x-min,2=y-min,3=z-min,4=y-max,5=z-max, 6=x-max,7=cylinder -- confirming, again, that this mesher's box face numbering is geometry-dependent and must be checked per case, never assumed. Solver: gmres+Schur-complement fieldsplit with lu sub-solves (never algo="direct" for this suite's mixed/saddle-point CFD systems -- standing instruction after algo="direct" stalled badly on larger systems earlier in this suite). SUPG stays off (supg_scale=0.0): CellDiameter hits an unsupported FFCx JIT path in this environment (see case_01), and Re=100 on a surface-refined P2/P1 mesh does not need convective stabilization to stay bounded (confirmed in the feasibility run below). A tiny geometric asymmetry (cylinder offset 0.02D off the channel centerline) is used to seed the shedding instability -- a standard, legitimate technique (not a hack): an impulsively-started, perfectly symmetric flow around a perfectly symmetric mesh can take an arbitrarily long, mesh/round-off-dependent time to spontaneously break symmetry, and the asymmetry only nudges *when* shedding starts, not the resulting frequency (a property of the flow, not the perturbation). Run: HWLOC_COMPONENTS=-gl python3 validation/cfd/laminar_flow/case_02_karman_vortex_street/run_validation.py """ import json import sys import time from pathlib import Path import numpy as np CASE_DIR = Path(__file__).resolve().parent sys.path.insert(0, str(CASE_DIR.parents[3] / "server")) from pytorsocae import TorsoCAESession D = 1.0 # m, cylinder diameter R = D / 2.0 DZ = 0.5 * D # quasi-2D thickness LX = 10.0 * D # channel length LY = 10.0 * D # channel height (blockage ratio D/LY = 0.1 -- doubled from the # first pass's 0.2, isolating blockage as the single changed # variable: that run reached a saturated (amplitude_stationary= # true) St=0.1863, still above the [0.16,0.17] unconfined band -- # see results_blockage_0.2.json/lift_history_blockage_0.2.png) XC = 3.5 * D # cylinder center, distance from inlet YC = LY / 2.0 + 0.02 * D # tiny offset off centerline to seed shedding RHO = 1.0 U = 1.0 RE = 100.0 MU = U * D * RHO / RE # = 0.01 ST_LOW, ST_HIGH = 0.16, 0.17 # literature band at Re=100 (Roshko / Williamson) SIZE_MAX = D / 2.5 # far-field background mesh size -- coarsened from D/5: at D/5 # (~45k nodes) this case's LU sub-solves ran the WSL host out # of memory partway through the transient run (no Python # traceback -- OOM-killed, confirmed via `free -h` showing only # ~10GB available). 3D unstructured LU fill-in scales much worse # than linearly with DOF count, so mesh size is the main lever. CYL_SIZE = D / 8.0 # cylinder surface local refinement (~pi*D/CYL_SIZE = 25 elements # around the circumference -- adequate for a bluff-body wake # instability, which is dominated by geometry-scale flow # structures, not boundary-layer-scale resolution) CYL_TRANSITION = 2.0 * SIZE_MAX # scales with background size (sphere-case lesson); doubled # here vs. that case -- a thin quasi-2D slab (DZ) needs more # gradation room between the refined cylinder size and DZ # itself, or TetGen hits the same "ScaledJac" quality crash # even with transition=size_max (hit once building this case). # n_cores=1 (serial), NOT the n_cores=4 used elsewhere in this suite: this # run checkpoints coordinate-tagged velocity/pressure state to disk between # batches (see fem_backends/checkpoint_state.py) so a host crash/restart # loses at most one batch instead of the whole run (see # run_transient_resumable below). The checkpoint format is partition- # independent (KDTree-matched by dof coordinate, not raw array order), so it # would in principle also survive an n_cores change between batches -- kept # at n_cores=1 here anyway because parallel laminar_flow currently has a # separate, unrelated known-bad-accuracy issue at n_cores>1 (tracked # separately, not a checkpoint mechanism limitation). The mesh is small # (7.5k nodes) so serial is still fast enough (~7s/step measured below). SOLVER_OPTIONS = dict( algo="gmres", precond="schur", tol=1e-6, max_iter=4000, n_cores=1, device="cpu", ksp_gmres_restart=250, fieldsplit_velocity_pc_type="bjacobi", fieldsplit_velocity_sub_pc_type="lu", fieldsplit_pressure_pc_type="bjacobi", fieldsplit_pressure_sub_pc_type="lu", ) BATCH_STEPS = 25 # ~1-2 min/batch -- a crash mid-batch loses at most this much progress CHECKPOINT_PATH = CASE_DIR / "checkpoint.npz" def _build_domain(session: TorsoCAESession) -> None: session.csg_make_shape("box", {"dx": LX, "dy": LY, "dz": DZ, "cx": 0.0, "cy": 0.0, "cz": 0.0}, name="Box", body_id="body_0") session.csg_make_shape("cylinder", {"r": R, "h": DZ, "cx": XC, "cy": YC, "cz": 0.0}, name="Cylinder", body_id="body_1") session.csg_boolean("cut", "body_0", "body_1", name="Domain", body_id="body_2") session.csg_select("body_2") def _mesh_domain(session: TorsoCAESession): # Pass 1: coarse probe, just to discover the cylinder's surface tag. probe = session.mesh(mesh_id="mesh_probe", algo_id="delaunay", size_factor=1.0, size_min=0, size_max=SIZE_MAX, order=2, dim=3) cyl_tag = min(probe["surfaces"], key=lambda s: s["elements"])["tag"] # Pass 2: re-mesh with the cylinder surface locally refined. mesh_res = session.mesh(mesh_id="mesh_0", algo_id="delaunay", size_factor=1.0, size_min=0, size_max=SIZE_MAX, order=2, dim=3, surface_sizes={cyl_tag: {"size": CYL_SIZE, "transition": CYL_TRANSITION}}) return mesh_res, cyl_tag def _set_physics(session: TorsoCAESession, cyl_tag: int, solver_options: dict) -> None: session.solid("Domain").material(mu=MU, rho=RHO, supg_scale=0.0) # Verified via getCenterOfMass on the actual STEP file for this exact # geometry (see module docstring) -- not assumed. session.surface(1, scope_id="body_2").bc("inlet_velocity", values=[U, 0, 0]) # x-min session.surface(6, scope_id="body_2").bc("outlet_pressure", values=[0.0]) # x-max session.surface(2, scope_id="body_2").bc("symmetry", values=[0.0], components=["y"]) # y-min session.surface(4, scope_id="body_2").bc("symmetry", values=[0.0], components=["y"]) # y-max session.surface(3, scope_id="body_2").bc("symmetry", values=[0.0], components=["z"]) # z-min session.surface(5, scope_id="body_2").bc("symmetry", values=[0.0], components=["z"]) # z-max # cyl_tag left unassigned -> CFD driver default-applies no-slip 'wall'. session.set_physics("cfd", submodel="laminar_flow", backend="dolfinx") session.set_solver_options(force_tags=[cyl_tag], **solver_options) def run_transient(num_steps: int, dt: float, tag: str) -> dict: """Single-shot (non-resumable) run -- used only for the short feasibility /timing probes, never for the long production run (see run_transient_resumable).""" session = TorsoCAESession() _build_domain(session) mesh_res, cyl_tag = _mesh_domain(session) _set_physics(session, cyl_tag, dict(SOLVER_OPTIONS, num_steps=num_steps, dt=dt)) t0 = time.time() result = session.compute(mesh_ids=["mesh_0"]) elapsed = time.time() - t0 npz = np.load(result["npz_path"]) out = { "tag": tag, "n_nodes": mesh_res["nodes"], "cyl_tag": cyl_tag, "cyl_elements": next(s["elements"] for s in mesh_res["surfaces"] if s["tag"] == cyl_tag), "num_steps": num_steps, "dt": dt, "elapsed_s": elapsed, "elapsed_s_per_step": elapsed / num_steps, "force_time": npz["force_time"].tolist(), "force_values": npz["force_values"].tolist(), # shape (n_steps, 1, gdim) } (CASE_DIR / f"probe_{tag}.json").write_text(json.dumps(out, indent=2)) return out def run_transient_resumable(total_steps: int, dt: float) -> dict: """Runs the transient solve in small batches, checkpointing the mixed velocity/pressure state + accumulated force history to disk after every batch (see CHECKPOINT_PATH). If this process is killed (this case has been interrupted by host restarts three times while building it), re-running this function picks up from the last completed batch instead of starting over -- at most BATCH_STEPS of work is lost, not the whole run.""" session = TorsoCAESession() _build_domain(session) mesh_res, cyl_tag = _mesh_domain(session) _set_physics(session, cyl_tag, dict(SOLVER_OPTIONS)) if CHECKPOINT_PATH.exists(): ck = np.load(CHECKPOINT_PATH) state = { "u_dof_coords": ck["u_dof_coords"], "velocity": ck["velocity"], "p_dof_coords": ck["p_dof_coords"], "pressure": ck["pressure"], } done_steps = int(ck["done_steps"]) force_time = ck["force_time"].tolist() force_values = ck["force_values"].tolist() elapsed_total = float(ck["elapsed_total"]) print(f"resuming from checkpoint: {done_steps}/{total_steps} steps already done, " f"{elapsed_total:.1f}s elapsed so far", flush=True) else: state = None done_steps = 0 force_time, force_values = [], [] elapsed_total = 0.0 while done_steps < total_steps: batch = min(BATCH_STEPS, total_steps - done_steps) session.set_solver_options(num_steps=batch, dt=dt, time=done_steps * dt, force_tags=[cyl_tag]) t0 = time.time() compute_kwargs = {"initial_state": state} if state is not None else {} result = session.compute(mesh_ids=["mesh_0"], **compute_kwargs) elapsed_total += time.time() - t0 npz = np.load(result["npz_path"]) state = { "u_dof_coords": npz["u_dof_coords"], "velocity": npz["velocity"], "p_dof_coords": npz["p_dof_coords"], "pressure": npz["pressure"], } force_time.extend(npz["force_time"].tolist()) force_values.extend(npz["force_values"].tolist()) done_steps += batch np.savez(CHECKPOINT_PATH, done_steps=done_steps, force_time=np.array(force_time), force_values=np.array(force_values), elapsed_total=elapsed_total, **state) print(f"batch done: {done_steps}/{total_steps} steps, " f"elapsed_total={elapsed_total:.1f}s ({elapsed_total / done_steps:.2f}s/step)", flush=True) return { "n_nodes": mesh_res["nodes"], "cyl_tag": cyl_tag, "cyl_elements": next(s["elements"] for s in mesh_res["surfaces"] if s["tag"] == cyl_tag), "num_steps": total_steps, "dt": dt, "elapsed_s": elapsed_total, "force_time": force_time, "force_values": force_values, } NUM_STEPS = 1200 # bumped from 450: first run's tail amplitude was still growing ~2x within # the tail window (not saturated), so St from that window isn't trustworthy. # Continuing (via checkpoint) to see whether St settles once shedding saturates. DT = 0.2 # Discard the impulsive-start transient before measuring frequency -- the # feasibility run showed Fx dropping from a startup spike (-4.7) to a settled # ~-0.41 within 5 steps, but shedding itself (an unstable-mode GROWTH # process) takes much longer to reach a periodic limit cycle. Keep only the # back half of the run for frequency extraction. TAIL_FRACTION = 0.5 def _period_from_zero_crossings(t: np.ndarray, y: np.ndarray) -> float | None: """Average period from positive-going zero crossings of the detrended (mean-removed) signal -- simple, robust to a non-integer number of periods in the window (unlike a raw FFT bin, which quantizes frequency to n_periods/window_length).""" y0 = y - y.mean() signs = np.sign(y0) crossings = np.where((signs[:-1] < 0) & (signs[1:] >= 0))[0] if len(crossings) < 2: return None # linear-interpolate the exact zero-crossing time within each bracketing step t_cross = [] for i in crossings: t0, t1 = t[i], t[i + 1] y0_, y1_ = y0[i], y0[i + 1] frac = -y0_ / (y1_ - y0_) if (y1_ - y0_) != 0 else 0.0 t_cross.append(t0 + frac * (t1 - t0)) periods = np.diff(t_cross) return float(np.mean(periods)) def main(): out = run_transient_resumable(total_steps=NUM_STEPS, dt=DT) t = np.array(out["force_time"]) fx = np.array(out["force_values"])[:, 0, 0] fy = np.array(out["force_values"])[:, 0, 1] n_tail = int(len(t) * TAIL_FRACTION) t_tail, fy_tail = t[-n_tail:], fy[-n_tail:] period = _period_from_zero_crossings(t_tail, fy_tail) amp_tail = float(np.std(fy_tail)) amp_first_half_of_tail = float(np.std(fy_tail[: len(fy_tail) // 2])) amp_second_half_of_tail = float(np.std(fy_tail[len(fy_tail) // 2:])) # Amplitude should have stopped growing appreciably by the tail window -- # otherwise the run is still in the growth phase, not a periodic limit # cycle, and any measured "period" is not trustworthy. amplitude_stationary = ( amp_second_half_of_tail > 0 and 0.5 < (amp_first_half_of_tail / amp_second_half_of_tail) < 2.0 ) st = (D / (period * U)) if period else None st_in_band = st is not None and (ST_LOW <= st <= ST_HIGH) try: import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt fig, ax = plt.subplots(figsize=(9, 4)) ax.plot(t, fy, lw=1.0, color="#2563eb") ax.axvspan(t[-n_tail], t[-1], color="#2563eb", alpha=0.08, label="tail window (frequency extraction)") ax.set_xlabel("time"); ax.set_ylabel("F_y on cylinder (N)") ax.set_title(f"Karman vortex street: lift oscillation, Re={RE:.0f}" + (f" -- St={st:.4f}" if st else "")) ax.legend(loc="upper right") fig.tight_layout() fig.savefig(CASE_DIR / "lift_history.png", dpi=140) except Exception as exc: print(f"plot failed (non-fatal): {exc}") summary = { "case": "cfd/laminar_flow / Karman vortex street (2D cylinder crossflow)", "citation": "St = f*D/U at Re=100 ~ 0.16-0.17 (Roshko 1954; Williamson 1996 review)", "geometry": {"D_m": D, "LX_m": LX, "LY_m": LY, "DZ_m": DZ, "blockage_ratio": D / LY}, "material": {"mu": MU, "rho": RHO}, "U": U, "Re": RE, "mesh": {"n_nodes": out["n_nodes"], "cyl_tag": out["cyl_tag"], "cyl_elements": out["cyl_elements"], "size_max": SIZE_MAX, "cyl_size": CYL_SIZE, "cyl_transition": CYL_TRANSITION}, "num_steps": NUM_STEPS, "dt": DT, "elapsed_s": out["elapsed_s"], "tail_fraction": TAIL_FRACTION, "period": period, "amplitude_first_half_of_tail": amp_first_half_of_tail, "amplitude_second_half_of_tail": amp_second_half_of_tail, "amplitude_stationary": bool(amplitude_stationary), "St": st, "St_reference_band": [ST_LOW, ST_HIGH], "St_in_band": bool(st_in_band), "Fx_settled_mean": float(np.mean(fx[-n_tail:])), "pass": bool(st_in_band and amplitude_stationary), } (CASE_DIR / "results.json").write_text(json.dumps(summary, indent=2)) print(json.dumps(summary, indent=2), flush=True) # NOTE: the checkpoint is deliberately kept even on a completed run -- a # "pass=false, needs more time" verdict is a normal outcome (see this # case's own history), and deleting already-computed progress here once # cost a full restart-from-zero recovery. Only NUM_STEPS needs to grow; # run_transient_resumable() will pick up from here and run just the delta. if __name__ == "__main__": main()