diff --git a/src/CelerisLab/body/jets.py b/src/CelerisLab/body/jets.py index 5223cd3..0efe681 100644 --- a/src/CelerisLab/body/jets.py +++ b/src/CelerisLab/body/jets.py @@ -11,6 +11,10 @@ One physical link (receiver, direction) may be driven by at most one jet. Edge rule: the lateral bound |s| < 1 is exclusive (the profile vanishes there for parabolic/quartic); Sector r_min/r_max and Tube 0 <= along <= length are inclusive. Jets are immutable; an amplitude change replaces the Jet. + +``surface_link_velocity`` is the single host-only assignment/validation point +for legacy apertures plus jets (no PyCUDA import), shared by dense and fine +levels. """ from __future__ import annotations @@ -181,3 +185,50 @@ def fine_hits_to_edge_frame(hits_local, fine_origin, dx): def edge_to_node_frame(points): """Coarse cell-edge frame -> Simulation integer-node frame (shift by -1/2).""" return np.asarray(points, np.float64) - .5 + + +def surface_link_velocity(objects, cfg, flags, fluid_idx, dir, body_id, rx, ry, + jets=(), hits=None): + """Per-link prescribed surface velocity (n, 2) FP32 and aperture mask. + + Host-only and storage-neutral: rows follow the given physical cut-link + order. Returns (None, None) when neither an object carries surface + actuation nor a jet is given. Jets need true FP64 ``hits`` (n, 2) and may + not share a link with each other or with a legacy aperture. + """ + enabled = [obj for obj in objects if obj.surface_actuation is not None] + if not enabled and not jets: + return None, None + from ..lbm.descriptors import FLUID, SOLID + if np.any((flags[fluid_idx] & FLUID) == 0) or np.any(flags[fluid_idx] & SOLID): + raise ValueError("Surface cut links must have active fluid receivers") + if len(np.unique(fluid_idx.astype(np.uint64)*cfg.nq+dir)) != len(dir): + raise ValueError("Surface boundary requires unique cut-link writers") + surface_velocity = np.zeros((len(dir), 2), np.float32) + surface_aperture = np.zeros(len(dir), bool) + for obj in enabled: + obj.validate_surface_actuation(cfg) + selected = body_id == obj.obj_id + velocity, aperture = obj.surface_actuation.project(rx[selected], ry[selected]) + if not aperture.any(): + raise ValueError("Surface aperture hits no resolved cut links") + surface_velocity[selected] = velocity + surface_aperture[selected] = aperture + if jets: + if hits is None or len(hits) != len(dir): + raise ValueError("Jets require true hit points row-aligned with the cut links") + if (cfg is None or not cfg.is_d2q9 or cfg.store_precision != "FP32" + or cfg.streaming != "double_buffer" or cfg.ddf_shifting): + raise ValueError("Jets require D2Q9/FP32/unshifted/double_buffer") + shape, jet_id = jet_link_shapes(jets, hits) + driven = jet_id >= 0 + owners = {obj.obj_id: obj for obj in objects} + for bid in np.unique(body_id[driven]): + state = vars(owners[int(bid)].state) + if any(state.get(k, 0.) != 0 for k in ("vx", "vy", "omega")): + raise ValueError("Jets may only drive links of stationary bodies") + if np.any(surface_aperture & driven): + raise ValueError("Jet overlaps a SurfaceActuation aperture") + surface_velocity[driven] = jet_link_velocity(jets, shape, jet_id)[driven] + surface_aperture |= driven + return surface_velocity, surface_aperture diff --git a/src/CelerisLab/lbm/curved_links.py b/src/CelerisLab/lbm/curved_links.py index ad33085..a95f21b 100644 --- a/src/CelerisLab/lbm/curved_links.py +++ b/src/CelerisLab/lbm/curved_links.py @@ -17,6 +17,8 @@ from typing import Optional import numpy as np import pycuda.driver as cuda +from ..body.jets import surface_link_velocity # host-only; re-exported for callers + def _free_alloc(ptr: Optional[cuda.DeviceAllocation]) -> None: if ptr is not None: @@ -42,54 +44,6 @@ def _upload_columns(columns, stream=None): stream.synchronize() -def surface_link_velocity(objects, cfg, flags, fluid_idx, dir, body_id, rx, ry, - jets=(), hits=None): - """Per-link prescribed surface velocity (n, 2) FP32 and aperture mask. - - Host-only and storage-neutral: rows follow the given physical cut-link - order. Returns (None, None) when neither an object carries surface - actuation nor a jet is given. Jets need true FP64 ``hits`` (n, 2) and may - not share a link with each other or with a legacy aperture. - """ - enabled = [obj for obj in objects if obj.surface_actuation is not None] - if not enabled and not jets: - return None, None - from .descriptors import FLUID, SOLID - if np.any((flags[fluid_idx] & FLUID) == 0) or np.any(flags[fluid_idx] & SOLID): - raise ValueError("Surface cut links must have active fluid receivers") - if len(np.unique(fluid_idx.astype(np.uint64)*cfg.nq+dir)) != len(dir): - raise ValueError("Surface boundary requires unique cut-link writers") - surface_velocity = np.zeros((len(dir), 2), np.float32) - surface_aperture = np.zeros(len(dir), bool) - for obj in enabled: - obj.validate_surface_actuation(cfg) - selected = body_id == obj.obj_id - velocity, aperture = obj.surface_actuation.project(rx[selected], ry[selected]) - if not aperture.any(): - raise ValueError("Surface aperture hits no resolved cut links") - surface_velocity[selected] = velocity - surface_aperture[selected] = aperture - if jets: - from ..body.jets import jet_link_shapes, jet_link_velocity - if hits is None or len(hits) != len(dir): - raise ValueError("Jets require true hit points row-aligned with the cut links") - if (cfg is None or not cfg.is_d2q9 or cfg.store_precision != "FP32" - or cfg.streaming != "double_buffer" or cfg.ddf_shifting): - raise ValueError("Jets require D2Q9/FP32/unshifted/double_buffer") - shape, jet_id = jet_link_shapes(jets, hits) - driven = jet_id >= 0 - owners = {obj.obj_id: obj for obj in objects} - for bid in np.unique(body_id[driven]): - state = vars(owners[int(bid)].state) - if any(state.get(k, 0.) != 0 for k in ("vx", "vy", "omega")): - raise ValueError("Jets may only drive links of stationary bodies") - if np.any(surface_aperture & driven): - raise ValueError("Jet overlaps a SurfaceActuation aperture") - surface_velocity[driven] = jet_link_velocity(jets, shape, jet_id)[driven] - surface_aperture |= driven - return surface_velocity, surface_aperture - - @dataclass class CurvedLinkSoA: """Cut-link columns on host with matching ``*_gpu`` device buffers.""" diff --git a/src/CelerisLab/lbm/refinement.py b/src/CelerisLab/lbm/refinement.py index fb25211..34f658b 100644 --- a/src/CelerisLab/lbm/refinement.py +++ b/src/CelerisLab/lbm/refinement.py @@ -292,13 +292,20 @@ def prepare_fine_boundary(cfg, patch, objects, *, wall_aware=False, jets=()): # coalescence; rows 2/ny-3 are ordinary bulk cells only read as donors. if wall_aware and (y <= 1 or y >= cfg.ny-2): raise ValueError("Transfer footprint touches the coarse wall reconstruction stencil") - from ..body.jets import fine_hits_to_edge_frame, jet_link_shapes + from ..body.jets import fine_hits_to_edge_frame, surface_link_velocity hits = fine_hits_to_edge_frame([(link.hit_x, link.hit_y) for link in links], patch.fine_origin, dx) + packed = pack_cut_links_to_soa(links) if jets: - jet_link_shapes(jets, hits) # empty selection/overlap rejected before CUDA + # Same assignment checks as the stepper: empty/overlapping jets and + # jet-vs-SurfaceActuation overlap are rejected before CUDA. + jet_cfg = replace(fine_cfg) + jet_cfg.validate() + fluid_idx, dirs, _, body_id, rx, ry, _, _ = packed + surface_link_velocity(local, jet_cfg, flags, fluid_idx, dirs, body_id, rx, ry, + jets=jets, hits=hits) return FineBoundary2D(tuple(local), _readonly(flags, np.uint16), - tuple(_readonly(a, a.dtype) for a in pack_cut_links_to_soa(links)), + tuple(_readonly(a, a.dtype) for a in packed), _readonly(sorted(boundary), np.uint32), _readonly(sorted(fine), np.uint32), _readonly(sorted(coarse), np.uint32), reads, diff --git a/tests/unit/test_jets.py b/tests/unit/test_jets.py index 52e3590..85e39d9 100644 --- a/tests/unit/test_jets.py +++ b/tests/unit/test_jets.py @@ -276,3 +276,10 @@ def test_fine_boundary_keeps_edge_frame_hits_and_checks_jets(): refinement.prepare_fine_boundary(cfg, patch, [body], jets=[jet, jet]) with pytest.raises(ValueError, match="need fine boundary objects"): refinement.prepare_fine_boundary(cfg, patch, [], jets=[jet]) + # Jet vs legacy SurfaceActuation overlap on the fine level (root's case). + legacy = SimObject(0, geometry=CircleGeometry(*C, R), center=C, radius=R, + surface_actuation=SurfaceActuation(0., .5, .002)) + with pytest.raises(ValueError, match="SurfaceActuation"): + refinement.prepare_fine_boundary(cfg, patch, [legacy], jets=[jet]) + disjoint = Jet(Sector(C, math.pi, .5, R + .01), Radial(C, inward=False)) + refinement.prepare_fine_boundary(cfg, patch, [legacy], jets=[disjoint])