Energy-Dependent, Self-Adaptive Mesh h(p)-Refinement of a Constraint-Based Continuous Bubnov-Galerkin Isogeometric Analysis Spatial Discretization of the Multi-Group Neutron Diffusion Equation with Dual-Weighted Residual Error Measures

Abstract

Energy-dependent self-adaptive mesh refinement algorithms are developed for a continuous Bubnov-Galerkin spatial discretization of the multi-group neutron diffusion equation using NURBS-based isogeometric analysis (IGA). The spatially self-adaptive algorithms employ both mesh (h) and polynomial degree (p) refinement. Constraint-based equations are established across irregular interfaces with hanging-nodes; they are based upon master-slave relationships and the conservative interpolation between surface meshes. A similar Galerkin projection is employed in the conservative interpolation between volume meshes to evaluate group-to-group source terms over energy-dependent meshes; and to evaluate interpolation-based error measures. Enforcing continuity over an irregular mesh does introduce discretization errors. However, local mesh refinement allows for a better allocation of computational resources; and thus, more accuracy per degree of freedom. Two a posteriori interpolation-based error measures are proposed. The first heuristically minimizes local contributions to the discretization error, which becomes competitive for global quantities of interest (QoIs). However, for localized QoIs, over energy-dependent meshes, certain multi-group components may become under-resolved. The second employs duality arguments to minimize important error contributions, which consistently and reliably reduces the error in the QoI.

Keywords:

• Multi-group neutron diffusion equation
• continuous Bubnov-Galerkin
• isogeometric analysis
• energy-dependent self-adaptive mesh h(p)-refinement
• constraint-based local refinement
• dual-weighted residual error measures

1. Introduction

The modeling and simulation (M&S) of modern-generation nuclear power plants relies upon the solution of a complex transport equation (Cummins and Matzie 2018). This is the neutron transport equation (NTE), which describes the non-equilibrium averaged statistical behavior of neutrons within a host material. The NTE is a linearized form of the Boltzmann transport equation, which describes the non-equilibrium average statistical behavior of a gas of atoms, or molecules (Boltzmann 1872). The NTE can be computationally challenging to solve for large scale, general geometry, nuclear reactor physics and shielding problems. Therefore, a mathematical approximation to the NTE is often used in whole-core nuclear reactor physics simulations; this is the steady-state neutron diffusion equation (NDE), which is an elliptic partial-integro-differential equation (PIDE). The NDE is a function of spatial position (x, y, z); and energy (E).

There exist high-fidelity, geometry conforming, spatial discretization techniques that can be employed to solve the NTE, or NDE. The finite element method (FEM) was first applied to nuclear reactor physics problems in the 1970s (Semenza, Lewis, and Rossow 1972; Hansen and Kang 1973; Martin and Duderstadt 1977). The key advantage of FEM is the ability to calculate solution fields over curvilinear geometries (Semenza, Lewis, and Rossow 1972; Hansen and Kang 1973; Martin and Duderstadt 1977). The FEM consists of decomposing the problem domain into basic geometric primitives; e.g. triangles and quadrilaterals in 2D (Becker, Carey, and Oden 1981). The finite-element (FE) mesh is the union of the finite elements (FEs); and the geometry of each and the solution field within each are represented by a set of polynomial basis functions. The lowest-order FEM approximation employs linear FEs that represent curvilinear domains using straight-sided line segments in 2D. A number of different mathematical and computational FEM-based refinement strategies have been devised, ranging from heuristic-based approaches to the derivation of more-rigorous a posteriori error estimates (Wang and Ragusa 2009a; Wang, Bangerth, and Ragusa 2009; Wang and Ragusa 2011).

A uniform refinement strategy may be defined whereby the characteristic mesh spacing, h, is decreased across the whole mesh; i.e., $h\to 0;$ or the polynomial degree of the basis, p, is increased; i.e., $p\to \infty ;$ or both, in tandem. Such uniform refinement strategies are useful for verifying consistency in the numerical code since theoretical a priori rates of convergence are well established (Salari and Knupp 2000; Wang and Ragusa 2009b). However, they do not generally allow for an optimal allocation of computational resources. In the pursuit of numerical accuracy and computational efficiency, adaptive mesh refinement (AMR) strategies were developed. These AMR algorithms improve numerical accuracy per degree of freedom (DoF) by refining the computational mesh where necessary (Ainsworth and Oden 1997).

Adaptive h-refinement techniques were developed to rectify automatically any deficiencies in the coarsest mesh; for example, to eliminate any spurious numerical artifacts about discontinuities, or sharp gradients, in the solution field that result from any discontinuous coefficients, boundary data or source terms (Bonito, Devore, and Nochetto 2013; Mitchell and McClain 2014). Both uniform and local h-refinement strategies can attain only algebraic rates of convergence, which are restricted by the polynomial degree of the basis, p. The convergence of the former is also limited by the regularity of the solution, whereas that of the latter is not (Demkowicz 2007; Demkowicz et al. 2008). Whilst separately, both AMR-h and AMR-p strategies may perform better than their corresponding uniform counterparts, a combined AMR-hp strategy may address both the discontinuous and the smooth parts of the numerical solution. That is, regardless of the underlying regularity of the solution, an AMR-hp strategy could attain (quasi-)optimal exponential rates of convergence (Babuška and Guo 1992; Guo and Babuška 1986a, 1986b).

Such AMR strategies require the numerical solution obtained over the current mesh, as well as some higher-fidelity solution obtained over a more-refined reference mesh that provides a better approximation to the exact solution than the current (coarser) mesh. The difference between the two solutions provides a measure of the error between the two meshes that are used within the AMR algorithm. The use of two meshes of differing refinement enables computationally efficient FE-based error estimators to be determined in some appropriate norm. These error estimators can then used to discern those FEs that contribute the most to the overall error in the solution; and thus, these are flagged for further h-refinement or p-refinement. In the case of elliptic problems, for which a solution is sought in the $H^1$-space, it is natural to measure the spatial discretization error in the energy norm; or equivalently, the $H^1$-(semi)norm.

Since the local error estimators are evaluated a posteriori, it is possible to minimize the overall discretization error in an automated succession of computations performed on locally-refined meshes. This algorithm can be terminated once the solution has converged to some prescribed tolerance (Ainsworth and Oden 2000). Otherwise, the next optimized mesh is determined by maximizing the rate of reducing the local error contributions per increase in the DoFs by selecting some optimal local-refinement strategy. Although simple enough to implement in 1D, complexities arise when constructing non-conforming, or irregular, meshes in multiple dimensions. For example, the optimal refinement of a given FE may directly conflict with that of its neighbor. Therefore, to facilitate their usage, and to reduce complexity in the strategy itself, heuristics such as the minimum rule and the 1-irregularity rule were conceived (Wang and Ragusa 2009a; Demkowicz et al. 2000).

Wang and Ragusa developed such an AMR-h(p) technique for the 1D multi-group NDE, by taking their reference mesh as that resulting from one-level of uniform hp-refinement of the coarser current mesh (Wang and Ragusa 2005). The explicit computation over both was deemed to be computationally inefficient in multiple dimensions. Therefore, the AMR-hp strategy was redeveloped, this time employing interpolation-based error estimates for problems in 1D; and then, more generally, in multiple dimensions (Wang and Ragusa 2009a, 2007). This approach had been applied by Rachowicz and Demkowicz et al. to other elliptic systems (Demkowicz 2007; Rachowicz, Oden, and Demkowicz 1989; Demkowicz, Rachowicz, and Devloo 2002).

Indeed, the ultimate objective of the numerical method is to determine some QoI to some prescribed numerical accuracy. Generally, one is interested in some functional of the QoI over some prescribed region of the problem domain may be of interest; e.g. volumetric capture or fission rates. Rannacher et al. presented the original mathematical and computational framework for goal-based AMR strategies (Bangerth and Rannacher 2003). They utilized local contributions to the error, or the residual, weighted by a mathematically rigorous adjoint-based measure of their importance to a functional of a prescribed QoI (Bangerth and Rannacher 2003). Based upon this research, Lathouwers optimized the ${\mathrm{k}}_{\text{eff}}$ for nuclear reactor physics problems (Lathouwers 2011a). For spatial regions of the computational domain, large local errors may be ignored if they do not contribute significantly to the error in the prescribed QoI. Conversely, those spatial regions with small local errors may be flagged for refinement if they have a strong influence on the prescribed goal. The importance is computed as the solution to the appropriately formulated physical adjoint (dual) of the direct (primal) governing equations. Therefore, such refinement strategies are said to employ dual-weighted residual (DWR) error measures.

For multi-group nuclear reactor physics problems, the notion of optimality is extended to the use of energy-dependent meshes. Owing to the significant variation of nuclear-material properties with respect to space and, more importantly, energy the smoothness of the multi-group scalar neutron flux components can vary significantly between the energy groups. The fast and epithermal groups tend to exhibit smoother flux profiles than the thermal ones. Accordingly, any mesh adapted to the higher-energy groups would likely be too coarse for the lower-energy groups. Therefore, any effective AMR strategy should be able to resolve effectively each energy group, to some prescribed numerical accuracy.

Wang and Ragusa developed adaptive hp-refinement strategies for a FEM-based discretization of the multi-group NDE that included energy-dependent computational meshes in multiple dimensions source terms (Wang and Ragusa 2009a). For a FEM-based spatial discretization and spherical harmonic (P ${}_{\mathrm{N}}$) angular discretization of the first-order form of the neutron transport (FONT) equation, Gofin et al. achieved the transfer of information between energy-dependent meshes as a conservative Galerkin projection between volume meshes (Goffin et al. 2013). This approach is similar to that employed to calculate the interpolation-based error estimates. The information being exchanged between any two different meshes must be interpolated onto a super-mesh that is formed as the intersection of two implicated meshes (Farrell and Maddison 2011).

Constructing the energy-dependent meshes independently of one another yielded a significant increase in the computational times when compared to similar analyses performed on a single energy-independent mesh. Ragusa, and separately Wang and Ragusa, mitigated this issue by limiting their adaptive h-refinement algorithm to a sequence of adapted meshes belonging to a hierarchy of nested Cartesian meshes, forming a tree-structure (Wang, Bangerth, and Ragusa 2009; Ragusa 2008). The lowest level of the tree being the coarsest mesh, which is common to all energy groups. Any local refinement of the computational mesh, consisting of the bisection of a FE in each parametric direction, forms four new FEs in 2D. Although the numerical integration procedure was exact, the problem geometries were restricted to those that could be exactly meshed by straight-sided FEs at the coarsest level.

With regard to nuclear reactor physics simulations, the state-of-the-art is to use high-fidelity, multi-scale and multi-physics, heterogeneous models of nuclear reactors using software such as MOOSE, CASL VERA and ENRICO (Gaston et al. 2009; Mylonakis et al. 2014; Turner et al. 2016; Martineau et al. 2020). An important aspect of high-fidelity models is the accurate representation of curvilinear geometries such as nuclear fuel pins, control rods and burnable absorbers. However, the FEM can only approximate conic sections and quadric surfaces using straight-sided linear, or curvilinear higher-order, FEs (Hansen and Kang 1973). It is important to note that conic sections and quadric surfaces can only be represented exactly using rational polynomials; e.g. non-uniform rational B-splines (NURBS) (Welch et al. 2017b). In general, even high-fidelity numerical models exhibit some geometric error (Welch et al. 2017b). In which case, an ancillary mesh generator must be (re)visited anytime a change to the refinement of the computational mesh incurs a change in its representation of the underlying geometry. This process of mesh (re-)generation is computationally demanding (Cottrell, Hughes, and Bazilevs 2009). Moreover, one must ensure that both the fissile mass (volume) and the neutron leakage (surface area) are simultaneously preserved.

Isogeometric analysis (IGA) methods are able to perform the computational analysis on the exact, underlying geometry, which can be described by some computer-aided geometric design (CAGD) software. The NURBS basis functions are utilized to represent both the geometry and the solution field within the problem domain. Rhino3D, SolidWorks, Siemens NX, CATIA and OpenSCAD are examples of such CAGD software that employ NURBS basis functions. In one sense, IGA reverses the classical isoparametric concept. The numerical analysis is performed on a mesh that exactly represents the underlying CAGD geometry at the coarsest-level of refinement, as well as at any subsequent level of refinement of that mesh. The IGA meshes are hierarchical in structure and any AMR may be performed within the IGA algorithm at runtime. Therefore, unlike the FEM, IGA-based spatial discretizations do not require the use of ancillary mesh generators. Since most CAGD software provide the mathematical description of a given geometry using NURBS basis functions, they can be passed directly to the IGA-based software to serve in the coarsest-computational-mesh description. This unifies and streamlines the CAGD and computer-aided engineering (CAE) analysis pipeline (Hughes, Cottrell, and Bazilevs 2005). That being said, generating trivariate NURBS-based computational meshes for highly complex 3D geometrical domains is challenging; similar problems are encountered when generating hexahedral-dominant FE computational meshes (Akhras et al. 2017).

The IGA method was first applied to the mono-energetic (1G), or one-speed, NDE for (two-dimensional) 2D curvilinear geometries by Hall et al. using a continuous Bubnov-Galerkin (CBG-)IGA spatial discretization (Hall, Eaton, and Williams 2012). This was extended to the multi-group NDE by Welch et al. who noted that $C^0$-continuous NURBS-based IGA performed better than a mass-preserved polygonal standard FEM and a non-mass-preserved higher-order FEM in 2D (Welch et al. 2017b). Moreover, the IGA method, compared to the FEM, could achieve improved accuracy per DoF by utilizing the higher-order continuity within a NURBS patch within the problem domain (Welch et al. 2017b). IGA has also been applied to the second-order even-parity (EP) and the self-adjoint angular flux (SAAF) forms of the NTE in 2D; and again, when compared to the FEM, the IGA method yielded higher accuracy per DoF (Welch et al. 2017a; Latimer et al. 2020).

Welch et al. briefly investigated a heuristic-based adaptive h-refinement strategy for both the NDE and EP-P${}_{\mathrm{N}}$ equations (Welch 2018; Welch et al. 2017a). In this research, the IGA meshes were adapted to minimize the global L₂-error over energy-dependent meshes. Owens et al. investigated a heuristic-based adaptive h-refinement strategy using a non-conforming, or hanging-node, sweep-based discontinuous Galerkin (DG) spatial discretization of the discrete ordinates (S${}_{\mathrm{N}}$) first-order form of the neutron transport (FONT-S${}_{\mathrm{N}}$) equation (Owens et al. 2017). Subsequently, this research was extended to develop goal-oriented adaptive h-refinement (DWR-AMR-h) algorithms over energy-dependent meshes for the FONT-S${}_{\mathrm{N}}$ equation (Owens, Kópházi, and Eaton 2017). The conservative mesh-to-mesh interpolation between the energy-dependent meshes was based upon a geometric hierarchy of computational meshes.

Local-refinement strategies inevitably lead to non-conforming, or hanging-node, meshes, whereby the patches may no longer be mathematically connected in a parametric sense. Although this does not pose any issue for CAGD, it does present one for the spatial discretization of the governing equations over NURBS-based computational meshes. The NURBS patches must be coupled explicitly in order to produce a well-posed and solvable system of linear algebraic equations. There are a number of approaches to develop non-conforming spatial discretizations of elliptic PIDEs. One approach may be to employ interior penalty (IP) schemes for DG spatial discretizations, where a penalty parameter is used to enforce continuity in the solution field between adjacent NURBS patches (Wilson et al. 2018). However, such an approach would, at least, double the DoFs.

For CBG spatial discretizations of elliptic PDEs, continuity must be enforced in a strong sense (Hughes, Cottrell, and Bazilevs 2005). For an IGA spatial discretization, Cottrell et al. developed a constraint-based approach that couples adjacent NURBS patches of different levels of spatial refinement via an operator that encapsulates the process of (knot-insertion) refinement (Cottrell, Hughes, and Bazilevs 2009; Cottrell, Hughes, and Reali 2007). Such an approach necessitates an ability to determine, at a given non-conforming interface, the more-refined (slave-)patch and the less-refined (master-)patch in order to constrain exactly the former with respect to the latter. This is the approach that both Welch and Latimer et al. adopted in their separate studies of AMR-h for second-order forms of the NTE (Welch et al. 2017a; Latimer et al. 2020). Although such a strategy is restricted to $C^0$-continuity between the NURBS patches themselves, it does exploit the potential of higher-order continuity of the NURBS basis functions within a patch. However, the refinement is not strictly local since it must propagate within a given NURBS patch, incurring additional and unnecessary DoFs. Moreover, it was not demonstrated whether it was possible to incorporate local p-refinement capabilities. Generalizing the constraint-based approach proposed by Cottrell et al. Coox et al. weakly enforced higher-order continuity across adjacent non-conforming NURBS-patch interfaces in their method of Virtual Uncommon-Knot-Inserted Master-Slave (VUKIMS) (Coox et al. 2017).

This paper considers both NRG-AMR-h(p) and DWR-AMR-h(p) algorithms for a CBG-IGA spatial discretization of the multi-group NDE. The work presented is based upon the author’s PhD thesis (Wilson 2021).

The basic notions of NURBS-based IGA and energy-dependent meshes are introduced in Section 2. The AMR-h(p) strategies and accompanying error measures are devised in Section 3. The direct (primal) multi-group NDE, along with its physical adjoint (dual) equation, are presented in Section 4; and, as per a CBG-IGA spatial discretization, their weak-forms are derived in Section 5. Higher-order continuity within a NURBS patch is sacrificed in order to facilitate the wholly-local nature of the AMR-h(p) strategies; and $C^0$-continuity is strongly enforced across the solution space. Any hanging-nodes are addressed by establishing constraint equations across irregular patch interfaces, which are based upon master-slave relationships and the conservative interpolation between surface meshes. A similar Galerkin projection is employed in the conservative interpolation between volume meshes when evaluating the interpolation-based error measures and the group-to-group neutron source terms over energy-dependent meshes. Finally, numerical results for reactor physics benchmark test cases are presented in Section 6 to verify the computational performance and accuracy of the AMR algorithms. The method of manufactured solutions (MMS) is used to verify the numerical convergence of the AMR strategies.

2. NURBS-based Isogeometric analysis

The geometries of most engineering problems are produced from a basis of NURBS, or some other related spline-based description, used to generate complex surfaces bounding solid objects. The simplest description of the geometry is typically based upon a decomposition of the problem into its constituent material regions. Regardless, the physical domain is exactly represented by a tessellation of NURBS patches, the union of which forms the physical mesh, $\mathcal{T}(V).$ The geometrical mapping is a parametrisation of $\overrightarrow{r}:\widehat{V}\to V,$ which takes a point in the parametric domain, ${\overrightarrow{r}}^{\text{−}1}\in \widehat{V}\subset {\mathbb{R}}^{\widehat{d}},$ where the NURBS basis is defined, and maps it to a point in Euclidean space, $\overrightarrow{r}\in V\subset {\mathbb{R}}^d,$ for $d,\widehat{d}=1,2\text{\hspace{0.5em}or\hspace{0.5em}}3.$ A NURBS basis function defined over the real physical space, $R(\overrightarrow{r}),$ is the composition of that function, $R({\overrightarrow{r}}^{\text{−}1}),$ with the inverse of the geometrical mapping; i.e., $R(\overrightarrow{r})=\{R({\overrightarrow{r}}^{\text{−}1})°{\overrightarrow{r}}^{\text{−}1}:\forall {\overrightarrow{r}}^{\text{−}1}\in \widehat{V}\}.$ The parametric domain is local to a given patch. Consequently, a NURBS basis function has finite local support on a single given patch, in both $\widehat{V}$ and V, as per the bijective mapping; i.e.,: (1) $$\text{supp}(R)=\{\overrightarrow{r}\in V_e:R(\overrightarrow{r})\ne 0,\forall {\overrightarrow{r}}^{\text{−}1}\in {\widehat{V}}_e\},\hspace{1em}\forall V_e\in \mathcal{T}(V).$$(1)

Multi-dimensional, $d>1,$ NURBS basis functions are constructed in a tensor-product fashion from univariate B-splines, which can be expressed via the Cox-de Boor recursion formula (Cox 1972; de Boor 1972). A description of the parametric domain, denoted by some parameter $0\le \xi \le 1$ in 1D, begins with an open knot-vector, ${\Xi }_{\xi },$ which is a set of non-decreasing coordinates: ${\Xi }_{\xi }={\left\{{\xi }_i\right\}}_{i=1}^{n+p+1},$ where ${\xi }_i$ are the knots; p is the polynomial degree of the basis; and n is the total number of basis functions that span ξ. Knots are points of reduced continuity that partition the parametric domain into knot-spans; the $i^{\mathit{\text{th}}}$ knot-span is defined by the half-closed interval $[{\xi }_i,{\xi }_{i+1}).$ A knot-value is said to have a multiplicity of $\dot{m}$ knots if $\dot{m}$ number of knots take on the same value. The continuity of the basis over any given knot-value is $C^{p-\dot{m}}.$ Within knot-spans of non-zero length, the basis is infinitely differentiable, $C^{\infty }$-continuous; i.e., for ${\xi }_i\ne {\xi }_{i+1}.$

Increasing the multiplicity of a knot to $\dot{m}=p$ renders the NURBS basis interpolatory at that given knot-value. Effectively, a new patch boundary is formed and $C^0$-continuity is naturally enforced across the new boundary.

The NURBS basis functions inherit many of their properties from the B-splines; e.g. the convex hull property and the variation diminishing property (Piegl and Tiller 1997). The support of B-splines, ${\left\{N_i^p\right\}}_{i=1}^N,$ is compact and extends over $p+1$ knot-spans; i.e., for $i=1,\dots ,N:$ (2) $$\text{supp}(N_i^p)=\{\xi \in \widehat{V}\subset {\mathbb{R}}^1:N_i^p(\xi )\ne 0,\forall {\xi }_i\le \xi <{\xi }_{i+p+1}\}\text{;}$$(2) and $p+1$ functions always have support within any given knot-span.

