IMPROVED ANISOTROPIC LINEAR SOURCE FORMULATION FOR MULTIPHYSICS PROBLEMS

This work seeks to extend an existing formulation of the method of characteristics with linear source approximation for problems with dynamic cross sections. The previous formulation eliminated cross section dependence of precomputed coefﬁcients for systems with an isotropic source. The method is extended to include a formulation for spatially ﬂat anisotropic scattering that eliminates cross section dependence of precomputed coefﬁcients without adding additional operations; increasing efﬁciency in multiphysics simulations where cross sections can be subject to change. The new formulation is implemented in the MPACT code and tested on two problems: 3D transport assembly calculations using MPACT’s 2D/1D method and a 3D assembly with T/H feedback using MPACT’s 2D/1D method coupled with COBRA-TF. This work demonstrates that the new linear source formulation allows for the number of mesh ele-ments to be signiﬁcantly reduced while maintaining accuracy, leading to shorter run-times for 3D cases with ﬁxed cross sections, and substantial reduction of memory usage for 3D cases with ﬁxed cross sections. The multiphysics calculations show similar runtimes for the same accuracy with signiﬁcant reduction of memory. For similar accuracy, the method proved effective in reducing memory usage by, on average, 30% for 3D problems and 21% for multiphysics problems.


INTRODUCTION
The method of characteristics (MOC) is commonly used in reactor physics computations as a means to solve the neutron transport equation. The flat source (FS) approximation of the MOC equations is the most common formulation of this method. In this formulation, the source is assumed to be spatially flat in each spatial region, though it is not assumed to be isotropic. In [1] a linear source (LS) formulation based on the computation of track-based spatial moments over source regions to obtain the LS expansion coefficients was proposed. It was demonstrated that for typical reactor physics problems, spatial expansion of the anisotropic portion of the flux is unnecessary and provides little improvement over treating it as flat. This led to the proposal and testing of a source that was flat for the anisotropic portion and linear in the source regions for the isotropic portion. Originally called the LS-P0 source, it will be referred to as the LIFA (linear isotropic, flat anisotropic) source in this work. The isotropic LS and LIFA formulations involved a number of cross section dependent coefficients that could be precomputed for use in each MOC sweep and source iteration. However, MPACT [2] is primarily intended for multiphysics problems. In these problems cross sections are dynamic (i.e. not constant) over the course of the calculation. This means that the precomputed coefficients must be recomputed for every multiphysics iteration. With this in mind, the original formulation was improved [3] for isotropic scattering in order to allow for a cheaper computation each multiphysics iteration without relying on precomputed cross section dependent coefficients. This formulation is here extended to spatially flat anisotropic scattering.
This work involves the introduction of this extension of the improved formulation for spatially flat anisotropic scattering, the implementation of this improved LIFA formulation in MPACT, and analysis of results using this formulation. Results show that the LIFA source noticeably improves the calculation of both eigenvalue and peaking factor on a coarse mesh, and that even with a finer reference mesh for comparison the FS and LIFA methods, we observe that the LIFA method still often improves the calculation for comparable run times and substantial reduction in memory.

FORMULATION
Beginning with the characteristic transport equation for a 2D system: where ψ g mki is the the angular flux in group g along track k in direction m in computational cell i. Σ g t,i is the total cross section, and s m is the variable segment length in direction m. The index m is a combined index of of polar and azimuthal indices p and a, respectively. Let the source be a linear function of space in computational cell i, so that q g mi (s m ) = q g mi (x, y). Assuming that the linear variation of the source is isotropic, but the flat portion is anisotropic, then the source becomes: Next we define the inner product: where V i is the radial area of the computational cell and r = (x, y). Note that these integrals are discrete summations in practice. The flat portion of the source can then be expressed as: where R r l are the real spherical harmonics defined in [4]. R r l ψ g i are the angular moments of the flux such that R 0 0 ψ g i = ψ g i = φ g i , and φ g i is the integrated scalar flux in computational cell i. Notice that the flat portion of the source is generally anisotropic. The linear portions of the source are then expressed as the following: and similarly for y. Here the linear portions of the source are isotropic. The source coefficients φ g i,x and φ g i,y are solved by inverting the following 2 × 2 matrix: where M i,xy = xy i 4π and M i,xx and M i,yy are defined similarly. The MOC equation is now: where the track centered portion of the source is defined as: with (x c mki , y c mki ) being the local coordinate of the track center. The coefficient to express the source variation along a track is: The angular moments of the flux can then be computed as: where ∆ψ g mki is the flux difference between the start and end of the track, and w m is the discrete ordinate weight multiplied by the sin of the polar angle and the ray width. The spatial moments of the flux are computed as: for the x moment, and similarly for the y moment. Here ψ g,out mki is the outgoing track flux and x in mki is the incoming track x coordinate.

