UFDC Home  myUFDC Home  Help 



Full Text  
ADVANCED QUADRATURE SETS, ACCELERATION AND PRECONDITIONING TECHNIQUES FOR THE DISCRETE ORDINATES METHOD IN PARALLEL COMPUTING ENVIRONMENTS By GIANLUCA LONGONI A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2004 Copyright 2004 by GIANLUCA LONGONI I dedicate this research work to Rossana and I thank her for the support and affection demonstrated to me during these years in college. This work is dedicated also to my family, and especially to my father Giancarlo, who always shared my dreams and encouraged me in pursuing them. "Only he who can see the invisible, can do the impossible." By Frank Gaines ACKNOWLEDGMENTS The accomplishments achieved in this research work would have not been possible without the guidance of a mentor such as Prof. Alireza Haghighat; he has been my inspiration for achieving what nobody else has done before in the radiation transport area. I wish to thank Prof Glenn E. Sjoden for his endless help and moral support in my formation as a scientist. I also express my gratitude to Dr. Alan D. George, for providing me with the excellent computational facility at the HighPerformance Computing and Simulation Research Lab. I am thankful to Prof. Edward T. Dugan for his useful suggestions and comments, as well as to Dr. Ray G. Gamino from Lockheed Martin  KAPL and Dr. Joseph Glover, for being part of my Ph.D. committee. I also thank UFTTG for the interesting conversations regarding radiation transport physics and computer science. I am also grateful to the U.S. DOE Nuclear Education Engineering Research (NEER) program, the College of Engineering, and the Nuclear and Radiological Engineering Department at the University of Florida, for supporting the development of this research work. TABLE OF CONTENTS page A C K N O W L E D G M E N T S ................................................................................................. iv LIST OF TABLES ....................................................... ............ ....... ....... ix LIST OF FIGURES ......... ........................................... ............ xi ABSTRACT ........ .............. ............. ...... ...................... xvi CHAPTER 1 IN TR OD U CTION ............................................... .. ......................... .. 1.1 O v erv iew ............... ..................................................... ................ 1 1.2 The Linear B oltzm ann Equation................................ ....................................1 1.3 Advanced Angular Quadrature Sets for the Discrete Ordinates Method...........2 1.4 Advanced Acceleration Algorithms for the SN Method on Parallel Computing E n v iro n m en ts .................................................................... .. 5 1.5 The EvenParity Sim plified SN M ethod .............................................................6 1.6 A New Synthetic Acceleration Algorithm Based on the EPSSN Method ......10 1.7 An Automatic Preconditioning Algorithm for the SN Method: FAST (Flux Acceleration Simplified Transport) ...................................... ............... 12 1.8 Outline ............... ...... ... .......... ...................................13 2 THE DISCRETE ORDINATES METHOD..................................................14 2.1 D iscrete Ordinates M ethod (SN) ........................................... ............... 14 2.1.1 Discretization of the Energy Variable............... .... ...............14 2.1.2 Discretization of the Angular Variable............................................ 17 2.1.3 Discretization of the Spatial Variable............. ...................20 2.1.4 D ifferencing Schem es ...................................... ....................... ......... 21 2.1.4.1 LinearDiamond Scheme (LD).......................... ..............22 2.1.4.2 Directional ThetaWeighted Scheme (DTW) ............................23 2.1.4.3 Exponential DirectionalWeighted Scheme (EDW) ...................24 2.1.5 The Flux M om ents ................ ........... .....................................25 2.1.6 B oundary C conditions ........................................ ......................... 25 2.2 Source Iteration M ethod ............................................................................ 25 2.3 Pow er Iteration M ethod .................................... .... ........................ 26 2.4 Acceleration Algorithms for the SN M ethod.............................................. 27 3 ADVANCED QUADRATURE SETS FOR THE SN METHOD ...........................29 3.1 Legendre EqualWeight (PNEW) Quadrature Set .......................................30 3.2 LegendreChebyshev (PNTN) Quadrature Set..............................................31 3.3 The Regional Angular Refinement (RAR) Technique .................................33 3.4 Analysis of the Accuracy of the PNEW and PNTN Quadrature Sets..............34 3.5 Testing the Effectiveness of the New Quadrature Sets.................................38 3.5.1 K obayashi Benchm ark Problem 3 ..................................... ................... 38 3.5.2 CTScan Device for Medical/Industrial Imaging Applications.............43 4 DERIVATION OF THE EVENPARITY SIMPLIFIED SN EQUATIONS ............47 4.1 Derivation of the Simplified Spherical Harmonics (SPN) Equations..............48 4.2 Derivation of the EvenParity Simplified SN (EPSSN) Equations..................51 4.2.1 Boundary Conditions for the EPSSN Equations ..................................55 4.2.2 Fourier Analysis of the EPSSN Equations .......................................... 56 4.2.3 A New Formulation of the EPSSN Equations for Improving the Convergence Rate of the Source Iteration Method..............................59 4.3 Comparison of the Pi Spherical Harmonics and SP1 Equations...................60 5 NUMERICAL METHODS FOR SOLVING THE EPSSN EQUATIONS ...............65 5.1 Discretization of the EPSSN Equations Using the FiniteVolume Method ....65 5.2 Numerical Treatment of the Boundary Conditions........................................72 5.3 The Compressed Diagonal Storage Method ......................................... 74 5.4 Coarse M esh Interface Projection Algorithm ............. .................................. 75 5.5 Krylov Subspace Iterative Solvers..................... ............ ...... ............. 80 5.5.1 The Conjugate Gradient (CG) M ethod ..................................................82 5.5.2 The BiConjugate Gradient M ethod ............... ................... ............83 5.5.3 Preconditioners for Krylov Subspace M ethods ......................................84 6 DEVELOPMENT AND BENCHMARKING OF THE PENSSn CODE..................86 6.1 Development of the PENSSn (Parallel Environment Neutralparticle Sim plified Sn) Code..................... ........ ..................................... 87 6.2 Numerical Analysis of Krylov Subspace Methods........................................92 6.2.1 Coarse M esh Partitioning of the M odel................................................92 6.2.2 B oundary C conditions ........................................ ........................ 95 6.2.3 M material H eterogeneities....................... .................. .....................96 6.2.4 Convergence Behavior of Higher EPSSN Order Methods...................97 6.3 Testing the Incomplete Cholesky Conjugate Gradient (ICCG) Algorithm .....99 6.4 Testing the Accuracy of the EPSSN M ethod .............. ....... ....................100 6 .4 .1 S catterin g R atio .......................................................... ..................... 10 0 6.4.2 Spatial Truncation Error ............................................ ............... 103 6.4.3 Low D ensity M aterials...................................... ........................ 104 6.4.4 M material Discontinuities................... ..... ........................ 108 6.4.5 Anisotropic Scattering .............. ... ............. ............ ............... 111 6.4.6 Small Light Water Reactor (LWR) Criticality Benchmark Problem.... 115 6.4.7 Small Fast Breeder Reactor (FBR) Criticality Benchmark Problem.... 120 6.4.8 The MOX 2D Fuel Assembly Benchmark Problem.........................124 7 PARALLEL ALGORITHMS FOR SOLVING THE EPSSN EQUATIONS ON DISTRIBUTED MEMORY ARCHITECTURES ............................................. 128 7.1 Parallel Algorithm s for the PEN SSn Code...................................................128 7.2 D om ain D ecom position Strategies ..................................... ............... ..130 7.2.1 Angular D om ain D ecom position........................................................ 130 7.2.2 Spatial D om ain D ecom position............. ...................... ... ...............131 7.2.3 Hybrid Domain Decomposition ............. .................. ....................131 7.3 Parallel Performance of the PENSSn Code .................................. .............. 132 7.4 Parallel Performance of PENSSn Applied to the MOX 2D Fuel Assembly B enchm ark Problem ........................................................... ............... 139 8 DEVELOPMENT OF A NEW SYNTHETIC ACCELERATION METHOD BASED ON THE EPSSN EQUATION S ........................................ ......................... 140 8.1 The EPSSN Synthetic Acceleration M ethod ................................................141 8.2 Spectral Analysis of the EPSSN Synthetic Acceleration Method...............145 8.3 Analysis of the Algorithm Stability Based on Spatial Mesh Size ...............148 8.3.1 Comparison of the EPSSN Synthetic Acceleration with the Simplified A ngular M ultigrid M ethod................................................ ............... 150 8.4 Limitations of the EPSSN Synthetic Acceleration Method ..........................153 9 FAST: FLUX ACCELERATION SIMPLIFIED TRANSPORT PRECONDITIONER BASED ON THE EPSSN METHOD .............................154 9.1 Development and Implementation of FAST .............................................154 9.2 Testing the Performance of the FASTc Preconditioning Algorithm........... 157 9.2.1 Criticality Eigenvalue Problem ................................. ............... 157 9.2.2 Fixed Source Problem ....................................................................... 159 9.3 The MOX 3D Fuel Assembly Benchmark Problem................................161 9.3.1 MOX 3D Unrodded Configuration............................... ............... 162 9.3.2 MOX 3D RoddedA Configuration................... .................165 9.3.3 M OX 3D RoddedB Configuration....... ........................................ 167 10 SUMMARY, CONCLUSION, AND FUTURE WORK ........................................171 APPENDIX A EXPANSION OF THE SCATTERING TERM IN SPHERICAL HARMONICS..175 B PERFORMANCE OF THE NEW EPSSN FORMULATION ...............................177 L IST O F R E F E R E N C E S ......... .. ............... ................. .............................................. 180 BIOGRAPH ICAL SKETCH .............................................................. ............... 185 LIST OF TABLES Table pge 31. Evenmoments obtained with a PNEW S30 quadrature set. .....................................35 32. Evenmoments obtained with a PNTN S30 quadrature set .......................................36 33. CPU time and total number of directions required for the CTScan simulation........45 61. Comparison of number of iterations required to converge for the CG and BiCG algorithm s ......... ... ......... .. ............................... ........ ..... .. ............ 93 62. Number of Krylov iterations required to converge for the CG and BiCG algorithms with different boundary conditions. .............................................. ............... 95 63. Number of Krylov iterations required to converge for the CG and BiCG algorithms for heterogeneous the box in a box problem. ........................... ................................ 96 64. Number of Krylov iterations required to converge for the CG and BiCG algorithms for the EPSSs equations. ............................................... .............................. 97 65. Number of iterations for the ICCG and CG algorithms. ........................ ..........99 66. Two groups crosssections and fission spectrum. ..................................................106 67. Comparison of keff obtained with the EPSSN method using DFM versus PENTRAN* S6 (Note that DFM=1.0 implies no crosssections scaling). ................................106 68. Balance tables for the EPSSN and S16 methods and relative differences versus the S 16 so lu tio n ........................................... ................... ....... ........ 1 1 1 69. Integral boundary leakage for the EPSSN and S16 methods and relative differences versus the S16 solution. .................. .......... .. ........ .. ............ .............. 111 610. Fixed source energy spectrum and energy range .................................................112 611. Maximum and minimum relative differences versus the S8 method for energy group 1 and 2....................................................................... ...... ....... 113 612. Twogroup crosssections for the small LWR problem .............. ... .................116 612. Twogroup crosssections for the small LWR problem (Continued).....................116 613. Fission spectrum and energy range for the small LWR problem...........................117 614. Criticality eigenvalues calculated with different EPSSN orders and relative error compared to M onte Carlo predictions ....................................... ...... ...........117 615. CRWs estimated with the EPSSN method ......................................................118 616. Criticality eigenvalues for the small FBR model. .............................................122 617. CRWs estimated with the EPSSN and Monte Carlo methods............................122 618. Criticality eigenvalues and relative errors for the MOX 2D benchmark problem.125 71. Data relative to the load imbalance generated by the Krylov solver........................136 72. Parallel performance data obtained on PCPENII Cluster............... ... .................138 73. Parallel performance data obtained on Kappa Cluster. ...........................................138 74. Parallel performance data for the 2D MOX Fuel Assembly Benchmark problem (P C PE N II C luster)............ .......................................................... ................... 139 81. Spectral radius for the different iterative methods. .............................................147 82. Comparison of the number of inner iteration between EPSSN synthetic methods and unaccelerated transport ........................................................................... ..... .... 149 91. Criticality eigenvalues obtained with the preconditioned PENTRANSSn code for different E P SSN orders................................................. ............................. 159 92. Results obtained for the MOX 3D in the Unrodded configuration.........................162 93. Results obtained for the MOX 3D RoddedA configuration. ...............................165 94. Results obtained for the MOX 3D RoddedB configuration................................167 B1. Performance data for the standard EPSSN formulation.......................................177 B2. Performance data for the new EPSSN formulation..........................................177 B3. Ratio between Krylov iterations and inner iterations. ...........................................178 B4. Inner iterations and time ratios for different SSN orders............... .... ..............178 LIST OF FIGURES Figure pge 21. Cartesian spaceangle coordinates system in 3D geometry. ....................................16 22. Point weight arrangement for a S8 levelsymmetric quadrature set .........................19 23. S20 LQN quadrature set. ......................................................................20 31. S28 PNEW quadrature set. ............................................... .. .. ................................ 31 32. S28 PNTN quadrature set ................................................... .............................. 33 33. PNTN quadrature set (S16) refined with the RAR technique. ....................................34 34. Configuration of the test problem for the validation of the quadrature sets ..............37 35. Relative difference between the currents Jx and Jz for the test problem.....................37 3.6. Mesh distribution for the Kobayashi benchmark problem 3: A) Variable mesh distribution; B) Uniform mesh distribution.............. ............................................39 37. Comparison of S20 quadrature sets in zone 1 at x=5.0 cm and z=5.0 cm. .................40 38. Comparison of S20 quadrature sets in zone 2 aty=55.0 cm and z=5.0 cm. ................40 39. Comparison of PNEW quadrature sets for different SN orders in zone 1 at x=5.0 cm an d z= 5 .0 cm ...................................................... ................. 4 1 310. Comparison of PNEW quadrature sets for different SN orders in zone 2 aty=55.0 cm and z=5.0 cm ............. .... .......................... ...... .. ........... 41 311. Comparison of PNTN quadrature sets for different SN orders in zone 1 aty=5.0 cm and z= 5.0 cm ...................................................... ................. 42 312. Comparison of PNTN quadrature sets for different SN orders in zone 2 aty=55.0 cm and z= 5.0 cm ...................................................... ................. 42 313. Crosssectional view of the CTScan model on the xy plane.............................43 314. Scalar flux distribution on the xy plane obtained with an S20 levelsymmetric quadrature set. ........................................................................44 315. Scalar flux distribution on the xy plane obtained with an S50 PNTN quadrature set.44 316. Scalar flux distribution on the xy plane obtained with an S30 PNTN quadrature set biased with RAR. ............................. ..... ... ....... ....... .............. .. 45 3.17. Comparison of the scalar flux at detector position (x=72.0 cm).............................46 5.1. Fine mesh representation on a 3D Cartesian grid ................ ............. ...............67 5.2. View of a fine mesh along the xaxis ........ .. .... ............................. .... ........... 68 5.3. Representation of a coarse mesh interface................................. ......... ............. 76 5.4. Representation of the interface projection algorithm between two coarse meshes. ...79 61. D description of PEN SSn input file .............. ................................... ............... 88 62. Flowchart of the PEN SSn code ...................................................... ..................90 62. Flowchart of the PENSSn code (Continued)..........................................................91 63. Configuration of the 3D test problem ............................................. ............... 92 64. Convergence behavior of the CG algorithm for the nonpartitioned model ..............94 65. Heterogeneous configuration for the 3D test problem...................... ...............96 66. Distribution of eigenvalues for the EPSSs equations.........................................98 67. Configuration of the 2D criticality eigenvalue problem. ........................................101 68. Criticality eigenvalues as a function of the scattering ratio (c) for different methods. 101 69. Relative difference for criticality eigenvalues obtained with different EPSSN methods compared to the S16 solution (PENTRAN code). ...................................102 610. Plot of criticality eigenvalues for different mesh sizes........................................103 611. Plot of the relative difference of the EPSSN solutions versus transport S16 for different m esh sizes .................. ............................. .. .... ..... ............ 104 612. Uranium assembly test problem view on the xy plane ........................................105 613. Relative difference of physical quantities of interest calculated with EPSSN method compared to the S6 PENTRAN solution. .................................... .................107 614. Convergence behavior of the PENSSn with DFM=100.0 and PENTRAN S6 ...... 108 615. Geometric and material configuration for the 2D test problem ...........................109 616. Scalar flux distribution at material interface (y=4.84 cm)...................................109 617. Relative difference versus S16 calculations at material interface (y=4.84 cm).......110 618. Fraction of scalar flux values within different ranges of relative difference (R.D.) in en erg y g rou p 1 ................................................................... 1 12 619. Fraction of scalar flux values within different ranges of relative difference (R.D.) in en ergy g rou p 2 ....................................................................................... 1 13 620. Front view of the relative difference between the scalar fluxes obtained with the EPSSs and S8 methods in energy group 1................... ........................... 114 621. Rear view of the relative difference between the scalar fluxes obtained with the EP SS8 and S8 m ethods in energy group 1....................................... ...............115 622. Model view on the xy plane. A) view of the model from z=0.0 cm to 15.0 cm, B) view of the model from z=15.0 cm to z=25.0 cm .................................................. 115 623. M odel view on the xz plane .................................................... .................. 116 624. Normalized scalar flux for case 1, in group 1 along the xaxis at y=2.5 cm and z=7.5 c m ................... ........................................................ ................ 1 1 8 625. Scalar flux distributions. A) Case 1 energy group 1, B) Case 2 energy group 1, C) Case 1 energy group 2, D) Case 2 energy group 2 ..............................................119 626. View on the xy plane of the small FBR model................... ...................................120 627. View on the xz plane of the small FBR model................... ...................................121 628. Scalar flux distribution in energy group 1: A) Case 1; B) Case 2..........................123 629. Scalar flux distribution in energy group 4: A) Case 1; B) Case 2..........................123 630. Mesh distribution of the MOX 2D Fuel Assembly Benchmark problem. ............124 631. Scalar flux distribution for the 2D MOX Fuel Assembly benchmark problem (EP SS4): A) Energy group 1; B) Energy group 2; C) Energy group 3; D) Energy group 4; E) Energy group 5; F) Energy group 6; G) Energy group 7. ..........................126 632. Normalized pin power distribution for the 2D MOX Fuel Assembly benchmark problem (EPSS4): A) 2D view; B) 3D view..................................................127 71. Hybrid decomposition for an EPSS6 calculation (3 directions) for a system partitioned with 4 coarse meshes on 6 processors............ ................................. 132 72. Speedup obtained by running PENSSn on the Kappa and PCPENII Clusters........134 73. Parallel efficiency obtained by running PENSSn on the Kappa and PCPENII C lu ste rs ...................................... ................................................. 1 3 5 74. Angular domain decomposition based on the automatic load balancing algorithm. 137 75. Parallel fraction obtained with the PENSSn code. ................................................137 81. Spectrum of eigenvalues for the Source Iteration and Synthetic Methods based on different SSN orders............. .... ........................................................ .... .... .... 147 82. Number of inner iterations required by each acceleration method as a function of the m esh size. ..........................................................................149 83. Number of inner iterations as a function of the scattering ratio (DZ differencing scheme e)....................................................... ................... .. .... ... ... . 150 84. Number of inner iterations as a function of the scattering ratio (DTW differencing scheme e)....................................................... ................... .. .... .. ... . 15 1 85. Number of inner iterations as a function of the scattering ratio (EDW differencing scheme e)....................................................... ................... .. .... .. ... . 15 1 86. Number of inner iterations for the EPSS2 synthetic method obtained with DZ, DTW, and EDW differencing schemes. .................................................. ......... ...... 152 91. Card required in PENTRANSSn input deck to initiate SSN preconditioning.........155 92. Flowchart of the PENTRANSSn Code System............... ...... .................156 93. Ratio of total number of inner iterations required to solve the problem with preconditioned PENTRANSSn and nonpreconditioned PENTRAN .................158 94. Relative change in flux value in group 1 ...................................... ............... 160 95. Relative change in flux value in group 2.............. ......................... ...... ............ 160 96. Behavior of the criticality eigenvalue as a function of the outer iterations............163 97. Convergence behavior of the criticality eigenvalue. ..............................................164 98. Preconditioning and transport calculation phases with relative computation time requ ired ........................................................ ............................. 164 99. Behavior of the criticality eigenvalue as a function of the outer iterations..............165 910. Convergence behavior of the criticality eigenvalue. ............................................166 911. Preconditioning and transport calculation phases with relative computation time requ ired ........................................................ ............................. 167 912. Behavior of the criticality eigenvalue as a function of the outer iterations ..........168 913. Convergence behavior of the criticality eigenvalue. ............................................169 914. Preconditioning and transport calculation phases with relative computation time re q u ire d ......... .............. ...................... .................................................16 9 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy ADVANCED QUADRATURE SETS, ACCELERATION AND PRECONDITIONING TECHNIQUES FOR THE DISCRETE ORDINATES METHOD IN PARALLEL COMPUTING ENVIRONMENTS By Gianluca Longoni December 2004 Chair: Alireza Haghighat Major Department: Nuclear and Radiological Engineering In the nuclear science and engineering field, radiation transport calculations play a keyrole in the design and optimization of nuclear devices. The linear Boltzmann equation describes the angular, energy and spatial variations of the particle or radiation distribution. The discrete ordinates method (SN) is the most widely used technique for solving the linear Boltzmann equation. However, for realistic problems, the memory and computing time require the use of supercomputers. This research is devoted to the development of new formulations for the SN method, especially for highly angular dependent problems, in parallel environments. The present research work addresses two main issues affecting the accuracy and performance of SN transport theory methods: quadrature sets and acceleration techniques. New advanced quadrature techniques which allow for large numbers of angles with a capability for local angular refinement have been developed. These techniques have been integrated into the 3D SN PENTRAN (Parallel Environment Neutralparticle TRANsport) code and applied to highly angular dependent problems, such as CTScan devices, that are widely used to obtain detailed 3D images for industrial/medical applications. In addition, the accurate simulation of core physics and shielding problems with strong heterogeneities and transport effects requires the numerical solution of the transport equation. In general, the convergence rate of the solution methods for the transport equation is reduced for large problems with optically thick regions and scattering ratios approaching unity. To remedy this situation, new acceleration algorithms based on the EvenParity Simplified SN (EPSSN) method have been developed. A new standalone code system, PENSSn (Parallel Environment Neutralparticle Simplified Sn), has been developed based on the EPSSN method. The code is designed for parallel computing environments with spatial, angular and hybrid (spatial/angular) domain decomposition strategies. The accuracy and performance of PENSSn has been tested for both criticality eigenvalue and fixed source problems. PENSSn has been coupled as a preconditioner and accelerator for the SN method using the PENTRAN code. This work has culminated in the development of the Flux Acceleration Simplified Transport (FASTc) preconditioning algorithm, which constitutes a completely automated system for preconditioning radiation transport calculations in parallel computing environments. CHAPTER 1 INTRODUCTION 1.1 Overview In the nuclear engineering field, particle transport calculations play a keyrole in the design and optimization of nuclear devices. The Linear Boltzmann Equation (LBE) is used to describe the angular, energy and spatial variations of the particle distribution, i.e., the particle angular flux.1 Due to the integrodifferential nature of this equation, an analytical solution cannot be obtained, except for very simple problems. For real applications, the LBE must be solved numerically via an iterative process. To solve large, realworld problems, significant memory and computational requirements can be handled using parallel computing environments, enabling memory partitioning and multitasking. The objective of this dissertation is to investigate new techniques for improving the efficiency of the of the SN method for solving problems with highly angular dependent sources and fluxes in parallel environments. In order to achieve this goal, I have investigated two major areas: 1. Advanced quadrature sets for problems characterized with highly angular dependent fluxes and/or sources. 2. Advanced acceleration/preconditioning algorithms. 1.2 The Linear Boltzmann Equation The LBE is an integropartial differential equation, which describes the behavior of neutral particle transport. The Boltzmann equation, together with the appropriate boundary conditions, constitutes a mathematically wellposed problem having a unique solution. The solution is the distribution of particles throughout the phase space, i.e., space, energy, and angle. The distribution of particles is, in general, a function of seven independent variables: three spatial, two angular, one energy, and one time variable. The timeindependent LBE in its general integrodifferential form1 is given by Eq. 1.1. Q. V 7 (F, E, Q) + o, (F, E)VC(F, E, Q) = qe (r, E, Q) + JdE' JQ'o, (F,' E' ,Q'Q)(FE',Q') (1.1) 0 4;T + XE dE' JdQ r(F, E')(F, E', Q'). 47r 0 4;T In Eq. 1.1, I have defined the following quantities: V.: Angular flux [particles/cm2/MeV/sterad/sec] . F: Particle position in a 3D space [cm]. E: Particle energy [MeV]. 0: Particle direction unit vector. a,: Macroscopic total crosssection [1/cm/MeV]. qe,,: External independent source [particles/cm3/MeV/sterad/sec] . ao: Macroscopic doubledifferential scattering crosssection [1/cm/sterad/MeV]. S: Fission spectrum [1/MeV]. v : Average number of neutrons generated per fission. 7f : Macroscopic fission crosssection [1/cm/MeV]. 1.3 Advanced Angular Quadrature Sets for the Discrete Ordinates Method The discrete ordinates method (SN) is widely used in nuclear engineering to obtain a numerical solution of the integrodifferential form of the Boltzmann transport equation. The method solves the LBE for a discrete number of ordinates (directions) with associated weights.2 The combination of discrete directions and weights is referred to as quadrature set.3 The major drawback of the SN method is the limited number of directions involved, which, in certain situations, may lead to the so called rayeffects, which appear as unphysical flux oscillations. In general, this behavior occurs for problems with highly angular dependent sources and/or fluxes, or when the source is localized in a small region of space, in low density or highly absorbent media. In the past, several remedies for rayeffects have been proposed; the most obvious one is to increase the number of directions of the quadrature set, or equivalently the SN order. However, this approach can lead to significant memory requirements and longer computational times. Carlson and Lathrop proposed a number of remedies45 for ray effects based on specialized quadrature sets, which satisfy higher order moments of the direction cosines. Remedies based on firstcollision approximations have also been investigated.3 The source of particles generated from firstcollision processes is often less localized than the true source; hence, the flux due to this source is usually less prone to rayeffects than the flux from the original source. If the true source is simple enough, analytic expressions can be obtained for the uncollided flux and used to produce a first collision source; however, for general sources and deeppenetration problems, this method can be very timeconsuming. An alternative approach6 is to expand the angular flux in terms of spherical harmonics (PN). The PN method does not suffer from ray effects, because the angular dependency in the angular flux is treated using continuous polynomial functions. However, the PN method has found limited applicability due to its computational intensive requirements. One of the most widely used techniques for generating quadrature sets is the level symmetric3 (LQN) method; however, the LQN method yields negative weights beyond order S20. In problems with large regions of air or highly absorbent materials, higher order (>20) quadrature sets are needed in order to mitigate rayeffects; therefore, it is necessary to develop techniques which allow for higher order or biased quadrature sets. In the past, different techniques have been investigated to remedy this issue. The equal weight quadrature set (EQN), developed by Carlson,7 is characterized by positive weights for any SN order. Other quadrature sets have been derived, by relaxing the constraints imposed by the LQN method;2 for this purpose, the Gauss quadrature formula and Chebyshev polynomials have been used for onedimensional cylindrical or two dimensional rectangular geometries.4 In a recent study, the uniform positiveweight quadrature sets8 (UEN and UGN) have been derived. The UEN quadrature set is derived by uniformly partitioning the unit sphere into the number of directions defined by the SN order while the UGN quadrature set selects the ordinates along the zaxis as roots of Legendre polynomials. A new biasing technique, named Ordinate Splitting (OS), has been developed9 for the Equal Weight (EW) quadrature set; the OS technique is a method which is suitable to solve problems in which the particle angular flux and/or source are peaked along certain directions of the unit sphere. The idea is to select a direction of flight of the neutron and split it into a certain number of directions of equal weights, while conserving the original weight. This new biasing technique has been implemented in the PENTRAN code10 and it has been proven very effective for medical physics applications such as CTScan devices. 1115 In this research work, I have developed new quadrature sets1112 based on Legendre (PN) and Chebyshev (TN) polynomials. The LegendreChebyshev (PNTN) quadrature set is derived by setting the polar angles equal to the roots of the Legendre polynomial of order N, and the azimuthal angles are calculated by finding the roots of the Chebyshev polynomials. The Legendre EqualWeight (PNEW) quadrature set is derived by choosing the polar angles as the roots of the Legendre polynomial of order N, while the azimuthal angles are calculated by equally partitioning a 90 degree angle. The set of directions is then arranged on the octant of the unit sphere similar to the levelsymmetric triangular pattern. The main advantage of these new quadrature sets is the absence of negative weights for any SN order, and their superior accuracy compared to other positive schemes such as equal weight quadrature sets. Moreover, I have developed a new refinement technique, referred to as Regional Angular Refinement (RAR), which leads to a biased angular quadrature set.1315 The RAR technique consists of fitting an auxiliary quadrature set in a region of the unit sphere, where refinement is needed. These quadrature sets have been applied successfully to large problems, such as a CTScan device used for medical/industrial imaging9 and a TimeofFlight (TOF) experiment simulation.16 The benefit of using biased quadrature sets is to achieve accurate solutions for highly angular dependent problems with reduced computational time. 1.4 Advanced Acceleration Algorithms for the SN Method on Parallel Computing Environments Radiation transport calculations for realistic systems involve the solution of the discretized SN equations. A typical 3D transport model requires the discretization of the SN equations in 300,000 spatial meshes, 47 energy groups, and considering an Ss calculation, a total of 80 directions on the unit sphere. These figures yield approximately 1.13 billion unknowns. In terms of computer memory, this number translates into 9 GBytes of RAM for storage (in single precision) of only the angular fluxes. It is clear that an efficient solution for such a problem is out of the scope of a regular workstation available with current technology. Therefore, it is necessary to develop new algorithms capable of harnessing the computational capabilities of supercomputers. Based on this philosophy, in the late 1990s G. Sjoden and A. Haghighat have developed a new 3D parallel radiation transport code: PENTRAN (Parallel Environment Neutralparticle Transport).10 However, besides the size and complexity of the problem being solved, other aspects come into play, especially when dealing with criticality eigenvalue problems. Because of the physics of these problems, the convergence rate of the currently used iterative methods is quite poor. For realistic problems, such as wholecore reactor calculations performed in a 3D geometry, the solution of the SN equations may become impractical if proper acceleration methods17 are not employed. The main philosophy behind the novel acceleration algorithms developed in this work is to employ a simplified mathematical model which closely approximates the SN equations, yet can be solved efficiently. The new acceleration/preconditioning algorithms have been developed during the course of this research in three major phases: 1. Development of the PENSSn code based on the EvenParity Simplified SN (EP SSN) method. 2. Investigation of a new synthetic acceleration algorithm based on the EPSSN method. 3. Development of an automated acceleration system for the SN method on parallel environments: FAST (Flux Acceleration Simplified Transport). 1.5 The EvenParity Simplified SN Method The EvenParity Simplified SN (EPSSN) method is developed based on the philosophy considered in the PN and Simplified Spherical Harmonics (SPN) methods.18 The spherical harmonics (PN) approximation to the transport equation is obtained by expanding the angular flux using spherical harmonics functions truncated to order N, where Nis an odd number; these functions form a complete basis in the limit of the truncation error. In the limit of N oc, the exact transport solution is recovered.1 In 3D geometries, the number of PN equations grows as (N + 1)2. The PN equations can be reformulated in a secondorder form, cast as (N + 1)2 / 2 diffusionlike equations, characterized by an elliptic operator. However, the number of these equations is very large, and the coupling involves not only angular moments, but also mixed spatial derivatives of these moments.6 Because of these issues, to reduce the computing time in the early 1960s, Gelbard et al. proposed the Simplified Spherical Harmonics18 or SPN method. The Gelbard procedure consists of extending the spatial variable to 3D by substituting the second order derivatives in the 1D PN equations with the 3D Laplacian operator. As a result, coupling of spatial derivatives is eliminated, yielding only (N+ 1) equations as compared to (N + 1)2. Further, since the SPN equations can be reformulated as secondorder elliptic equations, effective iterative techniques such as Krylov subspace1920 methods, can be employed. The main disadvantage of the SPN equations is that the exact solution to the transport equation is not recovered as N oc, due to terms that are inherently omitted in replacing a 1D leakage term with a simplified 3D formulation. However, for idealized systems characterized by homogeneous materials and isotropic scattering, the SPN and the SN equations yield the same solution, given the same quadrature set and spatial discretization is used for both methods. Despite this fact, the SPN equations yield improved solutions2122 compared to the currently used diffusion equation. The theoretical basis for the SPN equations has been provided by many authors using variational methods and asymptotic analysis.22 It has been shown that these equations are higherorder asymptotic solutions of the transport equation. Moreover, Pomraning has demonstrated that the SPN equations, for odd N, are a variational approximation to the onegroup evenparity transport equation with isotropic scattering in an infinite homogeneous medium. Recently, the SPN formulation has received renewed interest, especially in reactor physics applications. For applications such as the MOX fuel assemblies2223 or for reactor problems with strong transport effects,24 diffusion theory does not provide accurate results, while the SPN equations improve the accuracy of the solution within a reasonable amount of computation time. Initially, I derived the SP3 equations starting from the 1D P3 equations;25 however, for developing a general order algorithm, I derived a general formulation using the even parity form of the SN transport equation.26 The EvenParity Simplified SN (EPSSN) formulation has some interesting properties that make it suitable to develop algorithms of any arbitrary order. Chapter 4 is dedicated to this issue. To make the method more effective, the convergence properties of the EPSSN equations were improved by modifying the scattering term; it will be shown that this improved derivation is problem dependent but can reduce the number of iterations significantly. I have developed a general 3D parallel code,23 PENSSn (Parallel Environment Neutralparticle Simplified SN), based on the EPSSN equations. The EPSSN equations are discretized with a finitevolume approach, and the spatial domain is partitioned into coarse meshes with variable fine grid density in each coarse mesh emulating the PENTRAN grid structure. A projection algorithm is used to couple the coarse meshes based on the values of the evenparity angular fluxes on the interfaces. PENSSn uses iterative solvers based on Krylov subspace19 methods: the Conjugate Gradient (CG) and the BiConjugate Gradient (BiCG) solvers. In addition, I have developed a preconditioner based on the Incomplete Cholesky factorization20 for the CG method. The finitevolume discretization of the EPSSN equations in a 3D Cartesian geometry yields sparse matrices with a 7diagonal sparse structure. Therefore, I optimized the memory management of PENSSn by using a Compressed Diagonal Storage (CDS) approach, where only the nonzero elements on the diagonals are stored in memory. The PENSSn code is designed for parallel computing environments with angular, spatial and hybrid (angular/spatial) domain decomposition algorithms.23 The space decomposition algorithm partitions the 3D Cartesian space into coarse meshes which are then distributed among the processors while the angular decomposition algorithm partitions the directions among the processors. The hybrid decomposition algorithm is a combination of the two algorithms discussed above. Note that the hybrid decomposition combines the benefits of memory partitioning offered by the spatial decomposition algorithm, and the speedup offered by the angular decomposition algorithm. The code is written in Fortran 90, and for seamless parallelization, the MPI (Message Passing Interface) library27 is used. I have tested the PENSSn code for problems characterized by strong transport effects, and have shown that the improvements over the diffusion equation can be significant.26 The solutions obtained with the EPSSN method are in good agreement with SN and Monte Carlo methods; however, the computation time is significantly lower. Hence, these results indicate that the EPSSN method is an ideal candidate for the development of an effective acceleration or preconditioning algorithm for radiation transport calculations. 1.6 A New Synthetic Acceleration Algorithm Based on the EPSSN Method As mentioned earlier in this chapter, the solution of the linear Boltzmann equation is obtained numerically via an iterative process. The most widely used technique to iteratively solve the transport equation is the Source Iteration (SI) method3 or Richardson iteration. The convergence properties of this method are related to the spectral radius of the transport operator. It can be shown that for an infinite, homogeneous medium, the spectral radius of the transport operator is dominated by the scattering ratio c, given by c = (1.2) U t where c, is the macroscopic scattering crosssection and is the macroscopic total cross section. Note that Eq. 1.2 can be obtained by Fourier analysis in an infinite homogeneous medium. The asymptotic convergence rate v. is defined as v, log(c). (1.3) Eq. 1.3 indicates that for problems with scattering ratios close to unity, the unaccelerated SI method is ineffective, because the asymptotic convergence rate tends toward zero. Hence, the use of an acceleration scheme is necessary. In the past, many acceleration techniques have been proposed25 for solving the steadystate transport equation. The synthetic methods have emerged as effective techniques to speedup the convergence of the SI iterative process.28 In the synthetic acceleration process a lowerorder approximation of the transport equation (e.g., diffusion theory) is corrected using the transport equation at each iteration. In this way the spectral radius of the accelerated transport operator is reduced with consequent speed up of the iteration process. Two categories of synthetic methods have been investigated so far, the Diffusion Synthetic Acceleration (DSA) and the Transport Synthetic Acceleration (TSA).2930 The first method has been proven to be very effective for 1D problems and for certain classes of multidimensional problems. However, recently it has been shown that for multi dimensional problems with significant material heterogeneities, the DSA method fails to converge efficiently.3132 The same behavior, along with possible divergence, has been reported also for TSA.30 I have developed and tested a new synthetic acceleration algorithm33 based on the simplified form of the EvenParity transport equations (EPSSN). I tested the EPSSN algorithm for simple 3D problems and I concluded that it is affected by instability problems. These instabilities are similar in nature to what has been reported for DSA by Warsa, Wareing, and Morel.3132 The main problems affecting the stability of the synthetic methods are material heterogeneities and the inconsistent discretization of the lowerorder operator with the transport operator, which leads to divergence if the mesh size is greater than 1 mean free path. Moreover, because the synthetic acceleration method has been implemented into the PENTRAN code, another consideration comes into play. The spatial differencing in the PENTRAN code system is based on an Adaptive Differencing Strategy1o (ADS); with this method the code automatically selects the most appropriate differencing scheme based on the physics of the problem. Hence, the discretization of the lowerorder operator should be consistent with every differencing scheme present in the code. This task is feasible if we consider only the lineardiamond (LD) differencing scheme; however, the complexity increases if we consider the Directional Theta Weighted (DTW) or the Exponential Directional Weighted (EDW) differencing schemes.34 Moreover, it has been shown that even a fully consistent discretization of the lowerorder operator does not guarantee the convergence of the synthetic method. Due to the issues discussed above, I have developed a system that utilizes the EP SSN method as preconditioner for the SN method. 1.7 An Automatic Preconditioning Algorithm for the SN Method: FAST, (Flux Acceleration Simplified Transport) The last phase of the development of an effective acceleration algorithm for the SN method has culminated in the development of the FASTc system (Flux Acceleration Simplified Transport). The FAST algorithm is based on the kernel of the PENSSn code. The main philosophy followed in the development of the system is completely antithetic to the synthetic acceleration approach. The idea is to quickly obtain a relatively accurate solution with the EPSSN method and to use it a preconditioned initial guess for the SN transport calculation. This approach has the main advantage of decoupling the two methods, hence avoiding all the stability issues discussed above. The FASTc system is currently implemented into the PENTRANSSn Code System for distributed memory environments, and it is completely automated from a user point of view. Currently, the system has successfully accelerated large 3D criticality eigenvalue problems, speedingup the calculations by a factor of 3 to 5 times, and hence reducing significantly the spectral radius. 1.8 Outline The remainder of this dissertation is organized as follows. Chapter 2 provides the theory for the discrete ordinates method. The discretization of the phase space variables in the transport equation will be discussed, along with the proper boundary conditions. Chapter 3 discusses the theoretical development of the advanced and biased quadrature sets for the discrete ordinates method. It also presents the application of the new quadrature sets for the simulation of a CTScan device and for the Kobayashi benchmark problem 3. Chapter 4 discusses the derivation of the EPSSN equations. Chapter 5 discusses the numerical methods for the solution of the EPSSN equations, along with the iterative solvers based on Krylov subspace methods. Chapter 6 addresses the numerics and accuracy of the EPSSN equations. Chapter 7 presents the parallel algorithms implemented in the PENSSn code for distributed memory architectures. Chapter 8 describes the development of a new synthetic acceleration algorithm based on the EPSSN method and its limitations. Chapter 9 focuses on the development of the FASTC preconditioner; the performance of the algorithm is measured with two test problems and a large 3D wholecore criticality eigenvalue calculation. Chapter 10 will draw the conclusions on the objectives accomplished and it will point out some aspects for future development. CHAPTER 2 THE DISCRETE ORDINATES METHOD In this chapter, the Discrete Ordinates Method (SN) will be discussed in detail. The discretized form of the transport equation is formulated in a 3D Cartesian geometry. I also address the iterative techniques and acceleration methods used to solve the SN transport equations. 2.1 Discrete Ordinates Method (SN) The Discrete Ordinates Method (SN) is widely used to obtain numerical solutions of the linear Boltzmann equation. In the SN method, all of the independent variables (space, energy and angle) are discretized as discussed below. 2.1.1 Discretization of the Energy Variable The energy variable of the transport equation is discretized using the multigroup approach.3 The energy domain is partitioned into a number of discrete intervals (g=1... G), starting with the highest energy particles (g=l), and ending with the lowest (g=G). The particles in energy group g are those with energies between Eg1 and Eg. The multigroup crosssections for a generic reaction process x are defined as E1 dE dOQa (, E)yf(F, E, Q) g(F) = H 4;T (2.1) f'g l dE4 dQ^(F,E,,Q) Based on the definition given in Eq. 2.1, the group constants are defined in Eqs. 2.2, 2.3, and 2.4, for the "total," "fission" and "scattering" processes, respectively. t() g dE ) (2.2) gdEf dQy (F,E,Q) dEg1 dE 4OQf (F, E)V(F, E, Q) Of(, () ,= (2.3) g 1dE d (F, E,Q) JE 1 dE'S dn'z (F, E' E, E'.)y(F, E, Q') ( sgg(7)= 4 ' (2.4) E' l dE'I dV (F, E, ) fE, 4;T With the group constants defined above, the multigroup formulation of the transport equation is written as G Q V /i (F, Q) + Ug, (F)Vg (F, Q) = JI'Jgg (F, YQ'Q)g, (F, Q') g V'= 14 (2.5) 1 G + x _Cv f,g,(rF)Og,(rF) + q9( ,), g'=l for g 1, G, where the angular flux in group g is defined as g,(rF,) = J ldEj (F,E,Q ). (2.6) fg In Eq. 2.5, q (F, Q) is the angular dependent fixed source; in general, for criticality eigenvalue problems, this term is set to zero. The scalar flux in Eq. 2.5 is defined as g (F) = JdfQ (Fg,). (2.7) 4;T In a 3D Cartesian geometry, the "streaming" term can be expressed as a a a .V=, + +7 (2.8) Ox Qy Oz where the direction cosines are defined as u=a Qi, r = j, 5 =Qk. Figure 21 shows the Cartesian spaceangle system of coordinates in three dimensions. (2.9) Figure 21. Cartesian spaceangle coordinates system in 3D geometry. The multigroup transport equation, with the scattering kernel expanded in terms of Legendre polynomials and the angular flux in terms of spherical harmonics is given by Eq. 2.10. The complete derivation of the scattering kernel expansion in spherical harmonics shown in Eq. 2.10 is given in Appendix A. P +7 + ) +JCg (x,Y,z)JVif (x,y,z,,u(P)= G L 1 (2/ + 1)o g (x, yz){ ()q0' (, y, z) + g'= ll= (2.10) 21 kP" ( g)[ (X, y, Z) cos(kp) + q~,g (X, y, z) sin(ko)] } k=1 ( + k)! + g f v g,g (x, y, yz)O,g (x, y, z) + q (F, ), k g1 where u : direction cosine along the xaxis r : direction cosine along the yaxis S: direction cosine along the zaxis ag: total macroscopic crosssection (: azimuthal angle, i.e. arctan  Vyg (x, y,, z, ), (): angular flux in energy group g ,sl,g' g: 1 t moment of the macroscopic transfer crosssection P, (u) : /th Legendre Polynomial Sg (//): /th flux moment Pk (/) : associated th, kth Legendre Polynomial kc,g (/)): cosine associated th, kth flux moment u,g (u): sine associated th, kth flux moment ,g: group fission spectrum k: criticality eigenvalue va ,g : fission neutron generation crosssection 2.1.2 Discretization of the Angular Variable The angular variable, 0, in the transport equation is discretized by considering a finite number of directions, and the angular flux is evaluated only along these directions. Each discrete direction can be visualized as a point on the surface of a unit sphere with an associated surface area which mathematically corresponds to the weight of the integration scheme. The combination of the discrete directions and the corresponding weights is 18 referred to as quadrature set. In general, quadrature sets should satisfy the following properties:3 * The associated weights must be positive and normalized to a constant, usually chosen to be one M .m =1.0. (2.11) m=l * The quadrature set is usually chosen to be symmetric over the unit sphere, so the solution will be invariant with respect to a 90degree axis rotation and reflection. This condition results in the oddmoments of the direction cosines having the following property ,n = Wmr7m = C',m =0.0, for n odd. (2.12) m=l m=l m=l * The quadrature set must lead to accurate values for moments of the angular flux (i.e., scalar flux, currents); this requirement is satisfied by the following conditions on the evenmoments of the direction cosines as follows M M M 1 IWm/ = Wmr,' = Wm ,for n even. (2.13) m=l m=l m=l + 1 A widely used method for generating a quadrature set is the levelsymmetric technique (LQN). In this technique, the directions are ordered on each octant of the unit sphere along the zaxis () on N/2 distinct levels. The number of directions on each level N is equal to +1, for i=1, N/2. It is worth noting that in 3D geometries, the total 2 number of directions is M=N(N+2), where N is the order of the SN method. 22 N Considering + + 2 + = 1 and i + j +k = + 2, where N refers to the number 2 of levels and i, j, k are the indices of the direction cosines, we derive a formulation for determining the directions as follows A2 = ,2 +(i1)A, (2.14) where 2(1 3 2)2 N 1 = 2(1 2) for 2< i< ,and 0< 1,2 < (2.15) (N 2) 2 3 In Eq. 2.14 the choice of ul determines the distribution of directions on the octant. If the value of/ul is small, the ordinates will be clustered near the poles of the sphere; alternatively, if the value of kl is large, the ordinates will be placed far from the poles. The weight associated to each direction, called a point weight, is then evaluated with another set of equations. For example, in the case of an Ss levelsymmetric quadrature set, this condition can be formulated as follows 2p, + 2p2 =i, (2.16a) 2p2 +3 w2, (2.16b) 2p2 w3, (2.16c) lp1 = w4, (2.16d) wherepi, p2 andp3 are point weights associated with each direction, and w1, w2, w3, W4 are the weights associated with the levels, as shown in Figure 22. point weights (p) Figure 22. Point weight arrangement for a Ss levelsymmetric quadrature set. Level weights w4  w3 __ w2 __ wi As an example, Figure 23 shows the S20 LQN quadrature set for one octant of the unit sphere. Figure 23. S20 LQN quadrature set. Note that, this quadrature set is limited by unphysical negative weights beyond order S20. Therefore, if a higher order quadrature set is needed beyond S20, another formulation has to be developed, which satisfies the even and oddmoments conditions. To address this issue, I have developed new quadrature techniques based on the Gauss Legendre quadrature formula and on the Chebyshev polynomials. 2.1.3 Discretization of the Spatial Variable The linear Boltzmann equation, given in Eq. 2.10, can be rewritten in an abbreviated form as Pm + 7my + + m C (F) ,f,(F) =Qmg(F(), (2.17) form = 1,M and g=1,G. The angle and energy dependence are denoted by the indices m and g, respectively. The right hand side of Eq. 2.17 represents the sum of the scattering, fission and external sources. The spatial domain is partitioned into computational cells, bounded by x1/2, X3/2,..., xI+1i2; yi/2, y3/2, ..., yJ+~/2; z1/2, z3/2, ..., ZK+12. The crosssections are assumed to be constant within each cell and they are denoted by cr,,,k. Eq. 2.17 is then integrated over the cell volume V,k A Ayj Azk, and then divided by the cell volume to obtain the volume and surfaceaveraged fluxes, therefore reducing to P" (w +/1 2,j,k,m,g T 1 2,j,k,m,g ( 1,j+1 2,k,m,g T 1 j 2,k,m,g) J(2.18) + jM 1,k+1 2,m,g ,j,k11 2,m,g ) 7 i,j,k i,j,k,m,g i,j ,k,m,g Azk In Eq. 2.18, the indices i, j, k represent the cellcenter values, while i+1/2,j+1/2, k+1/2 refer to the surface values. 2.1.4 Differencing Schemes For the SN method, different classes of differencing schemes are available. Low order differencing schemes require only the angular fluxes, and the average values at the cell boundaries to be related at the cell average value. Various forms of Weighted Difference (WD) schemes belong to this class. Highorder differencing schemes require higher order moments, and may be linear or nonlinear. Discontinuous, characteristic, and exponential schemes are examples of highorder differencing schemes.35 The solution of the SN equations is obtained by marching along the discrete directions generated in each octant of the unit sphere; this process is usually referred to as a transport sweep.3 For each computational cell, the angular fluxes on the three incoming surfaces are already known, from a previous calculation or boundary conditions. The cell center fluxes and the fluxes on the three outgoing surfaces must be calculated, hence additional relationships are needed. The additional relationships are referred to as the "differencing schemes ". The general form of WD schemes can be expressed as T ,j,k,mg al,J,k,m,g T +1 /2,j,k,,g ,,,k,m,g)R 1/2,j,k,m,g, (2.19a) ~ ,j,k,m,g ,,k,g b,k,, +1/2,k,m,g +(1 b,,j,k,m,g ) ,J 1/2,k,m,g, (2.19b) Y,,j,k,m,g = C,j,k,m,g Tj,,k+1/2,m,g + ( C,,J,k,m,g )Tj,k1/2,m,g (2.19c) The values a,,j,km,g, bi,j,k,m,g, and c,,k,m,g are determined based on the type of weighted scheme employed. 2.1.4.1 LinearDiamond Scheme (LD) In the LD scheme, the cellaverage flux is an arithmetic average of any two opposite boundary fluxes; hence the weights are set to constant values a=b=c=1/2. ,,j,k,m,g ( z +1/2,j,k,m,g z 1/2,j,k,m,g) 2 i ,J+l/2,k,m,g ,J1/2,k,m,g2.20) I (q ,j,k+1/2,m,g + j7,k1/2m,g) For example, in the direction /, > 0, 77 > 0, > 0 the outgoing fluxes are obtained as follows S1+1/2,j,k,m,g = 2T,,,k,m,g ~ 1/2,],k,m,g, (2.21a) S,j+1/2,k,m,g = 2 ,km ,j1/2,k,m,g (2.21b) T,],k+1/2,m,g = 2i,],k,m,g ,],k1/2,m,g 2 (2.21c) We can then eliminate the fluxes on the outgoing surfaces in Eq. 2.18 and obtain the centercell angular flux r"k Y + 12m g +"iQJ + Qek 2u 1/2,j,k,m,g + 2r Tl,j1/2,k,m,g + 2 T ,j,,k1 /2,m,g [ i,j,k ,,m,g (2.22) j, +2m 2u 27 2( Ax, Ayj Azk The outgoing fluxes are then evaluated using Eqs. 2.21a, b, and c. The LD differencing scheme may yield negative angular fluxes in regions where the flux gradient is large, even if the incoming fluxes and scattering source are positive. In this case, one approach is to set negative fluxes equal to zero, and then the cellaverage flux is recalculated to preserve the balance of particles. This approach is referred to as negative flux fixup (DZ). The DZ scheme performs better than the LD in practical applications, but the linearity and accuracy of the LD equations is not preserved. 2.1.4.2 Directional ThetaWeighted Scheme (DTW) This scheme uses a directionbased parameter to obtain the weighting factors a, b, and c which are restricted to the range 0.5 and 1.0. The DTW scheme uses a direction based parameter to obtain an angular flux weighting factor, which ensures positivity of the angular flux and removes the oscillations due to the spatial and angular discretization.36 The DTW average cell angular flux is given by ql,j,k + aI/A + Vlny + V a,, ,mgx b y c Az V ,j,,mg a,],k,m,g i,],k,m,g Cj,k,m,gAz (2.23) U M n17 m n n a,,km,gAx b,,k,m,g C,,,k,m,gA In the DTW scheme, the weights (a,,km,g, by,km,g, c.,,km,,g) are restricted to the range between 0.5 and 1.0, approaching second order accuracy when all weights are equal to 0.5, which is in this case is equivalent to the LD scheme. 2.1.4.3 Exponential DirectionalWeighted Scheme (EDW) The Exponential Directional Weighted (EDW) differencing scheme34, implemented in PENTRAN, is a predictorcorrector scheme, which utilizes the DTW to predict a solution that is then corrected using an exponential fit. The EDW is an inherently positive scheme, and the auxiliary equations derived for this method are given in Eq. 2.24. Vm (x,y, z) = aexp(AP, (x)/ /pm )exp(A P,(y)/ r7, )exp(AkPl (z) /n m). (2.24) The DTW scheme is used to calculate the angular fluxes (f) needed for the estimation of the coefficients A,, A and Ak given in Eqs. 2.25a, b, and c, respectively. A, t,x Vmx )Pm, (2.25a) 2 4 A ( outjy Vny)1m (2.25b) 2 A Ak out,z fin,z )m (2.25c) 2R, where the subscripts in and out refers to the incoming and outgoing surface averaged angular fluxes. The cellaverage angular flux formulated with the EDW scheme is given in Eq. 2.26. 2A 2A 2e A V4 = exp( )1 exp( ) lexp(Ak) 1 S+ m Vx + i + m V (2.26) SAhx y Az where fl is calculated using the coefficients given in Eqs. 2.25a, b, and c.37 2.1.5 The Flux Moments The flux moments are obtained from the angular fluxes using the following formulation ,,Jk,l = Wm ,,kP (m). (2.27) m=l Note that Eq. 2.27 for /=0 yields the scalar flux. The Associated Legendre moments are calculated using Eqs. 2.28 and 2.29. M ijk,C, = Wm T,,j,kPm() COS(I. .), (2.28) m=l M ",,k,Sl = m ,,,k P" ( m ) i n( .) (2.29) m=l 2.1.6 Boundary Conditions Three major boundary conditions can be expressed with a general formula as YT( ,m) = aT(, i'), (2.30) where 'm n = 'm n. Depending on the value of coefficient a, the three boundary conditions are: a=l, reflective boundary condition. a=0, vacuum or nonreentrant boundary condition. a=P, albedo boundary condition. 2.2 Source Iteration Method Due to the integrodifferential nature of the transport equation, the solution of the multigroup discrete ordinates equations is obtained by means of an iterative process, named the source iteration.3 The source iteration method consists of guessing a source (i.e., ingroup scattering source), then sweeping through the angular, spatial and energy domains of the discretized system with the appropriate boundary conditions. When the sweep is completed, integral quantities such as scalar flux and flux moments are obtained from the angular fluxes, and then a new ingroup scattering source is calculated, and the iteration process continues until a convergence criterion, shown in Eq. 2.31, is satisfied. Typical tolerances for fixed source calculations range in the order of 1.0e3 to 1.0e4. < e. (2.31) If fission and/or upscattering processes are present, outer iterations are performed on the fission and transfer scattering sources. Acceleration techniques may be applied between source iterations to speedup the convergence rate by determining a better guess for the flux moments and the source. 2.3 Power Iteration Method Criticality eigenvalue problems are solved using the method of power iteration.38 For this method, it is assumed that the eigenvalue problem has a largest positive eigenvalue, k>O, with an associated fission distribution F(F) that is nonnegative. Hence by considering ko>O and F (F) > 0 as initial guesses, we calculate the eigenvalue at iteration i as follows JdE dfF' 1(F) k1[ = k' fdEf (2.32) fdEfdrF' (F) where F' (F) = voa (F, E)q~ (F, E) is the fission source distribution at iteration i. The iterative process is continued until the desired convergence is reached, as shown in Eq. 2.33. k' k'1 1 < E. (2.33) k"' Generally, the tolerance required for criticality calculation is 1. Oe4 to 1.0e6. 2.4 Acceleration Algorithms for the SN Method Many acceleration methods have been proposed to speedup the convergence of the iterative methods used to solve the steadystate transport equation.17 There are a number of problems where standard nonaccelerated iterative methods converge too slowly to be practical. Most of these problems are characterized by optically thick regions with scattering ratio near unity. The three major acceleration approaches are the Coarse Mesh Rebalance (CMR), Multigrid, and Synthetic methods. The CMR approach is based on the fact that the converged solution must satisfy the particle balance equation.3 By imposing this balance condition on the unconverged solution over coarse regions of the problem domain, it is possible to obtain an iteration procedure that usually converges more rapidly to the correct solution. However, this method is highly susceptible to the choice of the coarse mesh structure and can be unstable. The multigrid approach has been used to accelerate the SN calculations; the basic principle of the method is to solve the equations on a coarse grid and to project the solution onto a finer grid. Different types of multigrid approaches exist, such as "/" Slash cycle, VCycle and WCycle, and the Simplified Angular Multigrid39 (SAM). In Chapter 8, I will compare the results obtained with the EPSSN synthetic acceleration and the SAM method. The synthetic acceleration approach is based on using a lowerorder operator, generally diffusion theory, as a means to accelerate the numerical solution of the transport equation.17 In the late 1960s, Gelbard and Hageman developed a synthetic acceleration method based on the diffusion and the S4 equations.28 Later, Reed independently derived a similar synthetic acceleration scheme40 and pointed out some limitations of the method derived by Gelbard and Hageman. The synthetic method developed by Reed has the advantage of being very effective for small mesh sizes, but it is unstable for mesh sizes greater than 1 mfp (mean free path). Later, Alcouffe independently derived the Diffusion Synthetic Acceleration (DSA) method.29 He addressed the issue of stability of the method and derived an unconditionally stable DSA algorithm. Alcouffe pointed out that in order to obtain an unconditionally stable method, the diffusion equation must be derived consistently from the discretized version of the transport equation. In this way, the consistency between the two operators is preserved. Recently, the DSA method has been found to be ineffective3032 for multidimensional problems with strong heterogeneities, even when a consistent discretization of the lower order and transport operators is performed. CHAPTER 3 ADVANCED QUADRATURE SETS FOR THE SN METHOD This chapter covers the development of advanced quadrature sets for solving the neutron transport equation via the Discrete Ordinates (SN) method. The levelsymmetric (LQN) quadrature set is the standard quadrature set for SN calculations; although, as discussed in the previous chapter, this quadrature set is limited to order 20. The Equal Weight (EW) quadrature set has been proposed to resolve the issue of negative weights; this quadrature set is generated by partitioning the unit sphere into M directions, where M = N(N + 2) and by assigning an equal weight to each direction w, = The EW M quadrature set yields positive weights for any SN order; however, it does not completely satisfy the evenmoment conditions given in Eqs. 213. I have developed and tested new quadrature sets based on the Legendre (PN) and Chebyshev (TN) polynomials. In this chapter, I discuss the Legendre EqualWeight (PN EW), the LegendreChebyshev (PNTN) quadrature sets, and the Regional Angular Refinement (RAR) technique for local angular refinements. The PNEW and PNTN have no limitations on the number of directions, and the RAR technique is an alternative to the Ordinate Splitting (OS) technique.1l The OS technique has been developed to refine the directions of a standard quadrature set using equalspaced and equalweight directions, while the RAR utilizes the PNTN quadrature set over a subset of ordinates. The main difference between the OS and RAR techniques is in the refining methodology; the OS technique focuses on the refinement of each single direction, while the RAR considers a sector of the unit sphere. Also the biased quadrature set generated with RAR satisfies the conditions on the odd and evenmoments of the direction cosines (Eq. 212 and 213). To examine the effectiveness of the new techniques for angular quadrature generation, each technique has been implemented into the PENTRAN code10 (Parallel Environment NeutralParticle TRANsport), and utilized for a number of problems of practical interest. 3.1 Legendre EqualWeight (PNEW) Quadrature Set In order to develop a quadrature set which is not limited to order S20, I have investigated the GaussLegendre quadrature technique.4 This quadrature set is characterized by the same arrangement of directions as the LQN, but the directions and weights are evaluated differently. Given the SN order for the discrete set of directions, we apply the GaussLegendre quadrature formula using the following recursive formulation (j + 1)P,+1 = (2j + 1)4P7 jPl, forj =0, N, (3.1) where 1< <1, P_() = 0, and P() = 1. (3.2) The levels or polar angles, along the zaxis are set equal to the roots of Eq. 3.1. The values represent the levels of the quadrature set. Once we have evaluated the ordinates along the zaxis, we obtain the weights associated with each level using the following recursive formulation 2 N w2 ,for i = 1, (3.3) [l d J In order to complete the definition of each discrete direction, the azimuthal angle is evaluated on each level by equally partitioning a 90 degree angle into N  i + 2 angular 2 intervals, where i=1, N/2. Hence, the weight associated with each direction is given by w , for i = ... I (3.4) In Eq. 3.4, j =1... i +1 is the number of directions with equal weights on the ith 2 level. Figure 31 shows the directions and the associated weights for an S28 PNEW quadrature set on one octant of the unit sphere. wdej@t 00106670 0D01 066?3 0010AE68 o o 1CO 0010068 O OC9W33 oOON6648M 0009A6437 00092M388 00o c634 O0tJ3M292 000366243 DOCWAMi 0006i 47 GO DEM88 Figure 31. S28 PNEW quadrature set. Note that in Figure 31, all directions on the same level have the same weight, as indicated by the color. 3.2 LegendreChebyshev (PNTN) Quadrature Set In the PNTN quadrature set, similar to PNEW, we choose the levels on the zaxis equal to the roots of the GaussLegendre quadrature formula given in Eqs. 3.1 and 3.2; however, the azimuthal angles on each level are set equal to the roots of the Chebyshev (TN) polynomials of the first kind. Chebyshev polynomials of the first kind are formulated as follows T, [cos(o)] cos(/o). (3.5) The Chebyshev polynomials are orthogonal and satisfy the following conditions 0, I# k dyT()(y y)( y2) 1/2 =z,l= k = 0 Izr/2,1 =k 0 y =cos(m) (3.6) Again, using the ordering of the LQN quadrature set, we define the azimuthal angles on each level using the following formulation 0J = 2i2j+l1) (3.7) where o), e 0, 2) and i= 1, N/2. In Eq. 3.7, i is the level number and j = 1... i +1. The level and point weights are 2 generated in the same way as for the PNEW. Note that both PNEW and PNTN quadrature sets do not present negative weights for SN orders higher than 20. Figure 32 shows the directions and the associated weights for an S28 PNTN quadrature set on one octant of the unit sphere. Figure 32. S40 PNTN quadrature set over the unit sphere. 3.3 The Regional Angular Refinement (RAR) Technique The RAR method is developed for solving problems with highly peaked angular fluxes and/or sources. The approach consists of two steps. In the first step, we derive a PNTN quadrature set of arbitrary order on one octant of the unit sphere, as described in Section 3.2. In the second step, we define the area (angular segment) of the octant to be refined along with the order (N') of the PN'TN' quadrature set to be used in this region. The area to be refined is characterized by the polar range (min, max) and the azimuthal range (pmin, pmax). Generally, the biased region is selected based on the physical properties of the model. For example, if a directional source is forward peaked along the xaxis, the quadrature set will be refined on the pole along the xaxis. The levels of the PN'TN' quadrature set are calculated using Eqs. 3.1 and 3.2; therefore, they are mapped onto the (min, max) subdomain using the following formula: max 2 + max + mm for i=1, N'2. (3.8) 2 2 Hence, the azimuthal angles are evaluated using the Chebyshev polynomials as follows (3.9) where co, i 0o, i=1, N'/2 andj=l, '/2. S2 (N,')2 The number of directions in the refined region is equal to The weights in the S2 refined region are renormalized to preserve the overall normalization on the unit sphere. Figure 33 shows the RAR technique applied to an S16 PNTN quadrature set. In the biased region, which extends from =0.0 to =0.2 along the zaxis, and from qp=0 to (p=10 on the azimuthal plane, an Slo PNTN quadrature set is fitted. Figure 33. PNTN quadrature set (S16) refined with the RAR technique. 3.4 Analysis of the Accuracy of the PNEW and PNTN Quadrature Sets The main advantage of the LQN technique is the fact that it preserves moments of the direction cosines, thereby leading to an accurate solution. The PNEW and PNTN 0)i ) N'2j + 1 (Pmx (Pm +( max + (3,mm 2 N' 2 2 quadrature sets attempt to preserve these quantities independently along the a, fl, and 5 axes. Therefore, I verified the capability of the new quadrature sets in preserving the even moments of the direction cosines, which are directly related to the accuracy of the quadrature set. Table 31 compares the evenmoments of the direction cosines calculated with an S30 PNEW quadrature set with the exact value. Table 31. Evenmoments obtained with a PNEW S30 quadrature set. M 1=1 0.333333333 0.194318587 0.135907964 0.103874123 0.083691837 0.06983611 0.059749561 0.052087133 0.046074419 0.041234348 0.037257143 0.033932989 0.031114781 0.028696358 0.026599209 0.333333333 0.194318587 0.135907964 0.103874123 0.083691837 0.06983611 0.059749561 0.052087133 0.046074419 0.041234348 0.037257143 0.033932989 0.031114781 0.028696358 0.026599209 0.333333333 0.2 0.142857143 0.111111111 0.090909091 0.076923077 0.066666667 0.058823529 0.052631579 0.047619048 0.043478261 0.043478261 0.037037037 0.034482759 0.032258065 Exact value 1+n 0.333333333 0.2 0.142857143 0.111111111 0.090909091 0.076923077 0.066666667 0.058823529 0.052631579 0.047619048 0.043478261 0.043478261 0.037037037 0.034482759 0.032258065 As expected, the PNEW preserves exactly the evenmoments conditions along the axis, while on the a, rlaxes these conditions are only partially preserved; the maximum relative difference between the evenmoments calculated with PNEW and the exact solution is 17.0%. Table 32 shows the comparison of the evenmoments evaluated with an S30 PNTN quadrature set and the exact value. Moment Order (n even) 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 Table 32. Evenmoments obtained with a PNTN S30 quadrature set. M Moment Order (n even) 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 0.333333333 0.199999962 0.142857143 0.111111111 0.090909091 0.076923077 0.066666667 0.058823529 0.052631579 0.047619048 0.043478261 0.043478261 0.037037037 0.034482759 0.032258065 0.333333333 0.2 0.142857143 0.111111111 0.090909091 0.076923077 0.066666667 0.058823529 0.052631579 0.047619048 0.043478261 0.043478261 0.037037037 0.034482759 0.032258065 Exact value 1+n 0.333333333 0.2 0.142857143 0.111111111 0.090909091 0.076923077 0.066666667 0.058823529 0.052631579 0.047619048 0.043478261 0.043478261 0.037037037 0.034482759 0.032258065 It is clear from Table 32 that the PNTN quadrature set completely satisfies the evenmoment conditions. This is possible because both roots of Legendre and Chebyshev polynomials satisfy the evenmoment conditions given by Eqs. 2.13. In order to further verify the accuracy of the PNEW and PNTN quadrature sets, and to check their accuracy, I used a simple test problem, consisting of a homogeneous parallelepiped, where an isotropic source is placed in its lower left corner as shown in Figure 34. Due to the symmetry of the problem, it is expected that the particle currents flowing out of regions A and B (see Figure 34) have the same value. 0.333333333 0.199999962 0.142857143 0.111111111 0.090909091 0.076923077 0.066666667 0.058823529 0.052631579 0.047619048 0.043478261 0.043478261 0.037037037 0.034482759 0.032258065 HUMh a 9Z 25JU im UuIEN ar~flEBD LIC. Y k.9 QK EL Y 9161 (Wago ROeIKkW B.C. I rnl V3CLPJM B Figure 34. Configuration of the test problem for the validation of the quadrature sets. Figure 35 shows the relative difference between the particle current in regions A and B for the EW, PNEW, and PNTN as compared to the LQN technique. 12 S10 8 2  5 4 So   8 10 12 14 16 18 20 Sn Order EW PnTn  PnEW Figure 35. Relative difference between the currents Jx and Jz for the test problem. Levelsymmetric is considered as the reference because it preserves moments of both azimuthal and polar direction cosines. The PNTN yields almost perfect symmetry, while the PNEW and Equal Weight show maximum relative differences of 4% and 10%, respectively. It is worth noting that the loss in accuracy of the EW quadrature set is attributed to the fact that the evenmoment conditions are not satisfied. The PNEW yields higher accuracy compared to EW, because the evenmoment conditions are satisfied along the zaxis. 3.5 Testing the Effectiveness of the New Quadrature Sets In this section, the effectiveness of the new quadrature sets is examined by simulating two test problems: Kobayashi benchmark problem 3 and a CTScan device for industrial/medical imaging applications. 3.5.1 Kobayashi Benchmark Problem 3 To examine the effectiveness of the new quadrature sets, I have used the first axial slice of the Kobayashi43 3D benchmark problem 3 with pure absorber. Figures 36 show two different mesh distributions: Figure 36A is obtained from a previous study42, where an appropriate variable mesh was developed; Figure 36B shows a uniform mesh distribution that I have developed for the current study. The uniform mesh is used in order to separate the effects of the angular discretization from the spatial discretization. The reference semianalytical solutions are evaluated in two spatial zones shown in Figures 36A and 36B (zone 1 along yaxis, at every 10.0 cm intervals between 5.0 and 95.0 cm; zone 2 along xaxis, y = 55.0 cm, every 10.0 cm, between 5.0 cm and 55.0 cm). A B 5  t 5 5 55 5 25 4 5 55 x Figur sh distribution for the obayashi benhark probe ar e distribution shown in Figure 36a. In the previous study, by taking advantage of the l  ""  Figuriable 3.6. Mesh distribution for the obayashi benchmark problevelsymmetric quadraturiable mesh presented a maximum relative error of 6% in zone 1. In the current study, the solution disobtained with the levelsymmetric quadrature set and uniform spatial mesh yields attribution. maximum relative 37 shows the ratio of10% in zone calculated to the exact solution (C/E) for the level symmetric, PEW and PNEW and PTN quadrature sets undeof order 20 for zonestimate 1. Also, in this figure,x by 51.9%I present and 8.5%, solution34 obtained in a previous st point of zone 1; the variables due to the fact that the distribution shown in Figure 36a. In the previous study, byaround tweaking advantage ofmpared to the TNvariable mesh distribution, the solution obtained with the levelsymmetric quadrature sets presented a maximum relative error of 6% in zone 1. In the current study, the solution obtained with the levelsymmetric quadrature set and uniform spatial mesh yields a maximum relative error of10% in zone 1. In zone 1, the PNEW and PNTN quadrature sets underestimate the scalar flux by 51.9% and 8.5%, respectively, on the last point of zone 1; this is due to the fact that the PNEW quadrature set has fewer directions clustered around the yaxis as compared to the PNTN and levelsymmetric quadrature sets. 1.20E+00 1.00E+00  .2 8.00E01 S 6.00E01 O 4.00E01 2.00E01 0.00E+00 5 15 25 35 45 55 65 75 85 95 Y (cm) LQn S20 (Uniform mesh) LQn S20 (Variable mesh)  PnTn S20 X PnEW S20 Figure 37. Comparison of S20 quadrature sets in zone 1 at x=5.0 cm and z=5.0 cm. Figure 38 compares the scalar flux obtained in zone 2 of the benchmark problem. While using a uniform spatial mesh, the PNTN quadrature set yields slightly better results compared to levelsymmetric. 1.40E+00 1.20E+00 0 1.OOE+00 utj 6.00E01 0 4.00E01 2.00E01 2.00E01  0.OOE+00 5 15 25 35 45 55 X (cm) LQn S20 (Uniform mesh) LQn S20 (Variable mesh)  PnTn S20 PnEW S20 Figure 38. Comparison of S20 quadrature sets in zone 2 aty=55.0 cm and z=5.0 cm. In zone 2, the maximum relative error obtained with PNTN is 18.6%, while for levelsymmetric it is 21.9% using the uniform mesh distribution, and 6% using variable meshing. However, an error of 28.2% is observed for the PNEW quadrature set. Figures 39 and 310, show the solutions obtained with the PNEW quadrature set for different SN orders compared to levelsymmetric S20, in zone 1 and 2 respectively. 1.10E+00 1.00E+00  9.00E01 8.00E01 ,u 7.00E01 6.00E01 5.00E01 4.00E01 5 15 25 35 45 55 65 75 85 95 Y (cm)  LQn S20 I PnEW S20  PnEW S22  PnEW S24  PnEW S26 Figure 39. Comparison of PNEW quadrature sets for different SN orders in zone 1 at x=5.0 cm and z=5.0 cm. 1.40E+00 1.30E+00 1.20E+00 1.10E+00        % 1.10E+00 ,u 1.00E+00 9.00E01 8.00E01 7.00E01 5 15 25 35 45 55 X (cm) 4 LQn S20  PnEW S20  PnEW S22 PnEW S24  PnEW S26 Figure 310. Comparison of PNEW quadrature sets for different SN orders in zone 2 at y=55.0 cm and z=5.0 cm. In zone 1 (Figure 39), the PNEW is not as accurate as levelsymmetric, because fewer directions are clustered near the yaxis; however in zone 1, the solution improves somewhat by increasing the SN order. In zone 2 (Figure 310) the PNEW set yields inaccurate results, with a maximum relative error of 36% for the S20 case. Figure 311 compares the ratios of different computed solutions to the exact solution; the computed solutions wered obtained with the PNTN quadrature set for orders S20, S22, S24, S26 and with the S20 levelsymmetric quadrature set. It appears that the increase in the quadrature order does not have a noticeable effect in improving the accuracy. However, this behavior can be attributed to the fact we have retained the same spatial mesh discretization. 1.15E+00 1.10E+00  S1.05E+00 S1.00E+00 ,u 9.50E01 9.00E01 8.50E01 8.00E01 5 15 25 35 45 55 65 75 85 95 Y (cm) 4 LQn S20  PnTn S20  PnTn S22 X PnTn S24  PnTn S26 Figure 311. Comparison of PNTN quadrature sets for different SN orders in zone 1 at y=5.0 cm and z=5.0 cm. In zone 2 (Figure 312), the solution obtained with an S22 PNTN quadrature set is more accurate than what obtained with levelsymmetric. The S22 PNTN yields a maximum relative error of 9% compared to 22% from levelsymmetric. Again, in zone 2, the accuracy somewhat decreases as the SN order increases, because the spatial mesh is not consistently refined. 1.30E+00 1.20E+00  2 1.10E+00 W 1.00E+00 3 9.00E01 8.00E01 7.00E01 5 15 25 35 45 55 X (cm) LQn S20 PnTn S20 PnTn S22  PnTn S24  PnTn S26 Figure 312. Comparison of PNTN quadrature sets for different SN orders in zone 2 at y=55.0 cm and z=5.0 cm. 3.5.2 CTScan Device for Medical/Industrial Imaging Applications The model of a CTScan device used for medical/industrial applications is used in this section to verify the accuracy and performance of the RAR technique. A CTScan device utilizes a collimated xray source (fanbeam) to scan an object or a patient. The main components of a CTScan device are an xray source mounted on a rotating gantry and an array of sensors. The patient is positioned on a sliding bed that is moved inside the CTScan. The mesh distribution for this model, is shown in Figure 313. 20.0 cM 740 0m Aimay ofsewtrs Figure 313. Crosssectional view of the CTScan model on the xy plane. Figure 313 shows the simplified PENTRAN model which represents the xray directional source ("fan" beam), a large region of air and an array of sensors. The size of this model is 74 cm along the yaxis and 20 cm along the xaxis. The array of detectors is located at 72 cm from the source along the xaxis. The materials are described using onegroup crosssections from the 20group gamma of the BUGLE96 crosssections library. The group corresponds to an xray source emitting photons in an energy range of 100 KeV to 200 KeV. The crosssections were prepared using a P3 expansion for the scattering kernel. Because of the presence of large void regions and a directional source, the solution of the transport equation is significantly affected by the rayeffects.3 One remedy is to use high order quadrature sets with biasing, such as RAR. We compared the solutions obtained with an Ss0 PNTN quadrature set. The RAR technique has been applied to an S30 PNTN quadrature set; the biased region on the positive octant extends from z=0.0 cm to z=0.3 cm and the azimuthal angle extends from 0.0 to 5.0 degrees. In the biased region an Slo PNTN quadrature set is used. The PNTN quadrature set biased with RAR resulted in 142 directions per octant. The unbiased Ss0 PNTN quadrature set yielded 325 directions per octant. The S20 levelsymmetric quadrature set yielded 55 directions per octant. Figures 314, 315, and 316 show the flux distributions in the xy plane, obtained with the levelsymmetric S20, Ss0 PNTN and S30 PNTN with RAR, respectively. Fluxtpkm;/E) 1.54E18 1.22Ef4 g69E07 7,6 E407 20 26E4107 4 82E07 3 2E407 234EIO7 S10 1.93E407 1.91EH7 11 E407 1.19E~07 9.47EO6 7.93E411B 5 94E10 10 20 30 40 50 60 70 71 EI x (Ic) 3 73E4~ 2 .6EiO 1 51E.O5 Figure 314. Scalar flux distribution on the xy plane obtained with an S20 level symmetric quadrature set. Flux (p.ra1lE) 8.19E4B6 6.34' E" .D2E46O 2 3.23E20 30 40 50 60 70 X(Cm) 2.9E45G 2. 1956E10 1.24E4DG 5.D2E*05 10 20 30 40 50 60 70 16C4e0 x (cm) 2.SDEID t.24E4DS Figure 315. Scalar flux distribution on the xy plane obtained with an S50 PNTN quadrature set. I 10 >,,il 10 20 30 40 50 60 70 x (cm) FIIjh i..Ir,T, F I S6 EtCS 4.41 E+0 3A9E 06 2.77 E' 2 ,19 EW 1.74E+46 138E+06 1JOSE4c6 5 j5 E+05 5.43 E.+5 I4.31E105 3.41EB5 3.70E+OM 2.14E+D5 8.146E+05 .70 E+04 Figure 316. Scalar flux distribution on the xy plane obtained with an S30 PNTN quadrature set biased with RAR. The above results indicate that the levelsymmetric quadrature set exhibits significant rayeffects, while S50 PNTN and S30 PNTN with RAR quadrature sets, yield similar solutions without any rayeffects. The main advantage of using a biased quadrature set is the significant reduction in computational cost and memory requirement. Table 33 compares the CPU time and memory requirements for the three calculations presented above. Table 33. CPU time and total number of directions required for the CTScan simulation. Quadrature Set Directions CPU Memory ratio Time ratio Time(sec) S50 PNTN 2600 166.4 1.0 1.0 S30 PNTN RAR (S0o) 1136 79.4 0.51 0.47 S20 LQN 440 33.3 0.2 0.2 a memory and time ratio are referred to the S50 PNTN quadrature set. The RAR technique lessens the ray effect in the flux distribution and greatly reduces the computational time by more than a factor of 2 compared to S50. The new quadrature sets biased with the OS rather than the RAR technique have also been examined based on the CTScan model.9 Figure 317 compares the results of PENTRAN with a reference Monte Carlo solution. For all cases, the first direction of the lowest level in quadrature set is split in 9 or 25 directions; for example, PNTN 22255 corresponds to PNTN S22 with direction 55 split in 9 directions. All the quadrature sets biased with the OS technique yield accurate results within the statistical uncertainty of the Monte Carlo predictions. Due to the significant rayeffects, the levelsymmetric S20 quadrature set without ordinate splitting yields poor accuracy. Note that, even by using high order quadrature sets, such as PNTN S28 (840 directions), the solution at detector position is under predicted by 21%. 9.00E+08 8.50E+08  8.00E+08 7.50E+08 x 7.00E+08  > 6.50E+08 >c 6.00E+08 \ 5.50E+08  5.00E+08 S20LS (No OS) 4.50E+08 4.00E+08 0 2 4 6 8 10 12 14 16 18 20 YAxis Position (cm) PnTn22255 PnTn20346 MCNP Tally Flux S20LS346 S20EW546 Figure 3.17. Comparison of the scalar flux at detector position (x=72.0 cm). CHAPTER 4 DERIVATION OF THE EVENPARITY SIMPLIFIED SN EQUATIONS This chapter presents the initial derivation of the Simplified Spherical Harmonics (SPN) equations starting from the PN equations in 1D geometry, and it discusses the issues related to the coupling of the SPN moments on the vacuum boundary conditions. Because of this peculiarity, the implementation of the general SPN equations into a computer code proved to be cumbersome. However, I will present the initial derivation of the SP3 equations, successively implemented into a new computer code named PENSP3 (Parallel Environment Neutralparticles SP3). To overcome the difficulties related to the coupling of the SPN moments in the vacuum boundary conditions, I adopted a different formulation based on the EvenParity Simplified SN (EPSSN) equations. These equations are derived starting from the 1D SN equations, and using the same assumptions made for the derivation of the SPN equations; however, the main advantage of this formulation is the natural decoupling of the even parity angular fluxes for the vacuum boundary conditions. Therefore, a Fourier analysis of the EPSSN equations will follow, along with the derivation of a new formulation to accelerate the convergence of the source iteration method applied to the EPSSN equations. This chapter is concluded with the derivation of the 3D P1 equations. I will compare the P1 equations with the SP1 equations, and I will describe the assumptions made in the derivation of the SP1 equations and the relation with the spherical harmonics P1 formulation. 4.1 Derivation of the Simplified Spherical Harmonics (SPN) Equations The SPN equations were initially proposed by Gelbard18 in the early 1960s. However, they did not receive much attention due to the weak theoretical support. Recently, the theoretical foundations of the SPN equations have been significantly strengthened using a variational analysis approach.2122 I derive the multigroup SPN equations starting from the 1D multigroup PN equations and by applying the procedure originally outlined by Gelbard. The multigroup 1D PN equations are given by n+1 a,^ (X) n QO i (X) r 1 +, + [1, g(X) Cs, (x)(., (x) = q,,,g(), (4.1) 2n +1 x 2n +1 x for n=0, N and g=1, G, where G L 2+1 s(+1 qg (x) P 4 (/)og, (x)Og, (X) + qf,g (x) + (x) (4.2) g'=1 1=0 4i 7 g' g In Eqs. 4.1 and 4.2, I defined the following quantities: o,g (x), total macroscopic crosssection in group g. os,,g (x), Legendre moment of the ingroup macroscopic scattering crosssection of order n. ,sg',,g (x), Legendre moment of the group transfer macroscopic scattering crosssection of order 1. ng (x), Legendre moment of the angular flux of order n. q,g(x), fission source. S, (x), Legendre moment of the inhomogeneous source of order n. G, total number of energy groups and L, order of the Legendre expansion of the macroscopic scattering crosssection (L In Eq. 4.1, the angular flux is expanded in terms of Legendre polynomials as L21+1 ,g (x, P) = 21 (p)1,g (x) (4.3) 0o 4;r In Eq. 4.1, the term g is defined to be identically zero when n=N. This dx assumption closes the PN equations, yielding N+ equations with N+ unknowns. The procedure prescribed by Gelbard to obtain the 3D SPN equations from the 1D PN, consists of the following steps: 4. Replace the partial derivative operator in Eq. 4.1 for even n with the divergence operator (V). 5. Replace the partial derivative operator in Eq. 4.1 for odd n with the gradient operator (V). By applying this procedure to Eq. 4.1, it reduces to n+l n 4.4a) V + (F)+ V (r) + tg () s g(), (F ,g (F) = q, (F) (4.4a) 2n +1 2n +1 forn=0,2,...,N1, g=1, G, and re V, n+l n Vn+ I, () + g, (V ) + tg () agg(F)ng () = q (r), (4.4b) 2n +1 2n +1 forn=1,3,...,N, g=1, G, and F V. The SPN equations can be reformulated in terms of a secondorder elliptic operator, if one solves for the oddparity moments using Eqs. 4.4b i.e., 1 n+1 1 n qg (F) (+) 1= n V(1 lg() nI 1g()) + i (4.5) ran,g 2n+ an,g 2n+1 Can,g where an,g = tt,g (F) sng,gg () (4.6) and then substitute Eqs. 4.5 into Eqs. 4.4a to obtain n+l n+2 1 n+l n+l 1 2n + 1 2n + 3 an+1 g 2n + 1 2n + 3 an+lg n n 1 n nIl 1( V Vng n2,g (4.7) 2n+1 2n1 oan,,g 2n+12n 1 oan,,g +n+1 q ng n qy , S g2 + 1 Can+,g 2n +1 Can,g for n=0,2,...,N1 andg=1, G. Note that for simplicity, the spatial dependency of q has been eliminated. The SPN equations yield a system of (N+1)/2 coupled partial differential equations that can be solved using standard iterative methods, such as GaussSeidel or Krylov subspace methods. The main disadvantage of this formulation is that it yields Marshaklike vacuum boundary conditions, coupled through the moments. Because of this issue the implementation of this formulation into a computer code becomes cumbersome. However, to study the effectiveness of the SPN method, I developed a 3D parallel SP3 code,25 PENSP3 (Parallel Environment Neutralparticles SP3). The PENSP3 code is based on Eqs. 4.8a and 4.8b, which are derived assuming isotropic scattering and isotropic inhomogeneous source. DI, (F)VF(F)+ t,,g ()F(F)= 2,,t (F02g (r)+ ,,g g ()00,o (r) + 1 GCoG C( (4.8a) + o^ (F' g' (F) + k g(F)' (F) + S (4.8a) g1=1 k 9 1 g' g D2,g 0'42,g 0 + g (F)02,g 0 5 g' gfor g 1, G, forg=1, G, where 1 9 (F) 202,g )+00g(F), Dg (F) = and D2,g(F) 3 g (4.9) The SP3 Marshaklike vacuum boundary conditions for Eqs. 4.8a and 4.8b are given by 3 (f )+ D,) 0( V)F (f)= i 2/+/b ( 2 f, p, (p)d/d (4.10a) 2 8 2,;g ()+ D ( V) () = F ()+2P (, )d .(4.10 Ob) 40 40 g Eqs. 4.10a and 4.10b are termed Marshaklike boundary conditions, because in 1D geometry they reduce to standard Marshak boundary conditions. Implementation of these formulations into a computer code is difficult because of the coupling of the SPN moments. The reflective boundary condition is represented by setting the oddmoments equal to zero on the boundary, i.e., i O (Fb)= 0, for n odd, (4.11) where Fb e OV and h is the normal to the surface considered. 4.2 Derivation of the EvenParity Simplified SN (EPSSN) Equations I have derived the EvenParity Simplified SN (EPSSN) equations starting from the 1D SN equations given by um a (X1 ) + (x)(x, )VX,l)= Q(X, um), form 1, N, (4.12) where L L Q(x, P/m)= (2/ +1)P ( ) (m (x),j (x) + (2/ + i)P1 (,m)SI (x) + qf (x), (4.13) 1l0 1/0 1 1 N qf(x) = v, (x)O (x); (x) = M WP (/Pm)(x,/ m). (4.14) k 2 m= A GaussLegendre symmetric quadrature set (PN) is considered, where /m e (1,1), M wm = 2.0, and M N(N+2). In Eq. 4.13, L is the order of the Legendre expansion for m=1 both the macroscopic scattering crosssection and the inhomogeneous source (L Therefore, the even and oddparity angular fluxes are defined by 1M = 1 [(x, m) + i(x,m)] (even), (4.15a) 2 w = [y(x, m) (x,/m )] (odd). (4.15b) 2 To reformulate the 1D SN equations in terms of the even and oddparity angular fluxes, I rewrite Eq. 4.12 for /m as m ) + r, (x)y(x,'m) = Q(x,/m), for m 1, N. (4.16) ox Then, I add Eqs. 4.12 and Eq. 4.16 to obtain am [(xm) ( m + ot (X)[),(X, Pm)+ m(X,m)]= Q(x, m) + Q(x,m), (4.17) and use the definitions of even and oddparity angular fluxes (given by Eqs. 4.15a and 4.15b), to obtain 2/1m O (x)+ 2c, (x), '(x) = (21+ l)[P~(/pm)+ P,(p/m ) ,,(x)O,(x)+ 1 0 (4.18) (21 + 1)[PI (pm) + P, (/m)]S, (x) + 2qf(x). 10 Consider the following identities for the Legendre polynomials.3 P, (/) = P, (p), for I even, (4.19a) 53 P( (P) = P (p/), for I odd. (4.19b) Eq. 4.18 can be rewritten as m m (x)+ C, (x)yi (x)= 1 (2/1+)P (m )[ (x)A (x) + +S (x)] + qf (x). (4.20) 1=0,2 even Similarly, by subtracting Eq. 4.16 from Eq. 4.12, I obtain m [(x, () X)] + )[y(x, ) (x,, )] Q(x, ,m) Q(x,/i), (4.21) and by using the definitions of even and oddparity angular fluxes, I obtain aE= L 2/,m l'(x) + 2 (x),y (x) = (2/+1)[Pl (,m)P (,m )] ,(x) l(x)+ L10 (4.22) Z(2 + 1)[PI(/m) P(/ ) (). 10 Following the use of the Legendre polynomial identities (Eqs. 4.19a and 4.19b), Eq. 4.22 reduces to P A (x) +o a() (x)= f (2 +1) (1 (pm)[a (x)o (x) + S (x)]. (4.23) x 1=1,3 odd Now, the oddparity angular fluxes are then obtained from Eq. 4.23 as Vm (x)= V y (x)+ 1 (2/+l)P, (p,,) )(x)cl (x)()+S,(x)]. (4.24) C (X) xt () 1t1,3 odd Then, using Eq. 4.24 in Eq. 4.20, I obtain Ia k /a (X) + ct (x)/ (X) = (2/ + 1)P (pm, )[a (x)x (x) + S (x)] x o (x) x 1=0,2 evn (4.25) a 'm (2/+ 1)PI (p,)[, (x)q (x) + S,(x)]+ qf (x). S( dd) 1=1,3 odd Finally, the EPSSN equations in 3D Cartesian geometry with anisotropic scattering kernel and anisotropic inhomogeneous source of arbitrary order L, are obtained by applying the procedure outlined by Gelbard (i.e., substitution of first order partial differential operators with the gradient operator) to Eq. 4.25. 2 L1 V Vl m(r)2 + Et(r))m (r) / (2/ + I)P, (km )[, (F)q1 (F) + S, (F)] ( ) 1=0,2 even S\ (4.26) () 01,3 odd for m 1, N/2, where (F) = m v (r)+ 1 (2/+l)P, (m) J (F) (F)+ (F) (4.27) ,t (r) m (F) 1,3 odd Due to the symmetry of the GaussLegendre quadrature set, the EPSSN equations only need to be solved on half of the angular domain, e.g. / e (0,1). The moments of the even and oddparity angular fluxes are evaluated by N/2 S(r)= ,WmP (,m )Vf (), for / even, (4.28a) m=l and N/2 ((r) = ZwmP (m )V (0), for odd. (4.28b) m=l The multigroup form of Eqs. 4.26 and 4.27 with anisotropic scattering and source are written as S VgE g) + (r)g, (F) = f (21 + 1)P, (m )c,,,g,'g (F)q,,g (r) even G L L1 V. g 1 (21+1)P,(,m) ,(),g,(F) + (2/ + )P,(,m,)S,,(F) tg (F )g'=1=1,3 1=0,2 odd even G L  (2/+l)P, (),u),, (F) +qf ,g() ut,g ( )g' =1 =1 3 (4.29) and t ) t,g ( ) g'=1=1,3 odd (4.30) form 11, N2 and g 1, G. 4.2.1 Boundary Conditions for the EPSSN Equations The boundary conditions for the EPSSN equations are based on the assumption that the angular flux on the boundary surface is azimuthally symmetric about the surface normal vector. The EPSSN boundary conditions follow directly from the 1D evenparity SN boundary conditions. Hence, by considering the positive half of the angular domain for the EPSSN equations, the 1D source boundary condition at the right boundary face is given by w (x,, P/ ) n (x,,/ )=V(x,/m), (4.31) and the corresponding condition in 3D is given by (F,,P) n ( 0 /(F Pr) = Vy( ,p) (4.32) Note that in 3D geometry, the EPSSN formulation requires the incoming boundary flux to be azimuthally symmetric about the surface normal vector. The 3D albedo boundary condition is given by n.(F,/m.) = (/m) (4.33) 1+a Note that in Eq. 4.33, the vacuum boundary condition is obtained by setting a = 0 while the reflective boundary condition is obtained by setting a = 1. Note that the main advantage of the EPSSN formulation compared to SPN is the decoupling of the evenparity angular fluxes for the vacuum boundary conditions. 4.2.2 Fourier Analysis of the EPSSN Equations The EPSSN equations are solved iteratively using the source iteration method. This method is based on performing iterative cycles on the scattering source; moreover, the method has a clear physical interpretation that allows one to predict classes of problems where it should yield fast convergence. The source iteration method for the EPSSN equations is defined by HL,m,g (F)m (r) = qm,gg +q', for m 1, N/2 and g 1, G, (4.34) where, HL,mg is the EPSSN leakage plus collision operator, q ,m,gg is the ingroup scattering source, and q' is a fixed source term that includes scattering transfers from energy groups other than g, the external source and fission sources. The iterative method begins by assuming a flux guess; then, Eq. 4.34 is solved for / '+1 and the ingroup scattering source is updated. This process continues until a certain convergence criterion is satisfied. The convergence rate of any iterative method is characterized by the spectral radius. For the source iteration method, in an infinite homogeneous medium, it is well known that the spectral radius is equal to the scattering ratio (c), given by c, = gC g (4.35) 'tg The scattering ratio is bounded between 0 and 1; hence, a cratio approaching one means that the problem will converge slowly, while, oppositely, a cratio close to zero, indicates a fast converging problem. Fourier analysis is the tool of choice to analyze the convergence behavior of iterative methods. For simplicity, I will consider the 1D EPSSN equations with isotropic scattering and source, given by a fU2 a S(X) /2 (x, ) + U (x)Ql'1 2(x, ) = co (x)ol(x) + Sex(x), (4.36) ox o (x) ox which following division by ua (x), reduces to a f2+a (x, S) + 1 ) = c S(x)+ x(x) (4.37) x C 2(x) x / + V+ o 2 (x) where c is the scattering ratio defined by Eq. 4.35. The EPSSN equations can be rewritten in terms of the error between two consecutive iterations as a /12 aU E 1 /2 (.,O) (4.38) t (x (x) 2(xx 2 where E/2 ( P) = IY /z(X, P) 1/ 2(x, P), (4.39) o,1 (x) = 0 (x) 0"_,1 (x). (4.40) The error terms are then expanded in terms of the Fourier modes, considering an infinite homogeneous medium, the Fourier ansatz is defined as follows 11/2 (, p) = f(p) exp(iAx), and 5o,1 (x) = exp(iix), (4.41) where i = and A ( o, o). By substituting the above relations into Eq. 4.38, I obtain the EPSSN equations mapped onto the frequency domain, resulting in function f(/) given by f(u) = .(4.42) 1+ p Therefore, the spectrum of eigenvalues is obtained by observing that the error in the scalar flux at iteration 1+1 can be written as so, (x) = d 1 i/2(x, ) =exp() exdP ) dJ (/A) =exp(iAx)Jd (4.43) 1+ P By performing the integration over the angular variable in Eq. 4.43, the spectrum of eigenvalues is found to be equal to arctan  W(A) = arctan (4.44) Ut The result obtained in Eq. 4.44 is similar to what is obtained for the SN equations. The spectral radius is found to be equal top = max[)(A)] = c However, the convergence behavior of the EPSSN equations is also affected by the value of the total scattering crosssection. For optically thin media, where the total crosssection assumes small values, Eq. 4.44 suggests that the convergence should be very fast, and in the limit as ot > 0, the spectral radius will tend to zero. 4.2.3 A New Formulation of the EPSSN Equations for Improving the Convergence Rate of the Source Iteration Method As discussed in the previous paragraph, the performance of the source iteration method applied to the EPSSN equations is similar to the SN equations. However, I have derived a new formulation of the EP SSN equations which reduces the spectral radius for the source iteration method. Appendix B addresses the performance of the new formulation for a criticality eigenvalue benchmark problem; note that the new formulation is a key aspect for the successful implementation of an acceleration method for the SN equations. The main idea behind the new formulation is to remove the ingroup component of the scattering kernel for each direction. In order to reformulate the EP SSN equations, we note that the evenmoments in the ingroup portion of the scattering kernel can be expanded as follows L1 L1 N/2 1 (2/+ 1)P, (Im )a ,g (F)qg (F) = (2/ + I)P1 (um )Ci,^ (F) ] WmP, (/m ) ,g (F). 1=0,2 1=0,2 m=l even even (4.45) The term V 'g (F) is consistently removed from the ingroup portion of the scattering kernel and from the collision term on the lefthand side of the EP SSN equations. 2 L1 N/2 V. / VV/,() + ,(F)R (F)= 1 (2/+l) I(/(m)cl,gg (F) 1 Wm,P, (m )Vf,g (F) tg ) ) 1=0,2 m'=1 even m' m G L1 G L + 0 (21+ I)P, (km,) g'g ()1,g ( m 1 (2/ + )P, (m ),g'>g 1,g g'=l 1=0,2 ut,g g'=11=1,3 g' even odd L1 m G L + 1(21+1)P (,m)Sg (F). 1 (2/+l)P(),g() +f,g() 1=0,2tg (F '11,3 even odd for m 1, N/2 and g 1, G. (4.46) Note that in Eq. 4.46, the total crosssection is replaced, in the collision term, with a directiondependent removal crosssection as follows L1 ,g (F)= ctg (F) (2/ + 1)P (/km)2 mQ ,gg(F). (4.47) 1=0,2 even In this new formulation, the main idea is to remove a "degree of freedom" from the iteration process in order to reduce the iterations on the component ,g (rF) This modification leads to a drastic reduction of the spectral radius. 4.3 Comparison of the P1 Spherical Harmonics and SP1 Equations In order to understand the assumptions on which the SPN and the EPSSN equations are based, it is useful to examine the 3D P1 spherical harmonics equations.6 The expansion in spherical harmonics of the angular flux can be written as follows V(~, )= 1 ( (2/ +1)Pm (cos )[ 1m (F) cos(m)o) + /y, (F) sin(mp)l, (4.48) 1=0 m=0 where 0 < 9 < zn andO < ( < 2i . In the following discussion, for simplicity I will assume a Pi expansion of the angular flux in spherical harmonics, given by , ) =PO (cos 9)o (F) + 3P0 (cos 9)/0oF (F) + 3P' (cos 9)[i/f, () cos ( + 72,(F) sin P]. (4.49) By substituting the definitions of the Associated Legendre polynomials6 in Eq. 4.49, I obtain the Pi expansion for the angular flux: ,(, ) = Voo (F) + 3p/o (F) 3 sin 9[,, (F) cos o + Yii (F) sin ], (4.50) where u = cos 9. The derivation of the SPN equations outlined by Gelbard, assumes implicitly that the angular flux be azimuthally independent, and hence symmetric with respect to the azimuthal variable. By introducing this assumption on the P1 expansion of the angular flux in Eq. 4.50, I obtain (r, )= ) /, fdp = t oo (F) + 3V'o (F) 3 sin .9 [y9 (F)cos p + 7, (F)sin (pip. (4.51) Therefore, by performing the integration on Eq. 4.51, I obtain (F,/ ) = 00 (F) + 3/ /0 (F). (4.52) It is evident that the angular flux obtained in Eq. 4.52 is equivalent to the SPi angular flux where, y,'o is the scalar flux and y/10 is the total current. The general formulation of the multigroup PN equations6, with anisotropic scattering and source, is obtained by substituting Eq. 4.48 into the linear Boltzmann equation and deriving a set of coupled partial differential equations for the moments , (F) and 7 (F). /I I 2(1 + m + 1) + 2(1 m) + Oz Oz Ox +(l + m +2)(l + m + 1V) ' m ax + 2(21+ 1)lo,gV I (1I 1 + Yg+l\m 9y 1y V'l + 0) 1g+l m )( lm) l 1m + c 'x 2Slm,g , 2(1+ m + 1) +l +2(1 Oz + (1 + m + 2)(1 + m + 1) +2(21+ 1)ol,g = 2S' nm) IlIm Oz 'Ylg+ l+1 + ay + Ov 0/g+ lm a IIm + aym M1 ax ay (1 m 1)( m) g g / O lm 1 ay (4.53b) forg= 1, G, where ,g = tg sl,gg * Therefore, the P1 equations are obtained by evaluating Eqs. 4.53a and 4.53b for /=0, 1 and m=0, 1, as follows (1=0, m=0) 2 8z 2 c " z +2K + 2 Ox ay ar1" dr +6a +2 x + 2 Vf & dy2110 2Soo g, (4.54a) (4.54b) + + 2o,g /0 cYx 6"21 910 S+ 67 VI 0 2S1o g 9<;g (4.54c) (4.53a) + aiy } 1 a (1 1, m=0) k 4 2 Oz 4 20 +2 020 + 6 8z 8z ay 21 + 601,Y = 2S'o g. I y O a y ax)J 20g V20 ay +0 +12K 9y ++ 60,g vl11 x a, ) /20 12 22 a 22 + 6 71,gY 2SIg, (4.54e) 2S',g. (4.54f) The terms with /1> and m> are dropped from Eqs. 4.54c through f, yielding the following relationships 1 1 301,g Og0 S10,gO az 30g 1, g 1 70 S',g 710 + 3017 az 307g 01 0g+ 700 S g /g 30 lg a Oy .3.g 30 301g (4.55a) (4.55b) (4.55c) 1 1 17 o  30,I Sy Y0go + S 0'g ax J 301 Then, by substituting Eqs. 4.55a, c and d in Eq. 4.54a, I obtain ,1 a 1 30g OX300 + 30,7, gx 3(7,g c 1 y 301,g Oy~ + O,gy00 = S00, + Sl, 8x where a S, S~S al (S' S 3 0 ,g x3 g aJ g31, J 1,g Z 30,~, 8x (30,~, B Analogously, by using Eqs. 4.55b, c and d in Eq. 4.54b, I obtain (=1, m=1) 6 a~ f21 2 a az y (4.54d) 6 21g 8z (4.55d) (4.56) (4.56a) 1 1 8' 8 1 8gf . V0 00 + 0+ og + g00 SOg + Sg (4.57) 3cg cx 3c g y 3cl,g ax where SSI I'og SI'g + Slg (4.57a) Ozg az 30,,g ) 33,g) y 3,g" Eqs. 4.56 and 4.57 constitute a coupled system of partial differential equations for Vog and yog, which must be solved iteratively. Recall that the assumption made in the SPN methodology is that the angular flux is azimuthally symmetric; therefore, to obtain the SPI equations (Eq. 4.58 or 4.59), terms such as yg are dropped from Eqs. 4.56 and 4.57, as follows V. g v +C00 =SOV. g (4.58) 3cl, g 0 3c lg) or V. + ,gq0 = S, V. g (4.59) 3cg Ig 3c g ) Here, I can also conclude that in the case of a homogeneous medium, with isotropic scattering, the P1 and the SPi equations yield the same solution, because the azimuthal dependency on the angular flux is removed. Note that this result can also be generalized to the SPN equations. CHAPTER 5 NUMERICAL METHODS FOR SOLVING THE EPSSN EQUATIONS This chapter addresses the numerical techniques utilized to solve the EPSSN equations; I will describe the discretization of the EPSSN equations in a 3D Cartesian geometry using the finitevolume method, along with the matrix operator formulation utilized and the boundary conditions. I will also introduce the Compressed Diagonal Storage (CDS) method, which is fundamental for reducing the memory requirements and the computational complexity of the iterative solvers. Further, a new coarse mesh based projection algorithm for elliptictype partial differential equations will be presented. Finally, I will describe a class of iterative solvers based on the Krylov subspace methods, such as the Conjugate Gradient (CG) and the BiConjugate Gradient methods (BiCG). The CG and BiCG methods have been implemented to solve the linear systems of equations arising from the finitevolume discretization of the EPSSN equations. Furthermore, the issue of preconditioning of the CG methodology will be discussed. 5.1 Discretization of the EPSSN Equations Using the FiniteVolume Method The EPSSN equations derived in Chapter 4 are discretized using the finitevolume approach. For this purpose, I consider a general volume Vin a 3D Cartesian geometry. The volume Vis then partitioned into nonoverlapping subdomains Vj, called coarse meshes. Note that, the coarse mesh subdomains are generally defined along the boundaries of material regions. As I will discuss in Chapter 7, the main purpose of this approach is to partition the problem for parallel processing. The discretization of the spatial domain is completed by defining a finemesh grid onto each coarse mesh. I have derived a formulation of the discretized EPSSN equations which allows for variable fine mesh density on different regions of the problem; this approach is very effective to generate an effective mesh distribution, because it allows a finer refinement only in those regions where higher accuracy is needed. The finitevolume discretization of the multigroup EPSSN equations (Eqs. 4.29) is obtained by performing a triple integration on a finite volume, dr dxdydz, as follows I , ()dr = [Q ,m()+Qextg, m(F)Q Pr, (5.1) v vtgQ where G L1 Q,,,m (F)= (2 1+ )P, (/m ,),,,g',, (F) ,,g' (F) g'= 1 0,2 even (5.1a) m Z (21+1)P, (m)csl.g' >g (f),g' (g) tg ( ) g'=1l=1,3 odd L1 L Qextgm(F)= (21+ )P (/,m)Sg() 'm (21+l)(P, um),g, (F), (5.1b) 1=0,2 ut,g (r)=1,3 even odd and Qfg () = VCf, ()0 (F) (5.1c) For this derivation, I consider a central finitedifference scheme for generic mesh element with coordinates x,, y, and Zk; an example of a fine mesh element and its neighbor points is shown in Figure 5.1. (i, j, k+1) S (i,j+1, k) (i1,j,k) (i, j, k (i+1,j,k) z (i,j1, k) (i,j, k1) x Figure 5.1. Fine mesh representation on a 3D Cartesian grid. The generic fine mesh element is defined by the discretization step sizes, Axc, Ayc, and Az, along the x, y and zaxis, respectively. Note that the discretization steps are constant within each coarse mesh; hence, a nonuniform mesh distribution is not allowed. The discretization steps are defined as follows Lc Lc Lc Ax = Ay = Az = and Avc = AxcAyAzc, (5.2) Nc' N y Nc force 1, Ncm where, Ncm is the total number of coarse meshes; L L and Lc are the dimensions of the coarse mesh (c), along the x, y and zaxis, respectively; and N Ny, and NC refer to the number of fine meshes along the x, y and zaxis, respectively. Note that, Eq. 5.1 is numerically integrated on a generic finite volume Ave. I will first consider the integration of the elliptic or leakage operator (first term in Eq. 5.1) as follows 68 J,2 r1: dx/ l/ dy k+1:d z: J ) \ (r) S2 neihbo p Yins lon t fh k Vxa x IL km 2 8x / 2 L Kd1/ ^ d12 1 a iV wm' (r), = L[A Vtg(rg) ) gg(F) Vg ()JJ 11 2 1 + 2 7 A X CL (rF) Oyz V1M' g 2F(kY 1) Jm ,_12 For simplicity, I will derive the discretized operator along the xaxis; the treatment is analogous along they and zaxis. Figure 5.2 represents the view of a fine mesh and its neighbor points along the xaxis. Ax, Ax,/2 Ax,/2 X;1 X;1/2 X X + 1/2 X,+1 x Figure 5.2. View of a fine mesh along the xaxis. In Figure 5.2, au represents a generic macroscopic crosssection (e.g., total, fission, etc.) which is constant within the fine mesh. In Eq. 5.3, I evaluate the rightside and leftside partial derivatives along the xaxis at x,+112. fE()2 E( )(Xl/ g 2,y, Zk) mg( ,yj,,zk) Jtg(XI,Y,,Zk) )j/2 fE2( ) (X m mgl fJm,g )mg/22Z ;,Zk)m= ,g Z(5.4) '(E+) +2 ZVk I m,g( +1 ,YJ,Zk) ~m', g)( +1/2,Yj,Zk) 1'g (,+l/2,yj,Zk)= mtg(Xl YJIZk) (5.5) t,g (X+1 ,, Zk) x, c/2 In order for the elliptic operator to be defined, the function ,/ g (x, y, z) must be continuous along with its first derivative f, (x, y, z) and second derivative, which translates into the fact that the evenparity angular flux belongs to a C2 functional space, or rg (x, y, z) e C2. Therefore, the following relationships hold true Vg (X,+l/2 ,YjZk ) mg )(X/29 1,Yj,Zk)= m ,g+(X 1/2,YJ,Zk), (5.6) and fg(X+/2 1 YZk) fmg )(X1/2,Y,Zk)= (1/2 ,Zk). (5.7) Therefore, I eliminate the value ofV/lg ( +1/2, Y, Z) in Eqs. 5.4, obtaining the second order, centralfinite differencing formula for the evenparity angular flux: dx E dx E E d+l1,],kmg/+1,jkmg + ,djkmg/,Jkmg +l/2,,j,m,g (5.8) d I+l,j,k,m,g d ,j,k,m,g and the evenparity current density z+l/2,j,k,m,g dx x +l,j,k,m,g V ,jk,m,g 5.9) dj,k,m,g +d+l,j,k,m,g In Eqs. 5.8 and 5.9, I have defined the pseudodiffusion coefficients along the x axis, as 2 2 2 dx m dx km dx m (5.10) i,],k,m,g z ,j,k,m,g i l,]j,k,m,g 1 7t,c,j,k,g ct, +l,j,k,g c (t,,1,,k,g c Analogously, the expression for fE (x 1/2, k) is obtained as follows E 2dX ,k,m,gd ,kg ( E E 1/2,j,k,mg d + x (f ,],k,m,g i1,j,k,m,g) (5.11) Ij,k,m,g d l,j,k,m,g The partial derivatives along the y and zaxis are discretized in a similar fashion, yielding the finitevolume discretized elliptic operator given by i d xJy Jl/Zdy J;k+lPdz K Q. 12m V 1,2( 1 Xmg 2 k2dj j Z 12m lk E g S1/2 i 2 Y k1/2 d g () , ,E H H, H. A Az, + [ ,m,g E 1,,k,n,g l ,k,nm,g ,,n,g ( ,k,,g k g (512) 2dJk m,gdx1J,k,m g 2dJ,k,m gd 1,],k,m g Y = d1 +d ,ak 1mg dx d, (5.13a) , Ax A [r j,k m,g +1,jkm,g kj,k,m,g ( 1,j l,,k,gg where 2d dx 2d d2xk d,j,k,m,g +l,J ,k,m,g ,j,k ,m,g ,k,m,g 1 ) ]'i,+l,m,g d +d l 'd'lmg +d m (5.13) ,jikmg +l,j,k,m,g i,j,k,m,g Il,j ,k,m,g 2dy dY 2cY dY S i, ,k,m,g i,j+l,k,m,g i, ,k,m,g d ,jI,k,m,g i,],k,m,g j+l,,,g i,],k,m,g Ij ij,k,m,g 2Ykdk,m,g z ,k, z ,Yk,k,m,g ,k, ,gd (5.13 i,,j,k,m,g + d,j,k+l,m,g i d j,k,m,g d, ,,k1,m,g and 2 2 2 S/m /m jf /m54 dy dy = dd = (5.14) ,,t,,j,k,g lc ,t,,j+l1,k,g (7c t,i,J1,k,g y c 2 2 2 dz = mI z = (5.15) I,j,k,m,g n di,j,k+l,m,g d ,j,k,m,g n (5.15) tl,],k,g z c t,l,],k+l,g c t,l,],k1,g Z c Finally, by integrating the remaining terms of the EPSSN equations, I obtain the complete multigroup EPSSN formulation with anisotropic scattering and source as follows  Ay Az,[a,+,m, l,,k,m,g V,k,m,g ,zi,m,g ( i,k,m,g  Ax cAz [,l+l,m,g Vm(V,J +l,k,m,g V ijk,m,g ) ,jl,m,g ,(Vi,k,m,g  AX Ay 1,/k ,g jk+l,m,g I,],k,m,g )Yk,k,m,g ( Ijk,m,g G Li + t,gj kV ,,k,m,g c = Av (2/+1)P ( )g,g'>g,,j,k,g',,j g'= 1=0,2 even AG L Ay Az, m Z(2/+l)P(m)csl,g'>g,,,k (1,g',+l/i2,,k 't,g,i,j,k g'= 1=1,3 odd G L Ax Az, m Z(21+ l)P, (Um)gsl,g'>g,,,,k (1,g',,j+l/2,k 't,g,z,j,k g'= 1=1,3 odd G L SAx c (21+y )PI (Am)Crsl,g'> g',,.klg,,k+l/ 2 t,g,z,j,k g'= 1=1,3 odd L1 + 1(21+1)P (/m)S,g,,,j,kAv + 1 0,2 even G L Ay c m i:(21+1)1P(,m)sl,g'g,,J,k (Sl,g',+l/2,J,k t,g,z,j,k g'= 1=1,3 odd G L AxcAz, Z :(2/+ l)P(m)'sl,g'>g,,,k (,gj+ /2,k 0 t,g,z,],k g'= 1=1,3 odd Axc, m Z (2+l)P (P sl,g'g,,,,k (S,g' +/2 0t,g,z,],k g'= 1=1,3 odd Qf,g,,j,k Av , V1,j ,k,m,g ,Av  k c l,g',1/2,],k) l,g',z,j1/2,k ,g',,,k1/2 ) S I,g',1/2,jk Sl,g',,j 1/2,k ) Sl,g',,j,k 1/2)+ for c=1, Ncm, m=1, N/2, L=0, N1, g=1, G. The EPSSN equations discretized with the finitevolume method can be expressed in a matrix operator form characterized by a 7diagonal banded structure. (5.16) 72 Dm g U Ug Um g LX D Ux UY U m,g m,g mg mg mg Sm,g m,g mg mg m,g g m,g g m,g L LX D Ux m g mg nm,g mg Lz L Lx D Ux Lg L LX D mg mg m,g m,g force 1, Nc,; m 1, N/2; g 1, G, where Dx = AycAzc(azn, +a1 ,,n, )+AxcZ c j,nim,g + ~,~1,g +AcAyc c(k,k+1mg + k ,kim,g Ug c =AycAzca,+,, LX = AycAzca, mg mg mg C mg UY, =AxcAzc ,,,,g LY, = Ax AzyP , g = c ycYk,k+l,m,g, g g = AcycYk,k1,m,g for i=2, N 1, j=2, Ny , k=2, N 1, c=1, Ncm, m=1, N/2, g=1, G. 5.2 Numerical Treatment of the Boundary Conditions The boundary conditions for the EPSSN equations are discretized as well using the finitevolume method. In general, the BCs can be prescribed at back (Xb), front (+xb), left (yb), right (+yb), bottom (Zb), and top (+Zb). The reflective boundary conditions are simply derived as follows: Xb) f,g,2,,k = 0 +Xb) ',g,N 1 2,,k =0, (5.17a) Yb) ,g,, /2,k = 0 +yb) m,g,,N +1 = 0, (5.17b) +Zb) )i,,2,,, N+1+2 = 0 . (5.17c) Zb) V,g,,,12 = 0, The vacuum boundary conditions are obtained from Eq. 4.32, by setting a = 0. Hence, the vacuum boundary conditions along the x, y and zaxis are given below: Front side vacuum boundary condition x = +Xb o aN,m,g E VN,+1/2,],k,m,g Nx ,j,k,m,g S+aN,,m,g aNmg 1 (2/+1)Pl (m ))[s,g'>g,Njk g,N,+1/2,j,k +Sl,g',NJ+1/2,j,k +aN,,m,g t,g,Nx,j,k g'=1=1,3 odd (5.18a) Back side vacuum boundary condition x = Xb o a,m,g E V1/2,j,k,m,g Vl ,],k,m,g + 1 l,m,g aMg G L (5.18b) + al,m,g 1 (21+1) P( Um)[ sl,g'>g,l,J,k I ,g',3/2,j,k +Sl,g',3/2,j,k 1+al,m,g (t,gl,j,k g'=11=1,3 odd Right side vacuum boundary condition y = +Yb 0 bN,,m,g E 1,Ny+l/2,k,m,g I+b ,m,g Ny,k,m,g bNymm,,g 1b G +_: (21+1)PI(um) sl,g'>g,iNyk 1,g',iNy+1/2,k +SI,g',iN +1/2,k + Ny,m,g 't,g,l,Ny,k g'=l1=1,3 odd (5.19a) Left side vacuum boundary condition y = yb o bmg E Vz,l1/2,k,m,g 1 +b 1,l,k,m,g l,m,g +,m,g 1 ,(2/1+l)P(/m)[sg'>g,,l,kO,g',i,3/2,k +SI,g',,3/2,kl +b 1,m,g t,g,i,l,k g'=11=1,3 odd 74 Bottom side vacuum boundary condition z = Zb o CN,m,g E Vl,J,N,+1 2,m,g = ,],N ,m,g S NN,mn,g Nmg 1 (2 +1 )P,) ( o12 +S, 121 + N,m,g gt,g,i,j,N, g'= =13 odd (5.20a) Top side vacuum boundary condition z = +Z O Cl,m,g E i,],j1 2 ,m,g + g V ',],,,g Simg G L (5.20b) + ,m,g 1 Z Z (21+ 1)P,(, m)[r,g g,z,j,o,g',',,3/2 +S,,,,,3/2 1+ ',mg ut,g,c,j,1 g'=1=1,3 odd where 2 ,m 2 ,m 1,m,g A N,,m,g A , c c t,g,l,1,k A c t,gz,Ny,k 2,tm 2,tm Cimg CN,,m,g ,g c g, 1 zc t,g,jN and Nc Ncm Ncm N, = N ,, N, Y X: = Y Nc c=l c=l c=1l 5.3 The Compressed Diagonal Storage Method Due to the sparse structure of the matrices involved, I have adopted the Compressed Diagonal Storage (CDS) method in order to efficiently store the matrix operators. The CDS method stores only the nonzero elements of the coefficient matrix and it uses an auxiliary vector to identify the column position of each element. Due to the banded structure of the coefficients matrix, a mapping algorithm is easily defined for a generic square matrix as follows: a(i, j) e A,J a(i, d) E Ad i= 1,I i= 1, jcol(i,d). (5.21) j = 1, J d = 3,3 The algorithm defined in Eq. 5.21, maps the full structure of the matrix A into a compressed diagonal structure, where for each element on row i, there is an associated diagonal index ranging from 3 to 3, with index 0 being the main diagonal, and an auxiliary vectorjcol, which stores the column position of each element. If we consider a 360x360 full matrix in single precision, with a total of 129600 elements, the memory required for allocating the matrix is roughly 2.1 MB. However, if the CDS method is used, the total number of nonzero elements to be stored is only 2520, for a total memory requirement of 42 KB, which is a reduction of a factor of 50 compared to the full matrix storage. Moreover, since the CDS method stores only nonzero elements, I have also obtained a reduction in the number of operations involved in the matrixvector multiplication algorithms. 5.4 Coarse Mesh Interface Projection Algorithm The partitioning of the spatial domain into nonoverlapping coarse meshes leads to a situation in which the EPSSN equations have to be discretized independently for each coarse mesh. Therefore, each coarse mesh is considered as an independent transport problem; however, to obtain the solution on the whole domain, an interface projection algorithm has to be used in conjunction with an iterative method. The matrix operators have to be modified on the interfaces in order to couple the equations on each coarse mesh. For explanatory purposes, consider Figure 5.3, which shows the interface region between two coarse meshes. Coarse mesh 1 Coarse mesh 2 I I AxN A, XN XN1/2 XN X1' X3/2' X2' XN+1/2 X1/2' Figure 5.3. Representation of a coarse mesh interface The coordinates xN+1/2 and xl/2, represent the interface on coarse mesh 1 and 2, respectively. As shown in Figure 5.3, the discretization of the elliptic operator for coarse mesh 1, using the central finite difference method, would require the values of the even parity angular flux at points XN1, xN, and xlr. Similarly, in coarse mesh 2, the discretization would involve the value of the EP angular flux at points xN, xr,, and x2'. However, the point xr, is located on coarse mesh 2 and point xN is located on coarse mesh 1; hence this term does not appear explicitly in the matrix operator for both coarse meshes. In order to couple the equations on the interface, I have reformulated the discretized equations by bringing the unknown points on the right side of the equations. The numerical discretization of the EPSSN equations in coarse mesh 1 would yield [ aNIg(m^g ,) NN ia(1,m,g 1mg im,g )]+ Et g, H, Ax( = AxQ Ax + ,g ENn,g N N= gN Next,g,mN NQf,g,N (5.22) where 2dx dlmg 2dx mgdXlmg (.3 aN,l',m,g 2d 'g and aNN lmg dXNg N,g (5.23) N,mg l',m,g N,mg N 1,m,g The coefficient d,,,,g depends on the material properties and fine mesh discretization of coarse mesh 2, and it is calculated a priori; however, in Eq. 5.22, the term /,,,g is unknown, and hence has to be evaluated iteratively by placing it in the source term, as shown in Eq. 5.24. E E E E A, aN,l,m,g N,m,g + aN,N,m,g N,m,g N,N,m,g N 1,m,g t,g,N Nmg NQs,g,,N +N XNQext,g,mN + Nf,g,N+ N,1', mg g * (5.24) A similar equation can be formulated for coarse mesh 2, as follows a,2' ( g i 1,N,m,g (,,g lm,g)] (5.25) + ag ,vj Rg1x = AQs,g,m,1 + extg,mi, + ZhlQf g,1, or Eaig^/mg +ai g g + a,N,m,g Vmg + g ,',m,g =  l',2',m,g f2',m,g ',2',m,g l',m,g V l',Hmg(l',m) (5.26) lQs,g,m,l' + AX'Qext,gml' + a 'Qf,g, + al',N,m,g N,m,g, where 2dX dx 2dX dx aI2 g 2'"mg lm and a Nmg (5.27) 1',2,mg d +dx d d 1',m,g 2',m,g l',m,g Nmg Therefore, Eq. 5.24 and 5.26 are coupled through the value of the EP angular fluxes ',m,g and N,m,g The EPSSN equations are solved iteratively starting in coarse mesh 1, and assuming an initial guess for Vim,g. Once the calculation is completed the value of N,,g in Eq. 5.26, is set equal to /N,g Hence, once the calculation is completed on coarse mesh 2, the value obtained for /1 mg is used in Eq. 5.24, to update the value of mg; this procedure continues until a convergence criterion is satisfied. In a 3D Cartesian geometry the coupling on the coarse mesh interfaces is achieved exactly as described above; however, in this case the coarse meshes can be discretized with different fine mesh grid densities. The variable grid density requires a projection algorithm in order to map the EP angular fluxes and the pseudodiffusion coefficients among different grids. As stated earlier in this chapter, the variable density grid approach is very effective to refine only those regions of the model where a higher accuracy is needed; note that the main constraint on the fine mesh grid is the mesh size being smaller than the mean free path for that particular material region. The main philosophy behind the projection algorithm is derived from the multigrid method, where a prolongation/injection operator is used to map a vector onto grids with different discretizations. Figure 54 shows the application of the projection algorithm along theyaxis on the interface between two coarse meshes. Coarse mesh 1 (Finer grid) Coarse mesh 2 (Coarser grid) Grid 1 Grid 2 F3 F4 G5F G6F G7F G8F F iF,t F. G1F G2F G3F G4F x x Figure 5.4. Representation of the interface projection algorithm between two coarse meshes. For simplicity, I will consider the projection of a vector between two coarse meshes, along the yaxis, as shown in Figure 5.4. The finetocoarse projection of a vector is obtained by collapsing the values as follows 4 FIc Z W'FGQF, 1=1 wF=A for i=l, 4 Ac (5.28) (5.29) In Eq. 5.29, AF and Ac, are the areas associated with the finemesh and coarse mesh grid, respectively. Conversely, the coarsetofine projection is obtained as follows G1F = wIFF, (5.30a) G2F = W2FF (5.30b) where G3F 3FFc, (5.30c) G4F = 4FFC (5.30d) In general, the finetocoarse mesh projection is obtained with the following formulation NF Fc =Z WFGJF, (5.31) J=1 where AJF WF A (5.32) AC The weights in Eq. 5.32 are the ratios of the areas of the fine meshes intercepted by the coarse meshes on which the values are being mapped. Similarly, the coarsetofine mesh projection algorithm is defined as follows N, GF JC FJc, (5.33) j=1 where AF wc (5.34) A,, By using the above formulations, the evenparity angular fluxes and the pseudo diffusion coefficients are projected among coarse meshes with different grid densities. Note that the projected pseudodiffusion coefficients need to be calculated only one time at the beginning of the calculation, while, the projections for the EP angular fluxes have to be updated at every iteration. 5.5 Krylov Subspace Iterative Solvers Due to the size and sparse structure of the matrix operators obtained from the discretization of the EPSSN equations, direct solution methods such as LU decomposition and Gaussian elimination do not perform effectively both in terms of computation time and memory requirements. In contrast, the Krylov subspace iterative methods, such as Conjugate Gradient (CG), are specifically designed to efficiently solve large linear systems of equations characterized by sparse matrix operators. Note that in many engineering applications, the matrix operators resulting from a finitedifference discretization is usually positivedefinite and diagonally dominant. These conditions are fundamental in ensuring the existence of a unique solution. A matrix is positivedefinite if it satisfies the following condition X' Ax > 0, for every vector # 0. (5.35) Moreover, a matrix is defined to be diagonally dominant if the following condition holds true. a,, > a, for i 1, n. (5.36) J1 The CG algorithm is based on the fact that the solution of the linear system Ax = b is equivalent to finding the minimum of a quadratic form given by f(2)= 1 T A b T+c. (5.37) 2 The minimum of the quadratic form of Eq. 5.37 is evaluated by calculating its gradient as follows 8x, f '(x) = (5.38) The gradient of a function is a vector field, and for a given point x, points in the direction of the greatest increase of f(i)). Because the matrix A is positivedefinite, the surface defined by the function f(x) presents a paraboloid shape, which ensures the existence of a global minimum. Moreover, the diagonal dominance of the matrix A ensures the existence of a unique solution. By applying Eqs. 5.37 and Eq. 5.38, we derive the formulation for the gradient of the function f(2), given by f'(2) = AI + A A b. (5.39) 2 2 If the matrix A is symmetric, Eq. 5.39 reduces to f'(x) = Ax b. (5.40) Therefore, by setting f' (2) in Eq. 5.40 equal to zero, we find the initial problem that we wish to solve. 5.5.1 The Conjugate Gradient (CG) Method The CG method is based on finding the minimum of the function f(x) using a line search method. The calculation begins by guessing a first set of search directions do using the residual as follows: d, = :0 = b Axo. (5.41) The multiplier a for the search directions is calculated as follows FTF a, (5.42) d, Ad where i is the iteration index. The multiplier a is chosen such that the function f(x) is minimized along the search direction. Therefore, the solution and the residuals are updated using Eqs. 5.43 and 5.44. x+1 = X, +ad, (5.43) ft+ = a,Ad,. (5.44) The GramSchmidt orthogonalization method is used to update the search directions by requiring the residuals to be orthogonal at two consecutive iterations. The orthogonalization method consists of calculating the search directions S= +, (5.45) where the coefficients / are given by T + 1,+1 (5.46) r,r, Note that Eq. 5.44 indicates that the new residuals are a linear combination of the residual at the previous iteration and Ad It follows that the new search directions are produced by a successive application of the matrix operator A on the directions at a previous iteration d,. The successive application of the matrix operator A on the search directions d, generates a vector space called Krylov subspace, represented by K, = span ..4 .. .4 d 0,..., A''do }. (5.47) This iterative procedure is terminated when the residuals satisfy the following convergence criterion MAXr )< E, (5.48) where E is the value of the tolerance, which is usually set to 1.0e6. 5.5.2 The BiConjugate Gradient Method The BiConjugate Gradient (BiCG) has been developed for solving nonsymmetric linear systems. The update relations for the residuals are similar to the CG method; 