An adaptive reduced order model for the angular discretization of the Boltzmann transport equation using independent basis sets over a partitioning of the space‐angle domain

This article presents a new reduced order model (ROM) for the angular discretization of the Boltzmann transport equation. The angular ROM is built over a partitioning of the space‐angle phase‐space, by generating independent, optimized angular basis function sets for each partition. The advantage is that each basis function set is optimized to represent the neutron flux distribution in a particular partition of space and angle, rather than being optimized for the entire domain. This serves to reduce the total number of basis functions required, and therefore the solve time. Two numerical examples are presented to demonstrate the efficacy of the methods—a dog‐leg duct problem, involving particle streaming, and the Watanabe–Maynard problem, which includes significant particle scattering. In both cases, it is shown that the method reduces the angular flux error, for a given basis size or solve time, by around an order of magnitude in comparison to other, similar ROM methods. An adaptive version of the method is also presented, whereby the number of basis functions within each space‐angle partition can vary independently. It is shown to potentially provide further significant reductions in error.

example, reactor cores are constructed from many tens of thousands of fuel and control rods, so large spatial meshes are required to resolve their full domain. Furthermore, predicting the distribution of neutrons within cores accurately involves solving the Boltzmann transport equation (BTE), which has a seven-dimensional phase-space.
Over the years, many articles have been published which focus on the formulation of efficient solutions to the BTE, such as the use of self-adaptive methods [4][5][6] which reduce the total degrees of freedom by distributing resolution efficiently-that is, only where it is beneficial. However, over the last decade reduced order models (ROMs) have gained popularity, as they are capable of forming and applying basis functions which are optimized for a particular class of problem. They have been demonstrated to reduce problem sizes by orders of magnitude, resulting in a similar scale of reductions to solve times.
ROMs have been applied extensively across many fields of engineering. Recent developments include the application of proper orthogonal decomposition (POD) to transient heat conduction, 7 turbulent supersonic jets, 8 acoustic waves, 9 incompressible magnetohydrodynamics, 10 and many other problems. POD-based methods have been combined with neural networks to model plasticity, 11 unsteady flows in a combustion problem, 12 the viscous Burgers equation, 13 and structural damage. 14 Purely neural network-based ROMs have also been implemented, with applications including the modeling of biomass fast pyrolysis 15 and pH reactors. 16 Many other model order reduction strategies have been developed, including the Volterra 17 and Fourier 18 series expansions, space mapping, 19 the Kriging 20 and harmonic balance 21 methods, and radial basis functions. 22 Related to the work of this article, ROMs have also been developed that use and apply basis functions over a domain decomposition, or localized regions of a problem. These have been used to solve advection-diffusion problems, 23 Maxwell's equations, 24 the Stokes equations, 25 and elliptic multi-scale problems. Mixed approaches have also been described using both high-order methods and ROMs over a domain decomposition, applying ROMs over regions that they can resolve well, with high order methods used elsewhere. Applications include solutions to the Laplace equations 26 and, more recently, the compressible Euler equations. 27 The field of nuclear engineering and general radiation transport simulation has made use of ROMs to solve a variety of problems. Early examples of ROMs in the field include Reference 28, which used several ROM techniques to analyze the transient dynamics of reactor-driven systems; Reference 29, which applied them to one dimensional transient radiation problems; and Reference 30, which applied POD to the spatial dimensions of eigenvalue problems in the context of reactor physics. Further developments include the application of POD to the angular dimension of the BTE, [31][32][33] space-angle ROMs for radiative heat transfer, 34 and the use of range-finding algorithms for the linear transformation of parameters in multiphysics problems. 35 More recently, POD has been used to model fuel burnup 36 and reactor power distributions. 37 Many articles have also considered the problem of control rod movement using various ROMs, including early work based on proper generalized decomposition, 38 and recent publications which employed POD, 39 the empirical interpolation method, 40 and neural networks. 41 Recently, a new adaptive POD method known as discontinuous POD (DPOD) was developed, which applied angular POD functions to optimally resolve the angular dimension of the BTE over distinct partitions of the angular domain. 42 The method was demonstrated to create stable basis sets which provided efficient solutions to multi-dimensional transport problems. In this article, a simple but highly effective new modification is presented to extend this capability, in order to improve the efficiency of the method by further reducing the dimensionality of problems. Known here as regional discontinuous POD (RDPOD), the new model forms and applies independent sets of angular basis functions, which are each optimized to resolve different partitions of the space-angle phase-space of the BTE. That is, the method forms separate angular POD functions over distinct regions of space, rather than using a single basis set for all regions of space as developed in the previous method. 42 While this is a simple modification, it is of particular importance as the neutron flux distribution often varies substantially over space and angle. Capturing the characteristics of the neutron flux using a single set of basis functions can therefore place high demands on the original ROM formulation. Partitioning the spatial and angular dimensions and creating separate angular ROMs for each partition can help to overcome this issue. As the variation in neutron flux distributions within each partition is considerably smaller than the variation across the entire problem, the number of basis functions required to resolve each distribution is reduced substantially. This can lead to smaller systems of equations, and therefore reduce solve times.
While there has been a large amount of research interest in developing ROMs which are localized over decomposed domains, this is the first to do so using angular basis functions to resolve the BTE. Furthermore, this article presents an implementation of angular adaptivity for the RDPOD basis sets, adaptive RDPOD (ARDPOD), which allows the number of RDPOD basis functions to appropriately vary between space-angle partitions. In addition, this article presents a method which enables the efficient transfer of information between regions with different angular basis sets.
The capabilities of these ROMs are demonstrated with two problems-one involving only advection and absorption, and another which also includes scattering. It is shown that the reduction in basis sizes and solve times for a given error can exceed an order of magnitude in comparison to previous work on angular ROMs. 31,42 Further significant improvements were achieved through the use of adaptive resolution, which reduced error by up to 2 orders of magnitude compared to the non-adaptive method in some cases. While these reductions are already significant, perhaps their benefits would be greater still when applied to large-scale reactor physics problems with far greater levels of variation in their neutron flux distributions.
The sections of this article are arranged as follows. Section 2 presents the BTE together with its general angular and spatial discretizations. Section 3 presents the development of the ROM, and describes the communication between different basis sets, as well as the adaptive angular resolution procedure. Section 4 presents two numerical examples which demonstrate the capabilities of each method. Section 5 completes the article by presenting its conclusions.

