Files
CelerisLab/tests/validation/multigrid/spacetime_union_plan.py
T

213 lines
12 KiB
Python

"""Host-only experimental union event contract; no solver or runtime admission.
Consumes existing level-wide union geometry, never its legacy transfer tables.
The frozen 707afe6 rectangle oracle remains independent and unchanged. Topology
acceptance proves packet accounting, not interpolation accuracy or GPU support.
"""
from types import SimpleNamespace
import numpy as np
import spacetime_reference as ref
class UnresolvedTopology(ValueError):
"""An experimental coupling limitation, not a production geometry rejection."""
def __init__(self, code, **detail):
self.code, self.detail = code, detail
super().__init__(f"experimental space-time topology unresolved: {code}: {detail}")
def owned_slots(layout, directions, cells):
"""Map dense logical reads/writes to existing owner-core storage, not halos."""
cells = np.asarray(cells)
if not np.issubdtype(cells.dtype, np.integer) or np.any((cells < 0) | (cells >= np.prod(layout.shape))):
raise ValueError('invalid dense cell index')
owner = layout.owner.ravel()[cells]
if np.any(owner == layout.hole):
raise UnresolvedTopology('missing_block_owner', cells=cells[owner == layout.hole].tolist())
return layout.slot(owner, directions)
def _first_entry(covered, origin, child, c):
"""First continuous union entry of one characteristic, in exact quarter times.
Probe the cells immediately before/after each crossed grid line using integer
eighth-cell coordinates. Return only actual exterior faces of the UNION.
"""
start4 = 4*origin + 2*child + 1
for t4 in (1, 3):
before = (2*start4 + c*(2*t4-1))//8
after = (2*start4 + c*(2*t4+1))//8
if not covered[tuple(before)] and covered[tuple(after)]:
axes = []
for j in range(2):
if before[j] != after[j]:
neighbor = after.copy(); neighbor[j] = before[j]
if not covered[tuple(neighbor)]: axes.append(j)
return t4, axes
return None
class UnionPlan:
"""Two-level, level-wide square geometry; all coordinates are level-local.
Assumes a complete parent domain; does not validate nested parent coverage.
``pair`` is existing FixedPatchSet2D geometry, even for one rectangle. Scalar
window origins and nested state coverage need their own contract; they are
explicitly unadmitted here. No global rectangle-gap rule is imposed.
"""
phase = 'arrival_f'
runtime_supported = False
interpolation_accuracy = 'unqualified_for_new_unions'
def __init__(self, pair):
ny, nx = pair.coarse_shape
self.origin = tuple(pair.fine_origin)
if ny != nx or pair.fine_shape != (2*ny, 2*nx) or self.origin != (0., 0.):
raise UnresolvedTopology('requires_level_wide_square_pair',
coarse_shape=pair.coarse_shape, fine_shape=pair.fine_shape, origin=self.origin)
self.n, self.nf = nx, 2*nx
covered = np.array(pair.covered, dtype=bool, copy=True)
active = np.kron(covered, np.ones((2, 2), bool))
if not np.array_equal(active, pair.fine_kind == 0):
raise UnresolvedTopology('non_pair_active_coverage')
self.maps = SimpleNamespace(coarse_active=~covered, coarse_covered=covered, fine_active=active)
self.outgoing = [[], []]
# Check transit FIRST: a fine exit may not be silently routed into a
# covered coarse cell. Returning fine packets need a separately derived rule.
for q in range(1, 9):
c = np.array((ref.cpu.CY[q], ref.cpu.CX[q]))
boundary = active & ~np.roll(active, tuple(-c), axis=(0, 1))
for y, x in zip(*np.nonzero(boundary)):
p = np.array((y, x))
for m in (0, 1):
target = (p+(2-m)*c)//2
if covered[tuple(target)]:
raise UnresolvedTopology('fine_exit_reentry', q=q, substep=m,
fine_path_yx=[(p+j*c).tolist() for j in range(3-m)], coarse_target_yx=target.tolist())
self.outgoing[m].append((q, int(y*self.nf+x), int(target[0]*nx+target[1])))
self.outgoing = [np.array(x, dtype=np.int64).reshape(-1, 3) for x in self.outgoing]
self.groups = []
for iy, ix in zip(*np.nonzero(~covered)):
origin = np.array((iy, ix))
# Two-cell production margins make all nearby probes in bounds.
if not covered[max(0,iy-1):iy+2,max(0,ix-1):ix+2].any(): continue
for q in range(1, 9):
c = np.array((ref.cpu.CY[q], ref.cpu.CX[q]))
children = [np.array((dy, dx)) for dy in (0, 1) for dx in (0, 1)]
starts = [2*origin+child for child in children]
events = []
for p in starts:
event = None
for m in (0, 1):
old, new = p+m*c, p+(m+1)*c
if not active[tuple(old)] and active[tuple(new)]:
event = (m, int(old[0]*self.nf+old[1])); break
events.append(event)
if all(e is None for e in events): continue
entries = [_first_entry(covered, origin, child, c) for child in children]
axes = sorted({j for entry in entries if entry is not None for j in entry[1]})
if not axes:
raise UnresolvedTopology('missing_exterior_support', source_yx=origin.tolist(), q=q)
points = []
for p, entry, event in zip(starts, entries, events):
position = (p+.5)/2
planes = origin+(c > 0)
times = [(planes[j]-position[j])/c[j] for j in axes]
time = max(times)
if entry is not None and (4*time != entry[0] or event is None):
raise UnresolvedTopology('nonconvex_entry_support', source_yx=origin.tolist(), q=q,
quarter_time=entry[0], support_quarter_time=float(4*time))
offset = position+c*time-(origin+.5+c*.5)
points.append((float(offset[1]), float(offset[0]), float(time-.5)) if len(axes)==2
else (float(offset[1-axes[0]]), float(time-.5)))
kind = 'corner' if len(axes)==2 else 'face'; points = tuple(points)
samples, _, evaluation, coefficients = ref.interpolation(kind, points)
indices = []
for row in samples:
shift = np.array((row[1],row[0]), int) if kind=='corner' else np.eye(2,dtype=int)[1-axes[0]]*int(row[0])
pos = (origin+shift) % nx; indices.append(int(pos[0]*nx+pos[1]))
target = (origin+c) % nx
self.groups.append(ref.Group(q,int(iy*nx+ix),int(target[0]*nx+target[1]),kind,
points,np.array(indices),evaluation@coefficients,tuple(events)))
receivers = {(int(q),int(target)) for rows in self.outgoing for q,_,target in rows}
receivers |= {(g.q,g.receiver) for g in self.groups if None in g.events}
self.receivers = np.array(sorted(receivers),dtype=np.int64).reshape(-1,2)
self.proof = validate_ownership(self)
def validate_ownership(plan):
"""Independent mask-based event coverage and four-quarter receiver accounting."""
n, nf, active = plan.n, plan.nf, plan.maps.fine_active
actual = [(m,g.q,slot) for g in plan.groups for e in g.events if e is not None for m,slot in [e]]
expected = {(m,q,int((y-ref.cpu.CY[q])*nf+x-ref.cpu.CX[q]))
for m in (0,1) for q in range(1,9) for y,x in zip(*np.nonzero(active))
if not active[y-ref.cpu.CY[q],x-ref.cpu.CX[q]]}
if len(actual)!=len(set(actual)) or set(actual)!=expected:
raise ValueError('incoming writer uniqueness/coverage failed')
packets = {tuple(map(int,row)):[] for row in plan.receivers}
if len(packets)!=len(plan.receivers): raise ValueError('duplicate coarse receiver')
anchors = set()
for i,g in enumerate(plan.groups):
if not 1<=g.q<9 or not 0<=g.source<n*n or not plan.maps.coarse_active.ravel()[g.source]:
raise ValueError('invalid coarse source anchor')
if (g.q,g.source) in anchors: raise ValueError('duplicate coarse source group')
anchors.add((g.q,g.source))
origin=np.array(divmod(g.source,n)); c=np.array((ref.cpu.CY[g.q],ref.cpu.CX[g.q]))
target=(origin+c)%n
if g.receiver!=target[0]*n+target[1]: raise ValueError('group receiver provenance failed')
ordered=[]; entries=[]
for dy,dx in ((0,0),(0,1),(1,0),(1,1)):
p=2*origin+(dy,dx); event=None
for m in (0,1):
old,new=p+m*c,p+(m+1)*c
if not active[tuple(old)] and active[tuple(new)]:
event=(m,int(old[0]*nf+old[1])); break
ordered.append(event)
entries.append(_first_entry(plan.maps.coarse_covered,origin,np.array((dy,dx)),c))
if tuple(ordered)!=g.events or all(e is None for e in ordered):
raise ValueError('ordered child characteristic provenance failed')
axes=sorted({j for e in entries if e is not None for j in e[1]})
if not axes or g.kind!=('corner' if len(axes)==2 else 'face'):
raise ValueError('group supporting axes differ')
samples,_,evaluation,coefficients=ref.interpolation(g.kind,g.points)
indices=[]
for row in samples:
shift=np.array((row[1],row[0]),int) if g.kind=='corner' else np.eye(2,dtype=int)[1-axes[0]]*int(row[0])
pos=(origin+shift)%n; indices.append(int(pos[0]*n+pos[1]))
if not np.array_equal(g.samples,indices): raise ValueError('sample anchor provenance failed')
for k,(dy,dx) in enumerate(((0,0),(0,1),(1,0),(1,1))):
position=origin+(np.array((dy,dx))+.5)/2
time=max((origin[j]+int(c[j]>0)-position[j])/c[j] for j in axes)
offset=position+c*time-(origin+.5+c*.5)
point=(float(offset[1]),float(offset[0]),float(time-.5)) if g.kind=='corner' else (float(offset[1-axes[0]]),float(time-.5))
if g.points[k]!=point: raise ValueError('evaluation point provenance failed')
source=np.zeros(len(g.samples)); source[0 if g.kind=='corner' else 1]=1
if (not np.allclose(g.weights,evaluation@coefficients,rtol=0,atol=2e-14)
or not np.allclose(g.weights.mean(axis=0),source,rtol=0,atol=2e-14)):
raise ValueError('group coefficients/directional budget failed')
for k,event in enumerate(g.events):
if event is None:
if (g.q,g.receiver) not in packets: raise ValueError('missing leftover receiver')
packets[(g.q,g.receiver)].append(('leftover',i,k))
for m,rows in enumerate(plan.outgoing):
expected_out = {(q,int(y*nf+x),int(((y+(2-m)*ref.cpu.CY[q])//2)*n+(x+(2-m)*ref.cpu.CX[q])//2))
for q in range(1,9) for y,x in zip(*np.nonzero(active))
if not active[y+ref.cpu.CY[q],x+ref.cpu.CX[q]]}
actual_out = [tuple(map(int,row)) for row in rows]
if len(actual_out)!=len(set(actual_out)) or set(actual_out)!=expected_out:
raise ValueError('outgoing uniqueness/coverage failed')
for q,slot,target in actual_out: packets[(q,target)].append(('outgoing',m,q,slot))
expected_receivers = {(q,int(y*n+x)) for y,x in zip(*np.nonzero(plan.maps.coarse_active))
for q in range(1,9) if any(active[2*y+dy-ref.cpu.CY[q],2*x+dx-ref.cpu.CX[q]]
for dy in (0,1) for dx in (0,1)
if 0<=2*y+dy-ref.cpu.CY[q]<nf and 0<=2*x+dx-ref.cpu.CX[q]<nf)}
if set(packets)!=expected_receivers: raise ValueError('coarse receiver coverage failed')
tokens=[token for row in packets.values() for token in row]
if any(len(row)!=4 for row in packets.values()) or len(tokens)!=len(set(tokens)):
raise ValueError('receiver needs four distinct quarter-packets')
return dict(topology='supported',incoming_events=len(actual),outgoing_events=sum(map(len,plan.outgoing)),
groups=len(plan.groups),receivers=len(packets),receiver_quarter_packets=len(tokens),
interpolation_accuracy=plan.interpolation_accuracy,runtime_supported=False,
scope='pair-local ownership/coefficients; complete parent domain assumed, nested coverage unvalidated')