La lecture en ligne est gratuite
Le téléchargement nécessite un accès à la bibliothèque YouScribe
Tout savoir sur nos offres
Télécharger Lire

High performance computing for stability problems [Elektronische Ressource] : applications to hydrodynamic stability and neutron transport criticality / von Chandramowli Subramanian

De
99 pages
High Performance Computingfor Stability ProblemsApplications to Hydrodynamic Stability andNeutron Transport CriticalityZur Erlangung des akademischen Grades einesDOKTOR DER NATURWISSENSCHAFTENvon der Fakultät für Mathematik derUniversität Karlsruhe (TH)genehmigteDISSERTATIONvonDipl.-Math. Chandramowli SubramanianausFilderstadtTag der mündlichen Prüfung: 17. Februar 2011Referent: Prof. Dr. Vincent HeuvelineKorreferent: Prof. Dr. Götz AlefeldAbstractIn this work we examine two kinds of applications in terms of stability and perform nu-merical evaluations and benchmarks on parallel platforms. We consider the applicabilityof pseudospectra in the field of hydrodynamic stability to obtain more information than atraditional linear stability analysis can provide. We present parallel computational tech-niques as well as evaluations of pseudospectra of complex flow processes. In the field ofnuclear engineering we treat the criticality problem which describes the crucial criteria forthe stability of a nuclear reactor. For the resulting eigenvalue problem we highlight theDavidson method as an attractive alternative to the so far widely used power method.ZusammenfassungIn dieser Arbeit untersuchen wir zwei Klassen von Anwendungen hinsichtlich ihres Stabi-litätsverhaltens und führen entsprechende numerische Berechnungen und Benchmarks aufparallelenPlattformendurch.
Voir plus Voir moins

HighforPStabiliterformanceyProblemsComputing

ApplicationstoHydrodynamicStabilityand
NeutronTransportCriticality

ZurErlangungdesakademischenGradeseines

DOKTORDERNATURWISSENSCHAFTEN

vonderFakultätfürMathematikder
(TH)KarlsruheersitätUnivgenehmigte

TIONATDISSER

onvDipl.-Math.ChandramowliSubramanian
ausFilderstadt

TagdermündlichenPrüfung:17.Februar2011

KReferenorreferent:t:Prof.Prof.Dr.Dr.GötzVincentAlefeldHeuveline

AbstractInthisworkweexaminetwokindsofapplicationsintermsofstabilityandperformnu-
mericalevaluationsandbenchmarksonparallelplatforms.Weconsidertheapplicability
ofpseudospectrainthefieldofhydrodynamicstabilitytoobtainmoreinformationthana
traditionallinearstabilityanalysiscanprovide.Wepresentparallelcomputationaltech-
niquesaswellasevaluationsofpseudospectraofcomplexflowprocesses.Inthefieldof
nuclearengineeringwetreatthecriticalityproblemwhichdescribesthecrucialcriteriafor
thestabilityofanuclearreactor.Fortheresultingeigenvalueproblemwehighlightthe
Davidsonmethodasanattractivealternativetothesofarwidelyusedpowermethod.

ZusammenfassungIndieserArbeituntersuchenwirzweiKlassenvonAnwendungenhinsichtlichihresStabi-
litätsverhaltensundführenentsprechendenumerischeBerechnungenundBenchmarksauf
parallelenPlattformendurch.ZumeinenwirdderEinsatzvonPseudospektrenzurUntersu-
chungderhydrodynamischeStabilitätdiskutiertumstärkereAussagenalsinderklassischen
linearenStabilitätstheoriezuerhalten.WirpräsentierensowohleffizienteparalleleLöser
alsauchnumerischeErgebnissevonPseudospektrenvonkomplexenStrömungsprozessen.
DesWeiterenumfasstdieseArbeitdasKritikalitätsproblem,welchesdasausschlaggebende
KriteriumfürdieStabilitätvonKernreaktorenbeschreibt.FürdasresultierendeEigen-
wertproblemzeigenwirdasDavidsonVerfahrenalsattraktiveAlternativederweitläufig
verwendetenPotenzmethodeauf.

tswledgmenknocAIamespeciallygratefultoProf.Dr.VincentHeuvelinewhomadethisworkpossibleat
all.Iwanttothankhimforhisenduringsupportandmanyfruitfuldiscussions.Also,I
ammuchobligedtoProf.Dr.GötzAlefeldwhoawokemyinterestsinnumericalanalysis
inthefirstplace.Iwouldliketothankhimforcarefulreadingandhiskindsupport.I
amindebtedtoPDDr.GudrunThäterforcarefulreadingandallherhelpfulhints.Many
thanksgoalsotoDr.SergeVanCriekingenfortheexcellentcollaborationwhichgaveme
theopportunitytoworkinthefieldofneutrontransporttheory.Finally,Iamgratefulto
mycolleaguesinthegroupsNumericalSimulation,OptimizationandHighPerformance
ComputingandInstituteforAppliedandNumericalMathematics1forallthediscussions
andtheexcellentworkingatmosphere.

Contsten

ductiontroIn11.1FunctionSpacesandNotations.........................

yStabilitdynamicHydro22.1ModelingFluidFlowProcesses.........................
2.1.1ConservationofMass...........................
2.1.2BalanceofMomentum..........................
2.1.3ConstitutiveEquations..........................
2.1.4TheIncompressibleNavier-StokesEquations..............
2.1.5HeatConduction.............................
2.1.6TheOberbeck-BoussinesqEquations..................
2.2LinearStabilityTheory.............................
2.2.1StabilityoftheZeroSolution......................
2.2.2LinearStabilityofaSteadyFlow....................
2.2.3TheBénardProblem...........................
2.2.4BibliographicalRemarks.........................
2.3Pseudospectra...................................
2.3.1PseudospectraofMatrices........................
2.3.2PseudospectraofLinearOperators...................
2.3.3BoundsonMatrixandOperatorExponentials.............
2.3.4PseudospectraofMatrixPencils....................

3NeutronTransportCriticality
3.1NeutronTransport................................
3.1.1GeneralAssumptions...........................
3.1.2CrossSectionDefinitions.........................
3.1.3TheLinearBoltzmannEquation....................
3.2Criticality.....................................
3.2.1TheαEigenvalue.............................
3.2.2ThekEigenvalue.............................

4DiscretizationoftheEigenvalueProblems
4.1GalerkinFiniteElementSpectralApproximation...............
4.1.1ProblemFormulationforEllipticEigenvalueProblems........
4.1.2SpectralApproximationofCompactOperators............
4.1.3APrioriErrorEstimatesfortheFiniteElementApproximation...
4.1.4BibliographicalRemarks.........................

12

56778991010111314151516811923

2525262729313113

333333354045

5

6

7

8

A

4.2DiscretizationoftheNeutronTransportEquation...............
4.2.1EnergyDiscretization..........................
4.2.2SecondOrderEvenParityFormulation.................
4.2.3AngularandSpatialDiscretization...................

EigenvalueSolversandParallelization
5.1TheDavidsonMethod..............................
5.2ComputationofPseudospectra.........................
5.3Parallelization...................................
5.3.1SparseMatrixVectorMultiplicationforFiniteElementMethods..
5.3.2ParallelComputationofPseudospectra.................

PseudospectrainHydrodynamicStability
6.1IncompressibleFlowoveraBackwardFacingStep...............
6.2NaturalConvectioninaHorizontalAnnulus..................

TheoryReactorNuclearinApplications7.1TheTakeda1Benchmark............................
7.2TheNEAC5G7Benchmark...........................

SummaryokOutloand

EstimatesandCalculationsforthePoincaréConstant
A.1BoundsforthePoincaréConstant........................
A.2EvaluationofthePoincaréConstant......................
A.3ApplicationinNaturalConvection.......................

yBibliograph

45454748

515154565656

599567

717275

79

81818458

87

1Chapter

troInduction

Thequestionofstabilityiscrucialinavarietyofdisciplines,suchasengineering,control
theory,physics,andmathematics.Thiscanbeseeninalargenumberofapplications:
designingbridgesorbuildings,controllingnuclearreactors,stabilizingflowprocessesin
pipeaircrafts.lines,Tdevoelopcopeingwithelectronicthesekindsstabilityofconprobletrolformsvisehicles,anoraenormousvoidingchallengeturbulencesduetoaroundthe
highcomplexityofthewholesolutionprocess.First,anadequatemathematicalmodel
describingdiscretizingthetecunderlyhniquesingablephtoysicalapproproximatecessestheneedmostodelbeaccurateconstilytuted.havetoThen,beestablappropriished.ate
Inordertosolvetheresultinglarge-scaledproblemsinareasonabletime,onefurtherneeds
toemployefficientnumericalmethodscapableofexploitingparallelplatforms.
Inthisworkwetreattwodifferentstabilityproblems.First,wedevoteourselvesto
thestabilityoffluidflowprocesses,whichisamajorfieldinthetheoryoffluiddynamics,
namelythehydrodynamicstabilitytheory.Inthisrespect,theexaminationisbasedon
exertingaperturbationonalaminarflow,andthentostudytheevolutionoftheperturbed
orquangrotitwsy.inBytime.meansWofefoeigencusvalueontheproblems,linearwestabilitydetermineanalysiswhetherwhichtheisperturdirectlybationlinkeddecatoys
theanalysis,spaectrumpurofelythelinearlinearizedstabilityoperatoranalysiscannotdescribingguarathenfloteew.theUnlikstabeilitayofnonlineartheunderlyingstability
process.Thisworkdiscussestheapplicabilityofpseudospectratotacklethismatter.
Pseudospectraarealsobasedonthelinearizedproblem.Despiteofthequiteexpensive
evaluationofpseudospectracomparedtospectra,webelievethatsuchaninvestigationis
avvaluable,ailable,asandinbythemeansadvenoftofpseudosphighpectra,erformanceonemaycomputinggainadeepmoreerinsighcomputattinionalthepowstabiliteryis
behaviorofphysicalsystems.
Basedonresultsalreadyavailabletoapproximateeigenvaluesofellipticoperators,we
derivapproacethehallowsusmathematicaltotreatfoundationcomplexflotowproapprocessesximateappropseudoximatedspbectrayfiniteaswell.elemenOurtmethogeneralds.
Asforthecomputation,whichconsistsofevaluatingsingularvalues,wefocusonthe
Davidsonmethodandproposeanefficientcomputationalschemeforparallelhardware
hitectures.arcThesecondaspecttreatedinthismanuscriptisthecriticalityproblemarisinginthe
fieldofneutrontransporttheory.Itisconcernedwiththestabilityofnuclearfissionchain-
criticalitreactions.yTheproblemtypicalmayalapplsobeicationformistheulatedmobydelinmeansgandofanconeigentrollingvalueofnprouclearblem.Thereactors.largestThe

1

2

ODUCTIONINTR

eigenvalueoftheresultingproblemindicateshowfarthefissionchain-reactionisfromthe
desiredstateofequilibrium,wheretheamountofneutronsemittedequalstheamountof
ed.absorbneutronsSofar,thepowermethodhasbeenthemethodofchoicetosolvethecriticalityeigenvalue
problem.Somemoreelaboratedmethods,suchastheArnoldimethodortheJacobi-
Davidsonmethod,havealsobeensuccessfullyappliedtothisproblem.However,tothe
bestsimpleroftheapproauthors’ximationknoscwledge,hemes(PtheseandmethoSdshamethovedsonlytobeenresolveusedtheintheangularcondeptextofendence,some
seeSection4.2).1N
formIntulationhiswtoorkstwateefothecusoncriticalithetykproblem.eigenvalueWeprapplyoblemthewhicDahisvidsononeofmethothedtowidelysolveusedthis
generalizedeigenvalueproblemintheframeworkofamoresophisticatedapproximation
scnativhemeeto(PNthemethoratherd).simpleOurnpowerumericalmethoresultsdbyshospweedingthatupthisourmethodiscalculationsapromisingsignificanalter-tly.
Furthermore,theDavidsonmethodshowsmoreflexibilityandrobustnessthanthepower
d.methoThisthesisisorganizedasfollows.InthefollowingSection1.1weestablishsomebasic
definitionsandnotationsneededthroughoutthismanuscript.InChapter2weexplain
theconceptsofhydrodynamicstability.Westartbyestablishingtheequationsdescrib-
ingstabilityincompressibletheoryandflothewandutilizationnaturalofconvpseudospectionecprotraincesses.thisrespAfterwect.ards,ThewesubfojcusectofonlinChap-ear
terform3isulatetoeigendescribvealuetheprlinearoblemsBoltzdescribingmanntheequationcriticalitmoy.delingChapterthe4neutronisconcernedtransport,withandtheto
theoreticalbackgroundofthediscretizationmethodswhichareusedtoapproximatethe
modelsarisinginthehydrodynamicstabilitytheoryandthecriticalityproblem.Inboth
thiscapplicationshapterweillustratemploesythethebasicDavidsonmechanismmethodofwhithechisappliedtreatedinparallelizatiChapterontec5.hnMoreoiques.vIner,
Chaptercompressible6weflowpresenovertabacpseudospkwardectrfaacingoftwstepoasdifferenwelltasfloawpronaturalcesses:convWectioneproconsidercessaninanin-
annulus.NumericalresultsforcriticalityproblemsareshowninChapter7.InChapter8
wediscussourresultsandgiveanoutlooktopotentialfutureresearchfields.Addition-
itsally,inqualityAppnendixumericallyAwe.derivMoreoeaver,bwoundeapforplythethePoincobtainedaréconstantresultsintoaanannnaturulusalconandvcectionheck
problem.

1.1FunctionSpacesandNotations

WegiveashortintroductiontoLebesgueandSobolevspaces,asweneedthesedefinitions
throughoutthiswork,especiallyforthefiniteelementframework.Thesedefinitionsare
ratherstandardandcanforinstancebefoundin[21,50,93].
SupposeΩtobeaLebesguemeasurablesubsetofRd(usuallyd=2ord=3)with
non-emptyinterior.LetfbearealorcomplexvaluedfunctiononΩwhichisLebesgue
Bymeasurable.Ωf(x)dx

1.1.FunctionSpacesandNotations

3

wedenotetheLebesgueintegraloffoverΩ.For1≤p≤∞wedefinetheLebesguenorm:
1/p
|f(x)|pdx,1≤p<∞,
fLp(Ω)=Ω
Ω∈xesssup|f(x)|=inf{C≥0:|f(x)|<Calmosteverywhere},p=∞.
ItcanbeshownthattheLebesguespacesLp(Ω)={f:fLp(Ω)<∞}equippedwiththe
norm∙Lp(Ω)areBanachspacesforany1≤p≤∞.Inthiscontext,weidentifytwo
functionsfandgiff(x)=2g(x)almosteverywhere,i.e.f−gLp(Ω)=0.Forthespecial
casep=2,wehavethatL(Ω)isaHilbertspacewiththeinnerproduct
(f,g)0:=(f,g)L2(Ω)=f(x)g(x)dx.(1.1)
ΩLetLl1oc={f∈L1(K):Kcompact,K⊂Ω˚}denotethesetoflocallyintegrable
functions.∞Here,Ω˚referstotheinteriorofΩ.Fαurthermore,C0∞(Ω)isthesetoffunctions
inC(Ω)withcompactsupportinΩ.LetDbeshorthandforthepartialderivative
operator:∂|α|
αD:=∂x1α1∂x2α2∙∙∙∂xnαd.
Here,then-tupleα=(α1,α2,...,αn)isamulti-indexwithαi∈N0and|α|=id=1αi.
WedefineaweakderivativeDwαfofafunctionf∈Ll1oc(Ω)ifwehavethatDwαf∈
Ll1oc(Ω)and
Dwαf(x)ϕ(x)dx=(−1)|α|f(x)Dαϕ(x)dx
ΩΩholdsforanyϕ∈C0∞(Ω).Sinceforanyfunctionϕ∈C|α|(Ω)theclassicalderivative
coincideswiththeweakderivative,wedonotmakeanydistinctioninthenotationofD
.DandwEquippedwiththedefinitionofweakderivatives,wecansetuptheframeworkfor
Sobolevspaces.Foranyk∈N0wesettheSobolevnorm
1/p
DαfLpp(Ω),1≤p<∞,
pαfWpk(Ω)=|α|≤k
|αmax|≤kDfL∞(Ω),p=∞.
TheSobolevspacesdefinedbyWpk(Ω)={f∈Ll1oc(Ω):fWpk(Ω)<∞}endowedwiththe
correspondingSobolevnormareBanachspaces.Forthecasep=2,Hk(Ω):=W2k(Ω)are
Hilbertspaceswiththeinnerproduct
(f,g)k:=(f,g)Hk(Ω)=(Dαf,Dαg)L2(Ω).
k|≤α|Notethatfork=0,wehaveH0(Ω)=L2(Ω)withthesameinnerproduct(∙,∙)0asin
(1.1).Moreover,wedefineH0k(Ω)tobethecompletionofC0∞(Ω)withrespecttothenorm
∙Hk(Ω).

4

ODUCTIONINTR

ForavectorxwesetxH=xT(thecomplexconjugateofthetransposedvector).The
samenotationisusedforamatrixA,i.e.AH=AT.Finally,fortwovectorsx,yinRd
orCd,wedenotetheusualscalarproductbyx∙y=xHy.

Chapter2

yStabilitdynamicHydro

Stabilityofafluidflowisacrucialquestioninmanyapplicationsaffectingourdailylife:
Howstableisaflyingaircraftexposedtogustingwinds?Canparticlesinpipelines,trans-
portingoilorgas,causecongestions?Howcanturbulencesinaliquid,transportingthe
heatoutofanuclearreactor,beavoidedinordertomaintainafailure-freecoolingsystem?
Thesekindsofquestionsareverycomplex.Foratheoreticalanalysis,anappropriate
modelassumptionsdescribinghavethetobeunderlyingmadeinorderphenomenatogetneaedsformtoubelationconstituted.(e.g.partialInmostdifferencases,tialmanequa-y
chotions)ose,bjusothtpurfeweaspandecntsofumericalthecomplexanalysisflocanwbcopehaeviorwith.whicButhisnoobservmatteredarewhichreflectmodeledbwye
availabletheoreticalresults–mostlyforverysimpleconfigurations.
Theapproachofhydrodynamicstabilityistoinvestigatehowalaminarfluidflowbe-
haveswithrespecttoperturbations.Iftheperturbationdecaysintimeandtheflowreturns
toitsoriginalstateitissaidtobestable.Ontheotherhand,iftheperturbationcauses
thefloturbulence,wtocbuthaitngemainytoaalsodiffertakeenttheflostate,wtoitaissaiddifferentotbelamiinnarstable.state.InstabilitInstabilitymayytmeansriggerin
aproparticular,cedureaincomputedindustrialprosolutionductionmighistnotjustbeimpobservossible.ableinexperimentsorthecontrolof
stabilityInthetheoryfieldoisfhbasedydrodyonnamicexaminingstabilittheytherekineticaretenergywoofmainthefloapproacwbyhes.meansTheofinnonlinetegralar
inequalitytechniques.Thismethodisalsoreferredtoasenergymethod.Thelinearstability
theoryisconcernedwithalinearizedmodelofthefluidflowandestablishesstatementsby
meansofthespectrumofthelinearizedoperator.Ifalleigenvalueslieinthelefthalfof
inthethericomplexghthalfplane,ofthetheflocomplexwissaidplane,tobtheeflolinearwisstablinearle.Ifinstablethereis(andatleastthereforeoneeiinstablegenvalueas
see).willewprocInesses,Sectionwhere2.1wweegivconsidereatwbriefodifferenreviewtoftypthees.govFirst,erningwestudyequationsaninmodelingcompressiblefluidflofluidw
floderwtomodeledconsiderbyatheheatNaivdrivener-Stokfloesw.Thiequations.sphenomeAfterwnonards,iswformeextendulatedbthisymomeansdelofintheor-
Oberbeck-Boussinesqequations.SubjectofthesucceedingSection2.2istoestablishthe
basicdefinitionsoflinearstabilitytheoryintermsoftheNavier-Stokesequations.More-
over,relationwebpetoinwteenoutthestabilityandifferencesdeigebnetvwalueeenlinproblearems.andAfterwnonlinearards,stabilinitySectionands2.3howwethegiveclosean
outlinetoamoregeneralconceptthaneigenvalues,namelypseudospectra.Wedescribe

5

6

STYNAMICODHYDRABILITY

howpseudospectracanbeanassettogaindeeperinsightsinthestabilitybehaviorof
systems.ysicalph

2.1ModelingFluidFlowProcesses
Westartinthissectionbyderivingthebasicequationswhichdescribeflowprocessesof
incompressibleliquidsorgases,namelytheincompressibleNavier-Stokesequations,see
e.g.[7,28,82,114].Afterwards,wearealsoconcernedwithflowprocessesgovernedby
temperaturedifferencesandgravityforces,theso-calledOberbeck-Boussinesqequations.
Bothmodelsareconstitutedbytheconservationlawsofmassandmomentum,whereas
theOberbeck-Boussinesqequationsadditionallyrequirethebalanceofenergytoformulate
conduction.heattheLetΩ∈Rd(d=2ord=3)bearegionfilledwithafluidmovingduetoprevailing
internalandexternalforces.Inwhatfollows,weconsideraboundedtestvolumeV(t)∈Ω
whichsatisfiestherequirementstoapplythedivergencetheorem.Furthermore,wedenote
d-dimensionalvectorsbybold-facedletters,e.g.x∈ΩdenotesthespatialvariableinRd.
Foreachtimet,weassumethatthefluidhasawell-definedandcontinuousmassdensity
ρ(x,t)>0suchthatthemassm(t)ofthevolumeV(t)isgivenby
m(V(t))=V(t)ρ(x,t)dx.(2.1)
Thisconditionisalsoreferredtoascontinuumhypothesis.Thecontinuumhypothesis
allowsustospeakofphysicalpropertiesofaninfinitesimallysmallvolumeelementofthe
fluid,whichwecallmaterialparticleinwhatfollows.Wedefinethemovementofamaterial
particleη∈V(t)byafunctionx(η,t)suchthatthematerialparticleηattimetislocated
atpositionx(η,t)∈Ω.Weassumethefunctionxtobeinvertibleandfurthermoretobe
oth.smotlysufficienAflowmaybedescribedbythephysicalpropertiesofeachmaterialparticleasafunc-
tionoftime,whichiscalledLagrangeanformulation,i.e.foranyfixedmaterialparticleη
onefollowsitstrajectory.WhereasintheEulerianformulationtheflowisdescribedby
specifyingthephysicalpropertiesasafunctionoftimeforanyfixedpointofthedomain.
IntheframeworkoftheEulerianformulation,wedefinethevelocityofamaterialparticle
whichisatpositionx(η,t)bythed-dimensionalvectorfield
v(x,t)=∂tx(η,t).
Letf(x,t)representanyscalarphysicalpropertyofthefluidatpositionxandtimet.An
observationofthetemporalchangeoffatafixedpointinachosencoordinatesystem
isdescribedbythelocalderivative∂tf(x,t).Observingthetemporalvariationofffora
fixedparticleηbyf(x,t)=f(x(η,t),t)leadstothetimederivativeintheLagrangean
formulation,theso-calledmaterialderivative:
df=df(x(η,t),t)=v∙f+∂tf,
dtdtwherethenablaoperatoractsonlyonthespatialvariablex.Equippedwiththese
definitions,weareabletostatetheTransporttheorem.Therefore,weassumethatf,and
otherfunctionstobeintroducedlater,aresmoothenoughtoapplystandardoperations

2.1.ModelingFluidFlowProcesses

7

hem.tonTheorem2.1(Transporttheorem)Foranysmoothscalarfieldf(x,t)andanytest
volumeV(t)wehave
ddtV(t)f(x,t)dx=V(t)∂tf(x,t)+∙(f(x,t)v(x,t))dx.
Foraproof,seee.g.[28,114].
ThederivationoftheNavier-Stokesequationsisbasedonthefollowingconservation
laws,whicharediscussedinthefollowingsections.
(1)Conservationofmass:Massisneithercreatednordestroyed.
(2)Conservationofmomentum(Newton’ssecondlaw):Themomentumchangerateof
amaterialvolumeisequaltotheforceexerted.

MassofationConserv2.1.1Themassconservationlawstatesthatthemassm(V(t))ofanarbitrarytestvolumeis
constantintime,i.e.dd
dtm(V(t))=dtV(t)ρ(x,t)dx=0.
ApplyingtheTransporttheorem2.1,wededucethat
V(t)∂tρ(x,t)+∙(ρ(x,t)v(x,t))dx=0
holdsforanyV(t).Assumingtheintegrandtobecontinuous,wehavethatateachpoint
∂tρ+∙(ρv)=0,(2.2)
whichisknownasthecontinuityequation.

