VERIFICATION OF THE DIFFUSION AND TRANSPORT SOLVERS WITHIN DIF3D FOR 3D HEXAGONAL GEOMETRIES

The DIF3D code (DIFfusion 3D) has been a workhorse of fast reactor analysis work at Argonne National Laboratory for over 40 years. DIF3D was primarily built in the late 1970s as a three-dimensional multigroup diffusion equation solver operating on semi-structured grid geometries. In the mid-1990s, transport capabilities needed for high-leakage reactor configurations were added to DIF3D with the variational anisotropic nodal transport approach. Recent reactor design activities at Argonne are requiring that a thorough verification of the Argonne Reactor Computation (ARC) codes be performed. With DIF3D being central to the entire ARC system, the verification efforts are focused on the 3D Cartesian, 3D triangular, and 3D hexagonal core geometry options of DIF3D. Validation activities, while needed for the ongoing design activities at Argonne, are handled at a project-specific level. This paper summarizes the verification work so far on the forward and adjoint forms of the fixed source, inhomogeneous fixed source, and k-eigenvalue steady state transport and diffusion equations as implemented specifically for 3D triangular and hexagonal geometries in DIF3D. Since analytic solutions of the neutron diffusion and transport equations are either limited in scope or not possible, this verification required multiple tiers of problems unique to each solver and geometry type, each testing features independent and complementary arguments for why this separate testing of functionalities is acceptable. This separate testing was also supplemented with a high-level integral check of each the diffusion and transport capabilities and applicable geometries.


INTRODUCTION
The DIF3D code (DIFfusion 3D) was primarily built as a three-dimensional solver of the multigroup diffusion equation for semi-structured grid geometries using a finite difference spatial differencing methodology [1]. This capability still exists and is commonly referred to as DIF3D-FD. Later, in the early 1980s, a transverse integrated nodal methodology, referred to as DIF3D-Nodal [2], was built into DIF3D to improve the performance on large scale reactor problems. Later work, motivated by the research and development to support fast spectrum reactors that increase proliferation tolerance by eliminating blanket assemblies, the DIF3D-VARIANT solver [3] was added to DIF3D in the mid 1990s, and later updated [4], that extends the nodal concept to transport via a variational three-dimensional transport method applying spherical harmonic expansions in angle.
The present DIF3D capability can treat slab and cylindrical 1D domains, Cartesian, hexagonal, and R-Z 2D domains, and Cartesian, hexagonal-Z, triangular-Z, and R-Z-θ 3D domains. The DIF3D-FD diffusion solver can be applied to every geometry type in DIF3D except for hexagonal and hexagonal-Z geometries, where a triangular-Z mesh geometry is instead used to represent hexagonal reactor designs. The DIF3D-Nodal diffusion and DIF3D-VARIANT solvers only support 2D and 3D dimensional Cartesian and hexagonal geometries.
Recent reactor design activities at Argonne are requiring a thorough verification of the Argonne Reactor Computation (ARC) [5] codes (of which DIF3D is paramount) be performed to modern software quality assurance standards. In this context, software verification is defined as the act of confirming that the equations intended to be implemented within the software of interest were implemented correctly through use of primarily analytical solutions to the governing equations, with the usage of numerical or other methods when an analytical solution is not feasible. Note that the design work that requires ARC verification also required validation, however this validation work is not the subject of this paper.
More specifically, the required verification includes confirming: (1) DIF3D's ability to correctly translate the user's model in to DIF3D's preferred format; (2) that options planned for use have the desired effect; and (3) the correctness of the eigenvalue, fixed-source, forward, and adjoint solvers in DIF3D-FD and DIF3D-VARIANT for 3D Cartesian, triangular-Z and hexagonal-Z core geometry options. Note that the DIF3D-Nodal diffusion capability will not be discussed further.
This paper summarizes the verification tasks performed to verify the forward and adjoint forms of the fixed source, inhomogeneous fixed source, and eigenvalue forms of the transport and diffusion equation as implemented specifically hexagonal-Z and triangular-Z geometries in DIF3D. The paper will place emphasis on the most common use-case for DIF3D, the forward eigenvalue solver verification. The remainder of topics are required for the complete verification but will be discussed via separate work.

