diff --git a/src/CelerisLab/body/jets.py b/src/CelerisLab/body/jets.py index cda80b0..e75ecb7 100644 --- a/src/CelerisLab/body/jets.py +++ b/src/CelerisLab/body/jets.py @@ -216,11 +216,8 @@ def surface_link_velocity(objects, cfg, flags, fluid_idx, dir, body_id, rx, ry, 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") + from .preprocess.flag_overlay import validate_cut_link_ownership + validate_cut_link_ownership(flags, fluid_idx, dir, cfg.nq) surface_velocity = np.zeros((len(dir), 2), np.float32) surface_aperture = np.zeros(len(dir), bool) for obj in enabled: diff --git a/src/CelerisLab/body/manager.py b/src/CelerisLab/body/manager.py index feafb1f..819c42b 100644 --- a/src/CelerisLab/body/manager.py +++ b/src/CelerisLab/body/manager.py @@ -27,7 +27,7 @@ import pycuda.driver as cuda from .objects import SimObject from .registry import BodyRegistry -from .preprocess.flag_overlay import merge_flag_masks +from .preprocess.flag_overlay import merge_flag_masks, validate_cut_link_ownership from ..cuda.compiler_v2 import obs_layout from ..cuda.obs_layout import ObsLayout from ..lbm.descriptors import ( @@ -369,9 +369,9 @@ class ObjectManager: Args: objects: List of SimObject instances to process. - domain_flags: Optional uint16 flags for donor-cell classification. + domain_flags: Optional uint16 flags for receiver/donor classification. sensor_cell_counts_out: Int array to fill with per-object sensor - cell counts. If provided, ``fill(0)`` is called first. + cell counts after the full geometry has passed validation. hits_out: List that receives one FP64 (n, 2) array of true hit points row-aligned with the cut-link columns (for Jets). @@ -394,8 +394,8 @@ class ObjectManager: 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) + sensor_counts = (np.zeros_like(sensor_cell_counts_out) + if sensor_cell_counts_out is not None else None) ez = None if self.cfg and self.cfg.is_d3q19: @@ -430,9 +430,9 @@ class ObjectManager: if len(cells) > 0: s_cells.append(cells) s_ids.append(ids) - if (sensor_cell_counts_out is not None - and 0 <= obj.obj_id < sensor_cell_counts_out.size): - sensor_cell_counts_out[obj.obj_id] = int(len(cells)) + if (sensor_counts is not None + and 0 <= obj.obj_id < sensor_counts.size): + sensor_counts[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] @@ -461,8 +461,9 @@ class ObjectManager: r_fr_ids = np.concatenate(fr_ids) if fr_ids else np.zeros(0, dtype=np.int32) r_hit = np.concatenate(cl_hit) if cl_hit else np.zeros((0, 2)) - if hits_out is not None: - hits_out.append(r_hit) + # Passive walls need the same ownership guarantees as driven walls. + # Overlapping solids are not a supported union-boundary operation. + validate_cut_link_ownership(domain_flags, r_idx, r_dir, self.nq) if self.jets or 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 @@ -470,6 +471,11 @@ class ObjectManager: 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, jets=self.jets, hits=r_hit) + if sensor_cell_counts_out is not None: + sensor_cell_counts_out[:] = sensor_counts + if hits_out is not None: + hits_out.append(r_hit) + return ( r_idx, r_dir, r_q, r_bid, r_rx, r_ry, r_rz, r_fallback, @@ -487,8 +493,6 @@ class ObjectManager: the caller has already prepared and uploaded the final flag field for this geometry, such as during initialization. """ - self._telemetry_field = field - next_flags = self.build_flags(field.build_channel_flags()) if rebuild_flags else field.flag hits = [] @@ -499,6 +503,7 @@ class ObjectManager: fr_cells, fr_obj_id, ) = self.build_compact_lists(domain_flags=next_flags, hits_out=hits) + self._telemetry_field = field if rebuild_flags: field.flag = next_flags field.upload_flags() diff --git a/src/CelerisLab/body/preprocess/flag_overlay.py b/src/CelerisLab/body/preprocess/flag_overlay.py index 5951190..a3e0600 100644 --- a/src/CelerisLab/body/preprocess/flag_overlay.py +++ b/src/CelerisLab/body/preprocess/flag_overlay.py @@ -1,6 +1,6 @@ # CelerisLab/body/preprocess/flag_overlay.py """ -Flag overlay: merge per-object flag masks onto a clean domain base. +Flag overlay and cut-link ownership checks against the merged domain. """ from __future__ import annotations @@ -12,6 +12,15 @@ import numpy as np from ...lbm.descriptors import FLUID, SOLID, SENSOR_FLAG, FRC_REGION +def validate_cut_link_ownership(flags, fluid_idx, directions, nq): + """Reject solid receivers and multiple writers to one population.""" + if np.any((flags[fluid_idx] & FLUID) == 0) or np.any(flags[fluid_idx] & SOLID): + raise ValueError("Cut links must have active fluid receivers") + writers = fluid_idx.astype(np.uint64) * nq + directions + if len(np.unique(writers)) != len(directions): + raise ValueError("Boundary requires unique cut-link writers") + + def merge_flag_masks( base: np.ndarray, masks: List[np.ndarray], diff --git a/src/CelerisLab/simulation.py b/src/CelerisLab/simulation.py index e141131..b4f60f4 100644 --- a/src/CelerisLab/simulation.py +++ b/src/CelerisLab/simulation.py @@ -336,13 +336,19 @@ class Simulation: def read_surface_exchange(self, id: int): """Latest step's wall/aperture (Fx,Fy,Tz) exchange, not total nozzle force. + Supports SurfaceActuation and Jet-driven bodies, including zero-amplitude + apertures. Passive bodies without an aperture are rejected. 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: + obj = self.bodies.get(id) + curved = self.field.curved + has_aperture = (curved.aperture is not None + and np.any(curved.aperture & (curved.body_id == id))) + if obj.surface_actuation is None and not has_aperture: raise ValueError("Body has no surface actuation") self.stream.synchronize() - return self.field.curved.surface_exchange(id) + return curved.surface_exchange(id) @context_bound def read_body(self, id: int, *, @@ -451,18 +457,21 @@ class Simulation: @context_bound def initialize(self): """Initialize flow field and sync objects to GPU.""" - # Recompile if objects were added after construction. self._assert_allocation_contract() + # Validate the complete geometry before flags, populations or modules + # change. Compact-list construction publishes no outputs on failure. + next_flags = self.bodies.build_flags(self.field.build_channel_flags()) + if self.bodies.jets and self.bodies.count == 0: + raise ValueError("Jets need bodies with cut links") + self.bodies.build_compact_lists_for(self.bodies.objects, next_flags) + + # Recompile if objects were added after construction. if self._read_generated_n_objects() != self.bodies.count: self.recompile() - self.field.flag = self.field.build_channel_flags() - if self.bodies.count > 0: - self.field.flag = self.bodies.build_flags(self.field.flag) + self.field.flag = next_flags self.field.upload_flags() - if self.bodies.jets and self.bodies.count == 0: - raise ValueError("Jets need bodies with cut links") self.stepper.initialize() if self.bodies.count > 0: self.bodies.sync_to_gpu(self.field, rebuild_flags=False) diff --git a/tests/unit/test_jets.py b/tests/unit/test_jets.py index 9285f2e..947a70a 100644 --- a/tests/unit/test_jets.py +++ b/tests/unit/test_jets.py @@ -188,6 +188,55 @@ def jet_case(device, jets): return field, bodies +def test_public_surface_exchange_for_jet_preserves_aperture_at_zero_and_reverse(device): + from CelerisLab.simulation import Simulation + jet = Jet(Sector(CENTER, 0., .6, 3.), Radial(CENTER, inward=False), "quartic", 2e-3) + field, bodies = jet_case(device, [jet]) + sim = SimpleNamespace(bodies=bodies, field=field, _initialized=True, + stream=SimpleNamespace(synchronize=device.synchronize)) + cl = field.curved + assert bodies.get(0).surface_actuation is None + aperture, velocity = cl.aperture, cl.surface_velocity.copy() + assert aperture.any() and (~aperture).any() + # Synthetic per-link exchanges test readback/splitting, not kernel physics. + values = np.arange(cl.count * 3, dtype=np.float32).reshape(-1, 3) + device.htod(cl.link_exchange_gpu, values) + for amplitude in (2e-3, 0., -1e-3): + Simulation.set_jet_amplitude(sim, 0, amplitude) + result = Simulation.read_surface_exchange(sim, 0) + assert cl.aperture is aperture + np.testing.assert_array_equal(aperture, cl.jet_id >= 0) + np.testing.assert_allclose(cl.surface_velocity, velocity * (amplitude / 2e-3)) + for name, mask in (("wall", ~aperture), ("aperture", aperture)): + np.testing.assert_array_equal(result[name], values[mask].sum(axis=0, dtype=np.float64)) + + +def test_public_surface_exchange_keeps_legacy_and_passive_guards(device): + from CelerisLab.simulation import Simulation + field = make_field(device) + bodies = manager(field, [circle(SurfaceActuation(0., .6, .002), radius=2.)]) + bodies.sync_to_gpu(field) + sim = SimpleNamespace(bodies=bodies, field=field, + stream=SimpleNamespace(synchronize=device.synchronize)) + for value in Simulation.read_surface_exchange(sim, 0).values(): + np.testing.assert_array_equal(value, 0.) + field.curved.free_gpu_columns() + with pytest.raises(RuntimeError, match="diagnostic is allocated"): + Simulation.read_surface_exchange(sim, 0) + + # An exchange buffer elsewhere in the scene must not admit a passive body. + for jets in ([], [Jet(Sector(CENTER, 0., .6, 3.), Radial(CENTER))]): + field, bodies = jet_case(device, jets) + passive = bodies.add(circle(center=(2., 3.), radius=.6)) + bodies.sync_to_gpu(field) + sim.bodies, sim.field = bodies, field + with pytest.raises(ValueError, match="Body has no surface actuation"): + Simulation.read_surface_exchange(sim, passive) + if not jets: + with pytest.raises(ValueError, match="Body has no surface actuation"): + Simulation.read_surface_exchange(sim, 0) + + def test_manager_passes_true_hits_and_reuploads_only_velocity(device): # CPU memory double only: checks host/device bookkeeping, not GPU numerics. jet = Jet(Sector(CENTER, 0., .6, 3.), Radial(CENTER, inward=False), "quartic", 2e-3) diff --git a/tests/unit/test_state_regressions.py b/tests/unit/test_state_regressions.py index 2d3ad7b..60a5ce8 100644 --- a/tests/unit/test_state_regressions.py +++ b/tests/unit/test_state_regressions.py @@ -238,6 +238,81 @@ def test_overlapping_force_regions_rejected(device): bodies.build_compact_lists(bodies.build_flags(field.flag)) +@pytest.mark.parametrize("second_x,error", [(14., "unique cut-link writers"), + (17., "active fluid receivers")]) +def test_passive_overlap_rejected_before_publishing_lists(second_x, error): + cfg = LBMConfig(nx=32, ny=24) + cfg.validate() + bodies = ObjectManager(cfg.nx, cfg.ny, cfg.nz, cfg.nq, cfg) + for x in (14., second_x): + bodies.add(obj(x=x, y=12., r=4.3)) + bodies.add(obj(x=4., y=4., r=2., is_sensor=True)) + counts = np.full(bodies.count, 91, np.int32) + hits = [] + flags = bodies.build_flags(np.full(cfg.nx * cfg.ny, FLUID, np.uint16)) + for domain_flags in (flags, None): + with pytest.raises(ValueError, match=error): + bodies.build_compact_lists_for(bodies.objects, domain_flags, + sensor_cell_counts_out=counts, hits_out=hits) + np.testing.assert_array_equal(counts, 91) + assert hits == [] + + +@pytest.mark.parametrize("centers", [((10., 12.), (22., 12.)), + ((10.5, 12.), (19.5, 12.))]) +def test_separated_passive_bodies_keep_valid_writers(centers): + # The second case has a subcell physical gap with a resolved fluid column. + cfg = LBMConfig(nx=32, ny=24) + cfg.validate() + bodies = ObjectManager(cfg.nx, cfg.ny, cfg.nz, cfg.nq, cfg) + for x, y in centers: + bodies.add(obj(x=x, y=y, r=4.3)) + flags = bodies.build_flags(np.full(cfg.nx * cfg.ny, FLUID, np.uint16)) + indices, directions, _, ids, *_ = bodies.build_compact_lists(flags) + assert set(ids) == {0, 1} + assert np.all(flags[indices] & FLUID) and not np.any(flags[indices] & SOLID) + assert len(np.unique(indices.astype(np.uint64) * cfg.nq + directions)) == len(indices) + + +@pytest.mark.parametrize("entry", ["sync", "initialize", "runtime_edit"]) +@pytest.mark.parametrize("second_x", [5., 6.]) +def test_passive_overlap_failure_keeps_live_state(device, entry, second_x): + from CelerisLab.simulation import Simulation + field = make_field(device) + bodies = manager(field, [obj(), obj(x=9., y=3., r=1., is_sensor=True)]) + bodies.sync_to_gpu(field) + overlapping = obj(x=second_x) + if entry == "runtime_edit": + bodies.stage_add(overlapping) + else: + bodies.add(overlapping) + objects = bodies.objects + flags, curved, counts = field.flag.copy(), field.curved.fluid_idx.copy(), bodies.sensor_cell_counts.copy() + action = bodies.action.copy() + buffers = [(block, bytes(block.data), block.freed) for block in device.blocks] + writes = device.writes + sim = SimpleNamespace(bodies=bodies, field=field, lbm_cfg=field.cfg, + _initialized=entry == "runtime_edit", + _assert_allocation_contract=lambda: None, + stream=SimpleNamespace(synchronize=device.synchronize)) + with pytest.raises(ValueError, match="active fluid receivers|unique cut-link writers"): + if entry == "sync": + bodies.sync_to_gpu(field) + elif entry == "initialize": + Simulation.initialize(sim) + else: + Simulation.sync_bodies(sim) + assert bodies.objects == objects + assert bodies.has_pending_edit() == (entry == "runtime_edit") + assert sim._initialized == (entry == "runtime_edit") + assert device.writes == writes and len(device.blocks) == len(buffers) + assert all(bytes(block.data) == data and block.freed == freed for block, data, freed in buffers) + np.testing.assert_array_equal(field.flag, flags) + np.testing.assert_array_equal(field.curved.fluid_idx, curved) + np.testing.assert_array_equal(bodies.sensor_cell_counts, counts) + np.testing.assert_array_equal(bodies.action, action) + + def test_action_upload_and_scratch_are_pinned(device): field = make_field(device, streaming="esopull") bodies = manager(field, [obj(is_force_region=True)])