Efficient and robust Schur complement approximations in the augmented Lagrangian preconditioner for the incompressible laminar flows

Efficient and robust Schur complement approximations in the augmented Lagrangian

preconditioner for the incompressible laminar flows

He, Xin; Vuik, Cornelis

DOI

10.1016/j.jcp.2020.109286

Publication date

2020

Document Version

Final published version

Published in

Journal of Computational Physics

Citation (APA)

He, X., & Vuik, C. (2020). Efficient and robust Schur complement approximations in the augmented

Lagrangian preconditioner for the incompressible laminar flows. Journal of Computational Physics, 408,

[109286]. https://doi.org/10.1016/j.jcp.2020.109286

Important note

To cite this publication, please use the final published version (if applicable).

Please check the document version above.

Copyright

Other than for strictly personal use, it is not permitted to download, forward or distribute the text or part of it, without the consent of the author(s) and/or copyright holder(s), unless the work is under an open content license such as Creative Commons. Takedown policy

Please contact us and provide details if you believe this document breaches copyrights. We will remove access to the work immediately and investigate your claim.

This work is downloaded from Delft University of Technology.

Green Open Access added to TU Delft Institutional Repository

‘You share, we take care!’ – Taverne project

https://www.openaccess.nl/en/you-share-we-take-care

Otherwise as indicated in the copyright section: the publisher

is the copyright holder of this work and the author uses the

Dutch legislation to make this work public.

Contents lists available atScienceDirect

Journal

of

Computational

Physics

www.elsevier.com/locate/jcp

Efficient

and

robust

Schur

complement

approximations

in

the

augmented

Lagrangian

preconditioner

for

the

incompressible

laminar

flows

Xin He

a

,

Cornelis Vuik

b

,

∗

a_{StateKeyLaboratoryofComputerArchitecture,InstituteofComputingTechnology,ChineseAcademyofSciences,No.6KexueyuanSouth}

RoadZhongguancunHaidianDistrictBeijing,100190,P.R.China

b_{DelftInstituteofAppliedMathematics,DelftUniversityofTechnology,VanMourikBroekmanweg6,2628XE,Delft,theNetherlands}

a

r

t

i

c

l

e

i

n

f

o

a

b

s

t

r

a

c

t

Articlehistory:

Received12November2018

Receivedinrevisedform17January2020 Accepted25January2020

Availableonline7February2020 Keywords:

Navier-Stokesequations Stabilizedfiniteelementmethod Blockstructuredpreconditioners Schurcomplementapproximations AugmentedLagrangianpreconditioner

Thispaper introducesthree new Schur complementapproximations for theaugmented Lagrangian preconditioner. The incompressible Navier-Stokes equations discretized by a stabilized finite element method are utilized to evaluate thesenew approximations of the Schur complement. A widerangeof numerical experiments inthe laminar context determines the most efficient Schur complement approximation and investigates the effectof the Reynoldsnumber, meshanisotropy and refinementonthe optimal choice. Furthermore, the advantage over the traditional Schur complement approximation is exhibited.

©2020ElsevierInc.Allrightsreserved.

1. Introduction

Inthispaperweconsiderthenumericalsolutionofthesteady,laminarandincompressibleNavier-Stokes equationsas follows

−

ν



u

+ (

u

· ∇)

u

+ ∇

p

=

f on

,

∇ ·

u

=

0 on

.

(1)

Here u is thevelocity, p isthe pressure,the positive coefficient

ν

isthe kinematicviscosityand f isa givenforce field.



isa2Dor3Dboundedandconnecteddomainwiththeboundary

∂

.Ontheboundariesofthecomputationaldomain, eithertheDirichletboundarycondition u

=

g or Neumannboundarycondition

ν

∂u

∂n

−

np

=

0 isimposed,wheren denotes

theoutward-pointingunitnormaltotheboundary.

AfterthePicardlinearizationandFEMdiscretization[1], theincompressibleNavier-Stokesequationsconverttothe fol-lowinglinearsysteminsaddle-pointform



A BT B

−

C

 

u p



=



f g



with

A :=



A BT B

−

C



,

(2)

*

Correspondingauthor.

E-mailaddresses:hexin2016@ict.ac.cn(X. He),c.vuik@tudelft.nl(C. Vuik). URLs:http://english.ict.cas.cn/(X. He),http://www.ewi.tudelft.nl/(C. Vuik).

https://doi.org/10.1016/j.jcp.2020.109286

wherethematricesB andBT correspondtothedivergenceandgradientoperators,respectively.Picardlinearizationleadsto thematrix A inblock diagonalstructure,andeachdiagonalblockcorrespondstotheconvection-diffusionoperator.Dueto thepresenceoftheconvectiveterm, A isnotsymmetric.ForthefiniteelementdiscretizationsatisfyingtheLBB(‘inf-sup’) stability condition [1], no pressure stabilization is required and C

=

0 is taken. When LBB unstable finite elements are applied,thenonzeromatrixC correspondstoastabilizationoperator.

Blockstructuredpreconditioners[1–3] areoftenutilized toacceleratetheconvergenceoftheKrylovsubspacesolversfor saddlepointsystemsas(2).Theyarebasedontheblock

LDU

decompositionofthecoefficientmatrixgivenby

A = LDU =



A BT B

−

C



=



I1 O B A−1 I2

 

A O O S

 

I1 A−1BT O I2



,

(3)

where S

= −(

C

+

B A−1BT

)

is theso-calledSchur complement. A combinationofthisblock factorizationwitha suitable approximation oftheSchur complementisutilized tosuccessfullydesignthe blockstructured preconditioners,whichare givenasfollows

P

F

=



A O B



S

 

I1



A−1BT O I2



,

(4)

P

L

=



A O B



S



,

P

U

=



A BT O



S



.

(5)

Multiplying the

LD

and

DU

factors of(3) results inthe block lower- andupper-triangularpreconditioners

P

L and

P

U,

respectively. Preconditioner

P

F isbasedonthe multiplicationofthe

LDU

factors.Theterm



A−1 denotessome

approxi-mation oftheinverse actionof A,which isgiveneitherinan explicitformorimplicitlydefinedviaan iterativesolution methodwithaproperstoppingtolerance.

It is notpractical to explicitlyformthe exactSchur complementdueto theaction of A−1, typicallywhen thesize is large.Thisimpliesthatthemostchallenging taskistofindthespectrallyequivalentandnumericallycheapapproximation oftheSchur complement,whichisdenotedby



S in(4) and(5).Thereexistseveralstate-of-the-artapproximationsofthe Schurcomplement,e.g.theleast-squarecommutator(LSC)[4,5],pressureconvection-diffusion(PCD)operator[6,7] andthe approximations fromthe SIMPLE(R)[8–10] andaugmented Lagrangian (AL) preconditioner[11,12] etc. Among them, the ALpreconditionerexhibitsattractivefeatureswithstablefiniteelementmethods(FEM)usedforthediscretization,e.g.the purelyalgebraicandsimpleconstructionoftheSchurcomplementapproximationandrobustnesswithrespecttothemesh refinementandReynoldsnumber,atleastforacademicbenchmarks.Motivatedbytheseadvantages,thefurtherextensionto thecontextoffinitevolumemethod(FVM)[13] andthemodifiedvariant[14] withreducedcomputationalcomplexitiesare promoted.Recently,theauthorsofthispaperproposeanewvariantoftheALpreconditioner[15] fortheReynolds-Averaged Navier-Stokes(RANS)equationsdiscretizedbyastabilizedFVM,whicharewidelyusedtomodelturbulentflowsinindustrial computationalfluiddynamic(CFD)applications.

The roleoftheALtermforpreconditioningisverysimple:by varyingparameter

γ

itputs moreweightoneitherthe (1,1)orthe(2,2)blockoftheALpreconditioner.Ifonecannotaffordlargervaluesof

γ

,thenfindingasuitable(more compli-cated)preconditionerforSγ becomesimportantagain,whereSγ denotestheSchurcomplementfortheaugmentedsystem. Morediscussions on Sγ andtheinvolvedparameter

γ

aregiveninSection 3ofthispaper.Knownrepresentationsfor Sγ