THE BTE AND SPACE-ANGLE DISCRETIZATION
The steady-state, mono-energetic BTE governs the angular flux (⃗ r, Ω) in direction Ω at position ⃗ r, The first term in Equation (1) is the advection operator; the second is the removal term, which accounts for losses from both absorption and scattering; the third denotes the isotropic external source; and the fourth represents the scattering into the angle Ω from all angles Ω ′ . The macroscopic scattering and total removal cross sections are given by Σ s and Σ t , respectively. The symbols in brackets indicate the continuous dimensions associated with each variable-⃗ r and Ω for spatial and angular continuity, respectively.

Angular discretization of the BTE
As the process of discretization closely follows that of Reference 42, only a summary is presented here. The Galerkin method is applied to discretize the angular dimension of Equation (1). The angular flux Ψ(⃗ r, Ω) is approximated by a sum of N a angular basis functions  j (Ω) multiplied by the coefficients j (⃗ r), The approximation in Equation (2) is inserted into Equation (1), which is then weighted and integrated over all angles. The basis functions  are also used as weighting functions. Thus Equation (1) is expressed in the angularly discretized form, where A is a vector of matrices (A x , A y , A z ), each containing the components of the unit vectors Ω in the corresponding Cartesian axis. The removal and scattering terms are grouped into the matrix H. The coefficients of the angular expansion are contained in the vector Ψ(⃗ r), and the angularly discretized source term is represented by the vector Q(⃗ r). All matrices are of size N a × N a , and all vectors besides A are of size N a . The components of each matrix and vector in Equation (3) are given in Reference 42. Each class of problem is resolved using the discrete ordinates (S N ) method for n p cases with varying conditions, such as perturbations to the nuclear material cross-sections. The vectors of S N angular coefficients are partitioned into sets according to the angular octant Ω q of each angle, where Ψ q (⃗ r) is a vector containing the N a ∕8 coefficients with associated directions within the octant Ω q .