Equation (3)(3) $$R_{i,j}^{p,q}(\xi ,\eta )=\frac{N_i^p(\xi )M_j^q(\eta )W_{i,j}}{\sum _{i\prime =1}^N\sum _{j\prime =1}^M{N}_{i\prime }^p(\xi )M_{j\prime }^q(\eta )W_{i\prime ,j\prime }}.$$(3) expresses the tensor-product fashion in which bi-variate NURBS basis functions, $R_{i,j}^{p,q}(\xi ,\eta ),$ are formed from sets of univariate B-splines, ${\{N_i^p(\xi )\}}_{i=1}^N$ and ${\{M_j^q(\eta )\}}_{j=1}^M.$ Respectively, the polynomial degrees of each set of univariate B-splines are denoted by p and q; and the knot-vectors are denoted by ${\Xi }_{\xi }={\left\{{\xi }_i\right\}}_{i=1}^{N+p+1}$ and ${\Xi }_{\eta }={\left\{{\eta }_j\right\}}_{j=1}^{M+q+1}.$

Each NURBS basis function corresponds to a control-point, a vector-valued coefficient. The set of control-points serves as the scaffolding for the problem geometry; that is the control-net, $\mathcal{B}={\left\{{\overrightarrow{B}}_{i,j}\right\}}_{i,j=1}^{N,M}\subset {\mathbb{R}}^2$ in 2D. A NURBS object, $V\in {\mathbb{R}}^d,$ is the piecewise projective transformation of a B-spline object in ${\mathbb{R}}^{d+1}.$ The latter is defined as a set of projective control-points, to each of which is attributed a weight; i.e., $\mathcal{W}={\left\{W_{i,j}\right\}}_{i,j=1}^{N,M}\subset \mathbb{R}$ in 2D. The control-points in the real physical space are obtained by dividing each of the projective control-points by its respective weight. These weights convey a sense of affinity for the control-points in the real physical space. This same transformation is applied to every point in the projective space. The NURBS basis functions are rendered rational by the division of the polynomial weighting function, the denominator of Equation (3)(3) $$R_{i,j}^{p,q}(\xi ,\eta )=\frac{N_i^p(\xi )M_j^q(\eta )W_{i,j}}{\sum _{i\prime =1}^N\sum _{j\prime =1}^M{N}_{i\prime }^p(\xi )M_{j\prime }^q(\eta )W_{i\prime ,j\prime }}.$$(3) ; i.e., for $i,j=1,\dots ,N,M:$ (3) $$R_{i,j}^{p,q}(\xi ,\eta )=\frac{N_i^p(\xi )M_j^q(\eta )W_{i,j}}{\sum _{i\prime =1}^N\sum _{j\prime =1}^M{N}_{i\prime }^p(\xi )M_{j\prime }^q(\eta )W_{i\prime ,j\prime }}.$$(3)

The NURBS basis functions are point-wise non-negative for non-negative weights.

The parametrisation, $\overrightarrow{r}:\widehat{V}\to V,$ can now be defined as the linear combination of NURBS basis functions and control-points that define a surface, S: (4) $$S(\xi ,\eta )=\sum _{i=1}^N\sum _{j=1}^M{\overrightarrow{B}}_{i,j}{R}_{i,j}^{p,q}(\xi ,\eta ),\hspace{1em}\forall \xi ,\eta \in \widehat{V}.$$(4)

The gradient over the domain with respect to the parametric space, $\widehat{\nabla },$ can be expressed in terms of partial derivatives of the B-splines; but the gradient over the domain with respect to the real physical space, $\overrightarrow{\nabla },$ must be expressed as: (5) $${\underset{¯}{\underset{¯}{\mathbf{J}}}}_V·\overrightarrow{\nabla }{R}_i^{p,q}(\xi ,\eta )=\widehat{\nabla }{R}_i^{p,q}(\xi ,\eta ),\hspace{1em}\text{for\hspace{0.5em}}i=1,\dots ,N\times M.$$(5)

The convention of the gradient operator is retained with the unit-vectors defined parallel to each of the axes in a rectangular Cartesian coordinate system; i.e., for ${\overrightarrow{e}}_{\xi }={\overrightarrow{e}}_x={[1,0,0]}^T$ and ${\overrightarrow{e}}_{\eta }={\overrightarrow{e}}_y={[0,1,0]}^T:$ (6) $$\widehat{\nabla }\equiv {\overrightarrow{e}}_{\xi }\frac{\partial }{\partial \xi }+{\overrightarrow{e}}_{\eta }\frac{\partial }{\partial \eta }.$$(6)

The volumetric Jacobian $(2\times 2)$-matrix, is defined below: (7) $${\underset{¯}{\underset{¯}{\mathbf{J}}}}_V:=\left(\begin{array}{cc}\partial x/\partial \xi & \partial y/\partial \xi \\ \partial x/\partial \eta & \partial y/\partial \eta \end{array}\right).$$(7)

Stability of the discretization necessitates the mapping to be invertible and nonsingular, $\mathrm{\text{det}}{\underset{¯}{\underset{¯}{\mathbf{J}}}}_V>0,$ almost everywhere, such that the NURBS basis functions and their derivatives are square-integrable.

For the example of a pin-cell constructed from five bi-quadratic NURBS patches, the schematic in Figure 1 presents the domain of a single NURBS patch as that represented in the three classical spaces of IGA; viz. the real physical, the parametric and the parent spaces. The exact representation of the underlying geometry is maintained as the mesh is refined because the geometrical mapping and the Jacobian are invariant under refinement (Piegl and Tiller 1997).

Figure 1. An example of a 2D pin-cell represented by bi-quadratic NURBS patches (macro-elements) in the real physical space, $V,$ as denoted by the parameters x and y. The domain, represented in the parametric space, $\widehat{V},$ is formed from the tensor-product of two 1D spaces, as denoted by the parameters $\widehat{\xi }$ and $\widehat{\eta }.$ The parametric space has been decomposed into non-vanishing knot-spans (micro-elements), over which, numerical integration is performed by an affine mapping from the parent domain, $\tilde{V},$ as denoted by the parameters $\tilde{\xi }$ and $\tilde{\eta }.$ The patch boundaries are highlighted in red; whereas the non-vanishing knot-spans within a patch are separated by black lines. The image of given non-vanishing knot-span, shaded in dark grey, can be represented in all three spaces. (V. the web-based version for reference to color.)

One way in which to enrich the parametric domain, as suggested by Figure 1, is to form additional non-zero-knot-spans. Knots can be inserted with a multiplicity of one, KI-1-refinement, up to $p+1,$ KI-h-refinement. Both of these types of KI-refinement are demonstrated in Figure 2c and d, respectively, for the example of a quadratic NURBS curve. One clarifies that the latter is analogous to h-refinement in a standard FEM because the basis is reduced to $C^0$-continuity across the newly-formed element boundaries. Another way in which to enrich the basis is by degree-elevation (DE), which is analagous to p-refinement in a standard FEM. This is demonstrated in Figure 2b. One notes that the refinement processes of KI and DE are not commutative (Cottrell, Hughes, and Reali 2007).

Figure 2. The different mechanisms for refinement for the example of a quadratic NURBS curve, whose coarsest description is given in Figure 2a. DE-refinement is presented in Figure 2b; and the insertion of a single knot and p-number of knots at $\xi =0.5$ is presented in Figure 2c and Figure 2d, respectively. The blue or black segments represent the image of a knot-span in (x, y); and the red circles represent the control-points. (V. the web-based version for reference to color.)

2.1. Energy-dependent meshes

The microscopic neutron cross-section data can vary significantly across the range of possible neutron energies. Therefore, the spatial variation in the scalar neutron flux can also vary significantly across the different neutron energy groups. Faster neutrons tend to have longer mean free paths (MFPs) and typically exhibit smoother spatial distributions of the scalar neutron flux. Conversely, slower neutron energy groups tend to have shorter MFPs and typically exhibit more significant variations in the spatial distribution of the scalar neutron flux. This is especially the case near material interfaces and boundaries of the domain. Optimality could not be achieved over a single energy-independent mesh, $\mathcal{T}(V),$ and the associated finite dimensional function space, ${\widehat{\mathcal{V}}}_{h,p},$ as characterized by the mesh spacing, h, and the polynomial degree of the basis, p. One should expect to be able to maximize the accuracy per total DoF by solving the multi-group NTE over separate energy-dependent meshes, ${\mathcal{T}}^g(V),$ and their respective function spaces, ${\widehat{\mathcal{V}}}_{h,p}^g,$ which may be refined independently; i.e., $G\times \mathrm{\text{dim}}{\widehat{\mathcal{V}}}_{h,p}>\sum {}_g\mathrm{\text{dim}}{\widehat{\mathcal{V}}}_{h,p}^g.$

One designates the interpolation $\prod {}_{h,p}^t{v}^d:{\mathcal{V}}^d(V)\to {\mathcal{V}}^t(V)$ as the projection of some function, $v^d\in {\mathcal{V}}_{h,p}^d,$ defined on a donor mesh, ${\mathcal{T}}^d(V),$ onto a target mesh, ${\mathcal{T}}^t(V).$ The optimum target interpolant, $v^t\in {\mathcal{V}}_{h,p}^t,$ is approximated such that: $$\Vert v^t-v^d\Vert =\underset{v\in {\mathcal{V}}_{h,p}^t}{\mathrm{\text{min}}}\Vert v-v^d\Vert .$$

In order to preserve locality and continuity, the residual, as measured in some norm, is made to be negligibly small by solving the appropriate Sobolev minimization problem, from which conservation follows naturally (Jiao and Heath 2004). The functions are those taken from the same spaces spanned by the NURBS basis functions used to define each of the donor-mesh and the target-mesh, denoted by the superscripts d and t, respectively. By letting $N^i=\mathrm{\text{dim}}{\widehat{\mathcal{V}}}_{h,p}^i,$ the NURBS basis functions defined over a given mesh are also denoted by a similar superscript; i.e., for $i=d,t$ and for $j=1,\dots ,N^i:$ $$R_j^i(\overrightarrow{r})={(R(\overrightarrow{r}){|}_{{\mathcal{T}}^i})}_j,\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\overrightarrow{r}\in V.$$

The weighted-residual formulation is chosen to be optimal in the $H^1$-norm and the residual is minimized if the first variation in the control-variables of the target function, $v_i^t,$ is negligibly small for $i=1,\dots ,N^t.$ This is analogous to an orthogonal Galerkin projection of a discrete field onto the function space spanned by the target basis, where the projection is the function that minimizes the interpolation error in the desired norm (Farrell and Maddison 2011). This can be expressed in matrix form as: (8a) $$[{\underset{¯}{\underset{¯}{S}}}_V^{\mathit{\text{tt}}}+{\underset{¯}{\underset{¯}{M}}}_V^{\mathit{\text{tt}}}]·{\underset{¯}{v}}^t=[{\underset{¯}{\underset{¯}{S}}}_V^{\mathit{\text{td}}}+{\underset{¯}{\underset{¯}{M}}}_V^{\mathit{\text{td}}}]·{\underset{¯}{v}}^d\text{;}$$(8a) (8b) $$\Rightarrow {\underset{¯}{\underset{¯}{K}}}_V^{\mathit{\text{tt}}}·{\underset{¯}{v}}^t={\underset{¯}{\underset{¯}{K}}}_V^{\mathit{\text{td}}}·{\underset{¯}{v}}^d\text{;}$$(8b) where the consistent mass and stiffness matrices, ${\underset{¯}{\underset{¯}{M}}}_V^{\mathit{\text{tt}}}$ and ${\underset{¯}{\underset{¯}{S}}}_V^{\mathit{\text{tt}}},$ respectively, are formed by integrating the product of those functions, defined over ${\mathcal{T}}^t,$ that belong to the target function space, ${\widehat{\mathcal{V}}}_{h,p}^t.$ Whereas the projective mass and stiffness matrices, ${\underset{¯}{\underset{¯}{M}}}_V^{\mathit{\text{td}}}$ and ${\underset{¯}{\underset{¯}{S}}}_V^{\mathit{\text{td}}},$ respectively, require that in the same integrand are evaluated the product of functions defined over different spatial meshes, ${\mathcal{T}}^t$ and ${\mathcal{T}}^d(V),$ that belong to different function spaces, ${\widehat{\mathcal{V}}}_{h,p}^t$ and ${\widehat{\mathcal{V}}}_{h,p}^d,$ respectively. This is achieved by integrating over an implied super-mesh formed by the intersection of the target-mesh and the donor-mesh, ${\mathcal{T}}^g\cap {\mathcal{T}}^{g\prime },$ which decomposes the parametric domain into subdomains appropriate for numerical integration; v. Figure 3. Although it has not yet been detailed, an ability to perform local refinement has been assumed.

Figure 3. An example of the implied parametric super-mesh, ${\mathcal{T}}^{1\&2},$ formed as the intersection of two volume meshes, ${\mathcal{T}}^1$ and ${\mathcal{T}}^2,$ for the conservative interpolation of a scalar field from one to the other. The crosses represent the appropriate number of quadrature points required to properly integrate that element of the super-mesh. (V. the web-based version for reference to color.)

Looping over the elements formed by the intersection of two meshes, the bases in each parametric direction are sampled at a number of quadrature points, N_q, which is determined by: (9) $$N_q=\underset{i=t,d}{\mathrm{\text{max}}}{N}_q{|}_{{\mathcal{T}}^i},\hspace{1em}\forall V_{\widehat{e}}\in \underset{i}{\cap }{\mathcal{T}}^i(V).$$(9)

It was found that employing $p+1$ quadrature points in each parametric direction, over the canonical representation of a given element in the parent domain, provides sufficient accuracy.

Complexity is largely circumvented if the energy-dependent meshes are derived from a common coarsest-mesh description because any subsequent refinement results in a hierarchical set of nested computational meshes. The intersection of an element from one mesh with another from a different mesh always forms either an empty (disjoint) set, a shared interface, $V_{i\prime }^g\subseteq V_{j\prime }^{g\prime }$ or $V_{i\prime }^g\supseteq V_{j\prime }^g.$ The tensor-product structure of the parametric space of NURBS-based IGA lends itself to such an approach. Since the coarsest-mesh description is exactly the original CAGD description, each subsequent mesh nested within a hierarchy of meshes would also exactly represent the underlying geometry.

Some function may undergo a change of basis transformation by pre-multiplying the control-variables, ${\underset{¯}{v}}^d,$ by the combined transfer matrix, ${\underset{¯}{\underset{¯}{K}}}_V^{t\leftarrow d},$ for which the interpolation of the function from a donor mesh onto a target mesh, is optimal with respect to the $H^1$-norm. Although ${\underset{¯}{\underset{¯}{S}}}_V^{\mathit{\text{tt}}}$ is singular, ${\underset{¯}{\underset{¯}{K}}}_V^{\mathit{\text{tt}}}$ is symmetric-positive definite; and therefore, it is readily invertible using standard iterative smoothers: (10a) $${\underset{¯}{v}}^t={\left[{\underset{¯}{\underset{¯}{K}}}_V^{\mathit{\text{tt}}}\right]}^{-1}·{\underset{¯}{\underset{¯}{K}}}_V^{\mathit{\text{td}}}·{\underset{¯}{v}}^d,$$(10a) (10b) $$=\underset{(N^t\times N^d)(N^d\times 1)}{{\underset{¯}{\underset{¯}{K}}}_V^{t\leftarrow d}·{\underset{¯}{v}}^d.}$$(10b)

For continuous Galerkin discretizations, it becomes necessary to (re)use the sparse matrix storage structures and efficient smoothers to perform the interpolation (Balay et al. 2018).

Group-to-group scattering and fission contributions may not exist for some ordered pairs of energy groups, $(g\prime ,g),$ such that it would be unnecessary to perform the interpolation $\prod {}_{h,p}^g{\widehat{\varphi }}_{h,p}^{g\prime }:{\widehat{\mathcal{V}}}_{h,p}^{g\prime }\to {\widehat{\mathcal{V}}}_{h,p}^g;$ this is particularly pertinent to those problems that exhibit mostly down-scatter for $G\gg 1.$ For those ordered pairs in which information is transferred in both directions, additional savings can be made by ${\underset{¯}{\underset{¯}{K}}}_V^{g\leftarrow g\prime }={[{\underset{¯}{\underset{¯}{K}}}_V^{g\prime \leftarrow g}]}^T;$ in which case, only ${\underset{¯}{\underset{¯}{K}}}_V^{g\prime \to g}$ must be stored for $g\prime \ge g$ (Owens, Kópházi, and Eaton 2017).

The interpolation of a discrete scalar field between meshes is not unique in its application to energy-dependent meshes. It is also employed in the next section when devising interpolation-based error measures to drive the adaptive refinement strategies; and again as a means of enforcing $C^0$-continuity across irregular meshes.

3. The adaptive mesh refinement (AMR) strategies

IGA provides a good framework within which to implement AMR-h(p) strategies because the parametrisation does not vary under refinement; and thus, the underlying geometry is always exactly represented; cf Figure 2. The research presented in this paper employs only standard tensor-product NURBS patches in the proposed AMR algorithms. Therefore, these algorithms can be used with existing CAGD software.

Increased smoothness in the numerical solution is something to be desired when solving elliptic PDE problems; especially as higher-order continuity, $C^{p-1},$ between elements is achievable with NURBS-based IGA. However, this property is not investigated here. At most, $C^0$-continuity is considered to facilitate a wholly-local hp-refinement strategy. In fact, reduced continuity in the solution may help to reduce any spurious oscillations due to the presence of strong material heterogeneity (Reed 1971).

In the limit of a $C^0$-continuous basis, one can exploit the efficient and well-established pre-conditioners that have been developed for low-order continuity FEMs, which are not generally effective for bases that exhibit higher-order continuity (Gahalaut, Kraus, and Tomar 2013; Sangalli and Tani 2016).

Moreover, the tensor-product structure of the NURBS patches lends itself to the cell-based approach adopted by Wang and Ragusa (Wang, Bangerth, and Ragusa 2009). In this approach, any local isotropic h-refinement consists of bisecting an element in each parametric direction. Therefore, the AMR procedure can proceed as a sequence of adapted meshes that belong to a hierarchy of nested meshes.

The proposed AMR-h(p) strategy comprises an iterative mesh refinement procedure that requires all of the group-to-group sources to be fully converged, to some prescribed tolerance, before refining the spatial mesh(es) for the next AMR iteration; v. Figure A1 in the Appendix. Within a given AMR iteration, the multi-group components of the numerical solution must be known in order to evaluate the global error in the QoI, ${ϵ}^{(i)},$ as measured in an appropriate norm: (11) $${ϵ}^{(i)}=\Vert \mathit{\text{QoI}}(\varphi )-\mathit{\text{QoI}}({\widehat{\varphi }}_{h,p}^{(i)})\Vert .$$(11)

The QoI represents a physical quantity, such as the scalar neutron flux, at a point in the domain, or averaged over a region. Where ${\widehat{\varphi }}_{h,p}^{(i)}$ is the discrete solution obtained from the numerical analysis performed on the mesh(es), ${\mathcal{T}}^{(i)}(V),$ characterized by some combination of h and p, at the $i^{\mathit{\text{th}}}$-iteration of the AMR procedure, as denoted by the parenthesized index $(i).$ The relative global error, ${\overline{ϵ}}^{(i)},$ is defined as: (12) $${\overline{ϵ}}^{(i)}=\frac{\Vert \mathit{\text{QoI}}(\varphi )-\mathit{\text{QoI}}({\widehat{\varphi }}_{h,p}^{(i)})\Vert }{\Vert \mathit{\text{QoI}}(\varphi )\Vert }=\frac{{ϵ}^{(i)}}{\Vert \mathit{\text{QoI}}(\varphi )\Vert }\text{;}$$(12) where the AMR-iteration is terminated once ${\overline{ϵ}}^{(i)}<\tau ,$ for some prescribed tolerance, τ; or once the maximum number of iterations has been attained.

The proposed AMR strategy aims to maximize the decrease in the global error, $\Delta {ϵ}^{(i)},$ between successive iterations, per increase in the number of unknowns, $\Delta \mathit{\text{Do}}{F}^{(i)},$ by the selection of some combination of $h^{*},p^{*}$ (Ceze and Fidkowski 2013): (13a) $$\mathrm{\text{max}}\frac{\Delta {ϵ}^{(i)}}{\Delta \mathit{\text{Do}}{F}^{(i)}}=\mathrm{\text{max}}\frac{{ϵ}^{(i)}-{ϵ}^{(i+1)}}{|\mathit{\text{Do}}{F}^{(i)}-\mathit{\text{Do}}{F}^{(i+1)}|},$$(13a) (13b) $$=\underset{h^{*},p^{*}}{\mathrm{\text{max}}}\frac{\Vert \mathit{\text{QoI}}(\varphi )-\mathit{\text{QoI}}({\widehat{\varphi }}_{h,p}^{(i)})\Vert -\Vert \mathit{\text{QoI}}(\varphi )-\mathit{\text{QoI}}({\widehat{\varphi }}_{h^{*},p^{*}}^{(i+1)})\Vert }{|\mathit{\text{Do}}{F}_{h,p}^{(i)}-\mathit{\text{Do}}{F}_{h^{*},p^{*}}^{(i+1)}|},$$(13b) (13c) $$=\underset{h^{*},p^{*}}{\mathrm{\text{max}}}\overline{\Delta }{ϵ}_{h^{*},p^{*}}^{(i)}.$$(13c)

The mechanisms for spatial refinement are restricted to isotropic DE-1-refinement and isotropic KI-h-refinement; i.e., $h^{*},p^{*}=\{h/2,p\text{or}h,p+1\}.$ This is comparable to a standard AMR-h(p) strategy for a classical FEM. It is important to note from Equation (13) that the cost of performing KI-h-refinement becomes computationally more demanding as the polynomial degree of the basis, p, increases. This is due to the fact that the number of unknowns is increased for every knot inserted; and knots must be inserted with a multiplicity of $\dot{m}=p$ in order to render the basis $C^0$-continuous.

Two different AMR strategies are introduced in the following subsections. The first aims to minimize, almost everywhere, the discretization error in the solution, as measured in an appropriate norm. This is called the energy-driven (NRG-)AMR strategy, although the error does not have to be measured in the energy-norm. The second strategy aims to minimize the error in some specific goal. This is called the goal-based (DWR-)AMR strategy. Neither strategy considers coarsening of the computational mesh.

Whilst the multi-group NDE has yet to be formally presented and Discretized, it is sufficient to assume that the weak problem can be posed as finding $\widehat{\varphi }\in \widehat{\mathcal{V}},$ such that: (14) $$a(\widehat{\varphi },\widehat{w})=b(\widehat{w}),\hspace{1em}\forall \widehat{w}\in \widehat{\mathcal{V}}.$$(14)

