Compare commits

...
Author SHA1 Message Date
Frank14f 7b652b2423 test: construct shear remap fixture through per-level API 2026-10-28 10:21:14 +08:00
Frank14f 7482d95a42 Update channel CPU shim include and retain frozen projected driver 2026-10-28 09:57:15 +08:00
Frank14f 55067652e4 Migrate current callers to level-indexed refinement API 2026-10-28 09:53:20 +08:00
ClaudeandClaude Opus 5.5 b4eacf7526 FixedPatchStepper: level-indexed API, drop 15 fine/child compatibility views
- Hierarchy documents level/pair/shape/frame conventions; edge_frame(l)
  replaces the stepper's inline origin/spacing loop.
- Public per-level read-only state: layouts, boundaries, jets, leaf.
- set_jet_amplitude/read_raw_exchange/read_exchange_window: level= required.
- remap_child(child_objects=) -> remap_finest(objects=).
- _kernel_configs, _upload_boundary, _observed_boundary/_summarize split
  construction and observation duplication; behavior unchanged.
- fine_level_origin -> fine_level_offset (a (y, x) index, not an origin).
- docs: migration map in docs/development/src_coherence.md.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
2026-10-28 01:39:48 +00:00
ClaudeandClaude Opus 5.5 0c64f821a4 Blocked steps: one block_collide extension point; drop unused HaloPlan, collide_trt_no_force, MIN_PATCH_GAP
SRT keeps its two explicit arithmetic variants (level: unclamped no-force,
shifted direct-momentum feq; channel: uniform OneStep arithmetic), selected
at compile time by BLOCK_ROOT_BC; text moved verbatim, PTX parity pending.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
2026-10-28 01:37:38 +00:00
Frank14f 914906efd4 Merge periodic stepper alias cleanup into final stage candidate 2026-10-28 08:30:35 +08:00
Frank14f 55b0b34fc5 Merge direct Simulation Jet GPU regression 2026-10-28 08:29:52 +08:00
ClaudeandClaude Opus 5.5 ddc4a1a6c4 lifecycle: guard Simulation.stream after a failed jet update; scope docs
A set_jet_amplitude failure while creating the internal stream marked the
Simulation failed, but the public stream property still retried
cuda.Stream(). Guard it like the other device-touching methods and add a
regression test (no retry, close still retryable). Lifecycle wording now
names the guarded device-touching methods; cached readers and host staging
stay callable. compatibility.md records the caller-visible fail-closed change.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
2026-10-28 00:23:56 +00:00
Frank14f a9aa12e65d refactor: remove redundant periodic patch stepper alias 2026-10-28 08:03:29 +08:00
53 changed files with 510 additions and 408 deletions
+7 -7
View File
@@ -117,8 +117,8 @@ with CudaContext(0) as ctx, FixedPatchStepper(cfg, context=ctx, refinements=[(8,
level_objects=[(), (body,)]) as run: level_objects=[(), (body,)]) as run:
rest = np.asarray(compute_feq_d2q9(1.0, 0.0, 0.0), np.float32)[:, None, None] rest = np.asarray(compute_feq_d2q9(1.0, 0.0, 0.0), np.float32)[:, None, None]
# (9, ny, nx) postcollision populations per level; non-leaf cells are discarded. # (9, ny, nx) postcollision populations per level; non-leaf cells are discarded.
run.initialize(np.broadcast_to(rest, (9, *run.patch.coarse_shape)), run.initialize(np.broadcast_to(rest, (9, *run.hierarchy.pairs[0].coarse_shape)),
np.broadcast_to(rest, (9, *run.patch.fine_shape))) np.broadcast_to(rest, (9, *run.hierarchy.pairs[0].fine_shape)))
run.step(128) # coarse steps; fine level takes two per step run.step(128) # coarse steps; fine level takes two per step
coarse_g, fine_g = run.download() # NaN on covered, ghost and solid cells coarse_g, fine_g = run.download() # NaN on covered, ghost and solid cells
``` ```
@@ -143,7 +143,7 @@ Scope and limits:
from their union; adjacent, overlapping and concave unions are accepted when from their union; adjacent, overlapping and concave unions are accepted when
their actual dependencies close. There is no fixed rectangle-gap rule. their actual dependencies close. There is no fixed rectangle-gap rule.
The fine level is then one level-wide grid `(2*ny, 2*nx)`: The fine level is then one level-wide grid `(2*ny, 2*nx)`:
`run.patch.fine_shape` gives it, `initialize`/`download` use it, and only the `run.hierarchy.pairs[0].fine_shape` gives it, `initialize`/`download` use it, and only the
grid-aligned tiles intersecting the union plus its four-cell ghost band are stored. grid-aligned tiles intersecting the union plus its four-cell ghost band are stored.
Block sizes must divide the level-wide grid; rectangle edges need not align Block sizes must divide the level-wide grid; rectangle edges need not align
with tiles. Order of rectangles does not matter. L1 bodies/Jets (`level_objects[1]`/`level_jets[1]`) work with tiles. Order of rectangles does not matter. L1 bodies/Jets (`level_objects[1]`/`level_jets[1]`) work
@@ -159,8 +159,8 @@ Scope and limits:
with CudaContext(0) as ctx, FixedPatchStepper( with CudaContext(0) as ctx, FixedPatchStepper(
cfg, context=ctx, refinements=[((4, 4, 12, 12), (16, 4, 24, 12))], cfg, context=ctx, refinements=[((4, 4, 12, 12), (16, 4, 24, 12))],
level_blocks=[(4, 4), (4, 4)]) as run: level_blocks=[(4, 4), (4, 4)]) as run:
run.initialize(np.broadcast_to(rest, (9, *run.patch.coarse_shape)), run.initialize(np.broadcast_to(rest, (9, *run.hierarchy.pairs[0].coarse_shape)),
np.broadcast_to(rest, (9, *run.patch.fine_shape))) # (9, 32, 64) np.broadcast_to(rest, (9, *run.hierarchy.pairs[0].fine_shape))) # (9, 32, 64)
run.step(32) run.step(32)
coarse_g, fine_g = run.download() # coarse: NaN on every covered cell; coarse_g, fine_g = run.download() # coarse: NaN on every covered cell;
# fine: NaN on every non-ACTIVE cell # fine: NaN on every non-ACTIVE cell
@@ -193,10 +193,10 @@ Scope and limits:
cfg, context=ctx, refinements=[(4, 4, 12, 12), (8, 8, 16, 16)], cfg, context=ctx, refinements=[(4, 4, 12, 12), (8, 8, 16, 16)],
level_blocks=[(4, 4)] * 3, level_objects=[(), (), (body,)], level_blocks=[(4, 4)] * 3, level_objects=[(), (), (body,)],
level_jets=[(), (), (jet,)]) as run: level_jets=[(), (), (jet,)]) as run:
shapes = (run.patch.coarse_shape, run.patch.fine_shape, run.child_patch.fine_shape) shapes = (run.hierarchy.pairs[0].coarse_shape, run.hierarchy.pairs[0].fine_shape, run.hierarchy.pairs[1].fine_shape)
run.initialize(*(np.broadcast_to(rest, (9, *s)) for s in shapes)) # 16², 24², 24² run.initialize(*(np.broadcast_to(rest, (9, *s)) for s in shapes)) # 16², 24², 24²
run.step(4) run.step(4)
run.set_jet_amplitude(0, -.002, level=2) # default level=1 still addresses L1 run.set_jet_amplitude(0, -.002, level=2) # explicit level selects the driven level
run.step(4) run.step(4)
coarse_g, l1_g, l2_g = run.download() # solids/covered/nonleaves are NaN coarse_g, l1_g, l2_g = run.download() # solids/covered/nonleaves are NaN
``` ```
+22
View File
@@ -3,6 +3,13 @@
Changes that can make a previously accepted setup or call behave differently. Changes that can make a previously accepted setup or call behave differently.
Each entry names the commit, the affected entry points and the observable result. Each entry names the commit, the affected entry points and the observable result.
## Redundant periodic stepper alias removed
`PeriodicPatchStepper` is no longer exported from `lbm.refinement_stepper`.
Use `FixedPatchStepper` with the same arguments, including explicit
`level_blocks`. The old name was an alias for that exact class; periodic
boundaries remain the default. Existing imports of the old name must be updated.
## 87c3a31: cut-link ownership is validated for every body ## 87c3a31: cut-link ownership is validated for every body
`validate_cut_link_ownership` (`body/preprocess/flag_overlay.py`) now runs for `validate_cut_link_ownership` (`body/preprocess/flag_overlay.py`) now runs for
@@ -33,3 +40,18 @@ or negative amplitude (the aperture test is geometric). Bodies with neither
`SurfaceActuation` nor aperture links still raise `SurfaceActuation` nor aperture links still raise
`ValueError("Body has no surface actuation")`. The value is unchanged: the latest `ValueError("Body has no surface actuation")`. The value is unchanged: the latest
step's wall/aperture exchange, not total nozzle force. step's wall/aperture exchange, not total nozzle force.
## 7afb5ad: `Simulation` fails closed after a Jet amplitude device-write failure
`Simulation.set_jet_amplitude` previously only documented that a device error
leaves the per-link velocity column undefined; the instance stayed usable.
Now any error after input validation (internal stream creation or the upload)
marks the Simulation failed. Every later device-touching public method,
including the `stream` property (guard added in the follow-up repair commit on
`claude/jet-update-consistency`), raises
`RuntimeError("Simulation failed (...); close() and recreate it")`. The host
amplitude keeps its previous value. `close()` still works and stays retryable;
cached host readers and host staging calls are not blocked. Invalid input
(`ValueError`/`IndexError`) is unchanged: it raises before any change and the
instance stays usable. `FixedPatchStepper.set_jet_amplitude` already failed the
instance (status `failed`) and is unchanged.
+1 -1
View File
@@ -4,7 +4,7 @@
G-to-G explosion reconstruction. The omitted/default `"constant"` selector retains G-to-G explosion reconstruction. The omitted/default `"constant"` selector retains
the same kernel, launch arguments and device allocations. Projected mode adds two the same kernel, launch arguments and device allocations. Projected mode adds two
compact tables per adjacent pair; recursion, collision, BFL, restriction, buffer compact tables per adjacent pair; recursion, collision, BFL, restriction, buffer
ownership and synchronized publication are unchanged. `remap_child` preserves it. ownership and synchronized publication are unchanged. `remap_finest` preserves it.
Runtime eligibility is periodic-root, fully blocked D2Q9/SRT/FP32, body/Jet-free, Runtime eligibility is periodic-root, fully blocked D2Q9/SRT/FP32, body/Jet-free,
without LES; the fixed-patch backend has no forcing interface. Shifted and unshifted without LES; the fixed-patch backend has no forcing interface. Shifted and unshifted
+62
View File
@@ -0,0 +1,62 @@
# src coherence stage (from main 914906e)
Owner direction 2026-10-28: unified style over historical compatibility, inside `src`.
No numerical change is intended; kernel text was moved verbatim (PTX/bit parity is
pending a GPU run).
## Conventions (defined once, in `lbm.hierarchy.Hierarchy`)
- Level `l = 0 .. depth-1` (0 = root). Pair `hierarchy.pairs[l-1]` builds level `l`.
- Array shapes are `(ny, nx)`; rectangles, origins and centers are `(x, y)`.
- FixedPatchStepper bodies/Jets use the root EDGE frame (root cell `i` = `[i, i+1)`);
uniform `Simulation` objects use the integer-NODE frame. Nothing converts between them.
- `Hierarchy.edge_frame(l)` -> `((ox, oy), h)`: cell `(i, j)` of level `l`'s allocated
array (ghost band and holes included) covers `[o + h*i, o + h*(i+1))`. A scalar
FixedPatch2D window carries the half-coarse-cell ghost shift `bounds - HALO/2` in
`fine_origin`; a set's level-wide grid has origin 0. Level 1 boundaries are built from
the root config (the wall-aware fields live there); deeper ones from the parent frame.
- `FixedPatch2D.fine_level_offset` (was `fine_level_origin`) is an array INDEX offset
`(y, x)`, not an `(x, y)` origin.
## FixedPatchStepper migration map
| removed / changed | replacement |
|---|---|
| `patch`, `child_patch` | `hierarchy.pairs[0]`, `hierarchy.pairs[1]` |
| `base_layout`, `fine_layout`, `child_layout` (stepper properties) | `layouts[0]`, `layouts[1]`, `layouts[2]` (`FixedPatchSet2D.fine_layout(block)` is unchanged) |
| `fine_boundary`, `child_boundary` | `boundaries[1]`, `boundaries[2]` (host geometry; `None` = no bodies/Jets) |
| `_jets` | `jets[1]` |
| `_boundaries[l]` (dict from L1, level keys) | `_level_boundaries[l]` (list from L0, `None` at L0 unless single-level static circles) |
| `_fine_tables`, `_actuated`, `_jet_cache` | `_level_boundaries[1].tables` / `.actuated` / `.jet_cache` |
| `_fine_leaf`, `_child_leaf` | `leaf[1]`, `leaf[2]` (copies) |
| `_middle_kind` | `hierarchy.kinds[1]` |
| `_boundary_hosts` (private) | `_boundaries` (private list behind `boundaries`) |
| `remap_child(rects, child_objects=...)` | `remap_finest(rects, objects=...)` |
| `set_jet_amplitude(i, a)` (implicit level 1) | `set_jet_amplitude(i, a, level=l)` (required) |
| `read_raw_exchange(b)` / `read_exchange_window(b)` (implicit level 1) | `read_raw_exchange(b, level=l)` / `read_exchange_window(b, level=l)` (required) |
| `_level_name(l)` "base/fine/child" messages | `"level l"` |
| `refinement.MIN_PATCH_GAP` | deleted (no gap policy exists) |
| `params.cuh HaloPlan`, `collide_trt_no_force` | deleted (zero references) |
Kept: `step_count`, constructor keywords (`refinements`, `level_blocks`, `level_objects`,
`level_jets`, `coarse_boundary`, `explosion`), `Simulation.set_jet_amplitude(i, a)` (one level),
`RemapTopology.from_config(..., child_objects=, child_jets=)` (internal remap helper; "child" =
the replaced level relative to its parent).
## Responsibilities
- Construction: `_check_levels` (inputs), `_kernel_configs` (per-level kernel configs),
`Hierarchy.edge_frame` (geometry frames), `_build_level` / `_build_pair` / `_upload_boundary`
(device). Observation reads share `_observed_boundary` + `_summarize`.
- Blocked kernels: `operators/block_collision.cuh block_collide` is the one collision extension
point of `BlockStep` and `BlockChannelStep`. MRT shares one path; SRT keeps its two explicit
arithmetic variants (level: unclamped, no-force, shifted direct-momentum feq; channel: uniform
OneStep arithmetic with omega clamp), selected by `BLOCK_ROOT_BC`.
- `explosion="projected_linear"` stays an optional capability; `"constant"` stays default.
## Not done / open
- No single-`Simulation` facade over FixedPatchStepper: the two own different state models
(uniform node frame, body sync, forcing, checkpoint vs. immutable blocked hierarchy, edge
frame, remap-as-new-instance). A facade would only forward calls; no concrete need found.
- Tests, examples and validation callers are migrated separately (pinned runners stay frozen).
+1 -1
View File
@@ -61,7 +61,7 @@ Stationary whole-circle L2 bodies and Jets (`level_objects[2]`, `level_jets[2]`)
coordinates and a single L2 boundary object. Body IDs are physical object-list coordinates and a single L2 boundary object. Body IDs are physical object-list
positions, not region IDs or caller `obj_id`. Complete boundary stencils must positions, not region IDs or caller `obj_id`. Complete boundary stencils must
avoid the transfer footprint. Change a child Jet with avoid the transfer footprint. Change a child Jet with
`set_jet_amplitude(index, value, level=2)`; L1 bodies/Jets remain unsupported. `set_jet_amplitude(index, value, level=2)` (`level` is required); L1 bodies/Jets remain unsupported.
## Verification ## Verification
+2 -2
View File
@@ -6,8 +6,8 @@ it is; it adds no framework. Line references are to candidate 0f8fdbc3.
| Class | Contents | How it changes | Owner rule | | Class | Contents | How it changes | Owner rule |
|---|---|---|---| |---|---|---|---|
| 1. Module and its globals | Compile-time (visible to the compiler, fixed per module): `config.h` macros incl. `NX`/`NY`, `COLLISION_MODEL`, `BLOCK_ROOT_BC`, `WALL_*`. Runtime per-module globals (not visible to the compiler; written by host code): `__constant__` `d_params`, `d_w`, `d_cx`, `d_cy`, `d_block_channel`, `d_block_wall`; plus bound kernel handles | Globals are written once by the owner right after `load_module`. Uniform grid: `Simulation._recompile` (construction, body-topology sync) loads a new module and re-uploads; `update_runtime_params` rewrites `d_params` (omega, force, n_objects) between runs. Steppers: never rewritten after construction, so they are immutable for the instance's lifetime although they are runtime values. The blocked slot/reservoir root reads its slot from the `d_block_wall` constant (the former dense route used `WALL_*` macros; both were bit-identical in the parity gates). | Every owner loads its own module (`compiler.load_module` = `module_from_file`, no cache). A loaded module is never shared between live owners: globals are per module, so sharing would share mutable `d_params`. | | 1. Module and its globals | Compile-time (visible to the compiler, fixed per module): `config.h` macros incl. `NX`/`NY`, `COLLISION_MODEL`, `BLOCK_ROOT_BC`, `WALL_*`. Runtime per-module globals (not visible to the compiler; written by host code): `__constant__` `d_params`, `d_w`, `d_cx`, `d_cy`, `d_block_channel`, `d_block_wall`; plus bound kernel handles | Globals are written once by the owner right after `load_module`. Uniform grid: `Simulation._recompile` (construction, body-topology sync) loads a new module and re-uploads; `update_runtime_params` rewrites `d_params` (omega, force, n_objects) between runs. Steppers: never rewritten after construction, so they are immutable for the instance's lifetime although they are runtime values. The blocked slot/reservoir root reads its slot from the `d_block_wall` constant (the former dense route used `WALL_*` macros; both were bit-identical in the parity gates). | Every owner loads its own module (`compiler.load_module` = `module_from_file`, no cache). A loaded module is never shared between live owners: globals are per module, so sharing would share mutable `d_params`. |
| 2. Geometry, layout, tables | flags, curved-link tables, block layouts and halo pairs, transfer (explode/coalesce) tables, block kinds | Immutable for the life of an instance. Uniform grid: rebuilt by body sync. `FixedPatchStepper`: a geometry change is a new instance (`remap_child`), committed only after a successful transfer. | No in-place table edits. | | 2. Geometry, layout, tables | flags, curved-link tables, block layouts and halo pairs, transfer (explode/coalesce) tables, block kinds | Immutable for the life of an instance. Uniform grid: rebuilt by body sync. `FixedPatchStepper`: a geometry change is a new instance (`remap_finest`), committed only after a successful transfer. | No in-place table edits. |
| 3. Controls and dependent boundary columns | Jet amplitude -> per-link `surface_velocity`; body actions (uniform grid) | Between steps, without a table or module rebuild (`set_jet_amplitude`). | Input (index, level, amplitude) is validated before any device write; invalid input changes nothing and the instance stays usable. The host copy commits only after the upload returned. Any error after validation leaves the device column undefined and fails the instance: `FixedPatchStepper` status `failed`, `Simulation._failed`; every later call except `close()` raises, `close()` still releases resources and stays retryable. | | 3. Controls and dependent boundary columns | Jet amplitude -> per-link `surface_velocity`; body actions (uniform grid) | Between steps, without a table or module rebuild (`set_jet_amplitude`). | Input (index, level, amplitude) is validated before any device write; invalid input changes nothing and the instance stays usable. The host copy commits only after the upload returned. Any error after validation leaves the device column undefined and fails the instance: `FixedPatchStepper` status `failed`, `Simulation._failed`; every later device-touching public method (`Simulation`: `stream`, `set_jet_amplitude`, `initialize`, `sync_bodies`, `run`/`step`, `read_surface_exchange`, `read_body(s)`, `get_macroscopic`, `get_ddf`, `update_runtime_params`, `snapshot`, `restore`, save/load checkpoint) raises; cached host readers and host-side staging (e.g. `get_flags`, `add_body`/`remove_body` before a sync) do not touch the device and stay callable. `close()` still releases resources and stays retryable. |
| 4. Per-step state | population double buffers and their parity, raw-exchange ledgers, observations | Device-only during `step(n)` / `run(n)`; buffer swaps are host-handle swaps. | Valid on the host only after the single wait at the end of the call. | | 4. Per-step state | population double buffers and their parity, raw-exchange ledgers, observations | Device-only during `step(n)` / `run(n)`; buffer swaps are host-handle swaps. | Valid on the host only after the single wait at the end of the call. |
Storage routes are separate concerns: Storage routes are separate concerns:
+2 -2
View File
@@ -35,7 +35,7 @@ Tests: `tests/unit/test_refinement_mrt.py`, `tests/validation/multigrid/test_ref
| `slot_reservoir` | on-site velocity bottom slot, reservoir top, velocity 0, 2 <= slot width < nx, patch rows 4..ny-4 | `lbm/refinement_stepper.py` `FixedPatchStepper.__init__`; `lbm/refinement.py` `prepare_fine_boundary` | | `slot_reservoir` | on-site velocity bottom slot, reservoir top, velocity 0, 2 <= slot width < nx, patch rows 4..ny-4 | `lbm/refinement_stepper.py` `FixedPatchStepper.__init__`; `lbm/refinement.py` `prepare_fine_boundary` |
| `channel` | all levels blocked; Zou-He inlet, pressure-NEQ outlet, stationary bounce-back walls; transfers away from the closure | `lbm/channel.py` `validate_channel_config`, `validate_channel_transfers` | | `channel` | all levels blocked; Zou-He inlet, pressure-NEQ outlet, stationary bounce-back walls; transfers away from the closure | `lbm/channel.py` `validate_channel_config`, `validate_channel_transfers` |
| Explosion | `constant` (default); `projected_linear` only periodic, all-blocked SRT, no bodies/Jets, no LES (experimental) | `lbm/refinement_stepper.py` `FixedPatchStepper.__init__` | | Explosion | `constant` (default); `projected_linear` only periodic, all-blocked SRT, no bodies/Jets, no LES (experimental) | `lbm/refinement_stepper.py` `FixedPatchStepper.__init__` |
| L2 runtime remap (`remap_child`) | periodic three-level L2 sets only, no L0/L1 bodies, identical pinned boundary | `lbm/refinement_stepper.py` `FixedPatchStepper.remap_child` | | L2 runtime remap (`remap_finest`) | periodic three-level L2 sets only, no L0/L1 bodies, identical pinned boundary | `lbm/refinement_stepper.py` `FixedPatchStepper.remap_finest` |
Tests: `tests/unit/test_refinement_three_level_multipatch.py`, `tests/validation/multigrid/test_refinement_multipatch.py`, Tests: `tests/unit/test_refinement_three_level_multipatch.py`, `tests/validation/multigrid/test_refinement_multipatch.py`,
`test_refinement_three_level.py`, `test_refinement_three_level_multipatch_gpu.py`, `test_refinement_blocked_base.py`, `test_refinement_three_level.py`, `test_refinement_three_level_multipatch_gpu.py`, `test_refinement_blocked_base.py`,
@@ -84,6 +84,6 @@ Tests: `tests/unit/test_jets.py`, `test_jets_review.py`, `test_shifted_jets.py`,
## Not supported in this release ## Not supported in this release
3D, dynamic/automatic AMR (tagging, time-varying hierarchies beyond `remap_child`), 3D, dynamic/automatic AMR (tagging, time-varying hierarchies beyond `remap_finest`),
moving bodies on refined levels, pressure Jets, Jet checkpoints, LES or TRT on refined moving bodies on refined levels, pressure Jets, Jet checkpoints, LES or TRT on refined
levels, partial (non-dividing) tiles. levels, partial (non-dividing) tiles.
+5 -5
View File
@@ -33,7 +33,7 @@ origin is retained. Physical lower edges are `origin + spacing*(x,y)` on L0,
`origin + spacing*(parent.fine_origin + (x,y)/4)` on L2. `origin + spacing*(parent.fine_origin + (x,y)/4)` on L2.
Factory options `coarse_boundary` and `objects` must declare `"periodic"` and Factory options `coarse_boundary` and `objects` must declare `"periodic"` and
an empty collection; L2 bodies use `child_objects` as described below. an empty collection; L2 bodies use `objects` as described below.
LBMConfig's channel settings do not select this operator's boundary mode. LBMConfig's channel settings do not select this operator's boundary mode.
Geometry metadata is factory-created; Geometry metadata is factory-created;
manual mutation/forging of metadata is unsupported. Caller certifies the manual mutation/forging of metadata is unsupported. Caller certifies the
@@ -77,7 +77,7 @@ PYTHONPATH=src python -m unittest discover -s tests/unit -p 'test_synchronized_r
## Explicit runtime continuation ## Explicit runtime continuation
`replacement = current.remap_child(rectangles)` returns a new ready `replacement = current.remap_finest(rectangles)` returns a new ready
`FixedPatchStepper` at the **same absolute step_count**. It supports only `FixedPatchStepper` at the **same absolute step_count**. It supports only
periodic three-level runs with no L0/L1 bodies and an existing L2 periodic three-level runs with no L0/L1 bodies and an existing L2
`FixedPatchSet2D`; stationary L2 bodies/Jets follow the support contract below. `FixedPatchSet2D`; stationary L2 bodies/Jets follow the support contract below.
@@ -99,7 +99,7 @@ resources; an actual old-device synchronization/read failure marks the old
stepper failed, consistent with its existing failure semantics. stepper failed, consistent with its existing failure semantics.
```python ```python
replacement = current.remap_child([(20, 16, 36, 40)]) replacement = current.remap_finest([(20, 16, 36, 40)])
# Explicitly choose the returned continuation only after success. # Explicitly choose the returned continuation only after success.
current.close() current.close()
current = replacement current = replacement
@@ -119,7 +119,7 @@ Per-population conservation is not an advancement target under collision.
### Bounded pinned-body extension ### Bounded pinned-body extension
`replacement = current.remap_child(rectangles, child_objects=original_objects)` `replacement = current.remap_finest(rectangles, objects=original_objects)`
also accepts stationary Circles inside the L2 union. Each body's **physical also accepts stationary Circles inside the L2 union. Each body's **physical
boundary, complete BFL support and current per-link velocity must remain boundary, complete BFL support and current per-link velocity must remain
identical**. Rectangles are not pinned: they may be split, merged or resized, identical**. Rectangles are not pinned: they may be split, merged or resized,
@@ -130,7 +130,7 @@ coverage and two-hop BFL transfer-footprint checks still apply. This is manual
remapping, not adaptive body tracking or runtime body editing. remapping, not adaptive body tracking or runtime body editing.
Supply all original root edge-frame definitions in their construction order. Supply all original root edge-frame definitions in their construction order.
Do not pass `current.child_boundary.objects`: those are projected fine-local copies. Do not pass `current.boundaries[2].objects`: those are projected fine-local copies.
There is no retained mutable original-object registry. Before reading device There is no retained mutable original-object registry. Before reading device
state, remapping reprojects the supplied definition onto the old topology and state, remapping reprojects the supplied definition onto the old topology and
compares exact ordered circle centers/radii, flags, physical cut links, body compares exact ordered circle centers/radii, flags, physical cut links, body
+3 -3
View File
@@ -40,13 +40,13 @@ def config(mode):
def leaf_masks(run, mode): def leaf_masks(run, mode):
"""Known fluid leaves: uncovered coarse cells (minus wall rows in slot mode) """Known fluid leaves: uncovered coarse cells (minus wall rows in slot mode)
and fine ACTIVE cells outside the body.""" and fine ACTIVE cells outside the body."""
p = run.patch p = run.hierarchy.pairs[0]
x0, y0, x1, y1 = p.bounds x0, y0, x1, y1 = p.bounds
coarse = np.ones(p.coarse_shape, bool) coarse = np.ones(p.coarse_shape, bool)
coarse[y0:y1, x0:x1] = False coarse[y0:y1, x0:x1] = False
if mode == "slot_reservoir": if mode == "slot_reservoir":
coarse[[0, -1]] = False coarse[[0, -1]] = False
solid = (run.fine_boundary.flags.reshape(p.fine_shape) & SOLID) != 0 solid = (run.boundaries[1].flags.reshape(p.fine_shape) & SOLID) != 0
return coarse, (p.fine_kind == ACTIVE) & ~solid return coarse, (p.fine_kind == ACTIVE) & ~solid
@@ -94,7 +94,7 @@ def main():
cfg, context=context, coarse_boundary=args.mode, refinements=[BOUNDS], cfg, context=context, coarse_boundary=args.mode, refinements=[BOUNDS],
level_blocks=[tuple(args.base_block), tuple(args.fine_block)], level_blocks=[tuple(args.base_block), tuple(args.fine_block)],
level_objects=[(), (body,)]) as run: level_objects=[(), (body,)]) as run:
p = run.patch p = run.hierarchy.pairs[0]
# Rest equilibrium is already postcollision G, so no initial collision is # Rest equilibrium is already postcollision G, so no initial collision is
# needed. Non-leaf entries (covered, ghost, solid) are discarded by initialize. # needed. Non-leaf entries (covered, ghost, solid) are discarded by initialize.
rest = np.asarray(compute_feq_d2q9(1.0, 0.0, 0.0), np.float32)[:, None, None] rest = np.asarray(compute_feq_d2q9(1.0, 0.0, 0.0), np.float32)[:, None, None]
+1 -1
View File
@@ -153,7 +153,7 @@ def main():
with CudaContext(0) as context, FixedPatchStepper(cfg, context=context, with CudaContext(0) as context, FixedPatchStepper(cfg, context=context,
coarse_boundary='slot_reservoir', refinements=[BOUNDS], coarse_boundary='slot_reservoir', refinements=[BOUNDS],
level_blocks=LEVEL_BLOCKS) as run: level_blocks=LEVEL_BLOCKS) as run:
p = run.patch p = run.hierarchy.pairs[0]
masks = [np.ones(p.coarse_shape, bool), p.fine_kind == ACTIVE] masks = [np.ones(p.coarse_shape, bool), p.fine_kind == ACTIVE]
masks[0][[0, -1]] = False masks[0][[0, -1]] = False
masks[0][4:28, 16:32] = False masks[0][4:28, 16:32] = False
+3 -3
View File
@@ -16,7 +16,7 @@ from CelerisLab.config import load_lbm_config
from CelerisLab.cuda.context import CudaContext from CelerisLab.cuda.context import CudaContext
from CelerisLab.lbm.equilibrium import compute_feq_d2q9, compute_macro_from_ddf from CelerisLab.lbm.equilibrium import compute_feq_d2q9, compute_macro_from_ddf
from CelerisLab.lbm.refinement import ACTIVE from CelerisLab.lbm.refinement import ACTIVE
from CelerisLab.lbm.refinement_stepper import PeriodicPatchStepper from CelerisLab.lbm.refinement_stepper import FixedPatchStepper
def shear_equilibrium(shape, origin, dx, period): def shear_equilibrium(shape, origin, dx, period):
@@ -48,9 +48,9 @@ def main():
nx=16, ny=16, viscosity=.1, velocity=1e-3, rho=1.0, nx=16, ny=16, viscosity=.1, velocity=1e-3, rho=1.0,
lattice_model="D2Q9", collision="SRT", streaming="double_buffer", lattice_model="D2Q9", collision="SRT", streaming="double_buffer",
store_precision="FP32", ddf_shifting=False, les_enabled=False) store_precision="FP32", ddf_shifting=False, les_enabled=False)
with CudaContext(0) as context, PeriodicPatchStepper(cfg, context=context, refinements=[(4, 4, 12, 12)], with CudaContext(0) as context, FixedPatchStepper(cfg, context=context, refinements=[(4, 4, 12, 12)],
level_blocks=[(16, 16), (24, 24)]) as run: level_blocks=[(16, 16), (24, 24)]) as run:
p = run.patch p = run.hierarchy.pairs[0]
coarse = shear_equilibrium(p.coarse_shape, (0, 0), 1.0, cfg.ny) coarse = shear_equilibrium(p.coarse_shape, (0, 0), 1.0, cfg.ny)
fine = shear_equilibrium(p.fine_shape, p.fine_origin, .5, cfg.ny) fine = shear_equilibrium(p.fine_shape, p.fine_origin, .5, cfg.ny)
coarse_leaf = np.ones(p.coarse_shape, dtype=bool) coarse_leaf = np.ones(p.coarse_shape, dtype=bool)
+25 -3
View File
@@ -21,7 +21,17 @@ from .refinement import (ACTIVE, OUTER_GHOST, FixedPatch2D, FixedPatchSet2D,
@dataclass(frozen=True, eq=False) @dataclass(frozen=True, eq=False)
class Hierarchy: class Hierarchy:
"""configs[l]: validated LBMConfig of level l; pairs[l]: tables of (l, l+1); """configs[l]: validated LBMConfig of level l; pairs[l]: tables of (l, l+1);
kinds[l]: composed per-cell kind (child coverage OUTER on every non-finest level).""" kinds[l]: composed per-cell kind (child coverage OUTER on every non-finest level).
Index conventions (the single place they are defined):
* level l = 0 (root) .. depth-1; pair index p = l-1 builds level l from level p;
* array shapes are (ny, nx); rectangles, origins and centers are (x, y);
* the root EDGE frame puts root cell i at [i, i+1), centers at i+.5 (bodies and
Jets of FixedPatchStepper); uniform Simulation objects use the integer-NODE frame;
* cell (i, j) of level l's allocated array (ghost band and holes included) covers
[o + h*i, o + h*(i+1)) with (o, h) = edge_frame(l); level 1 is built from the root
config directly, deeper levels from their parent frame (see edge_frame).
"""
configs: tuple configs: tuple
pairs: tuple pairs: tuple
kinds: tuple kinds: tuple
@@ -30,6 +40,18 @@ class Hierarchy:
def depth(self): def depth(self):
return len(self.kinds) return len(self.kinds)
def edge_frame(self, level):
"""(origin (x, y), spacing) of level ``level``'s array in the root edge frame.
Each pair's fine_origin is in parent cells: zero for a set (level-wide grid),
the shifted window origin for a scalar FixedPatch2D.
"""
origin, spacing = (0.0, 0.0), 1.0
for pair in self.pairs[:level]:
origin = tuple(o + spacing * f for o, f in zip(origin, pair.fine_origin))
spacing /= 2
return origin, spacing
def is_rectangle(value): def is_rectangle(value):
value = tuple(value) value = tuple(value)
@@ -86,14 +108,14 @@ def _check_chain(coarse, fine, pair):
def leaf_partition(h): def leaf_partition(h):
"""Base-cell fraction published as ACTIVE over all levels (exact: powers of two). """Base-cell fraction published as ACTIVE over all levels (exact: powers of two).
Window (scalar) levels are offset by their fine_level_origin first. Window (scalar) levels are offset by their fine_level_offset first.
""" """
base = np.zeros(h.kinds[0].shape) base = np.zeros(h.kinds[0].shape)
oy = ox = 0 oy = ox = 0
for level, kind in enumerate(h.kinds): for level, kind in enumerate(h.kinds):
if level: if level:
pair = h.pairs[level - 1] pair = h.pairs[level - 1]
py, px = (pair.fine_level_origin if isinstance(pair, FixedPatch2D) else (0, 0)) py, px = (pair.fine_level_offset if isinstance(pair, FixedPatch2D) else (0, 0))
oy, ox = 2 * oy + py, 2 * ox + px oy, ox = 2 * oy + py, 2 * ox + px
y, x = np.nonzero(kind == ACTIVE) y, x = np.nonzero(kind == ACTIVE)
np.add.at(base, ((y + oy) >> level, (x + ox) >> level), 0.25 ** level) np.add.at(base, ((y + oy) >> level, (x + ox) >> level), 0.25 ** level)
-2
View File
@@ -34,9 +34,7 @@
#if BLOCK_ROOT_BC == 1 #if BLOCK_ROOT_BC == 1
#include "boundary/block_wall.cuh" #include "boundary/block_wall.cuh"
#endif #endif
#if BLOCK_ROOT_BC != 2
#include "operators/block_collision.cuh" #include "operators/block_collision.cuh"
#endif
// Kernels, in definition order: transfers, boundary staging, then the one full // Kernels, in definition order: transfers, boundary staging, then the one full
// step selected for this module (BlockStep, or BlockChannelStep for a channel root). // step selected for this module (BlockStep, or BlockChannelStep for a channel root).
@@ -44,16 +44,4 @@ struct RigidBodyControl2D {
int control_mode; // 0=force, 1=velocity, 2=displacement int control_mode; // 0=force, 1=velocity, 2=displacement
}; };
// ============================================================================
// Halo plan (multi-GPU / AMR 预留)
// ============================================================================
struct HaloPlan {
int face_id;
int peer_gpu;
int level;
size_t count;
int* send_idx;
int* recv_idx;
};
#endif // CELERIS_CORE_PARAMS_CUH #endif // CELERIS_CORE_PARAMS_CUH
@@ -1,4 +1,10 @@
// Collision used by periodic/fine and slot/reservoir BlockStep. // The one collision extension point of the blocked steps (BlockStep, BlockChannelStep).
// Compile-time selection: MRT shares one path; SRT keeps two explicit arithmetic
// variants, unchanged from their original kernels:
// level (BLOCK_ROOT_BC 0/1): unclamped omega, no forcing term; shifted storage
// builds feq from the retained drho and direct momentum (numerical variant);
// channel (BLOCK_ROOT_BC 2): matches uniform OneStep SRT arithmetic (omega clamp,
// zero Guo forcing term, compute_feq).
#ifndef CELERIS_OPERATORS_BLOCK_COLLISION_CUH #ifndef CELERIS_OPERATORS_BLOCK_COLLISION_CUH
#define CELERIS_OPERATORS_BLOCK_COLLISION_CUH #define CELERIS_OPERATORS_BLOCK_COLLISION_CUH
@@ -11,7 +17,13 @@
__device__ __forceinline__ void block_collide(float* f) { __device__ __forceinline__ void block_collide(float* f) {
#if COLLISION_MODEL == 2 #if COLLISION_MODEL == 2
collide_mrt_unforced(f); // SRT path below unchanged. collide_mrt_unforced(f);
#elif BLOCK_ROOT_BC == 2
float rho, ux, uy, feq[NQ], forcing[NQ] = {};
compute_rho_u(f, rho, ux, uy);
compute_feq(rho, ux, uy, feq);
const float omega = fminf(OMEGA_COLLISION_MAX, fmaxf(OMEGA_COLLISION_MIN, d_params.omega));
collide_srt(f, feq, forcing, omega);
#else #else
float rho, ux, uy, feq[NQ]; float rho, ux, uy, feq[NQ];
compute_rho_u(f, rho, ux, uy); compute_rho_u(f, rho, ux, uy);
@@ -57,25 +57,4 @@ __device__ __forceinline__ void collide_trt(float* __restrict__ f,
} }
} }
__device__ __forceinline__ void collide_trt_no_force(float* __restrict__ f,
const float* __restrict__ feq,
float omega)
{
const float wp = omega;
const float wm = compute_omega_minus(wp);
f[0] = f[0] - wp * (f[0] - feq[0]);
#pragma unroll
for (int i = 1; i < NQ; i += 2) {
const int ib = i + 1;
const float fi = f[i], fb = f[ib];
const float feqi = feq[i], feqb = feq[ib];
const float delta_s = 0.5f * ((fi - feqi) + (fb - feqb));
const float delta_a = 0.5f * ((fi - feqi) - (fb - feqb));
f[i] = fi - wp * delta_s - wm * delta_a;
f[ib] = fb - wp * delta_s + wm * delta_a;
}
}
#endif // CELERIS_OPERATORS_COLLISION_TRT_CUH #endif // CELERIS_OPERATORS_COLLISION_TRT_CUH
@@ -48,18 +48,7 @@ extern "C" __global__ void BlockChannelStep(
} else { } else {
apply_wall_bb_d2q9_rows(gy, d_block_channel.ny - 2, f, in + base, k); apply_wall_bb_d2q9_rows(gy, d_block_channel.ny - 2, f, in + base, k);
} }
if (!wall && !outlet) { if (!wall && !outlet) block_collide(f);
#if COLLISION_MODEL == 2
collide_mrt_unforced(f);
#else
// Match uniform SRT arithmetic, including its shifted equilibrium path.
float rho, ux, uy, feq[NQ], forcing[NQ] = {};
compute_rho_u(f, rho, ux, uy);
compute_feq(rho, ux, uy, feq);
const float omega = fminf(OMEGA_COLLISION_MAX, fmaxf(OMEGA_COLLISION_MIN, d_params.omega));
collide_srt(f, feq, forcing, omega);
#endif
}
} }
stream_pull_store(k, f, out + base); stream_pull_store(k, f, out + base);
} }
+3 -3
View File
@@ -115,8 +115,9 @@ class FixedPatch2D:
return self.bounds[0] - HALO / 2, self.bounds[1] - HALO / 2 return self.bounds[0] - HALO / 2, self.bounds[1] - HALO / 2
@property @property
def fine_level_origin(self) -> tuple[int, int]: def fine_level_offset(self) -> tuple[int, int]:
"""(y, x) of this window's first cell on the level-wide 2:1 fine grid.""" """Array index (y, x) of this window's first cell on the level-wide 2:1 fine grid
(an index offset, not an (x, y) origin)."""
return 2 * self.bounds[1] - HALO, 2 * self.bounds[0] - HALO return 2 * self.bounds[1] - HALO, 2 * self.bounds[0] - HALO
@property @property
@@ -151,7 +152,6 @@ class FixedPatch2D:
raise ValueError("levels are not synchronized: fine_steps must equal 2*coarse_steps") raise ValueError("levels are not synchronized: fine_steps must equal 2*coarse_steps")
MIN_PATCH_GAP = 4 # Legacy exported constant; union transfers impose no gap policy.
@dataclass(frozen=True, eq=False) @dataclass(frozen=True, eq=False)
+113 -125
View File
@@ -1,10 +1,11 @@
# CelerisLab/lbm/refinement_stepper.py # CelerisLab/lbm/refinement_stepper.py
"""Opt-in fixed patch backend with optional stationary fine-interior circles. """Blocked fixed 2:1 refinement hierarchy (any depth) with optional stationary circles/Jets.
The caller owns the CudaContext and must close this backend before that context. The caller owns the CudaContext and must close this backend before that context.
Only synchronized postcollision_g leaf populations are public. No checkpoint, Only synchronized postcollision_g leaf populations are public. No checkpoint,
forcing or runtime parameter mutation is offered. Bounded manual periodic forcing or runtime parameter mutation besides Jet amplitudes is offered. A bounded
L2 set replacement is available at completed base intervals via remap_child. manual periodic level-2 set replacement is available at completed base intervals
via remap_finest. Level/pair/shape/frame conventions: lbm.hierarchy.Hierarchy.
""" """
from dataclasses import dataclass, replace from dataclasses import dataclass, replace
from numbers import Integral from numbers import Integral
@@ -62,8 +63,22 @@ class _Level:
boundary: object = None boundary: object = None
def _level_name(level): def _kernel_configs(cfg, pairs, wall_aware):
return ("base", "fine", "child")[level] if level < 3 else f"level {level}" """Per-level kernel LBMConfig (window shape, viscosity doubling per level).
Geometry configs live in Hierarchy.configs; level 1 of a slot_reservoir root
replaces the root wall fields by the internal bounce-back/velocity ghost setup.
"""
configs = [replace(cfg)]
for index, pair in enumerate(pairs):
parent = configs[-1]
level = replace(parent, nx=pair.fine_shape[1], ny=pair.fine_shape[0],
viscosity=2 * parent.viscosity)
if index == 0 and wall_aware:
level = replace(level, y_wall_bc="bounce_back", top_wall_bc="velocity",
wall_slot_width=0, wall_normal_velocity=0.0, velocity=0.0)
configs.append(level)
return configs
MAX_GRID_DIM = 65535 # CUDA gridDim.y/z limit (rows, blocks) MAX_GRID_DIM = 65535 # CUDA gridDim.y/z limit (rows, blocks)
@@ -163,11 +178,18 @@ class FixedPatchStepper:
lattice velocity); every whole circle and its BFL stencil must fit strictly lattice velocity); every whole circle and its BFL stencil must fit strictly
away from that level's transfer footprint. away from that level's transfer footprint.
``explosion="projected_linear"`` selects experimental transverse linear ``explosion`` selects the parent->child reconstruction of every pair:
reconstruction on each adjacent pair. Default "constant" is unchanged. ``"constant"`` (default) or ``"projected_linear"``, an optional transverse
Projected mode requires periodic root, all-blocked SRT levels without linear reconstruction (Schornbaum-Ruede/Chen style) that is GPU-confirmed
bodies/Jets or LES; this backend supplies no forcing. Geometry is generic, more accurate than constant only on periodic two-level single-rectangle
but union/multilevel GPU and matched-physics qualification remain pending. Taylor-Green cases; union/multilevel GPU and matched-physics qualification
remain pending. It requires a periodic root and SRT levels without
bodies/Jets or LES; this backend supplies no forcing.
Per-level read-only state: ``hierarchy.pairs`` (pair l-1 builds level l),
``layouts``, ``boundaries`` (host geometry), ``jets`` and ``leaf`` masks,
each indexed by level. ``set_jet_amplitude``, ``read_raw_exchange`` and
``read_exchange_window`` take an explicit ``level=``.
One recursive hop advances every level: a child explodes from its parent's One recursive hop advances every level: a child explodes from its parent's
current hop source, takes two hops, then the parent gathers halos, stages current hop source, takes two hops, then the parent gathers halos, stages
@@ -233,42 +255,26 @@ class FixedPatchStepper:
layout = BlockLayout2D.build(*pairs[level - 1].fine_shape, *block, layout = BlockLayout2D.build(*pairs[level - 1].fine_shape, *block,
periodic=(False, False)) periodic=(False, False))
if layout.nblocks > MAX_GRID_DIM or layout.local_shape[0] > MAX_GRID_DIM: if layout.nblocks > MAX_GRID_DIM or layout.local_shape[0] > MAX_GRID_DIM:
raise ValueError(f"{_level_name(level)} blocks exceed launch limits") raise ValueError(f"level {level} blocks exceed launch limits")
self._layouts.append(layout) self._layouts.append(layout)
if self._channel: if self._channel:
base = blocked_base_tables(pairs[0], self._layouts[0], np.argsort(pairs[0].coalesce_coarse)) base = blocked_base_tables(pairs[0], self._layouts[0], np.argsort(pairs[0].coalesce_coarse))
validate_channel_transfers(pairs[0], self._layouts[0], base[1]) validate_channel_transfers(pairs[0], self._layouts[0], base[1])
# Kernel configurations per level (geometry configs live in self.hierarchy). self._level_cfgs = _kernel_configs(cfg, pairs, self._wall_aware)
self._level_cfgs = [replace(cfg)] # Host boundary geometry per level in the root edge frame (Hierarchy.edge_frame).
if pairs: self._boundaries = [None]
fine_cfg = replace(cfg, nx=pairs[0].fine_shape[1], ny=pairs[0].fine_shape[0],
viscosity=2 * cfg.viscosity)
if self._wall_aware:
fine_cfg = replace(fine_cfg, y_wall_bc="bounce_back", top_wall_bc="velocity",
wall_slot_width=0, wall_normal_velocity=0.0, velocity=0.0)
self._level_cfgs.append(fine_cfg)
for pair in pairs[1:]:
parent = self._level_cfgs[-1]
self._level_cfgs.append(replace(parent, nx=pair.fine_shape[1], ny=pair.fine_shape[0],
viscosity=2 * parent.viscosity))
# Per-level boundary geometry in the root edge frame (level l spacing 2**-l).
self._boundary_hosts = [None]
origin, spacing = (0.0, 0.0), 1.0
spacings = [spacing]
for level in range(1, depth): for level in range(1, depth):
pair = pairs[level - 1] pair = pairs[level - 1]
if level == 1: if level == 1: # the root config carries the wall-aware boundary fields
host = prepare_fine_boundary(cfg, pair, level_objects[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) _shifted_blocked=self._ddf_shifting)
else: else:
origin, spacing = self.hierarchy.edge_frame(level - 1)
host = prepare_fine_boundary(self._level_cfgs[level - 1], pair, level_objects[level], host = prepare_fine_boundary(self._level_cfgs[level - 1], pair, level_objects[level],
jets=level_jets[level], parent_origin=origin, jets=level_jets[level], parent_origin=origin,
parent_spacing=spacing, _shifted_blocked=self._ddf_shifting) parent_spacing=spacing, _shifted_blocked=self._ddf_shifting)
self._boundary_hosts.append(host) self._boundaries.append(host)
origin = tuple(o + spacing * f for o, f in zip(origin, pair.fine_origin))
spacing /= 2
spacings.append(spacing)
if not isinstance(observe_raw_exchange, bool): if not isinstance(observe_raw_exchange, bool):
raise ValueError("observe_raw_exchange must be bool") raise ValueError("observe_raw_exchange must be bool")
self._observe_raw_exchange = observe_raw_exchange self._observe_raw_exchange = observe_raw_exchange
@@ -280,7 +286,7 @@ class FixedPatchStepper:
self._window_valid = False self._window_valid = False
self._window_start = None self._window_start = None
self._exchange_valid = False self._exchange_valid = False
self._exchange_spacing = tuple(spacings) self._exchange_spacing = tuple(self.hierarchy.edge_frame(l)[1] for l in range(depth))
self._exchange_centers = tuple(tuple(tuple(o.geometry.center) for o in objects) self._exchange_centers = tuple(tuple(tuple(o.geometry.center) for o in objects)
for objects in level_objects) for objects in level_objects)
self._level_jets = level_jets self._level_jets = level_jets
@@ -289,7 +295,7 @@ class FixedPatchStepper:
tables = static_circle_tables(self._layouts[0], level_objects[0]) tables = static_circle_tables(self._layouts[0], level_objects[0])
root = _BlockedBoundary(None, tables, (), False) root = _BlockedBoundary(None, tables, (), False)
self._level_boundaries = [root] + [ self._level_boundaries = [root] + [
_prepare_blocked_boundary(self._boundary_hosts[level], self._layouts[level], _prepare_blocked_boundary(self._boundaries[level], self._layouts[level],
self._level_cfgs[level], level_jets[level], self._level_cfgs[level], level_jets[level],
_shifted_blocked=self._ddf_shifting) _shifted_blocked=self._ddf_shifting)
for level in range(1, depth)] for level in range(1, depth)]
@@ -297,13 +303,13 @@ class FixedPatchStepper:
if root is not None: if root is not None:
self._leaf[0] &= ~root.tables["solid"] self._leaf[0] &= ~root.tables["solid"]
for level in range(1, depth): for level in range(1, depth):
host = self._boundary_hosts[level] host = self._boundaries[level]
if host is not None: if host is not None:
self._leaf[level] &= (host.flags.reshape(pairs[level - 1].fine_shape) & SOLID) == 0 self._leaf[level] &= (host.flags.reshape(pairs[level - 1].fine_shape) & SOLID) == 0
self._state = [mask.copy() for mask in self._leaf] self._state = [mask.copy() for mask in self._leaf]
if self._channel: if self._channel:
self._leaf[0] &= channel_masks((cfg.ny, cfg.nx))[0] self._leaf[0] &= channel_masks((cfg.ny, cfg.nx))[0]
# Private constructor inputs for the bounded manual L2 continuation. # Private constructor inputs for remap_finest.
self._remap_cfg = replace(cfg) self._remap_cfg = replace(cfg)
self._level_blocks = level_blocks self._level_blocks = level_blocks
self.ctx = context self.ctx = context
@@ -325,21 +331,8 @@ class FixedPatchStepper:
for index, pair in enumerate(pairs): for index, pair in enumerate(pairs):
self._pairs.append(self._build_pair(index, pair)) self._pairs.append(self._build_pair(index, pair))
for level in self._levels: for level in self._levels:
b = level.boundary if level.boundary is not None:
if b is not None: self._upload_boundary(level)
b.flags = self._upload(b.tables["flags"])
b.device = {k: self._upload(v) for k, v in b.tables.items()
if k in BLOCKED_BOUNDARY_DEVICE_COLUMNS}
if self._observe_raw_exchange:
b.device.update(rx=self._upload(b.tables["rx"]),
ry=self._upload(b.tables["ry"]))
count = len(b.tables["cell"])
b.last_raw = self._allocate(count * 4 * 4)
b.interval_integral = self._allocate(count * 4 * 8)
if self._observe_window:
b.window_integral = self._allocate(count * 4 * 8)
b.accumulate = level.module.get_function("BlockAccumulateExchange")
b.kernel = level.module.get_function("BlockCurvedBoundary")
except BaseException: except BaseException:
try: try:
self.close() self.close()
@@ -347,6 +340,22 @@ class FixedPatchStepper:
pass # Preserve the original construction failure. pass # Preserve the original construction failure.
raise raise
def _upload_boundary(self, level):
"""Device columns, BFL kernel and (when observed) the exchange ledgers of one level."""
b = level.boundary
b.flags = self._upload(b.tables["flags"])
b.device = {k: self._upload(v) for k, v in b.tables.items()
if k in BLOCKED_BOUNDARY_DEVICE_COLUMNS}
if self._observe_raw_exchange:
b.device.update(rx=self._upload(b.tables["rx"]), ry=self._upload(b.tables["ry"]))
count = len(b.tables["cell"])
b.last_raw = self._allocate(count * 4 * 4)
b.interval_integral = self._allocate(count * 4 * 8)
if self._observe_window:
b.window_integral = self._allocate(count * 4 * 8)
b.accumulate = level.module.get_function("BlockAccumulateExchange")
b.kernel = level.module.get_function("BlockCurvedBoundary")
def _planned_bytes(self): def _planned_bytes(self):
"""Both G buffers, halo pairs and boundary columns of every level (transfer tables excluded). """Both G buffers, halo pairs and boundary columns of every level (transfer tables excluded).
@@ -516,14 +525,15 @@ class FixedPatchStepper:
self._status = "failed" self._status = "failed"
raise raise
def remap_child(self, rectangles, *, child_objects=()): def remap_finest(self, rectangles, *, objects=()):
"""Return a NEW ready stepper with an L2 set replaced at this base time. """Return a NEW ready stepper with the level-2 set replaced at this base time.
Periodic three-level set children only; L0/L1 and block sizes stay fixed. Periodic three-level hierarchies whose level 2 is a set only; levels 0/1 and
A body run requires its original root-edge-frame ``child_objects``: block sizes stay fixed. A body run requires its original root-edge-frame
level-2 ``objects``:
stationary circles with identical physical boundary and BFL support. stationary circles with identical physical boundary and BFL support.
Other fluid coverage may change, even in the same connected union. Other fluid coverage may change, even in the same connected union.
Current L2 Jets (including updated amplitudes) Current level-2 Jets (including updated amplitudes)
are preserved. Local projected boundary objects are not valid inputs. are preserved. Local projected boundary objects are not valid inputs.
The lock and stream synchronization certify The lock and stream synchronization certify
current postcollision G and the completed 1:2:4 counts internally. current postcollision G and the completed 1:2:4 counts internally.
@@ -541,19 +551,20 @@ class FixedPatchStepper:
pairs = self.hierarchy.pairs pairs = self.hierarchy.pairs
if (self._wall_aware or self._channel or self.hierarchy.depth != 3 if (self._wall_aware or self._channel or self.hierarchy.depth != 3
or not isinstance(pairs[1], FixedPatchSet2D) or not isinstance(pairs[1], FixedPatchSet2D)
or self._boundary_hosts[1] is not None or self._boundaries[1] is not None
or any(b is not None for b in self._level_boundaries[:2])): or any(b is not None for b in self._level_boundaries[:2])):
raise ValueError("remap_child requires periodic three-level L2 sets with no L0/L1 bodies") raise ValueError("remap_finest requires a periodic three-level hierarchy with a "
"level-2 set and no level-0/1 bodies")
cfg = self._remap_cfg cfg = self._remap_cfg
child_objects = tuple(child_objects) objects = tuple(objects)
actual = self._level_boundaries[2] actual = self._level_boundaries[2]
jets = () if actual is None else actual.jets jets = () if actual is None else actual.jets
old = RemapTopology.from_config(cfg, pairs[0], pairs[1].bounds, old = RemapTopology.from_config(cfg, pairs[0], pairs[1].bounds,
child_objects=child_objects, child_jets=jets, child_objects=objects, child_jets=jets,
_shifted_blocked=self._ddf_shifting) _shifted_blocked=self._ddf_shifting)
if (not same_physical_boundary(old.boundary, self._boundary_hosts[2]) if (not same_physical_boundary(old.boundary, self._boundaries[2])
or (self._boundary_hosts[2] is not None and actual is None)): or (self._boundaries[2] is not None and actual is None)):
raise ValueError("child_objects must reproduce the original root-edge-frame boundary") raise ValueError("objects must reproduce the original root-edge-frame level-2 boundary")
if actual is not None: if actual is not None:
velocity = actual.tables.get("surface_velocity") velocity = actual.tables.get("surface_velocity")
if velocity is None: if velocity is None:
@@ -561,18 +572,17 @@ class FixedPatchStepper:
_shifted_blocked=self._ddf_shifting) _shifted_blocked=self._ddf_shifting)
if (not same_physical_boundary(old.boundary, actual.host) if (not same_physical_boundary(old.boundary, actual.host)
or not np.array_equal(old.surface_velocity, velocity)): or not np.array_equal(old.surface_velocity, velocity)):
raise ValueError("child_objects changed boundary ownership or surface velocity") raise ValueError("objects changed boundary ownership or surface velocity")
new = RemapTopology.from_config(cfg, pairs[0], rectangles, new = RemapTopology.from_config(cfg, pairs[0], rectangles,
child_objects=child_objects, child_jets=jets, child_objects=objects, child_jets=jets,
_shifted_blocked=self._ddf_shifting) _shifted_blocked=self._ddf_shifting)
if (not same_physical_boundary(old.boundary, new.boundary) if (not same_physical_boundary(old.boundary, new.boundary)
or not np.array_equal(old.surface_velocity, new.surface_velocity)): or not np.array_equal(old.surface_velocity, new.surface_velocity)):
raise ValueError("remap requires the identical pinned physical boundary and BFL support") raise ValueError("remap requires the identical pinned physical boundary and BFL support")
base_block, fine_block, child_block = self._level_blocks
# Check sparse block alignment/launch size before any device read/allocation. # Check sparse block alignment/launch size before any device read/allocation.
layout = new.child.fine_layout(child_block) layout = new.child.fine_layout(self._level_blocks[2])
if layout.nblocks > MAX_GRID_DIM or layout.local_shape[0] > MAX_GRID_DIM: if layout.nblocks > MAX_GRID_DIM or layout.local_shape[0] > MAX_GRID_DIM:
raise ValueError("child blocks exceed launch limits") raise ValueError("level 2 blocks exceed launch limits")
with self.ctx.activate(): with self.ctx.activate():
try: try:
self._stream.synchronize() self._stream.synchronize()
@@ -588,8 +598,8 @@ class FixedPatchStepper:
phase="postcollision_g", **options) phase="postcollision_g", **options)
candidate = type(self)( candidate = type(self)(
cfg, context=self.ctx, refinements=[pairs[0].bounds, new.child.bounds], cfg, context=self.ctx, refinements=[pairs[0].bounds, new.child.bounds],
level_blocks=[base_block, fine_block, child_block], level_blocks=list(self._level_blocks),
level_objects=[(), (), child_objects], level_jets=[(), (), jets], level_objects=[(), (), objects], level_jets=[(), (), jets],
observe_raw_exchange=self._observe_raw_exchange, observe_raw_exchange=self._observe_raw_exchange,
observe_window_exchange=self._observe_window, explosion=self._explosion) observe_window_exchange=self._observe_window, explosion=self._explosion)
try: try:
@@ -606,10 +616,9 @@ class FixedPatchStepper:
raise raise
return candidate return candidate
def set_jet_amplitude(self, index, amplitude, *, level=1): def set_jet_amplitude(self, index, amplitude, *, level):
"""Replace one jet's amplitude between steps; no geometry/table rebuild. """Replace Jet ``index`` of ``level_jets[level]`` between steps; no geometry/table rebuild.
``level=1`` preserves the L1 default; use ``level=2`` for child_jets.
Invalid input raises before any change. A device/sync error marks the Invalid input raises before any change. A device/sync error marks the
instance failed (as step/initialize do); host state commits only after instance failed (as step/initialize do); host state commits only after
a successful upload. a successful upload.
@@ -618,7 +627,7 @@ class FixedPatchStepper:
self._require("new", "ready") self._require("new", "ready")
if (isinstance(level, bool) or not isinstance(level, Integral) if (isinstance(level, bool) or not isinstance(level, Integral)
or not 1 <= level < self.hierarchy.depth): or not 1 <= level < self.hierarchy.depth):
raise ValueError("Jet level must name an existing refined level (1 = first)") raise ValueError("Jet level must name an existing refined level (1 .. depth-1)")
boundary = self._level_boundaries[level] boundary = self._level_boundaries[level]
jets = list(self._level_jets[level]) jets = list(self._level_jets[level])
jets[index] = jets[index].with_amplitude(amplitude) # validates first jets[index] = jets[index].with_amplitude(amplitude) # validates first
@@ -758,7 +767,23 @@ class FixedPatchStepper:
self._window_start = self._step_count self._window_start = self._step_count
self._window_valid = True self._window_valid = True
def read_exchange_window(self, body_id, *, level=1): def _observed_boundary(self, level, body_id):
if (isinstance(level, bool) or not isinstance(level, Integral)
or not 0 <= level < len(self._levels)):
raise ValueError("level must name an existing level")
b = self._level_boundaries[level]
if (b is None or isinstance(body_id, bool) or not isinstance(body_id, Integral)
or not 0 <= body_id < len(b.host.objects)):
raise ValueError("body_id must name a body on this level")
return b
def _summarize(self, b, rows, level, body_id, window, **kw):
return summarize_exchange(rows, b.tables["body_id"], b.tables["aperture"],
body_id=body_id, level=level, h=self._exchange_spacing[level],
time_end=self._step_count, window=window,
reference_center=self._exchange_centers[level][body_id], **kw)
def read_exchange_window(self, body_id, *, level):
"""Non-destructive read of [window start, current step count] in base-lattice units.""" """Non-destructive read of [window start, current step count] in base-lattice units."""
with self._lock, self.ctx.activate(): with self._lock, self.ctx.activate():
self._require("ready") self._require("ready")
@@ -769,26 +794,16 @@ class FixedPatchStepper:
duration = self._step_count - self._window_start duration = self._step_count - self._window_start
if duration < 1: if duration < 1:
raise RuntimeError("exchange window requires at least one completed base interval") raise RuntimeError("exchange window requires at least one completed base interval")
if (isinstance(level, bool) or not isinstance(level, Integral) b = self._observed_boundary(level, body_id)
or not 0 <= level < len(self._levels)):
raise ValueError("level must name an existing level")
b = self._level_boundaries[level]
if (b is None or isinstance(body_id, bool) or not isinstance(body_id, Integral)
or not 0 <= body_id < len(b.host.objects)):
raise ValueError("body_id must name a body on this level")
rows = np.empty((len(b.tables["cell"]), 4), np.float64) rows = np.empty((len(b.tables["cell"]), 4), np.float64)
try: try:
cuda.memcpy_dtoh(rows, b.window_integral) cuda.memcpy_dtoh(rows, b.window_integral)
except BaseException: except BaseException:
self._status = "failed" self._status = "failed"
raise raise
return summarize_exchange(rows, b.tables["body_id"], b.tables["aperture"], return self._summarize(b, rows, level, body_id, "window", duration=duration)
body_id=body_id, level=level, h=self._exchange_spacing[level],
time_end=self._step_count, window="window",
reference_center=self._exchange_centers[level][body_id],
duration=duration)
def read_raw_exchange(self, body_id, *, level=1, window="last_base_step"): def read_raw_exchange(self, body_id, *, level, window="last_base_step"):
"""Read raw boundary reaction about its geometric center in the root edge frame. """Read raw boundary reaction about its geometric center in the root edge frame.
Repeated reads and step(0) preserve the window. last_substep samples Repeated reads and step(0) preserve the window. last_substep samples
@@ -802,13 +817,7 @@ class FixedPatchStepper:
raise RuntimeError("raw exchange observation is disabled") raise RuntimeError("raw exchange observation is disabled")
if not self._exchange_valid: if not self._exchange_valid:
raise RuntimeError("raw exchange requires a completed base interval") raise RuntimeError("raw exchange requires a completed base interval")
if (isinstance(level, bool) or not isinstance(level, Integral) b = self._observed_boundary(level, body_id)
or not 0 <= level < len(self._levels)):
raise ValueError("level must name an existing level")
b = self._level_boundaries[level]
if (b is None or isinstance(body_id, bool) or not isinstance(body_id, Integral)
or not 0 <= body_id < len(b.host.objects)):
raise ValueError("body_id must name a body on this level")
if window not in ("last_base_step", "last_substep"): if window not in ("last_base_step", "last_substep"):
raise ValueError("window must be last_base_step or last_substep") raise ValueError("window must be last_base_step or last_substep")
count = len(b.tables["cell"]) count = len(b.tables["cell"])
@@ -818,10 +827,7 @@ class FixedPatchStepper:
except BaseException: except BaseException:
self._status = "failed" self._status = "failed"
raise raise
return summarize_exchange(rows, b.tables["body_id"], b.tables["aperture"], return self._summarize(b, rows, level, body_id, window)
body_id=body_id, level=level, h=self._exchange_spacing[level],
time_end=self._step_count, window=window,
reference_center=self._exchange_centers[level][body_id])
def download(self): def download(self):
"""Return physical G (FP64 for shifted storage), NaN on non-leaves.""" """Return physical G (FP64 for shifted storage), NaN on non-leaves."""
@@ -889,28 +895,10 @@ class FixedPatchStepper:
if exc[0] is None: if exc[0] is None:
raise raise
# Read-only compatibility views over the level lists, kept on purpose for the existing # Read-only per-level state, indexed by level (0 = root). Pairs: hierarchy.pairs[l-1].
# tests and docs. Internal code must not use them; it reads hierarchy/_layouts/_level_*. layouts = property(lambda self: tuple(self._layouts))
patch = property(lambda self: self.hierarchy.pairs[0]) boundaries = property(lambda self: tuple(self._boundaries),
child_patch = property(lambda self: self.hierarchy.pairs[1] if self.hierarchy.depth > 2 else None) doc="Host boundary geometry per level (None: no bodies/Jets).")
base_layout = property(lambda self: self._layouts[0]) jets = property(lambda self: tuple(self._level_jets))
fine_layout = property(lambda self: self._layouts[1]) leaf = property(lambda self: tuple(m.copy() for m in self._leaf),
child_layout = property(lambda self: self._layouts[2] if len(self._layouts) > 2 else None) doc="Per-level (ny, nx) masks of public (non-NaN) download cells.")
fine_boundary = property(lambda self: self._boundary_hosts[1])
child_boundary = property(lambda self: self._boundary_hosts[2] if len(self._boundary_hosts) > 2
else None)
_jets = property(lambda self: self._level_jets[1])
_boundaries = property(lambda self: {i: b for i, b in enumerate(self._level_boundaries) if i})
_fine_tables = property(lambda self: None if self._level_boundaries[1] is None
else self._level_boundaries[1].tables)
_actuated = property(lambda self: self._level_boundaries[1] is not None
and self._level_boundaries[1].actuated)
_jet_cache = property(lambda self: None if self._level_boundaries[1] is None
else self._level_boundaries[1].jet_cache)
_fine_leaf = property(lambda self: self._leaf[1])
_child_leaf = property(lambda self: self._leaf[2] if len(self._leaf) > 2 else None)
_middle_kind = property(lambda self: self.hierarchy.kinds[1] if self.hierarchy.depth > 2 else None)
# Preserve the original entry point and its periodic default.
PeriodicPatchStepper = FixedPatchStepper
+5 -3
View File
@@ -90,6 +90,7 @@ class Simulation:
@context_bound @context_bound
def stream(self) -> cuda.Stream: def stream(self) -> cuda.Stream:
"""Internal CUDA stream for async DTOH coordination.""" """Internal CUDA stream for async DTOH coordination."""
self._require_usable()
if self._stream is None: if self._stream is None:
self._stream = cuda.Stream() self._stream = cuda.Stream()
return self._stream return self._stream
@@ -174,9 +175,10 @@ class Simulation:
"""Change one jet's signed amplitude; re-uploads per-link velocity only. """Change one jet's signed amplitude; re-uploads per-link velocity only.
Invalid input raises before any change. Any error after validation Invalid input raises before any change. Any error after validation
leaves the host amplitude unchanged but the device column undefined, (stream creation or upload) leaves the host amplitude unchanged but the
so the Simulation is marked failed: every later call except close() device column undefined, so the Simulation is marked failed: every later
raises; close and recreate it. device-touching public method (incl. ``stream``) raises; close() still
works. Close and recreate it.
""" """
self._require_usable() self._require_usable()
jet = self.bodies.jets[index].with_amplitude(amplitude) jet = self.bodies.jets[index].with_amplitude(amplitude)
+1
View File
@@ -42,6 +42,7 @@ inline float __fdiv_rn(float a,float b){return a/b;}
#include "step/one_step_double.cu" #include "step/one_step_double.cu"
#else #else
#include "core/cell_kind.cuh" #include "core/cell_kind.cuh"
#include "operators/block_collision.cuh"
#include "step/block_channel.cu" #include "step/block_channel.cu"
#endif #endif
int main(){ int main(){
+29 -4
View File
@@ -10,6 +10,7 @@ import numpy as np
import pytest import pytest
pytest.importorskip("pycuda.driver") pytest.importorskip("pycuda.driver")
from CelerisLab import simulation as module
from CelerisLab.body.jets import Jet, Radial, Sector from CelerisLab.body.jets import Jet, Radial, Sector
from tests.unit.test_pinned_body_remap import ready from tests.unit.test_pinned_body_remap import ready
from tests.unit.test_simulation_lifecycle import make_simulation from tests.unit.test_simulation_lifecycle import make_simulation
@@ -68,6 +69,30 @@ def test_simulation_upload_failure_keeps_host_and_blocks_reuse(monkeypatch):
assert sim._closed and sim.ctx.is_closed assert sim._closed and sim.ctx.is_closed
def test_simulation_stream_creation_failure_fails_closed(monkeypatch):
sim, jet = simulation(monkeypatch)
created = []
def stream():
created.append(1)
raise RuntimeError("stream creation failed")
monkeypatch.setattr(module.cuda, "Stream", stream, raising=False)
sim._stream = None
with pytest.raises(RuntimeError, match="stream creation failed"):
sim.set_jet_amplitude(0, 2e-3)
assert sim.bodies.jets[0] is jet and sim._failed is not None
assert sim.field.curved.writes == [] and created == [1]
for call in (lambda: sim.stream, lambda: sim.set_jet_amplitude(0, 2e-3)):
with pytest.raises(RuntimeError, match="Simulation failed"):
call()
assert created == [1] # the public stream property does not retry creation
sim.ctx.detach_failure = True
with pytest.raises(RuntimeError, match="detach failed"):
sim.close()
assert not sim._closed
sim.close() # close stays retryable after the failure
assert sim._closed and sim.ctx.is_closed
def test_simulation_before_initialize_commits_host_only(monkeypatch): def test_simulation_before_initialize_commits_host_only(monkeypatch):
sim, _ = simulation(monkeypatch) sim, _ = simulation(monkeypatch)
sim._initialized = False sim._initialized = False
@@ -78,7 +103,7 @@ def test_simulation_before_initialize_commits_host_only(monkeypatch):
def stepper_state(run): def stepper_state(run):
boundary = run._level_boundaries[2] boundary = run._level_boundaries[2]
return boundary.tables["surface_velocity"].copy(), boundary.jets, run._level_jets[2] return boundary.tables["surface_velocity"].copy(), boundary.jets, run.jets[2]
def test_stepper_invalid_input_writes_nothing_and_stays_ready(): def test_stepper_invalid_input_writes_nothing_and_stays_ready():
@@ -92,10 +117,10 @@ def test_stepper_invalid_input_writes_nothing_and_stays_ready():
run.set_jet_amplitude(*args, **kw) run.set_jet_amplitude(*args, **kw)
module.cuda.memcpy_htod.assert_not_called() module.cuda.memcpy_htod.assert_not_called()
np.testing.assert_array_equal(run._level_boundaries[2].tables["surface_velocity"], velocity) np.testing.assert_array_equal(run._level_boundaries[2].tables["surface_velocity"], velocity)
assert (run._level_boundaries[2].jets, run._level_jets[2]) == (jets, level_jets) assert (run._level_boundaries[2].jets, run.jets[2]) == (jets, level_jets)
assert run._status == "ready" assert run._status == "ready"
run.set_jet_amplitude(0, .017, level=2) run.set_jet_amplitude(0, .017, level=2)
assert run._level_jets[2][0].amplitude == .017 assert run.jets[2][0].amplitude == .017
finally: finally:
run.close() run.close()
@@ -112,7 +137,7 @@ def test_stepper_device_failure_keeps_host_and_blocks_reuse(where):
with pytest.raises(RuntimeError, match="upload failed"): with pytest.raises(RuntimeError, match="upload failed"):
run.set_jet_amplitude(0, .017, level=2) run.set_jet_amplitude(0, .017, level=2)
np.testing.assert_array_equal(run._level_boundaries[2].tables["surface_velocity"], velocity) np.testing.assert_array_equal(run._level_boundaries[2].tables["surface_velocity"], velocity)
assert (run._level_boundaries[2].jets, run._level_jets[2]) == (jets, level_jets) assert (run._level_boundaries[2].jets, run.jets[2]) == (jets, level_jets)
assert run._status == "failed" assert run._status == "failed"
module.cuda.memcpy_htod.side_effect = None module.cuda.memcpy_htod.side_effect = None
run._stream.synchronize.side_effect = None run._stream.synchronize.side_effect = None
+8 -8
View File
@@ -94,13 +94,13 @@ class MultiplePinnedRemap(unittest.TestCase):
run.set_jet_amplitude(0, .017, level=2) run.set_jet_amplitude(0, .017, level=2)
run.set_jet_amplitude(1, -.011, level=2) run.set_jet_amplitude(1, -.011, level=2)
before = run.download() before = run.download()
clone = run.remap_child((PINS+MOVED+ADDED)[::-1], child_objects=bodies) clone = run.remap_finest((PINS+MOVED+ADDED)[::-1], objects=bodies)
try: try:
self.assertEqual(tuple(j.amplitude for j in clone._boundaries[2].jets), (.017, -.011)) self.assertEqual(tuple(j.amplitude for j in clone._level_boundaries[2].jets), (.017, -.011))
np.testing.assert_array_equal(clone._boundaries[2].tables['surface_velocity'], np.testing.assert_array_equal(clone._level_boundaries[2].tables['surface_velocity'],
run._boundaries[2].tables['surface_velocity']) run._level_boundaries[2].tables['surface_velocity'])
clone.set_jet_amplitude(1, .005, level=2) clone.set_jet_amplitude(1, .005, level=2)
self.assertEqual(run._boundaries[2].jets[1].amplitude, -.011) self.assertEqual(run._level_boundaries[2].jets[1].amplitude, -.011)
clone.step() clone.step()
self.assertEqual((run.step_count, clone.step_count), (7, 8)) self.assertEqual((run.step_count, clone.step_count), (7, 8))
for got, want in zip(run.download(), before): for got, want in zip(run.download(), before):
@@ -116,7 +116,7 @@ class MultiplePinnedRemap(unittest.TestCase):
start = len(allocations) start = len(allocations)
module.cuda.memcpy_dtoh.reset_mock() module.cuda.memcpy_dtoh.reset_mock()
with self.assertRaisesRegex(ValueError, 'boundary'): with self.assertRaisesRegex(ValueError, 'boundary'):
run.remap_child(PINS+FLUID, child_objects=changed) run.remap_finest(PINS+FLUID, objects=changed)
self.assertEqual(len(allocations), start) self.assertEqual(len(allocations), start)
module.cuda.memcpy_dtoh.assert_not_called() module.cuda.memcpy_dtoh.assert_not_called()
with self.assertRaisesRegex(ValueError, 'pinned'): with self.assertRaisesRegex(ValueError, 'pinned'):
@@ -128,7 +128,7 @@ class MultiplePinnedRemap(unittest.TestCase):
cases = [(PINS+FLUID, bodies[:1]), (PINS+FLUID, bodies[::-1]), cases = [(PINS+FLUID, bodies[:1]), (PINS+FLUID, bodies[::-1]),
(PINS+FLUID, (bodies[0], body(center=(32.1, 10.)))), (PINS+FLUID, (bodies[0], body(center=(32.1, 10.)))),
(PINS+FLUID, (bodies[0], body(SurfaceActuation(0., .7, 0.), (32., 10.)))), (PINS+FLUID, (bodies[0], body(SurfaceActuation(0., .7, 0.), (32., 10.)))),
(PINS+FLUID, run.child_boundary.objects)] (PINS+FLUID, run.boundaries[2].objects)]
for index in range(2): for index in range(2):
changed = list(PINS+FLUID) changed = list(PINS+FLUID)
x0, y0, x1, y1 = changed[index] x0, y0, x1, y1 = changed[index]
@@ -141,7 +141,7 @@ class MultiplePinnedRemap(unittest.TestCase):
run._stream.synchronize.reset_mock() run._stream.synchronize.reset_mock()
for rects, supplied in cases: for rects, supplied in cases:
with self.subTest(rects=rects), self.assertRaises(ValueError): with self.subTest(rects=rects), self.assertRaises(ValueError):
run.remap_child(rects, child_objects=supplied) run.remap_finest(rects, objects=supplied)
self.assertEqual(len(allocations), start) self.assertEqual(len(allocations), start)
module.cuda.memcpy_dtoh.assert_not_called() module.cuda.memcpy_dtoh.assert_not_called()
run._stream.synchronize.assert_not_called() run._stream.synchronize.assert_not_called()
+13 -13
View File
@@ -61,11 +61,11 @@ class PinnedBodyRemap(unittest.TestCase):
run,module,allocations,storage,objects=ready() run,module,allocations,storage,objects=ready()
run.set_jet_amplitude(0,.017,level=2) run.set_jet_amplitude(0,.017,level=2)
before=run.download() before=run.download()
clone=run.remap_child(B,child_objects=objects) clone=run.remap_finest(B,objects=objects)
self.assertEqual(clone._boundaries[2].jets[0].amplitude,.017) self.assertEqual(clone._level_boundaries[2].jets[0].amplitude,.017)
np.testing.assert_array_equal(clone._boundaries[2].tables['surface_velocity'], np.testing.assert_array_equal(clone._level_boundaries[2].tables['surface_velocity'],
run._boundaries[2].tables['surface_velocity']) run._level_boundaries[2].tables['surface_velocity'])
self.assertTrue(same_physical_boundary(run.child_boundary,clone.child_boundary)) self.assertTrue(same_physical_boundary(run.boundaries[2],clone.boundaries[2]))
clone.step() clone.step()
self.assertEqual((run.step_count,clone.step_count),(7,8)) self.assertEqual((run.step_count,clone.step_count),(7,8))
for got,want in zip(run.download(),before): np.testing.assert_array_equal(got,want) for got,want in zip(run.download(),before): np.testing.assert_array_equal(got,want)
@@ -73,22 +73,22 @@ class PinnedBodyRemap(unittest.TestCase):
def test_static_body_without_actuation(self): def test_static_body_without_actuation(self):
run,module,allocations,storage,objects=ready(jets=()) run,module,allocations,storage,objects=ready(jets=())
clone=run.remap_child(A,child_objects=objects) clone=run.remap_finest(A,objects=objects)
self.assertFalse(clone._boundaries[2].actuated) self.assertFalse(clone._level_boundaries[2].actuated)
self.assertNotIn('surface_velocity',clone._boundaries[2].tables) self.assertNotIn('surface_velocity',clone._level_boundaries[2].tables)
for a,b in zip(run.download(),clone.download()): np.testing.assert_array_equal(a,b) for a,b in zip(run.download(),clone.download()): np.testing.assert_array_equal(a,b)
clone.close();run.close() clone.close();run.close()
def test_reject_before_device_reads_or_allocation(self): def test_reject_before_device_reads_or_allocation(self):
run,module,allocations,storage,objects=ready() run,module,allocations,storage,objects=ready()
start=len(allocations) start=len(allocations)
cases=[(B,()),(B,run.child_boundary.objects),(B,(body(center=(12.1,14.)),)), cases=[(B,()),(B,run.boundaries[2].objects),(B,(body(center=(12.1,14.)),)),
(((24,16,36,40),B[1]),objects),(((16,28,32,40),B[1]),objects), (((24,16,36,40),B[1]),objects),(((16,28,32,40),B[1]),objects),
((B[1],),objects)] ((B[1],),objects)]
module.cuda.memcpy_dtoh.reset_mock() module.cuda.memcpy_dtoh.reset_mock()
for rects,objs in cases: for rects,objs in cases:
with self.subTest(rects=rects),self.assertRaises(ValueError): with self.subTest(rects=rects),self.assertRaises(ValueError):
run.remap_child(rects,child_objects=objs) run.remap_finest(rects,objects=objs)
self.assertEqual(len(allocations),start) self.assertEqual(len(allocations),start)
module.cuda.memcpy_dtoh.assert_not_called() module.cuda.memcpy_dtoh.assert_not_called()
self.assertEqual((run._status,run.step_count),('ready',7)) self.assertEqual((run._status,run.step_count),('ready',7))
@@ -100,9 +100,9 @@ class PinnedBodyRemap(unittest.TestCase):
(body(SurfaceActuation(0.,.7,speed)),),()) (body(SurfaceActuation(0.,.7,speed)),),())
for spec in (SurfaceActuation(0.,.7,.003),SurfaceActuation(2.,.7,speed)): for spec in (SurfaceActuation(0.,.7,.003),SurfaceActuation(2.,.7,speed)):
module.cuda.memcpy_dtoh.reset_mock() module.cuda.memcpy_dtoh.reset_mock()
with self.assertRaises(ValueError): run.remap_child(B,child_objects=(body(spec),)) with self.assertRaises(ValueError): run.remap_finest(B,objects=(body(spec),))
module.cuda.memcpy_dtoh.assert_not_called() module.cuda.memcpy_dtoh.assert_not_called()
clone=run.remap_child(B,child_objects=objects) clone=run.remap_finest(B,objects=objects)
clone.close();run.close() clone.close();run.close()
def test_allocation_failure_keeps_old_ready_and_count(self): def test_allocation_failure_keeps_old_ready_and_count(self):
@@ -117,7 +117,7 @@ class PinnedBodyRemap(unittest.TestCase):
return allocate(size) return allocate(size)
with patch.object(module.cuda,'mem_alloc',side_effect=fail): with patch.object(module.cuda,'mem_alloc',side_effect=fail):
with self.assertRaisesRegex(RuntimeError,'allocation'): with self.assertRaisesRegex(RuntimeError,'allocation'):
run.remap_child(B,child_objects=objects) run.remap_finest(B,objects=objects)
self.assertEqual((run._status,run.step_count),('ready',7)) self.assertEqual((run._status,run.step_count),('ready',7))
for got,want in zip(run.download(),before): np.testing.assert_array_equal(got,want) for got,want in zip(run.download(),before): np.testing.assert_array_equal(got,want)
for a in saved: a.free.assert_not_called() for a in saved: a.free.assert_not_called()
+2 -2
View File
@@ -164,7 +164,7 @@ class ProjectedPairDispatch(unittest.TestCase):
run, host, _, _, storage = backend(refinements=FOUR[:2], explosion=explosion) run, host, _, _, storage = backend(refinements=FOUR[:2], explosion=explosion)
run._step_count = 7 run._step_count = 7
old = run.download() old = run.download()
clone = run.remap_child(FOUR[1]) clone = run.remap_finest(FOUR[1])
self.assertEqual(clone._explosion, explosion) self.assertEqual(clone._explosion, explosion)
self.assertEqual(clone.step_count, 7) self.assertEqual(clone.step_count, 7)
for pair in clone._pairs: for pair in clone._pairs:
@@ -182,7 +182,7 @@ class ProjectedPairDispatch(unittest.TestCase):
refinements=FOUR[:2], explosion='projected_linear') refinements=FOUR[:2], explosion='projected_linear')
run._step_count = 3 run._step_count = 3
before = run.download() before = run.download()
clone = run.remap_child(changed) clone = run.remap_finest(changed)
self.assertEqual(clone._explosion, 'projected_linear') self.assertEqual(clone._explosion, 'projected_linear')
self.assertEqual(clone.hierarchy.pairs[1].bounds, tuple(changed)) self.assertEqual(clone.hierarchy.pairs[1].bounds, tuple(changed))
self.assertNotEqual(clone._pairs[1].n_ghosts, run._pairs[1].n_ghosts) self.assertNotEqual(clone._pairs[1].n_ghosts, run._pairs[1].n_ghosts)
+5
View File
@@ -524,6 +524,11 @@ def test_public_jet_reversal_boundary_and_windows():
assert r.jet_interval(handle, "reverse", switch + 80, 80, diameter) == (80, None) assert r.jet_interval(handle, "reverse", switch + 80, 80, diameter) == (80, None)
assert r.jet_interval(handle, "fixed", switch, 80, diameter) == (80, None) assert r.jet_interval(handle, "fixed", switch, 80, diameter) == (80, None)
assert not calls assert not calls
refined_calls = []
refined = SimpleNamespace(set_jet_amplitude=lambda i, a, *, level: refined_calls.append((i, a, level)))
n, event = r.jet_interval(refined, "reverse", switch, 80, diameter, level=1)
assert n == 80 and event["step"] == switch and event["tstar"] == 300
assert refined_calls == [(0, -.001, 1), (1, .001, 1)]
rows = [dict(tstar=t, tstar_force=t - .5, cd=2.0, cl=-.1) for t in range(251, 306)] rows = [dict(tstar=t, tstar_force=t - .5, cd=2.0, cl=-.1) for t in range(251, 306)]
windows = r.jet_diagnostics(rows, "reverse") windows = r.jet_diagnostics(rows, "reverse")
assert len(windows) == 6 and all(w["complete"] for w in windows) assert len(windows) == 6 and all(w["complete"] for w in windows)
+2 -2
View File
@@ -177,7 +177,7 @@ class FixedPatchSetTests(unittest.TestCase):
self.assertEqual((fy, fx), (32, 48)) self.assertEqual((fy, fx), (32, 48))
rows = cells = 0 rows = cells = 0
for p in s.patches: for p in s.patches:
y0, x0 = p.fine_level_origin y0, x0 = p.fine_level_offset
ly, lx = p.fine_shape ly, lx = p.fine_shape
np.testing.assert_array_equal(s.fine_kind[y0:y0 + ly, x0:x0 + lx], p.fine_kind) np.testing.assert_array_equal(s.fine_kind[y0:y0 + ly, x0:x0 + lx], p.fine_kind)
gy, gx = np.divmod(p.ghost_fine.astype(np.int64), lx) gy, gx = np.divmod(p.ghost_fine.astype(np.int64), lx)
@@ -224,7 +224,7 @@ class FixedPatchSetTests(unittest.TestCase):
layout = s.fine_layout((4, 4)) layout = s.fine_layout((4, 4))
support = np.zeros(s.fine_shape, bool) support = np.zeros(s.fine_shape, bool)
for p in s.patches: for p in s.patches:
y0, x0 = p.fine_level_origin y0, x0 = p.fine_level_offset
support[y0:y0 + p.fine_shape[0], x0:x0 + p.fine_shape[1]] = True support[y0:y0 + p.fine_shape[0], x0:x0 + p.fine_shape[1]] = True
np.testing.assert_array_equal(layout.owned, support) np.testing.assert_array_equal(layout.owned, support)
self.assertTrue(layout.owned.ravel()[s.ghost_fine.astype(np.intp)].all()) self.assertTrue(layout.owned.ravel()[s.ghost_fine.astype(np.intp)].all())
+5 -2
View File
@@ -31,7 +31,10 @@ class RefinementMRTHost(unittest.TestCase):
"compute_feq(rho, ux, uy, feq);", "compute_feq(rho, ux, uy, feq);",
"collide_srt_no_force(f, feq, d_params.omega);"): "collide_srt_no_force(f, feq, d_params.omega);"):
self.assertIn(text, block_collision) self.assertIn(text, block_collision)
self.assertIn("collide_srt(f, feq, forcing, omega);", channel) self.assertIn("collide_srt(f, feq, forcing, omega);", block_collision)
self.assertIn("#elif BLOCK_ROOT_BC == 2", block_collision)
self.assertEqual(channel.count("block_collide(f);"), 1)
self.assertNotIn("collide_mrt_unforced(f);", channel)
for source in (blocks,): for source in (blocks,):
self.assertIn("(COLLISION_MODEL != 0 && COLLISION_MODEL != 2)", source) self.assertIn("(COLLISION_MODEL != 0 && COLLISION_MODEL != 2)", source)
self.assertIn("USE_LES", source) self.assertIn("USE_LES", source)
@@ -40,7 +43,7 @@ class RefinementMRTHost(unittest.TestCase):
# One shared unforced MRT body; each blocked/patch path calls it exactly once. # One shared unforced MRT body; each blocked/patch path calls it exactly once.
mrt = (KERNELS / "operators/collision_mrt.cuh").read_text() mrt = (KERNELS / "operators/collision_mrt.cuh").read_text()
self.assertEqual(mrt.count("void collide_mrt_unforced("), 1) self.assertEqual(mrt.count("void collide_mrt_unforced("), 1)
for source in (blocks + "\n" + block_collision, channel): for source in (blocks + "\n" + block_collision, channel + "\n" + block_collision):
self.assertEqual(source.count("collide_mrt_unforced(f);"), 1) self.assertEqual(source.count("collide_mrt_unforced(f);"), 1)
self.assertNotIn("collide_mrt_no_force(", source) self.assertNotIn("collide_mrt_no_force(", source)
self.assertNotIn("compute_rho_u(f, rho, ux, uy);\n collide_mrt", source) self.assertNotIn("compute_rho_u(f, rho, ux, uy);\n collide_mrt", source)
+10 -10
View File
@@ -29,7 +29,7 @@ class RuntimeChildRemap(unittest.TestCase):
saved = {key: value.copy() for key, value in storage.items()} saved = {key: value.copy() for key, value in storage.items()}
old_allocations = list(allocations) old_allocations = list(allocations)
module.cuda.memcpy_dtoh.reset_mock() module.cuda.memcpy_dtoh.reset_mock()
clone = run.remap_child(A) clone = run.remap_finest(A)
self.assertEqual(module.cuda.memcpy_dtoh.call_count, 3) self.assertEqual(module.cuda.memcpy_dtoh.call_count, 3)
for a, b in zip(clone.download(), before): for a, b in zip(clone.download(), before):
np.testing.assert_array_equal(a, b) np.testing.assert_array_equal(a, b)
@@ -48,7 +48,7 @@ class RuntimeChildRemap(unittest.TestCase):
run, module, allocations, storage = ready_backend() run, module, allocations, storage = ready_backend()
before = run.download() before = run.download()
expected = remap(topology(A), topology(B), before).g expected = remap(topology(A), topology(B), before).g
clone = run.remap_child(B) clone = run.remap_finest(B)
for level, want, got in zip(clone._levels, expected, clone.download()): for level, want, got in zip(clone._levels, expected, clone.download()):
np.testing.assert_array_equal(got, want) np.testing.assert_array_equal(got, want)
for ptr in (level.current, level.spare): for ptr in (level.current, level.spare):
@@ -65,17 +65,17 @@ class RuntimeChildRemap(unittest.TestCase):
module.cuda.memcpy_dtoh.reset_mock() module.cuda.memcpy_dtoh.reset_mock()
for status in ('new', 'failed', 'closed', 'advancing'): for status in ('new', 'failed', 'closed', 'advancing'):
run._status = status run._status = status
with self.assertRaises(RuntimeError): run.remap_child(B) with self.assertRaises(RuntimeError): run.remap_finest(B)
run._status = 'ready' run._status = 'ready'
for bad in ([], [(40,16,48,40)], [(16,16,32,24),(16,27,32,55)]): for bad in ([], [(40,16,48,40)], [(16,16,32,24),(16,27,32,55)]):
with self.assertRaises(ValueError): run.remap_child(bad) with self.assertRaises(ValueError): run.remap_finest(bad)
hosts, pairs = run._boundary_hosts, run.hierarchy.pairs hosts, pairs = run._boundaries, run.hierarchy.pairs
for attribute, value in (('_wall_aware', True), for attribute, value in (('_wall_aware', True),
('_boundary_hosts', [*hosts[:2], object()]), # child boundary ('_boundaries', [*hosts[:2], object()]), # child boundary
('_boundary_hosts', [hosts[0], object(), hosts[2]]), # fine boundary ('_boundaries', [hosts[0], object(), hosts[2]]), # fine boundary
('hierarchy', replace(run.hierarchy, pairs=(pairs[0], None)))): ('hierarchy', replace(run.hierarchy, pairs=(pairs[0], None)))):
with patch.object(run, attribute, value): with patch.object(run, attribute, value):
with self.assertRaises(ValueError): run.remap_child(B) with self.assertRaises(ValueError): run.remap_finest(B)
self.assertEqual(len(allocations), start) self.assertEqual(len(allocations), start)
module.cuda.memcpy_dtoh.assert_not_called() module.cuda.memcpy_dtoh.assert_not_called()
run.close() run.close()
@@ -102,7 +102,7 @@ class RuntimeChildRemap(unittest.TestCase):
initialize(candidate, *arrays) initialize(candidate, *arrays)
operation = patch.object(type(run), 'initialize', fail) operation = patch.object(type(run), 'initialize', fail)
with operation, self.assertRaisesRegex(RuntimeError, 'candidate'): with operation, self.assertRaisesRegex(RuntimeError, 'candidate'):
run.remap_child(B) run.remap_finest(B)
self.assertEqual((run._status, run.step_count), ('ready', 7)) self.assertEqual((run._status, run.step_count), ('ready', 7))
for a, b in zip(run.download(), before): np.testing.assert_array_equal(a, b) for a, b in zip(run.download(), before): np.testing.assert_array_equal(a, b)
for allocation in old_allocations: allocation.free.assert_not_called() for allocation in old_allocations: allocation.free.assert_not_called()
@@ -116,7 +116,7 @@ class RuntimeChildRemap(unittest.TestCase):
target = run._stream.synchronize if operation == 'sync' else module.cuda.memcpy_dtoh target = run._stream.synchronize if operation == 'sync' else module.cuda.memcpy_dtoh
target.side_effect = RuntimeError('old device failure') target.side_effect = RuntimeError('old device failure')
with self.assertRaisesRegex(RuntimeError, 'old device'): with self.assertRaisesRegex(RuntimeError, 'old device'):
run.remap_child(B) run.remap_finest(B)
self.assertEqual(run._status, 'failed') self.assertEqual(run._status, 'failed')
self.assertEqual(len(allocations), start) self.assertEqual(len(allocations), start)
target.side_effect = None target.side_effect = None
+3 -3
View File
@@ -77,13 +77,13 @@ class ShiftedJets(unittest.TestCase):
np.testing.assert_array_equal(boundary.tables['surface_velocity'],-original) np.testing.assert_array_equal(boundary.tables['surface_velocity'],-original)
self.assertEqual(run.step_count,7) self.assertEqual(run.step_count,7)
for a,b in zip(before.values,run._download_raw().values):self.assertEqual(a.tobytes(),b.tobytes()) 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) clone=run.remap_finest(B,objects=objects)
self.assertEqual(clone._level_jets[2][0].amplitude,-jet().amplitude) self.assertEqual(clone.jets[2][0].amplitude,-jet().amplitude)
np.testing.assert_array_equal(clone._level_boundaries[2].tables['surface_velocity'], np.testing.assert_array_equal(clone._level_boundaries[2].tables['surface_velocity'],
boundary.tables['surface_velocity']) boundary.tables['surface_velocity'])
self.assertEqual(clone.step_count,7) self.assertEqual(clone.step_count,7)
saved=run._download_raw() saved=run._download_raw()
with self.assertRaises(ValueError):run.remap_child(B,child_objects=(body(center=(12.1,14.)),)) with self.assertRaises(ValueError):run.remap_finest(B,objects=(body(center=(12.1,14.)),))
for a,b in zip(saved.values,run._download_raw().values):self.assertEqual(a.tobytes(),b.tobytes()) for a,b in zip(saved.values,run._download_raw().values):self.assertEqual(a.tobytes(),b.tobytes())
finally: finally:
if clone is not None:clone.close() if clone is not None:clone.close()
+2 -2
View File
@@ -95,7 +95,7 @@ class ShiftedStorage(unittest.TestCase):
before = run._download_raw() before = run._download_raw()
expected = remap(topology(A),topology(B),before.values) expected = remap(topology(A),topology(B),before.values)
with patch.object(run,'download',side_effect=AssertionError('public decode used')): with patch.object(run,'download',side_effect=AssertionError('public decode used')):
clone = run.remap_child(B) clone = run.remap_finest(B)
self.assertEqual((run.step_count,clone.step_count),(7,7)) self.assertEqual((run.step_count,clone.step_count),(7,7))
self.assertFalse(clone._exchange_valid) self.assertFalse(clone._exchange_valid)
for level, want, got in zip(clone._levels,expected.g,clone._download_raw().values): for level, want, got in zip(clone._levels,expected.g,clone._download_raw().values):
@@ -114,7 +114,7 @@ class ShiftedStorage(unittest.TestCase):
physical = tuple(np.broadcast_to(_WEIGHTS, level.shape).copy() for level in run._levels) physical = tuple(np.broadcast_to(_WEIGHTS, level.shape).copy() for level in run._levels)
for a in physical: a[1] += 2**-35 for a in physical: a[1] += 2**-35
run.initialize(*physical) run.initialize(*physical)
for a,mask in zip(run._download_raw().values,run._leaf): for a,mask in zip(run._download_raw().values,run.leaf):
self.assertTrue((a[1,mask] == np.float32(2**-35)).all()) self.assertTrue((a[1,mask] == np.float32(2**-35)).all())
self.assertEqual(run.step_count,0) self.assertEqual(run.step_count,0)
module.cuda.memcpy_htod.reset_mock() module.cuda.memcpy_htod.reset_mock()
+3 -3
View File
@@ -35,10 +35,10 @@ class UnionBodySupport(unittest.TestCase):
run,module,allocations,storage,objects=ready() run,module,allocations,storage,objects=ready()
try: try:
run.set_jet_amplitude(0,.017,level=2) run.set_jet_amplitude(0,.017,level=2)
clone=run.remap_child(((16,16,24,40),(24,16,32,40),A[1]),child_objects=objects) clone=run.remap_finest(((16,16,24,40),(24,16,32,40),A[1]),objects=objects)
try: try:
self.assertEqual(clone._boundaries[2].jets[0].amplitude,.017) self.assertEqual(clone._level_boundaries[2].jets[0].amplitude,.017)
self.assertTrue(same_physical_boundary(run.child_boundary,clone.child_boundary)) self.assertTrue(same_physical_boundary(run.boundaries[2],clone.boundaries[2]))
for a,b in zip(run.download(),clone.download()):np.testing.assert_array_equal(a,b) for a,b in zip(run.download(),clone.download()):np.testing.assert_array_equal(a,b)
finally:clone.close() finally:clone.close()
finally:run.close() finally:run.close()
+1 -1
View File
@@ -40,7 +40,7 @@ class Equivalence(unittest.TestCase):
table_sets(FixedPatchSet2D.from_config(CFG, [(8, 8, 20, 22)]))) table_sets(FixedPatchSet2D.from_config(CFG, [(8, 8, 20, 22)])))
def test_separated_rectangles_reproduce_current_tables(self): def test_separated_rectangles_reproduce_current_tables(self):
rects = [(4, 4, 10, 10), (16, 16, 26, 26)] # gap 6 >= MIN_PATCH_GAP rects = [(4, 4, 10, 10), (16, 16, 26, 26)] # separated rectangles; union dependencies determine legality
self.assertEqual(table_sets(union(*rects)), self.assertEqual(table_sets(union(*rects)),
table_sets(FixedPatchSet2D.from_config(CFG, rects))) table_sets(FixedPatchSet2D.from_config(CFG, rects)))
@@ -113,7 +113,7 @@ class ChannelRootGPU(unittest.TestCase):
arrays = [state(*level.shape[1:]) for level in backend._levels] arrays = [state(*level.shape[1:]) for level in backend._levels]
backend.initialize(*arrays) backend.initialize(*arrays)
backend.step(4) backend.step(4)
if bodies: backend.set_jet_amplitude(0, -.001) if bodies: backend.set_jet_amplitude(0, -.001, level=1)
backend.step(4) backend.step(4)
raw = backend._download_raw() raw = backend._download_raw()
physical = backend.download() physical = backend.download()
@@ -128,7 +128,7 @@ class ChannelRootGPU(unittest.TestCase):
else: else:
for a,b in zip(reference, raw.values): np.testing.assert_array_equal(a,b) for a,b in zip(reference, raw.values): np.testing.assert_array_equal(a,b)
with self.assertRaisesRegex(ValueError, "periodic"): with self.assertRaisesRegex(ValueError, "periodic"):
backend.remap_child([]) backend.remap_finest([])
if __name__ == "__main__": unittest.main() if __name__ == "__main__": unittest.main()
@@ -66,31 +66,31 @@ class ExchangeWindowCPU(unittest.TestCase):
arrays = [np.full(level.shape, 1 / 9, np.float32) for level in backend._levels] arrays = [np.full(level.shape, 1 / 9, np.float32) for level in backend._levels]
backend.initialize(*arrays) backend.initialize(*arrays)
with self.assertRaisesRegex(RuntimeError, "no valid exchange window"): with self.assertRaisesRegex(RuntimeError, "no valid exchange window"):
backend.read_exchange_window(0) backend.read_exchange_window(0, level=1)
backend.step(1) backend.step(1)
self.assertEqual(len(accumulations(m)), 2) # closed window: interval ledger only self.assertEqual(len(accumulations(m)), 2) # closed window: interval ledger only
backend.start_exchange_window() backend.start_exchange_window()
with self.assertRaisesRegex(RuntimeError, "at least one completed"): with self.assertRaisesRegex(RuntimeError, "at least one completed"):
backend.read_exchange_window(0) backend.read_exchange_window(0, level=1)
m["calls"].clear() m["calls"].clear()
backend.step(2) backend.step(2)
calls = accumulations(m) calls = accumulations(m)
self.assertEqual(len(calls), 8) # two fine hops x two ledgers x two base steps self.assertEqual(len(calls), 8) # two fine hops x two ledgers x two base steps
self.assertEqual({id(c[1][1]) for c in calls}, self.assertEqual({id(c[1][1]) for c in calls},
{id(b.interval_integral), id(b.window_integral)}) {id(b.interval_integral), id(b.window_integral)})
out = backend.read_exchange_window(0) out = backend.read_exchange_window(0, level=1)
self.assertEqual((out["time_start"], out["time_end"], out["sample_count"]), (1., 3, 4)) self.assertEqual((out["time_start"], out["time_end"], out["sample_count"]), (1., 3, 4))
self.assertEqual(backend.read_exchange_window(0)["time_end"], 3) # non-destructive self.assertEqual(backend.read_exchange_window(0, level=1)["time_end"], 3) # non-destructive
backend.initialize(*arrays) # invalidates backend.initialize(*arrays) # invalidates
with self.assertRaisesRegex(RuntimeError, "no valid exchange window"): with self.assertRaisesRegex(RuntimeError, "no valid exchange window"):
backend.read_exchange_window(0) backend.read_exchange_window(0, level=1)
backend.start_exchange_window() backend.start_exchange_window()
with patch.object(backend, "_hop", side_effect=RuntimeError("injected")): with patch.object(backend, "_hop", side_effect=RuntimeError("injected")):
with self.assertRaises(RuntimeError): with self.assertRaises(RuntimeError):
backend.step(1) backend.step(1)
self.assertFalse(backend._window_valid) self.assertFalse(backend._window_valid)
with self.assertRaises(RuntimeError): with self.assertRaises(RuntimeError):
backend.read_exchange_window(0) # failed status backend.read_exchange_window(0, level=1) # failed status
def test_disabled_default_has_no_window_cost(self): def test_disabled_default_has_no_window_cost(self):
with mocked_backend() as m, ExitStack() as stack: with mocked_backend() as m, ExitStack() as stack:
@@ -140,12 +140,12 @@ class ExchangeWindowGPU(unittest.TestCase):
summed = {part: [0., 0., 0., 0.] for part in ("total", "wall", "aperture")} summed = {part: [0., 0., 0., 0.] for part in ("total", "wall", "aperture")}
steps = 0 steps = 0
for speed in (0., .002, .002, -.001, 0., 0.): # amplitude updates stay in-window for speed in (0., .002, .002, -.001, 0., 0.): # amplitude updates stay in-window
windowed.set_jet_amplitude(0, speed) windowed.set_jet_amplitude(0, speed, level=1)
plain.set_jet_amplitude(0, speed) plain.set_jet_amplitude(0, speed, level=1)
windowed.step(1) windowed.step(1)
plain.step(1) plain.step(1)
steps += 1 steps += 1
obs = windowed.read_raw_exchange(0) obs = windowed.read_raw_exchange(0, level=1)
for part in summed: for part in summed:
o = obs[part] o = obs[part]
summed[part] = [summed[part][0] + o["impulse"][0], summed[part][1] + o["impulse"][1], summed[part] = [summed[part][0] + o["impulse"][0], summed[part][1] + o["impulse"][1],
@@ -153,7 +153,7 @@ class ExchangeWindowGPU(unittest.TestCase):
summed[part][3] + o["mass_into_fluid"]] summed[part][3] + o["mass_into_fluid"]]
for a, e in zip(windowed._download_raw().values, plain._download_raw().values): for a, e in zip(windowed._download_raw().values, plain._download_raw().values):
same_bits(self, a, e) same_bits(self, a, e)
out = windowed.read_exchange_window(0) out = windowed.read_exchange_window(0, level=1)
self.assertEqual((out["time_end"] - out["time_start"], out["sample_count"]), self.assertEqual((out["time_end"] - out["time_start"], out["sample_count"]),
(steps, 2 * steps)) (steps, 2 * steps))
self.assertEqual(substeps, 2 * steps) self.assertEqual(substeps, 2 * steps)
@@ -165,11 +165,11 @@ class ExchangeWindowGPU(unittest.TestCase):
diff = np.abs(got - np.array(summed[part])) diff = np.abs(got - np.array(summed[part]))
self.assertTrue((diff <= bound).all(), f"{part}: {diff} > {bound}") self.assertTrue((diff <= bound).all(), f"{part}: {diff} > {bound}")
np.testing.assert_array_equal(o["force"], o["impulse"] / steps) np.testing.assert_array_equal(o["force"], o["impulse"] / steps)
again = windowed.read_exchange_window(0) again = windowed.read_exchange_window(0, level=1)
np.testing.assert_array_equal(again["total"]["impulse"], out["total"]["impulse"]) np.testing.assert_array_equal(again["total"]["impulse"], out["total"]["impulse"])
windowed.initialize(*arrays) windowed.initialize(*arrays)
with self.assertRaises(RuntimeError): with self.assertRaises(RuntimeError):
windowed.read_exchange_window(0) windowed.read_exchange_window(0, level=1)
def audit_substeps(stepper, b, count, record): def audit_substeps(stepper, b, count, record):
@@ -220,7 +220,7 @@ class ExchangeWindowShiftedMRTGPU(unittest.TestCase):
for _ in range(steps): for _ in range(steps):
windowed.step(1) windowed.step(1)
plain.step(1) plain.step(1)
summed += components(windowed.read_raw_exchange(0)) summed += components(windowed.read_raw_exchange(0, level=1))
for a, e in zip(windowed._download_raw().values, plain._download_raw().values): for a, e in zip(windowed._download_raw().values, plain._download_raw().values):
same_bits(self, a, e) same_bits(self, a, e)
split.step(steps) # one call: same ledger arithmetic as five step(1) calls split.step(steps) # one call: same ledger arithmetic as five step(1) calls
@@ -231,28 +231,28 @@ class ExchangeWindowShiftedMRTGPU(unittest.TestCase):
cuda.memcpy_dtoh(rows[1], split._level_boundaries[1].window_integral) cuda.memcpy_dtoh(rows[1], split._level_boundaries[1].window_integral)
np.testing.assert_array_equal(rows[0].view(np.uint64), rows[1].view(np.uint64)) np.testing.assert_array_equal(rows[0].view(np.uint64), rows[1].view(np.uint64))
out = windowed.read_exchange_window(0) out = windowed.read_exchange_window(0, level=1)
self.assertEqual(record["substeps"], 2 * steps) self.assertEqual(record["substeps"], 2 * steps)
self.assertEqual((out["time_end"] - out["time_start"], out["sample_count"]), self.assertEqual((out["time_end"] - out["time_start"], out["sample_count"]),
(steps, 2 * steps)) (steps, 2 * steps))
bound = 4 * (count * record["substeps"]) * EPS * record["scale"] + 1e-300 bound = 4 * (count * record["substeps"]) * EPS * record["scale"] + 1e-300
diff = np.abs(components(out) - summed) diff = np.abs(components(out) - summed)
self.assertTrue((diff <= bound).all(), f"{diff} > {bound}") self.assertTrue((diff <= bound).all(), f"{diff} > {bound}")
np.testing.assert_array_equal(split.read_exchange_window(0)["total"]["impulse"], np.testing.assert_array_equal(split.read_exchange_window(0, level=1)["total"]["impulse"],
out["total"]["impulse"]) out["total"]["impulse"])
windowed.step(0) # no interval: window and its duration unchanged windowed.step(0) # no interval: window and its duration unchanged
again = windowed.read_exchange_window(0) again = windowed.read_exchange_window(0, level=1)
self.assertEqual(again["time_end"], out["time_end"]) self.assertEqual(again["time_end"], out["time_end"])
np.testing.assert_array_equal(again["total"]["impulse"], out["total"]["impulse"]) np.testing.assert_array_equal(again["total"]["impulse"], out["total"]["impulse"])
windowed.start_exchange_window() # reset after a window windowed.start_exchange_window() # reset after a window
with self.assertRaisesRegex(RuntimeError, "at least one completed"): with self.assertRaisesRegex(RuntimeError, "at least one completed"):
windowed.read_exchange_window(0) windowed.read_exchange_window(0, level=1)
record.update(scale=np.zeros(4), substeps=0) record.update(scale=np.zeros(4), substeps=0)
windowed.step(1) windowed.step(1)
single = components(windowed.read_raw_exchange(0)) single = components(windowed.read_raw_exchange(0, level=1))
reset = windowed.read_exchange_window(0) reset = windowed.read_exchange_window(0, level=1)
self.assertEqual(reset["time_end"] - reset["time_start"], 1) self.assertEqual(reset["time_end"] - reset["time_start"], 1)
bound = 4 * (count * record["substeps"]) * EPS * record["scale"] + 1e-300 bound = 4 * (count * record["substeps"]) * EPS * record["scale"] + 1e-300
self.assertTrue((np.abs(components(reset) - single) <= bound).all()) self.assertTrue((np.abs(components(reset) - single) <= bound).all())
@@ -81,7 +81,7 @@ def run_blocked_patch(speed):
level_objects=[(), (source,)], level_objects=[(), (source,)],
observe_raw_exchange=True) as backend: observe_raw_exchange=True) as backend:
c, f = backend._levels c, f = backend._levels
flags = backend.fine_boundary.flags.reshape(40, 40) flags = backend.boundaries[1].flags.reshape(40, 40)
active = np.zeros((40, 40), bool) active = np.zeros((40, 40), bool)
active[4:-4, 4:-4] = True active[4:-4, 4:-4] = True
active &= ((flags & FLUID) != 0) & ((flags & SOLID) == 0) active &= ((flags & FLUID) != 0) & ((flags & SOLID) == 0)
@@ -101,7 +101,7 @@ def run_blocked_patch(speed):
with patch.object(f, "kernel", fine): with patch.object(f, "kernel", fine):
backend.step(2) backend.step(2)
assert len(traces) == 4 and backend.step_count == 2 assert len(traces) == 4 and backend.step_count == 2
mass = backend.read_raw_exchange(0)["total"]["mass_into_fluid"] mass = backend.read_raw_exchange(0, level=1)["total"]["mass_into_fluid"]
return traces, roi, backend.download(), float(mass) return traces, roi, backend.download(), float(mass)
@@ -56,8 +56,8 @@ def make(context):
def topology(run, rectangles, objects): def topology(run, rectangles, objects):
return RemapTopology.from_config( return RemapTopology.from_config(
CFG, run.patch, rectangles, child_objects=objects, CFG, run.hierarchy.pairs[0], rectangles, child_objects=objects,
child_jets=run._boundaries[2].jets, child_jets=run._level_boundaries[2].jets,
) )
@@ -74,8 +74,8 @@ class MultiplePinnedRemapGPU(unittest.TestCase):
def change_amplitudes(self, run, amplitudes): def change_amplitudes(self, run, amplitudes):
# Each jet selects a different body. Updating one must preserve the # Each jet selects a different body. Updating one must preserve the
# other body's current per-link velocities, including its driven links. # other body's current per-link velocities, including its driven links.
boundary = run._boundaries[2] boundary = run._level_boundaries[2]
body_ids = run.child_boundary.curved_host[3] body_ids = run.boundaries[2].curved_host[3]
for index, amplitude in enumerate(amplitudes): for index, amplitude in enumerate(amplitudes):
before = boundary.tables['surface_velocity'].copy() before = boundary.tables['surface_velocity'].copy()
run.set_jet_amplitude(index, amplitude, level=2) run.set_jet_amplitude(index, amplitude, level=2)
@@ -87,11 +87,11 @@ class MultiplePinnedRemapGPU(unittest.TestCase):
def check_boundary(self, old, new, mesh, amplitudes): def check_boundary(self, old, new, mesh, amplitudes):
self.assertEqual(mesh.pinned_rectangles, PINS) self.assertEqual(mesh.pinned_rectangles, PINS)
self.assertTrue(same_physical_boundary(old.child_boundary, new.child_boundary)) self.assertTrue(same_physical_boundary(old.boundaries[2], new.boundaries[2]))
self.assertTrue(same_physical_boundary(new.child_boundary, mesh.boundary)) self.assertTrue(same_physical_boundary(new.boundaries[2], mesh.boundary))
# Compare physical cut-link rows, never layout-dependent blocked addresses. # Compare physical cut-link rows, never layout-dependent blocked addresses.
for run in (old, new): for run in (old, new):
boundary = run._boundaries[2] boundary = run._level_boundaries[2]
self.assertEqual(tuple(jet.amplitude for jet in boundary.jets), amplitudes) self.assertEqual(tuple(jet.amplitude for jet in boundary.jets), amplitudes)
exact(self, boundary.tables['surface_velocity'], mesh.surface_velocity, exact(self, boundary.tables['surface_velocity'], mesh.surface_velocity,
'latest physical per-link velocities') 'latest physical per-link velocities')
@@ -107,7 +107,7 @@ class MultiplePinnedRemapGPU(unittest.TestCase):
before, count = old.download(), old.step_count before, count = old.download(), old.step_count
mesh = topology(old, START, objects) mesh = topology(old, START, objects)
# Input order is irrelevant even with several pins and fluid regions. # Input order is irrelevant even with several pins and fluid regions.
new = owned.enter_context(old.remap_child(START[::-1], child_objects=objects)) new = owned.enter_context(old.remap_finest(START[::-1], objects=objects))
self.assertIsNot(old, new) self.assertIsNot(old, new)
self.assertEqual((old.step_count, new.step_count), (count, count)) self.assertEqual((old.step_count, new.step_count), (count, count))
self.check_boundary(old, new, mesh, amplitudes) self.check_boundary(old, new, mesh, amplitudes)
@@ -133,7 +133,7 @@ class MultiplePinnedRemapGPU(unittest.TestCase):
previous = current previous = current
self.change_amplitudes(previous, amplitudes) self.change_amplitudes(previous, amplitudes)
before, count = previous.download(), previous.step_count before, count = previous.download(), previous.step_count
old_mesh = topology(previous, previous.child_patch.bounds, objects) old_mesh = topology(previous, previous.hierarchy.pairs[1].bounds, objects)
new_mesh = topology(previous, rectangles, objects) new_mesh = topology(previous, rectangles, objects)
expected = remap_synchronized( expected = remap_synchronized(
old_mesh, new_mesh, before, base_count=count, old_mesh, new_mesh, before, base_count=count,
@@ -141,7 +141,7 @@ class MultiplePinnedRemapGPU(unittest.TestCase):
phase='postcollision_g', phase='postcollision_g',
).g ).g
current = owned.enter_context( current = owned.enter_context(
previous.remap_child(rectangles, child_objects=objects)) previous.remap_finest(rectangles, objects=objects))
self.assertIsNot(previous, current) self.assertIsNot(previous, current)
self.assertEqual((previous.step_count, current.step_count), (count, count)) self.assertEqual((previous.step_count, current.step_count), (count, count))
self.check_boundary(previous, current, new_mesh, amplitudes) self.check_boundary(previous, current, new_mesh, amplitudes)
@@ -36,12 +36,12 @@ class PinnedBodyRemapGPU(unittest.TestCase):
def check_boundary(self, before, after): def check_boundary(self, before, after):
for name in ('flags','boundary_cells','curved_hits'): for name in ('flags','boundary_cells','curved_hits'):
exact(self,getattr(after.child_boundary,name),getattr(before.child_boundary,name),name) exact(self,getattr(after.boundaries[2],name),getattr(before.boundaries[2],name),name)
for a,b in zip(after.child_boundary.curved_host,before.child_boundary.curved_host): for a,b in zip(after.boundaries[2].curved_host,before.boundaries[2].curved_host):
exact(self,a,b,'physical cut links') exact(self,a,b,'physical cut links')
exact(self,after._boundaries[2].tables['surface_velocity'], exact(self,after._level_boundaries[2].tables['surface_velocity'],
before._boundaries[2].tables['surface_velocity'],'current per-link velocities') before._level_boundaries[2].tables['surface_velocity'],'current per-link velocities')
self.assertEqual(after._boundaries[2].jets[0].amplitude,.017) self.assertEqual(after._level_boundaries[2].jets[0].amplitude,.017)
def test_noop_next_step_bitwise(self): def test_noop_next_step_bitwise(self):
with CudaContext(0) as ctx, ExitStack() as owned: with CudaContext(0) as ctx, ExitStack() as owned:
@@ -49,7 +49,7 @@ class PinnedBodyRemapGPU(unittest.TestCase):
old.initialize(*initial(old,False));old.step(2) old.initialize(*initial(old,False));old.step(2)
old.set_jet_amplitude(0,.017,level=2) old.set_jet_amplitude(0,.017,level=2)
before=old.download() before=old.download()
new=owned.enter_context(old.remap_child(A,child_objects=objects)) new=owned.enter_context(old.remap_finest(A,objects=objects))
self.check_boundary(old,new);self.check_reset(new,before) self.check_boundary(old,new);self.check_reset(new,before)
old.step();new.step() old.step();new.step()
for a,b in zip(new.download(),old.download()): exact(self,a,b,'no-op next step bitwise') for a,b in zip(new.download(),old.download()): exact(self,a,b,'no-op next step bitwise')
@@ -63,14 +63,14 @@ class PinnedBodyRemapGPU(unittest.TestCase):
current.set_jet_amplitude(0,.017,level=2) current.set_jet_amplitude(0,.017,level=2)
for rects in (B,A): for rects in (B,A):
previous=current;before=previous.download();count=previous.step_count previous=current;before=previous.download();count=previous.step_count
jets=previous._boundaries[2].jets jets=previous._level_boundaries[2].jets
old=RemapTopology.from_config(CFG,previous.patch,previous.child_patch.bounds, old=RemapTopology.from_config(CFG,previous.hierarchy.pairs[0],previous.hierarchy.pairs[1].bounds,
child_objects=objects,child_jets=jets) child_objects=objects,child_jets=jets)
new=RemapTopology.from_config(CFG,previous.patch,rects, new=RemapTopology.from_config(CFG,previous.hierarchy.pairs[0],rects,
child_objects=objects,child_jets=jets) child_objects=objects,child_jets=jets)
expected=remap_synchronized(old,new,before,base_count=count, expected=remap_synchronized(old,new,before,base_count=count,
level_steps=(count,2*count,4*count),phase='postcollision_g').g level_steps=(count,2*count,4*count),phase='postcollision_g').g
current=owned.enter_context(previous.remap_child(rects,child_objects=objects)) current=owned.enter_context(previous.remap_finest(rects,objects=objects))
self.check_boundary(previous,current);self.check_reset(current,expected) self.check_boundary(previous,current);self.check_reset(current,expected)
for a,b,om,nm in zip(current.download(),before,old.leaf_masks,new.leaf_masks): for a,b,om,nm in zip(current.download(),before,old.leaf_masks,new.leaf_masks):
exact(self,a[:,om & nm],b[:,om & nm],'physical-fluid overlap') exact(self,a[:,om & nm],b[:,om & nm],'physical-fluid overlap')
@@ -61,7 +61,7 @@ def census(run, index, buffer):
layout, b = level.layout, level.boundary layout, b = level.layout, level.boundary
g = download(buffer, layout.blocked_shape) g = download(buffer, layout.blocked_shape)
dense = layout.to_dense(g) dense = layout.to_dense(g)
leaf = run._leaf[index] leaf = run.leaf[index]
ny, nx = leaf.shape ny, nx = leaf.shape
solid = (b.host.flags.reshape(ny, nx) & SOLID) != 0 solid = (b.host.flags.reshape(ny, nx) & SOLID) != 0
y, x = np.nonzero(leaf) y, x = np.nonzero(leaf)
@@ -119,7 +119,7 @@ class RawExchangeGPU(unittest.TestCase):
measured._stream.synchronize() measured._stream.synchronize()
after = measured._levels[1].layout.to_dense( after = measured._levels[1].layout.to_dense(
download(dst, measured._levels[1].layout.blocked_shape)) download(dst, measured._levels[1].layout.blocked_shape))
change = totals(after, measured._leaf[1])-before change = totals(after, measured.leaf[1])-before
expected = external + np.array([-raw[:, 0].astype(float).sum(), expected = external + np.array([-raw[:, 0].astype(float).sum(),
-raw[:, 1].astype(float).sum(), -raw[:, 1].astype(float).sum(),
raw[:, 3].astype(float).sum()]) raw[:, 3].astype(float).sum()])
@@ -127,8 +127,8 @@ class RawExchangeGPU(unittest.TestCase):
samples.append(raw.copy()) samples.append(raw.copy())
measured._levels[1].kernel = audited_step measured._levels[1].kernel = audited_step
for speed in (0., .002, -.001, 0.): for speed in (0., .002, -.001, 0.):
measured.set_jet_amplitude(0, speed) measured.set_jet_amplitude(0, speed, level=1)
plain.set_jet_amplitude(0, speed) plain.set_jet_amplitude(0, speed, level=1)
samples.clear() samples.clear()
measured.step() measured.step()
plain.step() plain.step()
@@ -139,11 +139,11 @@ class RawExchangeGPU(unittest.TestCase):
np.testing.assert_array_equal(integral, expected) np.testing.assert_array_equal(integral, expected)
last = download(b.last_raw, (count, 4)) last = download(b.last_raw, (count, 4))
np.testing.assert_array_equal(last, samples[-1]) np.testing.assert_array_equal(last, samples[-1])
obs = measured.read_raw_exchange(0) obs = measured.read_raw_exchange(0, level=1)
assert obs['sample_count'] == 2 assert obs['sample_count'] == 2
assert b.tables['aperture'].any() assert b.tables['aperture'].any()
measured.step(0) measured.step(0)
np.testing.assert_array_equal(measured.read_raw_exchange(0)['total']['impulse'], np.testing.assert_array_equal(measured.read_raw_exchange(0, level=1)['total']['impulse'],
obs['total']['impulse']) obs['total']['impulse'])
def test_depth_four_interval_storage_and_sampling(self): def test_depth_four_interval_storage_and_sampling(self):
@@ -171,7 +171,7 @@ class RawExchangeGPU(unittest.TestCase):
run._stream.synchronize() run._stream.synchronize()
after = run._levels[3].layout.to_dense( after = run._levels[3].layout.to_dense(
download(dst, run._levels[3].layout.blocked_shape)) download(dst, run._levels[3].layout.blocked_shape))
change = totals(after, run._leaf[3])-before change = totals(after, run.leaf[3])-before
expected = external + np.array([-independent[:, 0].sum(), expected = external + np.array([-independent[:, 0].sum(),
-independent[:, 1].sum(), independent[:, 3].sum()]) -independent[:, 1].sum(), independent[:, 3].sum()])
np.testing.assert_allclose(change, expected, rtol=2e-4, atol=8e-5) np.testing.assert_allclose(change, expected, rtol=2e-4, atol=8e-5)
@@ -234,7 +234,7 @@ class RawExchangeGPU(unittest.TestCase):
samples.append(download(raw, (count, 4))) samples.append(download(raw, (count, 4)))
root.kernel = measured_boundary root.kernel = measured_boundary
for speed in (0., .002, -.001, 0.): for speed in (0., .002, -.001, 0.):
amr.set_jet_amplitude(0, speed) amr.set_jet_amplitude(0, speed, level=1)
cuda.memcpy_htod(velocity, b.tables['surface_velocity']) cuda.memcpy_htod(velocity, b.tables['surface_velocity'])
samples.clear() samples.clear()
amr.step(); uniform.step(2) amr.step(); uniform.step(2)
@@ -204,7 +204,7 @@ class BlockedBaseCPU(unittest.TestCase):
level_blocks=[(3, 4), FINE_BLOCK], level_blocks=[(3, 4), FINE_BLOCK],
level_objects=[(), (circle(actuation(.002)),)]) level_objects=[(), (circle(actuation(.002)),)])
try: try:
backend.initialize(state(NY, NX, 1), state(*backend.patch.fine_shape, 2)) backend.initialize(state(NY, NX, 1), state(*backend.hierarchy.pairs[0].fine_shape, 2))
coarse, fine = backend._levels coarse, fine = backend._levels
current, spare, fine_current = coarse.current, coarse.spare, fine.current current, spare, fine_current = coarse.current, coarse.spare, fine.current
backend.step(2) backend.step(2)
@@ -255,7 +255,7 @@ class BlockedBaseCPU(unittest.TestCase):
with mocked_backend() as m: with mocked_backend() as m:
backend = FixedPatchStepper(config("periodic"), context=m["context"], refinements=[BOUNDS], backend = FixedPatchStepper(config("periodic"), context=m["context"], refinements=[BOUNDS],
level_blocks=[(3, 4), FINE_BLOCK]) level_blocks=[(3, 4), FINE_BLOCK])
backend.initialize(state(NY, NX, 1), state(*backend.patch.fine_shape, 2)) backend.initialize(state(NY, NX, 1), state(*backend.hierarchy.pairs[0].fine_shape, 2))
m["functions"]["BlockCoalesce"].side_effect = RuntimeError("injected parent stage failure") m["functions"]["BlockCoalesce"].side_effect = RuntimeError("injected parent stage failure")
current = backend._levels[0].current current = backend._levels[0].current
with self.assertRaisesRegex(RuntimeError, "injected parent stage"): with self.assertRaisesRegex(RuntimeError, "injected parent stage"):
@@ -125,7 +125,7 @@ class BlockedFineCPU(unittest.TestCase):
patch.object(cuda, "mem_alloc", side_effect=alloc), patch.object(cuda, "memcpy_htod"): patch.object(cuda, "mem_alloc", side_effect=alloc), patch.object(cuda, "memcpy_htod"):
backend = FixedPatchStepper(config("periodic"), context=context, refinements=[BOUNDS], backend = FixedPatchStepper(config("periodic"), context=context, refinements=[BOUNDS],
level_blocks=[(12, 8), FINE_BLOCKS[0]]) level_blocks=[(12, 8), FINE_BLOCKS[0]])
self.assertEqual(len(backend.fine_layout.halo_dst), 0) self.assertEqual(len(backend.layouts[1].halo_dst), 0)
self.assertIsNone(backend._levels[1].halo_dst) self.assertIsNone(backend._levels[1].halo_dst)
backend._gather_halo(backend._levels[1], Mock()) # no launch without pairs backend._gather_halo(backend._levels[1], Mock()) # no launch without pairs
module.get_function.return_value.assert_not_called() module.get_function.return_value.assert_not_called()
@@ -178,7 +178,7 @@ class BlockedFineCPU(unittest.TestCase):
backend = FixedPatchStepper(config("periodic"), context=context, refinements=[BOUNDS], backend = FixedPatchStepper(config("periodic"), context=context, refinements=[BOUNDS],
level_blocks=[(12, 8), FINE_BLOCKS[1]], level_blocks=[(12, 8), FINE_BLOCKS[1]],
level_objects=[(), tuple(objects)]) level_objects=[(), tuple(objects)])
backend.initialize(state(NY, NX, 1), state(*backend.patch.fine_shape, 2)) backend.initialize(state(NY, NX, 1), state(*backend.hierarchy.pairs[0].fine_shape, 2))
return backend, calls return backend, calls
def test_same_buffer_halo_bfl_step_launch_order(self): def test_same_buffer_halo_bfl_step_launch_order(self):
@@ -212,7 +212,7 @@ class BlockedFineCPU(unittest.TestCase):
def test_actuation_without_velocity_buffer_fails_before_bfl(self): def test_actuation_without_velocity_buffer_fails_before_bfl(self):
for speed in (0.0, .002): # any spec, even zero amplitude for speed in (0.0, .002): # any spec, even zero amplitude
backend, calls = self.mocked_backend([circle(actuation(speed))]) backend, calls = self.mocked_backend([circle(actuation(speed))])
self.assertIsNotNone(backend._fine_tables["surface_velocity"]) self.assertIsNotNone(backend._level_boundaries[1].tables["surface_velocity"])
del backend._level_boundaries[1].device["surface_velocity"] del backend._level_boundaries[1].device["surface_velocity"]
with self.assertRaisesRegex(RuntimeError, "Surface velocity device buffer is missing"): with self.assertRaisesRegex(RuntimeError, "Surface velocity device buffer is missing"):
backend.step(1) backend.step(1)
@@ -20,7 +20,7 @@ from CelerisLab.lbm import refinement_stepper as implementation
from CelerisLab.lbm.blocks import BlockLayout2D from CelerisLab.lbm.blocks import BlockLayout2D
from test_refinement_gpu import (ROOT, POP_ATOL, MOMENT_ATOL, ref, postcollision, from test_refinement_gpu import (ROOT, POP_ATOL, MOMENT_ATOL, ref, postcollision,
totals, LBMConfig, CudaContext, compiler, totals, LBMConfig, CudaContext, compiler,
PeriodicPatchStepper) FixedPatchStepper)
def scalar_interval(coarse, fine, bounds, slot=None): def scalar_interval(coarse, fine, bounds, slot=None):
@@ -123,7 +123,7 @@ class FaultGuardsCPU(unittest.TestCase):
refinements=[bounds], level_blocks=[(1, 1)] * 2) refinements=[bounds], level_blocks=[(1, 1)] * 2)
def fake_backend(self): def fake_backend(self):
backend = PeriodicPatchStepper.__new__(PeriodicPatchStepper) backend = FixedPatchStepper.__new__(FixedPatchStepper)
layout = BlockLayout2D.build(1, 1, 1, 1, periodic=(False, False)) layout = BlockLayout2D.build(1, 1, 1, 1, periodic=(False, False))
backend._layouts = [layout, layout] # blocked storage on both levels backend._layouts = [layout, layout] # blocked storage on both levels
backend._ddf_shifting = False # unshifted on every level, as _check_storage expects backend._ddf_shifting = False # unshifted on every level, as _check_storage expects
@@ -175,7 +175,7 @@ class FaultGuardsCPU(unittest.TestCase):
context = Mock(activate=nullcontext, sm_arch="sm_70") context = Mock(activate=nullcontext, sm_arch="sm_70")
with patch.object(implementation, "HALO", 2): with patch.object(implementation, "HALO", 2):
with self.assertRaisesRegex(RuntimeError, "four-cell"): with self.assertRaisesRegex(RuntimeError, "four-cell"):
PeriodicPatchStepper(cfg, context=context, refinements=[(4, 4, 12, 12)], FixedPatchStepper(cfg, context=context, refinements=[(4, 4, 12, 12)],
level_blocks=[(16, 16), (24, 24)]) level_blocks=[(16, 16), (24, 24)])
p = implementation.FixedPatch2D.from_config(cfg, (4, 4, 12, 12)) p = implementation.FixedPatch2D.from_config(cfg, (4, 4, 12, 12))
from CelerisLab.lbm import hierarchy from CelerisLab.lbm import hierarchy
@@ -187,11 +187,11 @@ class FaultGuardsCPU(unittest.TestCase):
with patch.object(cuda, "Stream"), patch.object(cuda, "mem_get_info", return_value=(1 << 40, 1 << 40)), bad, patch.object(cuda, "mem_alloc") as alloc, \ with patch.object(cuda, "Stream"), patch.object(cuda, "mem_get_info", return_value=(1 << 40, 1 << 40)), bad, patch.object(cuda, "mem_alloc") as alloc, \
patch.object(cuda, "memcpy_htod") as upload: patch.object(cuda, "memcpy_htod") as upload:
with self.assertRaisesRegex(RuntimeError, "config chain is inconsistent"): with self.assertRaisesRegex(RuntimeError, "config chain is inconsistent"):
PeriodicPatchStepper(cfg, context=context, refinements=[p.bounds], FixedPatchStepper(cfg, context=context, refinements=[p.bounds],
level_blocks=[(16, 16), (24, 24)]) level_blocks=[(16, 16), (24, 24)])
with patch.object(hierarchy, "_check_chain"), \ with patch.object(hierarchy, "_check_chain"), \
self.assertRaisesRegex(RuntimeError, "omega/tau"): self.assertRaisesRegex(RuntimeError, "omega/tau"):
PeriodicPatchStepper(cfg, context=context, refinements=[p.bounds], FixedPatchStepper(cfg, context=context, refinements=[p.bounds],
level_blocks=[(16, 16), (24, 24)]) level_blocks=[(16, 16), (24, 24)])
alloc.assert_not_called() alloc.assert_not_called()
upload.assert_not_called() upload.assert_not_called()
@@ -202,7 +202,7 @@ class FaultGuardsCPU(unittest.TestCase):
patch.object(compiler, "load_module", return_value=module), \ patch.object(compiler, "load_module", return_value=module), \
patch.object(cuda, "memcpy_htod") as upload: patch.object(cuda, "memcpy_htod") as upload:
with self.assertRaisesRegex(RuntimeError, "layout mismatch"): with self.assertRaisesRegex(RuntimeError, "layout mismatch"):
PeriodicPatchStepper(cfg, context=context, refinements=[p.bounds], FixedPatchStepper(cfg, context=context, refinements=[p.bounds],
level_blocks=[(16, 16), (24, 24)]) level_blocks=[(16, 16), (24, 24)])
upload.assert_not_called() upload.assert_not_called()
@@ -304,7 +304,7 @@ class NormalBuildGPU(unittest.TestCase):
if not (x0 <= sx < x1 and y0 <= sy < y1): if not (x0 <= sx < x1 and y0 <= sy < y1):
corners.add(target) corners.add(target)
self.assertEqual(len(corners), 8) self.assertEqual(len(corners), 8)
self.assertEqual(set(backend.patch.coalesce_coarse.tolist()), targets) self.assertEqual(set(backend.hierarchy.pairs[0].coalesce_coarse.tolist()), targets)
backend.initialize(*initial) backend.initialize(*initial)
before = discrete_mass(initial) before = discrete_mass(initial)
for step in range(1, 9): for step in range(1, 9):
@@ -348,7 +348,7 @@ class NormalBuildGPU(unittest.TestCase):
state = ref.State(*(a * (1 + 1e-3*rng.standard_normal(a.shape)) state = ref.State(*(a * (1 + 1e-3*rng.standard_normal(a.shape))
for a in (initial.coarse_f, initial.fine_f))) for a in (initial.coarse_f, initial.fine_f)))
max_pop = max_drift = max_moment = 0.0 max_pop = max_drift = max_moment = 0.0
with CudaContext(0) as context, PeriodicPatchStepper( with CudaContext(0) as context, FixedPatchStepper(
LBMConfig(nx=16, ny=16, viscosity=.1), context=context, LBMConfig(nx=16, ny=16, viscosity=.1), context=context,
refinements=[(4, 4, 12, 12)], level_blocks=[(16, 16), (24, 24)]) as backend: refinements=[(4, 4, 12, 12)], level_blocks=[(16, 16), (24, 24)]) as backend:
backend.initialize(*postcollision(state, spec, maps)) backend.initialize(*postcollision(state, spec, maps))
@@ -381,10 +381,10 @@ class NormalBuildGPU(unittest.TestCase):
fine_mask[4:-4, 4:-4] = True fine_mask[4:-4, 4:-4] = True
initial[1][:, ~fine_mask] = np.nan initial[1][:, ~fine_mask] = np.nan
expected = scalar_interval(*[a.astype(np.float64) for a in initial], bounds) expected = scalar_interval(*[a.astype(np.float64) for a in initial], bounds)
with CudaContext(0) as context, PeriodicPatchStepper(cfg, context=context, refinements=[bounds], with CudaContext(0) as context, FixedPatchStepper(cfg, context=context, refinements=[bounds],
level_blocks=[(13, 17), (14, 16)]) as backend: level_blocks=[(13, 17), (14, 16)]) as backend:
self.assertEqual(backend.patch.fine_shape, (14, 16)) self.assertEqual(backend.hierarchy.pairs[0].fine_shape, (14, 16))
p = backend.patch p = backend.hierarchy.pairs[0]
allocated = 2 * 4 * sum(int(np.prod(level.layout.blocked_shape)) for level in backend._levels) \ allocated = 2 * 4 * sum(int(np.prod(level.layout.blocked_shape)) for level in backend._levels) \
+ sum(a.nbytes for a in (p.ghost_fine, p.ghost_coarse, p.coalesce_coarse, p.coalesce_fine)) + sum(a.nbytes for a in (p.ghost_fine, p.ghost_coarse, p.coalesce_coarse, p.coalesce_fine))
self.assertLess(allocated, 128*1024) self.assertLess(allocated, 128*1024)
@@ -22,7 +22,7 @@ import CelerisLab
from CelerisLab.config import LBMConfig from CelerisLab.config import LBMConfig
from CelerisLab.cuda import compiler from CelerisLab.cuda import compiler
from CelerisLab.cuda.context import CudaContext from CelerisLab.cuda.context import CudaContext
from CelerisLab.lbm.refinement_stepper import FixedPatchStepper, PeriodicPatchStepper from CelerisLab.lbm.refinement_stepper import FixedPatchStepper
ROOT = Path(__file__).resolve().parents[3] ROOT = Path(__file__).resolve().parents[3]
POP_ATOL = 3e-6 POP_ATOL = 3e-6
@@ -91,8 +91,8 @@ class CoupledGPU(unittest.TestCase):
ddf_shifting=shifted) ddf_shifting=shifted)
with FixedPatchStepper(cfg, context=context, refinements=[(4, 4, 12, 12)], with FixedPatchStepper(cfg, context=context, refinements=[(4, 4, 12, 12)],
**storage) as backend: **storage) as backend:
self.assertAlmostEqual(backend.patch.tau_coarse, spec.tau_coarse) self.assertAlmostEqual(backend.hierarchy.pairs[0].tau_coarse, spec.tau_coarse)
self.assertAlmostEqual(backend.patch.tau_fine, spec.tau_fine) self.assertAlmostEqual(backend.hierarchy.pairs[0].tau_fine, spec.tau_fine)
self.assertEqual(tuple(level.shape for level in backend._levels), self.assertEqual(tuple(level.shape for level in backend._levels),
((9, 16, 16), (9, 24, 24))) ((9, 16, 16), (9, 24, 24)))
# Physical FP64 g is encoded by initialize before FP32 rounding. # Physical FP64 g is encoded by initialize before FP32 rounding.
@@ -101,7 +101,7 @@ class CoupledGPU(unittest.TestCase):
backend.step(1) backend.step(1)
actual = backend.download() actual = backend.download()
self.assertEqual(backend.step_count, 1) self.assertEqual(backend.step_count, 1)
backend.patch.check_synchronization( backend.hierarchy.pairs[0].check_synchronization(
1, next_state.tick, phase="postcollision_g") 1, next_state.tick, phase="postcollision_g")
for got, want, arrival in zip(actual, expected, arrivals): for got, want, arrival in zip(actual, expected, arrivals):
np.testing.assert_array_equal(np.isnan(got), np.isnan(want)) np.testing.assert_array_equal(np.isnan(got), np.isnan(want))
@@ -124,7 +124,7 @@ class CoupledGPU(unittest.TestCase):
"wrong-phase controls PASS", flush=True) "wrong-phase controls PASS", flush=True)
def test_cleanup_attempts_every_allocation(self): def test_cleanup_attempts_every_allocation(self):
backend = PeriodicPatchStepper.__new__(PeriodicPatchStepper) backend = FixedPatchStepper.__new__(FixedPatchStepper)
backend.ctx = Mock(activate=nullcontext) backend.ctx = Mock(activate=nullcontext)
backend._lock = RLock() backend._lock = RLock()
backend._status = "failed" backend._status = "failed"
@@ -146,7 +146,7 @@ class CoupledGPU(unittest.TestCase):
spec = ref.make_spec(collision="srt") spec = ref.make_spec(collision="srt")
maps = ref.build_maps(spec) maps = ref.build_maps(spec)
with CudaContext(0) as context: with CudaContext(0) as context:
with PeriodicPatchStepper(LBMConfig(nx=16, ny=16, viscosity=.1), with FixedPatchStepper(LBMConfig(nx=16, ny=16, viscosity=.1),
context=context, refinements=[(4, 4, 12, 12)], context=context, refinements=[(4, 4, 12, 12)],
level_blocks=[(16, 16), (24, 24)]) as backend: level_blocks=[(16, 16), (24, 24)]) as backend:
backend.initialize(*postcollision(ref.make_state(spec), spec, maps)) backend.initialize(*postcollision(ref.make_state(spec), spec, maps))
@@ -165,13 +165,13 @@ class CoupledGPU(unittest.TestCase):
with self.assertRaisesRegex(RuntimeError, "failed"): with self.assertRaisesRegex(RuntimeError, "failed"):
operation() operation()
raw = level_storage(fine) raw = level_storage(fine)
self.assertTrue(np.isfinite(raw.ravel()[backend.patch.coalesce_fine]).all()) self.assertTrue(np.isfinite(raw.ravel()[backend.hierarchy.pairs[0].coalesce_fine]).all())
self.assertFalse(context.is_closed) self.assertFalse(context.is_closed)
print("partial failure: two fine hops completed, no coarse step published; cleanup PASS", flush=True) print("partial failure: two fine hops completed, no coarse step published; cleanup PASS", flush=True)
def eight_steps(self, backend, spec, maps): def eight_steps(self, backend, spec, maps):
"""Rest and seeded nonequilibrium, 8 steps against the scalar reference.""" """Rest and seeded nonequilibrium, 8 steps against the scalar reference."""
p = backend.patch p = backend.hierarchy.pairs[0]
for case in ("rest", "nonequilibrium_seed1"): for case in ("rest", "nonequilibrium_seed1"):
with self.subTest(case=case): with self.subTest(case=case):
state = ref.make_state(spec, velocity=(0, 0) if case == "rest" else (.03, -.02)) state = ref.make_state(spec, velocity=(0, 0) if case == "rest" else (.03, -.02))
@@ -223,7 +223,7 @@ class CoupledGPU(unittest.TestCase):
maps = ref.build_maps(spec) maps = ref.build_maps(spec)
cfg = LBMConfig(nx=16, ny=16, viscosity=.1) cfg = LBMConfig(nx=16, ny=16, viscosity=.1)
with CudaContext(0) as context: with CudaContext(0) as context:
with PeriodicPatchStepper(cfg, context=context, refinements=[(4, 4, 12, 12)], with FixedPatchStepper(cfg, context=context, refinements=[(4, 4, 12, 12)],
level_blocks=[(16, 16), (24, 24)]) as backend: level_blocks=[(16, 16), (24, 24)]) as backend:
self.eight_steps(backend, spec, maps) self.eight_steps(backend, spec, maps)
with patch.object(backend._pairs[0], "explode", side_effect=RuntimeError("injected launch failure")): with patch.object(backend._pairs[0], "explode", side_effect=RuntimeError("injected launch failure")):
@@ -62,8 +62,8 @@ class JetStepperCPU(unittest.TestCase):
C = (16.0, 12.0) C = (16.0, 12.0)
whole = Jet(Sector(C, .1, math.pi, 3.), Radial(C, inward=False), amplitude=0.) whole = Jet(Sector(C, .1, math.pi, 3.), Radial(C, inward=False), amplitude=0.)
backend, _ = mocked([circle()], [whole], FINE_BLOCKS[1]) backend, _ = mocked([circle()], [whole], FINE_BLOCKS[1])
_, jet_id = backend._jet_cache _, jet_id = backend._level_boundaries[1].jet_cache
self.assertGreater(len(np.unique(backend._fine_tables["block"][jet_id == 0])), 1) self.assertGreater(len(np.unique(backend._level_boundaries[1].tables["block"][jet_id == 0])), 1)
backend.close() backend.close()
def test_fine_ownership_uses_remapped_physical_owner(self): def test_fine_ownership_uses_remapped_physical_owner(self):
@@ -76,10 +76,10 @@ class JetStepperCPU(unittest.TestCase):
source = SimObject(99, geometry=CircleGeometry(*C, 2.15), center=C, radius=2.15) source = SimObject(99, geometry=CircleGeometry(*C, 2.15), center=C, radius=2.15)
whole = Jet(Sector(C, .1, math.pi, 3.), Radial(C, inward=False), amplitude=0.) whole = Jet(Sector(C, .1, math.pi, 3.), Radial(C, inward=False), amplitude=0.)
backend, _ = mocked([source], [whole], FINE_BLOCKS[1]) backend, _ = mocked([source], [whole], FINE_BLOCKS[1])
_, jet_id = backend._jet_cache _, jet_id = backend._level_boundaries[1].jet_cache
owners = backend.fine_boundary.curved_host[3][jet_id == 0] owners = backend.boundaries[1].curved_host[3][jet_id == 0]
self.assertEqual(np.unique(owners).tolist(), [0]) self.assertEqual(np.unique(owners).tolist(), [0])
self.assertGreater(len(np.unique(backend._fine_tables["block"][jet_id == 0])), 1) self.assertGreater(len(np.unique(backend._level_boundaries[1].tables["block"][jet_id == 0])), 1)
self.assertEqual(source.obj_id, 99) self.assertEqual(source.obj_id, 99)
backend.close() backend.close()
cfg = config("periodic") cfg = config("periodic")
@@ -97,8 +97,8 @@ class JetStepperCPU(unittest.TestCase):
def test_blocked_jet_velocity_matches_legacy_rows(self): def test_blocked_jet_velocity_matches_legacy_rows(self):
legacy, _ = mocked([circle(actuation(SPEED))], (), FINE_BLOCKS[1]) legacy, _ = mocked([circle(actuation(SPEED))], (), FINE_BLOCKS[1])
driven, _ = mocked([circle()], [jet()], FINE_BLOCKS[1]) driven, _ = mocked([circle()], [jet()], FINE_BLOCKS[1])
a = legacy._fine_tables["surface_velocity"] a = legacy._level_boundaries[1].tables["surface_velocity"]
b = driven._fine_tables["surface_velocity"] b = driven._level_boundaries[1].tables["surface_velocity"]
diff = np.abs(a - b).max() diff = np.abs(a - b).max()
print(f"fine jet vs legacy velocity rows: max |dU| = {diff:.3e} (exact={diff == 0})") print(f"fine jet vs legacy velocity rows: max |dU| = {diff:.3e} (exact={diff == 0})")
self.assertLessEqual(diff, 4 * np.finfo(np.float32).eps * SPEED) self.assertLessEqual(diff, 4 * np.finfo(np.float32).eps * SPEED)
@@ -106,20 +106,20 @@ class JetStepperCPU(unittest.TestCase):
def test_blocked_amplitude_update_rewrites_only_velocity_column(self): def test_blocked_amplitude_update_rewrites_only_velocity_column(self):
backend, htod = mocked([circle()], [jet()], FINE_BLOCKS[1]) backend, htod = mocked([circle()], [jet()], FINE_BLOCKS[1])
before = backend._fine_tables["surface_velocity"].copy() before = backend._level_boundaries[1].tables["surface_velocity"].copy()
device = backend._level_boundaries[1].device["surface_velocity"] device = backend._level_boundaries[1].device["surface_velocity"]
with patch.object(cuda, "memcpy_htod", htod): with patch.object(cuda, "memcpy_htod", htod):
htod.reset_mock() htod.reset_mock()
backend.set_jet_amplitude(0, -SPEED) backend.set_jet_amplitude(0, -SPEED, level=1)
(ptr, sent), _ = htod.call_args (ptr, sent), _ = htod.call_args
self.assertEqual(htod.call_count, 1) self.assertEqual(htod.call_count, 1)
self.assertIs(ptr, device) self.assertIs(ptr, device)
np.testing.assert_allclose(sent, -before, rtol=1e-6, atol=0) np.testing.assert_allclose(sent, -before, rtol=1e-6, atol=0)
np.testing.assert_array_equal(backend._fine_tables["surface_velocity"], sent) np.testing.assert_array_equal(backend._level_boundaries[1].tables["surface_velocity"], sent)
with self.assertRaises(ValueError): with self.assertRaises(ValueError):
backend.set_jet_amplitude(0, 1.5) backend.set_jet_amplitude(0, 1.5, level=1)
self.assertEqual(htod.call_count, 1) self.assertEqual(htod.call_count, 1)
self.assertEqual(backend._jets[0].amplitude, -SPEED) self.assertEqual(backend.jets[1][0].amplitude, -SPEED)
backend.close() backend.close()
def test_copy_failure_marks_failed_and_keeps_host_amplitude(self): def test_copy_failure_marks_failed_and_keeps_host_amplitude(self):
@@ -128,24 +128,24 @@ class JetStepperCPU(unittest.TestCase):
status = backend._status status = backend._status
for bad, error in ((1.5, ValueError), (float("nan"), ValueError)): for bad, error in ((1.5, ValueError), (float("nan"), ValueError)):
with self.assertRaises(error): with self.assertRaises(error):
backend.set_jet_amplitude(0, bad) backend.set_jet_amplitude(0, bad, level=1)
self.assertEqual(backend._status, status) # invalid input stays usable self.assertEqual(backend._status, status) # invalid input stays usable
with self.assertRaises(IndexError): with self.assertRaises(IndexError):
backend.set_jet_amplitude(1, -SPEED) backend.set_jet_amplitude(1, -SPEED, level=1)
self.assertEqual(backend._status, status) self.assertEqual(backend._status, status)
before = backend._fine_tables["surface_velocity"].copy() before = backend._level_boundaries[1].tables["surface_velocity"].copy()
boom = Mock(side_effect=RuntimeError("copy failed")) boom = Mock(side_effect=RuntimeError("copy failed"))
with patch.object(cuda, "memcpy_htod", boom), patch.object(cuda, "memcpy_htod_async", boom), \ with patch.object(cuda, "memcpy_htod", boom), patch.object(cuda, "memcpy_htod_async", boom), \
patch.object(cuda, "pagelocked_empty", side_effect=lambda shape, dtype: np.empty(shape, dtype)): patch.object(cuda, "pagelocked_empty", side_effect=lambda shape, dtype: np.empty(shape, dtype)):
with self.assertRaisesRegex(RuntimeError, "copy failed"): with self.assertRaisesRegex(RuntimeError, "copy failed"):
backend.set_jet_amplitude(0, -SPEED) backend.set_jet_amplitude(0, -SPEED, level=1)
self.assertEqual(backend._status, "failed", fine_block) self.assertEqual(backend._status, "failed", fine_block)
self.assertEqual(backend._jets[0].amplitude, SPEED) self.assertEqual(backend.jets[1][0].amplitude, SPEED)
np.testing.assert_array_equal(backend._fine_tables["surface_velocity"], before) np.testing.assert_array_equal(backend._level_boundaries[1].tables["surface_velocity"], before)
with self.assertRaisesRegex(RuntimeError, "failed"): with self.assertRaisesRegex(RuntimeError, "failed"):
backend.step(1) backend.step(1)
with self.assertRaisesRegex(RuntimeError, "failed"): with self.assertRaisesRegex(RuntimeError, "failed"):
backend.set_jet_amplitude(0, -SPEED) backend.set_jet_amplitude(0, -SPEED, level=1)
backend.close() backend.close()
def test_jet_overlapping_legacy_aperture_rejected_before_cuda(self): def test_jet_overlapping_legacy_aperture_rejected_before_cuda(self):
@@ -170,7 +170,7 @@ def run(mode, fine_block, *, objects, jets, switch=None):
backend.step(STEPS) backend.step(STEPS)
else: else:
backend.step(STEPS // 2) backend.step(STEPS // 2)
backend.set_jet_amplitude(0, switch) backend.set_jet_amplitude(0, switch, level=1)
backend.step(STEPS - STEPS // 2) backend.step(STEPS - STEPS // 2)
return backend.download(), raw(backend) return backend.download(), raw(backend)
@@ -36,14 +36,14 @@ def level_fine(patch_set, windows):
"""Level-wide fine array with each patch window filled from `windows`.""" """Level-wide fine array with each patch window filled from `windows`."""
out = np.full((9, *patch_set.fine_shape), np.nan, np.float32) out = np.full((9, *patch_set.fine_shape), np.nan, np.float32)
for p, w in zip(patch_set.patches, windows): for p, w in zip(patch_set.patches, windows):
y0, x0 = p.fine_level_origin y0, x0 = p.fine_level_offset
out[:, y0:y0 + p.fine_shape[0], x0:x0 + p.fine_shape[1]] = w out[:, y0:y0 + p.fine_shape[0], x0:x0 + p.fine_shape[1]] = w
return out return out
def window(patch_set, index, fine): def window(patch_set, index, fine):
p = patch_set.patches[index] p = patch_set.patches[index]
y0, x0 = p.fine_level_origin y0, x0 = p.fine_level_offset
return fine[:, y0:y0 + p.fine_shape[0], x0:x0 + p.fine_shape[1]] return fine[:, y0:y0 + p.fine_shape[0], x0:x0 + p.fine_shape[1]]
@@ -175,7 +175,7 @@ class MultipatchBodiesGPU(unittest.TestCase):
for step in self.STEPS: for step in self.STEPS:
if switch is not None and step == 8: if switch is not None and step == 8:
for index, amplitude in switch: for index, amplitude in switch:
backend.set_jet_amplitude(index, amplitude) backend.set_jet_amplitude(index, amplitude, level=1)
backend.step(step - backend.step_count) backend.step(step - backend.step_count)
out.append(backend.download() + raw(backend)) out.append(backend.download() + raw(backend))
return out return out
@@ -128,8 +128,10 @@ class ThreeLevelCPU(unittest.TestCase):
run = mocked_child() run = mocked_child()
b = run._levels[2].boundary b = run._levels[2].boundary
before = b.tables["surface_velocity"].copy() before = b.tables["surface_velocity"].copy()
with self.assertRaises(TypeError):
run.set_jet_amplitude(0, -.002) # level is required
with self.assertRaises(IndexError): with self.assertRaises(IndexError):
run.set_jet_amplitude(0, -.002) # default stays L1 run.set_jet_amplitude(0, -.002, level=1) # level 1 has no Jets
for level in (0, 3, True): for level in (0, 3, True):
with self.assertRaises(ValueError): with self.assertRaises(ValueError):
run.set_jet_amplitude(0, -.002, level=level) run.set_jet_amplitude(0, -.002, level=level)
@@ -138,7 +140,7 @@ class ThreeLevelCPU(unittest.TestCase):
self.assertEqual(upload.call_count, 1) self.assertEqual(upload.call_count, 1)
self.assertIs(upload.call_args.args[0], b.device["surface_velocity"]) self.assertIs(upload.call_args.args[0], b.device["surface_velocity"])
np.testing.assert_array_equal(b.tables["surface_velocity"], -before) np.testing.assert_array_equal(b.tables["surface_velocity"], -before)
self.assertEqual(run._jets, ()) self.assertEqual(run.jets[1], ())
with patch.object(cuda, "memcpy_htod", side_effect=RuntimeError("copy failed")): with patch.object(cuda, "memcpy_htod", side_effect=RuntimeError("copy failed")):
with self.assertRaisesRegex(RuntimeError, "copy failed"): with self.assertRaisesRegex(RuntimeError, "copy failed"):
run.set_jet_amplitude(0, .002, level=2) run.set_jet_amplitude(0, .002, level=2)
@@ -150,17 +152,17 @@ class ThreeLevelCPU(unittest.TestCase):
def test_child_public_masks_include_solids(self): def test_child_public_masks_include_solids(self):
from CelerisLab.lbm.descriptors import SOLID from CelerisLab.lbm.descriptors import SOLID
run = mocked_child() run = mocked_child()
solid = run.child_boundary.flags.reshape(run.child_patch.fine_shape) & SOLID != 0 solid = run.boundaries[2].flags.reshape(run.hierarchy.pairs[1].fine_shape) & SOLID != 0
self.assertTrue(solid.any()) self.assertTrue(solid.any())
self.assertFalse(run._child_leaf[solid].any()) self.assertFalse(run.leaf[2][solid].any())
arrays = [np.ones(level.shape, np.float32) for level in run._levels] arrays = [np.ones(level.shape, np.float32) for level in run._levels]
run._mask_nonleaves(*arrays) run._mask_nonleaves(*arrays)
self.assertTrue(np.isnan(arrays[2][:, solid]).all()) self.assertTrue(np.isnan(arrays[2][:, solid]).all())
with patch.object(cuda, "memcpy_htod") as upload: with patch.object(cuda, "memcpy_htod") as upload:
run.initialize(*arrays) run.initialize(*arrays)
sent = upload.call_args_list[-1].args[1] sent = upload.call_args_list[-1].args[1]
dense = run.child_layout.to_dense(sent) dense = run.layouts[2].to_dense(sent)
self.assertTrue(np.isnan(dense[:, ~run._child_leaf]).all()) self.assertTrue(np.isnan(dense[:, ~run.leaf[2]]).all())
run.close() run.close()
@@ -169,7 +171,7 @@ class ThreeLevelGPU(unittest.TestCase):
from CelerisLab.cuda.context import CudaContext from CelerisLab.cuda.context import CudaContext
feq = np.asarray(compute_feq_d2q9(1.0, .03, -.02), np.float32)[:, None, None] feq = np.asarray(compute_feq_d2q9(1.0, .03, -.02), np.float32)[:, None, None]
with CudaContext(0) as ctx, FixedPatchStepper(CFG, context=ctx, **KW) as run: with CudaContext(0) as ctx, FixedPatchStepper(CFG, context=ctx, **KW) as run:
shapes = ((16, 16), run.patch.fine_shape, run.child_patch.fine_shape) shapes = ((16, 16), run.hierarchy.pairs[0].fine_shape, run.hierarchy.pairs[1].fine_shape)
run.initialize(*(np.broadcast_to(feq, (9, *s)).copy() for s in shapes)) run.initialize(*(np.broadcast_to(feq, (9, *s)).copy() for s in shapes))
run.step(8) run.step(8)
self.assertEqual(run.step_count, 8) self.assertEqual(run.step_count, 8)
@@ -182,7 +184,7 @@ class ThreeLevelGPU(unittest.TestCase):
def test_nonequilibrium_masks_and_fresh_middle_feeds(self): def test_nonequilibrium_masks_and_fresh_middle_feeds(self):
from CelerisLab.cuda.context import CudaContext from CelerisLab.cuda.context import CudaContext
with CudaContext(0) as ctx, FixedPatchStepper(CFG, context=ctx, **KW) as run: with CudaContext(0) as ctx, FixedPatchStepper(CFG, context=ctx, **KW) as run:
p1, p2 = run.patch, run.child_patch p1, p2 = run.hierarchy.pairs[0], run.hierarchy.pairs[1]
run.initialize(state(16, 16, 1), state(*p1.fine_shape, 2), state(*p2.fine_shape, 3)) run.initialize(state(16, 16, 1), state(*p1.fine_shape, 2), state(*p2.fine_shape, 3))
f = run._levels[1] f = run._levels[1]
feeds, original = [], run._advance feeds, original = [], run._advance
@@ -195,7 +197,7 @@ class ThreeLevelGPU(unittest.TestCase):
run.step(1) run.step(1)
self.assertEqual(feeds, [cur, spare]) # current G, then the new G of hop 1 self.assertEqual(feeds, [cur, spare]) # current G, then the new G of hop 1
run.step(7) run.step(7)
middle = run._middle_kind middle = run.hierarchy.kinds[1]
for g, leaf in zip(run.download(), (~p1.covered, middle == ACTIVE, p2.fine_kind == ACTIVE)): for g, leaf in zip(run.download(), (~p1.covered, middle == ACTIVE, p2.fine_kind == ACTIVE)):
self.assertTrue(np.isfinite(g[:, leaf]).all()) self.assertTrue(np.isfinite(g[:, leaf]).all())
self.assertTrue(np.isnan(g[:, ~leaf]).all()) self.assertTrue(np.isnan(g[:, ~leaf]).all())
@@ -227,15 +229,15 @@ class ThreeLevelBodyGPU(unittest.TestCase):
if stop == 4 and switch: if stop == 4 and switch:
run.set_jet_amplitude(0, -.002, level=2) run.set_jet_amplitude(0, -.002, level=2)
public, stored = run.download(), raw(run) public, stored = run.download(), raw(run)
leaves = [~run.patch.covered, run._fine_leaf, run._child_leaf] leaves = [~run.hierarchy.pairs[0].covered, run.leaf[1], run.leaf[2]]
if mode != "periodic": if mode != "periodic":
leaves[0][[0, -1]] = False leaves[0][[0, -1]] = False
for g, mask in zip(public, leaves): for g, mask in zip(public, leaves):
self.assertTrue(np.isfinite(g[:, mask]).all()) self.assertTrue(np.isfinite(g[:, mask]).all())
self.assertTrue(np.isnan(g[:, ~mask]).all()) self.assertTrue(np.isnan(g[:, ~mask]).all())
invalid = [run.patch.covered.copy(), run._middle_kind == OUTER_GHOST, invalid = [run.hierarchy.pairs[0].covered.copy(), run.hierarchy.kinds[1] == OUTER_GHOST,
(run.child_patch.fine_kind == OUTER_GHOST) | (run.hierarchy.pairs[1].fine_kind == OUTER_GHOST) |
((run.child_boundary.flags.reshape(run.child_patch.fine_shape) & SOLID) != 0)] ((run.boundaries[2].flags.reshape(run.hierarchy.pairs[1].fine_shape) & SOLID) != 0)]
if mode != "periodic": if mode != "periodic":
invalid[0][[0, -1]] = True invalid[0][[0, -1]] = True
for g, mask in zip(stored, invalid): for g, mask in zip(stored, invalid):
@@ -54,9 +54,9 @@ def jets(speed):
def frames(run): def frames(run):
middle_origin = np.asarray(run.patch.fine_origin) middle_origin = np.asarray(run.hierarchy.pairs[0].fine_origin)
return ((np.zeros(2), 1.), (middle_origin, .5), return ((np.zeros(2), 1.), (middle_origin, .5),
(middle_origin + .5*np.asarray(run.child_patch.fine_origin), .25)) (middle_origin + .5*np.asarray(run.hierarchy.pairs[1].fine_origin), .25))
def physical_state(shape, origin, spacing): def physical_state(shape, origin, spacing):
@@ -92,15 +92,15 @@ class ThreeLevelMultipatchGPU(unittest.TestCase):
for index in range(len(drive)): for index in range(len(drive)):
run.set_jet_amplitude(index, -.002, level=2) run.set_jet_amplitude(index, -.002, level=2)
public, stored = run.download(), raw(run) public, stored = run.download(), raw(run)
masks = [~run.patch.covered, run._fine_leaf, run._child_leaf] masks = [~run.hierarchy.pairs[0].covered, run.leaf[1], run.leaf[2]]
if mode != 'periodic': if mode != 'periodic':
masks[0][[0, -1]] = False masks[0][[0, -1]] = False
invalid = [run.patch.covered.copy(), run._middle_kind == OUTER_GHOST, invalid = [run.hierarchy.pairs[0].covered.copy(), run.hierarchy.kinds[1] == OUTER_GHOST,
run.child_patch.fine_kind == OUTER_GHOST] run.hierarchy.pairs[1].fine_kind == OUTER_GHOST]
if mode != 'periodic': if mode != 'periodic':
invalid[0][[0, -1]] = True invalid[0][[0, -1]] = True
if run.child_boundary is not None: if run.boundaries[2] is not None:
invalid[2] |= (run.child_boundary.flags.reshape(run.child_patch.fine_shape) & SOLID) != 0 invalid[2] |= (run.boundaries[2].flags.reshape(run.hierarchy.pairs[1].fine_shape) & SOLID) != 0
for level, (g, mask, dense, bad) in enumerate(zip(public, masks, stored, invalid)): for level, (g, mask, dense, bad) in enumerate(zip(public, masks, stored, invalid)):
label = f'{mode} stop={stop} level={level}' label = f'{mode} stop={stop} level={level}'
self.assertTrue(np.isfinite(g[:, mask]).all(), label) self.assertTrue(np.isfinite(g[:, mask]).all(), label)
@@ -39,7 +39,7 @@ def initial(run, uniform):
class RuntimeChildRemapGPU(unittest.TestCase): class RuntimeChildRemapGPU(unittest.TestCase):
def check_public(self, run): def check_public(self, run):
for array, leaf in zip(run.download(), (~run.patch.covered, run._fine_leaf, run._child_leaf)): for array, leaf in zip(run.download(), (~run.hierarchy.pairs[0].covered, run.leaf[1], run.leaf[2])):
self.assertTrue(np.isfinite(array[:, leaf]).all()) self.assertTrue(np.isfinite(array[:, leaf]).all())
self.assertTrue(np.isnan(array[:, ~leaf]).all()) self.assertTrue(np.isnan(array[:, ~leaf]).all())
@@ -57,7 +57,7 @@ class RuntimeChildRemapGPU(unittest.TestCase):
old.initialize(*initial(old, False)) old.initialize(*initial(old, False))
old.step(2) old.step(2)
before = old.download() before = old.download()
new = owned.enter_context(old.remap_child(A)) new = owned.enter_context(old.remap_finest(A))
self.assertEqual((old.step_count, new.step_count), (2, 2)) self.assertEqual((old.step_count, new.step_count), (2, 2))
self.check_reset(new, before) self.check_reset(new, before)
for got, want in zip(old.download(), before): exact(self, got, want, 'old unchanged') for got, want in zip(old.download(), before): exact(self, got, want, 'old unchanged')
@@ -75,12 +75,12 @@ class RuntimeChildRemapGPU(unittest.TestCase):
current.step(2) current.step(2)
for rectangles in (B, A): for rectangles in (B, A):
before, count = current.download(), current.step_count before, count = current.download(), current.step_count
old_mesh = RemapTopology.from_config(CFG, current.patch, current.child_patch.bounds) old_mesh = RemapTopology.from_config(CFG, current.hierarchy.pairs[0], current.hierarchy.pairs[1].bounds)
new_mesh = RemapTopology.from_config(CFG, current.patch, rectangles) new_mesh = RemapTopology.from_config(CFG, current.hierarchy.pairs[0], rectangles)
expected = remap_synchronized(old_mesh, new_mesh, before, base_count=count, expected = remap_synchronized(old_mesh, new_mesh, before, base_count=count,
level_steps=(count, 2*count, 4*count), phase='postcollision_g').g level_steps=(count, 2*count, 4*count), phase='postcollision_g').g
previous = current previous = current
current = owned.enter_context(previous.remap_child(rectangles)) current = owned.enter_context(previous.remap_finest(rectangles))
self.assertEqual(current.step_count, count) self.assertEqual(current.step_count, count)
self.check_reset(current, expected) self.check_reset(current, expected)
for got, want, old_g, old_leaf, new_leaf in zip( for got, want, old_g, old_leaf, new_leaf in zip(
@@ -8,14 +8,17 @@ from contextlib import ExitStack
import unittest import unittest
import numpy as np import numpy as np
from CelerisLab.cuda.context import CudaContext from CelerisLab.cuda.context import CudaContext
from research.shear_refinement.example import (CHILDREN, make, initialize, observe, from CelerisLab.lbm.refinement_stepper import FixedPatchStepper
topology, objects) from research.shear_refinement.example import (CFG, PARENTS, CHILDREN, initialize,
observe, topology, objects)
class ShearRefinementGPU(unittest.TestCase): class ShearRefinementGPU(unittest.TestCase):
def test_measured_growth_matches_explicit_manual_continuation(self): def test_measured_growth_matches_explicit_manual_continuation(self):
with CudaContext(0) as ctx, ExitStack() as owned: with CudaContext(0) as ctx, ExitStack() as owned:
old = owned.enter_context(make(ctx)) old = owned.enter_context(FixedPatchStepper(
CFG, context=ctx, refinements=[PARENTS, CHILDREN],
level_blocks=[(4, 4)] * 3, level_objects=[(), (), objects()]))
initialize(old) initialize(old)
old.step(2) old.step(2)
before = old.download() before = old.download()
@@ -25,10 +28,10 @@ class ShearRefinementGPU(unittest.TestCase):
self.assertGreater(report['requested'], 0) self.assertGreater(report['requested'], 0)
self.assertEqual(report['corridor_unknown'], 0) self.assertEqual(report['corridor_unknown'], 0)
self.assertEqual(report['corridor_unavailable'], 0) self.assertEqual(report['corridor_unavailable'], 0)
auto = owned.enter_context(old.remap_child(chosen, child_objects=objects())) auto = owned.enter_context(old.remap_finest(chosen, objects=objects()))
explicit = owned.enter_context(old.remap_child(manual, child_objects=objects())) explicit = owned.enter_context(old.remap_finest(manual, objects=objects()))
for run in (auto, explicit): for run in (auto, explicit):
self.assertEqual(run.child_patch.bounds[0], CHILDREN[0]) self.assertEqual(run.hierarchy.pairs[1].bounds[0], CHILDREN[0])
self.assertEqual(run.step_count, 2) self.assertEqual(run.step_count, 2)
run.step(2) run.step(2)
for a, b, leaf in zip(auto.download(), explicit.download(), topology(manual).leaf_masks): for a, b, leaf in zip(auto.download(), explicit.download(), topology(manual).leaf_masks):
@@ -58,7 +58,7 @@ def census(run,index,buffer):
"""Independent physical link/face census of actual staged receiver windows.""" """Independent physical link/face census of actual staged receiver windows."""
level=run._levels[index];layout=level.layout;b=level.boundary level=run._levels[index];layout=level.layout;b=level.boundary
stored=download(buffer,layout.blocked_shape);dense=layout.to_dense(stored) stored=download(buffer,layout.blocked_shape);dense=layout.to_dense(stored)
leaf=run._leaf[index];ny,nx=leaf.shape leaf=run.leaf[index];ny,nx=leaf.shape
solid=(b.host.flags.reshape(ny,nx)&SOLID)!=0 solid=(b.host.flags.reshape(ny,nx)&SOLID)!=0
y,x=np.nonzero(leaf);owner=layout.owner[y,x].astype(np.intp) 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) cells=int(np.prod(layout.local_shape));block,cell=np.divmod(owner,cells)
@@ -93,7 +93,7 @@ class Ledger:
self.max_row=max(self.max_row,float(abs(raw.astype(float)-independent).max())) self.max_row=max(self.max_row,float(abs(raw.astype(float)-independent).max()))
self.original(src,dst,*args,**kwargs);run._stream.synchronize() self.original(src,dst,*args,**kwargs);run._stream.synchronize()
after=level.layout.to_dense(download(dst,level.layout.blocked_shape)) after=level.layout.to_dense(download(dst,level.layout.blocked_shape))
change=raw_totals(after,run._leaf[index])-before change=raw_totals(after,run.leaf[index])-before
expected=external+np.array([-independent[:,0].sum(),-independent[:,1].sum(),independent[:,3].sum()]) 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) 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.max_balance=max(self.max_balance,float(abs(change-expected).max()))
@@ -113,7 +113,7 @@ class Ledger:
for name in ('impulse','angular_impulse','mass_into_fluid'): 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) 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) self.test.assertEqual(obs['sample_count'],2**self.index)
for a,leaf in zip(run._download_raw().values,run._leaf): for a,leaf in zip(run._download_raw().values,run.leaf):
self.test.assertTrue(np.isfinite(a[:,leaf]).all()) self.test.assertTrue(np.isfinite(a[:,leaf]).all())
self.test.assertTrue((1+a[:,leaf].sum(axis=0,dtype=np.float64)>0).all()) self.test.assertTrue((1+a[:,leaf].sum(axis=0,dtype=np.float64)>0).all())
return obs return obs
@@ -143,7 +143,7 @@ class ShiftedJetsGPU(unittest.TestCase):
self.assertTrue(run._level_boundaries[1].tables['aperture'].all()) self.assertTrue(run._level_boundaries[1].tables['aperture'].all())
for _ in range(4):ledger.step() for _ in range(4):ledger.step()
error=max(float(abs(a[:,leaf]-b[:,leaf]).max()) for a,b,leaf in error=max(float(abs(a[:,leaf]-b[:,leaf]).max()) for a,b,leaf in
zip(run._download_raw().values,expected.values,run._leaf)) zip(run._download_raw().values,expected.values,run.leaf))
self.assertLessEqual(error,1e-6) # inherited uniform-state cap self.assertLessEqual(error,1e-6) # inherited uniform-state cap
print(f'matched whole-circle rho={rho}: max_abs_raw_state={error:.9g}',flush=True) print(f'matched whole-circle rho={rho}: max_abs_raw_state={error:.9g}',flush=True)
ledger.report(f'matched whole-circle rho={rho}') ledger.report(f'matched whole-circle rho={rho}')
@@ -161,7 +161,7 @@ class ShiftedJetsGPU(unittest.TestCase):
ledger=Ledger(self,measured,1) ledger=Ledger(self,measured,1)
for speed in SPEEDS: for speed in SPEEDS:
before=measured._download_raw();count=measured.step_count before=measured._download_raw();count=measured.step_count
measured.set_jet_amplitude(0,speed);plain.set_jet_amplitude(0,speed) measured.set_jet_amplitude(0,speed, level=1);plain.set_jet_amplitude(0,speed, level=1)
self.exact_raw(before,measured._download_raw());self.assertEqual(measured.step_count,count) self.exact_raw(before,measured._download_raw());self.assertEqual(measured.step_count,count)
ledger.step();plain.step() ledger.step();plain.step()
self.exact_raw(measured._download_raw(),plain._download_raw()) self.exact_raw(measured._download_raw(),plain._download_raw())
@@ -193,14 +193,14 @@ class ShiftedJetsGPU(unittest.TestCase):
old._initialize_raw(initial(old,.99));ledger=Ledger(self,old,2) 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() for speed in (.002,-.001):old.set_jet_amplitude(0,speed,level=2);ledger.step()
before=old._download_raw();count=old.step_count before=old._download_raw();count=old.step_count
topology=lambda rects:RemapTopology.from_config(cfg,old.patch,rects,child_objects=(obj,), topology=lambda rects:RemapTopology.from_config(cfg,old.hierarchy.pairs[0],rects,child_objects=(obj,),
child_jets=old._level_jets[2],_shifted_blocked=True) child_jets=old.jets[2],_shifted_blocked=True)
expected=_remap_stored(topology(a),topology(b),before,shifted=True,base_count=count, expected=_remap_stored(topology(a),topology(b),before,shifted=True,base_count=count,
level_steps=(count,2*count,4*count),phase='postcollision_g') level_steps=(count,2*count,4*count),phase='postcollision_g')
new=owned.enter_context(old.remap_child(b,child_objects=(obj,))) new=owned.enter_context(old.remap_finest(b,objects=(obj,)))
self.exact_raw(before,old._download_raw()) self.exact_raw(before,old._download_raw())
self.exact_raw(_StoredPopulations(expected.g,True),new._download_raw()) self.exact_raw(_StoredPopulations(expected.g,True),new._download_raw())
self.assertEqual(new._level_jets[2][0].amplitude,-.001) self.assertEqual(new.jets[2][0].amplitude,-.001)
np.testing.assert_array_equal(new._level_boundaries[2].tables['surface_velocity'],old._level_boundaries[2].tables['surface_velocity']) 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() newledger=Ledger(self,new,2);ledger.step();newledger.step()
self.assertEqual((old.step_count,new.step_count),(count+1,count+1)) self.assertEqual((old.step_count,new.step_count),(count+1,count+1))
@@ -61,7 +61,7 @@ class ShiftedRemapGPU(unittest.TestCase):
def check_finite(self,run): def check_finite(self,run):
raw=run._download_raw() raw=run._download_raw()
self.assertTrue(raw.shifted) self.assertTrue(raw.shifted)
for a,leaf in zip(raw.values,run._leaf): for a,leaf in zip(raw.values,run.leaf):
self.assertTrue(np.isfinite(a[:,leaf]).all()) self.assertTrue(np.isfinite(a[:,leaf]).all())
self.assertTrue(np.isnan(a[:,~leaf]).all()) self.assertTrue(np.isnan(a[:,~leaf]).all())
self.assertTrue((1+a[:,leaf].sum(axis=0,dtype=np.float64)>0).all()) self.assertTrue((1+a[:,leaf].sum(axis=0,dtype=np.float64)>0).all())
@@ -71,7 +71,7 @@ class ShiftedRemapGPU(unittest.TestCase):
old=owned.enter_context(make(ctx)) old=owned.enter_context(make(ctx))
old._initialize_raw(initial(old));old.step(2) old._initialize_raw(initial(old));old.step(2)
before=old._download_raw() before=old._download_raw()
new=owned.enter_context(old.remap_child(A)) new=owned.enter_context(old.remap_finest(A))
self.check_reset(new,before.values) self.check_reset(new,before.values)
self.snapshot_equal(old._download_raw(),before,'old unchanged') self.snapshot_equal(old._download_raw(),before,'old unchanged')
self.snapshot_equal(new._download_raw(),before,'identity raw copy') self.snapshot_equal(new._download_raw(),before,'identity raw copy')
@@ -85,11 +85,11 @@ class ShiftedRemapGPU(unittest.TestCase):
old=owned.enter_context(make(ctx)) old=owned.enter_context(make(ctx))
old._initialize_raw(initial(old));old.step(2) old._initialize_raw(initial(old));old.step(2)
before=old._download_raw();count=old.step_count before=old._download_raw();count=old.step_count
old_mesh=RemapTopology.from_config(CFG,old.patch,A) old_mesh=RemapTopology.from_config(CFG,old.hierarchy.pairs[0],A)
new_mesh=RemapTopology.from_config(CFG,old.patch,B) new_mesh=RemapTopology.from_config(CFG,old.hierarchy.pairs[0],B)
expected=_remap_stored(old_mesh,new_mesh,before,shifted=True, expected=_remap_stored(old_mesh,new_mesh,before,shifted=True,
base_count=count,level_steps=(count,2*count,4*count),phase='postcollision_g') base_count=count,level_steps=(count,2*count,4*count),phase='postcollision_g')
new=owned.enter_context(old.remap_child(B)) new=owned.enter_context(old.remap_finest(B))
self.check_reset(new,expected.g) self.check_reset(new,expected.g)
self.snapshot_equal(new._download_raw(),_StoredPopulations(expected.g,True),'host raw fill') self.snapshot_equal(new._download_raw(),_StoredPopulations(expected.g,True),'host raw fill')
for a,b,om,nm in zip(before.values,new._download_raw().values,old_mesh.leaf_masks,new_mesh.leaf_masks): for a,b,om,nm in zip(before.values,new._download_raw().values,old_mesh.leaf_masks,new_mesh.leaf_masks):
@@ -103,7 +103,7 @@ class ShiftedRemapGPU(unittest.TestCase):
with CudaContext(0) as ctx,make(ctx) as run: with CudaContext(0) as ctx,make(ctx) as run:
run._initialize_raw(initial(run));run.step(2) run._initialize_raw(initial(run));run.step(2)
before=run._download_raw() before=run._download_raw()
with self.assertRaises(ValueError):run.remap_child([]) with self.assertRaises(ValueError):run.remap_finest([])
self.snapshot_equal(run._download_raw(),before,'rejection preserves source') self.snapshot_equal(run._download_raw(),before,'rejection preserves source')
self.assertEqual(run.step_count,2) self.assertEqual(run.step_count,2)
run.step();self.check_finite(run) run.step();self.check_finite(run)
+5 -4
View File
@@ -81,7 +81,7 @@ def public_jets(center, radius, scenario):
) )
def jet_interval(handle, scenario, begin, n, diameter): def jet_interval(handle, scenario, begin, n, diameter, *, level=None):
"""Apply a reversal only at a completed interval boundary, before advancing.""" """Apply a reversal only at a completed interval boundary, before advancing."""
if scenario != "none": if scenario != "none":
boundaries = [round(t * diameter / SPEED) for t in (250, 260, 270, 280, 290, 300, 305)] boundaries = [round(t * diameter / SPEED) for t in (250, 260, 270, 280, 290, 300, 305)]
@@ -92,8 +92,9 @@ def jet_interval(handle, scenario, begin, n, diameter):
if begin < switch: if begin < switch:
return min(n, switch - begin), None return min(n, switch - begin), None
if begin == switch: if begin == switch:
handle.set_jet_amplitude(0, -JET_PEAK) target = {} if level is None else dict(level=level)
handle.set_jet_amplitude(1, JET_PEAK) handle.set_jet_amplitude(0, -JET_PEAK, **target)
handle.set_jet_amplitude(1, JET_PEAK, **target)
return n, dict(step=begin, tstar=300.0, amplitudes=[-JET_PEAK, JET_PEAK]) return n, dict(step=begin, tstar=300.0, amplitudes=[-JET_PEAK, JET_PEAK])
return n, None return n, None
@@ -448,7 +449,7 @@ def run_case(a, h, periodic, case):
) )
n = min(max(1, round(width * p["diameter"] / SPEED)), max_steps - begin) n = min(max(1, round(width * p["diameter"] / SPEED)), max_steps - begin)
n, change = jet_interval(backend if refined else sim, a.jet_scenario, n, change = jet_interval(backend if refined else sim, a.jet_scenario,
begin, n, p["diameter"]) begin, n, p["diameter"], level=1 if refined else None)
if change is not None: if change is not None:
amplitude_changes.append(change) amplitude_changes.append(change)
impulse = np.zeros(2) impulse = np.zeros(2)
+1 -1
View File
@@ -73,7 +73,7 @@ def run_variant(backend, variant, n_intervals, deadline):
else: else:
for _ in range(INTERVAL): for _ in range(INTERVAL):
backend.step(1) backend.step(1)
totals.append(components(backend.read_exchange_window(0))) totals.append(components(backend.read_exchange_window(0, level=1)))
t.check_deadline(deadline) t.check_deadline(deadline)
return np.array(totals), (np.array(scales) if scales else None) return np.array(totals), (np.array(scales) if scales else None)