VERIFICATION STRATEGY
No non-trivial analytic solution of the neutron diffusion and transport equations can be generated for the triangular-Z and hexagonal-Z geometries for multigroup, multi-region, fixed source and eigenvalue reference solutions. The verification of DIF3D therefore requires multiple tiers of problems unique to each solver and geometry type, each testing features independently and present complementary arguments for why that separate testing of functionality is acceptable. The parameters in the verification tests are chosen so the reference solution is only reasonably obtained through a correct implementation of the governing equations and algorithms in DIF3D. For example, problems with anisotropic scattering will use cross sections such that the expected solution is significantly different if the anisotropic scattering data is ignored by DIF3D. To add confidence in the separate testing as well as to provide confidence that the chosen problems did not fortuitously yield a correct solution, a final integral check of each the diffusion and transport capabilities and applicable geometries are also incorporated in the plan. These final integral problems utilize previously published benchmark problems with externally published solutions.
The tiers of cases used for this for the DIF3D-FD triangular-Z and DIF3D-VARIANT hexagonal-Z are shown in Table I. These tasks show the high-level type of problem and the source of the reference solution that was used for this tiered verification. The solution types include analytical/hand calculations, the method of manufactured solutions [6], comparison to published results for equivalent models, and comparison with a higher-fidelity Monte Carlo neutron transport solver, MCNP Version 6.2 [7].
As this is verification of the DIF3D solver and not the multigroup cross section generation process, the MCNP reference solutions were multigroup MCNP computations utilizing identical cross sections to the DIF3D cases. Using MCNP as the reference solution source naturally requires that the relevant MCNP6.2 capabilities be verified. Therefore, the MCNP models utilized only the basic features, such as only using planar surfaces to describe the geometry, avoiding nested universes, and not taking advantage of weight windows. The MCNP verification is not discussed in this paper.

Integral Verification
Published Results

Multi-Composition Boundary Condition Tests
Comparison to Full-Core

Integral Verification
Published Results, MCNP6.2 * Bold items are those discussed in this paper

VERIFICATION CASES AND RESULTS
The specific tasks for the verification of the DIF3D-FD Triangular-Z and DIF3D-VARIANT Hexagonal-Z solvers are discussed separately in Sections 3.1 and 3.2 below. Note that, for brevity, the only tasks that are discussed are those denoted in Table I with bold text. Each of these sections will begin with an overview of the relevant solver details needed for the verification work. For example, a brief discussion is provided on how DIF3D-FD uses regular triangular meshing to model a hexagonal reactor geometry before discussing any results. After this overview, the verification tasks discussed above are examined and results compared to the reference.
Before proceeding, it should be noted that DIF3D nomenclature refers to spatial regions assigned a specific set of isotopes as "zones". For clarity, this paper instead refers to these as "compositions" consistent with the literature.

DIF3D-FD Triangular-Z Verification
The DIF3D-FD solver is a finite-difference based diffusion solver. The user input for this geometry assumes that the input is a hexagonal lattice that is broken into a user selected number of triangles (6, 24, 54, etc.). The input for triangular-Z geometries is thus limited to specifying the ring and position in the hexagonal lattice along with the axial range. This geometry is converted internally into a triangular-Z mesh and stored in the geometry description (GEODST) file, a standardized file format established by the Committee on Computer Code Coordination [8].
The triangular mesh used to describe the reactor can have outer bounds that are rhombic with the axes at 120 degrees (third-core symmetry) or 60 degrees (sixth-core symmetry), or rectangular (full-core, quartercore, and half-core symmetry). The quarter-and half-core symmetry options are not included in this verification effort. While the algorithm for solving each of these triangular meshes is the same, the coefficients used in creating the solution matrices are unique to each of these triangular mesh types and therefore each must be tested. Regardless of the type of radial domain boundary, each axial mesh must have an identical radial domain.
The left hand side of Figure 1 provides an example of a rectangular domain for a single plane of a threering reactor. The solid lines denote the triangle mesh lines and the red, blue, green, and yellow simply denote different composition types. The right hand side provides a mesh of the same reactor assuming a 60 degree rhombic domain and thus sixth-core symmetry. In either of these meshes the white regions are mesh elements that are not filled when defining the three ring reactor such as this one, though they are still present in the computational mesh. By default, DIF3D treats "background regions" in a manner equivalent to a pure black absorber with no neutron return at the re-entrant corners. The user can also opt to fill all background regions with a single composition of their choice.

