Two-level preconditioned conjugate gradient methods with applications to bubbly flow problems

Gradient Methods with Appli ations to

Bubbly Flow Problems

PROEFSCHRIFT

ter verkrijging vande graadvando tor

aande Te hnis heUniversiteit Delft,

opgezag vande Re torMagni us prof.dr.ir. J.T.Fokkema,

voorzittervan hetCollege voorPromoties,

in hetopenbaarteverdedigen op

maandag8september 2008om15:00 uur

door

Jok Man TANG

wiskundig ingenieur

Samenstellingpromotie ommissie:

Re tor Magni us, voorzitter

Prof.dr.ir. C.Vuik, Te hnis he Universiteit Delft,promotor

Prof.dr.ir. H. Bijl, Te hnis he Universiteit Delft

Prof.dr.ir. B.J.Boersma, Te hnis he Universiteit Delft

Prof.dr. R.Nabben, Te hnis he UniversitätBerlin, Duitsland

Prof.dr.ir. C.W.Oosterlee, Te hnis he Universiteit Delft

Prof.dr.ir. S. Vandewalle, KatholiekeUniversiteit Leuven,België

Prof.dr.ir. P. Wesseling, Te hnis he Universiteit Delft

Keywords: onjugate gradient method, two-level pre onditioners, deation, domain

de omposition,multigrid, bubbly ows, Poisson equationwitha dis ontinuous

oe- ient, singularsymmetri positivesemi-denite matri es.

The work des ribed in this dissertation was arried out in the se tion of Numeri al

Analysisatthe Department ofApplied Mathemati s,DelftInstituteofApplied

Math-emati s,Fa ultyof Ele tri alEngineering,Mathemati sandComputer S ien e,Delft

UniversityofTe hnology, TheNetherlands.

Part of the resear h des ribed in this dissertation has been funded by the Dut h

BSIK/BRICKSproje t.

Two-Level Pre onditioned Conjugate Gradient Methods with Appli ations to Bubbly

Flow Problems.

DissertationatDelft UniversityofTe hnology.

Copyright

Tomyparents,

mymother LisaTang-Lam,

Two-Level Pre onditioned Conjugate Gradient Methods

with Appli ations to Bubbly Flow Problems

Jok M.Tang

The Pre onditioned Conjugate Gradient (PCG) method is one of the most popular

iterative methodsforsolvinglargelinearsystems witha symmetri andpositive

semi-denite oe ient matrix. However, if the pre onditioned oe ient matrix is

ill- onditioned, the onvergen e of the PCG method typi ally deteriorates. Instead, a

two-levelPCGmethod anbeused. The orrespondingtwo-levelpre onditionerusually

treatsunfavorableeigenvaluesofthe oe ientmatrixee tively,sothatthetwo-level

PCG method is expe ted to onverge faster than the original PCG method. Many

two-level pre onditioners are known in the elds of deation, multigrid and domain

de ompositionmethods. Severalof them are dis ussedin this thesis,where the main

fo usison the deationmethod.

Weshow some theoreti al properties of the deation method, whi hgive insights

into the ee tiveness of this method. A ru ial omponent of the deation

pre on-ditioner is the hoi e of proje tion ve tors. Several hoi es are dis ussed and

exam-ined. Weadvo atethatsubdomainproje tionve tors, whi hare basedondisjointand

pie ewise- onstantve tors, areamongthe best hoi esfor a lassofproblems.

Subsequently,weexaminetheappli ationofthedeationmethodtolinearsystems

with singular oe ient matri es. Several mathemati ally equivalent variants of the

original deation method are proposed to deal with the possible singularity of this

oe ientmatrix. Inaddition,twoapproa hesaredis ussed inorder tohandle oarse

linear systems with a Galerkin matrix, whi h are involved in ea h iteration of the

deation method. After the dis ussionof the implementation ande ien y issuesof

the deation method, it is demonstrated that this method is usually faster than the

originalPCG method.

Moreover, we presenta omparison between thedeation method and other

well-knowntwo-levelPCG methods,amongthemthe balan ing-Neumann-Neumann,

nonsym-metri V- y les. Asthe parameters ofthe orrespondingtwo-levelpre onditionersare

abstra t, we show that these methods are strongly onne ted to ea h other. The

omparisonisalsodonewherethedierenttwo-levelPCGmethodsadopttheirtypi al

and optimized set of parameters. Numeri al experiments show that some multigrid

methodsare attra tivein additiontothedeation method.

The major appli ation of this thesis is the Poisson equationwith a dis ontinuous

oe ient, whi h is derived from 2-D and 3-D bubbly ow problems. Most of the

performed numeri al experiments in thisthesis are based on thisequation. Both

sta-tionary andtime-dependent experiments are arried out to emphasize the theoreti al

results. Weshowthattwo-levelPCG methodsaresigni antlyfasterthantheoriginal

PCGmethod inalmost allexperiments. Hen e, omputationsinvolvedinbubblyows

Tweelaags Gepre onditioneerde Ge onjugeerde Gradiënten Methoden

met Toepassingen in Stromingsproblemen met Bellen

Jok M.Tang

Degepre onditioneerdege onjugeerdegradiënten(PCG)methodeiséénvandemeest

populaire iteratieve methoden voor het oplossenvan groots halige lineaire systemen,

waarbij de oë iëntenmatrix symmetris h en positief semi-deniet is. E hter, als

degepre onditioneerde oë iëntenmatrix sle htge onditioneerd is, danvertoont de

PCGmethodelangzame onvergentie. InplaatshiervankandetweelaagsePCG

meth-ode gebruikt worden die gebaseerd is op een tweelaagse pre onditioner. Deze

pre- onditioner elimineert de ee ten vande kleine en grote eigenwaarden van de

oë- iëntenmatrix,waardoorde tweelaagsePCG methode sneller onvergeert dande

oor-spronkelijke methode. Vele tweelaagse pre onditioners zijn bekend in de vakgebieden

van deatie, multirooster en domein de ompositie methoden. In dit proefs hrift

on-derzoeken wedeze pre onditioners nader, waar we onsvoornamelijk on entreren op

dedeatiemethode.

Welatentheoretis heeigens happenvandedeatiemethodezien,dieinzi htgeven

in deee tiviteit van dezemethode. Een ru iale omponentvan de deatie

pre on-ditioner is dekeuze vande proje tieve toren. Diversekeuzesworden beargumenteerd

en onderzo ht. Welaten zien datsubdomeinproje tieve toren, diegebaseerd zijnop

disjun te en stuksgewijs onstante ve toren, een van de beste keuzes zijn voor een

spe i iekeklassevan problemen.

Vervolgensonderzoeken wedetoepassingvande deatiemethode oplineaire

sys-temen waarbij de oë iëntenmatrix singulieris. Vers heidene wiskundig equivalente

variantenafgeleidvandeorigineledeatiemethodewordenbehandeld. Dezevarianten

zijnbestandtegendemogelijkesingulariteitvande oë iëntenmatrix. Verderworden

twee varianten bekeken die ges hikt zijn om kleinere lineaire systemen binnen de

de-atiemethodeoptelossenwaarbijdeGalerkinmatrixbetrokkenis. Nahetbehandelen

vandeimplementatieendee iëntievandedeatiemethode,latenweziendatdeze

Verder presenteren we een vergelijking tussen de deatie methode en andere

be-kende tweelaagse PCG methoden, waaronder de gebalan eerde Neumann-Neumann,

additief grof-rooster orre tie en multirooster methoden gebaseerd op symmetris he

en niet-symmetris he V- y li. Indien de parameters in de te bes houwen tweelaagse

pre onditionersgelijkzijn,kunnenweaantonendatdevers hillendemethodensterkaan

elkaargerelateerd zijn. Devergelijkingisverder ook uitgevoerd,waarbij detweelaagse

PCG methoden hun karakteristieke en geoptimaliseerde verzameling van parameters

aannemen. Numerieke experimenten laten zien dat sommige multirooster methoden

attra tief zijnnaast dedeatiemethode.

Debelangrijkstetoepassinginditproefs hriftisdePoissonvergelijkingmeteen

dis- ontinue oë iënt,hetgeenafgeleidisvan2-Den3-Dtwee-fasestromingsproblemen

metbellen. Demeeste vande uitgevoerdenumerieke experimenten zijngebaseerd op

deze vergelijking. Zowel stationaire alstijdsafhankelijkeexperimenten zijn uitgevoerd

omde theoretis he resultaten te onderbouwen. We latenzien datin bijnaalle

experi-menten de tweelaagse PCG methoden signi ant sneller onvergeren dande originele

PCG methode, waardoor de berekeningen voor twee-fase stromingen met bellen

After four years of intensive s ienti resear h, this do toral thesis is the absolute

rowningglory. Thesewordsofthanksare,in myopinion,a ru ialpartof thethesis,

asthisisthe onlywaytolet people whohave ontributeddire tly andindire tly ome

tothe frontandto thankthem ina deservingway. Onlywiththe help andsupportof

these very talented people has my PhD resear h su eeded and be ome the nished

produ tthat liesinfront ofyou now.

ThisPhD resear his partiallysponsoredbythe BSIK/BRICKSproje t, whi his a

six-years ienti resear hprogramthatisaninitiativeoftheDut hgovernment. Iam

thankfultoallofthe peoplebehind thisproje tfortheirnan ialsupport. Inaddition,

I would liketothank DelftUniversityof Te hnologyforsponsoringa largepart ofthe

thesisprinting ost.