A QUICK NOTE ON OPTIMIZATION
Traditionally, MPACT has accumulated the entire angular moment during each sweep. This means that the following calculation is done for each segment: A method of speeding up the anisotropic source calculations for MOC algorithms was proposed by [5] and allowed additional optimization of both the LIFA and FS methods. All calculations in this paper include this optimization. The fundamental concept of the optimization is to instead accumulate just the flux change ∆ψ g mki on each segment. So then at the end of the sweep only the following computation is required: Here the boxed term has already been accumulated during the sweep. This means that for a given phase region of angle m, group g, and region i, the multiplication by the spherical harmonics and other coefficients only has to be evaluated once compared to the number of segments in that discrete phase region with the original accumulation method. Because the LIFA method treats the angular moments as spatially flat, the optimization is identical in implementation between the FS and LIFA methods so the comparisons remain consistent.
For shared memory parallelism, it has been observed that this optimization works best for when the number of threads sweeping over the rays is low. For a large level of shared parallelism of the ray sweeps, it may be that most processes do not compute multiple segments for every source region, and indeed some processes may not compute any segments in some regions. To avoid this, a simple conditional can be implemented to ignore calculations for regions where a thread has stored zero integrated angular flux.

B&W Critical Experiments
The B&W 1484 Core 1 critical model from [6] served as the 2D/1D benchmark calculation. This model was essentially the 2D model extruded to a height of 72.45 cm with vacuum boundary conditions at the top and reflective boundary conditions at the bottom. All models were run with 64 azimuthal angles and 3 polar angles on the unit half sphere. All models were run with 0.05 cm ray spacing using MPACT's 51 group cross section library. There were a total of 3 different meshes. The fine reference mesh had 208 source regions in each fuel pin and 256 source regions for each moderator pin. The intermediate, or "mid", mesh had 48 regions in each fuel pin and 36 regions in each moderator pin. The coarse mesh had 32 regions in each fuel pin and 4 regions in each moderator pin. All cases were run with an axial mesh where each slice is 10.35 cm tall making for a total of 7 slices in the system. The axial mesh parameters were not varied since the linear source in this implementation is only linear in the radial plane and therefore will not have an affect on accuracy for varying axial meshes. All runs for these cases were done with 6 MPI processes each running 4 shared memory threads with OpenMP, so each of the available 24 cores were being fully utilized.

VERA progression problem #6
The VERA progression problem 6 [7] served as the 3D multiphysics benchmark. This is a steady state hot full power assembly problem with T/H feedback. The T/H portion of the simulation is computed using COBRA-TF [8]. The geometry is a single Westinghouse 17x17-type fuel assembly with a height of 406 cm with reflective BC in the plane and vacuum above and below the assembly.
All models were run with 64 azimuthal angles and 2 polar angles on the unit half sphere. All models were run with 0.03 cm ray spacing using MPACT's 51 group cross section library. There were a total of 3 different meshes. The fine reference mesh had 208 source regions in each fuel pin, 128 source regions for each guide tube pin, and 225 regions for each moderator pin. The intermediate, or "mid", mesh had 48 regions in each fuel pin and 40 regions in each moderator or guide tube pin. The coarse mesh had 8 regions in each fuel pin, 7 regions in each guide tube pin, and 9 regions in each moderator pin. All calculations for these cases were performed with 58 MPI processes running a single thread.

RESULTS
The results analyzed for this work focused on the global fission power peaking factor (F q ), and the system multiplication factor, k eff . These quantities represent perhaps the two most important global engineering quantities of interest of reactor systems, so their accuracy is paramount.