Figure 1. Rectangular Triangular-Z Mesh Examples
Triangular mesh refinement is also available to the user. Specifically, the user can request the default mesh (6 triangles per hexagon, as shown in Figure 1) be subdivided by decreasing the spacing between the mesh lines endpoints on the problem boundary. For example, Figure 2 shows the 60 degree rhombic mesh from Figure 1 on the left, and the same model with 2 triangle subdivisions on the right. Here we see that the mesh line spacing has been halved, thus yielding four triangles inside each original triangle, or 24 triangles per hexagon instead of the original 6 triangles per hexagon.

Figure 2. Rectangular Triangular-Z Mesh Subdivision Example
With this primer complete, the following sub-sections provide methods and results for each of the diffusion verification tasks identified earlier.

1-Composition, 1-Group Spatial Solution via the MMS
To verify the diffusion solution on triangular-Z geometries, a modified approach to the Method of Manufactured Solutions (MMS) was applied to verify the fixed source with fissile media (inhomogeneous fixed source) solver [6]. A similar approach was also used for the standard eigenvalue solver but is not presented for brevity.
In the modified methodology, an analytic solution to the diffusion equation is assumed which satisfies the boundary conditions used in DIF3D. The discrete approximation associated with a given DIF3D input is used to construct a fixed source from the diffusion equation that, when used by DIF3D, should result in the discrete flux solution. With sufficient mesh refinement, the DIF3D solution should match the analytic solution that was assumed which demonstrates that the solver methodology for the differential equation of study has been implemented correctly.
The typical one-group eigenvalue diffusion equation is provided in Equation (1). The reference solution for this was constructed by first assuming that the flux was zero on the boundaries of each of the hexagon. This flux solution, with a per-hexagon boundary condition function is provided in Equation (2).
The boundary condition function can be defined as shown in Equation (3). This boundary condition is usable on any size hexagonal geometry grid noting that it will display zero flux solutions within the domain and cause a negative fixed source at those points with respect to the diffusion equation. To introduce peaks in the domain, we use witch functions to define ‫,ݔ(݂‬ ‫,ݕ‬ ‫)ݖ‬ of the form provided in Equation (4). In this equation, the X and Y values are the maximum and minimum x and y values that occur for a given hexagonal grid which for N rings is determined according to ‫ݔ‬ = ܰܲ − ଶ A utility program was created to generate this source and write the appropriate fixed source input. The details of this methodology are not included here. The following values are used for the multigroup cross sections when generating the source: the diffusion coefficient has a value of 2/3 cm; the removal cross section is 0.25 1/cm; and the fission production cross section is 0.252 1/cm.
The overall manufactured solution is provided in Figure 3. This problem is a three-ring hexagonal system with a 40 cm height and 10 cm hexagonal pitch. This small size is driven by the need to reduce the computational and mesh burden due to the considerable mesh refinement necessary. The witch functions are set to have maximum values at axial positions 4 cm and 28 cm and on opposite sides of the domain. The planer slice shown at the right of Figure 3 shows how the two witch and sine functions combine to cause a complex axial flux shape.  The above problem was solved with DIF3D-FD using a series of progressively refined meshes. During the mesh refinement it was observed that as the reference flux solution approaches a local zero the DIF3D relative error grows significantly. This effect is naturally exacerbated as the mesh is refined because a higher fraction of any individual mesh element can contain a zero reference flux. Therefore the error metric used for this case is absolute, vice relative, errors. The following mesh refinements were investigated: 5 triangles per edge and 5 axial meshes, 10 triangles per edge and 10 axial meshes, 15 triangles per edge and 20 axial meshes, and 20 triangles per edge and 25 axial meshes. The absolute errors for the first (5/5) and last (20/25) cases are provided in Figure 4. In this case, one can see that the maximum error in any given flux magnitude is dropping as the mesh is refined. Further, to obtain relative errors, any flux value less than 0.001 was set to zero to eliminate the small trivial fluxes. In this case, the progressive mesh refinement discussed above yielded maximum errors of 38%, 34%, 4.8%, and 3.7%. The results presented here are sufficient evidence that DIF3D-FD is correctly solving the stated problem given the constraints of accuracy that can be applied with the fixed source input.  The preceding results demonstrate that DIF3D-FD solves the triangular-Z diffusion equation correctly. For completeness, a similar exercise to the above was also performed for the DIF3D-VARIANT solver with a 1 st order angular treatment (equivalent to the diffusion approximation). DIF3D-VARIANT was also found to correctly solve the hexagonal-Z diffusion equation. Finally, the above MMS approach is also used to verify the DIF3D peak flux and peak power edits though this is not discussed here.

