From 5dae81220733d2402f67403c7572ceb6501e24dd Mon Sep 17 00:00:00 2001 From: Ben Whewell Date: Wed, 14 Oct 2026 22:14:34 -0600 Subject: [PATCH 1/4] Fit unstructured nearby problems with a constrained spline per material region and add a one-shot adjoint k estimate. --- CHANGELOG.md | 17 + docs/theory/nearby_problems.md | 67 ++- docs/user_guide/verification.md | 35 +- src/ndiffusion/_spline_fit.py | 767 ++++++++++++++++++++++++++++++ src/ndiffusion/nearby.py | 95 +++- tests/test_nearby_unstructured.py | 77 +++ tests/test_spline_fit.py | 178 +++++++ 7 files changed, 1218 insertions(+), 18 deletions(-) create mode 100644 src/ndiffusion/_spline_fit.py create mode 100644 tests/test_spline_fit.py diff --git a/CHANGELOG.md b/CHANGELOG.md index b2f891b..cfb1b62 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,23 @@ All notable changes to ndiffusion are recorded here. The format follows ## [Unreleased] +### Added + +- A constrained spline fit for the method of nearby problems on the + unstructured mesh: one quintic spline per material region, fit to all regions + and groups together with flux and current continuity at the interfaces and the + boundary conditions on the boundary, plus corner singular functions where the + interfaces meet at a corner. `nearby_fixed_source` and `nearby_k_eigenvalue` + use it when given `bc=`; `fit="local"` keeps the per-cell polynomial fit. +- `nearby_k_eigenvalue(adjoint_solver=...)` returns a one-shot eigenvalue error + estimate from the adjoint flux (`k_error_adjoint`) and its per-cell + contributions (`k_error_cells`); `nearby=False` skips the nearby solve. + +### Fixed + +- The nearby fit treated a one-group `nusigf` as a fission matrix (its size + matches both layouts), so a one-group `nearby_k_eigenvalue` failed. + ## [1.0.2] - 2026-10-12 ### Fixed diff --git a/docs/theory/nearby_problems.md b/docs/theory/nearby_problems.md index 0756981..a313053 100644 --- a/docs/theory/nearby_problems.md +++ b/docs/theory/nearby_problems.md @@ -36,12 +36,44 @@ each material region, since the flux gradient jumps at a material interface: |---|---| | 1-D | quintic interpolating spline over each material block (cubic or linear on short blocks) | | 2-D structured | quintic splines along each direction, per material block | -| 2-D unstructured | least-squares polynomial of up to fourth degree over a stencil of same-material neighbors | +| 2-D unstructured, `bc=` given | one constrained quintic tensor spline per material region, all fit together (below) | +| 2-D unstructured, `fit="local"` | least-squares polynomial of up to fourth degree over a stencil of same-material neighbors | The residual is the continuous operator evaluated pointwise at the cell centers. Sources are per unit volume, so one definition serves all three geometries. +### The constrained spline fit + +Independent fits on either side of an interface disagree there by about the +scheme's own error, which is enough to spoil the estimate near interfaces. On +the unstructured mesh the fit therefore ties the regions together. Each +connected material region gets one quintic tensor-product B-spline on its +bounding box, with knots about 2-3 cells apart (2.5 shortest diffusion lengths, +clipped to that range). All regions and groups are fit at once by least squares +to the cell values, with heavily weighted rows at two points on every face: + +- interface faces: $\phi$ and $D\,\partial\phi/\partial n$ continuous; +- boundary faces: the face's Robin condition. + +Two variants are used, chosen from the mesh: + +- **smooth**: also the diffusion equation and its gradient on both sides of + every interface and boundary face. Accurate on layouts whose interfaces are + straight or gently curved and run to the boundary. +- **corner**: picked when the interfaces meet at a corner (turning more than 30 + degrees) or enclose a small region (under 8 knot spans across). There the + extra rows pose an ill-conditioned problem and the fit stops matching its + data, so the diffusion equation is imposed softly at every cell center + instead, and the singular functions $r^\alpha\Phi(\theta)$ of each material + corner are added to the regions around it. The exponents $\alpha$ come from + the corner angles and diffusion coefficients; each function is cut off + radially before the shortest interface out of the corner ends, so it adds no + kink where an interface stops. + +A fixed-source fit always uses the corner variant, whose interior rows carry +the source. + ## Eigenvalue problems For a k-eigenvalue problem the fit also defines an eigenvalue, @@ -56,11 +88,32 @@ $M\phi = F\phi/k + r$ is solved by a power iteration that reuses the fixed-source solver for each outer step, and $k_\text{nearby} - k_\text{fit}$ estimates the eigenvalue error. +The same estimate is available without the nearby solve. With the adjoint flux +$\phi^+$ (from a solver run on the adjoint cross sections), + +$$ +k_h - k \approx \sum_c \eta_c, \qquad +\eta_c = \frac{\phi^+_c\,(k_h M\phi_\text{fit} - F\phi_\text{fit})_c V_c} + {\sum_c \phi^+_c (M\phi_\text{fit})_c V_c}, +$$ + +and the per-cell contributions $\eta_c$ show where the eigenvalue error comes +from. + ## Limits -The estimate is only as good as the fit. Block-wise fits are independent on -either side of a material interface and need not agree there, so the estimate is -most reliable away from interfaces; on a strongly heterogeneous problem the -eigenvalue estimate is better read as an indicator of size than as a -correction. On a homogeneous problem it recovers the error closely - see -{doc}`../user_guide/verification` for an example. +The estimate is only as good as the fit. The 1-D and structured fits are +independent on either side of a material interface and need not agree there, +so on those meshes the estimate is most reliable away from interfaces. + +On layouts with material corners or small enclosed regions the eigenvalue error +is a small net of large contributions of both signs ($\sum_c|\eta_c|$ can be +hundreds of times $|\sum_c \eta_c|$). The local flux estimate stays good there, +but the eigenvalue estimate is better read as an indicator of size, with +$\sum_c|\eta_c|$ as a loose bound. + +The fit also needs the mesh to resolve the layout: a few cells per knot span and +several knot spans across each region. On a curved interface it estimates the +error of the scheme on the polygons the mesh uses, not the error of replacing +the curve by those polygons. On a homogeneous problem it recovers the error +closely - see {doc}`../user_guide/verification` for an example. diff --git a/docs/user_guide/verification.md b/docs/user_guide/verification.md index 46671c0..aa198df 100644 --- a/docs/user_guide/verification.md +++ b/docs/user_guide/verification.md @@ -69,10 +69,37 @@ is solved by a power iteration that reuses the fixed-source solver for each outer step. Along with the flux estimate it returns `k_curve_fit` and `k_nearby`, and `k_nearby - k_curve_fit` estimates the eigenvalue error. -The fits are independent in each material, so they need not agree at a material -interface. The flux error estimate is most reliable away from interfaces, and on -a strongly heterogeneous core the eigenvalue estimate is better read as an -indicator of size than as a correction. +Given an `adjoint_solver` - a k-eigenvalue solver built on +`make_adjoint_materials(mats)` with the same mesh and boundary conditions - it +also returns `k_error_adjoint`, a one-shot estimate of the same eigenvalue error, +and `k_error_cells`, its contribution from each cell. With `nearby=False` the +nearby power iteration is skipped and only the one-shot estimate is formed. + +On structured meshes the fits are independent in each material, so they need +not agree at a material interface, and the flux estimate is most reliable away +from interfaces. + +## Unstructured meshes + +On the unstructured mesh, pass the boundary conditions the solver was built +with as `bc=`. The curve fit is then one constrained spline per material region +that keeps the flux and current continuous across interfaces and meets the +boundary conditions (see {doc}`../theory/nearby_problems`): + +```python +keig = nd.KEigenSolverUnstructured2D(mats, mesh, bc, use_cg=True) +adj = nd.KEigenSolverUnstructured2D(nd.make_adjoint_materials(mats), mesh, bc, use_cg=True) +fixed = nd.FixedSourceSolverUnstructured2D(mats, mesh, bc) +r = nd.nearby_k_eigenvalue(keig, fixed, mats, mesh=mesh, bc=bc, adjoint_solver=adj) +``` + +Without `bc=` (or with `fit="local"`) the older per-cell polynomial fit is used. + +On a layout with material corners or small enclosed regions, such as a core +with rodded assemblies, the eigenvalue error is a small net of large +contributions of both signs. The flux estimate stays good there, but read the +eigenvalue estimate as an indicator of size, with `abs(r.k_error_cells).sum()` +as a loose bound. ## Adjoint problems diff --git a/src/ndiffusion/_spline_fit.py b/src/ndiffusion/_spline_fit.py new file mode 100644 index 0000000..4a4d888 --- /dev/null +++ b/src/ndiffusion/_spline_fit.py @@ -0,0 +1,767 @@ +"""Constrained global spline fit of an unstructured-mesh flux, for MNP. + +One smooth function per connected material region and group: a quintic +tensor-product B-spline on the region's bounding box, with knots about `stride` +cell sizes apart, so the fit does not depend on how the region is meshed. All +regions and groups are fit together by sparse least squares to the cell +centroid values, plus heavily weighted constraint rows at two Gauss points on +every interface and boundary face: + + interface: phi and D dphi/dn continuous across the face + boundary: the face's own Robin condition + "smooth": also the diffusion equation and its gradient on each side + +Two configurations are needed. "smooth" (hard PDE rows on interfaces and +boundaries) is accurate on layouts without material corners, but around a +corner, or around a small region enclosed by interfaces, the PDE rows and the +continuity rows together pose an ill-conditioned Cauchy problem and the fit +stops matching its data. "corner" drops those PDE rows, imposes the diffusion +equation softly at every centroid instead, and adds the corner singular +functions r^alpha Phi(theta) of every interior material corner to the regions +around it. "auto" picks "corner" when the mesh has a material corner turning +more than 30 degrees or a small enclosed inclusion, and "smooth" otherwise. + +The constraint sets are consistent but rank-deficient, so they enter as +weighted rows of one normal-equation system rather than as exact equalities. +""" + +from collections import defaultdict, namedtuple + +import numpy as np + +K = 5 # spline degree + +SplineFit = namedtuple("SplineFit", ["curve", "lap", "config", "evaluate", "info"]) + +CONFIGS = { + "smooth": dict(pde_orders=(0, 1), interior_pde=None, enrich=False), + "corner": dict(pde_orders=(), interior_pde=(1.0, (0, 1)), enrich=True), +} +CORNER_TURN = 30.0 # degrees; a polygonal curve turns a little at every vertex + + +def _scipy(): + try: + import scipy.sparse as sp + from scipy.interpolate import BSpline + from scipy.optimize import brentq + from scipy.sparse.linalg import spsolve + except ImportError as exc: # pragma: no cover - exercised only without scipy + raise ImportError( + "The method of nearby problems needs SciPy for the unstructured " + "spline fit. Install it with `pip install ndiffusion[nearby]` (or " + "`pip install scipy`)." + ) from exc + return sp, BSpline, brentq, spsolve + + +# --------------------------------------------------------------------------- +# Mesh topology +# --------------------------------------------------------------------------- + + +def _edge_owners(mesh): + offs = np.asarray(mesh.cell_offsets) + cv = np.asarray(mesh.cell_vertices) + owner = defaultdict(list) + for c in range(len(offs) - 1): + vs = cv[offs[c]:offs[c + 1]] + for a, b in zip(vs, np.roll(vs, -1)): + owner[(min(a, b), max(a, b))].append(c) + return owner + + +def mesh_faces(mesh, mat): + """Interface faces (interior, materials differ) and boundary faces, each a + dict with midpoint, unit tangent, length and cells; boundary faces also + carry their tag and outward normal.""" + vx, vy = np.asarray(mesh.vx), np.asarray(mesh.vy) + offs, cv = np.asarray(mesh.cell_offsets), np.asarray(mesh.cell_vertices) + owner = _edge_owners(mesh) + centre = np.array([[vx[cv[offs[c]:offs[c + 1]]].mean(), vy[cv[offs[c]:offs[c + 1]]].mean()] + for c in range(len(offs) - 1)]) + tag = {(min(a, b), max(a, b)): t + for a, b, t in zip(mesh.bface_v0, mesh.bface_v1, mesh.bface_bc_tag)} + inter, bound = [], [] + for (a, b), cs in owner.items(): + tx, ty = vx[b] - vx[a], vy[b] - vy[a] + ln = np.hypot(tx, ty) + mid = np.array([0.5 * (vx[a] + vx[b]), 0.5 * (vy[a] + vy[b])]) + f = dict(mid=mid, t=np.array([tx, ty]) / ln, cells=cs, length=ln) + if len(cs) == 2 and mat[cs[0]] != mat[cs[1]]: + inter.append(f) + elif len(cs) == 1: + nrm = np.array([ty, -tx]) / ln + if np.dot(nrm, mid - centre[cs[0]]) < 0: + nrm = -nrm + f.update(n=nrm, tag=tag[(a, b)]) + bound.append(f) + return inter, bound + + +def material_corners(inter, min_turn=0.0, tol=1e-6): + """Interface vertices where the interface turns by more than `min_turn` + degrees, or where three or more interface faces meet, as an (n, 2) array.""" + by_vertex = defaultdict(list) + for f in inter: + for s in (-0.5, 0.5): + p = f["mid"] + s * f["length"] * f["t"] + by_vertex[(round(p[0], 6), round(p[1], 6))].append(f["t"]) + out = [] + for key, ts in by_vertex.items(): + if len(ts) >= 3: + out.append(key) + elif len(ts) == 2: + c = abs(np.dot(ts[0], ts[1])) + if abs(c - 1.0) > tol and np.degrees(np.arccos(min(1.0, c))) > min_turn: + out.append(key) + return np.array(out, dtype=float).reshape(-1, 2) + + +def interface_breaks(inter, min_len, tol=1e-9): + """x of every vertical and y of every horizontal run of interface faces + whose total length reaches `min_len`. Used as knots, they put the + re-entrant edges of a non-convex region on knot lines instead of cutting a + knot span (a cut basis function makes the fit spike there).""" + runs = ({}, {}) + for f in inter: + tx, ty = f["t"] + for axis, along in ((0, abs(tx) < tol), (1, abs(ty) < tol)): + if along: + c = round(float(f["mid"][axis]), 6) + runs[axis][c] = runs[axis].get(c, 0.0) + f["length"] + return tuple(tuple(sorted(c for c, L in r.items() if L >= min_len)) for r in runs) + + +# --------------------------------------------------------------------------- +# Tensor B-spline regions +# --------------------------------------------------------------------------- + + +def _knots(a, b, hk, breaks=()): + """Clamped knot vector on [a, b] with spacing <= hk; every coordinate in + `breaks` strictly inside (a, b) becomes a simple knot.""" + eps = 1e-9 * max(1.0, b - a) + pts = [a] + sorted(p for p in set(breaks) if a + eps < p < b - eps) + [b] + inner = [] + for lo, hi in zip(pts[:-1], pts[1:]): + n = max(1, int(np.ceil((hi - lo) / hk - 1e-9))) + inner.extend(np.linspace(lo, hi, n + 1)[1:]) + return np.r_[[a] * (K + 1), inner[:-1], [b] * (K + 1)] + + +def _basis(x, t, nu=0): + """Dense (len(x), n_basis) matrix of the nu-th derivative of the B-splines.""" + _, BSpline, _, _ = _scipy() + x = np.clip(np.asarray(x, float), t[K], t[-K - 1]) + if nu == 0: + return BSpline.design_matrix(x, t, K).toarray() + return BSpline(t, np.eye(len(t) - K - 1), K).derivative(nu)(x) + + +class _Region: + def __init__(self, cells, material, box, hk, breaks): + x0, x1, y0, y1 = box + self.cells, self.material = cells, material + self.tx, self.ty = _knots(x0, x1, hk, breaks[0]), _knots(y0, y1, hk, breaks[1]) + self.nx, self.ny = len(self.tx) - K - 1, len(self.ty) - K - 1 + self.n = self.nx * self.ny + + def rows(self, x, y, dx=0, dy=0): + """Dense rows of d^dx/dx d^dy/dy of the tensor basis at the points.""" + Bx, By = _basis(x, self.tx, dx), _basis(y, self.ty, dy) + return (Bx[:, :, None] * By[:, None, :]).reshape(len(x), -1) + + def rows_sparse(self, x, y, dx=0, dy=0): + """The same rows, sparse: a tensor row has at most (K+1)^2 nonzeros.""" + sp = _scipy()[0] + Bx, By = _basis(x, self.tx, dx), _basis(y, self.ty, dy) + ix = np.argsort(-np.abs(Bx), axis=1)[:, :K + 1] + iy = np.argsort(-np.abs(By), axis=1)[:, :K + 1] + vx, vy = np.take_along_axis(Bx, ix, 1), np.take_along_axis(By, iy, 1) + cols = (ix[:, :, None] * self.ny + iy[:, None, :]).reshape(len(x), -1) + vals = (vx[:, :, None] * vy[:, None, :]).reshape(len(x), -1) + rr = np.repeat(np.arange(len(x)), cols.shape[1]) + return sp.csr_matrix((vals.ravel(), (rr, cols.ravel())), shape=(len(x), self.n)) + + +def regions(mesh, mat, area, stride, breaks=((), ())): + """Connected same-material cell sets with their spline bases, and the + region index of every cell.""" + vx, vy = np.asarray(mesh.vx), np.asarray(mesh.vy) + offs, cv = np.asarray(mesh.cell_offsets), np.asarray(mesh.cell_vertices) + nc = len(offs) - 1 + adj = defaultdict(list) + for cs in _edge_owners(mesh).values(): + if len(cs) == 2 and mat[cs[0]] == mat[cs[1]]: + adj[cs[0]].append(cs[1]) + adj[cs[1]].append(cs[0]) + comp = -np.ones(nc, int) + out = [] + for c0 in range(nc): + if comp[c0] >= 0: + continue + stack, members = [c0], [] + comp[c0] = len(out) + while stack: + c = stack.pop() + members.append(c) + for n in adj[c]: + if comp[n] < 0: + comp[n] = len(out) + stack.append(n) + members = np.array(sorted(members)) + vs = np.concatenate([cv[offs[c]:offs[c + 1]] for c in members]) + hk = stride * np.sqrt(np.median(area[members])) + box = (vx[vs].min(), vx[vs].max(), vy[vs].min(), vy[vs].max()) + out.append(_Region(members, mat[c0], box, hk, breaks)) + return out, comp + + +def small_inclusions(mesh, bc, n_groups, regs, comp, area, stride, max_spans=8.0): + """(material, spans) of every region that the mirrored domain encloses and + whose equal-area diameter is under `max_spans` knot spans (absorber rods, + small seeds). Enclosed: it touches no non-reflective boundary and does not + run across the whole domain in x or y (the mirrors would extend that one + into an infinite strip).""" + owner = _edge_owners(mesh) + vx, vy = np.asarray(mesh.vx), np.asarray(mesh.vy) + offs, cv = np.asarray(mesh.cell_offsets), np.asarray(mesh.cell_vertices) + tol = 1e-9 * max(np.ptp(vx), np.ptp(vy)) + + def reflective(t): + return all(bc[t * n_groups + g].A == 0.0 for g in range(n_groups)) + + open_regions = {comp[owner[(min(a, b), max(a, b))][0]] + for a, b, t in zip(mesh.bface_v0, mesh.bface_v1, mesh.bface_bc_tag) + if not reflective(t)} + out = [] + for i, R in enumerate(regs): + if i in open_regions: + continue + vs = np.concatenate([cv[offs[c]:offs[c + 1]] for c in R.cells]) + if (vx[vs].min() - vx.min() < tol and vx.max() - vx[vs].max() < tol) or \ + (vy[vs].min() - vy.min() < tol and vy.max() - vy[vs].max() < tol): + continue + hk = stride * np.sqrt(np.median(area[R.cells])) + spans = 2.0 * np.sqrt(area[R.cells].sum() / np.pi) / hk + if spans < max_spans: + out.append((int(R.material), float(spans))) + return out + + +def _diff_penalty(nx, ny, order): + """order-th coefficient differences along x and y (P-spline penalty) for a + coefficient grid stored ix * ny + iy.""" + sp = _scipy()[0] + + def d1(n): + D = sp.identity(n, format="csr") + for _ in range(order): + D = (sp.eye(D.shape[0] - 1, D.shape[0], 1) - sp.eye(D.shape[0] - 1, D.shape[0])) @ D + return D + + return sp.vstack([sp.kron(d1(nx), sp.identity(ny)), sp.kron(sp.identity(nx), d1(ny))]) + + +# --------------------------------------------------------------------------- +# Corner singular functions +# --------------------------------------------------------------------------- +# +# At an interior material corner the flux behaves like smooth + c r^alpha +# Phi(theta), with Phi = A_k cos(alpha t) + B_k sin(alpha t) in sector k. +# Continuity of Phi and D Phi' across the rays gives a transfer matrix per +# sector; an exponent is a root of trace(M) = 2 for the product M. Each mode +# is harmonic in every sector and continuous in value and current across the +# rays, so it can be added to every region around the corner. That only holds +# while the rays are interfaces: where an interface stops (IAEA assembly +# edges), the ray's extension crosses one material and an uncut mode kinks +# there, a line source the pointwise residual cannot see. So each mode is cut +# off radially, chi = 1 to R1 and 0 past R2, with R2 the shortest run of the +# same interface along the rays (infinite when every ray reaches the boundary). + + +def _transfer(alpha, w, D): + c, s = np.cos(alpha * w), np.sin(alpha * w) + return np.array([[c, s / (D * alpha)], [-D * alpha * s, c]]) + + +def _trace_minus_two(alpha, sectors): + M = np.eye(2) + for w, D in sectors: + M = _transfer(alpha, w, D) @ M + return np.trace(M) - 2.0 + + +def corner_exponents(sectors, amax=2.0, n=4000): + """Roots alpha in (0, amax] of trace(M) = 2 for sectors [(width, D), ...].""" + brentq = _scipy()[2] + a = np.linspace(1e-4, amax, n) + f = np.array([_trace_minus_two(x, sectors) for x in a]) + roots = [] + for i in range(n - 1): + if f[i] == 0.0: + roots.append(a[i]) + elif f[i] * f[i + 1] < 0: + roots.append(brentq(_trace_minus_two, a[i], a[i + 1], args=(sectors,))) + # double roots (alpha = 1 in a homogeneous medium) touch zero without a sign change + for i in range(1, n - 1): + if abs(f[i]) < 1e-6 and abs(f[i]) <= abs(f[i - 1]) and abs(f[i]) <= abs(f[i + 1]): + if not any(abs(a[i] - r) < 1e-3 for r in roots): + roots.append(a[i]) + return sorted(roots) + + +class CornerMode: + """chi(r) r^alpha Phi(theta) for one group at one corner.""" + + def __init__(self, x0, y0, alpha, group, starts, widths, AB, regions, R1=np.inf, R2=np.inf): + self.x0, self.y0, self.alpha, self.group = x0, y0, alpha, group + self.starts, self.widths, self.AB = np.asarray(starts), np.asarray(widths), np.asarray(AB) + self.regions = set(regions) + self.scale = 1.0 + self.R1, self.R2 = R1, R2 + + def _chi(self, r): + """Cutoff and its first three r-derivatives (C^3 septic step).""" + out = [np.ones_like(r), np.zeros_like(r), np.zeros_like(r), np.zeros_like(r)] + if not np.isfinite(self.R2): + return out + w = self.R2 - self.R1 + t = np.clip((r - self.R1) / w, 0.0, 1.0) + s = 35 * t**4 - 84 * t**5 + 70 * t**6 - 20 * t**7 + s1 = 140 * t**3 - 420 * t**4 + 420 * t**5 - 140 * t**6 + s2 = 420 * t**2 - 1680 * t**3 + 2100 * t**4 - 840 * t**5 + s3 = 840 * t - 5040 * t**2 + 8400 * t**3 - 4200 * t**4 + return [1 - s, -s1 / w, -s2 / w**2, -s3 / w**3] + + def _sector(self, theta): + rel = np.mod(theta - self.starts[0], 2 * np.pi) + edges = np.cumsum(np.r_[0.0, self.widths]) + k = np.clip(np.searchsorted(edges, rel, side="right") - 1, 0, len(self.widths) - 1) + return k, rel - edges[k] + + def _raw(self, x, y): + x = np.atleast_1d(np.asarray(x, float)) + y = np.atleast_1d(np.asarray(y, float)) + dx, dy = x - self.x0, y - self.y0 + r = np.maximum(np.hypot(dx, dy), 1e-12) + th = np.arctan2(dy, dx) + k, t = self._sector(th) + A, B, a = self.AB[k, 0], self.AB[k, 1], self.alpha + Phi = A * np.cos(a * t) + B * np.sin(a * t) + dPhi = a * (-A * np.sin(a * t) + B * np.cos(a * t)) + val = self.scale * r**a * Phi + dr = self.scale * a * r ** (a - 1) * Phi + dth = self.scale * r**a * dPhi + gx = np.cos(th) * dr - np.sin(th) / r * dth + gy = np.sin(th) * dr + np.cos(th) / r * dth + return r, th, val, gx, gy + + def eval(self, x, y): + """Value and gradient of chi m.""" + r, th, val, gx, gy = self._raw(x, y) + c, c1, _, _ = self._chi(r) + return c * val, c * gx + c1 * np.cos(th) * val, c * gy + c1 * np.sin(th) * val + + def lap(self, x, y): + """Laplacian of chi m and its gradient: m is harmonic, so + lap(chi m) = m (chi'' + (2 alpha + 1) chi' / r).""" + r, th, val, gx, gy = self._raw(x, y) + _, c1, c2, c3 = self._chi(r) + b = 2 * self.alpha + 1 + q = c2 + b * c1 / r + q1 = c3 + b * (c2 / r - c1 / r**2) + return q * val, q * gx + q1 * np.cos(th) * val, q * gy + q1 * np.sin(th) * val + + +def _ray_run(x0, y0, theta, inter, mat, pair, bverts, h_cell): + """Length of the straight interface between materials `pair` leaving + (x0, y0) at angle theta; inf if it reaches the domain boundary.""" + u = np.array([np.cos(theta), np.sin(theta)]) + tol = 1e-6 * max(1.0, h_cell) + segs = [] + for f in inter: + if abs(f["t"][0] * u[1] - f["t"][1] * u[0]) > 1e-6: + continue + d = f["mid"] - (x0, y0) + if abs(d[0] * u[1] - d[1] * u[0]) > tol or {mat[c] for c in f["cells"]} != pair: + continue + s = d @ u + segs.append((s - 0.5 * f["length"], s + 0.5 * f["length"])) + end = 0.0 + for a, b in sorted(segs): + if b <= tol: + continue + if a > end + tol: + break + end = max(end, b) + p = np.array([x0, y0]) + end * u + if len(bverts) and np.min(np.hypot(*(bverts - p).T)) < tol: + return np.inf + return end + + +def corner_modes(mesh, inter, comp, mats, cx, cy, h_cell, alpha_max=2.5, tol=1e-3, + min_turn=0.0, cutoff=0.5): + """Singular modes at every interior material corner, each cut off (chi = 1 + to cutoff * R2, 0 past R2) before its shortest interface ray ends.""" + G = mats.n_groups + D = np.asarray(mats.D).reshape(-1, G) + mat = np.asarray(mesh.material_id) + vx, vy = np.asarray(mesh.vx), np.asarray(mesh.vy) + bv = np.unique(np.r_[np.asarray(mesh.bface_v0), np.asarray(mesh.bface_v1)]) + bverts = np.column_stack([vx[bv], vy[bv]]) + modes = [] + for x0, y0 in material_corners(inter, min_turn=min_turn): + rays = [] + for f in inter: + a = f["mid"] - 0.5 * f["length"] * f["t"] + b = f["mid"] + 0.5 * f["length"] * f["t"] + for p, q in ((a, b), (b, a)): + if np.hypot(*(p - (x0, y0))) < 1e-6 * max(1.0, h_cell): + rays.append(np.arctan2(q[1] - y0, q[0] - x0)) + rays = np.sort(np.mod(rays, 2 * np.pi)) + if len(rays) < 2: + continue + rays = rays[np.r_[True, np.diff(rays) > 1e-9]] + widths = np.diff(np.r_[rays, rays[0] + 2 * np.pi]) + sec_mat, sec_reg = [], [] + for th0, w in zip(rays, widths): + th = th0 + 0.5 * w + px, py = x0 + 2.0 * h_cell * np.cos(th), y0 + 2.0 * h_cell * np.sin(th) + c = int(np.argmin(np.hypot(cx - px, cy - py))) + sec_mat.append(mat[c]) + sec_reg.append(comp[c]) + R2 = min(_ray_run(x0, y0, th, inter, mat, {sec_mat[i - 1], sec_mat[i]}, bverts, h_cell) + for i, th in enumerate(rays)) + for g in range(G): + sectors = [(w, D[m, g]) for w, m in zip(widths, sec_mat)] + for alpha in corner_exponents(sectors, amax=alpha_max): + if abs(alpha - round(alpha)) < tol: + continue + M = np.eye(2) + Ts = [] + for w, Dk in sectors: + T = _transfer(alpha, w, Dk) + Ts.append(T) + M = T @ M + s = np.linalg.svd(M - np.eye(2))[2][-1] # (Phi, D Phi') at the first ray + AB = [] + for (w, Dk), T in zip(sectors, Ts): + AB.append((s[0], s[1] / (Dk * alpha))) + s = T @ s + AB = np.array(AB) + peak = max(np.abs(A * np.cos(alpha * t) + B * np.sin(alpha * t)).max() + for (A, B), t in zip(AB, (np.linspace(0, w, 50) for w, _ in sectors))) + modes.append(CornerMode(x0, y0, alpha, g, rays, widths, AB / peak, sec_reg, + R1=cutoff * R2, R2=R2)) + return modes + + +# --------------------------------------------------------------------------- +# The fit +# --------------------------------------------------------------------------- + + +def _interior_rows(regs, cx, cy, off, N, N_spline, modes, A_op, D, G, orders): + """PDE rows at every centroid, in the cell's own region: + -D_g lap(phi_g) + sum_g' A[g,g'] phi_g' (order 0) and its gradient (order 1).""" + sp = _scipy()[0] + blocks = [] + for r, R in enumerate(regs): + x, y = cx[R.cells], cy[R.cells] + mm = R.material + A = A_op(mm) + for order in orders: + for ex, ey in ([(0, 0)] if order == 0 else [(1, 0), (0, 1)]): + val = R.rows_sparse(x, y, ex, ey) + lap = R.rows_sparse(x, y, 2 + ex, ey) + R.rows_sparse(x, y, ex, 2 + ey) + mv = [] + for e, md in enumerate(modes): + if r in md.regions: + v, gx, gy = md.eval(x, y) + lv, lgx, lgy = md.lap(x, y) + mv.append((e, md.group, v if order == 0 else (gx if ex else gy), + lv if order == 0 else (lgx if ex else lgy))) + for g in range(G): + ii, jj, vv = [], [], [] + for gp in range(G): + b = A[g, gp] * val + if gp == g: + b = b - D[mm, g] * lap + b = b.tocoo() + ii.append(b.row) + jj.append(off[(gp, r)] + b.col) + vv.append(b.data) + for e, ge, v, lv in mv: + col = A[g, ge] * v - (D[mm, g] * lv if ge == g else 0.0) + if np.any(col): + ii.append(np.arange(len(R.cells))) + jj.append(np.full(len(R.cells), N_spline + e)) + vv.append(col) + blocks.append(sp.csr_matrix((np.concatenate(vv), (np.concatenate(ii), np.concatenate(jj))), + shape=(len(R.cells), N))) + return sp.vstack(blocks, format="csr") + + +def choose_config(mats, mesh, bc, stride, inclusion_spans=8.0): + """"corner" if the mesh has a material corner turning more than 30 degrees + or a small enclosed inclusion, else "smooth"; also returns the inclusions.""" + from ndiffusion._core import cell_areas + + mat = np.asarray(mesh.material_id) + inter, _ = mesh_faces(mesh, mat) + if len(material_corners(inter, min_turn=CORNER_TURN)): + return "corner", [] + inclusions = [] + if inclusion_spans > 0: + area = np.asarray(cell_areas(mesh)) + regs, comp = regions(mesh, mat, area, stride) + inclusions = small_inclusions(mesh, bc, mats.n_groups, regs, comp, area, stride, inclusion_spans) + return ("corner" if inclusions else "smooth"), inclusions + + +def fit(mats, mesh, flux, k, bc, stride=2.0, config="auto", weight=1e3, ridge=1e-10, + penalty=None, source=None, inclusion_spans=8.0): + """Fit the cell-centroid flux on an unstructured mesh. + + `k` is the eigenvalue whose fission term enters the PDE rows (pass + ``np.inf`` for a fixed-source problem, with `source` the per-unit-volume + source at the centroids); `bc` is the solver's tag-indexed Robin list. + `stride` is the knot spacing in cell sizes. `penalty` = (order, lam) adds a + P-spline difference penalty, which a strongly graded mesh needs (knots are + spaced from the median cell, so large cells leave coefficients unsupported). + + Returns ``SplineFit(curve, lap, config, evaluate, info)``: the fit and its + Laplacian at the centroids, both (n_cells, G); the configuration used; and + ``evaluate(cells, x, y)``, the fit at points inside the given cells.""" + sp, _, _, spsolve = _scipy() + from ndiffusion._core import cell_areas, cell_centroids + + inclusions = [] + if config == "auto": + config, inclusions = choose_config(mats, mesh, bc, stride, inclusion_spans) + if config not in CONFIGS: + raise ValueError(f"config must be 'auto' or one of {sorted(CONFIGS)}, got {config!r}") + pde_orders, interior_pde, enrich = (CONFIGS[config][key] for key in ("pde_orders", "interior_pde", "enrich")) + if source is not None: + if pde_orders: + raise ValueError("a fixed-source fit needs config='corner' (the hard PDE rows carry no source)") + interior_pde = (interior_pde[0], (0,)) + + G = mats.n_groups + D = np.asarray(mats.D).reshape(-1, G) + Rm = np.asarray(mats.removal).reshape(-1, G) + S = np.asarray(mats.scatter).reshape(-1, G, G) + chi = np.asarray(mats.chi).reshape(-1, G) + nf = np.asarray(mats.nusigf).reshape(-1, G) + + def A_op(m): + return np.diag(Rm[m]) - S[m] - np.outer(chi[m], nf[m]) / k + + flux = np.asarray(flux, float).reshape(-1, G) + mat = np.asarray(mesh.material_id) + cx, cy = (np.asarray(a) for a in cell_centroids(mesh)) + area = np.asarray(cell_areas(mesh)) + inter, bound = mesh_faces(mesh, mat) + h_cell = np.sqrt(np.median(area)) + breaks = interface_breaks(inter, 3.0 * h_cell) + regs, comp = regions(mesh, mat, area, stride, breaks) + off, N = {}, 0 + for g in range(G): + for r, R in enumerate(regs): + off[(g, r)] = N + N += R.n + gscale = 1.0 / np.abs(flux).max(axis=0) + modes = [] + if enrich: + modes = corner_modes(mesh, inter, comp, mats, cx, cy, h_cell, min_turn=CORNER_TURN) + for md in modes: + cells = np.where(np.isin(comp, list(md.regions)))[0] + md.scale = 1.0 / max(np.abs(md.eval(cx[cells], cy[cells])[0]).max(), 1e-300) + N_spline = N + N = N_spline + len(modes) + + def mode_terms(r, px, py, deriv, toward=None): + """{column: (group, value)} of the modes reaching region r at one point: + deriv None for the value, (nx, ny) for a directional derivative, or + ("lap", None | (nx, ny)) for the Laplacian or its derivative. A mode's + gradient jumps on its rays, so `toward` (a point inside the region) + moves the evaluation a hair to that side.""" + out = {} + if toward is not None: + d = np.array([toward[0] - px[0], toward[1] - py[0]]) + d /= np.linalg.norm(d) + px, py = px + 1e-8 * d[0], py + 1e-8 * d[1] + for e, md in enumerate(modes): + if r not in md.regions: + continue + lap_ = deriv is not None and deriv[0] == "lap" + v, gx, gy = md.lap(px, py) if lap_ else md.eval(px, py) + d = deriv[1] if lap_ else deriv + out[N_spline + e] = (md.group, v[0] if d is None else d[0] * gx[0] + d[1] * gy[0]) + return out + + # data rows + di, dj, dv, rhs = [], [], [], [] + row = 0 + for r, R in enumerate(regs): + B = R.rows_sparse(cx[R.cells], cy[R.cells]).tocoo() + for g in range(G): + di.append(row + B.row) + dj.append(off[(g, r)] + B.col) + dv.append(B.data * gscale[g]) + for e, md in enumerate(modes): + if md.group == g and r in md.regions: + di.append(row + np.arange(len(R.cells))) + dj.append(np.full(len(R.cells), N_spline + e)) + dv.append(md.eval(cx[R.cells], cy[R.cells])[0] * gscale[g]) + rhs.append(flux[R.cells, g] * gscale[g]) + row += len(R.cells) + Adata = sp.csr_matrix((np.concatenate(dv), (np.concatenate(di), np.concatenate(dj))), shape=(row, N)) + y = np.concatenate(rhs) + + # constraint rows, each a {column offset: dense 1-row block} + crows = [] + + def pde(r, px, py, toward): + R = regs[r] + mm = R.material + A = A_op(mm) + for order in pde_orders: + for ex, ey in ([(0, 0)] if order == 0 else [(1, 0), (0, 1)]): + lapr = R.rows(px, py, 2 + ex, ey) + R.rows(px, py, ex, 2 + ey) + val = R.rows(px, py, ex, ey) + mt_ = mode_terms(r, px, py, None if order == 0 else (ex, ey), toward) + ml_ = mode_terms(r, px, py, ("lap", None if order == 0 else (ex, ey)), toward) + for g in range(G): + blocks = {} + for gp in range(G): + b = A[g, gp] * val + if gp == g: + b = b - D[mm, g] * lapr + if np.any(b): + blocks[off[(gp, r)]] = b + for col, (ge, v) in mt_.items(): + c = A[g, ge] * v - (D[mm, g] * ml_[col][1] if ge == g else 0.0) + if c != 0.0: + blocks[col] = np.array([[c]]) + crows.append(blocks) + + gp2 = np.array([-1, 1]) / np.sqrt(3) / 2 # two Gauss points per face + for f in inter: + c1, c2 = f["cells"] + r1, r2 = comp[c1], comp[c2] + nrm = np.array([f["t"][1], -f["t"][0]]) + m1, m2 = regs[r1].material, regs[r2].material + t1, t2 = (cx[c1], cy[c1]), (cx[c2], cy[c2]) + for s_ in gp2: + p = f["mid"] + s_ * f["length"] * f["t"] + px, py = p[:1], p[1:] + v1, v2 = regs[r1].rows(px, py), regs[r2].rows(px, py) + n1 = nrm[0] * regs[r1].rows(px, py, 1, 0) + nrm[1] * regs[r1].rows(px, py, 0, 1) + n2 = nrm[0] * regs[r2].rows(px, py, 1, 0) + nrm[1] * regs[r2].rows(px, py, 0, 1) + val1, val2 = mode_terms(r1, px, py, None, t1), mode_terms(r2, px, py, None, t2) + dn1, dn2 = mode_terms(r1, px, py, tuple(nrm), t1), mode_terms(r2, px, py, tuple(nrm), t2) + for g in range(G): + blk = {off[(g, r1)]: v1, off[(g, r2)]: -v2} + for col in set(val1) | set(val2): + if (val1.get(col) or val2.get(col))[0] == g: + blk[col] = np.array([[val1.get(col, (g, 0.0))[1] - val2.get(col, (g, 0.0))[1]]]) + crows.append(blk) + blk = {off[(g, r1)]: D[m1, g] * n1, off[(g, r2)]: -D[m2, g] * n2} + for col in set(dn1) | set(dn2): + if (dn1.get(col) or dn2.get(col))[0] == g: + blk[col] = np.array([[D[m1, g] * dn1.get(col, (g, 0.0))[1] + - D[m2, g] * dn2.get(col, (g, 0.0))[1]]]) + crows.append(blk) + pde(r1, px, py, t1) + pde(r2, px, py, t2) + for f in bound: + c0 = f["cells"][0] + r = comp[c0] + R = regs[r] + tc = (cx[c0], cy[c0]) + for s_ in gp2: + p = f["mid"] + s_ * f["length"] * f["t"] + px, py = p[:1], p[1:] + v = R.rows(px, py) + dn = f["n"][0] * R.rows(px, py, 1, 0) + f["n"][1] * R.rows(px, py, 0, 1) + mv, md_ = mode_terms(r, px, py, None, tc), mode_terms(r, px, py, tuple(f["n"]), tc) + for g in range(G): + b = bc[f["tag"] * G + g] + blk = {off[(g, r)]: b.A * v + b.B * dn} + for col, (ge, val_) in mv.items(): + if ge == g: + blk[col] = np.array([[b.A * val_ + b.B * md_[col][1]]]) + crows.append(blk) + pde(r, px, py, tc) + + ci, cj, cvv = [], [], [] + for i, blocks in enumerate(crows): + norm = np.sqrt(sum(np.sum(b**2) for b in blocks.values())) + for o, b in blocks.items(): + jj = np.nonzero(b[0])[0] + ci.append(np.full(len(jj), i)) + cj.append(o + jj) + cvv.append(b[0, jj] / norm) + C = (sp.csr_matrix((np.concatenate(cvv), (np.concatenate(ci), np.concatenate(cj))), shape=(len(crows), N)) + if crows else sp.csr_matrix((0, N))) + data_w = np.sqrt(Adata.multiply(Adata).sum(axis=1)).mean() + w = weight * data_w + Nm = Adata.T @ Adata + (w**2) * (C.T @ C) + ridge * sp.identity(N) + rhs_total = Adata.T @ y + if interior_pde is not None: + w_int, orders = interior_pde + Ci = _interior_rows(regs, cx, cy, off, N, N_spline, modes, A_op, D, G, orders) + rn = np.maximum(np.sqrt(np.asarray(Ci.multiply(Ci).sum(axis=1)).ravel()), 1e-300) + Ci = sp.diags(1.0 / rn) @ Ci + wi = w_int * data_w + Nm = Nm + (wi**2) * (Ci.T @ Ci) + if source is not None: + q = np.asarray(source, float).reshape(-1, G) + qrows = np.concatenate([q[R.cells, g] for R in regs for g in range(G)]) + rhs_total = rhs_total + (wi**2) * (Ci.T @ (qrows / rn)) + if penalty is not None: + order, lam = penalty + Dp = sp.block_diag([_diff_penalty(R.nx, R.ny, order) for g in range(G) for R in regs], format="csr") + if modes: + Dp = sp.hstack([Dp, sp.csr_matrix((Dp.shape[0], len(modes)))], format="csr") + Nm = Nm + lam * (Adata.multiply(Adata).sum() / N) * (Dp.T @ Dp) + coef = spsolve(Nm.tocsc(), rhs_total) + + curve = np.zeros_like(flux) + lap = np.zeros_like(flux) + for r, R in enumerate(regs): + x, yy = cx[R.cells], cy[R.cells] + V = R.rows_sparse(x, yy) + L = R.rows_sparse(x, yy, 2, 0) + R.rows_sparse(x, yy, 0, 2) + for g in range(G): + cg = coef[off[(g, r)]:off[(g, r)] + R.n] + curve[R.cells, g] = V @ cg + lap[R.cells, g] = L @ cg + for e, md in enumerate(modes): + if r in md.regions: + curve[R.cells, md.group] += coef[N_spline + e] * md.eval(x, yy)[0] + lap[R.cells, md.group] += coef[N_spline + e] * md.lap(x, yy)[0] + + def evaluate(cells, x, y): + """The fit at points (x, y), each inside the matching cell, as (n, G).""" + cells, x, y = np.asarray(cells), np.asarray(x, float), np.asarray(y, float) + out = np.zeros((len(x), G)) + for r in np.unique(comp[cells]): + sel = np.where(comp[cells] == r)[0] + R = regs[r] + V = R.rows_sparse(x[sel], y[sel]) + for g in range(G): + out[sel, g] = V @ coef[off[(g, r)]:off[(g, r)] + R.n] + for e, md in enumerate(modes): + if r in md.regions: + out[sel, md.group] += coef[N_spline + e] * md.eval(x[sel], y[sel])[0] + return out + + info = dict(n_unknowns=N, n_constraints=len(crows), regions=len(regs), inclusions=inclusions, + modes=[(md.x0, md.y0, md.alpha, md.group, md.R2) for md in modes], + constraint_residual=np.linalg.norm(C @ coef) / max(1.0, np.linalg.norm(coef))) + return SplineFit(curve, lap, config, evaluate, info) diff --git a/src/ndiffusion/nearby.py b/src/ndiffusion/nearby.py index ac43e11..53522bf 100644 --- a/src/ndiffusion/nearby.py +++ b/src/ndiffusion/nearby.py @@ -24,8 +24,16 @@ eigenvalue error via a curve-fit ``k`` and a nearby ``k``. The leakage term ``-div(D grad phi)`` needs the curve fit's second derivatives; -these come from quintic/cubic splines on structured meshes (via scipy) and from -a per-cell least-squares quadratic reconstruction on the unstructured mesh. +these come from quintic/cubic splines on structured meshes (via scipy). On the +unstructured mesh, given the boundary conditions, they come from one constrained +quintic tensor spline per material region, fit to all regions and groups at once +with flux and current continuity at the interfaces and the Robin conditions on +the boundary (``_spline_fit``); without them, from a per-cell least-squares +polynomial reconstruction (``fit="local"``). + +For a k-eigenvalue problem the eigenvalue error can also be estimated in one +shot from the adjoint flux, ``sum phi+ (k L phi_cf - F phi_cf) V / sum phi+ L +phi_cf V``, which skips the nearby power iteration. Everything here is orchestration in Python: the actual solves reuse the compiled ``FixedSourceSolver*`` classes, so no C++ changes are required. Sources are @@ -59,7 +67,10 @@ "k_nearby", "residual", "nearby_rate", + "k_error_adjoint", + "k_error_cells", ], + defaults=(None, None), ) @@ -97,7 +108,7 @@ def _xs_arrays(mats): chi = np.asarray(mats.chi, dtype=float).reshape(n_mat, G) scatter = np.asarray(mats.scatter, dtype=float).reshape(n_mat, G, G) nusigf = np.asarray(mats.nusigf, dtype=float) - if nusigf.size == n_mat * G * G: + if _fission_matrix_mode(mats): # by size alone, one group would look like a matrix nusigf = nusigf.reshape(n_mat, G, G) else: nusigf = nusigf.reshape(n_mat, G) @@ -461,11 +472,39 @@ def _curvelap_unstructured(mats, mesh, flux): return curve, lap, area, mat_ids -def _fit_dispatch(mats, flux, medium_map, edges_x, edges_y, geometry, mesh): +def _spline_stride(mats, area): + """Knot spacing in cell sizes: 2.5 of the shortest diffusion length + sqrt(D / Sigma_r), kept between 2 and 3 cells.""" + L_min = float(np.min(np.sqrt(np.asarray(mats.D) / np.asarray(mats.removal)))) + return float(np.clip(2.5 * L_min / np.sqrt(np.median(area)), 2.0, 3.0)) + + +def _curvelap_spline(mats, mesh, flux, bc, k, source): + from ndiffusion import _spline_fit + from ndiffusion._core import cell_areas + + area = np.asarray(cell_areas(mesh)) + config = "auto" if source is None else "corner" # only the corner rows carry a source + f = _spline_fit.fit(mats, mesh, flux, k, bc, stride=_spline_stride(mats, area), + config=config, source=source) + return f.curve, f.lap, area, np.asarray(mesh.material_id, dtype=int) + + +def _fit_dispatch(mats, flux, medium_map, edges_x, edges_y, geometry, mesh, + fit=None, bc=None, k=np.inf, source=None): """Pick the curve fitter from the supplied geometry descriptor.""" flux = np.asarray(flux, dtype=float) + if fit not in (None, "spline", "local"): + raise ValueError(f"fit must be None, 'spline' or 'local', got {fit!r}") if mesh is not None: + if fit == "spline" or (fit is None and bc is not None): + if bc is None: + raise ValueError("fit='spline' needs bc=, the tag-indexed boundary " + "conditions the solver was built with") + return _curvelap_spline(mats, mesh, flux, bc, k, source) return _curvelap_unstructured(mats, mesh, flux) + if fit == "spline": + raise ValueError("fit='spline' is for the unstructured mesh (mesh=)") if edges_y is not None: return _curvelap_2d(mats, medium_map, edges_x, edges_y, geometry, flux) if edges_x is not None: @@ -491,6 +530,8 @@ def nearby_fixed_source( edges_y=None, geometry=None, mesh=None, + bc=None, + fit=None, return_nearby=True, ): """Run the method of nearby problems for a fixed-source diffusion solve. @@ -509,6 +550,13 @@ def nearby_fixed_source( ``(medium_map, edges_x, geometry)`` for 1-D, ``(medium_map, edges_x, edges_y, geometry)`` for 2-D structured, or ``mesh=`` for the unstructured FVM mesh. + bc : list of BoundaryCondition, optional + Unstructured mesh only: the tag-indexed boundary conditions the solver + was built with. Given ``bc``, the curve fit is the constrained spline + fit; without it, the per-cell polynomial fit. + fit : {None, "spline", "local"} + Force the unstructured curve fit; ``None`` picks ``"spline"`` when + ``bc`` is given. return_nearby : bool, default True If True, also solve the nearby problem and return the error estimate. @@ -528,7 +576,8 @@ def nearby_fixed_source( flux = np.asarray(numerical.flux, dtype=float) curve, lap, _vol, mat_ids = _fit_dispatch( - mats, flux, medium_map, edges_x, edges_y, geometry, mesh + mats, flux, medium_map, edges_x, edges_y, geometry, mesh, + fit=fit, bc=bc, k=np.inf, source=source, ) loss = _assemble_loss(mats, mat_ids, curve, lap) # fixed-source: no fission residual = loss - source @@ -551,6 +600,10 @@ def nearby_k_eigenvalue( edges_y=None, geometry=None, mesh=None, + bc=None, + fit=None, + adjoint_solver=None, + nearby=True, tol=1e-8, max_iter=200, ): @@ -570,6 +623,17 @@ def nearby_k_eigenvalue( mats : Materials medium_map, edges_x, edges_y, geometry, mesh Geometry descriptor, as in :func:`nearby_fixed_source`. + bc, fit + The unstructured curve fit, as in :func:`nearby_fixed_source`. + adjoint_solver : matching KEigen* solver, optional + Built on :func:`make_adjoint_materials` of ``mats`` (same mesh and + boundary conditions). Its flux gives the one-shot eigenvalue error + estimate ``k_error_adjoint`` and its per-cell contributions + ``k_error_cells``. + nearby : bool, default True + Run the nearby power iteration. With ``nearby=False`` (and an + ``adjoint_solver``) only the fit and the one-shot estimate are formed; + ``nearby_flux`` and ``k_nearby`` are then None. tol : float Convergence tolerance on the flux change and ``|dk|``. max_iter : int @@ -582,13 +646,19 @@ def nearby_k_eigenvalue( ``nearby_flux``, ``k_nearby``, ``residual``, and ``nearby_rate`` (the curve-fit fission production integral). The arrays are shaped ``(n_cells, G)``, like the solver flux. - ``k_nearby - k_curve_fit`` estimates the eigenvalue discretization error. + ``k_nearby - k_curve_fit`` estimates the eigenvalue discretization error, + and so does ``k_error_adjoint`` when an ``adjoint_solver`` is given + (``k_error_cells`` sums to it). On layouts whose error contributions + cancel - material corners, small inclusions - ``k_error_cells`` and the + local error are more reliable than the net eigenvalue estimate; the sum + of ``abs(k_error_cells)`` is a (loose) bound. """ numerical = keig_solver.solve() flux = np.asarray(numerical.flux, dtype=float) curve, lap, vol, mat_ids = _fit_dispatch( - mats, flux, medium_map, edges_x, edges_y, geometry, mesh + mats, flux, medium_map, edges_x, edges_y, geometry, mesh, + fit=fit, bc=bc, k=numerical.keff, ) volG = vol[:, None] @@ -603,6 +673,15 @@ def nearby_k_eigenvalue( residual = loss - prod / k_cf # (n_cells, G), exact source of the nearby problem + k_error = cells = None + if adjoint_solver is not None: + adj = np.asarray(adjoint_solver.solve().flux, dtype=float).reshape(-1, mats.n_groups) + cells = (np.sum(adj * (numerical.keff * loss - prod) * volG, axis=1) + / float(np.sum(adj * loss * volG))) + k_error = float(cells.sum()) + if not nearby: + return NearbyKResult(numerical, curve, k_cf, None, None, residual, rate, k_error, cells) + # --- Nearby power iteration (reuses the fixed-source solver's operator M) --- phi = curve.copy() k = k_cf @@ -631,4 +710,6 @@ def nearby_k_eigenvalue( k_nearby=k, residual=residual, nearby_rate=rate, + k_error_adjoint=k_error, + k_error_cells=cells, ) diff --git a/tests/test_nearby_unstructured.py b/tests/test_nearby_unstructured.py index 06821ee..1886bbd 100644 --- a/tests/test_nearby_unstructured.py +++ b/tests/test_nearby_unstructured.py @@ -198,3 +198,80 @@ def test_error_estimate_tracks_true_error(self): _, te, est = _triangle_error(30, 24) ratio = np.linalg.norm(est) / np.linalg.norm(te) assert 0.5 < ratio < 3.0 + + +ZERO_BC = [nd.BoundaryCondition(A=1.0, B=0.0)] + + +def _solve_mms_spline(mesh): + cx, cy, _area, _neigh = nearby._unstructured_geometry(mesh) + phi_exact = np.sin(AX * cx) * np.sin(AY * cy) + source = (D * (AX**2 + AY**2) + SIG_A) * phi_exact + m = one_group_absorber() + solver = nd.FixedSourceSolverUnstructured2D( + m, mesh, ZERO_BC, epsilon=1e-12, max_inner=40000, omega=1.8, verbose=False) + result = nd.nearby_fixed_source(solver, m, source, mesh=mesh, bc=ZERO_BC) + te = result.numerical.flux[:, 0] - phi_exact + return np.linalg.norm(result.error_estimate[:, 0]) / np.linalg.norm(te) + + +class TestNearbyUnstructuredSpline: + """The constrained spline fit (picked when bc= is given) on the same + manufactured problem, over every cell, boundary rows included.""" + + def test_quad_estimate_is_exact(self): + assert _solve_mms_spline(quad_mesh(40, 32, LX, LY)) == pytest.approx(1.0, abs=0.01) + + def test_triangle_estimate_is_close(self): + # the local fit gives 1.55-2.07 here; measured 1.04-1.07 with the spline + assert 0.9 < _solve_mms_spline(triangle_mesh(30, 24, LX, LY)) < 1.2 + + def test_spline_needs_bc(self): + mesh = quad_mesh(8, 8, LX, LY) + m = one_group_absorber() + solver = nd.FixedSourceSolverUnstructured2D(m, mesh, ZERO_BC, epsilon=1e-10, verbose=False) + with pytest.raises(ValueError, match="bc="): + nd.nearby_fixed_source(solver, m, np.ones(64), mesh=mesh, fit="spline") + with pytest.raises(ValueError, match="fit must be"): + nd.nearby_fixed_source(solver, m, np.ones(64), mesh=mesh, bc=ZERO_BC, fit="global") + + +class TestNearbyKSpline: + """Bare quarter square, one group: phi = cos(pi x / 2L) cos(pi y / 2L) and + k = nu_f / (Sigma_a + D B^2) with B^2 = 2 (pi / 2L)^2.""" + + L = 10.0 + + def setup(self, n): + m = nd.Materials() + m.n_mat, m.n_groups = 1, 1 + m.D, m.removal, m.scatter, m.chi, m.nusigf = [1.0], [0.1], [0.0], [1.0], [0.13] + mesh = quad_mesh(n, n, self.L, self.L) + vx, vy = np.asarray(mesh.vx), np.asarray(mesh.vy) + far = self.L - 1e-9 # zero flux on x = L and y = L, reflective on the axes + mesh.bface_bc_tag = [int(min(vx[a], vx[b]) > far or min(vy[a], vy[b]) > far) + for a, b in zip(mesh.bface_v0, mesh.bface_v1)] + bc = [nd.BoundaryCondition(A=0.0, B=1.0), nd.BoundaryCondition(A=1.0, B=0.0)] + kw = dict(epsilon=1e-8, max_outer=10000, use_cg=True) + keig = nd.KEigenSolverUnstructured2D(m, mesh, bc, **kw) + adj = nd.KEigenSolverUnstructured2D(nd.make_adjoint_materials(m), mesh, bc, **kw) + fixed = nd.FixedSourceSolverUnstructured2D(m, mesh, bc, epsilon=1e-12, max_inner=40000, + omega=1.8, verbose=False) + k_exact = 0.13 / (0.1 + 2 * (np.pi / (2 * self.L)) ** 2) + return m, mesh, bc, keig, adj, fixed, k_exact + + def test_adjoint_and_nearby_estimates_match_the_true_error(self): + m, mesh, bc, keig, adj, fixed, k_exact = self.setup(20) + r = nd.nearby_k_eigenvalue(keig, fixed, m, mesh=mesh, bc=bc, adjoint_solver=adj, + tol=1e-10, max_iter=2000) + true = r.numerical.keff - k_exact + assert r.k_error_adjoint / true == pytest.approx(1.0, abs=0.05) + assert (r.k_nearby - r.k_curve_fit) / true == pytest.approx(1.0, abs=0.05) + assert r.k_error_cells.shape == (400,) + assert r.k_error_cells.sum() == pytest.approx(r.k_error_adjoint) + + def test_one_shot_only(self): + m, mesh, bc, keig, adj, fixed, _ = self.setup(10) + r = nd.nearby_k_eigenvalue(keig, fixed, m, mesh=mesh, bc=bc, adjoint_solver=adj, nearby=False) + assert r.nearby_flux is None and r.k_nearby is None + assert np.isfinite(r.k_error_adjoint) diff --git a/tests/test_spline_fit.py b/tests/test_spline_fit.py new file mode 100644 index 0000000..c9e4e62 --- /dev/null +++ b/tests/test_spline_fit.py @@ -0,0 +1,178 @@ +"""Tests for the constrained global spline fit used by MNP on unstructured meshes +(ndiffusion._spline_fit): corner exponents and modes, the choice between the +"smooth" and "corner" configurations, and the fit itself on a problem with an +analytic solution.""" + +import numpy as np +import pytest + +import ndiffusion as nd + +pytest.importorskip("scipy") +from ndiffusion import _spline_fit as sf # noqa: E402 + +REFLECT = nd.BoundaryCondition(A=0.0, B=1.0) +ZERO = nd.BoundaryCondition(A=1.0, B=0.0) + + +def quad_mesh(n, x0, x1, material, tag): + """n x n quads on [x0, x1]^2; material(x, y) per cell centre, tag(x, y) per + boundary face midpoint.""" + h = (x1 - x0) / n + vid = lambda i, j: i * (n + 1) + j # noqa: E731 + m = nd.UnstructuredMesh2D() + m.vx = [x0 + i * h for i in range(n + 1) for j in range(n + 1)] + m.vy = [x0 + j * h for i in range(n + 1) for j in range(n + 1)] + cv, offs, mats = [], [0], [] + for i in range(n): + for j in range(n): + cv += [vid(i, j), vid(i + 1, j), vid(i + 1, j + 1), vid(i, j + 1)] + offs.append(len(cv)) + mats.append(material(x0 + (i + 0.5) * h, x0 + (j + 0.5) * h)) + m.cell_vertices, m.cell_offsets, m.material_id = cv, offs, mats + bv0, bv1, tags = [], [], [] + for k in range(n): + for a, b in ((vid(k, 0), vid(k + 1, 0)), (vid(n, k), vid(n, k + 1)), + (vid(k + 1, n), vid(k, n)), (vid(0, k + 1), vid(0, k))): + bv0.append(a) + bv1.append(b) + tags.append(tag(0.5 * (m.vx[a] + m.vx[b]), 0.5 * (m.vy[a] + m.vy[b]))) + m.bface_v0, m.bface_v1, m.bface_bc_tag = bv0, bv1, tags + return m + + +def two_group(n_mat): + data = [dict(D=[1.4, 0.4], Siga=[0.01, 0.08 + 0.02 * i], Scat=[[0.0, 0.0], [0.02, 0.0]], + nuSigf=[0.0, 0.12], chi=[1.0, 0.0]) for i in range(n_mat)] + data[-1]["D"] = [1.0, 0.25] + return nd.make_materials(data, 2) + + +class TestCornerExponents: + @pytest.mark.parametrize("ratio", [0.1, 0.5, 10.0]) + def test_roots_satisfy_the_trace_condition(self, ratio): + q = np.pi / 2 + sectors = [(q, 1.0), (q, 1.0), (q, ratio), (q, 1.0)] + roots = sf.corner_exponents(sectors, amax=2.5) + assert roots and min(roots) < 1.0 # a material contrast makes the corner singular + for a in roots: + assert abs(sf._trace_minus_two(a, sectors)) < 1e-6 # double roots only to ~1e-6 + + def test_inverse_ratio_gives_the_same_exponents(self): + q = np.pi / 2 + a = sf.corner_exponents([(q, 1.0), (q, 1.0), (q, 0.2), (q, 1.0)]) + b = sf.corner_exponents([(q, 1.0), (q, 1.0), (q, 5.0), (q, 1.0)]) + np.testing.assert_allclose(a, b, atol=1e-8) + + +def corner_problem(): + """A 270/90-degree material corner at the origin of [-10, 10]^2.""" + mats = two_group(2) + mesh = quad_mesh(20, -10.0, 10.0, lambda x, y: int(x > 0 and y > 0), lambda x, y: 0) + inter, _ = sf.mesh_faces(mesh, np.asarray(mesh.material_id)) + area = np.asarray(nd.cell_areas(mesh)) + regs, comp = sf.regions(mesh, np.asarray(mesh.material_id), area, 2.0) + cx, cy = (np.asarray(a) for a in nd.cell_centroids(mesh)) + return mats, mesh, inter, comp, cx, cy + + +class TestCornerModes: + def test_modes_found_at_the_corner_only(self): + mats, mesh, inter, comp, cx, cy = corner_problem() + modes = sf.corner_modes(mesh, inter, comp, mats, cx, cy, 1.0, min_turn=sf.CORNER_TURN) + assert modes + assert all(abs(md.x0) < 1e-9 and abs(md.y0) < 1e-9 for md in modes) + assert all(np.isinf(md.R2) for md in modes) # both interfaces run to the boundary + + def test_value_and_current_continuous_across_the_rays(self): + mats, mesh, inter, comp, cx, cy = corner_problem() + D = np.asarray(mats.D).reshape(-1, 2) + for md in sf.corner_modes(mesh, inter, comp, mats, cx, cy, 1.0, min_turn=sf.CORNER_TURN): + for th in md.starts: + t = np.array([np.cos(th), np.sin(th)]) + n = np.array([-t[1], t[0]]) + p = 3.0 * t + (va, gxa, gya), (vb, gxb, gyb) = (md.eval(*(p + s * 1e-7 * n)) for s in (1, -1)) + Da = D[int((p + 0.5 * n)[0] > 0 and (p + 0.5 * n)[1] > 0), md.group] + Db = D[int((p - 0.5 * n)[0] > 0 and (p - 0.5 * n)[1] > 0), md.group] + assert va[0] == pytest.approx(vb[0], abs=1e-5) + ja = Da * (gxa[0] * n[0] + gya[0] * n[1]) + jb = Db * (gxb[0] * n[0] + gyb[0] * n[1]) + assert ja == pytest.approx(jb, rel=1e-4, abs=1e-6) + + def test_cut_off_laplacian_matches_finite_differences(self): + md = sf.CornerMode(0.0, 0.0, 0.73, 0, [0.0], [2 * np.pi], np.array([[1.0, 0.4]]), [0], R1=2.0, R2=5.0) + h = 1e-4 + for x, y in ((1.0, 0.7), (2.6, 1.9), (-3.1, 2.2), (0.4, -4.4)): + fd = (md.eval(x + h, y)[0] + md.eval(x - h, y)[0] + md.eval(x, y + h)[0] + + md.eval(x, y - h)[0] - 4 * md.eval(x, y)[0]) / h**2 + assert md.lap(x, y)[0][0] == pytest.approx(fd[0], abs=1e-4) + assert md.lap(1.0, 0.7)[0][0] == 0.0 # harmonic inside R1 + + +class TestConfig: + def test_material_corner_picks_corner(self): + mats, mesh, *_ = corner_problem() + assert sf.choose_config(mats, mesh, [REFLECT] * 2, 2.0)[0] == "corner" + + def test_straight_stripes_pick_smooth(self): + mats = two_group(2) + mesh = quad_mesh(20, 0.0, 20.0, lambda x, y: int(x > 10), lambda x, y: 0) + assert sf.choose_config(mats, mesh, [REFLECT] * 2, 2.0) == ("smooth", []) + + def test_small_enclosed_region_is_an_inclusion(self): + mats = two_group(2) + mesh = quad_mesh(30, 0.0, 30.0, lambda x, y: int(abs(x - 15) < 2 and abs(y - 15) < 2), + lambda x, y: 1 if x > 29.9 or y > 29.9 else 0) + bc = [REFLECT, REFLECT, ZERO, ZERO] + area = np.asarray(nd.cell_areas(mesh)) + mat = np.asarray(mesh.material_id) + regs, comp = sf.regions(mesh, mat, area, 2.0) + inc = sf.small_inclusions(mesh, bc, 2, regs, comp, area, 2.0) + assert [m for m, _ in inc] == [1] # the rod, not the open background + assert sf.small_inclusions(mesh, bc, 2, regs, comp, area, 2.0, max_spans=1.0) == [] + + +class TestFit: + L = 10.0 + + def solve(self, n): + m = nd.Materials() + m.n_mat, m.n_groups = 1, 1 + m.D, m.removal, m.scatter, m.chi, m.nusigf = [1.0], [0.1], [0.0], [1.0], [0.13] + mesh = quad_mesh(n, 0.0, self.L, lambda x, y: 0, lambda x, y: 1 if max(x, y) > self.L - 1e-9 else 0) + bc = [REFLECT, ZERO] + r = nd.KEigenSolverUnstructured2D(m, mesh, bc, epsilon=1e-8, max_outer=10000, use_cg=True).solve() + return m, mesh, bc, r + + def test_cosine_mode_curvature_and_k_estimate(self): + """Bare quarter square: phi = cos(pi x / 2L) cos(pi y / 2L), so + lap(phi) / phi = -2 (pi / 2L)^2 and k = nu_f / (Sigma_a + D B^2).""" + m, mesh, bc, r = self.solve(20) + phi = np.asarray(r.flux) + f = sf.fit(m, mesh, phi, r.keff, bc, stride=2.0) + assert f.config == "smooth" + assert np.abs(f.curve - phi).max() < 1e-3 * np.abs(phi).max() + B2 = 2 * (np.pi / (2 * self.L)) ** 2 + cx, cy = (np.asarray(a) for a in nd.cell_centroids(mesh)) + inner = (cx < 0.8 * self.L) & (cy < 0.8 * self.L) + np.testing.assert_allclose(f.lap[inner, 0] / f.curve[inner, 0], -B2, rtol=1e-2) + + # one-shot adjoint estimate of the k error (the operator is self-adjoint here) + area = np.asarray(nd.cell_areas(mesh)) + loss = m.D[0] * -f.lap[:, 0] + m.removal[0] * f.curve[:, 0] + prod = m.nusigf[0] * f.curve[:, 0] + eta = np.sum(phi[:, 0] * (r.keff * loss - prod) * area) / np.sum(phi[:, 0] * loss * area) + true = r.keff - m.nusigf[0] / (m.removal[0] + m.D[0] * B2) + assert eta / true == pytest.approx(1.0, abs=0.05) + + def test_fixed_source_needs_the_corner_config(self): + m, mesh, bc, r = self.solve(10) + with pytest.raises(ValueError, match="corner"): + sf.fit(m, mesh, np.asarray(r.flux), np.inf, bc, config="smooth", source=np.ones((100, 1))) + + def test_evaluate_reproduces_the_centroid_values(self): + m, mesh, bc, r = self.solve(10) + f = sf.fit(m, mesh, np.asarray(r.flux), r.keff, bc) + cx, cy = (np.asarray(a) for a in nd.cell_centroids(mesh)) + np.testing.assert_allclose(f.evaluate(np.arange(len(cx)), cx, cy), f.curve, rtol=1e-12, atol=1e-14) From 8b76e6d4428a5a6efadb0ef8779baaebc013099b Mon Sep 17 00:00:00 2001 From: Ben Whewell Date: Wed, 14 Oct 2026 23:11:19 -0600 Subject: [PATCH 2/4] Fix the ruff line-length and unused-variable errors in the spline fit. --- src/ndiffusion/_spline_fit.py | 35 +++++++++++++++++++++++------------ tests/test_spline_fit.py | 1 - 2 files changed, 23 insertions(+), 13 deletions(-) diff --git a/src/ndiffusion/_spline_fit.py b/src/ndiffusion/_spline_fit.py index 4a4d888..08ed6ad 100644 --- a/src/ndiffusion/_spline_fit.py +++ b/src/ndiffusion/_spline_fit.py @@ -448,7 +448,7 @@ def corner_modes(mesh, inter, comp, mats, cx, cy, h_cell, alpha_max=2.5, tol=1e- M = T @ M s = np.linalg.svd(M - np.eye(2))[2][-1] # (Phi, D Phi') at the first ray AB = [] - for (w, Dk), T in zip(sectors, Ts): + for (_, Dk), T in zip(sectors, Ts): AB.append((s[0], s[1] / (Dk * alpha))) s = T @ s AB = np.array(AB) @@ -500,8 +500,8 @@ def _interior_rows(regs, cx, cy, off, N, N_spline, modes, A_op, D, G, orders): ii.append(np.arange(len(R.cells))) jj.append(np.full(len(R.cells), N_spline + e)) vv.append(col) - blocks.append(sp.csr_matrix((np.concatenate(vv), (np.concatenate(ii), np.concatenate(jj))), - shape=(len(R.cells), N))) + ij = (np.concatenate(ii), np.concatenate(jj)) + blocks.append(sp.csr_matrix((np.concatenate(vv), ij), shape=(len(R.cells), N))) return sp.vstack(blocks, format="csr") @@ -518,7 +518,8 @@ def choose_config(mats, mesh, bc, stride, inclusion_spans=8.0): if inclusion_spans > 0: area = np.asarray(cell_areas(mesh)) regs, comp = regions(mesh, mat, area, stride) - inclusions = small_inclusions(mesh, bc, mats.n_groups, regs, comp, area, stride, inclusion_spans) + inclusions = small_inclusions(mesh, bc, mats.n_groups, regs, comp, area, stride, + inclusion_spans) return ("corner" if inclusions else "smooth"), inclusions @@ -544,10 +545,13 @@ def fit(mats, mesh, flux, k, bc, stride=2.0, config="auto", weight=1e3, ridge=1e config, inclusions = choose_config(mats, mesh, bc, stride, inclusion_spans) if config not in CONFIGS: raise ValueError(f"config must be 'auto' or one of {sorted(CONFIGS)}, got {config!r}") - pde_orders, interior_pde, enrich = (CONFIGS[config][key] for key in ("pde_orders", "interior_pde", "enrich")) + preset = CONFIGS[config] + pde_orders, interior_pde = preset["pde_orders"], preset["interior_pde"] + enrich = preset["enrich"] if source is not None: if pde_orders: - raise ValueError("a fixed-source fit needs config='corner' (the hard PDE rows carry no source)") + raise ValueError("a fixed-source fit needs config='corner' " + "(the hard PDE rows carry no source)") interior_pde = (interior_pde[0], (0,)) G = mats.n_groups @@ -619,7 +623,8 @@ def mode_terms(r, px, py, deriv, toward=None): dv.append(md.eval(cx[R.cells], cy[R.cells])[0] * gscale[g]) rhs.append(flux[R.cells, g] * gscale[g]) row += len(R.cells) - Adata = sp.csr_matrix((np.concatenate(dv), (np.concatenate(di), np.concatenate(dj))), shape=(row, N)) + Adata = sp.csr_matrix((np.concatenate(dv), (np.concatenate(di), np.concatenate(dj))), + shape=(row, N)) y = np.concatenate(rhs) # constraint rows, each a {column offset: dense 1-row block} @@ -663,12 +668,14 @@ def pde(r, px, py, toward): n1 = nrm[0] * regs[r1].rows(px, py, 1, 0) + nrm[1] * regs[r1].rows(px, py, 0, 1) n2 = nrm[0] * regs[r2].rows(px, py, 1, 0) + nrm[1] * regs[r2].rows(px, py, 0, 1) val1, val2 = mode_terms(r1, px, py, None, t1), mode_terms(r2, px, py, None, t2) - dn1, dn2 = mode_terms(r1, px, py, tuple(nrm), t1), mode_terms(r2, px, py, tuple(nrm), t2) + dn1 = mode_terms(r1, px, py, tuple(nrm), t1) + dn2 = mode_terms(r2, px, py, tuple(nrm), t2) for g in range(G): blk = {off[(g, r1)]: v1, off[(g, r2)]: -v2} for col in set(val1) | set(val2): if (val1.get(col) or val2.get(col))[0] == g: - blk[col] = np.array([[val1.get(col, (g, 0.0))[1] - val2.get(col, (g, 0.0))[1]]]) + jump = val1.get(col, (g, 0.0))[1] - val2.get(col, (g, 0.0))[1] + blk[col] = np.array([[jump]]) crows.append(blk) blk = {off[(g, r1)]: D[m1, g] * n1, off[(g, r2)]: -D[m2, g] * n2} for col in set(dn1) | set(dn2): @@ -706,8 +713,11 @@ def pde(r, px, py, toward): ci.append(np.full(len(jj), i)) cj.append(o + jj) cvv.append(b[0, jj] / norm) - C = (sp.csr_matrix((np.concatenate(cvv), (np.concatenate(ci), np.concatenate(cj))), shape=(len(crows), N)) - if crows else sp.csr_matrix((0, N))) + if crows: + ij = (np.concatenate(ci), np.concatenate(cj)) + C = sp.csr_matrix((np.concatenate(cvv), ij), shape=(len(crows), N)) + else: + C = sp.csr_matrix((0, N)) data_w = np.sqrt(Adata.multiply(Adata).sum(axis=1)).mean() w = weight * data_w Nm = Adata.T @ Adata + (w**2) * (C.T @ C) + ridge * sp.identity(N) @@ -725,7 +735,8 @@ def pde(r, px, py, toward): rhs_total = rhs_total + (wi**2) * (Ci.T @ (qrows / rn)) if penalty is not None: order, lam = penalty - Dp = sp.block_diag([_diff_penalty(R.nx, R.ny, order) for g in range(G) for R in regs], format="csr") + Dp = sp.block_diag([_diff_penalty(R.nx, R.ny, order) for _ in range(G) for R in regs], + format="csr") if modes: Dp = sp.hstack([Dp, sp.csr_matrix((Dp.shape[0], len(modes)))], format="csr") Nm = Nm + lam * (Adata.multiply(Adata).sum() / N) * (Dp.T @ Dp) diff --git a/tests/test_spline_fit.py b/tests/test_spline_fit.py index c9e4e62..858fe93 100644 --- a/tests/test_spline_fit.py +++ b/tests/test_spline_fit.py @@ -121,7 +121,6 @@ def test_straight_stripes_pick_smooth(self): assert sf.choose_config(mats, mesh, [REFLECT] * 2, 2.0) == ("smooth", []) def test_small_enclosed_region_is_an_inclusion(self): - mats = two_group(2) mesh = quad_mesh(30, 0.0, 30.0, lambda x, y: int(abs(x - 15) < 2 and abs(y - 15) < 2), lambda x, y: 1 if x > 29.9 or y > 29.9 else 0) bc = [REFLECT, REFLECT, ZERO, ZERO] From 6e6913f758c62033e1a08a179c0219d65627d642 Mon Sep 17 00:00:00 2001 From: Ben Whewell Date: Sat, 17 Oct 2026 09:11:08 -0600 Subject: [PATCH 3/4] Assemble the unstructured operator as a sparse matrix so the nearby solve and its adjoints can use the exact discrete transpose, and estimate errors in normalized quantities of interest. --- CHANGELOG.md | 11 ++ docs/api/python.md | 4 + docs/theory/nearby_problems.md | 59 +++++++- docs/user_guide/verification.md | 34 +++++ src/ndiffusion/__init__.py | 7 + src/ndiffusion/fv_operator.py | 223 ++++++++++++++++++++++++++++ src/ndiffusion/nearby.py | 251 ++++++++++++++++++++++++++++---- tests/test_fv_operator.py | 211 +++++++++++++++++++++++++++ 8 files changed, 774 insertions(+), 26 deletions(-) create mode 100644 src/ndiffusion/fv_operator.py create mode 100644 tests/test_fv_operator.py diff --git a/CHANGELOG.md b/CHANGELOG.md index 6c611cd..eecd767 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -17,6 +17,17 @@ All notable changes to ndiffusion are recorded here. The format follows - `nearby_k_eigenvalue(adjoint_solver=...)` returns a one-shot eigenvalue error estimate from the adjoint flux (`k_error_adjoint`) and its per-cell contributions (`k_error_cells`); `nearby=False` skips the nearby solve. +- `unstructured_operator(mats, mesh, bc)` assembles the unstructured solvers' + loss and production operators as sparse matrices, non-orthogonal correction + included, so a converged solve satisfies them to the solver's tolerance. +- `nearby_k_eigenvalue(method="direct")` factorizes that operator once for the + nearby solve and takes the one-shot estimate from the exact discrete adjoint. + A solver run on adjoint cross sections does not transpose the non-orthogonal + correction, which biases adjoint-weighted estimates by 1-2% on triangles and + hexagons. +- `nearby_qoi` estimates the error of normalized quantities of interest (ratios + of reaction-rate integrals, such as region powers) from a generalized adjoint, + with per-cell contributions for adaptation. ### Fixed diff --git a/docs/api/python.md b/docs/api/python.md index b636c7a..ab6d6c9 100644 --- a/docs/api/python.md +++ b/docs/api/python.md @@ -106,9 +106,13 @@ from C++; the rest is pure Python. nearby_fixed_source nearby_k_eigenvalue + nearby_qoi + unstructured_operator fission_source NearbyFixedResult NearbyKResult + NearbyQoIResult + UnstructuredOperator ``` ## Modules diff --git a/docs/theory/nearby_problems.md b/docs/theory/nearby_problems.md index a313053..f2dd404 100644 --- a/docs/theory/nearby_problems.md +++ b/docs/theory/nearby_problems.md @@ -98,7 +98,64 @@ k_h - k \approx \sum_c \eta_c, \qquad $$ and the per-cell contributions $\eta_c$ show where the eigenvalue error comes -from. +from. The two estimates agree closely on smooth layouts; where the error is a net +of large contributions of both signs (see below) they can differ by 5-15%. + +## The discrete operator and its adjoint + +The unstructured solvers never form a matrix, and they treat the non-orthogonal +part of the face flux as a deferred source that is iterated to convergence. What +a converged solve satisfies is nonetheless one linear system, and +`unstructured_operator` assembles it: the loss operator $\mathbf{M}$, including +the non-orthogonal correction $\mathbf{C}$ with the solver's least-squares +gradients, and the production operator $\mathbf{F}$, both volume-integrated. + +The explicit matrix gives the exact discrete adjoint. A solver run on adjoint +cross sections transposes the scattering and fission but not $\mathbf{C}$, +because the least-squares gradient correction is not symmetric; on triangles and +hexagons that leaves an error of 1-2% in an adjoint-weighted estimate that does +not shrink with the mesh size. On orthogonal meshes $\mathbf{C} = 0$ and the two +agree. + +With `method="direct"`, `nearby_k_eigenvalue` factorizes $\mathbf{M}$ once and +uses it for every step of the nearby power iteration, and takes the adjoint for +the one-shot estimate from the bordered system below. This is usually far faster +than the iterative route. + +## Quantities of interest + +For a normalized quantity of interest, a ratio of two reaction-rate integrals + +$$ +J(\phi) = \frac{\int a\cdot\phi\,dV}{\int b\cdot\phi\,dV} +$$ + +(a region power, or the absorption in an absorber per unit production), the +nearby estimate is $J(\phi_\text{nearby}) - J(\phi_\text{fit})$. `nearby_qoi` +also gives it without the nearby solve. Linearizing about the discrete solution, + +$$ +J(\phi_\text{nearby}) - J(\phi_\text{fit}) \approx +\Gamma^T\mathbf{V}r + J(\phi_h) - J(\phi_\text{fit}), +$$ + +where $\Gamma$ solves the generalized adjoint problem +$(\mathbf{M} - \mathbf{F}/k_h)^T\Gamma = \partial J/\partial\phi$ with +$\Gamma^T\mathbf{F}\phi_h = 0$. The ratio makes $J$ independent of the +eigenvector's normalization, so $\partial J/\partial\phi$ is orthogonal to +$\phi_h$, which is exactly the condition for this singular system to have a +solution. It is solved through the bordered matrix + +$$ +\begin{pmatrix}(\mathbf{M} - \mathbf{F}/k_h)^T & \mathbf{F}\phi_h\\ +(\mathbf{F}\phi_h)^T & 0\end{pmatrix}, +$$ + +one factorization for any number of quantities. $J(\phi_\text{fit})$ is +integrated with a degree-4 rule over each cell. The per-cell products +$\Gamma_c (\mathbf{V}r)_c$ form a goal-oriented error map that can drive mesh +refinement. The difference from the nearby-solve estimate is the curvature of +$J$; it shrinks with the mesh size. ## Limits diff --git a/docs/user_guide/verification.md b/docs/user_guide/verification.md index aa198df..3072970 100644 --- a/docs/user_guide/verification.md +++ b/docs/user_guide/verification.md @@ -95,6 +95,40 @@ r = nd.nearby_k_eigenvalue(keig, fixed, mats, mesh=mesh, bc=bc, adjoint_solver=a Without `bc=` (or with `fit="local"`) the older per-cell polynomial fit is used. +With `method="direct"` the solver's operator is assembled as a sparse matrix +(`unstructured_operator`) and factorized once. The nearby problem is then solved +with that factorization, the one-shot eigenvalue estimate uses the exact discrete +adjoint, and neither a fixed-source nor an adjoint solver is needed. On meshes of +a few tens of thousands of cells this is the faster route: + +```python +r = nd.nearby_k_eigenvalue(keig, None, mats, mesh=mesh, bc=bc, method="direct") +``` + +### Quantities of interest + +`nearby_qoi` estimates the error of a normalized quantity of interest, a ratio +of two reaction-rate integrals $\int a\cdot\phi\,dV / \int b\cdot\phi\,dV$ +with weights given per cell and group. For example, the absorption in a set of +rodded cells per unit fission production: + +```python +G = mats.n_groups +mat = np.asarray(mesh.material_id) +sig_a = ... # (n_mat, G) absorption cross sections +nusigf = np.asarray(mats.nusigf).reshape(-1, G) +num = sig_a[mat] * (mat == RODDED)[:, None] # (n_cells, G) +den = nusigf[mat] +q = nd.nearby_qoi(r, mats, mesh, bc, num, den) +q.value, q.error_adjoint, q.error_nearby # J(phi_h) and two estimates of its error +``` + +`error_adjoint` comes from a generalized adjoint solve on the exact discrete +operator, `error_nearby` from the nearby flux (when it was computed), and +`cells` holds the per-cell contributions to `error_adjoint`, a map of where the +error in that quantity comes from. Stacking the weights as +`(n_q, n_cells, G)` evaluates several quantities with one factorization. + On a layout with material corners or small enclosed regions, such as a core with rodded assemblies, the eigenvalue error is a small net of large contributions of both signs. The flux estimate stays good there, but read the diff --git a/src/ndiffusion/__init__.py b/src/ndiffusion/__init__.py index cb2ec61..9c593c6 100644 --- a/src/ndiffusion/__init__.py +++ b/src/ndiffusion/__init__.py @@ -27,6 +27,7 @@ ) from ndiffusion.adjoint import make_adjoint_materials from ndiffusion.create import boundary_conditions, make_materials, make_medium_map +from ndiffusion.fv_operator import UnstructuredOperator, unstructured_operator from ndiffusion.kinetics import ( DELAYED_U235_6GROUP, make_delayed_data, @@ -36,9 +37,11 @@ from ndiffusion.nearby import ( NearbyFixedResult, NearbyKResult, + NearbyQoIResult, fission_source, nearby_fixed_source, nearby_k_eigenvalue, + nearby_qoi, ) from ndiffusion.postprocess import ( KAPPA_U235, @@ -108,9 +111,13 @@ "make_adjoint_materials", "nearby_fixed_source", "nearby_k_eigenvalue", + "nearby_qoi", + "unstructured_operator", "fission_source", "NearbyFixedResult", "NearbyKResult", + "NearbyQoIResult", + "UnstructuredOperator", # post-processing "cell_volumes", "reaction_rate", diff --git a/src/ndiffusion/fv_operator.py b/src/ndiffusion/fv_operator.py new file mode 100644 index 0000000..59e72a8 --- /dev/null +++ b/src/ndiffusion/fv_operator.py @@ -0,0 +1,223 @@ +"""The unstructured finite-volume operator as a sparse matrix. + +The unstructured solvers are matrix-free and treat the non-orthogonal correction +as a deferred source, iterated to convergence. A converged fixed-source solve +therefore satisfies one linear system, + + M phi = V q, M = diag(a_diag) - offdiag(D_harm a_coef) - scatter V - C, + +where ``C phi`` is the face sum of ``D_f (w0 grad_0 + (1 - w0) grad_1) . T`` with +the solver's least-squares cell gradients, and a k-eigenvalue solve satisfies +``M phi = F phi / k`` with ``F`` the volume-integrated fission production. +:func:`unstructured_operator` assembles ``M`` and ``F`` explicitly, mirroring +``solver_2d_unstructured.cpp`` (``fill_face_geometry``, ``build_lsq_gradients``, +``build_diagonals``, ``cell_gradients``, ``non_orthogonal_correction``). + +The explicit matrix is what an exact discrete adjoint needs: running a solver on +adjoint cross sections (:func:`make_adjoint_materials`) transposes the scattering +and fission but not ``C``, because the least-squares gradient correction is not +symmetric, so on a non-orthogonal mesh it solves ``(M + C)^T - C`` rather than +``M^T``. On an orthogonal mesh ``C`` vanishes and the two agree. + +Unknowns are ordered like the solver flux, ``index = cell * G + g``. Periodic +meshes are not supported. Requires SciPy. +""" + +from collections import namedtuple + +import numpy as np + +UnstructuredOperator = namedtuple("UnstructuredOperator", + ["loss", "production", "correction", "area"]) +UnstructuredOperator.__doc__ = """Sparse unstructured FVM operator. + +``loss`` is ``M`` (CSR, volume-integrated, correction included), ``production`` is +``F``, ``correction`` is the non-orthogonal part ``C`` alone (``M = M_imp - C``), +and ``area`` holds the cell areas.""" + + +def _faces(mesh): + """Interior faces (c0, c1, v_lo, v_hi) and boundary faces (cell, v_lo, v_hi).""" + off = np.asarray(mesh.cell_offsets, dtype=np.int64) + cv = np.asarray(mesh.cell_vertices, dtype=np.int64) + n = len(off) - 1 + cell = np.repeat(np.arange(n), np.diff(off)) + nxt = np.arange(len(cv)) + 1 + nxt[off[1:] - 1] = off[:-1] + va, vb = cv, cv[nxt] + lo, hi = np.minimum(va, vb), np.maximum(va, vb) + key = lo * (int(cv.max()) + 1) + hi + order = np.argsort(key, kind="stable") + ks = key[order] + first = np.r_[True, ks[1:] != ks[:-1]] + start = np.flatnonzero(first) + count = np.diff(np.r_[start, len(ks)]) + if count.max() > 2: + raise ValueError("mesh has an edge shared by more than two cells") + pair = count == 2 + i0, i1 = order[start[pair]], order[start[pair] + 1] + ib = order[start[~pair]] + return n, (cell[i0], cell[i1], lo[i0], hi[i0]), (cell[ib], lo[ib], hi[ib]) + + +def unstructured_operator(mats, mesh, bc): + """Assemble the sparse operator of the unstructured FVM solvers. + + Parameters + ---------- + mats : Materials + mesh : UnstructuredMesh2D + bc : list of BoundaryCondition + Tag-indexed, ``bc[tag * G + g]``, as passed to the solver. + + Returns + ------- + UnstructuredOperator + ``loss`` (M), ``production`` (F), ``correction`` (C), ``area``. A + fixed-source solution satisfies ``loss @ phi.ravel() = (q * area[:, None]).ravel()`` + and a k-eigenvalue one ``loss @ phi.ravel() = production @ phi.ravel() / k``, + both to the solver's tolerance. + """ + import scipy.sparse as sp + + from ndiffusion._core import cell_areas, cell_centroids + from ndiffusion.nearby import _fission_matrix_mode, _xs_arrays + + if len(mesh.periodic_a0): + raise ValueError("unstructured_operator does not support periodic meshes") + G = mats.n_groups + D, removal, chi, nusigf, scatter = _xs_arrays(mats) + mat = np.asarray(mesh.material_id, dtype=np.int64) + vx, vy = np.asarray(mesh.vx, dtype=float), np.asarray(mesh.vy, dtype=float) + cx, cy = (np.asarray(a, dtype=float) for a in cell_centroids(mesh)) + area = np.asarray(cell_areas(mesh), dtype=float) + n, (c0, c1, a, b), (bcell, ba, bb) = _faces(mesh) + + # Boundary tags (an edge missing from the boundary list gets tag 0, as in the solver). + bv0, bv1 = np.asarray(mesh.bface_v0), np.asarray(mesh.bface_v1) + tags = np.asarray(mesh.bface_bc_tag) if len(mesh.bface_bc_tag) else np.zeros(len(bv0), int) + tagmap = {(min(p, q), max(p, q)): t + for p, q, t in zip(bv0.tolist(), bv1.tolist(), tags.tolist())} + btag = np.array([tagmap.get((p, q), 0) for p, q in zip(ba.tolist(), bb.tolist())], + dtype=np.int64) + n_tags = len(bc) // G + A_bc = np.zeros((len(btag), G)) + B_bc = np.zeros((len(btag), G)) + valid = (btag >= 0) & (btag < n_tags) + for g in range(G): + A_bc[valid, g] = [bc[t * G + g].A for t in btag[valid]] + B_bc[valid, g] = [bc[t * G + g].B for t in btag[valid]] + + # Interior faces: surface vector, over-relaxed split S = E + T, interpolation weight. + ex, ey = cx[c1] - cx[c0], cy[c1] - cy[c0] + sx, sy = vy[b] - vy[a], -(vx[b] - vx[a]) + flip = sx * ex + sy * ey < 0.0 + sx[flip] *= -1.0 + sy[flip] *= -1.0 + length = np.hypot(sx, sy) + dist = np.hypot(ex, ey) + eux, euy = ex / dist, ey / dist + e_dot_s = eux * sx + euy * sy + s_dot_s = sx * sx + sy * sy + tangential = e_dot_s <= 1e-12 * np.sqrt(s_dot_s) + e_mag = np.where(tangential, 0.0, s_dot_s / np.where(tangential, 1.0, e_dot_s)) + a_coef = np.where(tangential, length / dist, e_mag / dist) + tx = np.where(tangential, 0.0, sx - e_mag * eux) + ty = np.where(tangential, 0.0, sy - e_mag * euy) + mx, my = 0.5 * (vx[a] + vx[b]), 0.5 * (vy[a] + vy[b]) + d0 = np.hypot(cx[c0] - mx, cy[c0] - my) + d1 = np.hypot(cx[c1] - mx, cy[c1] - my) + w0 = np.where(d0 + d1 > 0.0, d1 / np.where(d0 + d1 > 0.0, d0 + d1, 1.0), 0.5) + + # Boundary faces: normal distance from the centroid to the face. + bsx, bsy = vy[bb] - vy[ba], -(vx[bb] - vx[ba]) + bex = 0.5 * (vx[ba] + vx[bb]) - cx[bcell] + bey = 0.5 * (vy[ba] + vy[bb]) - cy[bcell] + bflip = bsx * bex + bsy * bey < 0.0 + bsx[bflip] *= -1.0 + bsy[bflip] *= -1.0 + blen = np.hypot(bsx, bsy) + bdn = (bex * bsx + bey * bsy) / blen + pos = bdn > 0.0 + bdn_safe = np.where(pos, bdn, 1.0) + b_acoef = np.where(pos, blen / bdn_safe, 0.0) + + # Weighted least-squares gradient coefficients (boundary samples at the face midpoint). + odx = np.r_[ex, -ex, bsx / blen * bdn] + ody = np.r_[ey, -ey, bsy / blen * bdn] + oc = np.r_[c0, c1, bcell] + w = 1.0 / (odx**2 + ody**2) + Mxx = np.bincount(oc, w * odx * odx, n) + Mxy = np.bincount(oc, w * odx * ody, n) + Myy = np.bincount(oc, w * ody * ody, n) + det = Mxx * Myy - Mxy**2 + ok = np.abs(det) > 1e-14 * (Mxx + Myy) ** 2 + det = np.where(ok, det, 1.0) + ixx, ixy, iyy = (np.where(ok, v, 0.0) for v in (Myy / det, -Mxy / det, Mxx / det)) + lx = ixx[oc] * (w * odx) + ixy[oc] * (w * ody) + ly = ixy[oc] * (w * odx) + iyy[oc] * (w * ody) + nI = len(c0) + l0x, l1x, lbx = lx[:nI], lx[nI:2 * nI], lx[2 * nI:] + l0y, l1y, lby = ly[:nI], ly[nI:2 * nI], ly[2 * nI:] + + cells = np.arange(n) + rows, cols, vals = [], [], [] + crows, ccols, cvals = [], [], [] + nonorth = np.flatnonzero(np.hypot(tx, ty) > 1e-10 * length) + nf = len(nonorth) + for g in range(G): + Dg = D[mat, g] + a_ij = 2.0 * Dg[c0] * Dg[c1] / (Dg[c0] + Dg[c1]) * a_coef + diag = np.bincount(c0, a_ij, n) + np.bincount(c1, a_ij, n) + A, B = A_bc[:, g], B_bc[:, g] + den = A + np.where(pos, B / bdn_safe, 0.0) + good = valid & (np.abs(den) > 1e-30) + den_safe = np.where(good, den, 1.0) + diag += np.bincount(bcell, np.where(good, Dg[bcell] * b_acoef * A / den_safe, 0.0), n) + diag += removal[mat, g] * area + rows += [cells * G + g, c0 * G + g, c1 * G + g] + cols += [cells * G + g, c1 * G + g, c0 * G + g] + vals += [diag, -a_ij, -a_ij] + for gp in range(G): + if gp != g: + rows.append(cells * G + g) + cols.append(cells * G + gp) + vals.append(-scatter[mat, g, gp] * area) + if nf == 0: + continue + # Cell gradients as n x n operators, grad_c = sum_f lsq_f (phi_nbr - phi_c); a boundary + # face samples the Robin surface value, phi_s - phi_c = delta phi_c. + delta = np.where(good, -A / den_safe, 0.0) + gr = np.r_[c0, c0, c1, c1, bcell] + gc = np.r_[c1, c0, c0, c1, bcell] + Gx = sp.csr_matrix((np.r_[l0x, -l0x, l1x, -l1x, lbx * delta], (gr, gc)), shape=(n, n)) + Gy = sp.csr_matrix((np.r_[l0y, -l0y, l1y, -l1y, lby * delta], (gr, gc)), shape=(n, n)) + f = nonorth + face = np.r_[np.arange(nf), np.arange(nf)] + W = sp.csr_matrix((np.r_[w0[f], 1.0 - w0[f]], (face, np.r_[c0[f], c1[f]])), shape=(nf, n)) + Df = 2.0 * Dg[c0[f]] * Dg[c1[f]] / (Dg[c0[f]] + Dg[c1[f]]) + flux = sp.diags(Df * tx[f]) @ (W @ Gx) + sp.diags(Df * ty[f]) @ (W @ Gy) + S = sp.csr_matrix((np.r_[np.ones(nf), -np.ones(nf)], (np.r_[c0[f], c1[f]], face)), + shape=(n, nf)) + Cg = (S @ flux).tocoo() + crows.append(Cg.row * G + g) + ccols.append(Cg.col * G + g) + cvals.append(Cg.data) + + N = n * G + cat = lambda parts: np.concatenate(parts) if parts else np.zeros(0) # noqa: E731 + C = sp.csr_matrix((cat(cvals), (cat(crows).astype(np.int64), cat(ccols).astype(np.int64))), + shape=(N, N)) + M = sp.csr_matrix((cat(vals), (cat(rows), cat(cols))), shape=(N, N)) - C + + frows, fcols, fvals = [], [], [] + matrix_mode = _fission_matrix_mode(mats) + for g in range(G): + for gp in range(G): + coef = nusigf[mat, g, gp] if matrix_mode else chi[mat, g] * nusigf[mat, gp] + frows.append(cells * G + g) + fcols.append(cells * G + gp) + fvals.append(coef * area) + F = sp.csr_matrix((cat(fvals), (cat(frows), cat(fcols))), shape=(N, N)) + F.eliminate_zeros() + return UnstructuredOperator(M.tocsr(), F, C, area) diff --git a/src/ndiffusion/nearby.py b/src/ndiffusion/nearby.py index 53522bf..f27c40c 100644 --- a/src/ndiffusion/nearby.py +++ b/src/ndiffusion/nearby.py @@ -69,8 +69,15 @@ "nearby_rate", "k_error_adjoint", "k_error_cells", + "fit", + "operator", ], - defaults=(None, None), + defaults=(None, None, None, None), +) + +NearbyQoIResult = namedtuple( + "NearbyQoIResult", + ["value", "fit_value", "error_adjoint", "error_nearby", "cells"], ) @@ -487,12 +494,15 @@ def _curvelap_spline(mats, mesh, flux, bc, k, source): config = "auto" if source is None else "corner" # only the corner rows carry a source f = _spline_fit.fit(mats, mesh, flux, k, bc, stride=_spline_stride(mats, area), config=config, source=source) - return f.curve, f.lap, area, np.asarray(mesh.material_id, dtype=int) + return f.curve, f.lap, area, np.asarray(mesh.material_id, dtype=int), f def _fit_dispatch(mats, flux, medium_map, edges_x, edges_y, geometry, mesh, fit=None, bc=None, k=np.inf, source=None): - """Pick the curve fitter from the supplied geometry descriptor.""" + """Pick the curve fitter from the supplied geometry descriptor. + + Returns (curve, lap, volumes, material ids, fit), where fit is the + ``SplineFit`` for the constrained spline and None otherwise.""" flux = np.asarray(flux, dtype=float) if fit not in (None, "spline", "local"): raise ValueError(f"fit must be None, 'spline' or 'local', got {fit!r}") @@ -502,13 +512,13 @@ def _fit_dispatch(mats, flux, medium_map, edges_x, edges_y, geometry, mesh, raise ValueError("fit='spline' needs bc=, the tag-indexed boundary " "conditions the solver was built with") return _curvelap_spline(mats, mesh, flux, bc, k, source) - return _curvelap_unstructured(mats, mesh, flux) + return (*_curvelap_unstructured(mats, mesh, flux), None) if fit == "spline": raise ValueError("fit='spline' is for the unstructured mesh (mesh=)") if edges_y is not None: - return _curvelap_2d(mats, medium_map, edges_x, edges_y, geometry, flux) + return (*_curvelap_2d(mats, medium_map, edges_x, edges_y, geometry, flux), None) if edges_x is not None: - return _curvelap_1d(mats, medium_map, edges_x, geometry, flux) + return (*_curvelap_1d(mats, medium_map, edges_x, geometry, flux), None) raise ValueError( "Provide a geometry descriptor: (medium_map, edges_x[, edges_y], " "geometry) for structured meshes, or mesh= for the unstructured mesh." @@ -575,7 +585,7 @@ def nearby_fixed_source( numerical = solver.solve(source) flux = np.asarray(numerical.flux, dtype=float) - curve, lap, _vol, mat_ids = _fit_dispatch( + curve, lap, _vol, mat_ids, _fit = _fit_dispatch( mats, flux, medium_map, edges_x, edges_y, geometry, mesh, fit=fit, bc=bc, k=np.inf, source=source, ) @@ -604,21 +614,24 @@ def nearby_k_eigenvalue( fit=None, adjoint_solver=None, nearby=True, + method="iterative", tol=1e-8, max_iter=200, ): """Run the method of nearby problems for a k-eigenvalue diffusion solve. - The nearby problem is solved by a Python power iteration that reuses - ``fixed_solver`` (which inverts the same loss+scatter operator ``M`` as the - k-eigenvalue solver): each outer step solves ``M phi = (1/k) F phi + r`` and - updates ``k`` from the production integral ratio. + The nearby problem ``M phi = (1/k) F phi + r`` is solved by a power + iteration. With ``method="iterative"`` each outer step is a + ``fixed_solver`` solve (it inverts the same loss+scatter operator ``M`` as + the k-eigenvalue solver). With ``method="direct"`` (unstructured mesh and + ``bc=`` only) ``M`` is assembled by :func:`unstructured_operator` and + factorized once, which is much faster, and ``fixed_solver`` may be None. Parameters ---------- keig_solver : KEigenSolver / KEigenSolver2D / KEigenSolverUnstructured2D Solver whose ``.solve()`` returns a ``DiffusionResult`` (flux + keff). - fixed_solver : matching FixedSource* solver + fixed_solver : matching FixedSource* solver, or None with ``method="direct"`` Built with the *same* materials, mesh, and boundary conditions. mats : Materials medium_map, edges_x, edges_y, geometry, mesh @@ -629,11 +642,14 @@ def nearby_k_eigenvalue( Built on :func:`make_adjoint_materials` of ``mats`` (same mesh and boundary conditions). Its flux gives the one-shot eigenvalue error estimate ``k_error_adjoint`` and its per-cell contributions - ``k_error_cells``. + ``k_error_cells``. With ``method="direct"`` the adjoint is instead the + exact discrete one, from the assembled operator, and this is ignored. nearby : bool, default True - Run the nearby power iteration. With ``nearby=False`` (and an - ``adjoint_solver``) only the fit and the one-shot estimate are formed; - ``nearby_flux`` and ``k_nearby`` are then None. + Run the nearby power iteration. With ``nearby=False`` only the fit and + the one-shot estimate are formed; ``nearby_flux`` and ``k_nearby`` are + then None. + method : {"iterative", "direct"} + How the nearby problem and the adjoint are solved (see above). tol : float Convergence tolerance on the flux change and ``|dk|``. max_iter : int @@ -647,16 +663,27 @@ def nearby_k_eigenvalue( curve-fit fission production integral). The arrays are shaped ``(n_cells, G)``, like the solver flux. ``k_nearby - k_curve_fit`` estimates the eigenvalue discretization error, - and so does ``k_error_adjoint`` when an ``adjoint_solver`` is given + and so does ``k_error_adjoint`` when an adjoint is available (``k_error_cells`` sums to it). On layouts whose error contributions cancel - material corners, small inclusions - ``k_error_cells`` and the local error are more reliable than the net eigenvalue estimate; the sum - of ``abs(k_error_cells)`` is a (loose) bound. + of ``abs(k_error_cells)`` is a (loose) bound. ``fit`` is the + constrained spline fit (``SplineFit``) when one was used and + ``operator`` the assembled operator with ``method="direct"``; both feed + :func:`nearby_qoi`. """ + if method not in ("iterative", "direct"): + raise ValueError(f"method must be 'iterative' or 'direct', got {method!r}") + if method == "direct" and (mesh is None or bc is None): + raise ValueError("method='direct' needs the unstructured mesh= and bc=") + if method == "iterative" and nearby and fixed_solver is None: + raise ValueError("the iterative nearby solve needs fixed_solver") + numerical = keig_solver.solve() flux = np.asarray(numerical.flux, dtype=float) + G = mats.n_groups - curve, lap, vol, mat_ids = _fit_dispatch( + curve, lap, vol, mat_ids, spline = _fit_dispatch( mats, flux, medium_map, edges_x, edges_y, geometry, mesh, fit=fit, bc=bc, k=numerical.keff, ) @@ -673,23 +700,43 @@ def nearby_k_eigenvalue( residual = loss - prod / k_cf # (n_cells, G), exact source of the nearby problem + op = lu = None + if method == "direct": + from scipy.sparse.linalg import splu + + from ndiffusion.fv_operator import unstructured_operator + + op = unstructured_operator(mats, mesh, bc) + lu = splu(op.loss.tocsc()) + k_error = cells = None - if adjoint_solver is not None: - adj = np.asarray(adjoint_solver.solve().flux, dtype=float).reshape(-1, mats.n_groups) + adj = None + if method == "direct": + adj = _discrete_adjoint(op, flux, numerical.keff) + elif adjoint_solver is not None: + adj = np.asarray(adjoint_solver.solve().flux, dtype=float).reshape(-1, G) + if adj is not None: cells = (np.sum(adj * (numerical.keff * loss - prod) * volG, axis=1) / float(np.sum(adj * loss * volG))) k_error = float(cells.sum()) if not nearby: - return NearbyKResult(numerical, curve, k_cf, None, None, residual, rate, k_error, cells) + return NearbyKResult(numerical, curve, k_cf, None, None, residual, rate, k_error, cells, + spline, op) + + # --- Nearby power iteration: M phi = F phi / k + r --- + if lu is not None: + def solve_m(q): + return lu.solve((q * volG).ravel()).reshape(-1, G) + else: + def solve_m(q): + return np.asarray(fixed_solver.solve(q).flux, dtype=float) - # --- Nearby power iteration (reuses the fixed-source solver's operator M) --- phi = curve.copy() k = k_cf for _ in range(max_iter): Fphi = _fission_apply(mats, mat_ids, phi) prod_old = float(np.sum(Fphi * volG)) - q = Fphi / k + residual - phi_new = fixed_solver.solve(q).flux + phi_new = solve_m(Fphi / k + residual) prod_new = float(np.sum(_fission_apply(mats, mat_ids, phi_new) * volG)) k_new = k * prod_new / prod_old if prod_old != 0.0 else k @@ -712,4 +759,158 @@ def nearby_k_eigenvalue( nearby_rate=rate, k_error_adjoint=k_error, k_error_cells=cells, + fit=spline, + operator=op, ) + + +def _bordered_adjoint_lu(op, flux, keff): + """LU of [[(M - F/k)^T, F phi], [(F phi)^T, 0]], the generalized adjoint system. + + The appended row and column make the singular adjoint operator invertible: + the column absorbs the small inconsistency of an iteratively converged + (phi, k), and the row fixes the free multiple of the adjoint eigenvector.""" + import scipy.sparse as sp + from scipy.sparse.linalg import splu + + Fphi = op.production @ np.asarray(flux, dtype=float).ravel() + col = sp.csc_matrix(Fphi[:, None]) + A = (op.loss - op.production / keff).T + return splu(sp.bmat([[A, col], [col.T, None]], format="csc")), Fphi + + +def _discrete_adjoint(op, flux, keff): + """The exact discrete adjoint eigenvector, (n_cells, G), normalized to = 1.""" + lu, Fphi = _bordered_adjoint_lu(op, flux, keff) + N = len(Fphi) + rhs = np.zeros(N + 1) + rhs[N] = 1.0 + # (M - F/k)^T psi + mu F phi = 0 with = 1: mu vanishes and psi is the adjoint. + return lu.solve(rhs)[:N].reshape(np.shape(flux)) + + +# Dunavant degree-4 rule on a triangle: (barycentric a, weight) orbits. +_D4 = ((0.445948490915965, 0.223381589678011), (0.091576213509771, 0.109951743655322)) + + +def _cell_quadrature(mesh, cx, cy): + """Degree-4 quadrature on every cell, by fanning it into triangles from its + centroid. Returns (cell, x, y, weight) arrays; the weights sum to the areas.""" + off = np.asarray(mesh.cell_offsets, dtype=np.int64) + cv = np.asarray(mesh.cell_vertices, dtype=np.int64) + vx, vy = np.asarray(mesh.vx, dtype=float), np.asarray(mesh.vy, dtype=float) + cell = np.repeat(np.arange(len(off) - 1), np.diff(off)) + nxt = np.arange(len(cv)) + 1 + nxt[off[1:] - 1] = off[:-1] + P0 = np.column_stack([cx[cell], cy[cell]]) + P1 = np.column_stack([vx[cv], vy[cv]]) + P2 = np.column_stack([vx[cv[nxt]], vy[cv[nxt]]]) + tri = 0.5 * np.abs((P1[:, 0] - P0[:, 0]) * (P2[:, 1] - P0[:, 1]) + - (P2[:, 0] - P0[:, 0]) * (P1[:, 1] - P0[:, 1])) + cells, xs, ys, ws = [], [], [], [] + for s_, w in _D4: + for lam in ((s_, s_, 1 - 2 * s_), (s_, 1 - 2 * s_, s_), (1 - 2 * s_, s_, s_)): + p = lam[0] * P0 + lam[1] * P1 + lam[2] * P2 + cells.append(cell) + xs.append(p[:, 0]) + ys.append(p[:, 1]) + ws.append(w * tri) + return np.concatenate(cells), np.concatenate(xs), np.concatenate(ys), np.concatenate(ws) + + +def _ratio(num, den, phi, vol): + """J = sum num phi V / sum den phi V for a stack of weights (n_q, n, G).""" + top = np.einsum("qcg,cg,c->q", num, phi, vol) + bot = np.einsum("qcg,cg,c->q", den, phi, vol) + return top / bot + + +def nearby_qoi(result, mats, mesh, bc, numerator, denominator): + """Discretization error of normalized quantities of interest. + + A quantity of interest here is a ratio of two reaction-rate integrals, + + J(phi) = int a . phi dV / int b . phi dV, + + with per-unit-volume weights ``a`` (numerator) and ``b`` (denominator) per + cell and group, e.g. the absorption cross section in one region over the + production cross section everywhere (a normalized region power or + absorption rate). The ratio makes ``J`` independent of the eigenvector's + arbitrary normalization. + + The one-shot estimate solves the generalized adjoint problem + ``(M - F/k)^T Gamma = dJ/dphi`` with `` = 0`` on the exact + discrete operator (:func:`unstructured_operator`) and returns + `` + J(phi_h) - J(fit)``, the linearization of the nearby + estimate ``J(phi_nearby) - J(fit)``. It needs the exact transpose: a solver + run on adjoint cross sections does not transpose the non-orthogonal + correction and leaves a bias of 1-2% on triangles and hexagons. + + Parameters + ---------- + result : NearbyKResult + From :func:`nearby_k_eigenvalue` with the constrained spline fit + (``bc=`` given). Its ``operator`` is reused when present + (``method="direct"``). + mats, mesh, bc + As passed to :func:`nearby_k_eigenvalue`. + numerator, denominator : array-like + Weights ``a`` and ``b``, shaped ``(n_cells, G)``, or ``(n_q, n_cells, G)`` + to evaluate several quantities with one factorization. + + Returns + ------- + NearbyQoIResult + ``value`` J(phi_h), ``fit_value`` J(fit) (by degree-4 quadrature of the + fit over each cell), ``error_adjoint`` the one-shot estimate of + ``J(phi_h) - J(phi_exact)``, ``error_nearby`` the nearby-solve estimate + ``J(phi_nearby) - J(fit)`` (None without a nearby solve), and ``cells``, + the per-cell contributions ``Gamma . V r`` (n_q, n_cells), a + goal-oriented error map for adaptation. Arrays have a leading ``n_q`` + axis only when the weights do. + """ + from ndiffusion._core import cell_centroids + from ndiffusion.fv_operator import unstructured_operator + + if result.fit is None: + raise ValueError("nearby_qoi needs a result from the constrained spline fit (bc= given)") + G = mats.n_groups + phi = np.asarray(result.numerical.flux, dtype=float).reshape(-1, G) + n = phi.shape[0] + num = np.asarray(numerator, dtype=float) + den = np.asarray(denominator, dtype=float) + single = num.ndim == 2 + num = num.reshape(-1, n, G) + den = np.broadcast_to(den.reshape(-1, n, G), num.shape) + + op = result.operator if result.operator is not None else unstructured_operator(mats, mesh, bc) + area = op.area + keff = result.numerical.keff + J_h = _ratio(num, den, phi, area) + + cx, cy = (np.asarray(a, dtype=float) for a in cell_centroids(mesh)) + qc, qx, qy, qw = _cell_quadrature(mesh, cx, cy) + fvals = result.fit.evaluate(qc, qx, qy) + J_fit = (np.einsum("qkg,kg,k->q", num[:, qc], fvals, qw) + / np.einsum("qkg,kg,k->q", den[:, qc], fvals, qw)) + + # dJ/dphi in the volume-integrated layout: (a V - J b V) / int b phi dV + bot = np.einsum("qcg,cg,c->q", den, phi, area) + grads = (num - J_h[:, None, None] * den) * area[None, :, None] / bot[:, None, None] + lu, _ = _bordered_adjoint_lu(op, phi, keff) + N = n * G + rhs = np.zeros((N + 1, len(J_h))) + rhs[:N] = grads.reshape(len(J_h), N).T + Gam = lu.solve(rhs)[:N].T.reshape(len(J_h), n, G) + Vr = np.asarray(result.residual, dtype=float).reshape(n, G) * area[:, None] + cells = np.einsum("qcg,cg->qc", Gam, Vr) + est = cells.sum(axis=1) + J_h - J_fit + near = None + if result.nearby_flux is not None: + phi_n = np.asarray(result.nearby_flux, dtype=float).reshape(n, G) + near = _ratio(num, den, phi_n, area) - J_fit + + if single: + return NearbyQoIResult(float(J_h[0]), float(J_fit[0]), float(est[0]), + None if near is None else float(near[0]), cells[0]) + return NearbyQoIResult(J_h, J_fit, est, near, cells) diff --git a/tests/test_fv_operator.py b/tests/test_fv_operator.py new file mode 100644 index 0000000..de52454 --- /dev/null +++ b/tests/test_fv_operator.py @@ -0,0 +1,211 @@ +"""The assembled unstructured operator reproduces the matrix-free solvers.""" + +import numpy as np +import pytest + +import ndiffusion as nd +from ndiffusion.materials import two_group + +pytest.importorskip("scipy") + +from test_nearby_unstructured import quad_mesh, triangle_mesh # noqa: E402 + +# fuel and reflector, IAEA-like two-group constants +ROWS = [(1.5, 0.4, 0.010, 0.080, 0.020, 0.0, 0.135), + (2.0, 0.3, 0.000, 0.010, 0.040, 0.0, 0.0)] + + +def problem(kind, n=12, L=20.0, jitter=0.0, seed=0): + """Fuel in x < L/2, reflector beyond; reflective on x = 0 and y = 0, vacuum-like + Robin elsewhere. ``jitter`` moves interior vertices off the interface by up to that + fraction of a cell (n must be even so the interface falls on vertices).""" + mesh = (triangle_mesh if kind == "tri" else quad_mesh)(n, n, L, L) + vx, vy = np.asarray(mesh.vx, float), np.asarray(mesh.vy, float) + if jitter: + rng = np.random.default_rng(seed) + inner = (vx > 1e-9) & (vx < L - 1e-9) & (vy > 1e-9) & (vy < L - 1e-9) + inner &= np.abs(vx - L / 2) > 1e-9 # keep the material interface straight + h = L / n + vx[inner] += jitter * h * rng.uniform(-1, 1, inner.sum()) + vy[inner] += jitter * h * rng.uniform(-1, 1, inner.sum()) + mesh.vx, mesh.vy = vx, vy + far = L - 1e-9 + mesh.bface_bc_tag = [int(min(vx[a], vx[b]) > far or min(vy[a], vy[b]) > far) + for a, b in zip(mesh.bface_v0, mesh.bface_v1)] + cx, _cy = (np.asarray(a) for a in nd.cell_centroids(mesh)) + mesh.material_id = (cx > L / 2).astype(int) + mats = two_group(ROWS) + bc = [nd.BoundaryCondition(A=0.0, B=1.0), nd.BoundaryCondition(A=0.0, B=1.0), + nd.BoundaryCondition(A=0.4692, B=1.5), nd.BoundaryCondition(A=0.4692, B=0.4)] + return mats, mesh, bc + + +def keig(mats, mesh, bc): + return nd.KEigenSolverUnstructured2D(mats, mesh, bc, epsilon=1e-12, max_outer=20000, + max_inner=400, use_cg=True).solve() + + +class TestUnstructuredOperator: + @pytest.mark.parametrize("kind,jitter", [("quad", 0.0), ("tri", 0.0), ("tri", 0.2), ("quad", 0.2)]) + def test_fixed_source_solution_satisfies_the_matrix(self, kind, jitter): + mats, mesh, bc = problem(kind, jitter=jitter) + op = nd.unstructured_operator(mats, mesh, bc) + q = np.random.default_rng(1).random((len(op.area), 2)) + fs = nd.FixedSourceSolverUnstructured2D(mats, mesh, bc, epsilon=1e-13, max_inner=100000, omega=1.5) + phi = np.asarray(fs.solve(q).flux) + rhs = (q * op.area[:, None]).ravel() + assert np.linalg.norm(op.loss @ phi.ravel() - rhs) / np.linalg.norm(rhs) < 1e-10 + + def test_correction_vanishes_only_on_orthogonal_meshes(self): + assert nd.unstructured_operator(*problem("quad")).correction.nnz == 0 + assert nd.unstructured_operator(*problem("tri")).correction.nnz > 0 + + def test_eigenvalue_solution_satisfies_the_matrix(self): + mats, mesh, bc = problem("tri", jitter=0.2) + op = nd.unstructured_operator(mats, mesh, bc) + r = keig(mats, mesh, bc) + phi = np.asarray(r.flux).ravel() + Fphi = op.production @ phi / r.keff + assert np.linalg.norm(op.loss @ phi - Fphi) / np.linalg.norm(Fphi) < 1e-8 + + def test_adjoint_materials_do_not_transpose_the_correction(self): + mats, mesh, bc = problem("tri", jitter=0.2) + op = nd.unstructured_operator(mats, mesh, bc) + q = np.random.default_rng(2).random((len(op.area), 2)) + fs = nd.FixedSourceSolverUnstructured2D(nd.make_adjoint_materials(mats), mesh, bc, + epsilon=1e-13, max_inner=100000, omega=1.5) + x = np.asarray(fs.solve(q).flux).ravel() + rhs = (q * op.area[:, None]).ravel() + solver_op = (op.loss + op.correction).T - op.correction + assert np.linalg.norm(solver_op @ x - rhs) / np.linalg.norm(rhs) < 1e-10 + assert np.linalg.norm(op.loss.T @ x - rhs) / np.linalg.norm(rhs) > 1e-3 + + def test_fission_matrix_mode(self): + mats, mesh, bc = problem("quad", n=4) + F = nd.unstructured_operator(mats, mesh, bc).production + chi = np.asarray(mats.chi).reshape(-1, 2) + nsf = np.asarray(mats.nusigf).reshape(-1, 2) + mats.nusigf = np.einsum("mg,mh->mgh", chi, nsf).ravel() + mats.chi = np.zeros_like(chi).ravel() + assert abs(nd.unstructured_operator(mats, mesh, bc).production - F).max() < 1e-15 + + def test_periodic_mesh_is_rejected(self): + mats, mesh, bc = problem("quad", n=4) + mesh.periodic_a0, mesh.periodic_a1 = [0], [1] + mesh.periodic_b0, mesh.periodic_b1 = [20], [21] + with pytest.raises(ValueError, match="periodic"): + nd.unstructured_operator(mats, mesh, bc) + + +def run(kind="tri", method="direct", nearby=True, n=30, **kw): + """Resolved enough for the fit: knots about 2 cm apart, five spans across each region.""" + mats, mesh, bc = problem(kind, n=n, jitter=0.2) + ks = nd.KEigenSolverUnstructured2D(mats, mesh, bc, epsilon=1e-12, max_outer=20000, + max_inner=400, use_cg=True) + fixed = None + if method == "iterative": + fixed = nd.FixedSourceSolverUnstructured2D(mats, mesh, bc, epsilon=1e-13, + max_inner=100000, omega=1.5) + r = nd.nearby_k_eigenvalue(ks, fixed, mats, mesh=mesh, bc=bc, method=method, nearby=nearby, + tol=1e-12, max_iter=20000, **kw) + return mats, mesh, bc, r + + +def region_weights(mats, mesh): + """Absorption in the fuel's outer half (x > 5) over total production.""" + G = 2 + mat = np.asarray(mesh.material_id) + cx, _ = (np.asarray(a) for a in nd.cell_centroids(mesh)) + rem = np.asarray(mats.removal).reshape(-1, G) + sc = np.asarray(mats.scatter).reshape(-1, G, G) + sig_a = rem - (sc.sum(axis=1) - np.einsum("mgg->mg", sc)) + nsf = np.asarray(mats.nusigf).reshape(-1, G) + num = sig_a[mat] * ((mat == 0) & (cx > 5.0))[:, None] + return num, nsf[mat] + + +class TestNearbyDirect: + def test_direct_and_iterative_nearby_solves_agree(self): + _, _, _, d = run(method="direct") + _, _, _, i = run(method="iterative") + assert d.k_nearby == pytest.approx(i.k_nearby, abs=1e-9) + assert np.abs(d.nearby_flux - i.nearby_flux).max() < 1e-6 * np.abs(i.nearby_flux).max() + assert d.operator is not None and i.operator is None + + def test_exact_adjoint_matches_adjoint_materials_on_quads(self): + mats, mesh, bc = problem("quad") + ks = nd.KEigenSolverUnstructured2D(mats, mesh, bc, epsilon=1e-12, max_outer=20000, + max_inner=400, use_cg=True) + adj = nd.KEigenSolverUnstructured2D(nd.make_adjoint_materials(mats), mesh, bc, + epsilon=1e-12, max_outer=20000, max_inner=400, use_cg=True) + d = nd.nearby_k_eigenvalue(ks, None, mats, mesh=mesh, bc=bc, method="direct", nearby=False) + i = nd.nearby_k_eigenvalue(ks, None, mats, mesh=mesh, bc=bc, adjoint_solver=adj, nearby=False) + assert d.k_error_adjoint == pytest.approx(i.k_error_adjoint, rel=1e-6) + + def test_direct_needs_mesh_and_bc(self): + mats, mesh, bc = problem("quad", n=4) + ks = nd.KEigenSolverUnstructured2D(mats, mesh, bc, epsilon=1e-10, use_cg=True) + with pytest.raises(ValueError, match="direct"): + nd.nearby_k_eigenvalue(ks, None, mats, mesh=mesh, method="direct") + + +class TestNearbyQoI: + def test_one_shot_is_the_exact_linearization(self): + """ equals dJ/d(eps) of the problem with source eps V r, by central + differences, on a non-orthogonal mesh.""" + from scipy.sparse.linalg import splu + + mats, mesh, bc, r = run() + num, den = region_weights(mats, mesh) + q = nd.nearby_qoi(r, mats, mesh, bc, num, den) + linear = q.cells.sum() + + op = r.operator + lu = splu(op.loss.tocsc()) + phi0 = np.asarray(r.numerical.flux).ravel() + Vr = (r.residual * op.area[:, None]).ravel() + + def J(eps): + phi, k = phi0.copy(), r.numerical.keff + for _ in range(20000): + Fp = op.production @ phi + new = lu.solve(Fp / k + eps * Vr) + k_new = k * (op.production @ new).sum() / Fp.sum() + done = abs(k_new - k) < 1e-15 and np.abs(new - phi).max() < 1e-13 * np.abs(phi).max() + phi, k = new, k_new + if done: + break + p = phi.reshape(-1, 2) + return np.sum(num * p * op.area[:, None]) / np.sum(den * p * op.area[:, None]) + + eps = 0.01 + central = (J(eps) - J(-eps)) / (2 * eps) + assert linear == pytest.approx(central, rel=1e-4) + + def test_one_shot_tracks_the_nearby_estimate(self): + mats, mesh, bc, r = run() + num, den = region_weights(mats, mesh) + q = nd.nearby_qoi(r, mats, mesh, bc, num, den) + assert q.error_nearby == pytest.approx(q.error_adjoint, rel=0.02) + assert q.value == pytest.approx( + np.sum(num * r.numerical.flux * r.operator.area[:, None]) + / np.sum(den * r.numerical.flux * r.operator.area[:, None])) + + def test_stacked_weights(self): + mats, mesh, bc, r = run(nearby=False) + num, den = region_weights(mats, mesh) + one = nd.nearby_qoi(r, mats, mesh, bc, num, den) + both = nd.nearby_qoi(r, mats, mesh, bc, np.stack([num, 2 * num]), den) + assert both.error_adjoint.shape == (2,) and both.error_nearby is None + assert both.error_adjoint[0] == pytest.approx(one.error_adjoint) + assert both.error_adjoint[1] == pytest.approx(2 * one.error_adjoint) + assert both.cells.shape == (2, len(r.operator.area)) + + def test_needs_the_spline_fit(self): + mats, mesh, bc = problem("quad", n=6) + ks = nd.KEigenSolverUnstructured2D(mats, mesh, bc, epsilon=1e-10, use_cg=True) + fixed = nd.FixedSourceSolverUnstructured2D(mats, mesh, bc, epsilon=1e-12, omega=1.5) + r = nd.nearby_k_eigenvalue(ks, fixed, mats, mesh=mesh, fit="local", nearby=False) + num, den = region_weights(mats, mesh) + with pytest.raises(ValueError, match="spline"): + nd.nearby_qoi(r, mats, mesh, bc, num, den) From 302771d7f7656c1dabf925288e73b6e0d16c7f27 Mon Sep 17 00:00:00 2001 From: Ben Whewell Date: Sat, 17 Oct 2026 12:03:30 -0600 Subject: [PATCH 4/4] Import the shared test meshes relatively so the operator tests also collect from the repository root. --- tests/test_fv_operator.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_fv_operator.py b/tests/test_fv_operator.py index de52454..8713be4 100644 --- a/tests/test_fv_operator.py +++ b/tests/test_fv_operator.py @@ -8,7 +8,7 @@ pytest.importorskip("scipy") -from test_nearby_unstructured import quad_mesh, triangle_mesh # noqa: E402 +from .test_nearby_unstructured import quad_mesh, triangle_mesh # noqa: E402 # fuel and reflector, IAEA-like two-group constants ROWS = [(1.5, 0.4, 0.010, 0.080, 0.020, 0.0, 0.135),