tumMomenofBalance2.1.2Foranycontinuum,therearetwotypesofforcesactingonapieceofmaterialV(t).Body
(orexternal)forces,suchasgravityorelectromagneticforces,canberegardedasacting
throughoutavolume.Thesearequantifiedby
Fext(V(t))=ρ(x,t)f(x,t)dx,
)t(Vwiththevectorvaluedfunctionf(x,t)representinganexternalforcedensity.Whereas
contact(orinternal)forcesactthroughthesurfaceofV(t)andaredescribedbythestress
vectors(x,t)actingperunitarea:
Fint(V(t))=∂V(t)s(x,t)ds.

8

STYNAMICODHYDRABILITY

unitHere,nor∂V(mal.t)Itreferscantobetheshobwnoundarythats(ofx,tthe)=testn∙vσ,olumewhereVσ(t)=.σijLetisnthedenotestressthetensorout,wardsee
e.g.[7,82].
dThelinearmomentumconservationlawstatesthat
dtV(t)ρ(x,t)v(x,t)dx=V(t)ρ(x,t)f(x,t)dx+∂V(t)s(x,t)ds(2.3)
holdsforanytestvolumeV(t).Applyingthetransporttheoremandthedivergencetheo-
rem,and,moreover,usingthecontinuityequation(2.2)yields
ρ∂tv+ρ(v∙)v=ρf+∙σ.(2.4)
Thisresultingequation(2.4)iscallednon-conservativeformofthemomentumequation.
tensorUsingisthesymmetric,principlei.e.ofσTconserv=σ,ationseeofe.g.[7,angular82].momentum,onecanshowthatthestress

ationsEqueConstitutiv2.1.3Sofar,thederivedequationsaregenericandappropriateforallkindsoffluids.However,
thenumberofunknowns
•velocityvectorv∈Rd,
•densityρ∈R,
•stresstensorσ∈Rd×d
doesnotmatchthenumberofconditionsimposedbythefollowingequations
•continuityequation:∂tρ+∙(ρv)=0,
•momentumequation:ρ∂tv+ρ(v∙)v=ρf+∙σ,
•conservationofangularmomentum:σ=σT.
Inordertocompletethesystem,weneedtorelatethestresstensortothevelocityand
density.Theserelationsarederivedempiricallyinphysicalexperiments.Thefirstassump-
tionwemake,isthatthefluidisStokesian,i.e.thestresstensorissphericallysymmetric
whenthefluidisatrest:
σ|v=0=−pI,
wherep=p(x,t)denotesthethermodynamicpressure.Thismeans,thattheshearstress
componentsofσvanish,whilethenormalstressesareequaltothepressure.Apparently,
forafluidinmotionthisisnotthecase.Here,weset
σ|v=0=−pI+τ,(2.5)
definewhereτtheisadeformationsymmetric(orrtensorateofstrdescribingain)tensortheshasearD=1stresses.v+To(getv)Ta.Wrelationefassumeorτthatwe
2τ=F(D),whereFisanappropriatecontinuousfunctionmappingsymmetrictensorson

2.1.ModelingFluidFlowProcesses

9

symmetrictensors.IfwefurthermoreassumeFtobelinear(i.e.weconsiderNewtonian
fluids)andisotropic,equation(2.5)yields
σ=(−p+λ∙v)I+2µD,
see[7,28,82].Here,λisthevolumeviscosityandµtheshearviscosity.

2.1.4TheIncompressibleNavier-StokesEquations
Incompressibilitymeansthatthedensityisconstantinspaceandtime.Thus,byvirtueof
thecontinuityequation(2.2)thevelocityissolenoidal,i.e.∙v=0.Thisimplies
∙σ=−p+(λ∙v)+2µ∙D
=−p+µ(Δv+(∙v))
=−p+µΔv.
obtainSummingtheup,incomprwith(essibl2.4)eandNavier(2.2-S),tokesandedefininquationsg:thekinematicviscosityasν=µ/ρ,we
∂tv−νΔv+(v∙)v+1p=f,
ρ(2.6)∙v=0,
wheresuitableinitialandboundaryconditionsonvhavetobeimposed.

ConductionHeat2.1.5Uptonow,wehavestudiedonlyisothermalflowprocesses,i.e.weassumedaconstanttem-
perature,whereasinthissectionweintroduceflowprocesseswithheatvariation.Typical
applicationsincludecoolingsystemsandinsulationtechnologies.
Thermalconvections,wheretheconvectionisonlydrivenbysomeexternalforces(e.g.by
afan),andtheeffectsofbuoyancyareneglected,areknownasforcedconvections.How-
ever,nonewflowphenomenaarise,sincethemotionofthefluidisunaffectedbythe
temperatureandcanbedeterminedinthewayasintheprecedingSection2.1.4.Onthe
otherhand,thecasewheretheconvectionisonlygovernedbybuoyancyforcesiscalled
naturalconvection.Inthiscasethefluidwouldbeatrestifnotemperaturevariation
wasintroduced.Inmanyapplicationsthemotionofthefluidisinducedbytemperature
differencebetweenboundaries.Mixedconvectionisthecombinationofforcedandnatural
convection,i.e.buoyancyforcesandexternalforcesarebothconsidered.
Inthiswork,weconfineourselvestonaturalconvection,wherethefluidismovingdue
totemperaturedifferencesandgravity.Thisphenomenonisformulatedbymeansofthe
Oberbeck-Boussinesqequations(orBoussinesqapproximation).Theseequationsarebased
onthecontinuityequation,themomentumequation,andanadditionalequationforthe
conduction.heatSofar,wehavededucedtheequationsofcontinuityandmotionbyusingthelawsof
conservationofmassandmomentumrespectively.Toformulatetheheatconductionina
fluid,weneedtoexploittheconservationofenergyaswellwhichisexpressedbythefirst

10

STYNAMICODHYDRABILITY

lawofthermodynamics.Theequationofheatconductionreads
ρ∂t(cVT)+ρv∙(cVT)=∙(kT)−p∙v+Φ,(2.7)
see[26,105,114].Here,T=T(x,t)denotesthetemperature.Furthermore,cVisthe
specificheatatconstantvolume,kisthethermalconductivity.TheenergydissipationΦ
perunitvolumebyviscosityreadsinthecaseofNewtonianfluids
Φ=2µD:D−32(∙v)2,
wherethescalarproductoftwotensorsisdefinedbyA:B=AijBij.
ij2.1.6TheOberbeck-BoussinesqEquations
Boussinesq[18]andOberbeck[73]recognizedthatifthetemperaturevariationsaresmall,
thedynamicsofafluidcanbeapproximatedbyassumingaconstantdensityeverywhere,
exceptinthebuoyancy.Togiveabriefsketchofthisargumentation,supposethatthe
densitycanbequantifiedbythelineardependence
ρ=ρ0[1+α(T0−T)],
whereT0denotesthetemperatureforwhichρ=ρ0andαisthevolumetricexpansion
coefficient.Usually,αisintherangeof10−3to10−4K−1(e.g.perfectgas:α≈3∙10−3K−1;
liquidswhicharemostlyusedinexperiments:α≈5∙10−4K−1).Ifwehavesmall
temperaturedifferences,sayupto10K,thevariationsinthedensityareatmostoneper
cent.Thatiswhythedensityvariationsmaybeignoredinalltermsof(2.2),(2.4)and
(2.7),withoneexception:thetermdescribingthegravityforce,see[26,38,105].This
externalforceisgivenbyfg=ρg,wheregistheaccelerationduetogravity.
Next,wehavetosimplifyequation(2.7).Sincewehavereplacedthecontinuityequation
by∙v=0,wecanignorethetermp∙v.Furthermore,wecantreatcVandkas
constants.Finally,weneglecttheviscousdissipationΦsinceitismostcommonlyofvery
loworder.Summingup,theOberbeck-Boussinesqequationsread
∂tv−νΔv+(v∙)v+1p=g[1+α(T0−T)],
ρ∙v=0,(2.8)
∂tT+v∙T−κΔT=0,
wherewehavereplacedρ0byρandκ=k/(ρcV)denotesthethermaldiffusity.Forboth,
thevelocityvandthetemperaturedistributionT,suitableinitialandboundaryconditions
havetobeprescribed.Amoredetailedphysicalbackgrounddescribingtheapproximations
weusedinthissectioncanbefoundin[105].

TheoryyStabilitLinear2.2InlemsthisforsethectionweincompressibleshowtheNaclosevier-Stokrelationesbetwequations.eenstabilitFurtheyanalysisrmore,weandbrieigeneflyvaluediscussprob-the

TheoryyStabilitLinear2.2.

11

stabilityoftheBénardproblem.Laterinthiswork,inSection2.3.3,thelinearstability
theoryisputinamoreabstractframeworkwhichallowsustoestablishboundsbymeans
ofspectraandpseudospectra.
Forsimplicity,wescalethepressurepbysettingρ=1.Then,theviscousfluidflowwe
considerisgovernedby

∂tv−νΔv+(v∙)v+p=f,
∙v=0(2.9)
inaboundeddomainΩ,withboundaryconditionsv=hon∂Ωandaninitialvalue
v(x,0)=v0(x).
Beforewegodeeperintothelinearstabilitytheory,westartwithanexampletoshow
thedifferentapproachesoflinearandthenonlinearstability.Inwhatfollows,weassume
vandptobesufficientlysmoothandthedomainΩunderconsiderationtofulfillthe
requirementsofthedivergencetheorem.

2.2.1StabilityoftheZeroSolution
Inthisexampleweexaminethestabilityofthezerosolutionof(2.9),wheretheexternal
v0force(x)fforisxgiv∈enΩ,byandapzerootenbtialfoundary=v−Φalues,.i.e.Moreov(vx,er,t)w=e0impforoseanyanx∈∂initialΩ.vPhalueysicav(xlly,,0)this=
setupcorrespondstoamodelofacontainerfilledwithaliquidwhichisarbitrarilymoved
toaroundthegrafort<vitational0,andforcethen−helΦ.dfixedEmpiricallyatt,=w0e.expForectt≥the0ptheconerturbationtainertoisdecaonlyyexpandosedthe
liquidtoreturntoitsinitialstate.

yStabilitNonlinearThenonlinearstabilityconsidersthekineticenergyE(t)whichisgivenby
E(t)=1v02=1(v,v)0.
22Here,(∙,∙)0istheinnerproductinL2(Ω)and∙0thecorrespondingL2(Ω)norm,cf.Sec-
tion1.1.DifferentiatingE(t)withrespecttotyields
dtdE(t)=(∂tv,v)0.
Next,thefirstequationof(2.9)ismultipliedbyvandwetaketheintegraloverΩ.Then
weintegratebypartstoobtain
dE(t)=−νv02+1(v∙v,∙v)0+(p+Φ,∙v)0,
2dtsee[94].Byassumptionvisdivergencefreewhichimplies
dtdE(t)=−νv02.

12

ABILITYSTYNAMICODHYDR

FromthePoincaréinequality(A.1)wededuce
dtdE(t)≤−νkp2v02=−2νkp2E(t),
henceE(t)≤e−2νkp2tE(0).
Sincekp>0,wehavethatv0→0ast→∞,whichmeansthatallinitialperturbations
dieaway.NotethatthekinematicviscosityνaswellasthePoincaréconstantkpdetermine
theevolutionofthedisturbance.

yStabilitLinearThelinearstabilitytheoryfollowsadifferentapproach.Weassumetheperturbationtobe
smallinordertoneglectproductsofhigherordercontainingv.Thismeans,weconsider
problemlinearthe∂tv−νΔv+p=f,
0=v∙inΩ,withthesameinitialandboundaryconditionsasbefore.Supposethatvandpare
comprisedofsuperpositionsofnormalmodes
v(x,t)=v˜(x)eλt,p(x,t)=p˜(x)eλt,
wheretheeigenvaluesλandthecorrespondingeigenfunctions(v˜,p˜)satisfy
λv˜=νΔv˜−p,˜(2.10)
0=∙v˜
inΩ,withzeroDirichletboundaryconditionsforv˜.Clearly,anymodewithReλ<0
decaysintime,whereasamodewithReλ>0growsexponentially.Supposewehave
countablymanyeigenvalueswithnoaccumulationpointat0.Thiscanbeshowninthe
frameworkofellipticoperatorsassumingaboundeddomainΩwithasufficientlysmooth
boundary∂Ω,seeSection4.1.Wesaythatthezerosolutionislinearstableifallnormal
modesdecay,i.e.Reλ<0foralleigenvalues.Ontheotherhand,ifthereexistsone
eigenvalueλwithReλ>0,thezerosolutionissaidtobelinearinstable,see[37,38].
Inourcase,itiseasilyshownthatReλ<0holdsforallmodes.Wemultiplythefirst
equationof(2.10)byv˜andintegrateoverΩ.Thenintegratingbypartsunderconsideration
of∙v˜=0leadsto
λv˜02=−νv˜02.
AsaconsequenceofthePoincaréinequality(A.1)weobtain
λ≤−νkp2<0
foralleigenvaluesλ.Notethat,asinthecaseofnonlinearstability,thePoincaréconstant
andthekinematicviscosityareessentialforthestability.Bothquantitiesprescribethe
evolutionofthenormalmodes.
Aswehaveseenexemplarilyforthestabilityofthezerosolution,bothapproaches,
linearandnonlinearstability,arecloselyrelatedtoeigenvalueproblems.Theenergy

TheoryyStabilitLinear2.2.

13

methodemploysthePoincaréconstant,whichcanbedeterminedbymeansofaneigenvalue
theproblemspectrum(forbofoundedtheopderatoromains),inseeaddition.AppendixTheA,nonlinearwhereasthestabilitlineyartheorystabilitycantheoryguaranusestee
ifthstabiliteyupp,sierncebitoundprotendsvidesantouppinfiniteryb,nooundforstatementhetevinolutiontermsofofthestabilitpyerturbation.canbeHomade.wever,On
theeigenvotheraluehalyingnd,inthethelinearrighthalfstabilitofythetheorycomplexcanplane.assureButinstabilitlinearyasstabilitsoonyasdoeswenohatveimplyone
general.inystabilitnonlinear

2.2.2LinearStabilityofaSteadyFlow
Toinvestigatethestabilityofanarbitrarysteadyflow(V(x),P(x))governedby(2.9),we
disturbthesteadysolutionatt=0byafunction(v0,p0).Theselectedbasicflow(V,P)
isassumedtobeknown,eitheranalyticallyoronlybymeansofnumericalcomputations.In
ordertostudytheevolutionoftheperturbation(v(x,t),p(x,t)),weinserttheperturbed
quantities(V+v,P+p)in(2.9)yielding
∂tv−νΔv+(v∙)V+(V∙)v+(v∙)v+p=0,
∙v=0(2.11)
inΩ,withzeroboundaryconditionsv=0on∂Ωandinitialconditionsv(x,0)=v0,
p(x,0)=p0.
Wefollowthestabilitytheoryofordinarydifferentialequationsbydefiningstabilityin
thesenseofLyapunov,see[4,37,51].Therefore,weassumethefunctionsvandpunder
considerationtobeelementsofaBanachspace,where,forbrevity,wedenoteallnormsby
∙.Abasicflow(V,P)issaidtobestableifforallε>0thereexistsaδ=δ(ε)such
thatv(x,0)<δandp(x,0)<δ
impliesv(x,t)<εandp(x,t)<εforallt>0.
Otherwiseitissaidtobeinstable.Thismeans,astableflowensuresthatallinitially
smallperturbationsremainsmallforalltime.If,additionally,theperturbationsdecay
asymptotically,wesaythatthebasicflowisasymptoticallystable,i.e.itisstableand
v(x,t)→0andp(x,t)→0ast→∞.
Thelinearstabilitytheoryassumestheperturbationstobesmall.Sowemayneglect
productsoftheperturbedquantitiesin(2.11)yieldingthelinearproblem
∂tv−νΔv+(v∙)V+(V∙)v+p=0,(2.12)
∙v=0,
withthesameinitialandboundaryconditions.Sincethebasicflowissteady,thecoeffi-
cientsin(2.12)areindependentoft.Supposethatasolutionoftheinitialvalueproblem
(2.12)maybeseparated,i.e.itisalinearsuperpositionofnormalmodes
v(x,t)=v˜(x)eλt,p(x,t)=p˜(x)eλt,

14

ABILITYSTYNAMICODHYDR

wheretheeigenvaluesλandcorrespondingeigenfunctions(v˜,p˜)satisfy
λv˜=νΔv˜−(v˜∙)V−(V∙)v˜−p,˜(2.13)
0=∙v˜
inΩandv˜=0on∂Ω.Alleigenvaluesλof(2.13)areeitherrealoroccurincomplex
conjugatepairs.Again,ifReλ<0,thecorrespondingmodediesoutintime.Whereasa
modewithReλ>0resultsclearlyininstability.Finally,amodewithReλ=0iscalled
neutrallystableandmaytriggernonlinearinstability.
Ifwehavecountablymanyeigenvalueswithnoaccumulationpointat0,weconclude
thatthebasicflowofthelinearizedproblemisstableifallnormalmodesarestable,
i.e.Reλ<0foralleigenvalues.IfReλ>0foratleastoneeigenvalue,thebasicflowis
instable,see[37,38].Notethatlinearinstabilitymeansinparticularnonlinearinstability.
Physically,weexpectinstabilitiessuchasturbulencestooccurforfastmovinginviscid
fluids.Thesecharacteristicsareinherentinadimensionlessquantity,namelytheReynolds
number.ThisnumberisdefinedasRe=VL/νwithcharacteristicvelocityVandchar-
acteristiclengthL.Typically,aflowbecomesinstableastheReynoldspassesacertain
threshold,whichiscalledcriticalReynoldsnumber.
Hence,thecriticalReynoldsnumberRecisdefinedasthesmallestnumbersuchthat
thebasicflowunderconsiderationisstableforallRe≤Rec,andbecomesinstablefora
Re>Rec.Intermsoflinearstabilitytheorythismeansthatalleigenvaluesof(2.13)have
nonnegativerealpartforRe≤Rec,andthereisatleastoneeigenvaluewithpositivereal
partforaRe>Rec.

ProblemBénardThe2.2.3Weconsideralayerofafluidconfinedbetweentwoparallelplaneswhichisheatedfrom
below.Ifthetemperaturegradientissufficientlylargetoovercomethegravitationalforce,
atessellatedpatternofcellularmotionmaybeobserved.ThisphenomenoniscalledBénard
(orRayleigh-Bénard)convection,see[15,26,94].
Letthecoordinatesofthespatialvariablexbedenotedbyx=(x,y,z)inthethree-
dimensionalcaseandbyx=(x,z)inthetwo-dimensionalcase.Inordertodescribethe
naturalconvectionprocess,supposethatthefluidisinaninfinitelayerz∈(0,l),and
thatwehavefixedtemperaturesT0atz=0andTlatz=lwithT0>Tl.Forthe
Oberbeck-Boussinesqequations(2.8),thereexistsasolutionatrestwithalineartemper-
distributionaturev=0,T=−βz+T0,
wherethetemperaturegradientβisgivenby
β=T0−Tl.
dThestabilityofthissystemisdependentonthedimensionlessRayleighnumberdefinedby
Ra=αgβl4,
κνwheregisgivenbythegravitationalvectorg=(0,0,−g)T.
Asforthelinearstabilityanalysis,thereexistsasmallestnumberRaclinsuchthatfor

ectraPseudosp2.3.

15

anyRa>Ralcinthebasicstateatrestislinearinstableandhenceinstable.Ingeneral,this
doesnotimplythatthesolutionisstableforallRa<Ralin.Consideringthenonlinear
stability,thereexistsacriticalRayleighnumberRacnlsuchcthatthesolutionisstablefor
anyRa<Racnl.ForthesolutionatrestonecanshowthatRaclin=Racnl,seee.g.[94].In
otherwords,thissolutionislinearstableifandonlyifitisnonlinearstable,andthatitis
linearinstableifandonlyifitisnonlinearinstable.
AstheRayleighnumberexceedsthiscriticalvalue,adifferentsteadysolution–namely
aBénardconvection–canberealized.ThecriticalRayleighnumberforthisprimary
bifurinstancecationbefodepundendsin[on37,the38].bByoundaryfurthersetupincreasingchosenthe(rigidRaoryleighfreenumbber,oundaries)morebifandcaurcationsnfor
inlaboratoryexperimentscanberealized,wherealsotime-dependentbasicstatesare
observed.Unlikefortheprimarybifurcation,thesedependonthePrandtlnumber
ν=rPκ

aswell,seee.g.[37].
Inthesetupwestudylater(cf.Section6.2andAppendixA)noexactsolutionofthe
Oberbeck-Boussinesqequations(2.8)isknown.Inparticular,thereexistsnosolutionat
restasfortheBénardproblemdescribedabove.Therefore,undercertainconditions,a
simplifiedmodelisusedforwhichanalyticalsolutionsareknown.InAppendixAwe
evaluatetherelativeerroremergingfromthissimplification.Inthisframework,weshow
thesignificanceofthePoincaréconstantforthesetypeofevaluations.
ThePoincaréconstantisanimportanttoolforthequalitativeandquantitativedescrip-
tionoffluiddynamicsingeneral.Especially,itisessentialforthestabilitybehavioraswe
haveseenexemplarilyforthestabilityofthezerosolutionoftheNavier-Stokesequations.
InAppendixAwederiveaboundforthePoincaréconstantinaspecialsetupandperform
numericalcomputationswhichshowthattheestablishedboundisalmostsharp.

RemarksBibliographical2.2.4Therequotedaresofarmanwyemenreferencestiononhereh[26ydro,37dy]asnamicbooksstabilitfoycusedoanalysis.ntheInlinearadditionstabilittoythetheory.referencesThe
energymethodistreatedin[94]withanemphasisonthermalconvection.Bothlinearand
innonlineartermsofstathebilityNaavrecoier-Stokveredesin[equa38,tions90].Acanrelatibeonfoundbetwinee[n89].linearStabiandlitynonlineaofrsolutionsstabilitofy
ordinarysolutionsindifferenBanactialhspacesequationswearerefertodescrib[31].edin[4,8].Forageneraltreatiseofstabilityof

ectraPseudosp2.3Assystem.wehavHoewseeven,er,eigeninvexpalueserimenareats,vaeryfluidimpfloortanwisttoolobservtoedstudytobtheecomestabilitturbulyofenatalthougdynamicalh
anmosteigenvnotablyaluetheanalysiscaseforindicatesstronglylinearnonstabilinormalty,seeproblems.[104]Oneandremedyreferencesistointherevin.estigateThistheis
nonlinearstabilitywhichcanassurethestabilityofafluidmotion.However,tohandle
nonlinearproblemsmoreelaboratedtechniquesthanforlinearproblemsarenecessary.
Sincewepreferalineartechnique,weconsideranotherencouragingapproachbymeansof

16

ODHYDRABILITYSTYNAMIC

pseudospectra,whichareamoregeneraltoolthanspectra.
Weconsideradynamicalsystemdescribingtheevolutionofaperturbationuby
dAu,=udtwhereweassumeAtobealinearoperator.UndercertaintAconditions,thesolutioncan
beinitialexpresseddisturbanceinformsu(0).ofIfathematrixsporectrumopoferatorAliesexpinonenthelefttialehalfuof0,thewherecomplexu0represenplane,tsetAtheu
0tendstozerofort→∞andnoinstabilityisdetected.Consideringthepseudospectraof
AonemaydiscoverthatetAu0becomesarbitrarilylargeforfinitetandhencemaytrigger
instabilities.InthissectionweestablishlowerandupperboundsonetAbymeansofpseudospectra.
Westartbyintroducingsomebasicpropertiesofpseudospectraofmatricesandlinear
operators.Ourdiscussionfollows[103].

MatricesofectraPseudosp2.3.1Webeginwiththreeequivalentdefinitionsforpseudospectraofmatrices.Inwhatfollows
wewrite(z−A)instead−1of(zI−A),whereIistheidentityonCn.Moroever,weusethe
convention(z−A)=∞foranyz∈σ(A),whereσ(A)denotesthespectrumofA.
Definition2.2LetA∈Cn×nandε>0.Thesetsσε(A)definedby
σε(A)={z∈C:(z−A)−1>ε−1},(2.14)
σε(A)={z∈C:z∈σ(A+E)forsomeE∈Cn×nwithE<ε},(2.15)
σε(A)={z∈C:(z−A)v<εforsomev∈Cnwithv=1}(2.16)
arecalledε-pseudospectrumofA.
Noteinequalitiesthatathesthesedefinitionsleadtodepequivendalenontthedefininorm.tionsWforehaclosedvecophoseneratorsindefinitionsBanachusingspacesstrictin
Section2.3.2,see[25].
Theorem2.3Theabovedefinitions(2.14),(2.15)and(2.16)areequivalent.
ProTof.oFproorvez(∈2.15σ()A)⇒(the2.16),equivlet(alenceA+isE)obv=vious.zvSoforwesomeassumeE∈zC∈/nσ×(nA).withE<εand
v=1.Thenwehave(z−A)v=Ev≤E<ε,hence(2.16).
Toprove(−2.161)⇒(2.14),let−1(z−A)1v=su1withu=v=1ands<ε.Thus,we
haveTo(zpro−veA)(2.14)≥⇒((z2.15−),A)suppuose=s(vz−=As)−>1ε,>whic1.hmeansThen,(2.14there).exists<εand
ε11−u,vNext,withweushow=thatv=there1existssatisfyingE∈(zC−n×An)withu=Esv.=sHence,andwEevha=vesuzvwhic−hAv=impliessu.
(Amatrix+E)ofvthe=zvform−Esu=+ssuuw=HzvwithandvHwhence=1.(By2.15).virtueofThereforthee,wecHahn-BanachoosehEastheoremarank-1-there
isalinearfunctionalLonCnwithLvn=v=1andHL=1.DuetotheHrepresentation
theoremofRiesz,thereexistsaw∈CwithLv=vw.Thisimpliesvw=Lv=1.
FromL=1weobtainmaxx=1xHw=1.Thus,E=smaxx=1u(wHx)=
smaxx=1wHx=s,whichcompletestheproof.✷

Pseudosp2.3.ectra

17

Figure2.1:Pseudospectraforlogε∈{0,−0.5,−1,...−4}withrespecttothespectral
normofthediscretizedconvectiondiffusionoperator−0.05Δu+0.3ux+0.7uydefinedon
aunitsquarewithzeroDirichletboundaryconditions.

Thenexttheorempresentssomebasicpropertiesofpseudospectra.Inthisrespectwe
definethesumoftwosetsC,DbyC+D={c+d:c∈C,d∈D}.
Theorem2.4LetA∈Cn×nandε>0.
(1)σε(A)isnonempty,open,andbounded.
(2)If0<ε1≤ε2,thenσε1(A)⊆σε2(A).
(3)ε>0σε(A)=σ(A).
(4)σε(A)⊇σ(A)+B(ε),whereB(ε)={z∈C:|z|<ε}denotestheε-ballaroundthe
origin.(5)σε(A)consistsofatmostnconnectedcomponents,eachcontainingatleastoneeigen-
.AofvalueIfwechoosethespectralnorm,i.e.∙=∙2,wecanderivefurtherproperties.First,
notethatthespectralnormofamatrixAcorrespondstoitslargestsingularvalue,i.e.the
largesteigenvalueofAHA,seee.g.[98].ThisimpliesthespectralnormoftheinverseA−1
tobetheinverseofthesmallestsingularvalue.Thus,bydenotingthesmallestsingular
valueof(z−A)bysmin(z−A)wehave
1(z−A)−12=smin(z−A).
Thisleadstoafurtherrepresentationofpseudospectra,namely,
σε(A)={z∈C:smin(z−A)<ε}.(2.17)

Theorem2.5LetA∈Cn×nand∙=∙2.Thenforallε>0,
(1)σε(AH)=σε(A),

(2.17)

18

ODHYDRABILITYSTYNAMIC

(2)σε(A)=σ(A)+B(ε)ifandonlyifAisnormal.
Prthatoof.AisThenormal.firstThen,assertionbyfollothewsspectraldirectlyftheorem,rom(2.17there).Fexistsoratheseunitarycondmatrixassertion,Usuchsuppthatose
HofAA=.UΛSinceU,forwhereanyΛunitarydenotesmatrixtheUdiagonwealhavematrix(z−withUAUelHemen)−1ts2=equalUt(oz−theA)−eigen1UvH2alues=
(z−A)−12,weconclude
(z−A)−12=(z−UAUH)−12=(z−Λ)−12
11=minλi|(z−λi)|=dist(z,σ(A)),
andhenceσε(A)=σ(A)+B(ε).Foraproofoftheconversestatement,wereferto[103].
✷Feigenurthervaluespropandertiesofpseudosppseudospectraofectramatricescanbisefounddiscussedin[in103[4].0].Thecloserelationbetween

2.3.2PseudospectraofLinearOperators
Thedefinitionsofthelastsectioncanbegeneralizedtolinearoperatorsactingonan
arbitraryBanachspaceX.WeconsiderclosedlinearoperatorsA:D(A)→X,where
D(A)⊆XdenotesthedomainofA.ThesetofclosedlinearoperatorsmappingfromX
toXisdenotedbyC(X).Thus,bydefinition,A∈C(X)ifforanysequenceukinD(A)
convergingtoalimitu∈XsuchthatthesequenceAukisconvergingtoalimitv∈X
impliesu∈D(A)andAu=v.ThesetofboundedlinearoperatorsmappingXintoitself
isdenotedbyB(X).Againwewrite(z−A)insteadof(zI−A),whereIdenotesthe
identityontheconsideredspace.
Thespectrumσ(A)ofaclosedoperatorisdefinedbyallz∈Csuchthateither(z−A)
isnotinvertibleor(z−A)−1isnotboundedonX,seeDefinition4.5.
Asshownin[103],ifz∈/σ(A),then(z−A)−1≥dist(z,σ(A))−1.Therefore,weuse
thesameconventionasbeforebysetting(z−A)−1=∞forz∈σ(A)andproceedwith
theons.definitianalogousDefinition2.6LetA∈C(X)andε>0.Thesetsσε(A)definedby
σε(A)={z∈C:(z−A)−1>ε−1},(2.18)
σε(A)={z∈C:z∈σ(A+E)forsomeE∈C(X)withE<ε},(2.19)
σε(A)={z∈C:(z−A)v<εforsomev∈D(A)withv=1}(2.20)
arecalledε-pseudospectrumofA.
Theorem2.7Theabovedefinitions(2.18),(2.19)and(2.20)areequivalent.
Proof.See[103].✷

Theorem2.8LetA∈C(X)andε>0.
(1)σε(A)isnonemptyandopen.

ectraPseudosp2.3.

19

(2)If0<ε1≤ε2,thenσε1(A)⊆σε2(A).
(3)ε>0σε(A)=σ(A).
(4)σε(A)⊇σ(A)+B(ε).
(5)Anyboundedcomponentofσε(A)hasanonemptyintersectionwithσ(A).
Asstatedintheprecedingtheorem,somebasicpropertiesaresimilarasinthefinite
dimensionalcase(cf.Theorem2.4).However,inthegeneralcasepseudospectracanbecome
unbounded,andonlyforboundedcomponentsitcanbeassuredthattheycontainelements
ofthespectrum.

2.3.3BoundsonMatrixandOperatorExponentials
Instudyinglinearstability,equationsoftheform
d(2.21)Au=udtthearise.generaHere,lAsolutionisaoflinear(2.21op)eisratorgivenactingbyuon(t)a=eBanactAuh,spacewithanX.IfinitialAisvbalueounudedandonXthe,
00exponentialfunctiondefinedby
∞etA=1tnAn.(2.22)
!k=0kandTheopanalyticeratorwithexpresponenecttialtocontvinergestheforcompleanyxplane,complexseen[um59b].erIft.AisMoreounbver,ounded,itisbdefinitioundedon
(2.22)isnotapplicablebecausethedomainofAnmaybecomenarrowerforincreasingn.
Inthiscasetheoperatorexponentialcanbegeneralizedbymeansofsemigroups,see[79].
InthissectionwepresentlowerandupperboundsforetA,whereweassumeAtobe
vabalidforoundedclosopederatoropeonratorsaBanadefiningchaspacCeXsemigroup.withnormFor∙this.Somediscussionofwtheseereferresulttos[79are,103also].
0Inthesequel,weassumeXtobeaBanachspace.
Definition2.9LetA∈B(X).ThespectralabscissaofAisdefinedby
α(A)=z∈σsup(A)Rez.
Moreover,theε-pseudospectralabscissaofAisdefinedby
αε(A)=z∈σεsup(A)Rez.
Theorem2.10LetAbeaboundedoperator.Thenthereexistconstantsω∈Rand
M≥1,suchthat
etA≤Meωt(2.23)

20HYDRODYNAMICSTABILITY
holdsforallt≥0.Foranyz∈CwithRez>ωwehavez∈/σ(A),and,moreover,
∞(z−A)−1=e−tzetAdt,(2.24)
0and1(z−A)−1≤Rez−ω.(2.25)
Furthermore,foranyclosedcontourΓenclosingσ(A)initsinteriorwehave
1etA=2πiΓetz(z−A)−1dz.(2.26)
Proof.Theinequality(2.23)followsdirectlyfrom(2.22).Fortheproofsof(2.24),(2.25)
and(2.26)wereferto[79].✷
Theorem2.11Letε>0andLεdenotethearclengthoftheboundaryofσε(A)orofits
Thenl.hulonvexcetA≤Lεetαε(A)(2.27)
επ2holdsforallt≥0.
Proof.SupposethatLε<∞,otherwise(2.27)obviouslyholds.Then,by(2.26)wehave
foranyt≥0andε>0
etA≤21π|etz|(z−A)−1dz
1∂σε(A)
≤2πε∂σε(A)etRezdz
≤1etαε(A)Lε.
επ2Clearly,theseestimatesapplyfortheconvexhullofσε(A)aswell.✷
loBywerbreplacingoundsσdε(Aealing)bywithitspsconvexeudosphullectra,onewemaypresenreducetsometheresultsconstantinvLε.olvingBeforethewspeectralstate
abscissa.Theorem2.12LetA∈B(X).Thenwehave
etA≥etα(A)
forallt≥0and,moreover,forthestrictLyapunovexponent
1tlim→∞tlogetA=α(A).

ectraPseudosp2.3.

21

(2.28)

Proof.Aproofcanbefoundin[103].✷
Nowweareabletostatethemainresultofthischapterwhichestablisheslowerbounds
onthenormoftheoperatorexponential.
Theorem2.13LetAbeaboundedoperatoractingonaBanachspace.
(1)Theε-pseudospectralabscissaαε(A)isfiniteforallε>0.
(2)Ifz∈CwithRez>0,then
supetA≥Rez(z−A)−1.(2.28)
0≥thaveeW(3)t≥sup0etA≥K(A),(2.29)
ewherK(A):=supαε(A)=sup(Rez)(z−A)−1
ε≥0εRez>0
denotestheKreissconstant.
(4)Ifa=Rez>0andL=Rez(z−A)−1,then
aτsupetA≥eτa1+e−1(2.30)
Lτ≤<t0

(2.29)

(2.30)

forallτ>0.
(5)Leta=RezandL=Rez(z−A)−1.SupposethatetA≤Mforallt≥0and
L/M∈(−∞,1].Thenwehave
eτA≥eτa−eτa−1=1−(eτa−1)(1−L/M),(2.31)
ML/ML/andLτeτA≥1−(z−A)−1.(2.32)
Proof.WesetM=supt≥0etAandsupposeM<∞.Thus,wecanchooseω=0in
(2.23).Byvirtueof(2.24)wehaveforanyzwithRez>0
∞(z−A)−1≤MetRezdt=M,
zRe0whichimplies(2.28)andhence(2.29).
Toprove(2.30),wesetMτ=sup0<t≤τeτAandsupposeMτ<∞.Thus,etA≤Mτ
for0<t≤τ,etA≤Mτ2forτ<t≤2τ,andsoon.By(2.24)thisimplies
∞∞(z−A)−1≤(j+1)τe−taMτj+1dt=τe−tadte−τajMτj+1,
j=0jτ0j=0

22

ABILITYSTYNAMICODHYDR

whereinthelaststeptwassubstitutedbyt+jτ.IfMτ≥eτa,wehaveτa>0,L>1
andhence(2.30)holds.SoweassumeMτ<eτaandwemaysumupthegeometricseries
toobtain−11−τaMτeτa−1
(z−A)≤a(1−e)1−Mτe−τa=a(eτa/Mτ−1).
Invertingthisexpressionleadsto
a1a(eτa/Mτ−1)
L=(z−A)−1≥eτa−1,
andconsequentlyeτaeτa−1
Mτ−1≤L,
whichproves(2.30).
Using(2.30)weshowthatαε(A)isfiniteforeachε>0.Supposethecontrary,that
a=Rezbecomesarbitrarilylargeforsomevalueof(z−A)−1=ε−1.Settingτ=c/a
in(2.30)forsomec>0yields
csupeτa≥ec1+(e−1)ε.
ac/a≤<t0Takinga→∞showsthateτamustbearbitrarilylargeforarbitrarilysmallt,contra-
).2.23(dictingToprove(2.31)wesetP=eτAandfor0≤t≤τusing(2.23)weconclude
etA≤M,e(τ+t)A≤PM,e(2τ+t)A≤P2M,
andsoon.IfP≥eτawehaveL/M≥1,aL≥0whichimplies(2.31).ForP<eτaweget
inthesamemannerasintheproofof(2.30)
∞τL≤e−tadte−τajPjM=1(1−e−τa)M−τa,
a0j=0a1−Pe
andhenceL1−e−τa
M≤1−Pe−τa.
us,ThPe−τa≥1−1−e−τa,
ML/whichimplies(2.31).
Finallytakinga→0andL→0usingl’Hôpital’sruleproves(2.32).✷
Formatrices,wecanalsoformulateanupperboundbymeansoftheKreissconstant,
].103[seeTheorem2.14(KreissMatrixtheorem)ForanyA∈Cn×n,
tAK(A)≤t≥sup0e≤enK(A)

ectraPseudosp2.3.

holds,whereK(A)istheKreissconstantofAasdefinedinTheorem2.13.

23

2.3.4PseudospectraofMatrixPencils
AsitisshowninSection4.1,thediscretizationofaneigenvalueproblembymeansoffinite
elementmethodsleadstoageneralizedmatrixeigenvalueproblemoftheform
(2.33)u.λM=AuOnecandefinethepseudospectrumasσε(M−1A)asin[85].Thismaybethenatural
definition,however,besidesassumingMtobenonsingular,onemightpreferthegeneral-
izedformforcomputationalreasons.Adifferentdefinitionisestablishedin[84],wherethe
matrixMisassumedtobeHermitianpositivedefinite.Inthiscasetheε-pseudospectrum
canbedefinedbymeansoftheCholeskyfactorizationM=FHFastheε-curveof
(z−F−HAF−1)−1,whereF−Hisshorthandfor(FH)−1.Itiseasytoshowthatas-
sumingM=FHFandchoosing∙=∙2,thelattertwodefinitionsareequivalent,see
].103[However,Mmaybesingularasforinstancebydiscretizingtheincompressibilitycon-
straintoftheNavier-Stokesequations.Thisleadsustothedefinitionsfollowedin[43],
wherethepseudospectrumisdefinedbyperturbingbothAandMindependently.
Definition2.15Letε>0andα,µ>0.ThenwedefineforA,M∈Cn×ntheε-
pseudospectrumofthematrixpencil(2.33)equivalentlyby
σε(A,M)={z∈C:(zM−A)−1>(ε(α+|z|µ))−1},
σε(A,M)={z∈C:z∈σ(A+Aδ,M+Mδ)forsomeAδ,Mδ∈Cn×nwith
Aδ<εαandMδ<εµ},
σε(A,M)={z∈C:(zM−A)u<ε(α+|z|µ)forsomeu∈Cnwithu=1}.
Sincewearemainlyinterestedinε-pseudospectraaroundtheorigin,wesetα=1and
µ=0,whichisalsothedefinitionusedin[108,109].Inthiscase,Definition2.15is
applicableifwereplacethesecondidentityby
σε(A,M)={z∈C:z∈σ(A+Aδ,M)forsomeAδ∈Cn×nwithAδ<εα}.
FornonsingularmatricesM,thespectrumofthematrixpencilσ(A,M)coincideswith
thespectrumσ(M−1A).However,thispropertydoesnotholdforpseudospectraingeneral,
butthefollowinginclusionsapply
σε/M(M−1A)⊂σε(A,M)⊂σεM−1(M−1A).(2.34)
Fromthelastinclusionweknowthatσε(A,M)isboundedfornonsingularM.Onthe
contrary,thisisnotvalidforsingularmatricesMingeneralasshowninthenexttheorem,
].108[seeTheorem2.16LetA,M∈Cn×n.
(1)IfMisnonsingular,thenσε(A,M)isboundedforanyε>0.
Au(2)IfMissingular,thenσε(A,M)=Cforε>ε∗=Mu=0min,u=0u.

24

HYDRABILITYSTYNAMICOD

(3)IfMisnotthenullmatrix,thenσε(A,M)=∅forallε>0.
Proof.Thefirststatementisadirectconsequenceof(2.34).Thenextstatement(2)follows
.2.15DefinitionfromToprove(3),supposeσ(A,M)=∅(otherwise,forz∈σ(A,M)wehave(zM−A)−1=
∞byconvention).LeteHidenotethe−1ithunitvectorofCn.For1≤i,j≤nwedefine
thefunctionsϕi,j(z)=ei(zM−A)ej,whichareanalyticinC.Asaconsequenceof
Liouville’stheorem,eachofthesefunctionsϕi,jareeitherunboundedorconstant.Ifat
leastoneϕi,jisunbounded,then(zM−A)−1≥ϕi,j(z)→∞forsome|z|→∞and
thereforeσ(A,M)ε=∅−for1eachε>0.Ontheotherhand,ifallϕi,jareconstant,then
theentriesof(zM−A)areindependentofz,andhenceMisthenullmatrix.✷
WehaveseeninTheorem2.4thatσε(A)containsthediskswithradiusεaroundthe
eigenvalues.Furthermore,ifAisnormal,σε(A)isevenidenticaltotheunionofthese
disks.Formatrixpencils,thisstatementcanbegeneralizedby

σε(A,M)⊇σ(A,M)+B(ε/M),
providedthatMisnotthenullmatrix.Incontrasttopseudospectraofnormalmatri-
ces,pseudospectraofapencil(2.33)withAandMnormalcannotbedeterminedbyits
eigenvaluesalone.Infact,σε(A,M)canbemuchlargerthantheunionofdisksaround
σ(A,M)withradiusε/M,see[108].

3Chapter

NeutronTransportCriticality

Thesubjectofutmostimportanceforthesafetyofnuclearreactorsisthecriticalityprob-
lem.Itexaminestheevolutionofthefissionchainreactionofneutrons.Ifthisreaction
isinthedesirableself-sustainedstate,thenuclearreactorissaidtobecritical.Inasub-
criticalstate,thepopulationoffreeneutrondecays,whereasinasupercriticalstate,the
chainreactiongrowsexponentiallyleadingtoanuncontrolledexplosion.Thepopulation
ofneutronsinanuclearreactorismodeledbymeansoftheneutrontransportequation
whichfindsitsapplicationsinmedicalphysicsandnuclearradiationshieldingtechnology
ell.wasInthefollowingSection3.1wereviewthelinearBoltzmannequationformodeling
neutrontransportandexplainthebasicmechanismofanuclearreactor.Afterwards,in
Section3.2,westatethecriticalityproblembymeansofeigenvalueproblems.

TNeutron3.1ortransp

ThemodelofneutrontransportgoesbacktoLudwigBoltzmann[17]whoestablished
anintegro-differentialequationforthestudyofdilutegases.Incontrasttothefluidflow
modelofSection2.1,itisbasedonstatisticaldistributionsofparticles.Ithasnotonlybeen
appliedsuccessfullyondescribingdilutegasesbutalsoonmodelingradiativetransportin
planetaryandstellaratmospheresandneutrontransportinnuclearreactors,see[24].
Thebasicmechanismofthecurrentgenerationofnuclearreactors(typicallyusingthe
U235isotopeofuraniumasfuel)isdepictedinFigure3.1.Weconsiderthethreemost
importanttypesofinteractionsbetweenneutronsandnuclei.Freeneutronscausefission
reactionswiththenucleicontainedinthefuelrodsreleasingenergyaswellasfission
neutrons.Theseneutronsaresloweddownbythemoderator(scatteringreaction).This
enablesthechainreaction,becausesloweddownneutronsaremorelikelytocausefission
reactions.Inorderforthechainreactionnottorunoutofcontrol,controlrodscanbe
insertedinthenuclearreactorcatchingfreeneutrons(capturingreaction).Byhandling
thesecontrolrods,thenuclearreactorcanbekeptinacriticalstate.
WegiveabriefintroductionofthelinearBoltzmanntransportequationusedinreactor
physicstounderstanditsbasicmechanisms.Foramoredetaileddiscussionwereferto
[13,66].Atreatiseofnuclearphysicswithamoreintroductorycharactercane.g.befound
].115[in

25

26

barrierprotection

dsrotrolcon

CRITICALITYTTRANSPORONNEUTR

Figure3.1:Nuclearreactorprinciple.

rofuelds

ratordemo

AssumptionsGeneral3.1.1Forthemodelusedinthiswork,westatethefollowingassumptions.Inthiscontext,a
particledenoteseitheraneutronoraphoton(gammaray).

(1)Particlesmaybeconsideredaspoints,i.e.theycanentirelybedescribedbytheir
positionandvelocity.

(2)areCollisionsemittedmaybeimmediatelyconsidered.Formoinstandelstaneous,withdei.e.layedafteraneutrcollonsisionseethe[13,66emerging].particles

(3)Onlywhichtheareexpsmallectedinvaluecomparisonofthewithparticletheadensitverageyistakparticleenindtoensitaccouny,tare,i.e.neglectedfluctuat.ions,

(4)Particlestravelinstraightlinesbetweenpointcollisions.

(5)paredSincetotheatomicparticledensities,densitieswinenmaucylearneglectreactorsandparticle-particleotherinapplicationsteractions.aresmallcom-

(6)Materialpropertiesareassumedtobeisotropic.

Foradetaileddiscussionandphysicaljustificationsoftheseassumptions,wereferto[13,
].66

ortranspTNeutron3.1.

u

uδFigure3.2:Abeamofneutronstransmittingthroughaslab,cf.[66].

27

3.1.2CrossSectionDefinitions
Atthelengthscaleconsideredhere,wecannotdeterministicallypredictthataparticlehits
atargetnucleuslikeonabilliardtable.Thesekindsofinteractionscanonlybeformulated
bymeansofprobabilitieswhicharequantifiedbycrosssections.
energyWeEconsiderandadirectionbeamuofimpingingparticlesperpwithinendiculartensitytoIa(i.e.targetIslabparticlescpomprisedersecofond),atomswithof
asingleisotopeasinFigure3.2.Themicroscopiccrosssectionσ˜(E)istheeffectivecross
−sectional242areapernucleusseenbytheparticlesandisusuallymeasuredinbarns(1barn=
10cm).Itexpressestheprobabilitythataparticleinteractswithatargetnucleus.Let
ndenotethenumberofnucleiperunitvolumeofthemedium.Thenthebeamintensity
isgovernedby
I(u+δu)=I(u)[1−nσ˜(E)δu],
whereweconsideronlyparticleswhichhavenotmadeacollisionwithnuclei.By
σ(E)=nσ˜(E)
weprobabilitdenoteytheofmacrparticleoscinopicteractioncrosspsecterion,unitwhichdistancehasofunitspartiofcleinvtraersevel.lengNoteth.Itthatexpressescomparedthe
tothemicroscopiccrosssection,wheretheprobabilityisbasedononetargetnucleus,the
macroscopiccrosssectionquantifiestheprobabilitybasedononevolumeunit.
Ifmorethanoneisotopeisconsidered,theatomdensitiesmaynotbeuniform,and,
hence,themacroscopiccrosssectionisingeneraldependentonthespatialvariabler.In
thiscasethetotalcrosssectionreads
σ(r,E)=ni(r)σ˜i(E),
iwhereweusetheindexitorepresenteachtypeofnucleus.Thetotalcrosssectionmaybe
dividedintoparticularcrosssectionsfordifferenttypesofparticlereactions.Bydenoting

28

ONNEUTRCRITICALITYTTRANSPOR

Figure3.3:TotalneutroncrosssectionforU238asafunctionofincidentneutronenergy
measuredinelectronvolt(1eV≈1.6×10−19joule);adaptedfrom[66].

themicroscopiccrosssectionforareactiontyperforanisotopeiasσ˜ri(E),wehave
σr(r,E)=ni(r)σ˜ri(E).
iWeconsiderabsorptionandscatteringreactionsseparatelybysetting
σ˜(E)=σ˜a(E)+σ˜s(E).
Foradiscussionofthesekindsofcrosssectionswedistinguishbetweenneutronandgamma
sections.crossyra

SectionsCrossNeutronInthecaseofneutronabsorption,asumofvarioustypesofphysicalreactionsoccur.In
nuclearreactorapplications,themostimportantarecaptureandfissionreactions.In
acapturereaction,thereisonlyonecapturegammarayemitted,whereasinafission
reactionameanvalueofν(E)neutronsandadditionallygammaraysareemitted.Thus,
considerewσ˜a(E)=σ˜c(E)+σ˜f(E),
whereσ˜cisthecaptureandσ˜fthefissioncrosssection.
Thescatteringcrosssectioniscomprisedofelasticσ˜n(E)andinelasticσ˜n(E)scattering
sections:crossσ˜s(E)=˜σn(E)+σ˜n(E).

ortranspTNeutron3.1.

29

Inelasticscatteringinteractionsmomentumandkineticenergyoftheincidentparticle
areWhereasconservined,inelai.e.sticthescatterinimpinginggthereparticleisalossandofthekinettargeticrenergyemainofinthetrinsicallyneutron,uncwhilehanged.the
energystateofthenucleusiselevated.
Figure3.3showsthecomplexbehaviorofneutroncrosssectionsexemplarilyforUranium
238has.tobThisedatadetermincannotedbeempiricallydetermineasadbyfunctiomeansnofofenergypropfoerrtieseachofntheuclidenanduclides.foreacRather,htypite
.teractionniof

SectionsCrossyRaGammabeUnlikedeterminedneutronbycrossfirstsectprinciplions,estheoresignificanstimatedtcompfromonentsanalyticalofgammaapproraximatycrossions.Wesectionssetthecan
absorptioncrosssectionequaltothephotoelectriccrosssection,whichinvolvesabsorption
ofaphotononreleasinganelectronfromoneoftheorbitalshells:
σ˜a(E)=σ˜pe(E).
WeassumethegammarayscatteringtobecomprisedofComptonscatteringandpair
scattering:ductionoprσ˜s(E)=σ˜cs(E)+σ˜pp(E).
Comptonscatteringconsistsofscatteringoflowerenergyphotonsbyfreeelectrons.Pair
productioninvolvesthematerializationofaphotonintoanelectron-positronpair,which
.energyphotonthesharegammaSinceraygammacrossraysectionsinterainctionsstatingdothenotneutroninfluencethetranspportopulationequationofinneutrothens,followewingcansection.neglect

3.1.3TheLinearBoltzmannEquation
Inthesequel,letΩdenotetheneutrondirectionoftravelasitiscommonnotationin
neutrontransporttheory.ThereshouldbenoconfusionwithΩdenotingadomaininfluid
problems.wfloForathree-dimensionalproblemweneedsevenindependentvariablestodescribethe
distributionofneutrons:threespatialcoordinatesr,twoanglesfortheneutrondirection
oftravelΩ,theneutronenergyE,andthetimet.Notethattheenergyofaneutroncan
alsobeexpressedbyitsvelocity.
Inordertostatetheneutrontransportequation,wefirstneedtointroducesomenota-
tions.Thescatteringcrosssectionσs(r,E→E,Ω∙Ω)isdefinedsuchthatσs(r,E→
E,Ω∙Ω)dEdΩquantifiestheprobabilityperunitdistanceoftravelthataneutronat
positionrwithenergyEtravelingindirectionΩwillscatterintoanenergyintervaldE
aboutEintoasolidangleintervaldΩaboutΩ.
Letν(E)denotethemeannumberofneutronsproducedinafissionbyaneutronwith
energyE.Furthermore,letχ(E)bedefinedsuchthatχ(E)dEisequaltotheprobability
thataneutronproducedinafissionwillhaveanenergywithinanenergyintervaldEabout
.ETheneutronangulardensity,denotedbyN(r,Ω,E,t),isdefinedastheexpectednum-
berofneutronsatpositionrwithdirectionΩandenergyEattimetperunitvolumeper

30

CRITICALITYTTRANSPORONNEUTR

unitsolidangleperunitenergy.Theneutronangularfluxisdefinedby
ψ(r,Ω,E,t)=vN(r,Ω,E,t),
wherev=v2istheparticlespeed.
Thetransportequationexpressedintermsoftheangularfluxψreads
v1∂∂t+Ω∙+σ(r,E)ψ(r,Ω,E,t)=q(r,Ω,E,t),(3.1)
wherethenablaoperatoroperatesonlyonthespatialvariablerandqdenotesthe
emissiondensity,cf.[13,66].Theemissiondensityqiscomprisedofcontributionsof
externalsourcesqex,scatteredneutronsqs,andfissionneutronsqf,i.e.
q=qex+qs+qf.
Externalsourcesareconsideredtobeknownandtobeindependentofthefluxψ.The
emissiondensityforscatteredneutronsisgivenby
∞qs(r,Ω,E,t)=σs(r,E→E,Ω∙Ω)ψ(r,Ω,E,t)dΩdE,
S0whereSdenotestheunitsphere.Inourcaseweneglectdelayedneutrons,seeSection3.1.1,
andthus,theemissiondensityforfissionreads
∞qf(r,Ω,E,t)=χ(E)ν(E)σf(r,E)ψ(r,Ω,E,t)dΩdE.
S0Insertingqinequation(3.1)implies
v1∂∂t+Ω∙+σ(r,E)ψ(r,Ω,E,t)=qex(r,Ω,E,t)
+∞σs(r,E→E,Ω∙Ω)ψ(r,Ω,E,t)dΩdE(3.2)
S0∞+χ(E)ν(E)σf(r,E)ψ(r,Ω,E,t)dΩdE.
S0LetVdenotethedomainwithinwhichtheneutrontransportproblemistobesolvedand
letΓ=∂Vbeitsboundary.Tosolvethetransportequation(3.2),onehastoimposean
initialconditionψ(r,Ω,E,0)att=0aswellasappropriateboundaryconditions.Letn
denotetheunitvectornormaltoΓpointingoutwards.Thentheincomingfluxisspecified
ybψ(r,Ω,E,t)=ψin(r,Ω,E,t),n∙Ω<0,r∈Γ.
Thecommoncasewhereψin=0isreferredtoasvacuumboundarycondition.Reflective
boundaryconditions,whichprescribethatalloutgoingneutronsarereflectedback,are
ybharacterizedcψ(r,Ω,E,t)=ψ(r,Ω,E,t),n∙Ω=−n∙Ω,Ω×Ω∙n=0,r∈Γ,
whereΩdenotestheincidentdirectionandΩthereflectiondirection.

yCriticalit3.2.

31

yCriticalit3.2Thestabilityofasystemcontainingfissilenuclidesischaracterizedbyitspopulationof
freereactionneutrodiesns.outIfinthentime,umbertheoffreesystemisneutronsaidistobdecaesubyingcriticinaltime,.Onwhicthehothemeansrthehand,itfissionis
calledsupercriticalifthepopulationoffreeneutronsisgrowing.Finally,itisdefinedto
becriticalifthenumberoffreeneutronreachesatimeindependentequilibriuminthe
absenceofexternalsourcesofneutrons,i.e.ifthereexistsasteadynonnegativesolutionof
thesourcefreetransportequation(3.2)
∞Ω∙+σ(r,E)ψ(r,Ω,E)=0Sσs(r,E→E,Ω∙Ω)ψ(r,Ω,E)dΩdE
(3.3)∞+χ(E)0Sν(E)σf(r,E)ψ(r,Ω,E)dΩdE,
withTheseacriticalitppropriateycboundaryharacterizationscondicantionsbe(e.g.reformvacuumulatedoraseigenreflectivvealueboundaryproblems.conditions).

3.2.1TheαEigenvalue
Supposethatwehaveasymptoticsolutionsoftheform
ψ(r,Ω,E,t)=ψα(r,Ω,E)eαt
satisfyingtheimposedinitialandboundaryconditions.Insertingthistothesourcefree
formulationof(3.2)yieldstheeigenvalueproblem
∞α+Ω∙+σ(r,E)ψα(r,Ω,E)=σs(r,E→E,Ω∙Ω)ψα(r,Ω,E)dΩdE
v0S∞
+χ(E)ν(E)σf(r,E)ψα(r,Ω,E)dΩdE.
S0Assumethatthereexistsanexpansionofthesolutionψαineigenfunctionsψi.Letα0
denotetheeigenvaluehavingthelargestrealpartandψ0anassociatedeigenfunction.
Forlarget,weexpectthatthesolutionoftheinitialvalueproblemisproportionalto
ψ0eα0t.Furthermore,weassumeforphysicalreasonsthatα0isreal,otherwisenegative
orimaginarydensitiescouldoccur.Thus,wemakethefollowingdistinctionbythesignof
α0forcharacterizingthecriticality:
>0:supercritical,
<0:subcritical.
α0=0:critical,

3.2.2ThekEigenvalue
Toadjustingderivethethenkumbeigenerofvalueneutronsform,weemittedassbyumefission.thattheThatsystemmeanswcanebcanemadereplaceνcriticalbyνb/ky

32

TRANSPORONNEUTRCRITICALITYT

(3.4)

in(3.3)resultingintheeigenvalueproblem
Ω∙+σ(r,E)ψ(r,Ω,E)=
∞σs(r,E→E,Ω∙Ω)ψ(r,Ω,E)dΩdE(3.4)
S0+χ(E)∞ν(E)σf(r,E)ψ(r,Ω,E)dΩdE.
kS0Ifthelargesteigenvaluekisequaltoone,thecriticalitycondition(3.3)isclearlysatisfied.
Thecasek<1meansthatthenumberofneutronsperfissiontomakethesystemcritical
islargerthantheactualν,i.e.thesystemissubcritical.Conversely,itissupercriticalfor
k>1.Wrappingup,wehave:

critical,:1=k>1:supercritical,
<1:subcritical.
Ifthesystemiscritical,i.e.α0=0andk=1,thecorrespondingeigenfunctionsare
identical.However,foranysubcriticalorsupercriticalsystem,theeigenfunctionsmay
differ.Intheαeigenvalueproblemwehaveanadditionaltermα0/vwhichisreferredto
astimeabsorption.Sincethistermmaycausedifficultiesinnumericalcomputations,one
treatsthecriticalityproblembyevaluatingthekeigenvalueratherthantheαeigenvalue
inmanyapplications,see[13,66].

4Chapter

DiscretizationoftheEigenvalue
Problems

Inordertosolvetheeigenvalueproblemsarisinginthecontextofhydrodynamicstabil-
ity(2.13)andcriticality(3.4)numerically,weneedappropriatediscretizationtechniques.
InthefollowingSection4.1wederiveapriorierrorestimatesforthespectrumandthe
pseudospectrumintermsofafiniteelementapproximation.Afterwards,inSection4.2we
treatthediscretizationofthekeigenvaluecriticalityproblem.

4.1GalerkinFiniteElementSpectralApproximation
Thespectralapproximationtheorywepresentisputintheframeworkofellipticdifferential
operatorswhichariseinthefieldoffluiddynamics.Itisformulatedbymeansofabounded
operatorwithacompactinverseinordertoapplytheapproximationtheoryofcompact
operators.Thispresentationfollows[16,53].

4.1.1ProblemFormulationforEllipticEigenvalueProblems
VariationalFormulation
LetAbealineardifferentialoperatordefinedonaboundeddomainΩ⊂Rdoforder2m:
Au=(−1)|β|∂βaαβ(x)∂αu.(4.1)
|α|≤m|β|≤m
WeassumetheoperatorAtobeuniformlyelliptic,i.e.thereexistsaconstantε>0such
thataαβ(x)ξα+β≥ε|ξ|2m∀x∈Ω,∀ξ∈Rd.
|α|,|β|=m
Theclassicalformulationoftheeigenvalueproblemforthisoperatorreads
Findλ∈C,u=0suchthatAu=λuinΩ,
Bju=0on∂Ω,j=1,...,m,(4.2)
whereBjdenoteappropriateboundaryoperators.NotethatinhomogeneousDirichlet
boundaryconditionsarenotallowedsinceweneedthespaceoffunctionsfulfillingthe

33

34

DISCRETIZATIONOFTHEEIGENVALUEPROBLEMS

boundaryconditionstobeavectorspace.
Theassociatedvariationalformulation(orweakformulation)of(4.2)isgivenby:
Findλ∈C,u∈Vsuchthata(u,ϕ)=λ(u,ϕ)0∀ϕ∈V,u0=1,(4.3)
wherea:V×V→Cisthesesquilinearformgeneratedbythe2operatorA.WeassumeV
tobeanappropriate(complex)SobolevspacesuchthatV⊂L(Ω)⊂VbuildsaGelfand
triple,i.e.V⊂L2(Ω)iscontinuouslyanddenslyembedded(forinstanceV=H0m(Ω)).
Here,VdenotesthecontinuousdualspaceofV.
InthecaseofhomogeneousDirichletboundaryconditions,thecorrespondingsesquilin-
readsformeara(u,ϕ)=aαβ(x)∂αu∂βϕdµ.(4.4)
Ω|α|,|β|≤m
Theadjointeigenvalueproblemassociatedto(4.3)seeksλ∗∈Candu∗∈Vsuchthat
a(ϕ,u∗)=λ∗(u∗,ϕ)0∀ϕ∈V,u∗0=1.(4.5)
Clearly,theprimalandadjointeigenvaluesarerelatedtoeachotherbyλ∗=λ.
Further,weassumethesesquilinearforma(∙,∙)tobeV-coercive.Thismeansitis
bounded(orcontinuous),i.e.
|a(u,v)|≤CBuVvV∀u,v∈V,
withaconstantCB>0,andthereexistCK∈RandCE>0suchthat
Rea(v,v)≥CEvV2−CKv20∀v∈V.(4.6)
ForV=H0m(Ω),assumingsufficientlysmoothcoefficientsaαβ,thiscoercivitypropertyis
derivedbytheuniformellipticityofAusingGårding’stheorem,seee.g.[21,120].
Thesesquilinearforma˜:V×V→Cdefinedbya˜(u,v)=a(u,v)+CK(u,v)0possesses
thesameeigenfunctionsasa(∙,∙),andtheeigenvaluesofa˜(∙,∙)arethesameeigenvalues
ofa(∙,∙)shiftedbytheconstantCK.Therefore,withoutlossofgenerality,wecanassume
CK=0in(4.6),i.e.a(∙,∙)isV-elliptic.
Furthermore,weassumetheembeddingV⊂L2(Ω)tobecompact.ForV=H0m(Ω)or
V=Hm(Ω)andassumingtheboundary∂Ωtobesufficientlysmooth,thisisaconsequence
oftheRellich-Kondrachovcompactnesstheorem,see[3,116].
LetA∈B(V,V)betheuniquelinearoperatorassociatedtothesesquilinearforma(∙,∙),
i.e.(Au)(v)=a(u,v).LetB(X,Y)denotethespaceofboundedlinearoperatorsfromX
toY.DuetotheV-ellipticityofa(∙,∙),wehavethatA−1∈B(V,V).
Definition4.1Withthecontinuous,denseand,compactembeddingsIV→L2(Ω),IL2(Ω)→V
andIV→V,wedefinetheoperators
T:L2(Ω)→L2(Ω),T=IV→L2(Ω)◦A−1◦IL2(Ω)→V,
T˜:V→V,T˜=A−1◦IV→V.
Lemma4.2TheoperatorsTandT˜arecompact.
Proof.TandT˜aredefinedascompositionsofcompactandboundedlinearoperatorsand
✷ct.compathereforeare

4.1.GalerkinFiniteElementSpectralApproximation

35

Notethatλisaneigenvalueofthevariationalformulation(4.3)ifandonlyifµ=λ1isan
eigenvalueofT.ThisholdsforT˜respectively.

DiscretizationtElemenFiniteGalerkinLetVh⊂Vbeafinite-dimensionalsubspaceendowedwiththenorm∙V.Theapprox-
imationsoftheeigenvalueproblems(4.3)and(4.5)byaGalerkinfiniteelementmethod
read:Findλh∈C,uh∈Vhsuchthata(uh,ϕh)=λh(uh,ϕh)0∀ϕh∈Vh,uh0=1.(4.7)
Findλh∗∈C,uh∗∈Vhsuchthata(ϕh,uh∗)=λh∗(uh∗,ϕh)0∀ϕh∈Vh,uh∗0=1.(4.8)
Again,wehavethatλh∗=λh.Inalgebraicnotationthisresultsinthefollowinggeneralized
problemsalueveigenAhxh=λhMhxh,(4.9)
AhHxh∗=λh∗Mhxh∗,(4.10)
whereAh=a(ϕjh,ϕih)i,jisthestiffnessmatrixandMh=(ϕhj,ϕhi)0i,jisthesym-
metricandpositivedefinitemassmatrix.Here{ϕih}denotesabasisofVh.Ifwehavereal
coefficientsaαβandchoosethebasis{ϕi}tobereal,weobtainrealmatricesAhandMh.
Analogoustothecontinuouscase,wedefineAh−1:V→Vhastheboundedlinear
operatorsuchthatforf∈VwehaveAh−1f=uh,whereuhisthesolutionof
a(uh,ϕh)=(f,ϕh)0∀ϕh∈Vh.
NowweareabletodefinethediscretecounterpartofTandT˜.
Definition4.3WiththeembeddingsIVh→L2(Ω),IL2(Ω)→V,IVh→VandIV→V,wedefine
atorseroptheTh:L2(Ω)→L2(Ω),T=IVh→L2(Ω)◦Ah−1◦IL2(Ω)→V,
T˜h:V→V,T˜=IVh→V◦A−1◦IV→V.
Lemma4.4TheoperatorsThandT˜harecompact.
Proof.ThandT˜haredefinedascompositionsofcompactandboundedlinearoperators
✷compact.thereforeareandAgain,notethatλhisaneigenvaluesolutionofthediscreteformulation(4.7)ifandonly
ifµh=λ1hisaneigenvalueofTh.ThisholdsforT˜hrespectively.
4.1.2SpectralApproximationofCompactOperators
ThepreviouslydefineddifferentialoperatorsAandAhhavecompactinverses.Therefore,
wecanapplythespectralapproximationtheoryofcompactoperatorsbasedon[75].We
startwithsomemainresultsaboutcompactoperatorsincludingtheRiesz-Schaudertheory,
seee.g.[59,113,120].

36

DISCRETIZATIONOFTHEEIGENVALUEPROBLEMS

Inthesequel,wewrite(z−T)insteadof(zI−T),whereIdenotestheidentity.
Throughoutthissection,letXbeaBanachspace.
Definition4.5LetT∈C(X)beaclosedlinearoperator.
(1)Theresolventsetρ(T)ofTisdefinedby
ρ(T)={z∈C:(z−T)−1existsinB(X)}.
(2)Foranyz∈ρ(T)
Rz=Rz(T)=(z−T)−1
definestheresolventofT.
setThe(3)σ(T)=C\ρ(T)
iscalledthespectrumofT.
Anyµ∈σ(T)withN(µ−T)={0}iscalledeigenvalueofT,whereNdenotesthekernel.
Foraneigenvalueµ,thesetN(µ−T)iscalledeigenspaceofTassociatedtoµ,andany
elementu∈N(µ−T),u=0iscalledeigenvectororeigenfunctionofT.
Theorem4.6(Riesz-Schaudertheory)LetT:X→Xbeacompactlinearoperator.
(1)Thesetσ(T)\{0}consistsofatmostcountablymanyelements.
(2)Eachµ∈σ(T)\{0}isaneigenvalueofT.
(3)Foranyε>0theset{λ∈σ(T):|λ|≥ε}isfinite,i.e.thespectrumσ(T)possesses
atmostoneaccumulationpointat0.
(4)Foranyµ∈σ(T)\{0}theascentαµdefinedby
αµ=max{α∈N:N(µ−T)α−1=N((µ−T)α)}
isfinite.Moreover,thegeometricmultiplicitydefinedbydimN(µ−T)andthe
algebraicmultiplicitydefinedbydimN((µ−T)αµ)arefinite.Forµ∈σ(T)\{0},the
elementsofN((µ−T)αµ)arecalledthegeneralizedeigenfunctionsofTassociated
.µto(5)Foreachµ∈σ(T)\{0}wehavethedirectdecomposition
X=N((µ−T)αµ)⊕R((µ−T)αµ),
whereRdenotestherange.
(6)Ifµ∈σ(T)\{0},thentheidentityσ(T|R((µ−T)αµ))=σ(T)\{µ}holds.
Theorem4.7LetT∈C(X).
(1)Theresolventsetρ(T)isopeninC.
(2)ThefunctionCz→Rz(T)∈B(X)isanalyticineachconnectedcomponentof
.)T(ρ

4.1.GalerkinFiniteElementSpectralApproximation

37

(3)IfTiscompact,thenintheneighborhoodofaneigenvalueµwithascentαµthe
resolventcanbeexpandedinaLaurentseries:
∞1Rζ(T)
Rz(T)=Ak(z−µ)k,whereAk=2πi(ζ−µ)(k+1)dζ,
k=−αµΓ
andΓisaJordancurveenclosingµandcontainingorintersectingnootherelements
ofσ(T).
Definition4.8LetT∈B(X)becompact.Thenfor0=µ∈σ(T)theDunfordintegral
bygiven1E(µ,T)=2πiΓRz(T)dz,
whereΓisaJordancurveenclosingµandnotcontainingorintersectinganyotherelements
ofσ(T),iscalledspectralprojectorofT.ByCauchy’sintegraltheorem,thisdefinition
doesnotdependonthechoiceoftheJordancurveΓ.
NotethatE(µ,T)correspondstotheresidueofRz(T),i.e.E(µ,T)=A−1.Inthefollowing
theoremwepresentsomepropertiesofthespectralprojectorwhichjustifyitsdenotation.

Theorem4.9LetT∈B(X)becompact.
(1)Ifµ∈σ(T)\{0},thenE(µ,T)isaprojection.
(2)Ifµ1,µ2∈σ(T)\{0}andµ1=µ2,thenE(µ1,T)E(µ2,T)=0.
(3)Letαµdenotetheascentoftheeigenvalueµ∈σ(T)\{0}.ThentherangeR(E(µ,T))
ofthespectralprojectorcorrespondstothespaceofgeneralizedeigenfunctionsasso-
ciatedtoµ,i.e.
R(E(µ,T))=N((µ−T)αµ).
Inthefollowinglemmaweshowthatanapproximationofaboundedlinearoperator
approximatesthespectrumandtheresolventaswell.
Lemma4.10LetF∈B(X)andε>0.Thenthereexistsδ>0suchthatforallS∈B(X)
withF−S<δwehave
σ(S)⊂Bε(σ(F)),
Rz(F)−Rz(S)<εforanyz∈Bε(σ(F)),
whereBε(σ(F))denotestheε-diskaroundσ(F),i.e.
Bε(σ(F))={µ+ξ:µ∈σ(F),|ξ|<ε}.
Proof.i)First,weshowthatthemappingA→A−1iscontinuousforinvertiblebounded
islinearopen.opAnyerators.B∈B(AssumeX)withA∈BB(−X)AisinAv−1er<tible.1isTheinvsetertibleofinvandwertibleehaopveeratorsinB(X)
∞B−1=(A−(B−A))−1=A−1(I−(B−A)A−1)−1=A−1(B−A)A−1k.
=0k