Spatial discretization of the BTE
The discontinuous Galerkin finite element method is applied to discretize the spatial dimensions of the BTE. 42,43 The angularly discretized flux Ψ(⃗ r) in Equation (3) is approximated as a sum of the N s spatial basis functions  j (⃗ r) multiplied by the coefficients Ψ j , This approximation is inserted into Equation (3), which is then converted to its weak form by weighting with the set of basis functions  (⃗ r) and integrating over the volume of each element, V e . Green's theorem is applied to the advection term, splitting it into an integral over V e and another over the element boundary Γ e . Finally, the boundary term is split into inflow and outflow components, as first order upwinding is used to calculate the flow at element boundaries. 5,44 The result is a full order discretized formulation of the BTE, wheren is the unit vector normal to the element boundary; Ψ out j is the outflow, given by the angular flux vector of the element in question; and Ψ in j is the inflow, given by the angular flux vectors of the element's upstream neighbors. The matrices A in = (A in ⋅n) and A out = (A out ⋅n) are formed to pass the correct incoming and outgoing information, respectively, through the element's surface. For S N , they are produced by simply retaining the negative and positive diagonal elements, respectively, of the matrix (A ⋅n). In the general case, a Riemann approach can be employed. 45 Equation (6) gives the general form of the BTE after discretization in both space and angle, and can therefore be applied to any arbitrary angular approximation.

THE DISCONTINUOUS ANGULAR ROM
This section presents a new angular POD method that forms independent angular basis sets which are each optimized to resolve a specific partition of the space-angle phase-space of the BTE. The angular domain Ω and spatial domain V are partitioned into the two sets, where Ω q and V r represent subsets of the angular and spatial domains, respectively, as illustrated in Figure 1. The DPOD method described in Reference 42 corresponds to the use of a single spatial region spanning the domain, that is, n r =1. If a single angular partition spanning the full sphere is also used, then the standard angular POD method described in Reference 31 is recovered.
In line with the DPOD method, the angular domain is partitioned into eight octants. Each spatial partition may in principle contain any arbitrary set of n r elements in the spatial domain. A cross-product of the two partitions in Equation (7) forms the space-angle partition of the complete phase-space, which is given by, where the phase-space partition spanning the angular octant Ω q and spatial region V r is denoted Z q,r . The POD functions are formed for each partition via the method of snapshots. 46 The angular coefficients from the S N model are used to form the snapshot matrices for each partition. The angular coefficient vectors Ψ are partitioned by octant as shown in Equation (4), and each component Ψ q is assigned to a spatial region V r based on its position. Separate snapshot matrices S q,r are formed for each partition Z q,r , resulting in 8 × n r matrices.
The snapshot matrices S q,r also include vectors of angular coefficients from the elements adjacent to V r , as this helps to preserve information when mapping between regions. This is shown in Figure 1B, where the region enclosed in red takes snapshots from all elements shaded in blue. For each Z q,r , the associated snapshot set is therefore defined as, where each Ψ q,r,i is a vector of size N a ∕8 containing the angular coefficients of the ith snapshot associated with the partition Z q,r . The term n s , which may vary for each snapshot set, denotes the total number of snapshots in S q,r . This is given by n s = N r × N p , where N r is the number of FEM basis nodes in partition V r , including the neighboring nodes as shown in blue in Figure 1B, and N p is the number of problem variations used to train the ROM.
The RDPOD basis sets for each Z q,r can now be formed through the SVD of each snapshot matrix, where U q,r and V q,r are unitary matrices of sizes N a ∕8 × N a ∕8 and n s × n s , respectively. The column vectors of U q,r contain the optimized basis vectors that best represent the snapshot data, ordered such that the first n a columns form the optimal n a basis vectors in the Frobenius norm. Here, n a is the same for all q and r, though this is not required. A method for varying the number of basis functions by region and octant is described in Section 3.1. The RDPOD basis matrices  q,r are formed by truncating each snapshot matrix such that only the first n a columns are retained. The fraction of the information in U q,r which is retained in  q,r can be determined from the singular values, where I q,r varies from 0 to 1, with 1 being total capture of the snapshot information. The matrices  q,r can be used to map the angular coefficients between the full and ROMs through the relationship, where q,r contains a vector of n a coefficients of  q,r , for each node within Z q,r . The combined mapping over all angular partitions in each region V r can be expressed as, Equations (4) and (12) enable this to be compactly written as, which holds for any spatial position ⃗ r inside region r. Substituting Equation (14) into Equation (3) and premultiplying by  T r projects the angularly discretized equations onto the POD space. This projection is applied to each region V r separately, Equation (15) is then spatially discretized as described in Section 2.1, resulting in the fully discretized RDPOD formulation of the BTE. The spatially discretized forms of Equation (15), equivalent to projections of Equation (6) onto the RDPOD bases, are constructed separately for each spatial partition V r using their own optimized RDPOD angular basis sets. Communication between the elements within each partition, and between neighboring partitions, is implemented through the surface integrals of Equation (6). As mentioned, a Riemann approach can be employed to form the incoming and outgoing matrices inside the surface integrals for a general angular discretization. They can also be obtained through preand post-multiplication of the full order matrices A in and A out by the RDPOD mapping matrices  T r and  r . For adjacent elements within the same region V r , the incoming and outgoing surface matrices of Equation (6) are given by, respectively. To obtain an element's incoming surface information from an adjacent element belonging to a different region, say V r ′ , one must use the correct mapping to account for the fact that the incoming vector employs a different RDPOD basis. This is simply achieved by mapping the RDPOD coefficients r ′ from the incoming element to the full model space, applying the (S N ) incoming advection operator, then projecting the resulting vector onto the basis of region V r . That is, the incoming advection matrices in the RDPOD formulation are given by, The matrices A in r,r ′ are precomputed for each pair of regions r and r ′ with a common border. They can then be employed in Equation (6), where the incoming surface integral can be explicitly written as, As stated previously, in the case of a single spatial region spanning the domain, the DPOD method presented in Reference 42 is recovered. Using a single spatial region together with a single angular region spanning the sphere, the original angular POD method in Reference 31 is recovered.

