Nonlinear Krylov Acceleration Applied to a Discrete Ordinates Formulation of the k-Eigenvalue Problem

0.34 MB·English

by  Matthew T. Calef

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