The bulk of my gratitude goes out to my advisor, Prof. Kees Vuik. During my

Applied Mathemati s studies,Kees played aprominent rolebya ting asa motivating

andenthusiasti supervisorforvarioussubje tsandtheMS proje t. Iamverygrateful

tohim for onsequentlyoeringmeapositionofPhDstudent. Ienjoyedthepastfour

years,in whi hhe fun tionedasa terri supervisor, olleagueand ompanion. Ihave

learneda onsiderableamountfromhis knowledge,his onstru tive riti ism,and the

numerousdis ussionswehad. IamgratefultoKeesforthefa tthathegavemespa e

and freedom to do my do toral resear h, to attend ourses, tostart ooperations at

home and abroad, and to travel the world. It has been a privilege to attend several

onferen es together. In spiteofhis busy s hedule, hehas always been there, aswell

asgivenhisfull100per entinterestinmysu ess. IoweKeesforthefa tthathehas

alwaystakenthetimeandtroublereadthroughmypaperwork,ex hangemathemati al

ideas, helpget bugs outof program odes,answer mymany questionsandmore. His

enthusiasm, patien e and obliging manner led to various publi ations and a higher

standard for this thesis. Furthermore, he has inspired and motivated me in both the

easyanddi ulttimes,whi hensuredthatIbe ame amu hmore ompletes ientist.

Ihavehadatremendoustimewiththerenownedandprominentnumeri al analysis

groupundertheguidan eofProf. PietWesseling,andsubsequently Prof. Kees Vuik.

I feel honored to have been part of this group of talented olleagues, namely Fons

Prof. Kees Oosterlee, Jennifer Ryan, Guus Segal, Peter Sonneveld, Fred Vermolen,

and the many MS , PhD and postgraduate students of whom the list is too long to

re ount. I thank you all for the pleasant atmosphere in the department, the many

valuable onversations,thehilarious oee breaks,and,above all,thevariousamusing

outingsand dinners outside working hours. I wish ournumeri al analysis department

the verybest inthe future.

I would liketo name Diana Droog, our se retary ofthe numeri al analysis group,

separately. Sheis a true professional and extraordinarilygood at her job. These past

fewyears,shehas onstantly helpedme,and assistedmein variousan illary matters,

whi hensured that I ould fo us fullyon my resear hinstead ofon these matters of

lesserimportan e. Inaddition,sheis afantasti womanwhohasalwaysbrightenedup

ourdepartment. Iwanttotakethisopportunitytothankherforallherwork,support,

and heerful onversations.

During the whole of my studies and PhD work, I have had the fortune and the

privilegetoalwaysbeableto ountontwobuddies,namelySandervanVeldhuizenand

Jelle Hijmissen. In the past nine years, I have had mu h fun with and support from

them, bothon ana ademi and a so iallevel. Therefore, it ismore than self-evident

tometoaskthem toa t asmyparanymphs. I wanttothankthem up frontfortheir

assistan eandsupportat the timeI willbeupholdingmy do toral degree.

Iae tionatelythankthemembersofthedo toralexamination ommitteefor

read-ingthroughthisdo toral thesisaswellasprovidingvaluablesuggestionsandremarks.

I realize fullwell it is not easy toread this bulkywork in su h a short period of time.

Therefore,myutmostrespe tgoesouttothem. Moreover,IwouldliketothankS ott

Ma La hlanfor providingvariousparts of thethesis with riti al omments.

Further-more,Ihavere eivedusefultipsandvaluablehelp on erningthese a knowledgements

andthe propositions belonging tothis dissertation frommany people, espe iallyfrom

Merel Keijzer, JenniferRyan andFredVermolen. I am indebted toKar Fee, who has

ele troni allydevelopedasket hofthe overofthisthesisandwhoispreparedtoa t

asphotographerduringthePhD eremony. Finally,IthankOlof Borgwitfordesigning

the nalversion ofthe overof thedissertation.

In the past four years, I have had the privilege to spend two working visits at

TU Berlin in Prof. Reinhard Nabben's group. I am espe ially indebted to Reinhard

for oering me these opportunities. He is a very spe ial and enthusiasti s ientist I

have learned mu h from and withwhom I have had the honor of working on several

papers together. I thank him for the many pleasant onversations and for all the

time he invested in me. In addition, Reinhard was a fantasti host during my two

visits to Berlin. He made sure that I felt at home in both this metropolis and his

department, where his group be ame a real familyto me. For this, I have Veroni a

Twilling (se retary), Yogi Erlangga, Christian Mense, Elisabeth Ludwig and Sadegh

Jokar to thank as well. They have shown me that ondu ting resear h is also for a

greatpart onne ted tofriendship andhavingfun.

dis-several publi ations. Thanks to S ott, I have learned about the eld of multigrid as

wellastheglit hes inthe mathemati altheory inthis area. Inaddition,hehas helped

mesubstantiallywithvariouspapersbyprovidingmewith riti al ommentary onboth

the ontents and the use of English. Moreover, he has ontributed signi antly to

Chapters 7,9and 10 of this dissertation,for whi hI express mysin ere gratitudeto

him.

My former roommate and olleagueSander van der Pijl wasof great value and a

towerofstrength tomeinthestarting phaseofmyPhDresear h. Hefamiliarizedme

withthe world of bubblyowsand the urrent methods withsupplementary program

odes. ThankstoSander,Ihadasolidbasetoadvan emyworkinthespe i resear h

I was doing. Above all,he was a very so ial and ami able roommate with a dry wit

anda nete hnique in slime-volley.

During the lastyears, I had the pleasure ofsharing myo e withthe many MS

andPhDstudents. ApartfromSandervanderPijl,FangFangwastheonewithwhom

Idelightedinspendingmostofthetimewithintheo e. Irendermythanksforevery

ni e onversation,herhelpfulnessand patien eover the pastthree years.

It has been a great virtue to work together with several other professionals I will

name shortly. I would like tothank Prof. Bendiks Jan Boersma andEmil Coyajee for

theiruseful help and ourinteresting dis ussions. Furthermore, I thank them fortheir

lari ationsandformakingavailabletheirparallelprogram ode,theTe plotsoftware

pa kage, aswell asfor usingtheir omputer lusterfor al ulations. Moreover, I

ren-dermythankstoTijmenCollignon,PeterLu as,AllanGersborg-Hansen(DTU,

Den-mark), ChristopheVandeker khove (KULeuven, Belgium) andAuke Ditzel(MARIN,

Wageningen)forthe shortbutee tive ollaborationsinregardtotheemploymentof

two-levelPCG methods onvarious appli ationsrelevanttothem.

Ihavehadtheprivilegetoworkwithterri peopleonvariouspubli ations. Forthis,

InotonlytaketheopportunitytothankProf. Kees Vuik,Prof. ReinhardNabben,and

S ottMa La hlan,but IalsothankYogiErlanggaandElisabethLudwig. Additionally,

IthankProf. PietWesselingforhismanygoodsuggestionsand orre tionsinregardto

anumberofthesepapers. IalsoenjoyedthevariouspleasanttalksPietandIhadduring

boththe time of mygraduation and PhD. Lastly,I thank the editors and anonymous

referees who,withtheirsuggestionsandkeenobservations,have onsiderablyboosted

the standardof thepapers,and onsequently ofthisdissertation.

I owemu hgratitudetosystem administrators EefHartman and parti ularlyKees

Lemmens. They have always been there for me when omputer or network issues

arose andwhenIhadpressing questions on erning omputers andprogram odes. In

addition,theyhavetaughtmethefundamentalsofLinuxandmanysoftwarepa kages.

In the timeof my PhD, I had the opportunity to visit several onferen es in The

Netherlands (of whi h the Wouds hoten onferen es were the most splendid) and

abroad(su h asin China,Gree e, Italy,Poland andUSA). I havehad the privilegeto

meetmanyinspiringpeopleduringthesetrips. Ispe i allywanttomentionandthank

respe tively. Namely,IgivemythankstoProf. RobertBeauwens,whohandedmethe

beststudent paper prize. Indoing so,he motivated meto ontinue myPhD resear h

withmu henthusiasm andpassion.

The very best foreign adventure has without a doubt been my memorable visit

to the groupof Prof. Ni k Trefethen in Oxford. Thiswas on the invitation of Prof.

Gene Golub,who oered me the opportunity to broaden myhorizons in the world of

numeri al mathemati s. Gene was a brilliant and perfe t numeri al mathemati ian,

and,perhaps, thegreatest of thepast de ades. Therefore, he hasbeen mygreat idol

during the timeof my do toral resear h. It wasa privilege tospend aweek with him

in Oxford, during whi htime I dis overed thatGene is wonderful asa person as well.

Despitehisboundlessfameandrenown,heisoneofthemostself-ea ing,humorous,

involvedand ongenialpeopleIhavehadtheprivilegetomeetduringthelastfewyears.

The wayGene asso iated withme wasthe same ashe did withprominent s ientists.

He wasgenuinely interested in me, as a s ientistbut ertainly also asa person. The

way he motivated and stimulated people is truly brilliant. Furthermore, Gene was an

ex ellenthost, who reatedopportunitiesforfun a tivitiesbothwithin andoutsideof

working hours. He alsoshowed me that the world ofnumeri al mathemati sis wider

thansimple mathemati s,asdemonstratedbyourvisitstomusi alsandtohisfriends

inLondon. Gene'ssuddendeathat the endof 2007has notsurprisinglybeen aheavy

blow and di ult to grasp. The very greatest numeri al mathemati ian, and above