2D/1D results
The proceeding tables make use of the following terminology. Fine is the reference results with FS, Mid is the FS run on the intermediate mesh, Coarse is the coarse mesh run with FS, and LIFA is the coarse mesh run with the LIFA source. The N term refers to the order of anisotropic scattering. In all cases the Fine Mesh is considered to be the reference result as it represents a significantly fine spatial mesh. Different anisotropic scattering orders are not to be compared to each other when comparing LIFA and FS results, but rather serve as a demonstration of the fact that anisotropic scattering can be important for accurate core simulations in some circumstances.
Reference eigenvalues and the associated errors in pcm are shown in Tab. 1 for all cases. Lower errors for the LIFA source are observed for this test. The Mid Mesh is shown to have lower error than the Coarse Mesh for the FS. In all cases it is observed that the LIFA error is smaller, often by a substantial amount. Notice that the LIFA results in still lower errors even for higher order scattering up to P3. The memory usage for each case is shown in Tab. 4. These results show that LIFA takes a little less than twice as long as FS for an identical mesh but uses about the same memory. The intermediate mesh FS cases still showed worse errors than the LIFA method on the coarse mesh. The FS also took much longer for higher order scattering compared to LIFA and cost a little less than twice as much memory. One important thing to notice is that the FS, regardless of mesh, scales significantly worse in both time and memory as the order of anisotropy increases compared to the LIFA calculations. This is due to the fact that for each ray-segment the LIFA requires significantly more computational work compared to the FS since it is also computing linear moments. The accumulation of additional anisotropic moments is performed per-cell instead of per segment. The number of ray-segments is usually much larger than the number of cells, and thus the LIFA solver has much more work to do for the base (isotropic) case. Therefore, adding angular moments is relatively cheap compared to the base isotropic case. Since the FS calculation has less segment work to do, the per-cell angular moment calculation represents a larger fraction of the overall work.

Multiphysics results
These results make use of the same terminology as the 2D/1D results. Reference eigenvalues and the associated errors in pcm are shown in Tab. 5 for all cases. Lower errors for the LIFA source are observed for this test compared to the coarse mesh. The Mid Mesh is shown to have lower error than the Coarse Mesh for the FS. However, the difference between the LIFA error and the coarse mid mesh error is observed to be quite small to the point where the difference is likely cancellation of errors.

CONCLUSIONS
Results show that the improved formulation of the LIFA method can be used successfully for anisotropic scattering. In the 2D/1D results, the eigenvalues for LIFA were on average 100 pcm better, and the peaking factors were on average 49% better for the same mesh. Comparing the intermediate mesh to the coarse mesh LIFA in the 2D/1D calculations still showed better results from the LIFA calculation, and the coarse mesh LIFA calculation used on average 30% less memory and took on average 30% less runtime. For these 3D models, we observed much better scaling for LIFA for increased anisotropy compared to FS for a given mesh.
In the multiphysics results, the eigenvalues for LIFA were on average 60 pcm better, and the peaking factors were on average 73% better for the same mesh. Comparing the intermediate mesh to the coarse mesh LIFA in the multiphysics calculations again showed better results from the LIFA calculation. The coarse mesh LIFA calculation used on average 21% less memory and only took on average 4% more runtime. For these cases, the plurality of wall time is spent in COBRA-TF, and the CMFD calculations in MPACT, neither of which is effected by mesh size or linear source approximation, so comparing run times is difficult since even one additional iteration will likely result in a noticeably longer time. For these 3D models, we still observed much better scaling for LIFA for increased anisotropy compared to FS for a given mesh.
In conclusion, we see that for larger 3D and multiphysics problems, the LIFA method has the potential to decrease the costs of a calculation by coarsening the mesh, while maintaining comparable accuracy. As memory is often the limiting factor in many of these calculations, the consistent reduction of memory for the same accuracy results through the use of LIFA on a coarser mesh should allow for larger and more complex problems on a fixed set of resources. Now that the LIFA method formulation is implemented in MPACT, future VERA runs should be able to take advantage of coarser meshes for the neutronics portion of the multiphysics problems. The fact that it is implemented using the improved multiphysics formulation will additionally allow for more computational savings for multiphysics problems since all precomputed parameters are cross section independent.