38

DISCRETIZATIONOFTHEEIGENVALUEPROBLEMS

obtainewus,Th∞21−A−1−B−1≤A−1B−AA−1k=B−AA−1→0(B→A),
k=11−B−AA
whichimpliesthecontinuityoftheinversemappingforlinearboundedoperators.
ii)SinceI−z−1F→Iasz→∞andusingi)weobtain
zlim→∞Rz(F)=zlim→∞z1(I−z−1F)=0.
Hence,thereexistsNε≥0satisfyingRz(F)≤Nεforanyz∈C\Bε(σ(F)).
iii)Letz∈C\Bε(σ(F))andS−F<δ:=Nε−1.Thenwehave
∞Rz(F)[(S−F)Rz(F)]i=(z−F)−1(I−(S−F)(I−F)−1−1
(4.11)=0i=(z−S)−1=Rz(S),
whichimpliesthat(z−S)−1existsandisboundedandthereforez∈ρ(S)holds.Conse-
quentlyforS−F<δ1wehaveσ(S)⊂Bε(σ(F)).
iv)LetS−F≤δ:=εNε2+εNε−1<Nε−1.Using(4.11)wehave
∞iNε2S−F
Rz(S)−Rz(F)≤Rz(F)(S−F)Rz(F)≤1−S−FNε<ε,
=1iwhichcompletestheproof.✷
Inthesequel,letTand{Th}h>0becompactoperatorsdefinedonaHilbertspacesuch
thath→lim0T−Th=0.(4.12)
Thisassertionissatisfiedinthecontextoffiniteelementsandwillbeprovenlater(Sec-
tion4.1.3).Weproceedbyderivingaprioriconvergenceestimatesforeigenpairsofcompact
ators.eropLetΓbeaJordancurveenclosinganeigenvalueµ∈σ(T)\{0}asinDefinition4.8.If
{Th}h>0isasequenceofcompactoperatorssatisfying(4.12),Lemma4.10showsthatfor
hsmallenoughwehaveΓ⊂ρ(Th),andtherefore,analogouslytoDefinition4.8,wecan
definetheoperator1
E(µ,Th)=2πiΓRz(Th)dz.
Forbrevity,wesetE=E(µ,T)andEh=E(µ,Th)iftheargumentsusedareapparent
text.conthefrom

