diff --git a/pyproject.toml b/pyproject.toml index 55ecc56..672663f 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -7,7 +7,7 @@ name = "CelerisLab" version = "0.5.1" description = "GPU-accelerated Lattice Boltzmann Method (LBM) CFD solver using CUDA" readme = "README.md" -requires-python = ">=3.8" +requires-python = ">=3.10" license = {text = "MIT"} authors = [ {name = "Frank14f"} @@ -19,8 +19,6 @@ classifiers = [ "Topic :: Scientific/Engineering :: Physics", "License :: OSI Approved :: MIT License", "Programming Language :: Python :: 3", - "Programming Language :: Python :: 3.8", - "Programming Language :: Python :: 3.9", "Programming Language :: Python :: 3.10", "Programming Language :: Python :: 3.11", ] diff --git a/setup.py b/setup.py index 996a3d4..df0b419 100644 --- a/setup.py +++ b/setup.py @@ -26,15 +26,13 @@ setup( 'numpy>=1.19.0', 'scipy>=1.5.0', ], - python_requires='>=3.8', + python_requires='>=3.10', classifiers=[ 'Development Status :: 3 - Alpha', 'Intended Audience :: Science/Research', 'Topic :: Scientific/Engineering :: Physics', 'License :: OSI Approved :: MIT License', 'Programming Language :: Python :: 3', - 'Programming Language :: Python :: 3.8', - 'Programming Language :: Python :: 3.9', 'Programming Language :: Python :: 3.10', 'Programming Language :: Python :: 3.11', ], diff --git a/src/CelerisLab/__init__.py b/src/CelerisLab/__init__.py index 519caa0..71a207b 100644 --- a/src/CelerisLab/__init__.py +++ b/src/CelerisLab/__init__.py @@ -1,28 +1,16 @@ -# CelerisLab/__init__.py -""" -CelerisLab: GPU-Accelerated Computational Physics Simulation Library +"""GPU LBM library. CPU utilities do not initialize or import CUDA eagerly.""" +from importlib import import_module -Usage:: +__version__ = "0.5.1" +__all__ = ["Simulation", "lbm", "body", "common", "cuda", "config", "__version__"] - from CelerisLab import Simulation - sim = Simulation("configs/config_lbm.json") - sim.add_body("circle", center=(100, 50), radius=10) - sim.initialize() - sim.run(1000) - data = sim.get_macroscopic() - force = sim.read_force(0) -""" - -__version__ = "0.5.0" - -from . import common, cuda, lbm, body, config - -try: - from .simulation import Simulation - __all__ = ["Simulation", "lbm", "body", "common", "cuda", "config", - "__version__"] -except ImportError as e: - import warnings - warnings.warn(f"Simulation not available: {e}", ImportWarning) - __all__ = ["lbm", "body", "common", "cuda", "config", "__version__"] +def __getattr__(name): + if name == "Simulation": + value = import_module(".simulation", __name__).Simulation + elif name in ("lbm", "body", "common", "cuda", "config"): + value = import_module("." + name, __name__) + else: + raise AttributeError(name) + globals()[name] = value + return value diff --git a/src/CelerisLab/body/__init__.py b/src/CelerisLab/body/__init__.py index 79296de..adc81e8 100644 --- a/src/CelerisLab/body/__init__.py +++ b/src/CelerisLab/body/__init__.py @@ -2,8 +2,8 @@ """ Body / object management for immersed and rigid objects. """ -from .objects import SimObject, ObjectState, ObjectControl -from .manager import ObjectManager +from .objects import SimObject, ObjectState, ObjectControl, SurfaceActuation +from importlib import import_module from .registry import BodyRegistry from .sync_plan import BodySyncPlan from .geometry.base import CutLink, SensorCell, Geometry @@ -11,7 +11,7 @@ from .geometry.circle import CircleGeometry from .coupling.soa_packer import pack_cut_links_to_soa, pack_sensor_to_soa __all__ = [ - "SimObject", "ObjectState", "ObjectControl", + "SimObject", "ObjectState", "ObjectControl", "SurfaceActuation", "ObjectManager", "BodyRegistry", "BodySyncPlan", @@ -19,3 +19,11 @@ __all__ = [ "CircleGeometry", "pack_cut_links_to_soa", "pack_sensor_to_soa", ] + + +def __getattr__(name): + if name != "ObjectManager": + raise AttributeError(name) + value = import_module(".manager", __name__).ObjectManager + globals()[name] = value + return value diff --git a/src/CelerisLab/body/ddf_patch.py b/src/CelerisLab/body/ddf_patch.py index c145345..0b25cd7 100644 --- a/src/CelerisLab/body/ddf_patch.py +++ b/src/CelerisLab/body/ddf_patch.py @@ -56,7 +56,7 @@ def patch_ddf_for_body_sync( Args: field: LBMField instance (``field.ddf`` is modified in-place). - old_ddf: Physical DDF snapshot before topology change (n*NQ, float32). + old_ddf: Q-major physical DDF snapshot (NQ*n, float32). old_flags: Flag array before topology change. new_flags: Flag array after topology change. added_solid_mask: Bool array -- True where fluid became solid. @@ -78,31 +78,34 @@ def upload_patched_ddf( ) -> None: """Upload only the DDF cells that changed during patching. - Uses per-cell ``cuda.memcpy_htod`` for both ``ddf_gpu`` and - ``temp_gpu``. Leaves all other cells (especially boundary cells + Copies Q-major contiguous cell runs to both ``ddf_gpu`` and + ``temp_gpu`` using the field storage encoding. Leaves all other cells (especially boundary cells with NaN) untouched on the GPU, preserving their streaming- compatible values from before the sync. """ changed_indices: set[int] = set() changed_indices.update(np.where(added_solid_mask)[0].tolist()) changed_indices.update(np.where(released_fluid_mask)[0].tolist()) - if curved_fluid_indices is not None: + if curved_fluid_indices is not None and np.any(added_solid_mask): changed_indices.update(curved_fluid_indices.tolist()) if not changed_indices: field.invalidate_host_ddf_cache() field._host_ddf_step = None return - nq = field.nq - cell_bytes = nq * 4 - ddf_ptr = int(field.ddf_gpu) - temp_ptr = int(field.temp_gpu) - - for idx in sorted(changed_indices): - offset = idx * cell_bytes - cell_data = field.ddf[idx * nq:(idx + 1) * nq].tobytes() - cuda.memcpy_htod(ddf_ptr + offset, cell_data) - cuda.memcpy_htod(temp_ptr + offset, cell_data) + # Bulk fields are Q-major; the compact EsoPull scratch is cell-major. + # Encode only selected values with the same shift/precision path as a + # whole-field upload, then copy contiguous cell runs in each direction. + indices = np.array(sorted(changed_indices), dtype=np.intp) + encoded = field._encode_ddf(field.ddf.reshape(field.nq, field.n)[:, indices]) + breaks = np.flatnonzero(np.diff(indices) != 1) + 1 + for positions in np.split(np.arange(indices.size), breaks): + start = int(indices[positions[0]]) + for direction in range(field.nq): + offset = (direction * field.n + start) * field.store_bytes + values = np.ascontiguousarray(encoded[direction, positions]) + cuda.memcpy_htod(int(field.ddf_gpu) + offset, values) + cuda.memcpy_htod(int(field.temp_gpu) + offset, values) field.invalidate_host_ddf_cache() field._host_ddf_step = None @@ -119,7 +122,7 @@ def _patch_added_solid(field, added_solid_mask: np.ndarray) -> None: indices = np.where(added_solid_mask)[0] nq = field.nq for idx in indices: - field.ddf[idx * nq:(idx + 1) * nq] = feq + field.ddf.reshape(nq, field.n)[:, idx] = feq def _patch_boundary_fluid( @@ -135,13 +138,13 @@ def _patch_boundary_fluid( nq = field.nq feq_fallback = compute_feq_d2q9(field.cfg.rho, 0.0, 0.0) for idx in curved_fluid_indices: - f = field.ddf[idx * nq:(idx + 1) * nq] + f = field.ddf.reshape(nq, field.n)[:, idx] if np.any(np.isnan(f)): - field.ddf[idx * nq:(idx + 1) * nq] = feq_fallback + field.ddf.reshape(nq, field.n)[:, idx] = feq_fallback continue rho, ux, uy = compute_macro_from_ddf(f) feq = compute_feq_d2q9(float(rho), float(ux), float(uy)) - field.ddf[idx * nq:(idx + 1) * nq] = feq + field.ddf.reshape(nq, field.n)[:, idx] = feq def _patch_released_fluid( @@ -173,7 +176,7 @@ def _patch_released_fluid( old_macro: dict[int, tuple[float, float, float]] = {} fluid_indices = np.where(old_is_fluid)[0] for idx in fluid_indices: - f = old_ddf[idx * nq:(idx + 1) * nq] + f = old_ddf.reshape(nq, field.n)[:, idx] if np.any(np.isnan(f)): continue rho, ux, uy = compute_macro_from_ddf(f) @@ -199,7 +202,7 @@ def _patch_released_fluid( ux_avg = np.mean([old_macro[n][1] for n in fluid_nbrs]) uy_avg = np.mean([old_macro[n][2] for n in fluid_nbrs]) feq = compute_feq_d2q9(float(rho_avg), float(ux_avg), float(uy_avg)) - field.ddf[idx * nq:(idx + 1) * nq] = feq + field.ddf.reshape(nq, field.n)[:, idx] = feq filled[idx] = True prev_layer_filled[idx] = ( float(rho_avg), float(ux_avg), float(uy_avg)) @@ -221,7 +224,7 @@ def _patch_released_fluid( uy_avg = np.mean([prev_layer_filled[n][2] for n in filled_nbrs]) feq = compute_feq_d2q9( float(rho_avg), float(ux_avg), float(uy_avg)) - field.ddf[idx * nq:(idx + 1) * nq] = feq + field.ddf.reshape(nq, field.n)[:, idx] = feq filled[idx] = True current_layer[idx] = ( float(rho_avg), float(ux_avg), float(uy_avg)) @@ -236,4 +239,4 @@ def _patch_released_fluid( u_inlet * 4.0 * y * (ny - 1 - y) / ((ny - 1) ** 2) if ny > 2 else u_inlet) feq = compute_feq_d2q9(field.cfg.rho, float(ux_fallback), 0.0) - field.ddf[idx * nq:(idx + 1) * nq] = feq + field.ddf.reshape(nq, field.n)[:, idx] = feq diff --git a/src/CelerisLab/body/manager.py b/src/CelerisLab/body/manager.py index dcb7d2d..46a5592 100644 --- a/src/CelerisLab/body/manager.py +++ b/src/CelerisLab/body/manager.py @@ -31,6 +31,7 @@ from .preprocess.flag_overlay import merge_flag_masks from ..cuda.compiler_v2 import obs_layout from ..cuda.obs_layout import ObsLayout from ..lbm.descriptors import ( + FLUID, SOLID, D2Q9_EX, D2Q9_EY, D3Q19_EX, @@ -76,6 +77,7 @@ class ObjectManager: self._obs_alloc_nbytes: int = 0 self.action = np.zeros(0, dtype=np.float32) + self._action_pinned: Optional[np.ndarray] = None self.obs_pinned: Optional[np.ndarray] = None self.obs_force_view: Optional[np.ndarray] = None self.obs_torque_view: Optional[np.ndarray] = None @@ -106,13 +108,23 @@ class ObjectManager: # -- Registry delegation ------------------------------------------------- def add(self, obj: SimObject) -> int: """Register an object and return its id.""" + obj.validate_surface_actuation(self.cfg) body_id = self._registry.add(obj) self._resize_buffers() return body_id def remove(self, obj_id: int): + """Remove immediately and renumber remaining IDs, as runtime sync does. + + Previously returned numeric IDs are invalidated by this operation. + """ self._registry.remove(obj_id) + remaining = self.objects + self._registry.clear() + for obj in remaining: + self._registry.add(obj) self._resize_buffers() + self._refresh_action_from_objects() def get(self, obj_id: int) -> SimObject: return self._registry.get(obj_id) @@ -149,6 +161,7 @@ class ObjectManager: The object is not added to the formal registry until ``commit_pending()`` is called. This implicitly activates the edit window. """ + obj.validate_surface_actuation(self.cfg) self._edit_active = True self._pending_add.append(obj) @@ -190,7 +203,8 @@ class ObjectManager: new = SimObject(obj_id=oid, geometry=obj.geometry, center=obj.center, radius=obj.radius, is_sensor=obj.is_sensor, - is_force_region=obj.is_force_region) + is_force_region=obj.is_force_region, + surface_actuation=obj.surface_actuation) new.state = obj.state new.control = obj.control return new @@ -295,6 +309,7 @@ class ObjectManager: cl_host = plan.curved_host field.curved.assign_host(*cl_host) + field.curved.assign_surface(plan.next_objects, self.cfg, plan.next_flags) s_cells, s_ids = plan.sensor_host field.sensors.assign_host(s_cells, s_ids) fr_cells, fr_ids = plan.force_region_host @@ -363,6 +378,13 @@ class ObjectManager: cl_rx, cl_ry, cl_rz, cl_fallback = [], [], [], [] s_cells, s_ids = [], [] fr_cells, fr_ids = [], [] + if domain_flags is None: + domain_flags = self.build_flags_for( + objects, np.full(self.nx * self.ny * self.nz, FLUID, dtype=np.uint16), + nx=self.nx, ny=self.ny, nz=self.nz, + ) + active_fluid = ((domain_flags & FLUID) != 0) & ((domain_flags & SOLID) == 0) + force_owners: set[int] = set() if sensor_cell_counts_out is not None: sensor_cell_counts_out.fill(0) @@ -374,6 +396,7 @@ class ObjectManager: ex, ey = D2Q9_EX, D2Q9_EY for obj in objects: + obj.validate_surface_actuation(self.cfg) if (not obj.is_sensor) and (not obj.is_force_region) and hasattr(obj, 'get_curved_list'): ( fluid_idx, dirs, q_vals, body_ids, @@ -394,6 +417,8 @@ class ObjectManager: cl_fallback.append(fallback_vals) if obj.is_sensor and hasattr(obj, 'get_sensor_list'): cells, ids = obj.get_sensor_list(self.nx, self.ny, self.nz) + keep = active_fluid[cells] + cells, ids = cells[keep], ids[keep] if len(cells) > 0: s_cells.append(cells) s_ids.append(ids) @@ -402,7 +427,12 @@ class ObjectManager: sensor_cell_counts_out[obj.obj_id] = int(len(cells)) if obj.is_force_region and hasattr(obj, 'get_sensor_list'): cells, ids = obj.get_sensor_list(self.nx, self.ny, self.nz) + keep = active_fluid[cells] + cells, ids = cells[keep], ids[keep] if len(cells) > 0: + if force_owners.intersection(map(int, cells)): + raise ValueError("Overlapping force_region fluid cells are unsupported") + force_owners.update(map(int, cells)) fr_cells.append(cells) fr_ids.append(ids) @@ -422,6 +452,13 @@ class ObjectManager: r_fr_cells = np.concatenate(fr_cells) if fr_cells else np.zeros(0, dtype=np.uint32) r_fr_ids = np.concatenate(fr_ids) if fr_ids else np.zeros(0, dtype=np.int32) + if any(obj.surface_actuation is not None for obj in objects): + # Reject unresolved apertures/ownership before sync or restart mutates state. + from ..lbm.curved_links import CurvedLinkSoA + links = CurvedLinkSoA() + links.assign_host(r_idx, r_dir, r_q, r_bid, r_rx, r_ry, r_rz, r_fallback) + links.assign_surface(objects, self.cfg, domain_flags) + return ( r_idx, r_dir, r_q, r_bid, r_rx, r_ry, r_rz, r_fallback, @@ -441,21 +478,24 @@ class ObjectManager: """ self._telemetry_field = field - if rebuild_flags: - field.flag = self.build_flags(field.build_channel_flags()) - field.upload_flags() + next_flags = self.build_flags(field.build_channel_flags()) if rebuild_flags else field.flag ( cl_fluid_idx, cl_dir, cl_q, cl_body_id, cl_rx, cl_ry, cl_rz, cl_fallback_class, sensor_cells, sensor_obj_id, fr_cells, fr_obj_id, - ) = self.build_compact_lists(domain_flags=field.flag) + ) = self.build_compact_lists(domain_flags=next_flags) + + if rebuild_flags: + field.flag = next_flags + field.upload_flags() field.curved.assign_host( cl_fluid_idx, cl_dir, cl_q, cl_body_id, cl_rx, cl_ry, cl_rz, cl_fallback_class, ) + field.curved.assign_surface(self.objects, self.cfg, field.flag) field.sensors.assign_host(sensor_cells, sensor_obj_id) field.force_regions.assign_host(fr_cells, fr_obj_id) @@ -570,9 +610,16 @@ class ObjectManager: def _upload_action_async(self, stream: cuda.Stream) -> None: """Async H2D of host action array to action_gpu.""" + for obj in self.objects: + obj.validate_surface_actuation(self.cfg) + if obj.surface_actuation is not None and self.action[(obj.obj_id+1)*3*self.cfg.dim-1] != 0: + raise ValueError("Surface actuation cannot share the rotation action") if self.action_gpu is None or self.action.size == 0: return - cuda.memcpy_htod_async(self.action_gpu, self.action, stream) + if self._action_pinned is None or self._action_pinned.size != self.action.size: + self._action_pinned = cuda.pagelocked_empty(self.action.size, np.float32) + self._action_pinned[:] = self.action + cuda.memcpy_htod_async(self.action_gpu, self._action_pinned, stream) @property def obs_n_slots(self) -> int: @@ -684,9 +731,12 @@ class ObjectManager: The GPU action buffer is updated at the next ``run()`` call. """ self._validate_body_id(body_id) + obj = self.get(body_id) + obj.validate_surface_actuation(self.cfg, dict(vars(obj.state), omega=omega)) dim = self.cfg.dim if self.cfg else 2 slot = 3 * dim base = body_id * slot + self.get(body_id).state.omega = float(omega) self.action[base + slot - 1] = np.float32(omega) def set_force_state(self, body_id: int, fx: float = 0.0, fy: float = 0.0) -> None: @@ -698,6 +748,14 @@ class ObjectManager: dim = self.cfg.dim if self.cfg else 2 slot = 3 * dim base = body_id * slot + obj = self.get(body_id) + if not obj.is_force_region: + raise ValueError("set_force_state requires a force_region object") + # vx/vy are the existing force-region action aliases. Keep their + # authoritative host state in sync so edits to other bodies cannot + # discard this control during a full action refresh or topology sync. + obj.state.vx = float(fx) + obj.state.vy = float(fy) self.action[base + 0] = np.float32(fx) self.action[base + 1] = np.float32(fy) diff --git a/src/CelerisLab/body/objects.py b/src/CelerisLab/body/objects.py index 02a45fc..de3afbc 100644 --- a/src/CelerisLab/body/objects.py +++ b/src/CelerisLab/body/objects.py @@ -15,6 +15,7 @@ Design: from __future__ import annotations from dataclasses import dataclass +import math from typing import Optional, Tuple import numpy as np @@ -53,6 +54,37 @@ class ObjectControl: mode: int = 0 # 0=force, 1=velocity, 2=displacement +@dataclass(frozen=True) +class SurfaceActuation: + """Fixed circular aperture: angles in radians, velocity in lattice units. + + ``angle`` is measured counterclockwise from +x; ``half_width`` is the arc + half-angle. Positive normal_velocity blows from solid into fluid. The C1 + profile (1-s*s)^2 uses the true hit angle, not a fluid-node aperture mask. + This immutable specification is set before initialization; no moving nozzle. + """ + angle: float + half_width: float + normal_velocity: float + + def __post_init__(self): + if not all(math.isfinite(v) for v in (self.angle, self.half_width, self.normal_velocity)): + raise ValueError("Surface actuation values must be finite") + if not 0 < self.half_width <= math.pi or abs(self.normal_velocity) >= 1: + raise ValueError("Require 0 < half_width <= pi and abs(normal_velocity) < 1") + + def project(self, rx, ry): + """Project stored true hit-center vectors; returns velocity and aperture mask.""" + radius = np.hypot(rx, ry) + if np.any(radius <= 0): + raise ValueError("Surface normal requires nonzero hit-center distance") + delta = np.arctan2(np.sin(np.arctan2(ry, rx)-self.angle), + np.cos(np.arctan2(ry, rx)-self.angle)) + aperture = np.abs(delta) < self.half_width + speed = self.normal_velocity*np.maximum(0., 1-(delta/self.half_width)**2)**2 + return np.stack((speed*rx/radius, speed*ry/radius), axis=-1).astype(np.float32), aperture + + # --------------------------------------------------------------------------- # SimObject container # --------------------------------------------------------------------------- @@ -77,6 +109,7 @@ class SimObject: radius: float = 0.0, is_sensor: bool = False, is_force_region: bool = False, + surface_actuation: Optional[SurfaceActuation] = None, ) -> None: self.obj_id = obj_id self.center = center if center is not None else (0.0, 0.0) @@ -84,6 +117,10 @@ class SimObject: self._geometry = geometry self._is_sensor = is_sensor self._is_force_region = is_force_region + self._surface_actuation = (SurfaceActuation(**surface_actuation) + if isinstance(surface_actuation, dict) else surface_actuation) + if self._surface_actuation is not None and not isinstance(self._surface_actuation, SurfaceActuation): + raise TypeError("surface_actuation must be a SurfaceActuation or dict") c = self.center z0 = float(c[2]) if len(c) > 2 else 0.0 @@ -92,6 +129,22 @@ class SimObject: ) self.control = ObjectControl() + @property + def surface_actuation(self): + return self._surface_actuation + + def validate_surface_actuation(self, cfg, state=None): + if self.surface_actuation is None: + return + from .geometry.circle import CircleGeometry + if (type(self.geometry) is not CircleGeometry or self.is_sensor or self.is_force_region + or cfg is None or not cfg.is_d2q9 or cfg.store_precision != "FP32" + or cfg.streaming != "double_buffer" or cfg.ddf_shifting): + raise ValueError("Surface actuation requires a circle and D2Q9/FP32/unshifted/double_buffer") + values = vars(self.state) if state is None else state + if any(values.get(k, 0.) != 0 for k in ("vx", "vy", "omega")): + raise ValueError("Surface actuation requires a stationary circle") + @property def geometry(self) -> Geometry: if self._geometry is None: diff --git a/src/CelerisLab/body/preprocess/flag_overlay.py b/src/CelerisLab/body/preprocess/flag_overlay.py index 80fc385..5951190 100644 --- a/src/CelerisLab/body/preprocess/flag_overlay.py +++ b/src/CelerisLab/body/preprocess/flag_overlay.py @@ -9,6 +9,8 @@ from typing import List import numpy as np +from ...lbm.descriptors import FLUID, SOLID, SENSOR_FLAG, FRC_REGION + def merge_flag_masks( base: np.ndarray, @@ -16,9 +18,9 @@ def merge_flag_masks( ) -> np.ndarray: """Merge object flag masks onto a clean domain base. - For each mask in *masks*, non-zero entries overwrite the corresponding - positions in *base*. Callers should pass a fresh channel-layout flag - array to prevent stale obstacle bits from surviving after geometry edits. + Apply solid occupancy first, then sensor/force-region annotations only + to fluid cells. Object order cannot turn an obstacle or channel wall into + fluid. Callers supply a fresh channel-layout base after geometry edits. Args: base: uint16 domain flag array (e.g. from build_channel_flags()). @@ -29,6 +31,10 @@ def merge_flag_masks( """ merged = np.array(base, copy=True) for mask in masks: - nonzero = mask != 0 - merged[nonzero] = mask[nonzero] + solid = (mask & SOLID) != 0 + merged[solid] = mask[solid] + fluid = ((merged & FLUID) != 0) & ((merged & SOLID) == 0) + for mask in masks: + annotations = mask & (SENSOR_FLAG | FRC_REGION) + merged[fluid] |= annotations[fluid] return merged diff --git a/src/CelerisLab/common/__init__.py b/src/CelerisLab/common/__init__.py index fdb21d1..18e1b47 100644 --- a/src/CelerisLab/common/__init__.py +++ b/src/CelerisLab/common/__init__.py @@ -1,8 +1,12 @@ -# CelerisLab/common/__init__.py -""" -Common utilities and preprocessing functions. -""" - -from . import preprocess, pathline, checkpoint, render, streakline +"""Shared host utilities; optional GPU/plotting dependencies load on demand.""" +from importlib import import_module __all__ = ["preprocess", "pathline", "checkpoint", "render", "streakline"] + + +def __getattr__(name): + if name not in __all__: + raise AttributeError(name) + value = import_module("." + name, __name__) + globals()[name] = value + return value diff --git a/src/CelerisLab/common/checkpoint.py b/src/CelerisLab/common/checkpoint.py index cd3b4f4..b4fcd73 100644 --- a/src/CelerisLab/common/checkpoint.py +++ b/src/CelerisLab/common/checkpoint.py @@ -12,15 +12,19 @@ Usage:: ``lattice_model``, ``nx``, ``ny``, ``nz``, and ``store_precision``. """ +import hashlib import json +from dataclasses import asdict import os import time import numpy as np import pycuda.driver as cuda +from ..lbm.descriptors import FLUID, SOLID, OBSTACLE, SENSOR_FLAG, FRC_REGION + # File format version — bump on incompatible layout changes. -CKPT_VERSION = 1 +CKPT_VERSION = 4 _SUPPORTED_STORE_PRECISIONS = ("FP32", "FP16S") @@ -33,6 +37,66 @@ def _require_supported_precision(store_precision: str, source: str) -> None: ) +def _domain_flags(field, bodies, flags): + """Keep actual domain boundaries, removing only registered object overlays. + + A solid object's hidden underlay is unknowable from merged flags: use FLUID + there canonically. This base is for fixed-geometry restart, not body removal. + """ + base = flags & np.uint16(0xffff ^ (SENSOR_FLAG | FRC_REGION)) + for obj in bodies.objects: + mask = obj.get_flag_mask(field.nx, field.ny, field.nz) + base[(mask & SOLID) != 0] = FLUID + if np.any(base & OBSTACLE): + raise ValueError("Checkpoint flags contain an unregistered obstacle") + return base + + +def _geometry_records(field, bodies, base_flags=None): + """Describe actual geometry/projection, not just the number of objects.""" + records = [] + flags = bodies.build_flags(field.build_channel_flags() if base_flags is None else base_flags) + for obj in bodies.objects: + digest = hashlib.sha256() + digest.update(obj.get_flag_mask(field.nx, field.ny, field.nz).tobytes()) + if obj.is_sensor or obj.is_force_region: + arrays = obj.get_sensor_list(field.nx, field.ny, field.nz) + else: + arrays = obj.get_curved_list(field.nx, field.ny, field.nq, + nz=field.nz, domain_flags=flags) + for array in arrays: + digest.update(np.ascontiguousarray(array).tobytes()) + geometry = obj.geometry + records.append({ + "obj_id": obj.obj_id, + "geometry": type(geometry).__module__ + "." + type(geometry).__qualname__, + "center": list(obj.center), "radius": float(obj.radius), + "is_sensor": obj.is_sensor, "is_force_region": obj.is_force_region, + "projection_sha256": digest.hexdigest(), + }) + return records, flags + + +def _validate_legacy_objects(obj_data, bodies): + """v1 can safely restart built-in circles/sensors with unchanged geometry. + + It did not save authoritative force-region controls. Refuse those rather + than silently restoring zero force. Unknown/custom geometry is not inferred. + """ + from ..body.geometry.circle import CircleGeometry + + for od, obj in zip(obj_data, bodies.objects): + if obj.is_force_region: + raise ValueError("Legacy v1 checkpoint cannot restore force_region controls") + if type(obj.geometry) is not CircleGeometry: + raise ValueError("Legacy v1 checkpoint only supports built-in circle/sensor geometry") + expected = (float(od["x"]), float(od["y"]), float(od.get("z", 0.0))) + actual = (*obj.geometry.center, float(obj.state.z)) + if (expected != actual or float(od["radius"]) != obj.geometry._radius + or od.get("type") != obj.obj_type): + raise ValueError("Checkpoint object geometry/type mismatch") + + def save_checkpoint(field, stepper, lbm_cfg, bodies, path=None): """Save full simulation state to an HDF5 file. @@ -61,10 +125,25 @@ def save_checkpoint(field, stepper, lbm_cfg, bodies, path=None): cuda.memcpy_dtoh(tmp_buf, field.temp_gpu) flag_buf = np.empty(field.n, dtype=np.uint16) cuda.memcpy_dtoh(flag_buf, field.flag_gpu) + domain_flags = _domain_flags(field, bodies, flag_buf) + geometry, expected_flags = _geometry_records(field, bodies, domain_flags) + if not np.array_equal(flag_buf, expected_flags): + raise ValueError("Checkpoint flags do not match current object/domain geometry") + bodies.build_compact_lists_for(bodies.objects, domain_flags=expected_flags) + applied_action = None + if bodies.action_gpu is not None: + applied_action = np.empty_like(bodies.action) + cuda.memcpy_dtoh(applied_action, bodies.action_gpu) + for obj in bodies.objects: + if (obj.surface_actuation is not None + and applied_action[(obj.obj_id + 1) * 3 * lbm_cfg.dim - 1] != 0): + raise ValueError("Surface actuation requires zero applied omega") with h5py.File(path, "w") as hf: # Version and metadata - hf.attrs["ckpt_version"] = CKPT_VERSION + # Older readers must not silently restart a driven surface as a wall. + hf.attrs["ckpt_version"] = (CKPT_VERSION if any( + obj.surface_actuation is not None for obj in bodies.objects) else 3) hf.attrs["celerislab_version"] = __import__( "CelerisLab").__version__ hf.attrs["timestamp"] = time.strftime("%Y-%m-%dT%H:%M:%S") @@ -74,25 +153,7 @@ def save_checkpoint(field, stepper, lbm_cfg, bodies, path=None): hf.attrs["n_objects"] = bodies.count # Config as JSON string (human-readable in HDF5 viewers) - cfg = lbm_cfg - config_dict = { - "lattice_model": cfg.lattice_model, - "nx": cfg.nx, "ny": cfg.ny, "nz": cfg.nz, - "data_type": cfg.data_type, - "viscosity": cfg.viscosity, "velocity": cfg.velocity, - "rho": cfg.rho, - "collision": cfg.collision, "streaming": cfg.streaming, - "store_precision": cfg.store_precision, - "ddf_shifting": cfg.ddf_shifting, - "les_enabled": cfg.les_enabled, "les_cs": cfg.les_cs, - "les_closed_form": cfg.les_closed_form, - "trt_magic_param": cfg.trt_magic_param, - "inlet_profile": cfg.inlet_profile, - "inlet_scheme": cfg.inlet_scheme, - "outlet_mode": cfg.outlet_mode, - "omega_min": cfg.omega_min, "omega_max": cfg.omega_max, - } - hf.attrs["config_json"] = json.dumps(config_dict) + hf.attrs["config_json"] = json.dumps(asdict(lbm_cfg)) # GPU arrays — gzip compression for portable, compact files hf.create_dataset("ddf", data=ddf_buf, compression="gzip", @@ -101,8 +162,24 @@ def save_checkpoint(field, stepper, lbm_cfg, bodies, path=None): compression_opts=4) hf.create_dataset("flag", data=flag_buf, compression="gzip", compression_opts=4) + hf.create_dataset("domain_flag", data=domain_flags, compression="gzip", + compression_opts=4) - # Object states + hf.attrs["geometry_json"] = json.dumps(geometry) + hf.attrs["surface_actuation_json"] = json.dumps([ + asdict(obj.surface_actuation) if obj.surface_actuation is not None else None + for obj in bodies.objects + ]) + hf.attrs["runtime_params_json"] = json.dumps(field._runtime_params) + hf.attrs["obs_accum_steps"] = bodies._obs_accum_steps + if applied_action is not None: + hf.create_dataset("applied_action", data=applied_action) + if bodies.obs_gpu is not None: + obs = np.empty(bodies.obs_total_floats, dtype=np.float32) + cuda.memcpy_dtoh(obs, bodies.obs_gpu) + hf.create_dataset("obs", data=obs) + + # Object states and host controls (next run uploads the restored action). obj_states = [] for obj in bodies.objects: s = obj.state @@ -112,6 +189,7 @@ def save_checkpoint(field, stepper, lbm_cfg, bodies, path=None): "x": s.x, "y": s.y, "z": s.z, "theta": s.theta, "vx": s.vx, "vy": s.vy, "omega": s.omega, "radius": obj.radius, + "control": asdict(obj.control), }) hf.attrs["objects_json"] = json.dumps(obj_states) @@ -124,6 +202,9 @@ def load_checkpoint(path, field, stepper, lbm_cfg, bodies): The current config's lattice_model, nx, ny, nz, and store_precision must match the checkpoint. Other parameters (collision, LES, …) may differ for parameter-study restarts. + v3 restores the saved domain boundaries; v1/v2 retain their default-channel + validation. v4 also requires matching static surface actuation. Object + geometry and roles must match in every version. Args: path: HDF5 checkpoint file path. @@ -138,7 +219,7 @@ def load_checkpoint(path, field, stepper, lbm_cfg, bodies): with h5py.File(path, "r") as hf: # Version check ver = hf.attrs.get("ckpt_version", 0) - if ver != CKPT_VERSION: + if ver not in (1, 2, 3, CKPT_VERSION): raise ValueError( f"Checkpoint version mismatch: file={ver}, " f"expected={CKPT_VERSION}. " @@ -146,6 +227,16 @@ def load_checkpoint(path, field, stepper, lbm_cfg, bodies): # Config compatibility check saved_cfg = json.loads(hf.attrs["config_json"]) + if "on_site_velocity" in (saved_cfg.get("y_wall_bc"), lbm_cfg.y_wall_bc): + for key in ("y_wall_bc", "velocity", "wall_normal_velocity"): + if saved_cfg.get(key) != getattr(lbm_cfg, key): + raise ValueError(f"Checkpoint on-site boundary mismatch on '{key}'") + if saved_cfg.get("wall_slot_width", 0) != lbm_cfg.wall_slot_width: + raise ValueError("Checkpoint on-site boundary mismatch on 'wall_slot_width'") + if saved_cfg.get("top_wall_bc", "velocity") != lbm_cfg.top_wall_bc: + raise ValueError("Checkpoint on-site boundary mismatch on 'top_wall_bc'") + if lbm_cfg.top_wall_bc == "reservoir" and saved_cfg.get("rho") != lbm_cfg.rho: + raise ValueError("Checkpoint reservoir boundary mismatch on 'rho'") for key in ("lattice_model", "nx", "ny", "nz", "store_precision", "streaming", "ddf_shifting"): saved_val = saved_cfg.get(key) @@ -170,31 +261,112 @@ def load_checkpoint(path, field, stepper, lbm_cfg, bodies): "the object count before loading, or rebuild the " "simulation with the correct number of objects.") - # Restore GPU arrays + # Read and validate everything before modifying host/device state. + obj_data = json.loads(hf.attrs.get("objects_json", "[]")) + if (len(obj_data) != bodies.count + or [od.get("obj_id") for od in obj_data] != list(range(bodies.count))): + raise ValueError("Checkpoint object IDs must match the dense current registry") + surface_specs = [asdict(obj.surface_actuation) if obj.surface_actuation is not None else None + for obj in bodies.objects] + if ver >= 4 and "surface_actuation_json" not in hf.attrs: + raise ValueError("Checkpoint is missing surface actuation specification") + saved_surface_specs = json.loads(hf.attrs.get( + "surface_actuation_json", json.dumps([None] * bodies.count))) + if saved_surface_specs != surface_specs: + raise ValueError("Checkpoint surface actuation mismatch") + domain_flags = None + if ver >= 3: + if "domain_flag" not in hf: + raise ValueError("Checkpoint is missing domain_flag") + domain_flags = hf["domain_flag"][:] + if domain_flags.shape != (field.n,) or domain_flags.dtype != np.uint16: + raise ValueError("Invalid checkpoint domain_flag shape/dtype") + geometry, expected_flags = _geometry_records(field, bodies, domain_flags) + if ver == 1: + _validate_legacy_objects(obj_data, bodies) + elif json.loads(hf.attrs["geometry_json"]) != geometry: + raise ValueError("Checkpoint object geometry/role mismatch; recreate the matching case") + + dtype = np.dtype(np.float16 if field.store_bytes == 2 else np.float32) ddf_buf = hf["ddf"][:] tmp_buf = hf["temp"][:] flag_buf = hf["flag"][:] + for name, array, shape, expected_dtype in ( + ("ddf", ddf_buf, (field.n * field.nq,), dtype), + ("temp", tmp_buf, (field.n * field.nq,), dtype), + ("flag", flag_buf, (field.n,), np.dtype(np.uint16)), + ): + if array.shape != shape or array.dtype != expected_dtype: + raise ValueError(f"Invalid checkpoint {name} shape/dtype") + if not np.array_equal(flag_buf, expected_flags): + raise ValueError("Checkpoint flags do not match current object/domain geometry") + if ver >= 3 and not np.array_equal(domain_flags, _domain_flags(field, bodies, flag_buf)): + raise ValueError("Checkpoint domain_flag is not a canonical object-free domain") + restored_step = int(hf.attrs["step_count"]) + if restored_step < 0: + raise ValueError("Checkpoint step_count must be nonnegative") + raw_params = json.loads(hf.attrs.get("runtime_params_json", "{}")) + params = {key: float(raw_params[key]) for key in ("omega", "fx", "fy", "fz") + if key in raw_params} + if not all(np.isfinite(value) for value in params.values()): + raise ValueError("Checkpoint runtime parameters must be finite") + if lbm_cfg.y_wall_bc == "on_site_velocity" and ( + bodies.count or any(params.get(k, 0.0) != 0 for k in ("fx", "fy", "fz"))): + raise ValueError("on_site_velocity does not support body forces or objects") + obs_accum_steps = int(hf.attrs.get("obs_accum_steps", 0)) + if obs_accum_steps < 0: + raise ValueError("Checkpoint obs_accum_steps must be nonnegative") + # Validate compact ownership before restoring any state or GPU arrays. + bodies.build_compact_lists_for(bodies.objects, domain_flags=expected_flags) + # A changed viscosity is an intentional parameter-study restart. + if saved_cfg.get("viscosity") != lbm_cfg.viscosity: + params["omega"] = float(lbm_cfg.omega) + params["n_objects"] = bodies.count + applied_action = hf["applied_action"][:] if "applied_action" in hf and ver >= 2 else None + if applied_action is not None and ( + applied_action.shape != bodies.action.shape or applied_action.dtype != np.float32): + raise ValueError("Invalid checkpoint applied_action shape/dtype") + obs = hf["obs"][:] if "obs" in hf else None + if obs is not None: + from ..cuda.obs_layout import obs_layout + size = obs_layout(lbm_cfg.dim, bodies.count).total_floats + if obs.shape != (size,) or obs.dtype != np.float32: + raise ValueError("Invalid checkpoint obs shape/dtype") + state_fields = ("x", "y", "z", "theta", "vx", "vy", "omega") + states = [{key: float(od.get(key, 0.0)) for key in state_fields} + for od in obj_data] + controls = [{key: float(od.get("control", {}).get(key, 0.0)) + for key in ("ax", "ay", "alpha")} for od in obj_data] + modes = [int(od.get("control", {}).get("mode", 0)) for od in obj_data] + if not all(np.isfinite(value) for record in states + controls for value in record.values()): + raise ValueError("Checkpoint object state/control must be finite") + for obj, state in zip(bodies.objects, states): + obj.validate_surface_actuation(lbm_cfg, state) + if (obj.surface_actuation is not None and applied_action is not None + and applied_action[(obj.obj_id + 1) * 3 * lbm_cfg.dim - 1] != 0): + raise ValueError("Checkpoint surface actuation requires zero applied omega") + for obj, state, control, mode in zip(bodies.objects, states, controls, modes): + for key, value in state.items(): + setattr(obj.state, key, value) + for key, value in control.items(): + setattr(obj.control, key, value) + obj.control.mode = mode + + field.flag = flag_buf.copy() + # Rebuild compact lists and action from the validated matching geometry. + # This also supports loading into a constructed but uninitialized case. + bodies.sync_to_gpu(field, rebuild_flags=False) + cuda.memcpy_htod(bodies.action_gpu, + bodies.action if applied_action is None else applied_action) cuda.memcpy_htod(field.ddf_gpu, ddf_buf) cuda.memcpy_htod(field.temp_gpu, tmp_buf) cuda.memcpy_htod(field.flag_gpu, flag_buf) - field.flag = flag_buf.copy() - - # Restore step count - stepper._step_count = int(hf.attrs["step_count"]) - - # Restore object states - obj_data = json.loads(hf.attrs.get("objects_json", "[]")) - for od in obj_data: - oid = od["obj_id"] - try: - obj = bodies.get(oid) - obj.state.x = od["x"] - obj.state.y = od["y"] - obj.state.z = float(od.get("z", 0.0)) - obj.state.theta = od["theta"] - obj.state.vx = od["vx"] - obj.state.vy = od["vy"] - obj.state.omega = od["omega"] - except KeyError: - pass # object not in current sim — skip + field.update_params(**params) + stepper._step_count = restored_step + field.invalidate_host_ddf_cache() + if obs is not None: + bodies.obs_pinned[:] = obs + cuda.memcpy_htod(bodies.obs_gpu, obs) + bodies._obs_accum_steps = obs_accum_steps + # v1 lacks observations; sync_to_gpu has explicitly reset the interval. diff --git a/src/CelerisLab/config.py b/src/CelerisLab/config.py index 4edb505..bff1c45 100644 --- a/src/CelerisLab/config.py +++ b/src/CelerisLab/config.py @@ -12,6 +12,7 @@ macro dicts consumed by `cuda/compiler_v2.py` for header generation. import json import os +import math from dataclasses import dataclass, field from typing import Any, Dict, List, Optional @@ -33,10 +34,16 @@ INLET_SCHEME_MAP = { "regularized": 3, } OUTLET_MAP = {"neq_extrap": 0, "zero_gradient": 1, "blended": 2} -Y_WALL_BC_MAP = {"bounce_back": 0, "free_slip": 1} +Y_WALL_BC_MAP = {"bounce_back": 0, "free_slip": 1, "on_site_velocity": 2} +TOP_WALL_BC_MAP = {"velocity": 0, "reservoir": 1} DTYPE_MAP = {"FP32": "float", "FP64": "double"} +def _float_literal(value: float) -> str: + """Emit a round-trip numeric literal; display precision must not alter physics.""" + return repr(float(value)) + "f" + + # --------------------------------------------------------------------------- # LBM config # --------------------------------------------------------------------------- @@ -80,6 +87,9 @@ class LBMConfig: inlet_regularized_neq_damp: float = 0.5 outlet_srt_neq_damp: float = 0.5 y_wall_bc: str = "bounce_back" + wall_normal_velocity: float = 0.0 + wall_slot_width: int = 0 + top_wall_bc: str = "velocity" omega_min: float = 0.01 omega_max: float = 1.99 @@ -101,21 +111,38 @@ class LBMConfig: return self.lattice_model == "D3Q19" def validate(self): - assert self.lattice_model in LATTICE_MODEL_MAP, ( - f"Unknown lattice_model: {self.lattice_model}. " - f"Supported: {list(LATTICE_MODEL_MAP.keys())}") + for name, allowed in ( + ("lattice_model", LATTICE_MODEL_MAP), ("collision", COLLISION_MAP), + ("streaming", STREAMING_MAP), ("y_wall_bc", Y_WALL_BC_MAP), + ("top_wall_bc", TOP_WALL_BC_MAP), + ("data_type", DTYPE_MAP), ("store_precision", PRECISION_MAP), + ("inlet_profile", INLET_MAP), ("inlet_scheme", INLET_SCHEME_MAP), + ("outlet_mode", OUTLET_MAP), + ): + if getattr(self, name) not in allowed: + raise ValueError(f"Unknown {name}: {getattr(self, name)!r}") lm = LATTICE_MODEL_MAP[self.lattice_model] - self.dim = lm["dim"] - self.nq = lm["nq"] - self.lattice_model_id = lm["id"] - assert self.collision in COLLISION_MAP, f"Unknown collision: {self.collision}" - assert self.streaming in STREAMING_MAP, f"Unknown streaming: {self.streaming}" - assert self.y_wall_bc in Y_WALL_BC_MAP, f"Unknown y_wall_bc: {self.y_wall_bc}" - assert self.data_type in DTYPE_MAP, f"Unknown data_type: {self.data_type}" - assert self.store_precision in PRECISION_MAP, ( - f"Unknown store_precision: {self.store_precision}") - assert self.inlet_scheme in INLET_SCHEME_MAP, ( - f"Unknown inlet_scheme: {self.inlet_scheme}") + self.dim, self.nq, self.lattice_model_id = lm["dim"], lm["nq"], lm["id"] + for name in ("nx", "ny", "nz", "threads_per_block"): + value = getattr(self, name) + if isinstance(value, bool) or not isinstance(value, int) or value <= 0: + raise ValueError(f"{name} must be a positive integer.") + if self.threads_per_block > 1024: + raise ValueError("threads_per_block must not exceed 1024.") + for name in ("viscosity", "velocity", "rho", "les_cs", "trt_magic_param", + "outlet_blend_alpha", "inlet_trt_neq_damp", + "inlet_regularized_neq_damp", "outlet_srt_neq_damp", + "omega_min", "omega_max", "wall_normal_velocity"): + if not math.isfinite(getattr(self, name)): + raise ValueError(f"{name} must be finite.") + if self.viscosity <= 0 or self.rho <= 0: + raise ValueError("viscosity and rho must be positive.") + if self.les_cs < 0 or self.trt_magic_param <= 0: + raise ValueError("les_cs must be nonnegative and trt_magic_param positive.") + if not 0 <= self.outlet_blend_alpha <= 1: + raise ValueError("outlet_blend_alpha must lie in [0, 1].") + if self.omega_min <= 0: + raise ValueError("omega_min must be positive.") if self.store_precision == "FP16C": raise ValueError( "store_precision='FP16C' is not supported in the current runtime path. " @@ -126,9 +153,22 @@ class LBMConfig: "Only data_type='FP32' is supported: CUDA kernels and LBMField " "paths use float32; LBtype is not wired through device code yet." ) - assert self.nx > 0 and self.ny > 0 and self.nz > 0 - if self.dim == 2: - assert self.nz == 1, "nz must be 1 for 2D" + if self.dim == 2 and self.nz != 1: + raise ValueError("nz must be 1 for 2D") + if (isinstance(self.wall_slot_width, bool) or not isinstance(self.wall_slot_width, int) + or self.wall_slot_width != 0 and not 2 <= self.wall_slot_width <= self.nx): + raise ValueError("wall_slot_width must be 0 (disabled) or an integer in [2, nx]") + if self.wall_slot_width and (self.y_wall_bc != "on_site_velocity" or self.velocity != 0): + raise ValueError("wall_slot_width requires on_site_velocity and velocity=0") + if self.top_wall_bc == "reservoir" and (self.y_wall_bc != "on_site_velocity" or self.velocity != 0): + raise ValueError("reservoir requires on_site_velocity and velocity=0") + if self.y_wall_bc == "on_site_velocity": + if (not self.is_d2q9 or self.collision != "SRT" + or self.streaming != "double_buffer" or self.store_precision != "FP32" + or self.ddf_shifting or self.les_enabled): + raise ValueError("on_site_velocity requires D2Q9/SRT/FP32/unshifted/double_buffer without LES") + if abs(self.wall_normal_velocity) >= 1: + raise ValueError("on_site_velocity requires abs(wall_normal_velocity) < 1") if not (0.0 <= self.inlet_trt_neq_damp <= 1.0): raise ValueError("method.inlet.trt_neq_damp must lie in [0, 1].") if not (0.0 <= self.inlet_regularized_neq_damp <= 1.0): @@ -150,6 +190,36 @@ class LBMConfig: f"[{self.omega_min:.6f}, {self.omega_max:.6f}]." ) + def validate_runtime(self): + """Validate the currently complete Simulation path, not local operators.""" + self.validate() + if self.is_d3q19: + raise NotImplementedError( + "D3Q19 Simulation is not supported yet: boundary/auxiliary kernels " + "and object projection are not dimension-complete. D3Q19 local " + "descriptors and collision operators remain available for research." + ) + if self.nx < 3 or self.ny < 4: + raise ValueError("Simulation requires nx >= 3 and ny >= 4: x ghost layers and distinct wall-adjacent rows.") + + def allocation_signature(self) -> tuple: + """Parameters that cannot change while field allocations/state are retained.""" + self.validate() + return (self.lattice_model, self.nx, self.ny, self.nz, self.data_type, + self.store_precision, bool(self.ddf_shifting), self.streaming, + (self.velocity, self.wall_normal_velocity, self.wall_slot_width, + self.top_wall_bc, self.rho if self.top_wall_bc == "reservoir" else None) + if self.y_wall_bc == "on_site_velocity" else None) + + @property + def top_wall_normal_velocity(self) -> float: + """Volume-balanced top velocity target, or reservoir initialization guess.""" + if self.wall_slot_width in (0, self.nx): + return self.wall_normal_velocity + weights = (max(0.0, 1.0-((2*x-(self.nx-1))/self.wall_slot_width)**2)**2 + for x in range(self.nx)) + return self.wall_normal_velocity * math.fsum(weights) / self.nx + def to_macros(self) -> Dict[str, Any]: """Return flat dict of macro_name → value for header generation.""" self.validate() @@ -161,28 +231,32 @@ class LBMConfig: "NZ": self.nz, "LATTICE_MODEL": self.lattice_model_id, "LBtype": DTYPE_MAP[self.data_type], - "VIS": f"{self.viscosity:.10f}", - "RHO": f"{self.rho:.1f}", - "U0": f"{self.velocity}", + "VIS": repr(float(self.viscosity)), + "RHO": repr(float(self.rho)), + "U0": repr(float(self.velocity)), "COLLISION_MODEL": COLLISION_MAP[self.collision], "STREAMING_MODEL": STREAMING_MAP[self.streaming], "STORE_PRECISION": PRECISION_MAP[self.store_precision], "USE_DDF_SHIFTING": int(self.ddf_shifting), "USE_LES": int(self.les_enabled), - "LES_CS": f"{self.les_cs:.6f}f", + "LES_CS": _float_literal(self.les_cs), "LES_CLOSED_FORM": int(self.les_closed_form), "INLET_PROFILE": INLET_MAP[self.inlet_profile], "INLET_SCHEME": INLET_SCHEME_MAP[self.inlet_scheme], "OUTLET_MODE": OUTLET_MAP[self.outlet_mode], - "OUTLET_BLEND_ALPHA": f"{self.outlet_blend_alpha:.3f}f", + "OUTLET_BLEND_ALPHA": _float_literal(self.outlet_blend_alpha), "OUTLET_BACKFLOW_CLAMP": int(self.outlet_backflow_clamp), "Y_WALL_BC": Y_WALL_BC_MAP[self.y_wall_bc], - "OMEGA_COLLISION_MIN": f"{self.omega_min:.2f}f", - "OMEGA_COLLISION_MAX": f"{self.omega_max:.3f}f", - "TRT_MAGIC_PARAM": f"{self.trt_magic_param:.6f}f", - "INLET_TRT_NEQ_DAMP": f"{self.inlet_trt_neq_damp:.4f}f", - "INLET_REGULARIZED_NEQ_DAMP": f"{self.inlet_regularized_neq_damp:.4f}f", - "OUTLET_SRT_NEQ_DAMP": f"{self.outlet_srt_neq_damp:.4f}f", + "WALL_NORMAL_VELOCITY": _float_literal(self.wall_normal_velocity), + "WALL_SLOT_WIDTH": self.wall_slot_width, + "TOP_WALL_NORMAL_VELOCITY": _float_literal(self.top_wall_normal_velocity), + "TOP_WALL_BC": TOP_WALL_BC_MAP[self.top_wall_bc], + "OMEGA_COLLISION_MIN": _float_literal(self.omega_min), + "OMEGA_COLLISION_MAX": _float_literal(self.omega_max), + "TRT_MAGIC_PARAM": _float_literal(self.trt_magic_param), + "INLET_TRT_NEQ_DAMP": _float_literal(self.inlet_trt_neq_damp), + "INLET_REGULARIZED_NEQ_DAMP": _float_literal(self.inlet_regularized_neq_damp), + "OUTLET_SRT_NEQ_DAMP": _float_literal(self.outlet_srt_neq_damp), } @@ -200,8 +274,11 @@ class BodyConfig: def _find_config(filename: str, explicit_path: Optional[str] = None) -> str: """Search for config file in standard locations.""" candidates = [] - if explicit_path: - candidates.append(explicit_path) + if explicit_path is not None: + explicit_path = os.fspath(explicit_path) + if not os.path.isfile(explicit_path): + raise FileNotFoundError(f"Explicit config path does not exist: {explicit_path}") + return os.path.abspath(explicit_path) env_dir = os.environ.get("CELERISLAB_CONFIG_DIR") if env_dir: candidates.append(os.path.join(env_dir, filename)) @@ -242,6 +319,9 @@ def load_lbm_config(path: Optional[str] = None) -> LBMConfig: inlet_regularized_neq_damp=m["inlet"].get("regularized_neq_damp", 0.5), outlet_srt_neq_damp=m["outlet"].get("srt_neq_damp", 0.5), y_wall_bc=m.get("y_wall_bc", "bounce_back"), + wall_normal_velocity=m.get("wall_normal_velocity", 0.0), + wall_slot_width=m.get("wall_slot_width", 0), + top_wall_bc=m.get("top_wall_bc", "velocity"), omega_min=m["omega_guard"]["min"], omega_max=m["omega_guard"]["max"], threads_per_block=c["threads_per_block"], @@ -255,4 +335,4 @@ def load_body_config(path: Optional[str] = None) -> BodyConfig: fpath = _find_config("config_body.json", path) with open(fpath) as f: d = json.load(f) - return BodyConfig(objects=d.get("objects", [])) \ No newline at end of file + return BodyConfig(objects=d.get("objects", [])) diff --git a/src/CelerisLab/configs/CONFIG.md b/src/CelerisLab/configs/CONFIG.md index e4fa6ef..c50411e 100644 --- a/src/CelerisLab/configs/CONFIG.md +++ b/src/CelerisLab/configs/CONFIG.md @@ -1,7 +1,7 @@ # CelerisLab 配置说明 **唯一配置入口：`config_lbm.json`** -Python `config.py` 只负责读取和校验，不是配置位置。 +Python `config.py` 只负责读取和校验，不是配置位置。显式配置路径不存在时会报错，不会退回默认配置；仅未指定路径时按环境目录、工作目录、包内配置顺序查找。 --- @@ -9,7 +9,7 @@ Python `config.py` 只负责读取和校验，不是配置位置。 | 字段 | 类型 | 默认 | 允许值 | 说明 | |------|------|------|--------|------| -| `lattice_model` | string | `"D2Q9"` | `D2Q9`, `D3Q19` | 格子模型，决定维度和速度数量 | +| `lattice_model` | string | `"D2Q9"` | `D2Q9`, `D3Q19` | 格子模型，决定维度和速度数量；`Simulation` 当前仅支持 D2Q9。D3Q19 保留供 descriptor / 局部算子研究，完整运行时在创建 GPU 上下文前明确拒绝 | | `nx` | int | 512 | >0 | x 方向格点数（流向），建议整除 `threads_per_block` 以减少尾块浪费 | | `ny` | int | 256 | >0 | y 方向格点数 | | `nz` | int | 1 | >0 | z 方向格点数，D2Q9 时必须为 1 | @@ -58,16 +58,16 @@ Python `config.py` 只负责读取和校验，不是配置位置。 ## 配置 → 代码映射 -``` **`lbm/kernels/config/*.h` 由编译器根据 `LBMConfig` 自动生成，请勿手改。** -仓库中可能检入过历史版本的 `config/*.h`，与当前默认 JSON **不一定一致**；若以 Python 入口运行仿真，`Simulation` / `compiler_v2.generate_config` 会在编译前重写这些头文件。若只用 `nvcc` 单独编译 `kernel_v2.cu` 且跳过 `generate_config`，则可能读到陈旧宏，属不受支持用法。 +仓库中检入的 `config/*.h` 可能是历史版本。`compiler_v2.generate_config` 会复制 kernel 源码到独立临时构建目录，并在副本中生成配置；不会改写包内头文件或 PTX。返回目录须显式传给 `compile_kernel(build_dir=...)`，编译结果与配置摘要均保存在该目录。`Simulation` 负责目录清理；直接使用编译器的调用者须调用 `cleanup_build_dir()`。跳过生成配置直接编译仓库中的 `kernel_v2.cu` 属不受支持用法。 +``` config_lbm.json ↓ load_lbm_config() config.py:LBMConfig (读取 + 校验) ↓ to_macros() -cuda/compiler_v2.py (生成 config/*.h) +cuda/compiler_v2.py (独立构建目录中生成 config/*.h) ↓ config_grid.h → NX, NY, NZ, NT, LATTICE_MODEL； DIM / NQ 由 LATTICE_MODEL 推导写入（与 grid.lattice_model 一致，非 JSON 顶层字段） @@ -100,4 +100,66 @@ step/one_step_*.cu → kernel 编排 D2Q9/D3Q19 的 flag 扩展位 `0x0800` 定义为 `FRC_REGION`，用于标记 force_region 类型物体的格子。 - ForceRegionKernel 不依赖 flag 位执行，flag 只用于后处理可视化。 -- `merge_flag_masks` 会正确合并 force_region 的 `FLUID | FRC_REGION` 标记。 \ No newline at end of file +- `merge_flag_masks` 会正确合并 force_region 的 `FLUID | FRC_REGION` 标记。 +# J1 on-site velocity planes (limited reference case) + +Set `method.y_wall_bc` to `"on_site_velocity"` and +`method.wall_normal_velocity` to V (default 0). This selects an x-periodic +domain with physical planes through fluid rows y=1 and y=NY−2, H=NY−3. +The bottom velocity is (0,V); the top is (`physics.velocity`,V). +These are constant **Cartesian** velocities: positive V injects through the +bottom and leaves through the top. The selected rows do not also receive +halfway reflection. Initialization uses a linear tangential profile and V. + +Supported only for D2Q9/SRT/FP32/unshifted/double_buffer without LES, objects +or body forces. Boundary targets and the mode are fixed at construction; +checkpoint restart requires matching targets. With the slot disabled this is an +entire-plane Couette-throughflow reference. The default +`bounce_back` mode and existing configurations retain their previous path. + +### Static centered bottom slot + +Optional `method.wall_slot_width` is an integer W: 0 (default) keeps the above +uniform planes; 2≤W≤NX enables a centered bottom slot and requires +`physics.velocity=0`. The same support restrictions apply. Width and amplitude +are fixed at construction and must match on checkpoint restart; older v3 +checkpoints without this field mean W=0. + +For W str: return os.path.join(_KERNEL_DIR, filename) @@ -37,18 +48,60 @@ _HEADER = "// AUTO-GENERATED by CelerisLab compiler – DO NOT EDIT MANUALLY\n" def _write(path: str, content: str): - os.makedirs(os.path.dirname(path), exist_ok=True) - tmp_path = f"{path}.tmp" - with open(tmp_path, "w") as f: - f.write(content) - os.replace(tmp_path, path) + # Headers and metadata belong to this caller's private build directory. + Path(path).parent.mkdir(parents=True, exist_ok=True) + Path(path).write_text(content, encoding="utf-8") -def generate_config(cfg: LBMConfig, n_objects: int = 0): - """Render all config/*.h files from *cfg*.""" - m = cfg.to_macros() - cdir = os.path.join(_KERNEL_DIR, "config") +def generate_config(cfg: LBMConfig, n_objects: int = 0, + *, build_root: Optional[str] = None) -> str: + """Snapshot kernels and render headers in a fresh private directory. + Returns its path. Pass it explicitly to compile_kernel for interleaved + builds. Legacy generate_config(); compile_kernel() remains context-local. + The caller owns cleanup_build_dir() once the loaded module is no longer + needed for recompilation/debugging. No package files are changed. + """ + macros = cfg.to_macros() + if isinstance(n_objects, bool) or not isinstance(n_objects, int) or n_objects < 0: + raise ValueError("n_objects must be a nonnegative integer") + if cfg.y_wall_bc == "on_site_velocity" and n_objects: + raise ValueError("on_site_velocity currently requires an object-free domain") + build_dir = tempfile.mkdtemp(prefix="celeris-build-", dir=build_root) + try: + shutil.copytree(_KERNEL_DIR, build_dir, dirs_exist_ok=True, + ignore=shutil.ignore_patterns("*.ptx", "*.cubin", "*.tmp")) + _render_config(macros, n_objects, os.path.join(build_dir, "config")) + digest = hashlib.sha256() + for path in sorted(Path(build_dir).rglob("*")): + if path.is_file(): + digest.update(str(path.relative_to(build_dir)).encode()) + digest.update(path.read_bytes()) + metadata = {"source_sha256": digest.hexdigest(), "config": asdict(cfg), + "n_objects": n_objects} + _write(os.path.join(build_dir, ".celeris-build"), json.dumps(metadata, sort_keys=True)) + except BaseException: + shutil.rmtree(build_dir) + raise + _last_build_dir.set(build_dir) + _last_ptx_path.set(None) + return build_dir + + +def cleanup_build_dir(build_dir: Optional[str]) -> None: + """Remove only a directory created by generate_config; safe to repeat.""" + if build_dir is None or not os.path.exists(build_dir): + return + path = Path(build_dir).resolve() + if path == Path(_KERNEL_DIR).resolve() or not (path / ".celeris-build").is_file(): + raise ValueError("Refusing to remove a directory not owned by the compiler") + shutil.rmtree(path) + if _last_build_dir.get() == str(build_dir): + _last_build_dir.set(None) + _last_ptx_path.set(None) + + +def _render_config(m: dict, n_objects: int, cdir: str): _write(os.path.join(cdir, "config_grid.h"), f"""{_HEADER}\ // Layer: Global / Grid #ifndef CELERIS_CONFIG_GRID_H @@ -122,6 +175,10 @@ def generate_config(cfg: LBMConfig, n_objects: int = 0): #define OUTLET_BLEND_ALPHA {m['OUTLET_BLEND_ALPHA']} #define OUTLET_BACKFLOW_CLAMP {m['OUTLET_BACKFLOW_CLAMP']} #define Y_WALL_BC {m['Y_WALL_BC']} +#define WALL_NORMAL_VELOCITY {m['WALL_NORMAL_VELOCITY']} +#define WALL_SLOT_WIDTH {m['WALL_SLOT_WIDTH']} +#define TOP_WALL_NORMAL_VELOCITY {m['TOP_WALL_NORMAL_VELOCITY']} +#define TOP_WALL_BC {m['TOP_WALL_BC']} #define OMEGA_COLLISION_MIN {m['OMEGA_COLLISION_MIN']} #define OMEGA_COLLISION_MAX {m['OMEGA_COLLISION_MAX']} @@ -184,36 +241,47 @@ def generate_config(cfg: LBMConfig, n_objects: int = 0): # --------------------------------------------------------------------------- # Compilation # --------------------------------------------------------------------------- -def compile_kernel(arch: str = "sm_70", - source: str = "kernel_v2.cu", - output: Optional[str] = None) -> str: - """Compile *source* to PTX. Returns path to the .ptx file. - - The translation unit includes ``config/*.h``. Those headers must match the - intended ``LBMConfig`` — call :func:`generate_config` before compiling. - Invoking ``nvcc`` on ``kernel_v2.cu`` alone can pick up stale checked-in - headers from ``lbm/kernels/config/``. - """ - src = kernel_path(source) - if output is None: - output = source.replace(".cu", ".ptx") - dst = kernel_path(output) - import shutil +def compile_kernel(arch: str = "sm_70", source: str = "kernel_v2.cu", + output: Optional[str] = None, *, + build_dir: Optional[str] = None) -> str: + """Compile a generated private source snapshot; never use package headers.""" + # A failed attempt must not make implicit load_module() reuse older code. + _last_ptx_path.set(None) + build_dir = build_dir or _last_build_dir.get() + if not build_dir or not (Path(build_dir) / ".celeris-build").is_file(): + raise RuntimeError("Call generate_config() and pass its build_dir before compiling") + root = Path(build_dir).resolve() + src = (root / source).resolve() + dst = (root / (output or source.replace(".cu", ".ptx"))).resolve() + if root not in src.parents or root not in dst.parents: + raise ValueError("Compiler source and output must stay inside the build directory") + if dst == src or dst.suffix != ".ptx": + raise ValueError("Compiler output must be a separate .ptx file") + dst.parent.mkdir(parents=True, exist_ok=True) nvcc = shutil.which("nvcc") or "/usr/local/cuda/bin/nvcc" - subprocess.run( - [nvcc, "-ptx", f"-arch={arch}", "--use_fast_math", src, "-o", dst], - check=True, - ) - return dst + fd, tmp = tempfile.mkstemp(prefix=dst.stem + "-", suffix=".ptx", dir=dst.parent) + os.close(fd) + command = [nvcc, "-ptx", f"-arch={arch}", "--use_fast_math", str(src), "-o", tmp] + try: + subprocess.run(command, check=True) + os.replace(tmp, dst) + metadata = json.loads((root / ".celeris-build").read_text()) + metadata.update(arch=arch, flags=["--use_fast_math"], + ptx_sha256=hashlib.sha256(dst.read_bytes()).hexdigest()) + _write(str(dst) + ".json", json.dumps(metadata, sort_keys=True)) + finally: + if os.path.exists(tmp): + os.unlink(tmp) + _last_ptx_path.set(str(dst)) + return str(dst) -# --------------------------------------------------------------------------- -# Module loading -# --------------------------------------------------------------------------- -def load_module(ptx_path: Optional[str] = None) -> cuda.Module: - """Load a compiled PTX as a PyCUDA module.""" - if ptx_path is None: - ptx_path = kernel_path("kernel_v2.ptx") +def load_module(ptx_path: Optional[str] = None): + """Load explicit PTX or the last PTX compiled in this Python context.""" + ptx_path = ptx_path or _last_ptx_path.get() + if not ptx_path: + raise RuntimeError("No compiled PTX; call compile_kernel() first") + import pycuda.driver as cuda return cuda.module_from_file(ptx_path) diff --git a/src/CelerisLab/cuda/context.py b/src/CelerisLab/cuda/context.py index 199f3c3..1354beb 100644 --- a/src/CelerisLab/cuda/context.py +++ b/src/CelerisLab/cuda/context.py @@ -1,41 +1,77 @@ -# CelerisLab/cuda/context.py -""" -CUDA device selection, context lifecycle, and device info queries. -""" +"""CUDA context ownership with balanced, operation-scoped activation.""" +from __future__ import annotations -import subprocess -from typing import Optional +from contextlib import contextmanager +from functools import wraps +from threading import RLock import pycuda.driver as cuda -class CudaContext: - """Manages a single CUDA device context. +def context_bound(function): + """Activate an object's ctx for one GPU operation (nested calls are safe). - Usage:: - - ctx = CudaContext(device_id=0) - # ... use PyCUDA ... - ctx.close() - - Or as a context-manager:: - - with CudaContext(0) as ctx: - ... + Low-level fields constructed under an explicitly current context may omit + ctx. Simulation attaches its owner before exposing the field to callers. """ + @wraps(function) + def wrapped(self, *args, **kwargs): + ctx = getattr(self, "ctx", None) + if ctx is None: + return function(self, *args, **kwargs) + with ctx.activate(): + return function(self, *args, **kwargs) + return wrapped + +class CudaContext: + """Own a context without leaving it on the caller's context stack. + + Use ``with ctx.activate():`` for GPU work. Independent instances may be + closed in any order, provided no activation of the closing owner is active. + ``with CudaContext(device_id) as ctx:`` activates and closes one owner. + """ def __init__(self, device_id: int = 0): + self._ctx = None + self._lock = RLock() + self._active = 0 + self._scopes = [] cuda.init() if device_id < 0 or device_id >= cuda.Device.count(): - raise ValueError( - f"device_id {device_id} invalid; " - f"{cuda.Device.count()} device(s) available." - ) + raise ValueError(f"device_id {device_id} invalid; {cuda.Device.count()} device(s) available.") self.device_id = device_id self.device = cuda.Device(device_id) self._ctx = self.device.make_context() + # make_context pushes immediately. Only balanced activate scopes may + # leave this owner current; close must never pop another instance. + self._ctx.pop() + + @contextmanager + def activate(self): + with self._lock: + if self._ctx is None: + raise RuntimeError("CUDA context is closed") + context = self._ctx + context.push() + self._active += 1 + try: + yield self + finally: + try: + context.pop() + finally: + with self._lock: + self._active -= 1 + + @property + def is_active(self) -> bool: + with self._lock: + return self._active > 0 + + @property + def is_closed(self) -> bool: + return self._ctx is None - # -- properties ---------------------------------------------------------- @property def name(self) -> str: return self.device.name() @@ -47,7 +83,6 @@ class CudaContext: @property def sm_arch(self) -> str: - """Return 'sm_XX' string for nvcc -arch.""" cc = self.device.compute_capability() return f"sm_{cc[0]}{cc[1]}" @@ -55,25 +90,36 @@ class CudaContext: def total_memory(self) -> int: return self.device.total_memory() - # -- lifecycle ----------------------------------------------------------- def close(self): - if self._ctx is not None: - self._ctx.pop() + with self._lock: + if self._ctx is None: + return + if self._active: + raise RuntimeError("Cannot close a CUDA context during an active operation") + self._ctx.detach() self._ctx = None def __enter__(self): + scope = self.activate() + scope.__enter__() + self._scopes.append(scope) return self def __exit__(self, *exc): - self.close() + self._scopes.pop().__exit__(*exc) + if not self._scopes: + self.close() def __del__(self): - self.close() + try: + self.close() + except Exception: + # Explicit close/context-manager exit is the supported lifecycle. + pass def detect_compute_capability(device_id: int = 0) -> str: - """Return 'sm_XX' for the given device without creating a persistent context.""" + """Query architecture without creating a persistent current context.""" cuda.init() - dev = cuda.Device(device_id) - cc = dev.compute_capability() + cc = cuda.Device(device_id).compute_capability() return f"sm_{cc[0]}{cc[1]}" diff --git a/src/CelerisLab/lbm/__init__.py b/src/CelerisLab/lbm/__init__.py index 25d25de..23d5225 100644 --- a/src/CelerisLab/lbm/__init__.py +++ b/src/CelerisLab/lbm/__init__.py @@ -1,16 +1,17 @@ -# CelerisLab/lbm/__init__.py -""" -Lattice Boltzmann Method module for fluid simulation. -""" +"""Lattice Boltzmann operators and lazily loaded GPU resources.""" +from importlib import import_module +from .descriptors import * # preserve existing flag/lattice exports -try: - from .descriptors import * # noqa: flags, lattice vectors - from .field import LBMField - from .stepper import LBMStepper - from .initializers import add_vortex +__all__ = ["LBMField", "LBMStepper", "add_vortex"] +_LAZY = {"LBMField": ("field", "LBMField"), + "LBMStepper": ("stepper", "LBMStepper"), + "add_vortex": ("initializers", "add_vortex")} - __all__ = ["LBMField", "LBMStepper", "add_vortex"] -except ImportError as e: - import warnings - warnings.warn(f"LBM module not fully available: {e}", ImportWarning) - __all__ = [] \ No newline at end of file + +def __getattr__(name): + if name not in _LAZY: + raise AttributeError(name) + module, symbol = _LAZY[name] + value = getattr(import_module("." + module, __name__), symbol) + globals()[name] = value + return value diff --git a/src/CelerisLab/lbm/curved_links.py b/src/CelerisLab/lbm/curved_links.py index 5fbcee3..139de98 100644 --- a/src/CelerisLab/lbm/curved_links.py +++ b/src/CelerisLab/lbm/curved_links.py @@ -23,6 +23,25 @@ def _free_alloc(ptr: Optional[cuda.DeviceAllocation]) -> None: ptr.free() +def _upload_columns(columns, stream=None): + """Copy low-frequency topology data; keep pinned staging alive to completion.""" + if stream is None: + for device, host in columns: + cuda.memcpy_htod(device, host) + return + staging = [] + try: + for device, host in columns: + pinned = cuda.pagelocked_empty(host.shape, host.dtype) + pinned[:] = host + staging.append(pinned) + cuda.memcpy_htod_async(device, pinned, stream) + finally: + # Returning pageable/temporary storage while DMA is active is invalid. + # Topology uploads are infrequent and intentionally complete on return. + stream.synchronize() + + @dataclass class CurvedLinkSoA: """Cut-link columns on host with matching ``*_gpu`` device buffers.""" @@ -56,6 +75,48 @@ class CurvedLinkSoA: count: int = 0 capacity: int = 0 + # Optional fixed surface velocity derived from true hit-center columns. + surface_velocity: Optional[np.ndarray] = None + aperture: Optional[np.ndarray] = None + surface_velocity_gpu: Optional[cuda.DeviceAllocation] = None + link_exchange_gpu: Optional[cuda.DeviceAllocation] = None + + def assign_surface(self, objects, cfg, flags): + """Keep the original eight-column geometry/checkpoint contract unchanged.""" + enabled = [obj for obj in objects if obj.surface_actuation is not None] + self.surface_velocity = self.aperture = None + if not enabled: + return + from .descriptors import FLUID, SOLID + if np.any((flags[self.fluid_idx] & FLUID) == 0) or np.any(flags[self.fluid_idx] & SOLID): + raise ValueError("Surface cut links must have active fluid receivers") + if len(np.unique(self.fluid_idx.astype(np.uint64)*cfg.nq+self.dir)) != len(self.dir): + raise ValueError("Surface boundary requires unique cut-link writers") + self.surface_velocity = np.zeros((len(self.dir), 2), np.float32) + self.aperture = np.zeros(len(self.dir), bool) + for obj in enabled: + obj.validate_surface_actuation(cfg) + selected = self.body_id == obj.obj_id + velocity, aperture = obj.surface_actuation.project(self.rx[selected], self.ry[selected]) + if not aperture.any(): + raise ValueError("Surface aperture hits no resolved cut links") + self.surface_velocity[selected] = velocity + self.aperture[selected] = aperture + + def surface_exchange(self, body_id): + """Latest boundary step's (Fx,Fy,Tz), split by aperture; not nozzle thrust. + + Caller synchronizes the CUDA stream/context. Rebuild/restart clears this + diagnostic; it is not the time-accumulated rigid-body observation. + """ + if self.link_exchange_gpu is None: + raise RuntimeError("No surface actuation diagnostic is allocated") + values = np.empty((self.count, 3), np.float32) + cuda.memcpy_dtoh(values, self.link_exchange_gpu) + body = self.body_id == body_id + return {name: values[body & mask].sum(axis=0, dtype=np.float64) + for name, mask in (("wall", ~self.aperture), ("aperture", self.aperture))} + def assign_host( self, fluid_idx: np.ndarray, @@ -125,16 +186,21 @@ class CurvedLinkSoA: raise ValueError( "CurvedLinkSoA column 'fallback_class' must use classes {0, 1}" ) + self.surface_velocity = self.aperture = None def upload(self, stream: Optional[cuda.Stream] = None) -> None: """Upload host columns to GPU with capacity-based allocation reuse. Args: - stream: If set, use async host-to-device copies on this stream. + stream: If set, use pinned copies ordered on this stream and wait + for completion before releasing staging memory. """ n = int(len(self.fluid_idx)) self.count = n if n == 0: + _free_alloc(self.surface_velocity_gpu) + _free_alloc(self.link_exchange_gpu) + self.surface_velocity_gpu = self.link_exchange_gpu = None return if n > self.capacity or self.fluid_idx_gpu is None: @@ -160,25 +226,24 @@ class CurvedLinkSoA: assert self.rz_gpu is not None assert self.fallback_class_gpu is not None - if stream is not None: - cuda.memcpy_htod_async(self.fluid_idx_gpu, self.fluid_idx, stream) - cuda.memcpy_htod_async(self.dir_gpu, self.dir, stream) - cuda.memcpy_htod_async(self.q_gpu, self.q, stream) - cuda.memcpy_htod_async(self.body_id_gpu, self.body_id, stream) - cuda.memcpy_htod_async(self.rx_gpu, self.rx, stream) - cuda.memcpy_htod_async(self.ry_gpu, self.ry, stream) - cuda.memcpy_htod_async(self.rz_gpu, self.rz, stream) - cuda.memcpy_htod_async( - self.fallback_class_gpu, self.fallback_class, stream) - else: - cuda.memcpy_htod(self.fluid_idx_gpu, self.fluid_idx) - cuda.memcpy_htod(self.dir_gpu, self.dir) - cuda.memcpy_htod(self.q_gpu, self.q) - cuda.memcpy_htod(self.body_id_gpu, self.body_id) - cuda.memcpy_htod(self.rx_gpu, self.rx) - cuda.memcpy_htod(self.ry_gpu, self.ry) - cuda.memcpy_htod(self.rz_gpu, self.rz) - cuda.memcpy_htod(self.fallback_class_gpu, self.fallback_class) + _upload_columns([ + (self.fluid_idx_gpu, self.fluid_idx), + (self.dir_gpu, self.dir), + (self.q_gpu, self.q), + (self.body_id_gpu, self.body_id), + (self.rx_gpu, self.rx), + (self.ry_gpu, self.ry), + (self.rz_gpu, self.rz), + (self.fallback_class_gpu, self.fallback_class), + ], stream) + _free_alloc(self.surface_velocity_gpu) + _free_alloc(self.link_exchange_gpu) + self.surface_velocity_gpu = self.link_exchange_gpu = None + if self.surface_velocity is not None: + self.surface_velocity_gpu = cuda.mem_alloc(self.surface_velocity.nbytes) + self.link_exchange_gpu = cuda.mem_alloc(n*3*4) + _upload_columns([(self.surface_velocity_gpu, self.surface_velocity), + (self.link_exchange_gpu, np.zeros((n, 3), np.float32))], stream) def free_gpu_columns(self) -> None: """Release device columns only (host arrays unchanged).""" @@ -190,6 +255,9 @@ class CurvedLinkSoA: _free_alloc(self.ry_gpu) _free_alloc(self.rz_gpu) _free_alloc(self.fallback_class_gpu) + _free_alloc(self.surface_velocity_gpu) + _free_alloc(self.link_exchange_gpu) + self.surface_velocity_gpu = self.link_exchange_gpu = None self.fluid_idx_gpu = None self.dir_gpu = None self.q_gpu = None @@ -229,7 +297,8 @@ class SensorSoA: """Reallocate GPU columns and upload. Args: - stream: If set, use async host-to-device copies on this stream. + stream: If set, use pinned copies ordered on this stream and wait + for completion before releasing staging memory. """ self.free_gpu_columns() n = int(len(self.cells)) @@ -239,12 +308,10 @@ class SensorSoA: self.cells_gpu = cuda.mem_alloc(int(self.cells.nbytes)) self.obj_id_gpu = cuda.mem_alloc(int(self.obj_id.nbytes)) - if stream is not None: - cuda.memcpy_htod_async(self.cells_gpu, self.cells, stream) - cuda.memcpy_htod_async(self.obj_id_gpu, self.obj_id, stream) - else: - cuda.memcpy_htod(self.cells_gpu, self.cells) - cuda.memcpy_htod(self.obj_id_gpu, self.obj_id) + _upload_columns([ + (self.cells_gpu, self.cells), + (self.obj_id_gpu, self.obj_id), + ], stream) def free_gpu_columns(self) -> None: _free_alloc(self.cells_gpu) @@ -281,7 +348,8 @@ class ForceRegionSoA: """Reallocate GPU columns and upload. Args: - stream: If set, use async host-to-device copies on this stream. + stream: If set, use pinned copies ordered on this stream and wait + for completion before releasing staging memory. """ self.free_gpu_columns() n = int(len(self.cells)) @@ -291,12 +359,10 @@ class ForceRegionSoA: self.cells_gpu = cuda.mem_alloc(int(self.cells.nbytes)) self.obj_id_gpu = cuda.mem_alloc(int(self.obj_id.nbytes)) - if stream is not None: - cuda.memcpy_htod_async(self.cells_gpu, self.cells, stream) - cuda.memcpy_htod_async(self.obj_id_gpu, self.obj_id, stream) - else: - cuda.memcpy_htod(self.cells_gpu, self.cells) - cuda.memcpy_htod(self.obj_id_gpu, self.obj_id) + _upload_columns([ + (self.cells_gpu, self.cells), + (self.obj_id_gpu, self.obj_id), + ], stream) def free_gpu_columns(self) -> None: _free_alloc(self.cells_gpu) @@ -341,7 +407,8 @@ class EsoPullCurvedBuffer: """Reallocate GPU columns and upload. Args: - stream: If set, use async host-to-device copies on this stream. + stream: If set, use pinned copies ordered on this stream and wait + for completion before releasing staging memory. """ self.free_gpu_columns() n = int(len(self.corr_value)) @@ -352,14 +419,11 @@ class EsoPullCurvedBuffer: self.corr_value_gpu = cuda.mem_alloc(int(self.corr_value.nbytes)) self.target_cell_gpu = cuda.mem_alloc(int(self.target_cell.nbytes)) self.target_slot_gpu = cuda.mem_alloc(int(self.target_slot.nbytes)) - if stream is not None: - cuda.memcpy_htod_async(self.corr_value_gpu, self.corr_value, stream) - cuda.memcpy_htod_async(self.target_cell_gpu, self.target_cell, stream) - cuda.memcpy_htod_async(self.target_slot_gpu, self.target_slot, stream) - else: - cuda.memcpy_htod(self.corr_value_gpu, self.corr_value) - cuda.memcpy_htod(self.target_cell_gpu, self.target_cell) - cuda.memcpy_htod(self.target_slot_gpu, self.target_slot) + _upload_columns([ + (self.corr_value_gpu, self.corr_value), + (self.target_cell_gpu, self.target_cell), + (self.target_slot_gpu, self.target_slot), + ], stream) def free_gpu_columns(self) -> None: _free_alloc(self.corr_value_gpu) diff --git a/src/CelerisLab/lbm/field.py b/src/CelerisLab/lbm/field.py index 9d64cfe..5a986ed 100644 --- a/src/CelerisLab/lbm/field.py +++ b/src/CelerisLab/lbm/field.py @@ -21,6 +21,7 @@ import numpy as np import pycuda.driver as cuda from ..config import LBMConfig +from ..cuda.context import context_bound from .curved_links import CurvedLinkSoA, SensorSoA, ForceRegionSoA, EsoPullCurvedBuffer from .descriptors import FLUID, SOLID, BC_WALL, BC_INLET, BC_OUTLET @@ -122,6 +123,8 @@ class LBMField: if self.cfg.is_d2q9: flag[0:self.nx] = SOLID | BC_WALL flag[(self.ny - 1) * self.nx:self.ny * self.nx] = SOLID | BC_WALL + if self.cfg.y_wall_bc == "on_site_velocity": + return flag # x periodic; physical y planes are the adjacent fluid rows. flag[::self.nx] = SOLID | BC_INLET flag[self.nx - 1::self.nx] = SOLID | BC_OUTLET flag[0] = SOLID | BC_WALL @@ -154,6 +157,7 @@ class LBMField: return self.force_regions.count # -- DDF download with precision decode ---------------------------------- + @context_bound 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.""" @@ -190,6 +194,7 @@ class LBMField: # See load_physical_node_esopull / store_physical_node_esopull. self._host_ddf_step = int(step_id) if step_id is not None else None + @context_bound 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: @@ -201,6 +206,7 @@ class LBMField: return self._cached_w # -- d_params upload ----------------------------------------------------- + @context_bound def _upload_params(self): """Pack LBMParams struct and upload to __constant__ d_params. @@ -221,11 +227,16 @@ class LBMField: ptr, size = self.module.get_global("d_params") cuda.memcpy_htod(ptr, data) + @context_bound def update_params(self, **kwargs): """Re-upload d_params after changing runtime-adjustable values. Accepted keys: omega, fx, fy, fz, n_objects. """ + if self.cfg.y_wall_bc == "on_site_velocity" and ( + any(float(kwargs.get(k, 0.0)) != 0 for k in ("fx", "fy", "fz")) + or int(kwargs.get("n_objects", 0)) != 0): + raise ValueError("on_site_velocity does not support body forces or objects") self._runtime_params.update({ "omega": float(kwargs.get("omega", self._runtime_params["omega"])), "fx": float(kwargs.get("fx", self._runtime_params["fx"])), @@ -234,46 +245,34 @@ class LBMField: "n_objects": int(kwargs.get("n_objects", self._runtime_params["n_objects"])), }) - p = self._runtime_params - fmt = "=f fff I" - data = struct.pack( - fmt, - p["omega"], - p["fx"], p["fy"], p["fz"], - p["n_objects"], - ) - ptr, _ = self.module.get_global("d_params") - cuda.memcpy_htod(ptr, data) + self._upload_params() # -- Host ↔ GPU transfers ------------------------------------------------ - def upload_ddf(self): - upload_data = self.ddf + def _encode_ddf(self, physical: np.ndarray) -> np.ndarray: + """Encode a Q-major physical array (full field or selected cells).""" + encoded = np.asarray(physical, dtype=np.float32).copy() if getattr(self.cfg, 'ddf_shifting', False): - # Host stores physical DDF; GPU stores shifted DDF (f - w). - w = self._read_lattice_weights() - upload_data = upload_data.copy() - ddf_2d = upload_data.reshape(self.nq, -1) - for i in range(self.nq): - ddf_2d[i] -= w[i] - + encoded.reshape(self.nq, -1)[:] -= self._read_lattice_weights()[:, None] if self.store_bytes == 2: if self.store_dtype != np.float16: - raise RuntimeError( - "Invalid 2-byte DDF storage dtype; expected float16 for FP16S." - ) - buf = (upload_data * 32768.0).astype(self.store_dtype) - cuda.memcpy_htod(self.ddf_gpu, buf) - cuda.memcpy_htod(self.temp_gpu, buf) - else: - cuda.memcpy_htod(self.ddf_gpu, upload_data) - cuda.memcpy_htod(self.temp_gpu, upload_data) + raise RuntimeError("Invalid 2-byte DDF storage dtype; expected float16 for FP16S.") + encoded *= np.float32(32768.0) + return np.ascontiguousarray(encoded, dtype=self.store_dtype) + + @context_bound + def upload_ddf(self): + upload_data = self._encode_ddf(self.ddf) + cuda.memcpy_htod(self.ddf_gpu, upload_data) + cuda.memcpy_htod(self.temp_gpu, upload_data) self.invalidate_host_ddf_cache() # The canonical download_ddf is defined above in the class body. + @context_bound def upload_flags(self): cuda.memcpy_htod(self.flag_gpu, self.flag) + @context_bound def upload_compact_lists(self, stream: Optional[cuda.Stream] = None) -> None: """Upload compact lists (curved, sensor, force_region) to GPU. @@ -294,6 +293,7 @@ class LBMField: self.esopull_curved.upload(stream=stream) # -- Read lattice descriptors from CUDA module --------------------------- + @context_bound 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: @@ -318,17 +318,20 @@ class LBMField: return self._cached_vectors_2d # -- Macroscopic field extraction ---------------------------------------- + @context_bound def get_macroscopic(self, *, step_id: int | None = None): """Download DDF and compute rho, ux, uy [, uz]. - For double_buffer streaming: host-side decode from raw DDF. - For esopull streaming: GPU-side MacroscopicEsoPullKernel. + Double-buffer returns moments of the current spatial postcollision + field. EsoPull returns moments of the completed-count next-pull + (precollision) state via MacroscopicEsoPullKernel. These are different + numerical phases, so pointwise equality across streaming modes is not + a general contract. .. note:: - In EsoPull mode, ``get_ddf()`` returns backing-layout raw storage, - NOT physical DDF. Only ``get_macroscopic()`` is semantically correct. - Verified: Kan99b K2 (200k steps, MRT) -- bit-identical to - double-buffer for all metrics (St, mean_CL, mean_CD, amp_CL, amp_CD). + In EsoPull mode, ``get_ddf()`` returns backing-layout values after + precision and shifting decode, not spatial physical DDF. Use this + macroscopic path for the documented next-pull interpretation. """ if getattr(self.cfg, 'streaming', 'double_buffer') == 'esopull': return self.get_macroscopic_esopull(step_id=step_id) @@ -358,6 +361,7 @@ class LBMField: raise ValueError(f"Unsupported lattice_model: {self.cfg.lattice_model}") + @context_bound def get_macroscopic_esopull(self, *, step_id: int | None = None) -> dict: """Compute macroscopic field from EsoPull backing layout on GPU. @@ -366,11 +370,10 @@ class LBMField: .. note:: Only supports D2Q9. D3Q19 raises NotImplementedError. - This is the ONLY correct path for macroscopic diagnostics in EsoPull - mode. Host-side raw DDF download (get_ddf) is NOT physical DDF - when streaming is esopull -- the backing layout scatters paired - directions across neighbor cells. - Verified: Kan99b K2 -- bit-identical to double-buffer. + At the completed step count this reads the next-pull/precollision + state. Host get_ddf decodes precision and shifting but retains the + scattered backing layout; it is not spatial physical DDF. This + phase differs from the double-buffer postcollision diagnostic. Args: step_id: Current step count used for EsoPull parity decoding. @@ -408,6 +411,7 @@ class LBMField: # LBMField. # -- EsoPull sync decode/encode (compact list) -------------------------- + @context_bound def esopull_sync_decode(self, patch_indices, step_count, stream=None): """Decode a subset of cells from EsoPull backing to physical DDF. @@ -420,9 +424,11 @@ class LBMField: (scratch_gpu, scratch_host, idx_gpu): scratch buffer handles. """ n = len(patch_indices) + if n == 0: + return None, np.empty(0, dtype=np.float32), None nq = self.nq scratch_gpu = cuda.mem_alloc(n * nq * 4) - scratch_host = np.empty(n * nq, dtype=np.float32) + scratch_host = cuda.pagelocked_empty(n * nq, np.float32) idx_gpu = cuda.mem_alloc(n * 4) cuda.memcpy_htod(idx_gpu, np.asarray(patch_indices, dtype=np.uint32)) @@ -433,20 +439,27 @@ class LBMField: block=(tpb, 1, 1), grid=((n + tpb - 1) // tpb, 1, 1), stream=stream) - cuda.memcpy_dtoh_async(scratch_host, scratch_gpu, stream) - stream.synchronize() + if stream is None: + cuda.memcpy_dtoh(scratch_host, scratch_gpu) + else: + cuda.memcpy_dtoh_async(scratch_host, scratch_gpu, stream) + stream.synchronize() return scratch_gpu, scratch_host, idx_gpu + @context_bound def esopull_sync_encode(self, scratch_gpu, scratch_host, idx_gpu, n, step_count, stream=None): - """Upload patched physical DDF, collide + encode to EsoPull layout. + """Restore the same-count next-pull physical state, without collision. - The backing layout stores post-collision values. This method - applies collision to the patched physical DDF before encoding, - matching the store_f_esopull convention at parity *step_count*. - Frees scratch_gpu and idx_gpu after encoding. + Scratch is cell-major physical FP32. The kernel applies shifting and + storage encoding and is the exact inverse of same-count decode. """ - cuda.memcpy_htod_async(scratch_gpu, scratch_host, stream) + if n == 0: + return + if stream is None: + cuda.memcpy_htod(scratch_gpu, scratch_host) + else: + cuda.memcpy_htod_async(scratch_gpu, scratch_host, stream) fn = self.module.get_function("EncodePhysicalToCells") tpb = self.cfg.threads_per_block @@ -455,15 +468,45 @@ class LBMField: block=(tpb, 1, 1), grid=((n + tpb - 1) // tpb, 1, 1), stream=stream) + # Keep host/device scratch alive until the async upload and kernel end. + if stream is None: + cuda.Context.synchronize() + else: + stream.synchronize() scratch_gpu.free() idx_gpu.free() + self.invalidate_host_ddf_cache() - def snapshot(self): - self.download_ddf(force=True) - self._ddf_snap = self.ddf.copy() + @context_bound + def snapshot(self, *, step_id: int | None = None): + """Save a field-only snapshot; simulation time and controls are unchanged.""" + if self.cfg.streaming == "esopull": + if step_id is None: + raise ValueError("EsoPull field snapshots require the current step_id") + indices = np.arange(self.n, dtype=np.uint32) + scratch, physical, idx_gpu = self.esopull_sync_decode(indices, step_id) + self._ddf_snap = physical.reshape(self.n, self.nq).T.copy().reshape(-1) + scratch.free() + idx_gpu.free() + else: + self.download_ddf(force=True) + self._ddf_snap = self.ddf.copy() - def restore(self): + @context_bound + def restore(self, *, step_id: int | None = None): + """Restore the saved field at the current clock, including EsoPull parity.""" if self._ddf_snap is None: raise RuntimeError("No snapshot to restore") + if self.cfg.streaming == "esopull" and step_id is None: + raise ValueError("EsoPull field snapshots require the current step_id") self.ddf = self._ddf_snap.copy() - self.upload_ddf() + if self.cfg.streaming == "esopull": + indices = np.arange(self.n, dtype=np.uint32) + scratch_gpu = cuda.mem_alloc(self.n * self.nq * 4) + physical = cuda.pagelocked_empty(self.n * self.nq, np.float32) + physical[:] = self.ddf.reshape(self.nq, self.n).T.reshape(-1) + idx_gpu = cuda.mem_alloc(indices.nbytes) + cuda.memcpy_htod(idx_gpu, indices) + self.esopull_sync_encode(scratch_gpu, physical, idx_gpu, self.n, step_id) + else: + self.upload_ddf() diff --git a/src/CelerisLab/lbm/initializers.py b/src/CelerisLab/lbm/initializers.py index 9c8fb2b..305f2ea 100644 --- a/src/CelerisLab/lbm/initializers.py +++ b/src/CelerisLab/lbm/initializers.py @@ -10,7 +10,7 @@ import numpy as np from scipy.special import jv, expi from typing import Tuple -import pycuda.driver as cuda +from .equilibrium import compute_feq_d2q9 def add_vortex(field, center: Tuple[float, float], @@ -22,6 +22,12 @@ def add_vortex(field, center: Tuple[float, float], """ if not field.cfg.is_d2q9: raise NotImplementedError("Vortex init only for D2Q9") + if field.cfg.streaming != "double_buffer": + raise NotImplementedError( + "Vortex initialization requires double_buffer physical DDF; " + "raw EsoPull storage cannot be edited without its time/parity.") + if not np.isfinite(radius) or radius <= 0: + raise ValueError("Vortex radius must be finite and positive") if vortex_type not in ("lamb", "oseen", "taylor"): raise ValueError(f"Unknown vortex type: {vortex_type}") @@ -49,7 +55,6 @@ def add_vortex(field, center: Tuple[float, float], # Read lattice descriptors from CUDA module (single source of truth) cx, cy = field._read_lattice_vectors() - w = field._read_lattice_weights() # Compute current macroscopic from DDF rho_old = np.sum(f, axis=0) @@ -68,10 +73,7 @@ def add_vortex(field, center: Tuple[float, float], for j in range(1, ny - 1): for i in range(1, nx - 1): u, v, p = float(u_new[j, i]), float(v_new[j, i]), float(p_new[j, i]) - for e in range(nq): - eu = cx[e] * u + cy[e] * v - u2 = u * u + v * v - f[e, j, i] = w[e] * (3*p + 3*eu + 4.5*eu*eu - 1.5*u2) + f[:, j, i] = compute_feq_d2q9(3*p, u, v) field.ddf = f.reshape(-1) field.upload_ddf() @@ -103,16 +105,23 @@ def _lamb_vortex(r, theta, radius, strength): def _oseen_vortex(r, theta, radius, strength): kappa = 2 * np.pi * radius**2 * strength - r_safe = np.where(r > 1e-12, r, 1e-12) - exp_term = 1 - np.exp(-4 * r**2 / radius**2) - u_vor = -kappa / (2*np.pi*r_safe) * exp_term * np.sin(theta) - v_vor = kappa / (2*np.pi*r_safe) * exp_term * np.cos(theta) zeta = 4 * r**2 / radius**2 - p_vor = ( - -kappa**2 / (8*np.pi**2 * r_safe**2) - * (-2*zeta*(expi(-zeta) - expi(-2*zeta)) - + (1 - np.exp(-zeta))**2) - ) + exp_term = -np.expm1(-zeta) + speed = kappa / (2*np.pi) * np.divide( + exp_term, r, out=np.zeros_like(r, dtype=float), where=r != 0) + u_vor = -speed * np.sin(theta) + v_vor = speed * np.cos(theta) + # Factor out zeta/r**2 before evaluating the pressure. At the center, + # Ei(-zeta)-Ei(-2*zeta) tends to -log(2), not NaN. The expansion also + # avoids cancellation for near-center samples and preserves the same model. + pressure_shape = np.empty_like(zeta, dtype=float) + small = zeta < 1e-5 + z = zeta[small] + pressure_shape[small] = 2*np.log(2.0) - z + z*z/2 - 7*z*z*z/36 + z = zeta[~small] + pressure_shape[~small] = ( + -2*(expi(-z) - expi(-2*z)) + np.expm1(-z)**2 / z) + p_vor = -kappa**2 / (2*np.pi**2*radius**2) * pressure_shape return u_vor, v_vor, p_vor @@ -127,4 +136,4 @@ def _taylor_vortex(r, theta, radius, strength, nu): * np.exp(-2 * r**2 / radius**2) / (np.pi**2 * radius**6) ) - return u_vor, v_vor, p_vor \ No newline at end of file + return u_vor, v_vor, p_vor diff --git a/src/CelerisLab/lbm/kernels/boundary/bounce_back.cuh b/src/CelerisLab/lbm/kernels/boundary/bounce_back.cuh index 99384f1..69210ba 100644 --- a/src/CelerisLab/lbm/kernels/boundary/bounce_back.cuh +++ b/src/CelerisLab/lbm/kernels/boundary/bounce_back.cuh @@ -5,10 +5,7 @@ // We replace those directions with the opposite-direction population at the // SAME node from the INPUT buffer (half-way BB). // -// In EsoPull mode (t > 0), reads via load_physical_dir_esopull to resolve -// the opposite direction's value from the EsoPull single-buffer layout. -// For Kan99b (cylinder far from walls), this is correct in practice. -// A more rigorous pre_streaming fix is pending further investigation. +// EsoPull stages its pre-streaming donors in PrepareEsoPullBoundaries. // ============================================================================ #ifndef CELERIS_BOUNDARY_BOUNCE_BACK_CUH @@ -21,33 +18,19 @@ __device__ inline void apply_wall_bb_d2q9_y_pull(unsigned int y, float* __restrict__ f, const fpxx* __restrict__ fi_in, - unsigned long k, - const unsigned long* j = nullptr, - unsigned long t = 0ul) + unsigned long k) { if (y == 1) { // Directions sourced from y=0 wall in pull step: +y, +x+y, -x+y. - if (t == 0ul) { - f[3] = load_ddf(fi_in, index_f(k, 4u)); - f[5] = load_ddf(fi_in, index_f(k, 6u)); - f[8] = load_ddf(fi_in, index_f(k, 7u)); - } else { - f[3] = load_physical_dir_esopull(k, 4u, fi_in, j, t); // N ← S - f[5] = load_physical_dir_esopull(k, 6u, fi_in, j, t); // NE ← SW - f[8] = load_physical_dir_esopull(k, 7u, fi_in, j, t); // NW ← SE - } + f[3] = load_ddf(fi_in, index_f(k, 4u)); + f[5] = load_ddf(fi_in, index_f(k, 6u)); + f[8] = load_ddf(fi_in, index_f(k, 7u)); } else if (y == (unsigned int)(NY - 2)) { // Directions sourced from y=NY-1 wall in pull step: -y, -x-y, +x-y. - if (t == 0ul) { - f[4] = load_ddf(fi_in, index_f(k, 3u)); - f[6] = load_ddf(fi_in, index_f(k, 5u)); - f[7] = load_ddf(fi_in, index_f(k, 8u)); - } else { - f[4] = load_physical_dir_esopull(k, 3u, fi_in, j, t); // S ← N - f[6] = load_physical_dir_esopull(k, 5u, fi_in, j, t); // SW ← NE - f[7] = load_physical_dir_esopull(k, 8u, fi_in, j, t); // SE ← NW - } + f[4] = load_ddf(fi_in, index_f(k, 3u)); + f[6] = load_ddf(fi_in, index_f(k, 5u)); + f[7] = load_ddf(fi_in, index_f(k, 8u)); } } @@ -62,31 +45,17 @@ __device__ inline void apply_wall_bb_d2q9_y_pull(unsigned int y, __device__ inline void apply_wall_freeslip_d2q9_y_pull(unsigned int y, float* __restrict__ f, const fpxx* __restrict__ fi_in, - unsigned long k, - const unsigned long* j = nullptr, - unsigned long t = 0ul) + unsigned long k) { if (y == 1) { - if (t == 0ul) { - f[3] = load_ddf(fi_in, index_f(k, 4u)); - f[5] = load_ddf(fi_in, index_f(k, 7u)); - f[8] = load_ddf(fi_in, index_f(k, 6u)); - } else { - f[3] = load_physical_dir_esopull(k, 4u, fi_in, j, t); // N ← S (normal) - f[5] = load_physical_dir_esopull(k, 7u, fi_in, j, t); // NE ← SE (tangential mirror) - f[8] = load_physical_dir_esopull(k, 6u, fi_in, j, t); // NW ← SW (tangential mirror) - } + f[3] = load_ddf(fi_in, index_f(k, 4u)); + f[5] = load_ddf(fi_in, index_f(k, 7u)); + f[8] = load_ddf(fi_in, index_f(k, 6u)); } else if (y == (unsigned int)(NY - 2)) { - if (t == 0ul) { - f[4] = load_ddf(fi_in, index_f(k, 3u)); - f[6] = load_ddf(fi_in, index_f(k, 8u)); - f[7] = load_ddf(fi_in, index_f(k, 5u)); - } else { - f[4] = load_physical_dir_esopull(k, 3u, fi_in, j, t); // S ← N (normal) - f[6] = load_physical_dir_esopull(k, 8u, fi_in, j, t); // SW ← NW (tangential mirror) - f[7] = load_physical_dir_esopull(k, 5u, fi_in, j, t); // SE ← NE (tangential mirror) - } + f[4] = load_ddf(fi_in, index_f(k, 3u)); + f[6] = load_ddf(fi_in, index_f(k, 8u)); + f[7] = load_ddf(fi_in, index_f(k, 5u)); } } diff --git a/src/CelerisLab/lbm/kernels/boundary/curved_boundary.cuh b/src/CelerisLab/lbm/kernels/boundary/curved_boundary.cuh index b287c5d..6b566f8 100644 --- a/src/CelerisLab/lbm/kernels/boundary/curved_boundary.cuh +++ b/src/CelerisLab/lbm/kernels/boundary/curved_boundary.cuh @@ -48,26 +48,27 @@ __device__ __forceinline__ float bouzidi_linear_moving_correction( { if (fallback_class != CURVED_FALLBACK_BOUZIDI) { // Fallback replaces the q < 1/2 branch when its donor is illegal. - // For D2Q9, alpha_i = 3 w_i, so the moving-wall addend is - // +2 * alpha_i * (c_i . u_w) in the q < 1/2 branch. - return 2.0f * alpha_ci_dot_uw; + // Here dir points from fluid to solid and alpha_i = 3 w_i rho. + // The opposite-direction reflected population receives + // -2 * alpha_i * (c_i . u_w) in the q < 1/2 branch. + return -2.0f * alpha_ci_dot_uw; } if (q < 0.5f) { - return 2.0f * alpha_ci_dot_uw; + return -2.0f * alpha_ci_dot_uw; } - return (alpha_ci_dot_uw / q); + return -(alpha_ci_dot_uw / q); } __device__ __forceinline__ float compute_moving_wall_addend( unsigned int dir, float q, unsigned char fallback_class, - float Uw, float Vw) + float Uw, float Vw, float rho = 1.0f) { #if DIM == 2 const float ci_dot_uw = (float)d_cx[dir] * Uw + (float)d_cy[dir] * Vw; #elif DIM == 3 const float ci_dot_uw = (float)d_cx[dir] * Uw + (float)d_cy[dir] * Vw + (float)d_cz[dir] * 0.0f; #endif - const float alpha_ci_dot_uw = 3.0f * d_w[dir] * ci_dot_uw; + const float alpha_ci_dot_uw = 3.0f * d_w[dir] * rho * ci_dot_uw; return bouzidi_linear_moving_correction(q, fallback_class, alpha_ci_dot_uw); } @@ -81,7 +82,7 @@ __device__ inline void apply_bouzidi_link( float Uw, float Vw, float rx, float ry, unsigned char fallback_class, - float* __restrict__ obs, int body_id) + float* __restrict__ obs, int body_id, float* link_exchange = nullptr) { const unsigned int dir_opp = (unsigned int)opp_dir((int)dir); @@ -108,13 +109,18 @@ __device__ inline void apply_bouzidi_link( } else { f_reflected = compute_bouzidi_reflection(f_toward, 0.0f, 0.0f, q, fallback_class); } - f_reflected += compute_moving_wall_addend(dir, q, fallback_class, Uw, Vw); + float rho_wall = USE_DDF_SHIFTING ? 1.0f : 0.0f; + if (Uw != 0.0f || Vw != 0.0f) { + for (unsigned int i = 0; i < NQ; ++i) + rho_wall += load_ddf(fi, index_f(k_f, i)); + } + f_reflected += compute_moving_wall_addend(dir, q, fallback_class, Uw, Vw, rho_wall); // Write to the solid source node so the subsequent pull step loads the // corrected incoming population into the adjacent fluid node. store_ddf(fi, index_f(k_s, dir_opp), f_reflected); - if (obs != nullptr) { + if (obs != nullptr || link_exchange != nullptr) { #if USE_DDF_SHIFTING const float w_i = (float)d_w[dir]; const float fx = (float)d_cx[dir] * ((f_toward + w_i) + (f_reflected + w_i)); @@ -123,10 +129,15 @@ __device__ inline void apply_bouzidi_link( const float fx = (float)d_cx[dir] * (f_toward + f_reflected); const float fy = (float)d_cy[dir] * (f_toward + f_reflected); #endif - atomicAdd(&obs[obs_force_index(body_id, 0)], fx); - atomicAdd(&obs[obs_force_index(body_id, 1)], fy); const float tz = rx * fy - ry * fx; - atomicAdd(&obs[obs_torque_index(body_id, 0)], tz); + if (link_exchange != nullptr) { + link_exchange[0] = fx; link_exchange[1] = fy; link_exchange[2] = tz; + } + if (obs != nullptr) { + atomicAdd(&obs[obs_force_index(body_id, 0)], fx); + atomicAdd(&obs[obs_force_index(body_id, 1)], fy); + atomicAdd(&obs[obs_torque_index(body_id, 0)], tz); + } } } #endif // DIM == 2 @@ -170,7 +181,12 @@ __device__ inline void apply_bouzidi_link( } else { f_reflected = compute_bouzidi_reflection(f_toward, 0.0f, 0.0f, q, fallback_class); } - f_reflected += compute_moving_wall_addend(dir, q, fallback_class, Uw, Vw); + float rho_wall = USE_DDF_SHIFTING ? 1.0f : 0.0f; + if (Uw != 0.0f || Vw != 0.0f) { + for (unsigned int i = 0; i < NQ; ++i) + rho_wall += load_ddf(fi, index_f(k_f, i)); + } + f_reflected += compute_moving_wall_addend(dir, q, fallback_class, Uw, Vw, rho_wall); store_ddf(fi, index_f(k_s, dir_opp), f_reflected); diff --git a/src/CelerisLab/lbm/kernels/boundary/inlet/common.cuh b/src/CelerisLab/lbm/kernels/boundary/inlet/common.cuh index edba38e..cd37f5d 100644 --- a/src/CelerisLab/lbm/kernels/boundary/inlet/common.cuh +++ b/src/CelerisLab/lbm/kernels/boundary/inlet/common.cuh @@ -46,11 +46,8 @@ __device__ __forceinline__ float west_velocity_rho_closure_d2q9( { float rho_num = f[0] + f[3] + f[4] + 2.0f * (f[2] + f[6] + f[8]); #if USE_DDF_SHIFTING - // f[i] = physical_f[i] - w[i]; restore the sum of weights for - // the directions appearing in the numerator (w[0]=4/9, w[3]=w[4]=1/9, - // w[2]=w[6]=w[8]=1/36, with f[2],f[6],f[8] double-counted): - // sum(w) = 4/9 + 1/9 + 1/9 + 2*(1/36+1/36+1/36) = 6/9 + 6/36 = 5/6 - rho_num += 5.0f / 6.0f; + // Sum of the descriptor weights in this closure (including doubles). + rho_num += 1.0f; #endif return rho_num / (1.0f - ux_target); } @@ -64,11 +61,8 @@ __device__ __forceinline__ float west_velocity_rho_closure_d3q19( float rho_num = f[0] + f[3] + f[4] + f[5] + f[6] + f[11] + f[12] + f[17] + f[18] + 2.0f * (f[2] + f[8] + f[10] + f[14] + f[16]); #if USE_DDF_SHIFTING - // f[i] = physical_f[i] - w[i]; restore the sum of weights for - // the directions appearing in the numerator (w[0]=1/3, w[3..6,11,12,17,18]=1/18, - // w[2,8,10,14,16]=1/36 double-counted): - // sum(w) = 1/3 + 8*(1/18) + 5*(2/36) = 12/36 + 16/36 + 10/36 = 38/36 = 19/18 - rho_num += 19.0f / 18.0f; + // Sum of the descriptor weights in this closure (including doubles). + rho_num += 1.0f; #endif return rho_num / (1.0f - ux_target); } diff --git a/src/CelerisLab/lbm/kernels/boundary/inlet/zou_he_local.cuh b/src/CelerisLab/lbm/kernels/boundary/inlet/zou_he_local.cuh index c2526e5..c8352ce 100644 --- a/src/CelerisLab/lbm/kernels/boundary/inlet/zou_he_local.cuh +++ b/src/CelerisLab/lbm/kernels/boundary/inlet/zou_he_local.cuh @@ -21,34 +21,32 @@ // west_velocity_rho_closure_d2q9() need clean known-direction values. The // free-slip wall interferes with these at the top/bottom inlet corners. // The legacy name is kept for now (not renamed yet). -__device__ inline void repair_zou_he_west_knowns_d2q9( - float* __restrict__ f, - const fpxx* __restrict__ fi_in, - unsigned int x, - unsigned int y, - const unsigned long* j_neb = nullptr, - unsigned long t = 0ul) +__device__ inline void repair_west_knowns_from_donor_d2q9( + float* __restrict__ f, const float* __restrict__ donor, + unsigned int x, unsigned int y) { if (x != 0u) return; - - const unsigned long k_int = linear_index(x + 1u, y); - if (y == 1u) { - if (t == 0ul || j_neb == nullptr) { - f[3] = load_ddf(fi_in, index_f(k_int, 3u)); - f[8] = load_ddf(fi_in, index_f(k_int, 8u)); - } else { - f[3] = load_physical_dir_pre_streaming(k_int, 3u, fi_in, j_neb, t); - f[8] = load_physical_dir_pre_streaming(k_int, 8u, fi_in, j_neb, t); - } + f[3] = donor[3]; f[8] = donor[8]; } else if (y == (unsigned int)(NY - 2)) { - if (t == 0ul || j_neb == nullptr) { - f[4] = load_ddf(fi_in, index_f(k_int, 4u)); - f[6] = load_ddf(fi_in, index_f(k_int, 6u)); - } else { - f[4] = load_physical_dir_pre_streaming(k_int, 4u, fi_in, j_neb, t); - f[6] = load_physical_dir_pre_streaming(k_int, 6u, fi_in, j_neb, t); - } + f[4] = donor[4]; f[6] = donor[6]; + } +} + +// Double-buffer donors are immutable during the pull kernel. EsoPull calls +// the helper above with the same spatial state staged before in-place writes. +__device__ inline void repair_zou_he_west_knowns_d2q9( + float* __restrict__ f, const fpxx* __restrict__ fi_in, + unsigned int x, unsigned int y) +{ + if (x != 0u || (y != 1u && y != (unsigned int)(NY - 2))) return; + const unsigned long k_int = linear_index(x + 1u, y); + if (y == 1u) { + f[3] = load_ddf(fi_in, index_f(k_int, 3u)); + f[8] = load_ddf(fi_in, index_f(k_int, 8u)); + } else { + f[4] = load_ddf(fi_in, index_f(k_int, 4u)); + f[6] = load_ddf(fi_in, index_f(k_int, 6u)); } } diff --git a/src/CelerisLab/lbm/kernels/boundary/inlet_outlet.cuh b/src/CelerisLab/lbm/kernels/boundary/inlet_outlet.cuh index cd9b812..0add5ba 100644 --- a/src/CelerisLab/lbm/kernels/boundary/inlet_outlet.cuh +++ b/src/CelerisLab/lbm/kernels/boundary/inlet_outlet.cuh @@ -131,29 +131,19 @@ __device__ __forceinline__ void apply_inlet_esopull_d2q9( float* __restrict__ f, unsigned int x, unsigned int y, - const fpxx* __restrict__ fi, - unsigned long t) + const float* __restrict__ f_neb) { const float u_target = inlet_target_u((float)y); const float v_target = 0.0f; // Free-slip y-walls: repair west inlet ghost-node known directions. #if Y_WALL_BC == 1 - { - unsigned long j_neb[NQ]; - compute_neighbors(x + 1u, y, j_neb); - repair_zou_he_west_knowns_d2q9(f, fi, x, y, j_neb, t); - } + repair_west_knowns_from_donor_d2q9(f, f_neb, x, y); #endif #if INLET_SCHEME == 0 apply_zou_he_left_velocity_inlet_d2q9(f, u_target, v_target); #elif INLET_SCHEME == 1 || INLET_SCHEME == 3 - const unsigned long k_neb = linear_index(x + 1u, y); - unsigned long j_neb[NQ]; - compute_neighbors(x + 1u, y, j_neb); - float f_neb[NQ]; - load_f_esopull(k_neb, f_neb, fi, j_neb, t); #if INLET_SCHEME == 1 apply_channel_stabilized_inlet_d2q9(f, f_neb, (float)y); #else @@ -167,17 +157,8 @@ __device__ __forceinline__ void apply_inlet_esopull_d2q9( } __device__ __forceinline__ void apply_outlet_esopull_d2q9( - float* __restrict__ f, - unsigned int x, - unsigned int y, - const fpxx* __restrict__ fi, - unsigned long t) + float* __restrict__ f, const float* __restrict__ f_neb) { - const unsigned long k_neb = linear_index(x - 1u, y); - unsigned long j_neb[NQ]; - compute_neighbors(x - 1u, y, j_neb); - float f_neb[NQ]; - load_f_esopull(k_neb, f_neb, fi, j_neb, t); apply_pressure_outlet_d2q9(f, f_neb); } #endif diff --git a/src/CelerisLab/lbm/kernels/boundary/planar_velocity.cuh b/src/CelerisLab/lbm/kernels/boundary/planar_velocity.cuh new file mode 100644 index 0000000..87db746 --- /dev/null +++ b/src/CelerisLab/lbm/kernels/boundary/planar_velocity.cuh @@ -0,0 +1,44 @@ +// On-site D2Q9 Zou-He velocity closure, paired lattice ordering. +// Input is poststream f at the arrival time. Only unknown populations change; +// the caller owns collision/store. No force, shifting or off-lattice walls. +#ifndef CELERIS_BOUNDARY_PLANAR_VELOCITY_CUH +#define CELERIS_BOUNDARY_PLANAR_VELOCITY_CUH + +__device__ __forceinline__ float bottom_wall_normal_velocity(unsigned int x) +{ +#if WALL_SLOT_WIDTH > 0 && WALL_SLOT_WIDTH < NX + // C1 compact profile centered at (NX-1)/2; amplitude is its continuous peak. + const float s = (2.0f*(float)x-(float)(NX-1))/(float)WALL_SLOT_WIDTH; + const float w = fmaxf(0.0f, 1.0f-s*s); + return WALL_NORMAL_VELOCITY*w*w; +#else + return WALL_NORMAL_VELOCITY; // Disabled or explicitly uniform full width. +#endif +} + +__device__ __forceinline__ void reconstruct_bottom_velocity(float* f, float u, float v) +{ + const float rho = (f[0] + f[1] + f[2] + 2.0f*(f[4] + f[6] + f[7])) / (1.0f-v); + f[3] = f[4] + (2.0f/3.0f)*rho*v; + f[5] = f[6] + 0.5f*(f[2]-f[1]) + 0.5f*rho*u + (1.0f/6.0f)*rho*v; + f[8] = f[7] + 0.5f*(f[1]-f[2]) - 0.5f*rho*u + (1.0f/6.0f)*rho*v; +} + +__device__ __forceinline__ void reconstruct_top_velocity(float* f, float u, float v) +{ + const float rho = (f[0] + f[1] + f[2] + 2.0f*(f[3] + f[5] + f[8])) / (1.0f+v); + f[4] = f[3] - (2.0f/3.0f)*rho*v; + f[6] = f[5] + 0.5f*(f[1]-f[2]) - 0.5f*rho*u - (1.0f/6.0f)*rho*v; + f[7] = f[8] + 0.5f*(f[2]-f[1]) + 0.5f*rho*u - (1.0f/6.0f)*rho*v; +} + +// Fixed top density and tangential u=0; positive Cartesian v exits upward. +// J_y=rho*v follows from known poststream populations (ordinary D2Q9). +__device__ __forceinline__ void reconstruct_top_reservoir(float* f, float rho) +{ + const float jy = f[0]+f[1]+f[2]+2.0f*(f[3]+f[5]+f[8])-rho; + f[4] = f[3] - (2.0f/3.0f)*jy; + f[6] = f[5] + 0.5f*(f[1]-f[2]) - (1.0f/6.0f)*jy; + f[7] = f[8] + 0.5f*(f[2]-f[1]) - (1.0f/6.0f)*jy; +} +#endif diff --git a/src/CelerisLab/lbm/kernels/kernel_v2.cu b/src/CelerisLab/lbm/kernels/kernel_v2.cu index 2e4b7bc..5637974 100644 --- a/src/CelerisLab/lbm/kernels/kernel_v2.cu +++ b/src/CelerisLab/lbm/kernels/kernel_v2.cu @@ -56,6 +56,9 @@ #include "boundary/bounce_back.cuh" #include "boundary/curved_boundary.cuh" #include "boundary/inlet_outlet.cuh" +#if Y_WALL_BC == 2 +#include "boundary/planar_velocity.cuh" +#endif #include "boundary/ibm_kernels.cuh" // ============================================================================ diff --git a/src/CelerisLab/lbm/kernels/operators/collision_mrt.cuh b/src/CelerisLab/lbm/kernels/operators/collision_mrt.cuh index bc934f1..efc59c1 100644 --- a/src/CelerisLab/lbm/kernels/operators/collision_mrt.cuh +++ b/src/CelerisLab/lbm/kernels/operators/collision_mrt.cuh @@ -16,21 +16,13 @@ // m[7] = pxx (stress, s₇ = s_nu = ω = 1/(3ν + 0.5)) // m[8] = pxy (stress, s₈ = s_nu) // -// NOTE on ddf_shifting (USE_DDF_SHIFTING): -// MRT is NOT shift-invariant in moment space because the moment transform -// of the weight vector w produces non-zero entries in m[1] (energy) and -// m[2] (energy^2). Unlike SRT/TRT where (f-feq) cancels weights pairwise, -// the MRT moment construction and inverse transform do not preserve the -// shifted-space semantics. +// D2Q9 transforms the stored populations directly. For shifted storage, +// M*w = (1,-2,1,0,0,0,0,0,0) is removed from the equilibrium moments. +// This avoids rounding small populations through an O(1) weight background. +// D3Q19 projects g-feq; both arrays use the same storage representation. // -// Therefore when USE_DDF_SHIFTING=1, this function temporarily unshifts -// the distributions to physical space before collision, then shifts them -// back after. This ensures the moment construction, relaxation, and -// inverse transform all operate on physically meaningful distributions. -// -// Forcing contract: -// Fin[i] is the raw Guo source term without the (1 - ω/2) prefactor. -// This operator applies the prefactor internally, matching collide_srt(). +// Fin is the raw Guo source. Applying half before and half after linear +// relaxation gives the model-specific M^-1 (I-S/2) M source correction. // ============================================================================ #ifndef CELERIS_OPERATORS_COLLISION_MRT_CUH @@ -43,14 +35,6 @@ __device__ __forceinline__ void collide_mrt(float* __restrict__ g, const float* __restrict__ Fin, float omega) { - // MRT is not shift-invariant in moment space. When USE_DDF_SHIFTING=1, - // temporarily unshift to physical space before collision. -#if USE_DDF_SHIFTING - g[0] += (float)d_w[0]; g[1] += (float)d_w[1]; g[2] += (float)d_w[2]; - g[3] += (float)d_w[3]; g[4] += (float)d_w[4]; g[5] += (float)d_w[5]; - g[6] += (float)d_w[6]; g[7] += (float)d_w[7]; g[8] += (float)d_w[8]; -#endif - const float s_rho = 0.0f; const float s_e = 1.2f; const float s_eps = 1.2f; @@ -58,9 +42,10 @@ __device__ __forceinline__ void collide_mrt(float* __restrict__ g, const float s_q = 1.2f; const float s_jy = 0.0f; const float s_nu = omega; - const float c_tau = 1.0f - 0.5f * omega; - const float p = rho * (1.0f / 3.0f); + // Trapezoidal source: collide (g + Fin/2), then add Fin/2. + // Equivalent to M^-1 (I-S/2) M Fin for this relaxation matrix. + for (int i = 0; i < 9; ++i) g[i] += 0.5f * Fin[i]; float m[9]; m[0] = g[0] + g[1] + g[2] + g[3] + g[4] + g[5] + g[6] + g[7] + g[8]; @@ -75,9 +60,17 @@ __device__ __forceinline__ void collide_mrt(float* __restrict__ g, const float u2 = ux * ux + uy * uy; float meq[9]; +#if USE_DDF_SHIFTING + const float drho = rho - 1.0f; + meq[0] = drho; + meq[1] = -2.0f * drho + 3.0f * rho * u2; + meq[2] = drho - 3.0f * rho * u2; +#else + const float p = rho * (1.0f / 3.0f); meq[0] = 3.0f * p; meq[1] = -6.0f * p + 3.0f * rho * u2; meq[2] = 3.0f * p - 3.0f * rho * u2; +#endif meq[3] = rho * ux; meq[4] = -rho * ux; meq[5] = rho * uy; @@ -108,14 +101,9 @@ __device__ __forceinline__ void collide_mrt(float* __restrict__ g, #pragma unroll for (int i = 0; i < 9; i++) { - g[i] += c_tau * Fin[i]; + g[i] += 0.5f * Fin[i]; } -#if USE_DDF_SHIFTING - g[0] -= (float)d_w[0]; g[1] -= (float)d_w[1]; g[2] -= (float)d_w[2]; - g[3] -= (float)d_w[3]; g[4] -= (float)d_w[4]; g[5] -= (float)d_w[5]; - g[6] -= (float)d_w[6]; g[7] -= (float)d_w[7]; g[8] -= (float)d_w[8]; -#endif } __device__ __forceinline__ void collide_mrt_no_force(float* __restrict__ g, @@ -133,11 +121,7 @@ __device__ __forceinline__ void collide_mrt(float* __restrict__ g, const float* __restrict__ Fin, float omega) { - // MRT is not shift-invariant. Temporarily unshift to physical space. -#if USE_DDF_SHIFTING - #pragma unroll - for (int i = 0; i < 19; i++) g[i] += (float)d_w[i]; -#endif + // g and compute_feq use the same shifted/unshifted representation. float feq[19]; compute_feq(rho, ux, uy, uz, feq); @@ -145,7 +129,7 @@ __device__ __forceinline__ void collide_mrt(float* __restrict__ g, const float s_nu = omega; const float s_bulk = 1.2f; const float s_high = 1.4f; - const float c_tau = 1.0f - 0.5f * omega; + const float one_minus_s_nu = 1.0f - s_nu; const float one_minus_s_bulk = 1.0f - s_bulk; @@ -154,7 +138,8 @@ __device__ __forceinline__ void collide_mrt(float* __restrict__ g, float neq[19]; #pragma unroll for (int i = 0; i < 19; i++) { - neq[i] = g[i] - feq[i]; + // Half-source before relaxation; the other half is added below. + neq[i] = g[i] - feq[i] + 0.5f * Fin[i]; } float pixx = 0.0f, piyy = 0.0f, pizz = 0.0f; @@ -214,13 +199,9 @@ __device__ __forceinline__ void collide_mrt(float* __restrict__ g, + one_minus_s_nu * neq_dev + one_minus_s_bulk * neq_iso + one_minus_s_high * neq_h - + c_tau * Fin[i]; + + 0.5f * Fin[i]; } -#if USE_DDF_SHIFTING - #pragma unroll - for (int i = 0; i < 19; i++) g[i] -= (float)d_w[i]; -#endif } #endif // NQ diff --git a/src/CelerisLab/lbm/kernels/operators/collision_trt.cuh b/src/CelerisLab/lbm/kernels/operators/collision_trt.cuh index 1538652..ed28496 100644 --- a/src/CelerisLab/lbm/kernels/operators/collision_trt.cuh +++ b/src/CelerisLab/lbm/kernels/operators/collision_trt.cuh @@ -11,7 +11,7 @@ // // Forcing contract: // Fin[i] is the raw Guo source term without the (1 - ω/2) prefactor. -// This operator applies the prefactor internally, matching collide_srt(). +// Even/odd source components use their own (1 - omega_plus/minus / 2). // ============================================================================ #ifndef CELERIS_OPERATORS_COLLISION_TRT_CUH @@ -33,6 +33,7 @@ __device__ __forceinline__ void collide_trt(float* __restrict__ f, const float wp = omega; const float wm = compute_omega_minus(wp); const float c_tau = 1.0f - 0.5f * wp; + const float c_odd = 1.0f - 0.5f * wm; f[0] = f[0] - wp * (f[0] - feq[0]) + c_tau * Fin[0]; @@ -49,8 +50,10 @@ __device__ __forceinline__ void collide_trt(float* __restrict__ f, const float delta_s = 0.5f * ((fi - feqi) + (fb - feqb)); const float delta_a = 0.5f * ((fi - feqi) - (fb - feqb)); - f[i] = fi - wp * delta_s - wm * delta_a + c_tau * Fin[i]; - f[ib] = fb - wp * delta_s + wm * delta_a + c_tau * Fin[ib]; + const float force_s = c_tau * 0.5f * (Fin[i] + Fin[ib]); + const float force_a = c_odd * 0.5f * (Fin[i] - Fin[ib]); + f[i] = fi - wp * delta_s - wm * delta_a + force_s + force_a; + f[ib] = fb - wp * delta_s + wm * delta_a + force_s - force_a; } } diff --git a/src/CelerisLab/lbm/kernels/operators/forcing_guo.cuh b/src/CelerisLab/lbm/kernels/operators/forcing_guo.cuh index 95d84c8..fe211d1 100644 --- a/src/CelerisLab/lbm/kernels/operators/forcing_guo.cuh +++ b/src/CelerisLab/lbm/kernels/operators/forcing_guo.cuh @@ -3,7 +3,8 @@ // // Velocity correction: u* = u + F·dt / (2ρ) // Forcing term: -// Fin[i] = (1 - ω/2) · w_i · [ (c_i - u*)/cs² + (c_i · u*)/cs⁴ · c_i ] · F +// Fin[i] = w_i · [ (c_i - u*)/cs² + (c_i · u*)/cs⁴ · c_i ] · F +// The collision operator applies its relaxation-dependent prefactor. // simplified: Fin[i] = 9·w_i · [ (c_i·F)(c_i·u* + 1/3) - u*·F/3 ] // ============================================================================ diff --git a/src/CelerisLab/lbm/kernels/refinement.cu b/src/CelerisLab/lbm/kernels/refinement.cu new file mode 100644 index 0000000..ab4bf20 --- /dev/null +++ b/src/CelerisLab/lbm/kernels/refinement.cu @@ -0,0 +1,120 @@ +// Fixed, unforced internal 2:1 patch. Separate compilation entry; +// default OneStep/channel kernels are deliberately not included. +#include +#include +#include "config.h" +#include "core/precision.cuh" +#include "core/descriptors.cuh" +#include "core/layout.cuh" +#include "core/params.cuh" +#include "core/flags.cuh" +#include "core/obs.cuh" +#include "operators/macro.cuh" +#include "operators/collision_srt.cuh" +#include "streaming/pull_double_buffer.cuh" +#include "boundary/planar_velocity.cuh" +#include "boundary/curved_boundary.cuh" + +extern "C" { +#include "step/curved_double.cu" +} + +#if DIM != 2 || NQ != 9 || COLLISION_MODEL != 0 || STORE_PRECISION != 0 || USE_DDF_SHIFTING || USE_LES || STREAMING_MODEL != 0 +#error "Refinement requires unshifted FP32 D2Q9 SRT double buffers without LES" +#endif + +// Enabled only in the validation build; no diagnostic arrays in normal builds. +#ifdef CELERIS_REFINEMENT_AUDIT +extern "C" { +__device__ unsigned int collision_counts[TOTAL_CELLS]; +__device__ unsigned int explosion_counts[TOTAL_CELLS]; +__device__ unsigned int incoming_counts[NQ * TOTAL_CELLS]; +} +#endif + +__device__ __forceinline__ void patch_collide(float* f, unsigned long k) { + float rho, ux, uy, feq[NQ]; + compute_rho_u(f, rho, ux, uy); + compute_feq(rho, ux, uy, feq); + collide_srt_no_force(f, feq, d_params.omega); +#ifdef CELERIS_REFINEMENT_AUDIT + ++collision_counts[k]; +#endif +} + +extern "C" __global__ void PatchExplode( + const float* coarse, float* fine, const uint32_t* targets, + const uint32_t* donors, unsigned int count, unsigned long coarse_cells) +{ + unsigned int row = blockIdx.x * blockDim.x + threadIdx.x; + if (row >= count) return; + for (unsigned int q = 0; q < NQ; ++q) + fine[index_f(targets[row], q)] = coarse[q * coarse_cells + donors[row]]; +#ifdef CELERIS_REFINEMENT_AUDIT + ++explosion_counts[targets[row]]; +#endif +} + +extern "C" __global__ void PatchFineStep(const float* in, float* out, const uint16_t* flags) { + unsigned int x, y; unsigned long k; + index_from_thread(x, y, k); + if (x >= NX || y >= NY) return; + float f[NQ]; + // Outer halo is invalid after EACH hop, including when reusing the buffer + // that held explosion data. Inner halo only streams; it never collides. + if (x < 2 || y < 2 || x >= NX - 2 || y >= NY - 2 + || (flags != nullptr && is_solid(flags[k]))) { + for (int q = 0; q < NQ; ++q) f[q] = NAN; + } else { + unsigned long j[NQ]; + compute_neighbors(x, y, j); // These interior targets never wrap. + stream_pull_load(k, f, in, j); + if (x >= 4 && y >= 4 && x < NX - 4 && y < NY - 4) + patch_collide(f, k); + } + stream_pull_store(k, f, out); +} + +extern "C" __global__ void PatchCoarseStep( + const float* in, float* out, const float* fine, + const uint64_t* targets, const uint64_t* children, unsigned int count, + unsigned int x0, unsigned int y0, unsigned int x1, unsigned int y1, + unsigned int wall_aware) +{ + unsigned int x, y; unsigned long k; + index_from_thread(x, y, k); + if (x >= NX || y >= NY) return; + // Both buffers were poisoned at initialization; covered cells stay invalid. + if (x >= x0 && x < x1 && y >= y0 && y < y1) return; + if (wall_aware && (y == 0 || y == NY - 1)) return; + unsigned long j[NQ]; + compute_neighbors(x, y, j); + float f[NQ]; + stream_pull_load(k, f, in, j); + if (x >= x0 - 1 && x <= x1 && y >= y0 - 1 && y <= y1) { + for (unsigned int q = 1; q < NQ; ++q) { + uint64_t slot = index_f(k, q); + unsigned int lo = 0, hi = count; + while (lo < hi) { + unsigned int mid = lo + (hi - lo) / 2; + if (targets[mid] < slot) lo = mid + 1; + else hi = mid; + } + if (lo < count && targets[lo] == slot) { + const uint64_t* s = children + 4ul * lo; + f[q] = 0.25f * (fine[s[0]] + fine[s[1]] + fine[s[2]] + fine[s[3]]); +#ifdef CELERIS_REFINEMENT_AUDIT + ++incoming_counts[slot]; +#endif + } + } + } + // Includes the eight grazing-corner links with UNCOVERED coarse donors. + // All replacement arrivals must be in place before computing rho/u. + if (wall_aware) { + if (y == 1) reconstruct_bottom_velocity(f, 0.0f, bottom_wall_normal_velocity(x)); + else if (y == NY - 2) reconstruct_top_reservoir(f, RHO); + } + patch_collide(f, k); + stream_pull_store(k, f, out); +} diff --git a/src/CelerisLab/lbm/kernels/step/aux_kernels.cu b/src/CelerisLab/lbm/kernels/step/aux_kernels.cu index 34f1266..2b985a3 100644 --- a/src/CelerisLab/lbm/kernels/step/aux_kernels.cu +++ b/src/CelerisLab/lbm/kernels/step/aux_kernels.cu @@ -15,15 +15,7 @@ // (pre-streaming semantics = matching double-buffer fi_in[k, dir]). // SensorKernel / ForceRegionKernel run after EsoPullStep, using // t = step_count+1 for post-streaming physical-value semantics. -// Wall boundary functions (bounce_back.cuh) use load_physical_dir_esopull -// inside EsoPullStep (post-streaming, matching load_f_esopull output). -// -// NOTE: The wall boundary helpers use post-streaming reads (unlike Bouzidi) -// because they operate inside EsoPullStep on f[] already reconstructed -// by load_f_esopull. Replacing post-streaming f[i] with post-streaming -// f[opp(i)] is the correct half-way bounce-back for EsoPull. -// -// Verified: Kan99b K2 (60k steps, MRT) — bit-identical to double-buffer. +// Y-wall donors are staged by PrepareEsoPullBoundaries before any backing writes. // ============================================================================ #ifndef CELERIS_STEP_AUX_KERNELS_CU @@ -32,54 +24,7 @@ #include "../streaming/esopull_semantic_helpers.cuh" #include "../boundary/esopull_curved_helpers.cuh" -__device__ __forceinline__ float action_omega(const float* action, int body_id) -{ - if (action == nullptr || body_id < 0) return 0.0f; - unsigned int base = (unsigned int)body_id * (unsigned int)OBS_BODY_SLOT_FLOATS; - return action[base + (unsigned int)(OBS_BODY_SLOT_FLOATS - 1)]; -} - -__global__ void CurvedBoundaryKernel( - fpxx* fi, - const unsigned int* cl_fluid_idx, - const unsigned char* cl_dir, - const float* cl_q, - const float* cl_rx, - const float* cl_ry, - const float* cl_rz, - const unsigned char* cl_fallback_class, - const int* cl_body_id, - const float* action, - float* obs, - unsigned int n_curved) -{ - unsigned int tid = threadIdx.x + blockIdx.x * blockDim.x; - if (tid >= n_curved) return; - - const unsigned long k_f = (unsigned long)cl_fluid_idx[tid]; - const unsigned int dir = (unsigned int)cl_dir[tid]; - const float q = cl_q[tid]; - const int bid = cl_body_id[tid]; - const float rx = cl_rx[tid]; - const float ry = cl_ry[tid]; - const unsigned char fallback_class = cl_fallback_class[tid]; - const float omega = action_omega(action, bid); - -#if DIM == 2 - const float Uw = -omega * ry; - const float Vw = omega * rx; - (void)cl_rz; - apply_bouzidi_link(dir, q, k_f, fi, Uw, Vw, rx, ry, fallback_class, obs, bid); -#elif DIM == 3 - const float rz = cl_rz[tid]; - const float Uw = -omega * ry; - const float Vw = omega * rx; - const float Ww = 0.0f; - apply_bouzidi_link(dir, q, k_f, fi, Uw, Vw, Ww, - rx, ry, rz, fallback_class, obs, bid); -#endif -} -// End of CurvedBoundaryKernel +#include "curved_double.cu" // --------------------------------------------------------------------------- // PrepareEsoPullCurvedKernel -- read-only curved prep for EsoPull @@ -146,7 +91,12 @@ __global__ void PrepareEsoPullCurvedKernel( } else { f_reflected = compute_bouzidi_reflection(f_toward, 0.0f, 0.0f, q, fallback); } - f_reflected += compute_moving_wall_addend(dir, q, fallback, Uw, Vw); + float rho_wall = USE_DDF_SHIFTING ? 1.0f : 0.0f; + if (Uw != 0.0f || Vw != 0.0f) { + for (unsigned int i = 0; i < NQ; ++i) + rho_wall += load_physical_dir_pre_streaming(k_f, i, fi, j_f, t); + } + f_reflected += compute_moving_wall_addend(dir, q, fallback, Uw, Vw, rho_wall); unsigned long dst_cell; unsigned int dst_slot; @@ -249,14 +199,13 @@ __global__ void SensorKernel( } // --------------------------------------------------------------------------- -// ForceRegionKernel — operator-split local Guo forcing correction +// ForceRegionKernel — operator-split exact-difference force-density impulse // -// Runs after the main collide+stream step. Applies a Guo forcing term to -// cells within a force_region footprint. This is an operator-split correction, -// not inline in OneStep's collide_dispatch. -// -// The Guo formula and c_tau = 1 - omega/2 conventions match collide_srt, -// collide_mrt, and collide_trt exactly. +// Runs after the main step. The raw source evaluated at u + F/(2 rho) +// equals feq(rho, u + F/rho) - feq(rho, u) for the quadratic equilibrium. +// No collision prefactor: collision has already happened. This preserves +// density and adds exactly F momentum, independently of viscosity. It is +// a split impulse, not the collision-integrated global Guo discretization. // --------------------------------------------------------------------------- __device__ __forceinline__ float action_force_x(const float* action, int body_id) { @@ -319,11 +268,9 @@ __global__ void ForceRegionKernel( float Fin[NQ]; compute_guo_forcing(ux, uy, fx, fy, Fin); - float omega = d_params.omega; - float c_tau = 1.0f - 0.5f * omega; #pragma unroll for (int i = 0; i < NQ; i++) - f[i] += c_tau * Fin[i]; + f[i] += Fin[i]; #elif DIM == 3 float rho_n, ux, uy, uz; compute_rho_u(f, rho_n, ux, uy, uz); @@ -334,11 +281,9 @@ __global__ void ForceRegionKernel( float Fin[NQ]; compute_guo_forcing(ux, uy, uz, fx, fy, 0.0f, Fin); - float omega = d_params.omega; - float c_tau = 1.0f - 0.5f * omega; #pragma unroll for (int i = 0; i < NQ; i++) - f[i] += c_tau * Fin[i]; + f[i] += Fin[i]; #endif // Write back with correct semantics diff --git a/src/CelerisLab/lbm/kernels/step/curved_double.cu b/src/CelerisLab/lbm/kernels/step/curved_double.cu new file mode 100644 index 0000000..1b16f87 --- /dev/null +++ b/src/CelerisLab/lbm/kernels/step/curved_double.cu @@ -0,0 +1,61 @@ +// Shared compact double-buffer curved boundary; caller supplies core/BFL headers. +#ifndef CELERIS_CURVED_DOUBLE_CU +#define CELERIS_CURVED_DOUBLE_CU + +__device__ __forceinline__ float action_omega(const float* action, int body_id) +{ + if (action == nullptr || body_id < 0) return 0.0f; + unsigned int base = (unsigned int)body_id * (unsigned int)OBS_BODY_SLOT_FLOATS; + return action[base + (unsigned int)(OBS_BODY_SLOT_FLOATS - 1)]; +} + +__global__ void CurvedBoundaryKernel( + fpxx* fi, + const unsigned int* cl_fluid_idx, + const unsigned char* cl_dir, + const float* cl_q, + const float* cl_rx, + const float* cl_ry, + const float* cl_rz, + const unsigned char* cl_fallback_class, + const int* cl_body_id, + const float* action, + float* obs, + unsigned int n_curved, + const float* surface_velocity = nullptr, float* link_exchange = nullptr) +{ + unsigned int tid = threadIdx.x + blockIdx.x * blockDim.x; + if (tid >= n_curved) return; + + const unsigned long k_f = (unsigned long)cl_fluid_idx[tid]; + const unsigned int dir = (unsigned int)cl_dir[tid]; + const float q = cl_q[tid]; + const int bid = cl_body_id[tid]; + const float rx = cl_rx[tid]; + const float ry = cl_ry[tid]; + const unsigned char fallback_class = cl_fallback_class[tid]; + const float omega = action_omega(action, bid); + +#if DIM == 2 + float Uw = -omega * ry; + float Vw = omega * rx; + if (surface_velocity != nullptr) { + // Independent prescribed fluid velocity on fixed geometry. + const float u = surface_velocity[2*tid], v = surface_velocity[2*tid+1]; + if (u != 0.0f || v != 0.0f) { Uw += u; Vw += v; } + } + (void)cl_rz; + apply_bouzidi_link(dir, q, k_f, fi, Uw, Vw, rx, ry, fallback_class, obs, bid, + link_exchange == nullptr ? nullptr : link_exchange + 3*tid); +#elif DIM == 3 + const float rz = cl_rz[tid]; + const float Uw = -omega * ry; + const float Vw = omega * rx; + const float Ww = 0.0f; + apply_bouzidi_link(dir, q, k_f, fi, Uw, Vw, Ww, + rx, ry, rz, fallback_class, obs, bid); +#endif +} +// End of CurvedBoundaryKernel + +#endif diff --git a/src/CelerisLab/lbm/kernels/step/esopull_sync.cu b/src/CelerisLab/lbm/kernels/step/esopull_sync.cu index 54a7dc8..30ce615 100644 --- a/src/CelerisLab/lbm/kernels/step/esopull_sync.cu +++ b/src/CelerisLab/lbm/kernels/step/esopull_sync.cu @@ -47,10 +47,8 @@ __global__ void DecodeCellsToPhysical( // --------------------------------------------------------------------------- // EncodePhysicalToCells: physical DDF -> backing layout (compact list) // -// Reads physical (pre-collision) DDF from *fi_phys*, applies collision -// to produce post-collision DDF, then encodes into the EsoPull backing -// layout on *fi*. The backing layout stores post-collision values, -// matching what store_f_esopull would produce at step *t*. +// Inverse of DecodeCellsToPhysical at the SAME completed step count. +// Conversion only: no collision or time advancement, including unchanged halo. // --------------------------------------------------------------------------- __global__ void EncodePhysicalToCells( fpxx* fi, @@ -70,10 +68,9 @@ __global__ void EncodePhysicalToCells( for (int i = 0; i < NQ; i++) f[i] = fi_phys[tid * NQ + i]; - // Collide: physical (pre-collision) -> post-collision - float rho_n, ux, uy; - compute_rho_u(f, rho_n, ux, uy); - collide_dispatch(f, rho_n, ux, uy); +#if USE_DDF_SHIFTING + for (int i = 0; i < NQ; ++i) f[i] -= (float)d_w[i]; +#endif store_physical_node_esopull(k, f, fi, j, t); } diff --git a/src/CelerisLab/lbm/kernels/step/init_flow.cu b/src/CelerisLab/lbm/kernels/step/init_flow.cu index a576883..a31072e 100644 --- a/src/CelerisLab/lbm/kernels/step/init_flow.cu +++ b/src/CelerisLab/lbm/kernels/step/init_flow.cu @@ -6,98 +6,41 @@ #ifndef CELERIS_STEP_INIT_FLOW_CU #define CELERIS_STEP_INIT_FLOW_CU -// ---- Helper: write equilibrium for a given velocity ---- +// Reuse the canonical equilibrium, including its cancellation-free shifted +// expression. Recomputing physical feq then subtracting w erases tiny speeds. __device__ __forceinline__ void write_equilibrium( fpxx* fi, unsigned long k, float u_init) { - for (int i = 0; i < NQ; i++) { - float cu = (float)d_cx[i] * u_init; - float val = d_w[i] * RHO * (1.0f + 3.0f*cu + 4.5f*cu*cu - - 1.5f*u_init*u_init); -#if USE_DDF_SHIFTING - val -= d_w[i]; + float f[NQ]; +#if DIM == 2 + compute_feq(RHO, u_init, 0.0f, f); +#else + compute_feq(RHO, u_init, 0.0f, 0.0f, f); #endif - store_ddf(fi, index_f(k, (unsigned int)i), val); - } + for (int i = 0; i < NQ; ++i) + store_ddf(fi, index_f(k, (unsigned int)i), f[i]); } -__device__ __forceinline__ void write_rest_equilibrium( - fpxx* fi, unsigned long k) +__device__ __forceinline__ void write_rest_equilibrium(fpxx* fi, unsigned long k) { - for (int i = 0; i < NQ; i++) { - float val = d_w[i] * RHO; -#if USE_DDF_SHIFTING - val -= d_w[i]; -#endif - store_ddf(fi, index_f(k, (unsigned int)i), val); - } + write_equilibrium(fi, k, 0.0f); } -// ---- EsoPull layout helpers: equilibrium in EsoPull backing layout at t=0 ---- -// In EsoPull t=0 (even step), store_f_esopull writes: -// f[i_odd] → fi[ j[i_odd], i_odd ] (neighbor, forward slot) -// f[i_even] → fi[ k, i_even ] (local, reverse slot) -// These helpers write equilibrium values directly in the correct EsoPull -// layout so that load_f_esopull(t=0) reconstructs the correct f_eq[]. +// Init is load-compatible at t=0. Same-time encoding is the inverse of that +// read, unlike store_f_esopull(0), which advances the streaming interpretation. #if DIM == 2 __device__ __forceinline__ void write_equilibrium_esopull( fpxx* fi, unsigned long k, unsigned int x, unsigned int y, float u_init) { - unsigned long j[NQ]; - compute_neighbors(x, y, j); - - // --- rest particle: always local, slot 0 --- - float val0 = d_w[0] * RHO * (1.0f - 1.5f * u_init * u_init); -#if USE_DDF_SHIFTING - val0 -= d_w[0]; -#endif - store_ddf(fi, index_f(k, 0u), val0); - - // --- paired directions (i=1,3,5,7): {E,W}, {N,S}, {NE,SW}, {SE,NW} --- - for (int i = 1; i < NQ; i += 2) { - float cu_odd = (float)d_cx[i] * u_init; - float val_odd = d_w[i] * RHO - * (1.0f + 3.0f * cu_odd + 4.5f * cu_odd * cu_odd - - 1.5f * u_init * u_init); - float cu_even = (float)d_cx[i + 1] * u_init; - float val_even = d_w[i + 1] * RHO - * (1.0f + 3.0f * cu_even + 4.5f * cu_even * cu_even - - 1.5f * u_init * u_init); -#if USE_DDF_SHIFTING - val_odd -= d_w[i]; - val_even -= d_w[i + 1]; -#endif - // t=0 even step (load_f_esopull): f[i] reads from local slot i+1, - // f[i+1] reads from neighbor slot i. - // Therefore: fi[n, i+1] = f_eq[i] (odd goes to local reverse), - // fi[j[i], i] = f_eq[i+1] (even goes to neighbor forward) - store_ddf(fi, index_f(k, (unsigned int)(i + 1)), val_odd); - store_ddf(fi, index_f(j[i], (unsigned int)i), val_even); - } + unsigned long j[NQ]; compute_neighbors(x, y, j); + float f[NQ]; compute_feq(RHO, u_init, 0.0f, f); + store_physical_node_esopull(k, f, fi, j, 0ul); } __device__ __forceinline__ void write_rest_equilibrium_esopull( fpxx* fi, unsigned long k, unsigned int x, unsigned int y) { - unsigned long j[NQ]; - compute_neighbors(x, y, j); - - float val0 = d_w[0] * RHO; -#if USE_DDF_SHIFTING - val0 -= d_w[0]; -#endif - store_ddf(fi, index_f(k, 0u), val0); - - for (int i = 1; i < NQ; i += 2) { - float val_odd = d_w[i] * RHO; - float val_even = d_w[i + 1] * RHO; -#if USE_DDF_SHIFTING - val_odd -= d_w[i]; - val_even -= d_w[i + 1]; -#endif - store_ddf(fi, index_f(k, (unsigned int)(i + 1)), val_odd); - store_ddf(fi, index_f(j[i], (unsigned int)i), val_even); - } + write_equilibrium_esopull(fi, k, x, y, 0.0f); } #endif // DIM == 2 @@ -106,10 +49,12 @@ __device__ __forceinline__ uint16_t channel_flag_from_coords( { if (y == 0u || y == (unsigned int)(NY - 1)) return (uint16_t)(FLAG_SOLID | FLAG_BC_WALL); +#if Y_WALL_BC != 2 if (x == 0u) return (uint16_t)(FLAG_SOLID | FLAG_BC_INLET); if (x == (unsigned int)(NX - 1)) return (uint16_t)(FLAG_SOLID | FLAG_BC_OUTLET); +#endif return FLAG_FLUID; } @@ -142,7 +87,21 @@ __global__ void InitTubeFlow_v2(uint16_t* flag, fpxx* fi) if (is_solid(fl)) { write_rest_equilibrium(fi, k); } else { +#if Y_WALL_BC == 2 + // Initial g(0)=feq: both boundary velocities are already satisfied. + // Start from linear Couette, not the sought exponential steady profile. + float f[NQ]; +#if WALL_SLOT_WIDTH > 0 && WALL_SLOT_WIDTH < NX + const float t = (float)(y-1u)/(float)(NY-3); + const float v = (1.0f-t)*bottom_wall_normal_velocity(x)+t*TOP_WALL_NORMAL_VELOCITY; + compute_feq(RHO, 0.0f, v, f); +#else + compute_feq(RHO, U0*(float)(y-1u)/(float)(NY-3), WALL_NORMAL_VELOCITY, f); +#endif + stream_pull_store(k, f, fi); +#else write_equilibrium(fi, k, inlet_target_u((float)y)); +#endif } #elif DIM == 3 unsigned int x, y, z; unsigned long k; @@ -195,4 +154,4 @@ __global__ void InitEsoPull(uint16_t* flag, fpxx* fi) #endif } -#endif // CELERIS_STEP_INIT_FLOW_CU \ No newline at end of file +#endif // CELERIS_STEP_INIT_FLOW_CU diff --git a/src/CelerisLab/lbm/kernels/step/one_step_double.cu b/src/CelerisLab/lbm/kernels/step/one_step_double.cu index 3188db6..ef582b9 100644 --- a/src/CelerisLab/lbm/kernels/step/one_step_double.cu +++ b/src/CelerisLab/lbm/kernels/step/one_step_double.cu @@ -88,7 +88,16 @@ void OneStep( bounce_back_swap(f); if (is_fluid(fl) && (y == 1u || y == (unsigned int)(NY - 2))) { -#if Y_WALL_BC == 1 +#if Y_WALL_BC == 2 + // Reconstruct f(t+dt) on the physical plane, then collide below. + // This replaces, rather than compounds, the halfway wall correction. + if (y == 1u) reconstruct_bottom_velocity(f, 0.0f, bottom_wall_normal_velocity(x)); +#if TOP_WALL_BC == 1 + else reconstruct_top_reservoir(f, RHO); +#else + else reconstruct_top_velocity(f, U0, TOP_WALL_NORMAL_VELOCITY); +#endif +#elif Y_WALL_BC == 1 apply_wall_freeslip_d2q9_y_pull(y, f, fi_in, k); #else apply_wall_bb_d2q9_y_pull(y, f, fi_in, k); diff --git a/src/CelerisLab/lbm/kernels/step/one_step_esopull.cu b/src/CelerisLab/lbm/kernels/step/one_step_esopull.cu index 807504e..b853382 100644 --- a/src/CelerisLab/lbm/kernels/step/one_step_esopull.cu +++ b/src/CelerisLab/lbm/kernels/step/one_step_esopull.cu @@ -13,7 +13,7 @@ #if DIM == 2 __device__ __forceinline__ void apply_boundary_esopull( float* f, uint16_t fl, unsigned int x, unsigned int y, - const fpxx* fi, unsigned long t) + const float* boundary_values) { // Curved-BC nodes have been pre-processed by PrepareEsoPullCurvedKernel // + ApplyEsoPullCurvedKernel, so is_curved() cells here have valid backed @@ -24,10 +24,10 @@ __device__ __forceinline__ void apply_boundary_esopull( bool interior_y = (y > 0u) && (y < (unsigned int)(NY - 1)); if (is_inlet(fl) && interior_y) { - apply_inlet_esopull_d2q9(f, x, y, fi, t); + apply_inlet_esopull_d2q9(f, x, y, boundary_values + 6u * NX + y * NQ); } else if (is_outlet(fl) && interior_y) { - apply_outlet_esopull_d2q9(f, x, y, fi, t); + apply_outlet_esopull_d2q9(f, boundary_values + 6u * NX + (NY + y) * NQ); } else { bounce_back_swap(f); @@ -56,6 +56,43 @@ __device__ __forceinline__ void apply_boundary_esopull_3d( } #endif +// One read-only preparation for all D2 planar boundary donors, before curved +// or in-place writes. [6*NX wall values][2*NY*NQ x-boundary donor values]. +// t=0 uses the initialized physical state; later t uses spatial postcollision +// values, matching the double-buffer donor reads. Slots retain DDF shifting. +__global__ void PrepareEsoPullBoundaries(const fpxx* fi, float* boundary_values, + unsigned long t) +{ +#if DIM == 2 + const unsigned int n = threadIdx.x + blockIdx.x * blockDim.x; + const unsigned int side = blockIdx.y; + if (side >= 2u) return; + if (n < NX) { + const unsigned int y = side == 0u ? 1u : NY - 2u; + const unsigned long k = linear_index(n, y); + unsigned long j[NQ]; compute_neighbors(n, y, j); + const unsigned int lower[3] = {4u, 6u, 7u}; + const unsigned int upper[3] = {3u, 5u, 8u}; + for (unsigned int a = 0; a < 3u; ++a) { + unsigned int donor = side == 0u ? lower[a] : upper[a]; +#if Y_WALL_BC == 1 + if (a > 0u) donor = side == 0u ? lower[3u-a] : upper[3u-a]; +#endif + boundary_values[(side * NX + n) * 3u + a] = + load_physical_dir_pre_streaming(k, donor, fi, j, t); + } + } + if (n < NY) { + const unsigned int x = side == 0u ? 1u : NX - 2u; + const unsigned long k = linear_index(x, n); + unsigned long j[NQ]; compute_neighbors(x, n, j); + float* donor = boundary_values + 6u * NX + (side * NY + n) * NQ; + for (unsigned int i = 0; i < NQ; ++i) + donor[i] = load_physical_dir_pre_streaming(k, i, fi, j, t); + } +#endif +} + // ============================================================================ // EsoPullStep — Esoteric-Pull single-buffer step kernel // ============================================================================ @@ -63,7 +100,7 @@ __global__ __launch_bounds__(NT, 2) void EsoPullStep( fpxx* fi, uint16_t* flag, float* action, float* obs, - unsigned long t) + const float* boundary_values, unsigned long t) { #if DIM == 2 unsigned int x, y; unsigned long k; @@ -87,18 +124,17 @@ void EsoPullStep( // Boundary nodes if (is_solid(fl)) - apply_boundary_esopull(f, fl, x, y, fi, t); + apply_boundary_esopull(f, fl, x, y, boundary_values); if (is_obstacle(fl) && !is_curved(fl)) bounce_back_swap(f); // ---- y-wall adjacent row correction ---- if (is_fluid(fl) && (y == 1u || y == (unsigned int)(NY - 2))) { -#if Y_WALL_BC == 1 - apply_wall_freeslip_d2q9_y_pull(y, f, fi, k, j, t); -#else - apply_wall_bb_d2q9_y_pull(y, f, fi, k, j, t); -#endif + const unsigned int side = y == 1u ? 0u : 1u; + const float* staged = boundary_values + (side * NX + x) * 3u; + if (side == 0u) { f[3] = staged[0]; f[5] = staged[1]; f[8] = staged[2]; } + else { f[4] = staged[0]; f[6] = staged[1]; f[7] = staged[2]; } } // Collision diff --git a/src/CelerisLab/lbm/kernels/streaming/esopull_semantic_helpers.cuh b/src/CelerisLab/lbm/kernels/streaming/esopull_semantic_helpers.cuh index f88c55e..a046e3d 100644 --- a/src/CelerisLab/lbm/kernels/streaming/esopull_semantic_helpers.cuh +++ b/src/CelerisLab/lbm/kernels/streaming/esopull_semantic_helpers.cuh @@ -14,14 +14,14 @@ // ForceRegionKernel(t): reads/writes fi with post-streaming semantics // SensorKernel(t): reads fi with post-streaming semantics // MacroscopicEsoPullKernel(t): reads fi with post-streaming semantics -// Wall boundary functions(t): reads fi with pre-streaming semantics +// PrepareEsoPullBoundaries(t): stages pre-streaming donors before writes // // The "t" parameter always represents the completed step count corresponding // to the current physical interpretation of the backing layout. // Double-buffer mode: t = 0, raw fi[k,i] = physical direction i. // -// Verified: Kan99b K2 (200k steps, MRT) — bit-identical to double-buffer for -// Strouhal, CD, CL, C'L, C'D metrics (FP32, D=20/30). +// The same-time read/write mapping is checked by actual-header regression. +// A physical node here is the next pull state, not the spatial postcollision node. // ============================================================================ #ifndef CELERIS_STREAMING_ESOPULL_SEMANTIC_HELPERS_CUH @@ -130,14 +130,8 @@ __device__ __forceinline__ void store_physical_node_esopull( fpxx* fi, const unsigned long* j, unsigned long t) { store_ddf(fi, index_f(k, 0u), f[0]); - for (unsigned int i = 1; i < NQ; i += 2) { - if (t & 1ul) { - store_ddf(fi, index_f(j[i], i + 1), f[i]); - store_ddf(fi, index_f(k, i), f[i + 1]); - } else { - store_ddf(fi, index_f(j[i], i), f[i]); - store_ddf(fi, index_f(k, i + 1), f[i + 1]); - } + for (unsigned int i = 1; i < NQ; ++i) { + store_ddf(fi, index_f(esopull_src_cell(k, i, j), esopull_slot(i, t)), f[i]); } } diff --git a/src/CelerisLab/lbm/refinement.py b/src/CelerisLab/lbm/refinement.py new file mode 100644 index 0000000..0069026 --- /dev/null +++ b/src/CelerisLab/lbm/refinement.py @@ -0,0 +1,292 @@ +"""CPU preparation for one fixed D2Q9 2:1 patch and optional interior circles. + +Cell indices use y*nx+x; population slots use q*(nx*ny)+cell, matching +core/layout.cuh. Bounds are half-open (x0, y0, x1, y1) in coarse cells. +The four fine halo cells implement the pinned CPUREF volumetric geometry. +Base topology is independent of boundaries; prepare_fine_boundary checks them. +""" +from dataclasses import dataclass, replace +from numbers import Integral + +import numpy as np + +from ..config import LBMConfig +from .descriptors import D2Q9_EX, D2Q9_EY, FLUID, SOLID + +ACTIVE, INNER_GHOST, OUTER_GHOST = 0, 1, 2 +HALO = 4 + + +def _readonly(values, dtype): + array = np.ascontiguousarray(values, dtype=dtype) + array.setflags(write=False) + return array + + +@dataclass(frozen=True, eq=False) +class FixedPatch2D: + """Layer metadata and compact exchange tables, constructed by from_config. + + Shapes are (ny, nx), including the fine halo. ``fine_origin`` is its lower + cell edge in coarse physical units, with coarse dx=dt=1 and fine dx=dt=.5. + ``ghost_fine`` and ``ghost_coarse`` pair fine cells with coarse donors. + ``coalesce_coarse`` contains incoming population slots; each corresponding + row of ``coalesce_fine`` contains four fine slots in (dy, dx) row order. + Weights are all 1/4. Closed corner links are included even when the coarse + donor is uncovered. These tables contain no DDF data and allocate no GPU. + + This is preparation, not permission to advance these levels with OneStep: + its channel boundaries and fused collision do not implement patch coupling. + """ + coarse_shape: tuple[int, int] + bounds: tuple[int, int, int, int] + tau_coarse: float + tau_fine: float + fine_kind: np.ndarray + ghost_fine: np.ndarray + ghost_coarse: np.ndarray + coalesce_coarse: np.ndarray + coalesce_fine: np.ndarray + + @classmethod + def from_config(cls, cfg: LBMConfig, bounds) -> "FixedPatch2D": + """Prepare maps without changing cfg or allocating coarse-sized tables. + + The initial backend contract is D2Q9/SRT/FP32/unshifted/double_buffer, + without LES. Two coarse cells of margin keep all halo donors in bounds; + this does NOT certify that those donors are fluid or periodic. + """ + checked = replace(cfg) + checked.validate() + if (not checked.is_d2q9 or checked.collision != "SRT" + or checked.streaming != "double_buffer" + or checked.store_precision != "FP32" + or checked.ddf_shifting or checked.les_enabled): + raise ValueError("FixedPatch2D requires D2Q9/SRT/FP32/unshifted/double_buffer without LES") + bounds = tuple(bounds) + if len(bounds) != 4 or any(isinstance(v, (bool, np.bool_)) + or not isinstance(v, Integral) for v in bounds): + raise ValueError("bounds must be four integers: (x0, y0, x1, y1)") + x0, y0, x1, y1 = bounds = tuple(map(int, bounds)) + nx, ny = checked.nx, checked.ny + if not (2 <= x0 < x1 <= nx - 2 and 2 <= y0 < y1 <= ny - 2): + raise ValueError("patch must be nonempty with a two-coarse-cell margin") + fx, fy = 2 * (x1 - x0) + 2 * HALO, 2 * (y1 - y0) + 2 * HALO + if max(nx * ny, fx * fy) > np.iinfo(np.uint32).max: + raise ValueError("layer cell counts exceed uint32 indexing") + tau_c = 3.0 * checked.viscosity + 0.5 + tau_f = 6.0 * checked.viscosity + 0.5 + if not checked.omega_min <= 1.0 / tau_f <= checked.omega_max: + raise ValueError("fine relaxation rate lies outside configured omega guards") + + kind = np.full((fy, fx), OUTER_GHOST, dtype=np.uint8) + kind[2:-2, 2:-2] = INNER_GHOST + kind[HALO:-HALO, HALO:-HALO] = ACTIVE + gy, gx = np.nonzero(kind != ACTIVE) + ghost_fine = gy * fx + gx + ghost_coarse = (y0 + (gy - HALO) // 2) * nx + x0 + (gx - HALO) // 2 + + targets, sources = [], [] + # Only the outside ring can receive fine-to-coarse populations. + # Integer doubled midpoints include all eight diagonal corner passes. + for y in range(y0 - 1, y1 + 1): + xs = range(x0 - 1, x1 + 1) if y in (y0 - 1, y1) else (x0 - 1, x1) + for x in xs: + for q in range(1, 9): + mx, my = 2 * x + 1 - D2Q9_EX[q], 2 * y + 1 - D2Q9_EY[q] + if 2 * x0 <= mx <= 2 * x1 and 2 * y0 <= my <= 2 * y1: + targets.append(q * nx * ny + y * nx + x) + a, b = HALO + 2 * (x - x0), HALO + 2 * (y - y0) + sources.append([q * fx * fy + (b + dy) * fx + a + dx + for dy in (0, 1) for dx in (0, 1)]) + return cls((ny, nx), bounds, tau_c, tau_f, + _readonly(kind, np.uint8), + _readonly(ghost_fine, np.uint32), + _readonly(ghost_coarse, np.uint32), + _readonly(targets, np.uint64), + _readonly(sources, np.uint64).reshape(-1, 4)) + + @property + def fine_shape(self) -> tuple[int, int]: + return self.fine_kind.shape + + @property + def fine_origin(self) -> tuple[float, float]: + return self.bounds[0] - HALO / 2, self.bounds[1] - HALO / 2 + + @property + def table_bytes(self) -> int: + return sum(a.nbytes for a in (self.fine_kind, self.ghost_fine, + self.ghost_coarse, self.coalesce_coarse, + self.coalesce_fine)) + + def check_synchronization(self, coarse_steps: int, fine_steps: int, *, phase: str) -> None: + """Check a proposed backend handoff; does not convert or collide data. + + Production storage is g=C(f), while CPUREF returns arrival_f. A future + g-to-g coarse step must collide coarse cells once and fine cells twice; + fine halo ghosts never collide. Coalescence must precede coarse collision. + Converting reference f to g for comparison is an observation operation, + not an additional production collision. Counts here are completed steps. + """ + if phase != "postcollision_g": + raise ValueError("patch backend handoff requires postcollision_g") + for value in (coarse_steps, fine_steps): + if isinstance(value, (bool, np.bool_)) or not isinstance(value, Integral) or value < 0: + raise ValueError("step counts must be nonnegative integers") + if fine_steps != 2 * coarse_steps: + raise ValueError("levels are not synchronized: fine_steps must equal 2*coarse_steps") + + +@dataclass(frozen=True, eq=False) +class FineBoundary2D: + """Fixed local projection plus auditable, finite transfer dependency sets. + + transfer_reads rows are (hop, direction, cell), traced from coalescence at + hop 2 back to exploded G at hop 0. No previous coarse interval is traversed. + Cell footprints also include the whole explosion write/donor tables. + """ + objects: tuple + flags: np.ndarray + curved_host: tuple + boundary_cells: np.ndarray + transfer_fine_cells: np.ndarray + transfer_coarse_cells: np.ndarray + transfer_reads: np.ndarray + + +def _transfer_footprint(patch): + ny, nx = patch.fine_shape + cells = nx * ny + pending = [(2, *divmod(int(slot), cells)) for slot in patch.coalesce_fine.ravel()] + reads = set() + while pending: + hop, q, k = pending.pop() + if (hop, q, k) in reads: + continue + reads.add((hop, q, k)) + if hop == 0: + continue + y, x = divmod(k, nx) + kind = patch.fine_kind[y, x] + if kind == OUTER_GHOST: + raise ValueError("Transfer reads an invalid post-hop outer ghost") + # A collided output depends on all nine arrivals; ghosts only stream. + for direction in (range(9) if kind == ACTIVE else (q,)): + sx, sy = x - D2Q9_EX[direction], y - D2Q9_EY[direction] + if not (0 <= sx < nx and 0 <= sy < ny): + raise ValueError("Transfer dependency leaves the fine allocation") + pending.append((hop - 1, direction, sy * nx + sx)) + fine = set(map(int, patch.ghost_fine)) | {k for _, _, k in reads} + coarse = set(map(int, patch.ghost_coarse)) | set(map( + int, patch.coalesce_coarse % np.prod(patch.coarse_shape))) + return fine, coarse, _readonly(sorted(reads), np.int64) + + +def prepare_fine_boundary(cfg, patch, objects, *, wall_aware=False): + """Project fixed circles and reject boundary/transfer overlap before CUDA. + + Input centers/radii are in coarse physical coordinates (dx=dt=1). + SurfaceActuation retains its lattice-velocity meaning: input coarse units + are converted into independent fine-local objects with dx=dt=.5. + Original objects are not retained or mutated. No runtime geometry updates. + """ + if not objects: + return None + checked = replace(cfg) + checked.validate() + if (patch.coarse_shape != (cfg.ny, cfg.nx) or not checked.is_d2q9 + or checked.collision != "SRT" or checked.streaming != "double_buffer" + or checked.store_precision != "FP32" or checked.ddf_shifting or checked.les_enabled + or patch.tau_coarse != 3*checked.viscosity+.5): + raise ValueError("Fine boundary requires the matching fixed-patch configuration") + from ..body.geometry.circle import CircleGeometry + from ..body.objects import SimObject + from ..body.coupling.soa_packer import pack_cut_links_to_soa + + ny, nx = patch.fine_shape + fine_cfg = replace(cfg, nx=nx, ny=ny, viscosity=2 * cfg.viscosity) + flags = np.full(nx * ny, FLUID, np.uint16) + active = patch.fine_kind.ravel() == ACTIVE + local, links = [], [] + x0, y0, x1, y1 = patch.bounds + dx = dt = .5 + for bid, obj in enumerate(objects): + if (type(obj.geometry) is not CircleGeometry or obj.is_sensor or obj.is_force_region + or any(getattr(obj.state, name) != 0 for name in ("vx", "vy", "omega")) + or any(vars(obj.control).values())): + raise ValueError("Fine boundaries require stationary circles without other roles") + cx, cy = obj.geometry.center + radius = obj.geometry._radius + if (not np.isfinite([cx, cy, radius]).all() or radius <= 0 + or cx-radius < x0 or cx+radius > x1 or cy-radius < y0 or cy+radius > y1): + raise ValueError("Entire fine circle must lie inside the physical patch") + center = tuple((v-o)/dx-.5 for v, o in zip((cx, cy), patch.fine_origin)) + spec = obj.surface_actuation + if spec is not None: + spec = replace(spec, normal_velocity=spec.normal_velocity * dt / dx) + item = SimObject(bid, geometry=CircleGeometry(*center, radius/dx), + center=center, radius=radius/dx, surface_actuation=spec) + item.validate_surface_actuation(fine_cfg) + mask = item.get_flag_mask(nx, ny) + solid = (mask & SOLID) != 0 + if not solid.any() or np.any(solid & ~active) or np.any(flags[solid] & SOLID): + raise ValueError("Fine solids must be resolved, disjoint and entirely active") + flags[solid] = mask[solid] + local.append(item) + for item in local: + cut = item.geometry.build_cut_links(nx, ny, domain_flags=flags) + if not cut: + raise ValueError("Fine circle has no resolved cut links") + if item.surface_actuation is not None: + _, aperture = item.surface_actuation.project( + np.array([link.rx for link in cut], np.float32), + np.array([link.ry for link in cut], np.float32)) + if not aperture.any(): + raise ValueError("Surface aperture hits no resolved cut links") + for link in cut: + link.body_id = item.obj_id + links.extend(cut) + + boundary, writers, receivers = set(), set(), set() + for link in links: + r, q = link.fluid_idx, link.dir + y, x = divmod(r, nx) + sx, sy = x + D2Q9_EX[q], y + D2Q9_EY[q] + s = sy * nx + sx + if not (active[r] and flags[r] & FLUID and active[s] and flags[s] & SOLID): + raise ValueError("Cut-link receiver/source must be active fluid/solid") + writer = (s, q+1 if q % 2 else q-1) + if writer in writers: + raise ValueError("Duplicate fine boundary writer") + writers.add(writer) + receivers.add((r, q)) + boundary.update((r, s)) + if link.q < .5: + dx_, dy_ = x - D2Q9_EX[q], y - D2Q9_EY[q] + if not (0 <= dx_ < nx and 0 <= dy_ < ny and active[dy_*nx+dx_]): + raise ValueError("q<.5 donor crosses the fine active interface; fallback cannot hide it") + boundary.add(dy_*nx+dx_) + # Reverse coverage: every fluid pull across a solid neighbor needs a writer. + required = set() + for r in np.flatnonzero(active & ((flags & FLUID) != 0)): + for q in range(1, 9): + if flags[r + D2Q9_EX[q] + nx*D2Q9_EY[q]] & SOLID: + required.add((int(r), q)) + if required != receivers: + raise ValueError("Fine solid pull coverage does not match the cut links") + fine, coarse, reads = _transfer_footprint(patch) + if fine & (boundary | set(map(int, np.flatnonzero(flags & SOLID)))): + raise ValueError("Boundary stencil overlaps the two-hop transfer footprint") + for k in coarse: + y, x = divmod(k, cfg.nx) + if x0 <= x < x1 and y0 <= y < y1: + raise ValueError("Transfer coarse donor/target is not an uncovered bulk cell") + # Coarse reconstruction targets y=1,ny-2 pull from y=0..2,ny-3..ny-1. + if wall_aware and (y <= 2 or y >= cfg.ny-3): + raise ValueError("Transfer footprint touches the coarse wall reconstruction stencil") + return FineBoundary2D(tuple(local), _readonly(flags, np.uint16), + tuple(_readonly(a, a.dtype) for a in pack_cut_links_to_soa(links)), + _readonly(sorted(boundary), np.uint32), + _readonly(sorted(fine), np.uint32), + _readonly(sorted(coarse), np.uint32), reads) diff --git a/src/CelerisLab/lbm/refinement_stepper.py b/src/CelerisLab/lbm/refinement_stepper.py new file mode 100644 index 0000000..f38aaa0 --- /dev/null +++ b/src/CelerisLab/lbm/refinement_stepper.py @@ -0,0 +1,302 @@ +# CelerisLab/lbm/refinement_stepper.py +"""Opt-in fixed patch backend with optional stationary fine-interior circles. + +The caller owns the CudaContext and must close this backend before that context. +Only synchronized postcollision_g leaf populations are public. No checkpoint, +snapshot, forcing, runtime parameter mutation or dynamic refinement is offered. +""" +from dataclasses import dataclass, replace +from numbers import Integral +import struct +from threading import RLock + +import numpy as np +import pycuda.driver as cuda + +from ..cuda import compiler_v2 as compiler +from .refinement import ACTIVE, HALO, FixedPatch2D, prepare_fine_boundary +from .descriptors import SOLID + + +@dataclass +class _Level: + shape: tuple + module: object + current: object + spare: object + kernel: object + grid: tuple + + +class FixedPatchStepper: + """D2Q9/SRT/FP32, one unforced internal 2:1 patch in an x-periodic box. + + ``coarse_boundary`` is ``periodic`` (default) or ``slot_reservoir``: + a bottom on-site slot and top fixed-density reservoir from ``cfg``. + Periodic mode does not apply the wall/slot boundary fields in ``cfg``. + Outer walls remain coarse; fine halo donors must be ordinary fluid. + ``fine_objects`` accepts fixed circles in coarse physical coordinates; + velocities are in coarse lattice units (dx=dt=1). Independent local copies + are projected at fine cell centers. Their complete stencil must be active + and disjoint from the actual two-hop transfer footprint. + These use the edge/cell-center frame (coarse centers I+.5), whereas ordinary Simulation objects use the integer-node frame. + + ``initialize(coarse_g, fine_g)`` accepts physical Q-major (9,ny,nx) arrays + already in postcollision phase. It does not collide or build channel flags. + Every ``step()`` publishes only after both levels complete on one stream. + On a device/launch failure the instance becomes unusable: close and recreate + it, rather than pretending the partially advanced fine buffer rolled back. + """ + def __init__(self, cfg, bounds, *, context, coarse_boundary="periodic", fine_objects=()): + if HALO != 4: + raise RuntimeError("refinement.cu requires a four-cell fine halo") + if coarse_boundary not in ("periodic", "slot_reservoir"): + raise ValueError("coarse_boundary must be periodic or slot_reservoir") + self.patch = FixedPatch2D.from_config(cfg, bounds) + self._wall_aware = coarse_boundary == "slot_reservoir" + if self._wall_aware: + if (cfg.y_wall_bc != "on_site_velocity" or cfg.top_wall_bc != "reservoir" + or cfg.velocity != 0 or not 2 <= cfg.wall_slot_width < cfg.nx): + raise ValueError("slot_reservoir requires on_site_velocity, reservoir, " + "velocity=0 and 2 <= wall_slot_width < nx") + if self.patch.bounds[1] < 4 or self.patch.bounds[3] > cfg.ny - 4: + raise ValueError("slot_reservoir requires y0 >= 4 and y1 <= ny-4") + self.fine_boundary = prepare_fine_boundary( + cfg, self.patch, tuple(fine_objects), wall_aware=self._wall_aware) + self._fine_leaf = self.patch.fine_kind == ACTIVE + if self.fine_boundary is not None: + self._fine_leaf &= (self.fine_boundary.flags.reshape(self.patch.fine_shape) & SOLID) == 0 + self.ctx = context + self._lock = RLock() + self._status = "new" + self._step_count = 0 + self._allocations = [] + self._levels = [] + self._stream = None + self._fine_flags = self._fine_curved = self._curved = None + self._block = (cfg.threads_per_block, 1, 1) + fine_cfg = replace(cfg, nx=self.patch.fine_shape[1], ny=self.patch.fine_shape[0], + viscosity=2 * cfg.viscosity) + if self._wall_aware: + fine_cfg = replace(fine_cfg, y_wall_bc="bounce_back", top_wall_bc="velocity", + wall_slot_width=0, wall_normal_velocity=0.0, velocity=0.0) + with self.ctx.activate(): + try: + self._stream = cuda.Stream() + for config, name, tau in ((replace(cfg), "PatchCoarseStep", self.patch.tau_coarse), + (fine_cfg, "PatchFineStep", self.patch.tau_fine)): + config.validate() + if not np.isclose(config.omega, 1 / tau, rtol=1e-14, atol=0): + raise RuntimeError("level omega/tau mismatch") + n_objects = (len(self.fine_boundary.objects) + if name == "PatchFineStep" and self.fine_boundary is not None else 0) + build = compiler.generate_config(config, n_objects=n_objects) + try: + ptx = compiler.compile_kernel(self.ctx.sm_arch, source="refinement.cu", + build_dir=build) + module = compiler.load_module(ptx) + finally: + compiler.cleanup_build_dir(build) + ptr, size = module.get_global("d_params") + if size != struct.calcsize("=ffffI"): + raise RuntimeError("LBMParams host/device layout mismatch") + cuda.memcpy_htod(ptr, struct.pack("=ffffI", config.omega, 0, 0, 0, n_objects)) + shape = (9, config.ny, config.nx) + size = int(np.prod(shape)) * 4 + self._levels.append(_Level(shape, module, self._allocate(size), + self._allocate(size), module.get_function(name), + ((config.nx + self._block[0] - 1) // self._block[0], + config.ny, 1))) + p = self.patch + order = np.argsort(p.coalesce_coarse) + self._ghosts = self._upload(p.ghost_fine) + self._donors = self._upload(p.ghost_coarse) + self._targets = self._upload(p.coalesce_coarse[order]) + self._children = self._upload(p.coalesce_fine[order]) + self._explode = self._levels[1].module.get_function("PatchExplode") + if self.fine_boundary is not None: + from .curved_links import CurvedLinkSoA + boundary = self.fine_boundary + self._fine_flags = self._upload(boundary.flags) + # SoA owns its allocations; never also register them in _allocations. + self._fine_curved = CurvedLinkSoA() + self._fine_curved.assign_host(*boundary.curved_host) + self._fine_curved.assign_surface(boundary.objects, fine_cfg, boundary.flags) + self._fine_curved.upload() + self._curved = self._levels[1].module.get_function("CurvedBoundaryKernel") + except BaseException: + try: + self.close() + except Exception: + pass # Preserve the original construction failure. + raise + + def _allocate(self, size): + pointer = cuda.mem_alloc(size) + self._allocations.append(pointer) + return pointer + + def _upload(self, array): + array = np.ascontiguousarray(array) + pointer = self._allocate(array.nbytes) + cuda.memcpy_htod(pointer, array) + return pointer + + def _require(self, *states): + if self._status not in states: + raise RuntimeError(f"Fixed patch is {self._status}; expected {' or '.join(states)}") + + @property + def step_count(self): + return self._step_count + + def _mask_nonleaves(self, coarse, fine): + x0, y0, x1, y1 = self.patch.bounds + coarse[:, y0:y1, x0:x1] = np.nan + if self._wall_aware: + coarse[:, (0, -1), :] = np.nan + fine[:, ~self._fine_leaf] = np.nan + + def initialize(self, coarse_g, fine_g): + """Upload an explicit G_0; invalid non-leaves are discarded.""" + with self._lock, self.ctx.activate(): + self._require("new", "ready") + arrays = [np.array(a, dtype=np.float32, order="C", copy=True) + for a in (coarse_g, fine_g)] + if any(a.shape != level.shape for a, level in zip(arrays, self._levels)): + raise ValueError("Expected Q-major (9,ny,nx) arrays for both levels") + x0, y0, x1, y1 = self.patch.bounds + leaf = np.ones(self.patch.coarse_shape, dtype=bool) + leaf[y0:y1, x0:x1] = False + if self._wall_aware: + leaf[[0, -1]] = False + for array, mask in zip(arrays, (leaf, self._fine_leaf)): + if not np.isfinite(array[:, mask]).all() or (array[:, mask].sum(axis=0) <= 0).any(): + raise ValueError("Leaf populations must be finite with positive density") + self._mask_nonleaves(*arrays) + try: + for level, array in zip(self._levels, arrays): + cuda.memcpy_htod(level.current, array) + cuda.memcpy_htod(level.spare, array) + if self._fine_curved is not None and self._fine_curved.link_exchange_gpu is not None: + cuda.memset_d8(self._fine_curved.link_exchange_gpu, 0, self._fine_curved.count * 12) + self._step_count = 0 + self._status = "ready" + except BaseException: + self._status = "failed" + raise + + def _apply_fine_boundary(self, buffer, **kw): + """Stage solid source slots; exchange buffer holds this fine substep only.""" + cl = self._fine_curved + if cl is None: + if self.fine_boundary is not None: + raise RuntimeError("Fine boundary device buffers are missing") + return + if self._curved is None or cl.fluid_idx_gpu is None: + raise RuntimeError("Fine boundary device buffers are missing") + driven = any(obj.surface_actuation is not None for obj in self.fine_boundary.objects) + if driven and (cl.surface_velocity_gpu is None or cl.link_exchange_gpu is None): + raise RuntimeError("Surface boundary device buffers are missing") + self._curved(buffer, cl.fluid_idx_gpu, cl.dir_gpu, cl.q_gpu, + cl.rx_gpu, cl.ry_gpu, cl.rz_gpu, cl.fallback_class_gpu, cl.body_id_gpu, + np.uintp(0), np.uintp(0), np.uint32(cl.count), + cl.surface_velocity_gpu if cl.surface_velocity_gpu is not None else np.uintp(0), + cl.link_exchange_gpu if cl.link_exchange_gpu is not None else np.uintp(0), + grid=((cl.count + self._block[0] - 1) // self._block[0], 1), **kw) + + def step(self, n=1): + """Advance G→G: one explosion, two fine collisions, one coarse collision.""" + if isinstance(n, bool) or not isinstance(n, Integral) or n < 0: + raise ValueError("n must be a nonnegative integer") + with self._lock, self.ctx.activate(): + self._require("ready") + c, f = self._levels + p = self.patch + kw = dict(block=self._block, stream=self._stream) + try: + self._status = "advancing" + for _ in range(n): + self._explode(c.current, f.current, self._ghosts, self._donors, + np.uint32(len(p.ghost_fine)), np.uint64(np.prod(p.coarse_shape)), + grid=((len(p.ghost_fine) + self._block[0] - 1) // self._block[0], 1), + **kw) + flags = self._fine_flags if self._fine_flags is not None else np.uintp(0) + if self.fine_boundary is not None and self._fine_flags is None: + raise RuntimeError("Fine boundary flags are missing") + self._apply_fine_boundary(f.current, **kw) + f.kernel(f.current, f.spare, flags, grid=f.grid, **kw) + self._apply_fine_boundary(f.spare, **kw) + f.kernel(f.spare, f.current, flags, grid=f.grid, **kw) + c.kernel(c.current, c.spare, f.current, self._targets, self._children, + np.uint32(len(p.coalesce_coarse)), + *(np.uint32(v) for v in p.bounds), np.uint32(self._wall_aware), + grid=c.grid, **kw) + self._stream.synchronize() + c.current, c.spare = c.spare, c.current + self._step_count += 1 + self._status = "ready" + except BaseException: + self._status = "failed" + raise + + def download(self): + """Return two synchronized G arrays, NaN on all non-leaves; no host cache.""" + with self._lock, self.ctx.activate(): + self._require("ready") + arrays = [np.empty(level.shape, np.float32) for level in self._levels] + try: + for level, array in zip(self._levels, arrays): + cuda.memcpy_dtoh(array, level.current) + except BaseException: + self._status = "failed" + raise + self._mask_nonleaves(*arrays) + return tuple(arrays) + + def close(self): + """Release owned resources; retain them for retry if context activation fails.""" + with self._lock: + if self._status == "closed": + return + self._status = "failed" + with self.ctx.activate(): + self._status = "closed" + error = None + try: + if self._stream is not None: + self._stream.synchronize() + except Exception as exc: + error = exc + self._levels.clear() + self._explode = None + if getattr(self, "_fine_curved", None) is not None: + try: + self._fine_curved.free() + except Exception as exc: + error = error or exc + self._fine_curved = None + self._curved = self._fine_flags = None + for pointer in self._allocations: + try: + pointer.free() + except Exception as exc: + error = error or exc + self._allocations.clear() + self._stream = None + if error is not None: + raise error + + def __enter__(self): + return self + + def __exit__(self, *exc): + try: + self.close() + except Exception: + if exc[0] is None: + raise + + +# Preserve the original entry point and its periodic default. +PeriodicPatchStepper = FixedPatchStepper diff --git a/src/CelerisLab/lbm/stepper.py b/src/CelerisLab/lbm/stepper.py index 347f48a..940e62d 100644 --- a/src/CelerisLab/lbm/stepper.py +++ b/src/CelerisLab/lbm/stepper.py @@ -24,6 +24,8 @@ class LBMStepper: if self._esopull: self.step_fn = module.get_function("EsoPullStep") self.init_fn = module.get_function("InitEsoPull") + self.prepare_boundaries_fn = module.get_function("PrepareEsoPullBoundaries") + self.boundary_values_gpu = cuda.mem_alloc((6 * cfg.nx + 2 * cfg.ny * cfg.nq) * 4) else: self.step_fn = module.get_function("OneStep") self.init_fn = module.get_function("InitTubeFlow_v2") @@ -41,6 +43,7 @@ class LBMStepper: else: self.grid = ((cfg.nx + tpb - 1) // tpb, cfg.ny, cfg.nz) + self.boundary_grid = ((max(cfg.nx, cfg.ny) + tpb - 1) // tpb, 2, 1) self._step_count = 0 def initialize(self): @@ -54,7 +57,9 @@ 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(step_id=self._step_count, force=True) + # Initialization publishes a new t=0 state, including EsoPull parity. + self._step_count = 0 + f.download_ddf(step_id=0, force=True) def step( self, @@ -71,12 +76,17 @@ class LBMStepper: launch_kw["stream"] = stream for _ in range(n): + if self._esopull: + self.prepare_boundaries_fn( + f.ddf_gpu, self.boundary_values_gpu, np.uint64(self._step_count), + block=self.block, grid=self.boundary_grid, **launch_kw, + ) self._launch_curved(action_gpu, obs_gpu, **launch_kw) if self._esopull: self.step_fn( f.ddf_gpu, f.flag_gpu, - action_gpu, obs_gpu, + action_gpu, obs_gpu, self.boundary_values_gpu, np.uint64(self._step_count), block=self.block, grid=self.grid, **launch_kw, @@ -137,12 +147,17 @@ class LBMStepper: ) else: assert cl.fluid_idx_gpu is not None + if cl.surface_velocity is not None and ( + cl.surface_velocity_gpu is None or cl.link_exchange_gpu is None): + raise RuntimeError("Surface boundary device buffers are missing") self.curved_fn( f.ddf_gpu, cl.fluid_idx_gpu, cl.dir_gpu, cl.q_gpu, cl.rx_gpu, cl.ry_gpu, cl.rz_gpu, cl.fallback_class_gpu, cl.body_id_gpu, action_gpu, obs_gpu, np.uint32(f.n_curved), + cl.surface_velocity_gpu if cl.surface_velocity_gpu is not None else np.uintp(0), + cl.link_exchange_gpu if cl.link_exchange_gpu is not None else np.uintp(0), block=(tpb, 1, 1), grid=grid_c, **launch_kw, ) diff --git a/src/CelerisLab/simulation.py b/src/CelerisLab/simulation.py index 6a1b3fd..64cf484 100644 --- a/src/CelerisLab/simulation.py +++ b/src/CelerisLab/simulation.py @@ -18,12 +18,13 @@ from __future__ import annotations from typing import Dict, Optional, Tuple, Any import os +import operator import numpy as np import pycuda.driver as cuda from .config import LBMConfig, BodyConfig, load_lbm_config, load_body_config -from .cuda.context import CudaContext +from .cuda.context import CudaContext, context_bound from .cuda import compiler_v2 as compiler from .lbm.field import LBMField from .lbm.stepper import LBMStepper @@ -45,39 +46,50 @@ class Simulation: device_id: int = 0): # Internal stream for async DTOH coordination self._stream: Optional[cuda.Stream] = None + self._build_dir = None + self._closed = False # Load configs self.lbm_cfg = load_lbm_config(lbm_config_path) + self.lbm_cfg.validate_runtime() self.body_cfg = load_body_config(body_config_path) # CUDA context self.ctx = CudaContext(device_id) - arch = self._resolve_compile_arch() + try: + with self.ctx.activate(): + arch = self._resolve_compile_arch() + self._build_dir = compiler.generate_config(self.lbm_cfg, n_objects=0) + self._ptx_path = compiler.compile_kernel(arch=arch, build_dir=self._build_dir) + self._module = compiler.load_module(self._ptx_path) + self._compiled_n_objects = 0 + self._allocation_signature = self.lbm_cfg.allocation_signature() + self._compiled_viscosity = self.lbm_cfg.viscosity - # Compile kernel - compiler.generate_config(self.lbm_cfg, n_objects=0) - self._compiled_n_objects: int = 0 - self._ptx_path = compiler.compile_kernel(arch=arch) - self._module = compiler.load_module(self._ptx_path) - - # LBM field & stepper - self.field = LBMField(self.lbm_cfg, self._module) - self.stepper = LBMStepper(self.field, self._module, self.lbm_cfg) - - # Body manager - self.bodies = ObjectManager( - self.lbm_cfg.nx, self.lbm_cfg.ny, self.lbm_cfg.nz, - self.lbm_cfg.nq, self.lbm_cfg, - ) - - self._initialized = False - self._load_objects_from_body_config() - if self.bodies.count > 0: - self.recompile() - self._assert_object_count_contract(expected_count=self.bodies.count) + self.field = LBMField(self.lbm_cfg, self._module) + self.field.ctx = self.ctx + self.stepper = LBMStepper(self.field, self._module, self.lbm_cfg) + self.bodies = ObjectManager( + self.lbm_cfg.nx, self.lbm_cfg.ny, self.lbm_cfg.nz, + self.lbm_cfg.nq, self.lbm_cfg, + ) + self._initialized = False + self._load_objects_from_body_config() + if self.bodies.count > 0: + self.recompile() + self._assert_object_count_contract(expected_count=self.bodies.count) + except Exception: + try: + self.close() + except Exception: + # Construction must report its original failure. Cleanup is + # best effort here; an explicit close failure remains retryable. + pass + raise # -- Stream access ------------------------------------------------------- @property + @context_bound def stream(self) -> cuda.Stream: """Internal CUDA stream for async DTOH coordination.""" if self._stream is None: @@ -86,7 +98,7 @@ class Simulation: # -- Object management --------------------------------------------------- def add_body(self, type: str = "circle", center=None, - radius: float = 0.0, **kwargs) -> int: + radius: float = 0.0, *, surface_actuation=None, **kwargs) -> int: """Add a simulation body. Returns body_id. Before ``initialize()`` the object is added to the formal registry @@ -98,12 +110,20 @@ class Simulation: (future: ``"polygon"``, ``"mesh"``). center: (x, y) center coordinates (lattice units). radius: body radius (lattice units). + surface_actuation: Fixed circle aperture dict/SurfaceActuation: + angle, half_width (radians), normal_velocity (positive outward). + D2Q9/FP32/unshifted/DB only; set before initialization or stage + a new body. No runtime motion or amplitude setter. **kwargs: passed to the geometry constructor. Returns: int body_id (formal id before initialize; -1 if staged). """ + if self.lbm_cfg.y_wall_bc == "on_site_velocity": + raise ValueError("on_site_velocity currently requires an object-free domain") type_lower = str(type).strip().lower() + if surface_actuation is not None and type_lower != "circle": + raise ValueError("Surface actuation is only supported on circles") if type_lower == "circle": from .body.geometry.circle import CircleGeometry @@ -111,7 +131,7 @@ class Simulation: geom = CircleGeometry(cx, cy, radius) obj = SimObject(obj_id=-1, geometry=geom, center=center, radius=radius, - is_sensor=False) + is_sensor=False, surface_actuation=surface_actuation) elif type_lower == "sensor": from .body.geometry.circle import CircleGeometry cx, cy = float(center[0]), float(center[1]) @@ -142,7 +162,8 @@ class Simulation: Before ``initialize()`` the object is removed from the formal registry immediately. After ``initialize()`` the removal is staged and will be - committed at the next ``sync_bodies()``. + committed at the next ``sync_bodies()``. Removal compacts numeric IDs; + obtain the remaining IDs again after immediate removal or runtime sync. Args: id: body_id from ``add_body()``. @@ -174,6 +195,8 @@ class Simulation: center = tuple(spec.get("center", ())) radius = float(spec.get("radius", 0.0)) omega = spec.get("omega", None) + if spec.get("surface_actuation") is not None and obj_type != "cylinder": + raise ValueError("Surface actuation is only supported on cylinders in body config") if obj_type == "": raise ValueError( @@ -185,9 +208,10 @@ class Simulation: ) if obj_type == "cylinder": - body_id = self.add_body("circle", center=center, radius=radius) + body_id = self.add_body("circle", center=center, radius=radius, + surface_actuation=spec.get("surface_actuation")) if omega is not None: - self.bodies.get(body_id).state.omega = float(omega) + self.bodies.set_body_state(body_id, omega=float(omega)) elif obj_type == "sensor": self.add_body("sensor", center=center, radius=radius) else: @@ -210,6 +234,10 @@ class Simulation: vy: translation velocity y (future). ``None`` = no change. """ obj = self.bodies.get(id) + proposed = dict(vars(obj.state)) + proposed.update({k: float(v) for k, v in (("omega", omega), ("vx", vx), ("vy", vy)) + if v is not None}) + obj.validate_surface_actuation(self.lbm_cfg, proposed) changed = False if omega is not None: obj.state.omega = float(omega) @@ -276,6 +304,19 @@ class Simulation: raise RuntimeError("No obs buffer. Call run() first.") return self.bodies.read_sensor(id, normalize=normalize) + @context_bound + def read_surface_exchange(self, id: int): + """Latest step's wall/aperture (Fx,Fy,Tz) exchange, not total nozzle force. + + Unlike read_body this is not accumulated/time-averaged. It resets after + geometry sync or checkpoint load; advance before interpreting it. + """ + if self.bodies.get(id).surface_actuation is None: + raise ValueError("Body has no surface actuation") + self.stream.synchronize() + return self.field.curved.surface_exchange(id) + + @context_bound def read_body(self, id: int, *, stream: cuda.Stream | None = None, normalize: bool = True): @@ -301,6 +342,7 @@ class Simulation: stream.synchronize() return self.bodies.read_body(id, normalize=normalize) + @context_bound def read_bodies(self, *, stream: cuda.Stream | None = None) -> np.ndarray: """Return all bodies' telemetry as a flat float32 array. @@ -320,31 +362,39 @@ class Simulation: return self.bodies._obs_array() # -- Compilation --------------------------------------------------------- + @context_bound def _recompile(self, expected_n_objects: int | None = None) -> None: """Re-generate config headers and recompile the kernel. If *expected_n_objects* is ``None``, uses ``self.bodies.count``. Preserves ``step_count`` and GPU DDF/flag buffers across the - module reload. Does NOT preserve ``action_gpu`` / ``obs_gpu`` - (those are recreated by ``sync_to_gpu`` or ``run()``). + module reload. Packed action/observation buffers are managed separately + by the topology synchronization caller. This is safe to call before or after ``initialize()``. """ if expected_n_objects is None: expected_n_objects = self.bodies.count + self.lbm_cfg.validate_runtime() + self._assert_allocation_contract() arch = self._resolve_compile_arch() - compiler.generate_config(self.lbm_cfg, n_objects=expected_n_objects) + next_dir = compiler.generate_config(self.lbm_cfg, n_objects=expected_n_objects) + try: + next_ptx = compiler.compile_kernel(arch=arch, build_dir=next_dir) + next_module = compiler.load_module(next_ptx) + except Exception: + compiler.cleanup_build_dir(next_dir) + raise + old_dir = self._build_dir + self._build_dir, self._ptx_path, self._module = next_dir, next_ptx, next_module self._compiled_n_objects = int(expected_n_objects) - # Remove stale PTX to prevent PyCUDA module handle conflicts. - import os as _os - if _os.path.exists(self._ptx_path): - _os.remove(self._ptx_path) - self._ptx_path = compiler.compile_kernel(arch=arch) - 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() + if self.lbm_cfg.viscosity != self._compiled_viscosity: + self.field._runtime_params["omega"] = self.lbm_cfg.omega + self._compiled_viscosity = self.lbm_cfg.viscosity self.field._upload_params() _prev_step_count = self.stepper._step_count @@ -354,6 +404,7 @@ class Simulation: self.stepper._step_count = _prev_step_count self._assert_object_count_contract(expected_count=expected_n_objects) + compiler.cleanup_build_dir(old_dir) def recompile(self): """Re-generate config headers and recompile kernel. @@ -369,11 +420,12 @@ class Simulation: self._recompile() # -- Initialization ------------------------------------------------------ + @context_bound def initialize(self): """Initialize flow field and sync objects to GPU.""" # Recompile if objects were added after construction. - if self.bodies.count > 0 and \ - self._read_generated_n_objects() != self.bodies.count: + self._assert_allocation_contract() + if self._read_generated_n_objects() != self.bodies.count: self.recompile() self.field.flag = self.field.build_channel_flags() @@ -390,6 +442,7 @@ class Simulation: self._initialized = True # -- Runtime body topology sync ------------------------------------------- + @context_bound def sync_bodies(self) -> None: """Apply pending body edits (add/remove) to a running simulation. @@ -400,6 +453,7 @@ class Simulation: if not self._initialized: raise RuntimeError("Call initialize() before sync_bodies()") + self._assert_allocation_contract() # 1. Synchronize GPU self.stream.synchronize() @@ -470,6 +524,7 @@ class Simulation: self._assert_runtime_contracts() # -- Stepping ------------------------------------------------------------ + @context_bound def run(self, steps: int, *, stream: cuda.Stream | None = None, upload_act: bool = True, @@ -479,18 +534,24 @@ class Simulation: """Advance simulation by *steps* time steps. Args: - steps: Number of LBM steps. - stream: CUDA stream for async operations. ``None`` = internal stream. + steps: Nonnegative integer number of LBM steps. + stream: CUDA stream for async operations. ``None`` = internal stream. + External streams must be created under ``with sim.ctx.activate():``. upload_act: If True, upload host action array to ``action_gpu`` before the step group. sync_obs: If True, download ``obs_gpu`` to host pinned buffer after the step group. zero_obs: If True (default), zero all obs segments (force, torque, sensor) on GPU and reset the step accumulator before stepping. - checkpoint_interval: If >0, save checkpoint every N steps. + checkpoint_interval: Nonnegative integer; if >0, save every N steps. """ if not self._initialized: raise RuntimeError("Call initialize() first") + steps = operator.index(steps) + checkpoint_interval = operator.index(checkpoint_interval) + if steps < 0 or checkpoint_interval < 0: + raise ValueError("steps and checkpoint_interval must be nonnegative") + self._assert_runtime_contracts() # Discard any uncommitted body edits before stepping. if self.bodies.has_pending_edit(): self.bodies.clear_pending_edits() @@ -507,7 +568,6 @@ class Simulation: if zero_obs: self.bodies.zero_obs_async(stream) - self._assert_runtime_contracts() if checkpoint_interval > 0: done = 0 while done < steps: @@ -518,8 +578,10 @@ class Simulation: obs_gpu=self.bodies.obs_gpu, stream=stream, ) + self.bodies.increment_obs_steps(batch) done += batch if done < steps or done % checkpoint_interval == 0: + stream.synchronize() self.save_checkpoint() else: self.stepper.step( @@ -529,9 +591,9 @@ class Simulation: stream=stream, ) - # Accumulate step count for obs time-normalisation - if steps > 0: - self.bodies.increment_obs_steps(steps) + # Count each batch before any checkpoint captures observations. + if steps > 0: + self.bodies.increment_obs_steps(steps) # Async download obs if sync_obs: @@ -544,10 +606,12 @@ class Simulation: self.run(n) # -- Data access --------------------------------------------------------- + @context_bound def get_macroscopic(self) -> Dict[str, np.ndarray]: """Download DDF and return rho, ux, uy [, uz].""" return self.field.get_macroscopic(step_id=self.stepper.step_count) + @context_bound def get_ddf(self) -> np.ndarray: self.field.download_ddf(step_id=self.stepper.step_count) return self.field.ddf.copy() @@ -556,6 +620,7 @@ class Simulation: return self.field.flag.copy() # -- Runtime parameter updates ------------------------------------------- + @context_bound def update_runtime_params(self, **kwargs): """Update __constant__ d_params without recompiling. @@ -564,29 +629,54 @@ class Simulation: self.field.update_params(**kwargs) # -- Snapshots ----------------------------------------------------------- + @context_bound def snapshot(self): - self.field.snapshot() + """Snapshot the flow field; restore keeps the current clock and controls.""" + self.field.snapshot(step_id=self.stepper.step_count) + @context_bound def restore(self): - self.field.restore() + """Restore the saved flow field at the current streaming parity.""" + self.field.restore(step_id=self.stepper.step_count) # -- Checkpoint / Restore (HDF5) ---------------------------------------- + @context_bound def save_checkpoint(self, path: str = None): """Save full simulation state to an HDF5 file.""" from .common.checkpoint import save_checkpoint return save_checkpoint(self.field, self.stepper, self.lbm_cfg, self.bodies, path) + @context_bound def load_checkpoint(self, path: str): """Restore simulation state from an HDF5 checkpoint.""" from .common.checkpoint import load_checkpoint + if not self._initialized: + self.initialize() load_checkpoint(path, self.field, self.stepper, self.lbm_cfg, self.bodies) self._initialized = True # -- Cleanup ------------------------------------------------------------- def close(self): - self.ctx.close() + if getattr(self, "_closed", False): + return + ctx = getattr(self, "ctx", None) + if ctx is not None and ctx.is_active: + raise RuntimeError("Cannot close Simulation during an active CUDA operation") + if ctx is not None and not ctx.is_closed: + with ctx.activate(): + cuda.Context.synchronize() + # Release owned GPU resources while their context is current. + self.stepper = None + self.bodies = None + self.field = None + self._module = None + self._stream = None + ctx.close() + compiler.cleanup_build_dir(getattr(self, "_build_dir", None)) + self._build_dir = None + self._closed = True def __enter__(self): return self @@ -658,9 +748,16 @@ class Simulation: def _assert_runtime_contracts(self) -> None: """Validate object-count and obs layout contracts before stepping.""" + self._assert_allocation_contract() self._assert_object_count_contract() self._assert_obs_layout_contract() + def _assert_allocation_contract(self) -> None: + if self.lbm_cfg.allocation_signature() != self._allocation_signature: + raise ValueError( + "Grid, lattice, streaming or storage configuration changed after " + "allocation. Construct a new Simulation for structural changes.") + # -- Module-level helpers ------------------------------------------------------- @@ -668,9 +765,9 @@ def _build_patch_set(added_solid_mask, released_fluid_mask, curved_fluid_idx, nx=None, ny=None): """Return sorted uint32 array of cell indices that need DDF patching. - Includes a 1-cell halo because ``store_physical_node_esopull`` writes - to neighbor slots (D2Q9 directions 1,3,5,7). Neighbors must be - decoded/encoded together to avoid corrupting the backing layout. + Includes a 1-cell halo for the old-fluid donors used by released-cell + reconstruction. EsoPull same-time encoding itself has disjoint ownership + and does not require a halo to avoid storage collisions. """ s = set() s.update(np.where(added_solid_mask)[0].tolist()) @@ -699,12 +796,14 @@ def _build_patch_set(added_solid_mask, released_fluid_mask, curved_fluid_idx, def _scatter_to_field(field, scratch_host, patch_indices): """Copy decoded physical DDF from scratch buffer into field.ddf.""" nq = field.nq + populations = field.ddf.reshape(nq, -1) for j, idx in enumerate(patch_indices): - field.ddf[idx * nq:(idx + 1) * nq] = scratch_host[j * nq:(j + 1) * nq] + populations[:, idx] = scratch_host[j * nq:(j + 1) * nq] def _gather_from_field(field, scratch_host, patch_indices): """Copy patched physical DDF from field.ddf back to scratch buffer.""" nq = field.nq + populations = field.ddf.reshape(nq, -1) for j, idx in enumerate(patch_indices): - scratch_host[j * nq:(j + 1) * nq] = field.ddf[idx * nq:(idx + 1) * nq] + scratch_host[j * nq:(j + 1) * nq] = populations[:, idx]