476 lines
21 KiB
Python
476 lines
21 KiB
Python
"""Isolated NumPy CPU reference: one static 2:1 rectangular patch, volumetric transfer.
|
|
|
|
API version CPUREF-0.3 (ClaudeWorker, 2026-10-13). NOT EXECUTED BY THE AUTHOR:
|
|
written without a shell; Grok/Horsea run it. A test-reference CANDIDATE, not a
|
|
certified truth and not a project implementation. v0.2 (Drive
|
|
1gfYrjdpzLxQ8LCZAH_pMMQRLSSuvFMIf, SHA256 b6adbad7...) is kept unchanged.
|
|
|
|
Scope of evidence
|
|
Checked for conservation and bookkeeping only. Interface accuracy/stress
|
|
transfer is unvalidated (for every tau, tau_coarse = 1 included).
|
|
|
|
Changes from v0.2 (jointly agreed small patch; numerical algorithm unchanged):
|
|
* explosion_count removed (it was a constant); explosion_writes is the evidence.
|
|
* Corner-pass links are derived from the coalescence map (rows whose coarse
|
|
source cell is not covered): all 8, each with a unique key "<corner>_q<i>",
|
|
e.g. "lower_left_q8", "lower_left_q7". corner_packets is keyed the same way.
|
|
* Diagnostic renames: fine_valid_after -> hop2_valid; budget closure_error ->
|
|
bookkeeping_check.
|
|
* Per-packet path records (coalesce_paths, origin_g copy) moved to the tests;
|
|
cheap in-code guards kept (invalid packet, path through an outer ghost).
|
|
* Coalescence moved into a _coalesce helper.
|
|
* state.tick must be a non-negative even integer (synchronised phase); other
|
|
values are rejected instead of int()-truncated.
|
|
|
|
Contract: AMR_Shared_CPU_API_2026-10-13.md
|
|
SHA256 ee8010e763d38937b95c34dadaf97415cdb8b6c23e3c79a9f1f32bd5d4e112f6
|
|
plus root amendment: make_spec(*, collision='identity'|'srt'), stored in Spec.
|
|
Snapshot: 20-file snapshot
|
|
SHA256 e45df47d0fabde7199d345ada0729bca9ee64a833f433cdd776a4774b6af0673
|
|
Dependencies: NumPy and the standard library only. No CelerisLab/CUDA imports.
|
|
|
|
Conventions
|
|
D2Q9 paired ordering cx=(0,1,-1,0,0,1,-1,1,-1), cy=(0,0,0,1,-1,1,-1,-1,1).
|
|
Arrays are Q-major [q, y, x], float64, unshifted physical populations.
|
|
Coarse cell (x, y) has centre (x+.5, y+.5) on the doubly periodic [0,16)^2.
|
|
The patch covers coarse cells [4,12)^2 (covered). Fine index a = 0..23 (halo 4)
|
|
has centre 4 + (a-4+.5)/2 and coarse parent 4 + (a-4)//2; active fine is
|
|
[:, 4:20, 4:20]. Fine cell volume is 1/4 of a coarse cell.
|
|
Invalid entries (covered coarse, fine ghosts at synchronisation) hold NaN.
|
|
|
|
Phase rule
|
|
advance_coarse_once maps arrival f(t) to arrival f(t+1). Production stores
|
|
post-collision g. Compare C(f_ref) with g_prod at a matched tick; never
|
|
compare f with g directly and never apply C to a stored g again.
|
|
|
|
Synchronised phase "arrival_f". One advance_coarse_once (Schornbaum & Ruede
|
|
2016, Sec. 4 order):
|
|
1 collide coarse active cells g_c = C_c(f_c(t)) tau_c
|
|
2 collide fine active cells g_f = C_f(f_f(t)) tau_f
|
|
3 explode: every fine ghost cell (inner and outer) receives g_c of its coarse
|
|
parent in all directions; once per coarse step, no rescaling
|
|
4 fine stream, hop 1; targets are active cells and inner ghosts (inner ghosts
|
|
stream but never collide)
|
|
5 collide fine active cells only; outer ghosts are poisoned (NaN, invalid)
|
|
6 fine stream, hop 2; active cells now hold f_f(t+1)
|
|
7 coarse periodic stream; incoming populations of covered cells are invalid
|
|
8 coalesce: for each non-covered coarse cell K and direction i whose link
|
|
midpoint (centre(K) - c_i/2) lies in the CLOSED patch, replace the coarse
|
|
incoming population by the mean of K's four inner-ghost children's slot i
|
|
(includes the 8 corner-pass links, e.g. (4,3)->(3,4), i=8).
|
|
|
|
Assumptions (stated, not validated by rest-state tests)
|
|
* Ghosts hold g_c(t) for the whole coarse step and inner ghosts never collide,
|
|
so ghost-sourced inflow into the first fine layer is an uncollided coarse
|
|
packet in both hops.
|
|
* The closed corner rule matters only through the mid-step fine collision:
|
|
with identity collision the open-corner variant is also conservative, so
|
|
identity-only tests cannot distinguish the two rules; use SRT with
|
|
non-equilibrium states.
|
|
* Budget: bookkeeping_check is self-consistency of this file's flux accounting
|
|
(it subtracts interface_residual), NOT conservation evidence; test
|
|
leaf_before/leaf_after and interface_residual independently.
|
|
"""
|
|
from __future__ import annotations
|
|
|
|
from dataclasses import dataclass
|
|
from typing import NamedTuple
|
|
|
|
import numpy as np
|
|
|
|
API_VERSION = "CPUREF-0.3"
|
|
CONTRACT_SHA256 = "ee8010e763d38937b95c34dadaf97415cdb8b6c23e3c79a9f1f32bd5d4e112f6"
|
|
SNAPSHOT_SHA256 = "e45df47d0fabde7199d345ada0729bca9ee64a833f433cdd776a4774b6af0673"
|
|
|
|
NQ = 9
|
|
CX = np.array([0, 1, -1, 0, 0, 1, -1, 1, -1], dtype=np.int64)
|
|
CY = np.array([0, 0, 0, 1, -1, 1, -1, -1, 1], dtype=np.int64)
|
|
OPP = np.array([0, 2, 1, 4, 3, 6, 5, 8, 7], dtype=np.int64)
|
|
W = np.array([4 / 9] + [1 / 9] * 4 + [1 / 36] * 4, dtype=np.float64)
|
|
FINE_VOLUME = 0.25
|
|
# Moment weights per direction: (mass, x-momentum, y-momentum).
|
|
MOMENTS = np.stack([np.ones(NQ), CX.astype(np.float64), CY.astype(np.float64)], axis=1)
|
|
|
|
ACTIVE, INNER_GHOST, OUTER_GHOST = 0, 1, 2
|
|
|
|
|
|
@dataclass(frozen=True)
|
|
class Spec:
|
|
patch_id: str = "P0"
|
|
n_coarse: int = 16
|
|
patch_lo: int = 4
|
|
patch_hi: int = 12
|
|
halo: int = 4
|
|
dt_coarse: float = 1.0
|
|
dt_fine: float = 0.5
|
|
tau_coarse: float = 0.8
|
|
tau_fine: float = 1.1
|
|
collision: str = "identity"
|
|
|
|
@property
|
|
def n_fine(self) -> int:
|
|
return 2 * (self.patch_hi - self.patch_lo) + 2 * self.halo
|
|
|
|
|
|
@dataclass(frozen=True)
|
|
class State:
|
|
coarse_f: np.ndarray
|
|
fine_f: np.ndarray
|
|
phase: str = "arrival_f"
|
|
tick: int = 0
|
|
|
|
|
|
class CornerLink(NamedTuple):
|
|
key: str # "<corner>_q<i>", unique over the 8 corner-pass links
|
|
i: int
|
|
source_xy: tuple # uncovered coarse source cell (x, y)
|
|
target_xy: tuple # coarse target cell K (x, y)
|
|
|
|
|
|
@dataclass(frozen=True)
|
|
class Maps:
|
|
coarse_active: np.ndarray # (n, n) bool [y, x]
|
|
coarse_covered: np.ndarray # (n, n) bool
|
|
fine_active: np.ndarray # (nf, nf) bool
|
|
fine_inner_ghost: np.ndarray # (nf, nf) bool
|
|
fine_outer_ghost: np.ndarray # (nf, nf) bool
|
|
fine_class: np.ndarray # (nf, nf) int: ACTIVE / INNER_GHOST / OUTER_GHOST
|
|
fine_parent: np.ndarray # (nf,) coarse parent index (same for x and y)
|
|
fine_centre: np.ndarray # (nf,) physical centre (same for x and y)
|
|
coalesce: np.ndarray # (M, 3) rows (i, ky, kx)
|
|
corner_links: tuple # 8 x CornerLink, sorted by key
|
|
|
|
|
|
def make_spec(*, collision: str = "identity") -> Spec:
|
|
if collision not in ("identity", "srt"):
|
|
raise ValueError("collision must be 'identity' or 'srt'")
|
|
spec = Spec(collision=collision)
|
|
if not np.isclose(spec.tau_fine - 0.5, 2.0 * (spec.tau_coarse - 0.5)):
|
|
raise ValueError("tau values violate acoustic 2:1 scaling")
|
|
if not np.isclose(spec.dt_fine, 0.5 * spec.dt_coarse):
|
|
raise ValueError("dt values violate 2:1 time refinement")
|
|
return spec
|
|
|
|
|
|
def _readonly(array) -> np.ndarray:
|
|
array = np.array(array)
|
|
array.setflags(write=False)
|
|
return array
|
|
|
|
|
|
def _children_start(spec: Spec, k: int) -> int:
|
|
"""First fine index of the two children of coarse index k."""
|
|
return spec.halo + 2 * (k - spec.patch_lo)
|
|
|
|
|
|
def build_maps(spec: Spec) -> Maps:
|
|
n, lo, hi, nf = spec.n_coarse, spec.patch_lo, spec.patch_hi, spec.n_fine
|
|
idx = np.arange(n)
|
|
inside = (idx >= lo) & (idx < hi)
|
|
covered = inside[:, None] & inside[None, :]
|
|
|
|
a = np.arange(nf)
|
|
parent = lo + (a - spec.halo) // 2
|
|
centre = lo + (a - spec.halo + 0.5) / 2
|
|
active_1d = (a >= spec.halo) & (a < nf - spec.halo)
|
|
ring_1d = (parent >= lo - 1) & (parent <= hi) # children of patch or first ring
|
|
fine_active = active_1d[:, None] & active_1d[None, :]
|
|
near = ring_1d[:, None] & ring_1d[None, :]
|
|
inner = near & ~fine_active
|
|
outer = ~near
|
|
fine_class = np.where(fine_active, ACTIVE, np.where(inner, INNER_GHOST, OUTER_GHOST))
|
|
|
|
rows = []
|
|
for ky in range(n):
|
|
for kx in range(n):
|
|
if covered[ky, kx]:
|
|
continue
|
|
for i in range(1, NQ):
|
|
mx = kx + 0.5 - CX[i] / 2.0
|
|
my = ky + 0.5 - CY[i] / 2.0
|
|
if lo <= mx <= hi and lo <= my <= hi:
|
|
rows.append((i, ky, kx))
|
|
coalesce = np.array(rows, dtype=np.int64).reshape(-1, 3)
|
|
for i, ky, kx in coalesce:
|
|
y0, x0 = _children_start(spec, ky), _children_start(spec, kx)
|
|
if not inner[y0:y0 + 2, x0:x0 + 2].all():
|
|
raise AssertionError("coalescence target children must be inner ghosts")
|
|
|
|
# Corner-pass links: coalesced rows whose coarse source cell is not covered.
|
|
links = []
|
|
for i, ky, kx in coalesce.tolist():
|
|
sx, sy = kx - int(CX[i]), ky - int(CY[i])
|
|
if covered[sy, sx]:
|
|
continue
|
|
mx, my = kx + 0.5 - CX[i] / 2.0, ky + 0.5 - CY[i] / 2.0
|
|
if mx not in (lo, hi) or my not in (lo, hi):
|
|
raise AssertionError("corner-pass link must cross a patch corner point")
|
|
name = ("lower" if my == lo else "upper") + ("_left" if mx == lo else "_right")
|
|
links.append(CornerLink(f"{name}_q{i}", i, (sx, sy), (kx, ky)))
|
|
links.sort(key=lambda link: link.key)
|
|
if len(links) != 8 or len({link.key for link in links}) != 8:
|
|
raise AssertionError("expected 8 uniquely keyed corner-pass links")
|
|
|
|
return Maps(
|
|
coarse_active=_readonly(~covered), coarse_covered=_readonly(covered),
|
|
fine_active=_readonly(fine_active), fine_inner_ghost=_readonly(inner),
|
|
fine_outer_ghost=_readonly(outer), fine_class=_readonly(fine_class),
|
|
fine_parent=_readonly(parent), fine_centre=_readonly(centre),
|
|
coalesce=_readonly(coalesce), corner_links=tuple(links))
|
|
|
|
|
|
def _equilibrium(rho, ux, uy) -> np.ndarray:
|
|
"""D2Q9 second-order equilibrium; inputs are 1-D arrays of equal length."""
|
|
cu = CX[:, None] * ux[None, :] + CY[:, None] * uy[None, :]
|
|
usq = ux * ux + uy * uy
|
|
return W[:, None] * rho[None, :] * (1.0 + 3.0 * cu + 4.5 * cu * cu - 1.5 * usq[None, :])
|
|
|
|
|
|
def make_state(spec: Spec, rho: float = 1.0, velocity=(0.0, 0.0)) -> State:
|
|
maps = build_maps(spec)
|
|
ux, uy = velocity
|
|
feq = _equilibrium(np.array([rho], float), np.array([ux], float), np.array([uy], float))[:, 0]
|
|
n, nf = spec.n_coarse, spec.n_fine
|
|
coarse = np.full((NQ, n, n), np.nan)
|
|
fine = np.full((NQ, nf, nf), np.nan)
|
|
coarse[:, maps.coarse_active] = feq[:, None]
|
|
fine[:, maps.fine_active] = feq[:, None]
|
|
return State(coarse_f=coarse, fine_f=fine, phase="arrival_f", tick=0)
|
|
|
|
|
|
def _collide(spec: Spec, f: np.ndarray, valid: np.ndarray, tau: float,
|
|
mask: np.ndarray, counts: np.ndarray) -> np.ndarray:
|
|
"""Collision on masked cells (counted here); other cells copied unchanged."""
|
|
if not valid[:, mask].all():
|
|
raise RuntimeError("collision touches a cell with an invalid direction")
|
|
g = f.copy()
|
|
if spec.collision == "srt":
|
|
fm = f[:, mask]
|
|
rho = fm.sum(axis=0)
|
|
ux = (CX[:, None] * fm).sum(axis=0) / rho
|
|
uy = (CY[:, None] * fm).sum(axis=0) / rho
|
|
g[:, mask] = fm - (fm - _equilibrium(rho, ux, uy)) / tau
|
|
counts[mask] += 1
|
|
return g
|
|
|
|
|
|
def _explode(gc, coarse_valid, gf, fine_valid, maps: Maps, writes: np.ndarray) -> None:
|
|
"""Copy each ghost's coarse parent into it (all directions); counts writes."""
|
|
ghosts = maps.fine_inner_ghost | maps.fine_outer_ghost
|
|
gy, gx = np.nonzero(ghosts)
|
|
py, px = maps.fine_parent[gy], maps.fine_parent[gx]
|
|
if not coarse_valid[:, py, px].all():
|
|
raise RuntimeError("explosion source is not a valid coarse cell")
|
|
gf[:, gy, gx] = gc[:, py, px]
|
|
fine_valid[:, gy, gx] = True
|
|
np.add.at(writes, (gy, gx), 1)
|
|
|
|
|
|
def _stream_periodic(g: np.ndarray) -> np.ndarray:
|
|
out = np.empty_like(g)
|
|
for i in range(NQ):
|
|
out[i] = np.roll(g[i], (CY[i], CX[i]), axis=(0, 1)) # out[y,x] = g[y-cy, x-cx]
|
|
return out
|
|
|
|
|
|
def _shift(a: np.ndarray, dy: int, dx: int, fill) -> np.ndarray:
|
|
"""out[y, x] = a[y-dy, x-dx]; `fill` where the source lies outside the array."""
|
|
out = np.full(a.shape, fill, dtype=a.dtype)
|
|
ny, nx = a.shape
|
|
out[max(dy, 0):ny + min(dy, 0), max(dx, 0):nx + min(dx, 0)] = \
|
|
a[max(-dy, 0):ny + min(-dy, 0), max(-dx, 0):nx + min(-dx, 0)]
|
|
return out
|
|
|
|
|
|
def _fine_hop(g: np.ndarray, valid: np.ndarray, maps: Maps):
|
|
"""One fine pull-stream to active cells and inner ghosts, with exchange records."""
|
|
targets = maps.fine_active | maps.fine_inner_ghost
|
|
src_class = maps.fine_class.astype(np.float64)
|
|
out = np.full_like(g, np.nan)
|
|
out_valid = np.zeros_like(valid)
|
|
inflow = np.zeros(3) # into active cells from ghost sources (fine-volume units)
|
|
outflow = np.zeros(3) # from active sources into inner-ghost targets
|
|
sources = {"active": 0, "inner_ghost": 0, "outer_ghost": 0, "outside": 0}
|
|
with np.errstate(invalid="ignore"):
|
|
for i in range(NQ):
|
|
dy, dx = int(CY[i]), int(CX[i])
|
|
moved = _shift(g[i], dy, dx, np.nan)
|
|
cls = _shift(src_class, dy, dx, np.nan)
|
|
out[i] = moved
|
|
out_valid[i] = _shift(valid[i], dy, dx, False) & targets
|
|
from_ghost = maps.fine_active & (cls >= INNER_GHOST)
|
|
to_ghost = maps.fine_inner_ghost & (cls == ACTIVE)
|
|
inflow += MOMENTS[i] * moved[from_ghost].sum()
|
|
outflow += MOMENTS[i] * moved[to_ghost].sum()
|
|
act_cls = cls[maps.fine_active]
|
|
sources["active"] += int(np.sum(act_cls == ACTIVE))
|
|
sources["inner_ghost"] += int(np.sum(act_cls == INNER_GHOST))
|
|
sources["outer_ghost"] += int(np.sum(act_cls == OUTER_GHOST))
|
|
sources["outside"] += int(np.sum(np.isnan(act_cls)))
|
|
out[:, ~targets] = np.nan
|
|
out[~out_valid] = np.nan
|
|
return out, out_valid, FINE_VOLUME * inflow, FINE_VOLUME * outflow, sources
|
|
|
|
|
|
def _coalesce(spec: Spec, maps: Maps, f2, hop2_valid, fc_new, new_valid, owner):
|
|
"""Overwrite each listed coarse incoming slot with the mean of its four
|
|
inner-ghost children (fc_new, new_valid, owner updated in place).
|
|
Returns (replaced, coalesced, writes, reads)."""
|
|
n, nf = spec.n_coarse, spec.n_fine
|
|
replaced = np.zeros(3)
|
|
coalesced = np.zeros(3)
|
|
writes = np.zeros((NQ, n, n), dtype=np.int64)
|
|
reads = np.zeros((NQ, nf, nf), dtype=np.int64)
|
|
for i, ky, kx in maps.coalesce.tolist():
|
|
y0, x0 = _children_start(spec, ky), _children_start(spec, kx)
|
|
kids = (slice(y0, y0 + 2), slice(x0, x0 + 2))
|
|
if not hop2_valid[i][kids].all():
|
|
raise RuntimeError(f"invalid coalescence packet (i={i}, K=({kx},{ky}))")
|
|
ym, xm = y0 - int(CY[i]), x0 - int(CX[i]) # hop-2 sources (intermediates)
|
|
if (maps.fine_class[ym:ym + 2, xm:xm + 2] == OUTER_GHOST).any():
|
|
raise RuntimeError("coalescence packet passed through an outer ghost")
|
|
value = f2[i][kids].mean()
|
|
if new_valid[i, ky, kx]:
|
|
replaced += MOMENTS[i] * fc_new[i, ky, kx]
|
|
fc_new[i, ky, kx] = value
|
|
new_valid[i, ky, kx] = True
|
|
owner[i, ky, kx] = 2
|
|
coalesced += MOMENTS[i] * value
|
|
writes[i, ky, kx] += 1
|
|
reads[i][kids] += 1
|
|
return replaced, coalesced, writes, reads
|
|
|
|
|
|
def _leaf_totals(coarse_f, fine_f, maps: Maps) -> np.ndarray:
|
|
mc = coarse_f[:, maps.coarse_active]
|
|
mf = fine_f[:, maps.fine_active]
|
|
return MOMENTS.T @ mc.sum(axis=1) + FINE_VOLUME * (MOMENTS.T @ mf.sum(axis=1))
|
|
|
|
|
|
def _moments_of_difference(after, before, mask) -> np.ndarray:
|
|
return MOMENTS.T @ (after[:, mask] - before[:, mask]).sum(axis=1)
|
|
|
|
|
|
def advance_coarse_once(spec: Spec, maps: Maps, state: State):
|
|
n, nf = spec.n_coarse, spec.n_fine
|
|
if state.phase != "arrival_f":
|
|
raise ValueError("state.phase must be 'arrival_f'")
|
|
tick = state.tick
|
|
if isinstance(tick, bool) or not isinstance(tick, (int, np.integer)) or tick < 0 or tick % 2:
|
|
raise ValueError("state.tick must be a non-negative even integer (synchronised phase)")
|
|
tick = int(tick)
|
|
if state.coarse_f.shape != (NQ, n, n) or state.fine_f.shape != (NQ, nf, nf):
|
|
raise ValueError("state arrays have the wrong shape")
|
|
fc = np.array(state.coarse_f, dtype=np.float64)
|
|
ff = np.array(state.fine_f, dtype=np.float64)
|
|
if not (np.isfinite(fc[:, maps.coarse_active]).all() and np.isfinite(ff[:, maps.fine_active]).all()):
|
|
raise ValueError("active populations must be finite")
|
|
fc[:, ~maps.coarse_active] = np.nan # invalid input entries are never trusted
|
|
ff[:, ~maps.fine_active] = np.nan
|
|
coarse_valid = np.broadcast_to(maps.coarse_active, (NQ, n, n)).copy()
|
|
fine_valid = np.broadcast_to(maps.fine_active, (NQ, nf, nf)).copy()
|
|
leaf_before = _leaf_totals(fc, ff, maps)
|
|
|
|
count_c = np.zeros((n, n), dtype=np.int64)
|
|
count_f = np.zeros((nf, nf), dtype=np.int64)
|
|
explosion_writes = np.zeros((nf, nf), dtype=np.int64)
|
|
|
|
# 1-2: collisions at t
|
|
gc = _collide(spec, fc, coarse_valid, spec.tau_coarse, maps.coarse_active, count_c)
|
|
gf = _collide(spec, ff, fine_valid, spec.tau_fine, maps.fine_active, count_f)
|
|
coarse_collision = _moments_of_difference(gc, fc, maps.coarse_active)
|
|
fine_collision = FINE_VOLUME * _moments_of_difference(gf, ff, maps.fine_active)
|
|
|
|
# 3: explosion, once per coarse step
|
|
_explode(gc, coarse_valid, gf, fine_valid, maps, explosion_writes)
|
|
|
|
# 4: hop 1
|
|
f1, valid1, in1, out1, src1 = _fine_hop(gf, fine_valid, maps)
|
|
# 5: mid-step fine collision, active only; outer ghosts poisoned
|
|
g1 = _collide(spec, f1, valid1, spec.tau_fine, maps.fine_active, count_f)
|
|
fine_collision += FINE_VOLUME * _moments_of_difference(g1, f1, maps.fine_active)
|
|
valid1 = valid1.copy()
|
|
g1[:, maps.fine_outer_ghost] = np.nan
|
|
valid1[:, maps.fine_outer_ghost] = False
|
|
# 6: hop 2
|
|
f2, hop2_valid, in2, out2, src2 = _fine_hop(g1, valid1, maps)
|
|
|
|
# 7: coarse stream
|
|
streamed = _stream_periodic(gc)
|
|
streamed_valid = np.stack([np.roll(coarse_valid[i], (CY[i], CX[i]), axis=(0, 1)) for i in range(NQ)])
|
|
to_covered = np.zeros(3)
|
|
for i in range(NQ):
|
|
to_covered += MOMENTS[i] * streamed[i][maps.coarse_covered & streamed_valid[i]].sum()
|
|
fc_new = streamed.copy()
|
|
new_valid = streamed_valid & maps.coarse_active[None, :, :]
|
|
fc_new[~new_valid] = np.nan
|
|
owner = np.zeros((NQ, n, n), dtype=np.int64)
|
|
owner[new_valid] = 1
|
|
|
|
# 8: coalescence of incoming directions from the fine region
|
|
replaced, coalesced, coalesce_writes, coalesce_reads = _coalesce(
|
|
spec, maps, f2, hop2_valid, fc_new, new_valid, owner)
|
|
|
|
if not new_valid[:, maps.coarse_active].all() or not np.isfinite(fc_new[:, maps.coarse_active]).all():
|
|
raise RuntimeError("coarse incoming population without a valid owner")
|
|
if not hop2_valid[:, maps.fine_active].all() or not np.isfinite(f2[:, maps.fine_active]).all():
|
|
raise RuntimeError("fine active population received an invalid value")
|
|
ff_new = np.full((NQ, nf, nf), np.nan)
|
|
ff_new[:, maps.fine_active] = f2[:, maps.fine_active]
|
|
leaf_after = _leaf_totals(fc_new, ff_new, maps)
|
|
|
|
fine_inflow, fine_outflow = in1 + in2, out1 + out2
|
|
interface = (coalesced - to_covered - replaced) + (fine_inflow - fine_outflow)
|
|
budget = {
|
|
"coarse_collision": coarse_collision,
|
|
"fine_collision": fine_collision,
|
|
"coarse_to_covered": to_covered,
|
|
"coarse_replaced": replaced,
|
|
"coalesced": coalesced,
|
|
"fine_inflow": fine_inflow,
|
|
"fine_outflow": fine_outflow,
|
|
"interface_residual": interface,
|
|
# Self-consistency of this file's accounting only; defined using interface_residual.
|
|
"bookkeeping_check": (leaf_after - leaf_before) - (coarse_collision + fine_collision + interface),
|
|
}
|
|
|
|
corner_packets = {}
|
|
for link in maps.corner_links:
|
|
kx, ky = link.target_xy
|
|
y0, x0 = _children_start(spec, ky), _children_start(spec, kx)
|
|
kids = [(y, x) for y in (y0, y0 + 1) for x in (x0, x0 + 1)]
|
|
corner_packets[link.key] = {
|
|
"direction": link.i,
|
|
"source_coarse_xy": link.source_xy,
|
|
"target_coarse_xy": link.target_xy,
|
|
"children_fine_yx": kids,
|
|
"children_centre_xy": [(float(maps.fine_centre[x]), float(maps.fine_centre[y])) for y, x in kids],
|
|
"packets": [float(f2[link.i, y, x]) for y, x in kids],
|
|
"weights": [0.25] * 4,
|
|
"coalesced": float(fc_new[link.i, ky, kx]),
|
|
}
|
|
|
|
diagnostics = {
|
|
"api_version": API_VERSION,
|
|
"collision": spec.collision,
|
|
"tick_before": tick,
|
|
"tick_after": tick + 2,
|
|
"collision_count_coarse": count_c,
|
|
"collision_count_fine": count_f,
|
|
"explosion_writes": explosion_writes,
|
|
"coarse_incoming_owner": owner, # 0 invalid (covered), 1 coarse stream, 2 coalescence
|
|
"coalesce_writes": coalesce_writes,
|
|
"coalesce_reads": coalesce_reads,
|
|
"fine_hop_sources": {"hop1": src1, "hop2": src2},
|
|
# Pre-cleanup hop-2 validity, including valid inner ghosts. The returned
|
|
# fine_f is valid exactly on maps.fine_active and NaN elsewhere.
|
|
"hop2_valid": hop2_valid,
|
|
"coarse_valid_after": new_valid,
|
|
"leaf_before": leaf_before, # (mass, momentum_x, momentum_y)
|
|
"leaf_after": leaf_after,
|
|
"budget": budget,
|
|
"corner_packets": corner_packets, # 8 entries keyed "<corner>_q<i>"
|
|
}
|
|
new_state = State(coarse_f=fc_new, fine_f=ff_new, phase="arrival_f", tick=tick + 2)
|
|
return new_state, diagnostics
|