Add body-independent Jet: true-hit spatial selection, fields, profiles

Jet(region, field, profile, amplitude) with Sector/Tube regions,
Radial/Parallel fields and uniform/parabolic/quartic profiles. Selection
uses true FP64 hit points over all scene cut links, no body association.
Overlap between jets or with a SurfaceActuation aperture is rejected.
CurvedLinkSoA caches unit-amplitude shapes; set_jet_amplitude rewrites
and re-uploads only the surface-velocity column. SurfaceActuation.as_jet
is the legacy adapter (equivalence measured in tests). Kernel and SoA
geometry columns unchanged. Manager/Simulation/fine-path wiring pending.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

src/CelerisLab/body/__init__.py

@@ -9,6 +9,7 @@ from .sync_plan import BodySyncPlan
 from .geometry.base import CutLink, SensorCell, Geometry
 from .geometry.circle import CircleGeometry
 from .coupling.soa_packer import pack_cut_links_to_soa, pack_sensor_to_soa
+from .jets import Jet, Sector, Tube, Radial, Parallel
 
 __all__ = [
     "SimObject", "ObjectState", "ObjectControl", "SurfaceActuation",
@@ -18,6 +19,7 @@ __all__ = [
     "CutLink", "SensorCell", "Geometry",
     "CircleGeometry",
     "pack_cut_links_to_soa", "pack_sensor_to_soa",
+    "Jet", "Sector", "Tube", "Radial", "Parallel",
 ]
 
 

src/CelerisLab/body/jets.py

@@ -0,0 +1,151 @@
+# CelerisLab/body/jets.py
+"""Body-independent jets: spatial region + direction field + profile + amplitude.
+
+A Jet selects cut links purely by the true wall hit point (FP64, before SoA
+packing), over every link in the scene; it carries no body association.
+Per-link velocity = amplitude * profile(s) * field(hit), where ``s`` in [-1, 1]
+is the region's lateral coordinate. Velocities are lattice units of the level
+whose links are evaluated; callers convert frames/units before evaluation.
+One physical link (receiver, direction) may be driven by at most one jet.
+"""
+
+from __future__ import annotations
+
+import math
+from dataclasses import dataclass
+
+import numpy as np
+
+PROFILES = {
+    "uniform": lambda s: np.ones_like(s),
+    "parabolic": lambda s: 1 - s*s,
+    "quartic": lambda s: (1 - s*s)**2,
+}
+
+
+def _finite(*values):
+    if not all(math.isfinite(v) for v in values):
+        raise ValueError("Jet parameters must be finite")
+
+
+@dataclass(frozen=True)
+class Sector:
+    """Annular sector about ``origin``: |angle offset| < half_width, r_min <= r <= r_max."""
+    origin: tuple
+    angle: float
+    half_width: float
+    r_max: float
+    r_min: float = 0.0
+
+    def __post_init__(self):
+        _finite(*self.origin, self.angle, self.half_width, self.r_max, self.r_min)
+        if not (0 < self.half_width <= math.pi and 0 <= self.r_min <= self.r_max):
+            raise ValueError("Sector requires 0 < half_width <= pi and 0 <= r_min <= r_max")
+
+    def select(self, hits):
+        d = hits - np.asarray(self.origin, np.float64)
+        r = np.hypot(d[:, 0], d[:, 1])
+        delta = np.angle(np.exp(1j*(np.arctan2(d[:, 1], d[:, 0]) - self.angle)))
+        return (np.abs(delta) < self.half_width) & (r >= self.r_min) & (r <= self.r_max), \
+            delta/self.half_width
+
+
+@dataclass(frozen=True)
+class Tube:
+    """Rectangle from ``origin`` along ``axis_angle`` for ``length``, |lateral| < half_width."""
+    origin: tuple
+    axis_angle: float
+    length: float
+    half_width: float
+
+    def __post_init__(self):
+        _finite(*self.origin, self.axis_angle, self.length, self.half_width)
+        if self.length <= 0 or self.half_width <= 0:
+            raise ValueError("Tube requires positive length and half_width")
+
+    def select(self, hits):
+        d = hits - np.asarray(self.origin, np.float64)
+        c, s = math.cos(self.axis_angle), math.sin(self.axis_angle)
+        along, lateral = d[:, 0]*c + d[:, 1]*s, -d[:, 0]*s + d[:, 1]*c
+        return (along >= 0) & (along <= self.length) & (np.abs(lateral) < self.half_width), \
+            lateral/self.half_width
+
+
+@dataclass(frozen=True)
+class Radial:
+    """Unit vector along hit - origin (outward) or its reverse (inward)."""
+    origin: tuple
+    inward: bool = True
+
+    def __post_init__(self):
+        _finite(*self.origin)
+
+    def direction(self, hits):
+        d = hits - np.asarray(self.origin, np.float64)
+        r = np.hypot(d[:, 0], d[:, 1])[:, None]
+        if np.any(r <= 0):
+            raise ValueError("Radial field undefined at its origin")
+        return (-d if self.inward else d)/r
+
+
+@dataclass(frozen=True)
+class Parallel:
+    """Uniform unit vector at ``angle`` (radians, counterclockwise from +x)."""
+    angle: float
+
+    def __post_init__(self):
+        _finite(self.angle)
+
+    def direction(self, hits):
+        return np.broadcast_to([math.cos(self.angle), math.sin(self.angle)], hits.shape)
+
+
+@dataclass
+class Jet:
+    """Region (Sector|Tube), field (Radial|Parallel), profile name, signed amplitude.
+
+    Only ``amplitude`` may change after registration; it never rebuilds links.
+    """
+    region: object
+    field: object
+    profile: str = "uniform"
+    amplitude: float = 0.0
+
+    def __post_init__(self):
+        if self.profile not in PROFILES:
+            raise ValueError(f"Jet profile must be one of {sorted(PROFILES)}")
+        _finite(self.amplitude)
+        if abs(self.amplitude) >= 1:
+            raise ValueError("Jet amplitude must satisfy abs(amplitude) < 1 (lattice units)")
+
+
+def jet_link_shapes(jets, hits):
+    """Unit-amplitude per-link vectors (n, 2) FP64 and jet index (n,) int32 (-1 = none).
+
+    ``hits`` are true FP64 wall hit points (n, 2) row-aligned with the packed
+    links. Raises if a jet selects nothing or two jets select the same link.
+    """
+    hits = np.asarray(hits, np.float64).reshape(-1, 2)
+    shape = np.zeros_like(hits)
+    jet_id = np.full(len(hits), -1, np.int32)
+    for j, jet in enumerate(jets):
+        selected, s = jet.region.select(hits)
+        if not selected.any():
+            raise ValueError(f"Jet {j} selects no resolved cut links")
+        if np.any(jet_id[selected] >= 0):
+            raise ValueError(f"Jet {j} overlaps links already driven by another jet")
+        jet_id[selected] = j
+        shape[selected] = (PROFILES[jet.profile](s[selected])[:, None]
+                           * jet.field.direction(hits[selected]))
+    return shape, jet_id
+
+
+def jet_link_velocity(jets, shape, jet_id):
+    """FP32 (n, 2) velocity = amplitude[jet_id] * shape; O(n), no geometry work."""
+    amplitude = np.array([jet.amplitude for jet in jets] + [0.0], np.float64)
+    return (amplitude[jet_id][:, None]*shape).astype(np.float32)
+
+
+def fine_hits_to_physical(hits_local, fine_origin, dx):
+    """Fine-local cell-index hits -> physical frame: origin + dx*(hit + 1/2)."""
+    return np.asarray(fine_origin, np.float64) + dx*(np.asarray(hits_local, np.float64) + .5)

src/CelerisLab/body/objects.py

@@ -84,6 +84,18 @@ class SurfaceActuation:
         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
 
+    def as_jet(self, center, radius, tol=1e-6):
+        """Legacy adapter: equivalent body-independent Jet on a circle of ``radius``.
+
+        The thin annulus replaces the old body filter; equality with ``project``
+        holds up to FP32 lever-arm rounding (measured in tests, not assumed).
+        """
+        from .jets import Jet, Radial, Sector
+        band = tol*max(1.0, radius)
+        return Jet(Sector(tuple(center), self.angle, self.half_width, radius + band,
+                          max(0.0, radius - band)),
+                   Radial(tuple(center), inward=False), "quartic", self.normal_velocity)
+
 
 # ---------------------------------------------------------------------------
 # SimObject container

src/CelerisLab/lbm/curved_links.py

@@ -42,14 +42,17 @@ def _upload_columns(columns, stream=None):
         stream.synchronize()
 
 
-def surface_link_velocity(objects, cfg, flags, fluid_idx, dir, body_id, rx, ry):
+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 no object carries surface actuation.
+    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:
+    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):
@@ -66,6 +69,16 @@ def surface_link_velocity(objects, cfg, flags, fluid_idx, dir, body_id, rx, ry):
             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")