[11,14] suggestwaystoutilizeearlierdevelopedpreconditionersforthenon-augmentedproblem.Thepaperisbuiltonthis simpleobservationandtheoriginalideaasgivenin[11].Thisobservationisalreadyexploited,forexample,in[14,16].The challengesencounteredintheturbulentcalculations[17–19] areinevitablefactorswhichcouldcausethebreakdownofthe AL preconditioner,includingthehighReynoldsnumber,high-aspectratiocellsnearthevery thinboundary layerandthe significant variation inthevalue ofviscositydueto thepresenceof theeddy-viscosity.Toovercome thesechallenges, an alternative methodtoapproximatethe SchurcomplementfortheALpreconditionerisintroduced in[15], whichleads to a new variantof theAL preconditioner. Thisnew methodapproximates the Schur complementthrough its inverseform andfacilitatestheutilizationoftheexistingSchurcomplementapproximations.Amongtheavailablecandidates,theSchur complementapproximationfromtheSIMPLEpreconditioner[8,20] ischosenandsubstitutedintotheinverseSchur comple-mentapproximationfortheALpreconditioner.ThischoiceismotivatedfromthenotionthatitreducestoascaledLaplacian matrix[8,20] withtheconsideredFVMandits promisingefficiencyontheturbulentapplicationsofthemaritimeindustry [8,21].Consequently,theso-arisingnewvariantoftheALpreconditionerreducesthenumberofKrylovsubspaceiterations byafactorupto36 comparedtotheoriginalone[15].

Since the new methodto approximate the Schur complementfor theAL preconditioner usethe existing Schur com-plement approximations, the following questions straightforwardly raise. Does the utilization of other existing Schur complement approximationsdeliver a better performance than that fromthe SIMPLEpreconditioner? Ifso, whichSchur complementapproximationisthemostefficientone?Doestheoptimalchoicedependonthetestproblemandparameters arising from the physics and discretization,e.g. the Reynolds number andgrid size? To answer these questions,in this paper we utilizethe existing Schur complementapproximations not onlyfromthe SIMPLE preconditionerbutalso from the LSCand PCD operators to construct the new Schur complement approximation in the AL preconditioner. Moreover, extensivecomparisonsbetweentheconsidered Schurcomplementapproximationsare carriedoutonawide rangeof nu-mericalexperimentstoevaluatetheeffectoftheReynoldsnumber,meshanisotropyandrefinementontheoptimalchoice.

Thesenumericalevaluationsareconsideredinthecontextoflaminarflows,whichismotivatedbytheexpectationthatthe obtainedresultscanprovideafundamentalguidelineforthemorecomplicatedturbulentflowcalculations.

The structure of this paper is given as follows. The utilized stabilization method and a brief survey of the existing approximationsoftheSchurcomplementareintroducedinSection2.Section3illustratesthemethodusingtheseexisting SchurcomplementapproximationstoconstructthenewapproximationoftheSchurcomplementintheALpreconditioner. Section4includesnumericalresultsonvaryinglaminarbenchmarks.ConclusionsandfutureworkareoutlinedinSection5.

2. StabilizationmethodandsurveyofSchurcomplementapproximations 2.1. Stabilization

Inthispaperwe usethemixedFEMwhichdoesnotuniformlysatisfya discreteinf-supcondition [1] todiscretizethe Navier-Stokesequationsgoverninglaminarflows,whichischosenbythefollowingconsiderations.Firstly,theexistingSchur complement approximationsare originally designedwith finite element methods used for discretization. Therefore,it is expectedtoapplythenewSchurcomplementapproximationfortheALpreconditionerintheFEMcontext.Inaddition,this closesagapintheapplicationofthenewSchurcomplementapproximation.Secondly,boththestabilizedFEM[1] andFVM [17] lead tosaddlepoint systemwithanonzero

(

2

,

2

)

block whicharisesfromthe pressurestabilization. Thanks to this similarity,aminoradaptionisrequiredtoextendthenewvariantoftheALpreconditionerfromthestabilizedFVMtothe stabilizedFEM.Finally,theutilizationofstabilizedFVMdegradesthegeneralitytosomeextentsincetheSchurcomplement approximationintheSIMPLEpreconditionerreducestoaspecialformation[8,20].However,thisspecialformationcannot beobtainedwithotherstabilizationanddiscretizationmethods.UsingstabilizedFEM,allSchurcomplementapproximations consideredinthispaperareexpressedintheirdefinedmanners,includingthatfromtheSIMPLEpreconditioner.Inthisway, aconvincingevaluationofthenovelSchurcomplementapproximationfortheALpreconditionercanbeexpected.

Basedontheabovemotivations,inthispaperweusethe Q1- Q1 mixedfiniteelementapproximationwheretheequal

first-orderdiscretevelocitiesandpressurearespecifiedonacommonsetofnodes.Amongtheavailablestabilization meth-ods[22–27] specifiedforthe Q1- Q1discretization,wechoosetheapproachintroducedin[25].Themainmotivationisthat

therearefewstabilizationparametersrequiredinthefollowingoperator

C(proj)

(

ph

,

qh

)

=

1

ν

(

ph

− 

0ph

,

qh

− 

0qh

),

(6)

where



0 is the L2 projectionfrom thepressure approximationspaceinto thespace P0 ofthepiecewise constant basis

function.Thisprojectionisdefinedlocally:



0phisaconstantfunctionineachelement

2

k

∈

Th.Itisdeterminedsimplyby

thefollowinglocalaveraging



0ph

|

2k

=

1

|2

k

|



2k ph

,

for all

2

k

∈

Th

,

(7)

where

|2

k

|

is thearea ofelement k. Due tothe locality asillustrated by equation (7), thestabilizationmatrix C can be

assembledfromthecontributionmatricesonmacroelementsinthesamewayasassemblingastandardfiniteelementmass matrix.Takingthe2Drectangulargridasanexample,the4

×

4 macroelementcontributionmatrixC(macro)_isgivenby

C(macro)

=

1

ν

(

M

(macro)

₋

_qqT

_|2

k

|),

(8)

where M(macro) isthe4

×

4 macroelementmassmatrixforthe bilineardiscretizationandq

= [

1

/

4

,

1

/

4

,

1

/

4

,

1

/

4

]

T isthe localaveragingoperator.The nullspaceofthemacroelementmatrixC(macro) _{andassembled stabilizationmatrixC consist}

ofconstantvector,see[1,25] formoredetails.

Contrary to other pressure stabilizationmethods [27,28] which utilize the viscosity and velocity fields to derive the scalingparameterinfrontofthestabilizationmatrix,thealternativeemployedinthispaperonlyinvolvestheviscosity co-efficient.Resultsinnumericalexperimentsectiondemonstratethattheutilizedstabilizationmethodresultsinareasonable andsmoothcalculationofthevelocityandpressureunknownsrangingfromthemoderatetolargeReynoldsnumbers.The assessment ofother pressure stabilizationmethods andtheir effects on theproposed preconditioning techniques bythis paperisincludedinfutureresearchplan.

2.2. SurveyofSchurcomplementapproximations

As follows we briefly introduce severalstate-of-the-art Schur complementapproximations which are utilized to con-structthenewapproximation ofthe SchurcomplementfortheALpreconditioner.Werefer formore detailsoftheSchur complementapproximationtothesurveys[2,3,29,30] andthebooks[1,31].

In the following illustration, we use the notation p to indicate the operators definedon the pressure space andthe notationu fortheoperatorsdefinedonthevelocityspace.

(1) Thepressureconvection-diffusionoperator



SP C D.

Thisapproximation,denotedby



SP C D,isproposedbyKayetal. [6] anddefinedas



SP C D

= −

LpA−p1Mp

,

(9)

where Mp isthe pressure massmatrix, and Ap and Lp are the discrete pressureconvection-diffusion andLaplacian

operators, respectively.AlthoughthePCDSchur complementapproximation(9) isoriginally proposedforstablefinite elementmethods,itisstraightforwardlyapplicableforthediscretizationsneedingastabilizationterm,e.g.the Q1- Q1

pair.Formoredetailsaboutthisextensionwereferto[1].Ontheotherhand,thisapproximationrequiresusersto pro-vide thediscreteoperators Ap andLp andpresetsomeartificialpressureboundaryconditionsonthem.Theboundary

conditionscouldstronglyeffecttheperformancesoappropriateonesshouldbecarefullyselectedbasedontheproblem characteristic[32,33].ApplyingthePCDSchurcomplementapproximationinvolvestheactionofaPoissonsolve,amass matrixsolveandamatrix-vectorproductwiththematrixAp.