allanamazing person, hasevidently left us. Gene, youhave o upied a spe ial pla e

in my heart and shown me the beauty of numeri al mathemati s and life. Despite

not always having spoken as mu h about the resear h, your energy and enthusiasm

have, to a great extent, led to the realization of this thesis. In memory of Gene,

weorganizedthe GeneGolub DCSEsymposium in Delft withinthe frameworkof the

worldwide 'Gene-around-the-world' proje t. With this, I also thank fellow organizers

DianaDroog, Martin vanGijzen, Kees Vuik, andMarielba Rojas, who made sure this

symposiumbe ame agreat su ess.

Asmyresear hallowedmetolosemyselfinit,sodidtheedu ationIwasallowedto

providetostudentsand the oursesI ouldattendin ordertobroaden myexperien e.

For the various edu ational tasks, I give thanks to the expertise and help of Henri

Corstens, Fons Daalderop, Martin van Gijzen, Kees Vuik and Peter Wilders. For the

ourses,Ispe i allythankMerelKeijzerandEvelynvande Veenfortheireortsand

extraordinary lasses,whi hhavebrought myEnglishtoasigni antlyhigher level.

The daily lun h break, whi h was mostly in the aula, was always a heerful aair

and the moment we ould settle down every day at the o e. For this, I not only

openlythankmy olleaguesfromthe numeri alanalysisgroup, butalsotheonesfrom

thefthoor,inparti ularJelleHijmissenandMirjamterBrake. Furthermore,Ithank

thevarious olleaguesofotherdepartmentsfortheir onversationsoutsideofthelun h

breaks. Pertaining to this,I primarily thinkof the enjoyable annual DIAM outings, at

whi h olleaguesof theentire mathemati sdepartment ouldgettoknowea h other

the annual weekend event where PhD students in numeri al mathemati s from The

Netherlands and Flanders assemble to hange thoughts and have a lot of fun. After

my parti ipation at the su essful PhDays 2006 in Bérismenil (Belgium), I had the

pleasure to organize both PhDays 2007 in Baars hot and PhDays 2008 in De Haan

(Belgium). Forthis,IhumblythankArthurvanDam,YvesFrederix,LiesbethVanherpe

and Sander van Veldhuizen for o-organizing PhDays 2007, and Tijmen Collignon,

KatrijnFrederix, Ri ardodaSilva,MariaUgryumova,JorisVanbiervlietandon emore

LiesbethVanherpefor o-organizingPhDays2008. It wasanhonorandverydelightful

toorganize these latest twosu essfuleditions.

In on lusion, I would like to thank my dearest friends, a quaintan es and family

for their support and friendship over the past few years. I would like to spe i ally

thankmyparents Lisaand Cheung,who havebelievedin meun onditionallyandhave

investedinmefrommy hildhoodon. Theirinexhaustibleloveandsupporthavegiven

me the energy and strength to omplete my PhD after years of eort. I am very

proud ofthem, whi hmakes it unsurprisingthat I dedi ate this thesis tomy parents.

Additionally,Ithank mysister Loraforhersupport andfaithduring myentirelife.

Last but not least, I owemany thanks to my beloved SiuMei, who has been my

greatest support and an hor over the years. She has onstantly en ouraged me to

ontinuemydo toralresear h. SiuMeihasalwaysbeenthere formeandhasstoodby

methe entire time. I haveonly beenable towrite this dissertationwithher patien e,

heerfulnessandlove. Isin erelythankSiuMeiforthis. Sheisdenitelythebestthing

thathas ever happenedtome.

Jok Tang

Summary iv

Samenvatting vii

A knowledgements ix

1 Introdu tion 1

1.1 Ba kground . . . 1

1.2 Two-levelPre onditioned ConjugateGradientMethods . . . 3

1.3 BubblyFlow Problems . . . 3

1.4 S opeof the Thesis . . . 6

1.5 Outline ofthe Thesis . . . 6

1.6 Notation . . . 9

2 Iterative Methods 11 2.1 Introdu tion . . . 11

2.2 Basi IterativeMethods . . . 11

2.3 Conjugate GradientMethod . . . 13

2.4 Pre onditionedConjugate GradientMethod . . . 15

2.5 Further Considerations . . . 19

2.5.1 Pre onditioning . . . 19

2.5.2 Starting Ve tors andTerminationCriteria . . . 21

2.6 Appli ation toBubblyFlows . . . 22

2.7 Con ludingRemarks. . . 23 3 DeationMethod 25 3.1 Introdu tion . . . 25 3.2 Preliminaries . . . 25 3.3 Deated CGMethod . . . 28 3.4 Deated PCGMethod . . . 30

3.5 Properties oftheDeation Method . . . 33

3.5.2 Results foranArbitrary DeationSubspa e . . . 35

3.5.3 TerminationCriteria . . . 38

3.6 Appli ation toBubblyFlowProblems . . . 39

3.6.1 Results ofNumeri al Experiments . . . 40

3.7 Con ludingRemarks. . . 41

4 Sele tionofDeation Ve tors 43 4.1 Introdu tion . . . 43

4.2 Choi esofDeation Ve tors . . . 44

4.2.1 Approximated Eigenve torDeation . . . 44

4.2.2 Re y lingDeation . . . 45

4.2.3 SubdomainDeation . . . 46

4.2.4 Multigridand MultilevelDeationVe tors . . . 48

4.2.5 Dis ussionofDierent Approa hes . . . 48

4.3 Appli ation toBubblyFlows . . . 49

4.3.1 Preliminaries . . . 49

4.3.2 Inexa tEigenve torDeation . . . 51

4.3.3 Level-SetDeation Ve tors . . . 55

4.3.4 SubdomainDeation . . . 56

4.3.5 Level-Set-SubdomainDeation . . . 58

4.3.6 Numeri al Experiments . . . 61

4.3.7 AnalysisofSmall Eigenvalues . . . 64

4.4 Con ludingRemarks. . . 67

5 SubdomainDeationappliedto SingularMatri es 69 5.1 Introdu tion . . . 69

5.2 Preliminaries . . . 71

5.3 DeationVariants . . . 74

5.4 Theoreti al Comparisonof DeationVariants . . . 75

5.4.1 On theConne tion ofthe SingularandInvertible Matrix . . . 76

5.4.2 Comparisonof theDeated SingularandInvertible Matrix . . 77

5.4.3 Comparison of the Pre onditioned Deated Singular and In-vertible Matrix . . . 80

5.4.4 Comparisonof thePre onditioned DeatedSingularMatri es 82 5.5 Appli ation toBubblyFlows . . . 83

5.5.1 Results ofICCGand DICCGwithVariant 5.2 . . . 83

5.5.2 Results ofthe Comparison between Variants5.1and 5.2 . . . 84

5.6 Con ludingRemarks. . . 86

6 Comparison ofTwo-Level PCG Methods Part I 87 6.1 Introdu tion . . . 87

6.2 Two-Level PCG Methods . . . 90

Multi-6.2.2 GeneralLinearSystems . . . 91

6.2.3 Denitionof the Two-LevelPCG Methods . . . 92

6.2.4 Aspe ts ofTwo-Level PCGMethods . . . 96

6.3 Theoreti al Comparison . . . 99

6.3.1 Spe tralAnalysisof theMethods . . . 99

6.3.2 Equivalen es betweenthe Methods . . . 103

6.4 Numeri al Comparison . . . 105

6.4.1 Setup ofthe Experiments . . . 105

6.4.2 Experiment usingStandard Parameters . . . 106

6.4.3 Experiment usingIna urate GalerkinSolves . . . 106

6.4.4 Experiment usingSevere TerminationToleran es . . . 108

6.4.5 Experiment usingPerturbed Starting Ve tors . . . 112

6.4.6 Further Dis ussion . . . 114

6.5 Con ludingRemarks. . . 115

7 Comparison ofTwo-Level PCG Methods Part II 117 7.1 Introdu tion . . . 117

7.2 Two-Level PCG Methods . . . 118

7.3 Spe tralProperties ofMG . . . 120

7.3.1 UnitEigenvaluesofthe MG-Pre onditionedMatrix . . . 120

7.3.2 Positive Denitenessof theMG pre onditioner . . . 122

7.4 Comparisonof aSpe ialCase ofMG andDEF . . . 125

7.5 Ee tof RelaxationParameters . . . 126

7.5.1 AnalysisofS alingRelaxation . . . 127

7.5.2 Optimal Choi eof . . . 128

7.6 Symmetrizingthe Smoother . . . 130

7.7 Numeri al Experiments . . . 132

7.7.1 1-DPoisson-likeProblem. . . 132

7.7.2 2-DBubblyFlow Problem . . . 134

7.8 Con ludingRemarks. . . 135

8 E ien yand Implementationofthe DeationMethod 137 8.1 Introdu tion . . . 137

8.2 Computationswiththe DeationMatrix . . . 139

8.2.1 Constru tion ofAZ . . . 139

8.2.2 Constru tion ofE . . . 139

8.2.3 Cal ulationofPy and P T y . . . 140

8.3 E ientSolution ofGalerkinSystems . . . 140

8.3.1 GalerkinSystems withinDICCG1 . . . 141

8.3.2 GalerkinSystems withinDICCG2 . . . 141

8.3.3 Comparisonof GalerkinMatri es. . . 144

8.3.4 DeationProperties fora SingularGalerkinMatrix . . . 144

8.5.1 Results forthe DeationMethodwith E ientImplementation 147