Lemma4.11TheoperatorE(µ,Th)isaspectralprojectorontothedirectsumofthe
generalizedeigenfunctionspacescorrespondingtotheeigenvaluesofThenclosedbyΓ.

4.1.GalerkinFiniteElementSpectralApproximation

39

Lemma4.12Forµ∈σ(T)wehave
E(µ,T)−E(µ,Th)→0,ash→0.
ThefollowingconvergenceresultisstrongerthanLemma4.10sinceitassurestheconver-
genceofeigenvaluesaccordingtotheiralgebraicmultiplicity.
Theorem4.13Foraneigenvalueµ∈σ(T)letσµdenoteitsalgebraicmultiplicity.Given
aJordancurveΓenclosingµandexcludingallothereigenvaluesofT,thereexistsh0>0
suchthatforall0<h<h0thecurveΓenclosesexactlyσµeigenvaluesofThcountedwith
multiplicities.aicalgebrtheirTheconsideredeigenvaluesofThintheprecedingTheorem4.13whichlieinsideΓare
denotedbyµ1(h),...µσµ(h).
Definition4.14LetM,NbetwosubspacesofaHilbertspaceX.Wedefine
δX(M,N)=x∈M,supx=1dist(x,N),
andthegapδˆX(M,N)betweenMandNas
δˆX(M,N)=max(δ(M,N),δ(N,M)).
Wewriteδ(M,N)andδˆ(M,N)ifthereisnoconfusionabouttheunderlyingHilbertspace
.XTheorem4.15Forhsmallenough,thereisaconstantCindependentofhsuchthat
δˆ(R(E),R(Eh))≤(E−Eh)|R(E)≤C(T−Th)|R(E).
Proof.See[75,Theorem1].✷
Theorem4.13showsthatalleigenvaluesµ1(h),...µσµ(h)convergetoµ.However,the
individualµi(h)mightberatherpoorapproximationstoµ.Neverthelesstheirarithmetic
meanisgenerallyabetterapproximation,cf.[20].Thus,wedefine
σµ1µˆ(h)=σµµi(h).
=1iTheorem4.16Forhsmallenough,thereisaconstantCindependentofhsuchthat
|µ−µˆ(h)|≤C(T−Th)|R(E).
Proof.See[75,Theorem2].✷
Arefinedestimationisgiveninthefollowingtheorem.
Theorem4.17Letϕ1,...,ϕσµbeanorthonormalbasisofR(E)andϕk∗=E∗ϕkfor
k=1,...,σµ.Forhsmallenough,thereisaconstantCindependentofhsuchthat
1σµ
|µ−µˆ(h)|≤σµ((T−Th)ϕj,ϕj∗)+C(T−Th)|R(E)(T∗−Th∗)|R(E∗).
=1j

40

DISCRETIZATIONOFTHEEIGENVALUEPROBLEMS

Nowwestatetheconvergenceofseparateeigenvaluesµk(h)insteadoftheirarithmetic
mean.Inthisview,theorderofconvergencedecreasesbymeansoftheascentαµofµ.
Theorem4.18Letϕ1,...,ϕσµbeanorthonormalbasisofR(E)andϕk∗=E∗ϕkfor
k=1,...,σµ.Forhsmallenough,thereisaconstantCindependentofhsuchthat
σµ
=1i,j|µ−µk(h)|αµ≤C((T−Th)ϕi,ϕj∗)+(T−Th)|R(E)(T∗−Th∗)|R(E∗).
Sofar,wehaveonlystatedthateigenvectorsofthediscreteproblemapproximategener-
alizedeigenfunctionsofthecontinuousproblem(Theorem4.15).Thenexttheoremshows
thattheyevenapproximateeigenfunctionsofthecontinuousproblem.
Theorem4.19Letµ(h)beaneigenvalueofThsuchthatlimh→0µ(h)=µ.Supposethat
forhsmallenough,whisaunitvectorsatisfying(µ(h)−Th)kwh=0,wherekisapositive
integerwithk≤αlµ.Then,foranyintegerlwithk≤l≤αµ,thereisavectoruwh∈R(E)
satisfying(µ−T)uwh=0and
uwh−wh≤C(T−Th)|R(E)(l−k+1)/αµ,
whereCisaconstantCindependentofhandl.
Fortheproofsofthelastthreetheorems,weagainreferto[75].

4.1.3APrioriErrorEstimatesfortheFiniteElementApproximation
ectraSpofximationApproWestartbyderivingestimatesofT−ThwhichallowustoapplytheresultsofSec-
tion4.1.2.Inthefollowing,weassumeThtobeaquasiuniformtriangulationThofΩ,
see[19].Furthermore,letthefiniteelementspaceVh⊂V⊂Hm(Ω)bedefinedsuchthat
therestrictionofanyvh∈VhoneachcellK∈Thisapolynomialofdegreeatmostequal
tor−1.Thenthefollowinglemmaholds,see[19].
Lemma4.20(Bramble-Hilbert)Lett≥2.Then,undertheassumptionsmadeabove,
thereexistsaconstantCsuchthat
χ∈infVhv−χj≤Cht−jvtforanyv∈Ht(Ω),0≤j≤t≤r.
Furthermore,weassumethattheoperatorA(seeSection4.1.1)isHs-regular,i.e.forany
f∈Hs−2mandu=A−1f,wehavethatu∈Hs,andthereisaconstantCs=C(Ω,s)
thathsucus≤Csfs−2m,(4.13)
where2mdenotestheorderoftheellipticoperatorA(see(4.1)).Inourcontext,inequality
(4.13)isassumedtoholdfors≤r.Ifwehavesufficientlyssmoothcoefficientsaαβanda
sufficientlysmoothboundary∂Ωaswell,onecansprovetheH-regularityofAundercertain
boundaryconditions,seee.g.[50].NotethatH-regularityofAimpliesinparticularthat
alleigenfunctionshaveHs-regularity,see[50,53].

4.1.GalerkinFiniteElementSpectralApproximation

41

Theorem4.21Givenf∈Ht(Ω)andϕ∈Hs(Ω)with0≤s,t≤r−2mwherer−1is
thepolynomialorderoftheconsideredfiniteelementdiscretization,wehave
|((T−Th)f,ϕ)0|≤Cht+s+2mftϕs.
Proof.Foranyχ∈Vhwehave
1(T−Th)f2m≤CE|a((T−Th)f,(T−Th)f)|(ellipticity)
=|a((T−Th)f,Tf−χ)|(Galerkinorthogonality)
C≤CEB(T−Th)fmTf−χm(continuity).
Hence,

CB(T−Th)fm≤CEχ∈infVhTf−χm
≤C˜ht+2m−mTft+2m(Bramble-Hilbert)
≤Cˆht+mft(regularity).(4.14)
Withadualityargument(alsoknownastheNitschetrick)wederiveanestimateinL2(Ω):
Foranyχ∈Vhwehave
|((T−Th)f,ϕ)0|=|a((T−Th)f,T∗ϕ)|
=|a((T−Th)f,T∗−χ)|(Galerkinorthogonality)
≤CB(T−Th)fmχ∈infVhT∗ϕ−χm(continuity)
≤CBCˆht+mftχ∈infVhT∗ϕ−χm(see(4.14))
≤Cht+mfths+mT∗ϕs+2m(Bramble-Hilbert)
≤Cht+s+2mftϕs(regularity),
whichcompletestheproof.✷
Theorem4.21nowallowsustoproveassumption(4.12)wemadeinSection4.1.2:
Corollary4.22ThconvergestoT,i.e.T−Th0→0ash→0.
Proof.Letg∈L2(Ω)andchooset=s=0inTheorem4.21.Thenwehave
m2(T−Th)g0=ϕ∈L2,supϕ0=1|(T−Th)g,ϕ)0|≤Chg0.

Theorem4.23Letr−1bethedegreeofthepolynomialsconsideredinthefiniteelement
discretization.Thenthefollowingestimates
(T−Th)|R(E)0≤Chrand(T∗−Th∗)|R(E)0≤Chr

42

DISCRETIZATIONOFTHEEIGENVALUEPROBLEMS

hold.m2−rProApplyingof.LetfTheorem∈R(E4.21).withThet=rregularit−2myandsassumption=0,we(ha4.13v)eforimpanliesyϕR(∈EL)2⊂(Ω)H(Ω).
|((T−Th)f,ϕ)0|≤C˜hrfr−2mϕ0,
hence,

(T−Th)|R(E)0=f∈R(E),supf0=1(T−Th)f0
=f∈R(E),supf0=1ϕ∈L2(Ω)sup,ϕ0=1|((T−Th)f,ϕ)0|
r≤C˜hf∈R(E),supf0=1fr−2m
r,hC≤wherethelastestimateisderivedbyusingthatallnormsareequivalentinthefinite-
dimensionalspacesR(E).
Thesecondassertioncanbeprovenanalogously.✷
Nowweareabletoprovethemaintheoremwhichstatesquantitativeconvergenceresults
theforalues.veigenTheorem4.24LetAbeaHr-regular,uniformlyellipticoperatoroforder2masdefined
in(4.1),wherer−1denotestheorderofthefiniteelementapproximation.Furthermore,
assumeλtobeaneigenvalueofAwithalgebraicmultiplicityσandascentα.Then,
forsufficientlysmallh,thereareexactlyσapproximatingeigenvalues{λh,i}i=1,...,σofthe
discreteproblem(4.9)countedaccordingtotheiralgebraicmultiplicitysuchthat
σ
σλ−1λh,i≤Cλh2(r−m),(4.15)
=1i|λ−λh,i|≤Cλh2(r−m)/αfori=1,...,σ,(4.16)
whereCλandCλareconstantsindependentofh.
Proof.Theorem4.21witht=s=r−2myields
σ((T−Th)ϕj,ϕj∗)0≤Ch2r−2m.(4.17)
=1jThefirstassertion(4.15)followsfromTheorem4.17combinedwiththeestimate(4.17)
andTheorem4.23.Thesecondassertion(4.16)followsanalogouslyfromTheorem4.18
togetherwiththeestimate(4.17)andTheorem4.23.✷
BycombiningTheorem4.15andTheorem4.23,wedirectlygainthefollowingresultabout
theconvergenceorderofthegapbetweenthegeneralizedeigenspacesandtheirapproxi-
mations.