(2) Theleast-squarecommutator



SL SC.

Elmanetal. [4] originallypropose thismethodforstablefiniteelementdiscretizations andthen extenditto alterna-tives[5] thatrequirestabilization. Forsystem(2) witha nonzerostabilizationoperatorC , theLSCSchurcomplement approximation



SL SC isdefinedas



SL SC

= −(

BM



u−1BT

+

C1

)(

BM



−u1AM



u−1BT

+

C2

)

−1

(

BM



−u1BT

+

C1

),

(10)

where M



u denotesthe diagonalapproximation ofthe velocity massmatrix Mu, i.e.M



u

=

diag

(

Mu

)

.Giventhe

stabi-lizationmatrixC assembledfromthemacroelementcontributionmatrixC(macro) _{(8),thecontributionmatricesC}(macro)

1

andC₂(macro) fortheassociatedstabilizationmatricesC1 andC2areintroducedby

C(₁macro)

=

ν

|2

k

|

·

C(macro)

,

C(₂macro)

=

ν

2

|2

k

|

2

·

C(macro)

,

(11)

where

ν

denotes theviscosity parameter. Forthe derivationof C₁(macro) and C₂(macro) we refer to [5]. The implemen-tationoftheLSCSchur complementapproximationdoesnotrequireanyartificialboundarycondition andconsistsof one matrix-vector product withthe middle term in (10) and two solves with the other term. When the LSC Schur complementapproximationisappliedtostablefiniteelementdiscretizations,thematricesC1andC2aresettozeroin

(10).

(3) Theapproximation



SS I M P L EfromtheSIMPLEpreconditioner.

SIMPLE(Semi-Implicit PressureLinkedEquation)isused byPatanker [18] asan iterativemethodtosolve the Navier-Stokes problem.Theschemebelongstotheclassofbasiciterativemethodsandexhibitsslowconvergence.Vuiketal. [9,10] useSIMPLEasapreconditionerinaKrylovsubspacemethod,achievinginthisway,amuchfasterconvergence. RegardingtheSchurcomplementS

