diff --git a/src/CelerisLab/lbm/block_stepper.py b/src/CelerisLab/lbm/block_stepper.py index 44a1109..13c1c9c 100644 --- a/src/CelerisLab/lbm/block_stepper.py +++ b/src/CelerisLab/lbm/block_stepper.py @@ -21,6 +21,7 @@ import pycuda.driver as cuda from ..cuda import compiler_v2 as compiler from .blocks import NQ, BlockLayout2D, blocked_link_tables from .descriptors import D2Q9_EX, D2Q9_EY, FLUID, SOLID +from .population_storage import _StoredPopulations, _encode_physical, _decode_stored MAX_GRID_DIM = 65535 # CUDA gridDim.y/z limit (rows, blocks) @@ -85,10 +86,10 @@ class PeriodicBlockStepper: checked.validate() if (not checked.is_d2q9 or checked.collision != "SRT" or checked.streaming != "double_buffer" or checked.store_precision != "FP32" - or checked.ddf_shifting or checked.les_enabled + or checked.les_enabled or (checked.ny, checked.nx) != layout.shape): raise ValueError("PeriodicBlockStepper requires matching D2Q9/SRT/FP32/" - "unshifted/double_buffer config without LES") + "double_buffer config without LES") ly, lx = layout.local_shape pairs = len(layout.halo_dst) tpb = checked.threads_per_block @@ -99,6 +100,7 @@ class PeriodicBlockStepper: links = 0 if self.boundary is None else len(self.boundary["cell"]) if links >= 2**31: raise ValueError("cut links exceed launch limits") + self._ddf_shifting = bool(checked.ddf_shifting) self.layout = layout self.ctx = context self._links = links @@ -177,14 +179,20 @@ class PeriodicBlockStepper: """Upload dense G_0; the spare buffer starts NaN so stale halos cannot hide.""" with self._lock, self.ctx.activate(): self._require("new", "ready") - dense = np.array(g, dtype=np.float32, order="C", copy=True) + dense = np.array(g, dtype=np.float64 if self._ddf_shifting else np.float32, + order="C", copy=True) if dense.shape != (NQ, *self.layout.shape): raise ValueError("Expected a Q-major (9, ny, nx) array") fluid = np.ones(self.layout.shape, bool) if self.boundary is None else ~self.boundary["solid"] if not np.isfinite(dense[:, fluid]).all() or (dense[:, fluid].sum(axis=0) <= 0).any(): raise ValueError("Fluid populations must be finite with positive density") dense[:, ~fluid] = np.nan # Solid cores are deliberately invalid. - blocked = self.layout.to_blocked(dense) + encoded = _encode_physical(dense, self._ddf_shifting) + values = encoded[:, fluid] + if (not np.isfinite(values).all() + or (values.sum(axis=0, dtype=np.float64) + int(self._ddf_shifting) <= 0).any()): + raise ValueError("Stored fluid populations must be finite with positive density") + blocked = self.layout.to_blocked(encoded) spare = np.full_like(blocked, np.nan) try: cuda.memcpy_htod(self._current, blocked) @@ -225,7 +233,12 @@ class PeriodicBlockStepper: raise def download(self): - """Return dense (9, ny, nx) FP32 core populations; no host cache.""" + """Return physical core populations (FP64 when shifted); no host cache.""" + snapshot = self._download_raw() + return _decode_stored(snapshot.values[0], snapshot.shifted) + + def _download_raw(self): + """Return one tagged dense FP32 population array without conversion.""" with self._lock, self.ctx.activate(): self._require("ready") blocked = np.empty(self.layout.blocked_shape, np.float32) @@ -234,7 +247,7 @@ class PeriodicBlockStepper: except BaseException: self._status = "failed" raise - return self.layout.to_dense(blocked) + return _StoredPopulations((self.layout.to_dense(blocked),), self._ddf_shifting) def close(self): """Release owned resources; retain them for retry if context activation fails.""" diff --git a/src/CelerisLab/lbm/population_storage.py b/src/CelerisLab/lbm/population_storage.py new file mode 100644 index 0000000..ffada5e --- /dev/null +++ b/src/CelerisLab/lbm/population_storage.py @@ -0,0 +1,24 @@ +"""Private D2Q9 FP32 storage boundary for the opt-in blocked steppers.""" +from dataclasses import dataclass +import numpy as np + +# Unit-density reference weights in FP64; sum(w)=1 matches rho=1+sum(g). +_WEIGHTS = np.array([4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36], + dtype=np.float64)[:, None, None] + + +@dataclass(frozen=True, eq=False) +class _StoredPopulations: + values: tuple + shifted: bool + + +def _encode_physical(values, shifted): + if not shifted: + return np.array(values, dtype=np.float32, order="C", copy=True) + return np.ascontiguousarray(np.asarray(values, dtype=np.float64) - _WEIGHTS, + dtype=np.float32) + + +def _decode_stored(values, shifted): + return values.astype(np.float64) + _WEIGHTS if shifted else values diff --git a/src/CelerisLab/lbm/refinement_stepper.py b/src/CelerisLab/lbm/refinement_stepper.py index 7e05e01..ef98b24 100644 --- a/src/CelerisLab/lbm/refinement_stepper.py +++ b/src/CelerisLab/lbm/refinement_stepper.py @@ -25,6 +25,7 @@ from .refinement_tables import (blocked_base_tables, blocked_fine_tables, # noq BLOCKED_BOUNDARY_DEVICE_COLUMNS) from .descriptors import SOLID from .raw_exchange import summarize_exchange +from .population_storage import _StoredPopulations, _encode_physical, _decode_stored @dataclass @@ -238,6 +239,9 @@ class FixedPatchStepper: level_blocks = [None if b is None else tuple(b) for b in level_blocks] _check_levels(refinements, level_blocks, level_objects, level_jets) depth = len(refinements) + 1 + self._ddf_shifting = bool(cfg.ddf_shifting) + if self._ddf_shifting and any(b is None for b in level_blocks): + raise ValueError("shifted populations require all-blocked levels; dense/mixed is unsupported") self._wall_aware = coarse_boundary == "slot_reservoir" if self._wall_aware and (cfg.y_wall_bc != "on_site_velocity" or cfg.top_wall_bc != "reservoir" or cfg.velocity != 0 or not 2 <= cfg.wall_slot_width < cfg.nx): @@ -469,23 +473,54 @@ class FixedPatchStepper: for array, leaf in zip(arrays, self._leaf): array[:, ~leaf] = np.nan - def initialize(self, *level_g, coarse_g=None, fine_g=None, child_g=None): - """Upload an explicit G_0, one Q-major array per level; non-leaves are discarded. + def _check_storage(self): + if any(bool(c.ddf_shifting) != self._ddf_shifting for c in self._level_cfgs): + raise ValueError("mixed population representations are unsupported") + if self._ddf_shifting and any(layout is None for layout in self._layouts): + raise ValueError("shifted populations require all-blocked levels") - Pass the level arrays positionally (L0 first), or use the original - two/three-level names coarse_g, fine_g[, child_g]; mixing is rejected. + def initialize(self, *level_g, coarse_g=None, fine_g=None, child_g=None): + """Upload physical postcollision G, encoding before any FP32 rounding. + + Shifted input should be physical FP64 to retain small perturbations. + Public initialization resets the step count; non-leaves are discarded. """ given = _initial_levels(level_g, coarse_g, fine_g, child_g) with self._lock, self.ctx.activate(): self._require("new", "ready") + self._check_storage() if len(given) != len(self._levels): raise ValueError("Expected one Q-major array per level") - arrays = [np.array(a, dtype=np.float32, order="C", copy=True) for a in given] + dtype = np.float64 if self._ddf_shifting else np.float32 + arrays = [np.array(a, dtype=dtype, order="C", copy=True) for a in given] if any(a.shape != level.shape for a, level in zip(arrays, self._levels)): raise ValueError("Expected Q-major (9,ny,nx) arrays for every level") for array, mask in zip(arrays, self._leaf): if not np.isfinite(array[:, mask]).all() or (array[:, mask].sum(axis=0) <= 0).any(): raise ValueError("Leaf populations must be finite with positive density") + encoded = tuple(_encode_physical(a, self._ddf_shifting) for a in arrays) + self._initialize_raw(_StoredPopulations(encoded, self._ddf_shifting)) + + def _initialize_raw(self, snapshot): + """Upload an internally tagged same-representation state without a codec.""" + with self._lock, self.ctx.activate(): + self._require("new", "ready") + self._check_storage() + if (not isinstance(snapshot, _StoredPopulations) + or type(snapshot.shifted) is not bool or snapshot.shifted != self._ddf_shifting): + raise ValueError("stored population representation mismatch") + if len(snapshot.values) != len(self._levels): + raise ValueError("Expected one stored array per level") + if any(not isinstance(a, np.ndarray) or a.dtype != np.float32 for a in snapshot.values): + raise ValueError("stored populations must be FP32") + arrays = [np.array(a, order="C", copy=True) for a in snapshot.values] + if any(a.shape != level.shape for a, level in zip(arrays, self._levels)): + raise ValueError("Expected Q-major (9,ny,nx) arrays for every level") + for array, mask in zip(arrays, self._leaf): + values = array[:, mask] + if (not np.isfinite(values).all() + or (values.sum(axis=0, dtype=np.float64) + int(self._ddf_shifting) <= 0).any()): + raise ValueError("Leaf populations must be finite with positive density") self._mask_nonleaves(*arrays) arrays = [a if level.layout is None else level.layout.to_blocked(a) for a, level in zip(arrays, self._levels)] @@ -518,7 +553,7 @@ class FixedPatchStepper: construction/upload errors never invalidate the old instance; failure to synchronize/read its device state does, as with step/download. """ - from .synchronized_remap import (RemapTopology, remap_synchronized, + from .synchronized_remap import (RemapTopology, remap_synchronized, _remap_stored, same_physical_boundary, boundary_velocity) with self._lock: @@ -558,14 +593,16 @@ class FixedPatchStepper: with self.ctx.activate(): try: self._stream.synchronize() - source = self.download() + source = self._download_raw() if self._ddf_shifting else self.download() except BaseException: self._status = "failed" raise count = self._step_count - transferred = remap_synchronized( + transfer = _remap_stored if self._ddf_shifting else remap_synchronized + options = dict(shifted=self._ddf_shifting) if self._ddf_shifting else {} + transferred = transfer( old, new, source, base_count=count, level_steps=(count, 2*count, 4*count), - phase="postcollision_g") + phase="postcollision_g", **options) candidate = type(self)( cfg, pairs[0].bounds, context=self.ctx, base_block=base_block, fine_block=fine_block, child_block=child_block, @@ -573,7 +610,10 @@ class FixedPatchStepper: child_objects=child_objects, child_jets=jets, observe_raw_exchange=self._observe_raw_exchange) try: - candidate.initialize(*transferred.g) + if self._ddf_shifting: + candidate._initialize_raw(_StoredPopulations(transferred.g, transferred.shifted)) + else: + candidate.initialize(*transferred.g) candidate._step_count = count except BaseException: try: @@ -772,9 +812,15 @@ class FixedPatchStepper: reference_center=self._exchange_centers[level][body_id]) def download(self): - """Return synchronized G arrays (one per level), NaN on all non-leaves; no host cache.""" + """Return physical G (FP64 for shifted storage), NaN on non-leaves.""" + snapshot = self._download_raw() + return tuple(_decode_stored(a, snapshot.shifted) for a in snapshot.values) + + def _download_raw(self): + """Current owned populations with a representation tag; no add/subtract w.""" with self._lock, self.ctx.activate(): self._require("ready") + self._check_storage() arrays = [np.empty(level.shape if level.layout is None else level.layout.blocked_shape, np.float32) for level in self._levels] try: @@ -786,7 +832,7 @@ class FixedPatchStepper: arrays = [a if level.layout is None else level.layout.to_dense(a) for a, level in zip(arrays, self._levels)] self._mask_nonleaves(*arrays) - return tuple(arrays) + return _StoredPopulations(tuple(arrays), self._ddf_shifting) def close(self): """Release owned resources; retain them for retry if context activation fails.""" diff --git a/src/CelerisLab/lbm/synchronized_remap.py b/src/CelerisLab/lbm/synchronized_remap.py index f39d7d2..d280e28 100644 --- a/src/CelerisLab/lbm/synchronized_remap.py +++ b/src/CelerisLab/lbm/synchronized_remap.py @@ -16,6 +16,7 @@ from .refinement import (ACTIVE, FixedPatch2D, FixedPatchSet2D, FineBoundary2D, nest_patch_set, prepare_fine_boundary) from .descriptors import SOLID from ..body.jets import surface_link_velocity +from .population_storage import _StoredPopulations @dataclass(frozen=True, eq=False) @@ -118,6 +119,7 @@ class RemapCandidate: phase: str origin: tuple spacing: float + shifted: bool = False def _same_parent(a, b): @@ -126,7 +128,7 @@ def _same_parent(a, b): and a.tau_coarse == b.tau_coarse and a.tau_fine == b.tau_fine) -def _read_leaves(topology, arrays, layouts): +def _read_leaves(topology, arrays, layouts, *, shifted=False): """Unpack only owned cores; nonleaf values are not evidence of validity.""" if len(arrays) != 3 or len(layouts) != 3: raise ValueError("expected three G arrays and three layouts") @@ -140,8 +142,10 @@ def _read_leaves(topology, arrays, layouts): a = layout.to_dense(a) if a.shape != (9, *leaf.shape): raise ValueError("G must have Q-major (9, ny, nx) shape") + if shifted and a.dtype != np.float32: + raise ValueError("shifted stored populations must be FP32") values = a[:, leaf] - if not np.isfinite(values).all() or (values.sum(axis=0, dtype=np.float64) <= 0).any(): + if not np.isfinite(values).all() or (values.sum(axis=0, dtype=np.float64) + int(shifted) <= 0).any(): raise ValueError("required leaf donors must be finite with positive density") result.append(a) if len({a.dtype for a in result}) != 1: @@ -166,6 +170,20 @@ def remap_synchronized(old, new, old_g, *, base_count, level_steps, phase, the source dtype. No equilibrium reconstruction or nonequilibrium scaling. All destination nonleaves are NaN. Inputs and their write flags are intact. """ + return _remap(old, new, old_g, base_count=base_count, level_steps=level_steps, + phase=phase, origin=origin, spacing=spacing, layouts=layouts) + + +def _remap_stored(old, new, snapshot, *, shifted, **kwargs): + """Private same-representation transfer; public remap remains physical.""" + if (not isinstance(snapshot, _StoredPopulations) or type(shifted) is not bool + or type(snapshot.shifted) is not bool or snapshot.shifted != shifted): + raise ValueError("stored population representation mismatch") + return _remap(old, new, snapshot.values, shifted=shifted, **kwargs) + + +def _remap(old, new, old_g, *, base_count, level_steps, phase, + origin=(0., 0.), spacing=1., layouts=None, shifted=False): if (not isinstance(old, RemapTopology) or not isinstance(new, RemapTopology) or not _same_parent(old.parent, new.parent) or old.child.tau_fine != new.child.tau_fine): @@ -189,7 +207,7 @@ def remap_synchronized(old, new, old_g, *, base_count, level_steps, phase, if (frame.shape != (2,) or not np.isfinite(frame).all() or not np.isscalar(spacing) or not np.isfinite(spacing) or spacing <= 0): raise ValueError("frame needs a finite (x,y) origin and positive spacing") - source = _read_leaves(old, old_g, (None,)*3 if layouts is None else layouts) + source = _read_leaves(old, old_g, (None,)*3 if layouts is None else layouts, shifted=shifted) candidate = [np.full((9, *mask.shape), np.nan, dtype=source[0].dtype) for mask in new.leaf_masks] for level, (before, after) in enumerate(zip(old.leaf_masks, new.leaf_masks)): @@ -210,6 +228,6 @@ def remap_synchronized(old, new, old_g, *, base_count, level_steps, phase, # Scale first to avoid an overflowing intermediate sum for finite FP64 G. mean = sum(source[2][:, cy, cx].astype(np.float64) * .25 for cy, cx in children) candidate[1][:, y, x] = mean - _read_leaves(new, candidate, (None,)*3) + _read_leaves(new, candidate, (None,)*3, shifted=shifted) return RemapCandidate(new, tuple(candidate), int(base_count), tuple(map(int, level_steps)), - phase, tuple(map(float, frame)), float(spacing)) + phase, tuple(map(float, frame)), float(spacing), shifted) diff --git a/tests/unit/test_refinement.py b/tests/unit/test_refinement.py index 48c84ee..1c72c79 100644 --- a/tests/unit/test_refinement.py +++ b/tests/unit/test_refinement.py @@ -143,7 +143,7 @@ class FixedPatchTests(unittest.TestCase): def test_unsupported_backend_and_fine_relaxation_rejected(self): for kwargs in ({"collision": "TRT"}, {"collision": "MRT"}, {"streaming": "esopull"}, - {"ddf_shifting": True}, {"les_enabled": True}, + {"les_enabled": True}, {"store_precision": "FP16S"}, {"lattice_model": "D3Q19"}, {"data_type": "FP64"}, {"omega_min": 1.0}): cfg = config(**kwargs) diff --git a/tests/unit/test_shifted_storage.py b/tests/unit/test_shifted_storage.py new file mode 100644 index 0000000..bcca95b --- /dev/null +++ b/tests/unit/test_shifted_storage.py @@ -0,0 +1,184 @@ +"""CPU host contracts; mocked CUDA executes real lifecycle, not numerical kernels.""" +import ast +import types +from pathlib import Path +from contextlib import nullcontext +from dataclasses import replace +import unittest +from unittest.mock import Mock, patch +import numpy as np +from CelerisLab.lbm.population_storage import (_StoredPopulations, _WEIGHTS, + _encode_physical, _decode_stored) +from CelerisLab.lbm.synchronized_remap import _remap_stored, remap_synchronized +from tests.unit.test_synchronized_remap import topology, A, B +from tests.unit import test_refinement_three_level_multipatch as backend + + +def shifted_fields(mesh): + rng = np.random.default_rng(129) + arrays = [] + for leaf in mesh.leaf_masks: + a = rng.normal(0, 1e-7, (9, *leaf.shape)).astype(np.float32) + a[:, ~leaf] = np.nan + y, x = np.nonzero(leaf) + a[0, y[:6], x[:6]] = np.array([0., -0., 2**-80, -2**-80, 2**-140, -2**-140], np.float32) + arrays.append(a) + return tuple(arrays) + + +def remap(old, new, arrays, shifted=True): + return _remap_stored(old, new, _StoredPopulations(arrays, shifted), shifted=True, + base_count=7, level_steps=(7, 14, 28), phase='postcollision_g') + + +def ready(): + with patch.object(backend, 'CFG', replace(backend.CFG, ddf_shifting=True)): + run, module, allocations = backend.mock_backend(children=A) + storage = {} + module.cuda.memcpy_htod.side_effect = lambda p, a: storage.__setitem__(id(p), np.array(a, copy=True)) + module.cuda.memcpy_dtoh.side_effect = lambda a, p: np.copyto(a, storage[id(p)]) + run._initialize_raw(_StoredPopulations(shifted_fields(topology(A)), True)) + run._step_count = 7 + return run, module, allocations, storage + + +class ShiftedStorage(unittest.TestCase): + def test_codec_preserves_sub_fp32_physical_perturbation(self): + physical = np.broadcast_to(_WEIGHTS, (9, 2, 2)).copy() + physical[1] += 2**-35 + encoded = _encode_physical(physical, True) + self.assertTrue((encoded[1] == np.float32(2**-35)).all()) + np.testing.assert_array_equal(_decode_stored(encoded, True), physical) + self.assertEqual(_decode_stored(encoded, True).dtype, np.float64) + np.testing.assert_array_equal(_encode_physical(_WEIGHTS, True), np.zeros((9,1,1), np.float32)) + self.assertEqual(_encode_physical(physical, False).dtype, np.float32) + + def test_actual_fill_loop_bitwise_copy_and_single_rounding_restriction(self): + old, new = topology(A), topology(B) + arrays = shifted_fields(old) + saved = [a.tobytes() for a in arrays] + out = remap(old, new, arrays) + self.assertTrue(out.shifted) + for a, b, om, nm, before in zip(arrays, out.g, old.leaf_masks, new.leaf_masks, saved): + self.assertEqual(a[:, om & nm].tobytes(), b[:, om & nm].tobytes()) + self.assertTrue(np.isnan(b[:, ~nm]).all()) + self.assertEqual(a.tobytes(), before) + self.assertFalse(np.shares_memory(a,b)) + y,x = np.nonzero(new.leaf_masks[2] & ~old.leaf_masks[2]) + self.assertEqual(out.g[2][:,y,x].tobytes(), arrays[1][:,y//2,x//2].tobytes()) + y,x = np.nonzero(new.leaf_masks[1] & ~old.leaf_masks[1]) + mean = sum(arrays[2][:,2*y+dy,2*x+dx].astype(np.float64)*.25 + for dy in (0,1) for dx in (0,1)) + self.assertEqual(out.g[1][:,y,x].tobytes(), mean.astype(np.float32).tobytes()) + def integral(mesh, values): + return sum(a[:,m].sum(axis=1,dtype=np.float64)/4**l + for l,(a,m) in enumerate(zip(values,mesh.leaf_masks))) + cast_error = abs(mean.astype(np.float32).astype(np.float64)-mean).sum(axis=1)/4 + self.assertTrue(np.all(abs(integral(old, arrays)-integral(new,out.g)) <= cast_error+2e-18)) + + def test_shifted_density_and_tag_validation(self): + old,new = topology(A),topology(B) + a = tuple(np.where(m[None], np.float32(-.001), np.nan).repeat(9,axis=0) for m in old.leaf_masks) + remap(old,new,a) + with self.assertRaisesRegex(ValueError,'representation'): remap(old,new,a,False) + with self.assertRaisesRegex(ValueError,'FP32'): remap(old,new,tuple(x.astype(np.float64) for x in a)) + a[0][:,old.leaf_masks[0]] = -1 + with self.assertRaisesRegex(ValueError,'positive density'): remap(old,new,a) + rest = tuple(np.where(m[None],np.float32(0),np.nan).repeat(9,axis=0) for m in old.leaf_masks) + remap(old,new,rest) + with self.assertRaisesRegex(ValueError,'positive density'): + remap_synchronized(old,new,rest,base_count=7,level_steps=(7,14,28),phase='postcollision_g') + + def test_runtime_remap_bypasses_public_codec_and_preserves_both_buffers(self): + run,module,allocations,storage = ready() + clone = None + try: + before = run._download_raw() + expected = remap(topology(A),topology(B),before.values) + with patch.object(run,'download',side_effect=AssertionError('public decode used')): + clone = run.remap_child(B) + self.assertEqual((run.step_count,clone.step_count),(7,7)) + self.assertFalse(clone._exchange_valid) + for level, want, got in zip(clone._levels,expected.g,clone._download_raw().values): + self.assertEqual(want.tobytes(),got.tobytes()) + for pointer in (level.current,level.spare): + np.testing.assert_array_equal(storage[id(pointer)],level.layout.to_blocked(want)) + for a,b in zip(before.values,run._download_raw().values): self.assertEqual(a.tobytes(),b.tobytes()) + for a in clone.download(): self.assertEqual(a.dtype,np.float64) + finally: + if clone is not None: clone.close() + run.close() + + def test_runtime_public_encode_and_raw_mismatch_fail_before_upload(self): + run,module,_,_ = ready() + try: + physical = tuple(np.broadcast_to(_WEIGHTS, level.shape).copy() for level in run._levels) + for a in physical: a[1] += 2**-35 + run.initialize(*physical) + for a,mask in zip(run._download_raw().values,run._leaf): + self.assertTrue((a[1,mask] == np.float32(2**-35)).all()) + self.assertEqual(run.step_count,0) + module.cuda.memcpy_htod.reset_mock() + with self.assertRaisesRegex(ValueError,'representation'): + run._initialize_raw(_StoredPopulations(run._download_raw().values,False)) + module.cuda.memcpy_htod.assert_not_called() + run._level_cfgs[0] = replace(run._level_cfgs[0],ddf_shifting=False) + with self.assertRaisesRegex(ValueError,'mixed'): run._download_raw() + finally: run.close() + + def test_dense_or_mixed_shifted_rejected_before_compilation(self): + module=backend.host_module() + cfg=replace(backend.CFG,ddf_shifting=True) + for kwargs in ({},dict(base_block=(4,4)),dict(fine_block=(4,4))): + with self.assertRaisesRegex(ValueError,'all-blocked|base_block requires fine_block'): + module.FixedPatchStepper(cfg,(4,4,20,28),context=None,**kwargs) + module.compiler.compile_kernel.assert_not_called() + + def test_uniform_block_codec_and_raw_snapshot(self): + path = Path(backend.r.__file__).with_name('block_stepper.py') + tree = ast.parse(path.read_text()) + tree.body = [n for n in tree.body if not ( + isinstance(n, ast.Import) and any(a.name == 'pycuda.driver' for a in n.names) + or isinstance(n, ast.ImportFrom) and n.module == 'cuda')] + module = types.ModuleType('CelerisLab.lbm._test_block_host') + module.__package__ = 'CelerisLab.lbm' + cuda, compiler = Mock(), Mock() + module.__dict__.update(cuda=cuda,compiler=compiler) + exec(compile(tree,str(path),'exec'),module.__dict__) + cuda.mem_get_info.return_value=(10**9,10**9) + cuda.mem_alloc.side_effect=lambda size: Mock() + gpu=Mock();gpu.get_global.return_value=(Mock(),20) + compiler.load_module.return_value=gpu + context=Mock(activate=nullcontext,sm_arch='sm_70') + layout=module.BlockLayout2D.build(16,16,4,4) + cfg=replace(backend.CFG,nx=16,ny=16,ddf_shifting=True) + run=module.PeriodicBlockStepper(cfg,layout,context=context) + storage={} + cuda.memcpy_htod.side_effect=lambda p,a: storage.__setitem__(id(p),np.array(a,copy=True)) + cuda.memcpy_dtoh.side_effect=lambda a,p: np.copyto(a,storage[id(p)]) + try: + physical=np.broadcast_to(_WEIGHTS,(9,16,16)).copy();physical[1]+=2**-35 + run.initialize(physical) + raw=run._download_raw() + self.assertTrue(raw.shifted) + self.assertEqual(len(raw.values),1) + self.assertTrue((raw.values[0][1]==np.float32(2**-35)).all()) + np.testing.assert_array_equal(run.download(),physical) + finally: run.close() + + def test_passive_shifted_body_raw_remap_preserves_support(self): + from tests.unit.test_pinned_body_remap import mesh, body, A as rects, B as moved + old,new=mesh(rects,(body(),)),mesh(moved,(body(),)) + result=remap(old,new,shifted_fields(old)) + self.assertTrue(result.shifted) + self.assertTrue(np.isnan(result.g[2][:,~new.leaf_masks[2]]).all()) + changed=mesh(moved,(body(center=(12.1,14.)),)) + with self.assertRaisesRegex(ValueError,'boundary'): + remap(old,changed,shifted_fields(old)) + + def test_fp64_physical_roundtrip_cannot_replace_raw(self): + tiny=np.full((9,1,1),np.float32(2**-80)) + self.assertFalse(np.array_equal(tiny,_encode_physical(_decode_stored(tiny,True),True))) + + +if __name__=='__main__': unittest.main() diff --git a/tests/unit/test_synchronized_remap.py b/tests/unit/test_synchronized_remap.py index 826600f..9f924b3 100644 --- a/tests/unit/test_synchronized_remap.py +++ b/tests/unit/test_synchronized_remap.py @@ -179,7 +179,7 @@ class SynchronizedRemapTests(unittest.TestCase): with self.assertRaises(ValueError): topology(rectangles) parent = self.old.parent - for changes in (dict(lattice_model='D3Q19'), dict(ddf_shifting=True), + for changes in (dict(lattice_model='D3Q19'), dict(collision='TRT'), dict(les_enabled=True), dict(nz=2)): with self.assertRaises(ValueError): RemapTopology.from_config(replace(CFG, **changes), parent, A)