8.5.2 Comparisonof DICCG1andDICCG2 . . . 150

8.5.3 Comparisonof DICCG2andADICCG2 . . . 151

8.6 Con ludingRemarks. . . 154

9 Comparison ofDeationand Multigridwith Typi al Parameters 157 9.1 Introdu tion . . . 157

9.2 Numeri al Methods . . . 159

9.2.1 DeationApproa h . . . 159

9.2.2 MultigridApproa hes . . . 160

9.3 Implementationand ComputationalCost. . . 166

9.3.1 Cost ofDeation . . . 166

9.3.2 Cost ofMultigrid . . . 167

9.3.3 SingularityofCoe ientMatrix . . . 169

9.3.4 Parallelization . . . 169

9.3.5 Implementation . . . 169

9.4 Numeri al Experiments . . . 170

9.4.1 VaryingDensityContrasts . . . 170

9.4.2 VaryingBubbly Radii . . . 171

9.4.3 VaryingNumber ofBubbles . . . 172

9.4.4 VaryingNumber ofGrid Points. . . 173

9.4.5 Di ultTestProblem . . . 174

9.5 Con ludingRemarks. . . 175

10 BubblyFlow Simulations 177 10.1 Introdu tion . . . 177

10.2 Mathemati al Modelofthe BubblyFlow . . . 177

10.3 BubblyFlow Simulations . . . 179

10.3.1 RisingAir Bubble inWater . . . 180

10.3.2 Falling Water Droplet in Air . . . 182

10.3.3 TwoRising andMergingAir Bubblesin Water . . . 183

10.4 Con ludingRemarks. . . 185

11 Con lusions 187 11.1 Con lusions . . . 187

11.2 Future Resear h . . . 189

A Basi Theoreti al Results 191 B DeterminationofBubbles fromthe Level-Set Fun tion 199 C More Insightsinto Deationapplied toSingularCoe ient Matri es 203 C.1 Theoreti al Results . . . 203

D E ient ImplementationofDeationOperations 209 D.1 E ientConstru tion ofS AZ andS E in 2-D . . . 209

D.1.1 Number ofNonzeroEntries in AZ . . . 210

D.1.2 Treatmentofthe Dierent Cases . . . 211

D.1.3 Constru tion ofS AZ . . . 212 D.1.4 Constru tion ofS E . . . 212 D.2 E ientConstru tion ofS AZ andS E in 3-D . . . 213

D.2.1 Number ofNonzeroEntries in AZ . . . 213

D.2.2 MatrixS AZ forEight Blo ks . . . 214

D.2.3 MatrixS AZ for27 Blo ks . . . 214

D.2.4 MatrixS AZ withVariable NumberofBlo ks. . . 216

D.2.5 Constru tion ofS E . . . 216

E FlopCounts forthe Deation Method 217 E.1 DeationOperations . . . 218

E.2 ICCG,DICCG1and DICCG2 . . . 220

F ParallelVersion ofthe DeationMethod 223 F.1 TraditionalParallelPre onditioners . . . 224

F.2 ParallelDeation . . . 224

G Two-Level PCG Methods applied toPorous-Media Flows 227 G.1 Problem Setting . . . 227

G.2 Experiment usingStandardParameters . . . 228

G.3 Experiment usingIna urate CoarseSolves . . . 230

G.4 Experiment usingSevere TerminationToleran es . . . 231

G.5 Experiment usingPerturbedStarting Ve tors . . . 234

H DICCG Variants appliedtoBubbly Flow Simulations 237 H.1 Simulation1: RisingAir Bubblein Water . . . 237

H.2 Simulation2: FallingWater Droplet inAir . . . 238

I Comparison ofDeationand Multigridfor aSpe ial Case 243

Bibliography 245

List ofPubli ations 261

1.1 A dropletsplash: anexampleofa two-phasebubbly ow problem. . . 4

1.2 Geometry of some stationary bubbly ows onsidered in this thesis

(m=numberof bubbles,s=radius ofthe bubbles). . . 4

1.3 A two-phasebubbly owwith thephases airand water. . . 5

3.1 Deation subdomains withk =3, whi hare hosen independently of

the densitygeometry ofthe bubblyow. . . 40

3.2 NormoftherelativeresidualsduringtheiterationsofICCGandDICCG k. 42

4.1 Representation ofthe2-Dsubdomainswithk =3inasquaredomain,

, onsistingof n=64 2

grid points. . . 47

4.2 A2-Dexampleofabubblyowproblem withm =5andtwodierent

situationsfor theperturbations,fÆ

i

g. . . 54

4.3 A2-Dbubblyow problemwithm=5,whi hillustratesthelevel-set,

subdomain andlevel-set-subdomaindeation te hnique.. . . 60

4.4 EigenvaluesofM 1 AandM 1 P W S

AforS-DICCG,appliedtothe

Pois-sonproblem with=1andn =16 2 .. . . 62 4.5 EigenvaluesofbothM 1 AandM 1 P W L A orrespondingtoL-DICCG 4,

for thePoisson problem with=10 6 and n=16 2 . . . 65 4.6 Eigenvalues of M 1 A and M 1 P W S

A orresponding to S-DICCG, for

the Poissonproblem with=10 6 andn=16 2 . . . 66 4.7 EigenvaluesofbothM 1 AandM 1 P W LS A orrespondingtoLS-DICCG,

for thePoisson problem with"=10 6

andn =16 2

. . . 67

5.1 Residuals of ICCG, DICCG 2 3

and DICCG 3 3

with Variant 5.2, for

the test asewithn =32 3

and=10 3

. . . 85

6.1 Relativeerrorsduringtheiterativepro ess,forthebubblyowproblem

withn =64 2

,and`standard' parameters. . . 107

6.2 Relativeerrorsduringtheiterativepro ess,forthebubblyowproblem

withparameters n=64 2

and k =8 2

,andaperturbed Galerkinmatrix

inverse,E

e

1

6.3 Relativeerrors duringtheiterativepro essforthebubblyowproblem

withparameters n =64 2

;k =8 2

, andvarioustermination riterion. . 111

6.4 Relativeerrors duringtheiterativepro essforthebubblyowproblem

withn =64 2

;k =8 2

,and perturbed starting ve tors. . . 113

7.1 Regionswhere  MG < DEF (Regions A 1 and A 2 )and  DEF < MG (Regions B 1 and B 2 ), for arbitrary M 1 =  M 1 , when Z onsists of eigenve torsofM 1

A. Thetwo onditionnumbersareequalalongthe

dotted and dotted-dashedlines. . . 127

7.2 Eigenvaluesasso iatedwithDEF andMG forthe test aseswithk =

20 aspresented inTable7.3. . . 134

8.1 Visualizationofthe resultsforthetest problem withm =27. . . 149

8.2 Results forthe test problem withm=27 for varyinggridsizes. ICCG

andDICCG1 withboth =5and =10 are presented. . . 149

8.3 Resultsforthetest problem withm=27 forvaryingdensity ontrast, .150

8.4 CPU timeofDICCG1and DICCG2for thetest problem withm =27

andvariousgrid sizes. . . 151

8.5 Convergen e ofthe residuals during aninner solve at oneiteration of

ADICCG2 25 3

(Test Problem 2). Theplots are similarfor the other

outeriterationsofthesametest ase,sin eoneappliestheina urate

solvestothe Galerkinmatrix, E. . . 153

8.6 Convergen eoftheresidualsoftheouteriterationsfromDICCG2 25 3

andADICCG2 25 3

(Test Problem2). . . 154

8.7 ResidualplotsoftheouteriterationsfromADICCG2 25 3

(Test

Prob-lem 1).. . . 155

10.1 Evolution of a risingbubble in water. Parameters l and t denote the

timestep anda tual time, respe tively. . . 181

10.2 ResultsforICCG,DICCG 20 3

andBoxMG-CGforthepressuresolves

during the real-lifesimulationwitha risingairbubble inwater. . . 181

10.3 Evolutionofa fallingdroplet inair. . . 182

10.4 ResultsforICCG,DICCG 20 3

andBoxMG-CGforthepressuresolves

during the real-lifesimulationofa fallingwater dropletin air.. . . 183

10.5 Evolutionoftworisingair bubblesinwater. . . 184

10.6 Results forDICCG 25 3

andBoxMG-CGforthe pressure solveduring

the real-life simulation with two rising air bubbles in water. ICCG is

omitted in these results, be ause it is not ompetitive withthe other

two methods. . . 184

B.1 A 2-D bubbly ow problem with m = 3 showing the appli ation of

D.1 Domain dividedinto ninesubdomains (k=9),so that ea h

subdo-main orrespondsto exa tlyonegroup. . . 210

D.2 Casesof gridpoints involvedinthe groups ofS

AZ

. . . 210

D.3 Cases of grid points involved in E :=Z T

AZ, denoted by E1, E2and

E3,whosevalues inAZ ontribute toe

1 ;e 2 ande 3 ,respe tively. . . . 213 D.4 Treatment ofBlo k 1. . . 214

G.1 Geometryof theproje tion ve tors andpie ewise- onstant oe ient

in theporous-mediaow. . . 228

G.2 Relativeerrors duringthe iterativepro ess,fortheporous-media

prob-lem withn =55 2

, k =7, and`standard' parameters. . . 230

G.3 Relative errors duringthe iterative pro ess forthe porous-media

prob-lem withn =55 2 ;k =7 andE

e

1 ,where a perturbation =10 8 is taken. . . 232