Adaptivity in angle
This section presents an adaptive algorithm using the RDPOD basis functions, known as adaptive RDPOD (ARDPOD). Instead of utilizing the same number of basis functions n a throughout the domain, each partition Z q,r has an associated number of basis functions n q,r , which can be modified independently. The adaptive algorithm uses an error metric to estimate the effect of increasing n q,r for each Z q,r , and adds basis functions where they are likely to minimize the total error. The contribution of each basis function to the solution is dependent only on the magnitudes of its coefficients. The basis functions form a hierarchical set, which implies that the coefficients of each successive function will tend toward zero as the approximation converges. It can thus be inferred that additional basis functions are likely to be most beneficial in partitions where the coefficient of the final basis function currently included is large. However, successive coefficients have been observed to oscillate between positive and negative about zero, and so a single small coefficient does not guarantee convergence. The error metric therefore considers the final two coefficients.
An initial solution is required in order to calculate the error metric and begin the adaptive process, and so the problem in question is first solved with, n q,r = 2, ∀q ∈ {1, 8}, ∀r ∈ {1, n r }.
Next, the relative contribution from the final two basis functions is calculated for each partition Z q,r , F q,r = ( q,r,n q,r −1 + q,r,n q,r ) where q,r,i denotes the ith angular coefficient in the partition Z q,r . A threshold value is set to some initial, relatively large value, typically 1. A variable known as the threshold divisor, , is set to some value greater than 1. In this article, = 2. Each iteration, basis functions are added according to the equation, The number of basis functions added in this process, n Δ , is counted and compared to the minimum number of basis functions to add per iteration, n + . If n Δ < n + , then = ∕ and Equation (21) is applied again with the new value of . The process repeats until n Δ > n + , at which point the next iteration can begin. Once n q,r has been adjusted, the process repeats until a desired number of basis functions is reached. The complete adaptive method is described by Algorithm 1.