+        shape, jet_id = jet_link_shapes(jets, hits)
+        driven = jet_id >= 0
+        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
 
 
@@ -107,12 +120,36 @@ class CurvedLinkSoA:
     aperture: Optional[np.ndarray] = None
     surface_velocity_gpu: Optional[cuda.DeviceAllocation] = None
     link_exchange_gpu: Optional[cuda.DeviceAllocation] = None
+    # Host-only jet cache (never uploaded or checkpointed).
+    jets: tuple = ()
+    jet_shape: Optional[np.ndarray] = None
+    jet_id: Optional[np.ndarray] = None
 
-    def assign_surface(self, objects, cfg, flags):
-        """Keep the original eight-column geometry/checkpoint contract unchanged."""
+    def assign_surface(self, objects, cfg, flags, jets=(), hits=None):
+        """Keep the original eight-column geometry/checkpoint contract unchanged.
+
+        Jet shapes are cached host-side so amplitude changes skip geometry work.
+        """
         self.surface_velocity = self.aperture = None
+        self.jets, self.jet_shape, self.jet_id = tuple(jets), None, None
         self.surface_velocity, self.aperture = surface_link_velocity(
-            objects, cfg, flags, self.fluid_idx, self.dir, self.body_id, self.rx, self.ry)
+            objects, cfg, flags, self.fluid_idx, self.dir, self.body_id, self.rx, self.ry,
+            jets=self.jets, hits=hits)
+        if self.jets:
+            from ..body.jets import jet_link_shapes
+            self.jet_shape, self.jet_id = jet_link_shapes(self.jets, hits)
+
+    def set_jet_amplitude(self, index, amplitude, stream=None):
+        """Rewrite jet rows from cached shapes; re-upload only the velocity column."""
+        from ..body.jets import Jet, jet_link_velocity
+        jet = self.jets[index]
+        Jet(jet.region, jet.field, jet.profile, amplitude)  # validate before mutation
+        jet.amplitude = float(amplitude)
+        driven = self.jet_id >= 0
+        self.surface_velocity[driven] = jet_link_velocity(
+            self.jets, self.jet_shape, self.jet_id)[driven]
+        if self.surface_velocity_gpu is not None:
+            _upload_columns([(self.surface_velocity_gpu, self.surface_velocity)], stream)
 
     def surface_exchange(self, body_id):
         """Latest boundary step's (Fx,Fy,Tz), split by aperture; not nozzle thrust.
@@ -198,6 +235,7 @@ class CurvedLinkSoA:
                     "CurvedLinkSoA column 'fallback_class' must use classes {0, 1}"
                 )
         self.surface_velocity = self.aperture = None
+        self.jets, self.jet_shape, self.jet_id = (), None, None
 
     def upload(self, stream: Optional[cuda.Stream] = None) -> None:
         """Upload host columns to GPU with capacity-based allocation reuse.