One assumes that the bilinear form, $a(·,·)=\sum {}_g{a}^g(·,·),$ has the properties of continuity and coercivity.

3.1. Energy-driven (NRG) error measures

The NRG-AMR strategy aims to minimize the residual, or rather the discretization error, $e=\varphi -{\widehat{\varphi }}_{h,p},$ over the entire computational domain by refining those regions that contribute significantly to e, as measured in the energy-norm; viz. from Equation (11)(11) $${ϵ}^{(i)}=\Vert \mathit{\text{QoI}}(\varphi )-\mathit{\text{QoI}}({\widehat{\varphi }}_{h,p}^{(i)})\Vert .$$(11) , $ϵ=\Vert e\Vert .$ By the equivalence of norms on the chosen function space and the assumption that the first-order derivative term is dominant, the discretization error can be measured in the $H^1$-seminorm, $|·{|}_{H^1};$ i.e.,: (15) $$|\varphi -{\widehat{\varphi }}_{h,p}{|}_{H^1}^2=|e{|}_{H^1}^2={\Vert \overrightarrow{\nabla }e\Vert }_{L_2}^2={\int }_V\sum _i^d[\frac{\partial }{\partial x_i}e(\overrightarrow{r}){]}^2\mathit{\text{ dV}}.$$(15)

A numerical reference solution, ${\widehat{\varphi }}_{\mathit{\text{ref}}},$ is proposed in lieu of the true solution, $\varphi .$ Although still an approximation, the numerical reference solution, calculated over some more-refined reference mesh(es), ${\mathcal{T}}_{\mathit{\text{ref}}},$ is a more-accurate approximation of the true solution. With respect to Equation (13), the reference solution must sufficiently resolve the local scales of either of the two proposed local-refinement options, KI-h-refinement and DE-1-refinement, in order to justify the use of one over the other in the next iteration of mesh refinement; i.e., ${\widehat{\mathcal{V}}}_{h^{*},p^{*}}(V_e)\subset {\widehat{\mathcal{V}}}_{\mathit{\text{ref}}}(V_e),$ (Demkowicz 2007). Accordingly, the reference mesh should be at least one-level of uniform KI-h-refinement and DE-1-refinement above that of the coarser current mesh; v. Figure 4, where ${\mathcal{T}}_{\mathit{\text{ref}}}={\mathcal{T}}_{\text{h/2},p\text{+1}}$ and therefore, ${\widehat{\varphi }}_{\mathit{\text{ref}}}={\widehat{\varphi }}_{\text{h/2},p\text{+1}}\in {\widehat{\mathcal{V}}}_{\text{h/2},p\text{+1}}.$ All measurements of error are now estimates made with respect to the reference solution; i.e., $\widehat{ϵ}=|\widehat{e}{|}_{H^1},$ where $\widehat{e}={\widehat{\varphi }}_{\mathit{\text{ref}}}-{\widehat{\varphi }}_{h,p}.$

Figure 4. A schematic of the interpolation between a more-refined reference mesh and a coarser current mesh; and the projections to prospective subspaces for the next iteration of the AMR procedure. (V. the web-based version for reference to color.)

For ${\widehat{\varphi }}_{\mathit{\text{ref}}}\in {\widehat{\mathcal{V}}}_{\mathit{\text{ref}}}$ and for ${\widehat{\varphi }}_{h,p}\in {\widehat{\mathcal{V}}}_{h,p}\subset {\widehat{\mathcal{V}}}_{\mathit{\text{ref}}}$ it can be shown that the global error estimate, $\widehat{ϵ},$ provides an upper-bound to the true global error, under the assumption that the reference solution is a much better approximation of the true solution than is ${\widehat{\varphi }}_{h,p}$ calculated on ${\mathcal{T}}_{h,p};$ i.e., $\Vert \varphi -{\widehat{\varphi }}_{\mathit{\text{ref}}}\Vert \ll \Vert {\widehat{\varphi }}_{\mathit{\text{ref}}}-{\widehat{\varphi }}_{h,p}\Vert$ (Wilson 2021).

The performance of a posteriori error estimation is measured by the (reciprocal) effectivity index (Bangerth and Rannacher 2003): (16) $$\theta =\frac{\widehat{ϵ}}{ϵ}.$$(16)

Again, if the true solution is not available, this too must be approximated by some numerical reference solution, must be more accurate than ${\widehat{\varphi }}_{\text{h/2},p\text{+1}}^g.$ It is desirable that $\theta \to 1$ as $h\to 0$ and $p\to \infty .$

Owens et al. investigated the effectiveness of using error estimates based upon the explicit calculation of a current coarse solution, as well as a more-refined reference solution, to drive their AMR-h strategy for the DG-IGA discretization of the FONT-S ${}_{\mathrm{N}}$ equation in 2D. The reference mesh was taken to be that resulting from a uniform DE-1-refinement of the current coarser mesh (Owens, Kópházi, and Eaton 2017). It was observed that the accuracy in the estimation of the error degraded as the polynomial degree of the NURBS basis increased. As per the tensor-product structure of the NURBS basis, the ratio of DoFs per element in the reference mesh to that of the current mesh, ${(p+2)}^d/{(p+1)}^d,$ is a monotonically decreasing function for positive integer values of p. Therefore, quadratic solutions are significantly more accurate than linear ones; yet, cubic solutions are not much more accurate than quadratic ones.

It would be computationally expensive to perform two explicit calculations over two different levels of mesh refinement for each energy group. Fortunately, as per Céa’s Lemma for conforming and consistent discretizations, one can show that the chosen approximation space will yield a solution, ${\widehat{\varphi }}_{h,p}\in {\widehat{\mathcal{V}}}_{h,p},$ whose error is bounded by the quasi-best approximation error. The quasi-best approximation error is bounded by the difference in the exact solution and any function that belongs to the finite-dimensional hp function space. In particular, one can choose the projection of the exact solution onto this approximation space, $\prod {}_{h,p}\varphi .$ In doing so, the entire AMR procedure can be driven by controlling the interpolation error that bounds the quasi-best approximation error, which, in turn, bounds the true discretization error; i.e., (Demkowicz 2004): $${\Vert \varphi -{\widehat{\varphi }}_{h,p}\Vert }_{H^1}\le \frac{C_{\alpha }}{C_{\beta }}\underset{{\widehat{\varphi }}_{h,p}\in {\widehat{\mathcal{V}}}_{h,p}}{\mathrm{\text{inf}}}{\Vert \varphi -{\widehat{\varphi }}_{h,p}\Vert }_{H^1}\le \frac{C_{\alpha }}{C_{\beta }}{\Vert \varphi -\prod _{h,p}\varphi \Vert }_{H^1}\text{;}$$ where the constants of continuity, $C_{\alpha },$ and coercivity, $C_{\beta },$ depend upon the nature of the discretization.

Substitution of the exact solution for the reference solution yields an iterative procedure that requires the numerical solution over a more-refined reference mesh and its projection onto a coarser current mesh. This projection, $\prod {}_{h,p}:{\widehat{\mathcal{V}}}_{\text{h/2},p\text{+1}}\to {\widehat{\mathcal{V}}}_{h,p},$ is performed as a conservative interpolation between volume meshes. Wang and Ragusa did not observe any loss in accuracy when employing such an interpolation-based approach; and especially so, in especially in the asymptotic range of convergence (Wang and Ragusa 2009a).

An estimate of the global discretization error, based upon the interpolation of a reference solution, is summarized below: (17a) $${ϵ}^2=\sum _{g=1}^G{({ϵ}^g)}^2,$$(17a) (17b) $$=\sum _{g=1}^G{\Vert {\varphi }^g-{\widehat{\varphi }}_{h,p}^g\Vert }_{H^1}^2,$$(17b) (17c) $$\le \sum _{g=1}^G{\Vert {\varphi }^g-\prod _{h,p}{\varphi }^g\Vert }_{H^1}^2,$$(17c) (17d) $$\lesssim \sum _{g=1}^G{\Vert {\widehat{\varphi }}_{\text{h/2},p\text{+1}}^g-\prod _{h,p}{\widehat{\varphi }}_{\text{h/2},p\text{+1}}^g\Vert }_{H^1}^2,$$(17d) (17e) $$\lesssim \sum _{g=1}^G\sum _{V_e\in {\mathcal{T}}_{h,p}^g}{\Vert {\widehat{\varphi }}_{\text{h/2},p\text{+1}}^g-\prod _{h,p}{\widehat{\varphi }}_{\text{h/2},p\text{+1}}^g\Vert }_{H^1(V_e)}^2,$$(17e) (17f) $$\lesssim \sum _{g=1}^G\sum _{V_e\in {\mathcal{T}}_{h,p}^g}{\mu }_e^g.$$(17f)

The global discretization error estimate, $\widehat{ϵ},$ is now a projection-based error estimate of the global discretization error, bounded by the sum of element-local interpolation-based error measures, ${\mu }_e^g;$ i.e.,: (18) $$\widehat{ϵ}={\left[\sum _{g=1}^G{({\widehat{ϵ}}^g)}^2\right]}^{\frac{1}{2}}={\left[\sum _{g=1}^G\sum _{V_e\in {\mathcal{T}}_{h,p}^g}{\mu }_e^g\right]}^{\frac{1}{2}}\text{;}$$(18) where by the Poincaré inequality, ${\mu }_e^g$ are defined in terms of the first-order Sobolev space seminorm, for $g=1,\dots ,G:$ (19) $${\mu }_e^g=|{\widehat{\varphi }}_{\text{h/2},p\text{+1}}^g-\prod _{h,p}{\widehat{\varphi }}_{\text{h/2},p\text{+1}}^g{|}_{H^1(V_e)}^2,\hspace{1em}\forall V_e\in {\mathcal{T}}_{h,p}^g(V).$$(19)

In the event that only a single energy-independent mesh is employed, the error measure is also made to be energy-independent; i.e., ${\mu }_e=\sum {}_g{\mu }_e^g.$

Driven by the residual in the interpolation-based error estimate, successive meshes at each iteration of the AMR procedure are obtained by minimizing local contributions to the global error. An element is flagged if its local contribution to the global error is greater than some threshold, which is taken as some fraction, ${\beta }_1\in [0,1],$ of the maximum error over a single element belonging to any of the current coarser meshes. For energy-dependent meshes, $g=1,\dots ,G,$ at a some given iteration count in the AMR procedure, (i), the first decision in the AMR procedure, D-1 in Figure A1, is thus, to determine those elements that satisfy: (20) $$w^{g(i)}{\mu }_e^{g(i)}\ge {\beta }_1\underset{g\prime =1,\dots ,G^{}}{\mathrm{\text{max}}}\underset{V_{e\prime }\in {\mathcal{T}}_{h,p}^{g\prime (i)}}{\mathrm{\text{max}}}{w}^{g\prime (i)}{\mu }_{e\prime }^{g\prime (i)},\hspace{1em}\forall V_e\in {\mathcal{T}}_{h,p}^{g(i)}(V).$$(20)

Although somewhat arbitrary, experience from the literature suggests that ${\beta }_1=\text{1/3}$ provides a good compromise of the number of elements refined at each AMR-iteration and the total number of AMR-iterations (Rachowicz, Oden, and Demkowicz 1989).

Non-smooth multi-group components of the solution are targeted by refining all of the energy-dependent meshes at once. One introduces group-importance factors, $w^{g(i)},$ in Equation (20)(20) $$w^{g(i)}{\mu }_e^{g(i)}\ge {\beta }_1\underset{g\prime =1,\dots ,G^{}}{\mathrm{\text{max}}}\underset{V_{e\prime }\in {\mathcal{T}}_{h,p}^{g\prime (i)}}{\mathrm{\text{max}}}{w}^{g\prime (i)}{\mu }_{e\prime }^{g\prime (i)},\hspace{1em}\forall V_e\in {\mathcal{T}}_{h,p}^{g(i)}(V).$$(20) to normalize the group error measures with respect to the respective multi-group component of the scalar neutron flux; i.e., for $g=1,\dots ,G:$ (21) $$w^g=\frac{1}{\Vert {\widehat{\varphi }}_{\text{h/2},p\text{+1}}^g\Vert }.$$(21)

For those elements that satisfy Equation (20)(20) $$w^{g(i)}{\mu }_e^{g(i)}\ge {\beta }_1\underset{g\prime =1,\dots ,G^{}}{\mathrm{\text{max}}}\underset{V_{e\prime }\in {\mathcal{T}}_{h,p}^{g\prime (i)}}{\mathrm{\text{max}}}{w}^{g\prime (i)}{\mu }_{e\prime }^{g\prime (i)},\hspace{1em}\forall V_e\in {\mathcal{T}}_{h,p}^{g(i)}(V).$$(20) , an optimal local-refinement strategy is determined by maximizing the rate at which the local contribution to the interpolation-based error estimate can be decreased, Equation (13)(10a) $${\underset{¯}{v}}^t={\left[{\underset{¯}{\underset{¯}{K}}}_V^{\mathit{\text{tt}}}\right]}^{-1}·{\underset{¯}{\underset{¯}{K}}}_V^{\mathit{\text{td}}}·{\underset{¯}{v}}^d,$$(10a) .

The options for local refinement are restricted to $1\times$ KI-h-refinement or $1\times$ DE-1-refinement. Since the reference mesh is chosen to be ${\mathcal{T}}_{\text{h/2},p\text{+1}}(V),$ it is possible to interpolate the reference solution onto the prospective subspaces, ${\widehat{\mathcal{V}}}_{h^{*},p^{*}}(V_e),$ and evaluate, $\overline{\Delta }{ϵ}_{h^{*},p^{*}}^{(i)}{|}_{V_e};$ v. Figure 4. For $h^{*},p^{*}=\{h/2,p\text{or}h,p+1\},$ the optimal combination of $h^{*},p^{*}$ on $V_e\in {\mathcal{T}}^g(V)$ is chosen to be that which maximizes the following: (22) $$\underset{h^{*},p^{*}}{\mathrm{\text{max}}}\overline{\Delta }{ϵ}_{h^{*},p^{*}}^{g(i)}{|}_{V_e}=\underset{h^{*},p^{*}}{\mathrm{\text{max}}}\frac{{\mu }_e^{g(i)}-|{\widehat{\varphi }}_{\text{h/2},p\text{+1}}^{g(i)}-\prod _{h^{*},p^{*}}{\widehat{\varphi }}_{\text{h/2},p\text{+1}}^{g(i)}{|}_{H^1(V_e)}^2}{|\mathit{\text{Do}}{F}_{h,p}^{g(i)}-\mathit{\text{Do}}{F}_{h^{*},p^{*}}^{g(i)}{|}_{V_e}}.$$(22)

Certain constraints are imposed upon the use of either DE-1-refinement or KI-h-refinement. The first corresponds to D-2 in Figure A1, where an upper-limit is placed upon the polynomial degree of the NURBS basis; i.e.,: (23) $$p{|}_{V_e}\le {\beta }_2.$$(23)

The second corresponds to D-3 in Figure A1. Since KI-h-refinement reduces continuity in the basis, one restricts its use to those elements in which the estimated local error, ${\mu }_e^{g(i)},$ can be reduced by a certain amount. For some parameter, ${\beta }_3\in [0,1],$ KI-h-refinement is only considered if: (24) $$\Delta {ϵ}_{\text{h/2},\mathrm{p}}^{g(i)}{|}_{V_e}\ge {\beta }_3{\mu }_e^{g(i)}\text{;}$$(24) where ${\beta }_3=0$ implies that no restriction is placed upon the use of KI-h-refinement; whereas ${\beta }_3=1$ implies that KI-h-refinement is strictly prohibited. Without being too restrictive, one employs ${\beta }_3=\text{1/2}.$