= −(

C

+

B A−1_BT

₎

_{ofsystem(2),theSIMPLEpreconditionerapproximates A byits}

diagonal,i.e.diag

(

A

)

,andobtainstheapproximation



SS I M P L E as



SS I M P L E

= −(

C

+

Bdiag

(

A

)

−1BT

).

(12)

Substituting



SS I M P L E and



A−1

=

diag

(

A

)

−1 into (4) leads tothe so-calledSIMPLEpreconditioner. Forstablefinite

el-ement discretizations, C

=

0 is set in system(2) and correspondingly in the Schur complementapproximation (12). The easyimplementationandpromisingperformanceonthe complicatedmaritimeproblems[8,21] maketheSIMPLE preconditioneranditsvariantsattractiveinrealworldapplications.

Themaingoalofthispaperistoutilize theabovementionedSchurcomplementapproximationstoconstructanew ap-proximationoftheSchurcomplementintheALpreconditioner,withmoredetailspresentedinthenextsection.Theoretical analysis andnumerical evaluationoftheabove Schur complementapproximationsfallout ofthe scopeofthiswork and we referto[1,3,34] formoreresults.Herewesummarizethekeydifferences.



SP C D requirestheconstructionofadditional

matrices onthe pressurespacewhile



SL SC and



SS I M P L E rely onmatriceswhich could be easily generatedorare readily

available. Asseenfrom



SL SC,thestabilizationterms C1(macro) andC

(macro)

2 areeasily obtainedbysubstitutingtheavailable

termC(macro)_{into(11).Ontheotherhand,}



_S

P C D easilyextendstothestabilizedelementsandaminoradaptionisrequired

by



SS I M P L E forthisextension.However,



SL SC doesnotimmediatelyapplyandneedsappropriatestabilizationtermsC1and C2.Wefurthernote that boundaryconditionsforthepressure unknowns,which havefewphysicalmeanings, havetobe

considered forLp and Ap in



SP C D.What boundaryconditionsworkbestwithaspecifictype ofproblemisusually based

onexperimentalknowledge[32,33].

3. AugmentedLagrangianpreconditioner

The focusofthissectionis thenewmethodtoapproximatetheSchur complementintheaugmentedLagrangian (AL) preconditioner.Inthefollowing,wefirstbrieflyrecalltheALpreconditionerandthenintroducethenewmethodfollowed byacomparisonwiththeoldone.

ThemotivationofapplyingtheALpreconditioneristocircumventthechallengeonfindingtheefficientapproximation oftheSchurcomplementS for theoriginal system(2), cf., [11,14].Toapply theALpreconditioner,theoriginalsystem(2) istransformedintoanequivalentonewiththesamesolution[13,14],whichisoftheform



Aγ BγT B

−

C

 

u p



=



fγ g



with

A

γ

:=



Aγ BTγ B

−

C



,

(13)

where Aγ

=

A

+

γ

BTW−1B,BT_γ

=

BT

−

γ

BTW−1C andfγ

=

f

+

γ

BTW−1g.Thistransformationisobtainedbymultiplying

γ

BTW−1 onboth sidesof thesecond row ofsystem(2) andaddingthe resulting equation tothe first one.Clearly, the transformedsystem(13) hasthesamesolutionassystem(2) foranyvalueof

γ

andanynon-singularmatrixW .TheSchur complementofthetransformedsystem(13) isSγ

= −(

C

+

B A−1

γ BγT

)

.

TheALpreconditionerisappliedfortheequivalentsystem(13),whichistobesolved.Usingtheblock

D

U decomposition

of

A

γ ,theidealALpreconditioner

P

I AL anditsvariant,i.e.themodifiedALpreconditioner

P

M AL,aregivenby

P

I AL

=



Aγ BTγ O



Sγ



and

P

M AL

=



Aγ BγT O



Sγ



,

(14)

where



Sγ and



Aγ denotetheapproximationsofSγ and Aγ ,respectively.

First we consider the approximation



Aγ . Given the original pivot matrix A

=



A1 O O A1



and the divergence matrix

B

=



B1 B2



inthe2Dcase,thetransformedpivotmatrixAγ

=

A

+

γ

BTW−1B canbewrittenas Aγ

=



A1

+

γ

BT₁W−1B1

γ

B₁TW−1B2

γ

BT₂W−1B1 A1

+

γ

B2TW−1B2



.

Contraryto

P

I AL,

P

M AL approximates Aγ byitsblockupper-triangularpart,i.e.



Aγ withazero(2,1)block,suchthatthe

difficulty ofsolvingthe systemswith Aγ is avoided[14].When applying

P

M AL oneneeds tosolve thesub-systems with

thediagonal blocksof Aγ ,i.e. A1

+

γ

B1TW−1B1 and A1

+

γ

BT2W−1B2, whichdonot containthe couplingbetweentwo

componentsofthevelocitysothat standardalgebraicmultigridmethodscanbeapplied[34].Thisadvantagemotivatesus tochoose

P

M AL inthispaperdespitetheobservationthattheperformanceof

P

M ALisdependentofthevaluesof

γ

,which

isseeninthenumericalexperimentsofthispaperandotherrelatedreferences[14].Theaboveadvantagealsomotivesto approximateAγ byitsblocklower-triangularpartwithazero(1,2)block.Numericalexperimentsdemonstratethatdifferent approximationsofAγ slightlyeffecttheperformanceofthemodifiedALpreconditionerfortheconsideredbenchmarks.For thisreason,inthispaperweonlyillustratetheresultsby applyingtheblock upper-triangularapproximationof Aγ inthe modifiedALpreconditioner.RegardingtheidealALpreconditioner

P

I AL,standardmultigridmethodsareineffectivetosolve

thesystemswith Aγ .Aspecializedmultigridalgorithmfor Aγ isbuiltin[11] andtheextension tothethreedimensional applicationsis recentlyproposed in[35]. Alternatively,previous work [12] suggests tosolve thesystems with Aγ bythe Krylovsubspacemethods,whichareacceleratedbytheapproximateinversepreconditionerbasedontheShermon-Morrison formula.Intherelatedwork[34],thecomparisonbetweenthemodifiedandidealALpreconditionersisrealizedbyapplying the direct solution method for the involved sub-systems. Although fewer Krylov iterations are needed by the ideal AL preconditioner,removing thedifficultytosolve thesub-systemswith Aγ makesthemodifiedALpreconditionerattractive inpractice.

3.1. NewSchurapproximationintheALpreconditioner

FindinganeffectiveapproximationoftheSchurcomplementSγ isthekeyfortheidealandmodifiedALpreconditioners. Thispaper ismeant tousethe available Schurapproximationsfor theoriginal system(2), asintroduced inSection 2,to constructanewapproximationofSγ .ThenewSchurcomplementapproximationisrealizedbyusingthefollowinglemma.

Lemma3.1.Assumingthatalltherelevantmatricesareinvertible,thentheinverseofSγ isgivenby

S−_γ1

=

S−1

(

I

−

γ

C W−1

)

−

γ

W−1

,

(15)

whereS

= −(

C

+

B A−1_BT

₎

_{denotestheSchurcomplementoftheoriginalsystem(2).}

Proof. Fortheproofwereferto[13,14].

2

Lemma3.1 is originally revealed by [14] and used to derive the old approximation of Sγ , which is discussed inthe next section. Here, Lemma 3.1is viewed from another side. Since Lemma 3.1 buildsthe connection between the Schur complements Sγ and S, the natural andsimple method to approximate Sγ is substituting the approximation of S into



S_γ−1_new

=

S−1

(

I

−

γ

C W−1

)

−

γ

W−1

,

(16)

where



S denotestheapproximationofS.

ThenovelapproachprovidesaframeworktousetheknownSchurcomplementapproximation



S fortheoriginalsystem (2) toconstruct



Sγ new inthe ALpreconditioner,whichisapplied tothetransformedsystem(13). SubstitutingtheSchur

complementapproximationsdemonstratedinSection2,i.e.



SP C D,



SL SC and



SS I M P L E intoexpression(16),threevariantsof



_Sγ newarederivedas

• 

S−_γ1_{P C D}

=

S−_{P C D}1

(

I

−

γ

C W−1

)

−

γ

W−1,

• 

S−_γ1_{L SC}

=

S−_{L SC}1

(

I

−

γ

C W−1

₎

₋

_γ

_W−1_,

• 

S−_γ1_{S I M P L E}

=

S−_{S I M P L E}1

(

I

−

γ

C W−1

)

−

γ

W−1.

Followingotherrelatedreferences[11,14],inthispaperwechoosethematrixparameterW tothediagonal approxima-tion ofthepressuremassmatrix,i.e.W

= 

Mp

=

diag

(

Mp

)

.It istrivialto obtaintheactionof W−1 inthetransformation

(13) andthenewSchurcomplementapproximation(16).ApplyingthenewSchurcomplementapproximation



Sγ new

con-verts tosolveasystemwithitandthechoiceof W

= 

Mp focusesthecomplexitymainlyonthesolveof



S.Thisimpliesa

limitedincrease ofthecomplexity whenimplementingthenewSchurcomplementapproximation



Sγ new comparedto



S.

Inaddition,theconsiderableeffortstooptimizetheapproximation



S canstraightforwardlyreducethecomputationaltime of



Sγ new.

WhenapplyingstabilizedFVM,theinverseofSγ isexpressedinasimilarmanner[15] asLemma3.1andthissimilarity facilitates theextension ofthenewSchur complementapproximationfromthestabilized FVMtothe stabilizedFEM. Re-gardingthenewSchurcomplementapproximation,therearetwomaindifferencesbetween[15] andthiswork.Firstly,only



_Sγ S I M P L E isconsideredin[15] and inthispaperweintroduce threevariants,i.e.



Sγ P C D,



Sγ L SC and



Sγ S I M P L E.Inthis

way,thecomparisonbetweenthemisexpectedtoanswerthequestionsraisedintheintroductionsectionandfindoutthe optimalchoice.Secondly,in[15] finitevolumediscretizationstabilizedbythepressure-weightedinterpolationmethod[36] is applied,which leads to



SS I M P L E in areducedform. The generality isdegraded since thisspecialform of



SS I M P L E can

notbeobtainedbyusingotherstabilizationanddiscretizationmethodsingeneral.Inthispaper,theapproximations



SP C D,



SL SC and



SS I M P L E are expressedintheir definedmannerssothata convincingassessmentofthenewSchurcomplement

approximationcanbeexpected.

Based onthe above approach,it is easy to see that there is no extra requirement on thevalue ofthe parameter

γ

. This advantage of the new Schur complement approximation can be more clearly seen in the next section, where the contradictoryrequirementsonthevaluesof

γ

intheoldapproacharepresented.

3.2. OriginalSchurapproximationintheALpreconditioner

The starting point to constructthe original approximation ofthe Schur complement inthe AL preconditioner is also Lemma3.1.However,thestrategyistotallydifferent.ChoosingW1

=

γ

C

+

Mp andsubstitutingW1intoexpression(15) we

have

S−_γ1

=

S−1

(

I

− (γ

C

+

Mp

−

Mp

)(

γ

C

+

Mp

)

−1

)

−

γ

(

γ

C

+

Mp

)

−1

=

S−1Mp

(

γ

C

+

Mp

)

−1

−

γ

(

γ

C

+

Mp

)

−1

= (γ

−1_S−1_M

p

−

I

)(

C

+

γ

−1Mp

)

−1

.

Forlargevaluesof

γ

suchthat



γ

−1_S−1_M

p



1,theterm

γ

−1S−1Mpcanbeneglectedsothatwehave



Sγ origasfollows



Sγorig

= −(

C

+

γ

−1Mp

).

(17)

As shownabove,thechoice ofW1

=

γ

C

+

Mp isusedtoderivetheoriginal Schurcomplementapproximation



Sγ orig.

However, thechoiceof W1

=

γ

C

+

Mp isnotpracticalsincetheactionofW1−1 isneededinthetransformedsystem(13).

Onepracticalchoiceistoomittheterm

γ

C inW1andreplaceMp byitsdiagonalapproximation,whichleadstoW

= 

Mp.

This modificationis onlyapplied tosimplifythe matrixparameter W and theoriginal Schur complementapproximation



_Sγ origremains thesameasgivenin(17).Insummary,thechoiceof W

= 

Mp and



Sγ origisusedinthispaperandother

relatedwork,forinstance[13,14] wherestabilizeddiscretizationsareemployed.

Thecontradictoryrequirementsintheaboveapproximationareshownasfollows.Theapproximation



Sγ origisobtained

ifandonlyifW1

=

γ



C

+

Mpandlargevaluesof

γ

arechosen.However, W

= 

MpisspectrallyequivalenttoW1

=

γ

C

+

Mp

onlywhen

γ

issmall.Thismeansthatitiscontradictorytotunethevalueof

γ

sothatW

= 

Mp and



Sγ orig couldbe

si-multaneously obtained. By contrast,thiscontradictory requirementsare avoided by applyingthe newSchur complement approximation asgiven inSection 3.1.This disadvantage ofthe original Schurcomplement approximation reflectsin the

Table 1

Summaryofthelinearsystemstobesolved,appliedpreconditionersandapproximationsof theSchurcomplementutilizedtherein.

Linear system Preconditioner Schur complement approximations Transformed system withAγ PM AL Sγ P C D,SγL SC,Sγ S I M P L E,Sγorig Original system withA PU SP C D,SL SC,SS I M P L E

Table 2

Pressuresub-system‘mass-p’withSγ inPM ALandS inPU,andthesystemsinvolvedtherein.

‘mass-p’withSγ new ‘mass-p’withS Systems involvedinS

Sγ P C D SP C D LpandMp

Sγ L SC SL SC (BM−u1BT+C1)twice

Sγ S I M P L E SS I M P L E C+Bdiag(A)−1BT

‘mass-p’withSγ orig – Systems involvedinSγorig

Sγ orig – C+γ−1Mp

slowerconvergencerateoftheKrylovsubspacesolvers comparedtothenewSchurcomplementapproximation.This con-clusionismadebasedonthefactthattheperformanceofthemodifiedALpreconditionerisevaluatedbyvaryingtheSchur complementapproximations.Seemoreresultsinthenumericalsection.

The applicationoftheoriginal Schurcomplementapproximation



Sγ orig involvesthesolution ofthesystemwith C

+

γ

−1_M

p.SincethecontributionstabilizationmatrixC(macro) onmacroelementsconsistsofthemacroelementpressuremass

matrix asillustrated in (8), the presence ofthe assembled pressure mass matrix Mp does not introduce morenon-zero

fill-ininthestabilizationmatrixC .

3.3. SummaryoftheSchurcomplementapproximations

AteachPicarditeration,wesolveeitherthetransformedsystem(13) withthecoefficientmatrix

A

γ ortheoriginal

sys-tem(2) withthecoefficientmatrix

A

.WeapplythemodifiedALpreconditioner

P

M AL (14) andtheblockupper-triangular

preconditioner

P

U (5) tothetransformedandoriginalsystems,respectively.TheSchurcomplementapproximationsapplied

in

P

M AL and

P

U aresummarizedinTable1.

Duetothesmallsizeoftestproblemsandthelackofcodeoptimization,thecomplexitycomparisonofpreconditioners

P

M AL and

P

U is done based on the followingcosts analysisin thispaper, instead of reportingthe computational time.

Firstly,we considerthecosts ofusingthemodified ALpreconditioner

P

M AL foraKrylovsubspacemethodthatsolvesthe

systemwith

A

γ . Thepreconditioner isapplied ateach Kryloviteration andthe modified ALpreconditioner involvesthe

solutionofthe momentumsub-system ‘mom-u’with



Aγ and thepressuresub-system ‘mass-p’with



Sγ .Furthermore,at eachKryloviterationadditionalcostsareexpressedintheproductofthecoefficientmatrix

A

γ withaKrylovresidualvector bres.Thus,thetotalcostsateachKryloviterationare

• P

M AL:mom-uwith



Aγ +mass-pwith



Sγ +

A

γ

×

bres.

Clearly, the difference ofcosts by applying

P

M AL arisesfromsolving thepressure sub-system ‘mass-p’with different

Schur complement approximations. Ifwe ignore the multiplications in the definitionof the newSchur complement ap-proximation



Sγ new, finding the solution of the pressure sub-system in

P

M AL withthree variants derived from



Sγ new,

i.e.,



Sγ P C D,



Sγ L SC and



Sγ S I M P L E isreducedtosolve thepressuresub-system in

P

U with



SP C D,



SL SC and



SS I M P L E,

re-spectively. Systems involved in



SP C D,



SL SC and



SS I M P L E are shownin Table 2. The costs of applying the original Schur

complementapproximation



Sγ origarealsoincludedinTable2foracomparisonwiththenewSchurcomplement

approxi-mation



Sγ new.Notethatallinvolvedsystemsareofthesamesize.Ifweassumeacomparablecomplexitytosolvedifferent

involvedsystems,theanalysisinTable2showsthatthecostsofusing

P

M ALwith



Sγ P C D and



Sγ L SC areroughlythesame

andtwotimesofthatwith



Sγ S I M P L E and



Sγ orig.

Secondly,we considerthecosts ofapplyingtheupperblock-triangularpreconditioner

P

U withdifferentSchur

comple-mentapproximations,whichareusedfortheoriginalsystem.Similartotheanalysisof

P

M AL,weobtainthetotalcosts at

everyKryloviterationas

• P

U:mom-uwith A +mass-pwith



S +

A

×

bres.

Also,varyingSchurcomplementapproximations



S resultsinthedifferenceofcostsby applying

P

U.Basedontheanalysis

inTable2andtheassumptionofacomparablesolutioncomplexity forall involvedsystems,wefindout thatthecostsof applying

P

U with



SP C D and



SL SC areroughlythesameandtwotimesofthatwith



SS I M P L E.

Lastly,wecomparethecostsbetween

P

M AL and

P

U.Asmentionedbefore,solvingthepressuresub-systemwiththenew

Schurcomplementapproximation



Sγ newin

P

M AL canbereducedtocalculatethesolutionofthepressuresub-systemwith



S,whichistheSchurcomplementapproximationusedin

P

U.Thus,thedifferenceofcosts between

P

M ALand

P

U focuses

on the solutionof themomentum sub-system andtheproduct ofthe coefficient matrixwith theKrylov residualvector. Morenon-zerofill-inin Aγ and

A

γ [13],comparedto A and

A

,resultsinaheaviermatrix-vectorproductwhenapplying

P

M AL ateachKryloviteration.However,theheaviercomplexityof

P

M AL couldbepaidoffbyareducednumberofKrylov

iterations. In this paper we obtain a faster convergence ratepreconditioned by

P

M AL with the new Schur complement

approximations, comparedto

P

U used forthe original system. The time advantage of

P

M AL needs a further assessment

whichisincludedinfutureresearchplan.

4. Numericalexperiments

Inthissection,wecarryoutnumericalexperimentsonthefollowing2Dlaminarbenchmarks:

(1) Flowoverafiniteflatplate(FP)

This example, known asBlasius flow, models a boundary layer flow over a flat plateon the domain



= (−

1

,

5

)

×

(

−

1

,

1

)

.Tomodelthisflow,theDirichletboundaryconditionux

=

1

,

uy

=

0 isimposedattheinflowboundary(x

= −

1;

−

1

≤

y

≤

1)andalsoonthetopandbottomofthechannel(

−

1

≤

x

≤

5; y

= ±

1),representingwallsmovingfromleft torightwithspeed unity.Theplateismodeledby imposingano-flowconditionontheinternalboundary(0

≤

x

≤

5;

y

=

0), andthe Neumanncondition is applied attheoutflowboundary (x

=

5;

−

1

<

y

<

1), i.e.,

ν

∂u

∂n

−

np

=

0. The

Reynoldsnumberisdefinedby Re

=

U L

/

ν

andthereferencevelocityandlengtharechosenasU

=

1 andL

=

5.Onthe FPflow,weconsiderfourReynoldsnumbersasRe

= {

102

_,

₁₀3

_,

₁₀4

_,

₁₀5

_}

_{,whichcorrespondtotheviscosityparameters}

ν

= {

5

·

10−2

,

5

·

10−3

,

5

·

10−4

,

5

·

10−5

}

,respectively.

Since stretched grid is typically needed to compute the flow accurately at large Reynolds numbers, stretched grid is generated based on the uniform Cartesian grid with 12

×

2n

·

2n cells. The stretching function is applied in the

y-directionwiththeparameterb

=

1

.

01 [cf. [8]]: y

=

(

b

+

1

)

− (

b

−

1

)

c

(

c

+

1

)

,

c

= (

b

+

1 b

−

1

)

1− ¯y

_,

_¯

_y

₌

₀

_,

₁

_/

_n

_,

₂

_/

_n

_{, ...}

₁

_.

₍₁₈₎

(2) Flowoverbackwardfacingstep(BFS)

The L-shapeddomain isknownasthebackward facingstep.APoisseuille flowprofile isimposedon theinflow(x

=

−

1

;

0

≤

y

≤

1).No-slipboundaryconditionsareimposedonthewalls.TheNeumannconditionisappliedattheoutflow (x

=

5

;

−

1

<

y

<

1)whichautomaticallysetstheoutflowpressuretozero.UsingthereferencevelocityandlengthU

=

1 andL

=

2 andtheviscosityparameters

ν

= {

2

·

10−2

_,

₂

_·

₁₀−3

_}

_{,thecorrespondingReynoldsnumbersare Re}

₌

_{U L}

_/

_ν

₌

{

102

,

103

}

.

TheBFSflowismorecomplicatedthantheflat-plateflowasitfeaturesseparation,afreeshear-layerandreattachment. OntheBFSflowwedonotconsidertheReynoldsnumberRe

>

103 _{sincetheincreaseoftheReynoldsnumberbyan}

orderofmagnitudewilltransferthe flowtobe turbulent.Onthiscasewe onlyconsideruniformCartesian gridwith 11

×

2n

_·

₂n _cells.

(3) Liddrivencavity(LDC)

Thisproblemsimulatestheflowinasquarecavity

(

−

1

,

1

)

2 withenclosedboundaryconditions.Alidmovingfromleft torightwithahorizontalvelocityas:

ux

=

1

−

x4 for

−

1

≤

x

≤

1 y

=

1

.

Inordertoaccuratelyresolvethesmallrecirculations,weconsiderstretchedgridaroundthefourcorners.Stretchedgrid isgeneratedbasedontheuniformCartesiangridwith2n

_·

₂n _{cells.Thestretchingfunctionisappliedinbothdirections}

withparametersa

=

0

.

5 andb

=

1

.

01 [8] x

=

(

b

+

2a

)

c

−

b

+

2a

(

2a

+

1

)(

1

+

c

)

,

c

= (

b

+

1 b

−

1

)

¯ x−a 1−a

_,

¯

_x

=

₀

_,

₁

_/

_n

_,

₂

_/

_n

_{, ....,}

₁

_.

₍₁₉₎

ThereferencevelocityandlengthU

=

1 andL

=

2 andtheviscosityparameters

ν

= {

2

·

10−2

_,

₂

_·

₁₀−3

_,

₂

_·

₁₀−4

_}

_resultin

thefollowingReynoldsnumbers Re

= {

102

,

103

,

104

}

.ForthesamereasonasBFS,alargerReynoldsnumber Re

>

104 isnotconsideredonthiscase.

Inordertoexploretheperformanceof

P

M AL and

P

U withvaryingSchurcomplementapproximationsassummarizedin

Table1andTable2,numericalevaluationsareclassifiedintofourcategoriesasfollows. (C1) OnsmallReynoldsnumberanduniformgrid

InthiscategoryweconsidertheFP,BFSandLDCcasesonthesmallReynoldsnumberRe

=

102 _{anduniformCartesian}

(C2) OnmoderateReynoldsnumberanduniformgrid

Inthiscategorywe applythemoderateReynoldsnumberRe

=

103 ontheFP,BFSandLDCcases.Similartothefirst classofexperiments,uniformCartesiangrid isusedheretocheckthe variationofperformance whenincreasingthe Reynoldsnumberbyanorderofmagnitude.

(C3) OnmoderateReynoldsnumbersandstretchedgrid

ThiscategorycontainsthetestscarriedoutontheFPandLDCcaseswithstretchedgrid.Thestretchingfunctionsfor theFPandLDCcasesare (18) and (19), respectively.Still,the moderateReynoldsnumber Re

=

103 _{isemployed for}

thetwotests.Comparingwiththesecondclassofexperiments,thiscategoryismeanttoinvestigatetheeffectofmesh anisotropy.

(C4) OnlargeReynoldsnumbersandstretchedgrid

The LDCcasewith Re

=

104 andFPcasewith Re

= {

104

,

105

}

are includedinthis classof teststoassess howthe Krylovsubspace solver behaves atrelatively large Reynolds numbers.Here stretched grid is employed to accurately resolvetheproblemcharacteristics.

Inthispaperallexperimentsarecarriedoutbasedontheblocks A,B,C ,C1,C2,Ap,Mp,Lp andMu andtheright-hand

side vector rhs,which are obtainedatthe middlestep ofthe wholenonlinear iterations. Numerical experimentsin [13] show that thenumber oflineariterations variesduring thenonlinear procedure.The motivation ofchoosing themiddle step of the nonlinear iterations to export the blocks and vector is that a representative number of linear iteration can be obtained, compared to the averaged numberof linear iterations through the whole nonlinearprocedure. The relative stopping tolerance to solve the linear systemby GMRES ischosen equal to 10−8_{. The restart functionality ofGMRES is}

not usedinthispaper. Sincethepreconditioners

P

M AL and

P

U involvevariousmomentum andpressuresub-systems,all

thesesub-systems are directlysolved inthis paperto avoidthesensitiveness ofiterative solverson the varyingsolution complexities.

AspointedoutinSection2,theapplicationoftheSchurcomplementapproximation



Sγ P C D needstopresetboundary

conditions forthepressure Laplacian Lp and convection-diffusion Ap operators. Inthis paper,we follow thesuggestions

of [32,33] to use Dirichletboundary conditionsalong inflow boundaries todefine Lp and Ap.This means that the rows

andcolumns of Lp and Ap corresponding tothe pressure nodes on an inflow boundary are treatedas though they are

associatedwithDirichletboundaryconditions.Fortheenclosedflow,wealgebraicallyaddh2I to Lp andAp tomakethem

non-singular,whereh denotes thegrid sizeand I is theidentity matrixofproper size.Such artificialpressureboundary conditionsare onlyimposed on thepreconditioner. Thecoefficient matrix andright-handside vector are notaffected by theseboundarynodemodifications.

4.1. OnsmallReynoldsnumberanduniformgrid

InthissubsectionwecarryoutexperimentsontheFP,BFSandLDCcaseswithuniformCartesiangridandsmallReynolds number Re

=

102_{. The numberof Krylov iterations to solve the transformed system preconditioned by the modified AL}

preconditioner

P

M AL isgiveninTable 3.The Schurcomplementapproximations



Sγ P C D,



Sγ L SC,



Sγ S I M P L E in

P

M AL are

derived fromthe newmethod



Sγ new (16) andtheapproximation



Sγ orig corresponds to theoriginal Schur complement

approximation(17).Inthispaper,thereportednumberofKryloviterationspreconditionedby

P

M AL isobtainedbyusingthe

optimalvalue of

γ

,whichresultsinthefastestconvergencerateoftheKrylovsubspacesolver.Thefollowingobservations aremadefromTable3.

Except



Sγ S I M P L E,weseethattheotherSchur complementapproximationsresultintheindependenceofKrylov

itera-tionsonthemeshrefinementatthethreetestcases.IntermsofthenumberofKryloviterations,



Sγ L SC issuperiortothe

otherSchurcomplementapproximationsontheFPandBFScasesbythereducednumberofiterationsandequallyefficient as



Sγ P C D and



Sγ orig onthe LDC case. Tounderstand thisadvantage, we take the FP caseasan exampleandplot the

eigenvalues ofthepreconditionedSchur complementmatrix



S−_{γ Sγ in}1 Fig.1.As canbe seen,



Sγ L SC leads tomore

clus-teredeigenvaluesandthesmallesteigenvaluefurtherawayfromzero.Suchadistributionofeigenvaluesisfavorableforthe Krylovsubspacesolverandafasterconvergenceratecanbeexpected. Weknowthattherecanbematriceswherethereis norelationbetweenthespectrumandtheconvergenceofGMRES[37],especiallyifthematrixisstronglynonnormal.We includethespectrumbecauseinourexamplesthepropertiesofthespectrumareinlinewiththeconvergenceproperties ofGMRES.Inaddition,thefield-of-valuestypeestimatesfortheaugmentedLagrangianpreconditionedmatrixareprovided by[38].

As analyzedin Section 3.3,ateach Krylov iterationthe costs ofapplying

P

M AL with



Sγ L SC are roughly thesame as



Sγ P C D and two timesof that using



Sγ S I M P L E and



Sγ orig. Ifwe assume the computational expenseof applying

P

M AL

with



Sγ orig tobe unit ateachiteration,the totalcostsby usingall Schurcomplementapproximationsonthefinestgrid

arepresentedinTable4andcalculatedbymultiplyingtheexpenseperiterationby thenumberofiterations.Intheother classesofevaluationswealsousethismethodtocalculatethetotalcomputationalcosts.

Results inTable4 show thatthe minimal computationalcosts are achievedby using



Sγ orig in

P

M AL.Althoughfewer

Table 3

Re=102_{and uniformgrid:thenumberofGMRESiterationstosolvethetransformedsystemwithA}

γ

precon-ditionedbyPM ALwithdifferentSchurcomplementapproximationsandtheoptimalvalueofγinparentheses.

Sγ P C D SγL SC Sγ S I M P L E Sγ orig FP case:

n=5 26(1.e-1) 17(8.e-2) 43(2.e-1) 38(2.e-1) n=6 25(1.e-1) 25(8.e-2) 67(2.e-1) 38(2.e-1) n=7 25(1.e-1) 26(8.e-2) 100(2.e-1) 38(2.e-1) BFS case:

n=5 34(2.e-2) 17(2.e-2) 42(1.e-1) 36(1.e-1) n=6 42(3.e-2) 21(2.e-2) 60(1.e-1) 36(1.e-1) n=7 45(3.e-2) 22(2.e-2) 87(1.e-1) 36(1.e-1) LDC case:

n=6 17(2.e-2) 17(2.e-2) 34(1.e-1) 19(1.e-1) n=7 18(2.e-2) 20(2.e-2) 48(1.e-1) 19(1.e-1) n=8 18(2.e-2) 22(2.e-2) 63(1.e-1) 19(1.e-1)

Table 4

Re=102_{and uniformgrid:thetotalcostsofapplyingP}

M AL withdifferentSchurcomplementapproximations

onthefinestuniformCartesiangrid.

Sγ P C D Sγ L SC Sγ S I M P L E Sγorig

FP case: 50 52 100 38

BFS case: 90 44 87 36

LDC case: 36 44 63 19

theheaviercostsof



Sγ L SC.Inthisclassofexperiments,itseemsthattheoriginalSchurcomplementapproximation



Sγ orig

ismoreefficientthantheotherapproximationsduetothefewercomputationalcostsintotal.

4.2. OnmoderateReynoldsnumberanduniformgrid

Inthissubsectionwe choosethemoderateReynoldsnumberRe

=

103 _{toevaluatethe performanceoftheSchur}

com-plement approximationsusedin themodified ALpreconditioner

P

M AL andcomparewiththeevaluations at Re

=

102 in

Section4.1.BasedonthenumberofKryloviterationspresentedinTable5,weseethattheindependenceofKryloviterations on themeshrefinement isachievedby usingtheSchur complementapproximations



Sγ P C D and



Sγ L SC in

P

M AL,which

is alsoobservedinSection 4.1.Contrarytothe observationsinSection 4.1,the originalSchur complementapproximation



_Sγ orig doesnotresultin themeshindependenceof Kryloviterations at Re

=

103.Withtheutilizationof



Sγ S I M P L E the

numberofKryloviterationsisdependentofthegridsizeatboth Re

=

102and103.

Results inTable5showthatthesmallestnumberofKryloviterations isobtainedbyusing



Sγ L SC in

P

M AL,whichalso

resultsintheminimaltotalcostsinTable6.ThetotalcostsarecalculatedbyusingthesamemethodasSection4.2.Taking themeshindependenceintoaccount,theutilizationof



Sγ L SC willleadtoa furtherreduction oftotalcosts onfinergrids

over



Sγ S I M P L E and



Sγ orig,whichrequiremoreiterationswithmeshrefinement.Comparedto



Sγ P C D whichalsoresults

inthemeshindependenceofKryloviterations,theapplicationof



Sγ L SC reducesthetotalcomputationalcostsatleasttwo

times onthe FPandBFScases,andthisreduction factorcanalso be expectedonfinergrids. Onthe LDCcase



Sγ L SC is

equallyefficientas



Sγ P C D.

For thetestsat Re

=

103 _{it showsthat}



_Sγ

L SC issuperior to theother Schur complementapproximationsby the

re-ductionofKryloviterations andtotalcomputationalcosts.IntheprevioustestswithRe

=

102,thesuperiorityof



Sγ origis

seen.ThisimpliesthattheoptimalSchurcomplementapproximationdependsontheReynoldsnumber.

4.3. OnmoderateReynoldsnumberandstretchedgrid

This subsectionismeant to investigatetheinfluence ofmesh anisotropyon theperformance of themodified AL pre-conditioner

P

M AL.TocomparewithSection4.2,weapplythestretched gridandmoderateReynoldsnumberRe

=

103 on

theFPandLDCcases.ThenumberofKryloviterationsandtotalcomputationalcostsarepresentedinTable7andTable8, respectively.From Table7wenotethatonly



Sγ P C D resultsinthemeshindependenceandtheminimalnumberofKrylov

iterations.Althoughthetotalcostsofapplying



Sγ P C D aremorethanthatbyusing



Sγ S I M P L Eand



Sγ origontheconsidered

finestgrid,asseenfromTable8,fewer costsintotalbyusing



Sγ P C D canbeexpectedonfinergridsduetothemesh

in-dependence.Therefore,wethinkthat



Sγ P C D issuperiortotheotherSchurcomplementapproximationsonthetestswith

Re

=

103 andstretchedgrid.

Note that on the FP and LDC cases with stretched grid,

P

M AL with



Sγ L SC is not mesh independent any more and

Fig. 1. FP and Re=102_{: plot of eigenvalues of the preconditioned matricesS}−1

γ Sγ at the uniform Cartesian grid with 12×25·25cells.

Table 5

Re=103_{and uniformgrid:thenumberofGMRESiterationstosolvethetransformedsystemwithA}

γ

precon-ditionedbyPM ALwithdifferentSchurcomplementapproximationsandtheoptimalvalueofγinparentheses.

Sγ P C D Sγ L SC Sγ S I M P L E Sγorig FP case:

n=5 54(8.e-3) 29(8.e-3) 34(2.e-2) 76(6.e-2) n=6 55(8.e-3) 18(8.e-3) 51(2.e-2) 90(6.e-2) n=7 56(8.e-3) 17(8.e-3) 99(2.e-2) 95(6.e-2) BFS case:

n=5 66(4.e-3) 45(3.e-3) 49(1.e-2) 71(3.e-2) n=6 63(4.e-3) 27(3.e-3) 77(1.e-2) 76(3.e-2) n=7 65(3.e-3) 29(3.e-3) 142(1.e-2) 84(3.e-2) LDC case:

n=6 30(4.e-3) 54(1.e-3) 66(7.e-3) 36(2.e-2) n=7 28(4.e-3) 29(4.e-3) 52(1.e-2) 42(2.e-2) n=8 29(4.e-3) 29(4.e-3) 85(1.e-2) 48(2.e-2)

FPcaseasanexample,onthefineststretchedgridofn

=

7 thenumberofKryloviterationspreconditionedby

P

M AL with



_Sγ L SC increasesbyafactorabout7comparedtothefinestuniformgrid.Thiscanbeseenbycomparingthecorresponding

resultsinTable5andTable7.Thelessefficiencyof

P

M AL with



Sγ L SC arisingfromthemeshanisotropyisalsoseenonthe

LDCcase.Onthe otherhand,thenumberofKryloviterations preconditionedby

P

M AL withtheotherSchur complement

Table 6

Re=103_{and uniformgrid:thetotalcostsofapplyingP}

M AL withdifferentSchurcomplementapproximations

onthefinestuniformCartesiangrid.

Sγ P C D Sγ L SC Sγ S I M P L E Sγorig

FP case: 112 34 99 95

BFS case: 130 58 142 84

LDC case: 58 58 85 48

Table 7

Re=103_{and stretchedgrid:thenumberofGMRESiterationstosolvethetransformedsystemwithA}

γ

precon-ditionedbyPM ALwithdifferentSchurcomplementapproximationsandtheoptimalvalueofγinparentheses.

Sγ P C D Sγ L SC Sγ S I M P L E Sγ orig FP case:

n=5 59(8.e-3) 90(7.e-3) 37(2.e-2) 69(6.e-2) n=6 66(8.e-3) 89(7.e-3) 63(2.e-2) 85(6.e-2) n=7 62(8.e-3) 117(6.e-3) 119(2.e-2) 92(6.e-2) LDC case:

n=6 65(2.e-3) 98(2.e-3) 57(7.e-3) 69(1.e-2) n=7 41(2.e-3) 58(2.e-3) 46(7.e-3) 40(1.e-2) n=8 38(2.e-3) 84(2.e-3) 75(7.e-3) 54(1.e-2)

Table 8

Re=103_{and stretchedgrid:thetotalcostsofapplyingP}

M ALwithdifferentSchurcomplementapproximations

onthefineststretchedgrid.

Sγ P C D Sγ L SC Sγ S I M P L E Sγorig

FP case: 124 234 119 92

LDC case: 76 168 75 54

Fig. 2. FP and Re=103_{: plot of eigenvalues of the preconditioned matrices}_S−1

γ L SCSγ at the uniform and stretched grids with 12×25·25cells.

The lessefficiencyof



Sγ L SC onthestretchedgridcan beexplainedbytheresultsinFig.2,whereweconsidertheFP

caseatRe

=

103 _{andplottheeigenvaluesofthepreconditionedmatrix}



_S−1

γ L SCSγ forbothuniformandstretched grids.As

seen fromFig.2,stretching thegrid considerablyspreadsthe distributionofthe eigenvaluesofthe preconditionedSchur complement



S−_γ1_{L SC}Sγ ,whichmakestheconvergenceoftheKrylovsubspacesolvermoredifficult.

4.4. OnlargeReynoldsnumberandstretchedgrid

In thissubsectionwe applylarge Reynoldsnumbers Re

≥

104 _{andstretched grids ontheLDCandFPcases.Results in}

Table9andTable10illustratethatthefastestconvergencerateoftheKrylovsubspacesolverandtheminimalcomputational costs intotal areachievedby using



Sγ S I M P L E in

P

M AL onthetwo tests. TakingtheFP caseat Re

=

104 asan example,

fromTable10weseethattheutilizationof



Sγ S I M P L E reducesthetotalcostsatleasttwotimeswithrespecttotheother

Table 9

Re=104_{and stretchedgrid:thenumberofGMRESiterationstosolvethetransformedsystemwithA}

γ

precon-ditionedbyPM ALwithdifferentSchurcomplementapproximationsandtheoptimalvalueofγinparentheses.

Sγ P C D Sγ L SC Sγ S I M P L E Sγ orig FP case:

n=5 363(8.e-4) 369(6.e-4) 35(2.e-3) 93(1.e-2) n=6 334(8.e-4) 336(6.e-4) 53(3.e-3) 128(2.e-2) n=7 346(8.e-4) 374(6.e-4) 83(4.e-3) 192(2.e-2) LDC case:

n=6 113(3.e-4) 97(2.e-4) 34(1.e-3) 46(5.e-3) n=7 143(3.e-4) 235(2.e-4) 45(1.e-3) 65(5.e-3) n=8 159(4.e-4) 309(2.e-4) 80(2.e-3) 106(5.e-3)

Table 10

Re=104_{and stretchedgrid:thetotalcostsofapplyingP}

M AL withdifferentSchurcomplementapproximations

onthefineststretchedgrid.

Sγ P C D Sγ L SC Sγ S I M P L E Sγ orig

FP case: 692 748 83 192

LDC case: 318 618 80 106

Table 11

FP andRe=105_{:the numberofGMRESiterationsandtotalcoststosolvethetransformedsystemwithA}

γ

preconditionedbyPM ALwithdifferentSchurcomplementapproximationsandtheoptimalvalueofγin

paren-theses.Thestretchedgridisapplied.

Sγ P C D SγL SC Sγ S I M P L E Sγ orig iterations: n=5 1000+ 1000+ 26(1.e-4) 136(1.e-3) n=6 1000+ 1000+ 35(2.e-4) 192(2.e-3) n=7 1000+ 1000+ 58(3.e-4) 310(2.e-3) total costs: n=7 2000+ 2000+ 58 310

onthe FPcase, whichis seenfromTable 11.Inthe contextoflarge Reynoldsnumbers, itappears that



Sγ S I M P L E isthe

optimalSchurcomplementapproximationinthemodifiedALpreconditioner

P

M AL.Incontrasttotheprevioustests,atlarge

ReynoldsnumbersnoneoftheconsideredSchurcomplementapproximationsleadtothemeshindependenceof

P

M AL.The

advantageof



Sγ S I M P L E onfinergridsneedsafurtherassessment,whichisincludedinfutureresearch.

To investigatethe effect ofthe Reynolds number, we take the FP caseasan exampleand in Fig.3 plot the number ofKrylov iterations preconditioned by

P

M AL atvarying Reynolds numbers.It appears that only



Sγ S I M P L E results inthe

robustnessof

P

M AL withrespecttotheReynoldsnumber.Tounderstandthereasons,wecomputetheextremaleigenvalues

ofthepreconditionedSchurcomplementmatrix



S−_{γ Sγ and}1 presenttheminTable12. Rmin andRmax denotethesmallest

andlargestrealpartsoftheeigenvaluesandImax correspondsthelargestimaginarypart.Theseextremalvaluescorrespond

totheboundariesoftherectangulardomaincontainingalleigenvalues.Regarding



S−_γ1_{S I M P L E}Sγ ,thevaluesofRmin slightly

decreaseandremainthesameorderofmagnitude.Togetherwiththedecreaseof Rmax

/

Rmin andImax,theeigenvaluesare

furtherclustered. However,fewerclusteredeigenvaluesare yieldedby usingtheotherSchurcomplementapproximations. Thisexplainstherobustnessof

P

M AL with



Sγ S I M P L E withrespecttotheReynoldsnumber.

Toinvestigatethe computedsolutionsatlargeReynoldsnumbers,we choosetheFPcase.Intheinviscidlimit Re

→ ∞

thesolutionissimplyux

=

1,uy

=

0 and p

=

constant.Sincetheshearboundarylayerisofwidthproportionalto

√

ν

and

within the layerthe horizontalvelocity increases rapidlyfromzero tounity, the plateseems “invisible”as Re

→ ∞

[1]. Tocheckthisfeature,inFig.4weillustratethecalculatedpressureandequallyspacedcontours ofthehorizontalvelocity between0 and0

.

95 atdifferentReynoldsnumbers.Thestretched gridwith12

×

26

_·

₂6 _{cells isutilized.At Re}

₌

₁₀3_,the

countersofthehorizontalvelocityshowtheevolutionoftheboundarylayerasthefluidpassingfromtheleadingedgeof the platetothe outflow.The parabolicshape ofthe velocity contours seemsconsistent withasymptotic theory[39] and thereportedresultsin[1].WhenincreasingtheReynoldsnumberstoRe

=

105_{,weseethattheplate“disappears”andthe}

differencebetweenthepressurevaluesdecreasesbyoneorderofmagnitudecomparedtothecaseofRe

=

103.Resultsin Fig.4demonstratethatthecomputedsolutions,rangingfromthemoderatetolargeReynoldsnumbers,seemreasonable.

4.5. SummaryoftheSchurcomplementapproximationsin

P

M AL

Based on the above fourclasses of numericalevaluations, inTable 13 we summarizethe optimal Schur complement approximationinthemodifiedALpreconditioner

P

M AL.ItshowsthattheoptimalSchurcomplementapproximation,which

Fig. 3. FP and stretched grid: plot of the number of GMRES iterations preconditioned byPM ALat varying Reynolds numbers.

Ateveryclassofevaluations,theoptimalSchurcomplementapproximationisproblemindependent.Numericalevaluations inthispapershowthat



Sγ origissuitableforthecalculationswithsmallReynoldsnumbersand



Sγ S I M P L E deliversabetter

performance forlarge Reynolds numbers due to its Reynolds robustness. In the context ofmoderate Reynolds numbers,



Sγ L SC ismoreefficientwithuniformgrids butsensitiveto meshanisotropy.Whenstretched gridsareemployed,



Sγ P C D

turns out to be the optimalchoice inthe moderateReynolds numbercontext. Exceptthe calculationsatsmall Reynolds numbers anduniformgrids,theoptimalSchurcomplementapproximationsonother classesoftestsarederived fromthe newmethod



Sγ newproposedinthispaper.Thisdemonstratestheadvantageofthenewapproachoverthetraditionalone



_{Sγ orig}.ThemeshindependenceofKryloviterationsisnotachievedbyusingtheoptimalSchurcomplementapproximation only forthe class of tests withlarge Reynolds numbers. The reason and possible improvement on this issueare to be consideredinfutureresearch.

4.6. Comparisonbetween

P

M ALand

P

U.

To apply the modified AL preconditioner

P

M AL,one needs totransform the original system(2) to an equivalent one

(13) withthecoefficient matrix

A

γ .Thistransformationconsumesadditionalcosts.Furthermore,ateachKryloviteration

extra costs arisefromthe productof

A

γ with aKrylovresidual vectordueto morefill-in in

A

γ [13]. Inthissense, the

heavier complexities of

P

M AL could be payed offonly by a reducednumberof Kryloviterations, compared to theblock

upper-triangularpreconditioner

P

U applied tothe originalsystem. In thissection,we considerthe comparisonsbetween

P

M AL and

P

U on theLDCandFP casesatthe largeReynolds number Re

=

104 and stretchedgrid whichrepresent stiff

testsontheconsideredpreconditioners.

ItisrevealedinSection 4.4that



Sγ S I M P L E turnsouttobethemostefficientSchurcomplementapproximationforthe

Fig. 4. FP and stretchedgrid:plotofthecalculatedpressureunknown(left)andcontoursofthehorizontalvelocitybetween0 and0.95 (right)atdifferent Reynoldsnumbers.