"""GPU gates for raw exchange. Run explicitly; requires PyCUDA and CUDA. Independent census reads actual staged incoming/outgoing populations, including non-body ACTIVE-boundary flux, then checks postcollision fluid momentum/mass. This is a discrete implementation gate, NOT a physical nozzle-load validation. """ from dataclasses import replace import unittest import numpy as np import pycuda.driver as cuda from CelerisLab.body.geometry.circle import CircleGeometry from CelerisLab.body.objects import SimObject from CelerisLab.body.jets import Jet, Tube, Parallel from CelerisLab.config import LBMConfig from CelerisLab.cuda.context import CudaContext from CelerisLab.lbm.refinement_stepper import FixedPatchStepper from CelerisLab.lbm.descriptors import D2Q9_EX, D2Q9_EY, SOLID CX, CY = np.asarray(D2Q9_EX), np.asarray(D2Q9_EY) OPP = np.array([np.flatnonzero((CX == -x) & (CY == -y))[0] for x, y in zip(CX, CY)]) W = np.array([4/9] + [1/9]*4 + [1/36]*4, np.float32) CFG = LBMConfig(nx=32, ny=24, viscosity=.1) BOUNDS = (8, 6, 24, 18) def circle(): return SimObject(17, geometry=CircleGeometry(16., 12., 2.15), center=(16., 12.), radius=2.15) def jet(): # Asymmetric aperture avoids hiding exchange errors by paired-jet cancellation. return Jet(Tube((16., 12.), 0., 3., .7), Parallel(0.), amplitude=0.) def state(shape, ux=.013): cu = CX[:, None, None]*ux values = W[:, None, None]*(1+3*cu+4.5*cu*cu-1.5*ux*ux) return np.broadcast_to(values, shape).astype(np.float32).copy() def download(pointer, shape, dtype=np.float32): result = np.empty(shape, dtype) cuda.memcpy_dtoh(result, pointer) return result def totals(g, leaf): values = g[:, leaf].astype(np.float64) return np.array([(values*CX[:, None]).sum(), (values*CY[:, None]).sum(), values.sum()]) def census(run, index, buffer): """Independent incoming/outgoing face census on ACTIVE fluid cells. Reads the receiver's actual block window, not a freshly gathered owner image (that would erase BFL/restriction staging). Classifies external flux solely from geometry, without consulting raw observation arrays. """ level = run._levels[index] layout, b = level.layout, level.boundary g = download(buffer, layout.blocked_shape) dense = layout.to_dense(g) leaf = run._leaf[index] ny, nx = leaf.shape solid = (b.host.flags.reshape(ny, nx) & SOLID) != 0 y, x = np.nonzero(leaf) owner = layout.owner[y, x].astype(np.intp) cells = int(np.prod(layout.local_shape)) block, cell = np.divmod(owner, cells) flat = g.reshape(layout.nblocks, 9, cells) external = np.zeros(3, np.float64) for q, (cx, cy) in enumerate(zip(CX, CY)): sy, sx = y-cy, x-cx incoming = flat[block, q, cell-cy*layout.local_shape[1]-cx].astype(np.float64) entering = ~leaf[sy, sx] & ~solid[sy, sx] external += np.array([cx, cy, 1])*incoming[entering].sum() dy, dx = y+cy, x+cx leaving = ~leaf[dy, dx] & ~solid[dy, dx] outgoing = flat[block, q, cell].astype(np.float64) external -= np.array([cx, cy, 1])*outgoing[leaving].sum() receiver, direction, _, _, rx, ry, _, _ = b.host.curved_host owner = layout.owner.ravel()[receiver].astype(np.intp) block, cell = np.divmod(owner, cells) a = flat[block, direction, cell].astype(np.float64) source = cell+CY[direction]*layout.local_shape[1]+CX[direction] incoming = flat[block, OPP[direction], source].astype(np.float64) jx, jy = CX[direction]*(a+incoming), CY[direction]*(a+incoming) independent = np.column_stack((jx, jy, rx*jy-ry*jx, incoming-a)) return totals(dense, leaf), external, independent class RawExchangeGPU(unittest.TestCase): def test_row_ledger_interval_jet_switch_and_no_population_effect(self): with CudaContext(0) as ctx: for block in ((32, 40), (16, 20), (8, 10)): with FixedPatchStepper(CFG, BOUNDS, context=ctx, fine_block=block, fine_objects=[circle()], jets=[jet()], observe_raw_exchange=True) as measured, \ FixedPatchStepper(CFG, BOUNDS, context=ctx, fine_block=block, fine_objects=[circle()], jets=[jet()]) as plain: arrays = [state(l.shape) for l in measured._levels] measured.initialize(*arrays) plain.initialize(*arrays) b = measured._level_boundaries[1] count = len(b.tables['cell']) kernel = measured._levels[1].kernel samples = [] def audited_step(src, dst, *args, **kw): measured._stream.synchronize() before, external, independent = census(measured, 1, src) raw = download(b.last_raw, (count, 4)) np.testing.assert_allclose(raw, independent, rtol=2e-6, atol=2e-7) kernel(src, dst, *args, **kw) measured._stream.synchronize() after = measured._levels[1].layout.to_dense( download(dst, measured._levels[1].layout.blocked_shape)) change = totals(after, measured._leaf[1])-before expected = external + np.array([-raw[:, 0].astype(float).sum(), -raw[:, 1].astype(float).sum(), raw[:, 3].astype(float).sum()]) np.testing.assert_allclose(change, expected, rtol=2e-4, atol=8e-5) samples.append(raw.copy()) measured._levels[1].kernel = audited_step for speed in (0., .002, -.001, 0.): measured.set_jet_amplitude(0, speed) plain.set_jet_amplitude(0, speed) samples.clear() measured.step() plain.step() for a, e in zip(measured.download(), plain.download()): np.testing.assert_array_equal(a, e) integral = download(b.interval_integral, (count, 4), np.float64) expected = sum(s.astype(float)*[.25, .25, .125, .25] for s in samples) np.testing.assert_array_equal(integral, expected) last = download(b.last_raw, (count, 4)) np.testing.assert_array_equal(last, samples[-1]) obs = measured.read_raw_exchange(0) assert obs['sample_count'] == 2 assert b.tables['aperture'].any() measured.step(0) np.testing.assert_array_equal(measured.read_raw_exchange(0)['total']['impulse'], obs['total']['impulse']) def test_dense_fine_storage_reference(self): """Same fine resolution/initial data with dense storage: a numerical reference for cut-link rows, not a uniform-whole-domain convergence study.""" with CudaContext(0) as ctx: with FixedPatchStepper(CFG, BOUNDS, context=ctx, fine_objects=[circle()], jets=[jet()]) as dense, \ FixedPatchStepper(CFG, BOUNDS, context=ctx, fine_block=(16, 20), fine_objects=[circle()], jets=[jet()], observe_raw_exchange=True) as blocked: arrays = [state(l.shape) for l in blocked._levels] dense.initialize(*arrays) blocked.initialize(*arrays) for speed in (0., .002, 0.): dense.set_jet_amplitude(0, speed) blocked.set_jet_amplitude(0, speed) dense.step(); blocked.step() for a, e in zip(blocked.download(), dense.download()): np.testing.assert_array_equal(a, e) count = dense._fine_curved.count expected = download(dense._fine_curved.link_exchange_gpu, (count, 3)) actual = download(blocked._level_boundaries[1].last_raw, (count, 4)) np.testing.assert_array_equal(actual[:, :3], expected) def test_depth_four_interval_storage_and_sampling(self): cfg = LBMConfig(nx=48, ny=48, viscosity=.02) obj = SimObject(0, geometry=CircleGeometry(16., 16., .5), center=(16., 16.), radius=.5) drive = Jet(Tube((16., 16.), 0., .75, .2), Parallel(0.), amplitude=0.) settings = dict(refinements=[[(6, 6, 26, 26)], [(20, 20, 44, 44)], [(56, 56, 72, 72)]], level_blocks=[(4, 4)]*4, level_objects=[(), (), (), (obj,)], level_jets=[(), (), (), (drive,)]) with CudaContext(0) as ctx: with FixedPatchStepper(cfg, context=ctx, observe_raw_exchange=True, **settings) as run, \ FixedPatchStepper(cfg, context=ctx, **settings) as plain: arrays = [state(l.shape, ux=0.) for l in run._levels] run.initialize(*arrays); plain.initialize(*arrays) b = run._level_boundaries[3] count = len(b.tables['cell']) kernel = run._levels[3].kernel samples, independent_samples = [], [] def audited_step(src, dst, *args, **kw): run._stream.synchronize() before, external, independent = census(run, 3, src) raw = download(b.last_raw, (count, 4)) np.testing.assert_allclose(raw, independent, rtol=2e-6, atol=2e-7) kernel(src, dst, *args, **kw) run._stream.synchronize() after = run._levels[3].layout.to_dense( download(dst, run._levels[3].layout.blocked_shape)) change = totals(after, run._leaf[3])-before expected = external + np.array([-independent[:, 0].sum(), -independent[:, 1].sum(), independent[:, 3].sum()]) np.testing.assert_allclose(change, expected, rtol=2e-4, atol=8e-5) samples.append(raw) independent_samples.append(independent) run._levels[3].kernel = audited_step for speed in (0., .002, 0.): run.set_jet_amplitude(0, speed, level=3) plain.set_jet_amplitude(0, speed, level=3) samples.clear(); independent_samples.clear() run.step(); plain.step() self.assertEqual(len(samples), 8) weights = np.array([1/64, 1/64, 1/512, 1/64]) integral = download(b.interval_integral, (count, 4), np.float64) np.testing.assert_array_equal(integral, sum(v.astype(float)*weights for v in samples)) np.testing.assert_allclose(integral, sum(v*weights for v in independent_samples), rtol=2e-6, atol=2e-7) np.testing.assert_array_equal(download(b.last_raw, (count, 4)), samples[-1]) for actual, expected in zip(run.download(), plain.download()): np.testing.assert_array_equal(actual, expected) result = run.read_raw_exchange(0, level=3) self.assertEqual(result['sample_count'], 8) self.assertEqual(result['time_start'], run.step_count-1) if run.step_count == 1: np.testing.assert_allclose(result['total']['force'], 0, atol=2e-7) self.assertAlmostEqual(result['total']['torque'], 0, delta=2e-7) def test_uniform_whole_domain_fine_short_time_reference(self): """Same root geometry, nu and Jet speed on a whole-domain fine mesh. Eight local hops are a short-time equivalence check, not convergence. The test adapter only enables the standalone kernel's optional inputs. """ from CelerisLab.lbm.blocks import BlockLayout2D from CelerisLab.lbm.block_stepper import PeriodicBlockStepper fine_cfg = replace(CFG, nx=64, ny=48, viscosity=.2) obj = SimObject(0, geometry=CircleGeometry(31.5, 23.5, 4.3), center=(31.5, 23.5), radius=4.3) layout = BlockLayout2D.build(48, 64, 8, 8) with CudaContext(0) as ctx: with PeriodicBlockStepper(fine_cfg, layout, context=ctx, objects=[obj]) as uniform, \ FixedPatchStepper(CFG, BOUNDS, context=ctx, fine_block=(16, 20), fine_objects=[circle()], jets=[jet()], observe_raw_exchange=True) as amr: amr.initialize(*(state(l.shape) for l in amr._levels)) uniform.initialize(state((9, 48, 64))) b = amr._level_boundaries[1] count = len(b.tables['cell']) links = obj.geometry.build_cut_links(64, 48) np.testing.assert_array_equal([v.dir for v in links], b.tables['dir']) np.testing.assert_allclose([v.q for v in links], b.tables['q'], atol=2e-7, rtol=0) # Same fine-grid force arms after translation to the same center. rx = uniform._upload(b.tables['rx']) ry = uniform._upload(b.tables['ry']) velocity = uniform._upload(b.tables['surface_velocity']) raw = uniform._allocate(count*16) original = uniform._curved samples = [] def measured_boundary(*args, **kw): original(*args[:7], velocity, rx, ry, raw, **kw) uniform._stream.synchronize() samples.append(download(raw, (count, 4))) uniform._curved = measured_boundary for speed in (0., .002, -.001, 0.): amr.set_jet_amplitude(0, speed) cuda.memcpy_htod(velocity, b.tables['surface_velocity']) samples.clear() amr.step(); uniform.step(2) expected = sum(v.astype(float)*[.25, .25, .125, .25] for v in samples) actual = download(b.interval_integral, (count, 4), np.float64) np.testing.assert_allclose(actual, expected, rtol=1e-5, atol=2e-6) if __name__ == '__main__': unittest.main()