From 9dddf5d8121d500aba2105f0642094adbfeca27f Mon Sep 17 00:00:00 2001 From: Codex Date: Sat, 17 Oct 2026 00:46:49 -0700 Subject: [PATCH] Gate shifted velocity Jets to all-blocked callers --- src/CelerisLab/body/jets.py | 7 +- src/CelerisLab/lbm/refinement.py | 4 +- src/CelerisLab/lbm/refinement_stepper.py | 18 +- src/CelerisLab/lbm/refinement_tables.py | 4 +- src/CelerisLab/lbm/synchronized_remap.py | 11 +- tests/unit/test_shifted_jets.py | 105 +++++++++ .../multigrid/test_shifted_jets_gpu.py | 208 ++++++++++++++++++ 7 files changed, 340 insertions(+), 17 deletions(-) create mode 100644 tests/unit/test_shifted_jets.py create mode 100644 tests/validation/multigrid/test_shifted_jets_gpu.py diff --git a/src/CelerisLab/body/jets.py b/src/CelerisLab/body/jets.py index 067c1c3..a967da7 100644 --- a/src/CelerisLab/body/jets.py +++ b/src/CelerisLab/body/jets.py @@ -189,13 +189,14 @@ def edge_to_node_frame(points): def surface_link_velocity(objects, cfg, flags, fluid_idx, dir, body_id, rx, ry, - jets=(), hits=None): + jets=(), hits=None, *, _shifted_blocked=False): """Per-link prescribed surface velocity (n, 2) FP32 and aperture mask. Host-only and storage-neutral: rows follow the given physical cut-link order. Returns (None, None) when neither an object carries surface actuation nor a jet is given. Jets need true FP64 ``hits`` (n, 2) and may not share a link with each other or with a legacy aperture. + Shifted storage is privately enabled only by the all-blocked stepper. """ enabled = [obj for obj in objects if obj.surface_actuation is not None] if not enabled and not jets: @@ -219,7 +220,9 @@ def surface_link_velocity(objects, cfg, flags, fluid_idx, dir, body_id, rx, ry, if hits is None or len(hits) != len(dir): raise ValueError("Jets require true hit points row-aligned with the cut links") if (cfg is None or not cfg.is_d2q9 or cfg.store_precision != "FP32" - or cfg.streaming != "double_buffer" or cfg.ddf_shifting): + or cfg.streaming != "double_buffer" + or (cfg.ddf_shifting and not (_shifted_blocked is True + and cfg.collision == "SRT" and not cfg.les_enabled))): raise ValueError("Jets require D2Q9/FP32/unshifted/double_buffer") shape, jet_id = jet_link_shapes(jets, hits) driven = jet_id >= 0 diff --git a/src/CelerisLab/lbm/refinement.py b/src/CelerisLab/lbm/refinement.py index 78b6a31..ab9a61c 100644 --- a/src/CelerisLab/lbm/refinement.py +++ b/src/CelerisLab/lbm/refinement.py @@ -422,7 +422,7 @@ def _circle_in_coverage(covered, cx, cy, radius): def prepare_fine_boundary(cfg, patch, objects, *, wall_aware=False, jets=(), - parent_origin=(0., 0.), parent_spacing=1.): + parent_origin=(0., 0.), parent_spacing=1., _shifted_blocked=False): """Project fixed circles and reject boundary/transfer overlap before CUDA. Input centers/radii and Jets are in the root edge frame. ``parent_origin`` @@ -555,7 +555,7 @@ def prepare_fine_boundary(cfg, patch, objects, *, wall_aware=False, jets=(), jet_cfg.validate() fluid_idx, dirs, _, body_id, rx, ry, _, _ = packed surface_link_velocity(local, jet_cfg, flags, fluid_idx, dirs, body_id, rx, ry, - jets=jets, hits=hits) + jets=jets, hits=hits, _shifted_blocked=_shifted_blocked) return FineBoundary2D(tuple(local), _readonly(flags, np.uint16), tuple(_readonly(a, a.dtype) for a in packed), _readonly(sorted(boundary), np.uint32), diff --git a/src/CelerisLab/lbm/refinement_stepper.py b/src/CelerisLab/lbm/refinement_stepper.py index ef98b24..0843349 100644 --- a/src/CelerisLab/lbm/refinement_stepper.py +++ b/src/CelerisLab/lbm/refinement_stepper.py @@ -287,11 +287,12 @@ class FixedPatchStepper: pair = pairs[level - 1] if level == 1: host = prepare_fine_boundary(cfg, pair, level_objects[1], - wall_aware=self._wall_aware, jets=level_jets[1]) + wall_aware=self._wall_aware, jets=level_jets[1], + _shifted_blocked=self._ddf_shifting) else: host = prepare_fine_boundary(self._level_cfgs[level - 1], pair, level_objects[level], jets=level_jets[level], parent_origin=origin, - parent_spacing=spacing) + parent_spacing=spacing, _shifted_blocked=self._ddf_shifting) self._boundary_hosts.append(host) origin = tuple(o + spacing * f for o, f in zip(origin, pair.fine_origin)) spacing /= 2 @@ -309,7 +310,8 @@ class FixedPatchStepper: self._level_jets = level_jets self._level_boundaries = [None] + [ _prepare_blocked_boundary(self._boundary_hosts[level], self._layouts[level], - self._level_cfgs[level], level_jets[level]) + self._level_cfgs[level], level_jets[level], + _shifted_blocked=self._ddf_shifting) for level in range(1, depth)] self._leaf = [kind == ACTIVE for kind in self.hierarchy.kinds] for level in range(1, depth): @@ -558,6 +560,7 @@ class FixedPatchStepper: with self._lock: self._require("ready") + self._check_storage() pairs = self.hierarchy.pairs if (self._wall_aware or self.hierarchy.depth != 3 or not isinstance(pairs[1], FixedPatchSet2D) @@ -569,19 +572,22 @@ class FixedPatchStepper: actual = self._level_boundaries[2] jets = () if actual is None else actual.jets old = RemapTopology.from_config(cfg, pairs[0], pairs[1].bounds, - child_objects=child_objects, child_jets=jets) + child_objects=child_objects, child_jets=jets, + _shifted_blocked=self._ddf_shifting) if (not same_physical_boundary(old.boundary, self._boundary_hosts[2]) or (self._boundary_hosts[2] is not None and actual is None)): raise ValueError("child_objects must reproduce the original root-edge-frame boundary") if actual is not None: velocity = actual.tables.get("surface_velocity") if velocity is None: - velocity = boundary_velocity(cfg, pairs[1], actual.host, jets) + velocity = boundary_velocity(cfg, pairs[1], actual.host, jets, + _shifted_blocked=self._ddf_shifting) if (not same_physical_boundary(old.boundary, actual.host) or not np.array_equal(old.surface_velocity, velocity)): raise ValueError("child_objects changed boundary ownership or surface velocity") new = RemapTopology.from_config(cfg, pairs[0], rectangles, - child_objects=child_objects, child_jets=jets) + child_objects=child_objects, child_jets=jets, + _shifted_blocked=self._ddf_shifting) if (not same_physical_boundary(old.boundary, new.boundary) or not np.array_equal(old.surface_velocity, new.surface_velocity)): raise ValueError("remap requires the identical pinned physical boundary and BFL support") diff --git a/src/CelerisLab/lbm/refinement_tables.py b/src/CelerisLab/lbm/refinement_tables.py index 65acd3a..b097f30 100644 --- a/src/CelerisLab/lbm/refinement_tables.py +++ b/src/CelerisLab/lbm/refinement_tables.py @@ -87,7 +87,7 @@ class _BlockedBoundary: accumulate: object = None -def _prepare_blocked_boundary(boundary, layout, cfg, jets): +def _prepare_blocked_boundary(boundary, layout, cfg, jets, *, _shifted_blocked=False): if boundary is None or layout is None: return None fluid_idx, dirs, q, body_id, rx, ry, _, fallback = boundary.curved_host @@ -102,6 +102,6 @@ def _prepare_blocked_boundary(boundary, layout, cfg, jets): cfg.validate() tables["surface_velocity"], tables["aperture"] = surface_link_velocity( boundary.objects, cfg, boundary.flags, fluid_idx, dirs, body_id, rx, ry, - jets=jets, hits=boundary.curved_hits) + jets=jets, hits=boundary.curved_hits, _shifted_blocked=_shifted_blocked) return _BlockedBoundary(boundary, tables, jets, actuated, jet_link_shapes(jets, boundary.curved_hits) if jets else None) diff --git a/src/CelerisLab/lbm/synchronized_remap.py b/src/CelerisLab/lbm/synchronized_remap.py index d280e28..c66d006 100644 --- a/src/CelerisLab/lbm/synchronized_remap.py +++ b/src/CelerisLab/lbm/synchronized_remap.py @@ -39,7 +39,7 @@ class RemapTopology: @classmethod def from_config(cls, cfg, parent, child_rectangles, *, coarse_boundary="periodic", objects=(), - child_objects=(), child_jets=()): + child_objects=(), child_jets=(), _shifted_blocked=False): if coarse_boundary != "periodic" or tuple(objects) or cfg.nz != 1: raise ValueError("remap requires periodic 2D state; use child_objects for a pinned L2 body") if type(parent) not in (FixedPatch2D, FixedPatchSet2D): @@ -51,7 +51,7 @@ class RemapTopology: viscosity=2 * cfg.viscosity) boundary = prepare_fine_boundary( fine_cfg, child, child_objects, jets=child_jets, - parent_origin=parent.fine_origin, parent_spacing=.5) + parent_origin=parent.fine_origin, parent_spacing=.5, _shifted_blocked=_shifted_blocked) pinned, velocity = (), None child_leaf = child.fine_kind == ACTIVE if boundary is not None: @@ -67,7 +67,8 @@ class RemapTopology: if match is not None: containing.add(match) pinned = tuple(sorted(containing)) - velocity = boundary_velocity(cfg, child, boundary, child_jets) + velocity = boundary_velocity(cfg, child, boundary, child_jets, + _shifted_blocked=_shifted_blocked) if velocity is not None: velocity.setflags(write=False) masks = (~parent.covered, middle_kind == ACTIVE, child_leaf) @@ -76,7 +77,7 @@ class RemapTopology: return cls(parent, child, masks, boundary, pinned, velocity) -def boundary_velocity(cfg, child, boundary, jets=()): +def boundary_velocity(cfg, child, boundary, jets=(), *, _shifted_blocked=False): """Physical cut-link order; never blocked addresses or stale Jet snapshots.""" local_cfg = replace(cfg, nx=child.fine_shape[1], ny=child.fine_shape[0], viscosity=4 * cfg.viscosity) @@ -84,7 +85,7 @@ def boundary_velocity(cfg, child, boundary, jets=()): fluid, dirs, _, ids, rx, ry, _, _ = boundary.curved_host velocity, _ = surface_link_velocity( boundary.objects, local_cfg, boundary.flags, fluid, dirs, ids, rx, ry, - jets=jets, hits=boundary.curved_hits) + jets=jets, hits=boundary.curved_hits, _shifted_blocked=_shifted_blocked) return velocity diff --git a/tests/unit/test_shifted_jets.py b/tests/unit/test_shifted_jets.py new file mode 100644 index 0000000..2b56794 --- /dev/null +++ b/tests/unit/test_shifted_jets.py @@ -0,0 +1,105 @@ +"""Shifted velocity Jets: bounded opt-in and host lifecycle, no CUDA numerics.""" +from dataclasses import replace +import unittest +from unittest.mock import patch +import numpy as np +from CelerisLab.body.jets import Jet, Tube, Parallel +from CelerisLab.body.objects import SurfaceActuation +from CelerisLab.lbm.synchronized_remap import RemapTopology +from CelerisLab.lbm.refinement import FixedPatchSet2D +from CelerisLab.lbm.population_storage import _StoredPopulations +from tests.unit import test_refinement_three_level_multipatch as backend +from tests.unit.test_pinned_body_remap import A, B, body, jet, mesh +from tests.unit.test_shifted_storage import shifted_fields + +CFG=replace(backend.CFG,ddf_shifting=True) + + +def topology(objects, jets, *, enabled=False, cfg=CFG): + return RemapTopology.from_config(cfg,FixedPatchSet2D.from_config(cfg,backend.PARENTS),A, + child_objects=objects,child_jets=jets,_shifted_blocked=enabled) + + +def ready(): + objects,jets=(body(),),(jet(),) + with patch.object(backend,'CFG',CFG),patch.object(backend,'bodies_and_jets',return_value=(objects,jets)): + run,module,allocations=backend.mock_backend(children=A,bodies=True) + 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(mesh(A,objects,jets)),True)) + run._step_count=7 + return run,module,allocations,storage,objects + + +class ShiftedJets(unittest.TestCase): + def test_default_guard_and_storage_neutral_host_rows(self): + with self.assertRaisesRegex(ValueError,'unshifted'): + topology((body(),),(jet(),)) + shifted=topology((body(),),(jet(),),enabled=True) + unshifted=topology((body(),),(jet(),),cfg=backend.CFG) + np.testing.assert_array_equal(shifted.surface_velocity,unshifted.surface_velocity) + self.assertEqual(shifted.surface_velocity.dtype,np.float32) + legacy=body(SurfaceActuation(.7,.85,.002)) + with self.assertRaisesRegex(ValueError,'unshifted'): + topology((legacy,),(),enabled=True) + + def test_rejections_before_device_allocation(self): + module=backend.host_module() + moving=body();moving.state.omega=.01 + cases=[((moving,),(jet(),)),((body(),),(jet(),jet())), + ((body(),),(Jet(Tube((100.,100.),0.,1.,.1),Parallel(0.)),)), + ((body(center=(8.,14.)),),(jet(),)), + ((body(),body(center=(34.,14.))), + (Jet(Tube((10.,14.),0.,26.,.3),Parallel(0.)),))] + for objects,jets in cases: + with self.subTest(objects=len(objects),jets=len(jets)),self.assertRaises(ValueError): + module.FixedPatchStepper(CFG,backend.PARENTS,context=None, + base_block=(4,4),fine_block=(4,4),child_block=(4,4), + child_rectangles_in_parent=A,child_objects=objects,child_jets=jets) + for blocks in ((None,None,None),((4,4),None,(4,4))): + with self.assertRaises(ValueError): + module.FixedPatchStepper(CFG,backend.PARENTS,context=None, + base_block=blocks[0],fine_block=blocks[1],child_block=blocks[2], + child_rectangles_in_parent=A,child_objects=(body(),),child_jets=(jet(),)) + module.cuda.mem_alloc.assert_not_called() + module.compiler.compile_kernel.assert_not_called() + + def test_switch_upload_only_then_raw_remap_keeps_live_velocity(self): + run,module,allocations,storage,objects=ready();clone=None + try: + boundary=run._level_boundaries[2] + before=run._download_raw() + original=boundary.tables['surface_velocity'].copy() + module.cuda.memcpy_htod.reset_mock() + run.set_jet_amplitude(0,-jet().amplitude,level=2) + module.cuda.memcpy_htod.assert_called_once() + self.assertIs(module.cuda.memcpy_htod.call_args.args[0],boundary.device['surface_velocity']) + np.testing.assert_array_equal(boundary.tables['surface_velocity'],-original) + self.assertEqual(run.step_count,7) + for a,b in zip(before.values,run._download_raw().values):self.assertEqual(a.tobytes(),b.tobytes()) + clone=run.remap_child(B,child_objects=objects) + self.assertEqual(clone._level_jets[2][0].amplitude,-jet().amplitude) + np.testing.assert_array_equal(clone._level_boundaries[2].tables['surface_velocity'], + boundary.tables['surface_velocity']) + self.assertEqual(clone.step_count,7) + saved=run._download_raw() + with self.assertRaises(ValueError):run.remap_child(B,child_objects=(body(center=(12.1,14.)),)) + for a,b in zip(saved.values,run._download_raw().values):self.assertEqual(a.tobytes(),b.tobytes()) + finally: + if clone is not None:clone.close() + run.close() + + def test_invalid_amplitude_does_not_upload_or_commit(self): + run,module,_,_,_=ready() + try: + before=run._level_boundaries[2].tables['surface_velocity'].copy() + module.cuda.memcpy_htod.reset_mock() + for amplitude in (float('nan'),float('inf'),1.): + with self.assertRaises(ValueError):run.set_jet_amplitude(0,amplitude,level=2) + module.cuda.memcpy_htod.assert_not_called() + np.testing.assert_array_equal(before,run._level_boundaries[2].tables['surface_velocity']) + finally:run.close() + + +if __name__=='__main__':unittest.main() diff --git a/tests/validation/multigrid/test_shifted_jets_gpu.py b/tests/validation/multigrid/test_shifted_jets_gpu.py new file mode 100644 index 0000000..48947b3 --- /dev/null +++ b/tests/validation/multigrid/test_shifted_jets_gpu.py @@ -0,0 +1,208 @@ +"""Shifted velocity-Jet GPU gates. Run only on an authorized device. + +PYTHONPATH=src python tests/validation/multigrid/test_shifted_jets_gpu.py -v +Uniform initial curved-body states rho=.99/1.01; amplitude 0,+.002,-.001,0; +four-level substep ledger; three-level live-Jet raw remap. No pressure actuator. +The census is an implementation-exchange ledger, not an independent physical +traction/nozzle-thrust model. Frozen inherited caps: raw rows rtol=2e-6/atol=2e-7; fluid ledger +rtol=2e-4/atol=8e-5. Observer on/off, interval accumulation and raw copies exact. +Report actual residuals too: passing old loose caps alone is not an accuracy gain. +""" +from contextlib import ExitStack +import unittest +import numpy as np +import pycuda.driver as cuda +from CelerisLab.config import LBMConfig +from CelerisLab.cuda.context import CudaContext +from CelerisLab.body.geometry.circle import CircleGeometry +from CelerisLab.body.objects import SimObject +from CelerisLab.body.jets import Jet, Tube, Parallel +from CelerisLab.lbm.descriptors import D2Q9_EX, D2Q9_EY, SOLID +from CelerisLab.lbm.refinement_stepper import FixedPatchStepper +from CelerisLab.lbm.population_storage import _StoredPopulations +from CelerisLab.lbm.synchronized_remap import RemapTopology, _remap_stored + +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.float64) +SPEEDS=(0.,.002,-.001,0.) + + +def circle(center,radius): + return SimObject(17,geometry=CircleGeometry(*center,radius),center=center,radius=radius) + + +def drive(center,radius): + return Jet(Tube(center,0.,radius*1.5,radius*.32),Parallel(0.),amplitude=0.) + + +def initial(run,rho,u=.013,v=-.007): + cu=3*(CX*u+CY*v) + g=W*((rho-1)+rho*(cu+.5*cu*cu-1.5*(u*u+v*v))) + return _StoredPopulations(tuple(np.broadcast_to(g[:,None,None],l.shape).astype(np.float32).copy() + for l in run._levels),True) + + +def download(pointer,shape,dtype=np.float32): + out=np.empty(shape,dtype);cuda.memcpy_dtoh(out,pointer);return out + + +def raw_totals(a,leaf): + g=a[:,leaf].astype(np.float64) + # Only differences at the same fixed leaf mask are used, so mass's unit + # background cancels. Momentum never requires adding weights. + return np.array([(g*CX[:,None]).sum(),(g*CY[:,None]).sum(),g.sum()]) + + +def census(run,index,buffer): + """Independent physical link/face census of actual staged receiver windows.""" + level=run._levels[index];layout=level.layout;b=level.boundary + stored=download(buffer,layout.blocked_shape);dense=layout.to_dense(stored) + 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=stored.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;dy,dx=y+cy,x+cx + incoming=flat[block,q,cell-cy*layout.local_shape[1]-cx].astype(np.float64)+W[q] + outgoing=flat[block,q,cell].astype(np.float64)+W[q] + entering=~leaf[sy,sx]&~solid[sy,sx];leaving=~leaf[dy,dx]&~solid[dy,dx] + external+=np.array([cx,cy,1])*(incoming[entering].sum()-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) + outgoing=flat[block,direction,cell].astype(np.float64) + source=cell+CY[direction]*layout.local_shape[1]+CX[direction] + reflected=flat[block,OPP[direction],source].astype(np.float64) + physical_sum=outgoing+reflected+2*W[direction] + jx,jy=CX[direction]*physical_sum,CY[direction]*physical_sum + rows=np.column_stack((jx,jy,rx*jy-ry*jx,reflected-outgoing)) + return raw_totals(dense,leaf),external,rows + + +class Ledger: + def __init__(self,test,run,index): + self.test,self.run,self.index=test,run,index + self.samples=[];self.max_row=0.;self.max_balance=0. + level=run._levels[index];self.original=level.kernel;b=level.boundary + def audited(src,dst,*args,**kwargs): + run._stream.synchronize() + before,external,independent=census(run,index,src) + raw=download(b.last_raw,(len(b.tables['cell']),4)) + np.testing.assert_allclose(raw,independent,rtol=2e-6,atol=2e-7) + self.max_row=max(self.max_row,float(abs(raw.astype(float)-independent).max())) + self.original(src,dst,*args,**kwargs);run._stream.synchronize() + after=level.layout.to_dense(download(dst,level.layout.blocked_shape)) + change=raw_totals(after,run._leaf[index])-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) + self.max_balance=max(self.max_balance,float(abs(change-expected).max())) + self.samples.append(raw.copy()) + level.kernel=audited + + def step(self): + run=self.run;self.samples.clear();run.step() + b=run._level_boundaries[self.index];h=run._exchange_spacing[self.index] + self.test.assertEqual(len(self.samples),2**self.index) + expected=np.zeros_like(self.samples[0],dtype=np.float64) + weights=np.array([h*h,h*h,h*h*h,h*h]) + for raw in self.samples:expected+=raw.astype(np.float64)*weights + actual=download(b.interval_integral,expected.shape,np.float64) + np.testing.assert_array_equal(actual,expected) + obs=run.read_raw_exchange(0,level=self.index) + for name in ('impulse','angular_impulse','mass_into_fluid'): + np.testing.assert_allclose(obs['total'][name],np.asarray(obs['wall'][name])+obs['aperture'][name],rtol=0,atol=2e-13) + self.test.assertEqual(obs['sample_count'],2**self.index) + for a,leaf in zip(run._download_raw().values,run._leaf): + self.test.assertTrue(np.isfinite(a[:,leaf]).all()) + self.test.assertTrue((1+a[:,leaf].sum(axis=0,dtype=np.float64)>0).all()) + return obs + + def report(self,label): + print(f'{label}: max_abs_raw_row={self.max_row:.9g}, max_abs_fluid_ledger={self.max_balance:.9g}',flush=True) + + +class ShiftedJetsGPU(unittest.TestCase): + def exact_raw(self,a,b): + self.assertTrue(a.shifted and b.shifted) + for x,y in zip(a.values,b.values):self.assertEqual(x.tobytes(),y.tobytes()) + + def test_matching_whole_circle_parallel_velocity_preserves_uniform_state(self): + cfg=LBMConfig(nx=32,ny=24,viscosity=.1,ddf_shifting=True) + center=(16.,12.);radius=2.15;obj=circle(center,radius) + # Rectangle contains every circle hit; prescribed surface velocity equals + # the uniform fluid velocity. Geometry itself remains fixed. + whole=Jet(Tube((center[0]-radius-1,center[1]),0.,2*radius+2,radius+1), + Parallel(0.),amplitude=.013) + for rho in (.99,1.01): + with self.subTest(rho=rho),CudaContext(0) as ctx: + with FixedPatchStepper(cfg,(8,6,24,18),context=ctx,base_block=(4,4),fine_block=(4,4), + fine_objects=(obj,),jets=(whole,),observe_raw_exchange=True) as run: + run._initialize_raw(initial(run,rho,u=.013,v=0.)) + expected=run._download_raw();ledger=Ledger(self,run,1) + self.assertTrue(run._level_boundaries[1].tables['aperture'].all()) + for _ in range(4):ledger.step() + error=max(float(abs(a[:,leaf]-b[:,leaf]).max()) for a,b,leaf in + zip(run._download_raw().values,expected.values,run._leaf)) + self.assertLessEqual(error,1e-6) # inherited uniform-state cap + print(f'matched whole-circle rho={rho}: max_abs_raw_state={error:.9g}',flush=True) + ledger.report(f'matched whole-circle rho={rho}') + + def test_uniform_initial_curved_rho_both_signs_switch_and_observer(self): + cfg=LBMConfig(nx=32,ny=24,viscosity=.1,ddf_shifting=True) + obj=circle((16.,12.),2.15);j=drive((16.,12.),2.15) + for rho in (.99,1.01): + with self.subTest(rho=rho),CudaContext(0) as ctx,ExitStack() as owned: + settings=dict(bounds=(8,6,24,18),context=ctx,base_block=(4,4),fine_block=(4,4),fine_objects=(obj,),jets=(j,)) + measured=owned.enter_context(FixedPatchStepper(cfg,observe_raw_exchange=True,**settings)) + plain=owned.enter_context(FixedPatchStepper(cfg,**settings)) + measured._initialize_raw(initial(measured,rho));plain._initialize_raw(initial(plain,rho)) + ledger=Ledger(self,measured,1) + for speed in SPEEDS: + before=measured._download_raw();count=measured.step_count + measured.set_jet_amplitude(0,speed);plain.set_jet_amplitude(0,speed) + self.exact_raw(before,measured._download_raw());self.assertEqual(measured.step_count,count) + ledger.step();plain.step() + self.exact_raw(measured._download_raw(),plain._download_raw()) + ledger.report(f'L1 rho={rho}') + + def test_four_level_interval_ledger(self): + cfg=LBMConfig(nx=48,ny=48,viscosity=.02,ddf_shifting=True) + obj=circle((16.,16.),.5);j=drive((16.,16.),.5) + with CudaContext(0) as ctx,ExitStack() as owned: + settings=dict(context=ctx,refinements=[[(6,6,26,26)],[(20,20,44,44)],[(56,56,72,72)]], + level_blocks=[(4,4)]*4,level_objects=[(),(),(),(obj,)],level_jets=[(),(),(),(j,)]) + measured=owned.enter_context(FixedPatchStepper(cfg,observe_raw_exchange=True,**settings)) + plain=owned.enter_context(FixedPatchStepper(cfg,**settings)) + measured._initialize_raw(initial(measured,1.01));plain._initialize_raw(initial(plain,1.01)) + ledger=Ledger(self,measured,3) + for speed in SPEEDS: + measured.set_jet_amplitude(0,speed,level=3);plain.set_jet_amplitude(0,speed,level=3) + ledger.step();plain.step();self.exact_raw(measured._download_raw(),plain._download_raw()) + ledger.report('L3 rho=1.01') + + def test_live_jet_remap_then_both_ledgers_continue(self): + cfg=LBMConfig(nx=48,ny=32,viscosity=.1,ddf_shifting=True) + parents=((4,4,20,28),(24,4,44,28));a=((16,16,32,40),(56,16,72,40));b=(a[0],(60,16,76,40)) + obj=circle((12.,14.),.8);j=drive((12.,14.),.8) + with CudaContext(0) as ctx,ExitStack() as owned: + old=owned.enter_context(FixedPatchStepper(cfg,parents,context=ctx,base_block=(4,4),fine_block=(4,4), + child_block=(4,4),child_rectangles_in_parent=a,child_objects=(obj,),child_jets=(j,),observe_raw_exchange=True)) + old._initialize_raw(initial(old,.99));ledger=Ledger(self,old,2) + for speed in (.002,-.001):old.set_jet_amplitude(0,speed,level=2);ledger.step() + before=old._download_raw();count=old.step_count + topology=lambda rects:RemapTopology.from_config(cfg,old.patch,rects,child_objects=(obj,), + child_jets=old._level_jets[2],_shifted_blocked=True) + expected=_remap_stored(topology(a),topology(b),before,shifted=True,base_count=count, + level_steps=(count,2*count,4*count),phase='postcollision_g') + new=owned.enter_context(old.remap_child(b,child_objects=(obj,))) + self.exact_raw(before,old._download_raw()) + self.exact_raw(_StoredPopulations(expected.g,True),new._download_raw()) + self.assertEqual(new._level_jets[2][0].amplitude,-.001) + np.testing.assert_array_equal(new._level_boundaries[2].tables['surface_velocity'],old._level_boundaries[2].tables['surface_velocity']) + newledger=Ledger(self,new,2);ledger.step();newledger.step() + self.assertEqual((old.step_count,new.step_count),(count+1,count+1)) + ledger.report('L2 source');newledger.report('L2 remapped') + + +if __name__=='__main__':unittest.main()