Constrained spline fit, exact discrete adjoint and QoI estimates for unstructured nearby problems

CHANGELOG.md

@@ -6,8 +6,33 @@ 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.
+- `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
 
+- 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.
 - A k-eigenvalue solve with an `epsilon` tighter than round-off allows no longer
   runs every inner iteration out to `max_inner`. The inner sweeps now stop once
   their change stops shrinking. This mostly hit the CG inner solver on fine

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

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,89 @@ $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. 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
 
-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.

docs/user_guide/verification.md

@@ -69,10 +69,71 @@ 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.
+
+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
+eigenvalue estimate as an indicator of size, with `abs(r.k_error_cells).sum()`
+as a loose bound.
 
 ## Adjoint problems
 

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",