Algorithm 1. Adaptive RDPOD
n q,r = 2 for all q and r. /* Iterate until the desired number of basis functions is reached. */ while sum(n q,r ) < maximumBasisFunctions do Solve the ROM. Calculate F q,r . n Δ = 0; while n Δ < n + do /* Loop through partitions and add basis functions. */ for q = 0 to 8 do for r = 0 to n r do if F q,r > then n q,r += 1; n Δ += 1; end end end /* Decrease the threshold value if necessary. */ if n Δ < n + then = ∕ ; end end end

NUMERICAL EXAMPLES
In this section, two numerical examples are presented in order to compare performance of DPOD, RDPOD, and ARD-POD. Uniform quadrilateral FEM spatial meshes are employed, using discontinuous bilinear basis functions. The full order method employs the S N discretization, with a sufficiently high angular resolution to ensure that the solutions have converged in angle. A sweep based solver is employed for all methods, which solves for the unknowns of each mesh element in turn, following the path of information flow. Standard sweeping is employed for the full order model, where each ray is resolved individually and elements visited in the order of information flow along the direction of the ordinate. 47,48 For the ROMs, as the angular coefficients are heavily coupled, all angular coefficients are solved simultaneously for each element. Four sweep directions are used for the two dimensional problems presented here when structured meshes are applied. Elements are swept from left to right, top to bottom (and then in reverse), and top to bottom, left to right (and then in reverse). However, the focus of this article is not on developing the most efficient solver technology to resolve the various angular models. In fact, there are several ways to increase efficiency, but each would require substantial development and warrant a separate article. So, while solver times are provided, these should be considered conservative estimates, and it is highly likely that reductions can and will be made in future. Thus the main purpose of this analysis is to demonstrate the methods' increased accuracy for a given basis size when compared with DPOD, which has previously been demonstrated to outperform both S N and standard angular POD by the same metric. 42 In the examples presented, the spatial partitions used to construct the ROMs are of regular structure, despite Equation (8) allowing for any arbitrary partitioning to be used. While the ideal implementation would minimize the variation in angular flux profiles within each partition, generating partitions in this manner would require significant additional work. This article will instead focus on demonstrating that even regularly shaped partitions can drastically improve modeling efficiency. This demonstration therefore forms a foundation for future work on the generation of optimal partitions for further improvements to efficiency.
Finally, it is sufficient to say that the ROM's computer memory requirements are relatively small. The memory usage increases linearly with the number of regions, n r , since each sub-region requires its own set of discretized angular matrices. However, as shown in the numerical examples, the angular size of the reduced system should be typically small, on the order of a few tens. Thus, taking a large angular ROM of size, say, 100, and using a reasonably large set of partitions, say, n r = 1000, a matrix free solver would require less than 0.6 GiB RAM. Clearly, if parts of the matrix resulting from the spatially discretized equations had to be stored then the memory requirements would increase substantially. While a matrix-free solver is most desirable to avoid large memory requirements, particularly for unstructured finite element meshes, this has yet to be developed. In this article, a uniform structured mesh is used, which enables a single element-wise discretized matrix (from the element-wise discretized equations 6) to be stored for each partition. Again using the example of 1000 partitions and 100 angular POD functions, only 1.5 GiB of memory is required for storage in this case.

