Compare commits
1
Commits
| Author | SHA1 | Date | |
|---|---|---|---|
|
|
f75d543586 |
@@ -123,7 +123,7 @@ with CudaContext(0) as ctx, FixedPatchStepper(cfg, (8, 6, 24, 18), context=ctx,
|
||||
```
|
||||
|
||||
Scope and limits:
|
||||
- D2Q9, SRT, FP32, double buffer, unshifted, no LES or forcing; fixed 2:1 patch(es), two levels (or one nested third).
|
||||
- D2Q9, SRT, FP32, double buffer, unshifted, no LES or forcing; fixed 2:1 patches, two or three levels with separated rectangles.
|
||||
- `coarse_boundary="periodic"` (default) or `"slot_reservoir"` (bottom slot, top
|
||||
reservoir; patch needs `y0 >= 4`, `y1 <= ny-4`). Periodic mode ignores the wall/slot config fields.
|
||||
- Object centers/radii use the coarse **cell-edge** frame (cell `I` spans `[I, I+1)`),
|
||||
@@ -144,7 +144,8 @@ Scope and limits:
|
||||
not matter. `fine_objects`/`jets` work as for one patch: each circle (same
|
||||
coarse cell-edge frame) must lie wholly inside one patch's active region with
|
||||
its boundary stencil clear of transfer cells; a jet may not select links of two
|
||||
different bodies, even across patches. No wall-touching patches or nesting yet.
|
||||
different bodies, even across patches. Wall-touching patches are unsupported.
|
||||
Three-level nesting uses the explicit child coordinates described below.
|
||||
|
||||
```python
|
||||
cfg = LBMConfig(nx=32, ny=16, viscosity=.1)
|
||||
@@ -164,7 +165,10 @@ Scope and limits:
|
||||
ACTIVE cells, away from L1's ghost band and transfer footprint (checked before
|
||||
any GPU allocation). `initialize` and `download` take/return three arrays; each
|
||||
level is NaN on its non-ACTIVE (and child-covered) cells. One coarse step = 2 L1
|
||||
steps = 4 L2 steps. One L1 patch only. Optional stationary circles and Jets
|
||||
steps = 4 L2 steps. This scalar `child=` form uses one L1 patch.
|
||||
For separated L1/L2 sets, use `child_rectangles_in_parent`; see
|
||||
[coordinate frames and supported layouts](docs/fixed_three_level_multipatch.md).
|
||||
Optional stationary circles and Jets
|
||||
live on L2 via `child_objects` and `child_jets`; `fine_objects`/`jets` retain
|
||||
their L1 meaning and must be empty in this mode. All object and Jet geometry
|
||||
uses the **L0 edge frame**, not the child bounds' L1 window frame. L2 has
|
||||
@@ -192,14 +196,26 @@ Scope and limits:
|
||||
run.step(4)
|
||||
coarse_g, l1_g, l2_g = run.download() # solids/covered/nonleaves are NaN
|
||||
```
|
||||
- Manual synchronized L2 replacement: `replacement = current.remap_child(rectangles,
|
||||
child_objects=original_objects)` preserves the absolute step count and live Jet
|
||||
amplitudes. L0/L1 and whole body-containing L2 rectangles stay fixed; separated
|
||||
fluid rectangles may change. The source remains independently owned. See
|
||||
[ownership, pinning and transfer rules](docs/synchronized_leaf_remap.md).
|
||||
- Runnable three-level example: `PYTHONPATH=src python examples/multigrid_remap.py`.
|
||||
It uses two separated L1 regions, an L2 stationary Circle/Jet and an explicit
|
||||
fluid-rectangle move, retaining original root-frame body definitions. Add
|
||||
`--check-geometry` for a CPU-only check or `--verify-jet` for an optional
|
||||
next-step comparison with an explicitly driven control. No flowfiles are saved.
|
||||
- Not available: per-link surface exchange/force readout in blocked mode,
|
||||
checkpoints, moving bodies, L1 bodies/Jets in a three-level chain,
|
||||
several patches per level with nesting, dynamic refinement, `Simulation`
|
||||
integration. Validation so far is storage/operator equivalence and short runs,
|
||||
not accuracy, convergence or long-time stability. The new L2-body GPU gates
|
||||
(layouts, masks, amplitude reversal, legacy equivalence) and independent
|
||||
every-substep numerical replay remain required acceptance checks; CPU tests
|
||||
alone do not establish GPU numerical correctness.
|
||||
checkpoints, moving bodies, L1 bodies/Jets in a three-level hierarchy,
|
||||
arbitrary depth, automatic AMR or `Simulation` integration.
|
||||
Storage/operator equivalence, remap and short-run gates do not establish
|
||||
convergence or long-time physical accuracy. The earlier shear-decay result is
|
||||
a research accuracy limitation, not a demonstrated solver bug. Pinned-body
|
||||
replay conditioned on actual L1 feeds is not an independent whole-body solver.
|
||||
See [three-level validation scope](docs/fixed_three_level_multipatch.md) and
|
||||
[remap validation scope](docs/synchronized_leaf_remap.md). The new usage example
|
||||
has CPU geometry/syntax checks only until its GPU command is run on Horsea.
|
||||
|
||||
## Installation
|
||||
|
||||
|
||||
@@ -1,73 +0,0 @@
|
||||
# Raw multilevel boundary exchange (opt-in candidate)
|
||||
|
||||
Use `FixedPatchStepper(..., observe_raw_exchange=True)` and, after a successful
|
||||
base step, `read_raw_exchange(body_id, level=l, window="last_base_step")`.
|
||||
The alternative window is `"last_substep"`. Body IDs are object positions within
|
||||
the named level, not the object's original user label. The torque reference is
|
||||
that circle's geometry.center (the center used to construct BFL arms), reported
|
||||
in the root edge frame. A differing SimObject.center does not change this reference. Current geometry
|
||||
restrictions are unchanged: fixed complete bodies on the finest level, isolated
|
||||
from the transfer footprint. New observation currently requires blocked body
|
||||
storage; dense-body configurations raise before device allocation. Existing
|
||||
dense diagnostics are unchanged.
|
||||
|
||||
Each result contains `wall`, `aperture`, `total`, with `force`, `torque`,
|
||||
`impulse`, `angular_impulse`, `mass_into_fluid`, and `mass_rate`. These are raw
|
||||
boundary momentum reaction, not validated stress or whole-nozzle mechanical
|
||||
load. No GI/advective correction or guessed thrust is added. Apertures retain
|
||||
their geometric classification when their prescribed velocity is zero.
|
||||
|
||||
All results are base-lattice units per unit out-of-plane depth. For local
|
||||
spacing and timestep h=dt, raw `(jx,jy,tz,dm)` converts to impulse
|
||||
`(h*dt*jx,h*dt*jy,h*h*dt*tz)` and mass `h*h*dm`.
|
||||
The accumulator intentionally accepts h and dt separately to make the spatial
|
||||
and temporal conversion explicit; the current acoustic hierarchy supplies dt=h.
|
||||
The stored rx/ry arms run from geometry.center to the actual wall hit, not to
|
||||
the fluid node. For this raw exchange j=c_i(a+b), moving the arm from the fluid
|
||||
node to the hit adds q*c_i, whose cross product with j is zero. Thus either arm
|
||||
gives the same raw torque (up to floating-point rounding); this equivalence
|
||||
must not be assumed for a future nonparallel stress/advective correction.
|
||||
|
||||
Force/torque are impulse/angular impulse divided by the reported window length.
|
||||
`dm=b-a` is positive into the exterior fluid, while j is the opposite of the
|
||||
fluid's momentum change across that cut link. Interpolated BFL can introduce
|
||||
mass drift; dm is an observation, not an assertion of impermeability.
|
||||
|
||||
`last_base_step` means [t-1,t] even after step(n); `last_substep` means [t-h,t]
|
||||
and is not a common instantaneous sample across levels. A Jet amplitude set
|
||||
before step(1) applies to all its substeps. Reads are repeatable; step(0) leaves
|
||||
them unchanged. Initialization and a remap candidate invalidate observation.
|
||||
Only successful stream synchronization publishes a window; failed instances
|
||||
reject reads rather than returning partially advanced measurements.
|
||||
|
||||
Storage is one FP32 [N,4] raw array and one FP64 [N,4] integral per measured
|
||||
level, plus row metadata. No substep history is stored. Same-stream sequencing
|
||||
is halo/restriction -> BFL -> independent row accumulation -> step. There are
|
||||
no accumulation atomics or body obs offsets (blocked n_objects remains zero).
|
||||
Disabled observation allocates neither observation buffer and launches no
|
||||
accumulation kernel. The standalone periodic block caller passes explicit
|
||||
null optional arguments.
|
||||
|
||||
## Gates
|
||||
|
||||
CPU (standard library unittest, plus existing NumPy dependency):
|
||||
|
||||
PYTHONPATH=src python -m unittest discover -s tests/unit -p test_raw_exchange.py -v
|
||||
|
||||
GPU, explicitly on a CUDA host with the existing environment:
|
||||
|
||||
PYTHONPATH=src python -m unittest discover -s tests/validation/multigrid -p test_raw_exchange_gpu.py -v
|
||||
|
||||
The GPU suite includes an independent actual-population cut-link census and
|
||||
ACTIVE-fluid momentum/mass ledger with non-body boundary flux, observation
|
||||
on/off population equivalence, dense/blocked row comparison, three block
|
||||
partitions, zero/step/reversed/disabled Jet amplitudes, depth-four eight-sample
|
||||
accumulation, and a short-time whole-domain uniform-fine comparison with the
|
||||
same physical geometry, viscosity and Jet velocity. The uniform reference
|
||||
uses a test-only adapter to the standalone kernel's optional inputs.
|
||||
|
||||
These are implementation checks. They do not establish physical Jet/stress
|
||||
accuracy, long-time AMR error, or convergence order. Resolution/long-time
|
||||
studies and independently validated continuum load benchmarks remain separate.
|
||||
GPU gates were prepared but not run in the authoring environment (no CUDA,
|
||||
PyCUDA or nvcc); run them before merging this candidate.
|
||||
@@ -1,8 +1,9 @@
|
||||
# Synchronized leaf G remap and manual continuation
|
||||
|
||||
This is a CPU state-transfer component, **not a runtime adaptive solver**.
|
||||
It changes only a separated L2 rectangle set in a three-level, body-free,
|
||||
periodic hierarchy. L0/L1 topology, relaxation parameters, direction order,
|
||||
It changes only a separated L2 rectangle set in a three-level periodic
|
||||
hierarchy, either body-free or with stationary Circles/Jets in unchanged pinned
|
||||
L2 rectangles (see below). L0/L1 topology, relaxation parameters, direction order,
|
||||
physical frame and completed time remain unchanged. Existing geometry checks
|
||||
validate child support/nesting and separation. The CPU helper itself adds no CUDA, scheduler, allocator, boundary handling,
|
||||
collision or ghost refill. A bounded runtime adapter is described below.
|
||||
@@ -29,11 +30,12 @@ origin is retained. Physical lower edges are `origin + spacing*(x,y)` on L0,
|
||||
`origin + spacing*(parent.fine_origin + (x,y)/4)` on L2.
|
||||
|
||||
Factory options `coarse_boundary` and `objects` must declare `"periodic"` and
|
||||
an empty body collection (all levels). LBMConfig's channel settings do not
|
||||
an empty L0/L1 body collection; use `child_objects`/`child_jets` for pinned L2
|
||||
bodies. LBMConfig's channel settings do not
|
||||
select this operator's boundary mode. Geometry metadata is factory-created;
|
||||
manual mutation/forging of metadata is unsupported. Caller certifies the
|
||||
actual state has no bodies/forcing, that all device work and old readers are
|
||||
quiescent, and that inputs are the **current owned-leaf postcollision G**.
|
||||
actual state has no undeclared bodies or forcing, that all device work and old
|
||||
readers are quiescent, and that inputs are the **current owned-leaf postcollision G**.
|
||||
Count/phase stamps check consistency, not freshness or active-buffer identity.
|
||||
|
||||
`current_g` is three NumPy Q-major arrays `(9,ny,nx)`. For blocked inputs,
|
||||
@@ -73,8 +75,9 @@ PYTHONPATH=src python -m unittest discover -s tests/unit -p 'test_synchronized_r
|
||||
|
||||
`replacement = current.remap_child(rectangles)` returns a new ready
|
||||
`FixedPatchStepper` at the **same absolute step_count**. It supports only
|
||||
body-free periodic three-level runs whose existing L2 uses
|
||||
`child_rectangles_in_parent`. Scalar or union L1 parents are supported; scalar
|
||||
periodic three-level runs whose existing L2 uses
|
||||
`child_rectangles_in_parent`, either body-free or under the pinned-body rules below.
|
||||
Scalar or union L1 parents are supported; scalar
|
||||
`child=` is deliberately rejected. L0/L1, relaxation and block sizes stay fixed.
|
||||
No tagger, scheduler, in-place replacement, or general rollback is provided.
|
||||
|
||||
@@ -105,9 +108,9 @@ and geometry, and candidate/source failure cleanup. GPU tests are predeclared
|
||||
in `tests/validation/multigrid/test_runtime_child_remap_gpu.py`: identity next
|
||||
step bitwise, shifted A→B→A with two intervals between changes, uniform moving
|
||||
and nonequilibrium states, exact CPU remap and overlap, both buffers reset,
|
||||
finite leaves, exact masks and absolute counts. They were not run on the CPU
|
||||
host. The separate independent FP64 capture/replay gate is max population
|
||||
error **3e-6**, with exact masks; remap mass/momentum allow FP32 roundoff.
|
||||
finite leaves, exact masks and absolute counts. These are GPU execution gates,
|
||||
not CPU-mock results. The separate body-free independent FP64 capture/replay
|
||||
gate is max population error **3e-6**, with exact masks; remap mass/momentum allow FP32 roundoff.
|
||||
Per-population conservation is not an advancement target under collision.
|
||||
|
||||
### Bounded pinned-body extension
|
||||
@@ -155,6 +158,28 @@ Multiple-body CPU coverage: `tests/unit/test_multiple_pinned_remap.py`.
|
||||
two independently updated Jets, variable fluid-region counts and add/move/remove
|
||||
immediate comparisons with CPU transfer. All nine solid-excluded area-weighted
|
||||
population integrals, exact overlap and nonleaf NaNs are checked.
|
||||
These new GPU cases have **not been run here**. They do not establish evolution
|
||||
accuracy. A 3e-6 FP32-vs-FP64 evolution gate requires a separate independent
|
||||
BFL reference; the body-free replay gate above does not establish that claim.
|
||||
These transfer/equivalence cases do not establish evolution accuracy. Historical
|
||||
pinned-body FP64 replay uses actual L1 feeds: it checks conditional L2 evolution,
|
||||
not an independent whole-body solver or whole-hierarchy accuracy. The body-free
|
||||
replay gate above does not establish that stronger claim either.
|
||||
|
||||
## Small runnable usage example
|
||||
|
||||
[`examples/multigrid_remap.py`](../examples/multigrid_remap.py) keeps the original
|
||||
root-frame Circle/Jet definitions, moves a fluid-only L2 rectangle, checks fluid
|
||||
leaves with geometry-derived masks and retires both old/new owners before their
|
||||
shared context. `--verify-jet` adds a temporary control solely to compare the
|
||||
next step against an explicitly set live amplitude; normal use needs no control.
|
||||
|
||||
```sh
|
||||
# Host-only: no CUDA import or allocation.
|
||||
PYTHONPATH=src python examples/multigrid_remap.py --check-geometry
|
||||
# Horsea: select an approved idle physical GPU; it becomes logical device 0.
|
||||
PYTHONPATH="$PWD/src" CUDA_VISIBLE_DEVICES=0 python examples/multigrid_remap.py --verify-jet
|
||||
```
|
||||
|
||||
The new example has passed host geometry and syntax/import checks only; its GPU
|
||||
run is pending Horsea. Expected scope: base counts 2 before/after remap, then 5;
|
||||
finite positive-density fluid leaves, NaN nonleaves, and exact next-step control
|
||||
agreement. This is a short public-API smoke run, not a physical accuracy test.
|
||||
No flowfiles are written.
|
||||
|
||||
@@ -0,0 +1,126 @@
|
||||
"""Three levels, separated regions and an explicit fluid-region remap.
|
||||
|
||||
Run on an approved idle GPU (logical device 0):
|
||||
PYTHONPATH=src CUDA_VISIBLE_DEVICES=0 python examples/multigrid_remap.py --verify-jet
|
||||
Host-only geometry check, without CUDA:
|
||||
PYTHONPATH=src python examples/multigrid_remap.py --check-geometry
|
||||
|
||||
A short public-API smoke example, not an accuracy/stability benchmark. No files
|
||||
are written. Edit the rectangles below; this is manual remapping, not auto AMR.
|
||||
"""
|
||||
import argparse
|
||||
from contextlib import ExitStack
|
||||
|
||||
import numpy as np
|
||||
|
||||
from CelerisLab.body.geometry.circle import CircleGeometry
|
||||
from CelerisLab.body.jets import Jet, Parallel, Tube
|
||||
from CelerisLab.body.objects import SimObject
|
||||
from CelerisLab.config import LBMConfig
|
||||
from CelerisLab.lbm.equilibrium import compute_feq_d2q9
|
||||
from CelerisLab.lbm.refinement import FixedPatchSet2D
|
||||
from CelerisLab.lbm.synchronized_remap import RemapTopology
|
||||
|
||||
CFG = LBMConfig(nx=48, ny=32, viscosity=.1) # D2Q9/SRT/FP32, periodic
|
||||
PARENTS = ((4, 4, 20, 28), (24, 4, 44, 28)) # L0 cell-edge rectangles
|
||||
PIN = (16, 16, 32, 24) # L1 allocated-grid indices; whole body region stays fixed
|
||||
BEFORE = (PIN, (56, 32, 72, 40))
|
||||
AFTER = (PIN, (60, 32, 76, 40)) # move only the separated fluid rectangle
|
||||
LIVE_AMPLITUDE = -.004
|
||||
|
||||
|
||||
def geometry():
|
||||
# Retain these original L0 edge-frame definitions for every remap; never
|
||||
# substitute the stepper's projected fine-local boundary objects.
|
||||
objects = (SimObject(0, geometry=CircleGeometry(12., 10., .8),
|
||||
center=(12., 10.), radius=.8),)
|
||||
jets = (Jet(Tube((10.5, 10.), 0., 3., .3), Parallel(0.), amplitude=.002),)
|
||||
parent = FixedPatchSet2D.from_config(CFG, PARENTS)
|
||||
meshes = tuple(RemapTopology.from_config(
|
||||
CFG, parent, rectangles, child_objects=objects, child_jets=jets)
|
||||
for rectangles in (BEFORE, AFTER))
|
||||
if any(mesh.pinned_rectangles != (PIN,) for mesh in meshes):
|
||||
raise RuntimeError("the body-containing rectangle must remain pinned")
|
||||
return objects, jets, meshes
|
||||
|
||||
|
||||
def diagnostics(run, mesh):
|
||||
"""Masks come from geometry, never from whether downloaded values are finite."""
|
||||
mass = 0.
|
||||
for level, (g, leaf) in enumerate(zip(run.download(), mesh.leaf_masks)):
|
||||
fluid = g[:, leaf].astype(np.float64)
|
||||
rho = fluid.sum(axis=0)
|
||||
if not np.isfinite(fluid).all() or np.any(rho <= 0):
|
||||
raise RuntimeError(f"invalid L{level} fluid populations/density")
|
||||
if not np.isnan(g[:, ~leaf]).all():
|
||||
raise RuntimeError(f"unexpected L{level} nonleaf values")
|
||||
mass += rho.sum() / 4**level # area weights: 1, 1/4, 1/16
|
||||
print(f"absolute base step {run.step_count}: finite fluid leaves, mass={mass:.9f}")
|
||||
# Driven boundary injects/extracts mass; mass need not stay constant in time.
|
||||
|
||||
|
||||
def verify_live_jet(current, objects):
|
||||
"""Optional public-API check; normal continuation needs no control solver."""
|
||||
# Identity remap copies the same state. Explicitly setting the known live
|
||||
# amplitude on the control detects a stale amplitude on the continuation.
|
||||
# No private fields or constructor-time Jet snapshots inspect the solver.
|
||||
with current.remap_child(AFTER, child_objects=objects) as control:
|
||||
control.set_jet_amplitude(0, LIVE_AMPLITUDE, level=2)
|
||||
current.step()
|
||||
control.step()
|
||||
for actual, expected in zip(current.download(), control.download()):
|
||||
if not np.array_equal(actual, expected, equal_nan=True):
|
||||
raise RuntimeError("live Jet continuation differs from explicitly set control")
|
||||
print("live Jet amplitude preserved: next step equals explicitly set control")
|
||||
|
||||
|
||||
def main():
|
||||
parser = argparse.ArgumentParser(description=__doc__.splitlines()[0])
|
||||
parser.add_argument("--check-geometry", action="store_true")
|
||||
parser.add_argument("--verify-jet", action="store_true",
|
||||
help="compare one step with an explicitly driven control")
|
||||
args = parser.parse_args()
|
||||
objects, jets, (old_mesh, new_mesh) = geometry()
|
||||
if args.check_geometry:
|
||||
print("geometry OK: three levels, two L1 regions, fixed body pin, moved fluid region")
|
||||
return
|
||||
|
||||
# Lazy imports keep --check-geometry useful on a CPU-only machine.
|
||||
from CelerisLab.cuda.context import CudaContext
|
||||
from CelerisLab.lbm.refinement_stepper import FixedPatchStepper
|
||||
|
||||
# ExitStack owns every old/new stepper, including on failure. All are closed
|
||||
# before the shared context. remap_child never consumes/closes its source.
|
||||
with CudaContext(0) as context, ExitStack() as owned:
|
||||
old = owned.enter_context(FixedPatchStepper(
|
||||
CFG, PARENTS, context=context, base_block=(4, 4), fine_block=(4, 4),
|
||||
child_block=(4, 4), child_rectangles_in_parent=BEFORE,
|
||||
child_objects=objects, child_jets=jets))
|
||||
rest = np.asarray(compute_feq_d2q9(1., 0., 0.), np.float32)[:, None, None]
|
||||
old.initialize(*(np.broadcast_to(rest, (9, *mask.shape))
|
||||
for mask in old_mesh.leaf_masks))
|
||||
old.step(2) # one base step advances L1 twice and L2 four times
|
||||
old.set_jet_amplitude(0, LIVE_AMPLITUDE, level=2)
|
||||
diagnostics(old, old_mesh)
|
||||
count = old.step_count
|
||||
current = owned.enter_context(old.remap_child(AFTER, child_objects=objects))
|
||||
if current is old or current.step_count != count:
|
||||
raise RuntimeError("remap must return a distinct continuation at the same time")
|
||||
diagnostics(current, new_mesh)
|
||||
|
||||
if args.verify_jet:
|
||||
verify_live_jet(current, objects) # optional control; advances one step
|
||||
else:
|
||||
current.step()
|
||||
if old.step_count != count or current.step_count != count + 1:
|
||||
raise RuntimeError("old/new absolute step ownership changed")
|
||||
old.close() # retire the old state only after successful replacement
|
||||
current.step(2)
|
||||
if current.step_count != count + 3:
|
||||
raise RuntimeError("continuation lost its absolute step count")
|
||||
diagnostics(current, new_mesh)
|
||||
print("smoke example only: no accuracy or long-time stability claim")
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
@@ -9,11 +9,6 @@ whose links are evaluated; callers convert frames/units before evaluation.
|
||||
One physical link (receiver, direction) may be driven by at most one jet, and
|
||||
one jet may drive the links of only one body (checked at assignment, any amplitude).
|
||||
|
||||
``condition`` names the prescribed boundary type. Only ``"velocity"`` (the
|
||||
default: amplitude is a signed lattice velocity, not a flow rate) is
|
||||
implemented; ``"pressure"`` is reserved and rejected at construction, before
|
||||
any amplitude validation, so it never reaches assignment or a kernel.
|
||||
|
||||
Edge rule: the lateral bound |s| < 1 is exclusive (the profile vanishes
|
||||
there for parabolic/quartic); Sector r_min/r_max and Tube 0 <= along <= length
|
||||
are inclusive. Jets are immutable; an amplitude change replaces the Jet.
|
||||
@@ -30,9 +25,6 @@ from dataclasses import dataclass, replace
|
||||
|
||||
import numpy as np
|
||||
|
||||
CONDITIONS = ("velocity",)
|
||||
RESERVED_CONDITIONS = ("pressure",)
|
||||
|
||||
PROFILES = {
|
||||
"uniform": lambda s: np.ones_like(s),
|
||||
"parabolic": lambda s: 1 - s*s,
|
||||
@@ -130,8 +122,7 @@ class Parallel:
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class Jet:
|
||||
"""Region (Sector|Tube), field (Radial|Parallel), profile name, signed amplitude,
|
||||
prescribed boundary condition (velocity only; pressure reserved).
|
||||
"""Region (Sector|Tube), field (Radial|Parallel), profile name, signed amplitude.
|
||||
|
||||
Frozen: change amplitude with ``with_amplitude`` (validated, new object);
|
||||
the cached per-link shapes stay valid because region/field/profile cannot change.
|
||||
@@ -140,14 +131,8 @@ class Jet:
|
||||
field: object
|
||||
profile: str = "uniform"
|
||||
amplitude: float = 0.0
|
||||
condition: str = "velocity"
|
||||
|
||||
def __post_init__(self):
|
||||
# The condition decides what amplitude and field mean: check it first.
|
||||
if self.condition in RESERVED_CONDITIONS:
|
||||
raise NotImplementedError(f"{self.condition} jets are reserved but not implemented")
|
||||
if self.condition not in CONDITIONS:
|
||||
raise ValueError(f"Jet condition must be one of {CONDITIONS + RESERVED_CONDITIONS}")
|
||||
if self.profile not in PROFILES:
|
||||
raise ValueError(f"Jet profile must be one of {sorted(PROFILES)}")
|
||||
_finite(self.amplitude)
|
||||
|
||||
@@ -211,7 +211,6 @@ class PeriodicBlockStepper:
|
||||
t = self._tables
|
||||
self._curved(self._current, t["block"], t["cell"], t["dir"], t["q"],
|
||||
t["fallback"], np.uint32(self._links), np.uintp(0), # no actuation
|
||||
np.uintp(0), np.uintp(0), np.uintp(0), # no observation
|
||||
grid=self._link_grid, **kw)
|
||||
self._step(self._current, self._spare, np.uint32(self.layout.halo),
|
||||
self._tables.get("flags", np.uintp(0)), np.uintp(0), # no kind
|
||||
|
||||
@@ -1,100 +0,0 @@
|
||||
# CelerisLab/lbm/hierarchy.py
|
||||
"""Host-only geometry of a fixed 2:1 refinement hierarchy of any depth (numpy, no CUDA).
|
||||
|
||||
Loops the existing adjacent-pair logic of refinement.py; it adds no transfer rule.
|
||||
``refinements[l]`` describes level l+1 in level-l integer cells: a scalar
|
||||
(x0, y0, x1, y1) is one FixedPatch2D window (its level keeps the window shape), a
|
||||
sequence of rectangles is a FixedPatchSet2D on the level-wide grid (origin 0).
|
||||
A scalar child is only accepted under a scalar parent. Every pair after the first
|
||||
is built by nest_patch/nest_patch_set, whose _nested_kind validates coverage,
|
||||
donors/receivers, surviving pulls and the parent's upward restriction footprint.
|
||||
"""
|
||||
from dataclasses import dataclass, replace
|
||||
from numbers import Integral
|
||||
|
||||
import numpy as np
|
||||
|
||||
from .refinement import (ACTIVE, OUTER_GHOST, FixedPatch2D, FixedPatchSet2D,
|
||||
_middle_config, nest_patch, nest_patch_set)
|
||||
|
||||
|
||||
@dataclass(frozen=True, eq=False)
|
||||
class Hierarchy:
|
||||
"""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)."""
|
||||
configs: tuple
|
||||
pairs: tuple
|
||||
kinds: tuple
|
||||
|
||||
@property
|
||||
def depth(self):
|
||||
return len(self.kinds)
|
||||
|
||||
|
||||
def is_rectangle(value):
|
||||
value = tuple(value)
|
||||
return len(value) == 4 and all(isinstance(v, (Integral, np.integer))
|
||||
and not isinstance(v, (bool, np.bool_)) for v in value)
|
||||
|
||||
|
||||
def build_hierarchy(cfg, refinements, *, wall_aware=False):
|
||||
"""Geometry for len(refinements) + 1 levels; wall_aware marks base rows 0 and -1 OUTER."""
|
||||
configs, pairs, kinds = [_validated(cfg)], [], []
|
||||
for level, entry in enumerate(refinements):
|
||||
scalar = is_rectangle(entry)
|
||||
if not scalar and not tuple(entry):
|
||||
raise ValueError("every refinement level needs at least one rectangle")
|
||||
if level == 0:
|
||||
pair = (FixedPatch2D if scalar else FixedPatchSet2D).from_config(configs[0], entry)
|
||||
kind = np.where(pair.covered, OUTER_GHOST, ACTIVE).astype(np.uint8)
|
||||
if wall_aware:
|
||||
kind[[0, -1]] = OUTER_GHOST
|
||||
elif scalar:
|
||||
if not isinstance(pairs[-1], FixedPatch2D):
|
||||
raise ValueError("a scalar child requires a scalar parent; use a rectangle list")
|
||||
pair, kind = nest_patch(configs[level - 1], pairs[-1], entry)
|
||||
else:
|
||||
pair, kind = nest_patch_set(configs[level - 1], pairs[-1], entry)
|
||||
kinds.append(np.asarray(kind))
|
||||
pairs.append(pair)
|
||||
configs.append(_validated(_middle_config(configs[level], pair)))
|
||||
_check_chain(configs[level], configs[level + 1], pair)
|
||||
if pairs:
|
||||
kinds.append(pairs[-1].fine_kind)
|
||||
else:
|
||||
kind = np.full((cfg.ny, cfg.nx), ACTIVE, np.uint8)
|
||||
if wall_aware:
|
||||
kind[[0, -1]] = OUTER_GHOST
|
||||
kinds.append(kind)
|
||||
return Hierarchy(tuple(configs), tuple(pairs), tuple(kinds))
|
||||
|
||||
|
||||
def _validated(cfg):
|
||||
"""LBMConfig.omega is computed only by validate(); replace() copies it stale."""
|
||||
cfg = replace(cfg)
|
||||
cfg.validate()
|
||||
return cfg
|
||||
|
||||
|
||||
def _check_chain(coarse, fine, pair):
|
||||
if not (np.isclose(pair.tau_coarse, 1 / coarse.omega, rtol=1e-12, atol=0)
|
||||
and np.isclose(pair.tau_fine, 1 / fine.omega, rtol=1e-12, atol=0)
|
||||
and tuple(pair.fine_shape) == (fine.ny, fine.nx)):
|
||||
raise RuntimeError("level config chain is inconsistent with the pair tables")
|
||||
|
||||
|
||||
def leaf_partition(h):
|
||||
"""Base-cell fraction published as ACTIVE over all levels (exact: powers of two).
|
||||
|
||||
Window (scalar) levels are offset by their fine_level_origin first.
|
||||
"""
|
||||
base = np.zeros(h.kinds[0].shape)
|
||||
oy = ox = 0
|
||||
for level, kind in enumerate(h.kinds):
|
||||
if level:
|
||||
pair = h.pairs[level - 1]
|
||||
py, px = (pair.fine_level_origin if isinstance(pair, FixedPatch2D) else (0, 0))
|
||||
oy, ox = 2 * oy + py, 2 * ox + px
|
||||
y, x = np.nonzero(kind == ACTIVE)
|
||||
np.add.at(base, ((y + oy) >> level, (x + ox) >> level), 0.25 ** level)
|
||||
return base
|
||||
@@ -76,12 +76,11 @@ extern "C" __global__ void BlockCoalesce(
|
||||
// cell, so the existing BFL writes the solid-source window image that this
|
||||
// block's next pull reads (no halo copy between this kernel and BlockStep).
|
||||
// surface_velocity (nullable, 2 floats per link): prescribed fluid velocity on
|
||||
// fixed geometry, same convention as CurvedBoundaryKernel. Optional raw rows
|
||||
// (jx,jy,torque,dm) are unique physical links, never body-observation slots.
|
||||
// fixed geometry, same convention as CurvedBoundaryKernel. No exchange output.
|
||||
extern "C" __global__ void BlockCurvedBoundary(
|
||||
fpxx* f, const uint32_t* block, const uint32_t* cell, const unsigned char* dir,
|
||||
const float* q, const unsigned char* fallback, unsigned int count,
|
||||
const float* surface_velocity, const float* rx, const float* ry, float* raw)
|
||||
const float* surface_velocity)
|
||||
{
|
||||
unsigned int row = blockIdx.x * blockDim.x + threadIdx.x;
|
||||
if (row >= count) return;
|
||||
@@ -91,24 +90,7 @@ extern "C" __global__ void BlockCurvedBoundary(
|
||||
if (u != 0.0f || v != 0.0f) { Uw = u; Vw = v; }
|
||||
}
|
||||
apply_bouzidi_link(dir[row], q[row], cell[row], f + (unsigned long)block[row] * NQ * TOTAL_CELLS,
|
||||
Uw, Vw, rx ? rx[row] : 0.0f, ry ? ry[row] : 0.0f,
|
||||
fallback[row], nullptr, -1, raw ? raw + 4ul * row : nullptr,
|
||||
raw ? raw + 4ul * row + 3 : nullptr);
|
||||
}
|
||||
|
||||
// Same-stream BFL -> this row-owned accumulation -> pull step. No atomics,
|
||||
// no cross-row reduction and no history allocation proportional to substeps.
|
||||
extern "C" __global__ void BlockAccumulateExchange(
|
||||
const float* raw, double* integral, unsigned int count, double h, double dt)
|
||||
{
|
||||
unsigned int row = blockIdx.x * blockDim.x + threadIdx.x;
|
||||
if (row >= count) return;
|
||||
const unsigned long k = 4ul * row;
|
||||
const double impulse = h * dt;
|
||||
integral[k] += impulse * (double)raw[k];
|
||||
integral[k + 1] += impulse * (double)raw[k + 1];
|
||||
integral[k + 2] += h * impulse * (double)raw[k + 2];
|
||||
integral[k + 3] += h * h * (double)raw[k + 3];
|
||||
Uw, Vw, 0.0f, 0.0f, fallback[row], nullptr, -1);
|
||||
}
|
||||
|
||||
// Pull/collide/store core cells only; halo cells of out are left untouched.
|
||||
|
||||
@@ -82,8 +82,7 @@ __device__ inline void apply_bouzidi_link(
|
||||
float Uw, float Vw,
|
||||
float rx, float ry,
|
||||
unsigned char fallback_class,
|
||||
float* __restrict__ obs, int body_id, float* link_exchange = nullptr,
|
||||
float* link_mass = nullptr)
|
||||
float* __restrict__ obs, int body_id, float* link_exchange = nullptr)
|
||||
{
|
||||
const unsigned int dir_opp = (unsigned int)opp_dir((int)dir);
|
||||
|
||||
@@ -121,10 +120,6 @@ __device__ inline void apply_bouzidi_link(
|
||||
// corrected incoming population into the adjacent fluid node.
|
||||
store_ddf(fi, index_f(k_s, dir_opp), f_reflected);
|
||||
|
||||
// Equal opposite weights cancel even with shifted DDFs.
|
||||
if (link_mass != nullptr)
|
||||
*link_mass = f_reflected - f_toward;
|
||||
|
||||
if (obs != nullptr || link_exchange != nullptr) {
|
||||
#if USE_DDF_SHIFTING
|
||||
const float w_i = (float)d_w[dir];
|
||||
|
||||
@@ -1,36 +0,0 @@
|
||||
"""Host-only reduction of raw cut-link observations; no stress/GI correction.
|
||||
|
||||
Rows are unique physical links. All quantities are 2D, per unit out-of-plane
|
||||
thickness; centers are in the root edge frame, torque is about the named body.
|
||||
"""
|
||||
import numpy as np
|
||||
|
||||
|
||||
def summarize_exchange(rows, body_ids, aperture, *, body_id, level, h,
|
||||
time_end, window, reference_center):
|
||||
"""Reduce FP64 interval integrals or convert FP32 last-substep raw rows.
|
||||
|
||||
dt_l=h_l under acoustic 2:1 refinement. Local momentum j -> base force
|
||||
h*j -> impulse h*dt*j; torque adds one h; mass is h**2*dm.
|
||||
"""
|
||||
if window not in ("last_base_step", "last_substep"):
|
||||
raise ValueError("unknown exchange window")
|
||||
values = np.array(rows, dtype=np.float64, copy=True)
|
||||
duration = 1.0 if window == "last_base_step" else h
|
||||
if window == "last_substep":
|
||||
values *= np.array([h*h, h*h, h*h*h, h*h])
|
||||
selected = np.asarray(body_ids) == body_id
|
||||
aperture = np.asarray(aperture, bool)
|
||||
|
||||
def reduce(mask):
|
||||
integral = values[mask].sum(axis=0, dtype=np.float64)
|
||||
return dict(impulse=integral[:2].copy(), angular_impulse=float(integral[2]),
|
||||
mass_into_fluid=float(integral[3]), force=integral[:2] / duration,
|
||||
torque=float(integral[2] / duration), mass_rate=float(integral[3] / duration))
|
||||
|
||||
return dict(wall=reduce(selected & ~aperture), aperture=reduce(selected & aperture),
|
||||
total=reduce(selected), body_id=int(body_id), level=int(level),
|
||||
time_start=time_end-duration, time_end=time_end, dt=h,
|
||||
sample_count=round(1/h) if window == "last_base_step" else 1,
|
||||
window=window, units="base_lattice_per_unit_depth",
|
||||
definition="raw_boundary_exchange", reference_center=tuple(reference_center), valid=True)
|
||||
File diff suppressed because it is too large
Load Diff
@@ -1,107 +0,0 @@
|
||||
# CelerisLab/lbm/refinement_tables.py
|
||||
"""Host-only (numpy) mapping of fixed-patch transfer and boundary tables to block storage.
|
||||
|
||||
No CUDA import: CPU tests import these directly. refinement_stepper re-exports the names.
|
||||
"""
|
||||
from dataclasses import dataclass, replace
|
||||
|
||||
import numpy as np
|
||||
|
||||
from .blocks import blocked_link_tables
|
||||
from ..body.jets import jet_link_shapes, surface_link_velocity
|
||||
from .descriptors import D2Q9_EX, D2Q9_EY
|
||||
from .refinement import ACTIVE, OUTER_GHOST
|
||||
|
||||
|
||||
def blocked_fine_tables(patch, layout, order):
|
||||
"""Remap fixed-patch fine tables to owner-core blocked storage (host only).
|
||||
|
||||
Returns (explosion target block cells, coalescence child slots in the same
|
||||
rows/child order as ``patch.coalesce_fine[order]``, per-block-cell kind with
|
||||
unused halo entries OUTER). Coarse donors/targets are unchanged.
|
||||
"""
|
||||
owner = layout.owner.ravel().astype(np.intp)
|
||||
cells = patch.fine_shape[0] * patch.fine_shape[1]
|
||||
q, cell = np.divmod(patch.coalesce_fine[order], np.uint64(cells))
|
||||
kind = blocked_kind(layout, patch.fine_kind)
|
||||
return (owner[patch.ghost_fine.astype(np.intp)].astype(np.uint64),
|
||||
layout.slot(owner[cell.astype(np.intp)], q), kind)
|
||||
|
||||
|
||||
def blocked_kind(layout, dense_kind):
|
||||
"""Per-block-cell kind; unused halo entries and unstored holes stay OUTER."""
|
||||
owner = layout.owner.ravel().astype(np.intp)
|
||||
owned = layout.owned.ravel() # sparse multipatch levels leave holes unstored
|
||||
kind = np.full(layout.nblocks * int(np.prod(layout.local_shape)), OUTER_GHOST, np.uint8)
|
||||
kind[owner[owned]] = np.asarray(dense_kind).ravel()[owned]
|
||||
return kind
|
||||
|
||||
def blocked_base_tables(patch, layout, order, wall_aware=False, *, composed_kind=None):
|
||||
"""Parent donor cells and receiver-local arrival slots, in transfer row order.
|
||||
|
||||
Stage arrivals only after child reads and parent halo gather, immediately
|
||||
before parent pull. Source *owner* slots are wrong at block seams.
|
||||
Nested pairs must pass composed_kind to retain inherited INNER/OUTER and
|
||||
sparse holes; donors and receivers must be owned, uncovered ACTIVE cells.
|
||||
"""
|
||||
owner = layout.owner.ravel()
|
||||
q, receiver = np.divmod(patch.coalesce_coarse[order], np.uint64(owner.size))
|
||||
if composed_kind is None:
|
||||
dense_kind = np.where(patch.covered, OUTER_GHOST, ACTIVE).astype(np.uint8)
|
||||
if wall_aware:
|
||||
dense_kind[[0, -1]] = OUTER_GHOST
|
||||
else:
|
||||
dense_kind = np.asarray(composed_kind)
|
||||
if dense_kind.shape != patch.coarse_shape:
|
||||
raise ValueError("composed parent kind has the wrong shape")
|
||||
donors = patch.ghost_coarse.astype(np.intp)
|
||||
receivers = receiver.astype(np.intp)
|
||||
legal = layout.owned.ravel() & (dense_kind.ravel() == ACTIVE) & ~patch.covered.ravel()
|
||||
if not legal[donors].all() or not legal[receivers].all():
|
||||
raise ValueError("child donors/receivers need owned uncovered parent ACTIVE cells")
|
||||
offset = (np.asarray(D2Q9_EY)[q.astype(np.intp)] * layout.local_shape[1]
|
||||
+ np.asarray(D2Q9_EX)[q.astype(np.intp)])
|
||||
source = owner[receivers].astype(np.int64) - offset
|
||||
return (owner[donors].astype(np.uint64), layout.slot(source, q),
|
||||
blocked_kind(layout, dense_kind))
|
||||
|
||||
|
||||
# Required device columns plus optional actuation; observation arms are uploaded
|
||||
# separately only when enabled. body_id/aperture remain host reduction metadata.
|
||||
BLOCKED_BOUNDARY_DEVICE_COLUMNS = ("block", "cell", "dir", "q", "fallback", "surface_velocity")
|
||||
|
||||
|
||||
@dataclass
|
||||
class _BlockedBoundary:
|
||||
"""One level's stationary BFL data and mutable Jet velocity column."""
|
||||
host: object
|
||||
tables: dict
|
||||
jets: tuple
|
||||
actuated: bool
|
||||
jet_cache: object = None
|
||||
flags: object = None
|
||||
device: object = None
|
||||
kernel: object = None
|
||||
last_raw: object = None
|
||||
interval_integral: object = None
|
||||
accumulate: object = None
|
||||
|
||||
|
||||
def _prepare_blocked_boundary(boundary, layout, cfg, jets):
|
||||
if boundary is None or layout is None:
|
||||
return None
|
||||
fluid_idx, dirs, q, body_id, rx, ry, _, fallback = boundary.curved_host
|
||||
tables, _ = blocked_link_tables(layout, boundary.flags, fluid_idx, dirs, q, fallback)
|
||||
# blocked_link_tables/map_cut_links already rejects repeated (fluid_idx, dir),
|
||||
# independently of body_id, and repeated local writers.
|
||||
tables.update(body_id=np.asarray(body_id, np.int32), rx=np.asarray(rx, np.float32),
|
||||
ry=np.asarray(ry, np.float32), aperture=np.zeros(len(dirs), bool))
|
||||
actuated = bool(jets) or any(o.surface_actuation is not None for o in boundary.objects)
|
||||
if actuated:
|
||||
cfg = replace(cfg)
|
||||
cfg.validate()
|
||||
tables["surface_velocity"], tables["aperture"] = surface_link_velocity(
|
||||
boundary.objects, cfg, boundary.flags, fluid_idx, dirs, body_id, rx, ry,
|
||||
jets=jets, hits=boundary.curved_hits)
|
||||
return _BlockedBoundary(boundary, tables, jets, actuated,
|
||||
jet_link_shapes(jets, boundary.curved_hits) if jets else None)
|
||||
@@ -1 +0,0 @@
|
||||
"""CelerisLab test package; keeps shared test helpers independent of discovery paths."""
|
||||
@@ -1,97 +0,0 @@
|
||||
"""CPU geometry tests for lbm.hierarchy (no CUDA). Geometry construction only, not LBM numerics."""
|
||||
import unittest
|
||||
|
||||
import numpy as np
|
||||
from CelerisLab.config import LBMConfig
|
||||
from CelerisLab.lbm import refinement as r
|
||||
from CelerisLab.lbm.hierarchy import build_hierarchy, leaf_partition
|
||||
|
||||
CFG = LBMConfig(nx=48, ny=48, viscosity=.02)
|
||||
L1 = [(6, 6, 26, 26), (32, 30, 42, 42)] # L0 cells; L1 ACTIVE x/y 12..52 and x 64..84, y 60..84
|
||||
L2 = [(20, 20, 44, 44), (72, 68, 76, 76)] # L1 cells; L2 ACTIVE 40..88 and x 144..152, y 136..152
|
||||
L3 = [(56, 56, 72, 72)] # L2 cells
|
||||
FOUR = [L1, L2, L3]
|
||||
SMALL = LBMConfig(nx=16, ny=16, viscosity=.1)
|
||||
|
||||
|
||||
class FourLevels(unittest.TestCase):
|
||||
@classmethod
|
||||
def setUpClass(cls):
|
||||
cls.h = build_hierarchy(CFG, FOUR)
|
||||
|
||||
def test_config_chain(self):
|
||||
self.assertEqual(self.h.depth, 4)
|
||||
for level, c in enumerate(self.h.configs):
|
||||
self.assertEqual((c.ny, c.nx), (48 << level, 48 << level))
|
||||
self.assertAlmostEqual(c.viscosity, .02 * 2 ** level, places=12)
|
||||
self.assertEqual(self.h.kinds[level].shape, (c.ny, c.nx))
|
||||
for level, p in enumerate(self.h.pairs):
|
||||
self.assertAlmostEqual(p.tau_coarse, 1 / self.h.configs[level].omega, places=12)
|
||||
self.assertAlmostEqual(p.tau_fine, 1 / self.h.configs[level + 1].omega, places=12)
|
||||
|
||||
def test_partition_and_composition(self):
|
||||
h = self.h
|
||||
np.testing.assert_array_equal(leaf_partition(h), np.ones((48, 48)))
|
||||
for level in (1, 2):
|
||||
covered, raw = h.pairs[level].covered, h.pairs[level - 1].fine_kind
|
||||
self.assertTrue((h.kinds[level][covered] == r.OUTER_GHOST).all())
|
||||
np.testing.assert_array_equal(h.kinds[level][~covered], raw[~covered])
|
||||
np.testing.assert_array_equal(h.kinds[3], h.pairs[2].fine_kind)
|
||||
|
||||
def test_depth_is_list_length(self):
|
||||
for depth in range(1, 5):
|
||||
h = build_hierarchy(CFG, FOUR[:depth - 1])
|
||||
self.assertEqual(h.depth, depth)
|
||||
np.testing.assert_array_equal(leaf_partition(h), np.ones((48, 48)))
|
||||
for level in range(depth - 1):
|
||||
np.testing.assert_array_equal(h.kinds[level], self.h.kinds[level])
|
||||
|
||||
def test_illegal_coverage_reports_first_guard(self):
|
||||
cases = [([L1, [(10, 20, 30, 44)]], "inherited ACTIVE"),
|
||||
([L1, L2[:1] + [(54, 54, 58, 58)]], "inherited ACTIVE"),
|
||||
([L1, L2, [(40, 40, 60, 60)]], "donors/receivers"),
|
||||
([L1, L2, [(41, 41, 60, 60)]], "donors/receivers"),
|
||||
([L1, L2, [(56, 56, 64, 64), (66, 56, 72, 64)]], "closer than")]
|
||||
for levels, reason in cases:
|
||||
with self.subTest(levels=levels), self.assertRaisesRegex(ValueError, reason):
|
||||
build_hierarchy(CFG, levels)
|
||||
build_hierarchy(CFG, [L1, L2, [(42, 42, 60, 60)]]) # minimal legal margin
|
||||
|
||||
|
||||
class LegacyForms(unittest.TestCase):
|
||||
def test_scalar_parent_and_child_match_existing_calls(self):
|
||||
local = r.FixedPatch2D.from_config(SMALL, (4, 4, 12, 12))
|
||||
child, kind = r.nest_patch(SMALL, local, (8, 8, 16, 16))
|
||||
h = build_hierarchy(SMALL, [(4, 4, 12, 12), (8, 8, 16, 16)])
|
||||
np.testing.assert_array_equal(h.kinds[0], np.where(local.covered, r.OUTER_GHOST, r.ACTIVE))
|
||||
np.testing.assert_array_equal(h.kinds[1], kind)
|
||||
np.testing.assert_array_equal(h.kinds[2], child.fine_kind)
|
||||
for name in ("ghost_fine", "ghost_coarse", "coalesce_coarse", "coalesce_fine"):
|
||||
np.testing.assert_array_equal(getattr(h.pairs[1], name), getattr(child, name))
|
||||
np.testing.assert_array_equal(leaf_partition(h), np.ones((16, 16)))
|
||||
|
||||
def test_set_child_under_scalar_parent(self):
|
||||
local = r.FixedPatch2D.from_config(SMALL, (4, 4, 12, 12))
|
||||
child, kind = r.nest_patch_set(SMALL, local, [(8, 8, 16, 16)])
|
||||
h = build_hierarchy(SMALL, [(4, 4, 12, 12), [(8, 8, 16, 16)]])
|
||||
np.testing.assert_array_equal(h.kinds[1], kind)
|
||||
np.testing.assert_array_equal(h.pairs[1].coalesce_fine, child.coalesce_fine)
|
||||
np.testing.assert_array_equal(leaf_partition(h), np.ones((16, 16)))
|
||||
|
||||
def test_scalar_child_under_set_parent_is_rejected(self):
|
||||
with self.assertRaisesRegex(ValueError, "scalar parent"):
|
||||
build_hierarchy(SMALL, [[(4, 4, 12, 12)], (12, 12, 20, 20)])
|
||||
|
||||
def test_wall_rows_and_single_level(self):
|
||||
h = build_hierarchy(SMALL, [(4, 4, 12, 12)], wall_aware=True)
|
||||
self.assertTrue((h.kinds[0][[0, -1]] == r.OUTER_GHOST).all())
|
||||
one = build_hierarchy(SMALL, [])
|
||||
self.assertEqual(one.depth, 1)
|
||||
self.assertTrue((one.kinds[0] == r.ACTIVE).all())
|
||||
walled = build_hierarchy(SMALL, [], wall_aware=True) # regression: no refinement
|
||||
self.assertTrue((walled.kinds[0][[0, -1]] == r.OUTER_GHOST).all())
|
||||
self.assertTrue((walled.kinds[0][1:-1] == r.ACTIVE).all())
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
unittest.main()
|
||||
@@ -1,32 +0,0 @@
|
||||
"""Jet boundary-condition hook: velocity implemented, pressure reserved (host only, no CUDA)."""
|
||||
import pytest
|
||||
|
||||
from CelerisLab.body.jets import Jet, Parallel, Tube
|
||||
|
||||
|
||||
def region():
|
||||
return Tube((16., 12.), 0., 3., .7)
|
||||
|
||||
|
||||
def test_default_is_velocity_and_survives_amplitude_change():
|
||||
jet = Jet(region(), Parallel(0.), amplitude=.01)
|
||||
assert jet.condition == "velocity"
|
||||
assert jet.with_amplitude(-.02).condition == "velocity"
|
||||
# positional calls keep their meaning; condition is the new last field
|
||||
assert (Jet(region(), Parallel(0.), "parabolic", .01, "velocity")
|
||||
== Jet(region(), Parallel(0.), "parabolic", .01))
|
||||
|
||||
|
||||
def test_pressure_is_rejected_before_amplitude_validation():
|
||||
# An amplitude that the velocity check would reject must not decide the error.
|
||||
for amplitude in (0., 5., float("nan")):
|
||||
with pytest.raises(NotImplementedError, match="pressure jets are reserved"):
|
||||
Jet(region(), Parallel(0.), amplitude=amplitude, condition="pressure")
|
||||
with pytest.raises(NotImplementedError):
|
||||
Jet(region(), None, profile="no-such-profile", condition="pressure")
|
||||
|
||||
|
||||
def test_unknown_condition_is_a_value_error():
|
||||
for condition in ("Velocity", "flow_rate", "", None):
|
||||
with pytest.raises(ValueError, match="condition"):
|
||||
Jet(region(), Parallel(0.), amplitude=.01, condition=condition)
|
||||
@@ -1,170 +0,0 @@
|
||||
"""CPU reduction and real host scheduling with CUDA doubles, not GPU physics."""
|
||||
from contextlib import nullcontext
|
||||
from dataclasses import replace
|
||||
from unittest.mock import Mock
|
||||
import numpy as np
|
||||
import unittest
|
||||
from CelerisLab.body.geometry.circle import CircleGeometry
|
||||
from CelerisLab.body.objects import SimObject, SurfaceActuation
|
||||
from CelerisLab.config import LBMConfig
|
||||
from CelerisLab.lbm.blocks import BlockLayout2D
|
||||
from CelerisLab.lbm.refinement import FixedPatch2D, prepare_fine_boundary
|
||||
from CelerisLab.lbm.refinement_tables import _prepare_blocked_boundary
|
||||
from CelerisLab.lbm.raw_exchange import summarize_exchange
|
||||
from tests.unit.test_refinement_three_level_multipatch import host_module
|
||||
|
||||
|
||||
def body(center=(16., 16.), radius=.5, speed=0.):
|
||||
return SimObject(91, geometry=CircleGeometry(*center, radius), center=center, radius=radius,
|
||||
surface_actuation=SurfaceActuation(.7, .85, speed))
|
||||
|
||||
|
||||
def backend(observe=True, obj=None):
|
||||
host = host_module()
|
||||
calls = []
|
||||
module = Mock()
|
||||
module.get_global.return_value = (Mock(), 20)
|
||||
module.get_function.side_effect = lambda name: Mock(side_effect=lambda *a, **k: calls.append((name, a, k)))
|
||||
host.compiler.load_module.return_value = module
|
||||
host.cuda.mem_alloc.side_effect = lambda n: Mock(size=n)
|
||||
host.cuda.memcpy_dtoh.side_effect = lambda array, pointer: array.fill(0)
|
||||
ctx = Mock(activate=nullcontext, sm_arch='sm_70')
|
||||
run = host.FixedPatchStepper(LBMConfig(nx=48, ny=48, viscosity=.02), context=ctx,
|
||||
refinements=[[(6, 6, 26, 26)], [(20, 20, 44, 44)], [(56, 56, 72, 72)]],
|
||||
level_blocks=[(4, 4)]*4, level_objects=[(), (), (), (body() if obj is None else obj,)],
|
||||
observe_raw_exchange=observe)
|
||||
arrays = [np.full(l.shape, 1/9, np.float32) for l in run._levels]
|
||||
run.initialize(*arrays)
|
||||
return run, host, calls, arrays
|
||||
|
||||
|
||||
class RawExchange(unittest.TestCase):
|
||||
def _check_weights_split_and_window(self, level):
|
||||
h = 2. ** -level
|
||||
rows = np.array([[2, -3, 5, 7], [-1, 4, -2, -6], [99, 99, 99, 99]], np.float32)
|
||||
kw = dict(body_id=0, level=level, h=h, time_end=3, reference_center=(7, 8))
|
||||
last = summarize_exchange(rows, [0, 0, 1], [False, True, False], window='last_substep', **kw)
|
||||
integral = rows.astype(float)*[h*h, h*h, h**3, h*h]/h
|
||||
full = summarize_exchange(integral, [0, 0, 1], [False, True, False], window='last_base_step', **kw)
|
||||
np.testing.assert_array_equal(last['total']['force'], np.array([1, 1])*h)
|
||||
self.assertEqual(last['total']['torque'], 3 * h * h)
|
||||
self.assertEqual(last['total']['mass_into_fluid'], h * h)
|
||||
for name in ('force', 'torque', 'mass_rate'):
|
||||
np.testing.assert_array_equal(last['total'][name], full['total'][name])
|
||||
for name in ('force', 'impulse', 'torque', 'angular_impulse', 'mass_into_fluid'):
|
||||
np.testing.assert_array_equal(last['total'][name],
|
||||
np.asarray(last['wall'][name])+last['aperture'][name])
|
||||
self.assertEqual(last['time_start'], 3 - h)
|
||||
self.assertEqual(full['time_start'], 2)
|
||||
self.assertEqual(last['sample_count'], 1)
|
||||
self.assertEqual(full['sample_count'], 2 ** level)
|
||||
self.assertEqual(full['reference_center'], (7, 8))
|
||||
|
||||
def test_final_substep_is_not_interval_mean(self):
|
||||
kw = dict(body_id=0, level=1, h=.5, time_end=1, reference_center=(0, 0))
|
||||
last = summarize_exchange([[4, 0, 0, 0]], [0], [False], window='last_substep', **kw)
|
||||
full = summarize_exchange([[1.5, 0, 0, 0]], [0], [False], window='last_base_step', **kw)
|
||||
self.assertEqual(last['total']['force'][0], 2)
|
||||
self.assertEqual(full['total']['force'][0], 1.5)
|
||||
|
||||
def test_host_metadata_zero_jet_and_rows(self):
|
||||
cfg = LBMConfig(nx=32, ny=32, viscosity=.1)
|
||||
pair = FixedPatch2D.from_config(cfg, (4, 4, 28, 28))
|
||||
boundary = prepare_fine_boundary(cfg, pair, [body(radius=2.15)])
|
||||
local = replace(cfg, nx=pair.fine_shape[1], ny=pair.fine_shape[0], viscosity=.2)
|
||||
local.validate()
|
||||
layout = BlockLayout2D.build(*pair.fine_shape, 8, 8, periodic=(False, False))
|
||||
b = _prepare_blocked_boundary(boundary, layout, local, ())
|
||||
for key, idx in [('body_id', 3), ('rx', 4), ('ry', 5)]:
|
||||
np.testing.assert_array_equal(b.tables[key], boundary.curved_host[idx])
|
||||
self.assertTrue(b.tables['aperture'].any())
|
||||
self.assertTrue((~b.tables['aperture']).any())
|
||||
self.assertFalse(b.tables['surface_velocity'].any())
|
||||
# Classification does not disappear when prescribed velocity is zero.
|
||||
self.assertEqual(len(b.tables['cell']), len(boundary.curved_host[0]))
|
||||
|
||||
def test_four_level_schedule_publication_and_fixed_storage(self):
|
||||
run, host, calls, arrays = backend()
|
||||
b = run._level_boundaries[3]
|
||||
rows = len(b.tables['cell'])
|
||||
self.assertEqual(b.last_raw.size, rows * 16)
|
||||
self.assertEqual(b.interval_integral.size, rows * 32)
|
||||
with self.assertRaisesRegex(RuntimeError, 'completed'):
|
||||
run.read_raw_exchange(0, level=3)
|
||||
run.step()
|
||||
names = [name for name, _, _ in calls]
|
||||
self.assertEqual(names.count('BlockAccumulateExchange'), 8)
|
||||
for i, (name, args, kw) in enumerate(calls):
|
||||
if name == 'BlockAccumulateExchange':
|
||||
self.assertEqual(names[i - 1], 'BlockCurvedBoundary')
|
||||
self.assertEqual(names[i + 1], 'BlockStep')
|
||||
self.assertIs(args[0], b.last_raw)
|
||||
self.assertIs(args[1], b.interval_integral)
|
||||
self.assertEqual(args[3:], (0.125, 0.125))
|
||||
self.assertIs(kw['stream'], run._stream)
|
||||
self.assertEqual(host.cuda.memset_d8_async.call_count, 1)
|
||||
value = run.read_raw_exchange(0, level=3)
|
||||
self.assertEqual(value['sample_count'], 8)
|
||||
self.assertEqual(value['time_start'], 0)
|
||||
run.step(0)
|
||||
self.assertEqual(run.read_raw_exchange(0, level=3)['time_end'], 1)
|
||||
self.assertEqual(run.read_raw_exchange(0, level=3, window='last_substep')['time_start'], 0.875)
|
||||
run.step(2)
|
||||
self.assertEqual(run.read_raw_exchange(0, level=3)['time_start'], 2)
|
||||
run.initialize(*arrays)
|
||||
with self.assertRaisesRegex(RuntimeError, 'completed'):
|
||||
run.read_raw_exchange(0, level=3)
|
||||
run.close()
|
||||
|
||||
def test_torque_reference_uses_geometry_center_not_optional_object_center(self):
|
||||
for center in (None, (3., 4.)):
|
||||
obj = SimObject(91, geometry=CircleGeometry(16., 16., .5), center=center, radius=.5)
|
||||
run, _, _, _ = backend(obj=obj)
|
||||
run.step()
|
||||
result = run.read_raw_exchange(0, level=3)
|
||||
self.assertEqual(result['reference_center'], (16.0, 16.0))
|
||||
self.assertNotEqual(result['reference_center'], tuple(obj.center))
|
||||
run.close()
|
||||
|
||||
def test_disabled_has_nulls_and_no_observation_allocation_or_launch(self):
|
||||
run, host, calls, _ = backend(False)
|
||||
b = run._level_boundaries[3]
|
||||
self.assertIs(b.last_raw, b.interval_integral)
|
||||
self.assertIs(b.interval_integral, None)
|
||||
run.step()
|
||||
for name, args, _ in calls:
|
||||
if name == 'BlockCurvedBoundary':
|
||||
self.assertEqual(tuple(args[8:]), (0, 0, 0))
|
||||
self.assertTrue(all((n != 'BlockAccumulateExchange' for n, _, _ in calls)))
|
||||
host.cuda.memset_d8_async.assert_not_called()
|
||||
with self.assertRaisesRegex(RuntimeError, 'disabled'):
|
||||
run.read_raw_exchange(0, level=3)
|
||||
run.close()
|
||||
|
||||
def test_failed_sync_never_publishes_partial_interval(self):
|
||||
run, host, calls, _ = backend()
|
||||
run.step()
|
||||
run._stream.synchronize.side_effect = RuntimeError('device failure')
|
||||
with self.assertRaisesRegex(RuntimeError, 'device failure'):
|
||||
run.step()
|
||||
self.assertFalse(run._exchange_valid)
|
||||
with self.assertRaisesRegex(RuntimeError, 'failed'):
|
||||
run.read_raw_exchange(0, level=3)
|
||||
run._stream.synchronize.side_effect = None
|
||||
run.close()
|
||||
|
||||
def test_dense_observation_rejected_before_device_access(self):
|
||||
host = host_module()
|
||||
ctx = Mock()
|
||||
with self.assertRaisesRegex(ValueError, 'blocked body'):
|
||||
host.FixedPatchStepper(LBMConfig(nx=32, ny=32, viscosity=.1), (4, 4, 28, 28),
|
||||
context=ctx, fine_objects=[body(radius=2.15)], observe_raw_exchange=True)
|
||||
ctx.activate.assert_not_called()
|
||||
|
||||
def test_all_depth_weights(self):
|
||||
for level in (1, 2, 3, 7):
|
||||
self._check_weights_split_and_window(level)
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
unittest.main()
|
||||
@@ -251,14 +251,14 @@ class ConstructorAndSchedule(unittest.TestCase):
|
||||
initial_c=c.current
|
||||
def record(name):
|
||||
return lambda source,*a,**kw: events.append((name,source))
|
||||
run._pairs[0].explode.side_effect=record('explode01')
|
||||
run._pairs[1].explode.side_effect=record('explode12')
|
||||
run._explode.side_effect=record('explode01')
|
||||
run._child.explode.side_effect=record('explode12')
|
||||
for i,level in enumerate(run._levels):
|
||||
level.kernel.side_effect=record(f'step{i}')
|
||||
if bodies:
|
||||
g.boundary.kernel.side_effect=record('bfl2')
|
||||
run._pairs[1].coalesce.side_effect=lambda src,dst,*a,**kw: events.append(('restrict21',dst))
|
||||
run._pairs[0].coalesce.side_effect=lambda src,dst,*a,**kw: events.append(('restrict10',dst))
|
||||
run._child.coalesce.side_effect=lambda src,dst,*a,**kw: events.append(('restrict21',dst))
|
||||
run._coalesce.side_effect=lambda src,dst,*a,**kw: events.append(('restrict10',dst))
|
||||
run._status='ready'
|
||||
expected=[('explode01',initial_c)]
|
||||
for middle in (f.current,f.spare):
|
||||
@@ -278,14 +278,13 @@ class ConstructorAndSchedule(unittest.TestCase):
|
||||
self.assertEqual([l.kernel.call_count for l in run._levels],[1,2,4])
|
||||
# A stale-feed implementation must fail the source-identity guard.
|
||||
feeds=[]
|
||||
advance=run._advance
|
||||
def checked_advance(index, buffer, **kw):
|
||||
if index == 2: # L1 -> L2 feed: must be the L1 hop source
|
||||
self.assertIs(buffer,(f.current,f.spare)[len(feeds)])
|
||||
feeds.append(buffer)
|
||||
advance(index,buffer,**kw)
|
||||
with patch.object(run,'_advance',side_effect=lambda index,buffer,**kw:
|
||||
checked_advance(index,f.current if index == 2 else buffer,**kw)):
|
||||
advance=run._advance_child
|
||||
def checked_advance(buffer, **kw):
|
||||
self.assertIs(buffer,(f.current,f.spare)[len(feeds)])
|
||||
feeds.append(buffer)
|
||||
advance(buffer,**kw)
|
||||
with patch.object(run,'_advance_child',side_effect=lambda buffer,**kw:
|
||||
checked_advance(f.current,**kw)):
|
||||
with self.assertRaises(AssertionError):
|
||||
run.step()
|
||||
# A per-region shared-parent advance must fail the 1:2:4 guard.
|
||||
@@ -311,24 +310,5 @@ class ConstructorAndSchedule(unittest.TestCase):
|
||||
mock_backend(**options)
|
||||
|
||||
|
||||
|
||||
|
||||
class InitialLevels(unittest.TestCase):
|
||||
"""initialize() entry adapter: positional per-level arrays or the legacy names, never both."""
|
||||
|
||||
def test_forms(self):
|
||||
a, b, c = object(), object(), object()
|
||||
f = host._initial_levels
|
||||
self.assertEqual(f((a, b), None, None, None), (a, b))
|
||||
self.assertEqual(f((a, b, c), None, None, None), (a, b, c))
|
||||
self.assertEqual(f((a, b), None, None, c), (a, b, c)) # old initialize(c, f, child_g=g)
|
||||
self.assertEqual(f((), a, b, None), (a, b))
|
||||
self.assertEqual(f((), a, b, c), (a, b, c))
|
||||
for args in (((a,), None, None, c), ((a, b, c), None, None, c), ((a,), a, b, None),
|
||||
((), a, None, None), ((), None, b, c)):
|
||||
with self.subTest(args=args), self.assertRaises(TypeError):
|
||||
f(*args)
|
||||
|
||||
|
||||
if __name__=='__main__':
|
||||
unittest.main()
|
||||
|
||||
@@ -1,5 +1,4 @@
|
||||
"""CPU mocks execute production remap lifecycle; no GPU numerical claim."""
|
||||
from dataclasses import replace
|
||||
import unittest
|
||||
from unittest.mock import patch
|
||||
|
||||
@@ -69,11 +68,8 @@ class RuntimeChildRemap(unittest.TestCase):
|
||||
run._status = 'ready'
|
||||
for bad in ([], [(40,16,48,40)], [(16,16,32,24),(16,27,32,40)]):
|
||||
with self.assertRaises(ValueError): run.remap_child(bad)
|
||||
hosts, pairs = run._boundary_hosts, run.hierarchy.pairs
|
||||
for attribute, value in (('_wall_aware', True),
|
||||
('_boundary_hosts', [*hosts[:2], object()]), # child boundary
|
||||
('_boundary_hosts', [hosts[0], object(), hosts[2]]), # fine boundary
|
||||
('hierarchy', replace(run.hierarchy, pairs=(pairs[0], None)))):
|
||||
for attribute, value in (('_wall_aware', True), ('child_boundary', object()),
|
||||
('fine_boundary', object()), ('child_patch', None)):
|
||||
with patch.object(run, attribute, value):
|
||||
with self.assertRaises(ValueError): run.remap_child(B)
|
||||
self.assertEqual(len(allocations), start)
|
||||
|
||||
@@ -1,264 +0,0 @@
|
||||
"""GPU gates for raw exchange. Run explicitly; requires PyCUDA and CUDA.
|
||||
|
||||
Independent census reads actual staged incoming/outgoing populations, including
|
||||
non-body ACTIVE-boundary flux, then checks postcollision fluid momentum/mass.
|
||||
This is a discrete implementation gate, NOT a physical nozzle-load validation.
|
||||
"""
|
||||
from dataclasses import replace
|
||||
import unittest
|
||||
import numpy as np
|
||||
import pycuda.driver as cuda
|
||||
|
||||
from CelerisLab.body.geometry.circle import CircleGeometry
|
||||
from CelerisLab.body.objects import SimObject
|
||||
from CelerisLab.body.jets import Jet, Tube, Parallel
|
||||
from CelerisLab.config import LBMConfig
|
||||
from CelerisLab.cuda.context import CudaContext
|
||||
from CelerisLab.lbm.refinement_stepper import FixedPatchStepper
|
||||
from CelerisLab.lbm.descriptors import D2Q9_EX, D2Q9_EY, SOLID
|
||||
|
||||
CX, CY = np.asarray(D2Q9_EX), np.asarray(D2Q9_EY)
|
||||
OPP = np.array([np.flatnonzero((CX == -x) & (CY == -y))[0] for x, y in zip(CX, CY)])
|
||||
W = np.array([4/9] + [1/9]*4 + [1/36]*4, np.float32)
|
||||
CFG = LBMConfig(nx=32, ny=24, viscosity=.1)
|
||||
BOUNDS = (8, 6, 24, 18)
|
||||
|
||||
|
||||
def circle():
|
||||
return SimObject(17, geometry=CircleGeometry(16., 12., 2.15), center=(16., 12.), radius=2.15)
|
||||
|
||||
|
||||
def jet():
|
||||
# Asymmetric aperture avoids hiding exchange errors by paired-jet cancellation.
|
||||
return Jet(Tube((16., 12.), 0., 3., .7), Parallel(0.), amplitude=0.)
|
||||
|
||||
|
||||
def state(shape, ux=.013):
|
||||
cu = CX[:, None, None]*ux
|
||||
values = W[:, None, None]*(1+3*cu+4.5*cu*cu-1.5*ux*ux)
|
||||
return np.broadcast_to(values, shape).astype(np.float32).copy()
|
||||
|
||||
|
||||
def download(pointer, shape, dtype=np.float32):
|
||||
result = np.empty(shape, dtype)
|
||||
cuda.memcpy_dtoh(result, pointer)
|
||||
return result
|
||||
|
||||
|
||||
def totals(g, leaf):
|
||||
values = g[:, leaf].astype(np.float64)
|
||||
return np.array([(values*CX[:, None]).sum(), (values*CY[:, None]).sum(), values.sum()])
|
||||
|
||||
|
||||
def census(run, index, buffer):
|
||||
"""Independent incoming/outgoing face census on ACTIVE fluid cells.
|
||||
|
||||
Reads the receiver's actual block window, not a freshly gathered owner image
|
||||
(that would erase BFL/restriction staging). Classifies external flux solely
|
||||
from geometry, without consulting raw observation arrays.
|
||||
"""
|
||||
level = run._levels[index]
|
||||
layout, b = level.layout, level.boundary
|
||||
g = download(buffer, layout.blocked_shape)
|
||||
dense = layout.to_dense(g)
|
||||
leaf = run._leaf[index]
|
||||
ny, nx = leaf.shape
|
||||
solid = (b.host.flags.reshape(ny, nx) & SOLID) != 0
|
||||
y, x = np.nonzero(leaf)
|
||||
owner = layout.owner[y, x].astype(np.intp)
|
||||
cells = int(np.prod(layout.local_shape))
|
||||
block, cell = np.divmod(owner, cells)
|
||||
flat = g.reshape(layout.nblocks, 9, cells)
|
||||
external = np.zeros(3, np.float64)
|
||||
for q, (cx, cy) in enumerate(zip(CX, CY)):
|
||||
sy, sx = y-cy, x-cx
|
||||
incoming = flat[block, q, cell-cy*layout.local_shape[1]-cx].astype(np.float64)
|
||||
entering = ~leaf[sy, sx] & ~solid[sy, sx]
|
||||
external += np.array([cx, cy, 1])*incoming[entering].sum()
|
||||
dy, dx = y+cy, x+cx
|
||||
leaving = ~leaf[dy, dx] & ~solid[dy, dx]
|
||||
outgoing = flat[block, q, cell].astype(np.float64)
|
||||
external -= np.array([cx, cy, 1])*outgoing[leaving].sum()
|
||||
receiver, direction, _, _, rx, ry, _, _ = b.host.curved_host
|
||||
owner = layout.owner.ravel()[receiver].astype(np.intp)
|
||||
block, cell = np.divmod(owner, cells)
|
||||
a = flat[block, direction, cell].astype(np.float64)
|
||||
source = cell+CY[direction]*layout.local_shape[1]+CX[direction]
|
||||
incoming = flat[block, OPP[direction], source].astype(np.float64)
|
||||
jx, jy = CX[direction]*(a+incoming), CY[direction]*(a+incoming)
|
||||
independent = np.column_stack((jx, jy, rx*jy-ry*jx, incoming-a))
|
||||
return totals(dense, leaf), external, independent
|
||||
|
||||
|
||||
class RawExchangeGPU(unittest.TestCase):
|
||||
def test_row_ledger_interval_jet_switch_and_no_population_effect(self):
|
||||
with CudaContext(0) as ctx:
|
||||
for block in ((32, 40), (16, 20), (8, 10)):
|
||||
with FixedPatchStepper(CFG, BOUNDS, context=ctx, fine_block=block,
|
||||
fine_objects=[circle()], jets=[jet()], observe_raw_exchange=True) as measured, \
|
||||
FixedPatchStepper(CFG, BOUNDS, context=ctx, fine_block=block,
|
||||
fine_objects=[circle()], jets=[jet()]) as plain:
|
||||
arrays = [state(l.shape) for l in measured._levels]
|
||||
measured.initialize(*arrays)
|
||||
plain.initialize(*arrays)
|
||||
b = measured._level_boundaries[1]
|
||||
count = len(b.tables['cell'])
|
||||
kernel = measured._levels[1].kernel
|
||||
samples = []
|
||||
|
||||
def audited_step(src, dst, *args, **kw):
|
||||
measured._stream.synchronize()
|
||||
before, external, independent = census(measured, 1, src)
|
||||
raw = download(b.last_raw, (count, 4))
|
||||
np.testing.assert_allclose(raw, independent, rtol=2e-6, atol=2e-7)
|
||||
kernel(src, dst, *args, **kw)
|
||||
measured._stream.synchronize()
|
||||
after = measured._levels[1].layout.to_dense(
|
||||
download(dst, measured._levels[1].layout.blocked_shape))
|
||||
change = totals(after, measured._leaf[1])-before
|
||||
expected = external + np.array([-raw[:, 0].astype(float).sum(),
|
||||
-raw[:, 1].astype(float).sum(),
|
||||
raw[:, 3].astype(float).sum()])
|
||||
np.testing.assert_allclose(change, expected, rtol=2e-4, atol=8e-5)
|
||||
samples.append(raw.copy())
|
||||
measured._levels[1].kernel = audited_step
|
||||
for speed in (0., .002, -.001, 0.):
|
||||
measured.set_jet_amplitude(0, speed)
|
||||
plain.set_jet_amplitude(0, speed)
|
||||
samples.clear()
|
||||
measured.step()
|
||||
plain.step()
|
||||
for a, e in zip(measured.download(), plain.download()):
|
||||
np.testing.assert_array_equal(a, e)
|
||||
integral = download(b.interval_integral, (count, 4), np.float64)
|
||||
expected = sum(s.astype(float)*[.25, .25, .125, .25] for s in samples)
|
||||
np.testing.assert_array_equal(integral, expected)
|
||||
last = download(b.last_raw, (count, 4))
|
||||
np.testing.assert_array_equal(last, samples[-1])
|
||||
obs = measured.read_raw_exchange(0)
|
||||
assert obs['sample_count'] == 2
|
||||
assert b.tables['aperture'].any()
|
||||
measured.step(0)
|
||||
np.testing.assert_array_equal(measured.read_raw_exchange(0)['total']['impulse'],
|
||||
obs['total']['impulse'])
|
||||
|
||||
def test_dense_fine_storage_reference(self):
|
||||
"""Same fine resolution/initial data with dense storage: a numerical
|
||||
reference for cut-link rows, not a uniform-whole-domain convergence study."""
|
||||
with CudaContext(0) as ctx:
|
||||
with FixedPatchStepper(CFG, BOUNDS, context=ctx, fine_objects=[circle()], jets=[jet()]) as dense, \
|
||||
FixedPatchStepper(CFG, BOUNDS, context=ctx, fine_block=(16, 20),
|
||||
fine_objects=[circle()], jets=[jet()], observe_raw_exchange=True) as blocked:
|
||||
arrays = [state(l.shape) for l in blocked._levels]
|
||||
dense.initialize(*arrays)
|
||||
blocked.initialize(*arrays)
|
||||
for speed in (0., .002, 0.):
|
||||
dense.set_jet_amplitude(0, speed)
|
||||
blocked.set_jet_amplitude(0, speed)
|
||||
dense.step(); blocked.step()
|
||||
for a, e in zip(blocked.download(), dense.download()):
|
||||
np.testing.assert_array_equal(a, e)
|
||||
count = dense._fine_curved.count
|
||||
expected = download(dense._fine_curved.link_exchange_gpu, (count, 3))
|
||||
actual = download(blocked._level_boundaries[1].last_raw, (count, 4))
|
||||
np.testing.assert_array_equal(actual[:, :3], expected)
|
||||
|
||||
def test_depth_four_interval_storage_and_sampling(self):
|
||||
cfg = LBMConfig(nx=48, ny=48, viscosity=.02)
|
||||
obj = SimObject(0, geometry=CircleGeometry(16., 16., .5), center=(16., 16.), radius=.5)
|
||||
drive = Jet(Tube((16., 16.), 0., .75, .2), Parallel(0.), amplitude=0.)
|
||||
settings = dict(refinements=[[(6, 6, 26, 26)], [(20, 20, 44, 44)], [(56, 56, 72, 72)]],
|
||||
level_blocks=[(4, 4)]*4, level_objects=[(), (), (), (obj,)],
|
||||
level_jets=[(), (), (), (drive,)])
|
||||
with CudaContext(0) as ctx:
|
||||
with FixedPatchStepper(cfg, context=ctx, observe_raw_exchange=True, **settings) as run, \
|
||||
FixedPatchStepper(cfg, context=ctx, **settings) as plain:
|
||||
arrays = [state(l.shape, ux=0.) for l in run._levels]
|
||||
run.initialize(*arrays); plain.initialize(*arrays)
|
||||
b = run._level_boundaries[3]
|
||||
count = len(b.tables['cell'])
|
||||
kernel = run._levels[3].kernel
|
||||
samples, independent_samples = [], []
|
||||
def audited_step(src, dst, *args, **kw):
|
||||
run._stream.synchronize()
|
||||
before, external, independent = census(run, 3, src)
|
||||
raw = download(b.last_raw, (count, 4))
|
||||
np.testing.assert_allclose(raw, independent, rtol=2e-6, atol=2e-7)
|
||||
kernel(src, dst, *args, **kw)
|
||||
run._stream.synchronize()
|
||||
after = run._levels[3].layout.to_dense(
|
||||
download(dst, run._levels[3].layout.blocked_shape))
|
||||
change = totals(after, run._leaf[3])-before
|
||||
expected = external + np.array([-independent[:, 0].sum(),
|
||||
-independent[:, 1].sum(), independent[:, 3].sum()])
|
||||
np.testing.assert_allclose(change, expected, rtol=2e-4, atol=8e-5)
|
||||
samples.append(raw)
|
||||
independent_samples.append(independent)
|
||||
run._levels[3].kernel = audited_step
|
||||
for speed in (0., .002, 0.):
|
||||
run.set_jet_amplitude(0, speed, level=3)
|
||||
plain.set_jet_amplitude(0, speed, level=3)
|
||||
samples.clear(); independent_samples.clear()
|
||||
run.step(); plain.step()
|
||||
self.assertEqual(len(samples), 8)
|
||||
weights = np.array([1/64, 1/64, 1/512, 1/64])
|
||||
integral = download(b.interval_integral, (count, 4), np.float64)
|
||||
np.testing.assert_array_equal(integral, sum(v.astype(float)*weights for v in samples))
|
||||
np.testing.assert_allclose(integral, sum(v*weights for v in independent_samples),
|
||||
rtol=2e-6, atol=2e-7)
|
||||
np.testing.assert_array_equal(download(b.last_raw, (count, 4)), samples[-1])
|
||||
for actual, expected in zip(run.download(), plain.download()):
|
||||
np.testing.assert_array_equal(actual, expected)
|
||||
result = run.read_raw_exchange(0, level=3)
|
||||
self.assertEqual(result['sample_count'], 8)
|
||||
self.assertEqual(result['time_start'], run.step_count-1)
|
||||
if run.step_count == 1:
|
||||
np.testing.assert_allclose(result['total']['force'], 0, atol=2e-7)
|
||||
self.assertAlmostEqual(result['total']['torque'], 0, delta=2e-7)
|
||||
|
||||
def test_uniform_whole_domain_fine_short_time_reference(self):
|
||||
"""Same root geometry, nu and Jet speed on a whole-domain fine mesh.
|
||||
Eight local hops are a short-time equivalence check, not convergence.
|
||||
The test adapter only enables the standalone kernel's optional inputs.
|
||||
"""
|
||||
from CelerisLab.lbm.blocks import BlockLayout2D
|
||||
from CelerisLab.lbm.block_stepper import PeriodicBlockStepper
|
||||
fine_cfg = replace(CFG, nx=64, ny=48, viscosity=.2)
|
||||
obj = SimObject(0, geometry=CircleGeometry(31.5, 23.5, 4.3),
|
||||
center=(31.5, 23.5), radius=4.3)
|
||||
layout = BlockLayout2D.build(48, 64, 8, 8)
|
||||
with CudaContext(0) as ctx:
|
||||
with PeriodicBlockStepper(fine_cfg, layout, context=ctx, objects=[obj]) as uniform, \
|
||||
FixedPatchStepper(CFG, BOUNDS, context=ctx, fine_block=(16, 20),
|
||||
fine_objects=[circle()], jets=[jet()], observe_raw_exchange=True) as amr:
|
||||
amr.initialize(*(state(l.shape) for l in amr._levels))
|
||||
uniform.initialize(state((9, 48, 64)))
|
||||
b = amr._level_boundaries[1]
|
||||
count = len(b.tables['cell'])
|
||||
links = obj.geometry.build_cut_links(64, 48)
|
||||
np.testing.assert_array_equal([v.dir for v in links], b.tables['dir'])
|
||||
np.testing.assert_allclose([v.q for v in links], b.tables['q'], atol=2e-7, rtol=0)
|
||||
# Same fine-grid force arms after translation to the same center.
|
||||
rx = uniform._upload(b.tables['rx'])
|
||||
ry = uniform._upload(b.tables['ry'])
|
||||
velocity = uniform._upload(b.tables['surface_velocity'])
|
||||
raw = uniform._allocate(count*16)
|
||||
original = uniform._curved
|
||||
samples = []
|
||||
def measured_boundary(*args, **kw):
|
||||
original(*args[:7], velocity, rx, ry, raw, **kw)
|
||||
uniform._stream.synchronize()
|
||||
samples.append(download(raw, (count, 4)))
|
||||
uniform._curved = measured_boundary
|
||||
for speed in (0., .002, -.001, 0.):
|
||||
amr.set_jet_amplitude(0, speed)
|
||||
cuda.memcpy_htod(velocity, b.tables['surface_velocity'])
|
||||
samples.clear()
|
||||
amr.step(); uniform.step(2)
|
||||
expected = sum(v.astype(float)*[.25, .25, .125, .25] for v in samples)
|
||||
actual = download(b.interval_integral, (count, 4), np.float64)
|
||||
np.testing.assert_allclose(actual, expected, rtol=1e-5, atol=2e-6)
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
@@ -194,7 +194,7 @@ class CoupledGPU(unittest.TestCase):
|
||||
f"max_corner_packet={max_corner:.9g}; ownership/counts PASS", flush=True)
|
||||
|
||||
# Fail before any additional time step; no stale synchronized export.
|
||||
with patch.object(backend._pairs[0], "explode", side_effect=RuntimeError("injected launch failure")):
|
||||
with patch.object(backend, "_explode", side_effect=RuntimeError("injected launch failure")):
|
||||
with self.assertRaisesRegex(RuntimeError, "injected"):
|
||||
backend.step()
|
||||
with self.assertRaisesRegex(RuntimeError, "failed"):
|
||||
|
||||
Reference in New Issue
Block a user