G.4 Relative errors duringthe iterative pro ess forthe porous-media

prob-lem withn =55 2

;k =7, andtermination toleran e Æ=10 16

.. . . . 233

G.5 Relative errors in the A norm during the iterative pro ess for the

porous-media problem with n = 55 2

;k = 7 2

, and perturbed starting

ve tors with =10 5

. Theplot ofthe relative errors in the2 norm

is omitted,sin e thetwoplotsare approximatelythe same. . . 235

H.1 Evolutionofthe risingbubblein water inthe rst 250timesteps. . . 238

H.2 Results for ICCG,DICCG1 10 3

andDICCG2 20 3

for the simulation

witha risingairbubble in water. . . 239

H.3 Evolutionofthe fallingdropletinair in therst 250 timesteps. . . . 240

H.4 Results for ICCG,DICCG1 10 3

andDICCG2 20 3

for the simulation

witha fallingwater dropletin air. . . 241

I.1 Fun tion  i (2  i )for i 2[0;2℄. . . 244

1.1 Notation forstandardmatri esand ve torswhere ;; 2N. . . 9

2.1 Results for ICCGapplied to 2-D bubbly ow problems. `# It' means

thenumberofrequirediterations,and`CPU'isthe orresponding

om-putationaltime inse onds. . . 22

2.2 Results forICCGappliedto3-Dbubbly ow problems. . . 23

3.1 Number of iterations for ICCG and DICCG k for various number of

bubbles,m,anddeation ve tors,k. . . 41

4.1 ResultsforthePoissonproblemwith=1andn=16 2

. `#It'means

the numberofrequired iterationsfor onvergen e. . . 61

4.2 Resultsforthe Poissonproblemwithm=5,=10 6

,andvaryinggrid

sizes,n. `#It' means thenumberof required iterations,and`CPU'is

the orresponding omputational timein se onds. . . 62

4.3 Results for the Poisson problem with m = 5, n = 64 2

, and varying

density ontrast,. . . 63

4.4 Results for the Poisson problem with  = 10 6

, n = 64 2

, and varying

numberof bubbles,m. . . 64

5.1 Corresponding matri esofthe proposeddeation variants. . . 75

5.2 Number of iterationsforICCG and DICCG(Variant 5.2) tosolve the

linearsystem 

A x =b with invertible 

A,forthe test asewithm=2 3

,

=10 3

, ands =0:05. . . 84

5.3 NumberofiterationsforICCGandDICCG(Variant5.2)tosolve 

Ax =

b,for thetest asewithm=3 3

,=10 3

,and s=0:05. . . 85

5.4 NumberofiterationsforDICCG (bothVariant5.1andVariant5.2)to

solveAx =b (withasingularA)and  Ax =b(with aninvertible  A)for m=2 3 . . . 85

6.1 List of methods that are ompared in this hapter. The operator of

ea hmethod anbeinterpreted asthepre onditionerP,givenin(6.1)

with A = A. Where possible, referen es to the methods and their

implementationsare presentedin the last olumn. . . 96

6.2 Choi es of parameters forea h method, used in the generalized

two-level PCGmethod asgiven inAlgorithm 7. . . 97

6.3 Extra omputational ost periteration of the two-levelPCG methods

ompared toPREC.IP =innerprodu ts, MVM=matrix-ve tor

mul-tipli ations, VU =ve torupdates and GSS=Galerkin systemsolves.

Note thatIP holdsforsparseZ andMVM holdsfordense Z. . . 99

6.4 Numberofrequired iterationsfor onvergen eandthe 2 normofthe

relativeerrorsofallmethods,forthebubblyowproblemwithn =64 2

,

and `standard' parameters. `NC' means no onvergen e within 250

iterations. . . 106

6.5 Computational ost within the iterations in terms ofnumber ofinner

produ ts(`IP'),ve torupdates(`VU'),Galerkinsystemsolves(`GSS'),

andpre onditioningstepwithM 1

(`PR'),forthebubblyowproblem

withn =64 2

,k =8 2

,and`standard'parameters. . . 108

6.6 Numberofrequired iterationsfor onvergen eandthe 2 normofthe

relativeerrorsofallmethodsforthebubblyowproblemwith

parame-tersn =64 2

andk =8 2

,andaperturbedGalerkinmatrixinverse,E

e

1

,

isusedwithvaryingperturbation . `NC'meansno onvergen ewithin

250iterations. . . 110

6.7 Numberofrequired iterationsfor onvergen eandthe 2 normofthe

relative errors of all methods, for the bubbly ow problem with

pa-rameters n =64 2

andk =8 2

. Various termination toleran es, Æ,are

tested. . . 110

6.8 Numberofrequired iterationsfor onvergen eandthe 2 normofthe

relativeerrors ofsomemethods,forthe bubblyow problemwithn =

64 2

;k = 8 2

, and perturbed starting ve tors. An asterisk (*) means

thatanextra uniquenessstepis appliedinthattest ase.. . . 114

7.1 Listof two-levelPCG methodsthat are ompared inthis hapter. . . 120

7.2 Test ases orresponding todierentregions aspresentedin Figure7.1.132

7.3 Resultsoftheexperimentwithtest asesaspresentedforthe

Poisson-likeprobleminTable7.2. Theresultsarepresentedintermsofnumber

ofiterations,# It., and onditionnumber,. . . 133

7.4 Numberofrequired iterationsfor onvergen eandthe 2 normofthe

relative errors of 2L-PCG methods, for the bubbly ow problem with

n =64 2 and  M 1 =M 1

. PRECrequires 137iterationsand leadsto

7.5 Numberofrequired iterationsfor onvergen eandthe 2 normofthe

relative errors of 2L-PCG methods, for the bubbly ow problem with

n =64 2 and M 1 =  M 1 +  M T  M T A  M 1 . PREC requires137

iterationsand leadstoarelative errorof4:610 7

. . . 135

8.1 Corresponding matri esofthe proposeddeation variantsin Chapter5.140

8.2 Convergen e resultsforICCG and DICCG1 for alltest problems with

n =100 3

and=10 3

. `#It' meansthenumberofrequired iterations

for onvergen e, and `CPU' means the total omputational time in

se onds. . . 148

8.3 Results for DICCG2 and ADICCG2 to solve Ax = b with n = 100 3

,

orresponding to Test Problem 1. ICCG requires 390 iterations and

37.0se ondsforTestProblem1and543iterationsand177.6se onds

forTestProblem2. `#It'=numberofiterationsoftheouterpro ess,

`CPU'=therequired omputingtime(inse onds)in ludingthesetup

timeof themethods, `NC'=no onvergen e within250iterations. . 152

9.1 Convergen e results for the experiment with n = 64 3

, m = 2 3

, s =

0:05,andvaryingthe ontrast,. `CPU',`#It.' and`RES'denotethe

total omputingtime, numberofiterationsor y les,andthea ura y

of the solution measured as the relative norm of the exa t residuals,

respe tively. . . 171

9.2 Convergen eresultsfortheexperimentwithn=64 3 ,m=2 3 ,=10 3 ,

andvaryingthe radiusofthe bubbles,s. . . 172

9.3 Convergen e results for the experiment with n = 64 3

, s =0:05,  =

10 3

,andvarying thenumber ofbubbles,m. . . 173

9.4 Convergen e results for the experiment with m = 2 3

, s = 0:05,  =

10 3

,andvarying thetotalnumber ofdegrees offreedom, n. . . 174

9.5 Convergen e results for the di ult test problem. The following

pa-rameters are kept onstant: n = 128 3 , m =3 3 , s =0:025,  =10 5 , andÆ =10 8 . . . 175

D.1 Explanationofthe variables. . . 211

D.2 Casesinvolvedin the blo ksforeight subdomainsin 3-D.. . . 215

D.3 Casesinvolvedin Blo ks19fork =27 inthe 3-D ase. . . 215

E.1 Results ofop ounts forstandardoperations. . . 217

G.1 Numberofrequirediterationsfor onvergen eofallproposedmethods,

fortheporous-mediaproblemwith`standard'parameters. The2 norm

G.2 Total omputational ost withinthe iterations in terms ofnumber of

inner produ ts (`IP'), ve tor updates (`VU'), Galerkin system solves

(`GSS'),pre onditioning step withM 1

(`PR'),forthe porous-media

problem withn=55 2

,k =7,and `standard'parameters. . . 229

G.3 Numberofrequired iterationsfor onvergen eandthe 2 normofthe

relative errors of all methods, for the porous-media problem with

pa-rametersn =55 2

andk =7. AperturbedGalerkinmatrixinverse,E

e

1

,

is usedwitha varyingperturbation, . . . 231

G.4 Numberofrequired iterationsfor onvergen eandthe 2 normofthe

relative errors of all methods, for the porous-media problem with

pa-rameters n = 55 2

and k =7. Various termination toleran es, Æ, are

tested. . . 234

G.5 Numberofrequired iterationsfor onvergen eandthe 2 normofthe

relative errors of some methods, for the porous-media problem with

parametersn=55 2

;k =7,andperturbedstartingve tors. Anasterisk

1 Conjugate Gradient(CG) solvingAx =b . . . 14

2 Pre onditionedCG (PCG)solving Ax=b (Original Variant) . . . 17

3 Pre onditionedCG (PCG)solving Ax=b (Pra ti al Variant) . . . . 17

4 Deated ConjugateGradient(DCG) solvingAx =b . . . 30

5 Deated PCG(DPCG) solvingAx =b (Original Variant). . . 31

6 Deated PCG(DPCG) solvingAx =b (Pra ti al Variant) . . . 32

7 GeneralizedTwo-LevelPCG Methodforsolving Ax=b. . . 96

8 Computationof Py . . . 138

Chapter

1

Introdu tion

1.1 Ba kground

The fo us of this thesis is on the numeri al solution of the linear partial dierential

equations(PDEs) resultingfromthe mathemati almodeling ofphysi al systemsand,

inparti ular, bubblyows. Weassumethatthese PDEshavealready beendis retized

in a sensible manner through the use of nite dieren es, nite volumes or nite

elements. Ourprimary fo usis one ient solutionof linearsystems,ofthe form

Ax =b; A2R nn

; n 2N; (1.1)

thatarisefromsu h dis retizations,wheren is thenumberofdegreesof freedomand

is alledthe dimensionof A. InEq.(1.1), the oe ient matrix,A,is assumed tobe

real,symmetri , andpositivesemi-denite(SPSD), i.e.,

A=A T ; y T Ay 08y 2R n ;

and has d zero eigenvalues with orresponding linearly independent eigenve tors. If

d >0, then A is singular. To guarantee that Eq. (1.1) is onsistent, the right-hand

side,b,is presumed to bein the rangeof A,i.e.,b2R(A) whereR(A):=fy 2R n

:

y =Aw forw 2R n

g. Thus, thenext assumptionholds throughout thisthesis.

Assumption 1.1. The oe ient matrix, A, is SPSD and has d zero eigenvalues.

Moreover, the linearsystem(1.1) is onsistent.

The nullspa e ofA is dened asN(A):=fw 2R n :Aw =0 n g,where 0 n is the

all-zero ve torwithn entries. Then, N(A) is theorthogonal omplement ofthe olumn

spa e ofA, i.e., N(A) =R(A) ?

. As a onsequen e, the linear system (1.1) is only

onsistentif b T

w =0is satisedfor allw 2N(A).

Linearsystem(1.1)is typi allylarge, sparse,and ill- onditioned. Thatmeansthat

urrent problems of interest involve millions of degrees of freedom, a xed number

of nonzero entries per row and olumn of A, and ondition number of A, denoted

in reases, respe tively. In this thesis, we denote by 

i

(B) (or, shortly, 

i

) the i-th

eigenvalue of anarbitrary symmetri matrix, B 2R nn

,where the set f

i

g is always

ordered in reasingly (unless otherwise stated), i.e., 

1   2  :::   n . This set, f i

g, is alled the spe trum of B and is denoted as (B). If B is SPSD, then its

(spe tral) ondition number is dened as the ratio of the largest and the smallest

nonzeroeigenvalues,i.e.,

(B):=  n  d+1 :

Thelinearsystem(1.1) an besolvedusingdire tmethods. Mostofthese solvers

generally involve expli it fa torization of (permutations of) A into a produ t of a

lowerandanuppertriangularmatrix. Important advantagesofdire t solversaretheir

robustness and generalappli ability. However, the bottlene k of dire t solvers is that

the matrix fa tor is often signi antly denser than A. On the one hand, this might

lead to an ex essive amount of omputations and, on the other hand, it might lead

to insu ient memory to form and store matrix fa tors. Therefore, dire t methods

are typi ally prohibitively expensive and in some ases impossible, even with the best

available omputing power.

Instead of dire t solvers, iterative methods are more attra tive touse to nd the

solutionof(1.1). Inthis ase,bothmemory requirementsand omputing time anbe

redu ed, espe iallyif A is large and sparse. Moreover, these methods are mandatory

for some numeri al dis retization methods, where A is not expli itly available. The

term `iterativemethod' refers toa widerange ofte hniques thatuseiterates, or

su - essive approximations, to obtain more a urate solutions toa linearsystem at ea h

iteration step. Krylovsubspa e iterative methods, espe ially the Conjugate Gradient

(CG)method ofHestenesandStiefel, areprominentiterativemethods tosolve(1.1).

Inthesemethods,the losestapproximationtothe solutionof(1.1)isfound ina

sub-spa e whosesize is iterativelyin reased. The onvergen e ofthese methods depends

highlyon (A), whi h again typi ally grows as the problem size in reases. To avoid

thein reaseiniterations,itis ommonpra ti etomodifytheKrylovsubspa emethod

inthehopesofredu ingthedi ultiesinsolvingthegivensystem. Ifthisisappliedto

CG,thentheresultingmethod is alledthepre onditionedConjugateGradient(PCG)

method. In this ase,(1.1) is multiplied by apre onditioner, M 1

, hosen to redu e

the ondition number of the iteration matrix from (A) to (M 1 2 AM 1 2 ), whi h is equivalent to (M 1

A). The resulting pre onditioned system that should be solved

reads M 1 Ax=M 1 b; (1.2)

whereM is assumedtobe symmetri andpositivedenite(SPD), i.e.,

M =M T ; y T My >08y 6=0 n :

PCG is more ee tive than original CG for many problems of interest. When M 1

y

M 1

A. Thatis,if(M 1

A)issigni antlylessthan (A),orifthe spe trumismore

lusteredthanthatoftheoriginalmatrix,we anexpe tsigni antlyfeweriterationsto

beneeded. However, even withsophisti atedpre onditioners,su h aspre onditioners

basedonin ompletefa torizations,(M 1

A)mightstillbe omelargerastheproblem

size or oe ient ratio in the original PDEs in reases. In this ase, PCG may suer

fromslow onvergen e due tothe presen eof unfavorableeigenvalues in(M 1

A).

1.2 Two-level Pre onditioned Conjugate Gradient Methods

Inadditiontoatraditionalpre onditioner,M 1

,ase ondkindofpre onditioner anbe

in orporatedtoimprovethe onditioningofthe oe ientmatrixevenfurther,sothat

theresultingapproa hee tivelytreats the ee tofallunfavorable eigenvalues. This

ombined pre onditioning is known as `two-level pre onditioning', and the resulting

iterative method is alled a `two-level PCG (2L-PCG) method'. In this ase, CG,

in ombination with a pre onditioner based on a multigrid (MG) method or domain

de ompositionmethod (DDM), an beregarded as a2L-PCG method, sin emostof

thesemethods relyonpre onditioningontwolevels. Thesepre onditionershavebeen

known fora longtime, dating ba kat least tothe 1930s.

The main fo us of this thesis is on the 2L-PCG method whosetwo-level

pre on-ditioner is based on a deation te hnique. The resulting method is often alled the

deation method and was introdu ed independently by Ni olaides and Dostal in the

1990s.

1.3 Bubbly Flow Problems

Themain appli ation ofthis thesis istwo-phasebubbly ows,asin Figure1.1.

Com-putationofthese owsis avery a tiveresear htopi in omputational uiddynami s

(CFD).Understandingthedynami sandintera tionofbubblesanddropletsinalarge

varietyofpro essesinnature,engineering,andindustryare ru ialfore onomi allyand

e ologi ally optimized design. Bubbly ows o ur, for example, in hemi al rea tors,

boiling,fuel inje tors, oating,and vol ani eruptions.

Two-phaseowsare ompli atedtosimulate,be ausethegeometryoftheproblem

typi allyvarieswithtime,andtheuidsinvolvedhaveverydierentmaterialproperties.

A simple example is thatof air bubbles in water, wherethe densities vary bya fa tor

of about 800. In this thesis, we onsider both stationary and time-dependent bubbly

ows,where the omputational domainisalwaysaunitsquareorunit ube lledwith

auidtoa ertain height. The bubblesanddroplets inthe domainare always hosen

su hthattheyare lo atedinastru turedwayandhaveequalradius,s,atthestarting

time. Typi al3-Dtest problems, onsideredin thisthesis, are depi tedin Figure1.2.

2-D test problems are always based on se tions of these 3-D domains. Throughout

thisthesis, lengthsare typi allygivenin entimeters ( m).

in-Figure1.1: Adropletsplash:anexampleofatwo-phasebubblyowproblem.

(a) m=1ands=0:1. (b) m=8ands=0:05. ( )m=27ands=0:025.

Figure 1.2: Geometry of some stationary bubbly ows onsidered in this thesis (m = number of

bubbles,s=radiusofthebubbles).

using operator-splitting te hniques. In these s hemes, equations for the velo ity and

pressure are solvedsequentially at ea h timestep. Inmany popular operator-splitting

methods, the pressure orre tion is formulated impli itly, requiring the solution of a

linearsystem(1.1) at ea htimestep. Thissystemtakesthe form ofa Poisson

equa-tion with dis ontinuous oe ients (also alled the `pressure(- orre tion) equation')

andNeumann boundary onditions, i.e.,

(

r



1 (x) rp(x)



= f(x); x2;  n p(x ) = g(x ); x2; (1.3)

where ;p;;x, and n denote the omputational domain, pressure, density, spatial

oordinates,andtheunitnormalve tortothe boundary,,respe tively. Right-hand

sidesf andg followexpli itlyfromtheoperator-splittingmethod,whereg issu hthat

massis onserved, leadingtoasingularbut ompatiblelinearsystem(1.1) 1 . We dene n x ;n y and n z

as the number of degrees of freedom in ea h spatial

dire tion, so that n = n x n y n z

. In this thesis, we perform the omputations on a

uniformCartesiangridwithn

x =n

y =n

z

. Furthermore,we onsidertwo-phasebubbly

1

owswith,for example,air andwater.

Inthis ase, ispie ewise onstant witharelatively large ontrast:

=

(

 0 =1; x2 0 ;  1 ="; x2 1 : (1.4)

Forowswithwaterandair,thedensity ontrast,denedas:=  0  1 =" 1 ,is10 3 ,

see Figure1.3. Inthis ase, 

0

is water, the main uid ofthe ow around the m air

bubbles,and

1

is theregion insidethe bubbles.

10

−3

10

−3

10

0

10

−3

10

−3

10

−3

Composition

water

air

air

air

air

air

Density

Figure1.3: Atwo-phasebubblyowwiththephasesairandwater.

Solving linear system (1.1), that is a dis retization of (1.3), within an

operator-splittingapproa hhaslongbeenre ognizedasa omputationalbottlene kinuid-ow

simulation,sin eittypi ally onsumesthebulkofthe omputingtime. Whether

nite-dieren e, nite-element,ornite-volumete hniques are usedtodis retize (1.3),the

resulting matrix is sparse, but with a bandwidth of n

x n

y

, in lexi ographi al ordering.

Foradis retizationona ubewith100degreesoffreedominea hdire tion,thismeans

that A is of dimension n = 10 6 , with bandwidth n x n y = 10 4 . It is well-known that

dire tsolutionte hniqueswithoutreorderingrequireanumberofoperationsthats ale

asn 7

3

forabandedCholeskyde omposition,O(10 14

)operationsintheexampleabove.

Thus,here,we onsiderthesolutionof (1.1)usingiterativete hniques. StandardPCG

methodsare notsuitable,sin etheyexhibitastrongsensitivitytothedensity ontrast

andgrid size. Thereis a realneed for two-levelpre onditioningin order toa elerate

the onvergen e of the iterative pro ess of PCG. Hen e, we apply 2L-PCG methods

tosolve(1.1).

Next,dene1

p

astheall-oneve torwithpentries. Then,thefollowingassumption

holds in our bubbly ow problem, whi h follows impli itly from the above problem

onstantatea hboundary;weusethefollowingboundary onditionsinthe3-D ase:

8

>

<

>

:

 n p(x )j x=0 =  n p(x )j x=1 = 1;  n p(x )j y=0 =  n p(x )j y=1 = 1;  n p(x )jz=0 =  n p(x )jz=1 = 1:

setting.

Assumption1.2. Inbubbly ow problems, we assumethat A is asingular M-matrix,

andthe equationsA1

n =0 n andb T 1 n =0aresatised.

A symmetri M-matrix is a square SPSD matrix whose o-diagonal entries are

less than or equal to zero. A ording to [15℄, d = 1 holds for a singular M-matrix,

A. In other words, we have rank A = n 1 and dim N(A) = 1, where rank B and

dim B denote the rank of matrix B and the dimension of subspa e B, respe tively.

Note that Assumption 1.1 follows immediately from Assumption 1.2, sin e (1.1) is

always onsistent. Forb =0

n

,this is trivial,and, forb 6=0

n

,b ?N(A)= spanf1

n g

resultingin thefa t thatb 2R(A). Therefore, althoughA is singular,(1.1) is always

onsistent andan innite number ofsolutions exists. Due tothe Neumann boundary

onditions, the solution, x, of(1.1) is xed up to a onstant, i.e., if x

1 is a solution thenx 1 + 1 n

is alsoa solutionof(1.1), where 2Ris any onstant. Thissituation

presents no real di ulty, sin e pressure is a relative variable, not an absoluteone in

the operator-splitting methods.

1.4 S ope of the Thesis

Thisthesis dealswitha elerationof PCG usingtwo-levelpre onditioningin order to

solve linear systems with an SPSD oe ient matrix. The main 2L-PCG method is

the deation method. In the literature, mu h is known about applying the deation

method to linear systems with invertible oe ient matri es and to problems with

xed and known density elds. In this thesis, we generalize it to linear systems with

singular oe ient matri esand to problems wherethe density eldvaries or annot

be des ribed expli itly. Moreover, we investigate the e ient implementation of the

deation method and the further improvements of the method. We also ompare

thedeationmethod withotherwell-known2L-PCGmethodsby onsideringboththe

abstra t variants and their optimal variants with theirtypi al parameters. Numeri al

experiments withbubbly owsare performed toillustratethe theoreti al results.

Remark 1.1. Manytheoreti al results presentedin thisthesis are generallyappli able

andarenot restri tedtoappli ationsofbubblyows,although these bubblyowsare

themain appli ationofthisthesis. Therefore,Assumption 1.2isnotdemanded inthe

generaldis ussion,butisonly requiredwhenthegeneraltheoreti al resultsare applied

tobubblyows. Inaddition,allresultsthatrequireAssumption1.2 analsobeapplied

toother elds wherethisassumption is fullled.

1.5 Outline of the Thesis

Chapter2: Iterative Methods. This hapter is devotedto theintrodu tion of

itera-tive methods, espe ially the CG and PCG methods. In most introdu tory books, CG

andPCGarederivedandanalyzedwhereAisassumedtobeinvertible,butwegivethe

methods and their on ise derivations for general SPSD oe ient matri es.

More-over,some propertiesof these methodsare presented,whi hare not fully learin the

literature. Finally, the drawba ks of PCG are illustrated using numeri al experiments

withbubblyows.

Chapter 3: Deation Method. In order to improve the onvergen e of the

stan-dardPCGmethod, thedeationmethod anditspre onditioned variantare introdu ed

in Chapter 3. We give their derivation in detail and present some new theoreti al

properties. Weshow, boththeoreti allyand numeri ally,that the deation method is

expe tedto bemore ee tive thanthe originalPCG method.

Chapter 4: Sele tion of Deation Ve tors. The su ess of the deation method

highlydependsonthe hoi eofthe so- alled`deationve tors'. Goodapproximations

of eigenve tors asso iated with unfavorable eigenvalues are often hosen, whi h are

usuallydenseand not straightforward to obtain. Inaddition,the densityeldis often

notknownexpli itlyinourbubbly ow appli ation,whi hmight leadtodi ultiesfor

approximating eigenve tors. We analyze this issue in more detail in Chapter 4 and

providesome strategies todetermine the best deation ve tors for bubblyow

prob-lems. We omeupwithseveralsuitable hoi es,whoseutilityisillustratedinnumeri al

experiments.

Chapter 5: Subdomain Deation applied to Singular Matri es. Theoreti al

re-sultsforthe deation method are well-knownif itis appliedtononsingular oe ient

matri es. Theappli ationofthismethod tosingular oe ientmatri esismore

om-pli ated and has not been widely onsidered in the literature. This issue is further

investigated in Chapter 5. We show equivalen es between deation methods applied

to singular and invertible oe ient matri es. We ome up with several

mathemati- allyequivalentvariantsofthetwo-levelpre onditioner orrespondingtothe deation

method. Numeri al experiments are used to show that these variants an be easily

appliedinpra ti e.

Chapter 6: Comparison of Two-level PCG Methods  Part I. The main fo usof

thisthesis is on the deation method, whereas other attra tive 2L-PCGmethods are

known in the literature. In Chapter 6, we ompare the deation method with some

prominent2L-PCGmethods omingfromtheelds ofdeation,DDMandMG.Both

atheoreti alandnumeri al omparisonareperformedusingtheabstra tformsofthese

methods. Weinvestigate theirspe tral properties, equivalen es, ee tiveness and

Chapter 7: Comparison of Two-Level PCG Methods  Part II. In Chapter 6, the

2L-PCG method based on the standard multigridV(1,1)- y le method is ex luded in

the omparison, sin e it has dierent spe tral properties and requiresa spe ial

theo-reti al treatment. Chapter 7 examines this method in more detail. We ompare the

2L-PCGmethods asdis ussed in Chapter 6,and show thatit dependson the hosen

parameters whi h2L-PCGmethod is the mostee tiveone.

Chapter8: E ien yandImplementationIssuesoftheDeationMethod. Inthe

previous hapters, we haveshown that the deation method is expe ted to onverge

fasterthanPCGintermsofiteration ounts. However,thedeationmethod needsto

beimplementede ientlyinorder toobtainafastmethod withrespe tto omputing

timeaswell. ThisissueisexaminedinChapter8,whereweshowhowea hstepofthe

deationalgorithm anbebestimplementedfora lassofproblems. Atea hiteration

ofthe deationmethod, oarselinearsystems should besolved,whi his usuallydone

by a dire t method. If the number of deation ve tors is relatively large, we show

that itis more attra tive to usean iterative method, so thatthe resultingmethod is

based onaninner-outer iterationpro ess. Thisis further detailedand illustratedwith

numeri al experiments in Chapter8.

Chapter 9: Comparison of Deation and Multigrid with Typi al Parameters. In

Chapter6and7,the2L-PCGmethodshavebeen omparedintheirabstra tforms. In

this ase,thedierentparameterswithinthesemethods anbearbitrary,butareequal

for ea h method, whi h allows us to perform a general omparison. The omparison

analsobe arriedoutwithtypi alparameters inthe methods. Ea h2L-PCGmethod

thentakesitsoptimizedsetofparameters thatistypi alintheeldwherethemethod

omes from. Chapter 9 is devoted to this omparison. The aim of this hapter is

to show whi h optimized 2L-PCG method is urrently the best one to apply for 3-D

bubblyow appli ations.

Chapter 10: Bubbly Flow Simulations. In the previous hapters, we have shown

that 2L-PCG methods are bene ial to use for stationary bubbly ow problems. In

Chapter 10, the exa t mathemati al model for the bubbly ows is formulated, so

that real-life time-dependent experiments an be performed. We show that 2L-PCG

methodsredu e signi antlythe omputations ofbubblyow simulationsand are less

sensitivetothe densityeld omparedwith standard PCGmethods.

Chapter 11: Con lusions. The main on lusions of the thesis and ideas for future

resear harepresented in Chapter 11.

Thisthesis isbasedon thete hni al reports[87,132,134,136,137,140,141,144℄,

thepro eedingpapers[138,139,142,147,172℄,and,espe ially,thejournalpapers[85,

1.6 Notation

Throughoutthis thesis,weusethe notation asgiveninTable1.1.

Notation Meaning

I identitymatrix withanappropriatedimension

e ( )



-th olumn ofI withdimension

e ( )

;

matrix with identi al olumnse ( )



1

;

matrix whoseentries are ones

1  olumnof1 ; 0 ;

matrix whoseentries are zeros

0



olumnof0

;

Chapter

2

Iterative Methods

2.1 Introdu tion

Re allthatthe mainfo usofthisthesisis onsolvingthelinearsystem(see Eq.(1.1))

Ax =b; A=[a

ij ℄2R

nn

; (2.1)

whereAisasparseandSPSD oe ientmatrix. Weaimatsolving(2.1)usingKrylov

iterativemethods, whi hare examinedin this hapter.

Westartthis hapterbyreviewingbasi iterativemethods. Thisisfollowedby

pre-sentingtheConjugateGradient(CG)method, whi his awell-knowniterative method

tosolve(2.1). Therateat whi hCGandgeneraliterative methods onvergedepends

greatlyonthespe trumofthe oe ientmatrix,A. Hen e,thesemethodsusually

in-volvease ondmatrixthattransforms Aintoonewithamorefavorablespe trum. The

resultingmethodisthen alledthepre onditionedConjugateGradient(PCG)method,

and is des ribed in Se tion 2.4. Some further onsiderations regarding

pre ondition-ing, starting ve tors and termination riteria of the iterative pro edure are dis ussed

inSe tion2.5. We on lude this hapterwiththe appli ationofthe solverstobubbly

ow problemsinorder toillustratethe performan e ofthe PCG methods.

2.2 Basi Iterative Methods

Iterative methods generate a sequen e of iterates, fx

j

g, that approximate the exa t

solution,x. These methodsessentiallyinvolvematrix Aonly inthe ontext of

matrix-ve tor multipli ations. The starting point of these methods is onsidering a splitting

ofA ofthe form

A=M N ; M ;N 2R nn

; (2.2)

whereMisassumedtobeinvertible. Ifthesplitting(2.2)issubstitutedintoEq.(2.1),

weobtain

From Eq.(2.3),a basi iterative method an be onstru ted asfollows: Mx j+1 =b+Nx j ; (2.4)

where the iterate, x

j+1

, in the (j+1)-th step an be determined from the previous

iterate,x

j

. Eq.(2.4) an be rewrittenas

x j+1 =x j +M 1 r j ; (2.5)

wherethe residualafter thej-th iteration isdenedas

r

j

:=b Ax

j ;

thatisameasureofthe dieren eoftheiterativeandtheexa tsolutionof(2.1). For

the iteration (2.5) to be pra ti al, it must be relatively easy to solve a linearsystem

with M as the oe ient matrix. For example, M = diag(A) is used in Ja obi

iter-ations,M onsists of the lower-triangular part of A in Gauss-Seidel iterations, and a

more ompli atedMisusedin(symmetri )su essiveover-relaxation((S)SOR)

itera-tions,that anbederivedfromGauss-Seideliterationsbyintrodu inganextrapolation

parameter. It has beenre ognized thatthese basi iterative methods are impra ti al,

be ause they onverge slowly,need good tuningofsome parameters, orrequire stri t

onditions for onvergen e. More basi iteration methods and their analysis an be

foundin [8,63,70,120,167℄.

When the rst iterationsof(2.5) are developed, oneobtains































x 0 ; x 1 = x 0 +M 1 r 0 ; x 2 = x 0 +2M 1 r 0 M 1 AM 1 r 0 ; x 3 = x 0 +3M 1 r 0 3M 1 AM 1 r 0 + M 1 A

2 M 1 r 0 ; . . . Thisyields x j+1 2x 0 + span



M 1 r j ;M 1 A(M 1 r j );:::;(M 1 A) j 1 (M 1 r j )

:

Subspa es ofthe form

K j (A;r 0 ):= span



r 0 ;Ar 0 ;A 2 r 0 ;:::;A j 1 r 0

are alled Krylov subspa es with dimension j, belonging to A and r

0

. Hen e, the

followingholds forbasi iterative methods:

x j+1 2x 0 +K j (M 1 A;M 1 r 0 ): (2.6)

Thesemethods are also alled Krylov(-subspa e)methods. FromEq.(2.6), itfollows

basisforK

j

,su hthatthe iterativemethod onvergesfastwithareasonablea ura y

ande ien y withrespe t tomemory storageand omputationaltime.

Krylovmethods anbedividedintostationaryandnonstationaryvariants. Methods

su h as Ja obi, Gauss-Seidel, (S)SOR iterations are stationary methods, sin e the

sameoperationsonthe urrentiterationve torsareperformedinea hiteration. They

are easy to understand and implement, but they are often not ee tive. On the

other hand, nonstationary methods, that have iteration-dependent oe ients, are

a relatively re ent development. Their analysis is ommonly harder to understand,

but they an be highly e ient. They rely on forming an orthogonal basis of the

Krylov sequen e



r 0 ;Ar 0 ;A 2 r 0 ;:::;A j 1 r 0

. The iterates are then onstru ted by

minimizing the residual over the subspa e formed. The prototypi al method in this

lassisthe (pre onditioned)ConjugateGradient((P)CG) method,whi his des ribed

inSe tion 2.3and2.4. Thisis a popularandee tive nonstationaryKrylovsolverfor

linearsystemswithanSPSD oe ientmatrix,asthestorageforonlyalimitednumber

ofve tors is required. For non-SPSD matri es, the Krylov solvers GMRES[121℄and

Bi-CGSTAB [160℄are popularmethods inuse,see[120,160℄.

2.3 Conjugate Gradient Method

TheConjugateGradient(CG)methodisprobablythemostprominentiterativemethod

for solving the SPSD linear system (2.1). It is dis overed independentlyby Hestenes

and Stiefel, and they jointly published the method in [72℄, whi h has be ome the

lassi al referen e on CG. We refer to [8,61,63,86,117,120,161℄ for more details

aboutthismethod.

The purpose of CG is to onstru t a sequen e, fx

j

g, that satises (2.6), with

M=I and thepropertythat

min x j 2K j (A;r0) jjx j xjj A (2.7) holds,where jjwjj A :=

p

(w;Aw)is usedforany w 2R n

. In otherwords,the erroris

minimized inthe A-semi-norm(that is oftenabbreviatedasthe A-norm,ifthere is no

ambiguity) at ea h iteration. Thisminimum is guaranteed to existin general, only if

A is SPSD. Moreover, CG requires thatsear h dire tion ve tors, fp

j

g, are onjugate

withrespe ttoA, i.e.,

(Ap

i ;p

j

)=0; i 6=j; (2.8)

hen ethename`ConjugateGradientmethod'. It anbeshownthat(2.8)isequivalent

tothe fa tthat theresiduals, fr

j

g, formanorthogonal set, i.e.,

(r

i ;r

j

)=0; i 6=j: (2.9)

Now,theCGmethod pro eedsasfollows. The(j+1)-thiterateisupdatedviathe

sear hdire tion:

where j 2R. Thisyields r j+1 =r j  j Ap j : (2.11)

It anbe shownthat

 j = (r j ;r j ) (Ap j ;p j ) (2.12) minimizes jjx j xjj A

over all possible hoi es of 

j

and ensures that Eq. (2.9) is

satised. Thesear hdire tions are updatedusingthe residuals:

p j+1 =r j+1 + j p j ; (2.13) where j 2Requal to  j = (r j+1 ;r j+1 ) (r j ;r j ) (2.14)

ensuresthatEq.(2.8)issatised. Infa t,it anbeshownthatEqs.(2.12)and(2.14)

make p

j+1

onjugate to all previous sear h dire tions, fp

i

: i = 1;:::;jg, and r

j+1

orthogonalto allpreviousresiduals,fr

i

:i=1;:::;jg.

The abovederivation leadstoAlgorithm1,seebelow. Thisis essentiallythe form

ofthe CGalgorithm thatappearedin [72℄.

Algorithm1Conjugate Gradient(CG)solvingAx =b

1: Sele t x 0 . Compute r 0 :=b Ax 0 andset p 0 :=r 0 .

2: for j :=0;1;:::;until onvergen e do

3: w j :=Ap j 4:  j := (r j ;r j ) (w j ;p j ) 5: x j+1 :=x j + j p j 6: r j+1 :=r j  j w j 7:  j := (r j+1 ;r j+1 ) (r j ;r j ) 8: p j+1 :=r j+1 + j p j 9: endfor 10: x it :=x j+1 Remark2.1.

 The iterative solutionof Ax = b is denoted byx

it

in Algorithm 1 todistinguish

itfromthe exa tsolution, x.

 It is straightforward to derive CG from the Lan zos algorithm for solving

sym-metri eigensystems andvi e versa. The relationship an beexploitedto obtain

relevant information about the eigensystem of A. We refer to [8,63,120℄ for

more details.