To recapitulate, the parameters employed to constrain the mechanisms of local refinement are: (25) $$\{\begin{array}{cc}{\beta }_1=\text{1/3}\hfill & \mathrm{(}\text{Equation\hspace{0.5em}20}\mathrm{)}\text{;}\hfill \\ {\beta }_2=p_{\mathrm{\text{max}}}\hfill & \mathrm{(}\text{Equation\hspace{0.5em}23}\mathrm{)}\text{;}\hfill \\ {\beta }_3=\text{1/2}\hfill & (\text{Equation\hspace{0.5em}}24).\hfill \end{array}$$(25)

3.2. Goal-based (DWR) error measures

The DWR-AMR strategy aims to minimize the residual in some specific QoI by refining those regions that contribute to the error in the prescribed goal. The research presented in this paper considers only two QoI that are commonly used in nuclear reactor physics applications: volumetric detector responses and the effective neutron multiplication factor, ${\mathrm{k}}_{\text{eff}}.$

For a consistent and similar discretization of both the (semi-)analytical primal and dual problems, the discrete variational bi-linear forms are adjoint-consistent, $a(\widehat{u},\widehat{v})-a^{†}(\widehat{v},\widehat{u})=0$ $\forall \widehat{u},\widehat{v}\in {\widehat{\mathcal{V}}}_{h,p};$ and (quasi-)optimal rates of convergence can be attained for the error measured in some prescribed goal (Hartmann 2007).

3.2.1. Volumetric detector response

The response of a detector uniformly distributed within some given detector material, denoted by the subscript d, is measured as some reaction rate over the subdomain, $V_d:$ (26) $$\mathit{\text{QoI}}(\varphi )=\sum _{g=1}^G{\int }_{V_d}{\Sigma }_d^g(\overrightarrow{r}){\varphi }^g(\overrightarrow{r})\mathit{\text{ dV}}.$$(26)

For example, the reaction rate could be that measured by fission chambers. More generally, the reactions are those characterized by the effective macroscopic cross-section of absorption; and ${\Sigma }_d^g$ is defined such that, for $g=1,\dots ,G:$ (27) $${\Sigma }_d^g(\overrightarrow{r})=\{\begin{array}{cc}{\Sigma }_a^g(\overrightarrow{r}),\hspace{1em}\hfill & \overrightarrow{r}\in V_d\text{;}\hfill \\ 0,\hspace{1em}\hfill & \overrightarrow{r}\notin V_d.\hfill \end{array}$$(27)

3.2.1.1. Fixed-source problems

By the arbitrariness of the adjoint fixed-source and letting $q_x^{†}={\Sigma }_d,$ the QoI is made equal to the adjoint linear functional for some test function $\varphi ;$ i.e., for ${\varphi }^{†}\in \mathcal{V}:$ (28a) $$a^{†}({\varphi }^{†},\varphi )=b^{†}(\varphi ),\hspace{1em}\forall \varphi \in \mathcal{V},$$(28a) (28b) $$=\mathit{\text{QoI}}(\varphi ),\hspace{1em}\forall \varphi \in \mathcal{V}.$$(28b)

The DWR-AMR strategy aims to minimize the absolute error in the goal, $|\mathit{\text{QoI}}(\varphi )-\mathit{\text{QoI}}({\widehat{\varphi }}_{h,p})|,$ Equation (11)(11) $${ϵ}^{(i)}=\Vert \mathit{\text{QoI}}(\varphi )-\mathit{\text{QoI}}({\widehat{\varphi }}_{h,p}^{(i)})\Vert .$$(11) . Again, one can establish an error estimate, $\widehat{ϵ},$ measured with respect to a reference solution, that provides an upper-bound to the true error in the QoI, ϵ (Wilson 2021).

For a given energy group, the reference solutions to the primal and the dual problems are solved on the same reference mesh, ${\mathcal{T}}_{\text{h/2},p\text{+1}}.$ Rather than explicitly calculate four different solutions over two different levels of mesh refinement, only the reference solutions are computed at P-1 in Figure A1. The DWR-AMR strategy is driven by estimating the error in the QoI by interpolation of the reference solutions onto the coarser current meshes; i.e.,: (29a) $${ϵ}^2=\sum _{g=1}^G{({ϵ}^g)}^2,$$(29a) (29b) $$\le \sum _{g=1}^G{\Vert {\varphi }^{g†}-{\widehat{\varphi }}_{h,p}^{g†}\Vert }_{H^1}^2{\Vert {\varphi }^g-{\widehat{\varphi }}_{h,p}^g\Vert }_{H^1}^2,$$(29b) (29c) $$\le \sum _{g=1}^G{\Vert {\varphi }^{g†}-\prod _{h,p}{\varphi }^{g†}\Vert }_{H^1}^2{\Vert {\varphi }^g-\prod _{h,p}{\varphi }^g\Vert }_{H^1}^2,$$(29c) (29d) $$\lesssim \sum _{g=1}^G{\Vert {\widehat{\varphi }}_{\text{h/2},p\text{+1}}^{g†}-\prod _{h,p}{\widehat{\varphi }}_{\text{h/2},p\text{+1}}^{g†}\Vert }_{H^1}^2{\Vert {\widehat{\varphi }}_{\text{h/2},p\text{+1}}^g-\prod _{h,p}{\widehat{\varphi }}_{\text{h/2},p\text{+1}}^g\Vert }_{H^1}^2,$$(29d) (29e) $$\lesssim \sum _{g=1}^G\sum _{V_e\in {\mathcal{T}}_{h,p}^g}{\Vert {\widehat{\varphi }}_{\text{h/2},p\text{+1}}^{g†}-\prod _{h,p}{\widehat{\varphi }}_{\text{h/2},p\text{+1}}^{g†}\Vert }_{H^1(V_e)}^2{\Vert {\widehat{\varphi }}_{\text{h/2},p\text{+1}}^g-\prod _{h,p}{\widehat{\varphi }}_{\text{h/2},p\text{+1}}^g\Vert }_{H^1(V_e)}^2,$$(29e) (29f) $$\lesssim \sum _{g=1}^G\sum _{V_e\in {\mathcal{T}}_{h,p}^g}{\mu }_e^g.$$(29f)

The global error estimate in the QoI, $\widehat{ϵ},$ is now a projection-based error estimate of $ϵ,$ bounded by the sum of element-local interpolation-based error measures, ${\mu }_e^g,$ Equation (18)(18) $$\widehat{ϵ}={\left[\sum _{g=1}^G{({\widehat{ϵ}}^g)}^2\right]}^{\frac{1}{2}}={\left[\sum _{g=1}^G\sum _{V_e\in {\mathcal{T}}_{h,p}^g}{\mu }_e^g\right]}^{\frac{1}{2}}\text{;}$$(18) ; i.e.,: for $g=1,\dots ,G:$ (30) $${\mu }_e^g=|{\widehat{\varphi }}_{\text{h/2},p\text{+1}}^g-\prod _{h,p}{\widehat{\varphi }}_{\text{h/2},p\text{+1}}^g{|}_{H^1(V_e)}^2|{\widehat{\varphi }}_{\text{h/2},p\text{+1}}^{g†}-\prod _{h,p}{\widehat{\varphi }}_{\text{h/2},p\text{+1}}^{g†}{|}_{H^1(V_e)}^2,\hspace{1em}\forall V_e\in {\mathcal{T}}_{h,p}^g(V).$$(30)

Those elements eligible for refinement at D-1 in Figure A1 are determined by Equation (20)(20) $$w^{g(i)}{\mu }_e^{g(i)}\ge {\beta }_1\underset{g\prime =1,\dots ,G^{}}{\mathrm{\text{max}}}\underset{V_{e\prime }\in {\mathcal{T}}_{h,p}^{g\prime (i)}}{\mathrm{\text{max}}}{w}^{g\prime (i)}{\mu }_{e\prime }^{g\prime (i)},\hspace{1em}\forall V_e\in {\mathcal{T}}_{h,p}^{g(i)}(V).$$(20) ; yet, one recognizes that the importance of each energy group is already accounted for by the error measures, Equation (30)(30) $${\mu }_e^g=|{\widehat{\varphi }}_{\text{h/2},p\text{+1}}^g-\prod _{h,p}{\widehat{\varphi }}_{\text{h/2},p\text{+1}}^g{|}_{H^1(V_e)}^2|{\widehat{\varphi }}_{\text{h/2},p\text{+1}}^{g†}-\prod _{h,p}{\widehat{\varphi }}_{\text{h/2},p\text{+1}}^{g†}{|}_{H^1(V_e)}^2,\hspace{1em}\forall V_e\in {\mathcal{T}}_{h,p}^g(V).$$(30) ; and therefore, $w^g=1$ for $g=1,\dots ,G.$ The same options for local refinement are available, as in the case of the NRG-AMR strategy; and the same constraints apply; qq.v. Equation (25)(25) $$\{\begin{array}{cc}{\beta }_1=\text{1/3}\hfill & \mathrm{(}\text{Equation\hspace{0.5em}20}\mathrm{)}\text{;}\hfill \\ {\beta }_2=p_{\mathrm{\text{max}}}\hfill & \mathrm{(}\text{Equation\hspace{0.5em}23}\mathrm{)}\text{;}\hfill \\ {\beta }_3=\text{1/2}\hfill & (\text{Equation\hspace{0.5em}}24).\hfill \end{array}$$(25) .

There does exist a subtle difference in choosing the refinement option that maximizes the rate at which the local contribution to the interpolation-based error estimate decreases, Equation (13): (31) $$\underset{h^{*},p^{*}}{\mathrm{\text{max}}}\overline{\Delta }{ϵ}_{h^{*},p^{*}}^{g(i)}{|}_{V_e}\mathrm{}=\mathrm{}\underset{h^{*},p^{*}}{\mathrm{\text{max}}}\frac{{\mu }_e^{g(i)}-|{\widehat{\varphi }}_{\text{h/2},p\text{+1}}^{g(i)}\mathrm{}-\mathrm{}\prod _{h^{*},p^{*}}\mathrm{}{\widehat{\varphi }}_{\text{h/2},p\text{+1}}^{g(i)}{|}_{H^1(V_e)}|{\widehat{\varphi }}_{\text{h/2},p\text{+1}}^{g†(i)}\mathrm{}-\mathrm{}\prod _{h^{*},p^{*}}\mathrm{}{\widehat{\varphi }}_{\text{h/2},p\text{+1}}^{g†(i)}{|}_{H^1(V_e)}}{|\mathit{\text{Do}}{F}_{h,p}^{g(i)}-\mathit{\text{Do}}{F}_{h^{*},p^{*}}^{g(i)}{|}_{V_e}}.$$(31)

3.2.1.2. Eigenvalue problems

For detector responses in fixed-source problems the QoI is a linear quantity. However, in eigenvalue problems, the QoI becomes a non-linear quantity and the fission sources must be known before the detector response can be determined. Therefore, one must first converge the dual-problem fission sources and the associated eigenvalue. Indeed, one may use the ${\mathrm{k}}_{\text{eff}}$ calculated from the primal problem to expedite the convergence of the dual problem (Wang, Bangerth, and Ragusa 2009; Goffin et al. 2013). Once these have been established, the dual problem can then be re-solved, as a fixed-source problem, with the QoI defined as per Equations (26)(26) $$\mathit{\text{QoI}}(\varphi )=\sum _{g=1}^G{\int }_{V_d}{\Sigma }_d^g(\overrightarrow{r}){\varphi }^g(\overrightarrow{r})\mathit{\text{ dV}}.$$(26) and (28) (Park 2006).

3.2.2. The effective neutron multiplication factor (${\mathrm{k}}_{\text{eff}}$)

For multiplying systems without any extraneous fixed-sources of neutrons, $q_x=0$ and hence $b=0,$ the weak primal problem is posed as finding $\widehat{\varphi }\in {\widehat{\mathcal{V}}}_{h,p}$ such that: (32) $$\tilde{a}(\widehat{\varphi },\widehat{w})=\frac{1}{{\mathrm{k}}_{\text{eff}}}{a}_F(\widehat{\varphi },\widehat{w}),\hspace{1em}\forall \widehat{w}\in {\widehat{\mathcal{V}}}_{h,p}\text{;}$$(32) and similarly, the weak dual problem is posed as finding ${\widehat{\varphi }}^{†}\in {\widehat{\mathcal{V}}}_{h,p}$ such that: (33) $${\tilde{a}}^{†}({\widehat{\varphi }}^{†},\widehat{w})=\frac{1}{{\mathrm{k}}_{\text{eff}}}{a}_F^{†}({\widehat{\varphi }}^{†},\widehat{w}),\hspace{1em}\forall \widehat{w}\in {\widehat{\mathcal{V}}}_{h,p}\text{;}$$(33) where $\tilde{a}$ and ${\tilde{a}}^{†}$ are the original bilinear forms in Equations (14)(14) $$a(\widehat{\varphi },\widehat{w})=b(\widehat{w}),\hspace{1em}\forall \widehat{w}\in \widehat{\mathcal{V}}.$$(14) and (28), respectively, less the fission-source terms; i.e.,: $\tilde{a}=a+a_F$ and ${\tilde{a}}^{†}=a^{†}+a_F^{†},$ respectively.

When solving eigenvalue problems, the solution should be normalized after each outer-iteration to avoid any problems of numerical underflow, or overflow. Moreover, the primal and the dual problems should not be normalized independently. Accordingly, one normalizes the direct solution such that $a_F(\varphi ,\varphi )=1;$ and one normalizes the adjoint solution such that $a_F^{†}({\widehat{\varphi }}^{†},\widehat{\varphi })=a_F(\widehat{\varphi },{\widehat{\varphi }}^{†})=1.$ Now, the QoI can be expressed as (Lathouwers 2011a,b): (34) $$\mathit{\text{QoI}}(\varphi )=\frac{1}{{\mathrm{k}}_{\text{eff}}}=\frac{1}{{\mathrm{k}}_{\text{eff}}}{a}_F^{†}({\widehat{\varphi }}^{†},\widehat{\varphi }).$$(34)

4. The multi-group neutron diffusion equation

The research presented in this paper considers the steady-state, energy-dependent, neutron diffusion equation (NDE). The spatial domain is Discretized using isogeometric analysis (IGA); and the energy domain is Discretized using the multi-group approximation.

The continuous energy domain, $(E_{\mathrm{\text{min}}},E_{\mathrm{\text{max}}}],$ is partitioned into distinct values, $E_{\mathrm{\text{max}}}=E_0>E_1>\dots >E_G=E_{\mathrm{\text{min}}},$ forming G number of energy groups; such that the $g^{\mathit{\text{th}}}$-group spans the half-closed interval $(E_g,E_{g-1}].$ Any integrals over the entire energy spectrum are expressed as a sum of individual integrals over each energy group; i.e.,: (35) $${\int }_0^{\infty }dE\approx {\int }_{E_G}^{E_0}dE\Rightarrow \sum _{g=1}^G{\int }_{E_g}^{E_{g-1}}\mathit{\text{dE}}.$$(35)

The defining equation is integrated over each energy group and the result is solved for the multi-group components of the scalar neutron flux, ${\varphi }^g;$ where: (36) $${\varphi }^g(\overrightarrow{r})={\int }_{E_g}^{E_{g-1}}\varphi (\overrightarrow{r},E)\mathit{\text{ dE}}.$$(36)

4.1. The direct (primal) equation

The direct (primal) multi-group NDE is posed on the closure of a connected and bounded domain, $V\subset {\mathbb{R}}^d,$ in d-dimensional Euclidean space for d = 2, 3. The strong form of the problem is cast as solving ${\varphi }^g:V\to \mathbb{R}$ for $g=1,\dots ,G:$ $$\{\begin{array}{ll}-\overrightarrow{\nabla }·D^g(\overrightarrow{r})\overrightarrow{\nabla }{\varphi }^g(\overrightarrow{r})+{\Sigma }_r^g(\overrightarrow{r}){\varphi }^g(\overrightarrow{r})=q^g(\overrightarrow{r}),\overrightarrow{r}\in V\text{;}& \mathrm{(}\text{37a}\mathrm{)}\hfill \\ {\varphi }^g(\overrightarrow{r})=g_D^g,\overrightarrow{r}\in \partial V_D\text{;}& \mathrm{(}\text{37b}\mathrm{)}\hfill \\ \overrightarrow{\nabla }{\varphi }^g(\overrightarrow{r})·\overrightarrow{n}(\overrightarrow{r})=g_N^g,\overrightarrow{r}\in \partial V_N\text{;}& \mathrm{(}\text{37c}\mathrm{)}\hfill \end{array}$$ where Dirichlet and Neumann-type BCs, Equations (37b)(25) $$\{\begin{array}{cc}{\beta }_1=\text{1/3}\hfill & \mathrm{(}\text{Equation\hspace{0.5em}20}\mathrm{)}\text{;}\hfill \\ {\beta }_2=p_{\mathrm{\text{max}}}\hfill & \mathrm{(}\text{Equation\hspace{0.5em}23}\mathrm{)}\text{;}\hfill \\ {\beta }_3=\text{1/2}\hfill & (\text{Equation\hspace{0.5em}}24).\hfill \end{array}$$(25) and (37c)(25) $$\{\begin{array}{cc}{\beta }_1=\text{1/3}\hfill & \mathrm{(}\text{Equation\hspace{0.5em}20}\mathrm{)}\text{;}\hfill \\ {\beta }_2=p_{\mathrm{\text{max}}}\hfill & \mathrm{(}\text{Equation\hspace{0.5em}23}\mathrm{)}\text{;}\hfill \\ {\beta }_3=\text{1/2}\hfill & (\text{Equation\hspace{0.5em}}24).\hfill \end{array}$$(25) , respectively, are prescribed on the external boundaries; i.e., $\partial V_{\mathit{\text{ext}}}=\partial V_D\cup \partial V_N,$ respectively; and the normal vector, $\overrightarrow{n},$ is outward-pointing. The defining equation can be recast in block-like group operator notation: (38a) $$L\varphi =q,$$(38a) (38b) $$=S\varphi +F\varphi +q_x\text{;}$$(38b) with the scalar neutron flux, $\varphi ={[{\varphi }^1,\dots ,{\varphi }^G]}^T;$ the extraneous fixed-sources, $q_x={[q_x^1,\dots ,q_x^G]}^T;$ the diagonal loss operator, L, which comprises within-group losses due to diffusive-leakage and total collision, $L^g,$ Equation (39a)(39a) $$L^g\varphi :=-\overrightarrow{\nabla }·D^g\overrightarrow{\nabla }{\varphi }^g+{\Sigma }_r^g{\varphi }^g\text{;}$$(39a) ; the scatter operator, S, which comprises group-to-group scatter terms, $S^g,$ Equation (39b)(39b) $$S^g\varphi :=\sum _{g\prime \ne g}{\Sigma }_s^{g\prime \to g}{\varphi }^{g\prime }\text{;}$$(39b) ; and the fission operator, F, which comprises multi-group fission terms, $F^g,$ Equation (39c)(39c) $$F^g\varphi :={\chi }^g\sum _{g\prime }{\nu }^{g\prime }{\Sigma }_f^{g\prime }{\varphi }^{g\prime }.$$(39c) . (39a) $$L^g\varphi :=-\overrightarrow{\nabla }·D^g\overrightarrow{\nabla }{\varphi }^g+{\Sigma }_r^g{\varphi }^g\text{;}$$(39a) (39b) $$S^g\varphi :=\sum _{g\prime \ne g}{\Sigma }_s^{g\prime \to g}{\varphi }^{g\prime }\text{;}$$(39b) (39c) $$F^g\varphi :={\chi }^g\sum _{g\prime }{\nu }^{g\prime }{\Sigma }_f^{g\prime }{\varphi }^{g\prime }.$$(39c)

Defining the operator, A, such that: (40) $$A\varphi =(L-S)\varphi ,$$(40) one can distinguish fixed-source problems: (41) $$(A-F)\varphi =q_x,$$(41) from eigenvalue ones: (42) $$A\varphi =\frac{1}{{\mathrm{k}}_{\text{eff}}}F\varphi .$$(42)

4.2. The adjoint (dual) Equation

Mathematically, the dual problem should be adjoint-consistent with the primal with respect to the discretization of the spatial and energy domains. By linearity of the inner-product, adjoint-consistency can be satisfied in each of the operators defined in Equation (39); i.e., (Lewis and Miller Jr 1984): (43) $$\sum _{g=1}^G[{\int }_V{\varphi }^{g†}(\overrightarrow{r})B^g{\varphi }^g(\overrightarrow{r})\mathit{\text{ dV}}]\equiv \sum _{g=1}^G[{\int }_V{\varphi }^g(\overrightarrow{r})B^{g†}{\varphi }^{g†}(\overrightarrow{r})\mathit{\text{ dV}}].$$(43)

The dual problem is posed in the same domain, $V\in {\mathbb{R}}^d,$ as the primal one; its strong form cast in multi-group operator notation as finding ${\varphi }^{g†}:V\to \mathbb{R}$ for $g=G,\dots ,1:$ (44a) $$L^{†}{\varphi }^{†}=q^{†},$$(44a) (44b) $$=S^{†}{\varphi }^{†}+F^{†}{\varphi }^{†}+q_x^{†}\text{;}$$(44b) with the scalar neutron importance, ${\varphi }^{†}={[{\varphi }^{1†},\dots ,{\varphi }^{G†}]}^T;$ the adjoint source, $q_x^{†}={[q_x^{1†},\dots ,q_x^{G†}]}^T;$ and the loss, scatter and fission adjoint operators, $L^{†},S^{†}$ and $F^{†},$ respectively: (45a) $$L^{g†}{\varphi }^{†}:=-\overrightarrow{\nabla }·D^g\overrightarrow{\nabla }{\varphi }^{g†}+{\Sigma }_r^g{\varphi }^{g†}\text{;}$$(45a) (45b) $$S^{g†}{\varphi }^{†}:=\sum _{g\prime \ne g}{\Sigma }_s^{g\to g\prime }{\varphi }^{g\prime }\text{;}$$(45b) (45c) $$F^{g†}{\varphi }^{†}:={\nu }^g{\Sigma }_f^g\sum _{g\prime }{\chi }^{g\prime }{\varphi }^{g\prime }.$$(45c)

The scatter and the fission operators are not self-adjoint; their adjoint operators are obtained by the transposition and the interchanging of arguments of the energy-integral kernels (Schunert et al. 2015). The loss operator is formally self-adjoint, $(L\varphi ,{\varphi }^{†})=(\varphi ,L{\varphi }^{†});$ however, the diffusive-leakage term necessitates an appropriate prescription of adjoint BCs. It simplifies the problem to restrict the general Dirichlet and Neumann boundary conditions (BCs), Equations (37b) and (37c), respectively, to the set of zero-flux, reflective and vacuum BCs, prescribed on the external boundaries, $\partial V_0,\partial V_R$ and $\partial V_V$ respectively: $$\{\begin{array}{ll}{\varphi }^g(\overrightarrow{r})=0,\overrightarrow{r}\in \partial V_0\text{;}& \mathrm{(}\text{46a}\mathrm{)}\hfill \\ -D^g(\overrightarrow{r})\overrightarrow{\nabla }{\varphi }^g(\overrightarrow{r})·\overrightarrow{n}=0,\overrightarrow{r}\in \partial V_R\text{;}& \mathrm{(}\text{46b}\mathrm{)}\hfill \\ -D^g(\overrightarrow{r})\overrightarrow{\nabla }{\varphi }^g(\overrightarrow{r})·\overrightarrow{n}=\frac{1}{2}{\varphi }^g(\overrightarrow{r}),\overrightarrow{r}\in \partial V_V.& (46c)\hfill \end{array}$$

The corresponding set of adjoint BCs are prescribed over the same boundaries: $$\{\begin{array}{ll}{\varphi }^{g†}(\overrightarrow{r})=0,\overrightarrow{r}\in \partial V_0\text{;}& \mathrm{(}\text{47a}\mathrm{)}\hfill \\ D^g(\overrightarrow{r})\overrightarrow{\nabla }{\varphi }^{g†}(\overrightarrow{r})·\overrightarrow{n}=0,\overrightarrow{r}\in \partial V_R\text{;}& \mathrm{(}\text{47b}\mathrm{)}\hfill \\ D^g(\overrightarrow{r})\overrightarrow{\nabla }{\varphi }^{g†}(\overrightarrow{r})·\overrightarrow{n}=\frac{-1}{2}{\varphi }^{g†}(\overrightarrow{r}),\overrightarrow{r}\in \partial V_V\text{;}& \mathrm{(}\text{47c}\mathrm{)}\hfill \end{array}$$ where the vacuum BC has been considered as a general Robin BC of the form: $\alpha \varphi +\overrightarrow{\nabla }\varphi ·\overrightarrow{n}=g,$ which assumes that the incoming partial current of neutrons, as well as the outgoing partial current of neutron importance, are both negligibly small.

In a similar manner to Equation (40)(40) $$A\varphi =(L-S)\varphi ,$$(40) for the primal problem, the operator $A^{†}$ is defined such that: (48) $$A^{†}{\varphi }^{†}=(L^{†}-S^{†}){\varphi }^{†}.$$(48)

Again, one distinguishes fixed-source problems: (49) $$(A^{†}-F^{†}){\varphi }^{†}=q_x^{†},$$(49) from eigenvalue ones: (50) $$A^{†}{\varphi }^{†}=\frac{1}{{\mathrm{k}}_{\text{eff}}^{†}}{F}^{†}{\varphi }^{†}.$$(50)

Under the assumption that the property of (semi-discrete) adjoint-consistency holds, it should be clear that ${\mathrm{k}}_{\text{eff}}^{†}={\mathrm{k}}_{\text{eff}};$ i.e., the Rayleigh quotient: (51) $${\mathrm{k}}_{\text{eff}}^{†}=\frac{(\varphi ,F^{†}{\varphi }^{†})}{(\varphi ,A^{†}{\varphi }^{†})}=\frac{(F\varphi ,{\varphi }^{†})}{(A\varphi ,{\varphi }^{†})}={\mathrm{k}}_{\text{eff}}\text{;}$$(51) where $|1/{\mathrm{k}}_{\text{eff}}|$ is determined after each ${(i)}^{\mathit{\text{th}}}$ power-iteration: (52a) $${\mathrm{k}}_{\text{eff}}^{(i)}=\frac{(F{\varphi }^{(i)},{\varphi }^{†})}{(A{\varphi }^{(i)},{\varphi }^{†})}=\frac{(F{\varphi }^{(i)},{\varphi }^{†})}{\frac{1}{k_{\text{eff}}^{(i+1)}}(F{\varphi }^{(i+1)},{\varphi }^{†})}\text{;}$$(52a) (52b) $$\Rightarrow \frac{1}{{\mathrm{k}}_{\text{eff}}^{(i+1)}}=\frac{1}{{\mathrm{k}}_{\text{eff}}^{(i)}}\frac{(F{\varphi }^{(i)},{\varphi }^{†})}{(F{\varphi }^{(i+1)},{\varphi }^{†})}\text{;}$$(52b) where ${\varphi }^{†}$ is not obliged to be the solution to the adjoint problem in the above iteration; and $\underset{(i)\to \infty }{\mathrm{\text{lim}}}{\mathrm{k}}_{\text{eff}}^{(i)}={\mathrm{k}}_{\text{eff}}.$

5. A continuous Bubnov-Galerkin isogeometric discretization

To start, one introduces some abstract function spaces. The Lebesgue space, $L_2(V),$ of some measurable functions, u, that are square-integrable on the bounded domain, V, is defined as: (53) $$L_2(V):=\{u\text{\hspace{0.5em}measureable}:{\Vert u\Vert }_{\mathit{\text{LE}}(V)}<\infty \}.$$(53)

This is a real linear Hilbert space furnished with an inner-product: (54) $${(u,v)}_{L_2(V)}:={\int }_{V}vu\mathit{\text{ dV;}}$$(54) and a norm: (55) $${\Vert u\Vert }_{L_2(V)}:={(u,u)}_{L_2(V)}^{\text{1/2}}.$$(55)

The s-order Sobolev space, $H^s(V)\subset L_2(V),$ is defined for those functions belonging to $L_2(V),$ that also admit all weak partial derivatives, $D^{\alpha },$ up to order s, belonging to $L_2(V);$ for $\overrightarrow{r}=(x_1,x_2,x_3)=(x,y,z)$ and $D^{\alpha }=\prod {}_{i=1}^d{D}_i^{{\alpha }_i}$ for $D_i^{{\alpha }_i}={\partial }^{{\alpha }_i}/{\partial }_{x_i}^{{\alpha }_i};$ and $\left|\alpha \right|=\sum {}_i^d{\alpha }_i$ for some multi-index $\alpha \in {\mathbb{R}}^d:$ (56) $$H^s(V):=\{u\in L_2(V):D^{\alpha }u\in L_2(V),\forall 0\le |\alpha |\le s\}.$$(56)

The Sobolev space, $H^s(V),$ is a Hilbert space equipped with an inner-product: (57) $${(u,v)}_{H^s(V)}:=\sum _{\left|\alpha \right|\le s}{(D^{\alpha }v,D^{\alpha }u)}_{L_2(V)}\text{;}$$(57) and norm: (58) $${\Vert u\Vert }_{H^s(V)}:={(u,u)}_{H^s(V)}^{\text{1/2}}.$$(58)

The restriction of such functions and their normal derivatives with support over the domain, V, are defined over the boundary, $\partial V,$ by trace operators, ${\gamma }_i:H^s(V)\to H^{s-i-1/2}(\partial V)$ $\forall 0\le i<s;$ viz. the Dirichlet trace, ${\gamma }_{0}u=u{|}_{\partial V},$ is well defined in $L_2(\partial V)$ for $s>\text{3/2}$ (Oden and Reddy 2011).

Now, the point of departure is the semi-discrete multi-group NDE, Equation (38), and its adjoint, Equation (44). Both the primal and the dual problem constitute systems of equations that are solved for the multi-group components ${\varphi }^g:V\to \mathbb{R}$ and ${\varphi }^{g†}:V\to \mathbb{R},$ respectively; with BCs, such as those in Equations (46) and (47), respectively, prescribed over separate sections of the external boundary.

Without loss of generality, an emphasis is placed upon manipulating the direct equation; and the same discretization is applied to the adjoint equation. The strong form of the defining equations is reformulated into a suitable abstract variational boundary value problem (VBVP), or weak form, which itself is obtained by multiplying Equation (38) by a test function, $w\in \mathcal{V}(V),$ and then by integrating over the entire domain in the real physical space, V. For $q_x\in L_2(V),$ rather than searching for a smooth solution in the space of continuous functions, a solution is sought in the s-order Sobolev space, Equation (56)(56) $$H^s(V):=\{u\in L_2(V):D^{\alpha }u\in L_2(V),\forall 0\le |\alpha |\le s\}.$$(56) , where derivatives are understood in a weak distributional sense. The (semi-)analytical solution is approximated by some infinite linear expansion of trial functions, $u(\overrightarrow{r})\in \mathcal{U}(V)\subseteq H^s(V),$ and control-variables, ${\varphi }_i:$ (59) $$\varphi (\overrightarrow{r})\approx \sum _{i=1}^{\infty }{\varphi }_i{u}_i(\overrightarrow{r}),\hspace{1em}\forall \overrightarrow{r}\in V.$$(59)

In the integral statement defined over V, one substitutes the approximation above into Equation (38); this yields a weighted-residual, which is made to be negligibly small everywhere. Thus, the residual is made orthogonal to every member of the set of test functions. The weighted-residual expression necessitates that the trial functions are strictly second-order differentiable; this is requirement is relaxed by application of the divergence theorem to the double-differential diffusion term. Now, the test and trial functions are required to be first-order differentiable; and as per a Bubnov-Galerkin approach, they are both taken from the same space and the residual is made orthogonal to the numerical solution; i.e., $\mathcal{U}(V),\mathcal{V}(V)\subseteq H^1(V)$ (Finlayson 1972).

The weak form naturally satisfies the reflective BC, Equation (46b); whereas the vacuum BC is implemented by substitution of Equation (46c); and Nitsche’s penalty method is used to weakly enforces the zero-flux BC, Equation (46a), at the Dirichlet boundaries via the following penalization, for ${\mu }_f>0$ (Nitsche 1971): (60) $$\sum _{\partial V_f\in \partial V_D}{\mu }_f{\int }_{\partial V_f}w(\overrightarrow{r})\varphi (\overrightarrow{r})\mathit{\text{ dS}}=\sum _{\partial V_f\in \partial V_D}{\mu }_f{\int }_{\partial V_f}w(\overrightarrow{r})g_D\mathit{\text{ dS}}.$$(60)

The domain in the real physical space, V, is decomposed into E number of non-degenerate, non-overlapping NURBS patches, $V_e;$ i.e., $V=\cup {}_{e=1}^E{V}_e.$ Within any given patch, the material is assumed to be homogeneous and have uniform properties; viz. the nuclear data is assigned as piecewise constants; e.g. ${\Sigma }_{r,e}={\Sigma }_r(\overrightarrow{r}){|}_{V_e},$ $\forall \overrightarrow{r}\in V_e\in \mathcal{T}.$

The union of macro-elements forms the general physical mesh, $\mathcal{T}(V)={\left\{V_e\right\}}_{e=1}^E.$ For simplicity of developing the spatial discretization, one assumes a single energy-independent mesh. The intersection of two patches belonging to the same physical mesh forms either an empty set (disjoint); an edge shared by indirect neighbors, which is actually a NURBS vertex in 2D; or a face shared by direct neighbors, which is actually a NURBS edge in 2D. A patch interface formed at the boundaries of two adjoining patches that share a common face belongs to the set of interior faces, $\partial V_{\mathit{\text{int}}};$ and therefore, an arbitrary patch face, $\partial V_f,$ belongs to the complete set $\partial V_f\in \partial V=\partial V_{\mathit{\text{int}}}\cup \partial V_{\mathit{\text{ext}}}.$

The integral terms over the entire domain are obtained as a summation of dense local contributions from over each patch, $V_e\in \mathcal{T},$ in a systematic assembly process. A virgin patch, with no internal refinement structure, can be decomposed into more non-vanishing knot-spans, $V_{\widehat{e}},$ via knot insertion; i.e., $V_e=\cup {}_{\widehat{e}}{V}_{\widehat{e}}.$ Each new micro-element retains its own homogeneity and well-defined material properties, as inherited from the macro-element; viz. $\Sigma {|}_{V_{\widehat{e}}}=\Sigma {|}_{V_e}$ $\forall V_{\widehat{e}}\in V_e.$ Integration is always performed over the smallest entity; so more rigorously, the assembly process loops over all non-vanishing knot-spans belonging to the mesh. Without loss of generality, however, the subscript e is retained, as pertaining to the element over which the integral quantity is evaluated; and the distinction between macro-element and the micro-element is made when necessary.

The weak form of Equation (41)(41) $$(A-F)\varphi =q_x,$$(41) is to find $\varphi \in \mathcal{V}(V)\subseteq H^1(V),$ such that: (61) $$a(\varphi ,w)=b(w),\hspace{1em}\forall w\in \mathcal{V}(V).$$(61)

The bilinear form, $a(·,·):\mathcal{V}(V)\times \mathcal{V}(V)\to \mathbb{R},$ retains its block-like structure and comprises a weak diffusive-leakage and group-removal term, $a_L,$ defined over volumes; a group-to-group scattering term, $a_S;$ a fission term, $a_F;$ as well as a term that incorporates the BCs prescribed on the external boundary, $a_{\mathit{\text{ext}}}:$ (62) $$a(\varphi ,w)=a_L(\varphi ,w)-a_S(\varphi ,w)-a_F(\varphi ,w)+a_{\mathit{\text{ext}}}(\varphi ,w).$$(62)

Suppressing any spatial dependence, each of these is defined below, respectively: (63a) $$\begin{array}{c}\begin{array}{c}{a}_L(\varphi ,w)=\sum _{g=1}^G[\sum _{V_e\in \mathcal{T}}{\int }_{V_e}\overrightarrow{\nabla }w·D_e^g\overrightarrow{\nabla }{\varphi }^g+w{\Sigma }_{r,e}^g{\varphi }^g\mathit{\text{ dV}}],\hfill \\ =\sum _{g=1}^G{a}_L^g({\varphi }^g,w)\text{;}\hfill \end{array}\\ \end{array}$$(63a) (64b) $$\begin{array}{c}\begin{array}{c}{a}_S(\varphi ,w)=\sum _{g=1}^G\left[\sum _{g\prime \ne g}^G\sum _{V_e\in \mathcal{T}}{\int }_{V_e}w{\Sigma }_{s,e}^{g\prime \to g}{\varphi }^{g\prime }\mathit{\text{ dV}}\right],\hfill \\ =\sum _{g=1}^G{a}_S^g(\varphi ,w)\text{;}\hfill \end{array}\\ \end{array}$$(64b) (63c) $$\begin{array}{c}\begin{array}{c}{a}_F(\varphi ,w)=\sum _{g=1}^G\left[\sum _{g\prime =1}^G\sum _{V_e\in \mathcal{T}}{\chi }_e^g{\int }_{V_e}w{\nu }_e^{g\prime }{\Sigma }_{f,e}^{g\prime }{\varphi }^{g\prime }\mathit{\text{ dV}}\right],\hfill \\ =\sum _{g=1}^G{a}_F^g(\varphi ,w)\text{;}\hfill \end{array}\\ \end{array}$$(63c) (63d) $$\begin{array}{c}{a}_{\mathit{\text{ext}}}(\varphi ,w)=\sum _{g=1}^G[\sum _{V_e\in \mathcal{T}}\sum _{\partial V_f\in V_e\cap \partial V_0}\mathrm{}\mathrm{}{\mu }_f{\int }_{\partial V_f}w{\varphi }^g\mathit{\text{ dS}}+\mathrm{}\mathrm{}\sum _{\partial V_f\in V_e\cap \partial V_V}\mathrm{}\frac{1}{2}{\int }_{\partial V_f}w{\varphi }^g\mathit{\text{ dS}}],\hfill \\ =\sum _{g=1}^G{a}_{\mathit{\text{ext}}}^g({\varphi }^g,w).\hfill \end{array}$$(63d)

The linear functional, $b(·):\mathcal{V}(V)\to \mathbb{R},$ comprises only fixed-sources, since $g_D=0:$ (64) $$\begin{array}{c}b(w)=\sum _{g=1}^G\left[\sum _{V_e\in \mathcal{T}}{\int }_{V_e}w{q}_{x,e}^g\mathit{\text{ dV}}\right],\hfill \\ =\sum _{g=1}^G{b}^g(w).\hfill \end{array}$$(64)

For multiplying systems, the discrete weak statement of Equation (42)(42) $$A\varphi =\frac{1}{{\mathrm{k}}_{\text{eff}}}F\varphi .$$(42) is to find $\varphi \in \mathcal{V}(V),$ such that: (65) $$\tilde{a}(\varphi ,w)=\frac{1}{{\mathrm{k}}_{\text{eff}}}{a}_f(\varphi ,w),\hspace{1em}\forall w\in \mathcal{V}(V)\text{;}$$(65) where $\tilde{a}$ is the original bilinear form in Equation (62)(62) $$a(\varphi ,w)=a_L(\varphi ,w)-a_S(\varphi ,w)-a_F(\varphi ,w)+a_{\mathit{\text{ext}}}(\varphi ,w).$$(62) , a, less the fission term; i.e.,: (66) $$\tilde{a}=a+a_F=a_L-a_S+a_{\mathit{\text{ext}}}.$$(66)

The discrete variational problem is formulated by restricting the function space to a finite-dimensional subset; i.e., ${\widehat{\mathcal{V}}}_{h,p}(V)\subset \mathcal{V}(V).$ As per the isoparametric concept, this space is chosen as that spanned by the same N number of NURBS basis functions used to describe the mesh, $\mathcal{T};$ where the NURBS basis functions, R, Equation (3)(3) $$R_{i,j}^{p,q}(\xi ,\eta )=\frac{N_i^p(\xi )M_j^q(\eta )W_{i,j}}{\sum _{i\prime =1}^N\sum _{j\prime =1}^M{N}_{i\prime }^p(\xi )M_{j\prime }^q(\eta )W_{i\prime ,j\prime }}.$$(3) , are characterized by the mesh spacing, h, as well as the polynomial degree of the basis, p: (67) $$\{\begin{array}{c}{[{\widehat{\mathcal{V}}}_{h,p}(V)]}^T=\text{span}{\hspace{0.5em}\{R_i^p\}}_{i=1}^N,\hfill \\ \text{with\hspace{0.5em}supp\hspace{0.5em}}(R^p)=\{\overrightarrow{r}\in V:R_i^p(\overrightarrow{r})\ne 0\text{\hspace{0.5em} for \hspace{0.5em}}i=1,\dots ,N\}.\hfill \end{array}$$(67)

Indeed, the tensor-product structure of the parametric domain permits the NURBS basis functions to be constructed from different polynomials in different parametric directions; however, for simplicity, they are characterized by a single p and unless necessary, the superscript is suppressed. Accordingly, Equation (59)(59) $$\varphi (\overrightarrow{r})\approx \sum _{i=1}^{\infty }{\varphi }_i{u}_i(\overrightarrow{r}),\hspace{1em}\forall \overrightarrow{r}\in V.$$(59) is further approximated by truncating the expansion to the $N^{\mathit{\text{th}}}$ term: (68) $$\varphi (\overrightarrow{r})\approx \widehat{\varphi }(\overrightarrow{r})=\sum _{i=1}^N{\varphi }_i{R}_i(\overrightarrow{r}).$$(68)

One chooses ${\widehat{\mathcal{V}}}_{h,p}(V)$ such that the weighted-residual becomes negligibly small with respect to variations in the control-variables, ${\varphi }_i,$ as the dimension of the finite-dimensional space increases; i.e., $\widehat{\varphi }\to \varphi$ as $N\to \infty .$ The weighted-residual statement must be satisfied for each member in the weighting set; and since the space of trial and test functions coincide, this forms the same number of algebraic equations as there are unknowns per energy group; i.e., $\mathit{\text{Do}}{F}^g=\mathrm{\text{dim}}{\widehat{\mathcal{V}}}_{h,p}(V)=N$ for $g=1,\dots ,G;$ and hence, $\mathit{\text{DoF}}=\sum {}_g\mathit{\text{Do}}{F}^g=G\times N.$ Consequently, for the analysis over a single energy-independent mesh, the discrete problem is posed as a system of ($\mathit{\text{GN}}\times \mathit{\text{GN}}$)-matrix equations.

The space ${\widehat{\mathcal{V}}}_{h,p}(V)$ is the collection of global functions whose restriction to a patch lives in the tensor-product parametric domain ${\widehat{V}}_e.$ Consequently, NURBS basis functions in the real physical space also have local support over a given patch; i.e., ${\widehat{\mathcal{V}}}_{h,p}^e(V_e)\subset {\widehat{\mathcal{V}}}_{h,p}(V)$ $\forall V_e\in \mathcal{T}.$ This piecewise definition assigns the functions a narrow base and results in the banded structure of the assembled matrices. It also introduces discontinuity into the basis; yet, the $H^1(V)$ space demands strict conformity of the numerical solution; viz. continuity across V. Dirichlet traces are well defined in $L_2(\partial V_e)$ for $s>\text{1/2};$ i.e., ${\gamma }_0:H^1(V_e)\to {\mathrm{H}}^{\text{1/2}}(\partial V_e);$ and the extension of functions defined up to the boundary of a patch are made continuous across an interface by virtually gluing them together. This is trivial in the case of regular meshes where the control-variables at patch interfaces exist in a one-to-one correspondence; i.e., $R(\overrightarrow{r}){|}_{V_e}=R\prime (\overrightarrow{r}){|}_{V_{e\prime }}$ $\forall \overrightarrow{r}\in \partial V_f\in V_e\cap V_{e\prime }$ and $\forall V_e,V_{e\prime }\in \mathcal{T}.$

The bilinear form, a, is said to possess the properties of continuity, or boundedness, and coercivity, or ellipticity. The latter is often guaranteed by nature of the problem; viz. physical non-negative nuclear cross-section data. As per the Lax-Milgram Lemma, these properties are sufficient conditions to assert well-posedness; viz. the admission of a unique and stable solution (Ern and Guermond 2004). Since ${\widehat{\mathcal{V}}}_{h,p}(V)\subset \mathcal{V},$ the discrete problem inherits well-posedness from the (semi-)analytical weak-problem (Oden and Reddy 2011). Céa’s Lemma of quasi-best approximation for conforming and consistent discretizations, Equation (69), states that up to a constant, as those defined by continuity and coercivity, $C_{\alpha }$ and $C_{\beta },$ respectively, the set ${\widehat{\mathcal{V}}}_{h,p}(V)$ is that chosen from all possible finite sets of function spaces that minimizes the error in the numerical solution, $\widehat{\varphi }\in {\widehat{\mathcal{V}}}_{h,p}(V),$ with respect to the true solution, $\varphi \in \mathcal{V}(V),$ as measured in the energy-norm; i.e., for some $\widehat{v}\in {\widehat{\mathcal{V}}}_{h,p}(V)$ (Donea and Huerta 2003): (69a) $$C_{\beta }{\Vert \varphi -\widehat{\varphi }\Vert }^2\le a(\varphi -\widehat{\varphi },\varphi -\widehat{\varphi }),$$(69a) (69b) $$\le a(\varphi -\widehat{\varphi },\varphi -\widehat{v})+a(\varphi -\widehat{\varphi },\widehat{v}-\widehat{\varphi }),$$(69b) (69c) $$\le a(\varphi -\widehat{\varphi },\varphi -\widehat{v}),$$(69c) (69d) $$\le C_{\alpha }\Vert \varphi -\widehat{v}\Vert \Vert \varphi -\widehat{\varphi }\Vert \text{;}$$(69d) (69e) $$\Rightarrow \Vert \varphi -\widehat{\varphi }\Vert \le \frac{C_{\alpha }}{C_{\beta }}\underset{\widehat{v}\in {\widehat{\mathcal{V}}}_{h,p}(V)}{\mathrm{\text{inf}}}\Vert \varphi -\widehat{v}\Vert \text{;}$$(69e) where the first line, Equation (69a)(69a) $$C_{\beta }{\Vert \varphi -\widehat{\varphi }\Vert }^2\le a(\varphi -\widehat{\varphi },\varphi -\widehat{\varphi }),$$(69a) , expresses in the bilinear form, a, the property of coercivity. The linearity of each argument of a is exhibited in the second line by introducing $\widehat{v};$ the second term is made to be negligibly small because the error in the numerical solution is orthogonal to the chosen subspace; i.e.,, since ${\widehat{\mathcal{V}}}_{h,p}(V)\subset \mathcal{V}(V):$ (70a) $$a(\varphi -\widehat{\varphi },w)=a(\varphi ,w)-a(\widehat{\varphi },w),\hspace{1em}\forall w\in \mathcal{V}(V),$$(70a) (70b) $$=b(w)-b(w),\hspace{1em}\forall w\in \mathcal{V}(V),$$(70b) (70c) $$=0.$$(70c)

The fourth line, Equation (69d)(69d) $$\le C_{\alpha }\Vert \varphi -\widehat{v}\Vert \Vert \varphi -\widehat{\varphi }\Vert \text{;}$$(69d) , is obtained by expression of continuity in the bilinear form in the third line. Finally, Céa’s Lemma is obtained, Equation (69e)(69e) $$\Rightarrow \Vert \varphi -\widehat{\varphi }\Vert \le \frac{C_{\alpha }}{C_{\beta }}\underset{\widehat{v}\in {\widehat{\mathcal{V}}}_{h,p}(V)}{\mathrm{\text{inf}}}\Vert \varphi -\widehat{v}\Vert \text{;}$$(69e) , by taking the infimum over all possible $\widehat{v}\in {\widehat{\mathcal{V}}}_{h,p}(V).$ A strong definition of consistency has been enforced; yet, Legendre-Gauss quadrature rules are employed to evaluate all integral quantities. Consequently, one expects the discretization to be asymptotically consistent.

Weakly coupled in energy, the within-group equations can be solved sequentially using the most up-to-date information for the group-to-group sources; whereby, the multi-group source terms can be treated as linear functionals that can be moved to the right-hand side of the equation. In the presence of extraneous fixed-sources, ${\widehat{\varphi }}^g\in {\widehat{\mathcal{V}}}_{h,p}(V)$ is sought for $g=1,\dots ,G:$ (71) $$a_L^g({\widehat{\varphi }}^g,\widehat{w})+a_{\mathit{\text{ext}}}^g({\widehat{\varphi }}^g,\widehat{w})=a_S^g(\widehat{\varphi },\widehat{w})+a_F^g(\widehat{\varphi },\widehat{w})+b^g(\widehat{w}),\hspace{1em}\forall \widehat{w}\in {\widehat{\mathcal{V}}}_{h,p}(V).$$(71)

Since the left-hand side, $a_L^g+a_{\mathit{\text{ext}}}^g,$ is SPD, one can employ efficient iterative smoothers to solve Equation (71)(71) $$a_L^g({\widehat{\varphi }}^g,\widehat{w})+a_{\mathit{\text{ext}}}^g({\widehat{\varphi }}^g,\widehat{w})=a_S^g(\widehat{\varphi },\widehat{w})+a_F^g(\widehat{\varphi },\widehat{w})+b^g(\widehat{w}),\hspace{1em}\forall \widehat{w}\in {\widehat{\mathcal{V}}}_{h,p}(V).$$(71) ; e.g. the preconditioned conjugate gradient (PCG) method.

The eigenvalue problem is cast as finding ${\widehat{\varphi }}^g\in {\widehat{\mathcal{V}}}_{h,p}(V)$ for $g=1,\dots ,G:$ (72) $$a_L^g({\widehat{\varphi }}^g,\widehat{w})+a_{\mathit{\text{ext}}}^g({\widehat{\varphi }}^g,\widehat{w})=a_S^g(\widehat{\varphi },\widehat{w})+\frac{1}{{\mathrm{k}}_{\text{eff}}}{a}_F^g(\widehat{\varphi },\widehat{w}),\hspace{1em}\forall \widehat{w}\in {\widehat{\mathcal{V}}}_{h,p}(V).$$(72)

5.1. Enforcing C⁰-continuity on irregular meshes

For a standard, conforming CBG-IGA discretization, any mesh refinement propagates throughout the domain in order to satisfy conditions of compatibility. This can become computationally demanding over large meshes. T-spline descriptions do overcome the limitation of the tensor-product refinement of NURBS bases. However, their usage becomes difficult to generalize in 3D, especially where additional DoFs may need to be introduced at T-junctions (Bazilevs et al. 2010; Buffa, Sangalli, and Vázquez 2014). Alternatively, seeking a solution in the broken Sobolev space does grant total flexibility in terms of local refinement. However, this comes at the price of significantly more DoFs. Moreover, in the case of an interior penalty scheme for a discontinuous Galerkin discretization, the stability of the analysis would depend upon the assignment of some penalty parameter (Wilson et al. 2018).

Figure 5 depicts a simple one-speed, or one-group (1G), problem, which has been solved for $\widehat{\varphi }\in {\widehat{\mathcal{V}}}_{h,p}(V),$ where $C^0$-continuity has been enforced in a strong sense; and separately for $\widehat{\varphi }\prime \in \widehat{\mathcal{V}}{\prime }_{h,p}(\mathcal{T}),$ where continuity has be enforced in a weak sense. In the case of the latter, the basis is discontinuous and continuity has been enforced across patch interfaces by means of an interior penalty scheme. In the case of the former, constraint equations have been applied over all irregular patch interfaces, $\partial V_{\mathit{\text{sm}}},$ formed between a more-refined (slave-)patch, $V_s,$ and a less-refined (master-)patch, $V_m;$ i.e.,: $$\partial V_{\mathit{\text{sm}}}=\{\partial V_f\in V_s\cap V_m\backslash \varnothing \subseteq \partial V_{\mathit{\text{int}}}:\mathrm{\text{dim}}{\widehat{\mathcal{V}}}_{h,p}(V_s)>\mathrm{\text{dim}}{\widehat{\mathcal{V}}}_{h,p}(V_m),\forall V_s,V_m\in \mathcal{T}(V)\}.$$

Figure 5. A comparison of the solutions (right) of a CBG-IGA discretization and its DG-IGA counterpart along the dashed line at $y=‐0.25$ over an irregular mesh (left) for the example of a a simple test case, which consists of a one-group, ($2\text{cm}\times 2\text{cm}$) homogeneous bare-square with isotropic unit-source, ${\Sigma }_t=1.00{\text{cm}}^{‐1},$ c = 0.7. Continuity in the former is enforced by constraining the hanging-nodes 6 & 8 via a Galerkin projection; whereas continuity in the latter is enforced by means of a penalty scheme. The global enumeration of the control-variables for the CBG discretization is given over the mesh on the left-hand side; where the patch boundaries are highlighted in red. (V. the web-based version for reference to color.)

To reiterate, any single intersection of patches that does not belong to an empty set comprises a maximum of two patches; i.e., for ${\cap }_e{V}_e\backslash \varnothing \ni \partial V_f\in \partial V_{\mathit{\text{int}}}$ then $\text{card}({\cap }_e{V}_e\backslash \varnothing )=2.$ Thus, the set of interfaces includes shared edges. A natural consequence of using open knot-vectors is that the NURBS basis is interpolatory at the NURBS-patch boundaries. As such, the boundary of a NURBS surface is a NURBS edge, which is influenced solely by the NURBS basis functions, and their associated control-points, at that given boundary. These control-points at the boundary form their own subset, ${\mathcal{B}}_f\subset \mathcal{B},$ separate from those that do not influence the boundary, ${\mathcal{B}}_n\subset \mathcal{B},$ such that $\mathcal{B}={\mathcal{B}}_n\cup {\mathcal{B}}_f.$ Accordingly, those functions that have support on the boundary of a slave-patch, $\partial V_f\in V_s\cap \partial V_{\mathit{\text{sm}}},$ span the proper subset ${\widehat{\mathcal{V}}}_{h,p}^s(\partial V_f)\subset {\widehat{\mathcal{V}}}_{h,p}(V_s).$ Similarly, those functions that have finite support on the boundary of a master-patch, $\partial V_f\in V_m\cap \partial V_{\mathit{\text{sm}}},$ span ${\widehat{\mathcal{V}}}_{h,p}^m(\partial V_f)\subset {\widehat{\mathcal{V}}}_{h,p}(V_m).$ One can express this more compactly for $e=s,m$ and $\forall V_e\in \mathcal{T}(V)$ and $\forall \partial V_f\in V_e\cap \partial V_{\mathit{\text{sm}}}:$ (73) $$\{\begin{array}{c}{[{\widehat{\mathcal{V}}}_{h,p}^e(\partial V_f)]}^T=\text{span}{\hspace{0.5em}\{{(R{|}_{\partial V_f})}_i\}}_{i=1}^{N_f^e},\hfill \\ \text{with \hspace{0.5em}supp\hspace{0.5em}}(R{|}_{\partial V_f})=\{\overrightarrow{r}\in \partial V_f:{(R(\overrightarrow{r}){|}_{\partial V_f})}_i\ne 0\text{\hspace{0.5em} for \hspace{0.5em}}i=1,\dots ,N_f^e\}\text{;}\hfill \end{array}$$(73) where the NURBS basis functions defined locally over a given element, $V_e\in \mathcal{T}(V),$ have also been denoted by a similar superscript; i.e.,, by letting $N^e=\mathrm{\text{dim}}{\widehat{\mathcal{V}}}_{h,p}(V_e),$ where ${\widehat{\mathcal{V}}}_{h,p}(V_e)\subset {\widehat{\mathcal{V}}}_{h,p}(V),$ then for $i=1,\dots ,N^e:$ $$R_i^e(\overrightarrow{r})={(R(\overrightarrow{r}){|}_{V_e})}_i,\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\overrightarrow{r}\in V_e.$$

Continuity in the geometry across an irregular patch interface is enforced by applying a change of basis to the linear combination of NURBS basis functions and control-points, such that the control-points on one side of the interface, $\left\{{\overrightarrow{B}}_f^s\right\},$ can be represented by $\left\{{\overrightarrow{B}}_f^m\right\}$ on an adjoining side. As per the isoparametric concept, continuity in the solution field at the irregular interface, $\left\{{\overrightarrow{\varphi }}_f\right\},$ can be constrained in a similar fashion; i.e.,: (74a) $${\underset{¯}{\overrightarrow{B}}}_f^s={\underset{¯}{\underset{¯}{K}}}_{\partial V}^{s\leftarrow m}·{\underset{¯}{\overrightarrow{B}}}_f^m\text{;}$$(74a) (74b) $${\underset{¯}{\varphi }}_f^s={\underset{¯}{\underset{¯}{K}}}_{\partial V}^{s\leftarrow m}·{\underset{¯}{\varphi }}_f^m.$$(74b)

The combined transfer matrix above, ${\underset{¯}{\underset{¯}{K}}}_{\partial V}^{s\leftarrow m},$ is that introduced for the conservative interpolation between volume meshes, Equation (10b)(10b) $$=\underset{(N^t\times N^d)(N^d\times 1)}{{\underset{¯}{\underset{¯}{K}}}_V^{t\leftarrow d}·{\underset{¯}{v}}^d.}$$(10b) . However, in this instance, the $H^1$-Sobolev minimization is performed over an irregular patch interface, $\forall \partial V_f\in \partial V_{\mathit{\text{sm}}},$ or a surface mesh. The interpolation $\prod {}_{h,p}^s{\widehat{v}}^m:{\widehat{\mathcal{V}}}^m(\partial V_f)\to {\widehat{\mathcal{V}}}^s(\partial V_f)$ designates the projection of some function, ${\widehat{v}}^m\in {\widehat{\mathcal{V}}}_{h,p}^m(\partial V_f),$ defined on the surface mesh, $\mathcal{T}(\partial V_f),$ of a master-patch onto that of a slave-patch. One notes that each patch may be characterized by some different combination of h and p. For some $\partial V_f\in \partial V_{\mathit{\text{sm}}},$ the optimum target interpolant, ${\widehat{v}}^s\in {\widehat{\mathcal{V}}}_{h,p}^s(\partial V_f),$ is approximated such that: $${\Vert {\widehat{v}}^s-{\widehat{v}}^m\Vert }_{H^1(\partial V_f)}=\underset{\widehat{v}\in {\widehat{\mathcal{V}}}_{h,p}^s(\partial V_f)}{\mathrm{\text{min}}}{\Vert \widehat{v}-{\widehat{v}}^m\Vert }_{H^1(\partial V_f)}.$$

The basis of the slave-patch should sufficiently represent any function that belongs to that of the master-patch, $\mathrm{\text{dim}}{\widehat{\mathcal{V}}}_{h,p}^s(\partial V_f)>\mathrm{\text{dim}}{\widehat{\mathcal{V}}}_{h,p}^m(\partial V_f).$ However, it is not required that ${\widehat{\mathcal{V}}}_{h,p}^s(\partial V_f)\supset {\widehat{\mathcal{V}}}_{h,p}^m(\partial V_f);$ in which case, the continuity would be enforced exactly. For clarity, the combined surface transfer matrix, for which the interpolation between surface meshes is optimal in the $H^1$-norm, is defined as: $$\underset{(N_f^s\times N_f^m)}{{\underset{¯}{\underset{¯}{K}}}_{\partial V}^{s\leftarrow m}}=\underset{(N_f^s\times N_f^s)}{{\left[{\underset{¯}{\underset{¯}{K}}}_{\partial V}^{\mathit{\text{ss}}}\right]}^{-1}}·\underset{(N_f^s\times N_f^m)}{{\underset{¯}{\underset{¯}{K}}}_{\partial V}^{\mathit{\text{sm}}}\text{;}}$$ with the consistent matrices formed over the slave-patch boundary, $\partial V_s,$ ${\underset{¯}{\underset{¯}{K}}}_{\partial V}^{\mathit{\text{ss}}}={\underset{¯}{\underset{¯}{S}}}_{\partial V}^{\mathit{\text{ss}}}+{\underset{¯}{\underset{¯}{M}}}_{\partial V}^{\mathit{\text{ss}}};$ and the projective matrices formed over the entire irregular interface, $\partial V_{\mathit{\text{sm}}},$ ${\underset{¯}{\underset{¯}{K}}}_{\partial V}^{\mathit{\text{sm}}}={\underset{¯}{\underset{¯}{S}}}_{\partial V}^{\mathit{\text{sm}}}+{\underset{¯}{\underset{¯}{M}}}_{\partial V}^{\mathit{\text{sm}}}.$

Although a surface mesh may exist in d-dimensions in the real physical space, it is mapped from a ($d-1$)-dimensional parametric domain. The gradient of NURBS basis functions over a boundary of the domain in the real physical space is expressed for 2D via: (75) $${\underset{¯}{\underset{¯}{\mathbf{J}}}}_{\partial V}·\overrightarrow{\nabla }{R}_i^p(\xi )=\widehat{\nabla }{R}_i^p(\xi ),\mathrm{\hspace{1em}}\text{for\hspace{0.5em}}i=1,\dots ,n\text{;}$$(75) and the surface Jacobian $(d-1\times d)$-matrix, is defined below: (76) $${\underset{¯}{\underset{¯}{\mathbf{J}}}}_{\partial V}:=\left(\begin{array}{cc}\partial x/\partial \xi & \partial y/\partial \xi \end{array}\right)\text{;}$$(76) where the non-square matrix is inverted by evaluating its pseudo-inverse (Jiao and Heath 2004): (77a) $${\left[{\underset{¯}{\underset{¯}{\mathbf{J}}}}_{\partial V}\right]}^{-1}={\left[{\underset{¯}{\underset{¯}{\mathbf{J}}}}_{\partial V}\right]}^T·{\left[{\underset{¯}{\underset{¯}{\mathbf{J}}}}_{\partial V}\right]}^{-T}·{\left[{\underset{¯}{\underset{¯}{\mathbf{J}}}}_{\partial V}\right]}^{-1},$$(77a) (77b) $$=\underset{(d\times d-1)}{{\left[{\underset{¯}{\underset{¯}{\mathbf{J}}}}_{\partial V}\right]}^T}·\underset{(d-1\times d-1)}{{[{\underset{¯}{\underset{¯}{\mathbf{J}}}}_{\partial V}·{\left[{\underset{¯}{\underset{¯}{\mathbf{J}}}}_{\partial V}\right]}^T]}^{-1}.}$$(77b)

The derivatives of the NURBS basis functions with respect to the boundary must be evaluable in the real physical space, Equation (75)(75) $${\underset{¯}{\underset{¯}{\mathbf{J}}}}_{\partial V}·\overrightarrow{\nabla }{R}_i^p(\xi )=\widehat{\nabla }{R}_i^p(\xi ),\mathrm{\hspace{1em}}\text{for\hspace{0.5em}}i=1,\dots ,n\text{;}$$(75) ; and so, the surface mesh must be non-degenerate and the surface Jacobian must be nonsingular.

Again, the notion of an implied super-mesh is employed to ensure that the integral quantities evaluated over an irregular interface are constructed using quadrature rules that sufficiently sample the bases either side of the interface. This is important for maintaining consistency; and thus, for conserving neutrons. However, the use of $p+1$ quadrature points over element boundaries that contain hanging-nodes would be inadequate because such an integrand is not generally a polynomial function. Nevertheless, the integrand over partial boundaries does remain a continuous function over smooth NURBS objects; and thus, additional quadrature points may still offer improved accuracy (Owens 2017).

The manner in which the constraint equations are incorporated into the globally assembled system is based upon the work of Cottrell et al. which builds upon that of Kagan et al. for restoring continuity of a basis in the presence of hanging-nodes on volumetric NURBS objects (Cottrell, Hughes, and Reali 2007; Kagan, Fischer, and Bar-Yoseph 2003). The coupling relationships depend only upon the mesh itself. Therefore, they can be generated in a pre-processing step along with the other dense element matrices.

This paper proposes a wholly-local AMR-h(p) strategy by considering, at most, C⁰-continuity across element interfaces. The roles of master and slave are assigned to those patches that intersect at an irregular patch interface and the constraint equations are enforced via a Galerkin projection, which is agnostic of the underlying refinement structure of either patch. In the assembly process, the local linear systems are constructed as integral quantities over the elements, which correspond to entire patches, $V_e.$ The discrete system of linear equations local to a given slave-patch, ${\underset{¯}{\underset{¯}{A}}}^s·{\underset{¯}{\varphi }}^s={\underset{¯}{b}}^s$ $\forall V_s\in \mathcal{T}(V),$ for any number of constrained sides, $\forall \partial V_f\in V_s\cap \partial V_{\mathit{\text{sm}}},$ is partitioned into those integral quantities that account for the behavior within the patch; and separately for any interaction with its neighbors: (78a) $${\underset{¯}{\underset{¯}{A}}}_{\mathit{\text{nn}}}^s·{\underset{¯}{\varphi }}_n^s+\sum _j[{\underset{¯}{\underset{¯}{A}}}_{n{f}_j}^s·{\underset{¯}{\varphi }}_{f_j}^s]={\underset{¯}{b}}_n^s\text{;}$$(78a) (78b) $$\sum _{\widehat{ı}}[{\underset{¯}{\underset{¯}{A}}}_{f_{\widehat{ı}n}}^s·{\underset{¯}{\varphi }}_n^s+\sum _{\widehat{\mathit{ȷ}}}[{\underset{¯}{\underset{¯}{A}}}_{f_{\widehat{ı}}{f}_{\widehat{\mathit{ȷ}}}}^s·{\underset{¯}{\varphi }}_{f_{\widehat{ı}}}^s]={\underset{¯}{b}}_{f_{\widehat{ı}}}^s]\text{;}$$(78b) where ${\underset{¯}{\underset{¯}{A}}}_{\mathit{\text{nn}}}^s$ is the left-hand-side partial matrix with coefficients that support only those internal functions that have no support at an irregular interface; ${\underset{¯}{\underset{¯}{A}}}_{\mathit{\text{ff}}}^s$ is the partial matrix with coefficients that support only those functions with support at an irregular interface; and ${\underset{¯}{\underset{¯}{A}}}_{\mathit{\text{nf}}}^e$ and ${\underset{¯}{\underset{¯}{A}}}_{\mathit{\text{fn}}}^e$ are the partial matrices with coefficients that support those functions both internal to the patch and at an irregular interface. The linear functional, ${\underset{¯}{b}}^s,$ and the solution vector, ${\underset{¯}{\varphi }}^s,$ are decomposed in a similar manner. However, the latter must also be constrained to the master-patch that adjoins the slave-patch at a given irregular interface, $\partial V_f\in V_m\cap \partial V_{\mathit{\text{sm}}},$ by Equation (74b)(74b) $${\underset{¯}{\varphi }}_f^s={\underset{¯}{\underset{¯}{K}}}_{\partial V}^{s\leftarrow m}·{\underset{¯}{\varphi }}_f^m.$$(74b) ; i.e.,: (79a) $${\underset{¯}{\underset{¯}{A}}}_{\mathit{\text{nn}}}^s·{\underset{¯}{\varphi }}_n^s+\sum _j[{\underset{¯}{\underset{¯}{A}}}_{n{f}_j}^s·{\underset{¯}{\underset{¯}{K}}}_{\partial V_j}^{s\leftarrow m_j}·{\underset{¯}{\varphi }}_{f_j}^{m_j}]={\underset{¯}{b}}_n^s\text{;}$$(79a) (79b) $$\sum _{\widehat{ı}}[{\underset{¯}{\underset{¯}{A}}}_{f_{\widehat{ı}}n}^s·{\underset{¯}{\varphi }}_n^s+\sum _{\widehat{\mathit{ȷ}}}[{\underset{¯}{\underset{¯}{A}}}_{f_{\widehat{ı}}{f}_{\widehat{\mathit{ȷ}}}}^s·{\underset{¯}{\underset{¯}{K}}}_{\partial V_{\widehat{\mathit{ȷ}}}}^{s\leftarrow m_{\widehat{\mathit{ȷ}}}}·{\underset{¯}{\varphi }}_{f_{\widehat{\mathit{ȷ}}}}^{m_{\widehat{\mathit{ȷ}}}}]={\underset{¯}{b}}_{f_{\widehat{ı}}}^s].$$(79b)

The above equations are tested against the larger space of the slave-patch. Yet, the trial functions that belong to the same space are restricted to act only in a linear combination, as per the Galerkin projection, that results in those functions existing in the function space as defined on the master-patch. Without over or under-constraining the DoFs at an irregular interface, so too must those test functions defined at the slave-patch boundary act only in a linear combination, such as to replicate those functions in the weighting space as defined on the master-patch. Therefore, those equations over irregular interfaces must be pre-multiplied by ${[{\underset{¯}{\underset{¯}{K}}}_{\partial V_f}^{s\leftarrow m}]}^T:$ (80a) $${\underset{¯}{\underset{¯}{A}}}_{\mathit{\text{nn}}}^s·{\underset{¯}{\varphi }}_n^s+\sum _j\text{}\underset{(N_n^s\times N_{f_j}^s)(N_{f_j}^s\times N_{f_j}^{m_j})(N_{f_j}^{m_j}\times 1)}{[{\underset{¯}{\underset{¯}{A}}}_{n{f}_j}^s·{\underset{¯}{\underset{¯}{K}}}_{\partial V_j}^{s\leftarrow m_j}·{\underset{¯}{\varphi }}_{f_j}^{m_j}]}\mathrm{}\mathrm{}=({\underset{¯}{b}}_n^s)\text{;}$$(80a) (80b) $$\sum _{\widehat{ı}}\text{}\underset{(N_{f_{\widehat{ı}}}^{m_{\widehat{ı}}}\mathrm{}\times N_{f_{\widehat{ı}}}^s)\mathrm{ }(N_{f_{\widehat{ı}}}^s\mathrm{}\times N_n^s)(N_n^s\mathrm{}\times 1)}{[{\left[{\underset{¯}{\underset{¯}{K}}}_{\partial V_{\widehat{ı}}}^{s\leftarrow m_{\widehat{ı}}}\right]}^T\mathrm{}·{\underset{¯}{\underset{¯}{A}}}_{f_{\widehat{ı}}n}^s·{\underset{¯}{\varphi }}_n^s}\text{}+\sum _{\widehat{\mathit{ȷ}}}\text{}\underset{(N_{f_{\widehat{ı}}}^{m_{\widehat{ı}}}\mathrm{}\times N_{f_{\widehat{ı}}}^s)\mathrm{ }(N_{f_{\widehat{ı}}}^s\mathrm{}\times N_{f_{\widehat{\mathit{ȷ}}}}^s)\mathrm{ }(N_{f_{\widehat{\mathit{ȷ}}}}^s\mathrm{}\times N_{f_{\widehat{\mathit{ȷ}}}}^{m_{\widehat{\mathit{ȷ}}}})\mathrm{ }(N_{f_{\widehat{\mathit{ȷ}}}}^{m_{\widehat{\mathit{ȷ}}}}\mathrm{}\times 1)}{[{\left[{\underset{¯}{\underset{¯}{K}}}_{\partial V_{\widehat{ı}}}^{s\leftarrow m_{\widehat{ı}}}\right]}^T\mathrm{}\mathrm{}·{\underset{¯}{\underset{¯}{A}}}_{f_{\widehat{ı}}{f}_{\widehat{\mathit{ȷ}}}}^s\mathrm{}·{\underset{¯}{\underset{¯}{K}}}_{\partial V_{\widehat{\mathit{ȷ}}}}^{s\leftarrow m_{\widehat{\mathit{ȷ}}}}\mathrm{}\mathrm{}·{\underset{¯}{\varphi }}_{f_{\widehat{\mathit{ȷ}}}}^{m_{\widehat{\mathit{ȷ}}}}]}\text{}=\underset{(N_{f_{\widehat{ı}}}^{m_{\widehat{ı}}}\mathrm{}\times N_{f_{\widehat{ı}}}^s)\mathrm{ }(N_{f_{\widehat{ı}}}^s\mathrm{}\times 1)}{{\left[{\underset{¯}{\underset{¯}{K}}}_{\partial V_{\widehat{ı}}}^{s\leftarrow m_{\widehat{ı}}}\right]}^T\mathrm{}·{\underset{¯}{b}}_{f_{\widehat{ı}}}^s].}$$(80b)

These equations can be assembled with those established on the master-patch adjoining the slave-patch at a given irregular interface, which, thus far, have remained unchanged, i.e., $\forall \forall \partial V_f\in V_m\cap \partial V_{\mathit{\text{sm}}}:$ (81a) $${\underset{¯}{\underset{¯}{A}}}_{\mathit{\text{nn}}}^m·{\underset{¯}{\varphi }}_n^m+{\underset{¯}{\underset{¯}{A}}}_{\mathit{\text{nf}}}^m·{\underset{¯}{\varphi }}_f^m={\underset{¯}{b}}_n^m\text{;}$$(81a) (81b) $${\underset{¯}{\underset{¯}{A}}}_{\mathit{\text{fn}}}^m·{\underset{¯}{\varphi }}_n^m+{\underset{¯}{\underset{¯}{A}}}_{\mathit{\text{ff}}}^m·{\underset{¯}{\varphi }}_f^m={\underset{¯}{b}}_f^m.$$(81b)

Referring back to the simple problem in Figure 5, $C^0$-continuity has been enforced by constraining those control-variables at the slave-patch boundary, in the bottom left-hand corner, on its west-side, ${\varphi }_{f_W}^s,$ and its north side, ${\varphi }_{f_N}^s,$ to those master-patches that adjoin the slave-patch on its west-side, $V_{m_W},$ and its north-side, $V_{m_N},$ respectively; i.e., $\prod {}_{h,p}^{m_W}{\varphi }_{f_W}^s:\{{\varphi }_5,{\varphi }_6,{\varphi }_9\}\to \{{\varphi }_5,{\varphi }_9\}$ and $\prod {}_{h,p}^{m_N}{\varphi }_{f_N}^s:\{{\varphi }_7,{\varphi }_8,{\varphi }_9\}\to \{{\varphi }_7,{\varphi }_9\},$ respectively. Again, those control-variables at regular and compatible boundaries exist in a one-to-one correspondence between all adjoining patches. As such, corner control-variables can be shared by several patches; e.g. ${\varphi }_9\in V_s\cap V_{m_W}\cap V_{m_N}.$ Consequently, a shared corner control-variable is constrained to itself across an irregular patch interface since it lies where the basis is interpolatory; e.g. ${\varphi }_9={\underset{¯}{\underset{¯}{K}}}^{s\leftarrow m_w}·{\varphi }_9.$ To ensure that such a corner control-variable does not become over-constrained, one must tally the number of times that it is constrained in order to correct the corresponding constraint equations; e.g.: $${\varphi }_9=\frac{1}{2}[{\underset{¯}{\underset{¯}{K}}}_{\partial V}^{s\leftarrow m_E}·{\varphi }_9]+\frac{1}{2}[{\underset{¯}{\underset{¯}{K}}}_{\partial V}^{s\leftarrow m_N}·{\varphi }_9]=\frac{1}{2}[{\varphi }_9+{\varphi }_9]={\varphi }_9.$$

A profound consequence of the proposed local-refinement technique is the proliferation of hanging-nodes at the irregular interfaces. More specifically, there may be instances where the constraints on the original set of (slave) control-variables are compounded, or multiply constrained (Demkowicz 2007). Multiply constrained control-variables complicate the numerical analysis because they create dependencies, which must be resolved to determine the order in which to apply the constraints and then recover the (slave) control-variables. To avoid such complications, one enforces a degree of mesh regularity.

Where corner control-variables can be shared by several patches, compounded constraints arise when the master to some slave-patch and that slave-patch are both constrained to some other common master-patch. For example, the left-hand side of Figure 6 presents such a dependency for the mesh of a simple problem. In this case, two patches in the lower quadrants have been h-refined; and consequently, ${\varphi }_8$ is multiply constrained; i.e., $\{{\varphi }_7,{\varphi }_8,{\varphi }_{13},{\varphi }_{14}\}\to \{{\varphi }_7,{\varphi }_{14}\},$ yet, $\{{\varphi }_4,{\varphi }_{10},{\varphi }_8\}\to \{{\varphi }_4,{\varphi }_8\}.$ Despite its constraints being compounded, ${\varphi }_{14},$ poses no problem because it is always constrained to itself; i.e., $\{{\varphi }_{14},{\varphi }_{19},{\varphi }_{20}\}\to \{{\varphi }_{14},{\varphi }_{20}\},$ yet, $\{{\varphi }_6,{\varphi }_{12},{\varphi }_{14}\}\to \{{\varphi }_6,{\varphi }_{14}\}.$ To alleviate the problem of multiply constraints, some regularity is enforced over the mesh. For example, on the right-hand side of Figure 6, an additional patch has been h-refined, as highlighted in blue, and now $\{{\varphi }_{11},{\varphi }_{12},{\varphi }_{13},{\varphi }_{17},{\varphi }_{18}\}\to \{{\varphi }_{11},{\varphi }_{18}\}.$

Figure 6. An example of forced regularity for a CBG-IGA discretization over an irregular mesh for the example of a simple problem; v. Figure 5. The initial mesh (left) presents a multiply constrained control-variable, ${\varphi }_8.$ Dependencies are resolved (right) by enforcing some degree of regularity; viz. the h-refinement highlighted in blue. For ease of viewing, only the enumeration of those control-variables that lie at an irregular patch interface are indicated. (V. the web-based version for reference to color.)

With regard to the interpolation-based nature of the proposed AMR strategy, an inherent bonus of defining the more-refined reference mesh, ${\mathcal{T}}_{\mathit{\text{ref}}}={\mathcal{T}}_{\text{h/2},p\text{+1}},$ as that resulting from one uniform hp-refinement of the coarser current mesh is that some regularity is implicitly enforced over the mesh where the numerical analysis is performed. Therefore, h-refinement resolves some of the dependencies and the elements are assured to have at least two compatible neighbors in 2D. Moreover, p-refinement ensures a basis of at least quadratic NURBS basis functions and that at least one central control-variable remains unconstrained since the support of interior functions become negligible at the boundary.

Figures 5 and 6 are useful for demonstrative purposes; however, without loss of generality, the proposed method for enforcing $C^0$-continuity via projection-based constraint equations is applicable to arbitrarily irregular meshes. Admittedly, the task of establishing the constraints and resolving any dependencies does become one of good bookkeeping.

6. The numerical results

This section presents some numerical results for the AMR-h(p) strategies developed in Section 2.1 and Section 3. Unless otherwise stated, all operations are performed to optimize the discretization error in the $H^1$-(semi)norm. The AMR strategies are compared against uniform refinement for a series of verification test cases using the method of manufactured solutions (MMS); and a couple of reactor physics benchmarks. Any data-plots that make reference to the DoFs of the Discretized system justly account for all unknowns that were solved for in the complete system.

The CBG-IGA spatial discretization schemes developed here have been implemented in a new modern Fortran code (Metcalf, Reid, and Cohen 2011). Adopting a distributed memory computing architecture, one makes use of the Open MPI message passing library; as well as the PETSc library and its capacity for parallelization of its sparse linear solvers and its compressed row storage (CRS) data structures for large sparse matrices (Balay et al. 2018). All jobs are submitted to the general purpose Imperial College London (ICL) HPC cluster computing service (CX1). The within-group matrix equations are converged to machine precision using PETSc’s conjugate gradient smoother, pre-conditioned by the BoomerAMG algebraic multi-grid method from the HYPRE library (Henson and Yang 2002). The scalar neutron flux is converged to a tolerance of $1.00\mathrm{E}-06$ between successive source iterations; and the eigenvalue is converged to a tolerance of $1.00\mathrm{E}-08$ between successive power iterations. An arbitrary upper-limit of ten is imposed upon the polynomial degree of the NURBS basis over the coarser current mesh; and therefore p_max = 11 over the reference mesh.

The number of outer-iterations required to attain convergence is reduced by bootstrapping the reference solution from the previous AMR-iteration, using it as an initial guess over the reference mesh in the current i^th-iteration (Wang 2009).

6.1. The method of manufactured solutions (MMS) verification

The method of manufactured solutions (MMS) is used to investigate the spatial discretizations of the multi-group NDE; cf. an overview by Salari et al. (Salari and Knupp 2000). A smooth manufactured solution is chosen to test against theoretical orders of accuracy (Rivière 2008; Warsa 2008): (82a) $$|e{|}_{H^1}\le C_1{h}^{\mathrm{\text{min}}\{p,r\}}\text{;}$$(82a) (82b) $${\Vert e\Vert }_{L_2}\le C_2{h}^{\mathrm{\text{min}}\{p+1,r\}}\text{;}$$(82b) where $e=\varphi -\widehat{\varphi }$ is the discretization error; h is the characteristic mesh-spacing, p is the polynomial degree of the underlying NURBS basis; r denotes the regularity index of the true solution; and $C_1$ and $C_2$ are constants independent of h and p. As per the tensor-product structure of the NURBS basis, Equation (3)(3) $$R_{i,j}^{p,q}(\xi ,\eta )=\frac{N_i^p(\xi )M_j^q(\eta )W_{i,j}}{\sum _{i\prime =1}^N\sum _{j\prime =1}^M{N}_{i\prime }^p(\xi )M_{j\prime }^q(\eta )W_{i\prime ,j\prime }}.$$(3) , the polynomial degree in the inequalities above, Equation (82), is actually $p={\mathrm{\text{min}}}_i\left\{p_i\right\}$ for $i=\xi ,\eta .$

For some manufactured solution, ${\varphi }^g,$ the (semi-)analytical manufactured fixed-source, $q_x^g,$ can be evaluated for $g=1,\dots ,G:$ (83) $$q_x^g=(L^g-S^g)\varphi .$$(83)

First, the manufactured solution and source are used to study the convergence rates of the uniform refinement of a Cartesian mesh and a more challenging pin-cell mesh. Second, one studies the NRG-AMR-h(p) strategy to verify all mesh-to-mesh interpolation procedures and the capacity to enforce continuity over irregular meshes that contain hanging-nodes.

6.1.1. MMS verification: a Cartesian mesh

The first MMS verification case consists of a Cartesian mesh defined over $V={[0,1]}^2.$ It comprises a homogeneous material whose properties are prescribed in Table 1. The coarsest mesh is composed of four equal-volume squares represented by bilinear NURBS patches. Assuming a one-speed, or one-group (1G), energy approximation, the manufactured solution is chosen as a standing wave equation: $$\{\begin{array}{ll}{\varphi }^n(\overrightarrow{r})=\hspace{0.17em}\mathrm{\text{cos}}\hspace{0.17em}(B_x^{n}x)\hspace{0.17em}\mathrm{\text{cos}}\hspace{0.17em}(B_y^{n}y),\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\overrightarrow{r}\in V\text{;}& \mathrm{(}\text{84a}\mathrm{)}\hfill \\ B_x^n,B_y^n:=(2n+1)\pi ,\hspace{1em}\forall n\ge 0\text{;}& \mathrm{(}\text{84b}\mathrm{)}\hfill \end{array}$$ where ${\varphi }^n$ denotes an odd harmonic mode as characterized by the geometric buckling, $B_x^n$ and $B_y^n.$ To allow for the maximum attainable convergence rates, the manufactured solution has been chosen to be smooth; yet, it cannot be represented exactly by arbitrarily high-order NURBS bases. Moreover, the homogeneous Neumann (reflective) BC is naturally satisfied.

Table 1. The nuclear data for the 2D 1G MMS verification.

Material	g	D^g(cm)	${\Sigma }_t^g$ (cm${}^{-1}$)	${\Sigma }_{s0}^{g\to 1}$ (cm${}^{-1}$)
(1) -	1	$3.333\mathrm{E}-01$	$1.000\mathrm{E}+00$	$5.000\mathrm{E}-01$

6.1.1.1. Uniform refinement

For the KI-refinement of linear, quadratic and cubic B-splines, the $H^1$-errors for the CBG-IGA spatial discretization are presented in Figure 7 for the manufactured solution chosen as the first, third and fifth harmonics of the standing wave equation; q.v. Equation (84). Qualitatively similar plots were obtained for the $L_2$-error; however, for brevity, these have not been presented. Under refinement, the continuous discretization schemes converge to the analytical solution; and asymptotically, they attain their theoretical orders of accuracy, as presented in Table 2.

Figure 7. The convergence plots, ${ϵ}_{H^1}$v h (left) and ${ϵ}_{H^1}$v DoFs (right), for the MMS verification of the uniform refinement of the CBG-IGA 1G NDE over a 2D Cartesian mesh. The manufactured solution is chosen as per Equation (84) for n = 0 (top row), n = 1 (middle row) and n = 2 (bottom row). (V. the web-based version for reference to color.)

Table 2. The asymptotic orders of accuracy, Equation (82), for the MMS verification of the uniform refinement of the CBG-IGA 1G NDE over a 2D Cartesian mesh.

		KI-1		KI-h
	p	${ϵ}_{H^1}$	${ϵ}_{L_2}$	${ϵ}_{H^1}$	${ϵ}_{L_2}$
(n = 0)	1	1.00	2.00	1.00	2.00
	2	2.00	3.00	2.00	3.00
	3	2.98	3.98	3.00	4.00
(n = 1)	1	0.99	2.00	0.99	2.00
	2	2.01	3.03	2.00	3.00
	3	3.01	4.03	3.00	4.00
(n = 2)	1	0.99	2.00	1.00	2.00
	2	2.04	3.08	2.00	2.99
	3	3.06	4.01	3.00	3.99

The convergence rates are evaluated as the slope of the last few points from each dataset.

Whilst the mesh spacing, h, does give some insight into how the discretization error varies under refinement, it does not reveal the true computational effort required; viz. the number of elements in the mesh does not vary under DE-refinement. Therefore, one should also consider the DoFs of the discrete system; this is a better metric by which to compare non-uniform refinement strategies.

For such a homogeneous problem with a smooth underlying solution, exploiting higher-order continuity in the basis yields more accuracy per DoF than that compared to a similar h-refinement. Moreover, one remarks that h-refinement becomes increasingly more expensive for increasingly large polynomial degrees.

An expected observation is that of the prolonged pre-asymptotic range of convergence, which can be seen on the left-hand side of Figure 7; it becomes increasingly large for increasingly-more oscillatory solutions. For Galerkin-based methods, unique solvability exists for certain conditions imposed upon the choice the discrete approximation space; viz. the minimum resolution for asymptotic convergence (Melenk, Parsania, and Sauter 2013).

6.1.1.2. Adaptive local refinement

The results of the NRG-AMR-h(p) strategy for chasing the global $H^1$-error in the CBG-IGA spatial discretization are presented in Figure 8 for the manufactured solution chosen as the first, third and fifth harmonics of the standing wave equation; q.v. Equation (84). These results are compared against those obtained for uniform h-refinement. For a coarsest mesh of linear B-splines, the NRG-AMR-hp strategy, denoted in the legend by P1-NRG-hp, exhibits exponential rates of convergence. Such a homogeneous and smooth problem favors p-refinement until p_max of the basis is attained, after which, only h-refinement is available. For coarsest meshes of linear and quadratic B-splines, the NRG-AMR-h strategy, denoted in the legend by P1-NRG-h and P2-NRG-h, respectively, yields algebraic rates of convergence, which are comparable with those of uniform refinement. One recalls that the computation is actually performed on the reference mesh, which is more refined than the current coarser mesh by one-level of uniform hp-refinement.

Figure 8. The convergence plots, ${ϵ}_{H^1}$v DoFs, for the MMS verification of the NRG-AMR-h(p) of the CBG-IGA 1G NDE over a 2D Cartesian mesh. The manufactured solution is chosen as per Equation (84) for n = 0 (top left), n = 1 (top right) and n = 2 (bottom). (V. the web-based version for reference to color.)

Inevitably, constraining a more-refined subspace to a less-refined one does introduce some numerical error at the irregular interface. Accordingly, it may appear that the local-refinement strategy does not offer much of a competitive edge for such a simple problem. However, these results do verify consistency in all mechanisms of the AMR strategy.

6.1.2. MMS verification: a pin-cell mesh

The second MMS verification case consists of a more challenging pin-cell mesh defined over $V={[0,1]}^2.$ It comprises a homogeneous material whose properties are prescribed in Table 3. The coarsest mesh is composed of nine equal-volume non-degenerate bi-quadratic NURBS patches. Assuming a two-group (2G) energy approximation, the manufactured solution is chosen as: $$\{\begin{array}{l}{\varphi }^g(\overrightarrow{r})=\hspace{0.17em}\mathrm{\text{sin}}\hspace{0.17em}(B_x^{g}x)\hspace{0.17em}\mathrm{\text{sin}}\hspace{0.17em}(B_y^{g}y),\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\overrightarrow{r}\in V\text{;\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.5em}}\mathrm{(}\text{85a}\mathrm{)}\\ B_x^g,B_y^g:=(2g+1)\pi ,\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\text{for \hspace{0.5em}}g=1,2\text{;\hspace{1em}\hspace{1em}\hspace{0.5em}}\mathrm{(}\text{85b}\mathrm{)}\end{array}$$ where ${\varphi }^g$ denotes an odd harmonic mode as characterized by the geometric buckling, $B_x^g$ and $B_y^g.$ Again, the manufactured solution has been chosen to be infinitely differentiable. However, now, the homogeneous Dirichlet (zero-flux) BC is naturally satisfied.

Table 3. The nuclear data for the 2D 2G MMS verification.

Material	g	D^g (cm)	${\Sigma }_t^g$ (cm${}^{-1}$)	${\Sigma }_{s0}^{g\to 1}$ (cm${}^{-1}$)	${\Sigma }_s^{g\to 2}$ (cm${}^{-1}$)
(1) -	1	$3.333\mathrm{E}-01$	$1.000\mathrm{E}+00$	$5.000\mathrm{E}-01$	$2.500\mathrm{E}-01$
(1) -	2	$3.333\mathrm{E}-01$	$1.000\mathrm{E}+00$	$0.000\mathrm{E}+00$	$5.000\mathrm{E}-01$

6.1.2.1. Uniform refinement

For the KI-refinement of quadratic, cubic and quartic NURBS bases, the $H^1$-errors for the CBG-IGA spatial discretization are presented in Figure 9 for the multi-group components of the manufactured solution chosen as the third and fifth harmonics of the standing wave equation; q.v. Equation (85). Qualitatively similar plots were obtained for the $L_2$-error; however, for brevity, these have not been presented. Under refinement, the continuous discretization schemes converge to the analytical solution; and asymptotically, they attain their theoretical orders of accuracy for KI-h-refinement, as presented in Table 4. However, for KI-1-refinement, one remarks a degree of super-convergence, which may be explained either by the higher-order continuity in the basis or the weak imposition of the Dirichlet BCs. Although not presented in this paper, for the same pin-cell mesh with naturally enforced reflective BCs and manufactured solutions chosen as per Equation (84), one did observe the expected theoretical rates of convergence.

Figure 9. The convergence plots, ${ϵ}_{H^1}$v h (left) and ${ϵ}_{H^1}$v DoFs (right), for the MMS verification of the uniform refinement of the CBG-IGA 2G NDE over a 2D pin-cell mesh. The manufactured solutions are chosen as per Equation (85) for g = 1, 2. (V. the web-based version for reference to color.)

Table 4. The asymptotic orders of accuracy, Equation (82), for the MMS verification of the uniform refinement of the CBG-IGA 2G NDE over a 2D pin-cell mesh.

p	KI-1		KI-h
	${ϵ}_{H^1}$	${ϵ}_{L_2}$	${ϵ}_{H^1}$	${ϵ}_{L_2}$
2	2.06	3.15	1.99	2.98
3	3.09	4.18	3.00	4.00
4	4.29	5.40	3.99	4.99

The convergence rates are evaluated as the slope of the last few points from each dataset.

6.1.2.2. Adaptive local refinement

The results of the NRG-AMR-h(p) strategy for chasing the global $H^1$-error in the CBG-IGA discretization are presented in Figure 10 for the multi-group components of the manufactured solution chosen as the third and fifth harmonics of the standing wave equation; q.v. Equation (85)(60) $$\sum _{\partial V_f\in \partial V_D}{\mu }_f{\int }_{\partial V_f}w(\overrightarrow{r})\varphi (\overrightarrow{r})\mathit{\text{ dS}}=\sum _{\partial V_f\in \partial V_D}{\mu }_f{\int }_{\partial V_f}w(\overrightarrow{r})g_D\mathit{\text{ dS}}.$$(60) . These results are compared against those obtained for uniform h-refinement. For a coarsest mesh of quadratic NURBS basis functions, the NRG-AMR-hp strategy, denoted in the legend by P2-NRG-hp, again, exhibits faster rates of convergence; although this time the mesh iteration is not dominated by p-refinement. For coarsest meshes of quadratic and cubic NURBS basis functions, the NRG-AMR-h strategy, denoted in the legend by P2-NRG-h and P3-NRG-h, respectively, again, yields algebraic rates of convergence, which are comparable with those of uniform refinement.

Figure 10. The convergence plots, ${ϵ}_{H^1}$v DoFs, for the MMS verification of the NRG-AMR-h(p) of the CBG-IGA 2G NDE over a 2D pin-cell mesh. The prefixes M1- and M2- denote the use of energy-independent and energy-dependent meshes, respectively. The manufactured solutions are chosen as per Equation (85) for g = 1, 2. (V. the web-based version for reference to color.)

Within a multi-group formalism, it is fitting to test the capacity for energy-dependent meshes, denoted in the legend by the prefix M2-, which yield more accuracy per DoF when compared to the refinement strategies that employ energy-independent meshes, denoted in the legend by the prefix M1-. The effectivity, ${\theta }_{H^1},$ of the interpolation-based error measures, as per Equation (16)(16) $$\theta =\frac{\widehat{ϵ}}{ϵ}.$$(16) , is measured and displayed in Figure 11. Since ${\theta }_{H^1}$ tends to values slightly larger than unity, one can approximate the remaining error in the discretization, which may be used to terminate the mesh iteration procedure.

Figure 11. The plots of the effectivity index, ${\theta }_{H^1}$v DoFs, for the MMS verification of the NRG-AMR-h(p) of the CBG-IGA 2G NDE over a 2D pin-cell mesh. The prefixes M1- and M2- denote the use of energy-independent and energy-dependent meshes, respectively. The manufactured solutions are chosen as per Equation (85) for g = 1, 2. (V. the web-based version for reference to color.)

Again, these results do verify consistency in all mechanisms of the AMR strategy over irregular meshes for rational geometries.

6.2. The 2D IAEA/ANL BSS-11 (LWR) benchmark

The 2D IAEA benchmark for the 2G NDE is the BSS-11 quarter-core LWR problem taken from the Argonne National Laboratory (ANL) Benchmark Problem Book, Supplement 2 (Argonne National Laboratory 1977). Taken at the mid-plane, $z=1.90\mathrm{E}+02$ cm, this 2D benchmark is adapted from the full 3D (xyz) one; an axial buckling, $B_z^2=8.00\mathrm{E}-05,$ approximates the effects of axial leakage (Theler 2013). Accordingly, the corrected macroscopic nuclear cross-section data are given in Table 5. Materials 1, 2 and 3 represent fueled regions, each with slightly different properties to mimic varying levels of enrichment or the presence of control-rods; whereas Material 4 represents a reflector region; v. Figure 12.

Figure 12. The geometry of the 2D 2G IAEA/ANL BSS-11 benchmark (Argonne National Laboratory 1977). The constituent materials are denoted by the numbers in parentheses; viz. fuel-1 (1)(1) $$\text{supp}(R)=\{\overrightarrow{r}\in V_e:R(\overrightarrow{r})\ne 0,\forall {\overrightarrow{r}}^{\text{−}1}\in {\widehat{V}}_e\},\hspace{1em}\forall V_e\in \mathcal{T}(V).$$(1) , fuel-2 (2)(2) $$\text{supp}(N_i^p)=\{\xi \in \widehat{V}\subset {\mathbb{R}}^1:N_i^p(\xi )\ne 0,\forall {\xi }_i\le \xi <{\xi }_{i+p+1}\}\text{;}$$(2) , fuel-2 & control-rod (3)(3) $$R_{i,j}^{p,q}(\xi ,\eta )=\frac{N_i^p(\xi )M_j^q(\eta )W_{i,j}}{\sum _{i\prime =1}^N\sum _{j\prime =1}^M{N}_{i\prime }^p(\xi )M_{j\prime }^q(\eta )W_{i\prime ,j\prime }}.$$(3) and reflector (4)(4) $$S(\xi ,\eta )=\sum _{i=1}^N\sum _{j=1}^M{\overrightarrow{B}}_{i,j}{R}_{i,j}^{p,q}(\xi ,\eta ),\hspace{1em}\forall \xi ,\eta \in \widehat{V}.$$(4) . (V. the web-based version for reference to color.)

Table 5. The nuclear data for the 2D 2G IAEA/ANL BSS-11 benchmark (Argonne National Laboratory 1977).

Material	g	D^g (cm)	${\Sigma }_r^{g*}$ (cm${}^{-1}$)	${\Sigma }_{s0}^{g\to 2}$ (cm${}^{-1}$)	${\nu }^g{\Sigma }_f^g$ (cm${}^{-1}$)
(1) Fuel-1	1	$1.50000\mathrm{E}+00$	$3.01200\mathrm{E}-02$	$2.00000\mathrm{E}-02$	$0.00000\mathrm{E}+00$
	2	$4.00000\mathrm{E}-01$	$8.00320\mathrm{E}-02$	$0.00000\mathrm{E}+00$	$1.35000\mathrm{E}-01$
(2) Fuel-2	1	$1.50000\mathrm{E}+00$	$3.01200\mathrm{E}-02$	$2.00000\mathrm{E}-02$	$0.00000\mathrm{E}+00$
	2	$4.00000\mathrm{E}-01$	$8.50320\mathrm{E}-02$	$0.00000\mathrm{E}+00$	$1.35000\mathrm{E}-01$
(3) Fuel-3 & Rod	1	$1.50000\mathrm{E}+00$	$3.01200\mathrm{E}-02$	$2.00000\mathrm{E}-02$	$0.00000\mathrm{E}+00$
	2	$4.00000\mathrm{E}-01$	$1.30032\mathrm{E}-01$	$0.00000\mathrm{E}+00$	$1.35000\mathrm{E}-01$
(4) Reflector	1	$2.00000\mathrm{E}+00$	$4.01600\mathrm{E}-02$	$4.00000\mathrm{E}-02$	$0.00000\mathrm{E}+00$
	2	$3.00000\mathrm{E}-01$	$1.00240\mathrm{E}-02$	$0.00000\mathrm{E}+00$	$0.00000\mathrm{E}+00$

To better approximate the effects of axial leakage, one applies an axial buckling, $B_z^2=8.00000\mathrm{E}-05.$ Accordingly, the effective macroscopic cross-sections of removal are corrected as ${\Sigma }_r^{g*}={\Sigma }_r^g+D^g{B}_z^2.$

One assumes no incoming current at the external boundaries, as modeled by (extrapolated) vacuum BCs; and as per the one-quarter symmetry in the problem, reflective BCs are prescribed at $x=y=0.00$ cm. The coarsest Cartesian mesh is constructed from 241 bilinear NURBS patches. A reference solution was obtained over the energy-independent mesh generated after the ${11}^{\mathit{\text{th}}}$ iteration of the DWR-AMR-hp strategy; viz. ${\mathrm{k}}_{\text{eff}}=1.0295886369,$ which is in good agreement with the literature (Argonne National Laboratory 1977). The spatial distribution of the multi-group components of both the scalar neutron flux and the scalar neutron importance are presented in Figure 13, for the QoI chosen as the ${\mathrm{k}}_{\text{eff}}.$

Figure 13. The spatial distributions of the multi-group components of the scalar neutron flux (left) and the scalar neutron importance (right) of the 2G NDE for the 2D IAEA/ANL BSS-11 benchmark. The QoI is the ${\mathrm{k}}_{\text{eff}}.$ (V. the web-based version for reference to color.)

Those meshes generated after the $9^{\mathit{\text{th}}}$ AMR-iteration of the AMR-h and AMR-hp strategies are presented in Figures 14 and 15, respectively.

Figure 14. The energy-independent meshes (top) and the energy-dependent meshes (middle and bottom) generated at the $9^{\mathit{\text{th}}}$ AMR-iteration of the NRG-AMR-h (left) and the DWR-AMR-h (right) of the CBG-IGA 2G NDE for the 2D IAEA/ANL BSS-11 benchmark. The QoI is the ${\mathrm{k}}_{\text{eff}}.$ (V. the web-based version for reference to color.)

Figure 15. The energy-independent meshes (top) and the energy-dependent meshes (middle and bottom) generated at the $9^{\mathit{\text{th}}}$ AMR-iteration of the NRG-AMR-hp (left) and the DWR-AMR-hp (right) of the CBG-IGA 2G NDE for the 2D IAEA/ANL BSS-11 benchmark. The QoI is the ${\mathrm{k}}_{\text{eff}}.$ (V. the web-based version for reference to color.)

At the cost of increasingly ill-conditioned matrices, the hp-refinement strategies performed better than the corresponding h-refinement ones; v. Figure 16. The computation times for the former were not prohibitive since the chosen algebraic multi-grid pre-conditioner worked effectively. Whilst imposing a certain degree of regularity to avoid any multiply constraints, the discretization developed in Section 5 does not appear to introduce any significant sources of error across the irregular interfaces of this heterogeneous problem.

Figure 16. The convergence plots, ϵ_QoI v DoFs, for the uniform refinement (black), the NRG-AMR (red) and the DWR-AMR (blue) of the CBG-IGA 2G NDE for the 2D IAEA/ANL BSS-11 benchmark. The prefixes M1- and M2- denote the use of energy-independent and energy-dependent meshes, respectively. The QoI is the ${\mathrm{k}}_{\text{eff}}.$ (V. the web-based version for reference to color.)

Both AMR strategies performed better than the uniform one. This is unsurprising as the ${\mathrm{k}}_{\text{eff}}$ is a global quantity and both AMR strategies chase the global error in the primal solution. However, at the cost of solving the dual problem, the DWR-AMR strategies do allow for a better allocation of resources; and better still when employing energy-dependent meshes.

6.3. The 2D BIBLIS (PWR) benchmark

The checkerboard fuel-loading pattern of the 2D BIBLIS benchmark for the 2G NDE better represents the core of an operating PWR. First proposed by Finemann and Wagner, Herbert provided an updated version to correspond to a rods-out configuration (Hébert 1985, 2019). This problem presents eight physically distinct regions, including a homogenized baffle-and-reflector region and seven smeared fueled regions of varying degrees of burn up; each element is a square of congruent sides of length $2.31226\mathrm{E}+01$ cm; v. Figure 17. The corresponding macroscopic nuclear cross-section data are given in Table 6; no axial buckling information is provided.

Figure 17. The geometry of the 2D 2G BIBLIS benchmark (Hébert 1985). The constituent materials are denoted by the numbers in parentheses; viz. fuel-1 (1)(1) $$\text{supp}(R)=\{\overrightarrow{r}\in V_e:R(\overrightarrow{r})\ne 0,\forall {\overrightarrow{r}}^{\text{−}1}\in {\widehat{V}}_e\},\hspace{1em}\forall V_e\in \mathcal{T}(V).$$(1) , fuel-2 (2)(2) $$\text{supp}(N_i^p)=\{\xi \in \widehat{V}\subset {\mathbb{R}}^1:N_i^p(\xi )\ne 0,\forall {\xi }_i\le \xi <{\xi }_{i+p+1}\}\text{;}$$(2) , baffle & reflector (3)(3) $$R_{i,j}^{p,q}(\xi ,\eta )=\frac{N_i^p(\xi )M_j^q(\eta )W_{i,j}}{\sum _{i\prime =1}^N\sum _{j\prime =1}^M{N}_{i\prime }^p(\xi )M_{j\prime }^q(\eta )W_{i\prime ,j\prime }}.$$(3) , fuel-3 (4)(4) $$S(\xi ,\eta )=\sum _{i=1}^N\sum _{j=1}^M{\overrightarrow{B}}_{i,j}{R}_{i,j}^{p,q}(\xi ,\eta ),\hspace{1em}\forall \xi ,\eta \in \widehat{V}.$$(4) , fuel-4 (5)(5) $${\underset{¯}{\underset{¯}{\mathbf{J}}}}_V·\overrightarrow{\nabla }{R}_i^{p,q}(\xi ,\eta )=\widehat{\nabla }{R}_i^{p,q}(\xi ,\eta ),\hspace{1em}\text{for\hspace{0.5em}}i=1,\dots ,N\times M.$$(5) , fuel-5 (5)(5) $${\underset{¯}{\underset{¯}{\mathbf{J}}}}_V·\overrightarrow{\nabla }{R}_i^{p,q}(\xi ,\eta )=\widehat{\nabla }{R}_i^{p,q}(\xi ,\eta ),\hspace{1em}\text{for\hspace{0.5em}}i=1,\dots ,N\times M.$$(5) , fuel-6 (7)(7) $${\underset{¯}{\underset{¯}{\mathbf{J}}}}_V:=\left(\begin{array}{cc}\partial x/\partial \xi & \partial y/\partial \xi \\ \partial x/\partial \eta & \partial y/\partial \eta \end{array}\right).$$(7) and fuel-7 (8). (V. the web-based version for reference to color.)

Table 6. The nuclear data for the 2D 2G BIBLIS benchmark (Hébert 1985).

Material	g	D^g (cm)	${\Sigma }_r^g$ (cm${}^{-1}$)	${\Sigma }_{s0}^{g\to 2}$ (cm${}^{-1}$)	${\nu }^g{\Sigma }_f^g$ (cm${}^{-1}$)
(1) -	1	$1.43600\mathrm{E}+00$	$2.72582\mathrm{E}-02$	$1.77540\mathrm{E}-02$	$5.87080\mathrm{E}-03$
	2	$3.63500\mathrm{E}-01$	$7.50580\mathrm{E}-02$	$0.00000\mathrm{E}+00$	$9.60670\mathrm{E}-02$
(2) -	1	$1.43660\mathrm{E}+00$	$2.72995\mathrm{E}-02$	$1.76210\mathrm{E}-02$	$6.19080\mathrm{E}-03$
	2	$3.63600\mathrm{E}-01$	$7.84360\mathrm{E}-02$	$0.00000\mathrm{E}+00$	$1.03580\mathrm{E}-01$
(3) -	1	$1.32000\mathrm{E}+00$	$2.57622\mathrm{E}-02$	$2.31060\mathrm{E}-02$	$0.00000\mathrm{E}+00$
	2	$2.77200\mathrm{E}-01$	$7.15960\mathrm{E}-02$	$0.00000\mathrm{E}+00$	$0.00000\mathrm{E}+00$
(4) -	1	$1.43890\mathrm{E}+00$	$2.74640\mathrm{E}-02$	$1.71010\mathrm{E}-02$	$7.45270\mathrm{E}-03$
	2	$3.63800\mathrm{E}-01$	$9.14080\mathrm{E}-02$	$0.00000\mathrm{E}+00$	$1.32360\mathrm{E}-01$
(5) -	1	$1.43810\mathrm{E}+00$	$2.72930\mathrm{E}-02$	$1.72900\mathrm{E}-02$	$6.19080\mathrm{E}-03$
	2	$3.66500\mathrm{E}-01$	$8.48280\mathrm{E}-02$	$0.00000\mathrm{E}+00$	$1.03580\mathrm{E}-01$
(6) -	1	$1.43850\mathrm{E}+00$	$2.73240\mathrm{E}-02$	$1.71920\mathrm{E}-02$	$6.42850\mathrm{E}-03$
	2	$3.66500\mathrm{E}-01$	$8.73140\mathrm{E}-02$	$0.00000\mathrm{E}+00$	$1.09110\mathrm{E}-01$
(7) -	1	$1.43890\mathrm{E}+00$	$2.72900\mathrm{E}-02$	$1.71250\mathrm{E}-02$	$6.19080\mathrm{E}-03$
	2	$3.67900\mathrm{E}-01$	$8.80240\mathrm{E}-02$	$0.00000\mathrm{E}+00$	$1.03580\mathrm{E}-01$
(8) -	1	$1.43930\mathrm{E}+00$	$2.73210\mathrm{E}-02$	$1.70270\mathrm{E}-02$	$6.42850\mathrm{E}-03$
	2	$3.68000\mathrm{E}-01$	$9.05100\mathrm{E}-02$	$0.00000\mathrm{E}+00$	$1.09110\mathrm{E}-01$

The average number of neutrons released per fission is $\nu =2.47$ for the fission spectrum $\chi ={[1.00,0.00]}^T.$

One assumes no incoming current at the external boundaries, as modeled by vacuum BCs; and as per the one-quarter symmetry in the problem, reflective BCs are prescribed at $x=y=0.00$ cm. The coarsest Cartesian mesh is constructed from 257 bilinear NURBS patches. Again, a reference was obtained over the energy-independent mesh generated after the ${11}^{\mathit{\text{th}}}$ iteration of the DWR-AMR-hp strategy. The spatial distribution of the multi-group components of both the scalar neutron flux and the scalar neutron importance are presented in Figure 18. From the literature, ${\mathrm{k}}_{\text{eff}}=1.0250103$ (Hébert 2019). However, in this instance, the QoI is chosen as the absorption in Material 5; v. Figure 17.

Figure 18. The spatial distributions of the multi-group components of the scalar neutron flux (left) and the scalar neutron importance (right) of the 2G NDE for the 2D BIBLIS benchmark. The QoI is the absorption in Material 5. (V. the web-based version for reference to color.)

The meshes generated after the $9^{\mathit{\text{th}}}$ AMR-iteration of the AMR-h and AMR-hp strategies are presented in Figures 19 and 20, respectively. Indicative timings for the different stages of analysis at the $9^{\mathit{\text{th}}}$ AMR-iteration are presented in Table 7. In any case, the relative time spent in the Assembly stage is negligle. However, the relative time spent in the Refinement stage became dominating for the AMR-h strategies, which includes the calculating of the interpolation-based error measures; and establishing the projection-based constraint equations over the more-refined reference meshes. Whereas, the relative time spent in the Solution stage became dominating for the AMR-hp strategies, which includes evaluating the group-to-group sources over energy-dependent meshes.

Figure 19. The energy-independent meshes (top) and the energy-dependent meshes (middle and bottom) generated at the $9^{\mathit{\text{th}}}$ AMR-iteration of the NRG-AMR-h (left) and the DWR-AMR-h (right) of the CBG-IGA 2G NDE for the 2D BIBLIS benchmark. The QoI is the absorption in Material 5. (V. the web-based version for reference to color.)

Figure 20. The energy-independent meshes (top) and the energy-dependent meshes (middle and bottom) generated at the $9^{\mathit{\text{th}}}$ AMR-iteration of the NRG-AMR-hp (left) and the DWR-AMR-hp (right) of the CBG-IGA 2G NDE for the 2D BIBLIS benchmark. The QoI is the absorption in Material 5. (V. the web-based version for reference to color.)

Table 7. The indicative timings for the different stages of analysis at the $9^{\mathit{\text{th}}}$ AMR-iteration of the NRG-AMR-h(p) and DWR-AMR-h(p) of the CBG-IGA 2G NDE over energy-dependent meshes for the 2D BIBLIS benchmark.

Stage	NRG-AMR-h	NRG-AMR-hp	DWR-AMR-h	DWR-AMR-hp
Assembly	1%	1%	¡1%	1%
Solution	42%	68%	23%	81%
Refinement	57%	31%	77%	18%
Total (s)	138	8	106	10

The Solution stage includes: evaluating the group-to-group sources over energy-dependent meshes. The Refinement stage includes: calculating the a posteriori interpolation-based error measures; and establishing the projection-based constraint equations over the more-refined reference meshes.

Due to the loading pattern and the relatively large dimensions of the fuel assemblies, this problem is challenging to solve with low-order spatial discretizations such as a linear FEM. One would not expect such difficulties for the proposed AMR strategies since the coarsest reference mesh is constructed from bi-quadratic NURBS patches. However, the fast and thermal components of the multi-group fluxes present very different spatial distributions.

Clearly, one single mesh could not suitably allow for an optimum allocation of resources. Accordingly, over energy-dependent meshes, almost twice as many DoFs were allocated to ${\mathcal{T}}^2$ than to ${\mathcal{T}}^1$ for the DWR-AMR strategy; whereas this was more than four times as many for the NRG-AMR strategy. The latter aims to distribute evenly the error in the solution in a global sense and thus, better approximate those sharp gradients in the thermal flux profile almost everywhere. This can be observed as the h-refinement about the material discontinuities in Figure 19. However, this does not contribute to a consistent reduction in the error in QoI; v. Figure 21. Accommodating the need to minimize those dominating errors in the thermal scalar neutron flux distribution, the NRG-AMR strategy over an energy-independent mesh outperforms that over energy-dependent ones. In the case of the former, ${\mathcal{T}}^1$ benefits from a refinement about Material 5, where it otherwise does not. Conversely, the DWR-AMR strategies consistently yield more accuracy per DoF by refining each mesh about the region of interest.

Figure 21. The convergence plots, ϵ_QoI v DoFs, for the uniform refinement (black), the NRG-AMR (red) and the DWR-AMR (blue) of the CBG-IGA 2G NDE for the 2D BIBLIS benchmark. The prefixes M1- and M2- denote the use of energy-independent and energy-dependent meshes, respectively. The QoI is the absorption in Material 5. (V. the web-based version for reference to color.)

7. Conclusion

This paper proposes self-adaptive NRG-AMR-h(p) and DWR-AMR-h(p) strategies for a CBG-IGA spatial discretization of the multi-group NDE, which can be readily extended to full 3D heterogeneous nuclear reactor physics problems. Higher-order continuity within a NURBS patch was sacrificed to accomodate wholly-local AMR algorithms. Where not trivially so, continuity in the solution field across irregular patch interfaces was enforced via constraint equations based upon a master-slave relationship and the conservative interpolation between surface meshes. A similar Galerkin projection was implemented for the conservative interpolation between volume meshes when evaluating group-to-group sources over energy-dependent meshes; and when calculating interpolation-based error measures. The spatial discretization results in large, sparse, SPD matrices that can benefit from PETSc’s capacity for distributed-memory CRS data structures; and the parallel implementation of its PCG smoothers. Moreover, the algebraic multi-grid pre-conditioner was effective for those problems posed over Cartesian meshes.

Whilst enforcing continuity over mesh irregularities does introduce some discretization error, the capacity for local refinement offers practical gains in a better allocation of computational resources. For global quantities such as the ${\mathrm{k}}_{\text{eff}},$ the NRG-AMR strategies for were competitive with the DWR-AMR ones. However, for localized QoI, the former becomes less performant than uniform h-refinement; and especially so over energy-dependent meshes, since certain multi-group components were under-resolved. Conversely, the DWR-AMR strategies consistently and reliably reduced the error in the QoI at each subsequent AMR-iteration.

Indeed, only the isotropic components of the numerical solution were used to evaluate the error measures in the $H^1$-norm. Instead, an anisotropic AMR algorithm could be driven by separating out the direction derivatives; and each parametric direction of a flagged element could be refined independently. Moreover, one could perform a sensitivity analysis on the floating parameters prescribed in Equation (25)(25) $$\{\begin{array}{cc}{\beta }_1=\text{1/3}\hfill & \mathrm{(}\text{Equation\hspace{0.5em}20}\mathrm{)}\text{;}\hfill \\ {\beta }_2=p_{\mathrm{\text{max}}}\hfill & \mathrm{(}\text{Equation\hspace{0.5em}23}\mathrm{)}\text{;}\hfill \\ {\beta }_3=\text{1/2}\hfill & (\text{Equation\hspace{0.5em}}24).\hfill \end{array}$$(25) , which are employed to constrain the mechanisms of local refinement.

Although not presented in this paper, the AMR-h(p) algorithms for the proposed CBG-IGA multi-group NDE did generate qualitatively similar meshes to those generated by a symmetric interior penalty scheme for a discontinuous Galerkin (SIP-DG-)IGA spatial discretization. Indeed, the latter incurs at least twice as many DoFs as the former for qualitatively similar meshes. Yet, the latter does offer a more-compact stencil and more flexibility in terms of local refinement. With regard to the mesh-to-mesh interpolation routines, the full transfer matrices could be formed locally and inverted as fully-discontinuous block matrices. Conversely, to avoid any loss of information for the globally-defined basis functions of the CBG-IGA discretization, it was necessary to form and inverted globally-assembled transfer matrices. That being said, it is difficult to compare directly the efficiency of the proposed CBG-IGA discretization and a SIP-DG-IGA one using only DoFs. Such a metric does not account for the computational effort attributed to the different stages of analysis; cf. Table 7.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data Statement

In accordance with EPSRC funding requirements, all supporting data used to create figures in this paper may be accessed via the following URL: https://doi.org/10.5281/zenodo.10120265.

Additional information

Funding

S. G. Wilson acknowledges the financial support of the Engineering and Physical Sciences Research Council (EPSRC) under their doctoral training partnership (DTP) scheme (EPSRC Grant No. EP/N509486/1). S. G. Wilson also acknowledges the financial support of Rolls-Royce for the funding of the CASE conversion of the EPSRC DTP PhD. Finally, the authors thank the high-performance computing (HPC) service support team at Imperial College London.

Appendix A.

AMR-procedure flowchart

Figure A1. The AMR-h(p) iteration procedure is enclosed in the red dotted line. The (black dashed) process-flow looping operations take priority. For NRG-AMR, the interpolation-based error-estimates are weighted by some appropriate energy-dependent function. For DWR-AMR, they are weighted by the error in the importance function; and therefore, for some prescribed QoI, Process P-1 must also resolve the adjoint (dual) problem. Which elements are eligible for further refinement is determined at Decision D-1; whereas constraints are placed on how they are refined at Decisions D-2 and D-3. (V. the web-based version for reference to color.)