N O D A L M E T H O D S : A N A L Y S I S , P E R F O R M A N C E A N D F A S T I T E R A T I V E S O L V E R S By J. David Moulton B. Eng. (Physics) McMaster University, 1988 M . Eng. (Physics) McMaster University, 1990 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF T H E REQUIREMENTS FOR T H E D E G R E E OF D O C T O R OF PHILOSOPHY in T H E FACULTY OF GRADUATE STUDIES INSTITUTE OF APPLIED MATHEMATICS DEPARTMENT OF MATHEMATICS We accept this thesis as conforming to the required standard T H E UNIVERSITY OF BRITISH COLUMBIA November 1996 Â© J . David Moulton, 1996 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of Br i t i sh Columbia , I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Mathematics The University of Br i t i sh Columbia 2075 Wesbrook Place Vancouver, Canada V 6 T 1Z1 Date: A b s t r a c t Nodal Methods have long been one of the most popular discretization techniques employed with in the reactor physics community, while remaining conspicuously absent from the mainstream numerical analysis literature. A fundamental reason for this anomaly is that the physical arguments which were used to develop and enhance these methods seemed at odds with more rigorous discretization techniques. To facilitate communication be-tween these distinct communities, a detailed chronological study of the lowest-order nodal methods for linear second order elliptic problems is presented. The presentation high-lights the underlying motivation of these methods and formalizes many of their renowned properties. In addition, various equivalence relations within this family of discretizations are demonstrated, and equivalences with specific non-conforming and mixed-hybrid finite element methods ( F E M s ) are established. Rigorous error bounds and stability properties follow immediately from these latter equivalence relations, corroborating the results of a more rudimentary truncation error analysis. A n inherent difficulty of reactor simulation is that the coefficients i n the mathematical model exhibit severe variations on two significantly different length scales. A s i n many other applications this is treated by defining an appropriate homogenization procedure which yields a simplified model with piecewise constant coefficients on a coarse scale suitable for efficient computation. Significant enhancements in accuracy are possible if the processes of homogenization and discretization are unified. We review the popular techniques that are based on this premise and rely on certain properties of the nodal discretization. In addition, we address the factors that contribute to their success i n reactor modelling and deter their generalization outside of the reactor physics community. i i A s an alternative to these highly specialized methods, we introduce a new multi- level homogenization technique which is readily applicable i n a general setting, and is shown to have many important attributes. Widespread acceptance of nodal methods has also been hindered by their use of non-standard unknowns, as this results in stencils that appear awkward and incompatible wi th sophisticated iterative solution techniques. Specifically, equivalence wi th certain mixed-hybrid F E M s reveals that the nodal discretizations result i n an indefinite system, which in two dimensions contains both cell-based and edge-based unknowns. Yet , in -herent in this structure is a natural partitioning of the system which may be exploited to define a hierarchy of reduced systems (i.e. Schur complements) that are symmetric positive definite. Unfortunately, the reduced systems that involve unknowns of only one type, and hence seem most compatible with sophisticated iterative methods, suffer a loss of sparsity. However, the structure inherent in this hierarchy may be uti l ized i n the con-struction of sparse approximate Schur complements for these systems. It is demonstrated that any one of these approximate operators, which are of either the standard 5 or 9-point family, may be uti l ized as an excellent preconditioner for conjugate gradient iterations. The efficiency of this approach is fully realized when the preconditioner is approximately inverted using only a single V - or W-cycle of a robust Black Box mult igr id solver. iit Table of Contents Abstract i i List of Tables ix List of Figures xi Acknowledgment xi i i 1 Introduction 1 1.1 Background Sketch 2 1.2 Outline 5 2 Background Theory 9 2.1 Neutron Transport Theory 9 2.1.1 The Neutron Transport Equation 9 2.1.2 Simplifying Assumptions 13 2.1.3 Transport Boundary Conditions . . . 14 2.2 Diffusion Theory . , 15 2.2.1 The Pi Approximat ion 15 2.2.2 Cr i t ica l i ty : The Reactor Eigenvalue 18 2.2.3 Approximate Boundary Conditions 19 3 Nodal Discretizations 21 3.1 Classical Approaches 23 iv 3.1.1 Fini te Volume and Nodal Balance 23 3.1.2 Standard Cel l Centered Finite Differences 25 3.2 Simulators 27 3.2.1 Fine Mesh Correction - Cross Coupling 29 3.2.2 F L A R E . 29 3.3 Transverse Integration 32 3.4 Polynomial Nodal Methods - N E M 33 3.4.1 Components of the N E M Discretization 34 A Polynomial Basis 34 Transverse Leakage Approximations 37 Weighted Residual Procedure 39 3.4.2 The Lowest-Order N E M 40 Interface Current Formulation 41 F l u x Formulation 43 Edge Based F l u x Formulation 44 3.5 Analy t i c Nodal Methods - N I M 45 3.5.1 The Nodal Integration Method 46 3.5.2 The constant-constant N I M for Conservative Diffusion 48 3.5.3 The constant-constant N I M for Nonconservative Diffusion 50 The F l u x Formulation . . 53 Edge Based F l u x Formulation 54 3.6 Truncation Error Analysis 55 3.6.1 Piecewise Constant Coefficients 57 3.6.2 Piecewise Smooth Coefficients 58 3.7 Numerical Experiments 60 v. 3.7.1 Example 1: A Simple Test 60 3.7.2 Example 2: The Iron-Water Benchmark 63 4 H o m o g e n i z a t i o n 69 4.1 Equivalence Techniques i n Reactor Model l ing 70 4.1.1 Coarse-Scale Properties 72 4.1.2 A One Dimensional Study 73 4.1.3 Mult idimensional Equivalence Techniques 76 F l u x Weighting 77 Generalized Equivalence Theory 78 4.2 Asymptot ic Homogenization Techniques 82 4.2.1 Mult idimensional Periodic Analysis 83 4.2.2 Mult idimensional Nearly Periodic Analysis 85 4.3 M u l t i g r i d Homogenization 86 4.3.1 Obtaining the Diffusion Tensor 87 4.3.2 A Numerical Homogenization Suite 90 Shape Dependence 91 Dependence on the Diffusivity 91 A Nonsymmetric Inhomogeneity 94 5 F i n i t e E l e m e n t M e t h o d s 95 5.1 P r i m a l Conforming F E M 96 5.1.1 Variational Formulation 96 5.1.2 Conforming Polynomial Basis 98 5.2 P r i m a l Nonconforming F E M 102 5.2.1 The Edge and Cel l Basis 102 v i One-Dimensional Basis 102 Two-Dimensional Basis 105 5.2.2 Polynomial Basis Functions 107 A One-Dimensional Study 107 Equivalence i n Two Dimensions 114 5.2.3 Analy t i c Basis Functions 116 A One-Dimensional Study â€¢ 116 Equivalences in Two Dimensions 121 5.3 M i x e d H y b r i d Fini te Element Methods 123 5.3.1 Variational Formulation 124 5.3.2 Polynomial Basis Functions 126 5.3.3 Analy t i c Basis Functions 128 6 R e d u c e d S y s t e m s a n d A p p r o x i m a t e Schur C o m p l e m e n t s 129 6.1 Lumping The Edge-Based Schur Complement 130 6.1.1 A Condit ion Bound 132 6.1.2 Bounding the elements of 136 6.2 A Two-Step Lumped Approximat ion 137 6.3 Lumping the Cell-Based Schur Complement 140 6.3.1 A Condit ion Bound 141 6.4 Approximate Mult i - level Inversion of the Preconditioner 144 6.5 Numerica l Experiments 145 6.5.1 A Toy Problem 146 6.5.2 Ground Water Flow 148 6.5.3 A Diffusive Checkerboard 151 vi i 7 C o n c l u s i o n s 154 7.1 Research Summary 154 7.2 Future Work 159 B i b l i o g r a p h y 163 vi i i L i s t of Tables 3.1 Convergence results for Example 1: the transverse averaged scalar flux on a 10 x 10 uniform mesh 61 3.2 Convergence results for Example 1 in the Loo-norm 62 3.3 Convergence results for Example 1: the cell averages 62 3.4 The diffusivity D, and the absorption cross-section E a , for Example 2 [25]. 64 3.5 Errors for Example 2 on a 10 x 10 computational grid 67 4.1 A comparison of the shape dependence of the diffusivity relative to the square (i.e. the percent relative difference - % R D ) is presented for the results of Bourgat [12] and the black box homogenization 92 5.1 Canonical Definitions in I D 103 5.2 Canonical Definitions i n 2D 105 6.1 Iteration counts for the preconditioners, S^u and for the toy problem with constant mesh spacing on the domain [0,1] x [0,1] 147 6.2 Iteration counts for the lumped preconditioner, <SM;U wi th constant mesh spacing on the domain [0,a] x [0,6]. {Note: a = 2/a,/3 = 2/6} 147 6.3 Iteration counts for the two-step lumped preconditioner, S^u<v wi th con-stant mesh spacing on the domain [0,a] x [0,6]. {Note: a = 2/a,/3 = 2/6} 147 6.4 Iteration counts for the cell-based preconditioner, S<f, wi th constant mesh spacing on the domain [0,a] x [0,6]. {Note: a = 2/a,/3 = 2/6} 148 6.5 Iteration Counts with single V ( l , l , l ) cycles 150 6.6 Iteration Timings, [s] for the flow problem with single V(l,1,1) cycles . . 150 6.7 Iteration counts for ")x = = 2/10 153 x L i s t of F i g u r e s 3.1 A n L x M tensor product mesh imposed on the domain 22 3.2 The location of the unknowns 25 3.3 Stencil weights of the flux formulation wi th D = 1 and E a = 0 44 3.4 The edge-based flux formulation wi th D = 1 and E a = 0 46 3.5 The transverse averaged scalar fluxes, {4>i+^,j-,4>i,j+^} obtained with the constant-constant N I M for this simple test case are shown in (a) while the corresponding node averages, </>tJ are plotted in (b) 62 3.6 Schematic of the iron-water benchmark problem [25] 64 3.7 The constant-constant N I M solution of Example 2 computed on a 10 x 10 grid: (a) transverse averaged scalar fluxes (b) cell averages . 67 3.8 Solution cross-sections for Example 2 are compared to a reference solution that was obtained using the finite volume discretization on a 640 x 640 uniform mesh: yj = 10.5 (a) ful l domain (b) zoom on layer 68 4.1 A compass based definition of an arbitrary 9-point stencil 87 4.2 Three inhomogeneities with the same area of 1/4, but different shape. . 92 4.3 A square inhomogeneity with diffusivity, A and an area of l/9th 93 4.4 A comparison of homogenized diffusivities 93 4.5 The homogenization of an L shaped inhomogeneity leads to a dense tensor. 94 5.1 The quadratic canonical cell and edge basis functions 109 6.1 Upper bound of the condition number, K (6.6) and of the decay rate A c o ( r ) .137 xi. 6.2 A schematic of Mose et al.'s [56] first example 149 6.3 Solution 4>(x,y) and streamlines of Mose et al.'s [56] first example (above) 150 6.4 A diffusive checkerboard configuration based on Dendy's [23] example. . . 152 6.5 The solution, <f)(x,y) of the diffusive checkerboard computed on a 24 x 24 uniform mesh 153 6.6 2D Exponential G r i d with the parameters 7^ = -jy â€” 2/10 153 7.1 Stream lines for the convection dominated flow of Smith and Hut ton [67]. 161 7.2 The convection dominated flow of [67] with Pe â€” [-D] - 1 = 500 computed using a constant-constant N I M : (a) cell averages (d) inflow and outflow profiles 161 xn A c k n o w l e d g m e n t I would like to express my sincere appreciation to my advisor, D r . U r i Ascher, for promoting and tolerating my diverse research interests. In addition I would like to thank h i m for his guidance, insight, patience and gregarious hospitality, al l of which were v i ta l factors in the completion of this thesis. For many enthusiastic, and often impromptu, discussions I convey my sincere grat-itude to D r . J i m Varah. I would also like to thank the other members of my advisory committee, Drs. Br ian Wetton and Michael Ward , both of whom have made positive contributions to this research. I have had the good fortune of visit ing Los Alamos several times over the past few years. I would like to express my gratitude to Dr . James H y m a n for his key role i n making my early visits possible and to acknowledge his incredibly contagious enthusiasm that continues to inspire my interest in homogenization well beyond its roots as a summer project. I also extend my thanks to Dr . J i m Morel for both my in i t ia l contact wi th nodal methods and the many encouraging and instructive discussions that we had during my latter visits. I am particularly grateful to D r . Joel Dendy for both his insightful observations regarding my research and his generous hospitality. To my parents and brother I convey my deepest thanks for their unconditional support and encouragement. The long and arduous journey that has been this thesis, could not have been completed without the continued support and encouragement of Tracey W a r d . I am especially indebted to her patience, sacrifice and understanding, which have been essential i n the endless final stage of this work. xiii Chapter 1 Introduction N o d a l m e t h o d s have l o n g been one of the m o s t p o p u l a r d i s c r e t i z a t i o n techniques e m -p l o y e d w i t h i n the reactor phys ics c o m m u n i t y to solve m u l t i - g r o u p d i f f u s i o n p r o b l e m s {e.g. [72, 34, 53, 68] }. A survey of these m e t h o d s c a n be f o u n d i n [52]. T h e i r success is a resul t of t h e i r e x c e p t i o n a l accuracy w h i c h m a y be a t t r i b u t e d to three d i s t i n c t aspects of the nodal ideology. Spec i f i ca l ly , a k i n to f in i te v o l u m e m e t h o d s , n o d a l m e t h o d s are p h y s i c a l l y m o t i v a t e d ce l l -based d i sc re t iza t ions t h a t , b y c o n s t r u c t i o n , r i g o r o u s l y enforce ce l l ba lance . T h e y u t i l i z e an i n t r i g u i n g choice of u n k n o w n s ( cons i s t ing , i n two d i m e n s i o n s for e x a m p l e , of ce l l a n d edge m o m e n t s ) w h i c h makes the n o d a l d i s c r e t i z a t i o n n a t u r a l l y c o m p a t i b l e w i t h the var ious h o m o g e n i z a t i o n techniques t h a t are c r u c i a l to reac tor m o d -e l l i n g . M o r e o v e r , the n e u t r o n current is o b t a i n e d a c c u r a t e l y a n d a u t o m a t i c a l l y f r o m the d i s c r e t i z a t i o n , thereby a v o i d i n g p r o b l e m a t i c f in i t e difference a p p r o x i m a t i o n s . D e s p i t e these d i s t i n g u i s h i n g character is t i cs a n d the appearance of d i f fus ive p h e n o m -e n a i n m a n y d i s c i p l i n e s , n o d a l m e t h o d s have rece ived o n l y l i m i t e d a t t e n t i o n o u t s i d e the reac tor phys ics c o m m u n i t y . T h i s a n o m a l y is due , i n p a r t , to the l i m i t e d ana lys i s of these d i s c r e t i z a t i o n s (see [39, 36, 25] ). In fact , the p h y s i c a l a r g u m e n t s w h i c h were used to deve lop a n d enhance these d i sc re t iza t ions m a k e i t d i f f i c u l t to e s t a b l i s h the s t r e n g t h of t h e i r m a t h e m a t i c a l f o u n d a t i o n . In a d d i t i o n , t h e i r successful a p p l i c a t i o n i n reac tor m o d e l l i n g relies o n h o m o g e n i z e d d i f fus ion a n d i n t e r a c t i o n coefficients t h a t are t y p i c a l l y a p p r o x i m a t e d t h r o u g h i d e a l i z e d l o c a l a u x i l i a r y c o m p u t a t i o n s . H e n c e the c o r r e s p o n d i n g d i s c r e t i z a t i o n error is d i f f i cu l t to a n a l y z e a n d the a p p l i c a t i o n of n o d a l m e t h o d s to genera l 1 Chapter 1. Introduction 2 diffusion problems seems nontrivial . Final ly, trie unusual choice of unknowns makes the resulting linear system difficult to solve as it appears incompatible wi th sophisticated iterative solution techniques. The objective of this research is to address the key issues, namely analysis and so-lut ion techniques, that have thus far restricted the application of these methods. In particular, we conduct a detailed study of the lowest-order nodal discretizations. This investigation suggests that the homogenization techniques employed i n reactor modelling are an integral part of the nodal methods dominance in this arena. The ult imate merit of nodal methods relative to other techniques is less clear in a general setting. However, for the related and important class of problems i n which piecewise constant coefficients are assumed, nodal discretizations are ideal. To enhance their competitiveness in these appli -cations we investigate a new general-purpose homogenization technique and we develop a new family of fast iterative solvers for the corresponding sparse linear systems. 1.1 Background Sketch The neutron distribution in a reactor core is one of the fundamental quantities required in the simulation and design of nuclear reactors. A n exact description of this distribution (i.e. the angular neutron density) is provided by the transport equation. Unfortunately, this integrodifferential equation is extremely difficult to solve and although its numerical treatment has improved significantly over the past 30 years, it is often replaced wi th approximate models [52]. The most popular of these is the PN approximation, which corresponds to an N T H order angular expansion of the angular neutron density, and is asymptotic to the transport equation [51]. This expansion leads to a system of ./V first order P D E s , and i n particular the Pi approximation leads to the first order form of the diffusion equation (i.e. a balance equation and Fick's law). In general, the energy Chapter 1. Introduction 3 dependence of the angular neutron density is discretized into groups such that a system of diffusion equations is obtained. These multi-group diffusion equations have been the backbone of reactor modelling for some 35 years and continue to play a significant role today. The diffusion model represents a significant simplification of the transport equation; however, the fine-scale variations in the coefficients in combination wi th the large phys-ical size of the reactor core pose a complex computational problem. In particular, a naive fine-mesh three-dimensional treatment of the multi-group diffusion equations was inconceivable i n the 60's [40] and remains intractable today [52]. Yet the coarse-scale properties of the solution (e.g. cell averages), that appear on a practical computational mesh, are significantly influenced by the fine-scale structure of the coefficients. Problems containing two or more significantly different length scales arise frequently i n physical applications and their treatment relies on various homogenization procedures. Specifi-cally, the fine-scale structure of the coefficients is averaged in some sense to define an homogenized diffusion problem in which the coefficients are piecewise smooth functions, typical ly constants, on the computational mesh. Unfortunately, the earliest treatments of the multigroup diffusion equations, which applied standard discretization techniques to this homogenized problem, proved inadequate. Efforts to bridge this gap between computational resources and sophisticated modelling have produced a powerful class of discretization techniques referred to as nodal methods. C o m m o n to all nodal methods is their rigorous enforcement of cell balance, a charac-teristic motivated by the underlying transport equation. In fact, the first nodal methods were based on a clever hybrid approach that combined transport physics and diffusion theory. However, the performance of these methods, now collectively known as nodal Chapter 1. Introduction 4 simulators, relied on the tuning of several reactor dependent free parameters [34]. More-over, these methods are not consistent with the multigroup diffusion equations, and as a result have largely been superseded by the modern, consistent family of nodal discretiza-tions [52]. These particular methods are based on a transverse integration procedure that reduces the multidimensional diffusion equation to a coupled set of O D E s involving transverse averaged quantities. This not only facilitates the use of unknowns compatible wi th the popular homogenization techniques (cell and edge moments i n two dimensions) but has contributed to their enhancement [65]. The treatment of these transverse inte-grated O D E s classifies the members of this family as either polynomial or analytic nodal methods. However, the influence of the enhanced homogenization techniques that were devel-oped specifically for reactor modelling is difficult to predict in a general setting. In ad-di t ion, the ut i l izat ion of transverse integration in conjunction wi th cell- and edge-based unknowns makes the discretizations themselves difficult to analyze. This uncertainty surrounding the mathematical foundation of nodal methods has significantly hindered their use in related applications outside of the reactor physics community. Recently, in the absence of homogenization, the equivalence of various nodal methods wi th cer-tain nonconforming and mixed-hybrid finite element methods ( F E M s ) was established [39, 37, 36]. Hence i n this restricted setting a solid theoretical foundation exists. The widespread popularity of nodal discretizations has also been hindered by the awkward linear systems which result. This is complicated by the fact that each equiva-lent discretization (mixed, mixed-hybrid and non conformal F E M s ) manifests this diffi-culty in a different way. For example, the mixed system is indefinite and although the el imination of certain unknowns leads to a positive definite system the relative sparsity is compromised. Various techniques have been employed to circumvent this l imitat ion Chapter 1. Introduction 5 {e.g. [61, 44, 2]}. C o n v e r s e l y , the n o n c o n f o r m a l F E M s are sparse b u t t y p i c a l l y c o n t a i n o n l y edge u n k n o w n s m a k i n g the d e f i n i t i o n of a su i tab le g r i d h i e r a r c h y a n d the i n t e r g r i d t ransfer operators d i f f i c u l t . Never the less , m u l t i l e v e l p r e c o n d i t i o n e r s have b e e n d e v e l o p e d for th i s s y s t e m a n d i n c l u d e the m u l t i l e v e l s u b s t r u c t u r i n g p r e c o n d i t i o n e r s d iscussed i n [29, 50]. H o w e v e r , the presence of s t rong ly a n i s o t r o p i c d i f fus ion m a y i n genera l degrade the p e r f o r m a n c e of these p r e c o n d i t i o n e r s . I n th is a n d m o s t o ther cases, solvers for the n o n c o n f o r m a l d i s c r e t i z a t i o n u t i l i z e the m i x e d - h y b r i d m e t h o d o n l y as a t o o l to es tab l i sh equiva lence w i t h the c o r r e s p o n d i n g m i x e d m e t h o d . W e c o n t e n d t h a t the m i x e d - h y b r i d m e t h o d is , i n fac t , c e n t r a l to the d e v e l o p m e n t of efficient solvers , as i t is c e n t r a l to the equiva lence of these d i s c r e t i z a t i o n s . 1.2 Outline T o e s t a b l i s h the se t t ing i n w h i c h n o d a l m e t h o d s o r i g i n a t e , a b r i e f i n t r o d u c t i o n to t r a n s -p o r t t h e o r y a n d the P i a p p r o x i m a t i o n is g i v e n i n C h a p t e r 2. T h e var ious t e r m s of the t r a n s p o r t e q u a t i o n are presented a long w i t h a d iscuss ion of the s i m p l i f y i n g a s s u m p t i o n s t h a t are t y p i c a l l y e m p l o y e d i n reactor m o d e l l i n g . T h e one-group d i f f u s i o n e q u a t i o n is de-r i v e d i n S e c t i o n 2.2.1 a n d the d i f fus ion a p p r o x i m a t i o n of the t y p i c a l t r a n s p o r t b o u n d a r y c o n d i t i o n s is addressed i n Sec t ion 2.2.3. A n h i s t o r i c a l o v e r v i e w of the m o s t p o p u l a r n o d a l m e t h o d s is presented i n C h a p t e r 3. W e res t r i c t our p r e s e n t a t i o n to the one-group d i f fus ion e q u a t i o n a n d b e g i n b y h i g h l i g h t i n g the k e y c o m p o n e n t s t h a t charac ter ize n o d a l m e t h o d s . P a r t i c u l a r l y , the c lass i ca l ce l l -based finite v o l u m e d i s c r e t i z a t i o n [Section 3.1] demonstra tes the r igorous t r e a t m e n t of ce l l b a l a n c e a n d serves to h i g h l i g h t the issues s u r r o u n d i n g i t . T h e m o t i v a t i o n of the n o d a l simulators is d iscussed i n Sec t ion 3.2 a long w i t h a d e t a i l e d r e v i e w of F L A R E [21]. I n a d d i t i o n to a r igorous t r e a t m e n t of ce l l ba lance , the m o d e r n , cons is tent , n o d a l Chapter 1. Introduction 6 methods rely on transverse integration. We review this technique and resulting transverse integrated O D E s in Section 3.3. The most popular variants of the polynomial and analytic modern, consistent, nodal methods are presented chronologically. The nodal expansion method ( N E M ) is discussed in Section 3.4 and the nodal integration method ( N I M ) appears i n Section 3.5. For the restricted case of conservative diffusion, a new equivalence is established between the lowest-order N E M and N I M discretizations. In Section 3.6 we begin to address the need for a strong mathematical foundation with the presentation of a new, rudimentary, analysis of these methods. The results of this analysis are verified i n the numerical examples of Section 3.7. Specifically, a contrived problem is used to highlight the second order convergence of the flux, while the superior accuracy of these methods, particularly N I M , is demonstrated in a more practical setting. The complex issues of homogenization are discussed i n Chapter 4. Specifically, i n Section 4.1 we review the equivalence-based techniques that are typical ly employed i n reactor modelling and are motivated by the desire to preserve certain coarse-scale proper-ties of the fine-scale solution. Unfortunately, closed-form expressions of the homogenized diffusion coefficient which satisfy exact equivalence relations are known in only a few cases, and thus, various approximations based on the numerical solution of local auxi l -iary problems have been developed. In Section 4.1.3 one of the most popular of these approximate techniques is discussed in detail and its crucial role is revealed. A more rig-orous mathematical approach to homogenization utilizes a two-scale asymptotic analysis and is well suited to periodic and nearly periodic problems. We review the important results of this analysis i n Section 4.2, and in particular, we identify the local auxiliary problems which define the homogenized diffusion tensor and deter its widespread appli -cation. In Section 4.3 we present a new multi-level technique [24] that does not rely on the solution of local auxiliary problems and demonstrates great promise for the periodic Chapter 1. Introduction 7 case. Moreover, this new approach is readily extended to nonperiodic problems. In Chapter 5 the equivalence of the lowest-order nodal discretizations wi th certain nonconforming and mixed-hybrid F E M s is studied. To provide an appropriate reference for comparison the standard bilinear conforming F E M is reviewed i n Section 5.1 and several relevant properties are noted. In Section 5.2 we explore the equivalence of these nodal methods wi th certain nonconforming F E M s presented i n [36, 37]. Specifically, we adopt a one-dimensional analysis that highlights the influence of the transverse integra-tion procedure and provides new insight into the characteristics and relationships of the nonconformal methods. Of particular importance is the more rigorous equivalence wi th N I M that this new analysis establishes. Final ly, the equivalence of mixed-hybrid F E M s w i t h the lowest-order N E M , that was established by Hennart and del Valle [39], is pre-sented in Section 5.3. In addition, we discuss the difficulty of extending this equivalence to N I M , but most importantly we demonstrate that the N I M system may be written i n the indefinite form associated with the mixed-hybrid F E M . This latter property is crucial to the success of our cell-based solver. In Chapter 6 we construct a new family of preconditioners that are based on the sparse approximation of various reduced systems. Specifically, we begin by el iminating the currents from the mixed-hybrid F E M system to obtain a reduced system for the scalar unknowns. Section 6.1 introduces the most common treatment of this reduced system, which eliminates the cell-based unknowns in favour of an exclusively edge-based system. This edge-based Schur complement is st i l l sparse and, moreover, is naturally partitioned by edge type (i.e. horizontal or vertical). However, a further reduction through the el imination of either edge type degrades the sparsity, and hence, a lumping procedure is introduced to approximate the Schur complement of this reduced system. The resulting approximation has a 9-point sparsity structure and is readily inverted wi th a standard Chapter 1. Introduction 8 mult igr id solver such as Black Box M u l t i g r i d [23]. We prove opt imal convergence rates for this edge-based preconditioner. Unfortunately, these rates depend on the aspect ra-tio of the grid, and hence, i n Section 6.2 we consider a two-step preconditioner which attempts to remove the directional bias of the lumping procedure. Ironically, this pre-conditioner is asymmetric, hence a more attractive possibility arises by returning to the reduced system for the scalar unknowns and proceeding to eliminate the edge variables [Section 6.3]. Again a loss of sparsity results and hence we apply a simple lumping to the original indefinite system that results in an approximate Schur complement which is the standard 5-point cell-centred discretization. Thus, efficient mult igr id inversion of this preconditioner is also possible. However, performing an exact mult igr id solve for each iteration of preconditioned conjugate gradient is undesirable, and hence in Section 6.4 we discuss the effectiveness of a single mult igr id cycle as an approximate inversion of the preconditioner. In Section 6.5 we present a numerical study of our solvers that includes some non-tr ivial problems. In all cases, our preconditioners perform extremely well : only a small number of iterations, independent of the grid size, is required. The robustness of the two-step edge-based and the cell-based preconditioners is evident i n the tests wi th high aspect ratio cells. F inal ly , conclusions and considerations for future research are offered i n Chapter 7. Chapter 2 Background Theory The foundation of reactor modelling is neutron transport theory. The corresponding equation, the neutron transport equation, is considered to be an essentially exact expres-sion governing the neutron distribution i n a reactor core. However, the generality and inherent complexity of this integrodifferential equation makes it exceedingly difficult to solve. Typical ly , simplifying approximations are made to reduce its complexity, and in cases such as diffusion theory, lead to a system of P D E s . Hence we present a concise introduction to transport theory which serves to highlight the natural hierarchy of the approximations that are inherent in any reactor model, and hence, any corresponding computation. A n important implicat ion of this hierarchy is that the terminology em-ployed i n vir tual ly al l approximate reactor models has its ancestry in transport theory. Thus, our overview of transport theory permits the introduction of this terminology i n its natural environment. 2.1 Neutron Transport Theory 2.1.1 The Neutron Transport Equation The fundamental quantity of interest in reactor modelling is the angular neutron density n(r, E, f i , i). It is actually a density in the statistical sense that n(r, D,Â£l,t)drdEdCl is the expected number of neutrons located in dr about r, within dE of the kinetic energy E, travelling within the solid angle dSl of f i , at a time in dt about t. Al though this 9 Chapter 2. Background Theory 10 definition may seem excessively detailed, its generality is fully uti l ized i n the derivation of its governing equation, the transport equation. The mechanisms by which neutrons interact with the various components of the re-actor naturally play a significant role in the development of this equation. Hence, these mechanisms are classified, with the two primary classes being scattering and absorp-tion. In addition, these interactions must be quantified, thus a macroscopic cross-section S a ( r , E) 1 is defined for each class that gives the probability per unit length that a neu-tron w i l l interact with a nucleus. Specifically, the probability per unit length of a neutron being absorbed at (r,E) is T,a(r,E). A n analogous definition holds for the macroscopic scattering cross-section S s (r , E). Interaction of any type combines these classes, and is described by the total macroscopic cross-section, E t ( r , Â£ ) = Â£ 8 ( r , Â£ ) + E 0 ( r , Â£ ) These general classes may be further subdivided. For example, absorption includes fission and radiative capture (among other mechanisms). Each of these mechanisms has a corresponding macroscopic cross-section, S/(r, E) and S 7 ( r , E) respectively. Al though a measure of the probability of neutron scatter has been defined, the angular dependence has not been considered. Moreover we have suppressed the dependence on the neutrons resultant energy E'. These dependencies are described by the double differential scattering cross-section, Â£ 5 ( r , Â£ E',n -> n') which gives the probability that a neutron init ial ly in the state (dE,dâ‚¬l) about (E,Q) is scattered into (dE',dQ ) about (E',Cl ). Integrating over energy and solid angle we x T o avoid confusion with the summation symbol, interaction coefficients wil l always appear with a subscript that denotes their classification. In definitions and statements which apply over all classes, the subscript a wi l l be used. Chapter 2. Background Theory 11 obtain, E s ( r , Â£ ) = E s ( r , Â£ , f i ) = J dft dE^Ss(r, E E', -)â€¢ )} clarifying the rather strange term: double differential cross-section. The dependence on incident direction is suppressed because the nuclei in any macroscopic sample are randomly oriented, thus averaging out any directional dependence. Having defined the neutron density and the means of interaction we can now de-fine the quantities that arise in the derivation and analysis of the transport equation. In particular, the speed of a neutron mult ipl ied by its interaction probability gives an interaction frequency (i.e. u E a ) . Hence, we define a reaction rate density by, f(r,E,to,t) = vXa(r,E)n(r,E,n,t) where the neutron speed is deduced from the kinetic energy, v = y2mE. The product vn(v, E,Cl,t) is often referred to as the angular neutron flux tp(r,E,H,t) = i m ( r , E , f i , i ) This unfortunate terminology is derived from the units, [speed] x [density] = [ cm] _ 2 [ s ] _ 1 , which are consistent wi th those of the term flux in many other disciplines (e.g. electro-magnetic theory). However, the term neutron flux is very misleading as it is a scalar quantity, not a vector, and hence characterizes the expected rate at which neutrons having a particular direction within dQ, about ft are passing through a given elemental surface at t ime i n dt about t regardless of the surface's orientation. 2 In contrast we define the angular current density, j ( r , Â£ , n , 0 = fi<p(r, 2 I n hopes of distinguishing it from any similar or related vector interpretation, and not merely to create an oxymoron, the flux is sometimes referred to as the scalar flux. Chapter 2. Background Theory 12 which assigns the direction to to the corresponding angular neutron flux and hence be-haves much more like a traditional flux. Specifically, if we consider an infinitesimal area element dA located at r then j(r, E,to,t)dAdEdto is the expected number of neutrons passing through the area dA per unit time with energy E in dE, direction to i n dto, at a t ime i n dt about t. Final ly , we address the question of neutron sources within the reactor. The primary mechanisms of neutron production are fission and various forms of decay. To describe the fission source we assume that v(E') is the average number of neutrons produced i n a fission event that was induced by a neutron of energy E', such that the total rate at which fission neutrons are created is, J dto dE'{v(E')Zf(E')(p(r,E,to,t)} The fission neutrons w i l l have an angular distribution that we assume to be isotropic, and a distribution i n energy which is described by the fission spectrum x(-E). The final consideration in defining the fission source is the time required for a fission event to release the fission neutrons. Amazingly, most (i.e. about 99%), of fission neutrons are emitted promptly (e.g. O(10~4) [s] in a thermal reactor [26]) while less than one percent are delayed, being emitted wi th a temporal distribution that depends on the exact fission chain and isotopes which produce them (e.g. C(10 _ 1 ) - O(102) [s]). If we approximate prompt as instantaneous then the prompt fission source may be written, sf(r, EAt) = *yp-J dto dE'{v(E')Xf(r, E')<p(r, E, to, t)} A detailed consideration of delayed neutrons, although fundamental i n reactor control, is beyond the scope of this thesis. Thus, we collect al l unspecified sources i n the generic term, q(r, E, to, t) and write the total source, s(r, E, to, t) = sf(r, E, to, t) + q(r, E, to, t) Chapter 2. Background Theory 13 The transport equation is commonly derived by considering the balance of neutrons wi th in an arbitrary elemental volume. Many authors have presented discussions of this derivation (e.g. [20, 16, 26, 9] ) which results in the following integrodifferential equation, ^ + n â€¢ Vif + Â£t(r, E)<f(v, E, ft, t) = f dQ I dE'\%s(r,E ^ E',n ^ Â£l)(p(r,E,Si,t)} + s(r,E,n,t) JA-K Jo ^ ' where the first term represents the rate of change with time; the second, the net change due to neutron flow; the third term, the flux lost to scatter and absorption; the fourth, the gain i n flux from scattered neutrons; and finally the source discussed above. 2.1.2 Simplifying Assumptions Introducing the transport equation (2.1) in such a general setting serves to convey the fundamental issues of neutronics modelling. However, some additional assumptions may now be introduced which simplify the transport equation without jeopardizing the insight it provides. Specifically we first introduce the one-speed approximation. Indeed this is a useful approximation in its own right, provided that the energy dependence of the interaction coefficients is removed carefully [26]. A less practical analogy i n neutron transport, although readily justified in other radiative settings such as monochromatic photon transport [57], is to assume that all scattering is elastic, Â£s(r, E fi ->â€¢ ft) = Â£s(r, E, ft n)5(E' - E) where S(E' â€” E) is the Dirac ^-function. This assumption decouples the energy depen-dence of the scattering source leaving all terms evaluated at the same energy, hence the explicit energy dependence may be suppressed. The second assumption is that the scattering is isotropic, Es(r, fY) = ^S s(r, E E') Chapter 2. Background Theory 14 Technically this is a relatively poor assumption. However, a more rigorous treatment is not justified for our purposes. F ina l ly we are interested in steady state problems so we w i l l assume that there is no t ime dependence. This requires that the fission source S/(r, E, fi, f), previously defined to be the source of prompt fission neutrons, be considered the source of all fission neutrons Sf(r,E,Q). Moreover it was assumed that the fission source is isotropic, thus we w i l l now assume the same of the generic source, q{r,E,n) = Â±Q{r,E) making the total source isotropic, s(r, E, Q) = - ^ 5 ( r , E) = i - ( S / ( r , E) + Q(r, E)) Combining al l of these assumptions results in the steady-state, one-speed transport equa-tion wi th isotropic scattering, ft- V V + St(rV(r,ft) = j <p(r,fi)dn +-^-S(r) (2.2) which, despite al l of these simplifications, is st i l l very difficult to solve. 2.1.3 Transport Boundary Conditions To complete the mathematical description of the steady state neutron transport problem (2.2), we must specify a domain V such that a boundary condition which prescribes the incident angular neutron flux may be imposed. The vacuum3 boundary condition is perhaps the most common in transport problems as it describes the interface between the medium of interest (i.e. the reactor core) and the surrounding vacuum. It s imply 3The vacuum boundary is also referred to as a free boundary and a purely absorbing boundary. Chapter 2. Background Theory 15 specifies that neutrons cannot return to the medium once they have left, and is expressed <p ( r , f i ) = 0 V (r G dV, ft n < 0) (2.3) where n is the outward unit normal to dV. Its complement is the reflecting boundary condition whose natural analogy as a perfect mirror suggests its most common application, the treatment of symmetry conditions. Mathematical ly the reflecting boundary condition is written V ? ( r , f i ) = </>(r,n r) V ( r e a y , f t - n < 0 ) (2.4) where ftr is the reflection defined such that ft â€¢ n = ftr â€¢ n and ft x ftr â€¢ n = 0. Fina l ly both of these conditions are included in the more general albedo condition, </?(r,ft) = av?(r , f t r ) V (r â‚¬ dV, ft â€¢ n < 0) (2.5) where the albedo, a = S s / S j Â£ [0,1]. In some applications it is necessary to generalize this condition to include an angular dependent albedo [57]. 2.2 D i f f u s i o n T h e o r y 2.2.1 T h e P i A p p r o x i m a t i o n M a n y of the difficulties encountered in attempting to solve the transport equation arise from the scattering integral. A technique commonly employed to approximate the influ-ence of this integral is the expansion of the angular flux i n a truncated series of spherical harmonics. Specifically we assume that ^ ( r , f t ) = J2 Â£ | ^ } V ( r ) y n m ( " ) (2-6) 71=0 m = â€” n where the expansion coefficients are denoted <pâ„¢(r) and y^ m (ft) are the spherical har-monics [16]. The uti l ization of this approximation in the transport equation yields a Chapter 2. Background Theory 16 system of N coupled first order partial differential equations that are referred to as the PN equations. We w i l l demonstrate that, under the simplifying assumptions of Section 2.1.2, the diffusion approximation arises naturally from the first order expansion (i.e. 7Y = 1) in spherical harmonics The validity of this truncated expansion has been studied through the asymptotic analysis of the transport equation. In fact, under certain scalings Larsen and Pomraning [51] have shown that, even with anisotropic scattering, PN theory is an asymptotic l imit of transport theory. Under the assumption of isotropic scatter the required scaling for Pi is such that the albedo a = 1 â€” O(e), thereby quantifying the more common stipulation that E t Â» Â£ a . Now expressing ft i n terms of spherical harmonics and uti l iz ing their orthogonality this may be written i n a more appealing form, <p(r,Q) 1 I 3 ^ (2.7) where the neutron flux <f>(r) is defined by J ( r ) Inserting this P i approximation (2.7) into the transport equation (2.2) yields, n â€¢ V[<f>(r) + 3ft â€¢ J ( r ) ] + E t ( r ) [cb(r) + 3ft â€¢ J ( r ) ] = E s ( r ) / [<f>{r) + 3n-J{r)]dn+XvZf[(t>(v) + 3Â£i- J ( r ) ] + Q ( r ) (2.8) Chapter 2. Background Theory 17 Next , integrating over ft one readily obtains V â€¢ J + E B ( r ) # r ) = i/xÂ£/(r)Â«Kr) + Q(r) (2.9) which is commonly referred to as the neutron continuity equation [27] or the neutron balance equation [41]. In fact, this balance equation is equivalently the zeroth angular moment of the transport equation, and hence, it is considered to be an exact expression. To develop a second equation we mult iply (2.8) by il and integrate over al l angles to obtain [16], i v # r ) + Â£ t ( r ) J ( r ) = 0 (2.10) which readily yields Fick 's Law, J ( r ) = - D ( r ) V < Â£ ( r ) (2.11) with the diffusion coefficient given by D(r) = [3Â£t(r)]_1. It is important to note that (2.11) does not correspond to the first angular moment of the transport equation. It is an approximate closure of the system made possible by the truncated angular expansion of ip(r, IT) [equation (2.7)]. In general an Nth order expansion w i l l lead to a system of N P D E s which is comprised of the first (./Vâ€” 1) angular moments of the transport equation and an approximate expression for the Nth angular moment that closes i t . We recognize that substitution of (2.11) into (2.9) gives the standard second order form of the diffusion equation - V â€¢ ( 0 ( r ) V # r ) ) + E 0 ( r ) # r ) = X ^ / ( r ) 0 ( r ) + Q(r) However, i n the following investigation of nodal discretizations it is important to preserve the distinction between the balance equation, which is an exact expression, and Fick 's Law, which is an approximation. In addition, not al l neutron transport problems arise Chapter 2. Background Theory 18 i n a m u l t i p l y i n g m e d i u m (i .e. S / ^ 0), a n d perhaps m o r e i m p o r t a n t l y , m a n y s igni f i cant p h e n o m e n a o u t s i d e of reactor phys ics are m o d e l l e d as di f fus ive processes. T h u s , we w i l l p r i m a r i l y inves t iga te the proper t ies of n o d a l m e t h o d s for general d i f f u s i o n p r o b l e m s . T o th is e n d i t is convenient to e x p l i c i t l y define two di f fus ion p r o b l e m s . T h e first, w h i c h we refer to as the one-group diffusion equation, is def ined by the f o l l o w i n g first-order s y s t e m , V-J(r) + Â£a(r)</>(r) = Q(r) (2.12a) J(r) = -D(r)V</)(r) (2.12b) a n d sub jec t to a general b o u n d a r y c o n d i t i o n of the f o r m , c(r)(f> + d(r)3 â€¢ n = g(r) (2.13) T h e l a t t e r p r o b l e m , w h i c h contains a fission source, w i l l be def ined i n the f o l l o w i n g sec t ion . 2.2.2 Criticality: The Reactor Eigenvalue I n a s teady-s ta te , one-group, source-free c a l c u l a t i o n we have the d i f f u s i o n e q u a t i o n , - V â€¢ (ZJ(r)V^(r)) + E â€ž ( r ) # r ) = A^ Â£,(r)</>(r) w h i c h , u n d e r homogeneous b o u n d a r y c o n d i t i o n s , m a y not have a s o l u t i o n . Y e t the exis -tence of a s o l u t i o n is prec ise ly the o b j e c t i v e of a m o d e l l i n g or c o n t r o l c o m p u t a t i o n as i t corresponds to a c o n f i g u r a t i o n w h i c h suppor ts a se l f -sustained c h a i n r e a c t i o n . T h u s , for p r a c t i c a l c o m p u t a t i o n s we i n t r o d u c e the multiplication constant fc, a n d w r i t e - V â€¢ (D(r)V</>(r)) + Ea(r)</>(r) = ^Xf(r)<f>(r) w h e r e k~l p e r f o r m s the role of an eigenvalue. In fact we are o n l y in teres ted i n the largest e igenvalue , as th is corresponds to the o n l y n o n n e g a t i v e e i g e n f u n c t i o n [26].This Chapter 2. Background Theory 19 transformation to an eigenvalue problem guarantees the existence of a solution, and moreover, allows one to characterize the criticality of the core: k > 1 supercritical k = 1 critical k < 1 subcritical suggesting i n rough terms the control which is needed to maintain a crit ical reactor core. Thus, we define the one-group, fission-source, diffusion problem, V â€¢ J(r) + E B(r)#r) = ^ ( r ^ r ) (2.14a) J = -Â£>(r)V</>(r) (2.14b) subject to the general (homogeneous) boundary condition (2.13), as the second problem. 2.2.3 A p p r o x i m a t e B o u n d a r y C o n d i t i o n s In practical computations the physical boundary conditions of a problem must be ap-proximated with in the context of the mathematical model. This is particularly evident wi th diffusion theory as the exact transport boundary conditions of Section 2.1.3 depend on ft, while the flux 0 (r ) and the net current J(r) are integrated quantities. Hence, it is not possible for diffusion theory to satisfy the transport boundary conditions exactly, and we must consider enforcing them in only the integral or weak sense. To this end we decompose the net current into partial currents, J - n = j â€ž + - j - (2-15) which are obtained by integrating over the respective half space ^ ' n ( r ) = / ( j ( r , f t )n }< i f t (2.16) The part ial currents give the expected number of neutrons passing through an elemental surface area, wi th directions characterized by either ft â€¢ n < 0 or ft â€¢ n > 0. To apply Chapter 2. Background Theory 20 these definitions within diffusion theory we must also employ the P i approximation of the angular flux (2.7), in W Â« ^ / { W') + 3 0 ' J ( r ) ] [ Â« â€¢ n ] } d O = ^ ( r ) Â± i j ( r ) â€¢ n (2.17) W i t h this expression i n hand, the approximation of the transport boundary conditions becomes transparent. In particular, the vacuum boundary condition is satisfied i n an integral sense by considering / l<p(rav,n)[n-n]\dÂ£l = 0 which is precisely the requirement that j~(r) = 0, and hence, from (2.17) we have the diffusion theory vacuum boundary condition, J n ( r ) = ^ ( r ) - ^ J ( r ) - n = 0 (2.18) Similar ly one can show that a reflecting boundary is defined by 3n(*)=ti(r) (2- 19a) and hence, corresponds to the following zero net current condition, J ( r ) - n = 0 (2.19b) In the general case of the albedo boundary condition, we obtain Jn(r) = aj+(r) (2.20a) which may be written, ^(1 - a)<f>(r) - + a ) J ( r ) â€¢ n = 0 (2.20b) We note that the l imit ing behaviour of (2.20b) i n the albedo is consistent wi th the transport conditions (i.e. a = 0 Jâ€ž ( r ) = 0 and a = 1 J ( r ) â€¢ n = 0). C h a p t e r 3 N o d a l D i s c r e t i z a t i o n s The combination of severe fine-scale variations i n the interaction coefficients and the large physical size of the reactor core made a naive, fine-mesh, three-dimensional com-putation inconceivable in the 60's [40] and continues to make it intractable today [52]. Nevertheless, the safe operation and continued development of nuclear reactors depend on sophisticated neutronics modelling. Research efforts to bridge the gap between com-putational resources and modelling demands has spanned some 35 years and is ongoing today. The most significant contribution of this research is a broad class of discretization techniques referred to as nodal methods. The distinct view of nodal ideology is that the diffusion equation, when expressed as a first order system of P D E s (2.12), is composed of an exact expression governing the neutron balance over an arbitrary volume, the bal-ance equation (2.12a), and an approximate relationship between the net current and the neutron flux, given by Fick's law (2.12b). This perspective differs significantly from that found in the numerical analysis community, which considers the entire system as an exact model , and it strongly influences the choices that are made throughout the development of the various nodal discretizations. To investigate these techniques we consider various nodal discretizations of the one-group diffusion equation (2.12) i n two dimensions (i.e. r = (x,y) and dr = dxdy) on a rectangular domain f i . The discretizations wi l l employ an L x M tensor product mesh having nodes fi.-j = [xi_^,xi+i] x [ J / J + I , Vj-^\- The cell centres are simply denoted (x;, yj), and it is convenient to define the mesh spacing Axi = xi+^ â€” Xi_^ and Ayj = yj+i â€”yj_^. 21 Chapter 3. Nodal Discretizations 22 ( X ; / 2 ' V M + ; / 2 ) (XL+l/2 > yM+l/2) 1 < t (Xl/2 A Xj (XL+l/2 Figure 3.1: A n L x M tensor product mesh imposed on the domain f2 Inherited naturally from the nodal perspective, and common to al l nodal discretiza-tions, is the choice of cell and edge-based unknowns. Generally these are taken to be moments up to some specified order, although we focus our discussion on the lowest-order case and hence employ simple averages. Specifically, the cell-based unknowns are averages defined by 4>i,j = A lA I + 2 / 3 + 2 (f>(x,y)dxdy (3.1) AxtAy3 Jx J ' - 1 3-? while the edge-based scalar unknowns, namely edge averages of the flux, are given by <f>i+Â±,j = (f)j(xi+i), 4>j{x) = -râ€”l <f>(x,y)dy (3.2a) faj+k = <f>i(yj+b)> &(y) = ^â€” 4>{x,y)dx (3.2b) LÂ±xi Jx. 1 Chapter 3. Nodal Discretizations 23 Similarly, the edge averaged currents are written 1 fyj+L l i m Jj(x), Jj(x) = â€” 2 [3(x,y) â€¢ nx]dy +*f. i ^Vj Jy. i 1-7 1 fX>+i r dx where = [1,0] T and n,, = [0,1] T . (3.3a) (3.3b) 3.1 C l a s s i c a l A p p r o a c h e s 3.1.1 F i n i t e V o l u m e a n d N o d a l B a l a n c e The nodal view of the first order system (2.12) suggests that a natural starting point is to integrate the exact balance equation (2.12a) over an arbitrary cell 0 4 ) J . Thus, i n the spirit of finite volume techniques we integrate equation (2.12a) over Clij and divide through by the cell volume Ar-j-j to obtain V â€¢ Jdr + 7-IL V â€¢ J* +2^7 IL SÂ°(RWR)RFR=^7IL Q{i)dr {3A) ArhJ J J A I ] A The first term of equation (3.4) is naturally handled by the divergence theorem, 1 ff 1 A r , V - J d r A r , J â€¢ nds where the line integral motivates the definition of the transverse averaged normal currents (3.3) and is thus discretized exactly by, 1 â€” / J â€¢ nds = â€”â€” J â€¢ 1 1 J - 1 + J â€¢ â€¢ 1 1 J â€¢ â€¢ 1 (3.5) We are not as fortunate with the second term, which is approximated by AT~ ff S a ( r ) ^ ( r ) d r = . ( E 0 ) (3.6) Chapter 3. Nodal Discretizations 24 where </>jj is the node average defined i n (3.1) and the cell-based absorption cross-section has been derived a priori wi th an homogenization procedure, = ftEa ( E 0 (a, y)) V r G fi,-j A discussion of the complex issues that surround this process w i l l be presented in Chap-ter 4. However, we note that if we continued to pursue an exact discretization of the balance equation, the homogenized absorption cross-section must i n fact be the mean value suggested by (3.6), fj EQ(r)</>(r)dr %)t<] = J-^j (3.7) / / <j>(r)dr which involves the unknown solution. Thus, the homogenized treatment of the second term is exact only if an iterative procedure is implemented, or i n the unlikely event that E a ( r ) is truly constant over For many, this approach may seem esoteric as it could be argued that (3.6) is a midpoint approximation of the integral, imply ing that ( E 0 ) i . and fyj are equivalently cell-centred point values. But the distinction is crit ical to nodal discretizations and typifies the desire to preserve the exact balance equation. F inal ly , we accommodate the last term by defining the node averaged source, such that combining this with (3.5) and (3.6) in equation (3.4) yields, Aâ€” \j- , â€¢-J+1]+-Jâ€” f j ~ , , - J+- J + ( S 0 ) . 4i j = Qi j (3.8) This is a discrete balance equation over (7,- ,^ and despite the aforementioned difficulty associated wi th defining the homogenized absorption cross-section, it is st i l l considered to be of greater physical significance than Fick's Law [cf. Section 2.2.1] Hence, the various Chapter 3. Nodal Discretizations 25 nodal methods have a common foundation in the balance equation, and are characterized by the techniques they employ to derive the additional equations necessary to close the discrete system. We note that this finite volume treatment of the one-group fission-source balance equation (2.14a) leads to the analogous discrete balance relation, 1 Ax J- - 7+ + l A y . J~J+i - + = f P ^ * " ( 3 ' 9 ) where the homogenized fission cross-section ( Â£ / ) is defined i n analogy wi th ( Â£ a ) (Xi-I/2 > 3^ -7/2) tyi+l/2Ji Ji+l/2,j Figure 3.2: The location of the unknowns 3.1.2 Standard Cell Centered Finite Differences One of the first techniques employed in reactor modelling was the uti l izat ion of finite dif-ference approximations of Fick 's law (2.12b) to eliminate the transverse averaged normal currents from the discrete balance equation (3.8) [40]. Specifically, to approximate the transverse averaged normal currents we substitute Fick 's law (2.12b) into their definition (3.3), to obtain, along wi th similar expressions for the other edges. We further assume that an homog-enized diffusion coefficient, whose evaluation w i l l be discussed i n Chapter 4, has been Chapter 3. Nodal Discretizations 26 denned on each cell [i.e. Dij = HD{D(X,I/)^] such that, J . ! . = l i m for which we use a one-sided difference approximation, 2Dit Similar arguments also yield, Thus, enforcing continuity of the normal current in the weak sense of the transverse average, we can eliminate i . j to obtain, Ax, '->i+l,j â€” (3.10) where A j ; + i = â€” Â£ ; ) and the interface diffusion coefficient D^+^j is given by the harmonic mean, 1 f A x ; A a : i + i 2DijD{+ijAx{+ DijAxi+i + 7J) i + 1 )jAa;j A x i + i ( 2 A j 2A+i , i Developing an analogous expression for Jij+i, followed by substitution into the discrete balance equation, we obtain a discrete system governing the node average, Ace Ace,--TJ Ax,-D; The convergence properties of this well known 5-point scheme have been investigated by a number of authors including Weiser and Wheeler [74] who have shown it is second-order under suitable assumptions. But the objective of nodal methods is to attain sufficient Chapter 3. Nodal Discretizations 27 accuracy on a very coarse mesh, approximately one cell per assembly. In this case the errors introduced by the approximations of J{+LJ and Jij+i can be severe [40]. Con-versely, the fine mesh required to reduce these errors, particularly i n three dimensions, is too costly in practice [52]. 3.2 S i m u l a t o r s The main objective of a neutronics simulator is the accurate and reliable prediction of i n -core neutron distributions that may, i n turn, be coupled to a global reactor model. Work on such simulators began in the early 60's and despite the severe memory restrictions of that era ( F L A R E [21] was developed for a system with 32Kb of R A M ) 3D simulations were st i l l targeted. The necessary compromise was the use of an extremely coarse mesh. However, a corresponding loss of accuracy could not be tolerated. This paradox was resolved wi th the development of a family of nodal methods, simulators, which are based on a clever hybrid of transport physics and diffusion theory. Simulators achieve impressive accuracy 1 once the associated free parameters are tuned to a particular reactor [34]. They have been used successfully i n a wide range of reactors including P W R s , L W R s , B W R s and even fast breeder reactors (see the review by G u p t a [34] and the references therein). However, from the perspective of truncation error analysis these methods are inconsistent wi th diffusion theory. Thus, it is primari ly for completeness that we provide a brief historical development of these methods. In simulator development, the nodal view of the first order system (2.14) is adopted. Since the discrete balance equation (3.9) is viewed as an exact relationship governing the average quantities it contains, it is the discretization of the net currents (3.10) which x In this context accuracy is relative to in core measurements and does not necessarily correspond to the associated discretization error. Chapter 3. Nodal Discretizations 28 is considered the primary source of error and hence must be improved. In fact, it was generally accepted that the underlying problem was the need to approximate the net current at an interface where it naturally embodied information from the neighbouring nodes. Thus, i n the following development due to Henry [40] and discussed by G u p t a [34] we consider rewriting the discrete balance equation (3.9) in terms of part ial currents (2.16). In particular we approximate the partial currents through the introduction of the part ial current coupling coefficients, = ^JyVj]HJ^X^y)}dy = 4 ! i ) ' 3 ) A x ^ (3- l la) Ji+^j = ~ty~. J {jnx(xi+b,y)}<ly = afill^Axi+^i+u (3.11b) wi th analogous definitions for the j/-oriented coupling coefficients aj^j 1 " 1 ' and a r j V ' ^ . The implicat ion of (3.11) is an approximation of the net current that may be writ ten, Ji+K = ~ Jty) = ~ L4+UA^+I<?W - a g J ^ A x , - ^ - ] (3.12) wi th the Ji,j+i defined similarly. Thus the partial current coupling coefficients represent a significant increase in the number of degrees of freedom. Final ly , ut i l iz ing (3.12) the balance equation (3.9) becomes, + A^+t A . A x , - ! uj) AyJ+1 { i J ) Aj / j - i (4+1^ + a?-}** + + 4*7*) + ( E â€ž ) . . - ^ ( E 7 ) . . V OJ) / V a/t,j fe V 3)i,3 which, provided we are able to evaluate the a's, is sti l l of the 5-point cell-centred variety. Unfortunately, this stencil w i l l be asymmetric i n general, making its solution significantly more difficult than its symmetric counterpart. Chapter 3. Nodal Discretizations 29 3.2.1 Fine Mesh Correction - Cross Coupling The Cross Coupling method of Wachspress [71] was perhaps the first neutronics sim-ulator and is based on the partial current coupling approach. It utilizes the diffusion theory expression for the partial currents [equation (2.17)] in conjunction wi th fine mesh calculations to compute the coupling coefficients. Specifically, a series of part ial core fine mesh calculations are performed to obtain a standard solution </>^m', d^l^p a n d Based on these fine mesh values the coupling coefficient is computed by first evaluating the part ial currents, e.g. where we acknowledge that the approximation of J â€¢ n employed here, although deemed inadequate on very coarse meshes, was considered reasonable for this fine mesh calcu-lation. W i t h the partial currents in hand, the coupling coefficients are readily obtained from equation (3.11). Unfortunately, the Cross Coupling method has the obvious weakness that if the prob-lem one desires to solve is not sufficiently close to the problem for which the fine mesh computations were performed, the results can be disastrous. It was for this reason, in conjunction wi th the additional cost of the fine mesh calculations, that other methods were sought. Nevertheless, attempts to utilize partial core or reduced dimension fine mesh computations to enhance the accuracy of the computed coarse mesh solution ap-pears i n many of the early methods, and w i l l appear again i n the more formal setting of Equivalence Theory [Section 4.1.3]. -4>W- -3.2.2 FLARE The first method that successfully approximated the coupling coefficients independently of diffusion theory and standard discretization techniques was Delp et al.'s F L A R E [21]. Chapter 3. Nodal Discretizations 30 Their approach first acknowledges that, the primary quantity of interest is the power produced i n each node, which is denoted Sij and defined by Sij = AxiAyj(Y,f).^4>itj Hence they introduce the modified coupling coefficients, _ _ * | _ a ( Â« ' + i . i ) _ k'ii+h,J J ^ ) - k Q.('-.i) _ k'ji+y where the substitution of from (3.11) has revealed that, w[f jj 1S the ratio of, the total number of neutrons leaving fi;j and entering 0 ;^ through their common edge, to the number of neutrons born i n In the F L A R E model the assumption is made that al l neutrons born i n O t J are absorbed either there, or i n one of its 4 nearest neighbours. This fundamental assumption implies that w^'1^ is equivalently the probability that a neutron born i n Qij is absorbed in its neighbour Clkti and facilitates the use of transport probabilities to estimate the coupling. Thus, the form of the balance equation employed i n most simulators is, â€” w) â€¢' n-OO ( h j ) ( i j ) ( h j ) ( h j ) A where the infinite medium neutron multiplication constant is given by, and = [ l - ^(-j)1"7' - w(ij)'') ~ - ^(i . i ) - ^] ' w n i c l 1 i s obtained directly from the balance equation, is also consistent with the definition of w^jy Chapter 3. Nodal Discretizations 31 This form of the balance equation is valid for all interior nodes wi th fuel bearing neighbours. Boundaries are treated with an albedo condition. For example at re = 0 we have, w; ' I " ( L i ) with w{()'J\ = 1 - wPA - (1 - a i 7 W(?'JJ - 1} ~ â„¢f!'-r0 for j = 2 , . . . , M - 1 where a i j is the albedo of the fictitious boundary cell and is considered a free parameter. A peculiarity of choosing S ; j as the unknowns is that non-fuel bearing cells must also be treated wi th albedo type boundary conditions because in the l imi t of i / ( S / ) . â€”> 0, S j â€” > 0 but wlk'1) â€”> oo. The absorption probability, w^'^ implemented in F L A R E is given by, (k,i) ,Mii3- M ? -= & ~ 9)~2A + where g is a free mixing parameter, M ? â€¢ is the migration area 2 and A is the mesh spacing (assumed constant) in the appropriate coordinate direction. Both terms arise from one-dimensional calculations, wi th the first based on transport theory and the second arising from a diffusion theory analysis. Although these approximations may seem crude, they are deceptively clever, chosen carefully to perform best on extremely coarse meshes. Nat-urally as simulators have evolved these kernels have improved, mult iple energy groups have been considered, and the strict nearest neighbour interaction assumption has been relaxed [34]. B u t despite these improvements and their many successes, simulators have not escaped the crit icism of being inconsistent methods which are laden with free param-eters that must be tuned. 2 T h e diffusion length Lij in the one-group diffusion model is defined by Lfj = (T,a). j / D i j , and is a measure of how far neutrons will diffuse before being absorbed. The migration area Mfj is an improved approximation of the true diffusion length (squared), that incorporates multigroup information. Chapter 3. Nodal Discretizations 32 3.3 Transverse Integration Mainta ining Henry's view [40] (i.e. Section 3.2) that the net current approximation is the primary source of error i n the discretization process, the use of part ial currents is considered a physically appealing proposition. However, the partial current coupling coefficients that are employed in nodal simulators lack any theoretical foundation. Thus researchers sought consistent discretization procedures which util ize the part ial currents directly. F innemann et al . [31] developed one of the first and most popular discretizations of this type, by introducing what has become a cornerstone of modern nodal methods, the transverse integration procedure. The transverse integration procedure first assumes that a homogenized diffusion co-efficient and absorption cross-section are defined on each cell. Equivalently, one may consider this an assumption that the problem itself has piecewise constant coefficients. In two dimensions transverse integration implies operating on an equation w i t h either of the following, - | ; { J ' ( Y ) } + ( S " ) . - > ' M = ~A~X~- K + i ( y ) - J?-M + Q t { y ) ( 3 1 3 b ) of the cell , while the transverse averaged sources Qj(x) and Qi(y) are defined i n analogy wi th the transverse averaged flux. The term pseudo-source is adopted because it w i l l prove useful to treat these transverse edge currents as sources despite the obvious fact Thus, transverse integrating the balance equation (2.12a) we obtain, Chapter 3. Nodal Discretizations 33 that they are fundamental unknowns. To emphasize this role the transverse leakage is commonly defined [52], %AX) = - (3-14a) cUy) =-^.{^(y) - ( 3 - 1 4 b ) Similar ly the transverse average of Fick's Law, which relies on the assumption of an homogenized diffusion coefficient, is simply, Jj(x) = - A j â€”{<7V(z)} x e [xZ_L,xi+i_] (3.15a) My) =-Dij-fiy-faiy)} y G [j/j-i>yj+*] (3.i5b) Inserting (3.14) and (3.15) into (3.13) we obtain a system of second order O D E s - A j ^ f a i O * ) } + ( E ^ X * ) = Â£'AX) + Qj(x) (3.16a) -D^3â€”{Uy)} + (Sa)f.^(y) = 4 , (y ) + &(y) (3-16b) There are two distinct applications of the transverse integrated equations. The first is the weighted residuals procedure of the nodal expansion method ( N E M ) which we outline below [Section 3.4.1] and the second is the Nodal Integration Method ( N I M ) which we discuss i n Section 3.5.1. 3.4 Polynomial Nodal Methods - NEM The development of modern consistent nodal discretizations began in the m i d 70's. These methods were based on local polynomial expansions. The first popular poly-nomial method was the nodal expansion method ( N E M ) of Finnemann et a l . [31]. In fact, although some variations and improvements have been considered [52], the N E M Chapter 3. Nodal Discretizations 34 ideology st i l l dominates the polynomial class of nodal methods. In its lowest-order form, N E M considers a quadratic expansion of the transverse averaged flux [i.e. <f>j{x) a n d 4>i(y)] on each cell. The expansion coefficients are determined by applying Fick 's Law i n combination wi th the discrete balance equation (3.8) and continuity of the normal cur-rent. Considerable effort has been made to utilize higher order polynomial expansions wi th in N E M . The difficulty this creates is centred around the evaluation of the higher order expansion coefficients. In particular, the weighted residual procedure that is typi -cally used relies on transverse integration of the balance equation (2.12a) and as a result an approximation of the transverse normal currents (i.e. the transverse leakage) is also required. In the following discussion we highlight the issues surrounding the weighted residual procedure and the approximation of the transverse leakage, although we focus on the lowest-order case. 3.4.1 C o m p o n e n t s of the N E M D i s c r e t i z a t i o n A P o l y n o m i a l Basis The N E M treatment of the transverse integrated O D E s (3.16) is based on a low order polynomial expansion of the transverse integrated flux. Specifically, we write, (3.17a) Ny (3.17b) k=i where VN Â£ span j f / 'o , . . . , IPN} is a polynomial of degree N [i.e. either Nx or Ny]. The basis functions employed i n [31] are chosen to satisfy certain properties at the expense of universal orthogonality. Chapter 3. Nodal Discretizations 35 D e f i n i t i o n 3.4.1. The first five N E M basis functions are given by, MO = i MO = f MO = ^ 2 ~ \ M0 = c(c-\){c + \) where Â£ = (a; â€” X{)/Axi in the expansion of <f>j(x) while Â£ = (y â€” yj)/Ayj for <f>i(y). Constructing the NEM basis. To construct the basis we first specify that '(/'/(rc) V/ > 1 be orthogonal to tpo{x). Hence mult iplying (3.17a) by tl>o(x) and integrating gives, a<CjÂ°l *[Mx)]2dx= / 1 4>3(x)Mx)dx J X 1 J X- l *-J ' - 2 It is natural to choose ^o(x) = 1 obtaining aj = As a consequence <4jÂ°' = </>;,j and thus these expansions are, by definition, consistent with the node average, = T / 4>j{x)dx = â€¢â€”- / </>i(y)^ y A x ; A y j J Now requiring that the linear and quadratic basis functions be orthogonal to ij)o(x) and {ipo(x),i()i(x)} respectively, we obtain a MO = Â« iÂ£ ( 3- 1 8 MO = ^[e-^\ (3.18b) where the normalization constants a l 5 ct2 are defined as follows. Recall that the definition of the transverse averaged flux implies that <?5>;Â±ij = 4lj{xiÂ±^)^ providing two constraints wi th which to define the normalization constants. Taking Nx = 2, substitution of x ; Â± i . Chapter 3. Nodal Discretizations 36 i n (3.18) results in the following linear system â€” \ot\ fiCt2 from which the coefficients are easily obtained, We take ai = 1 and a2 = 3 to obtain the basis functions T/>I(Â£) and 4>2(C)-This unusual construction of the quadratic expansion ensures the global continuity of the transverse averaged flux. However, we must prevent higher order terms from altering the edge averaged flux and thus higher order basis functions must satisfy, V > ; ( * , Â± i ) = o r = 3 , . . . , i v For example, for the cubic basis function, this constraint implies the form, so orthogonality with ipo(x) implies that CQ = 0. The normalization remains, and appears to be chosen arbitrarily, as setting c\ = 1 we obtain ^ 3 ( Â£ ) . S imilarly for the quartic basis functions we have, MO = [co + ClÂ£ + c2e] (t -\)(t+ \) and orthogonality gives c\ = 0 and c2 â€” â€”20c0. Again the normalization seems arbitrary, and taking c2 = 1 immediately yields V^ CO- '-' The consequence of this approach is that expansions of degree less than two are not considered while the expansion of degree two is completely determined by the edge and cell based fluxes. These quadratic expansions are given i n the following proposition. 0M) (*.2) hi (3.19) Chapter 3. Nodal Discretizations 37 Proposition 3.1. The coefficients of the lowest-order NEM expansion of the transverse averaged flux [i.e. NX = NY = 2 in equation (3.17)7 are given by (x'Â°) x Transverse Leakage Approximations In addition to the transverse averaged flux, the transverse integrated O D E s depend on the transverse leakages (3.14). Although these terms are treated as sources they are i n fact unknown, coupling the O D E s through Fick's law (3.15). Thus, in N E M the transverse leakages are also approximated with a polynomial expansion, 'h3 1=0 4S(y) = E6S,V(y) ^ [ M , M ] (3.20a) (3.20b) k=l where VM Â£ span{(/?o, â€¢ â€¢ â€¢, <PM} is a polynomial of degree M [i.e. either MX or MY]. In the following discussion the expansion of is considered in detail; analogous definitions and relations are assumed for d f ] . The quadratic expansion suggested in [31] utilizes the N E M basis [i.e. Definition 3.4.1], Hence, the expansion coefficients follow immediately from Definition 3.1, ,(37,0) _ r(x,0) ui,j ~ ' H i Â°i,3 ' J +L-i-hd Chapter 3. Nodal Discretizations 38 where we have adopted the notation (x,0) i r>+ C\'J(x)dx -7 ,(x). c7+u = l i n i C'JW l i m &Hx) C Unfortunately, the edge values of the transverse leakage are unknown and moreover are not considered with in the N E M discretization. Hence additional constraints are obtained by enforcing the continuity of the transverse leakage and its first derivative at cell edges. However, rather than incorporating these constraints into the discretization they are viewed in [31] as auxiliary equations, forming a tridiagonal system that relates the edge values j to the cell averages dfjÂ°\ This implies that (L + M) such systems must be inverted for each iteration of the global discretization. Unwi l l ing to accept the additional expense of these auxiliary calculations, most i m -plementations of N E M have employed the local technique suggested i n [10]. Specifically, a quadratic expansion in { l , Â£ , Â£ 2 } wi th Â£ = (x â€” X j ) / A x ; is adopted, and the expan-sion coefficients are determined such that the average transverse leakage over the central node fijj and its nearest neigbours ^;+i,j are reproduced. Thus the expansion coefficients of this local quadratic are characterized by the following system, (3.21) where (3l>k = AXl/Axk and 7 $ = + l ) , 7 $ = A , * ^ + 6A,fc + 3). A n alternative local scheme based on a piecewise linear expansion was also suggested in [31]. However, we are specifically interested i n the lowest-order N E M for which a A - 1 , Â» - f U a ) 1 (2) JL(*.O) r(x,0) 1 0 1 12 hj = r(x,0) A + M 2 7 j + M 1 (2) 12 7;+1,, r(x,0) Chapter 3. Nodal Discretizations 39 simple constant approximation, j Ax,-(3.22) is employed. We note that this zeroth moment approximation is consistent wi th the balance equation. Weighted Residual Procedure The development of the transverse integrated O D E s was motivated by the desire to util ize higher order polynomial expansions of the transverse averaged flux and transverse leakage. The transverse integrated O D E s facilitate this through a weighted residual procedure developed in [31]. Specifically, we consider a set of weight functions wn{^), n = 0 , . . . , Nw, such that operating on the transverse integrated O D E s (3.16) wi th / 1 wn(x)( â€¢ )dx , and / * wm(y)( â€¢ )dy Jx._i dy. y yields, f"* wm(y) \-D^2{Uy)}]dy + ( 2 a ) ,^g m ) = + f** rvm(y)^J(y)dy 3-2 J~2 where n = 0 , . . . , N^\ rn = 0 , . . . , N^ and cf)\1^ defines the weighted integral of the scalar flux, di*) _ h3 Jx. \ Jy-i * 7 J 7 x)wk(y)<f>(x,y)dxdy It is assumed that iwo(Â£) = 1, with x = Xi + Arc8-Â£ and y = y3 + AyjÂ£. We note that the most popular choice of weight functions is the N E M basis given in Definition 3.4.1. The resulting procedure is typically referred to as moment weighting [52]. Chapter 3. Nodal Discretizations 40 Final ly , to determine ivif' and Nw^ observe that with a quadratic expansion of the transverse average flux and a constant approximation of the transverse leakage no ad-ditional equations are required. Moreover, wi th itfo(Â£) = 1 both weighted transverse integrated O D E s yield the balance equation (3.8), implying that the quadratic/constant expansion is consistent wi th this zeroth order weighted residual procedure. Thus the required number of weighted residual equations in each spatial direction is Nffl = Nx + M x - 2 and Nly) = Ny + My - 2. 3.4.2 T h e L o w e s t - O r d e r N E M D e f i n i t i o n 3.4.2. The lowest-order NEM discretization of the one-group diffusion equa-t ion (2.12) is comprised of the cell balance equation (3.8) and the following cell-based normal currents, Jr+kJ = ^2~Kx-{^-^j ~ + 2&+h'i} (3.23a) J t h 3 = - 2 f ^ { - 2 & - * . ; + ~ &+bi) (3-23b) w i t h analogous expressions defining JfjÂ±L- The continuity of the normal current is i m -posed across al l interior edges, in the weak sense of the transverse average, J Z . i - = J t i ; J ~ ^ i = J + Â± i (3-24) i n conjunction wi th the discrete boundary equations to complete the discretization. Derivation of the Lowest-Order NEM. The lowest-order N E M assumes a constant leak-age approximation [i.e. Mx = My = 0] and hence the zeroth weighted residual of either transverse integrated O D E (3.16) is the balance equation (3.8). A quadratic flux expan-sion [i.e. Nx = Ny = 2] is also characteristic of,the lowest-order N E M . Substitution of this expansion (3.17a) into the transverse averaged Fick's Law (3.15a) immediately yields (3.23). The analogous treatment follows for the y oriented net currents. â€¢ Chapter 3. Nodal Discretizations 41 Interface Current Formulation Recal l the motivation was to develop a discretization that employed the part ial currents directly. This is now possible if we utilize the diffusion theory approximation of the partial currents (2.17) to write *Â»H = 20t+K + ^ j) (3.25a) Ji+y =Jt+^-37+y (3.25b) We note that enforcing the continuity of the normal current (3.24), i n conjunction wi th the continuity of the scalar flux ensures continuity of the partial currents. Definition 3.4.3. The interface current formulation of the lowest-order N E M [Defini-tion 3.4.2] is defined by the following expressions for the outgoing partial currents, where ^ = ( i + 8 f r f - ( 4 Â£ f with analogous definitions for the outgoing partial currents i n y (i.e. and and the coefficients r^'^ and T / J ' ^ . The partial current balance equation is writ ten i n the form, 1 2 ^ ( 1 + 4 ^ 1 + 1 2 M 1 + 4 ^ ) + (?.),. ,1 = <?,-Ax? I A x , J A y ] ( Ay,- I ' ^ + The discrete boundary conditions complete the discretization. Chapter 3. Nodal Discretizations 42 Derivation of the Interface Current Formulation. The objective is to eliminate the net currents and the edge averaged scalar flux from the N E M i n favour of the part ial currents. Substitution of (3.25) in (3.23) gives, { I + 8 Â£ ^ = { * ( 3 ' 2 8 a ) j1" 8S}^ = j1 + + 4 ^ &w+^w) + P - 2 8 " ) The interface current formulation now expresses the outgoing partial currents (i.e. j7_1 . and jf+L â€¢) i n terms of the incoming partial currents (i.e. â€¢ and ^ and the node average. Specifically, substituting the outgoing current J~_LJ defined by (3.28b) into (3.28a), collecting terms and simplifying yields (3.26a). Similarly we obtain (3.26b) and moreover, the rotated analogues hold in y. To express the balance equation in terms of incoming partial currents we employ (3.25b) to rewrite the average transverse leakage, 1 Ax,- Ax,- Oi+k,j Ji~h^ followed by substitution of the outgoing partial currents from (3.26), such that upon collecting terms and simplifying the coefficients we obtain, 1 A x , Thus equation (3.27) is obtained immediately with the substitution of this expression for the average transverse leakage in x, along with its y-oriented analogue into the balance equation (3.8). â€¢ A distinct advantage of the interface current formulation is that higher order approx-imations of the transverse leakage may be employed without significantly altering the sparsity structure of the system. This is not the case for the flux formulations which fol-low. Moreover, numerical studies have demonstrated that higher-order approximations Chapter 3. Nodal Discretizations 43 of the transverse leakage decrease the magnitude of the discretization error on a given mesh (e.g. [31, 52, 45]). Al though, to the best of our knowledge their influence on the order of the truncation error has not been studied. Flux Formulation Definition 3.4.4. The flux formulation of the lowest-order N E M [Definition 3.4.2] is defined by the following flux based balance equation, D 19/). . J.x:y) ~ 6 A x j ^ + 4 > ^ + AxtAy3 ' -6-^|(</>M_i + <f>iJ+h) =Q>,3 (3.29) w here (*,Â») _ AxjAyj ' i J ~ A x ? + A y $V) = 1 + AxiAyj Y1D. ^,3 and one-dimensional stencils that govern the transverse averaged flux 9 A j â€”1 Axi Ax; Ax i+1 Ji+h,3 | Di+l,3 A x i + 1 Di+1,3 + 2 D. A x i + 1 with the rotated analogue defining </>JJ+I [Figure 3.3]. (3.30) Derivation of the Flux Formulation. Ut i l iz ing the expressions for the one-sided currents given i n (3.23) we obtain, A x ; L along wi th an analogous expression for the y-oriented leakage. Substitution into the balance equation (3.8) immediately gives equation (3.29). The one-dimensional stencils that govern the transverse averaged flux [i.e. equation (3.30)] are obtained directly from imposing the continuity of normal currents (3.24). â€¢ Chapter 3. Nodal Discretizations 44 < 0 <|> -fi <8> -6 <2> 4 -0 ) Figure 3.3: Stencil weights of the flux formulation wi th D â€” 1 and S a = 0. Edge Based Flux Formulation Definition 3.4.5. The edge-based flux formulation of the lowest-order N E M [Definition 3.4.2] is defined by the following 7 point nearest neighbour flux based stencil, W ^ U . , - 4 - \f(x'y) 4 - f(x'y)] A,. .-. . J . , ^ 2 ^ ) ^ . , . â€¢ i _ 2 I , > â€¢ Â» ) â€ž ( * , ? / ) where 3(i-[/4rr)+i M A x , (3.31) (3.32a) (3.32b) and ti\xjy\ g\jV^ are given in Definition 3.4.4. The y-oriented analogue governs <f>ij+Â±-Chapter 3. Nodal Discretizations 45 Derivation of the Edge Based Flux Formulation. Solving (3.29) for the node average, h J 2 ,> .Â» ) 1 A Pi,3 ) + AxiAyj A y j ^ ' ^ t â€¢ ^ . j - i J â€¢ 6 D i j Hence substituting into (3.23a) and mult iplying through by Ai/j we obtain A y , Qi,j Ay,-Â«/., ! â€¢ 73 Ax,-2 - 3 A y , Q\T} Axi 4 - 3 A y , ^ h Ax t - ^(*.Â») A y 2 2 2 Pi,3 Qi,j Similar ly substitution into (3.23b) yields, A y , A y , Q%V) Axi Jx'y) jo . , A y ; A y 2 o^. , J , ) Thus we consider manipulating the coefficients to deduce, D. A y , l'3 A 2 - 3 A y , e%v) Axi .>â€¢Â») Pi,3 A , , A y , Ax,- 2 - 3 ^ r ] + 3 7 3 Pi,j (3.33) (3.34) (3.35) but 2 â€” 3[/u|j-'^] 1 = 3 ( l â€” [p\XjV^] *) â€” 1 and hence enforcing continuity of the normal current and ut i l iz ing the definitions of Cij'V^ a n d given i n (3.32), gives (3.60). â€¢ 3.5 Analytic Nodal Methods - N I M The second distinct class of nodal discretizations utilize analytic information that is ob-tained from the solution of the cell-based homogenized diffusion equation. Al though the analytic solution of the two-dimensional cell-based homogenized diffusion equation has been considered [58], typically analytic nodal methods are based on a particular treatment Chapter 3. Nodal Discretizations 46 Figure 3.4: The edge-based flux formulation wi th D = 1 and E a = 0. of the transverse integrated O D E s (3.16). Methods of this type include the Nodal Inte-gration Method ( N I M ) [32, 6] and the Nodal Green's Function Method ( N G F M ) [54]. In the lowest-order case these discretizations are i n fact equivalent [6], thus we focus on the more popular form, N I M . In contrast to the polynomial methods, the nodal integration method defines the transverse averaged flux [i.e. <f>j(x) and 4>i(y)} as the analytic solution of the corresponding transverse averaged O D E (3.16). This is facilitated by specifying cell-based boundary conditions i n conjunction with an approximation of the transverse edge currents [i.e. transverse leakage (3.14)] and the transverse averaged source. In the following we consider the assumption that these currents are approximately constant along their respective edges,and thus derive the constant-constant N I M discretization of the one-group diffusion equation (2.12). 3.5.1 The Nodal Integration Method Posing the cell-based transverse integrated O D E s (3.16) which govern the transverse average flux is the first step i n the N I M discretization. To this end we must specify the cell-based boundary conditions as well as the functional form of both the transverse edge currents and the transverse integrated source. Chapter 3. Nodal Discretizations 47 The definition of the edge average (3.2a) naturally yields Dirichlet boundary condi-tions for each cell [6], 4>j(xi-L) = fa-ij <f>j(xi+i) = (f>i+Ltj (3.36a) My,-*-) = <f>i,j-t = (3.36b) although other boundary conditions have been considered. In particular, Neumann boundary conditions were considered by Greenman et al . [33], but in the conservative case suffer the obvious problem of nonuniqueness. Fischer and Finnemann [32] employed part ial current boundary conditions as these naturally yield an interface current for-mulat ion of the discrete equations. In fact, this discretization is equivalent to the flux formulation that arises naturally with the Dirichlet boundary conditions (3.36). More-over, we are specifically interested i n the flux formulation, and thus restrict our discussion to the Dirichlet case. We assume that the edge-based normal currents are approximately constant, or equivalently we assume a constant transverse leakage. Various higher order approxi-mations are possible and were discussed i n Section 3.4.1 Final ly , the exact functional form of the source term Q(x, y) may be unknown, so we employ a piecewise constant mean value definition, Hence, the transverse integrated O D E s (3.16) become, ~Yx {Di*h*i{x)} + P-hhW = ~~Ky~j - J i - J + Q I * ( 3 - 3 7 A ) ~h { D i * h m } + ( S o ) Â« > ( y ) = ~ t {J+y - J - ^ l + Q i * ( 3 - 3 7 b ) Chapter 3. Nodal Discretizations 48 where x Â£ [ s 2 _ i , xi+i] and y Â£ [yj-^ respectively, subject to the Dirichlet B C ' s given i n (3.36). These O D E s are readily solved and thus, wi th expressions for 4>j{x) and 4>i(y) in hand, we outline the N I M discretization procedure. F irs t , we note that two independent definitions of the cell average are possible, (3.38) yielding 2 equations per node. Furthermore, ut i l iz ing the transverse averaged Fick's Law (3.15) we obtain expressions for the edge averaged normal currents, Jlu= l i r Â£ J^x) J 4 t * = l h ? Mv) (3-39) These 6 equations constitute only 5 linearly independent equations per node as the cell balance equation arises from two linear combinations. Imposing the continuity of the normal current i n the weak sense of the transverse average yields an equation for each interior edge [i.e. (L â€” 1)M + L(M â€” 1) equations] while the boundary conditions give rise to 2L + 2 M discrete boundary equations. Thus we have 7LM + L + M equations i n as many unknowns. 3.5.2 T h e constantâ€”constant N I M f o r C o n s e r v a t i v e D i f f u s i o n T h e o r e m 3 .1 . The constant-constant N I M discretization of the conservative one-group diffusion equation (2.12) on an L x M tensor product mesh is equivalent to the lowest-order NEM discretization given in Definition 3.4-2. Proof. We begin by solving the transverse integrated O D E s in a conservative medium. The solution of (3.37a) with (S a )^ . = 0 may be expressed as, (3.40) Chapter 3. Nodal Discretizations 49 hence substitution of (f>j(x) into Fick's Law (3.15a) gives ^ = - ^ t ^ - ^ Â± - " ^ } ( 3 - 4 1 a ) Similarly, the analogous expression for 4>i(y) gives, J U = -fS(Aj+* - tiFfa-k & w ~ J t^} ( 3 ' 4 1 b ) Evaluating the consistency expressions of the node average, and, (3.41c) (3.41d) These equations comprise the cell-based N I M discretization and we now proceed to prove their equivalence to the lowest-order N E M . Subtracting Jf_^ â€¢ from J~+l where both expressions are given i n (3.41a), and d i -viding through by Arc; we obtain, (3.42) which is the balance equation (3.8), and the first equation i n Definition 3.4.2. U t i l i z i n g the balance equation (3.42) in the expression for the node average (3.41c) we obtain, Jr+hJ - JU, = -6f^ (^ -+^ - + h-u) (3-43) Chapter 3. Nodal Discretizations 50 F ina l ly adding , . and J i + i j [equation (3.41a)] yields, (3.44) Hence adding (3.43) to (3.44) gives (3.23a) of Definition 3.4.2, while subtracting (3.43) from (3.44) we obtain (3.23b). Analogous expressions follow in y. Thus we have established the equivalence of the interior equations, and the conti-nuity of the normal current is enforced precisely as i n N E M . The equivalence of the discrete boundary conditions, follows as a consequence of the cell-based treatment that is employed i n both methods. â€¢ One implicat ion of Theorem 3.1 is that in the conservative case, the flux formulations of the constant-constant N I M are obtained as the conservative l imi t of those given i n Definition 3.4.4 and Definition 3.4.2. Specifically, with (^a) â€”> 0 we have p\Xj'y^ = 1. 3.5.3 The constantâ€”constant N I M for Nonconservative Diffusion Definition 3.5.1. The constant-constant NIM discretization of the one-group diffusion equation (2.12) is comprised of the cell balance equation (3.8) and the following cell-based normal currents, ^ = - ^ { ^ ~ ~ + ^ + ^ * i + t o } ( 3 - 4 5 a ) J?-iÂ»- = ~ A ^ H ^ + + 2 q ! ^ " ( a * " ^ * i + t o } ( 3 - 4 5 b ) where we have defined 1 (3.46) wi th , â€” \. .AT-. \W - - \ ^ ' (3.47) Chapter 3. Nodal Discretizations 51 and where analogous expressions define J^., , . The continuity of the normal current is imposed across al l interior edges, in the weak sense of the transverse average, T~ - T+ ' ' . â€¢ - u l (3.48) i n conjunction wi th the discrete boundary equations to complete the discretization. Deriving the constant-constant NIM Discretization. We begin by solving the transverse integrated O D E s . For example, consider equation (3.37a) whose solution may be written <f>j(x) s inh(Ai j [ r c i + i . - x]) sinh(Ag') s inh(Ai j [x â€” x^_i ) + Sinh(A>J) cosh [x â€” x,]) 1 _ cosh(Ag) { Qi,j (3.49) C (y.o) where Ajj and AJ^ are defined i n (3.47). A p p l y i n g the transverse integrated form of Fick 's Law (3.15a) we obtain, " ^ ^ ^ ( A i : ; ) ) ^ + i j " " + " Â£S0)} ( 3- 5Â° a ) sinh(Ag)) {<t>i+y - cosh(AgVJ_,J} - ^ ^ { Q ^ - Â£ & 0 ) } (3.50b) We continue by evaluating the consistency expressions of the node average, <f>i,j = i r-+* (pj(x)da = + + J - J - [ l - { Q * - A ^ } (3.51a) and, i - 4 ? " } (3.51b) Chapter 3. Nodal Discretizations 52 These equations comprise the cell-based N I M discretization and we now develop the expressions given in Definition 3.5.1. We begin by subtracting (3.50b) from (3.50a) to obtain, J~+h, ~ J t h , = -DijXij tanh(Aj J ) + , } + AxiSgfaj - C^f) (3.52) From (3.51a) we obtain an expression for (4>i+%,j + ^J '-^J) which upon substitution into (3.52) yields, after a l i t t le algebra, the balance equation, 4 } ' 0 ) + 4 f + ( S B ) , . X = (3.53) Now the balance equation (3.53) gives an expression for {Qi,j â€” Â£ ( j Â° ' } hence substitution into the node average (3.51a) gives, Jr+y - JU, = - ^ a $ { h - K - + <f>i+tj} (3-54) where a\Xj is defined in (3.46). A d d i n g (3.50b) to (3.50a) we find, wi th p\Xj defined i n (3.46). Now we simply add (3.54) to (3.55) to obtain (3.45a). Con-versely subtracting (3.54) from (3.55) yields (3.45b). â€¢ Chapter 3. Nodal Discretizations 53 The Flux Formulation Definition 3.5.2. The flux formulation of the constant-constant N I M discretization [Definition 3.5.1 ] is defined by the following flux based balance equation (3.56) where we have defined, J*,y) _ \j \i L1 ~ 1 fo c 7 \ and one-dimensional stencils of the form, OL. + A z - v " ~ t J ' ^ M y ' Ax- v "~ t + 1 ' J ' I (3.58) wi th a similar expression governing <f*i,j+%-Derivation of the Flux Formulation. The transverse leakage C\j'Â°^ is readily obtained from (3.54) and an analogous expression for may be similarly derived. Substi-tut ion of these expressions into the balance equation (3.8) gives (3.56), although the coefficient of <f>ij requires simplification. Specifically we have, Jx) aW ) f &w 6^ A , ? Ay] j ^ " ^ I [1 - [1 - 4''] â€¢ â€ž. = ^ â€¢ 4 [,_,<..] thus the definition of 7 ^ ' ^ . Imposing the continuity of the normal current w i t h (3.50a) and (3.50b) gives (3.58) â€¢ Chapter 3. Nodal Discretizations 54 Edge Based Flux Formulation Definition 3.5.3. The edge-based flux formulation of the constant-constant N I M given in Definition 3.5.1 is defined by the following 7-point nearest-neighbour flux based stencil, *,3 of] of] r , x ) 7i ,7 + v. . rv. . Ti',.?' 7i,j where A y ; - n . . ^ E i f [ i _ Â£ > ) ] i -Q{X)\ (3.61a) (3.61b) with and given in Definition 3.5.2. The y-oriented analogue governs (f>ij+L-Derivation of the Edge Based Flux Formulation. Solving the flux based balance equation (3.56) for the node average we obtain, 1 2D *,3 Hence substituting this expression into (3.45a) and mult iplying through by A y , gives, - 7 3 A x ; - O " ^ ) f A y ^ a j g ] 2 ] Ax,-1 v1'3 â€¢ ^ J Axt J**) ri+^ (3.62) (*) (2/) " i . i "1.1 A y ? a ! ? Chapter 3. Nodal Discretizations 55 Similarly from (3.45b) we obtain, (x) (y) A 2 (x) (z.y) [hd+t + hd-h) ~ (x-y)Â®^ Enforcing the continuity of the normal currents w i l l yield (3.60), however this is not apparent without further simplification of the coefficients. Specifically we note that, after a l i t t le algebraic manipulation, a\x}{ 1 -(oc) ~\ (x) (v) a\j Ay, I L M Ax; a) 'a) / hJ 1 ^ (x,y) Ax 1 - ^ and thus the definition of C | j ' ^ in (3.61a) and r?jj ^ in (3.61b) complete the derivation. â€¢ 3.6 Truncation Error Analysis The second-order convergence of the lowest-order N E M and N I M discretizations is rigor-ously proven through their equivalence to certain F E M s [see Chapter 5]. However, such equivalences have not been established for all nodal methods, and hence a more rudimen-tary treatment of consistency and stability would be required to prove convergence i n these cases. Unfortunately, a typical truncation analysis based on direct Taylor series ex-pansions is inappropriate for nodal discretizations. This is a consequence of the cell- and edge-based unknowns, which are not associated with a particular point i n their respective cell or along their respective edge. Hence a prudent treatment returns to the definitions of these unknowns, and considers the Taylor series expansion of the corresponding inte-grand. This technique is applicable to nodal methods in general, and is employed i n the following discussion to verify the second order consistency of the lowest-order N E M and Chapter 3. Nodal Discretizations 56 N I M , under the assumption of piecewise constant coefficients. It is then used to examine the consistency of these methods for piecewise smooth coefficients. We assume that the solution is sufficiently smooth on each cell fi,-j and write the Taylor series expansion about the cell centres (x ; ,y , ) in the form, J{x,y) â€¢ nx = p + px(x - x^ + py(y - y3) + â€¢â€¢â€¢ 3(x, y) â€¢ ny = q + qx(x - x{) + qy (y - y3) + â€¢ â€¢ â€¢ <j>(x,y) = <p + ipx(x - + <py(y -yj) + â€¢â€¢â€¢ where we have adopted a shorthand notation for the function and its derivatives evaluated at the centre of the cell (i.e. <f>(xi,yj) = ip). Thus, for example, the edge average < / > ; + i , j has the expansion, 1 Pi ^ JV-i-l ^ x^rX ^ Â»i+h>J = â€”l 1 n m 9{?,y))dy A y , = - T â€” / 2 \ <P + Â¥x{x ~ Xi) + ipy(y - yj) + â€¢ â€¢ â€¢ \dy Ax- A x 2 A y 2 = <P + -yV* + -cf-V** + -^^yy + 0{A3) (3.65) Although we assume that the source may be evaluated exactly, its analogous treatment is useful i n the following discussion. Once again we adopt a shorthand notation, Q(x,y) = 5 + sx(x - x^ + sy(y - y,) â€¢â€¢â€¢ such that the expansion of Qij may be written, A x 2 A y 2 Qi,i = s + ~ ^ s x x + - ^ s y y + 0 ( A 4 ) (3.66) We note again that the proposed truncation analysis relates only to consistency. To bridge the gap between consistency and convergence, the stability of the lowest-order Chapter 3. Nodal Discretizations 57 nodal discretizations must be established, and we accomplish this indirectly through their equivalence to certain F E M s [see Chapter 5]. Alternatively, a local mode analysis could be employed to prove the stability of these nodal discretizations directly. 3.6.1 Piecewise Constant Coefficients Using the approach outlined above we w i l l now investigate the truncation error of the expressions i n Definition 3.4.2 ( N E M ) and Definition 3.5.1 ( N I M ) . In particular, it was noted i n Section 3.1.1 that, under the assumption of piecewise constant coefficients and exact integration of the source, the balance equation (3.8) is an exact expression. Hence there is no truncation error, alance â€” 0 Next we consider the one-sided normal currents (3.23) of N E M , for which substitution of the necessary expansions gives A r c 2 T â€” 1 Jl 1 J T , â€” i o ^hjfxxx iÂ±^,i 12 with an analogous expression for the normal components in y. This is in contrast to the truncation error of the one-sided normal currents (3.45) of N I M , which are given by â€” Dij(pxxx + (Tia)..ipx Differentiating the balance equation with respect to x we see that, if the source was piecewise constant and <pyyx = 0 (i.e. ipyy =constant, a condition satisfied by an essentially one-dimensional problem) then this term would vanish, a property anticipated by the exact treatment of the transverse integrated O D E s . Moreover, i n situations where ipyyx and sx are small , the truncation error of N I M w i l l be significantly smaller than that of N E M . T 12 Chapter 3. Nodal Discretizations 58 Final ly , for both N E M and N I M continuity of normal current is imposed exactly i n the weak sense, but pointwise it is second order. Hence the lowest-order nodal methods are second order consistent on a non-uniform mesh. A similar approach may be used to treat the flux formulation and the edge-based flux formulation directly. 3.6.2 Piecewise Smooth Coefficients To evaluate the potential of nodal discretizations i n a more general setting and to sug-gest possible improvements we extend the previous analysis to include piecewise smooth coefficients. In particular, we assume that the coefficients may be expanded about the centre of the cell, T,a{x,y) = a + ax(x - xi) + ay(y - y3) + â€¢â€¢â€¢ D(x,y) = D + Dx(x - Xi) + Dy(y - y3) + This spatial dependence introduces an error in the balance equation, Tbalance = ^jyAx] G xip x + Ay](Ty(Py^ + 0(A3) B u t the truncation error of the one sided currents (both N E M and N I M ) becomes first order, TJ.+ 1 , = Â± ^ D ^ + 0 ( A 2 ) with analogous expressions for the normal currents in y. The continuity of normal cur-rents is unaffected. Thus the discretization is only first order consistent i n this more general setting. However, the first order term of the error is proportional to the deriva-tives of the coefficients and may be small i n practice. The restoration of second order consistency appears feasible but introduces addi-tional complexities into the discretization. Specifically, replacing the constant diffusion Chapter 3. Nodal Discretizations 59 coefficient with a diagonal tensor, D(x,y) D + Dx(x-Xi) 0 0 D + D y ( y - y j ) (3.67) would preserve the simple transverse average of Fick's law (3.15) while restoring second order truncation error in the one-sided currents. Note that the the use of N I M is st i l l possible, although it may not be possible to solve the transverse integrated O D E s with more general polynomial coefficients. Al though the truncation error in the balance equation is second order these gradients may be large, and hence on coarse meshes introduce unacceptable errors. E l iminat ing these second order terms within the nodal procedure is nontrivial . For example, assume that only linear terms are present in the expansion of the absorption coefficient, and consider the removal term of the balance equation, AxiAyj Jx Jy = Ax-Ay J J [o-+ crx(x - x^ + ay(y - yJ)\<f>(x,y)dxdy We have introduced the higher order cell-based moments of the flux (i.e. 4>tf\ 4>\Â°j^) t Â° accommodate the linear terms of E a (# , Â£/)â€¢ W i t h i n N E M a higher order expansion of the transverse averaged flux in conjunction with the weighted residual procedure [Section 3.4.1] may be used to derive the additional equations necessary to govern these new unknowns (e.g. [73]). Unfortunately, the application of N I M may not be feasible in this setting. Chapter 3. Nodal Discretizations 60 3.7 Numerical Experiments To confirm the theoretical prediction of the convergence rate a simple test problem wi th known solution was contrived. The results of these computations are presented i n Section 3.7.1. A more practical example from the reactor physics community, namely the iron-water benchmark problem [Dilber and Lewis [25]] is used to demonstrate the superior accuracy of nodal schemes, particularly N I M . 3.7.1 Example 1: A Simple Test To verify the results of the truncation error analysis [Section 3.6.1], consider the square domain [0,1] x [0,1], over which we take the diffusivity and the absorption to be constant. In particular, we choose D(x,y) = 1, E a (a ; ,y) = 1, the source term, Q(x,y) = 2 + 1 + (a2 + /3 2)TT 2 sin (2TTX) sin (2vry) and the Dirichlet boundary condition, (j>{x,y) = 2 V(x,y)e<9ft The exact solution is, 4>{x, y) = 2 + sin (airx) sin (2f3y) Thus, in the nodal discretizations we have Dij = ( S a ) . . = 1 wi th the cell-based source, ^ ~ T / o ?! f s in (7rAxi) ] f sin (7rAy,) ] . I n . . , n . Q , â€ž = 2 [1 + (aÂ»+ f y ] {-^f\ { ^ r | - câ„¢) - <*-Â») The results are displayed in Figure 3.5 where the transverse averages are displayed at their respective midpoints and the cell averages at the cell centres. B o t h surfaces Chapter 3. Nodal Discretizations 61 are i n good qualitative agreement with the exact solution. Al though this is encouraging we desire a quantitative estimate of the second order convergence of these nodal dis-cretizations. Unfortunately, the transverse averaged quantities associated wi th a given mesh, for example ( A x , A y ) do not exist on any finer mesh. This is also true of the cell averages. However, for certain fine mesh spacings, such as ( A x / r a , A y / n ) , where n is a positive integer, the coarse grid unknowns are, by definition, arithmetic means of the fine grid unknowns. Thus, the standard pointwise estimate of the convergence rate is complicated somewhat by the necessity to choose a coarse grid on which all errors w i l l be computed. In order to obtain convincing asymptotic convergence behaviour without considering excessively fine meshes, errors were computed on a 10 x 10 grid. The results for the transverse average quantities are displayed in Table 3.1, while the cell averages appear in Table 3.3. In both cases second order convergence of the constant-constant N I M discretization is observed. Similar results are readily obtained for the lowest-order N E M . It is also possible to obtain global convergence estimates by evaluating the error in an L p - n o r m . The results for the Loo-norm, which are shown in Table 3.2 and Table 3.3, also confirm the second order convergence of the constant-constant N I M discretization. Table 3.1: Convergence results for Example 1: the transverse averaged scalar flux on a 10 x 10 uniform mesh. Computat ional G r i d _i iAx.Av Ratio / iAx.Av Ratio 10 x 10 1.56802 x I O " 2 - 5.98930 x 10~ 3 -20 x 20 3.90134 x 10~ 3 4.019 1.49018 x I O " 3 4.019 40 x 40 9.74176 x 10~ 4 4.005 3.72102 x I O - 4 4.005 80 x 80 2.43472 x 10" 4 4.001 9.29979 x 10~ 5 4.001 160 x 160 6.08634 x I O " 5 4.000 2.32477 x I O - 5 4.000 Chapter 3. Nodal Discretizations 62 Table 3.2: Convergence results for Example 1 i n the L e o - n o r m . G r i d / lAx.Ay Ratio i iAx.Av ^H-hl+h oo Ratio 10 x 10 20 x 20 40 x 40 80 x 80 160 x 160 3.13604 x 10~ 2 8.20423 x I O " 3 2.07415 x I O " 3 5.19987 x I O " 4 1.30088 x I O " 4 3.822 3.955 3.989 3.997 3.13604 x 10~ 2 8.20423 x I O " 3 2.07415 x I O - 3 5.19987 x I O " 4 1.30088 x I O " 4 3.822 3.955 3.989 3.997 Table 3.3: Convergence results for Example 1: the cell averages. G r i d j. i&x,Ay 94,3 - <P4,3 Ratio oo Ratio 10 x 10 20 x 20 40 x 40 80 x 80 160 x 160 2.55358 x I O " 2 6.48514 x I O " 3 1.62752 x I O " 3 4.07269 x I O " 4 1.01841 x I O " 4 3.938 3.985 3.996 3.999 3.15640 x I O " 2 8.01607 x I O ' 3 2.06219 x I O " 3 5.19237 x I O " 4 1.30041 x I O " 4 3.938 3.887 3.972 3.993 0.0 (a) (b) Figure 3.5: The transverse averaged scalar fluxes, {<&+Â£j , <&j+Â£} obtained wi th the constant-constant N I M for this simple test case are shown in (a) while the corresponding node averages, </>,j are plotted in (b). Chapter 3. Nodal Discretizations 63 3.7.2 E x a m p l e 2: T h e I r o n - W a t e r B e n c h m a r k Since the inception of reactor modelling, organizations within the reactor physics com-munity such as the International Atomic Energy Agency ( I A E A ) [4] have been designing and cataloguing benchmark problems. M a n y of these benchmarks have been used by a number of authors (see review by Lawrence [52]) to demonstrate the superior accuracy of nodal methods. Unfortunately, the majority of these benchmarks are physically meaning-ful complex models of various reactor assemblies, and as such are beyond the scope of this thesis. A reasonable compromise is offered by the iron-water benchmark problem first employed by Dilber and Lewis [25]. Specifically, this benchmark requires relatively l i t t le background i n transport theory or reactor physics and yet provides compelling evidence of the nodal methods superiority, particularly highlighting the performance of N I M . The governing equation i n the iron-water benchmark is the one-group diffusion equa-tion (2.12). We recall that in this setting the diffusion coefficient D, and the macroscopic absorption cross-section Â£ a , are defined by, where Â£ t is the total macroscopic cross-section and Â£<, is the macroscopic scattering cross-section. These macroscopic properties of iron and water are given in Table 3.4 along wi th the corresponding values of the diffusion coefficient and macroscopic removal cross-section obtained with the above definitions. The domain of this shielding experiment, 0 = [0,30] x [0, 30], is divided into four regions which are shown schematically i n Figure 3.6. The source is piecewise constant and may be written, D = 1 E a â€” E j â€” S s 3S Chapter 3. Nodal Discretizations Table 3.4: The diffusivity D, and the absorption cross-section S a , for Example 2 [25] Material E s[cm] 1 D [cm] E a[cm] 1 Water 3.33 3.31002 0.10010 0.01998 Iron 1.33 1.10523 0.25063 0.22477 3 0 21 Clt- Water fi3- Iron Q 2 - Water ,Cl,-Water 0 12 15 21 3 0 x[cm] Figure 3.6: Schematic of the iron-water benchmark problem [25] Chapter 3. Nodal Discretizations 65 Final ly , to specify the boundary conditions it is convenient to define the following segments, = { ( 0 < z < l , y = 0 ) U ( s = 0 , 0 < y < l ) } r+,+ = { ( 0 < x < l , y = l ) U ( x = l , 0 < y < l ) } such that on the lower left segment we have a reflecting or zero net current boundary condition (2.19b), J n = 0 V ( x , y ) e r _ > _ whereas on the upper right segment a vacuum boundary is assumed which implies that the inward directed partial current is set to zero as in (2.18), _I^ + I j . n = 0 v ( * , y ) e r + , + The homogenized coefficients which appear i n the nodal discretizations are tr iv ia l ly de-fined from Table 3.4 as the grid is chosen such that both the diffusion coefficient D(x, y) and the macroscopic absorption cross-section E a ( x , y ) , which are piecewise constant over f l , are constant over each cell f2, j . For the purposes of comparison the three discretizations developed in this Chapter were used to compute the solution of the iron-water benchmark problem. The first, which w i l l be referred to as the finite volume discretization, is the classical 5-point cell-based scheme obtained in Section 3.1. The distinguishing feature of this discretization is the natural definition of the interface diffusion coefficient as the harmonic mean of its neigh-bours. The second discretization is the lowest-order N E M that was developed i n Section 3.4. The fundamental improvement of this N E M over the finite volume discretization is that the currents are continuous and second order convergent. The constant-constant N I M of Section 3.5 is the third method and incorporates analytic information from the Chapter 3. Nodal Discretizations 66 transverse integrated O D E s . In the nonconservative case this results i n a more costly stencil for which the truncation analysis [Section 3.6] suggests a significant improvement in accuracy may result. The N I M solution of the iron water benchmark problem on a 10 x 10 computational mesh is shown in Figure 3.7. The use of an efficient mult igrid solver [23] with the finite volume discretization en-abled the generation of the reference solution on a 640 x 640 uniform mesh. The error associated wi th the reference solution is bounded in the maximum-norm by 0.01. Hence the deviations given i n Table 3.5 provide compelling evidence of the superior accuracy which may be obtained with nodal schemes. Specifically, i n the Z/2-norm, the N E M er-ror is approximately 3 times less than that of the finite volume method, while the N I M error is nearly 7 times less than that of N E M . Similar reductions of the error in the Loo-norm are also shown. The exceptional performance of N I M is further highlighted in the cross-sections that are displayed in Figure 3.8, where the N I M cross-section is v ir tual ly indistinguishable from that of the reference solution. We also note that [Figure 3.8(b)] on this very coarse mesh the N E M discretization generates a solution wi th slightly negative, and therefore nonphysical, values near the layer. In contrast the N I M solution is positive everywhere. This observation is consistent with our forthcoming conjectures in Section 5.2.2 regarding N E M and in Section 5.2.3 regarding N I M . Underlying the above comparisons of accuracy, is the computational cost of each method. In particular, the sparsity structure of the lowest-order N E M and constant-constant N I M are identical, and hence, the computational cost in these cases is identical. However, the computational cost of the finite volume method, on a given mesh, is lower. Taking into account the flop counts of these methods, it is possible to show that for this example, the work required to obtain comparable accuracy using the finite volume method is approximately 5 times that of N I M . Moreover, this computational savings is Chapter 3. Nodal Discretizations 67 independent of the mesh-size that is used to obtain the N I M approximation. This is crit ical aspect of a fair comparison, and is facilitated by the new family of fast iterative solvers that are developed i n Chapter 6 and exhibit mesh-independent convergence. Table 3.5: Errors for Example 2 on a 10 x 10 computational grid. Discretization / IAX.ATJ fa-hj' 2 / lAx.Ay oo Fini te Volume lowest-order N E M constant-constant N I M 0.991 0.315 0.047 4.396 1.126 0.146 Chapter 3. Nodal Discretizations 68 0 1 0 2 0 3 0 1 5 2 0 2 5 3 0 X X (a) (b) Figure 3.8: Solution cross-sections for Example 2 are compared to a reference solution that was obtained using the finite volume discretization on a 640 x 640 uniform mesh: yj = 10.5 (a) full domain (b) zoom on layer. Chapter 4 Homogenizat ion In a practical reactor model, the coefficients of the multi-group diffusion equations exhibit strong variations on two significantly different length scales. Unfortunately, the compu-tational cost associated with the fine-scale mesh is prohibitive, and yet the fine-scale behaviour of the coefficients significantly influences the coarse-scale properties of the so-lut ion. Generally, an efficient computational treatment of such two-scale problems relies on averaging the differential operator (in some sense) over the microscopic or fine scale to obtain an homogenized problem governing the macroscopic or coarse-scale solution. It is specifically this class of problems for which nodal methods were developed [Chapter 3] and their excellent performance (cf. Section 3.7) is tightly coupled to the corresponding homogenized coefficients. We note that the strong multi-scale dependence of coefficients arises i n many applications outside of reactor modelling, including flow i n porous media, heat transfer, and electromagnetism. Hence, the process of homogenization has received considerable attention both by the reactor physics community (see review by Smith [65]) and by others (cf. [11, 62] ). The most popular techniques within the reactor physics community pursue a physi-cally motivated equivalence between certain coarse-scale properties (e.g. the coarse-scale cell average) of the solution of the two-scale problem, and the solution of the homogenized problem. A strict equivalence is not possible in general, and in Section 4.1 we review a popular approximation of Generalized Equivalence Theory. In other applications a rigor-ous mathematical technique based on a two-scale asymptotic analysis has been employed. 69 Chapter 4. Homogenization 70 This technique has been particularly effective for periodic and nearly (i.e. non-uniformly) periodic structures and we review the relevant results of this work i n Section 4.2. F inal ly , we present a new multilevel approximate homogenization algorithm in Section 4.3 that maintains many of the important properties of the rigorous asymptotic analysis and is considerably more efficient. 4.1 Equivalence Techniques in Reactor Modelling There are, in fact, two separate issues of homogenization i n a reactor core, and conse-quently, a two stage homogenization process. The first stage of homogenization utilizes various transport approximations and probabilistic methods to account for the break-down of the Pi approximation near interfaces, such as fuel pins, as well as the strong energy dependence of certain effects, such as absorption resonances [65]. These issues, while extremely important for the success of the model, are beyond the scope of the thesis. As our starting point we consider the output of this first stage homogenization procedure: the heterogeneous multi-group diffusion problem, which contains piecewise constant coefficients on the fine scale. In the interest of consistency and clarity we re-strict our discussion to the one-group form, 1 and explicit ly define the heterogeneous and homogenized one-group fission-source diffusion problems [cf. Definition 2.14]. Definition 4.1.1. The heterogeneous one-group diffusion problem is defined by, 1 T h e energy dependence as well as its potential homogenization, often referred to as condensation or collapse, play an important role in reactor modelling. This dependence is readily included in the analysis and appears in most of the methods discussed here. V - J + Â£ a 0 = ^ E / 0 J = -DV4> (4.1b) (4.1a) Chapter 4. Homogenization 71 where the parameters D(x,y), Â£ a ( x , y ) , E/(a; ,y) are prescribed functions of (x,y) that are piecewise constant on the fine-scale cells Fij and subject to particular boundary conditions. The corresponding heterogeneous solution is ((/>, J). Definition 4.1.2. The homogeneous one-group diffusion problem is defined by, V - J + S ^ = (4.2a) K J = -73v<Â£ (4.2b) where the parameters D(x, y) , E a ( x , y), E / ( x , y) are to be computed with a homogeniza-tion procedure (typically piecewise constant on the coarse-scale cells C ; j ) and ( 0 , J) is the corresponding solution. The objective set out by the coarse-mesh nodal methods is for these homogenized parameters to be piecewise constant functions (i.e. constant on an assembly size coarse-mesh cell Cij). In this situation we have shown that N E M and N I M lead to a discrete solution (bh that is a second order approximation of </>. In practice, errors of only 1-2% in average assembly powers are observed on meshes having only a single cell per assembly [65]. However, we demonstrated that for piecewise smooth coefficients these methods are only first order [Section 3.6] and in general, (f) â€” <f>= 0(AC) where A c is the assembly size coarse-mesh spacing. This leads to errors of approximately 5-10% between the average assembly powers of cbh and those of (b. Certainly a loss of fine-scale information is an inevitable consequence of the homogenization process. However, one would hope for better accuracy i n such coarse-scale properties as the assembly averaged powers. Thus, the essential objective of the homogenization procedure is to reduce the errors associated wi th various prescribed coarse-scale quantities. Chapter 4. Homogenization 72 4.1.1 Coarse-Scale Properties To motivate the homogenization techniques which follow we identify several integrated quantities that would ideally be preserved under homogenization. Studying the i m p l i -cation that this equivalence has in the definition of the homogenized parameters serves to demonstrate the extreme difficulties inherent i n this process. A natural quantity to consider is the assembly averaged reaction rate whose preservation would require (a = a or a = / ) which under the assumption of piecewise constant homogenized pa-rameters gives: Thus, we encounter the first well-known problem i n defining an homogenization proce-dure: the ideal homogenized interaction coefficients depend upon the unknown heteroge-neous solution as well as the homogeneous solution, the solution to the very problem we are working to define. Even if an efficient nonlinear iteration could be defined to solve for both homogenized quantities simultaneously, the fact remains that the heterogeneous solution is unknown. The edge-averaged net current on each edge of an assembly is another important quantity and its preservation provides insight into the additional problems associated wi th defining an ideal homogenized diffusion coefficient. Specifically we would require, I T,a(x,y)(/>(x,y)dxdy = / T,a(x,y)(j)(x,y)dxdy (4.3) (4.4) (4.5) where dC^ is the k edge of the coarse-mesh cell C,- j . Ut i l i z ing the homogenized Fick's Chapter 4. Homogenization 73 Law (4.2b), we could define the homogenized diffusion coefficient as, J (x , y) â€¢ nds lac Dij = â€”f (4.6) JdC% Therefore, overlooking the dependence on the unknown heterogeneous solution and the problem's own solution, each normal current equivalence (4.5) may result i n a different value of Dij. Hence, a constant diffusion coefficient on each coarse-mesh cell is insuf-ficient. This is not surprising: choosing J â€¢ n to be an edge-based quantity suggests, through Fick 's law, that D should also be edge-based, or at least be described by a cor-responding number of degrees of freedom. For example, one might consider the diagonal linear tensor (3.67) or alternatively a constant, ful l tensor. Based on such suggestions it is clear that procedures which may preserve or reduce the errors associated wi th some desir-able coarse-scale properties of the heterogeneous solution may necessitate modifications of the assumptions inherent in the derivation of standard nodal methods. Nevertheless, the most significant problem facing the definition of homogenized parameters is their dependence on the solutions of the heterogeneous and homogeneous problems. 4.1.2 A One Dimensional S t u d y To clarify the objectives of equivalence based homogenization and identify certain classes of problems for which equivalences may be determined independently of <f), we begin wi th a discussion of the one-dimensional case. Specifically, in perhaps the first equiva-lence based homogenization procedure Selengut [64] suggested that the replacement of a heterogeneous cell wi th its homogenized counterpart should not perturb the solution outside the cell i n any way. In one dimension this interpretation of the homogenization Chapter 4. Homogenization 74 of the cell C ; = , corresponds to the constraint that the boundary fluxes, 'i+k = 4>{xi+h) = 4>(xt+Â±) and the boundary currents, Ji-i = J(xi_i) = J(a?;_i), J,-+i = J(xi+i) = J(xi+ of the heterogeneous and homogenized problems be equal. As an example consider the homogenization of a layered structure composed of infinite slabs: the cell Cf, = [â€”6, b] for which the diffusion coefficient is given by D{x) = D0 \x\ < a Di a <\x\ <b To this end, we first note that assuming a planar source at large negative x wi th a constant diffusion coefficient the solution has the general form, cj>(x) = A + Bx where A and B are constants determined by the strength and exact location of the source. If we now consider inserting the cell Cf, and enforce continuity of the scalar flux, we obtain (x)={ A + B A + B (1 + e0v) {x + sgn(x)e 0a) (1 + e0v) A+.Bx \x\ < a a < \x\ < b \x\ > b where 'Di -D0 Chapter 4. Homogenization 75 and v = a/b is the volume fraction of material 0 i n the cell Cb- Thus the flux constraint becomes t r iv ia l , while the equivalence of net currents yields the well known result, D0D1 D v 1 â€” v - i {l-v)D0 + vD1 the harmonic average. A n interesting benefit of this approach is that, by considering an infinite cyl indrical geometry with azimuthal symmetry or an infinite spherical geometry, the corresponding two- or three-dimensional results are obtained. The general expression may be written, (1 + nv)D0 + n(l - v)Dl (1 - v)D0 + (n + v)Dx 1 " } where v = (a/b)n+1 and n + 1 is the number of dimensions [64]. A comparison of the 2D form of this result with commonly used averages, asymptotic analysis and the new mult igr id procedure is made i n Section 4.3. Difficulties in the extension of this analysis to relevant, truly multidimensional prob-lems was likely the reason that the impact of this research was primari ly its underlying equivalence ideology. In fact, it wasn't unt i l 1976 that Kollas and Henry [49] conducted a more detailed study of these ideas. To demonstrate their strategy we consider the constant coefficient source-free problem, ~(y3-^-<Â£(x)) = 0, a<x<b dx t ax > where, after a l i t t le algebra, one can deduce the relationship between the endpoint values, (4.8) wi th A = b â€” a. Now, consider a piecewise diffusion coefficient: 4>(a) J(a) 1 I 0 1 D(x) D0 a < x < XQ DI XQ < x < b Chapter 4. Homogenization 76 Using (4.8) we readily obtain, '1(b) J(b)_ where A i = XQ â€” a and A2 = b â€” x0. Thus, enforcing the equality of these boundary values yields, which once again is the well known discrete result that, in one dimension, the har-monic average is considered the exact homogenization. However, in their paper Henry and Kollas specifically considered the use of this equivalence analysis for the one-group, one-dimensional diffusion equation. For symmetric regions they were able to determine homogenized parameters, (D,T,a) that were independent of the boundary values. How-ever, more significantly they demonstrated that for a general region the homogenization process is overdetermined. We note that this distinction does not arise wi th the conser-vative, source-free, diffusion problem discussed above. 4.1.3 Multidimensional Equivalence Techniques In more than one dimension a weaker equivalence is sought, as was suggested by the discussion in Section 4.1.1. In particular, a number of popular techniques have been i m -plemented which endeavor to approximate the equivalences suggested by (4.3) and (4.5). The common thread among these methods is the use of locally defined auxil iary prob-lems whose solution may be computed very accurately on a fine mesh. The most common of these considers each assembly wi th homogeneous Neumann boundary conditions and thus for convenient reference we provide the following definition. \4>(a) J(a) "1 A L ' '1 Â£*1 0 1 0 1 j(b) 4i _L 42." â€¢1 Chapter 4. Homogenization 77 Definition 4.1.3. The heterogeneous assembly one-group diffusion problem is defined by, V-J + Zatp^xvXf? (4.10a) J = -DV4> (4.10b) for (x,y) â‚¬ Cij and is subject to the reflecting boundary condition [i.e. homogeneous Neumann B C s (2.19b)] J - n = 0 V(x,y)â‚¬dCi,j The parameters D(x,y), E a ( x , y ) , Y,j(x,y) are prescribed functions of (x,y) that are piecewise constant on the fine-scale cells F , j The assembly solution is denoted [J,(p). The local solution {J, ip) serves as an approximation of the heterogeneous solution i n the desired equivalence relation. In the following discussion we present the first and per-haps most common of these techniques, flux weighting, as well as a prototypical example of its popular successors, an approximation of generalized equivalence theory. Flux Weighting In the 1970's [41] and well into the 1980's [65] the most common approach to defin-ing homogenized parameters was flux weighting. Motivated by (4.4) the heterogeneous solution d>(x,y) on C , j is approximated by ip(x,y) [Definition 4.1.3]. In addit ion, it is assumed that the node averages for this auxiliary problem and the homogenized problem are equal; that is, / ${x,y)dxdy= I 4>(x,y)dxdy Chapter 4. Homogenization 78 Direct substitution into (4.4) yields the flux weighted cross-sections2 / Z,a(x,y)fi(x,y)dxdy ( Â£ C = ^ 7 ( 4 - n ) / if{x,y)dxdy Typica l ly the homogenized diffusion coefficient is taken from the transport definition [see Section 2.2.1], D = [ 3 E i ] _ 1 hence the flux weighting of Â£< gives, D\iw) = { . } (4.12) / ^p{x,y)da ixdy which interestingly is a weighted harmonic average. The natural cri t icism of this approach is regarding the error associated wi th the zero net current boundary conditions. Not only can it be large but, perhaps more importantly, it is difficult to predict. The unfortunate implicat ion is that none of the aforementioned integral quantities are preserved and the errors associated with them are unknown. More-over, i n the conservative case the solution (p(x,y) is unique only up to a constant, causing obvious problems with these definitions. G e n e r a l i z e d E q u i v a l e n c e T h e o r y Koebke [48] demonstrated that if the exact heterogeneous solution <p was known then relaxing the continuity assumption of the edge averaged scalar flux would facilitate the preservation of al l node-integrated reaction rates and edge averaged net currents. In addition, the coarse-scale edge average of the heterogeneous scalar flux, may be recovered from the corresponding jump conditions. Naturally, in practice <f> is not known and only an approximate equivalence is possible. Two very similar approximate techniques have 2 In his review Smith [65] refers to these as assembly homogenized cross sections (AXSs) Chapter 4. Homogenization 79 emerged based on this equivalence theory, and we present an overview of Smith et al.'s [66] approximation of generalized equivalence theory. We begin by writ ing the balance equation over a coarse-mesh cell C,-j for the hetero-geneous problem i n the form, - Ju,}+- â€¢%-*} (4.13) where we have employed the following definitions of the node integrated flux and inter-action coefficients, = A â€” A - / 2 / 2 4>{x,y)dxdy (4.14a) Ax,-Ay ; - J^. t ./â€ž. 1 (S^) = 1 S^(x,yMx,y)<fxdy (4.14b) (a = t or a = / ) . We note that this mean value definition of (T,a) also appears in (3.7), where it arose naturally during the derivation of the cell-balance equation. The edge averaged currents also follow the definitions of Chapter 3 and are written Jlh,3 = ^ j j + ^ ( x f Â± h , y ) - n x c 4y â€¢â€” 2 with analogous definitions of J^., , . Similarly, we can write the balance equation over the coarse-mesh cell for the homogenized problem + {ZtJij&jAxiAyj = â€” ( S / j . ^ . ^ j A x . - A y j where we note that we have util ized the assumption that the homogenized interaction coefficients E a are constant over each cell and al l other definitions are analogous to the heterogeneous case. Assuming that the equivalence of the edge averaged normal currents Chapter 4. Homogenization 80 [i.e. enforcing (4.5)] is satisfied and selecting ( E a ) t . . = ( E a ) ^ . , subtracting (4.15) from (4.13) reveals that the node-average has been preserved, fc= / (f>(x,y)dxdy = / 4>(x,y)dxdy = fa and hence we preserve the reaction rates (4.3). Thus, success of the equivalence process lies with the assumed equivalence of the edge integrated normal currents. Returning to the traditional nodal methodology of transverse integration we write the transverse integrated O D E s [cf. equation (3.16)] for the homogenized problem, ' l y { d ^ l { v ) } + { ^ Â» ~ \uX&h}*<W = %AV) (4.16b) where <f>j(x), 4>i(y) are the transverse averages defined in (3.2) and the pseudo source terms have been written in the form, It is at this point that we are forced to assume that the exact homogeneous solution (f>(x, y) is known such that these transverse leakages may be treated as known. Thus, al l that remains is the selection of the appropriate boundary conditions. Re-calling that the objective of the analysis is the preservation of the edge integrated normal currents we choose these as the boundary conditions, Chapter 4. Homogenization 81 This choice not only satisfies our objective but automatically guarantees the continuity of the edge integrated normal currents within the global homogenized problem. More-over, the homogenized diffusion coefficient is arbitrary, as we have satisfied al l desired relations between the heterogeneous and homogeneous problems. This indicates an un-fortunate inconsistency between the formulation of the continuous homogenized problem and the homogenized discretization. The manifestation of this inconsistency is the loss of continuity i n the edge averaged scalar flux of the discrete homogenized problem. Al lowing the discontinuity of the edge-averaged scalar flux, introduces the following j u m p conditions, Ax)- = ^ t - y Ax)+ = ^i+hJ Ji,j ~ T Ji,j â€” T with analogous conditions defining / / j ^ - Once again application of the equivalence technique requires an approximation of </>(rc,y), for which they chose the local assembly solution ip(x,y). This approximation implies that the equivalence theory cross-sections and diffusion coefficient are precisely the F W C s , The crucial difference is that the specific definition of D 4 J is irrelevant as the compu-tation of the discontinuity factors represent all necessary degrees of freedom. However, evaluation of the discontinuity factors requires an additional approximation. Specifically, we solve the local assembly problem [i.e. Definition 4.1.3] for the homogenized assembly and note that wi th constant cross sections and zero net current boundary conditions the solution (p(x,y) is spatially constant, Chapter 4. Homogenization 82 but we have (pij = ipij hence we define the assembly discontinuity factors ( A D F s ) as, Ax)Â± _ <Pi,j Ay)Â± _ &,j Ji,j ~=F Ji,3 [XT The use of A D F s i n G E T offers some improvement over the direct application of the F W C s i n that if the errors associated with the zero net current assumption are small reaction rates w i l l be very accurate. But clearly the same problems arise here as well : the errors in this assumption are difficult to predict and may lead to significant errors in the solution. 4.2 A s y m p t o t i c H o m o g e n i z a t i o n T e c h n i q u e s There are several asymptotic methods which consider the problem of homogenization i n periodic or nearly periodic composite structures (cf. [11, 62]). Perhaps the most popular of these is based on a two-scale treatment that is facilitated by the introduction of a fine-scale parameter, e. Specifically one considers the coarse-scale dependence r, and the fine-scale dependence p = r/e, to be independent variables. Then the second order form of the heterogeneous problem may be written, L Â£(&) = Q(r) cbt subject to suitable B C s where Le is uniformly elliptic in e. Note that to simplify our discussion we have made the common assumption that the source is independent of the fine scale. A n asymptotic expansion of <f)e is then introduced, cbe = 0 o ( r , p) + efar, p) + e2cf>2(r, p) + â€¢ â€¢ â€¢ such that one may investigate the behaviour of cj>t as e â€”> 0 + . For the fine-scale periodic and near periodic cases of linear second-order ell iptic operators, one finds that (in a weak Chapter 4. Homogenization 83 sense) <be converges to </>o as e â€”> 0 + , where (bo is the solution of the homogenized problem, L W = Q(r) 4>o subject to corresponding B C s and the second-order elliptic operator, L, is referred to as the homogenized operator. The coefficients of the homogenized problem vary on the coarse scale i n the nearly periodic case, and are in fact constant for fine-scale periodic problems. It is important to note that i n general, L is not related to the operator l_o that arises from evaluating the l imit e â€”> 0 + of the problem's coefficients. In the following discussion we review the two-scale asymptotic analysis of Bensoussan et al . [11] and, in particular, discuss the relevant results for the fine-scale periodic and nearly (i.e. non-uniformly) periodic cases. A n excellent introduction to these concepts may be found i n [42]. 4.2.1 M u l t i d i m e n s i o n a l P e r i o d i c A n a l y s i s We consider the fine-scale periodic family of linear differential operators, L Â£ = V r â€¢ (V Q V r ) + S a (4.17) where D(p), p = r/e is a diffusion tensor defined such that L e is uniformly ell iptic in e. In addition ut i l iz ing the fine-scale periodicity we introduce the following local fine-scale diffusion problem, - V p â€¢ V(p)Vp(j) = 0, i n F (4.18a) <f> periodic i n F (4.18b) The corresponding weak variational formulation (cf. Section 5.1.1) may be writ ten as: Chapter 4. Homogenization 84 find d> â‚¬ W(F) such that ap(cb, rb) = (D(p)Vpcb, V p t f ) = 0 W(F) (4.19) where the space of admissible functions is W(F) = {(f>\4> Â£ H1(F),<f> periodic} . The results of a two-scale asymptotic analysis [11] are summarized in the following theorem. Theorem 4.1. For the globally periodic problem with LÂ£ given in (4.17) the homogenized operator may be written, L = V r â€¢ (PV r) + f Â£ with the homogenized absorption cross-section given by, (4.20a) where F is a fine-scale cell and Ap its corresponding volume. In addition, the elements of the homogenized diffusion tensor are written, v{i'3) = Xpap^{p)" Pl' x'(p)" p i ) (4-20b) with xJ(p) ^ e solution of, ap(XJ(p)-Pj^)=0 y^eW(F) (4.21) and the bilinear form, ap(-,-) is given in (4.19). We note that (4.21) is a local ell iptic boundary value problem that must be solved to obtain the homogenized diffusion tensor. This situation is very similar to that found i n the equivalence analysis of the previous section. However, an apparent difference is that the homogenization of the absorption coefficient (4.20a) does not depend on the fine-scale solution (i.e. xJ m this case), and instead is simply the integral average. But Chapter 4. Homogenization 85 here a strict fine-scale periodicity has been assumed, yielding c/>n(r, p) = <^ >o(r)- Thus substitution of (f>t into (4.4) gives (Tia).. = Â£ a + C(e), indicating that these methods are consistent and that the strict equivalence approach impl ic i t ly captures higher order terms. Unfortunately, as we noted in Section 4.1.1 enforcing this strict equivalence is not possible i n practice. Moreover, the error associated wi th the asymptotic treatment is bounded, \\d>c - ^o||LÂ°Â°(n) < c e -4.2.2 Multidimensional Nearly Periodic Analysis In many applications a strict fine-scale periodicity does not exist, and hence we consider the extension to a fine-scale non-uniformly periodic family of linear differential operators, U = V r â€¢ (V (r, ^) V r ) + E a (r, (4.22) where D(r, p) is a diffusion tensor defined such that Le is uniformly elliptic i n e. Once again we summarize the relevant results of [11]. Theorem 4.2. For the globally periodic problem with Le given in (4.22) the homogenized operator may be written, L = V P ( Â£ ( r ) - V P ) +El(r) The homogenized absorption cross-section is given by, and the elements of the homogenized diffusion tensor are written, V^\v) = J _ a p ( X J ( r , p ) - P j , X \ r , p ) - Pi) with xJ(r?P) solution of, - V p â€¢ (2>(r, p)Vp) (^"(r, p) ~ Pi) = 0 (4.23) subject to periodic boundary conditions. Chapter 4. Homogenization 86 Recall that in the fine-scale periodic case the homogenized diffusion tensor is defined by only a single ell iptic boundary value problem, whereas in the nearly periodic case we have an elliptic boundary problem (4.23) for each point r, i n the global domain 0 . In general, a practical application of this result considers the solution of (4.23), and hence the computation of the homogenized diffusion tensor, at only a few points wi th in each coarse-mesh cell. A piecewise smooth function may then be defined to interpolate these values, once again suggesting the generalization of nodal discretizations to piecewise smooth coefficients [cf. Section 3.6]. 4.3 Multigrid Homogenization The mult igr id solution of second order ell iptic P D E s has been studied extensively since Brandt 's [13] seminal paper appeared i n 1977. Several robust mult igr id methods have been developed which, even in the presence of very rough coefficients, demonstrate opti-mal convergence rates in many practical situations. The efficiency of any mult igr id al-gorithm depends heavily on the coarse grid operators' approximation of the coarse-scale properties of the fine-grid discretization. This strongly suggests that a robust mult igr id algorithm for solving such heterogeneous problems [equation (4.1)] must employ coarse grid operators that are representative of the discretization of the corresponding homoge-nized problem [equation (4.2)]. We asked the following question i n [24]: what information, if any, regarding the homogenized diffusion coefficient (or possibly tensor) could be ex-tracted from these coarse grid operators? The results, summarized i n this section, are very promising. Recently Knapek [46, 47] addressed this question i n a complementary manner and we comment on his approach i n the following discussion. Chapter 4. Homogenization 87 4.3.1 O b t a i n i n g the D i f f u s i o n T e n s o r In the periodic case the asymptotic analysis of problems wi th rapidly varying coefficients demonstrated that the homogenization process leads to a ful l diffusion tensor. Intuitively, this is the case if we consider an essentially one-dimensional heterogeneous material (e.g. shale) which is not aligned with the coordinate axes. Clearly the pr imary axes of the homogenized diffusion problem should be aligned wi th those of the essentially one-dimensional structure. In general, this is only possible if the homogenization yields a fu l l diffusion tensor. Thus, our objective is to develop a local technique to compute a 2 x 2 constant tensor for each cell of a coarse grid, based solely on the coarse grid operator itself. To this end we adopt the compass based notation displayed i n Figure 4.1 (i.e. denotes the stencil weight in the direction negative signs are treated expl ic i t ly) , and generate the 9-point coarse grid operators wi th the operator-induced variational coarsening of Dendy's black box mult igr id [23]. Then observing that any diffusion stencil may be viewed as a superposition of currents, an analysis of the current passing through a cell centred reference frame yields the following results. r> N W, 1* CN, NE W C w co CSW Figure 4.1: A compass based definition of an arbitrary 9-point stencil Chapter 4. Homogenization 88 Theorem 4.3. Consider the conforming bilinear finite element stencil for conservative anisotropic diffusion subject to periodic boundary conditions and under the assumption of a piecewise constant diffusion tensor. A second-order approximation of the coarse-scale diffusion tensor T>ij on 0 2 j is given by, Ax (-QE , oNE , cNW \ ( ONE oNW \ Ay "+â€¢ >^l,k + Â°/+l,fc J Â°l+l,k) (oNE QNW \ Ay (-QN , QNE , cNW \ \Â°l,k ~ Â°l+l,k) A x " r ^^k -T~Â°l+l,kj (4.24) where (Ax , A y ) denotes the grid spacing, (l,k) = (i â€” |, j â€” |) and we have defined, $i,k = 2 si^k+i) -qN _ Ifc^v , ON \ For a strictly constant diffusion tensor (i.e. T>ij â€” V), equation (4.24) is an exact expression. Proof. A detailed proof may be found in [24]. â€¢ Corollary 4.1. Given a coarse grid operator generated by operator-induced coarsening of the conforming bilinear stencil described in Theorem 4-3, the diffusion tensor T>{j given in (4.24) is a second-order approximation of the multigrid homogenized diffusion tensor. Moreover, if the stencil weights are spatially constant then (4.24) gives the multigrid homogenized tensor exactly. Proof. For this purpose the two key properties of the 9-point stencil are that it is sym-metric and conservative. Operator-induced coarsening preserves these properties, and hence, the proof follows from that of Theorem 4.3. â€¢ Remark 4-1- We note that, in Theorem 4.3, assuming a conforming bilinear basis ensures â€”E â€”N that the only error is associated with Slk and Slk. However, in general a second-order error also arises with the approximation of first derivatives (i.e. for the current) at the Chapter 4. Homogenization 89 cell centre. These two errors manifest themselves very differently. Specifically, St k and S^k introduce a second order error in T>ij which vanishes when the stencil is spatially constant. In contrast the errors associated with the solution do not influence Z \ j but lead to an error in the resulting stencil. Numerical ly this error has been demonstrated to be second-order in [47], although a proof remains elusive. Based on Corollary 4.1 our mult igr id homogenization algorithm consists of first defin-ing the conforming bilinear discretization on a fine grid capable of resolving the fine-scale structure of the problem. Then, the operator-induced variational coarsening of black box mult igr id is applied unt i l a coarse-mesh spacing suitable for efficient computation is reached. From this coarse-mesh problem we compute the cell-based mult igr id homoge-nized diffusion tensors using the cell-based expression (4.24). A n alternative vertex-based approach is considered in [47] which inverts a 9 x 9 system that is defined over a group of 4 cells. These methods result i n equivalent homogenized diffusion tensors when the stencil is spatially constant, in which case it is natural to assume that 4 neighbouring cells w i l l have identical properties. However, this assumption may be too restrictive in a more general setting and hence we feel that the local technique is preferable. Remark 4-2. Note that we have not entirely escaped the dependence that the homoge-nization process has on the heterogeneous solution. Instead, the mult igr id homogeniza-tion approach shifts this information from the actual heterogeneous solution to the finite dimensional subspace of the discrete solution. Choosing a different fine grid discretiza-tion corresponds to choosing a different subspace and w i l l yield different results. This is evident if we consider a simple constant coefficient diffusion problem and employ the stan-dard 5-point vertex based discrete operator. The homogenized coefficients determined wi th mult igr id homogenization w i l l be i n error, approaching the correct diffusion coeffi-cient only i n the l imit of infinite coarsening. This suggests that for a given discretization Chapter 4. Homogenization 90 it may be necessary to utilize a compatible sequence of subspaces i n the definition of the coarse grid operators and hence the homogenization process. 4.3.2 A N u m e r i c a l H o m o g e n i z a t i o n Sui te Unfortunately, there are very few exact theoretical homogenization results, wi th the vast majority restricted to one-dimensional problems. However, based on the asymptotic analysis outlined i n Section 4.2, numerical computations may be performed for the case of periodic boundary conditions which converge to the homogenized diffusion tensor [12]. Thus, we have conducted a series of tests in the periodic case to gauge the potential of this new technique. To accomplish the homogenization of a representative cell we acknowledge that the smallest plausible computational grid is 3 x 3 and hence define the corresponding problem as a 3 x 3 t i l ing of representative cells. A fine grid is selected which adequately represents the fine-scale structure of the representative cell and the operator induced coarsening is applied. Final ly, following Corollary 4.1, the homogenized diffusion tensors are given by (4.24). Several rudimentary tests, which are not included here, were conducted to verify the proper functioning of the code [24]. These included, the demonstration that the algorithm reproduces a constant tensor and that the homogenization of essentially one-dimensional problems is exact. Unfortunately it was also revealed that this homogenization technique is not infallible: T>ij is simply the arithmetic average for a checkerboard problem. We believe that this pathology, which is l inked to the lumping inherent i n operator induced coarsening, is an isolated example that is overshadowed by the method's performance in other cases. In particular, the tests that we present in the following subsections, which characterize the dependence of the homogenized diffusion tensor for a representative cell f i r on the shape and diffusivity of an internal inhomogeneity f i i (i.e. fii Â£ Qr) [12], clearly highlight the strengths of this new homogenization technique. Chapter 4. Homogenization 91 S h a p e D e p e n d e n c e A n evaluation of the geometric dependence of the homogenized diffusion tensor is ac-complished wi th the set of representative cells shown i n Figure 4.2(a)-(c). The diffusion tensor on these representative cells is defined by, where O 0 = O r \ O i . In al l cases the area of Oi is 1/4, implying the standard averag-ing techniques (e.g. arithmetic, harmonic) would yield the same homogenized diffusion coefficient independently of the inhomogeneity's shape. Based on their symmetry we anticipate that the homogenized diffusion tensors w i l l in fact be scalar multiples of the identity and this was confirmed numerically. A comparison of the blackbox homoge-nized diffusion coefficient with Bourgat's asymptotic computation is provided i n Table 4.1, where percentage differences, relative to the square inhomogeneity, are also included. These results demonstrate that the relative sensitivity of blackbox homogenization to the shape of the heterogeneity is similar to the rigorous treatment of Bourgat. In a direct comparison we find that the blackbox results consistently overestimate the asymptotic value by approximately 2-3%. This is impressive when a commonly employed alternative such as the two-dimensional harmonic average underestimates the asymptotic value by approximately 17%. D e p e n d e n c e on the DifTusiv i ty The second case is simply a square inhomogeneity [Figure 4.3] defined by, to evaluate the dependence of the homogenized diffusion tensor on the parameter A. Symmetry once again guarantees that the homogenized diffusion tensor is s imply a scalar V ( z , y ) â‚¬ O 0 V ( x , y ) e O i V ( x , y ) G 0, V ( x , y ) â‚¬ 0 Chapter 4. Homogenization 92 fl, \ Â« 6 (a) (b) (c) Figure 4.2: Three inhomogeneities with the same area of 1/4, but different shape. Table 4.1: A comparison of the shape dependence of the diffusivity relative to the square (i.e. the percent relative difference - % R D ) is presented for the results of Bourgat [12] and the black box homogenization. Shape Bourgat % R D Black Box % R D Square 1.548 - 1.598 -Disk 1.516 -2 .06 1.563 - 2 . 1 9 Lozenge 1.573 +1.69 1.608 +0.63 mult iple of the identity. The mult igrid homogenization results are presented in Figure 4.4 along wi th the asymptotic results of Bourgat, the aforementioned 2D result of Selengut (4.7) and the commonly used means, ^ = / / V ) d x d y = l ( X + 8 ) h l rl [D(x,y)] 1 dxdy 9A (1 + 8A) We note the excellent agreement of the blackbox homogenized diffusion coefficient with the asymptotic results over 8 orders of magnitude in A. We also observe the catastrophic failure of the harmonic mean for A â€”> 0 + is contrasted by an overestimation of approx-imately 10% in the arithmetic mean. Whereas for A â€”Â» +oo the harmonic mean yields approximately a 10% underestimation, while the arithmetic mean grows linearly, display-ing an arbitrari ly large error. Having been derived in an infinite cyl indrical geometry the Chapter 4. Homogenization 93 result of Selengut (4.7) is very interesting to include because it demonstrates the effect of the geometry. Here we see that for A â€”>â€¢ 0 + the l imi t ing value is very insensitive to the geometry while for A â€”> + o o it is quite sensitive, underestimating the final result by approximately 4%. Nevertheless, it is the only closed form expression to mimic the behaviour of the numerical calculations over the entire range of A. 1/3 2/3 Figure 4.3: A square inhomogeneity with diffusivity, A and an area of 1/9 Ih. 2 . 0 1 .8 1 . 6 1 .4 1 . 2 0.0 o Arithmetic â€¢ Harmonic v Black Box * Bourgat o Selengut â€¢ o 10 10" 10 10 10" 1 0' Figure 4.4: A comparison of homogenized diffusivities. Chapter 4. Homogenization 94 A N o n s y m m e t r i c I n h o m o g e n e i t y h V ( x , y ) e O 0 10/ 2 V ( x , y ) â‚¬ f t i Consider the L shaped region shown in Figure 4.5 with the diffusion tensor, V(x,y) = from which the asymptotic computation yields [12], = Q v{as) = 1.915 -0.101 -0.101 1.915 2.016 0 0 1.814 Q1 where the matrix of eigenvectors Q is given by, and defines the principal axes of diffusion, in this case a rotation of 45 c homogenization also gives a full tensor Blackbox 1.959 -0.153 -0.153 1.959 Q 2.113 0 0 1.806 Q T Moreover, we remarkably obtain the exact principal axes of diffusion in this case. The only error is the scaling in each of these directions, approximately 5% and 0.4% respec-tively. 0 1/6 3/6 5/6 1 Figure 4.5: The homogenization of an L shaped inhomogeneity leads to a dense tensor. Chapter 5 Finite Element Methods The development of modern nodal methods is based almost entirely on physical argu-ments. However, traditionally no assurance was given of a correspondingly solid math-ematical foundation. Over the years there has been a continuing effort to remedy this situation wi th definite advances i n some areas. One such vein of research, pursued by Hennart [35, 36], has established the equivalence of certain nodal and nonconformal finite element methods ( F E M s ) . Equivalence results such as this form important contributions to the mathematical foundation of nodal methods, because the related F E M theoretical results are immediately applicable to the corresponding nodal methods. In particular, stability and optimality results for the corresponding nodal discretization may be proven in this way. We explore this equivalence for the lowest-order nodal methods that were introduced i n Chapter 3 and, in particular, we utilize a one-dimensional study to highlight informa-tion that is present within these nonconformal F E M s but was previously obscured. For example, in [36] a nonstandard Radau quadrature is suggested for the treatment of the mass matr ix , however we show that a simple lumping [69] is all that is required in the polynomial case [Section 5.2.2], and moreover, the exact evaluation is necessary i n the analytic case [Section 5.2.3]. The primary affect of this lumping procedure is to restore the cell-based balance relation. In addition, motivated by the well known equivalences between some nonconformal and mixed-hybr id finite element methods [14], Hennart et al . [39] have also studied the equivalence between mixed-hybr id F E M s and various nodal 95 Chapter 5. Finite Element Methods 96 methods. In Section 5.3 we demonstrate this equivalence for the polynomial case and discuss the difficulty associated with extending it to the analytic case. Specifically, we show that although it is possible to cast the const ant-const ant N I M in the indefinite form characteristic of mixed-hybr id F E M s , the appropriate basis functions are elusive. 5.1 P r i m a l C o n f o r m i n g F E M In its second-order form the finite element treatment of the one-group diffusion equation commonly employs a piecewise bilinear basis. The use of this s popular basis is motivated by two essential properties, namely, it is a low order polynomial that leads to a second or-der discretization and it is conforming. However, this basis performs poorly for problems wi th piecewise constant coefficients that exhibit severe variations. Thus, to highlight the strengths of the nodal discretizations, and hence the equivalent F E M s , we review this conforming case. 5.1.1 V a r i a t i o n a l F o r m u l a t i o n For completeness we briefly review the variational formulation and establish the notation that we employ i n the following sections. Consider the second order divergence form of the one-group diffusion equation in R n , - V - ( D ( r ) V 0 ) + E a ( r ) ^ = g(r ) Vr e ft with .D(r) > Dmin > 0 and E a ( r ) > 0. To avoid unnecessary complication we w i l l restrict this theoretical discussion to homogeneous Dirichlet boundary conditions, (b(r) = 0, Vr e dn It is well known (e.g. [15]) that the weak variational form of this diffusion problem may be written as: Chapter 5. Finite Element Methods 97 find d> Â£ V such that, a(<f>,v) = (DV<f>,Vv) + (Zad>,v) = (Q,v), Vu Â£ V (5.1) The space of admissible functions V = HQ(Q), H1^) = {v| \\v\\i < oo} #o(fi) = |u| u = 0 Vr Â£ dtt^f^H1^) where ||u||i = (||Vu||o+ IM|o)7 a n ( l IMIo = IMIz,2(fi)- The conformal approach to numer-ically approximating this variational problem is to replace V by some finite-dimensional subspace 14 (he. 14 C V). A basis for 14 is chosen, 14 = span{$i ,$ 2 , - - . ,$ jv} C V such that the approximate solution (f>h is given by, N 4>h(r) = EÂ«Â»*.-(r) i=l In the Galerkin approximation the corresponding linear system that determines the co-efficients a{ is obtained by direct substitution into (5.1), / N *. ^ j=i ' N = j = l,...,N (5.2) i=l which is a linear system AOL = b where the components of the load vector b are simply, while the matr ix A = A ^ + A ^ is decomposed into the stiffness matr ix A ^ and the mass matr ix A ^ wi th elements defined by Chapter 5. Finite Element Methods 98 5.1.2 Conforming Polynomial Basis To provide the necessary background for the one-dimensional discussions i n Sections 5.2.2 and 5.2.3 we recall the one-dimensional conforming piecewise linear F E M and its properties. The well-established extensions to a two-dimensional tensor product mesh are noted at the end of this discussion. We begin by stating the global piecewise linear basis functions, 1 Ax; x â€” x G ft; â€¢(x;+| â€” x) x e ft;+i (5-3) Axl+1 v ^ 0 x e ft/Aft; which leads to the following discretization. Definition 5.1.1. The piecewise linear FEM for one-dimensional one-group diffusion is defined by the following stencil, - C ; ( / > ; _ i + (di + d;+i)</>;+! - C ; + i 4 + 3 = / ; + fi+i (5.4a) where, Di r 2 r r ^ n 2 l . Di r 4 r - ^ i 2 l â€ž Axi Ci A and, \ _ . / ( S a ) ' , i \ ( - ) _ A . . A r . \W _ was defined i n (3.47). This is a well-known second-order discretization (see, e.g. [69]), for which the variational form implies that the discrete error eh(x) = 4>(x) - (f>h(x) Chapter 5. Finite Element Methods 99 is minimized over the tr ia l space Vh in the H1 (ft) norm. However, several relevant physical properties of the discretization are seldom referenced, hence we recount them here. To this end we first recall the definition of an M - m a t r i x given in Varga [70] (pg. 85). Definition 5.1.2. A real nxn matrix A = (ctij) with a; j < 0 for al l i ^ j is an M-matrix if A is nonsingular, and A~l > 0. Here the notation A - 1 > 0 indicates that the individual elements of the inverse are nonnegative, and thus, if A is an M - m a t r i x and the source Q(x) > 0 for al l x â‚¬ ft then a nonnegative (and hence physical) solution 4>h(x) > 0 for al l x Â£ ft is guaranteed. The following Corollaries of Varga [70] (pg. 85) are employed throughout our discus-sion to prove that certain global matrices A satisfy the above definition. Corollary 5.1. If A = (a;j) is a real, irreducibly diagonally dominant nxn matrix with < 0 for all i ^ j , and a;,; > 0 for all 1 < i < n then A~l > 0. Corollary 5.2. If A = ( a t J ) is a real, symmetric and nonsingular n x n irreducible matrix, where a,-^ - < 0 for all i ^ j , then A - 1 > 0 if and only if A is positive definite. From Corollary 5.1 it is t r iv ia l to conclude that A of Definition 5.1.1 is an M â€” matrix if [Aj-^] 2 < 2/3. In applications where the coefficients vary by orders of magnitude this restriction may be quite severe. Hence, even though Corollary 5.1 only gives a sufficient condition for A - 1 > 0, this breakdown may deter the use of this F E M discretization. The computation of the current is also required in a number of modelling problems. Typica l ly this is obtained through the application of Fick 's law to the finite element solution, (5.5) Chapter 5. Finite Element Methods 100 In this manner we find, JU ~ J+h = + (Â£: + i^t)^ -A^>+* and hence from Definition 5.1.1 it is apparent that the normal currents are discontinu-ous. Furthermore, the current is constant on each cell and hence it is not possible to satisfy a cell-based balance equation. Instead we find that the pr imal conforming F E M approximates a vertex-based balance relation. The common procedure of lumping the mass matrix [69] may be used to alter the properties of this discretization. In particular, lumping modifies the coefficients of (5.4a) to read, Axi Ax{ t J Ax{ Z 1 while the source term is unaffected. Then it is apparent from Corollary 5.1 that the global matr ix A is unconditionally an M-matr ix . Moreover, if the coefficients are constant on each cell , then the discretization satisfies the following vertex based balance equation, Ji+i - J i + ^ { A ^ ( Â£ A ) T + A x , - + i ( E 0 ) . + 1 } < k + i = ^ { A X , - Q , - + Axi+1Qi+1} (5.7) However, the lumping process also destroys the variational form, and hence we no longer minimize ||e/,,(x)||i over 14. We summarize these observations in the following proposition. Chapter 5. Finite Element Methods 101 Proposition 5.1. The piecewise linear FEM with the coefficients constant on each cell and exact integration of the stiffness and mass matrices is a second order discretization of the one-group diffusion equation that is characterized by the following properties 1. \\eh(x)\\i is minimized over the trial space Vh 2. cbh(x) > 0 for all x eft if [\\x)]2 < 2/3 for all i Â£ [1, L] 3- J + + h - J - H = 0(Ax) 4- A cell balance equation does not exist 5. A vertex based balance equation is approximated Lumping the mass matrix also produces a second order discretization which possesses the following properties 1. \\eh(x)\\i is no longer minimized over the trial space Vh 2. 4>h(x) > 0 for all x Â£ ft independent of the mesh 3. J + + h - J - + h = 0{Ax) 4- A cell balance equation does not exist 5. The vertex based balance equation (5.7) is satisfied exactly Thus, it is apparent that a lumping of the mass matrix (i.e. the lower-order terms) may be used to alter certain properties of the discretization. Specifically, we demonstrated that one may sacrifice the variational principle i n favour of certain physical properties, such as posit ivity of the F E M solution and a balance relation. Ideally, however, such a sacrifice should be unnecessary, and indeed all of the desired properties are obtained wi th the use of analytic basis functions (see Section 5.2.3). Final ly , we note that the extension of these observations to the two-dimensional case on a tensor product mesh is, apart from the vertex-based balance relation, quite natu-ral . The control volume for the lumped approximation of the bilinear finite elements is nonstandard and is discussed in some detail in [19]. Chapter 5. Finite Element Methods 102 5.2 P r i m a l N o n c o n f o r m i n g F E M In Section 5.1.2 it was demonstrated that the popular conforming piecewise bilinear basis is not necessarily the best choice, particularly if the continuous problem has piecewise constant coefficients which jump by orders of magnitude. In many situations it has been demonstrated that relaxing the conformity restriction permits the choice of more physically motivated basis functions which are ult imately more accurate. We w i l l demon-strate that this is precisely the case for the nodal discretizations whose corresponding basis functions are based on average quantities, and hence, enforce continuity i n only a weak sense. To simplify the development and to highlight the underlying influence of transverse integration we first derive the discretizations in I D , and then demonstrate their natural extension to a tensor product mesh in 2D. 5.2.1 T h e E d g e a n d C e l l Basis The underlying commonality between the various nodal methods is the choice of cell- and edge-based unknowns. To develop primal nonconforming finite element methods that are equivalent to specific nodal methods we first introduce a notation suitable for the cell and edge bases. O n e - D i m e n s i o n a l Bas is The distinction between cell and edge unknowns in ID simply states that for each cell, centred at a;,-, we have an associated unknown, </>;, which is defined as, while the edges are degenerate, being simply point values, 4>i+i = (f>(xi+i), i = 0 , . . . , L . Chapter 5. Finite Element Methods 103 The implicat ion of this degeneracy is that these methods are actually conforming i n I D . Considering the canonical element f2c = [â€”1,1] we adopt the following notation, where the superscript lk indicates the moment that normalizes this basis function and the subscript (k defines at which point in Qc the function is normalized to 1. In the event that either script does not apply to a particular basis function, an asterisk " * " is used. Moreover, al l other relevant moments and point values of a particular basis function are set to zero. Al though this notation may seem excessive in I D , it is a notation that extends well to higher dimensions and higher-order methods. The corresponding global Table 5.1: Canonical Definitions i n I D 4)(o basis functions are simply, otherwise otherwise Chapter 5. Finite Element Methods 104 where &(x) = 2(x â€” x ; ) / A x ; . ' Thus the global expansion of <f>h(x) is given by L - l (x) = ^ ^ + i* % w + E^.-0)(a;) i= l W i t h the definitions of the basis functions on the canonical element it is possible to perform an element by element construction of the stiffness and mass matrices (cf. [15]). Specifically we consider the contribution from the element $7; which is described by the following 3 x 3 block, ' s , - ( ^ J , ^ W ) ) s , - ( ^ Â° j , ^ } ) (5.8) where Similar ly we have the elemental mass matr ix , "m.-(<; ,<!) M,- = m ^ r ^ f - i ) ) m . -(^(i) .^(i)) (5.9) where, m Hence, the total elemental matrix is A ; = S; + M t . F inal ly the elemental contribution to the load vector is b4 = Ax t -( O ^ f i J ) which, under the assumption of a piecewise constant source gives, b, = [<5;,0,0 T Chapter 5. Finite Element Methods 105 Two-Dimensional Basis The two-dimensional cell- and edge-based unknowns were introduced i n Chapter 3, and provide 5 degrees of freedom on the canonical element flc = [â€”1,1] x [â€”1,1]. The corresponding extension of the one-dimensional basis function notation is summarized i n Table 5.2. In addition, the two-dimensional bases that we consider are unions Table 5.2: Canonica Definitions i n 2D ^ ( - 1 , 1 7 ) ^ = 1 ,(o,o), (*.*r \s. u \ \[j[^{i^)didr] = \ of one-dimensional bases. This property significantly simplifies the development and implementation of these methods and is stated in the following Corollary. Corollary 5.3. The lowest-order two-dimensional cell and edge basis functions that we consider span a union of one-dimensional spaces, and hence, may be written, = <Â°J(0+ < } ( Â» / ) - i (5.10a) = <J(0 (5.10b) = < \ ) ( 0 (5.10c) = 1>$(v) (5.10d) (5.10e) where the one-dimensional basis functions are characterized in Table 5.1 and span the corresponding one-dimensional space. Chapter 5. Finite Element Methods 106 Proof. A two-dimensional canonical basis function may be written and hence the correspondence with cell and edge unknowns yields the following under-determined system, a 0 1 0 0 1 0 0 1 " Â« i 0 1 0 0 0 0 1 a 2 0 0 1 0 0 0 1 b0 0 0 0 0 1 0 1 bi 0 0 0 0 0 1 1 b2 c 0 = ek where = [5k,i, &k,2, o~k,3, Â°~k,4, ^fe.s]T, and 5k,i is the Kronecker delta function. It is a simple matter to verify that (5.10) is a solution. â€¢ U t i l i z i n g (5.10) the global basis functions are simply, < i [&(*)] <1 [%(*/)] ^(-i)fe+i(y)] (5.11a) (5.11b) *S 0 ) (* ,y) (5.11c) where { Â£ i ( x ) , ^ ( y ) } define the affme mapping of onto Clc. The global expansion is L-l M L M - l L M M * , y) = E E y) + E E y) + E E ^ *S , 0 ) ( *> i=l j = l =1 i = l *'=1 J=l Chapter 5. Finite Element Methods 107 The elemental construction of the stiffness and mass matrices is now a t r iv ia l extension of the I D case, wi th the elements of the elemental stiffness matrix S,j given by, % Â«/ â€” I t / â€” 1 and similarly the elements of the elemental mass matrix are given by, We note that the evaluation of these integrals is greatly simplified wi th the relationships presented in (5.10) 5.2.2 P o l y n o m i a l Bas is F u n c t i o n s The objective is to demonstrate that a simple tensor product polynomial basis may be used to derive a nonconformal finite element method which is equivalent to the N E M dis-cretization. Recall that N E M and N I M are equivalent for conservative diffusion. Hence for this case we investigate the implicat ion of this additional equivalence on the interpre-tation of the nonconformal F E M . A O n e - D i m e n s i o n a l S t u d y Although nodal methods lead to conforming finite elements i n the simple I D case, the use of transverse integration in modern nodal methods effectively decouples the coordinate dependence and hence most treatments are essentially one-dimensional i n nature. Thus, a study of the one-dimensional case is an informative starting point, in which the basis functions and integrals are significantly less complicated. Chapter 5. Finite Element Methods 108 Definition 5.2.1. In one dimension the cell and edge quadratic basis functions on the canonical element Vlc may be written, (5.12a) (5.12b) (5.12c) and are displayed in Figure 5.1. The corresponding global quadratic basis functions are given by, 1 - - ^ { X i + L -x) + (Xi+l ~ 1 -Axi 4 Ax, +i â€¢ ( x - x i + i ) + Ax2 +i -(X x e f2;+i x e nh\{ni{jni+1} 6 Ax 0 2 V ^ ' + j X )(x-Xi_%) X^tti x e o^\Oj (5.13a) (5.13b) Remark 5.1. Basis functions such as $jÂ°^(a;) are often referred to as bubble functions as their support is restricted to a single cell and hence their appearance is bubble-l ike. Derivation of the ID Polynomial Basis. Wri t ing a general quadratic as the correspondence wi th cell and edge unknowns yields the following linear system gov-erning the expansion coefficients, 1 - 1 1 a0 1 1 1 = ek 1 0 1 3" 0.2 Chapter 5. Finite Element Methods 109 1.0 -0.5 0.0 0.5 (a) (b) (c) Figure 5.1: The quadratic canonical cell and edge basis functions. where e& = [<Â£fc,i, b~k,2, ^ , 3 ] â€¢ W i t h a l i tt le help from MAPLE [17] we obtain the canonical basis functions given in (5.12). If we note that, 2 A 2 'x-Xi) = l - - ^ ( x i H - x ) Ax. (x - X{_L) - 1 then the global basis functions (5.13) follow. (5.14) â€¢ Being low order polynomials the exact evaluation of the elemental stiffness and mass matrices (5.8) and (5.9) respectively, is simple and gives, S,- = 2 Di - 3 - 3 Axi Furthermore recalling that 2 1 1 2 M,-( S a ) . A x j 12 5 1 5 1 5 1 4 1 5 15 15 1 1 4 5 15 15 -(5.15) X i = 2-^-n(-)i 2 Axi we find that the total elemental matrix is simply, A,-Ax; 12(1 +fp!*)]2} - e { i + ^ [ A n 2 } -e{i+^[An2}" 4 { l + MA^] 2 } . - 6 { l + MA? 0] 2} 2 ( 1 - M W } 4{i+^[An 2}. (5.16) To highlight the properties of this discretization we present the corresponding stencils. Chapter 5. Finite Element Methods 110 Definition 5.2.2. The cell- and edge-based polynomial FEM for one-dimensional one-group diffusion is defined by the following local cell-based equation, and an edge-based stencil of the form, Definition 5.2.3. The reduced edge-based polynomial FEM for one-group diffusion is given by the following edge-based equation, + (di + di+^j4>i+L - ci+1(bi+i = fi + fi+i (5.19a) -Ci9i-L where. AxA 1 + | [ \ W ] 2 (5.19b) / i = 2 { l ^ P ' ( 5-1 9 d ) Derivation of the ID polynomial FEMs. The coefficients in equations (5.17) and (5.18) are simply read from the total elemental matrix (5.16). The support of ^ ^0\x) is 0 , and hence cbi may be eliminated locally. Solving for the cell average in (5.17) gives, + \ A x 2 h = A â€”-rf + 0 , - 1 + 7 ' , 2 1 Qi (5.20) 2 { 1 + I W ] } 1 2 A { l + I [A^] 2 } and hence, substitution into (5.18) yields (5.19a). â€¢ Chapter 5. Finite Element Methods 111 Lemma 5.1. The global coefficient matrix A of the ID reduced edge-based polynomial FEM (5.19a) is an M-matrix. Proof. The coefficients c,- are strictly positive, wi th the numerator attaining an absolute m i n i m u m at AJ3^ = \/2 of j | > 0, and hence the elements [A]itk < 0 V I ^ fc. In addition for \\x) > 0 we have di > Ci hence A is diagonally dominant. Thus by Corollary 5.2 it is an M - m a t r i x . â€¢ Obtaining an M - m a t r i x , independently of the off-diagonal terms that arise from the mass matr ix , appears to be a significant improvement over what is obtained wi th the piecewise linear basis. However, the piecewise quadratic basis functions are not nonneg-ative, and hence, we obtain unconditional positivity of only the point values, (f>h(%i+^) for i = 1 , . . . ,L â€” 1. In addition, upon evaluating the current (5.5) it is apparent that (5.17) is not the balance equation over $7;. Moreover, the current is not continuous at the cell edges. Hence, the quadratic basis does share some of the physical deficiencies that appeared previously with the piecewise linear basis functions. Ult imately , however, we seek to establish equivalences between certain nodal methods and certain F E M methods, and for this both the current continuity and the balance equation are essential. Thus, we consider the approximate treatment of the elemental mass matr ix . In particular, lumping the mass matr ix gives, 2 0 0 ' 0 0 0 0 0 0 x ~ ( S a ) . Axi Â£(M,-) = Mt = V %' Hence the total elemental matrix becomes, 1 2 { l + |[AP]2} A t- = S, + M i Axi 4 2 2 4 (5.21) Chapter 5. Finite Element Methods 112 such that the resulting discretization not only preserves the cell based balance equation but also restores the continuity of the current. However, the corresponding coefficients of the reduced edge formulation become, Di f l - l [ A t W ] 2 , 1 ) (5.22a) 1 Axai + ip i ' ) ] 2 and hence the lumped global matrix A is not an M - m a t r i x unless [ A | ^ ] 2 < 2/3. Thus we have exchanged the property that A is unconditionally an M - m a t r i x , for the adherence to cell balance and current continuity. This is in distinct contrast to the piecewise linear basis, for which lumping gives an M - m a t r i x unconditionally but produces an edge- as opposed to cell-based balance relation. We also note that this one-dimensional cell- and edge-based F E M (5.21) foreshadows the equivalence with N E M that we w i l l verify i n two dimensions. Two alternative treatments of the mass matrix have appeared i n the literature. Specif-ically, the nonstandard Radau quadrature suggested by Hennart [30, 36] also yields M;. The work of Arbogast and Chen [3] suggests that the required approximation can be viewed as a projection, m, (MO-MWa')) In this case the projection operator V () is simply the integral average. A treatment of the mass terms is entirely avoided with conservative diffusion, and moreover, the previously obscured correspondence of the cell- and edge-based F E M with N I M becomes evident. For completeness we first state the conservative forms of Definition 5.2.2 and Definition 5.2.3 (i.e. \\x) 0). Chapter 5. Finite Element Methods 113 C o r o l l a r y 5.4. The cell- and edge-based polynomial FEM for conservative diffusion yields the following flux based balance equation, - 6 - ^ _ L + 1 2 - ^ - & - 6 ^ ^ + i = AxiQi (5.23) Axi 2 Axi Axi 2 while at the edges x = x i + i we have, 2 - ^ - 6 - ^ - + + - 6 - ^ - ^ + 1 + 2 - ^ + } = 0 A x ; 2 Ao?; \ A a ; i A a j j + i / 2 Axi+i Axi+1 2 C o r o l l a r y 5.5. The reduced edge-based polynomial FEM for conservative one-group dif-fusion is given by the following edge-based equation, Di . f Di Di+l \ , A+i . 1 . ^ . 1 A x . - - " * + V A X T + Ax^Tj^+* " A ^ 7 * + * = 2 A X ^ + 2 * X i + l Q i + 1 ( 5 - 2 4 ) L e m m a 5.2. The finite element solution 4>h(x) obtained with the one-dimensional poly-nomial basis [Definition 5.2.1] is the exact solution of the one-dimensional, one-group, conservative diffusion equation if both the diffusion coefficient D(x) and the source Q(x) are constant over each cell. Proof. The exact solution is piecewise quadratic and is, therefore, included in the tr ia l space. Since the variational formulation chooses the optimal H1 (fi) element of the tr ial space, the F E M solution w i l l also be the exact solution. â€¢ Remark 5.2. It is straightforward to obtain the exact solution of the problem defined i n L e m m a 5.2 in explicit form: M x ) = ^ A s " * ^ - * + ^ A I ! " ^ + * + 2 ^ ( X " ~ X ) Q i ( 5 ' 2 5 ) where < x < This expression (5.25) clearly demonstrates that having Qi > 0 for i = 1 , . . . , L, and (f>h{xi+%) > 0 for ? = 1 , . . . , L â€” 1, are the necessary and sufficient conditions to have 4>h(x) > 0 V x â‚¬ fi. Therefore, having unconditionally that A is an M -matr ix is sufficient to ensure that for Q(x) > Owe obtain, unconditionally, a nonnegative F E M solution. Chapterd. Finite Element Methods 114 E q u i v a l e n c e i n T w o D i m e n s i o n s We extend the one-dimensional analysis and demonstrate the equivalence between the following cell- and edge-based nonconforming finite element method and the N E M dis-cretization. The two-dimensional cell and edge basis functions on the canonical element are defined by (5.10) using the one-dimensional basis functions of Definition 5.2.1. The corresponding global basis functions are defined according to (5.11). This relationship wi th the one- dimensional basis greatly simplifies the evaluation of the elemental stiffness and mass matrices, reducing most integrals to those of the previous section. It also gives a clear indication of the relevance of the one-dimensional study. Evaluating the elemental stiffness matrix we obtain, 6 Ax<y) - l 3 r Â« , i 3 r t j 3r,j 3 r * ' , i ri,j 0 0 â€”3r,-j 0 0 â€”3fjj 0 0 2f ; j ri,j - 3 r \ i 0 0 ?r- â€¢ (5.26) where g\Xj'v^ = AxiAyj/(Ax2 + Ay2) is given i n Definition 3.4.4, and we have defined, r 8 J - = Ayj/Axi and f 8 ) j = Axi/Ayj. Exact integration of the elemental mass terms gives, 14 1 1 1 1 5 5 5 5 5 1 4 1 0 0 5 15 15 1 1 4 0 0 5 15 15 1 0 0 4 1 5 15 15 1 0 0 1 4 5 15 15 -(5.27) Observe that there is no direct coupling of horizontal and vertical edges. In addition, the failure to satisfy the cell balance equation (3.8) and the continuity of normal currents first Chapter 5. Finite Element Methods 115 identified i n the I D case is also found here. Thus, equivalence wi th the nodal expansion method relies on an approximate treatment of the mass terms, and i n particular lumping gives 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 . o "0 o â€¢ 0 0 . (5.28) and leads to the following conclusion. Theorem 5.1. The cell and edge polynomial basis yields a nonconforming finite ele-ment discretization of the one-group diffusion equation that is equivalent to the NEM discretization [Definition 3-4-4] if the mass matrix is lumped. Proof. The coefficient of <pij i n the flux based balance equation is simply, \j ) x + -py>:v> )> = \2D Q Â« J (x,y) A l l other coefficients are read directly from (5.26) and thus, dividing through by A x ; A y j , equivalence wi th Definition 3.4.4 follows. â€¢ The properties of the exact and lumped discretizations parallel those of the one-dimensional case. However, there is one significant exception: A is not an M - m a t r i x . Yet these discretizations have been used extensively in reactor modell ing, an application i n which it is required that 4>h(x,y) > 0, wi th great success. Therefore, we investigated the inverse of a small constant coefficient conservative Dirichlet problem using M A P L E [17] and found that indeed A~* > 0. We conjecture that with the exact evaluation of the stiffness and mass matrices, A-1 > 0 unconditionally, while for the lumped case Chapter 5. Finite Element Methods 116 a condition similar to that found in I D exists. Note that nonhomogeneous boundary conditions could sti l l produce a nonphysical solution, but in many practical cases the B C s are homogeneous. 5.2.3 A n a l y t i c Bas is F u n c t i o n s The motivation of enforcing nodal balance through the naturally arising edge and cell based unknowns is taken one step further in analytic nodal methods. Specifically, we recall that, under the assumption of piecewise constant coefficients and source, the trans-verse integrated O D E s are solved exactly. In terms of the nonconformal finite element methods this suggests that the basis be spanned by the homogeneous and particular solutions of these transverse integrated O D E s [37]. Thus, we note that i n the case of conservative diffusion (i.e. S a = 0), the polynomial and analytic bases are equivalent independent of dimension. A O n e - D i m e n s i o n a l S t u d y Following the treatment of the polynomial F E M s we begin wi th an investigation of the analytic basis functions i n one dimension. D e f i n i t i o n 5.2.4. In one dimension the cell and edge analytic basis on the canonical element Q,c are given by, 1 <t } ls inh(Af>Q 1 2 [1-<*!*>] 2 s i n h ( A | * ) ) 2 1 8& l s i n h ( A t W Q 1 2 [1 - 2 s inh(Af >) 2 1 CO8h(A!*>0 (5.29a) 1 c o s h ( A ^ 0 (5.29b) (5.29c) w. here we recall 6^ = [\\x>]_1 t a n h ^ ) . Chapter 5. Finite Element Methods 117 Derivation of the ID Analytic Basis. We consider the basis {1, s i n h ( A ^ Â£ ) , cosh(A t-^Â£)}, which contains the homogeneous and particular solutions of the piecewise constant O D E on the canonical element Q,c. Thus a general canonical basis function is writ ten, rb(0 = a0 + a i smhiX^O + a2 cosh(A^Â£) yielding the system, 1 -sinh(Ai*>; 1 sinh(Aj 3 ; ) ) 1 0 A c o s h ^ ) a0 c o s h ( \ W ) ai t l - h ( ^ ) Â«2 W i t h a l i t t le help from M A P L E [17] we obtain the explicit formulae given i n (5.29). â€¢ Evaluating the integrals of the stiffness matrix we obtain, Axi (x,s) li 2 i 1 (x,s) â€” -co 2 i 2 i 1 (x,s) â€” â€”CX-2 i I | Q ( M + / ? ( M | where we have conveniently used the following definitions, {sinh(AfVsM\(V\W} (x,s) (s) <A - i i [ A tw ] 3 -{ \ W c o s h ( A J W ) - s i n h ( A f ) ) } ' { s i n h ( A 8 W ) c o s h ( A t W ) + A W } sinh2(Af>) Similarly, the elemental construction of the mass matrix may be expressed i n the form, Mi = ^ 2 (x.m) Ti 1 (x,m) â€” -a-2 2 1 (x,m) 2 i 1 (x,m) 2 i HÂ«f'm ) + #m )} HÂ« ! x ' m ) -# ' m ) } 1 (x,m) 2 i Chapter 5. Finite Element Methods 118 w here, OL (x,m) sinh(A W i [r ] 2 -s inh 2 (A^) V - ) L % j - \ I / {*,<*> c o s h ( A ^ ) - s i n h ^ ) } {*,<*> c o s h ( A ^ ) - s i n h ( A ^ ) } ' (x,m) _ smh(\\x)) cosh(AP) - A (*)\ \(x) T (x,m) A\<*>sinh2(A,<*>) 2A^ cosh(AP) + { [ A P ] 2 - A P c o s h ( A ^ ) s i n h ( A ^ ) } { A,<*> c o s h ( A ^ ) - s i n h ( A ^ ) } { A < * ) cosh(A<*>) - s i n l ^ ) } 1 hence the total elemental matrix may be written, 2i_ Ax; -2aS x ) -2Â«S*) -2aS x ) {Â«!X) + ^ } . (5.30) where, (Â«) = I f (*..) + [ A W a ^ ) = M W 0 . ri" = i ( : ^ , . ' ( x > s ) i [ \ ( x ) ~ \ 2 ( x , m ) \ T = 2 ^ â€¢+ iXi J Ti J = [1 -Although the elemental matrices provide a compact representation of the discretiza-t ion they obscure many of the interesting properties. Thus, in the following discussion we first present the corresponding stencils. Chapter 5. Finite Element Methods 119 D e f i n i t i o n 5.2.5. The cell- and edge-based analytic FEM for one-dimensional, one-group, diffusion is defined by the following local cell-based balance equation, Di -a D LA*). Axi"1 " 'Axi and edge-based stencils of the form Ax; -OL. (5.31) Di Axi (Â«!r) - /r) D. ^ ~ 2 A x LAX), â€¢OL, + (a) + # ') + T - â€” ( a , V i + A V i > A x ; A z i + i o ^.(g) A x i + i (5.32) Ci}+\(pi+l + Di+1 = 0 [Axi+i D e f i n i t i o n 5.2.6. The reduced edge-based analytic FEM for one-group diffusion is given by the following edge-based equation, i<f)i-L + (di + di+1^J<f)i+L - C i + i 0 i + | = /,â€¢ + fi+l where, Di_ Ax D Di A x , 8 inh(AJ a : ) ) Di A f ) c o s h ( A J - ) ) A x , s i n h ^ ) /,â€¢ = \8\x)AxiQi (5.33a) (5.33b) (5.33c) (5.33d) Derivation of the ID analytic FEMs. The coefficients i n equations (5.31) and (5.32) are simply read from the total elemental matrix (5.30). The support of $\Â°\x) is fi; and hence d>i may be eliminated locally. Solving for the cell average i n (5.31) gives, <k = 5 *\ + p|, ) ] 2 4 D i (5-34) and hence substitution into (5.32) yields (5.33a). The simplified form of the coefficients is obtained by applying the definitions of a\x\ (3^ and S\x^ along wi th a few hyperbolic trigonometric identities. â€¢ Chapter 5. Finite Element Methods 120 Similar ly to L e m m a 5.1 we obtain Lemma 5.3. The global coefficient matrix A of the ID reduced edge-based analytic FEM (5.33a) is an M-matrix independent ofX.\xK Without any special treatment of the mass terms the analytic basis functions uncon-dit ionally yield an M - m a t r i x . In addition, Fick's law (5.5) verifies that the cell balance equation is satisfied and that the currents are continuous over the entire domain. Thus, we have attained both the physical and mathematical objectives of the discretization. This is readily understood with the consideration of the following L e m m a . Lemma 5.4. The finite element solution 4>h(x) obtained with the one-dimensional an-alytic basis [Definition 5.2-4] z s the exact solution of the one-dimensional, one-group diffusion equation if the coefficients D(x), S a ( x ) and the source Q(x) are constant over each cell. Proof. A s i n L e m m a 5.2 the tr ial space includes the exact solution. â€¢ Remark 5.3. It is possible to show that the exact solution to the problem defined in L e m m a 5.4 may be written as smh(\\x)[xi+i - x]) smn(A) ') s'mh(X\x)[x - Xi_i}) f cosh ( A ^ 0 [a: -Xj])\ Q{ s inh(A 8 W ) â€¢ H I cosh(A J W ) )*i where < x < x i + L . This provides a straightforward verification that having Qi > 0 for i = 1,. . . , L, and </>/j(x8+i) > 0 for i = 1 , . . . , L â€” 1, are the necessary and sufficient conditions to have d>h(x) > 0 V x G f i . Therefore, having unconditionally that A is an M -matr ix is sufficient to ensure that for Q(x) > Owe obtain, unconditionally, a nonnegative F E M solution. Thus, we see that the analytic basis functions yield a discretization that <f>h{x) Chapter 5. Finite Element Methods 121 captures the desirable physical and mathematical properties, even i n the nonconservative case, while the polynomial basis does not. F ina l ly we examine the conservative l imit of the total elemental matr ix A to verify that it corresponds to the total elemental matrix of the polynomial F E M (5.16). We also demonstrate that a definite relationship exists between the analytic and polynomial F E M s . Specifically, performing the Taylor expansion of the elements of A; (5.30) we find, and hence truncating the expansions at second order in At- gives precisely the elemental matr ix obtained with the polynomial basis functions (5.16). Furthermore, evaluating the l imi t \\x) ->â€¢ 0 we obtain identical systems from the polynomial and analytic F E M s . Equivalences in Two Dimensions Following the polynomial case [Section 5.2.2] we now extend the one-dimensional analysis and demonstrate the equivalence between the following cell and edge-based nonconform-ing finite element method and the N I M discretization. The two-dimensional cell and edge analytic basis functions on the canonical element are defined by (5.10) using the one-dimensional basis given in Definition 5.2.4. Once again the corresponding global functions are defined according to (5.11). The relationship between the one- and two-dimensional basis functions that was first uti l ized in the polynomial case facilitates the exact evaluation of the elemental stiffness and mass matrices. Specifically, ut i l iz ing the !(Â«!â€¢> + flf")Â»{2 + 1 [*<â€¢)]' - J j p h V opi"]6} i(oi-Â» - fl") Â« {l - i [ A Â« ' > ] 2 - Â§-5[Xtr + 0[X<-Â»]'} Chapter 5. Finite Element Methods 122 I D results the stiffness matrix may be written, -x-r 2 (* i,3ai (a;,s) - i f 2 ' â€¢ a{y's) â€” i f - Q ^ ' S ) 2 ' _ I r . (a;,s) I r . 3^t,3 *) r-,(x,s) i^i,3 0 0 S- â€” 2D â€¢ _ I r . I r "(x> a) Cl(x<s) 0 0 2'* 0 0 l 4 i,3 _ I f . L 2 ' 1 0 0 1 4 (y>s) i,3 (5.36) where we have defined, 1 (s) (x,s) . - c\(*>s) (*i s) 1 /Q(*>S) A Â» J ~ ri>3~/i,j + ri<3^i,j 1 " i j â€” + J with * = x, or * = y. S imilarly the mass matrix is given by ct h3 ^h3 (Ej.^.Ax.-At/j A( m ) 1 (x.m) 1 (x.m) â€” â€” a - â€” - Q ! -2 J 2 t i â€ž ( v . m ) i,Jy>m) 2 i 2 t 1 (x,m) 2Uli,3 1 (x,m) 2ULi,3 lr\(x,m) l ^ ( i , m ) 4 U i , j 4 " â€” 1^(^,771) lr\(x,m) 4 l - J Â« , j 4 0 0 0 0 2 i 1 (y,m) 2 t 0 0 0 0 l A ( j , m ) 4 i,3 4^i,3 lz:(y,m) I A M 4 , _ , * J 4 W Â« , i (5.37) where we have defined, A ( m ) _ ^ . " 0 i,3 iiJ 1 . > > m ) 1 lt,J ,m) 1â€”(*,m) - A Chapter 5. Finite Element Methods 123 F ina l ly the total elemental matrix may be written, w here -2r,-jaJ*) - 2 r - o ^ " r- -Q{x) r . .=â€¢(*) 0 0 - 2 r i , 3 s a i , j r. .â€¢=â€¢(*) 0 0 (5. - 2 f â€¢ a{y) LI y U , 0 0 f- r. -(Â») 'l'3^i,3 0 0 =. .=â€¢(Â») â€¢=â€¢(*.Â»)' . ' 1 - lÂ»\.y which tr iv ia l ly yields the following theorem. Theorem 5.2. The cell and edge analytic basis yields a nonconformal finite element dis-cretization of the one-group diffusion equation that is equivalent to the constant-constant NIM discretization [Definition 3-4-4]-We note that this equivalence is stronger than that stated by Hennart [37] as it doesn't require the application of a nonstandard Radau quadrature to the mass matr ix , and hence, assures that ||e ,^(a;)||i is minimized over the tr ial space 14- In addition equivalence indicates that the continuity of the normal current and the cell-based balance equation follow. However, we do not have an M - m a t r i x . Similar to the conjecture in Section 5.2.2 our one-dimensional analysis suggests that A^1 > 0 unconditionally. Unfortunately, thus far we have been unsuccessful i n proving this conjecture. 5.3 Mixed Hybrid Finite Element Methods A n apparent weakness of the pr imal approach is that it abandons the rigorous treatment of cell balance which is central to the nodal view. In contrast a mixed F E M treats Chapter 5. Finite Element Methods 124 the first order form (2.12), viewing Fick's law (2.12b) as a constraint of the balance equation (2.12a) and hence explicit ly enforcing cell balance. Hybrid F E M s are a special class of mixed F E M s which temporarily relax certain inter-element continuity conditions, ult imately enforcing them in the weak variational sense. Thus, the hybridization of a mixed method captures the essence of nodal ideology and strongly suggests that an equivalence with N E M and N I M exists. Indeed such an equivalence has been established for the lowest order N E M by Hennart and del Valle [39] and we outline their results i n the following section. A discussion of the issues surrounding a potential equivalence wi th N I M is also presented. 5.3.1 Variational Formulation The typical (e.g. [5, 3]) introduction of mixed-hybrid F E M s is facilitated by first present-ing the mixed formulation followed by its hybridization. However, a stronger connection to nodal discretizations is apparent if we define the hybridization of the first-order system (2.12) and then proceed directly to the mixed-hybrid variational formulation. To this end we denote the set of edges E of the rectangles ftjj as Eh- Then the sets of boundary and interior edges are defined by Â£ f = { E Â£ E\ \ E C 3D,} and Â£Â° = Eh\EÂ®, respectively. The strong formulation of the corresponding hybrid one-group diffusion equation is given by, V nhl e n f c (5.39a) V nitj G tth (5.39b) V E â‚¬ El (5.39c) subject to boundary conditions. To avoid unnecessary complications in the following [D(v)] ' j + V ^ O V - J + Â£o(r)0 = Q J â€¢ n \ E + - J â€¢ n|s_ = 0 Chapter 5. Finite Element Methods 125 theoretical discussion we consider homogeneous Dirichlet boundary conditions, <f>(r) = 0, V r e dO We have adopted the standard ordering of these equations, wi th Fick 's law (5.39a) stated first, followed by the balance equation (5.39b), and finally the constraint that enforces the continuity of normal currents (5.39c). In addition, note that Fick 's law has been divided through by the diffusion coefficient so that the harmonic average w i l l arise from subsequent integration. In Chapter 4 it was demonstrated that i n one dimension this is the correct averaging, although no rigorous extension to two dimensions exists. More-over, this is the form of Fick's law [equation (2.10)] that was obtained with the Pi approximation of the transport equation [Section 2.2.1]. We define the following function spaces: . R-^vih) = {q| qe L2(n) x L2(n), v - q e L2{nhJ) v n t J e rth} Ro(tth) = {q |q â‚¬ R-i {Slh) , (q â€¢ n> Â£ + = (q â€¢ n) Â£_ VE G Â£ Â° h ) M_1(Â£h) = {u.\u.eL2(Â£h), p\E = 0 VEeÂ£dh] to facilitate the statement of the well known (cf. [14, 15]) weak mixed-hybrid variational formulation of (5.39) as: find (J h , 4>h, nh) e R{Qh) X U{Qh) X M(Â£l) such that [p-l3h,c{hj - V-q-fc)^. - {fih,<lh â€¢ n ) a n i j } = -(9><lh- n)- V qh e R(Slh) i,3 X>,Jfc-n)a n.. = 0 Vvhâ‚¬M(Â£h) where R C R-i, U C Â£/_i, M C M _ i are particular choices of finite dimensional sub-spaces. We further restrict the choice of subspaces such that the (weak) continuity Chapter 5. Finite Element Methods 126 constraint for the normal currents ensures that we find J/, Â£ Ro(^h)C H(div;Q), a conforming approximation of the corresponding mixed variational problem. Selecting a particular set of basis functions and performing the integration yields a system of the following form, A BT CT ' " Jh~ Qj' B -E 0 = â€” C 0 0 - Ph. 0 Like the mixed F E M , this system is indefinite; however, unlike the mixed F E M , A is block diagonal wi th each block corresponding to a cell f ^ j . The implications of this property on the sparsity of the reduced system are discussed i n detail i n Chapter 6. 5.3.2 P o l y n o m i a l Bas is F u n c t i o n s Recall the lowest-order Raviart-Thomas bases [60], which may be written RTl^ilH) = { q | ^ Â£ { { l , e } , { l ,77 } }Vn M Â£ n f c } f|/2_ ! uÂ°^(nh) = Â£ {i} v a j e n^n*7-* â€¢ AÂ£1(Slh) = {ii\Eeil} V F Â£ ^ 0 } f | M _ 1 The following theorem is due to Hennart [39]. T h e o r e m 5.3. The lowest-order Raviart-Thomas bases yield a mixed-hybrid finite ele-ment discretization of the one-group diffusion equation that is equivalent to the NEM discretization [Definition 3.4-4]-Proof. U t i l i z i n g the lowest-order Raviart-Thomas bases and assuming piecewise constant coefficients it is straightforward to perform the integration of the variational form exactly and obtain the following matrices. The matrix A is cell-based and hence has a block Chapter 5. Finite Element Methods 127 diagonal structure with each block itself being block diagonal. These blocks are given by, A i,3 AxiAy-j AD h3 48 0 A. (y) 4 W _ A(.V) â€” 2 1 1 2 (5.41) where the index k references the diagonal block of A. The matr ix B is also cell-based, and block diagonal, although the blocks are degenerate in the lowest-order case (i.e. [1 x 4]), Bi,3 k,k = [Ayj,-Ayj, A x , - , - A x , ] The last cell-based entry in (5.40) is E which in this lowest-order case is strictly diagonal wi th entries, E, 1,3 Final ly , the matr ix C is edge-based, and is composed of a two blocks corresponding to the horizontal and vertical edges, C â€” [ C ^ | C ^ ] T . Each of these blocks is bidiagonal, composed of the following cell-based 1 x 4 blocks, C g * + ) = [ - A y J 5 0, 0, 0] , C\p = [ 0, A y , , 0, 0] = [ 0 , 0 , - A x , , 0] , C%~] = [ 0 , 0 , 0, A x t ] where the Â± signs indicate which side of the edge lies within f i ; j . The block bidiagonal structure of C enforces the continuity of normal currents (3.24). Moreover, a direct evaluation of BJh â€” E<f>h = â€”Qh gives the balance equation (3.8). The definition of the one-sided currents (3.23) is obtained from AJh + BTcj)h + Cpn â€” 0 w i t h the diagonalization of A and the realization that the Lagrange multipliers are the edge unknowns of the N E M . â€¢ Chapter 5. Finite Element Methods 128 5.3.3 Analytic Basis Functions W h e n developing preconditioners for the nodal discretizations [Chapter 6] we make ex-tensive use of the indefinite formulation (5.40). Unfortunately, the equivalence of N I M with a particular mixed-hybrid F E M has not been established. Al though a strict equiv-alence would be ideal, for our purposes it is sufficient to cast the N I M system i n this indefinite form. Definition 5.3.1. The indefinite formulation of the constant-constant N I M discretiza-tion for nonconservative one-group diffusion [Definition 3.5.1] is given by (5.40) wi th , H 3 4 A , i 0 0 (i Â»>J J (5.42) where A: (x) 1,3 (x)n(x) (5.43) and A\yj is similarly defined. The matrices B, C and E are defined i n Section 5.3.2. t,3 Derivation of the Indefinite Formulation. We readily obtain the indefinite form of the constant-constant N I M discretization [Definition 3.5.1] from the one-sided net currents, (3.45). Specifically, we consider adding (a[xJ + B$) (3.45b) and (of] - B\XJ) (3.45a) to obtain, Thus mul t ip ly ing through by [As,-Aj/j] [4Z?jtj] gives the first row of A(x] defined m (5.43). S imilar ly adding (a\xJ - B\X)) (3.45b) and (a\x) + B\X}) (3.45a) yields the second row A\Xj given i n (5.43). Final ly , BJh â€” Ed>h = â€”Qh is simply the balance equation (3.8), while C J / , = 0 enforces the continuity of the normal currents. â€¢ C h a p t e r 6 R e d u c e d S y s t e m s a n d A p p r o x i m a t e Schur C o m p l e m e n t s A n important objective of this research is the development of fast iterative solvers for the nodal discretizations (equivalently the mixed-hybrid F E M discretization). One obstacle i n this pursuit has become apparent through the abstract formulation of the mixed-hybrid F E M s , namely the linear system is indefinite. Fortunately, this indefmiteness is readily circumvented with the block elimination of Jh in (5.40) which yields <t>h Ph '(4>.M) (6.1) with >(<M -SB BA-'C7 CA-'B? Sc BA-iQj - Q4 CA-'Qj where <5g = BA~lBT + E and Sc = CA~XCT. Thus far, the elementwise block diagonal structure of A , and hence A - 1 , has preserved the relative sparsity of the system. In fact, this is an abstract representation of the flux formulation that was given for N E M and N I M in Definitions 3.4.4 and 3.5.2 respectively. Thus, the exact sparsity structure of -5(0iM) is displayed in Figure 3.3. Perhaps more importantly, the Schur complement 5(0^) is symmetric positive definite indicating that it is an ideal candidate for iterative methods. In addition, the equivalence of this formulation wi th nonconformal F E M s was established in Section 5.2. However, the majority of reported works using the nodal, mixed-hybrid and nonconformal methods proceed to eliminate the cell-based unknowns {e.g. [38, 50, 18]}, so we consider this case first. 129 Chapter 6. Reduced Systems and Approximate Schur Complements 130 6.1 Lumping The Edge-Based Schur Complement Proceeding to eliminate the cell-based unknowns <ph from (6.1) leads to the Schur com-plement of the edge unknowns, S^h = Qs, (6-2) where SÂ» = C{A~1 - A~1BTSs1BA~l)CT Q 5 M = C(A-lBTSs1BA-1-A-1)QJ-CA'1BTS^Q4> Note that (6.2) is precisely the edge-based formulation that was defined in Definitions 3.4.5, 3.5.3 and shown schematically i n Figure 3.4. The relative sparsity of the original system (5.40) continues to be preserved in <S^ because SB is strictly diagonal (Figure 3.3). Yet , once again the most important property of this reduced system is that it is symmetric positive definite. In fact, it is this combination that has made (6.2) the most widely studied formulation of the mixed-hybrid finite element system {e.g. [50, 18]}. However, as we shall see, it is not without its disadvantages. To understand the inherent complexity of creating fast iterative solvers for (6.2) and to motivate our approach we partit ion the edge unknowns by orientation, vertical and horizontal, \ih = {u, v} such that (6.2) may be rewritten, ' A A u Qu Ayu Ayy V where Auu, Avv are tridiagonal and Auv is block bidiagonal, wi th each block itself bidiago-nal (Avu = A^v). A D I methods appear to be a natural choice and have been investigated by a number of authors [38, 8]. In addition, most standard preconditioners such as poly-nomial and incomplete factorizations have been studied [56]. However, none of these Chapter 6. Reduced Systems and Approximate Schur Complements 131 methods can attain the efficiency of multi-level solvers, the development of which has been hindered by the difficult hierarchy of grids and inter-grid transfer operators which must be defined for these edge-based stencils. Certainly multi- level solvers would be much easier to develop if only a single orientation of the edge unknowns remained. How-ever, the choice to eliminate either u or v is arbitrary and moreover it is apparent that the respective Schur Complements, AUU AyU (Ayy) AUV, Sy â€”~ Ayy Ayy (A yy ) Ayy suffer an equivalent loss of sparsity as a result of the appearance of (Ayy)-1 and respectively. Thus, posed purely as a problem in linear algebra, a suitable objective is to minimize the f i l l generated during the formation of Su (or <Sâ€ž). Further noting that each diagonal block of ( A u u ) _ 1 is pre and post mult ipl ied by bidiagonal matrices, the m i n i m a l f i l l attainable through approximations of ( A u u ) _ 1 results when it is diagonal. Determining the optimal diagonal approximation is a foreboding task. However, recalling that the best incomplete Cholesky factorizations preserve row sums it is reasonable to construct a diagonal approximation, Auu, which is composed simply of row sums of Auu. A simple inversion then gives (Auu)~l PS (Auu)~l. Hence the preconditioner may be written as S^-u â€” A â€¢ A A A (6.4) and the reduced system Sv = Avv â€” AVU(AUU)~1Auv is symmetric positive definite. In addition, we note that it is straightforward to prove that Sv is a second-order 9-point approximation of an ell iptic P D E . Thus, efficient mult igr id algorithms exist to "invert" it . In practice we use Dendy's black box mult igr id code [23], a robust solver which requires only the fine grid stencil as input. Chapter 6. Reduced Systems and Approximate Schur Complements 132 6.1.1 A Condition Bound We prove two results independently. The first is a proposition which makes the connection between the eigenvalues of the preconditioner and the difference between A U U and its approximation. Proposition 6.1. The eigenvalues of the'preconditioned system, X^JS^-u] 1S/J^j maybe written in the following form, X = 1, of multiplicity Nv \ = \(s*1su) = V+X(S^A:U) where Su = AUU â€” AUV (AVV) AVU and AUU = AUU â€” AUU. Proof. Observe that the inverse of the preconditioner <SM;u may be expressed as X = A! -1 A A A A - l \-vv A ( X T ' A u v . ] [S^u] -% = A â€” Auu. Hence the ). Since Su = Su + A*uu S-1 u -1 . r T â€” i - l C _ 1 A A~l \_Sv\ Avu \_AUU\ [<S^ ] uv. Block mult ipl icat ion and simplification give Su Su â€¢ Remark 6.1. Note that the preconditioner is symmetric positive definite, and hence may be factored: S^U = TTT. Thus, the preconditioned system is Q = JR~1SLXJR~T, and moreover, it is similar to [<SWU] 1SI1. Consequently, to bound K(Q) it is sufficient to bound the eigenvalues of [<SWâ€ž] 1SFI, and in fact, this is standard practice i n C G related estimates. Chapter 6. Reduced Systems and Approximate Schur Complements 133 Next we bound the condition number of the edge-based preconditioned system, i n -dependently of the mesh size, the diffusion coefficient and the absorption cross-section. We begin wi th the conservative case and consider the nonconservative case i n subsequent Corollaries. Theorem 6.1. Let r^j = A j / j / A x , - . Then the condition number of the preconditioned system with S a = 0 is bounded by K (K ] _ 1 S â€ž ) < m a x { ( l + r y , 3 } (6.5) With a uniform grid spacing r = rij, an even tighter bound is realized: ^[S^y'S,) < m a x | j ( l + r 2 ) , - ^ - ^ } (6.6) Proof. We employ a superelement analysis as proposed by Kuznetsov [50] to obtain a local bound on the condition number K. Consider first the elemental construction of the edge-based system by either uti l iz ing its equivalence with a nonconformal method or by elementwise el imination of the currents and the cell-based unknowns from the original system (5.40). O n a single cell fi;j we have (3 +a 13-a -(3 -(3 P - a P + a -(3 -j3 ~^p -p P + 7 p - 1 -P -P p - 1 P + j_ with a = Dijrij, p = 3Dijg\j'y\ 7 = A J ^ - J , fij = [r;j] 1 and where the unknowns on 0,-j have been ordered as [</>,_4> i+i j , <&j-i., 4>i,j+%] â€¢ The lumped approximation gives I" is straightforward to verify K, = ker(<S^) = ker( [S^-JS) = [1,1,1, l]Tâ€¢ This kernel is to be expected as <S^ essentially defines a pure Neumann problem on the cell , fi;j. F ina l ly 3 1 = As Al [AUU] Al As \AUU 2p 0 0 2p Chapter 6. Reduced Systems and Approximate Schur Complements 134 we consider the generalized eigenvalue problem, 3 > s = M S M ; U ] V B wi th ^ S J _ / C A p p l y i n g Theorem 6.3.1 of Rao and M i t r a [59] we find the that three "proper" eigenvalues are A s = [ l , 1, |(1 + r^)]. If nj > V2 A s <= [ l , |(1 + r?,-)] while if r i t i < y/2 A s G [|(1 + r ? - ) , l ] . Hence, if the mesh spacing is constant superelement analysis t r ivia l ly yields the bound given in (6.6). Whereas on a spatially dependent mesh we consider the following, ri,i > V2 V ( t , i ) ^ K < m a x | i ( l + r ? j ) } ritJ < \/2 V ( i , j) =^ K < rnaxl 3 1 < 3 r<.j l ( l + r i , , ) J I 1 + r1 1 ritj > V2 and rk,i < V2 K < max < )f - > < m a x { l + r2j] ri,j,rk,l I 1 -f- Tf_ j I ' Combining these bounds yields (6.5). â€¢ Remark 6.2. The form of the bound given in (6.5) is unchanged when the analysis is extended to a diagonal diffusion tensor. The crucial difference is the definition of r 2 J , . which is modified to read D*,3 Ax; (y) Hence we see that the implicat ion of poor conditioning on a spatially dependent high as-pect ratio grid is also indicative of poor conditioning in the presence of strong anisotropy. We now address the nonconservative case, for which the lowest-order N E M and constant-constant N I M are no longer equivalent. Corollary 6.1. The condition number of the preconditioned NEM system [Definition 3.4-5] for nonconservative diffusion satisfies the bound given in (6.5) of Theorem 6.1. Chapter 6. Reduced Systems and Approximate Schur Complements 135 Proof. Following the proof of Theorem 6.1 we consider a superelement analysis. The blocks of the elemental N E M system may be written, AL, = AL = 0 9 e e where 9 = ^Di,jQij'V^ [p\j'V^] 1 and Asvv is defined analogously. Ut i l i z ing the definitions of (ijV^ and we find that the lumped approximation gives, 6DiJriJ{l-\}i\'f]-1}+2e 0 0 6Dijrij{l-[i4'/]-1}+2d In contrast to the conservative case, the kernel is the null space, a property that follows from the well-posedness of the Neumann problem for E 0 > 0 and greatly simplifies the specification and treatment of the generalized eigenvalue problem. Specifically, we find that the eigenvalues may be written, 1 1 1 1 J 1 _ i _ rJA Thus the only non-unity eigenvalue is bounded by l < K < l ( l + rl) â€¢ and hence the bound in (6.5) follows Corollary 6.2. The condition number of the preconditioned NIM system for nonconser vative diffusion [Definition 3.5.3] is bounded by, f 1 0 , . . o x 10 Sr,u\ S^< m a x | â€” (1 + r . - J , â€” j (6.7) Proof. Following the proof of Theorem 6.1 we consider a superelement analysis. The blocks of the elemental N I M system may be written, " 9 9 As â€” 8jv) + o %Â°;v) + o AS â€”a* _ uv vu 9 9 Chapter 6. Reduced Systems and Approximate Schur Complements 136 where 9 = DijS^a^a^ 1 and Asvv is defined analogously. U t i l i z i n g the defini-tions of C,\Xj'v^ and -d^f-^ we find the lumped approximation gives, 2 ( A . r M [ l - ^ ) ] Â« S ) + ^) 0 0 2 ( A , ^ [ I - ^ V ! ? + ^ Hence the generalized eigenvalues are given by 1,1,1, 1 - 8{x)8{y) *,3 Â« , J To bound the non-unity eigenvalue we first note that X\~ \s(\\Xj , ritj). It is then straightforward to show that, ( j = and hence A s = ^ < A s (2,0) < A s < A s (0, r) = i (1 + r M ) and hence we obtain (6.7). â€¢ 6.1.2 B o u n d i n g the elements of Auu-The fut i l i ty of searching for diagonal approximations of Auu which eliminate the depen-dence on rij may be understood if we employ Theorem 2.4 of Demko et al . [22] which bounds the exponential decay rate of elements i n the inverse of a banded matr ix . Specifi-cally, we consider constant mesh spacing and constant diffusivity with Dirichlet boundary conditions for conservative diffusion such that the diagonal blocks of Auu are given by {Auu)33 = D ( 7 ^ 2 ) t r i [2 - r 2 , 8 + 2 r 2 , 2 - r 2 ] , j = 1 , . . . , M (6.8) and of dimension (L-l)x(L-l). A bound on the elemental decay rate may be written i n the form | { ( ^ ) - / } w | < C A ( r ) l ^ l (6.9) Chapter 6. Reduced Systems and Approximate Schur Complements 137 where the decay rate is itself bounded, A ( r ) < A 0 0 ( r ) = ^ ~ Q ; _ , a Â± = y/(4 + r 2 ) Â± |2 - r 2 | (6.10) We note that the bound is consistent with Theorem 6.1, clearly capturing its depen-dence on r . In particular, the exact nature of the approximation for r = y/2 is observed wi th A(\/2) = 0. But perhaps most importantly, Figure 6.1 reveals that A ^ r ) â€”> 1 as r â€”y oo, imply ing that the diagonal blocks of A~* are not only dense but al l the elements are of approximately the same magnitude. Thus, it is not surprising that a diagonal approximation is inadequate. Figure 6.1: Upper bound of the condition number, K (6.6) and of the decay rate Aoo(r). 6.2 A T w o - S t e p L u m p e d A p p r o x i m a t i o n El iminat ing u to obtain the Schur complement for v is an arbitrary choice and it may not always be the natural one. Given this preconditioner's dependence on the ordering of unknowns we postulate that it can be improved, especially for nonuniform grids, i n much the same manner that symmetric S O R (SSOR) improves upon SOR,( i .e . by composing the two possible orderings). To this end consider two splittings of S^, Sfi = Ei + Fi, Su, = E2 + F2 Chapter 6. Reduced Systems and Approximate Schur Complements 138 which are employed i n the following two stage iteration of the system S)Xph = Qg^-, E2K (fc+i) -Flfi[k) + QSfl - F 2 ^ ] + Qs^ (6.11a) (6.11b) Defining the residual i n the standard manner, = Q 5 â€” S^p^, allows (6.11a) to be rewritten i n the form Â»[k+Â» = ^ + E ^ (6.12) such that substitution of (6.12) into (6.11b) gives E ^ { k + 1 ) = - F 2 (â€ž[k) + E?rM) + QSli = E 2 ^ + (Ex - F2) ErV*> Hence the two-stage iteration is readily expressed as a single residual equation, {Ex (E, - F,)-1 E2} AÂ»lk+1) = r<*> (6.13) To apply this to the nodal equations we need only realize that a splitt ing exists for the lumped preconditioner. The splitt ing is described by Â£ 1 = A A A A Fi = 0 0 with Auu = Auu + A*uu. Similarly we introduce the lumping i n v writ ing Avv = Avv + A * Substitution into (6.13) yields a two step lumped preconditioner, A A A A -1 ' A A A A A 1 C A 'A" (6.14) Observe that S^y ^ <SWâ€ž and hence S ^ v w i l l be asymmetric i n general. This is unfortu-nate, as it precludes the use of preconditioned conjugate gradients used elsewhere i n this Chapter 6. Reduced Systems and Approximate Schur Complements 139 work. However, there are many alternative, transform free K r y l o v subspace solvers avail-able. In our numerical tests we employ preconditioned Orthores [75] wi th a truncation length of only five and observe excellent performance. The following result again relates the eigenvalues of the preconditioned system to the approximation difference i n the Schur complement, and like Proposition 3.1 can be easily proved. Proposition 6.2. The eigenvalues of the preconditioned system arising from the appli-cation of Sfj,;UiV are given by, A = I, of multiplicity Nv A = A - 1 1 - A AUV(^AVy) Ayy Sy AVU(^AUU) A â€” 1 , - 1 r an where Su â€” Auu Auv I Avv I A. - l Unfortunately, superelement analysis is not directly applicable to this preconditioner. Moreover, an easily obtainable bound such as A r < A r (<?u) SU - 1 A, is too crude to be of use. However, it is very encouraging that the eigenvalues may be expressed as a perturbation from unity, wi th the perturbation a product involving the errors associated with the lumped approximations. These errors possess a cell-based block diagonal structure hence for every cell at least one of A*u and A*vv is associated wi th a good approximation. Thus we anticipate S^UFV to be robust for high aspect ratio problems and we demonstrate this computationally in Section 6.5. Chapter 6. Reduced Systems and Approximate Schur Complements 140 6.3 Lumping the Cell-Based Schur Complement El iminat ing ph from (6.1) we obtain the Schur complement for the cell-based unknowns S*<PK = Qs, (6-15) where S$ = B(A-1 - A~1CTSc1CA-1)BT Qs, = Qt-BiA-i-A-WS^CA-^Qj The Schur complement S$ is symmetric positive definite and hence we consider the de-velopment of a preconditioner for it . Of particular importance to the robustness of the resulting solvers is that this reduced system (6.15) governs only the cell-based unknowns. However, the relative sparsity structure of the system has been compromised because Sc is block diagonal, wi th each block itself tridiagonal. Thus the matrix-vector mult ipl ica-t ion, S^d)h W 1 U involve the inversion of M , (L - 1) X (L - 1) and L, (M - 1) x ( M - 1) tridiagonal matrices. To develop a preconditioner for S^, and in particular one which is suitable for mul t igr id inversion, we return to the ful l system (5.40) and note that is positive definite if A is positive definite and range(i?) f ) range(C) = {0}. Hence, the primary objective is the minimizat ion of fill through symmetric positive definite approximations of A. To this end we continue the idea of the previous section and replace A , a block diagonal matr ix , wi th A a strictly diagonal lumped approximation. For N E M (and N I M w i t h S a = 0) lumping (5.41) we obtain, " 1 0 0 0 AxiAyi 0 1 0 0 0 0 1 0 . 0 0 0 1 (6.16) Chapter 6. Reduced Systems and Approximate Schur Complements 141 This simple approximation transforms S$ â€”> S$ where S$ is in fact the well known 5-point cell-centred discretization. In the nonconservative case lumping A j j of the N I M discretization (5.42) gives, Aij â€” r 0 0 0 AxiAyi 0 8$ hJ 0 0 0 0 0 0 0 0 (6.17) and hence we st i l l obtain a symmetric 5-point cell-centred discretization that is an M -matr ix , although it has nonstandard stencil weights. Thus we have derived an approxi-mate Schur complement that may be used as a preconditioner and which may again be efficiently inverted using a standard mult igr id method. 6.3.1 A Condition Bound The use of tridiagonal solves in the matrix-vector product S,pd>h makes a cell-based su-perelement analysis impossible. Certainly using superelements defined over entire rows and columns of the mesh is possible but obtaining analytic expressions for the corre-sponding eigenvalues seems improbable. It is, in fact, more promising to consider the ful l scalar system (i.e. the lumping of A in (5.40) transforms S^ifl) of (6.1) to S^^) for which we bound the condition number of the correspondingly preconditioned system by a constant independent of the mesh, the diffusion coefficient and the absorption cross-section. We begin wi th the conservative case, which defines a bound that w i l l be satisfied i n the nonconservative case. Theorem 6.2. For conservative diffusion the preconditioner <5(0,M) yields an effective condition number which is bounded by, - l >(<M < 3 (6.18) independently of the the mesh and the diffusion coefficient. Chapter 6. Reduced Systems and Approximate Schur Complements 142 Proof. Once again we employ superelement analysis to relate local bounds, which are readily obtained, to global ones. We begin by writ ing down the elemental problem, " 12(a + 7) â€” 6 a â€” 6 a â€” 6 7 â€” 6 7 â€” 6 a 4 a 2a 0 0 â€” 6 a 2 a 4 a 0 0 â€” 6 7 0 0 47 27 â€” 6 7 0 0 27 47 with a = DijAyj/Axi and 7 = DijAxi/Ayj and where the unknowns have been ordered, * s = [<l>i,3,^i-y,(l>i+y,^i,j-i,<l>i,j+i\T- Approximating A by A gives, " 4 (a + 7) - 2 a - 2 a - 2 7 - 2 7 - 2 a 2a 0 0 0 - 2 a 0 2 a 0 0 - 2 7 0 0 27 0 - 2 7 0 0 0 27 where it is apparent that we have preserved the kernel, iC = M5(^=M^)) = [1.1'1.1.1]T Solving the corresponding generalized eigenvalue problem, S ( S < M * s = kS^Qs W i t h * S Â± / C we obtain, A s = [1,1,3,3]. Since these eigenvalues are independent of the cell indices the global bound follows immediately. â€¢ Remark 6.3. The bound (6.18) is unchanged if we extend the analysis to a diagonal diffusion tensor. The implicat ion is that this preconditioner is robust w i t h respect to anisotropy as well as spatial variations in the diffusion tensor. Corollary 6.3. . For the lowest-order NEM discretization with nonconservative diffu-sion [Definition 3-4-4] l^e condition number of the preconditioned system satisfies the bound given in (6.18) of Theorem 6.2. Chapter 6. Reduced Systems and Approximate Schur Complements 143 Proof. In the case of N E M the addition of absorption requires only a t r iv ia l change in the proof of Theorem 6.2. Specifically the (1,1) elements of the superelement matrices are perturbed as follows, 12 (a + 7) - Â» 12 (a + 7) + e 4 ( a + 7) 4 ( a ' + 7) + e where e = (S a ) j . . A x t - A j / j . This perturbation makes the corresponding pure Neumann problem well posed and thus eliminates the kernel. The eigenvalues are easily evaluated with A s = [1,1,1,3,3] and hence we obtain the bound (6.18). â€¢ Corollary 6.4. For the constant-constant NIM discretization with nonconservative dif-fusion the condition number of the preconditioned system satisfies the bound given in (6.18) of Theorem 6.2. Proof. Following the proof of Theorem 6.2 we employ superelement analysis to obtain a global bound on the condition number. We begin by writ ing down the elemental problem, ' 4 ( a + 7) + e - 2 a - 2 a - 2 7 - 2 7 - 2 a (Â« + /?) ( a - / ? ) 0 0 - 2 a ( a - / 3 ) (a + 0) 0 0 - 2 7 0 0 (7 + 0 ( 7 - 0 - 2 7 0 0 ( 7 - 0 (7 + 0 with a = Dijrijafj, 7 = A j r ^ a - j , = D^r^Br) and c = DitJrhJj3\VJ. Approximat ing A by A (6.17) gives, 4(a + 7) -2/? -2(3 - 2 c - 2 c â€¢ -2(3 2/? 0 0 0 -2(3 0 2/? 0 0 - 2 c 0 0 2c 0 - 2 c 0 0 0 2c Solving the corresponding generalized eigenvalue problem, we obtain, Chapter 6. Reduced Systems and Approximate Schur Complements 144 Consider the eigenvalue, A W = 2J? = AJ-> t anh(Aff) ' (ft [Ag>coÂ«h(Ag>)-l] which is a symmetric function of \\j . It is possible to show that has a horizontal asymptote at A s = 1 and a local max imum of \{x\0) = 3. M oreover, it is monotonically decreasing for \\Xj > 0 implying that it satisfies the bound A ^ Â£ (1,3). The other non-unity eigenvalue may be bounded i n the same manner and thus we find that the condition number is bounded by (6.18) â€¢ Based on the equivalence of mixed and mixed-hybrid F E M s it is apparent that S$ is equivalently the reduced system obtained i n mixed F E M s . Generally, it is problematic to solve this system because of the loss of sparsity. Here we are able to take advantage of the restricted geometry, incurring only the additional cost of tridiagonal solves. A sensitivity to highly variable coefficients has also been a problem for some methods. A l l a n et al . [2] observed this and employed a lumped diagonal preconditioner to the mixed formulation, and in taking advantage of the assumed rectangular geometry obtained results comparable to ours. However, wi th the use of the mixed-hybrid discretization we have found a simpler approach to a robust cell-based preconditioner which i n the form t>(0lAt) has the potential to be extended to more general geometries. 6.4 Approximate Multi-level Inversion of the Preconditioner Using either of the aforementioned preconditioners would be clearly impract ical if we actually intended to "invert" them to some strict level of accuracy. However, this is unnecessary because the preconditioners are standard discretizations of (2.12) and hence standard mult igr id methods yield an error reduction per iteration which is independent of the grid size. This implies a spectral equivalence which allows us to employ single Chapter 6. Reduced Systems and Approximate Schur Complements 145 V - or W-cycle in lieu of a complete mult igrid solve. We state and prove the following L e m m a for completeness. L e m m a 6.1. Let A,B,C be n x n symmetric positive definite matrices, and assume p(I -C~1B)<c<l (6.19) where p denotes the spectral radius. Then Â« ( C " U ) < { ^ j ^ ( t f - 1 - 4 ) ( 6 - 2 Â° ) Proof. Let A be an eigenvalue of C~1B. Then (6.19) implies 1 â€” c < X < 1 + c and hence we readily obtain (6.20). â€¢ The practical implicat ion of Theorem 6.1 follows if we let A be the edge-based or the cell-based <S^, and let B be its symmetric preconditioner. B is further approximated by C, a single mult igr id cycle for which we have c independent of grid size, typical ly c ~ 0.1. Thus, < 1.3 implying that C is (almost) as effective as B. 6.5 N u m e r i c a l E x p e r i m e n t s A progressive test suite is presented which systematically highlights the strengths and weaknesses of each approximate Schur complement preconditioner. Beginning with a constant coefficient Dirichlet problem on the unit square [Section 6.5.1] we verify numer-ically that on a uniform mesh (r = 1) all three preconditioners exhibit mesh independent convergence. Continuing with this simple example we vary the aspect ratio of a constant mesh to confirm the breakdown of the simple edge based preconditioner S^-u and the robustness of and S<f,. In the next test we consider a groundwater flow problem [Section 6.5.2] with significant jumps in the diffusion coefficient. To isolate the influ-ence of this spatial dependence we solve this problem on a uniform grid, once again Chapter 6. Reduced Systems and Approximate Schur Complements 146 observing mesh independent convergence. To evaluate the relative cost of each precondi-tioner machine timings are also presented. Final ly , in Section 6.5.3 we present a diffusive checkerboard problem which combines significant jumps in D(x,y) wi th a spatially de-pendent grid containing high aspect ratio cells. The inevitable breakdown of <SWâ€ž is observed along with the robustness of S^u^v and S^. Note that the following relative residual convergence criteria, is employed in al l of the numerical tests that appear in this section. 6.5.1 A Toy Problem We begin with the Dirichlet problem of Section 3.7.1 in the conservative case (i.e. S a = 0). A comparison of the two edge-based preconditioners, and the cell-based preconditioner is presented i n Table 6.1 for the domain [0,1] x [0,1] with a = 3 = 2. It is clear from these results that all preconditioners generate an average residual reduction which is independent of the mesh size. Moreover, there is no difference in the iteration counts for those runs which inverted the preconditioner exactly (i.e. columns marked C) and those which only used a single V ( l , l , l ) cycle (i.e. columns marked V ) . Thus, these solvers offer an efficiency for the solution of the nodal equations which is comparable to standard mult igr id algorithms applied to standard discretizations of Poisson equations. Moreover, these results are consistent with the theoretical analysis performed in Chapter 6. Unfortunately, we must also demonstrate the vulnerability of 5 M ; u regarding high aspect ratio cells and hence we conducted several runs in which the size of the physical domain ([0, a] x [0, b]) was varied while the mesh size remained fixed. The results for the preconditioned edge-based solver S^u are summarized in Table 6.2. Here it is apparent that for r sufficiently close to 1, this solver demonstrates excellent efficiency. Performance Chapter 6. Reduced Systems and Approximate Schur Complements 147 is st i l l excellent as r approaches zero. However, just as the bound on K [Theorem 6.1, equation (6.6)] predicts, the efficiency is degraded as r increases. Conversely, the two step lumped preconditioner S^-U,v displays perfectly symmetric iteration counts on the 40 x 40 mesh [Table 6.3]. This symmetry in r is very encouraging, particularly i n light the preconditioner's asymmetry. Similarly, Table 6.4 clearly displays the r-independent convergence of S^,. Table 6.1: Iteration counts for the preconditioners, S^U and S$ for the toy problem wi th constant mesh spacing on the domain [0,1] X [0,1]. Mesh Size c V c V c V 20 x 20 6 6 4 4 11 11 40 x 40 6 6 4 4 11 11 80 x 80 6 6 4 4 11 11 Table 6.2: Iteration counts for the lumped preconditioner, S^;U wi th constant mesh spac-ing on the domain [0,a] x [0,6]. {Note: a = 2/a,fl = 2/6} Mesh Size Aspect ratio: r = b/a 1/8 1/4 1/2 1 2 4 8 10 x 10 9 9 9 6 7 15 19 20 x 20 9 10 9 6 7 16 30 40 x 40 9 10 9 6 7 15 29 Table 6.3: Iteration counts for the two-step lumped preconditioner, <SM;U)â€ž wi th constant mesh spacing on the domain [0,a] x [0,6]. {Note: a = 2/a,(3 = 2/6} Mesh Size Aspect ratio: r = b/a 1/8 1/4 1/2 1 2 4 8 10 x 10 7 7 5 4 5 6 7 20 x 20 9 7 5 4 5 6 8 40 x 40 9 7 5 4 5 7 9 Chapter 6. Reduced Systems and Approximate Schur Complements 148 Table 6.4: Iteration counts for the cell-based preconditioner, S$ wi th constant mesh spacing on the domain [0,a] x [0,6]. {Note: a = 2/a,3 = 2/6} Mesh Size Aspect ratio: r â€” b/a 1/8 1/4 1/2 1 2 4 8 . 10 x 10 9 10 10 10 10 10 9 20 x 20 10 10 11 11 11 . 11 10 40 x 40 11 11 11 11 11 11 11 6.5.2 G r o u n d W a t e r F l o w The saturated flow problem of Mose et al . [56] serves both as an excellent test of the preconditioners developed in Section 6 and as a showcase for the discretizations them-selves. The problem, shown in Figure 6.2 models the flow of a ground water through the channel formed by the impervious sides, J â€¢ n = 0 V O , y) Â£ T 7 = {x = 0, x = 100,0 < y < 100} and driven by the externally imposed gradient, <b(x, y) = 100 V(x , y) Â£ TT = {0 < x < 100, y = 100} (6.21) <f>(x, y) = 99 V ( x , y) Â£ TB = {0 < x < 100, y = 0} (6.22) The flow is impeded by two regions of extremely low hydraulic conductivity. Solutions computed on a very coarse mesh of 25 x 25 are shown i n Figure 6.3. Part icularly i m -pressive are the streamlines which are significantly more accurate than those obtained even wi th careful postprocessing of conforming methods [56]. Computations were per-formed on a uniform mesh with the resulting iteration counts presented i n Table 6.5. The columns marked CG(Sa) and CG(S$) record the performance of conjugate gradients ap-plied directly to the reduced systems and S<p, respectively. This provides an indication of the problem's conditioning. For the purposes of comparison iteration counts for diag-onal preconditioning of the edge-based system, denoted are also included. Consistent Chapter 6. Reduced Systems and Approximate Schur Complements 149 wi th the theoretical bounds, the entries for the three approximate Schur complement preconditioners are identical to the constant coefficient Dirichlet problem [Table 6.1]. In Table 6.6 we present machine timings for the three preconditioners 1 . Not surpris-ingly the simplest preconditioner, S^>u which in this case also yields the smallest bound on the condition of the system, is the fastest. However, we have established the lack of robustness i n S^u and it is encouraging to see that only approximately 50% more t ime is required by the robust preconditioners S^UfV and S^. It is also interesting to note that despite having approximately twice as many unknowns in the iterative process and needing OrthoRes in place of a simple conjugate gradient solver, S^UtV is competitive wi th S<f,. This is of particular interest to people with existing codes designed to solve the edge-based system. Final ly, we note that on the 200 x 200 problem the approximate Schur complement preconditioners are approximately 16 times faster than the simple diagonal preconditioner. 20 60 20 100 Dirichlet: <|)=100 Impervious â€” (Neumann: J-n=0) o Hydraulic Conductivity, D=l 4^ o 0 Dirichlet: <j)=99 0 100 Figure 6.2: A schematic of Mose et al.'s [56] first example. t i m i n g s were per formed on an H P A p o l l o 9000/735 (99 M H z P A - R I S C 7100) Chapter 6. Reduced Systems and Approximate Schur Complements 150 Figure 6.3: Solution <b(x,y) and streamlines of Mose et al.'s [56] first example (above) Table 6.5: Iteration Counts with single V ( l , l , l ) cycles Mesh Size CG(Sâ€ž) Â°Â» S,f, CG(S$) 50 x 50 795 342 6 4 11 667 100 x 100 1,675 600 6 4 11 1,709 200 x 200 3,049 1,027 6 4 11 4,486 Table 6.6: Iteration Timings, [s] for the flow problem with single V ( 1,1,1) cycles Mesh Size CG(Sâ€ž) s$ CG{St) 50 x 50 7.62 3.50 0.430 0.610 0.450 5.95 100 x 100 72.3 32.4 2.13 3.26 2.51 71.5 200 x 200 601 243 10.4 15.53 15.8 939 Chapter 6. Reduced Systems and Approximate Schur Complements 151 6.5.3 A Diffusive Checkerboard The final test of the preconditioners combines both huge jumps in the diffusion coefficient wi th a spatially dependent grid containing high aspect ratio cells. Shown schematically i n Figure 6.4 we consider a single checkerboard on the square domain [0, 24] x [0,24] wi th fix = {(0,12) x (0 ,12)} U { (12 ,24) x ( 1 2 , 2 4 ) } , 0 2 = {(0,12) x (0 ,12)} U { ( 1 2 , 2 4 ) ' x along TR = {x = 0,0 < y < 24} U {0 < x < 24,y = 0}. Vacuum boundary conditions are specified along Tv = {x = 24,0 < y < 24} U {0 < x < 24, y = 24}. In the constant coefficient case wi th a constant source this imposes as natural flow of neutrons along streamlines y â€” xâ€” constant. Introducing a piecewise constant definition of the diffusion coefficient and the source creates a significant perturbation of this idealized flow. Particularly, the relatively low diffusion coefficient in fl2 serves to channel the neutrons through a single point, (12,12) creating the large gradients observed i n Figure 6.5. Thus, this problem introduces a mixed boundary condition and a significantly worse singularity than was present i n the ground water flow problem Section 6.5.2. To further test the preconditioners and to accommodate internal layers near the interfaces of f i i and f i 2 we introduce an exponential (12,24)}. We begin by introducing reflective boundary conditions (2.19b), J n = 0 v ( x , y ) e r H (2.18), Chapter 6. Reduced Systems and Approximate Schur Complements 152 (Neumann: J>n=0) Figure 6.4: A diffusive checkerboard configuration based on Dendy's [23] example. mesh. Specifically on the region 0 < x < 12 we consider a grid spacing which is uniform in area such that 2 r12 I Ax = - e^xdx, x i + 1 = â€” log(e^x + l x Ax) L Jo Ix similar definitions are employed for 12 < x < 24 and for the grid i n y. A typical mesh generated with ^x = jy = 2/10 is displayed in Figure 6.6. The iteration counts for this grid are given in Table 6.7 where column r indicates the m a x i m u m aspect ratio present. We first note that the iteration counts for CG{Sll) and CG(S(f>) are comparable and extremely large, (16,055 and 19, 054 respectively on the 96 x 96 mesh) giving a clear indication of the conditioning of the problem. The influence of SD is quite impressive given its simplicity, achieving iteration counts which are only 25% higher than for the ground water flow problem. A s expected the iteration counts for <SWâ€ž are extremely poor, remaining consistent with the bound given i n Theorem 6.1. Conversely the results for both Â«? W t t ) t i and S<p remain excellent. We also note that a comparison of the iteration counts for S^-u^ in Table 6.7 wi th those for the constant coefficient Dirichlet problem in Table 6.3 with r = 8 (or r = 1/8) reveals a comparable Chapter 6. Reduced Systems and Approximate Schur Complements 153 X Figure 6.5: The solution, <f>(x,y) of Figure 6.6: 2D Exponential G r i d with the diffusive checkerboard computed the parameters 7^ = -yy = 2/10. on a 24 x 24 uniform mesh performance on this significantly more difficult problem. The equivalent comparison for Sj, leads to the same conclusion. Table 6.7: Iteration counts for jx â€” jy â€” 2/10. Mesh Size r CG(Sâ€ž) S$ CG{S+) 24 x 24 7.71 2,129 181 59 9 13 2,063 48 x 48 9.04 6,190 374 74 10 13 6,932 96 x 96 9.91 16,055 759 86 11 14 19,054 Chapter 7 Conclusions 7.1 Research Summary In this research we have addressed the fundamental issues that have restricted the prac-t ical use of nodal methods to the reactor physics community. Specifically, we have en-hanced the mathematical foundation of the lowest-order nodal methods by establishing new equivalences wi th certain F E M s and by providing a generally applicable truncation error analysis. In addition, we have developed a new family of preconditioners for the iterative solution of the resulting system of equations, which yield optimally fast solvers. F inal ly , we have introduced a new efficient multi-level homogenization procedure. Beginning in Section 3.6 we introduced a rudimentary truncation analysis that may be used to verify the consistency of nodal discretizations. The incompleteness of estab-lished equivalences between all nodal methods and corresponding F E M discretizations (e.g. ut i l izat ion of higher order transverse leakage approximations [Section 3.4.1]) makes this approach necessary. We demonstrated the ut i l i ty of this technique wi th the analysis of the N E M [Definition 3.4.2] and N I M [Definition 3.5.1] discretizations. Specifically, under the assumption of piecewise constant coefficients we verified the second-order con-sistency of these popular nodal methods. This analysis also served to highlight the exact treatment of both the cell balance relation, and the continuity (in the weak sense) of normal currents, which is instrumental i n their performance. In addition, the trunca-tion errors for the corresponding one-sided normal currents were used to characterize 154 Chapter 7. Conclusions 155 situations in which N I M is anticipated to be more accurate than N E M . This was demon-strated computationally in Section 3.7. Unfortunately, i n the general setting of piecewise smooth coefficients [i.e. Z},-,(r), ( E a ) . .(r)] this analysis revealed that the one-sided nor-mal currents of both N I M and N E M are only first order consistent, wi th the leading term depending on the derivative of the diffusion coefficient. Moreover, we determined that the balance equation is no longer exact, exhibiting a second-order error that depends on the first derivatives of the absorption coefficient. These properties are an asset i n applications which rely on homogenization procedures. Nevertheless, a discussion of the key issues surrounding the restoration of second-order consistency and exact cell balance was presented. This rudimentary analysis considered the nodal discretization in isolation, while ig-noring the influence of the homogenization procedure. However, the approximations of Equivalence Theory [Section 4.1.3], which are the most popular homogenization tech-niques employed i n reactor physics, introduce discontinuities i n the flux. Discontinuous F E M s [1] might offer a formal treatment of this scenario, but the jumps in the flux are determined through local auxiliary fine scale computations. This latter aspect of Equivalence Theory complicates the analysis and discourages its use i n more general settings. Moreover, the impressive accuracy that is obtained on assembly size meshes i n practical reactor computations reflects not only the aforementioned properties of the nodal discretization, but its integration with the approximations of Equivalence Theory. A n unfortunate consequence of this interdependence is that the discretization itself has received scant consideration outside the reactor physics community. Efforts to adopt a more rigorous treatment of homogenization, which may improve this situation, are con-sidered i n [76]. Unfortunately, even with an asymptotic treatment of the nearly periodic case, practical computations sti l l rely on the solution of local auxiliary problems. Chapter 7. Conclusions 156 In Section 4.3 we introduced a new mult igr id homogenization procedure which cir-cumvents many of these problems. In particular, the operator-induced coarsening of black box mult igr id is used to construct coarse-scale operators from a specified fine-scale dis-cretization. Since this coarsening procedure is a generalization of variational coarsening, a bilinear finite element discretization is used on the fine grid and assumed on the coarse grid. The diffusion tensor is then extracted from the coarse-grid operator according to Theorem 4.3. Thus we shift the assumption of a priori knowledge of the exact solution to a specification of the solution space. We believe that i n this respect our technique marks a considerable step forward i n practical homogenization. We demonstrated the excellent performance of our technique for a series of periodic test problems that included a computation involving nonzero off-diagonal terms. This abil ity is crucial for modelling layered structures that are not aligned with the coordinate axes. Nevertheless, there are many applications for which homogenization procedures al-ready exist, and frequently these do not address the possible need of off-diagonal entries i n the homogenized diffusion tensor. In this vein a few researchers have focused on the need for a stronger mathematical foundation to assess the applicabil ity of nodal dis-cretizations for general diffusion problems. This research established the equivalence of various nodal methods with certain nonconforming and mixed-hybrid F E M s [36, 39]. Our investigation of these equivalences in the lowest-order case produced several new results. In particular, we found that the nodal integration method ( N I M ) was equivalent to a pr imal nonconforming method, which employed an analytic basis, if exact integration is used to evaluate the elements of the stiffness and mass matrices [Theorem 5.2]. More-over, ut i l iz ing a quadratic polynomial basis in a pr imal nonconforming F E M is equivalent to the lowest-order nodal expansion method ( N E M ) if a special treatment of the mass Chapter 7. Conclusions 157 matr ix is employed. 1 We demonstrated that a simple mass lumping is sufficient, and cor-responds to a restoration of a cell balance relation and the continuity of normal currents. Recall ing that the pr imal nonconformal methods do not rigorously enforce cell balance, an apparent advantage of the analytic basis is that it intrinsically embodies this prop-erty. However, a rigorous treatment of cell balance does appear i n mixed-hybrid F E M s . In fact, Hennart and del Valle [39] have established the equivalence of the lowest-order N E M wi th the Raviart-Thomas mixed-hybrid methods [Section 5.3.2]. The one-dimensional analysis of Section 5.2 was instrumental in identifying the subtle aspects of the nonconformal equivalences. It has also provided important information regarding the one-dimensional properties of these discretizations that do not generalize completely to two dimensions. Specifically, in the one-dimensional case el imination of the cell unknowns yields an edge-based system which, i n certain cases, is an M - m a t r i x . In particular, wi th the analytic basis functions (i.e. corresponding to N I M ) we found that the edge-based system is unconditionally an M - m a t r i x . Similarly, the polynomial basis unconditionally yields an M-matr ix when the exact stiffness and mass matrices are used. However, once lumping is employed to establish a cell balance relation and continuity of the current (i.e. corresponding to N E M ) , we obtain an M - m a t r i x subject to a constraint on the mesh size that depends on ^ / ( E a ) ^ . / Z ) , ' j - Thus, the analytic basis yields a discretization that appears to capture the essential physical and mathematical properties of the differential problem while the polynomial basis does not. Moreover, the implicat ion of the edge-based system corresponding to an M - m a t r i x is that the global cell-and edge-based system satisfies a A~L > 0. We conjecture that although an M - m a t r i x does not arise i n two dimensions, the global cell- and edge-based system also satisfies A-1 > 0 in a manner consistent with the one-dimensional case. If proven, this would x Note that in the conservative case the analytic and polynomial bases (i.e. N I M and N E M ) are equivalent and hence quite fortuitously the polynomial basis inherits its advantageous properties. Chapter 7. Conclusions 158 establish an important advantage for N I M : i n difficult problems wi th severe variations i n the coefficients no restriction of mesh spacing would be necessary. Despite the many desirable properties that the families of nodal methods and mixed-hybr id finite element methods embody, their very design makes the solution of the re-sulting equations awkward and costly, thwarting many potential users. Yet , inherent i n this structure is a natural partitioning of the system which can be uti l ized i n the develop-ment of preconditioners. The existence of such a partit ioning suggests the investigation of various reduced systems (e.g. Schur complements) in conjunction wi th suitably sparse approximations. In particular, we presented three such approximations, <SM;â€ž, S^u,v and S4, and demonstrated that the preconditioning which resulted was optimal i n the sense that a fixed number of iterations, independent of the mesh size, was required to reduce the residual by a fixed amount. These solvers are competitive with multi- level methods because the preconditioner is only approximately inverted wi th a single V - or W - m u l t i g r i d cycle. Unfortunately, S^u which preconditions the most popular form of the discretization, is sensitive to high aspect ratio cells having r > 1. This shortcoming, an artifact of arbitrari ly approximating Auu as opposed to Avv, was alleviated by the two step pre-conditioner Sfj.;UtV. Conversely, an identical approximation is made in both coordinate directions during the development of S$ and hence its preconditioning is insensitive to high aspect ratio cells. Moreover, the bound obtained in Theorem 6.2 suggests that in a more general setting the system S^tfl) w i l l be of significant practical interest. In all cases the preconditioners were found to be robust with respect to spatial variations i n the diffusion coefficient. Chapter 7. Conclusions 159 7.2 Future Work Our analysis of nodal discretizations considered their most common treatment in reactor physics, assuming a tensor product mesh with piecewise constant homogenized coeffi-cients. In fact, this restricted setting represents a significant class of modelling problems. Nevertheless, these assumptions may be inadequate in more general settings. In Sec-tion 3.6 we discussed the issues surrounding the use of piecewise smooth coefficients and concluded that this is a feasible extension to pursue. In addition, the asymptotic ho-mogenization techniques indicated that a piecewise smooth diffusion tensor is necessary to accurately model general problems on a coarse mesh. Unfortunately, neither N E M nor N I M facilitates a straightforward extension in this both the tangential and normal components of the current are required at the cell interface. A similar situation arises wi th a logically rectangular mesh, particularly as the application of the transverse integration procedure becomes unclear. We feel that i n each of these cases examina-tion of the corresponding mixed-hybrid F E M offers an excellent starting point for the investigation. The incomplete analysis found with the lowest-order nodal discretizations is com-pounded for higher-order methods. M a n y methods are referred to in the literature as higher-order because a particular aspect of the discretization utilizes a higher-order ap-proximation, yet a correspondingly higher order truncation error is not necessarily at-tained [36]. Nevertheless, the construction of higher-order nodal methods that adhere to the principles outlined in this research, is a compelling topic for future consideration. The results of the new mult igrid homogenization technique presented i n Section 4.3 were restricted to the periodic case. The extension of the local result given in Theorem 4.3 to the nonperiodic case is straightforward. Specifically, no change is required for Neumann boundary conditions, as these are the natural B C s of the pr imal variational Chapter 7. Conclusions 160 formulation, while Dirichlet B C s may be incorporated directly. We note that if the stencil weights vary spatially on the coarse mesh this local result is a second order approxima-tion of the mult igr id homogenized diffusion tensor. Al though this may be unavoidable in the periodic case, for Neumann and Dirichlet B C s it seems possible to extract the exact mult igr id homogenized diffusion tensor through a series of one-dimensional least-squares problems. However, the exact mult igr id homogenized diffusion tensor is itself an approximation that reflects the operator induced coarsening procedure. A n analysis of this procedure is required to understand the errors that arise i n truly multidimensional problems, particularly in pathological cases such as the checkerboard. In addition to a theoretical analysis of the errors, it w i l l be very informative to apply the mult igr id homogenization technique to realistic modelling problems. The abil ity of the nodal integration method to bui ld analytic information into the discretization suggests that it may be a very powerful technique for discretizing other differential equations. Research in this area, while st i l l focused on aspects reactor simula-t ion, include applications of N I M to the transport equation [52] and problems i n f luid flow [28, 43, 7]. In a preliminary investigation we have applied N I M to the convection-diffusion equation (see also [55, 28]) wi th encouraging results. Specifically, the use of an analytic one-dimensional solution to the transverse inte-grated O D E s provides an intrinsic adaptive upwinding. Moreover, it may be viewed as a two-dimensional generalization of the well-known Scharfetter-Gummel discretization em-ployed in one-dimensional semiconductor modelling [63]. Our results for the circulating flow of Smith and Hut ton [67], a popular convection diffusion benchmark, are displayed i n Figure 7.1. Note the impressive outflow profile, which is exact i n the convective l imi t . However, there are also small oscillations observed near the layer [Figure 7.2]. Through further study of the convective l imit we identified the primary source of this instabil i ty: Chapter 7. Conclusions 161 the violation of the assumption that the transverse normal currents are approximately constant. The prospect of employing assumptions which are compatible with both the hy-perbolic and elliptic limits of the convection-diffusion equation represents an interesting approach for new discretizations. y -1 inflow 0 outflow X Figure 7.1: Stream lines for the convection dominated flow of Smith and Hutton [67]. (a) (b) Figure 7.2: The convection dominated flow of [67] with Pe = [D]"1 = 500 computed using a constant-constant N I M : (a) cell averages (d) inflow and outflow profiles. Chapter 7. Conclusions 162 The solvers that we considered i n Chapter 6 were based on preconditioning K r y l o v subspace methods. However, it is possible to view these preconditioners as coarse-grid op-erators for which the specification of intergrid transfer operators and a fine-grid smoother would yield a multi- level algorithm. Our investigation of this scenario for the cell-based coarse-grid operator S$ is very encouraging. Ul t imately this approach may play a sig-nificant role in the iterative solution of higher-order methods, since in this setting the lowest-order nodal discretizations may be viewed as coarse-grid operators themselves. Changes to the underlying assumptions of the discretization present new challenges for an iterative solver. It was noted i n Chapter 6 that the extension to a piecewise constant diagonal tensor was straightforward and that both <SM;Uit, and <S^ are expected to perform well . The treatment of piecewise smooth coefficients is unclear as it may necessitate the use of a higher-order discretization, i n which case we may appeal to the multi- level approach described above. In addition, both the extension to a ful l tensor and logically rectangular meshes pose a formidable challenge. Specifically, if the ult imate form of a nodal discretization i n these cases is assumed to be that of the corresponding mixed-hybrid F E M , it is apparent that A remains block diagonal, but the blocks are dense. Nevertheless we believe an extension of our technique is possible in these cases. B i b l i o g r a p h y [1] M . L . ADAMS AND W . R . MARTIN, Diffusion synthetic acceleration of discontin-uous finite element transport iterations, N u c l . Sci. Eng . , I l l (1992), pp. 145-167. [2] M . B . A L L E N , R . E . E w i N G , A N D P . L u , Well-conditioned iterative schemes for mixed finite-element models of porous-media flows, S L A M J . Sci. Stat. Comput . , 13 (1992), pp. 794-814. [3] T . ARBOGAST A N D Z . C H E N , On the implementation of mixed methods as non-conforming methods for second-order elliptic problems, M a t h . Comput . , 64 (1995), pp. 943-972. [4] A R G O N N E NATIONAL LABORATORY, Benchmark Problem Book: ANL-7416, Sup-plement 2, 1977. [5] D . N . A R N O L D AND F . BREZZI, Mixed and nonconforming finite element methods: Implementation, postprocessing and error estimates, M a t h . M o d e l . Numer. A n a l . , 19 (1985), pp. 7-32. [6] Y . Y . A Z M Y AND J . J . DORNING, A nodal integral approach to the numerical solu-tion of partial differential equations, in Advances in Reactor Computations, vol . II, Salt Lake Ci ty , U T , M a r c h 28-30, 1983, American Nuclear Society, A N S , pp. 893-909. [7] , Nodal integral method solutions of cavity thermal convection problems, in N u -merical Methods i n Laminar and Turbulent F low, IV, C . Taylor, M . D . Olson, P. M . 163 Bibliography 164 Gresho, and W . G . Habaski, eds., Swansea, Wales, U K , July 9-12 1985, Pineridge Press, U K , pp. 753-767. [8] Y . Y . A Z M Y A N D B . L . K l R K , Performance modeling of parallel algorithms for solving neutron diffusion problems, N u c l . Sci. Eng. , 120 (1995), pp. 1-17. [9] G . I. B E L L AND S. GLASSTONE, Nuclear Reactor Theory, Van Nostrand Reinhold Company, New York, 1970. [10] F . BENNEWITZ, H . F l N N E M A N N , AND M . R . W A G N E R , Higher order corrections in nodal reactor calculations, Trans. A m . N u c l . S o c , 22 (1975), pp. 250-251. [11] A . B E N S O U S S A N , J . - L . L I O N S , A N D G . P A P A N I C O L A O U , Asymptotic Analysis For Periodic Structures, vol . 5 of Studies in Mathematics and its Applicat ions, N o r t h -Holland, New York, 1978. [12] J . F . BOURGAT, Numerical experiments of the homogenization method for operators with periodic coefficients, in Computing Methods in A p p l i e d Science and Engineer-ing, I, R. Glowinski and J . -L . Lions, eds., Versailles, December 5-9 1977, Springer, pp. 330-356. [13] A . B R A N D T , Multi-level adaptive solutions to boundary-value problems, M a t h . Corn-put. , 31 (1977), pp. 333-390. [14] F . BREZZI A N D M . FORTIN, Mixed and Hybrid Finite Element Methods, no. 15 in Springer Series i n Computational Mathematics, Springer-Verlag, New York, 1991. [15] G . F . C A R E Y AND J . T . O D E N , Finite Elements: A Second Course, vol . II of The Texas F ini te Element Series, Prent ice -Hal l , Inc., Englewood Cliffs, N J , 1983. Bibliography 165 [16] K . M . C A S E A N D P . F . Z W E I F E L , Linear Transport Theory, Addison-Wesley, Reading, Massachusetts, 1967. [17] B . C H A R , K . G E D D E S , G . G O N N E T , B . L E O N G , M . M O N A G A N , A N D S . W A T T , Maple V: Language Reference Manual, Spinger-Verlag, New York, 1991. [18] Z . C H E N A N D D . K W A K , The analysis of multigrid algorithms for nonconforming and mixed methods for second order elliptic problems, I M A Preprint 1277, Institute for Mathematics and its Applications, Minneapolis, M N , 1994. [19] C . C O R D E S AND W . K l N Z E L B A C H , Continuous groundwater velocity fields and path lines in linear, bilinear and trilinear finite elements, Water Resour. Res., 28 (1992), pp. 2903-2911. [20] B . DAVISON, Neutron Transport Theory, Oxford University Press, London, 1957. [21] D . L . D E L P , D . L . F I S C H E R , J . M . H A R R I M A N , A N D M . J . S T I D W E L L , FLARE: A three-dimensional boiling reactor water simulator, A t o m i c Energy Commission Research and Development Report G E A P - 4 5 9 8 , A t o m i c Power Equipment Depart-ment, General Electric , July 1964. [22] S. D E M K O , W . F . MOSS, AND P . W . SMITH, Decay rates for inverses of band matrices, M a t h . Comput . , 43 (1984), pp. 491-499. [23] J . E . D E N D Y , Black box multigrid, J . Comput . Phys. , 48 (1982), pp. 366-386. [24] J . E . D E N D Y , J . M . H Y M A N , AND J . D . M O U L T O N , The black box multigrid nu-merical homogenization algorithm, Tech. Rep. L A U R - 9 6 - 3 5 8 8 , Theoretical Divis ion , Los Alamos National Laboratory, Los Alamos, N M , M a r c h 1994. Bibliography 166 [25] I. D l L B E R A N D E . E . L E W I S , Variational nodal methods for neutron transport, N u c l . Sci. Eng . , 91 (1985), pp. 132-142. [26] J . J . D U D E R S T A D T A N D L . J . H A M I L T O N , Nuclear Reactor Analysis, John W i l e y k Sons, Inc., New York, 1976. [27] J . J . D U D E R S T A D T A N D W . R . M A R T I N , Transport Theory, John W i l e y & Sons, Inc., New York, 1979. [28] P . D . E S S E R A N D R . J . W l T T , An upwind nodal integral method for incompressible fluid flow, N u c l . Sci. Eng. , 114 (1993), pp. 20-35. [29] R . E . E W I N G , Y . K U Z N E T S O V , R . D . L A Z A R O V , A N D S. M A L I A S S O V , Finite Element Modeling of Environmental Problems, John Wi ley & Sons, New York , 1995, ch. Substructure Preconditioning for Porous Flow Problems, pp. 303-329. [30] C . F E D O N - M A G N A U D , J . P . H E N N A R T , A N D J . J . L A U T A R D , On the relation-ship between some nodal schemes and the finite element method in static diffusion calculations, i n Proceedings of A Topical Meeting on Advances in Reactor Compu-tations, vol . 2, Salt Lake Ci ty , Utah , March 28-30, 1983, Amer ican Nuclear Society, Mathematics and Computat ion Divis ion, pp. 987-1000. [31] H . F l N N E M A N N , F . B E N N E W I T Z , A N D M . R . W A G N E R , Interface current techniques for multidimensional reactor calculations, Atomkernenergie, 30 (1977), pp. 123-128. [32] H . D . F I S C H E R A N D H . F l N N E M A N N , The nodal integration method - A diverse solver for neutron diffusion problems, Atomkernenergie-Kerntechnik, 39 (1981), pp. 229-236. Bibliography 167 [33] G . G R E E N M A N , K . S M I T H , A N D A . F . H E N R Y , Recent advances in an analytic nodal method for static and transient reactor analysis, in Computational Methods in Nuclear Engineering, vol. I, Williamsburg, Virginia, Apri l 23-25, 1979, American Nuclear Society, A N S , pp. 3-49-3-72. [34] N . K . G U P T A , Nodal methods for three-dimensional simulators, Progress in Nuclear Energy, 7 (1981), pp. 127-149. [35] J . P . H E N N A R T , A general finite element framework for nodal methods, in The Mathematics of Finite Elements and Applications V , J . R. Whiteman, ed., London, 1985, Academic Press, pp. 309-316. [36] , A general family of nodal schemes, S I A M J . Sci. Stat. Comput. , 7 (1986), pp. 264-287. [37] , On the numerical analysis of analytical nodal methods, Numerical Methods for Partial Differential Equations, 4 (1988), pp. 233-254. [38] J . P . H E N N A R T A N D E . D E L V A L L E , Fast elliptic solvers using mixed-hybrid nodal finite elements, in Advances in Numerical P D E ' s and Optimization, S. Gomez, J . P. Hennart, and R. A . Tapia, eds., Philadelphia, 1991, S I A M , pp. 171-184. [39] , On the relationship between nodal schemes and mixed-hybrid finite elements, Numerical Methods for Partial Differential Equations, 9 (1993), pp. 411-430. [40] A . F . H E N R Y , Refinements in accuracy of coarse-mesh finite-difference solutions of the group-diffusion equations, in Numerical Reactor Calculations, no. I A E A - S M -154/21 in Proceedings Series, Vienna, January 17-21, 1972, Internaltional Atomic Energy Agency, I A E A , pp. 447-471. Bibliography 168 [41] A . F . H E N R Y , ' Nuclear Reactor Analysis, The M I T Press, Cambridge, Mas-sachusetts, 1975. [42] M . H . H O L M E S , Introduction to Perturbation Methods, no. 20 i n Texts i n A p p l i e d Mathematics , Springer-Verlag, New York, 1995. [43] W . C . H O R A K A N D J . J . D O R N I N G , A nodal coarse-mesh method for the efficient numerical solution of laminar flow problems, J . Comput . Phys. , 59 (1985), pp. 405-440. [44] J . J . D O U G L A S , R . D U R A N , A N D P . P l E T R A , Numerical Approximation of Par-tial Differential Equations, North-Hol land, Amsterdam, 1987, ch. Formulation of Alternating-Direct ion Iterative Methods for M i x e d Methods i n Three Space, pp. 21-30. [45] H . K . J O O A N D C . H . K l M , A new approach to core-reflector boundary conditions for nodal reactor computations, N u c l . Sci. Eng . , 116 (1994), pp. 300-312. [46] S. K N A P E K , Matrix-dependent multigridrhomogenization for diffusion problems, in The Prel iminary Proceedings of the Copper Mountain Conference on Iterative M e t h -ods, T . Manteuffal and S. M c C o r m i c k , eds., vol . I, S I A M Special Interest Group on Linear Algebra, Cray Research, A p r i l 9-13 1996. [47] , Matrix-dependent multigrid-homogenization for diffusion problems, Preprint : submitted to S I A M J . Sci. Stat. Comput . , (1996). [48] K . K O E B K E , A new approach to homogenization and group condensation, i n Ho-mogenization Methods i n Reactor Physics, Lugano, Switzerland, November 13-15, 1978, International Atomic Energy Agency, I A E A , pp. 303-323. Bibliography 169 [49] J . G . K O L L A S AND A . F . H E N R Y , The determination of homogenized group diffu-sion theory parameters, N u c l . Sci. Eng . , 60 (1976), pp. 464-471. [50] Y . A . KUZNETSOV AND S. Y . MALIASSOV, Substructuring preconditioners for nonconforming finite element approximations of second order elliptic problems with anisotropy, Russ. J . Numer. A n a l . M o d e l . , 10 (1995), pp. 511-533. [51] E . W . L A R S E N AND G . C . POMRANING, The PN theory as an asymptotic limit of transport theory in planar geometry-I: Analysis, N u c l . Sci. Eng . , 109 (1991), pp. 49-75. [52] R . D . L A W R E N C E , Progress in nodal methods for the solution of the neutron diffu-sion and transport equations, Progress i n Nuclear Energy, 17 (1986), pp. 271-301. [53] , Three-dimensional nodal diffusion and transport methods for the analysis of fast-reactor critical experiments, Progress i n Nuclear Energy, 18 (1986), pp. 101-111. [54] R . D . L A W R E N C E AND J . J . DORNING, A nodal green's function method for multi-dimensional neutron diffusion calculations, N u c l . Sci. Eng . , 76 (1980), pp. 218-231. [55] E . P . E . M l C H E A L , J . J . DORNING, E . M . G E L B A R D , AND RlZWAN-UDDIN, Nodal integration method for the convection-diffusion heat equation, Trans. A m . N u c l . S o c , 69 (1993), pp. 239-241. [56] R . M O S E , P . S lEGEL, P . A C K E R E R , AND G . C H A V E N T , Application of the mixed hybrid finite element approximation in a groundwater flow model: Luxury or neces-sity, Water Resour. Res., 30 (1994), pp. 3001-3012. [57] J . D . M O U L T O N , Diffusion modelling of picosecond pulse propagation in turbid me-dia, Master's thesis, McMaster University, Hamil ton , Ontario, August 1990. Bibliography 170 [58] J . M . N O H A N D N . Z . C H O , A new approach of analytic basis function expansion to neutron diffusion nodal calculation, N u c l . Sci. Eng. , 116 (1994), pp. 165-180. [59] C . R . R A O A N D S. K . M l T R A , Generalized Inverse of Matrices and its Applications, John W i l e y & Sons, Inc., New York, 1971. [60] P . A . R A V I A R T A N D J . M . T H O M A S , A mixed finite element method for 2-nd order elliptic problems, in Mathematical Aspects of the Fini te Element Method , I. Gal l igani and E . Magenes, eds., vol . 606 of Lecture Notes i n Mathematics , Springer Verlag, 1977, pp. 292-315. [61] T . R U S T E N A N D R . W l N T H E R , A preconditioned iterative method for saddlepoint problems, S I A M J . M a t r i x A n a l . A p p l . , 13 (1992), pp. 887-904. [62] E . S A N C H E Z - P A L E N C I A , Non-Homogeneous Media and Vibration Theory, no. 127 i n Lecture Notes in Physics, Springer-Verlag, Ber l in , 1980. [63] D . L . S C H A R F E T T E R A N D H . K . G U M M E L , Large signal analysis of a silicon read diode oscillator, I E E E Trans. Electron Devices, 16 (1969), pp. 64-77. [64] D. S. S E L E N G U T , Diffusion coefficients for heterogenous systems, Nuclear Physics Research Quarterly Report HW-60220, Hanford Laboratories, General Electric , Richland, Washington, A p r i l 1959. [65] K . S . S M I T H , Assembly homogenization techniques for light water reactor analysis, Progress i n Nuclear Energy, 17 (1986), pp. 303-335. [66] K . S. S M I T H , A . F . H E N R Y , A N D R . A . L O R E T Z , The determination of homog-enized diffusion theory parameters for coarse mesh nodal analysis, in Advances i n Reactor Physics and Shielding, Sun Valley, September 14-19, 1980, American N u -clear Society, A N S , pp. 294-308. Bibliography 171 [67] R . M . SMITH AND A . G . H U T T O N , The numerical treatment of advection: A performance comparison of current methods, Numer. Heat Trarisf., 5 (1982), pp. 439-461. [68] J . W . S O N G AND J . K . K l M , An efficient nodal method for transient calculations in light water reactors, Nuclear Technology, 103 (1993), pp. 157-167. [69] G . S T R A N G AND G . J . F i x , An Analysis of the Finite Element Method, Series i n Automat ic Computat ion, Prent ice -Hal l , Inc., Englewood Cliffs, N J , 1973. [70] R . S. V A R G A , Matrix Iterative Analysis, Prent ice -Hal l , Inc., Englewood Cliffs, N J , 1962. [71] E . L . WACHSPRESS, R . D . BURGESS, AND S. B A R O N , Multichannel flux synthesis, N u c l . Sci. Eng . , 12 (1962), pp. 381-389. [72] M . R . W A G N E R , H . F INNEMANN , K . K O E B K E , AND H . - J . W I N T E R , Validation of the nodal expansion method and the depletion program MEDIUM-2 by benchmark calculations and direct comparison with experiment, Atomkernenergie, 30 (1977), pp. 129-135. [73] M . R . W A G N E R , K . K O E B K E , AND H . - J . WINTER, A nonlinear extension of the nodal expansion method, i n Advances in Mathematical Methods for the Solution of Nuclear Engineering Problems, vol . 2, M u n i c h , Germany, A p r i l 27-29 1981, A m e r i -can Nuclear Society, A N S , Mathematics and Computat ion Divis ion , pp. 43-58. [74] A . WEISER AND M . F . W H E E L E R , On convergence of block-centered finite differ-ences for elliptic problems, S I A M J . Numer. A n a l . , 25 (1988), pp. 351-375. [75] D . M . Y O U N G AND K . C . J E A , Generalized conjugate-gradient acceleration of nonsymmetrizable iterative methods, Linear Algebr. A p p l . , 34 (1980), pp. 159-194. Bibliography 172 [76] H . Z H A N G , R l Z W A N - U D D I N , A N D J . J . D O R N I N G , Systematic homogenization and self-consistent flux and pin power reconstruction for nodal diffusion methods - I: Diffusion equation-based theory, Nucl . Sci. Eng. , 121 (1995), pp. 226-244.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Nodal methods : analysis, performance and fast iterative...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Nodal methods : analysis, performance and fast iterative solvers Moulton, J. David 1996
pdf
Page Metadata
Item Metadata
Title | Nodal methods : analysis, performance and fast iterative solvers |
Creator |
Moulton, J. David |
Date Issued | 1996 |
Description | Nodal Methods have long been one of the most popular discretization techniques employed within the reactor physics community, while remaining conspicuously absent from the mainstream numerical analysis literature. A fundamental reason for this anomaly is that the physical arguments which were used to develop and enhance these methods seemed at odds with more rigorous discretization techniques. To facilitate communication between these distinct communities, a detailed chronological study of the lowest-order nodal methods for linear second order elliptic problems is presented. The presentation highlights the underlying motivation of these methods and formalizes many of their renowned properties. In addition, various equivalence relations within this family of discretizations are demonstrated, and equivalences with specific non-conforming and mixed-hybrid finite element methods (FEMs) are established. Rigorous error bounds and stability properties follow immediately from these latter equivalence relations, corroborating the results of a more rudimentary truncation error analysis. An inherent difficulty of reactor simulation is that the coefficients in the mathematical model exhibit severe variations on two significantly different length scales. As in many other applications this is treated by defining an appropriate homogenization procedure which yields a simplified model with piecewise constant coefficients on a coarse scale suitable for efficient computation. Significant enhancements in accuracy are possible if the processes of homogenization and discretization are unified. We review the popular techniques that are based on this premise and rely on certain properties of the nodal discretization. In addition, we address the factors that contribute to their success in reactor modelling and deter their generalization outside of the reactor physics community. As an alternative to these highly specialized methods, we introduce a new multi-level homogenization technique which is readily applicable in a general setting, and is shown to have many important attributes. Widespread acceptance of nodal methods has also been hindered by their use of nonstandard unknowns, as this results in stencils that appear awkward and incompatible with sophisticated iterative solution techniques. Specifically, equivalence with certain mixed-hybrid FEMs reveals that the nodal discretizations result in an indefinite system, which in two dimensions contains both cell-based and edge-based unknowns. Yet, inherent in this structure is a natural partitioning of the system which may be exploited to define a hierarchy of reduced systems (i.e. Schur complements) that are symmetric positive definite. Unfortunately, the reduced systems that involve unknowns of only one type, and hence seem most compatible with sophisticated iterative methods, suffer a loss of sparsity. However, the structure inherent in this hierarchy may be utilized in the construction of sparse approximate Schur complements for these systems. It is demonstrated that any one of these approximate operators, which are of either the standard 5 or 9-point family, may be utilized as an excellent preconditioner for conjugate gradient iterations. The efficiency of this approach is fully realized when the preconditioner is approximately inverted using only a single V - or W-cycle of a robust Black Box multigrid solver. |
Extent | 8985586 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-03-27 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0079927 |
URI | http://hdl.handle.net/2429/6598 |
Degree |
Doctor of Philosophy - PhD |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 1997-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_1997-196275.pdf [ 8.57MB ]
- Metadata
- JSON: 831-1.0079927.json
- JSON-LD: 831-1.0079927-ld.json
- RDF/XML (Pretty): 831-1.0079927-rdf.xml
- RDF/JSON: 831-1.0079927-rdf.json
- Turtle: 831-1.0079927-turtle.txt
- N-Triples: 831-1.0079927-rdf-ntriples.txt
- Original Record: 831-1.0079927-source.json
- Full Text
- 831-1.0079927-fulltext.txt
- Citation
- 831-1.0079927.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0079927/manifest