Grok
103f144404
Rows 2 and ny-3 are ordinary bulk cells that explosion only reads as donors; only rows 0/1 and ny-2/ny-1 are dropped or reconstructed after coalescence. Matches the existing y0>=4, y1<=ny-4 slot_reservoir contract. Adds margin-4 acceptance and margin-3 rejection for a central circle.
234 lines
11 KiB
Python
234 lines
11 KiB
Python
"""CPU-only fine circle projection and transfer/boundary isolation contracts."""
|
|
from dataclasses import asdict, replace
|
|
|
|
import numpy as np
|
|
import pytest
|
|
|
|
from CelerisLab.body.geometry.circle import CircleGeometry
|
|
from CelerisLab.body.objects import SimObject, SurfaceActuation
|
|
from CelerisLab.config import LBMConfig
|
|
from CelerisLab.lbm import refinement
|
|
from CelerisLab.lbm.descriptors import D2Q9_EX, D2Q9_EY, FLUID, SOLID
|
|
from CelerisLab.lbm.refinement import ACTIVE, INNER_GHOST, FixedPatch2D
|
|
|
|
|
|
def setup_patch():
|
|
cfg = LBMConfig(nx=24, ny=24, viscosity=.05)
|
|
cfg.validate()
|
|
return cfg, FixedPatch2D.from_config(cfg, (4, 4, 20, 20))
|
|
|
|
|
|
def circle(center=(12.25, 12.25), radius=2.15, speed=.002):
|
|
return SimObject(0, geometry=CircleGeometry(*center, radius), center=center,
|
|
radius=radius, surface_actuation=SurfaceActuation(.7, .85, speed))
|
|
|
|
|
|
@pytest.mark.parametrize('speed', [.002, -.002])
|
|
def test_fine_projection_uses_cell_center_bridge_and_preserves_source(speed):
|
|
cfg, patch = setup_patch()
|
|
source = circle(center=(12.25, 12.25 if speed > 0 else 11.75), speed=speed)
|
|
expected_center = (20., 20. if speed > 0 else 19.)
|
|
assert patch.fine_origin == (2., 2.)
|
|
before = (source.center, source.radius, source.geometry.center,
|
|
source.geometry._radius, asdict(source.state), asdict(source.surface_actuation),
|
|
vars(cfg).copy())
|
|
plan = refinement.prepare_fine_boundary(cfg, patch, [source])
|
|
local, = plan.objects
|
|
assert local is not source and local.geometry is not source.geometry
|
|
np.testing.assert_allclose(local.center, expected_center, rtol=0, atol=1e-14)
|
|
np.testing.assert_allclose(local.geometry.center, expected_center, rtol=0, atol=1e-14)
|
|
assert local.radius == pytest.approx(4.3)
|
|
assert local.geometry._radius == pytest.approx(4.3)
|
|
assert local.surface_actuation.normal_velocity == speed # dt/dx is unchanged.
|
|
assert (source.center, source.radius, source.geometry.center,
|
|
source.geometry._radius, asdict(source.state), asdict(source.surface_actuation),
|
|
vars(cfg).copy()) == before
|
|
_, _, _, _, rx, ry, _, _ = plan.curved_host
|
|
np.testing.assert_allclose(np.hypot(rx, ry), 4.3, rtol=0, atol=6e-7)
|
|
# These world points lie on the ORIGINAL physical circle, including the +.5.
|
|
world_x = patch.fine_origin[0]+.5*(local.center[0]+rx.astype(float)+.5)
|
|
world_y = patch.fine_origin[1]+.5*(local.center[1]+ry.astype(float)+.5)
|
|
np.testing.assert_allclose(np.hypot(world_x-source.center[0], world_y-source.center[1]),
|
|
source.radius, rtol=0, atol=3e-7)
|
|
|
|
|
|
def test_fine_boundary_receivers_sources_and_reverse_coverage():
|
|
cfg, patch = setup_patch()
|
|
plan = refinement.prepare_fine_boundary(cfg, patch, [circle()])
|
|
ny, nx = patch.fine_shape
|
|
flags = plan.flags
|
|
assert flags.shape == (nx*ny,) and flags.dtype == np.uint16
|
|
solid = (flags & SOLID) != 0
|
|
assert solid.any()
|
|
assert np.all(patch.fine_kind.ravel()[solid] == ACTIVE)
|
|
receivers, directions, q, _, _, _, _, fallback = plan.curved_host
|
|
cx, cy = np.asarray(D2Q9_EX), np.asarray(D2Q9_EY)
|
|
sources = receivers.astype(np.int64)+cx[directions]+nx*cy[directions]
|
|
assert np.all((flags[receivers] & FLUID) != 0)
|
|
assert not np.any(flags[receivers] & SOLID)
|
|
assert np.all(solid[sources])
|
|
assert np.all(patch.fine_kind.ravel()[receivers] == ACTIVE)
|
|
assert np.all(patch.fine_kind.ravel()[sources] == ACTIVE)
|
|
# A source population slot has a unique writer, not just a unique cell.
|
|
assert len(set(zip(sources.tolist(), directions.tolist()))) == len(receivers)
|
|
expected = set()
|
|
for cell in np.flatnonzero((patch.fine_kind.ravel() == ACTIVE) & ~solid):
|
|
y, x = divmod(int(cell), nx)
|
|
for direction in range(1, 9):
|
|
sx, sy = x+int(cx[direction]), y+int(cy[direction])
|
|
if 0 <= sx < nx and 0 <= sy < ny and solid[sy*nx+sx]:
|
|
expected.add((int(cell), direction))
|
|
assert set(zip(receivers.tolist(), directions.tolist())) == expected
|
|
donors = receivers[q < .5].astype(np.int64)-cx[directions[q < .5]]-nx*cy[directions[q < .5]]
|
|
assert np.all(patch.fine_kind.ravel()[donors] == ACTIVE)
|
|
assert not np.any(solid[donors])
|
|
assert not fallback.any()
|
|
stencil = set(receivers.tolist()) | set(sources.tolist()) | set(donors.tolist())
|
|
assert stencil == set(plan.boundary_cells.tolist())
|
|
assert not set(plan.boundary_cells.tolist()) & set(plan.transfer_fine_cells.tolist())
|
|
assert not set(np.flatnonzero(solid).tolist()) & set(plan.transfer_fine_cells.tolist())
|
|
|
|
|
|
def test_transfer_dependencies_follow_two_pulls_and_expand_active_collision():
|
|
cfg, patch = setup_patch()
|
|
plan = refinement.prepare_fine_boundary(cfg, patch, [circle()])
|
|
ny, nx = patch.fine_shape
|
|
cells = nx*ny
|
|
# Closed form for exactly two hops, independent of the production traversal:
|
|
# the final inner ghost only pulls; its predecessor may have collided.
|
|
offsets = np.asarray(D2Q9_EX)+nx*np.asarray(D2Q9_EY)
|
|
visited = set()
|
|
for slot in patch.coalesce_fine.ravel():
|
|
direction, target = divmod(int(slot), cells)
|
|
assert patch.fine_kind.ravel()[target] == INNER_GHOST
|
|
middle = target-int(offsets[direction])
|
|
visited.update(((2, direction, target), (1, direction, middle)))
|
|
if patch.fine_kind.ravel()[middle] == ACTIVE:
|
|
visited.update((0, q, middle-int(offsets[q])) for q in range(9))
|
|
else:
|
|
assert patch.fine_kind.ravel()[middle] == INNER_GHOST
|
|
visited.add((0, direction, middle-int(offsets[direction])))
|
|
actual = {tuple(map(int, row)) for row in plan.transfer_reads}
|
|
assert actual == visited
|
|
assert len(plan.transfer_reads) == len(visited)
|
|
expected_fine = set(patch.ghost_fine.tolist()) | {cell for _, _, cell in visited}
|
|
expected_coarse = set(patch.ghost_coarse.tolist()) | set(
|
|
(patch.coalesce_coarse % (cfg.nx*cfg.ny)).tolist())
|
|
assert set(plan.transfer_fine_cells.tolist()) == expected_fine
|
|
assert set(plan.transfer_coarse_cells.tolist()) == expected_coarse
|
|
for coarse_slot, children in zip(patch.coalesce_coarse, patch.coalesce_fine):
|
|
coarse_q, coarse_cell = divmod(int(coarse_slot), cfg.nx*cfg.ny)
|
|
coarse_y, coarse_x = divmod(coarse_cell, cfg.nx)
|
|
fine_q, fine_cell = np.divmod(children, cells)
|
|
np.testing.assert_array_equal(fine_q, coarse_q)
|
|
fine_y, fine_x = np.divmod(fine_cell, nx)
|
|
physical_centers = np.column_stack((
|
|
patch.fine_origin[0]+.5*(fine_x+.5),
|
|
patch.fine_origin[1]+.5*(fine_y+.5)))
|
|
np.testing.assert_array_equal(physical_centers.mean(axis=0),
|
|
(coarse_x+.5, coarse_y+.5))
|
|
# This case must exercise the collision fan-in, not only ghost straight lines.
|
|
active_middle = [(q, cell) for stage, q, cell in visited
|
|
if stage == 1 and patch.fine_kind.ravel()[cell] == ACTIVE]
|
|
assert active_middle
|
|
_, sample = active_middle[0]
|
|
for incoming in range(9):
|
|
source = sample-int(D2Q9_EX[incoming])-nx*int(D2Q9_EY[incoming])
|
|
assert (0, incoming, source) in actual
|
|
|
|
|
|
def test_fine_plan_arrays_are_readonly_and_contiguous():
|
|
cfg, patch = setup_patch()
|
|
plan = refinement.prepare_fine_boundary(cfg, patch, [circle()])
|
|
assert len(plan.curved_host) == 8
|
|
for array in (plan.flags, *plan.curved_host, plan.boundary_cells,
|
|
plan.transfer_fine_cells, plan.transfer_coarse_cells, plan.transfer_reads):
|
|
assert array.flags.c_contiguous and not array.flags.writeable
|
|
assert array.size
|
|
with pytest.raises(ValueError):
|
|
array.flat[0] = array.flat[0]
|
|
|
|
|
|
def test_empty_objects_keep_existing_refinement_path():
|
|
cfg, patch = setup_patch()
|
|
assert refinement.prepare_fine_boundary(cfg, patch, []) is None
|
|
|
|
|
|
def test_fine_rejects_cross_level_low_q_donor_without_fallback():
|
|
cfg, patch = setup_patch()
|
|
source = circle(center=(5.05, 12.25), radius=.7)
|
|
# In fine coordinates C=(5.6,20), R=1.4. Receiver (4,20), dir +x
|
|
# has q=.2, but its required reverse donor (3,20) is not active.
|
|
local = CircleGeometry(5.6, 20., 1.4)
|
|
links = local.build_cut_links(patch.fine_shape[1], patch.fine_shape[0])
|
|
target, = [link for link in links if link.fluid_idx == 20*40+4 and link.dir == 1]
|
|
assert target.q == pytest.approx(.2)
|
|
assert patch.fine_kind[20, 3] != ACTIVE
|
|
with pytest.raises(ValueError):
|
|
refinement.prepare_fine_boundary(cfg, patch, [source])
|
|
|
|
|
|
def test_fine_rejects_active_stencil_that_still_touches_transfer():
|
|
cfg, patch = setup_patch()
|
|
source = circle(center=(6.05, 12.25), radius=.7)
|
|
# Move the cross-level example two fine cells inward: all boundary reads
|
|
# are now active, but the q=.2 reverse donor is still a transfer dependency.
|
|
local = CircleGeometry(7.6, 20., 1.4)
|
|
nx = patch.fine_shape[1]
|
|
solid = (local.build_flag_mask(nx, patch.fine_shape[0]) & SOLID) != 0
|
|
assert np.all(patch.fine_kind.ravel()[solid] == ACTIVE)
|
|
links = local.build_cut_links(nx, patch.fine_shape[0])
|
|
stencil = set()
|
|
for link in links:
|
|
offset = int(D2Q9_EX[link.dir])+nx*int(D2Q9_EY[link.dir])
|
|
stencil.update((link.fluid_idx, link.fluid_idx+offset))
|
|
if link.q < .5:
|
|
stencil.add(link.fluid_idx-offset)
|
|
assert np.all(patch.fine_kind.ravel()[list(stencil)] == ACTIVE)
|
|
with pytest.raises(ValueError, match='transfer footprint'):
|
|
refinement.prepare_fine_boundary(cfg, patch, [source])
|
|
|
|
|
|
@pytest.mark.parametrize('change', [dict(nx=25), dict(collision='TRT'), dict(viscosity=.06)])
|
|
def test_fine_rejects_mismatched_patch_configuration(change):
|
|
cfg, patch = setup_patch()
|
|
with pytest.raises(ValueError, match='matching'):
|
|
refinement.prepare_fine_boundary(replace(cfg, **change), patch, [circle()])
|
|
|
|
|
|
def test_fine_rejects_nonzero_control_even_without_motion():
|
|
cfg, patch = setup_patch()
|
|
obj = circle()
|
|
obj.control.ax = .001
|
|
with pytest.raises(ValueError, match='stationary'):
|
|
refinement.prepare_fine_boundary(cfg, patch, [obj])
|
|
|
|
|
|
@pytest.mark.parametrize('objects', [
|
|
[circle(center=(4.25, 12.25), radius=1.)],
|
|
[circle(center=(12.5, 12.5), radius=.05)],
|
|
[circle(), circle(center=(13., 12.25))],
|
|
])
|
|
def test_fine_rejects_ghost_solids_unresolved_circle_and_overlap(objects):
|
|
cfg, patch = setup_patch()
|
|
with pytest.raises(ValueError):
|
|
refinement.prepare_fine_boundary(cfg, patch, objects)
|
|
|
|
|
|
@pytest.mark.parametrize('margin, ok', [(4, True), (3, False)])
|
|
def test_fine_wall_guard_allows_read_only_bulk_donor_rows(margin, ok):
|
|
cfg = LBMConfig(nx=24, ny=24, viscosity=.05)
|
|
cfg.validate()
|
|
patch = FixedPatch2D.from_config(cfg, (4, margin, 20, 24-margin))
|
|
_, coarse, _ = refinement._transfer_footprint(patch)
|
|
rows = {k // cfg.nx for k in coarse}
|
|
assert min(rows) == margin-2 and max(rows) == cfg.ny-margin+1
|
|
source = circle(center=(12.25, 12.25))
|
|
# Footprint never violates the plain bulk-cell checks; only the wall rule differs.
|
|
refinement.prepare_fine_boundary(cfg, patch, [source])
|
|
if ok:
|
|
refinement.prepare_fine_boundary(cfg, patch, [source], wall_aware=True)
|
|
else:
|
|
with pytest.raises(ValueError, match='wall reconstruction'):
|
|
refinement.prepare_fine_boundary(cfg, patch, [source], wall_aware=True)
|