Frank14f/CelerisLab
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Jets: one shared assignment check for dense and fine levels

Browse Source 
Move surface_link_velocity (pure host code) from lbm/curved_links.py to
body/jets.py, re-exported by curved_links, so refinement.py can reuse it
without importing PyCUDA. prepare_fine_boundary now runs the same
assignment validation as the steppers when jets are given, rejecting
jet-vs-SurfaceActuation overlap (root's C(12.25,12.25) R2.15 case) as
well as empty/overlapping jets. Test the overlap and a disjoint pair.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
This commit is contained in:
Claude
2026-10-15 08:35:12 +00:00
co-authored by Claude Opus 5.5
parent f7753be0e4
commit 038a759ff1
4 changed files with 70 additions and 51 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/CelerisLab/body/jets.py
+51
View File
@@ -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
src/CelerisLab/lbm/curved_links.py
+2 -48
View File
@@ -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."""
src/CelerisLab/lbm/refinement.py
+10 -3
View File
@@ -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,
tests/unit/test_jets.py
+7
View File
@@ -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])