4.1
The dog-leg duct problem The first example is a dog-leg duct problem. 49 Figure 2A shows a schematic of the domain. Region 1 is a source, region 2 is the duct, and region 3 is a heavy absorber. Vacuum boundary conditions are applied to the top and right boundaries, and reflective boundary conditions to the bottom and left boundaries. The spatial domain is discretized with a 140 × 180 mesh of discontinuous linear quadrilateral elements. The full order solutions used for snapshots and error calculations employed the S 50 angular discretization. Figure 2B shows the scalar flux distribution of the S 50 solution to the interpolation problem. Table 1 lists the material cross sections for the training and test problems. The snapshot matrix was formed from all four training solutions, and the resulting RDPOD bases were used to solve the test problems. The first test problem is referred to as the interpolation problem, as its material properties are within the range for which snapshots were produced, and the second is referred to as the extrapolation problem, as its material properties lie outside of this range. Figure 3 compares the L 2 -norms of the relative angular flux errors for the dog-leg duct interpolation and extrapolation problems using DPOD and RDPOD, with varying numbers of basis functions and regions. The results show that RDPOD consistently reduces error compared to DPOD for a given angular basis size. Furthermore, increasing the number of spatial partitions consistently reduces the error. For the largest set of spatial partitions, the error is reduced by approximately 1 order of magnitude with 16 basis functions. The reduction in error compared to DPOD continues to grow as the angular resolution is increased further.
F I G U R E 2 Schematic (A) and S 50 scalar flux solution (B) for the dog-leg duct interpolation problem. Region 1 is the source, region 2 is the duct, and region 3 is a highly absorbing material.

Problem
Material Source (cm −2 s −1 ) a (cm −1 ) s (cm −1 ) Note: Material 1 is the source, material 2 is the duct, and material 3 is a neutron absorbing material.  Figure 4 presents the L 2 -norm of the relative angular flux error against the solve time for both interpolation and extrapolation problems, for varying numbers of basis functions and regions. Similar trends are observed, with RDPOD reducing the solve time required to achieve a given error. For the largest spatial partition sets, the solve times are reduced by over an order of magnitude with a relative error of around 10 −5 . The trends in the graph indicate that greater efficiency will be achieved for lower error tolerances. Figure 5 presents the relative angular flux error versus the mean number of basis functions per node for DPOD, and for RDPOD and ARDPOD with varying numbers of regions. It is shown that ARDPOD significantly reduces the error for a given angular basis size. With just 12 basis functions, the error was reduced by approximately half an order of magnitude compared to RDPOD, and by more than an order of magnitude compared to DPOD. By 84 basis functions, the reductions in error for a given basis size had increased to almost 2 orders of magnitude compared to RDPOD and approximately 2.5 orders of magnitude compared to DPOD. Figure 6 presents the relative angular flux error against the solve time for DPOD, RDPOD, and ARDPOD. The adaptive solve times show the time taken to complete the final iteration of the adaptive process, starting from a zero solution. This gives an indication of the performance of ARDPOD once the adaptive stage is complete and an efficient basis function distribution has been generated. It is shown that the adaptive method drastically reduces the solve time required to reach a given level of error in comparison to both DPOD and RDPOD. At around 20 s of solve time, the ARDPOD errors are reduced by just over an order of magnitude in comparison to RDPOD, and by around 2 orders of magnitude when compared to DPOD. Figure 7 presents the relative angular flux error against the total solve times using ARDPOD. The total solve time is the time to complete the entire adaptive solution process. While the ARDPOD solver is far from optimized, the graphs still show that ARDPOD reduces the angular flux error by up to an order of magnitude for a given solve time when compared to DPOD, and that the error exhibits a modest increase in order of convergence. Figure 8 shows the number of basis functions associated with each octant in each spatial region, for the dog-leg duct extrapolation problem. The figures show that basis functions have predominantly been added to high flux regions, such as inside the duct and at its borders. This demonstrates that the adaptive method is successfully increasing the resolution in regions where additional basis functions are likely to be most beneficial, and thus explains the method's success in reducing the error for a given number of basis functions.