tests/unit/test_jets.py

@@ -0,0 +1,134 @@
+"""Body-independent Jet: true-hit spatial selection, fields, profiles, overlap, amplitude."""
+import math
+
+import numpy as np
+import pytest
+
+from CelerisLab.body.jets import (Jet, Parallel, Radial, Sector, Tube,
+                                  fine_hits_to_physical, jet_link_shapes)
+from CelerisLab.body.objects import SurfaceActuation
+from CelerisLab.lbm.curved_links import CurvedLinkSoA, surface_link_velocity
+from .test_surface_actuation import circle, projected
+
+CENTER, RADIUS = (6., 5.), 2.999
+
+
+def scene(spec=None):
+    """Single circle links plus their true FP64 hits, row-aligned with the SoA."""
+    obj = circle(spec)
+    cl, cfg, flags = projected(obj)
+    links = obj.geometry.build_cut_links(cfg.nx, cfg.ny, domain_flags=flags)
+    hits = np.array([(link.hit_x, link.hit_y) for link in links])
+    assert len(hits) == len(cl.dir)
+    return obj, cl, cfg, flags, hits
+
+
+def arc(hits, angle, half_width):
+    d = hits - CENTER
+    delta = np.angle(np.exp(1j*(np.arctan2(d[:, 1], d[:, 0]) - angle)))
+    return np.abs(delta) < half_width
+
+
+def test_concentric_sector_gives_inward_radial_and_tube_parallel_on_same_links():
+    *_, hits = scene()
+    half = math.radians(10)
+    expected = arc(hits, 0., half)
+    assert expected.sum() >= 2
+    sector = Jet(Sector(CENTER, 0., half, RADIUS + .1), Radial(CENTER), "uniform", 1e-3)
+    shape, jet_id = jet_link_shapes([sector], hits)
+    assert np.array_equal(jet_id >= 0, expected)
+    d = hits[expected] - CENTER
+    np.testing.assert_allclose(shape[expected], -d/np.hypot(*d.T)[:, None], atol=1e-15)
+    # Tube along +x through the same arc: lateral |y-5| < R sin(10 deg).
+    tube = Jet(Tube((CENTER[0], CENTER[1]), 0., RADIUS + .1, RADIUS*math.sin(half)),
+               Parallel(math.pi), "uniform", 1e-3)
+    shape_t, id_t = jet_link_shapes([tube], hits)
+    assert np.array_equal(id_t >= 0, expected)
+    np.testing.assert_allclose(shape_t[expected], [[-1., 0.]]*expected.sum(), atol=1e-15)
+
+
+def test_off_centre_radial_points_to_field_origin_not_body_centre():
+    *_, hits = scene()
+    origin = (12., 9.)
+    jet = Jet(Sector(CENTER, .1, math.pi, 10.), Radial(origin, inward=False))
+    shape, jet_id = jet_link_shapes([jet], hits)
+    assert np.all(jet_id == 0)
+    d = hits - origin
+    np.testing.assert_allclose(shape, d/np.hypot(*d.T)[:, None], atol=1e-15)
+    assert not np.allclose(shape, (hits - CENTER)/RADIUS, atol=1e-3)
+
+
+@pytest.mark.parametrize("profile,expected", [("uniform", 1.), ("parabolic", .75), ("quartic", .5625)])
+def test_profiles_on_lateral_coordinate(profile, expected):
+    hits = np.array([[1., .5]])
+    jet = Jet(Tube((0., 0.), 0., 2., 1.), Parallel(0.), profile)
+    shape, _ = jet_link_shapes([jet], hits)
+    np.testing.assert_allclose(shape[0], [expected, 0.], rtol=1e-15)
+
+
+def test_selection_uses_true_hit_not_clipped_q():
+    # Same link as test_aperture_uses_true_hit_not_fluid_node_or_clamped_q:
+    # receiver (9,5), direction (-1,+1); true t < .01 but packed q == .01.
+    t = (6-math.sqrt(36-8*(9-2.999**2)))/4
+    angle = math.atan2(t, 3-t)
+    _, cl, _, _, hits = scene()
+    row = np.flatnonzero((cl.fluid_idx == 5*12 + 9) & np.isclose(cl.q, .01))
+    assert len(row) == 1
+    clipped = np.array([9., 5.]) + .01*np.array([-1., 1.])
+    jet = Jet(Sector(CENTER, angle, 1e-4, RADIUS + .1), Radial(CENTER))
+    _, jet_id = jet_link_shapes([jet], hits)
+    assert jet_id[row[0]] == 0
+    with pytest.raises(ValueError, match="selects no"):
+        jet_link_shapes([jet], clipped[None, :])
+
+
+def test_fine_half_cell_frame():
+    np.testing.assert_allclose(fine_hits_to_physical([[0., 0.], [3., 1.5]], (4., 6.), .5),
+                               [[4.25, 6.25], [5.75, 7.]], rtol=0, atol=0)
+
+
+def test_links_are_distinct_by_receiver_and_direction_and_overlap_rejected():
+    # Two physical links (different receiver/direction) sharing one hit point
+    # are two rows, both driven; two jets on one link are rejected.
+    hits = np.array([[1., 1.], [1., 1.]])
+    a = Jet(Tube((0., 1.), 0., 2., .5), Parallel(0.))
+    shape, jet_id = jet_link_shapes([a], hits)
+    assert jet_id.tolist() == [0, 0]
+    b = Jet(Sector((1., 0.), math.pi/2, .1, 2.), Parallel(0.))
+    with pytest.raises(ValueError, match="overlaps"):
+        jet_link_shapes([a, b], hits)
+    obj, cl, cfg, flags, real = scene(SurfaceActuation(0., .3, .002))
+    with pytest.raises(ValueError, match="SurfaceActuation"):
+        surface_link_velocity([obj], cfg, flags, cl.fluid_idx, cl.dir, cl.body_id, cl.rx, cl.ry,
+                              jets=[Jet(Sector(CENTER, 0., .1, 4.), Radial(CENTER))], hits=real)
+    with pytest.raises(ValueError, match="row-aligned"):
+        surface_link_velocity([], cfg, flags, cl.fluid_idx, cl.dir, cl.body_id, cl.rx, cl.ry,
+                              jets=[a], hits=real[:-1])
+
+
+def test_amplitude_update_reuses_cached_shapes():
+    obj, cl, cfg, flags, hits = scene()
+    jet = Jet(Sector(CENTER, 0., .5, 4.), Radial(CENTER), "parabolic", 1e-3)
+    cl.assign_surface([obj], cfg, flags, jets=[jet], hits=hits)
+    shape, ids = cl.jet_shape, cl.jet_id
+    before = cl.surface_velocity.copy()
+    cl.set_jet_amplitude(0, -2e-3)
+    assert cl.jet_shape is shape and cl.jet_id is ids
+    np.testing.assert_allclose(cl.surface_velocity, -2*before, rtol=1e-6, atol=0)
+    assert np.array_equal(cl.aperture, ids >= 0)
+    with pytest.raises(ValueError):
+        cl.set_jet_amplitude(0, 1.5)
+    assert jet.amplitude == -2e-3
+
+
+@pytest.mark.parametrize("spec", [SurfaceActuation(0., .6, .002), SurfaceActuation(2., 1.1, -.003)])
+def test_legacy_adapter_measured_equivalence(spec):
+    obj, legacy, cfg, flags, hits = scene(spec)
+    cl = CurvedLinkSoA()
+    cl.assign_host(legacy.fluid_idx, legacy.dir, legacy.q, legacy.body_id,
+                   legacy.rx, legacy.ry, legacy.rz, legacy.fallback_class)
+    cl.assign_surface([circle()], cfg, flags, jets=[spec.as_jet(CENTER, RADIUS)], hits=hits)
+    assert np.array_equal(cl.aperture, legacy.aperture)
+    diff = np.abs(cl.surface_velocity - legacy.surface_velocity).max()
+    print(f"legacy adapter max |dU| = {diff:.3e} (exact={diff == 0})")
+    assert diff <= 4*np.finfo(np.float32).eps*abs(spec.normal_velocity)