Integral Benchmark Verification
The purpose of the integral verification test is to compare both the eigenvalue and flux distribution against externally obtained solutions. The VVER-440 model [9] was chosen as the comparison case given the available external solutions in the literature. The full-core geometry for this case is shown in Figure 5. This core has 30 degree symmetry and thus the 30 degree sector of interest is shown. Two group cross sections were provided in [9]. The hexagonal pitch is 14.7 cm. The total geometry height is 300 cm; the active core region is 250 cm and there are 25 cm of moderator regions above and below the active core. A portion of control rods are inserted 150 cm from the bottom.

Figure 5. VVER-440 Benchmark Model
The original reference solution was obtained in [9] by extrapolating the DIF3D-FD keff to the keff that would be obtained when using an infinitesimally small mesh element size. The result was a keff of 1.01132. A whole core power distribution assuming ten axial regions for each assembly was also provided. Table II provides the triangular mesh size (side length of the triangle), the axial mesh size, and the resulting DIF3D-FD keff result converged to an iterative error of 10 -8 . Because convergence with respect to axial mesh refinement is observed at an axial mesh size of 1.25 cm, extrapolation to the infinitesimally small mesh element size was performed based only on the radial mesh refinement. Using the 2.12 cm and 1.70 cm triangle meshing results with 1.0 cm axial meshing leads to an eigenvalue of 1.01155. Using the 1.70 cm and 1.41 cm triangle meshing results with 1.0 cm axial meshing leads to 1.01150. In addition to these linear extrapolations, two least squared approaches were also used. The first used a basis of (1, x, y), yielding an eigenvalue of 1.01135. The second used a basis of (1, x, y, x·y) and yielded an eigenvalue of 1.01168. These are all included in eigenvalue error of +/-20 pcm around the reference solution is assumed to be considered accurate given the variance in the above extrapolation schemes. The power distribution from DIF3D-FD was compared to the reference via a mesh-wise comparison for all axial meshes of all 37 assemblies. In Table III we provide four measures to describe the accuracy of the power distribution. The RMS error is the standard root-mean-square of the mesh-wise error in power density (all meshes have the same volume). The average error is the average of the absolute value of the mesh-wise error in power density. The peak powered mesh error is the error in the mesh with the maximum power density that occurs in axial mesh 4 of assembly 25 or ring 10 position 3. Finally, the maximum mesh error is the maximum error in any mesh from any assembly. Table III includes results from multiple sources and cases: DIF3D-FD with the 1.70 cm triangle and 1.0 axial mesh case from Table II; DIF3D-VARIANT with the P1 angular approximation (i.e., diffusion); a new DIF3D-FD solution with meshing that matches the meshing scheme used for the stated reference solution [9]; and results from references [9] and [10] that do not rely on DIF3D. From these results, it is clear that the updated DIF3D-FD results contain a non-trivial amount of error in it compared with the published reference solution. The DIF3D-VARIANT solution is slightly closer to the reference solution also having a maximum error well under 1%. The DIF3D-FD solution with a mesh that matches the reference interestingly contains non-zero errors in all mesh locations. This implies that either different convergence settings were used in the DIF3D-FD calculations (i.e. iterative error) or the stated reference is also an extrapolated result that as not stated in the literature. Regardless, the error measures are similar to DIF3D-VARIANT and accurate with respect to the reference solution. In reference [9], only a 1.28% average error is reported for ANC-H along with a keff error of 25 pcm. In reference [10], only the maximum mesh error is given for the coarsest mesh approximations although keff solutions for multiple meshes are provided. For the coarsest PARCS solution, the keff error is 45 pcm and the maximum mesh error was reported as 21.7%. The coarsest COREDAX solution had a keff error of 102 pcm and a maximum reported mesh error of 1.28%. For the finest mesh keff results reported, PARCS had a 45 pcm error and COREDAX had a 6 pcm error and thus one can assume that the power distribution errors were significantly reduced as well even though they were not provided. Excluding the PARCS result, these solutions are reasonably consistent with the DIF3D-FD and DIF3D-VARIANT solutions.

