ClaudeandClaude Opus 5.5
62c525e735
Sector/Tube/Radial store origin as a validated float 2-tuple, so caller lists/arrays cannot change a registered definition. Simulation.set_body rejects nonzero motion for bodies with Jet-driven links once initialized (before initialize the sync-time stationary check applies). prepare_fine_boundary accepts jets in the same coarse cell-edge frame as fine objects, records read-only FP64 edge-frame hits (FineBoundary2D .curved_hits, default None) and rejects empty or overlapping selections before CUDA. FixedPatchStepper wiring is not included. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
273 lines
13 KiB
Python
273 lines
13 KiB
Python
"""Body-independent Jet: true-hit spatial selection, fields, profiles, overlap, amplitude."""
|
|
import math
|
|
from types import SimpleNamespace
|
|
|
|
import numpy as np
|
|
import pytest
|
|
|
|
from CelerisLab.body.jets import (Jet, Parallel, Radial, Sector, Tube,
|
|
edge_to_node_frame, fine_hits_to_edge_frame,
|
|
jet_link_shapes)
|
|
from CelerisLab.body.objects import SurfaceActuation
|
|
from CelerisLab.lbm.curved_links import CurvedLinkSoA, surface_link_velocity
|
|
from .test_state_regressions import device, make_field, manager
|
|
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) & (cl.dir == 8))
|
|
assert len(row) == 1 and np.isclose(cl.q[row[0]], .01)
|
|
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():
|
|
edge = fine_hits_to_edge_frame([[0., 0.], [3., 1.5]], (4., 6.), .5)
|
|
np.testing.assert_array_equal(edge, [[4.25, 6.25], [5.75, 7.]])
|
|
np.testing.assert_array_equal(edge_to_node_frame(edge), [[3.75, 5.75], [5.25, 6.5]])
|
|
|
|
|
|
def test_same_physical_hit_in_simulation_and_patch_frames():
|
|
# One physical circle: Simulation node-frame centre C; patch edge-frame
|
|
# centre C + 1/2; fine-local centre (C + 1/2 - origin)/dx - 1/2 as in
|
|
# prepare_fine_boundary. Fine hits mapped edge -> node must lie on the
|
|
# Simulation circle, and one node-frame Jet selects the same arc on both.
|
|
from CelerisLab.body.geometry.circle import CircleGeometry
|
|
C, R, origin, dx = np.array([10.3, 9.6]), 2.15, np.array([4., 3.]), .5
|
|
local = (C + .5 - origin)/dx - .5
|
|
fine = CircleGeometry(*local, R/dx).build_cut_links(40, 40)
|
|
coarse = CircleGeometry(*C, R).build_cut_links(24, 20)
|
|
fine_hits = edge_to_node_frame(fine_hits_to_edge_frame(
|
|
[(l.hit_x, l.hit_y) for l in fine], origin, dx))
|
|
coarse_hits = np.array([(l.hit_x, l.hit_y) for l in coarse])
|
|
np.testing.assert_allclose(np.hypot(*(fine_hits - C).T), R, rtol=0, atol=1e-12)
|
|
np.testing.assert_allclose(np.hypot(*(coarse_hits - C).T), R, rtol=0, atol=1e-12)
|
|
jet = Jet(Sector(tuple(C), .4, .5, R + .01, R - .01), Radial(tuple(C), inward=False))
|
|
for hits in (fine_hits, coarse_hits):
|
|
shape, jet_id = jet_link_shapes([jet], hits)
|
|
d = hits - C
|
|
angle = np.arctan2(d[:, 1], d[:, 0])
|
|
assert np.array_equal(jet_id == 0, np.abs(angle - .4) < .5)
|
|
np.testing.assert_allclose(shape[jet_id == 0], d[jet_id == 0]/R, atol=1e-12)
|
|
# Relabelling instead of converting moves every hit off the circle.
|
|
wrong = fine_hits_to_edge_frame([(l.hit_x, l.hit_y) for l in fine], origin, dx)
|
|
assert np.abs(np.hypot(*(wrong - C).T) - R).max() > .3
|
|
|
|
|
|
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 cl.jets[0].amplitude == -2e-3 and jet.amplitude == 1e-3 # registered Jet untouched
|
|
np.testing.assert_allclose(cl.surface_velocity, -2*before, rtol=1e-6, atol=0)
|
|
|
|
|
|
def test_jet_definition_is_immutable_and_amplitude_validated():
|
|
jet = Jet(Sector(CENTER, 0., .5, 4.), Radial(CENTER), "parabolic", 1e-3)
|
|
for name, value in (("field", Parallel(0.)), ("region", Tube(CENTER, 0., 1., 1.)),
|
|
("profile", "uniform"), ("amplitude", float("nan"))):
|
|
with pytest.raises(AttributeError):
|
|
setattr(jet, name, value)
|
|
for bad in (float("nan"), float("inf"), 1.0):
|
|
with pytest.raises(ValueError):
|
|
jet.with_amplitude(bad)
|
|
assert jet.with_amplitude(-1e-3).field is jet.field and jet.amplitude == 1e-3
|
|
with pytest.raises(ValueError, match="finite"):
|
|
jet_link_shapes([jet], [[np.nan, 5.]])
|
|
|
|
|
|
@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)
|
|
|
|
|
|
def jet_case(device, jets):
|
|
field = make_field(device)
|
|
bodies = manager(field, [circle(radius=2.)])
|
|
bodies.jets.extend(jets)
|
|
bodies.sync_to_gpu(field)
|
|
return field, bodies
|
|
|
|
|
|
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)
|
|
field, bodies = jet_case(device, [jet])
|
|
cl = field.curved
|
|
assert cl.jet_id is not None and (cl.jet_id == 0).sum() >= 2
|
|
velocity_gpu = cl.surface_velocity_gpu
|
|
np.testing.assert_array_equal(np.frombuffer(velocity_gpu.data, np.float32).reshape(-1, 2),
|
|
cl.surface_velocity)
|
|
geometry = {name: bytes(getattr(cl, name + "_gpu").data)
|
|
for name in ("fluid_idx", "dir", "q", "body_id", "rx", "ry")}
|
|
writes, before = device.writes, cl.surface_velocity.copy()
|
|
cl.set_jet_amplitude(0, -1e-3, SimpleNamespace(synchronize=lambda: None))
|
|
assert device.writes == writes + 1 and cl.surface_velocity_gpu is velocity_gpu
|
|
np.testing.assert_array_equal(np.frombuffer(velocity_gpu.data, np.float32).reshape(-1, 2),
|
|
cl.surface_velocity)
|
|
assert all(bytes(getattr(cl, name + "_gpu").data) == data for name, data in geometry.items())
|
|
np.testing.assert_allclose(cl.surface_velocity, -.5*before, rtol=1e-6, atol=0)
|
|
np.testing.assert_array_equal(cl.surface_velocity[cl.jet_id < 0], 0.)
|
|
|
|
|
|
def test_manager_rejects_overlapping_jets_and_moving_bodies(device):
|
|
a = Jet(Sector(CENTER, 0., .6, 3.), Radial(CENTER))
|
|
b = Jet(Tube(CENTER, 0., 3., .5), Parallel(0.))
|
|
with pytest.raises(ValueError, match="overlaps"):
|
|
jet_case(device, [a, b])
|
|
field = make_field(device)
|
|
moving = circle(radius=2.)
|
|
moving.state.omega = 1e-3
|
|
bodies = manager(field, [moving])
|
|
bodies.jets.append(a)
|
|
with pytest.raises(ValueError, match="stationary"):
|
|
bodies.sync_to_gpu(field)
|
|
|
|
|
|
def test_origins_are_frozen_float_pairs():
|
|
point = [6., 5.]
|
|
array = np.array([6., 5.])
|
|
jets = [Jet(Sector(point, 0., .5, 4.), Radial(array)), Jet(Tube(array, 0., 1., 1.), Parallel(0.))]
|
|
point[0], array[0] = 100., 100.
|
|
assert jets[0].region.origin == (6., 5.) and jets[0].field.origin == (6., 5.)
|
|
assert jets[1].region.origin == (6., 5.)
|
|
assert all(type(v) is float for v in jets[0].region.origin)
|
|
for bad in ((1., 2., 3.), (1., float("nan"))):
|
|
with pytest.raises(ValueError):
|
|
Radial(bad)
|
|
|
|
|
|
def test_set_body_rejects_motion_of_jet_driven_body(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, lbm_cfg=field.cfg, _initialized=True)
|
|
with pytest.raises(ValueError, match="stationary"):
|
|
Simulation.set_body(sim, 0, omega=.01)
|
|
assert bodies.get(0).state.omega == 0
|
|
Simulation.set_body(sim, 0, omega=0.) # stationary updates stay allowed
|
|
|
|
|
|
def test_fine_boundary_keeps_edge_frame_hits_and_checks_jets():
|
|
from CelerisLab.body.objects import SimObject
|
|
from CelerisLab.body.geometry.circle import CircleGeometry
|
|
from CelerisLab.config import LBMConfig
|
|
from CelerisLab.lbm import refinement
|
|
cfg = LBMConfig(nx=24, ny=24, viscosity=.05)
|
|
cfg.validate()
|
|
patch = refinement.FixedPatch2D.from_config(cfg, (4, 4, 20, 20))
|
|
C, R = (12.25, 12.25), 2.15 # coarse cell-edge frame, as fine_objects
|
|
body = SimObject(0, geometry=CircleGeometry(*C, R), center=C, radius=R)
|
|
plan = refinement.prepare_fine_boundary(cfg, patch, [body])
|
|
assert plan.curved_hits.shape == (len(plan.curved_host[0]), 2)
|
|
assert not plan.curved_hits.flags.writeable
|
|
np.testing.assert_allclose(np.hypot(*(plan.curved_hits - C).T), R, rtol=0, atol=1e-12)
|
|
jet = Jet(Sector(C, 0., .5, R + .01), Radial(C, inward=False))
|
|
refinement.prepare_fine_boundary(cfg, patch, [body], jets=[jet])
|
|
away = Jet(Sector((2., 2.), 0., .1, .5), Parallel(0.))
|
|
with pytest.raises(ValueError, match="selects no"):
|
|
refinement.prepare_fine_boundary(cfg, patch, [body], jets=[away])
|
|
with pytest.raises(ValueError, match="overlaps"):
|
|
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])
|