4.1.GalerkinFiniteElementSpectralApproximation

43

Theorem4.25Forafiniteelementapproximationoforderr−1thereholdsforsufficiently
hlsmalδˆL2(Ω)(R(E),R(Eh))≤Chr,
whereCisaconstantindependentofh.
Finally,westatethequalitativeconvergenceofthegeneralizedeigenfunctions.
Theorem4.26LettheassumptionsofTheorem4.24hold.Furthermore,letAhdenotethe
linearboundedoperatorassociatedtothesesquilinearformin(4.7)andλhaneigenvalue
ofAhconvergingtoaneigenvalueλofA(i.e.of(4.3))withascentα.Supposethatfor
hsmallenough,wuhisaunitvectorsuchthat(λh−Ah)kwuh=0holdsforanyintegerk
with1≤k≤α.Then,foranyintegerlwithk≤l≤α,thereexistsavectoruhsatisfying
(λ−T)luh=0and
uh−wuh0≤Chr(l−k+1)/α,
whereCisaconstantindependentofhandl.
Proof.TheassertionfollowsdirectlyfromTheorem4.19andTheorem4.23.✷
ConvergenceresultscanalsobestatedintheHm-normbyconsideringT˜anditsdiscrete
counterpartT˜hforV=Hm(Ω)(cf.Definitions4.1and4.3).Weshowthisexemplarilyfor
theconvergenceorderofthegapbetweenR(E)andR(Eh).Usinginequality(4.14)we
thatevha(T−Th)fm≤Cˆht+mft.(4.18)
Fort=0weobtain(T−Th)fm≤Chmf0≤Chmf0andconsequently
h→lim0T˜−T˜hm=0.
Settingt=r−2min(4.18)andusinganalogousargumentsasintheproofofTheorem4.23,
obtainew(T−Th)|R(E)0≤Chr−m.
Thus,usingTheorem4.15weconcludethecorrespondingformulationofTheorem4.25:
δˆHm(Ω)(R(E),R(Eh))≤Chr−m,
whereE(andEh)isthespectralprojectorassociatedtoT˜(andT˜hrespectively).
ectraPseudospofximationApproThepreviouslystatederrorestimatesforeigenvaluesallowustoderiveapriorierror
estimateforpseudospectrawithrespecttothespectralnorm.Westartwithatheorem
quantifyingthespectralnormforcompactoperators.Inthisview,recallthatbyvirtueof
theRiesz-Schaudertheoryanycompactoperatorpossessesaneigenvaluewithmaximum
absolutevalue,cf.Theorem4.6.
Theorem4.27LetT:H→HbeacompactoperatoronarealorcomplexHilbertspace
H.Thenforthespectralnorm∙0wehave
T0=smax(T),

44

DISCRETIZATIONOFTHEEIGENVALUEPROBLEMS

wheresmax(T)denotesthelargestsingularvalueofT,i.e.thesquarerootofthelargest
eigenvalueofT∗T.
Proof.Weadapttheproofforthefinite-dimensionalcasewhichcanbefounde.g.in[98].
Wehave2∗
T0=x0sup=1(Tx,Tx)0=x0sup=1(TTx,x)0,
whereT∗Tisanormalcompactoperatorwithrealnonnegativeeigenvalues.Bythespec-
traltheoremforcompactoperators(seee.g.[113]),thereexistsanorthonormalsystem
{e1,e2,...}ofeigenvectorsassociatedtothenon-zeroeigenvalues{λ1,λ2...}ofT∗Tsuch
thatT∗Tx=λk(x,ek)0ek
kholdsforanyx∈H.Consequently,
22T0=x0sup=1kλk(x,ek)0(ek,x)0=x0sup=1kλk|(ek,x)|
≤λmax(T∗T)sup|(ek,x)|2,
x0=1k
andapplyingParseval’sidentityweconclude
T0≤smax(T).
Ontheotherhand,wederive
T02=sup(T∗Tx,x)0≥(T∗Temax,emax)0=λmax(T∗T),
=1x0whereemaxdenotesaneigenvectorassociatedtothelargesteigenvalueλmax(T∗T).This
✷of.prothecompletesTheorem4.27implies−11
A0=smin(A),
foraboundedoperatorAwithacompactinverse.
Foranyz∈C,shiftingtheellipticoperatorAin(4.1)toA−zIclearlyleavesthe
operatorelliptic.Forµ∈/σ(A)∩σ(Ah)wecandefinetheshiftedversionTµ(andThµ)of
thecompactoperatorT(andTh)analogouslytoDefinition4.1(andDefinition4.3).
Definition4.28Forµ∈/σ(A)∩σ(Ah)wedefine
Tµ:L2(Ω)→L2(Ω),T=IV→L2(Ω)◦(A−µIV→V)−1◦IL2(Ω)→V,
Thµ:L2(Ω)→L2(Ω),T=IVh→L2(Ω)◦(Ah−µIVh→V)−1◦IL2(Ω)→V
Nowweareabletostatethemaintheoremfortheconvergenceoftheresolventnorm.
Theorem4.29LetAbeaHr-regular,uniformlyellipticoperatoroforder2masdefined
in(4.1),wherer−1denotestheorderofthefiniteelementapproximation.Further,let

4.2.DiscretizationoftheNeutronTransportEquation

45

µµasc∈/entσ(ofA)∩theσ(larAhgest)andαeigenvaluedenoteof(theTµ)asc∗Tµent.ofThen,theforlarhgestsmallsingularenough,valueweςofhaveT,i.e.the
Tµ0−Thµ0≤Cςh2(r−m)/α,(4.19)
whereCςisaconstantindependentofh.
Prtorsoof.(Tµ)This∗Tµandassertion(Thµ)∗folloThµwstobyapprotheproximatofeoftheTheoremlargest4.24eigenvapplalueiedofon(Tµthe)∗Tµ.compactopera-✷

RemarksBibliographical4.1.4Therearemanypublicationsonthespectralapproximationtheoryforellipticoperators.
Inadditiontothereferencesalreadycitedwementionhere[9,20,50,76].Foraposteriori
errorestimateswereferto[53,56,57].Thecaseofnoncompactoperatorsisdescribed
in[34,35].Mixedandhybridmethodsaretreatedin[68].Withoutanyattemptof
onfinitecompleteness,elemenwtemenmethotionds[w2,e61]referasto[further19,21w,orks29,47on,sp81].ectralapproximation.Fortreatises

4.2DiscretizationoftheNeutronTransportEquation
Inthesetupofgeneralreactorsnoexactsolutionsofthetransportequationareknown.
Thisisespeciallyduetothevastamountofdetailedinformationaboutneutroncross
sectionswhichhastobeknown.Thus,sophisticatedapproximationmethodsandnumerical
methodsarenecessary.Thereexisttwogeneralapproaches:deterministicandMonteCarlo
methods.MonteCarlomethodsarestochasticmethodsthatsimulateafinitenumberof
particlehistories.Inthisworkweconfineourselvestodeterministicmethodsandrefer
to[13,66]foradiscussionofMonteCarlomethodswithrespecttoneutrontransport
problems.Thecriticalityprobleminvolvesthreeindependentvariables,namely,thepositionr,
thedirectionoftravelΩ,andtheenergyE.Intheapproachusedinthiswork,the
positiondependenceisapproximatedbymeansoffiniteelementmethods.Thedependence
onparticledirectionisexpandedasaseriesofsphericalharmonics,whereastheenergy
dependenceisapproximatedbypiecewiseconstantfunctions(multigroupapproximation).
Inthisworkwerestrictourselvestothekeigenvalueproblem(3.4).Thediscretizationof
thetimedependenttransportequationandtheαeigenvalueproblemaretreatedelsewhere,
see[66]and[63]respectively.

DiscretizationEnergy4.2.1Inordertodiscretizetheenergyvariable,theenergyintervalofinterestisdividedintoa
finitenumberofintervals(orgroups)(Eg,Eg−1),g=1,...,G.Thecrosssectionineach
groupisassumedtobeindependentoftheenergy.Thismaybeachievedforinstanceby
averagingwithineachinterval.Notethatforincreasinggtheenergyoftheassociated
decreasing.isgroup

46

DISCRETIZATIONOFTHEEIGENVALUEPROBLEMS

Weapproximatetheangularfluxwithineverygroupgbyafunctionψg(r,Ω)known
asthegroupangularflux:
ψg(r,Ω)=gψ(r,Ω,E)dE,
whereintegrationovergmeansintegrationovertheinterval(Eg,Eg−1).
Weassumefornowthattheenergydependenceisseparable,i.e.thattheangularflux
canbewrittenasaproductofaknownfunctionf(E)andthegroupangularflux:
ψ(r,Ω,E)=f(E)ψg(r,Ω).(4.20)
Notethatbydefinitionofthegroupangularfluxthefunctionf(E)isnormalizedinthe
thatsensef(E)dE=1
gholdsforanygroupg.
Integrating(3.4)betweenEgandEg−1,andsubstituting(4.20)intheresultingequa-
yieldstionsGΩ∙+σg(r)ψg(r,Ω)=σsg→g(r,Ω∙Ω)ψg(r,Ω)dΩ+1Fψ(4.21)
kS=1gforanyg=1,...,G,where
GFψ=χgνgσfg(r)ψg(r,Ω)dΩ,
kS=1gandψ=(ψ1,ψ2...,ψG)T.Themultigroupcrosssectionsin(4.21)aredefinedby
σg(r)=σ(r,E)f(E)dE,
gσsg→g(r,Ω∙Ω)=σs(r,E→E,Ω∙Ω)f(E)dEdE,
ggνgσfg(r)=gν(E)σf(r,E)f(E)dE,
χg=χ(E)dE.
andmoreoverwehaveset
gTheseparabilityassumption(4.20)isinfactnotneededtoderivethewithingroup
equations(4.21).Bydroppingthisassumptiononeevenobtainsimproveddefinitions
ofthemultigroupcrosssections.However,weskipthislengthyderivationandreferto
[13,14,66]fortheinterestedreader.
Wedefinethestreamingcollisionoperatorforagroupgas
Hg0gψg=[Ω∙+σg(r)]ψg(Ω,r)

4.2.DiscretizationoftheNeutronTransportEquation

andthegroup-to-groupscatteringoperatoras
Hg1gψg=σgs→g(r,Ω∙Ω)ψg(r,Ω)dΩ.
SMoreover,themultigrouptransportoperatorisdefinedby
Hgg=δggHg0g−Hg1g.
SettingtheblockmatrixHas
H11H12∙∙∙H1G
H21H22∙∙∙H2G
.H=.....,
HHGG1Gweobtaintheblockwiserepresentation
Hψ=1Fψ
k

ofthecriticalityproblem.

4.2.2SecondOrderEvenParityFormulation

47

(4.22)

Thesecondorderevenparityformulationallowsustoderiveself-adjointdiagonalblocks
HggintheblockmatrixH.Westartwiththewithin-groupequations(4.21),where
weassumeforsimplificationthescatteringmultigroupcrosssectiontobeisotropic.For
anisotropicscatteringthederivationbecomesmorecumbersomeandcanforinstancebe
foundin[65,106].Furthermore,forbrevityweassumetohaveonlyoneenergygroup.
Thismeanswecandropthegroupindices,andconsequentlyequation(4.21)reads
Ω∙+σ(r)ψ(r,Ω)=σs(r)ψ(r,Ω)dΩ+χνσf(r)ψ(r,Ω)dΩ.(4.23)
kSSWeproceedbysplittingtheangularfluxintoevenandoddangularparitycomponents
ψ(r,Ω)=ψ+(r,Ω)+ψ−(r,Ω),
wheretheevenparitycomponentψ+andtheoddparitycomponentψ−aredefinedby
1ψ±(r,Ω)=2[ψ(r,Ω)±ψ(r,−Ω)].
Then,weevaluate(4.23)at−Ωtogain
−Ω∙+σ(r)ψ(r,−Ω)=σs(r)ψ(r,Ω)dΩ+χνσf(r)ψ(r,Ω)dΩ.(4.24)
kSS

48

DISCRETIZATIONOFTHEEIGENVALUEPROBLEMS

Addingequations(4.23)and(4.24)yields
Ω∙ψ−(r,Ω)+σ(r)ψ+(r,Ω)=σs(r)ψ(r,Ω)dΩ+χνσf(r)ψ(r,Ω)dΩ.
kSSAnalogously,bysubtractingthesetwoequationswehave
Ω∙ψ+(r,Ω)+σ(r)ψ−(r,Ω)=0.
Finally,wecombinethetwolastequationstoeliminateψ−whichresultsinthesecond
orderevenparityformulationofthecriticalityproblem:
1−Ω∙σ(r)Ω∙ψ+(r,Ω)+σ(r)ψ+(r,Ω)=σs(r)Sψ+(r,Ω)dΩ
(4.25)+kχνσf(r)ψ+(r,Ω)dΩ.
SNotethatsinceψ+(r,Ω)=ψ+(r,−Ω),theevenparityformulation(4.25)onlyneedsto
besolvedoverhalfoftheangulardomain.
Thevacuumboundaryconditionfortheangularfluxcanbeshowntoresultinabound-
aryconditionfortheevenparitycomponent:
Ω∙ψ+(r,Ω)±σ(r)ψ+(r,Ω)=0,n∙Ω0,r∈Γ,
see[66].Thereflectingboundaryconditionsfortheevenparitycomponentaregivenby
ψ+(r,Ω)=ψ+(r,Ω),n∙Ω=−n∙Ω,Ω×Ω∙n=0,r∈Γ,
whereΩdenotestheincidentdirectionandΩthereflectiondirection.

DiscretizationSpatialandAngular4.2.3Theangularvariableistreatedmainlyusingthreetypesoftechniques.Integralmethods
arebasedonintegratingouttheangulardependencefromthetransportequation,but
theyleadtodensematrices.TheDiscreteordinate(orSN)collocationmethodconsists
inevaluatingtheneutrontransportequationalongafinitenumberofangulardirections.
Thismethodsissimpletousebutitsuffersfromanomaliesintheangularfluxdistribution
knownasrayeffects,see[66].
Inthisworkweconsiderthesphericalharmonic(orPN)methodwhichdoesnotsuffer
fromtherayeffectphenomenon.ThesphericalharmonicsarebasedonLegendrepolyno-
mialsandformanorthogonalbasisofthesquareintegrablefunctionsontheunitsphere.
Theideareliesonexpandingtheangulardependenceinaseriesofsphericalharmonics
.NorderattruncatedThenormalizedsphericalharmonicsaredefinedby
Ylm(Ω)=Ylm(ϑ,ϕ)=2l+1(l−m)!Plm(cosϑ)eimϕ,l=0,1,...,m=0,1,...,l,
4π(l+m)!
wherePlmaretheassociatedLegendrepolynomials,ϑ∈[0,π]denotesthepolarangle,and

4.2.DiscretizationoftheNeutronTransportEquation

49

ϕ∈[0,2π]theazimuthalangle.Thesefunctionsareorthonormalinthesensethat
2ππ
SYlm(Ω)Ylm(Ω)dΩ=00Ylm(ϑ,ϕ)Ylm(ϑ,ϕ)sinϑdϑdϕ=δllδmm.
Moreover,theysatisfytheparityproperty
Ylm(−Ω)=Ylm(π−ϑ,π+ϕ)=(−1)lYlm(ϑ,ϕ)=(−1)lYlm(Ω).
Hence,forexpandingtheevenparityangularflux,onlycomponentsYlmwithlevenare
needed.SplittingYlmintorealandimaginarypartsyields
Ylm(Ω)=Ylem(Ω)+iYlom(Ω),
whereYlemiseveninϕandYlomoddinϕ.
ThePNmethodforevenparityinΩusestherealandtheimaginarypartofYlmas
basisfunctions,wherelisevenwithl<N.Wedefinethevector
y(Ω)=Y00,Y20,Y21e,Y21o,Y22e,Yo22,Y40,...,YNe−1N−1,YNo−1N−1T
containingthebasisfunctionsfortheangulardependence.
Thedependenceonthespatialvariableistreatedbymeansoffiniteelementmethods.
Letfdenotethevectorcontainingafiniteelementbasisofdimensions,i.e.
f(r)=(f1(r),f2(r),...fs(r))T.
Then,thespatio-angulardiscretizationψh+ofψ+isgivenby
ψh+(r,Ω)=αijfi(r)yj(Ω),(4.26)
i,jwherethecoefficientsαijaretobedetermined.Here,yj(Ω)denotesthej-thcomponent
ofthevectory(Ω).
Inordertoderivetheweakmatrixformulationofthecriticalityproblem,oneplugsthe
spatio-angulardiscretization(4.26)intheevenparityformulationof(4.21),wheremost
commonlythescatteringcrosssectionisexpandedinLegendrepolynomials.Then,the
resultingequationistestedbytheexpansionfunctionsofψh+,i.e.wemultiplybyfunctions
whichareoftheformfi(r)yj(Ω)andintegrateoverspaceandangle.Thisresultsina
matrixformofageneralizedeigenvalueproblemseekingthelargesteigenvalueλ(which
correspondstothekeigenvaluetobeapproximated)andtheassociatedeigenvectoru
(whichholdsthecoefficientsαijofψh+foreachenergygroup)of
u.λH=uFThoughstraightforward,wedonotpresentthesetediouscalculationsandreferinsteadto
[106]foradetaileddescriptionoftheweakformanditsmatrixformulation.

50

TIONDISCRETIZA

OF

THE

ALUEEIGENV

PROBLEMS

5Chapter

EigenvalueSolversand
arallelizationP

Theeigenvalueproblemsarisinginthefieldsofhydrodynamicstabilityandcriticality
arehighlycomplex.Thisfactisespeciallyreflectedbythelargesizeandcomplexityof
theresultingdiscretizedproblems.Therefore,weneedpowerfuleigenvaluesolverswhich
inadditionallowanefficientparallelization.Inthisworkwefocusontheapplication
oftheDavidsonmethodcombinedwithdifferentpreconditioningtechniquestocompute
pseudospectraofellipticoperatorsandkeigenvaluesofthelinearBoltzmannequation.
InthefollowingsectionweestablishtheDavidsonmethodwithabriefreviewofeigen-
valueprojectionmethods.Thereafter,inSection5.2,wefocusonpseudospectracompu-
tationinthespectralnorm.InSection5.3wetreattheparallelizationschemesofbasic
linearalgebraroutinesintermsoffiniteelementdiscretizationsaswellasanadditional
parallelizationtechniqueforthecomputationofpseudospectra.

5.1TheDavidsonMethod
TheDavidsonmethodisaniterativealgorithmtocomputefewextremeigenvaluesand
thecorrespondingeigenvectorsoflargesparsematrices.Itwasoriginallydevelopedby
Davidson[32]forreal-symmetricmatricesinthefieldofquantumchemistry.Later,this
methodwasalsosuccessfullyappliedtononsymmetricproblems,see[69,88].Itisclosely
relatedtoKrylovsubspacemethodswhichareamongthemostimportantmethodsfor
solvinglargesparseeigenvalueproblems.TheconceptofKrylovmethodsreliesonthe
projectionoftheproblemontoasmallersubspace.
LetA∈Cn×nandK,LbetwosubspacesofCnofsamedimensionm.Weseekfor
λ∈Candu∈Cn\{0}suchthat
(5.1)λu.=AuTheconceptofprojectionmethodsistoapproximatetheexacteigenvectorubyavector
u˜belongingtothesubspaceKandthecorrespondingeigenvalueλbysomeλ˜suchthat
conditionov-GalerkinPetrtheAu˜−λ˜u˜⊥L(5.2)
issatisfied.LetV=[v1,v2,...,vm]denotethematrixwhosecolumnvectorsviforma
basisofKandW=[w1,w2,...,wm]thematrixwhosecolumnvectorswiformabasisof
L.SupposethatVandWarebiorthogonal,i.e.WHV=I,whereIisthem-dimensional

51

52

EIGENVALUESOLVERSANDPARALLELIZATION

Algorithm5.1Powermethod.
1:functionPower(A,k)
2:Chooseaninitialnormalizedvectorv0;
3:fork=1,2,...do
4:Computezk=Avk−1;
5:Computevk=zk/zk2;
6:ComputetheRayleighquotientαk=vkHAvk;
orfend7:ctionfunend8:

identitymatrix.Then,wecanreformulatetheprojectionmethodbyseekingλ˜∈Cand
y∈Cm\{0}satisfying
Hmy=˜λy,(5.3)
whereHm=WHAV∈Cm×mistheRayleighmatrix.Sincewehaveusuallymn,
solvingtheeigenvalueproblem(5.3)ismuchlessexpensivecomparedtosolvingtheoriginal
).5.1(problemInthecaseL=K,projectionmethodsarereferredtoasobliqueprojectionmethods,
whereasinthecaseL=K,theyarereferredtoasorthogonalprojectionmethods.General
convergenceresultsofprojectionmethodscanbefoundin[86].
Averysimpleprojectionmethodisobtainedbychoosing
L=K=span{Akv}
foranintegerk≥0.Theresultingprocedureisknownaspowermethodandgivenin
Algorithm5.1,seee.g.[49].Ateachiterationk,theeigenvalueproblemisprojected
ontotheone-dimensionalsubspacespannedbyAkv0.Underverymildconditionsthe
sequenceαkconvergestothelargesteigenvalueinmodulusλ1ofAandthesequencevk
tothecorrespondingeigenvector.Thisholds,ifλ1isdominant(i.e.thereisoneandonly
oneeigenvalueoflargestmodulus)andsemi-simple,andifv0hasacomponentofthe
eigenvectorassociatedwiththiseigenvalue,see[86].Theconvergencerateisdictatedby
thequotient|λ2|/|λ1|,whereλ2isthesecondlargesteigenvalueinmodulus.
MoresophisticatedprojectionmethodscanbederivedbyemployingKrylovsubspaces
whicharedefinedby
Km=Km(A,v)=span{v,Av,A2v,...,Am−1v}
foragivenv∈Cn.Krylovsubspaceshavedesirablepropertiesinthecontextofprojection
methods,seee.g.[86,87].Forinstance,ifwechoosebothLandKtobetheKrylov
subspaceKm(A,v),theprojectedmatrixHmisanupperHessenbergmatrix,i.e.ithas
zeroentriesbelowthefirstsubdiagonal.Inthiscase,theresultingprojectionmethodis
knownasArnoldimethod.IfwefurtherassumeAtobeHermitian,i.e.AH=A,the
RayleighmatricesHmaretridiagonal.ThentheArnoldialgorithmsimplifiestotheso-
calledLanczosmethod.
Tospeedupconvergenceonemayemploypreconditioningtechniques.Theconver-
gencepropertiesoftheArnoldimethodaredirectlylinkedtotheseparationpropertiesof
theeigenvalues,see[86].Abetterseparationofthedesiredeigenvalueleadstoafaster
convergence.Therefore,onemayreformulate(5.1)toaneigenvalueproblemwithbetter

5.1.TheDavidsonMethod

Algorithm5.2BlockversionofDavidsonmethod.
1:functionDavidson(A,m,l)
2:ChooseaninitialorthonormalmatrixV1=[v1,...,vl];
3:fork=1,2,...do
4:ComputetheRayleighmatrixHk=VkHAVk;
5:Computetheldesiredeigenpairs(λ˜k,i,yk,i)1≤i≤lofHk;
6:ComputetheRitzvectorsu˜k,i=Vkyk,ifori=1,...,l;
7:Computetheresidualsrk,i=Au˜k,i−λ˜k,iu˜k,ifori=1,...,l;
8:ifconvergencethenexit;
9:Computenewdirectionstk,i=Ck,irk,ifori=1,...,l;
10:ifdim(Vk)≤mthen
11:Vk+1=MGS(Vk,tk,1,...,tk,l);
else12:13:Vk+1=MGS(u˜k,1,...,u˜k,l,tk,1,...,tk,l);
ifend14:orfend15:ctionunfend16:

53

separationproperties.Forinstance,onemaytransform(5.1)to(A−σI)−1x=µx,where
µisclosetothedesiredeigenvalueλ.TheeigenvectorsofAand(A−σI)−1areidentical
andtheeigenvaluesλofAcanberecoveredbytheidentityµ=(λ−σ)−1.Thisisthe
conceptoftheshift-and-invertArnoldimethod.
AditionedmoreversiongeneraloftheapproacArnoldihiscmethohosendin,buttheinmethoorderdoftoDagainvidson.moreItisflexibilitbasedyitonaalloprwsecon-the
preconditioningtovaryateachiteration.NotethatinthiscasetheRayleighmatricesHm
dense.ecomebymaeigenAvbloaluesckvλiersionofaofmatrixtheDaAisvidsongivenmethoindtoAlgorithmcompute5.2,lseedesired[30,86(for].insThetanceintegerthemsmallest)refers
totherestartparameteranddeterminesthemaximumsizeofthebasisVjfromwhichthe
eigenvectorsarebuilt.Ateachiterationk,wecomputetheRayleighprojectionHkofthe
matrixA.Then,instep5,theprojectedeigenvalueproblemissolved.Thisproblemisof
desiredmaximumeigensizevaluesm+landandtheconeorrespcanondingapplyforeigenvinstanceectors,thseeeQRe.g.[49].algorithmInstepto11dete(andrminestepthe13
respectively)newdirectionstj,i=Cj,irj,iareincorp−1oratedinthebasisVj.Wechoosethe
preconditioningmatricesCj,itobeoftheformMi,whereMiisanapproximationof
(A−λiI).Inthiscontext,theabbreviationMGSstandsforthemodifiedGram-Schmidt
cedure.proionorthogonalizatNoticethatbesidessomebasicvectorroutines,onlysparsematrix-vectormultiplications
Avforagivenvectorvareneeded.SincethepreconditioningmatricesCmayvaryat
eachiterationk,theseneednottobebuiltactually.Instead,onemayusekiterativ,iemethods
inordertodeterminetk,i,i.e.onesolvesMitk,i=rk,iinstep9iteratively.Furthermore,
notethattheRayleighmatricesHkcanbecomputedonlybyupdatingthelastlrowsand
columnsateachiterationk.
Algorithm5.2isvalidforbothsymmetricandunsymmetriceigenvalueproblems.The
onlydifferenceis,thatintheunsymmetriccase,complexeigenvaluesmayoccur.Never-
theless,acomplexarithmeticiseasilyavoidedbysplittingcomplexvectorsintotworeal
vectors:oneholdingtherealpart,theothertheimaginarypart.Inordertoincorporatea

54

EIGENVALUESOLVERSANDPARALLELIZATION

complexvectortothebasisVj,wejustappendbothrealandimaginarypartastworeal
vectors.Thisseparationisnorestriction.Rather,itinducesmoreflexibility,sincewith
morevectorsinV,therearemorewaystocombinethem.
TheDavidsonjmethodcaneasilybeextendedtosolvegeneralizedeigenvalueproblems
u.λM=AuIfMissymmetricpositivedefinite,onemayuseanM-orthogonalizationprocedurein
step11andstep13inAlgorithm5.2.ForgeneralM,weprojectbothmatricesAandM
andsolveageneralizedeigenvalueprobleminstep5(forinstancewiththeQZalgorithm,
see˜[49]).Inbothcases,theresidualcomputationinstep7isreplacedbyrk,i=Au˜k,i−
Mλk,iu˜k,i.
Convergenceresultsforthesymmetriccasecanbefound[30,70,86].Thenonsym-
metricso-calledcaseJacisobi-Davidsonexaminedin[metho69,d.88].ItaddsAnadvanancemenadditionaltofobliquetheDaprovidsonjectionmethinodstepis7theof
thanAlgorithmtheDa5.2,vidsonsee[92meth].oHodwinevaller,thecasesaspJacobi-Daointedvidsonoutin[metho72].disnotproventobebetter

5.2ComputationofPseudospectra
Inthisworkweconsiderpseudospectrawithrespecttothespectralnorm.Then,pseu-
dospectracanbedeterminedbymeansofsmallestsingularvaluesofshiftedmatrices,
cf.Section2.3.1.Moreprecisely,wecomputethenormoftheresolvent(z−A)−12=
smin(z−A)fordifferentzonagridinthecomplexplane.
costForandlargehighmatrices,storagearequircompleteements.singularTherefvore,aluedifferendecomptositioncomputationalinduceshighmethodshacomputationalvebeen
developed,see[102,103]andreferencestherein.
Theconceptofthepoorman’spseudospectrumreliesontheequivalentdefinitionbased
σε(A)=σ(A+E),
ondisturbingthematrixA
<εEseeDefinition2.2.Theresultingalgorithmcomputesspectraofthedisturbedmatrices
A+TheoremEfor2.3,someitisrandomsufficientmatricthatesEEiswithofrankE1<whicε.hAswreduceeshavtheeseencostinofthescalingproofthisof
].84[seematrix,Especiallyforsmallermatricesthecomputationofpseudospectracanbeperformedby
meansofspectraldichotomymethods,see[48,64].First,ε-spectralspots,i.e.regionsin
thecomplexplainwhichcontaintheε-pseudospectrum,aredetermined.Then,spectral
thisprojectorsmethodassoprcoiatedvidestosptheseectralreprogionsjectorsarecandomputed.invarianIntadditionsubspacestotheasεwell.-pseudospHowevectrum,er,it
entailshighercomputationalcostsandisthereforepreferredforsmall-scaledproblems.
Anotherapproach,inparticularforlargematrices,istoprojectthematrixAontoa
subspaceandcomputethedesiredpseudospectraoftheprojectedmatrix.LetUdenotea
thematrixsubspacewhoseisccolumnshosentoformbeaninvarianorthonormalt,itiseasybasistoofshothewthatsubspaceunderconsideration.If
σε(UHAU)⊆σ(A),(5.4)

5.2.ectraPseudospofComputation

Algorithm5.3Computationofaspectralportrait.
1:functionx2−Drax1w_pory2trait−y1(A,(x1,x2),(y1,y2),nx,ny)
2:hx=nx−1,hy=ny−1;
3:z=x1+iy1;
4:forj=1,...,nxdo
5:fork=1,...,nydo
6://ComputeStartswithmin(thez−prA)eviouswithcomputeAlgorithmd5.2singular;subspace
7:z=z+ihy;//Nextline
orfend8:9:z=z−ihy+hx;//Nextcolumn
10:hy=−hy;//Changesweepdirectionalongthecolumn
11:12:endendffunctionor

55

uessee[of103the].Inpro[101jecte]adpromatrixjectionUHonAUtoconKrylovergvetosubspacestheeigenisvpropaluesosed.ofATcloseypicallyto,thethebeigenoundaryval-
ofthespectrum.However,Krylovsubspacesarenotinvariantingeneralandhenceinclu-
UsionHAU(5.4)onedomaesynotshowhold.thatNevtheertheless,resultingbpyappseudospendingectraoneareroconwtotainedtheinσHessenε(Ab),ergsee[10matrix1].
In(z−thisA)−case,1withtheε(z−A-pseudosp)−1ectrumdenotingoftheapseudoirectangunvlarersematrofix(zis−A).definedastheε-curveof
eigenInvthisaluewsolvorkers.wefoIncusthisonconthetextthecomputationDavidsonofmetpseudosphodwectasrasucbycessfullymeansofused,iterativsee[e55,sparse80].
speFctrolloalpwortringaittheofadefinitionmatrixAofbythepseudospplotectraoftheintromapducedinSection2.3.1,wedefinethe
z→sp(A)(z)=log10(z−A)−12=−log10[smin(z−A)],
wheresmin(z−A)denotesthesmallestsingularvalueofz−A.
ThecomputationofaspectralportraitofamatrixAisoutlinedinAlgorithm5.3.
Pleasenotethathereireferstotheimaginaryunit.First,wedefineagridinthedomain
ofinterestinthecomplexplane.Thegridsweconsideraregiveninarectangulardomain
[x1,x2]×[y1,y2]withnx×nypoints(nxinthehorizontal,Hnyinthev1/ert2ical).Then,
fortheanDaygridvidsonpmointethozwd.eFortcomputewonesighminb(zoring−A)grid=pλoinmints(zz1−andA)z2(,zw−eAexp)ectbtheymatrimeanscesof
(z1−A)and(z2−A)tohaveclosesingularvaluesandclosesingularvectors.Toimprove
performance,westarttheDavidsonalgorithmusingthesingularsubspacecomputedin
thecomputelast(zstep−Aas)Hin(zitial−A)guess,explicitlysee.lineOnly6inmatrixAlgorithmvector5.3pro.ductsNoteofthattheweformdoAvnotandneedAHtov
foragivenvectorvareneeded.Moreover,ifthematrixAisreal,thespectralportraitof
Aisclearlysymmetricwithrespecttotherealaxis.
origin.InhydroSincewedynamicexpectstabilittheyeigenwevarealuesofmainlyAtoinbeinterestedthisinregitheon,sptheectralpsmallestortraitsingulararoundvaluethe
inofzthe−ADaisvidsonlikelymethoclosedtotozero.beoftheTherefore,formwMe−c1,hoosewhereallMtheisaappropreconditioningximationofmatricAHesA.CkIt,i
iseasytochooseMtobepositivedefinitewhichisadesiredpropertyintheconvergence

56

EIGENVALUESOLVERSANDPARALLELIZATION

statementsoftheDavidsonmethod,see[23,80].
Foramatrixpencil(A,M),wedefinethespectralportraitas
z→gsp(A,M)(z)=log10(zM−A)−12=−log10[smin(zM−A)],
cf.Section2.3.4.Inthiscase,thecomputationofaspectralportraitisgivenbyAlgo-
rithm5.3,wherewehaveonlytoreplaceline6byapplyingtheDavidsonmethodon
(zM−A).
arallelizationP5.3Theproblemswehavetocopewithareusuallyofhighorder(∼105−107)andconsequently
entailhighcomputationalcosts.Tosolvethesekindofproblemsinreasonabletime,
parallelizationtechniquesarevirtuallymandatory.
Inthissectionwetreatparallelizationschemesonplatformswithdistributedmemory.
Sincethemosttimeconsumingpartiniterativesolversisthesparsematrixvectormulti-
plication,wefocusonthisissue.Wealsoshowhowthiscanbecombinedwithasecond
levelofparallelisminordertomanagetheveryexpensivecomputationofpseudospectra.

5.3.1SparseMatrixVectorMultiplicationforFiniteElementMethods
Inourframeworkassembledmatricesandvectorsaredistributedrow-wiseoverthepro-
cessesaccordingtothedistributionofthedegreesoffreedom(DoF)fromthefiniteelement
approach.Eachlocalsub-matrixisdividedintotwoblocks:adiagonalblockrepresenting
allcouplingsandinteractionswithinthesubdomain,andanoff-diagonalblockrepresenting
thecouplingsacrosssubdomaininterfaces.
InFigure5.1wefindadomainpartitioningintofoursubdomains.Inordertodetermine
thestructureoftheglobalsystemmatrices,first,eachsubdomainhastobeassociatedwith
asingleprocess.Then,eachprocessnotonlydetectscouplingswithinitsownsubdomain,
butalsocouplingstotheso-calledghostDoF,i.e.neighboringDoFwhichareownedbya
differentprocess.Intermsofaglobalsystemmatrix,DoFiandjinteractifthematrix
hasanon-zeroelementinrowiandcolumnj.
ThedistributedsparsematrixvectormultiplicationisgiveninAlgorithm5.4.While
everyprocessiscomputingitslocalcontributionofthematrixvectormultiplication,an
asynchronouscommunicationforexchangingtheghostvaluesisinitiated.Afterthiscom-
municationphasehasbeencompleted,thelocalcontributionsfromthecoupledghostDoF
areaddedaccordingly,see[6,58].
Forthecomputationofpseudospectra,weneedmultiplicationsofatransposedmatrix
withvectorsaswell.However,usingarow-wisedistributionresultsinanalltoallcommu-
nicationofthevectortobemultiplied,see[55].Thus,inordertoreducecommunication,
westorethematrixtwice:oncewitharow-wisedistributionandoncewithacolumn-wise
distribution.

5.3.2ParallelComputationofPseudospectra
Fordifferentgridpointsz,thecomputationofsingularvaluesofz−Acanbeper-
formedcompletelyindependentlywhichallowsaneasywaytoparallelizetheevaluations
of(z−A)−12.However,thisimpliesthatacopyofthematrixAisavailabletoeach

zationaralleliP5.3.

P0

P1

P2

P3

57

Figure5.1:Domainpartitioning:DoFofprocessP0aremarkedingreen(interiorDoFin
diagonalblock);theremainingDoFrepresentinter-processcouplingsforprocessP0(ghost
DoFinoff-diagonalblock).

Algorithm5.4Distributedmatrixvectormultiplication.
functionStartasynchrondistr_mvmuloustcomm(A,x,yunication:)
alues;vghosthangeExcyinSynct=Ahronizediagxcoint;mmunication;
yint=yint+Aoffdiagxghost;
functionend

••
••
P0•+P1P2P3•
••
diagonalblockinterioroffdiagonalblock
•ghost
Figure5.2:Distributedmatrixvectormultiplication.

58

EIGENVALUESOLVERSANDPARALLELIZATION

process.Inordertoavoidhighstoragecostsinthecaseoflargematrices,weutilizea
parallellinearalgebraaswell.
Weuseanapproachwhichisreferredtoashybridparallelism.Therefore,wepartition
themappeddomaintoeacofhgridofptheseointssubindoCmainintoofkgsubriddpoinomainstsofbuildingtheksamegroupssize.,proThenvidedpprothatcwessesehavaree
k×amongpprothecessesproacessvesailable.inorderWithintopeacherformgroupapthearallelsystemcomputationmatricesofandthevectorssmallestaresingularspread
vector,seeFigure5.3.

Ω

ΩΩ30

Im
Ω


Ω0Ω3
Re


Ω1Ω2


G0G1G2
Figure5.3:TheleftfigureshowsadecompositionofthephysicaldomainΩintofour
subdomains(p=4).Therightonedepictsthedistributionofgridpointsin[x1,x2]×
[y1,y2]⊂CamongthreegroupsGi(k=3).

ΩΩ21

6ChapterPseudospectraforApplicationsin
yStabilitdynamicHydroThesucceedingnumericalcomputationsofspectralportraitswerecarriedoutwiththe
3casewpredecessorehaveofusedthethefiniteparallelelementdirectsoftwsparsearelinearHiFlowsolv,ersee[1,MUMPS6].Inasathetwpreconditiono-dimensionalerfor
theDavidsonmethod.Thissolverisbasedonamultifrontalmethod,see[5].Although
thismethodmaynotbethebestchoiceintermsofscalability,wehaveachievedthe
bestperformanceresultsfortwo-dimensionalproblems.Forthree-dimensionalproblems,
iterativesolvers,suchastheCGmethod,havebeenturnedouttobesuccessfulaswell.
Theparallellinearalgebraemployedinthisworkwasimplementedbymeansofthelibrary
PETSc,see[10,11,12].
CentrTheenforumericalComputingbenc(SCC),hmarkswKarereplsruhe,erformedGermanony.theEacHPhofXtheC3000coresclusterontheattheCPUsSteinbuch(Quad-
coreIntelXeonprocessors)runsataclockspeedof2.53GHz.
6.1IncompressibleFlowoveraBackwardFacingStep
Weconsiderastationaryincompressiblefluidflowoverabackwardfacingstepasdepicted
inFigure6.1.Thissetupisoriginatedfromawell-knownoptimizationproblemwherethe
vortexbehindthestepistobereduced,seee.g.[27].ThefluidflowinthedomainΩis
governedbythesteady-stateNavier-Stokesequations
−R1eΔv+(v∙)v+p=0,(6.1)
∙v=0,
whereacteristicwevhaelovecitysetVρ=and1.cHereharacteristicRe=VlengthL/νL.denotesAsantheinflowReynoldsbnoundaryumbercondwithitionchar-we
setv=vinonΓin,
wherevinhasaparabolicprofile.OnΓrigidweimposeno-slipboundaryconditions,i.e.
v=0onΓrigid.
59

60

z

y

Γin21

x

PSEUDOSPECTRAINHYDRODYNAMICSTABILITY

Γgirid

Γidigr

20

Γdiigr14

Figure6.1:Geometryofthebackwardfacingstepbenchmark.

Γout

Figure6.2:Stationaryflowinthetwo-dimensionalbackwardfacingstepgeometrywith
Re=100(upper)andRe=500(lower).

Furthermore,weprescribefree-streamoutflowconditions(ordo-nothingconditions)
1Re∂nv−pn=0onΓout,
wherenreferstotheoutwardunitnormal.
InournumericalexampleswehavechosenVtobethemaximumvelocityofthe
parabolicinflowvin,andLtobethegapheight.StationarysolutionsforRe=100
andRe=500ofthetwo-dimensionalproblemareplottedinFigure6.2.
Westartwiththevariationalformulationofthesteadyproblem.Weset
H(Ω)={v∈(H1(Ω))d:v=0onΓin∪Γrigid}
ford=2ord=3.Moreover,letvinext∈(H1(Ω))dbeasolenoidalextensionoftheinflow
vinsuchthatvinext=0onΓrigid,seee.g.[47].Then,theweakformulationseeksforpairs
(V,P)∈H(Ω)×L2(Ω)+{vinext,0}suchthat
R1e(V,ϕ)0+((V∙)V,ϕ)0−(P,∙ϕ)0+(∙V,ψ)=0
holdsforany(ϕ,ψ)∈H(Ω)×L2(Ω).Wesolvethisequationnumericallywithacon-
formingGalerkinfiniteelementmethodbyreplacingthespacesH(Ω)andL2(Ω)byfinite-
dimensionalsubspacesHhandLh.Hence,wehavetofind{Vh,Ph}∈Hh×Lh+{vin,exth,0}

6.1.IncompressibleFlowoveraBackwardFacingStep61
satisfyingR1e(Vh,ϕh)0+((Vh∙)Vh,ϕh)0−(Ph,∙ϕh)0+(∙Vh,ψh)=0(6.2)
forany(ϕ,ψ)∈Hh×Lh⊂H(Ω)×L2(Ω).Herevin,exthisanadequateboundaryinterpolant
ofvinext,see[91].
WechooseatriangulationThofΩconsistingofquadrilaterals(ford=2)orhexahedrons
the(ford=reference3)suchelementhattKˆthe:=(0,obtained1)d.Thismeshismeansaffinefor,i.e.anyeacKh∈KT∈Ttherehisaexistsffineanequivaffinealentandto
horientationpreservingmappingFKsuchthatK=FK(Kˆ),seee.g.[41].
Weconsiderso-calledQkfiniteelementswhichconsistofcontinuous,piecewisepolyno-
mialfunctionsofdegreek,seee.g.[41].Let
Qk(Kˆ)=span{xiyj:0≤i,j≤k},ford=2,
Qk(Kˆ)=span{xiyjzl:0≤i,j,l≤k},ford=3,
bethepolynomialsuptotheorderkonthereferencecellKˆ.Then,thevectorspaceof
Qkelementsisgivenby
Xhk:={uh∈C(Ω):uh|K◦FK∈Qk(Kˆ)foranyK∈Th}.
InournumericalexperimentswehavechosenQelementsforthevelocityandQ
elementsforthepressure,i.e.Hh={vh∈(Xh2)d:vh2=0onΓin∪Γrigid}andLh=Xh11,
inordertosatisfytheinf-sup-condition
|(∙ϕ,ψ)0|
ψ∈infLhϕ∈supHhϕH1(Ω)ψL2(Ω)>0.
Forpairsuh={vh,ph}andξh={ϕh,ψh}∈Hh×Lhwedefinethesesquilinearforms
1a(uh,ξh)=Re(vh,ϕh)0+((vh∙)Vh,ϕh)0+((Vh∙)vh,ϕh)0
−(ph,∙ϕh)0+(∙vh,ψh)0,
m(uh,ξh)=(vh,ϕh)0,
aswheredescribwehaedveinderivSecteionda(2.2∙,.∙)Letjustξbhiy(i=1,linearization...,n)ofbe(a6.2)basisaroundoftheasteadyfinite-dimensionasolution{Vlh,Pspaceh}
Hh×Lh.WedenotebyAhthestiffnessmatrix
Ah=a(ξhj,ξhi)1≤i,j≤n,(6.3)
andbyMhthesymmetricpositivesemidefinitemassmatrix
ijMh=m(ξh,ξh)1≤i,j≤n.(6.4)
Thesmallestsingularvaluesςharedeterminedbyfindingthesmallesteigenvaluesςh2and

62

PSEUDOSPECTRAINHYDRODYNAMICSTABILITY

Algorithm6.1SolutionprocedureinHiFlow.
1:Readinitialmeshandperformh-refinement;
3:2:SolvDefinee(6.2FE)Mbyonthmeansecellsofa(e.g.NewtonQ1orQmetho2);d;
4:5:forAssemµonbleaAhgridandinCMhdo;
6:Computesmallestsingularvalueςhbymeansof(6.5)withtheDavidsonmethod;
forend7:8:Plotspectralportrait;

thecorrespondingeigenvectorsxhof
(µMh−Ah)H(µMh−Ah)xh=ςh2MhHMhxh,(6.5)
forsomeµonagridinthecomplexplane.Thecompletesolutionprocedureisoutlinedin
.6.1AlgorithmSpectralportraitsofthetwo-dimensionalproblemareplottedinFigure6.3withthe
Reynoldsnumberrangingfrom100to500asindicated.Here,wecomputedsingularvalues
for4symmetry,488gridalongpointhetsrealintheaxis.Thecomplex259,971domainunkno[−0.6wns,0.w2]ere×[0,0.distributed4]⊂Candamongexp4coresloitedandthe
theminutesoncomputation4×8corestimeonfortheeachHPofXtheC3000depictedclusterspdescriectralbedportabovraitse.Ftoorokaroincreasingund110−Reynolds135
numbersthepseudospectraprotrudemoreandmoreintotherighthalfofthecomplexplane
whichmayindicateinstability.
Inordertoteststrongscalabilityofourcomputationswehavesetupatwo-dimensional
b672enchmarksingularvproblemaluesonchoaosinggridRine[=−0.4100,0..2]This×[0b,0enc.4]⊂hmarkC.Theconsists259of,971theunknocomputationwnswereof
Sectiondistributed5.3.2.amongInT2iable,i=6.12,th..e.,8resultscoresintermsexploitingofthestronghybridscalabilityparallelismusingasonedescribgroupedarein
thegiven.unknoInTwnsableov6.2erw4eprofindthecesses,corrweesphaveondingobtainedresultsinthebcaseestofptwoerformance,groups.ByseeTdistriablebuting6.3.
Here,wespreadthegridpoints,dividedinverticalstripes,equallyamongeachgroupas
inachievFigureeda5.3sig.nificanAltoughtimprothisvseemsement,tosebeeaFigurequite6.4si.mpleloadbalancingtechnique,wehave
areplottedSolutionsinofFigurethestead6.5.yAsfloinwinthetthewo-dimensionalthree-dimensionalcase,casethevforortexReb=ehind100andtheRbace=kward500
facingstepenlargesastheReynoldsnumberincreases.Thistendencytoinstabilityis
ofreflectedthebycomplextheplane,pseudospseeectra,Figurewhic6.7h.areHere,protrudingthemdiscretizationoreandofmoretheintoproblemtherightyieldedhalf
as143,484intheunktwnowns,o-dimensiwhiconalhwcase,erei.e.spreadweaplottmonged8cores.singularWvealueshaveforchosen4,488thegridsamepointssetupin
[−0.6,0.2]×[0,0.4]⊂C.
tionTheoftheestimateperturb(e2.30d)inquantitTheoybyremmeans2.13ofstatesapseudospminimectra:umFgroorwthanyafactor=Reforzthe>0evoandlu-

6.1.IncompressibleFlowoveraBackwardFacingStep

200=eR(b)

400=eR(d)

(a)100=eR

300=eR(c)

500=eR(e)

63

Figure6.3:Spectralportraitsinthedomain[−0.6,0.2]×[−0.4,0.4]ofthetwo-dimensional
problemfordifferentReynoldsnumberswithcontourplotsfor{−7,−6,...,−1}.

64

PSEUDOSPECTRAINHYDRODYNAMICSTABILITY

groupsizenumbercorestime(s)speedupefficiency
1.1.665344168168426651851.561.280.390.64
0.292.30288932320.152.3927886464

Table6.1:Scalabilitytestusingonegroup.

groupsizenumbercorestime(s)speedupefficiency
1.1.4313841683216251331621.721.360.430.68
0.322.53170564320.172.65162812864

Table6.2:Scalabilitytestusingtwogroups.

numbercorestime(s)speedupefficiency
1.1.66534816431324471.542.720.680.77
643269113015.119.630.640.60
25612817735837.5018.570.590.58

Table6.3:Scalabilitytestwithgroupsizefour.

Figure6.4:Runningtimesfordifferentgroupsetups.

6.1.IncompressibleFlowoveraBackwardFacingStep

65

ReFigure=1006.5:(upper)StationaryandRfloew=in500the(lowthree-er).dimensionalbackwardfacingstepgeometrywith

Figure6.6:Plotoflowerbound(6.6)withrespecttothethree-dimensionalproblemwith
Reynoldsnumbersrangingfrom100to500asindicated.

66

PSEUDOSPECTRAINHYDRODYNAMICSTABILITY

200=eR(b)

400=eR(d)

R(a)100=e

(c)300=eR

500=eR(e)

Figure6.7:Spectralportraitsinthedomain[−0.6,0.2]×[−0.4,0.4]ofthethree-dimensional
problemfordifferentReynoldsnumberswithcontourplotsfor{−7,−6,...,−1}.

6.2.NaturalConvectioninaHorizontalAnnulus

67

L=a(z−A)−1wehave
aτsupetA≥eτa1+e−1.(6.6)
Lτ≤<t0Forthethree-dimensionalproblemwehaveplottedthisboundinFigure6.6,wherewehave
chosenawithinthecomplexdomainshowninFigure6.7suchthatLattainsitsmaximum
value.TheplotshowsthatforincreasingReynoldsnumbers,thelowerboundincreases
aswell.Thiscoincideswiththeobservationthatasteadyflowtendstoinstabilityfor
increasingReynoldsnumbers.

6.2NaturalConvectioninaHorizontalAnnulus
InthissectionweconsidernaturalconvectiveflowsofNewtonianfluidsinaninfinitelylong
gapbetweentwohorizontalcoaxialcylinders.Thisphysicalprocessappearsinavarietyof
applications,suchasdesigningenergystoragesystems,coolingofelectroniccomponents,
ormodelingaircraftcabininsulations.
WedenotetheradiusoftheinnercylinderbyRi,andbyRotheradiusoftheouter
cylinder(thus0<Ri<R0).Moreover,weprescribeconstantsurfacetemperaturesTi
ontheinnerandToontheouterjacketwithTo<Ti.Modelingthissystemwiththe
Oberbeck-Boussinseqequations(2.8),theflowpatternischaracterizedbythreedimen-
sionlessquantities.Wehavethegeometricparameter
R2iA=Ro−Ri
whichisknownasinverserelativegapwidth.Moreover,RadenotestheRayleighnumber
ybdefinedRa=αg(Ti−To)(Ro−Ri)3,
κνwiththevolumetricexpansioncoefficientα,gravityaccelerationg,kinematicviscosityν,
andthermaldiffusityκ.Finally,thePrandtlnumberisgivenby
ν.=rPκDependingonthesecharacteristics,theproblemisverycomplexandallowsverydifferent
possiblebehavioroftheflow.Uptonow,itwasmostlyinvestigatedbynumericalexper-
iments,seee.g.[118,119]andreferencestherein,butthereexistsometheoreticalresults
aswell,see[77,78]andreferencestherein.
Despitethefactthatspiralflowsmayoccur,weconfineourselvestoatwo-dimensional
model,whichisphysicallyreasonableforlowRayleighnumbersaspointedoutin[39].
Hence,thedomainunderconsiderationisgivenby
Ω={(r,ϕ):Ri<r<Ro,0≤ϕ<2π},
seeFigure6.8.
Asdescribedin[117],iftheRayleighnumberexceedsacriticalvalue,twosteadysolu-
tionsareobserved:Thedownwardflowconsistsoftwocounterrotatingvorticesineach
halfoftheannulus,whereastheupwardflowhasacrescentshapededdyoneachside,see

68

PSEUDOSPECTRAINHYDRODYNAMICSTABILITY

.6.9FigureThegoverningsteady-stateOberbeck-Boussinseqequationsread
−νΔv+(v∙)v+1p=g[1+α(T0−T)],
ρ∙v=0,
∂tT+v∙T−κΔT=0,
seeSection2.1.6.Weimposeno-slipboundaryconditionsforthevelocity,i.e.
v=0on∂Ω.
Moreover,wesetTi=1andTo=0whichresultsinthefollowingboundaryconditionsfor
erature:tempthe

T(x)=1forx2=Ri,
T(x)=0forx2=Ro.
ForthediscretizationwechooseQelementsforthevelocityandtemperatureandQ
elementsforthepressure.Thus,wit2hthenotationsHh={vh∈(X2)2:vh=0on∂Ω}1,
Lh={ph∈Xh1:Ωph=0}andWh={Th∈Xh2:Th=0onh∂Ω}andanappro-
priateboundaryinterpolantTextoftheimposedtemperatureontheinlet,weseekfor
(Vh,Ph,τh)∈Hh×Lh×Wh+i,h{0,0,Ti,hext}suchthat
ν(Vh,ϕh)0+((Vh∙)Vh,ϕh)0−(Ph,∙ϕh)0+gα(τh,ϕh)0
+(∙Vh,ψh)0(6.7)
+(Vh∙τh,ζ)0+κ(τh,ζh)0
=αTo(g,ϕh)0
holdsforany(ϕh,ψh,ζh)∈Hh×Lh×Wh.Forbrevity,wehavesetρ=1.Bylin-
earizationaroundasteadysolutionasintheprecedingSection6.1,weobtainfortriples
uh={vh,ph,Th}andξh={ϕh,ψh,ζh}∈(Xh2)2×Xh1×Xh2thesesquilinearform
a(uh,ξh)=ν(vh,ϕh)0+((vh∙)Vh,ϕh)0+((Vh∙)vh,ϕh)0
−(ph,∙ϕh)0+gα(Th,ϕh)0+(∙vh,ψh)0
+(vh∙τh,ζh)0+(Vh∙Th,ζh)0+κ(Th,ζh)0.
Moreover,wedefine
m(uh,ξh)=(vh,ϕh)0+(Th,ζh)0.
Then,usingthedefinitionsforthestiffnessmatrixandthemassmatrixasin(6.3)and
(6.4),weobtainthesameformulationasbefore:Weseekthesmallesteigenvalueς2ofthe
generalizedeigenvalueproblem
(µMh−Ah)H(µMh−Ah)xh=ςh2MhHMhxh.
inSpFigureectral6.10por,traitswherewwithehavrespesetecttAo=the2,Prsteady=0.7,solutionsRa=5de,000picted.ForineacFighureoft6.9heseareplotsplottedwe

6.2.NaturalConvectioninaHorizontalAnnulus

z

RorRϕi

x

Figure6.8:Two-dimensionalmodeloftheflowregion.

69

Figure6.9:Streamlinesandtemperaturedistributionofthedownwardflow(left)and
upwardflow(right)forA=2,Pr=0.7,Ra=5,000.

Figure6.10:Spectralportraitsofthedownwardflow(left)andtheupwardflow(right)in
thedomain[−0.4,0.2]×[−0.4,0.4]withcontourplotsfor{−4,−3.5,...,−1}.

70

PSEUDOSPECTRAINHYDRODYNAMICSTABILITY

havechosenadiscretizationresultingin429,568unknownsandevaluated2,624singular
valuesexploitingthesymmetryalongtherealaxis.
Thepseudospectraofthedownwardflowareprotrudingmoreintotherighthalfof
dothewnwcomplexardflowplanetobethancloserthetopseudospinstabilitectrayofthanthetheupwupardwardflow,flow.whicHohwmaevyer,iasndicateforthethe
beBénardinstableproblem,andtheiftdowownwsteadyardflowtosolutionsbearestable.Ifrealized,theitisReynoldsexpecntedumbtheerisupwardincreasedflowantoy
further,thedownwardflowisobservedtobecomeinstableaswell.

7Chapter

ReactorNuclearinApplicationsTheory

theUpktonoeigenw,vthealuepcowerriticalitmethoyd(problemcf.(Algori3.4)thmin5.1react)orhasbeenapplications,themethoseed[66of].cInhoiceorderforstoolvingfind
abetteralternativetothisratherbasicmethod,severalwayshavebeeninvestigated:
TheTheemploJacobi-DaymenvidtofsonthemethoimplicitdinthisrestartedcontextisArnoldidiscussedmethodinis[52,110examined],and,in[22more,111recen,112tly].,
theThesetecapplicationhniquesofhathevebeenJacobian-freeappliedtotheNewton-Krylodiffusiovnmethoequatiodnis(i.e.Pexamined1approin[45,ximation46,60for].
thesphericaldependence)ortothediscreteordinatetransportapproximation(i.e.SN
far.soximation)approInthischaptertheDavidsonmethod(cf.Algorithm5.2)iscomparedtothetraditional
pareowpererformedmethodinwithinthetheframewframeworkoforktheofPtheNsoftapprowareximation.Parafish,Theseen[107umerical].Thisevaluationsparallel
solverisbasedonanon-overlappingdomain-decomposition(DD)technique.Itusesthe
mspatialultigroupdiscretization,approachfoandrthesphericaenergylharmonicsdiscretization,fortheangularnon-conformingdiscrfiniteetizationelemenastsdescribfortheed
.4.2ectionSinfinallyThetheorderangle.ofAtdiscretizationtheenergyinPlevel,arafishweisuseasfaolloblocws:kvefirstrsiontheofentheergy,Gathenuss-Seidelthespace,methoandd
nogivenupscthebloatteringckiswisepresenrepresent,i.e.σtationsg→gof≡th0eformgulti<ggroup,blockdiscretization,triangularitseeyisSectionfullyac4.2.1hiev.edIf
andonlyoneGauss-Seideliterationissufficient.
Atthespatio-angularlevel,domain-decompositionisappliedyieldinganinteriorblock
seeshapeSectionsuchthat4.2.2.eachinTherefore,teriorthesediagonaldiablogonalckisblockssymmetriccanbeandctreatedorrespindepondstoendenonetlybydomain,dif-
vferenersiontprofothecesses.AGMREStthismethodspatio-angularbyfaclevtorizingel,weeachapplyinateriorblock-diadiagonalgonallyblockprebyanconditionedincom-
pletelibrariesLUSp(ILU)arseLib++[factorization.36]andTheSuperLUILU[33decomp].TheositiondifferenwastiterationsimplemenlotedopsbyofmeaPnsarafishoftheare
schematicallydisplayedinFigure7.1.
domain.MoreovIner,thisPwaarafishy,oneiscandesignedtakeinfullsucadvhaanwatageyofthattheoneproflexibilitcessyofcanthehandleresultingmoreprethancondi-one
tioningwithoutincreasingthenumberofprocesses.Furthermore,Parafishallowstheuser

71

72

Energy

Space

Angle

vidsonDaorerwoPiterations)(outer

APPLICATIONSINNUCLEARREACTORTHEORY

iterations)Gauss-Seidel(thermal

iterations)GMRESDomainositiondecomp(inner

ositiondecompILU

(wholeGlobaldomain)(eachLosubcaldomain)
Figure7.1:Parafishcalculationscheme,see[107].

tousedifferentfactorizationmethodsindifferentdomains.Notethattheparallelization
schemeusedinParafishdiffersfromtheonepresentedinSection5.3.Theinterfacenodes
areduplicatedamongneighboringdomains,see[107].

7.1TheTakeda1Benchmark
Thethree-dimensionalTakeda1benchmark[100]consistsofasmalllightwaterreactor
(LWR)coreofcubicshapeasdepictedinFigure7.2.Allunitsoflengtharegivenincm.
ThiscorrespondstothesimplifiedmodelofareactorasshowninFigure3.1withonecore
andtwocontrolrods,andassumingsymmetryalongtheplanesx=0,y=0andz=0
(ifwesettheorigininthecenterofthereactor).TheTakeda1benchmarkconsidersonly
oneoctantofthereactordescribingtheremainingonesbyreflectingboundaryconditions.
Hence,itprescribesreflectiveboundaryconditionsforx=0,y=0andz=0.Onthe
outerframe,i.e.forx=25,y=25andz=25,weimposevacuumboundaryconditions.
Furthermore,wehavetwoenergygroups,wherenoupscatteringispresent,i.e.σs2→1≡0.
Thecorrespondingcrosssectionvaluescanbefoundin[100].
Inthegeometrythreedifferentzonesaretobefound:thecorezone,themoderator
zoneandtheguidetubezone.Weconsiderbothinstancesofthisbenchmark:
•Case1:Thecontrolrodisintheguidetubezone.
•Case2:Thecontrolrodisoutofthecore,i.e.theguidetubezoneisvoided.
FortheangulardiscretizationweusethePNmethod,withNrangingfrom3to7as
indicated.Thespatialdependenceisdiscretizedbymeansofa25×25×25meshwith

7.1.TheTakeda1Benchmark

youndarybacuumv25

y25

acuumvyoundarb

73

215oundarybreflection5voundarybacuumoundarybreflection5oundarybacuumv
2133000152025x0152025x
reflectionboundaryreflectionboundary
(a)0<z<15(b)15<z<25

Figure7.2:Geometryofthe3DTakeda1benchmark.Zone1=core,zone2=moderator,
zone3=guidetube.

hexahedralNC6elements.Thiselementisathree-dimensionalgeneralizationofthetwo-
dimensionalrotatedQ1element,see[83].ItscapabilitytoaccuratelysolvetheTakeda1
benchmarkhasbeenassessedpreviouslyin[106].Thefiniteelementnodesaredistributed
among5×5×5domains,withdomaininterfacesparalleltothecubefaces.Notethat
thisdomaindecompositionimpliestheguidetubezonetobeentirelycomprisedwithin
domain.singleoneIncase1ofthebenchmark,wehaveoptimizedourcomputationaltimebyapplyingan
ILU(0)factorizationinalldomainsasapreconditionerfortheGMRESmethod.Incase2,
wehaveachievedthebestresultsbyapplyingavariantofanILUTP(incompleteLU
factorizationwiththresholdandpivoting)methodasdescribedin[67]inthevoidedguide
tubezoneandanILU(0)decompositioninallremainingdomains.Theresultsweobtained
in[96,97]aregivenintheTables7.1,7.2and7.3.Thesecomputationswereperformed
ontheHPXC3000clusterasdescribedinChapter6.Pleasenotethatwelistedonlythe
timesforsolvingtheeigenvalueproblemincludingthetimefortheILUdecomposition.
Weexcludedthetimeforassemblingthematrices.
Forcase1,weobtainedfortheP3(P5andP7)approximationthecriticaleigenvalue
keff=0.96122(0.96222and0.96241respectively)whichagreesdecentlywiththereference
valuesof[99],e.g.keff=0.9623±0.00048withtheMonte-Carlomethod.Satisfying
resultshavealsobeenachievedincase2,whereournumericalexperimentsyieldedforthe
P3(P5andP7)approximationthecriticaleigenvaluekeff=0.97247(0.97553and0.97633
respectively)whichiscomparablewiththereferencevalue0.9778±0.00046obtainedwith
theMonte-Carlomethodin[99].
Thecomputationofthecriticaleigenvalueandthecorrespondingeigenfunctionusing
theDavidsonmethodconsumesroughlyhalfofthetimewhichisneededbyapplyingthe
powermethodtogaincomparableresults.Noimprovementtothepowermethodresults
couldbemadebyutilizingaChebyshevacceleration.Accordingtothenumberofouter

74

APPLICATIONSINNUCLEARREACTORTHEORY

indro1:CaseDavidsonMethodPowerMethod
NumberCoresTime(s)SpeedupEfficiencyTime(s)SpeedupEfficiency
515.621.83.881.0.781.11.241.03.661.0.731.
251251.61.114.0320.020.560.163.11.233.2013.140.270.53

outdro2:CaseDavidsonMethodPowerMethod
NumberCoresTime(s)SpeedupEfficiencyTime(s)SpeedupEfficiency
132.71.1.58.01.1.
58.43.90.7815.43.80.75
253.39.80.395.810.00.40
1252.215.10.123.616.00.13
Table7.1:ResultswithP3approximation(675,000unknowns).

indro1:CaseDavidsonMethodPowerMethod
NumberCoresTime(s)SpeedupEfficiencyTime(s)SpeedupEfficiency
1598.524.11.4.080.821.46.67182.83.921.0.781.
125252.46.941.2514.230.330.574.012.115.1045.530.600.36

outdro2:CaseDavidsonMethodPowerMethod
NumberCoresTime(s)SpeedupEfficiencyTime(s)SpeedupEfficiency
1152.71.1.307.21.1.
538.24.10.8078.83.90.78
2517.38.80.3530.510.10.40
12512.612.10.1020.415.00.12
Table7.2:ResultswithP5approximation(1,687,500unknowns).

7.2.TheNEAC5G7Benchmark

indro1:CaseDavidsonMethodPowerMethod
NumberCoresTime(s)SpeedupEfficiencyTime(s)SpeedupEfficiency
15239.869.83.431.0.691.141.0479.73.401.0.681.
125255.019.247.5712.460.380.5035.19.513.6850.670.550.41

outdro2:CaseDavidsonMethodPowerMethod
NumberCoresTime(s)SpeedupEfficiencyTime(s)SpeedupEfficiency
1633.11.1.1117.61.1.
5161.63.90.78268.44.20.83
2590.77.00.28139.08.00.32
12577.88.10.07112.010.00.08
Table7.3:ResultswithP7approximation(3,150,000unknowns).

75

iterations,wehaveobtainedsimilarresultsasshowninTable7.4.Withregardstothe
speedup,boththeDavidsonandthepowermethodshowacomparablebehavior.Thisis
quiteevidentsinceinbothmethodsthesamebasiclinearalgebraroutinesareutilized.
TheDavidsonmethodappearstobearobustmethodtocomputeeigenvaluesandthe
correspondingeigenvectorsofageneralizedeigenvalueproblemsincetheresidualsofthe
formFv−λHvarecomputedaccurately(inexactarithmetic).Moreover,theDavidson
methodbuildsafullbasisofanadequateprojectionspace,whereasthepowermethod
projectsonlyontothespacespannedbythelastbuiltvectoroftheform(H−1F)jv.Hence,
theDavidsonmethodkeepstrackofmoreinformationandthereforeneedslessiterations.
BothmethodsrequiresolvinglinearsystemsoftheformHx=y.TheDavidsonmethod
turnedouttobemorerobustinpracticeinthesensethatthechoiceofthepreconditioner
islessrestrictivethaninthepowermethod.Thismeans,solvinglinearsystemsoftheform
Hx=ycanbedoneveryroughly,whilethepowermethodneedsamoreaccuratesolving
inordertoconverge.
ComparingthetwocasesoftheTakeda1benchmark,wenoticethattheperformance
andscalabilityachievedincase2(rodout)cannotkeepupwiththeresultsobtainedincase
1(rodin).Havingtherodoutofthecoreisthetoughercasebecauseofthelow-density
regionappearinginthevoidedguidetube.Asnoticedbefore,thevoidedguidetubeis
comprisedwithinonedomain,whichimpliesthatthecorrespondingworkloadwithinthis
domainislargercomparedtothecasewheretherodisinthecore.

7.2TheNEAC5G7Benchmark
Inthissectionweconsiderthetwo-dimensionalMOXfuelassemblybenchmarkissuedby
theNuclearEnergyAgency(NEA).ItusesC5MOXfueland7energygroups,henceits
“C5G7”nickname.AsdepictedinFigure7.3(again,allunitsoflengthareincm),its

76

APPLICATIONSINNUCLEARREACTORTHEORY

indro1:CaseDavidsonMethodPowerMethod
AngularNumberNumberNumberNumber
ApproximationOutersInnersOutersInners
P3621712549
PP75561942371413671617

outdro2:CaseDavidsonMethodPowerMethod
AngularNumberNumberNumberNumber
ApproximationOutersInnersOutersInners
P3727912644
P5720514813
P7733415918
Table7.4:Numberofouter(Davidsonorpower)andinner(GMRES)iterationsforthe
Takeda1benchmark.

geometryisaquartercorecontainingfourfuelassembliesandthesurroundingmoderator.
Eachfuelassemblyiscomprisedofa17×17latticeofsquarepincells.Eachofthese
pincellsismadeoutofacylindricalsectionsurroundedbymoderator.Thesecylindrical
sectionseithercontainafuel-cladmixorconstituteavailablespacefortheinsertionofan
absorberrod.ThepincellcompositionsaredetailedinFigure7.4.Thecorresponding
crosssectionvaluescanbefoundin[74].Notethat,unliketheTakeda1benchmark,the
C5G7benchmarkincludesupscattering.
TheParafishmodelusesafiniteelementmeshwithquadrilateralrotatedQ1elements
suchthateachpincellisdiscretizedinto14×14elements.Thisapproximationofthe
cylindricalsectionbyaCartesianmeshwasshownin[106]tobeappropriatebecauseit
preservesnotonlythevolumeratiobetweenthetwopincomponents,butalsothe“density”
ofthemeshinginboththesecomponents.Asfortheangulardiscretization,weconsider
heretheP3approximation.
WeoptimizedtherunningtimesofourbenchmarksbyapplyingtheILUTPprecondi-
tionerineverydomain.Asforthedecomposition,weobtainedgoodresultsusinga20×20
DD.Thegroupingwasdonesuchthatoneorseveralpin-cellsbelongtoonedomain,and,
withintheassemblies,adomainlimitisnecessarilyapin-celllimit(theconverseisnot
true).Theresulting400domainswerespreadamong32to400coressuchthatadecentload
balancingwasobtained.Thisloadbalancingwasperformedsemi-automatically,withsome
fine-tuningdone“byhand”.
ThecalculationsfortheC5G7benchmarkwererunontheJuRoPA(JülichResearch
onPetaflopArchitectures)supercomputerinstalledattheForschungszentrumJülich,Ger-
many.EachofthecoresontheCPUs(Quad-coreIntelXeonprocessors)runsataclock
GHz.2.93ofeedspWehavein[74]akeffectivereferencevalueofkeff=1.186550±0.00008obtainedbythe

7.2.

C5G7NEAThehmarkBenc

yvacuumboundary
26.64

atorderMo84.4217×1717×17
cellspincellspinoundarybreflection17×1717×17
42.21cellspincellspin0021.4242.84
oundarybreflection

oundarybacuumv

x26.64

Figure7.3:Geometryofthe2DC5G7benchmark.

cellPin(e.g.fuel)

Figure7.4:PincellscompositionintheC5G7benchmark,see[74].

77

78

APPLICATIONSINNUCLEARREACTORTHEORY

DavidsonMethodPowerMethod
NumNumbbererInnersOuters3,2871914,65137
NumberCoresTime(s)SpeedupEfficiencyTime(s)SpeedupEfficiency
6432961611.681.0.841.3255281.1.621.0.81
40025620308.005.410.640.68721007.385.280.590.66

Table7.5:Resultsforthe2DC5G7benchmarkwithP3approximation(14,962,584un-
knowns).TheinneriterationsrefertoGMRESiterationsinbothcases.Speed-upsand
efficienciesarecomputedwithrespecttothe32corecase.

Monte-Carlomethod.BothpowerandDavidsonmethodsyieldfortheP3approximationa
satisfyingkeffectiveof1.18319and1.18320respectively.TheDavidsonmethodwasagain
muchmoreefficienttoreachthissolution,asitcanbeseeninTable7.5,cf.[96].
Thecomputingtimesaswellasthenumberofiterationsaredrasticallyreducedbyusing
theDavidsonmethod.Again,applyingaChebyshevaccelerationtothepowermethod
couldnotimproveourresults.

8Chapter

okOutloandSummary

Inthisthesiswehaveconsideredtwoinvolvedclassesofapplicationsintermsofstability.
Wehaveillustratedthatanefficientimplementationexploitingparallelplatformsenables
ustoaddressmoreandmorecomplexproblems.
Withregardstothelinearstability,pseudospectrahavebeenshowntobeanencouraging
tooltofacilitatethestudiesofnonlineareffectswhichmaynotberevealedbyatraditional
splinearectralstpabilitortraitsyofanalysis.ellipticWehadifferenvetialestabliopshederatorsaintheoretitermscalofbacfiniteelekgroundmenttomethoapprods.ximateThe
anbacumericalkwardexamplesfacingstepcovaseredwellinasthisamannaturaluscriptconvincludectioneanprocessinincompressibleanannulus.fluidFfloorwbovother
problems,wehaveillustratedthattheDavidsonmethodiscapabletosolvetheinherent
singularvalueproblemsefficiently.
Aswehaveseen,pseudospectraprotrudingconsiderablyintotherighthalfofthecom-
plexplanemayindicateagrowthofperturbations.However,itstillremainsopentofind
aputationmoreofaccuratepseudospinectraterpretationonofparallelpseudospplatfectraorms.Wcanehabeveacshohievwned.thatNevapertheless,erformanatcom-more
sophisticatedloadbalancingtofurtherexploitthetrivialparallelismcanbedeveloped
inmentorderofrtoegionstuneinthethecoscalabilitmplexy.Fplaneordepinstance,endentoneonthemightwinclorkloadudeofsometheprodynamiccessgroups.assign-
Furthermore,usingparallelpreconditioningtechniques,suchasmultigridorSchurcom-
plementmethods,mayturnouttobevaluabletopreconditiontheinherenteigenvalue
computations.Asforthecriticalityproblem,ournumericalexperimentshaveshownthesuperiorityof
theDacomprehensiblevidsonmethoduedtointhefactcomparisthatonthewithDathesovidsonfarmethowidelydkusedeepspowtracerkmeofthomored.Thisinformationisfully
alongtheiterativeprocedure.Theextrastoragecostsandadditionalcomputationaleffort
nperumberiterationofofiterations.theDaFvidsonurthermore,methothdeareDaoutvidsonweighmethotedbdyappthesearstoignificanbetmorereductionrobustofthanthe
thepowermethod:First,theDavidsonmethodcontrolsthe(mathematically)exactresid-
ualsofthematrixpencil.Second,itneedsalessaccuratesolvingineachpreconditioning
stepthanthepowermethod.Thepowermethodappliedtoamatrixpencilmaynotneed
anmayexactnotconsolvingvergeofattheallifarithissinglinearsolvingissystemsdoneintopracticeroughly.asitisrequiredintheory,butit
mightTherebeisusefulstillpinotenthetialfieldtooffurtherneutrontunetrtheansportcomputatiotheoryninofcritigeneral.calIneigenvaluesparticular,whicforh

79

OUTLOOK

AND

YSUMMAR

80

for

eb

the

energy

endencedep

urthermore,Faluable.v

insteadellevenergyof

the

hemesc

toout

onsolving

dcouple

a

upscattering,

large

gincludin

problems

can

eb

ducedtroin

to

improev

ablock-wiseappliedGauss-Seidelmethodmayturn

especiallyforlarge-scaledproblems,aparallelization

.yscalabilit

AendixApp

PEstimatesoincaréandConstantinCalculationsanAnnforulusthe

nInthisumericalcevhapteraluatweionsderivwhiceahbgiveoundusforantheinsighPtaboincaréouttheconstanqualittinyofanthisannbulusound.andWpeerformapply
mothesedelresultsdescribingtomakenaturalquanconvtitativeection.statTheseementsofresultstheareappropublishedximationin[qualit54].yofansimplified

A.1BoundsforthePoincaréConstant

AsalreadyseeninSection2.2,thePoincaréconstantisaveryimportanttooltoinvestigate
fluidflows.Forinstance,itisneededtoproveexistenceofweaksolutionsandregularity
assumptionsinthecontextofellipticdifferentialequations,seee.g.[19].Butasseenin
Section2.2,itplaysamajorroleinthestabilitytheoryaswell.
LetΩ⊂Rnbeboundedinatleastonedirection.ThePoincaréconstantisdefinedas
thesmallestconstantkpsuchthat

u0≤kpu0(A.1)
holdsforallu∈H01(Ω).IfΩisbounded,onecanshowthatkp=λ1−1/2,whereλ1isthe
smallesteigenvalueoftheLaplaceproblem
−Δu=λuinΩ,u=0on∂Ω,(A.2)
seecane.beg.fo[95und].inThe[44fol].lowingtwoboundsweestablishholdforquitegeneraldomainsΩand

TheoremA.1IfΩ⊂{x∈Rn:−d/2<xn<d/2},thenwehave
du0≤πu0

forallu∈H01(Ω).
Proof.SinceC0∞(Ω)isdenseinH01(Ω),itissufficienttoshowtheassertionholdsforall

81

82ESTIMATESANDCALCULATIONSFORTHEPOINCARÉCONSTANT

u∈C0∞(Ω).Wedefinethefunction
)x(uU(x)=sin[π(xn+d/2)/d].
Clearly,Uisboundedandvanishesat−d/2andd/2.Wehave
d/2∂uππd2
0≤−d/2∂xn−ducotd(xn+2)dxn
d/2∂u2πd/2∂uπd
=−d/2∂xndxn−d−d/22∂xnucotd(xn+2)dxn
π2d/2πd2
+d2ucotd(xn+2)dxn.
2d/−Integratingthesecondtermbypartsyields
πd/2∂uπd
−d2∂xnucotd(xn+2)dxn
2d/−d/2d/2
=−πu2cotπ(xn+d)+πu2−πsin−2π(xn+d)dxn
dd2−d/2d−d/2dd2
d/22d/2
dd2−d/2d−d/2d2
=−πucosπ(xn+d)U−π2u2sin−2π(xn+d)dxn
=−π2u2sin−2π(xn+d)dxn.
2d/2
d−d/2d2
Therefore,wehave
d/2∂u2π2d/22−2πdπd
0≤−d/2∂xndxn−d2−d/2usind(xn+2)−cotd(xn+2)dxn.
Sincesin−2x−cot2x=1+cos2x=1+1−sin2x=2−1≥1,
sin2xsin2xsin2x
obtainew2d/2d2d/2∂u2
−d/2udxn≤π2−d/2∂xndxn.
LetHn−1denotethe(n−1)-dimensionalhyperplaneofthefirstn−1dimensions,i.e.
Hn−1=span{x1,x2,...,xn−1}.Then,wehave
2d/22d/2d2∂u2
u0=Hn−1−d/2udxnd(x1,...,xn−1)≤Hn−1−d/2π2∂xndxnd(x1,...,xn−1)
d2∂u2d2n∂u2d22
=π2Ω∂xndx≤π2Ω∂xidx=π2u0,
=1iwhichcompletestheproof.✷

A.1.BoundsforthePoincaréConstant

83

(A.3)

TheoremA.2LetΩ⊂{x∈Rn:−d/2<xi<d/2,i=1,...,n},then
du0≤π√nu0(A.3)
foranyu∈H01(Ω).
Proof.Again,itisenoughtoshowthetheoremforanyu∈C0∞(Ω).Fromtheproofofthe
lasttheoremweknowthat
dxdxud/22d2d/2∂u2
≤−d/2nπ2−d/2∂xii
holdsforanyi=1,...,n.Denotingthe(n−1)-dimensionalhyperplaneofthedimensions
(1,...,i−1,i+1,...,n)byHin−1,wededuce
nnu02=u2dx
Ω=1i2nd/2
=Hn−1−d/2udxid(x1,...,xi−1,xi+1,...,xn)
i=1ind2d/2∂u2
≤Hn−1π2−d/2∂xidxid(x1,...,xi−1,xi+1,...,xn)
i=1id2n∂u2d22
=π2Ω∂xidx=π2u0,
=1iwhichcompletestheproof.✷
Fromthelasttheorem,weobtainintermsoftheannulusΩA={(r,ϕ):A/2<r<
1+A/2,0≤ϕ<2π}theestimate
√2Akp≤π1+2.(A.4)
ThisboundisonlyusefulforsmallAbecausethegapwidthisnottakenintoaccount.If
wecombine(A.4)withtheboundderivedin[42],weobtain
12√2A
kp≤min21+A,π1+2.(A.5)

TheoremA.3IfΩ={(r,ϕ):Ri<r<Ro,0≤ϕ<2π},thenforthePoincaréconstant
kpwehavethebound
kp≤RoRo−Ri.
πRi

84ESTIMATESANDCALCULATIONSFORTHEPOINCARÉCONSTANT

Proof.Introducingpolarcoordinatesweobtain
2Ω|w|2d(x,y)Ωr|w˜|2d(r,ϕ)
Ωkp=w∈Hmax01(Ω)Ωw∙wd(x,y)=w˜∈Hmax10(Ω)r(∂rw˜)2+r−2(∂ϕw˜)2d(r,ϕ),(A.6)
wherew˜(r,ϕ)=w(x,y).From[71,Proposition1.1]weknowthattheargumentforwhich
(A.6)attainsitsmaximumvalueisaradialfunctionu˜1=u˜1(r).Hence,thereholds
2RRi0r|u˜1(r)|2dr
kp=R0r|u˜1(r)|2dr.
RiNotethatu˜1istheeigenfunctionassociatedtothesmallesteigenvalueoftheLaplace
problem(A.2).Further,weknowfromtheone-dimensionalLaplaceeigenvalueproblem
thatRo−Ri2RRi0|w˜(r)|2dr
π=w˜∈H01(maxRi,Ro)RRi0|w˜(r)|2dr,
whichincombinationwith(A.1)implies
k2≤RoRRi0|˜u1(r)|2dr≤RoRo−Ri2.
pRiRRi0|u˜1(r)|2drRiπ
✷ByvirtueofTheoremA.3wealreadyobtainanewboundforthePoincaréconstantwith
21respecttoΩA,namely
kp≤π1+A.
Togetherwith(A.4)wehave
√
kp≤min11+2,21+A.(A.7)
2ππAA.2EvaluationofthePoincaréConstant
InordertoperformanefficientevaluationofthePoincaréconstant,wewanttoexploitits
one-dimensionalcharacter.Therefore,wetransformtheeigenvalueproblem(A.2)topolar
coordinatesandusetheresultof[71]thatalleigenfunctionsassociatedtothesmallest
eigenvalueareradial.Hence,wehavethatkp=λ1−1/2,whereλ1isthesmallesteigenvalue
oftheone-dimensionalproblem
−u˜−1u˜=λu˜in(Ri,Ro),u˜(Ri)=u˜(Ro)=0.
rThesubsequentcomputationswereperformedintheframeworkofthefiniteelementsoft-
wareHiFlow,see[1].First,theproblemwasdiscretizedbymeansofafinitedifference
schemeofsecondorder.Then,thealgebraiceigenvalueproblemwassolvedbymeansof
theDavidsonmethod,seeSection5.1.Thenumericalresultsaswellthebounds(A.5)and

(A.7)

A.3.ApplicationinNaturalConvection

FigureA.1:aluationEvofoincaréPtsconstanwithesprectto.A

85

(A.7)areplottedinFigureA.1.Theperformedcomputationsshowthattheestablished
boundinTheoremA.3isalmostsharpforlargeA.

(A.8)

A.3ApplicationinNaturalConvection
Weconsiderthesteady-stateformulationoftheOberbeck-Boussinesqequations(2.8)asin
thesetupofSection6.2.Byusingpolarcoordinates(r,ϕ)andthedefineddimensionless
quantitiesA,Ra,Pr,weobtainthedimensionlessformulation
aR1Pr(v∙)v−Δv+p=Bsinϕer+Raτe3,
∙v=0,(A.8)
vrv∙τ−Δτ=rB
inΩA,withboundaryconditionsv=0,τ=0on∂ΩA,see[77].Here,τ=τ(r,ϕ)isthe
excesstemperatureandwehavesetB=ln(1+2/A)=ln(Ro/Ri).Furthermore,e3isthe
unitvectorintodirectionz,i.e.e3=sinϕer+cosϕeϕ.
Becausenoexactsolutionsareknown,simplifiedmodelsmaybemoreconvenient.For
instanceonemaydecouplethesystembyreplacingthefirstequationof(A.8)byalinear
yieldingequation−Δv+p=RBasin(ϕ)er,
∙v=0,(A.9)
v∙τ−Δτ=vr
rB

(A.9)

86ESTIMATESANDCALCULATIONSFORTHEPOINCARÉCONSTANT
FigureA.2:ContourplotsofE=0.01,0.05,0.1.
inΩAwithzeroboundaryconditionsforvandτon∂Ω.Notethatthethirdequation
isknostwn.illThenonlinear.modelF(orA.9)thisissystempreferredanforanalyticaldescribingsolutionnaturalincontermvesctionofinBesselcaseoffunctionsasmallis
A(usuallyA<2.8).
systemAn(approA.8)xiandmationthescsolutionhemetoofestiitsmatesimplitheficationrelativ(eA.9er)rorinbetwtermseenoftheasoluquantiontityofXtheunderfull
considerationisdescribedin[62].First,anupperboundf1onthenormoftheabsolute
errorδXisstated:
δX≤f1(kp,Ra,Pr,A).(A.10)
Second,alowerboundf2onthenormofthesolutionX0ofthesimplifiedproblemis
established:f2(Ra,Pr,A)≤X0.(A.11)
Finally,bycombining(A.10)and(A.11),weobtain
δXδXf1
E=X≤X0−δX≤f2−f1(A.12)
asanupperboundoftherelativeerrorE,providedthatf2>f1.Intermsofvandτ
thefunctionsf1andf2canbefoundin[62].Sincethesefunctionsarerathertechnical,
wedonotstatethemhere.Usingthederivedestimate(A.7)forthePoincaréconstant,we
havecomputedupperboundsfortherelativeerrorEbyevaluatingf1andf2.InFigureA.2
countourplotsoftheupperbound(A.12)asafunctionofRaandA,withfixedPrandtl
numberPr=0.7,intermsofvareshown.

yBibliograph

[1]HiFlow-ageneralfiniteelementtoolbox.http://www.hiflow.de,2011.
[2]Alonso,A.,andDelloRusso,A.Spectralapproximationofvariationally-posedeigen-
valueproblemsbynonconformingmethods.JournalofComputationalandAppliedMathe-
matics223,1(2009),177–197.
[3]Alt,H.W.LineareFunktionalanalysis,5thed.Springer,2006.
[4]Amann,H.GewöhnlicheDifferentialgleichungen,2nded.deGruyter,1995.
[5]Amestoy,P.R.,Duff,I.S.,Koster,J.,andL’Excellent,J.-Y.Afullyasynchronous
multifrontalsolverusingdistributeddynamicscheduling.SIAMJournalonMatrixAnalysis
andApplications23,1(2001),15–41.
[6]Anzt,H.,Augustin,W.,Baumann,M.,Bockelmann,H.,Gengenbach,T.,Hahn,
T.,Heuveline,V.,Ketelaer,E.,Lukarski,D.,Otzen,A.,Ritterbusch,S.,
Rocker,B.,Ronnås,S.,Schick,M.,Subramanian,C.,Weiß,J.-P.,andWil-
helm,F.HiFlow3-aflexibleandhardware-awareparallelfiniteelementpackage.EMCL
2010.eries,StPreprin[7]Aris,R.Vectors,tensors,andthebasicequationsoffluidmechanics.DoverPublications,
1989.[8]Aulbach,B.GewöhnlicheDifferentialgleichungen.SpektrumAkademischerVerlag,1997.
[9]Babuska,I.,andOsborn,J.Eigenvalueproblems.InFiniteElementMethods(Part1),
vol.2ofHandbookofNumericalAnalysis.Elsevier,1991,pp.641–787.
[10]Balay,S.,Buschelman,K.,Eijkhout,V.,Gropp,W.D.,Kaushik,D.,Knepley,
M.G.,McInnes,L.C.,Smith,B.F.,andZhang,H.PETScusersmanual.Tech.Rep.
ANL-95/11-Revision3.0.0,ArgonneNationalLaboratory,2008.
[11]Balay,S.,Buschelman,K.,Gropp,W.D.,Kaushik,D.,Knepley,M.G.,McInnes,
L.C.,Smith,B.F.,andZhang,H.PETScwebpage.http://www.mcs.anl.gov/petsc,
2011.[12]Balay,S.,Gropp,W.D.,McInnes,L.C.,andSmith,B.F.Efficientmanagement
ofparallelisminobjectorientednumericalsoftwarelibraries.InModernSoftwareToolsin
ScientificComputing(1997),E.Arge,A.M.Bruaset,andH.P.Langtangen,Eds.,Birkhäuser
163–202.pp.Press,[13]Bell,G.I.,andGlasstone,S.Nuclearreactortheory.VanNostrandReinhold,1970.
[14]Bell,G.I.,Hansen,G.E.,andSandmeier,H.A.Multitabletreatmentsofanisotropic
scatteringinSNmultigrouptransportcalculations.NuclearScienceandEngineering,28
(1967).[15]Bénard,H.Lestourbillonscellulairesdansunenappeliquide.RevuegénéraledesSciences
puresetappliquées11(1900).

87

88

BIBLIOGRAPHY

[16]Bertsch,C.Finite-Elemente-DiskretisierungundeinMehrgitterverfahrenfürunsymme-
trischeEigenwertprobleme.Diplomarbeit,UniversitätHeidelberg,1999.
[17]Boltzmann,L.WeitereStudienüberdasWärmegleichgewichtunterGasmolekülen.Wiener
275–370.(1872),66Berichte[18]Boussinesq,J.Théorieanalytiquedelachaleur.Gauthier-Villars,1903.
[19]Braess,D.FiniteElemente,4thed.Springer,2007.
[20]Bramble,J.H.,andOsborn,J.E.Rateofconvergenceestimatesfornonselfadjoint
eigenvalueapproximations.MathematicsofComputation27,123(1973),525–549.
[21]Brenner,S.C.,andScott,L.R.Themathematicaltheoryoffiniteelementmethods,
2008.pringer,Sed.3rd[22]CaPhysicso,Y.,ofRandeactorsLee,J.“NucleC.arAPowerKrylo:vsAubspaceSustainablemethoRdesourtoce”calculate(PHk-YSORandα-2008)eigen(Invaluesterlak.enIn,
Switzerland,September14–19,2008).
[23]Carpraux,J.F.,Erhel,J.,andSadkane,M.Spectralportraitfornonhermitianlarge
sparsematrices.Computing53(1994),301–310.
[24]Cercignani,C.TheBoltzmannequationanditsapplications.Appliedmathematicalsci-
1988.Springer,67.;ences[25]Chaitin-chatelin,F.,andHarrabi,A.Aboutdefinitionsofpseudospectraofclosed
operatorsinBanachspaces.Tech.Rep.TR/PA/98/08,CERFACS,1998.
[26]Chandrasekhar,S.Hydrodynamicandhydromagneticstability.DoverPublications,1981.
[27]Choi,H.,Hinze,M.,andKunisch,K.Instantaneouscontrolofbackward-facingstep
flows.AppliedNumericalMathematics31,2(1999),133–158.
[28]Chorin,A.J.,andMarsden,J.E.Amathematicalintroductiontofluidmechanics,
2000.pringer,Sed.3rd[29]Ciarlet,P.G.Thefiniteelementmethodforellipticproblems.Classicsinappliedmath-
ematics;40.SocietyforIndustrialandAppliedMathematics,2002.
[30]Crouzeix,M.,Philippe,B.,andSadkane,M.TheDavidsonmethod.SIAMJournal
onScientificComputing15(1994),62–76.
[31]Dalecki˘ı,J.L.,andKre˘ın,M.G.StabilityofsolutionsofdifferentialequationsinBanach
space.Translationsofmathematicalmonographs;43.AmericanMathematicalSociety,1974.
[32]Davidson,E.R.Theiterativecalculationofafewofthelowesteigenvaluesandcorre-
spondingeigenvectorsoflargereal-symmetricmatrices.JournalofComputationalPhysics
87–94.(1975),17[33]Demmel,J.W.,Eisenstat,S.C.,Gilbert,J.R.,Li,X.S.,andLiu,J.W.H.
Asupernodalapproachtosparsepartialpivoting.SIAMJournalonMatrixAnalysisand
Applications20,3(1999),720–755.
[34]Descloux,J.,Nassif,N.,andRappaz,J.Onspectralapproximation.Part1.The
problemofconvergence.RAIROAnalyseNumérique12,2(1978),97–112.
[35]Descloux,J.,Nassif,N.,andRappaz,J.Onspectralapproximation.Part2.Error
estimatesfortheGalerkinmethod.RAIROAnalyseNumérique12,2(1978),113–119.
[36]Dongarra,J.,Lumsdaine,A.,Niu,X.,Pozo,R.,andRemington,K.LAPACK
workingnote74:AsparsematrixlibraryinC++forhighperformancearchitectures.Tech.
1994.rep.,

BIBLIOGRAPHY

89

[37]Drazin,P.G.Introductiontohydrodynamicstability.Cambridgetextsinappliedmathe-
matics;32.CambridgeUniversityPress,2002.
[38]Drazin,P.G.,andReid,W.H.HydrodynamicStability,2nded.CambridgeUniversity
2004.Press,[39]Dyko,M.P.,Vafai,K.,andMojtabi,A.K.Anumericalandexperimentalinvesti-
gationofstabilityofnaturalconvectiveflowswithinahorizontalannulus.JournalofFluid
27–61.(1999),381chanicsMe[40]Embree,M.,andTrefethen,L.N.Generalizingeigenvaluetheoremstopseudospectra
theorems.SIAMJournalonScientificComputing23,2(2001),583–590.
[41]Ern,A.,andGuermond,J.-L.Theoryandpracticeoffiniteelements.Appliedmathe-
maticalsciences;159.Springer,2004.
[42]Ferrario,C.,Passerini,A.,andPiva,S.AStokes-likesystemfornaturalconvection
inahorizontalannulus.NonlinearAnalysis:RealWorldApplications9,2(2008),403–411.
[43]Fraysse,V.,Gueury,M.,Nicoud,F.,andToumazou,V.Spectralportraitsformatrix
1996.encils,p[44]vol.Galdi,1:G.LinearisedP.AnintrsteadyoductionproblemstotheofSpringermathematictralactstheinorynaturofalthephilosophyNavier-Stokes;38.equationsSpringer,,
1998.[45]Gill,D.F.,andAzmy,Y.Y.Jacobian-freeNewton-Krylovasanalternativetopower
iterationsforthek-eigenvaluetransportproblem.ANSTransactions92(2005),728–730.
[46]Gill,D.F.,andAzmy,Y.Y.AJacobian-freeNewton-Kryloviterativeschemeforcriti-
calitycalculationsbasedontheneutrondiffusionequation.InMathematics,Computational
Methods&ReactorPhysics(M&C2009)(SaratogaSprings,NewYork,USA,May3–7,
2009).[47]Girault,V.,andRaviart,P.-A.FiniteelementmethodsforNavier-Stokesequations.
Springerseriesincomputationalmathematics;5.Springer,1986.
[48]LineGodunoarAv,lgebrS.aK.,andanditsApplicSadkane,ationsM.279,1–3Computation(1998),of163–1pseudosp75.ectraviaspectralprojectors.
[49]Golub,G.H.,andVanLoan,C.F.Matrixcomputations,3rded.JohnsHopkins
1997.Press,yersitUniv[50]Hackbusch,W.Ellipticdifferentialequations.Springerseriesincomputationalmathemat-
2003.Springer,18.;ics[51]Hahn,W.Stabilityofmotion.DieGrundlehrendermathematischenWissenschaftenin
1967.Springer,138.;Einzeldarstellungen[52]Havet,M.SolutionofAlgebraicProblemsarisinginNuclearReactorSimulationsusing
Jacobi-DavidsonandMultigridMethods.PhDthesis,UniversitéLibredeBruxelles,Belgium,
2008.[53]Heuveline,V.Finiteelementapproximationsofeigenvalueproblemsforellipticpartial
differentialoperators.Habilitationsschrift,UniversitätHeidelberg,2002.
[54]Heuveline,V.,Passerini,A.,Subramanian,C.,andThäter,G.Estimatesand
calculationsforthePoincaréconstantinanannulus.EMCLPreprintSeries,toappear.
[55]Heuveline,V.,Philippe,B.,andSadkane,M.Parallelcomputationofspectralportrait
oflargematricesbyDavidsontypemethods.NumericalAlgorithms16,1(1997),55–75.

90

BIBLIOGRAPHY

[56]Heuveline,V.,andRannacher,R.Aposteriorierrorcontrolforfiniteelementap-
proximationsofellipticeigenvalueproblems.AdvancesinComputationalMathematics15,1
107–138.(2001),[57]Heuveline,V.,andRannacher,R.AdaptiveFEMforeigenvalueproblemswithappli-
cationinhydrodynamicstabilityanalysis.JournalofNumericalMathematics(2006),1–32.
[58]Heuveline,V.,Subramanian,C.,Lukarski,D.,andWeiß,J.-P.Amulti-platform
linearalgebratoolboxforfiniteelementsolversonheterogeneousclusters.InIEEECluster
2010,WorkshoponParallelProgrammingandApplicationsonAcceleratorClusters(PPAAC
2010)(Heraklion,Greece,September20–24,2010).
[59]Kato,T.Perturbationtheoryforlinearoperators.DieGrundlehrendermathematischen
WissenschafteninEinzeldarstellungen;132.Springer,1966.
[60]Knoll,D.A.,Park,H.,andSmith,K.ApplicationoftheJacobian-freeNewton-
Krylovmethodincomputationalreactorphysics.InMathematics,ComputationalMethods
&ReactorPhysics(M&C2009)(SaratogaSprings,NewYork,USA,May3–7,2009).
[61]Kolata,W.G.Approximationinvariationallyposedeigenvalueproblems.Numerische
159–171.(1978),29Mathematik[62]Lamacz,A.,Passerini,A.,andThäter,G.Naturalconvectioninhorizontalannuli:
evaluationoftheerrorfortwoapproximations.Submittedforpublication.
[63]Lathouwers,D.Computingtime-eigenvaluesusingtheeven-paritytransportform.Annals
ofNuclearEnergy33,10(2006),941–943.
[64]Lavallée,P.-F.,andSadkane,M.Computationofpseudospectrabyspectraldichotomy
methodsinaparallelenvironment.NumericalAlgorithms33,1(2003),343–355.
[65]Lewis,E.E.Second-orderneutrontransportmethods.InNuclearComputationalScience,
Y.AzmyandE.Sartori,Eds.Springer,2010,pp.85–115.
[66]Lewis,E.E.,andMiller,W.F.Computationalmethodsofneutrontransport.Wiley,
1984.[67]Li,partialX.pivS.,oting.andAShaCMo,TrM.ansactionsAsuponernoMdalathematicapproacalhtoSoftwarincome37plete,4(20LU10).factorizationwith
[68]byMermixedcier,andB.,hybridOsborn,methoJ.,ds.Rappaz,MathematicsJ.,andofRaviarComputationt,P.36A.,154,Eigenvalue427–453.approximation
[69]Morgan,R.B.GeneralizationsofDavidson’smethodforcomputingeigenvaluesoflarge
nonsymmetricmatrices.JournalofComputationalPhysics101,2(1992),287–291.
[70]ingMoreigengan,valuesR.B.,ofandsparseScott,symmetricD.S.matrices.GeneralizationsSIAMofJournaDalonvidson’sScientificmethodandforStatisticcomput-al
Computing7,3(1986),817–825.
[71]Nazarinequalitoyv,inA.I.sphericalTheandplane.one-dimensionalJournalcofharacterofMathematicanalextremSciencumespoin102t,of5the(2000),Friedric4473–hs
4486.[72]Notay,Y.IsJacobi-DavidsonfasterthanDavidson?SIAMJournalonMatrixAnalysis
andApplications26,2(2005),522–543.
[73]Oberbeck,A.ÜberdieWärmeleitungderFlüssigkeitenbeiBerücksichtigungderStrö-
mungeninfolgevonTemperaturdifferenzen.AnnalenderPhysik243(1879),271–292.
[74]OECD/NEA.Benchmarkondeterministictransportcalculationswithoutspatialhomogeni-
sation.A2-D/3-DMOXfuelassemblybenchmark.InNEA/NSC/DOC(2003)16(2005),
NuclearScience,OECD,NuclearEnergyAgency.

BIBLIOGRAPHY

91

[75]Osborn,J.E.Spectralapproximationforcompactoperators.MathematicsofComputation
712–725.(1975),131,29[76]Osborn,J.E.Approximationoftheeigenvaluesofanonselfadjointoperatorarisinginthe
studyofthestabilityofstationarysolutionsoftheNavier-Stokesequations.SIAMJournal
onNumericalAnalysis13,2(1976),185–197.
[77]Passerini,A.,Ferrario,C.,Růžička,M.,andThäter,G.Theoreticalresultson
steadyconvectiveflowsbetweenhorizontalcoaxialcylinders.SIAMJournalonApplied
(submitted).Mathematics[78]Passerini,A.,Růžička,M.,andThäter,G.Naturalconvectionbetweentwohorizontal
coaxialcylinders.ZeitschriftfürAngewandteMathematikundMechanik89,5(2009),399–
413.[79]Pazy,A.Semigroupsoflinearoperatorsandapplicationstopartialdifferentialequations,
2nded.Appliedmathematicalsciences;44.Springer,1992.
[80]largePhilippe,matrix.B.,LineandarAlgebrSadkane,aandM.itsApplicComputationationsof257the(1997),fundamen77–104.talsingularsubspaceofa
[81]Quarteroni,A.,andValli,A.Numericalapproximationofpartialdifferentialequations.
Springerseriesincomputationalmathematics;23.Springer,1994.
[82]Rannacher,R.NumerischeMathematik3-NumerikvonProblemenderKontinu-
umsmechanik.Vorlesungsskriptum,2006.
[83]RannaNumericalcher,MethoR.,dsandforPartialTurek,DifferS.AentialsimpleEquationsnonconforming8(1992),97–111.quadrilateralStokeselement.
[84]Riedel,K.S.Generalizedepsilon-pseudospectra.SIAMJournalonNumericalAnalysis
31,4(1994),1219–1225.
[85]Ruhe,A.TherationalKrylovalgorithmforlargenonsymmetriceigenvalues–mappingthe
resolventnorms(pseudospectrum),1994.
[86]Saad,Y.Numericalmethodsforlargeeigenvalueproblems.ManchesterUniversityPress,
1992.[87]Saad,Y.Iterativemethodsforsparselinearsystems,2nded.SocietyforIndustrialand
2003.,MathematicsApplied[88]Sadkane,M.Block-ArnoldiandDavidsonmethodsforunsymmetriclargeeigenvalueprob-
lems.JournalofNumericalMathematics64(1993),195–211.
[89]Sattinger,D.Themathematicalproblemofhydrodynamicstability.JournalofMathe-
maticsandMechanics19,9(1970),797–817.
[90]Schmid,P.J.,andHenningson,D.S.Stabilityandtransitioninshearflows.Applied
mathematicalsciences;142.Springer,2001.
[91]Scott,R.Interpolatedboundaryconditionsinthefiniteelementmethod.SIAMJournal
onNumericalAnalysis12,3(1975),404–427.
[92]Sleijpen,G.L.G.,andVanderVorst,H.A.AJacobi-Davidsoniterationmethod
forlineareigenvalueproblems.SIAMReview42,2(2000),267–293.
[93]Sobolev,S.L.Applicationsoffunctionalanalysisinmathematicalphysics.Translations
ofmathematicalmonographs;7.AmericanMathematicalSociety,1963.
[94]Straughan,B.Theenergymethod,stability,andnonlinearconvection.Appliedmathe-
maticalsciences;91.Springer,1992.
[95]Strauss,W.A.Partialdifferentialequations.Wiley,1992.

92

BIBLIOGRAPHY

[96]Subramanian,C.,VanCriekingen,S.,Heuveline,V.,Nataf,F.,andHavé,P.The
Davidsonmethodasanalternativetopoweriterationsforcriticalitycalculations.Submitted
publication.for[97]Subramanian,C.,VanCriekingen,S.,Heuveline,V.,Nataf,F.,andHavé,P.
TheDavidsonmethodasanalternativetopoweriterationsforcriticalitycalculations.In
MathematicsandComputationalMethodsappliedtoNuclearScienceandEngineering(MC
2011)(RiodeJaneiro,Brazil,May8–12,2011,submitted).
[98]Süli,E.,andMayers,D.F.Anintroductiontonumericalanalysis.CambridgeUniversity
2003.Press,[99]Takeda,T.,andIkeda,H.Finalreporton3Dneutrontransportbenchmark.Department
ofNuclearEngineering,OsakaUniversity,Japan,1991.
[100]Takeda,T.,Tamitani,M.,andUnesaki,H.Proposalof3Dneutrontransportbench-
mark.NEACRPA-953REV.1(1989).
[101]Toh,SIAMJouK.-C.,rnalandonScientificTrefethen,L.ComputingN.17,Calculation1(1996),of1–15.pseudospectrabytheArnoldiiteration.
[102]Trefethen,L.N.Computationofpseudospectra.ActaNumerica8(1999),247–295.
[103]Trefethen,L.N.,andEmbree,M.Spectraandpseudospectra.PrincetonUniversity
2005.Press,[104]Trefethen,L.N.,Trefethen,A.E.,Reddy,S.C.,andDriscoll,T.A.Hydrody-
namicstabilitywithouteigenvalues.Science261,5121(1993),578–584.
[105]Tritton,D.J.Physicalfluiddynamics.Themodernuniversityphysicsseries.VanNostrand
1979.Reinhold,[106]VanCriekingen,S.A2D/3Dcartesiangeometrynon-conformingsphericalharmonic
neutrontransportsolver.AnnalsofNuclearEnergy34,3(2007),177–187.
[107]VanCriekingen,S.,Nataf,F.,andHavé,P.Parafish:aparallelFE−PNneutron
transportsolverbasedondomaindecomposition.AnnalsofNuclearEnergy38,1(2011),
145–150.[108]vanDorsselaer,J.L.M.Pseudospectraformatrixpencilsandstabilityofequilibria.
BITNumericalMathematics37(1997),833–845.
[109]vanDorsselaer,J.L.M.Severalconceptstoinvestigatestronglynonnormaleigenvalue
problems.SIAMJournalonScientificComputing24,3(2002),1031–1053.
[110]Verdú,G.,Ginestar,D.,Miró,R.,andVidal,V.UsingtheJacobi-Davidsonmethod
toobtainthedominantlambdamodesofanuclearpowerreactor.AnnalsofNuclearEnergy
1274–1296.(2005),32[111]methoVerdú,d,anG.,efficienMiró,tR.,alternativGinestetoar,solvDe.,theandneutronVidal,diffusionV.Theequation.implicitAnnalsrestartedofNucleArnoldiar
Energy26(1999),579–593.
[112]WKryloarsa,vsuJ.Sbspace.,Witerationsareing,Tfor.A.,deterministicMorel,J.kE.,-eigenvalueMcGhee,J.calculations.M.,andNuclearLehoucq,ScienceR.andB.
Engineering147,1(2004),26–42.
[113]Werner,D.Funktionalanalysis,6thed.Springer,2007.
[114]Wesseling,P.Principlesofcomputationalfluiddynamics.Springerseriesincomputational
2009.Springer,29.;mathematics[115]Williams,W.S.C.Nuclearandparticlephysics.Oxfordsciencepublications.Clarendon
1992.Press,

BIBLIOGRAPHY

[116]

[117]

[118]

[119]

[120]

Wloka,J.Partialdifferentialequations.CambridgeUniversityPress,1992.

93

ders.Yoo,J.-S.InternationalDualsteadyJournalsolutionsofHeatinandnaturalFluidconvFlowection17,6bet(19ween96),horizon587–593.talconcentriccylin-

Yoo,nationalJ.-S.JournalNaturalofHeconatvandectionMassinaTrnarroansferw41,horizon20tal(1998),cylindrical3055–3073ann.ulus:Pr≤0.3.Inter-

Yoo,J.-S.Transitionandmultiplicityofflowsinnaturalconvectioninanarrowhorizontal
cylindricalannulus:Pr=0.4.InternationalJournalofHeatandMassTransfer42,4(1999),
709–722.

Yosida,K.Functionalanalysis,4thed.DieGrundlehrendermathematischenWissenschaf-
teninEinzeldarstellungen;123.Springer,1974.