DIF3D-VARIANT Hexagonal-Z Verification
DIF3D-VARIANT is a transport solver based on variational nodal methods on 3D Cartesian and hexagonal-Z geometries. Like the diffusion solver, DIF3D internally converts the user's hexagonal reactor description in to a hexagonal mesh of the reactor with one hexagon mesh per ring position from the user input. This hexagonal mesh is extruded with user-set axial mesh spacing and composition assignments.
Given that this is a variational nodal method solver, the DIF3D-VARIANT allows the user to individually set the polynomial orders used for the within node source, within node flux, and node surface leakage. The user also must specify the angular (spherical harmonics) order for both the within node flux and the node surface leakage. Since the spatial functions are represented with polynomials, models with strong spatial gradients will require many spatial orders to accurately represent. Similarly, models with high anisotropy of the flux or leakage (especially near an inhomogeneous source or a vacuum boundary) will require high angular orders to represent the physical angular flux. With finite computing resources, fully resolved spaceangle approximations are unrealistic. This aspect is further exacerbated by DIF3D's memory management system that constrains the problem size below the theoretical limit of the computing resource. As an example, only small and/or coarsely-meshed models can allocate enough memory to combine a 4 th order spatial source, 6 th order spatial flux, 1 st order leakage, and 13 th order spherical harmonics expansion for the angular flux and leakage. This memory limitation is the ultimate challenge in obtaining reference solution answers with the DIF3D-VARIANT solver, however this work aims to show that the solution is in fact converging to the reference before the memory limitation is reached. Table I provided the verification tasks to be performed for DIF3D-VARIANT. For brevity, this paper will focus only on the integral transport verification task. The purpose of this task is to provide a final check for the reasonability of DIF3D solutions for more realistic problems. To do this, the 4 th series of models of the Takeda and Ikeda benchmark set [11] was used including the provided four-group cross sections. The fullcore geometry for this case is shown in Figure 6. These models have three configurations termed "Rods-Out", "Rods-Half", and "Rods-In." In each case, six hexagonal positions are modified (only one is depicted due to symmetry) to simulate a control rod change. It is worth noting that the axial definitions of the control rod geometry are considerably simplified for the purposes of simplifying the benchmark model. The discussion that follows will focus only on the "Rods-In" case since this one led to the worst DIF3D to reference biases. The performance for other cases is similar although better than for "Rods-In". The Monte Carlo solutions published with the benchmark contained 50 to 60 pcm of statistical uncertainty for the reported eigenvalue; these were computed with less than three million histories due to computational limitations of the time. This uncertainty was considered too large for this work and therefore MCNP6.2 was used to recreate the reference solutions. In this case, the goal is to show that the DIF3D results are converging to the MCNP6.2-based reference solution.

