diff --git a/src/CelerisLab/body/__init__.py b/src/CelerisLab/body/__init__.py index adc81e8..7d678ab 100644 --- a/src/CelerisLab/body/__init__.py +++ b/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", ] diff --git a/src/CelerisLab/body/jets.py b/src/CelerisLab/body/jets.py new file mode 100644 index 0000000..7bb5016 --- /dev/null +++ b/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) diff --git a/src/CelerisLab/body/objects.py b/src/CelerisLab/body/objects.py index de3afbc..c3fb101 100644 --- a/src/CelerisLab/body/objects.py +++ b/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 diff --git a/src/CelerisLab/lbm/curved_links.py b/src/CelerisLab/lbm/curved_links.py index 7f15e89..3fdf0e4 100644 --- a/src/CelerisLab/lbm/curved_links.py +++ b/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. diff --git a/tests/unit/test_jets.py b/tests/unit/test_jets.py new file mode 100644 index 0000000..4ddde23 --- /dev/null +++ b/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)