"""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 def fine(block): """Two-level circle+jet construction on BOUNDS, blocked with a (12, 8) root tile.""" return dict(refinements=[BOUNDS], level_blocks=[(12, 8), block], level_objects=[(), (circle(),)], level_jets=[(), (jet(),)]) 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, context=ctx, **fine(block), observe_raw_exchange=True) as measured, \ FixedPatchStepper(CFG, context=ctx, **fine(block)) 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, level=1) plain.set_jet_amplitude(0, speed, level=1) 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, level=1) assert obs['sample_count'] == 2 assert b.tables['aperture'].any() measured.step(0) np.testing.assert_array_equal(measured.read_raw_exchange(0, level=1)['total']['impulse'], obs['total']['impulse']) 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. """ 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) with CudaContext(0) as ctx: with FixedPatchStepper(fine_cfg, context=ctx, refinements=[], level_blocks=[(8, 8)], level_objects=[(obj,)]) as uniform, \ FixedPatchStepper(CFG, context=ctx, **fine((16, 20)), 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: the # single level's static BFL gets the AMR Jet column and force arms; the raw # output pointer is substituted at launch (no stepper-side accumulation). root = uniform._level_boundaries[0] velocity = uniform._upload(b.tables['surface_velocity']) root.device.update(surface_velocity=velocity, rx=uniform._upload(b.tables['rx']), ry=uniform._upload(b.tables['ry'])) raw = uniform._allocate(count*16) original = root.kernel samples = [] def measured_boundary(*args, **kw): original(*args[:10], raw, **kw) # BlockCurvedBoundary raw-exchange output uniform._stream.synchronize() samples.append(download(raw, (count, 4))) root.kernel = measured_boundary for speed in (0., .002, -.001, 0.): amr.set_jet_amplitude(0, speed, level=1) 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()