A Nonlinear Krylov Accelerator for the Boltzmann k-Eigenvalue Problem MatthewT.Calef∗,ErinD.Fichtl,JamesS.Warsa,MarkusBerndt,NeilN.Carlson ComputationalPhysicsandMethods,LosAlamosNationalLaboratory,LosAlamos,NM87545-0001 Abstract 2 We compare variants of Anderson Mixing with the Jacobian-Free Newton-Krylov and Broyden methods applied to 1 0 the k-eigenvalue formulation of the linear Boltzmann transport equation. We present evidence that one variant of 2 AndersonMixingfindssolutionsinthefewestnumberofiterations. Weexamineandstrengthentheoreticalresultsof AndersonMixingappliedtolinearproblems. n a J 6 1. Introduction ] h p Thek-eigenvalueformulationofthelinearBoltzmanntransportequationiswidelyusedtocharacterizethecriti- - calityoffissioningsystems[1,2]. Physically,thelargesteigenvalue,generallydenotedbyk,istheeffectiveneutron p multiplication factor that, in the equation, scales the fission production term to achieve a steady-state solution. The m correspondingeigenmodedescribestheneutronfluxprofileforthatsteady-state(i.e.,critical)systemand,whenthe o systemisclosetocriticality,providesusefulinformationaboutthedistributionoftheneutronpopulationinspaceand c velocity[1].Mathematically,theequationisastandardeigenproblemforwhichpoweriterationiswell-suitedbecause . s theeigenmodeofinterestismostcommonlythatwiththelargestmagnitude[3]. Forthedeterministick-eigenvalue c i problemeachstepofatruepoweriterationincursaheavycomputationalcostduetotheexpenseoffullyinvertingthe s y transportoperator,thereforeanonlinearfixedpointiterationisgenerallyemployedinwhichanapproximateinversion h ofthisoperatorisperformedateachiteration. Inadditiontopoweriteration,theImplicitlyRestartedArnoldiMethod p hasbeenappliedtothisproblemandhastheadvantageofbeingabletocomputeadditionaleigenmodes[4].However, [ thetransportmatrixmustbefullyinvertedateachiteration,diminishingitscomputationalefficiencyandattractiveness 3 whenonlythedominanteigenmodeisdesired. v Recently,moresophisticatednonlineariterationmethodsemployingapproximateinversion,predominantlyJacobian- 8 FreeNewton-Krylov(JFNK),havebeenappliedwithgreatsuccess[5,6,7]. However,therehasnotyetbeenacom- 6 prehensivecomparisonofnonlinearsolvers. Thispaperpresentssuchacomparison, examiningtheperformanceof 5 threenonlinearsolvers—JFNK,Broyden’sMethodandAndersonMixing—appliedtothek-eigenvalueproblem. A 3 . variantofAndersonMixing[8],firstdescribedin[9],isofparticularinterestbecause,intheexperienceoftheauthors, 2 itisfrequentlycomputationallymoreefficientthanJFNKandBroyden’smethod. 1 1 JFNK is an inexact Newton’s method in which the inversion of the Jacobian is performed to arbitrary precision 1 using a Krylov method (most commonly GMRES) and the Jacobian itself is never formed, but rather its action is : approximatedusingfinitedifferencesofarbitrarilyclosestatedata. JFNKcanbeexpectedtoconvergequadratically v i in a neighborhood containing the solution (cf. [10] and the references therein). Each iteration of JFNK requires a X nested ‘inner’ iteration and the bulk of the computational effort is expended in this ‘inner’ Krylov inversion of the r Jacobianateach‘outer’Newtonstep. Attheendofeachinversion,theaccumulatedKrylovspaceisdiscardedeven a thoughtheJacobianisexpectedtochangeminimallyduringthefinalNewtonstepswhenasimilarspacewillberebuilt inthenextNewtoniteration. Ineffect,attheendofeachiteration,JFNKdiscardsinformationthatmaybeofusein successiveiterations. ∗Correspondingauthor Emailaddresses:[email protected](MatthewT.Calef),[email protected](ErinD.Fichtl),[email protected](JamesS.Warsa), [email protected](MarkusBerndt),[email protected](NeilN.Carlson) PreprintsubmittedtoElsevier January9,2012 In its standard formulation Broyden’s method (cf. [11]), like low memory BFGS (cf. [12]), uses differences in state from successive iterations to make low rank updates to the Jacobian. The Sherman-Morrison-Woodbury up- date rule is then used to compute the action of the inverse of the Jacobian after such an update. While Broyden’s methodisrestrictedtolow-rankupdates,itprovidesanexplicitrepresentationoftheJacobianallowingonetoemploy the Dennis-More´ condition [13] to show that it converges super-linearly in a neighborhood containing the solution. Further, it has been shown to solve linear problems of size N in at most 2N iterations (cf. [11] and the references therein.) Anderson Mixing [8] uses differences in state from successive iterations to infer information about the inverse of the Jacobian, which is assumed to be roughly constant in a neighborhood containing all the iterates. Unlike Broyden’smethod, theupdatescanbeofarbitraryrank. RecentresultsbyWalkerandNi[14]showthat, withmild assumptions, Anderson Mixing applied to a linear problem performs as well as the generalized minimum residual method (GMRES) [15]. In this regard, Anderson Mixing may be thought of as a nonlinear version of GMRES. In independent work, Carlson and Miller formulated a so-called nonlinear Krylov acceleration [9] method which we show to be a variant of Anderson Mixing. Further, we examine the hypotheses of a central theorem presented by WalkerandNiandarguethattheywill, withhighprobability, bemetandthattheycanbeomittedentirelyifoneis willingtoacceptasmallperformancepenalty. WhileAndersonMixingperformsbestforournumericalexperiments, thereisnotheorythattheauthorsofthispaperknowofthatcancharacterizeitsperformancefornonlinearproblems. Therestofthispaperisorganizedasfollows: InSection2wediscusssomeoftheanalysisofAndersonMixing, motivateanimplementationchoiceandelaborateontheresultsofWalkerandNiin[14]. InSection3wedescribethe k-eigenvalueproblem. OurnumericalresultscanbefoundinSection4andweconcludewithanoverview. 2. BackgroundofAndersonMixingandNonlinearKrylovAcceleration 2.1. NonlinearKrylovAcceleration In [9, 16] Carlson and Miller outlined an iterative method, dubbed nonlinear Krylov acceleration or NKA, for acceleratingconvergenceoffixed-pointiterationbyusinginformationgainedoversuccessiveiterations. Theproblem theyconsideristofindarootofthefunctionRN (cid:51)x→ f(x)∈RN. Oneapproachistoapplyafixed-pointiterationof theform xn+1 =xn− f(xn). (1) Fixed-point iterations can be viewed in the context of an approximate Newton iteration where the derivative of f, denotedDf,isreplacedbytheidentityI.ThemotivationbehindNKAisinsteadtoapproximateDf usinginformation from previous iterates, improving that approximation over successive iterations, and in cases where no applicable approximationisavailable,reverttoafixed-pointiterationwhereDf isreplacedwithI. NKArequiresaninitialguessx0andatthenthinvocationprovidesanupdatevn+1thatisusedtoderivethen+1st iteratefromthenth. Thismethodmaybewrittenas vn+1 = NKA[f(xn),...] xn+1 = xn−vn+1, where NKA[f(x ),...]istheupdatecomputedbythenonlinearKrylovaccelerator. Weusethebracketsandellipsis n toindicatethatNKAisstatefulanddrawsonpreviousinformationwhencomputinganupdate. Onitsfirstinvocation(n = 0)NKAhasnoinformationandsimplyreturns f(x ). Atiterationn > 0ithasaccess 0 totheMvectorsofdifferencesforsomenaturalnumberM ∈(0,n), v =x −x and w = f(x )− f(x) fori=n−M+1,...,n, i i−1 i i i−1 i and where, for convenience, we shall define W to be the span of the w vectors available at iteration n. If f has a n i constantandinvertiblederivative,Df,thenwewouldhave Dfv =w and Df−1w =v. (2) i i i i 2 WedenotebyP theoperatorthatprojectsorthogonallyontothesubspaceW andwritetheidentity Wn n f(x )=P f(x )+(f(x )−P f(x )). n Wn n n Wn n Notethat f(x )−P f(x )isorthogonaltoW . Ifthew vectorsarelinearlyindependent,thenthereisauniqueset n Wn n n i ofcoefficientsz(n) :=(z(n),z(n),...,z(n))∈Rnsothat 1 2 n (cid:88)n P f(x )= z(n)w, Wn n i i i=n−M+1 andhencebyEq.(2) (cid:88)n (cid:88)n Df−1P f(x )= Df−1 z(n)w = z(n)v. Wn n i i i i i=n−M+1 i=n−M+1 TheideamotivatingCarlsonandMilleristoproject f(x )ontoW (thespacewheretheactionof Df−1 isknown), n n computethatactionontheprojectionand,forlackofinformation,applyafixed-pointupdategiveninEq.(1)onthe portionof f(xn)thatisorthogonaltoWn. Theresultingformulaforxn+1is    vn+1 = (cid:88)n z(in)vi+f(xn)− (cid:88)n z(in)wi, i=n−M+1 i=n−M+1 xn+1 =xn−vn+1, wherethevectorofcoefficientsintheorthogonalprojection,z(n),isthesolutiontotheprojection,alternativelymini- mization,problem (cid:13) (cid:13) (cid:13)(cid:13) (cid:88)n (cid:13)(cid:13) z(n) =argminy∈RM(cid:13)(cid:13)(cid:13)f(xn)− yiwi(cid:13)(cid:13)(cid:13) . (cid:13) i=n−M+1 (cid:13)2 2.2. NKAasAndersonMixing CarlsonandMillerhadessentiallyrediscovered,albeitinaslightlydifferentform,theiterativemethodpresented muchearlierbyAndersonin[8],whichisnowcommonlyreferredtoasAndersonMixingorAndersonAcceleration. Anderson was studying the problem of finding a fixed point of some function RN (cid:51) x → G(x) ∈ RN. In Section 4 of[8]hedefinesafixedpointresidualforiterationias r =G(x)−x, i i i whichcanbeusedtomeasurehowx failstobeafixedpoint. WiththisAndersonproposedanupdatingschemeof i theform1    xn+1 =xn−(cid:88)M z˜(in)(xn−xn−i)−βnrn−(cid:88)M z˜(in)[rn−rn−i] i=1 i=1 wherethevectorofcoefficients,z˜(n) ∈RM,ischosentominimizethequantity (cid:13) (cid:13) (cid:13)(cid:13) (cid:88)M (cid:13)(cid:13) (cid:13)(cid:13)(cid:13)rn− z˜(in)[rn−rn−i](cid:13)(cid:13)(cid:13) , (cid:13) i=1 (cid:13)2 andwhereAndersonrequiresthatβ >0. HereAndersonconsidersadepthofMdifference-vectors. n 1Andersonusesθntodenotetheithupdatecoefficientofthenthiterate,wherewehavewrittenz˜(n)inkeepingwiththepresentationofNKA. i i 3 Theproblemoffindingarootof f mayberecastasafixedpointproblembydefiningG(x) = x− f(x),andthen r =−f(x). If,foreachn,onedefinesthevectors i i v˜(n) =x −x and w˜(n) = f(x )− f(x ) fori=1,...,M, i n n−i i n n−i thentheAndersonMixingupdatebecomes    xn+1 =xn−(cid:88)M z˜(in)v˜(in)+βnf(xn)−(cid:88)M z˜(in)w˜(in), i=1 i=1 where (cid:13) (cid:13) (cid:13)(cid:13) (cid:88)M (cid:13)(cid:13) z˜(n) =argminy∈RM(cid:13)(cid:13)(cid:13)f(xn)− yiw˜(in)(cid:13)(cid:13)(cid:13) . (cid:13) i=1 (cid:13)2 Note that span{w˜(1n),...,w˜(Mn)} = span{wn−M+1,...,wn} = Wn. In particular the update coefficients of NKA and AndersonMixingarerelatedbytheformula (cid:88)M z(n) =− z˜(n). (3) i j j=n−i+1 Further,foraconstantandinvertiblederivative, Df−1w˜(n) = v˜(n). Withthisitisclearthat,forlinearproblems, NKA i i is equivalent to Anderson Mixing when β = 1; the difference is only in the choice of basis vectors that span W n n resultinginachangeofvariablesgivenbyEq.(3). OneadvantageoftheNKAformulationisthatitusesdifferences ofsuccessivefunctionevaluationstoformabasisforW . Thesedifferencescanbeusedfor Msuccessiveiterations. n In contrast, Anderson’s formulation forms a basis for W using differences of the most recent function evaluation n withallpreviousevaluations,whichmustberecomputedateveryiteration. 2.3. AnalyticExaminationsofAndersonMixingAppliedtoaLinearProblem In [14], Walker and Ni consider the behavior of a generalized Anderson Mixing algorithm when it is applied to linear problems and when all previous vectors are stored, i.e. when M = n. They prove that, under modest assumptions,aclassofvariantsofAndersonMixingareequivalenttoGMRES.TheyconsiderAndersonMixingasa meanstofindafixed-pointofG. LikeAndersontheycomputer =G(x)−x ateachstep. i i i IntheirpresentationofAndersonMixing,though,theycomputetheupdatecoefficients2. (cid:13) (cid:13) (cid:13)(cid:13)(cid:88)n (cid:13)(cid:13) (cid:88)n z˜(n) =argminy∈Rn+1(cid:13)(cid:13)(cid:13) yiri(cid:13)(cid:13)(cid:13) subjectto z˜(in) =1, (cid:13)i=0 (cid:13)2 i=0 whicharethenusedtoformtheupdate (cid:88)n xn+1 = z˜(in)G(xi). i=0 AgainchoosingG(x)=x− f(x)(andhencer =−f(x))andnotingthattheconstraintrequiresthat i i (cid:88)n−1 z˜(n) =1− z˜(n), n i i=0 2WalkerandNiused fjandα(nj)wherewe,forconsistency,areusingrjandz˜(jn). 4 onehas,asnotedbyWalkerandNi,thattheaboveisequivalenttoAnderson’soriginalformulation (cid:13) (cid:13) (cid:13)(cid:13) (cid:88)n−1 (cid:13)(cid:13) z˜(n) =argminy∈Rn(cid:13)(cid:13)(cid:13)f(xn)− yi(f(xn)− f(xi))(cid:13)(cid:13)(cid:13) (cid:13) i=0 (cid:13)2 and    xn+1 =xn−(cid:88)n−1z˜(in)(xn−xi)+f(xn)−(cid:88)n−1z˜(in)(f(xn)− f(xi)). i=0 i=0 Because Anderson Mixing and NKA are equivalent to each other and to the Walker and Ni formulation, we may restateoneofWalkerandNi’stheoremsasfollows: Theorem 2.1 (Walker and Ni [14]). Suppose that we search for a root of the function f(x) = Ax−b starting with x using either the NKA update or Anderson’s update with β = 1. Suppose further that A is non-singular. Let 0 n xGMRES denote the nth iterate generated by GMRES applied to the problem Ax = b with starting point x and let n 0 rGMRES = AxGMRES −bdenotetheassociatedresidual. Ifforsomen > 0,rGMRES (cid:44) 0and(cid:107)rGMRES(cid:107) < (cid:107)rGMRES(cid:107) n n n−1 j 2 j−1 2 forall0< j<n,thenthen+1iterategeneratedbyeitherupdateruleisgivenbyxn+1 =(I−A)xGnMRES +b. It should be noted that the original presentation of the theorem given in [14] applies to a class of orthogonal projectionmethodswherewehavepresentedthetheoremasappliedtotwomembersofthisclass:Anderson’soriginal methodandtheNKAvariant. ItshouldalsobenotedthatthistheoremofWalkerandNishowsthat W =AK , n n whereK =span{r ,Ar ,...,An−1r }. n 0 0 0 Animmediatecorollaryis Corollary 2.2. Let r = Ax −b denote the residual associated with the nth iterate generated by either Anderson n n MixingorNKA.UndertheassumptionsofTheorem2.1 (cid:107)rn+1(cid:107)2 ≤(cid:107)I−A(cid:107)(cid:107)rGnMRES(cid:107)2, where(cid:107)·(cid:107)istheL2operatornorm. Proof. rn+1 = Axn+1−b (cid:104) (cid:105) = A (I−A)xGMRES +b −b n = AxGMRES −b−A(AxGMRES −b) n n = (I−A)(AxGMRES −b), n fromwhichtheclaimfollows. (cid:3) ThevalueofCorollary2.2isthat,inthelinearcase,convergenceofthenon-truncatedversionsofbothNKAand AndersonMixinghavethesamecharacterizationsasGMRES,i.e. whenconsideringanormalmatrixA,theresidual iscontrolledbythespectrumofA,providedasGMRESdoesn’tstagnate. 2.4. Non-stagnationofGMRES When considering a linear problem, the coefficients z˜(n) (Anderson Mixing) and z(n) (NKA) do two things. Through a minimization process, they select the best approximation within a given subspace, and they also select thesubspaceforthenextiteration. Thefailuremodeisthatthebestapproximationwithinagivensubspacedoesn’t use the information from the most recent iteration, in which case the next subspace will be the same as the current oneandAndersonMixingandNKAwillbecometrapped. Andersonbrieflydiscussesthisinhispresentationofthe 5 method3.Lemma2.1.5fromNi’sdissertation[17]addressesthismoredirectly.WhatWalkerandNirecognizein[14] isthatthisfailuremodecorrespondstothestagnationofGMRES. There have been several examinations of when GMRES stagnates. Zavorian, O’Leary and Elman in [18] and Greenbaum, Ptak and Strakos in [19] present examples where GMRES does not decrease the norm of the residual on several successive iterations, i.e. it stagnates. While GMRES converges for such cases, Anderson Mixing and NKA will not. Greenbaum and Strakos show in [20] that the residual for GMRES applied to a matrix A is strictly decreasing,i.e. nostagnation,if(cid:104)r,Ar(cid:105) (cid:44) 0forallrsatisfying(cid:107)r(cid:107) (cid:44) 0. Thisinturnwillensuretheconvergenceof 2 AndersonMixingandNKA. Itisimportanttobearinmindthat,inpractice,thesubspaceW islikelytohavemodestdimension,say10,and n issittinginRN where N ismuchlarger, oftenontheorderofthousandsormillions. Theslightestperturbationofa vectorwill,withprobabilityone,pushitoffofthislowdimensionalspace. Evenso,thereisasimplechangetothe update step that provides a theoretical guarantee that NKA will still converge even when GMRES stagnates. This guaranteecomesataslightperformancecostandisaccomplishedbymodifyingtheupdatecoefficientstoensurethat W (cid:40)W . n−1 n WeconsideramodificationtotheNKAupdaterule4:    xn+1 =xn−(cid:88)n viz(in)+f(xn)−(cid:88)n wiz(in), i=1 i=1 wherez(n)ischosentominimize (cid:13) (cid:13) (cid:13)(cid:13) (cid:88)n (cid:13)(cid:13) (cid:13)(cid:13)(cid:13)f(xn)− wiz(in)(cid:13)(cid:13)(cid:13) . (cid:13) i=1 (cid:13)2 Wenowaddthefollowingsafetycheck: if(cid:107)w (cid:107) (cid:44)0,thenperformthefollowingupdate n 2 z(n) ←z(n)±ε(cid:13)(cid:13)(cid:13)f(xn)−(cid:80)ni=1wiz(in)(cid:13)(cid:13)(cid:13)2 n n (cid:107)w (cid:107) n 2 forsomepositiveε,whereonechoosestoaddorsubtractbasedonwhichwillmakethemodifiedversionof|z(n)+1| n furtherfromzero. AswillbeshownintheproofofTheorem2.3thenumeratorofthemodificationiszeroonlywhen NKAhasfoundtheroot. WiththiswecanstrengthenCorollary2.2asitappliestoNKAasfollows: Theorem2.3. LetAbenon-singularsquarematrixandsupposethatwesearchforarootofthefunction f(x)=Ax−b startingwithx usingthemodifiedversionofNKA.LetrGMRES = AxGMRES −bdenotethenth GMRESresidualand 0 n n letr denotethenthNKAresidual,then n (cid:107)rn+1(cid:107)2 ≤(1+ε)(cid:107)I−A(cid:107)(cid:107)rGnMRES(cid:107)2. TheproofcanbefoundinAppendix A.Thelimitationofthisresultisthattheorthogonalprojectionwillbecome morepoorlyconditionedasεdecreases. NotethatforthesamereasonbothAndersonMixingandNKAcandevelop arbitrarilypoorlyconditionedprojectionproblems. 2.5. RelationshipBetweenAndersonMixingandBroyden FangandSaadin[21]provideacomparisonbetweenseveralmethodsincludingAndersonMixingandBroyden. Asnotedabove,thecentraldifferencebetweenthetwoisthatBroydenusessecantinformationtocomputelowrank updates to the approximation of Df, where Anderson Mixing uses the same information to compute arbitrary rank updatesto Df−1. FangandSaadclarifythisrelationship andplaceavariantof BroydenandAndersonMixingin a broaderclassofmulti-secantupdatemethods. 3Thisneedtoexpandthespaceoverwhichoneissearchingnecessitatesβn(cid:44)0. 4ThefollowingcouldbeadaptedtomanyformulationsofAndersonMixing 6 Fang’sandSaad’spresentationofAndersonMixing(cf. Eq. (24)in[21])inthenotationadoptedinthispaperis    xn+1 =xn− (cid:88)n viz(in)−β˜f(xn)− (cid:88)n wiz(in), i=n−M+1 i=n−M+1 where (cid:13) (cid:13) (cid:13)(cid:13) (cid:88)n (cid:13)(cid:13) z(n) =argminy∈RM(cid:13)(cid:13)(cid:13)f(xn)− wiyi(cid:13)(cid:13)(cid:13) , (cid:13) i=n−M+1 (cid:13)2 andwhereβ˜ ischosentobepositive. It is worth noting that Fang and Saad have presented a version of Anderson Mixing that is more closely related toCarlson’sandMiller’sNKAmethod,butwhereNKAwouldrequirethat,intheaboveformulation,β˜ = −1. This leadstothequestionwhatisthebestvalueofβ˜ foragivenproblem,orinthecaseofAnderson’spresentationwhat isthebestchoiceforβ = −β˜. Wereportnumericalresultsforβ˜ = −1(theNKAformulation)andβ˜ = 1(denoted n NKA ). −1 2.6. NKAinPractice ForourexperimentsweshallusetheNKAformulationofAndersonMixing. Weshallletβ˜ = 1andβ˜ = −1,as definedbyFangandSaad. Eachw vectorisrescaledtohaveunit-norm,andtheassociatedv vectorisscaledbythe i i same factor. In our implementation we use a Cholesky Factorization to solve the least-squares problem associated withtheorthogonalprojection. ThisrequiresthattheGramianofthew vectorsbepositivedefinite. Weenforcethis i bychoosingalinearlyindependentsubsetofthew vectorsasfollows: Atthebeginningofiterationnwestartwith i w andconsiderw . Iftheanglebetweenw andw islessthansometolerancewediscardw . Weiteratein n n−1 n n−1 n−1 thismanneroversuccessivelyoldervectors,keepingthemiftheangletheymakewiththespacespannedbythekept w vectorsisgreaterthanthetolerance. ThisensuresthatwemayuseCholeskyfactorizationandthatthecondition i numberoftheproblemisbounded. Memoryconstraintsrequirethat M <n,forcingustouseNKA,NKA andBroydeninasettingforwhichthere −1 arenotheoreticalresultsthattheauthorsofthispaperknowof. OneimportantdistinctionbetweenBroydenandNKA isthatforeachiterationNKAstoresapairofstatevectors,whileBroydenonlystoresone. Consequentlythememory requirementsforBroydenarehalfthatofNKAforthesamedepthofstoredvectors. DependingonthetolerancechosenateachlinearsolveforJFNK,onecanreasonablyexpecttheoreticalguarantees ofquadraticconvergencetoholdinsomeneighborhoodofthesolution,howeveridentifyingthatneighborhoodisoften considerablyharderthansolvingtheoriginalproblem. Insummary,forthenumericalexperimentswepresent,there islittletheoryregardingperformance,makingthefollowingnumericalresultsofvalue. 3. Thek-EigenvalueFormulationoftheBoltzmannTransportEquation Theproblemthatweusetocomparetheefficiencyofthesevariousnonlinearmethodsisthek-eigenvalueformu- lationofthelinearizedBoltzmannneutrontransportequation. (cid:90) (cid:90) (cid:90) (cid:104) (cid:105) 1 Ωˆ ·∇(cid:126) +Σ((cid:126)r,E) ψ((cid:126)r,Ωˆ,E)= dE(cid:48) dΩ(cid:48)Σ ((cid:126)r,E(cid:48) → E,Ωˆ(cid:48)·Ωˆ)ψ((cid:126)r,Ωˆ,E(cid:48))+ dE(cid:48)χ(E(cid:48) → E)ν¯Σ ((cid:126)r,E(cid:48))φ((cid:126)r,E(cid:48)). t s k f (4) Theunknownψdescribestheneutronfluxinspace,(cid:126)r ∈ R3,angleΩˆ ∈ S2 andenergyE. Σ isthetotalcrosssection, t Σ isthescatteringcrosssectionandΣ isthefissioncrosssection. Thequantitiesν¯ andχcharacterizetherateand s f energy spectrum of neutrons emitted in the fission process. Integrating the flux over angle gives φ, the scalar flux at a given position and energy. For a thorough examination of the mathematical models of fission, including the Boltzmanntransportequation,cf.[1,2]. DiscretizationofEq.(4)isaccomplishedusing 7 1. S (discreteordinates)inangle: AnS levelsymmetricquadraturesetischosenandthesolutioniscomputed N 4 attheabscissasofthatsetandthenintegratedoverangleusingthequadratureweights. 2. Multigroupinenergy: Crosssectionscanbeverynoisyanditisthereforenotpracticaltodiscretizetheenergy variable as one would a continuous function. Instead, the energy space is divided up into groups and it is assumedthattheenergycomponentforacrosssectionσcanbeseparatedout: σ((cid:126)r,E) ≈ f(E)σ ((cid:126)r), E < E ≤ E g g g−1 (cid:82)Eg−1dEσ((cid:126)r,E) σ = Eg g (cid:82)Eg−1dEf(E) Eg 3. Finiteelementordifferenceinspace: Thespatialvariablecanbetreatedinavarietyofdifferentways. Herewe exploreresultsfromtwodifferentdiscretizationstrategiesasemployedbytheLosAlamosNationalLaboratory productiontransportcodesPARTISN[22]andCapsaicin. PARTISNisusedtogenerateresultsonastructured mesh with diamond (central) difference and Capsaicin to generate results on an unstructured polygonal mesh usingdiscontinuousfiniteelements. ApplyingtheS andmultigroupapproximations,thelinearizedBoltzmannequationtakestheform N (cid:16)Ωˆ ·∇+σ ((cid:126)r)(cid:17)ψ ((cid:126)r)= 1 (cid:88)G σ ((cid:126)r)φ ((cid:126)r)+ 1 (cid:88)G ν¯σ ((cid:126)r)χg(cid:48)→g((cid:126)r)φ ((cid:126)r). (5) m t,g g,m 4π s,g(cid:48)→g g(cid:48) k f,g(cid:48) 4π g(cid:48) g(cid:48)=1 g(cid:48)=1 Here, • ψ represents the angular flux in direction Ωˆ in energy group g, which is the number of neutrons passing g,m (cid:16) m (cid:17) throughsomeunitareaperunittime # cm2·Mev·ster·sec • φ isthescalarflux,orangle-integratedangularflux. TheS quadratureusedtodefineangularabscissas(Ωˆ ) cagnbeusedtointegrateψoverallangle: φ =(cid:80)M w ψ Nwherew arethequadratureweights(cid:16) # m(cid:17) g m=1 m g,m m cm2·Mev·sec (cid:16) (cid:17) • σ isthetotalcrosssection,orinteractionprobabilityperarea,forgroupg # t,g cm2 • σ isthe‘inscatter’crosssection,whichistheprobabilityperareathataneutronwillscatterfromgroupg(cid:48) s,g(cid:48)→g (cid:16) (cid:17) intogroupg # cm2 (cid:16) (cid:17) • σ isthefissioncrosssectionforgroupg # f,g cm2 • ν¯ isthemeannumberofneutronsproducedperfissionevent • χ describestheenergyspectrumofemittedfissionneutronsingroupgproducedbyneutronsabsorbedfrom g(cid:48)→g groupg(cid:48) Applicationofthespatialdiscretizationthenyieldsamatrixequation.Forconvenience,thisequationcanbeexpressed inoperatornotationas 1 Lψ=MSDψ+ MFDψ. (6) k whereListhestreamingandremovaloperator,Sisthescatteringoperator,Fisthefissionoperator,andMandDare themoment-to-discreteanddiscrete–to–momentoperators,respectively,andDψ=φ. Ifthefluxesarearrangedbyenergygroup,theformofLisblockdiagonalanditcanbeinvertedontoaconstant right-hand-sideexactlybyperformingaso-calledsweepforeachenergygroupandangularabscissa. Thesweepisan exactapplicationofL−1andcanbethoughtofastrackingthemovementofparticlesthroughthemeshalongasingle directionbeginningattheincidentboundaries,whichvarybyangle,sothatallnecessaryinformationaboutneutrons streaming into a cell from other ‘upwind’ cells is known before it is necessary to make calculations for that cell. 8 Mosttransportcodescontainthemechanismtoperformthesweep,butcouldnotapplyLdirectlywithoutsignificant modification,thereforeEq.(6)isgenerallythoughtofasastandardeigenproblemoftheform kφ=(cid:16)I−DL−1MS(cid:17)−1DL−1MFφ. Becauseweseekthemodewiththelargestmagnitudeeigenvalue,poweriterationisonepossibleiterativetechnique: φz+1 = (cid:16)I−DL−1MS(cid:17)−1DL−1MFφz kz+1 = WWTTFFφφz+1 z (cid:16) (cid:17) Full inversion of I−DL−1MS is expensive, however, so it is generally more efficient to use a nonlinear fixed (cid:16) (cid:17) point iteration (FPI) in which the operator I−DL−1MS is inverted only approximately using a nested inner-outer iterationschemetoconvergethescatteringterm. Commonly, oneormore‘outeriterations’areconductedinwhich thefissionsourceislagged. Eachouteriterationconsistsofaseriesofinner‘within-group’iterations(oneormore perenergygroup)inwhichtheinscatterislaggedsothatthewithin-groupscatteringcanbeconverged. Ingeneral, theiterationcanbewrittenas φz+1 = P(kz)φz (8a) kz+1 = kzWWTTFFφφz+1 (8b) z whereW isavectorofweightsandgenerallyrepresentsavolumeintegral. Thek-eigenvaluecanbethoughtofasthe ratioofneutronsinsuccessivegenerationsandnewneutronsarecreatedinfissionevents,thereforeitmakessenseto adjust k using the ratio of fission sources. If P(k ) = (cid:16)I−DL−1MS(cid:17)−1DL−1MF we recover a true power iteration, z buttherearenumerouspossibleformsforthisoperator. Itistypicallymoreefficienttolimitthenumberofsweeps, however. ThemostbasicapproachistoapplyL−1 viaasweep. Inlightofthefactthatwearelaggingthescattering altogether,wealsoconsideranotherupdatefortheeigenvaluethataccountsforthechangeinthescatteringaswell. (cid:32) (cid:33) 1 φz+1 = DL−1M S+ k F φz (9a) z kz+1 = (cid:16) WTFφz+1 (cid:17). (9b) WT k1zFφz−S(φz+1−φz) Clearly,ifweassumethatthescatteringisconverged,thenthesecondterminthedenominatorofEq.(9b)iszeroand werecoverEq.(8b). Foraconvergedsystem,thistermwillindeedbezeroandtheratioofthefissionsourceswillgo toonebecausethefissionsourcehasbeensuitablyadjustedby 1 sothatthenetneutronproductioniszero. k FromthiswehavethatafixedpointoftheiterationinEq.(9)isalsoarootoftheresidualfunction (cid:32) (cid:33)   F φk = 1−PW(TkF)φΦ(φ)k , (10) WTFφ where (cid:32) (cid:32) (cid:33)(cid:33) 1 P(k)= I−DL−1M S+ F k In order to initialize the iteration, a single sweep is performed on a vector of ones, E = [11...1]T, and the scaled two-normofthefluxisthennormalizedtoone: DL−1M(S+F)E φ0 = (cid:13)(cid:13)(cid:13)DL−1M(S+F)E(cid:13)(cid:13)(cid:13) 2,s k = 1 0 9 Here (cid:107)·(cid:107) indicates the L2-norm normalized by the square root of the number of degrees of freedom. Once the 2,s residual is formulated, it is possible to apply any of the nonlinear solvers discussed to seek a root of F given in Eq. (10). The following section contains a comparison of JFNK, Broyden, NKA and NKA for the k-eigenvalue −1 problem. 4. Results Resultsaregivenforacylinderwitha3.5cmradiusandaheightof9cmmodeledintwo-dimensionalcylindrical coordinates,similartothecylindricalproblemthatwasstudiedin[4]inthree-dimensionalCartesiancoordinates. The problemconsistsofacentral5-cmlayerofBoron-10with1cmthickwaterlayersoneithersideand1cmlayersof highlyenricheduraniumontheends. Thetop,bottomandradialboundariesareeitherallvacuumorallreflective(the innerradialboundaryisasymmetryconditionincylindricalcoordinates). Thisa‘difficult’problem(i.e.,onewitha highdominanceratio)becausetheproblemissymmetric, havingtwofissileregionsthatareeffectively‘decoupled’ becauseofthereflector(water)andabsorber(boron)betweenthem,whichmeansthesecondeigenmodeisclosetothe fundamentalmode;reflectionfurtherincreasesthedifficultyoftheproblem. A16-groupHansen-Roachcrosssection data set is used to generate the results presented here (the 16 group data was collapsed to 5 groups in the original paper[4]). A 175 (r-axis) by 450 (z-axis) mesh of equally sized squares is used in PARTISN for both the reflected and unreflected problems. The Capsaicin results are computed on an unstructured mesh comprising roughly the same numberofcellsintherandzaxesasthePARTISNcomputation,foratotalof79855(possiblynon-convex)polygons, eachofwhichhas3to6vertexesonacell.Thesecodesusesimilarmethods,butthereareseveraldifferencesthatlead todifferentiterativebehaviorandslightlydifferentresults. Firstly,PARTISNisanorthogonalmeshcodethatutilizes a diamond (central) difference (DD) spatial discretization, requiring the storage of only one cell-centered value per energygroup[2].Capsaicinusesanunstructuredmesh,whichhastheadvantageofbeingabletomodelgeometrically complexproblemswithhigherfidelity,butitcomesatacost. Findinganefficientsweepscheduleonanunstructured mesh that minimizes the latency in parallel computations is a difficult problem in itself [23, 24, 25]. In contrast, a parallel sweep schedule that is nearly optimal is easy to implement for structured meshes [26]. Furthermore, a discontinuousfiniteelementmethod(DFEM)spatialdiscretizationisemployedinCapsaicinsuchthatthenumberof unknownspercellisequaltothenumberofnodes[27]. WhiletheDFEMhasbetteraccuracythanDD,potentially enablingtheuseofcommensuratelyfewermeshcellsforthesamesolutionaccuracy,itismorecostlytocomputethe solutiontotheDFEMequationsthantheDDequations,andtheDFEMisthuslessefficientthanDDwhencalculations areperformedonmeshesofsimilarsize. Theseconddifferencebetweenthetwocodesisthatreflectingboundaryconditionscancauseinstabilitieswiththe DDspatialdiscretizationandtheangulardifferencingusedinPARTISNso,inordertoachievestability,thereflected fluxis‘relaxed’. Thisisdoneusingalinearcombinationofthenewandoldreflectedfluxesastheboundarysource forthenextiteration: ψz+1 =rψz +(1−r)ψz+1. relaxed refl refl Therelaxationparameter,r,is,bydefault, 1 forthisproblem. ThisrelaxationisnotnecessaryinCapsaicinbecause 2 it uses a more accurate spatial and angular discretization, including the use of starting directions for calculations in cylindrical coordinates [2], as well the use of reflected starting directions when reflection conditions are specified on the radial boundary. The eigenvalues computed by PARTISN are k = 0.1923165 and k = 0.8144675 for the unreflected and reflected cases, respectively. Because the codes employ different discretizations, the eigenvalues calculated by Capsaicin were k = 0.191714 for the unreflected case and k = 0.814461 for the reflected case. The Capsaicindiscretizationresultsinroughly50milliondegreesoffreedomwhilethePARTISNdiscretizationresultsin alittleover5milliondegreesoffreedom. ThethreenonlinearsolversthatweconsiderinthispaperareimplementedintheNOXnonlinearsolverpackage thatispartofthelargersoftwarepackageTrilinos10.6[28].WeuseJFNKasitisimplementedinNOXandhavewrit- tenimplementationsofbothNKAandBroyden’smethodintheNOXframeworkasusersupplieddirectionclasses. These two direction classes are available for download on Sourceforge at http://sourceforge.net/projects/ 10