The Watanabe-Maynard problem
The second example is the Watanabe-Maynard problem. 50 Figure 9A depicts a schematic of the spatial domain. Region 1 is the source, region 2 is a void, and region 3 is a moderate absorber and scatterer. Table 2 shows the material properties of each region for the training and test solutions. Reflective boundary conditions are applied to the bottom and left boundaries, and vacuum boundary conditions are applied to the top and right. The spatial domain is discretized with a 160 × 160 mesh of discontinuous linear quadrilateral elements. The S 50 angular discretization was employed to produce full order angular solutions. Figure 9B depicts the scalar flux distribution of the S 50 solution to the interpolation problem.    Table 2 lists the material cross sections for the training and test problems. As previously, the snapshot matrix was formed from all four training solutions, and the resulting bases were used to solve both an interpolation and an extrapolation problem. Figure 10A,B depicts the L 2 -norm of the angular flux error for the Watanabe-Maynard problems, for both DPOD and RDPOD with varying numbers of regions. As the figure shows, increasing the number of RDPOD regions consistently reduces the angular flux error for a given basis size, excluding 4 basis functions. With just 8 basis functions, the maximum reduction in error is an order of magnitude, which increases steadily to 2 orders of magnitude by 84 basis functions. Figure 11 plots the relative angular flux error against the solve time for DPOD and RDPOD. It is again shown that increasing the number of spatial regions consistently decreases the error for a given solve time, with a maximum reduction of more than 2 orders of magnitude.
In some cases, moving from 4 to 8 basis functions decreased both the error and the solve time, which was not expected. This is likely because the iterative solver converged slowly for the bases with 4 functions, so adding a basis function to each octant not only reduced the error as expected, but also allowed the solver to converge in fewer iterations and thereby reduced the solve time. Figure 12 plots the relative angular flux error against the mean number of basis functions per node for all three methods with varying numbers of regions. As the figure shows, ARDPOD reduces the error compared to RDPOD in all cases, with a peak reduction of approximately an order of magnitude. Compared to DPOD, ARDPOD reduces the error by approximately 2.5 orders of magnitude. The figure also shows, once again, that increasing the number of regions consistently reduces error for both RDPOD and ARDPOD. Figure 13 compares the relative angular flux error to the solve times for DPOD, RDPOD, and ARDPOD. As explained in the discussion of Figure 6, the solve times given are for a single iteration, and are intended to be indicative. The cumulative solve time is not considered, as the solver has not yet been optimized for this purpose. The figure shows that ARDPOD drastically decreases the error for a given solve time, with a peak reduction of more than 3 orders of magnitude. This demonstrates that, while RDPOD and ARDPOD both reduce error compared to DPOD for a given basis size, they do not significantly affect the solve time. Similar anomalies to Figure 11 are seen with 4 basis functions per node, likely for the same reason. Figure 14 presents the relative angular flux error against the cumulative solve time for ARDPOD. As mentioned in the discussion of Figure 7, the adaptive method is not yet fully optimized, and significant improvements in this metric are likely possible. Despite this, the graphs show that ARDPOD performs at least as well as DPOD in the case of 2 × 2 regions, and significantly better with more regions. At best, it offers over an order of magnitude reduction in error with the same solve time as DPOD.

CONCLUSIONS
This article has developed upon the method of DPOD proposed in a recent article. 42 A new ROM for the angular dimension of the BTE, known as RDPOD, has been described. The novelty of RDPOD lies in its separation of the spatial domain into multiple regions, each of which has its own optimized DPOD basis set. A method of projecting flux between each reduced order basis without full order calculations has also been derived and implemented. Finally, an adaptive algorithm based on the RDPOD bases has been presented. The RDPOD method is shown to consistently decrease the relative angular flux error for a given solve time and basis size when compared to DPOD, by an amount proportional to the number of spatial regions. This was the expected result, as increasing the number of regions reduces the number of elements per region, and therefore allows each basis function to be better optimized for the angular flux profiles it represents. The relative angular flux error for a given number of basis functions was reduced by up to an order of magnitude for the dog-leg duct problem, and 2 orders of magnitude for the Watanabe-Maynard problem, demonstrating that the method can benefit both advection and scattering problems. Similar results were observed when comparing the error to the solve time, as RDPOD did not significantly affect the solve time for a given basis size compared to DPOD.
The adaptive method, known as ARDPOD, was demonstrated to further reduce error for a given basis size compared to RDPOD. The effect of adaptivity was significant-for both numerical examples, the error was reduced by up to 3 orders of magnitude compared to DPOD. Compared to RDPOD with the same number of regions, ARDPOD typically reduced the error by up to an order of magnitude, and more in some cases. The same was true when comparing the error to the solve time for each iteration. However, as previously mentioned, the cumulative time to reach a given adaptive stage was not optimized. Additional optimization of the adaptive method could significantly reduce the cumulative solve time, but this is left for future work.