Integral Benchmark Verification
The MCNP reference solution was generated here with one million histories per batch, 100 inactive batches, and 16,000 active batches for a total of 16 billion active histories. The resultant value of keff predicted by MCNP was 0.88105 +/-0.00001. The relative errors for all tally bins are less than 0.4%.

Rods-In DIF3D Results
The DIF3D-VARIANT solver was used to solve the same reactor and multigroup cross sections that was analyzed by MCNP. Since an exact analytical match is not expected, this verification focuses on ensuring that the DIF3D solution converges to the MCNP solution as the spatial and angular orders are increased.
The DIF3D "Rods In" eigenvalue errors are provided in Table IV. This table is organized similar to the  previous case with  The spatial flux order (8-3-5, 10-3-5, and 12-3-5) reduces the eigenvalue error by a few tens of pcm in some cases. The spatial leakage approximation (10-3-5, 10-4-5, and 10-5-5) is slightly closer to convergence with less than a 5 pcm change between the 10-4-5 and 10-5-5 cases. Combined, this type of behavior indicates that additional spatial flux and leakage refinement may be useful, if it could be performed. However, given that the flux order and leakage order approximations have reverse trends (increasing the flux order increases the eigenvalue while increasing the leakage order decreases the eigenvalue), it is sufficient to assume that the 10-4-5 approximation is closest to the fully spatially refined result.
The last aspect of refinement to consider is angular refinement. The two 13 th angular order cases that completed show a few pcm change compared to the equivalent P11 cases. Therefore it is possible that increasing the angular order to 15 and above may be worth another several pcm. Overall, these results are sufficiently accurate given the inability to further refine the space-angle approximations.
The peak flux and flux root-mean-squared-errors followed similar trends to the eigenvalue errors shown in Table IV. Therefore, in lieu of repeating these, we will simply examine the 10-4-5 P7 case because, as stated, this case contains offsetting errors from the flux and leakage spatial approximations. For this case, -0.05%, -0.03%, -0.02%, and -0.03% differences between the peak fluxes in the DIF3D and MCNP cases were observed; these peak fluxes were verified to exist in the same, or symmetrically equivalent locations. The overall root-mean-squared-error is 1.1x10 -6 , 3.1x10 -7 , 9.2x10 -8 , and 7.0x10 -8 for groups 1-4.
The region-wise flux comparisons for DIF3D and MCNP for the 10-4-5 P7 case are provided in Figure 7. This figure provides two plots, each showing the ratio of the DIF3D-to-MCNP group 1 fluxes. The left figure shows this ratio in the bottom axial plane while the right figure provides it for the central plane. The groups 2, 3, and 4 behaviors exhibit the same general trends, though smaller in magnitude. These results show overall very close agreement for the majority of the regions, with the maximal errors being observed in the outer ring next to the vacuum boundary. These outer ring locations are typically where DIF3D is known to perform the weakest due to the difficulty in representing vacuum boundary conditions with spherical harmonics. The MCNP relative error for these outer ring locations was between 0.1 and 0.3% for the bottom-most plane, and 0.03 to 0.07% for the central axial plane.
Ratios in the Bottommost Axial Plane Ratios in the Central Axial Plane Figure 7. Ratio of DIF3D to MCNP Fluxes in Group 1 These results indicate that the DIF3D can closely match the MCNP reference, converging with increasing spatial and angular orders, and therefore can be used to provide confidence in the verification of DIF3D.

CONCLUSIONS
This paper provides an overview of the verification tasks performed to verify the forward and adjoint forms of the fixed source, inhomogeneous fixed source, and k-eigenvalue forms of the transport and diffusion equation as implemented specifically for 3D hexagonal geometries in DIF3D. Since analytic solutions of the neutron diffusion and transport equations are either limited in scope or not possible , this verification required multiple tiers of problems unique to each solver and geometry type, each testing features independent and complementary arguments for why this separate testing of functionalities is acceptable. This separate testing was also supplemented with a high-level integral check of each the diffusion and transport capabilities and applicable geometries. The results of this verification work were provided primarily for the forward eigenvalue cases and showed how the tiers of problems were used to provide the verification basis for the software. The cases not discussed herein will be the subject of future publications.