IFISS, a Matlab toolbox for modelling incompressible flow Elman, Howard C. and Ramage, Alison ...

19 Pages·2017·0.23 MB·English

Algorithm 866: IFISS, a Matlab toolbox for
modelling incompressible flow
Elman, Howard C. and Ramage, Alison and Silvester,
David J.
2007
MIMS EPrint: 2005.45
Manchester Institute for Mathematical Sciences
School of Mathematics
The University of Manchester
Reports available from: http://eprints.maths.manchester.ac.uk/
And by contacting: The MIMS Secretary
School of Mathematics
The University of Manchester
Manchester, M13 9PL, UK
ISSN 1749-9097
Algorithm 866: IFISS, A Matlab Toolbox for
Modelling Incompressible Flow
HOWARDC.ELMAN
UniversityofMaryland
ALISONRAMAGE
UniversityofStrathclyde
and
DAVIDJ.SILVESTER
UniversityofManchester
IFISSisagraphicalMatlabpackagefortheinteractivenumericalstudyofincompressibleflow
problems.Itincludesalgorithmsfordiscretizationbymixedfiniteelementmethodsandaposteriori
error estimation of the computed solutions. The package can also be used as a computational
laboratoryforexperimentingwithstate-of-the-artpreconditionediterativesolversforthediscrete
linearequationsystemsthatariseinincompressibleflowmodelling.Auniquefeatureofthepackage
isitscomprehensivenature;foreachproblemaddressed,itenablesthestudyofbothdiscretization
anditerativesolutionalgorithmsaswellastheinteractionbetweenthetwoandtheresultingeffect
onoverallefficiency.
Categories and Subject Descriptors: G.1.3 [Numerical Analysis]: Numerical Linear Algebra—
Linearsystems(directanditerativemethods);G.1.8[NumericalAnalysis]:PartialDifferential
Equations—Ellipticequations,finiteelementmethods,iterativesolutiontechniques,multigridand
multilevelmethods;G.4[MathematicalSoftware];J.2[PhysicalSciencesandEngineering]
GeneralTerms:Algorithms,Design,Documentation
AdditionalKeyWordsandPhrases:Matlab,finiteelements,incompressibleflow,iterativesolvers,
stabilization
The work of H. C. Elman was supported by the U.S. National Science Foundation under grant
DMS0208015andbytheU.S.DepartmentofEnergyundergrantDOEG0204ER25619.
Authors’addresses:H.C.Elman,Dept.ofComputerScience,UniversityofMaryland,CollegePark,
MD20742;email:[email protected];A.Ramage,Dept.ofMathematics,UniversityofStrathclyde,
26RichmondStreet,GlasgowG11XH,UK;email:[email protected];D.J.Silvester,School
of Mathematics, University of Manchester, Sackville Street, Manchester M60 1QD, UK; email:
[email protected].
Permissiontomakedigitalorhardcopiesofpartorallofthisworkforpersonalorclassroomuseis
grantedwithoutfeeprovidedthatcopiesarenotmadeordistributedforprofitordirectcommercial
advantageandthatcopiesshowthisnoticeonthefirstpageorinitialscreenofadisplayalong
withthefullcitation.CopyrightsforcomponentsofthisworkownedbyothersthanACMmustbe
honored.Abstractingwithcreditispermitted.Tocopyotherwise,torepublish,topostonservers,to
redistributetolists,ortouseanycomponentofthisworkinotherworksrequirespriorspecificper-
missionand/orafee.PermissionsmayberequestedfromthePublicationsDept.,ACM,Inc.,2Penn
Plaza,Suite701,NewYork,NY10121-0701USA,fax[phone],[email protected].
(cid:2)C 2007ACM0098-3500/2007/06-ART14$5.00DOI10.1145/1236463.1236469http://doi.acm.org/
10.1145/1236463.1236469
ACMTransactionsonMathematicalSoftware,Vol.33,No.2,Article14,Publicationdate:June2007.
2 • H.C.Elmanetal.
ACMReferenceFormat:
Elman,H.C.,Ramage,A.,andSilvester,D.J.2007.Algorithm866:IFISS,amatlabtoolboxfor
modellingincompressibleflow.ACMTrans.Math.Softw.,33,2,Article14(June2007),18pages.
DOI=10.1145/1236463.1236469 http://doi.acm.org/10.1145/1236463.1236469
1. INTRODUCTION
ThisarticledescribestheIncompressibleFlowIterativeSolutionSoftware(IFISS)
package, a collection of over 290 matlab functions and m-files organised as a
matlabtoolbox.Thetoolboxcanbeusedtocomputenumericalsolutionsofpar-
tialdifferentialequations(PDEs)thatareusedtomodelsteadyincompressible
fluid flow. It includes algorithms for discretization by finite element methods,
fast iterative solution of the algebraic systems that arise from discretization,
andaposteriorierroranalysisofthecomputeddiscretesolutions.Thepackage
wasproducedinconjunctionwiththemonographFiniteElementsandFastIter-
ativeSolversbyElmanetal.[2005b],anditwasusedtoperformthecomputa-
tionalexperimentsdescribedtherein.Itisstructuredasastand-alonepackage
forstudyingdiscretisationalgorithmsforPDEsandforexploringanddevelop-
ingalgorithmsinnumericallinearandnonlinearalgebraforsolvingtheasso-
ciated discrete systems. It can also be used as a pedagogical tool for studying
theseissues,ormoreelementaryonessuchasthepropertiesofKrylovsubspace
iterativemethods.
Four PDEs are treated: the Poisson equation, steady-state versions of the
convection-diffusion equation, Stokes equations, and Navier-Stokes equations.
Thefirsttwoequationsareubiquitousinscientificcomputing,seeforexample
Ockendonetal.[1999],Milleretal.[1995],Morton[1996]orRoosetal.[1996].
The latter two PDEs constitute the basis for computational modelling of the
flowofanincompressibleNewtonianfluid.
The first main feature of the package concerns problem specification and
finiteelementdiscretization.Foreachoftheequationslisted,IFISSoffersachoice
of two-dimensional domains on which the problem can be posed, along with
boundary conditions and other aspects of the problem, and a choice of finite
element discretizations on a quadrilateral element mesh. The package allows
thestudyof
—accuracyoffiniteelementsolutions,
—differentchoicesofelements,
—aposteriorierroranalysis.
In addition, special features associated with individual problems can be ex-
plored.Theseincludetheeffectsofboundarylayersonsolutionqualityforthe
convection-diffusion equation, and the effects of discrete inf-sup stability on
accuracyfortheStokesandNavier-Stokesequations.
The second main feature of the package concerns iterative solution of the
discrete algebraic systems, with emphasis on preconditioned Krylov subspace
methodsforlinearsystemsofequations.TheKrylovsubspacemethodsarecho-
sentomatcheachproblem.Forexample,thediscretePoissonequation,which
ACMTransactionsonMathematicalSoftware,Vol.33,No.2,Article14,Publicationdate:June2007.
Algorithm866:IFISS,AMatlabToolboxforModellingIncompressibleFlow • 3
hasasymmetricpositivedefinitecoefficientmatrix,canbetreatedbythecon-
jugate gradient method (CG) [Hestenes and Stiefel 1952] whereas the discrete
convection-diffusion and Navier-Stokes equations require a method such as
the generalized minimum residual method (GMRES) [Saad and Schultz 1986],
whichisdesignedfornonsymmetricsystems(see§3.1fordetails).Thekeyfor
fastsolutionliesinthechoiceofeffectivepreconditioningstrategies.Thepack-
age offers a range of options, including algebraic methods such as incomplete
LUfactorizations,aswellasmoresophisticatedandstate-of-the-artmultigrid
methodsdesignedtotakeadvantageofthestructureofthediscretelinearized
Navier-Stokesequations.Inaddition,thereisachoiceofiterativestrategies,Pi-
carditerationorNewton’smethod,forsolvingthenonlinearalgebraicsystems
arisingfromthelatterproblem.
AuniquefeatureoftheIFISSpackageisitscomprehensivenature,wherefor
eachproblemitaddresses,itenablesthestudyofbothdiscretizationanditer-
ative solution algorithms as well as the interaction between the two and the
resultingeffectonsolutioncost.Itisorganizedinamodularfashion,withsep-
arate program units (typically matlab functions) for individual tasks such as
discretization,erroranalysis,iterativesolution,andotheraspectssuchasprob-
lem definition and visualization tools grouped in separate subdirectories. The
packageprovidesaconvenientstartingpointforresearchprojectsthatrequire
theconstructionofnewproblemclasses,discretizationsorsolutionalgorithms.
A detailed description of how new problems can be incorporated into IFISS is
giveninElmanetal.[2005b].
Theremainderofthearticleisorganizedasfollows.Section2identifiesthe
PDEstreated,describesthefiniteelementmethodsusedtoapproximatethem,
and outlines the strategies used for a posteriori error estimation. Section 3
describestheiterativesolutionalgorithmsprovidedforthelinearsystemsthat
ariseforeachofthebenchmarkproblems.Finally,Section4describessomeof
thedesignfeaturesofthepackage.
2. APPROXIMATION
A key feature of the IFISS package is the facility to construct various finite el-
ement solutions for a range of PDE problems and to carry out a posteriori
erroranalysisofthecomputeddiscretesolutions.Thetheoryoffiniteelement
methodsiswellestablished(seeforexampleBraess[1997];BrennerandScott
[1994]; Ciarlet [1978]; Elman et al. [2005b]; Gunzburger [1989]; and Johnson
[1987]).AnoutlineofthespecificapproximationmethodsimplementedinIFISS
ispresentedbelowforeachPDEinturn.
2.1 ThePoissonEquation
Theproblemconsideredhereis
−∇2u = f in (cid:2), (1)
∂u
u= g on∂(cid:2) , =0on∂(cid:2) ,
D D ∂n N
ACMTransactionsonMathematicalSoftware,Vol.33,No.2,Article14,Publicationdate:June2007.
4 • H.C.Elmanetal.
where (cid:2) is a two-dimensional domain with boundary ∂(cid:2) consisting of two
nonoverlappingsegments∂(cid:2) ∪∂(cid:2) .Thefunction f isagivensource(orload)
D N
functionand g isgivenDirichletboundarydata.Inaddition, (cid:5)nistheoutward-
D
pointing normal to the boundary and ∂u denotes the directional derivative in
∂n
thenormaldirection.TheweakformulationofthePDEproblemhasasolution
u in the Sobolev space H1((cid:2)) that satisfies the Dirichlet boundary condition
togetherwiththeconditionthat
(∇u,∇v)=(f,v) ∀v∈ V :={v∈H1((cid:2))|v=0on∂(cid:2) },
D
where (·,·) denotes the L inner product on (cid:2). Given some finite-dimensional
2
subspace V ⊂ V,theassociateddiscretesolutionu satisfiesu = g on∂(cid:2)
h h h D D
and
(∇u ,∇v )=(f ,v ) ∀v ∈ V , (2)
h h h h h h
where f isthe L ((cid:2))orthogonalprojectionof f into V .Forastandardfinite
h 2 h
elementgridwithnuinteriornodesandn∂ Dirichletboundarynodes,wechoose
a set of basis functions φ , i = 1,...,n for V , and look for an approximate
i u h
solutionoftheform
(cid:2)nu n(cid:2)u+n∂
u (x, y)= u φ (x, y)+ u φ (x, y), (3)
h i i i i
i=1 i=nu+1
wher(cid:3)etheboundarycoefficientsui,i =nu+1,...,nu+n∂ arechosentoensure
that nu+n∂ u φ (x, y)interpolatestheboundarydata g on∂(cid:2) .Thisleads
toanni=×nu+n1 siystiemoflinearequations D D
u u
Au=f, (4)
wheretheentriesofthesolutionuaretheunknowncoefficientsu in(3).
i
Two standard finite element approximation methods are implemented in
IFISS:
—Bilinear quadrilateral Q : The isoparametrically mapped square element
1
with bilinear basis functions of the form (ax + b)(cy + d). There are four
unknownsperelement.
—Biquadratic quadrilateral Q : The bilinearly mapped square element with
2
biquadraticbasisfunctionsoftheform(ax2+bx+c)(dy2+ey + f).There
arenineunknownsperelement,correspondingtotheterms1, x, y, x2, xy,
y2, xy2, x2y and x2y2.
These finite element approximations typically have a discontinuous normal
derivative across interelement boundaries. Consequently, it is convenient to
definethefluxjumpacrossedge E adjoiningelementsT and S as
(cid:4)(cid:4) (cid:5)(cid:5)
∂v
:=(∇v| −∇v| )· (cid:5)n =(∇v| −∇v| )· (cid:5)n , (5)
∂n T S E,T S T E,S
where (cid:5)n is the outward normal with respect to E and ∇u · (cid:5)n is the
E,T h E,T
discrete (outward-pointing) normal flux. An a posteriori estimate of the dis-
cretizationerrore=u−u maybecomputedfromtheequidistributedinterior
h
ACMTransactionsonMathematicalSoftware,Vol.33,No.2,Article14,Publicationdate:June2007.
Algorithm866:IFISS,AMatlabToolboxforModellingIncompressibleFlow • 5
edgefluxjump R := 1[[∂uh]], togetherwiththeinteriorelementresidual
E 2 ∂n
R :={f +∇2u }| .
T h T
Thisisaconsequenceofthefactthattheerrorsatisfies
(cid:6) (cid:7)
(cid:2) (cid:2) (cid:2)
(∇e,∇v) = (R ,v) − (cid:9)R ,v(cid:10) (6)
T T T E E
T∈Th T∈Th E∈E(T)
for all v ∈ V, where E(T) is the set of edges associated with element T; (cid:9)·,·(cid:10)
E
denotesthe L innerproducton E;andT isthesetofallelements.SeeElman
2 h
etal.[2005b,p.50]forfurtherdetails.
TheIFISSsoftwarecomputesanaposteriorierrorestimatewheneverQ1 ap-
proximationisused.Thisisparticularlystraightforwardtoimplement,inthat
R ispiecewiselinearand R = f| canbeapproximatedbyapiecewisecon-
E T T
stantfunctionR0 byevaluating f attheelementcentroid.Considerthehigher
T
ordercorrectionspace
Q = Q ⊕B ,
T T T
sothat Q isthespacespannedbybiquadraticedgebubblesψ ,and B isthe
T E T
spacespannedbyinteriorbiquadraticbubblesφ .Theenergynormoftheerror
T
isestimatedbysolvinga5×5localPoissonproblemposedovereachelement
ofthegrid1:
(cid:8) (cid:9) (cid:2)
(∇e ,∇v) = R0,v − (cid:9)R ,v(cid:10) ∀v∈Q .
T T T T E E T
E∈E(T)
Oncee hasbeencomputed,thelocal,elementwise,errorestimatorisgivenby
T
η =(cid:12)∇e (cid:12) ,andtheglobalerrorestimatorisgivenby
T T T
(cid:10) (cid:11)
(cid:2) 1/2
η:= η2 ≈(cid:12)∇(u−u )(cid:12).
T h
T∈Th
This approach for error estimation was originally introduced in Bank and
Weiser[1985].Itbelongstoaclassofmethodsreferredtoasimplicitestimators,
seeAinsworthandOden[2000].Ananalysisofitseffectivenessasanestimator
canbefoundinElmanetal.[2005b].
2.2 TheConvection-DiffusionEquation
Theequation
−(cid:7)∇2u+ w(cid:5) ·∇u=0 in (cid:2), (7)
(with(cid:7) >0)arisesinnumerousmodelsofflowsandotherphysicalphenomena.
Inatypicalapplication,theunknownfunctionurepresentstheconcentration
of a pollutant being transported (or ‘convected’) along a stream moving at ve-
locity, w(cid:5) andalsosubjecttodiffusiveeffects.Theconvection-diffusionequation
1Thelocalproblemdefinitionneedstobeslightlymodifiedforthoseelementshavingoneormore
edgesontheDirichletboundary∂(cid:2)D.
ACMTransactionsonMathematicalSoftware,Vol.33,No.2,Article14,Publicationdate:June2007.
6 • H.C.Elmanetal.
considered in IFISS is (7) posed on a two-dimensional domain (cid:2), together with
theboundaryconditions
∂u
u= g on∂(cid:2) , =0on∂(cid:2) . (8)
D D ∂n N
Wewillrefertothevelocityvector w(cid:5) asthewind.
ItiswellknownthatapplyingtheGalerkinfiniteelementmethodto(7)and
(8)willoftenresultinadiscretesolutionthatexhibitsnonphysicaloscillations.
This unwanted effect is minimized in IFISS through use of the streamline dif-
fusionmethodofHughesandBrooks[1979]:thediscreteproblemisstabilized
by adding diffusion in the streamline direction. The resulting discrete formu-
lation, which is the analogue of (2) for the convection-diffusion equation with
discretestabilization,is
(cid:2)
(cid:7)(∇u ,∇v )+(w(cid:5) ·∇u ,v ) + δ (w(cid:5) ·∇u , w(cid:5) ·∇v ) =0 ∀v ∈ V , (9)
h h h h k h h k h h
k
wherethesumistakenoverallelements inthegrid.
k
To implement (9), we need a way of choosing locally defined parameters δ .
k
OnewayofdoingthisissuggestedbytheanalysisinFischeretal.[1999]and
ElmanandRamage[2002].Specifically,wechoose
⎧ (cid:15) (cid:16)
⎨ hk 1− 1 ifPk >1
δ = 2|w(cid:5)k| Phk h (10)
k ⎩
0 ifPk ≤1.
h
Here, |w(cid:5) | is the Euclidean norm of the wind at the element centroid, h is
k k
a measure of the element length in the direction of the wind,2 and Pk :=
h
|w(cid:5) |h /(2(cid:7))istheso-calledelementPecletnumber.
k k
For all convection-diffusion problems solved in IFISS, the discretization is
done using Q approximation on either uniform or stretched grids with sta-
1
bilization parameter (10). There is also scope for a user-defined stabilization
parameter so that the effect of adding streamline diffusion can be studied in
detail. Further information specific to finite element modelling of convection-
diffusionproblemscanbefoundinElmanetal.[2005b],GreshoandSani[1998],
Morton[1996],QuarteroniandValli[1997]andRoosetal.[1996].Theaposteri-
orierrorestimationstrategybuiltintotheIFISSsoftwareisanaturalextension
ofthatusedforthePoissonequation.Giventhediscretesolutionu ,theerror
h
e=u−u ∈ V isestimatedlocallybysolvingthe5×5localPoissonproblem
h (cid:8) (cid:9) (cid:2)
(cid:7)(∇e ,∇v) = R0,v −(cid:7) (cid:9)R ,v(cid:10) ∀v∈Q
T T T T E E T
E∈E(T)
posedovereachelementofthegrid,wheretherighthandsidedataisgivenby
the equidistributed interior edge flux jumps R := 1[[∂uh]], together with the
E 2 ∂n
piecewiseconstantinteriorresiduals
R0 :={f − w(cid:5) ·∇u }| .
T h T
2Onarectanglewithsidesoflengthhx,hy inwhichthewindformsanangleθ =arctan(|wy/wx|)
atthecentroid,hk isgivenbymin(hx/cosθ,hy/sinθ).
ACMTransactionsonMathematicalSoftware,Vol.33,No.2,Article14,Publicationdate:June2007.
Algorithm866:IFISS,AMatlabToolboxforModellingIncompressibleFlow • 7
In this case, the local, elementwise, error estimator is given by η = (cid:12)∇e (cid:12) ,
T T T
andtheglobalerrorestimatorisagaingivenby
(cid:10) (cid:11)
(cid:2) 1/2 (cid:17) (cid:17)
η:= η2 ≈(cid:17)∇(u−u )(cid:17).
T h
T∈Th
This approach was originally suggested by Kay and Silvester [2001] and is a
simplifiedversionoftheimplicitestimatorintroducedandanalyzedinVerfu¨rth
[1998].
2.3 TheStokesEquations
TheStokesequations,givenby
−∇2(cid:5)u+∇p = (cid:5)0, (11)
∇· (cid:5)u = 0, (12)
togetherwithboundaryconditions(analogousto(8))
(cid:5)u= w(cid:5) on∂(cid:2) , ∂(cid:5)u − (cid:5)np= (cid:5)0on∂(cid:2) , (13)
D ∂n N
represent a fundamental model of viscous flow. The variable (cid:5)u is a vector-
valued function representing the velocity of the fluid, and the scalar function
p represents the pressure. The equation (11) represents conservation of the
momentum of the fluid (the momentum equation), and equation (12) enforces
conservation of mass (the incompressibility constraint). The crucial modelling
assumption is that the flow is ‘low speed,’ so that convection effects can be
neglected.
The formal aspects of finite element discretization of the Stokes equations
are the same as for the Poisson equation. However, two different types of test
functionareneededinthiscase,namelyvector-valuedvelocityfunctions (cid:5)v =
h
[v ,v ]T and scalar pressure functions q . By extending this notation to the
x y h
discretecaseinanobviousway,theweakformulationcanbewrittenas
(∇(cid:5)u ,∇(cid:5)v )−(p ,∇· (cid:5)v ) = 0 ∀(cid:5)v ∈X ,
h h h h h h
(q ,∇· (cid:5)u ) = 0 ∀q ∈ M ,
h h h h
where X and M are appropriate finite dimensional subspaces of H1((cid:2)) and
h h
L ((cid:2)),respectively.Thefactthatthevelocityandpressurespacesareapproxi-
2
matedindependentlyleadstothenomenclature,mixedapproximation.Imple-
mentationagainentailsdefiningappropriatebasesforthevelocityandpressure
finiteelementspacesandconstructingtheassociatedfiniteelementcoefficient
matrix (see Elman et al. [2005b, section 5.3]). In matrix form, the resulting
linearsystemis
(cid:4) (cid:5)(cid:4) (cid:5) (cid:4) (cid:5)
A BT u f
= . (14)
B O p g
This is a system of dimension n +n , where n and n are the numbers of
u p u p
velocityandpressurebasisfunctions,respectively.ThematrixAisa2×2block
ACMTransactionsonMathematicalSoftware,Vol.33,No.2,Article14,Publicationdate:June2007.
8 • H.C.Elmanetal.
diagonal matrix with scalar Laplacian matrices defining the diagonal blocks,
andthematrix Bisann ×n rectangularmatrix.
p u
If(14)istoproperlyrepresentacontinuousStokesproblem,thenthefinite
dimensional approximation spaces need to be chosen carefully. The theory of
howtochoosebasisfunctionsappropriatelyisdealtwithbytheinf-supstability
condition. We omit a discussion of this analysis; full details can be found in
Brezzi and Fortin [1991], and an accessible discussion is available in Elman
etal.[2005b,Section5.3.1].Twoinf-supstablemixedmethodsareimplemented
inIFISS:
—Q –Q .Theso-calledTaylor-HoodmethodusesQ approximationforvelocity
2 1 2
andcontinuousQ approximationforpressure.
1
—Q2–P−1. This uses Q2 approximation for velocity together with a discontin-
uous linear pressure. Specifically, the P−1 element has a central node with
three associated degrees of freedom, the pressure at the centroid, and its x
and y derivatives.
It is an unfortunate fact of life that the simplest mixed approximations,
such as the case when velocities and pressures are defined on the same set
of grid points, are not inf-sup stable. In such cases however, it is possible to
use stabilization techniques to circumvent the inf-sup condition. Details can
befoundinElmanetal.[2005b,section5.3.2].Inpractice,stabilisationofthe
lowestordermethodsleadstoaStokessystemoftheform
(cid:4) (cid:5)(cid:4) (cid:5) (cid:4) (cid:5)
A BT u f
= (15)
B −βC p g
forsomen ×n symmetricpositivesemidefinitematrixC andastabilization
p p
parameterβ >0.TwostabilizedelementpairsareimplementedinIFISS:
—Q –P .Thisisthemostfamousexampleofanunstableelementpair,using
1 0
bilinear approximation for velocity and a constant approximation for the
pressure.
—Q –Q .Herethevelocityandpressuredegreesoffreedomaredefinedatthe
1 1
samesetofgridpoints.Thisisveryappealingfromthepointofviewofease
ofprogrammingandcomputationalefficiency.
StokesflowtestproblemscanbesolvedwithIFISSusinganyofthefourdis-
cretizations listed above, with stabilization applied where necessary. An opti-
malvalueofthestabilizationparameterisdeterminedautomatically,although
the parameter value can also be changed by a user who wishes to investigate
theeffectsofstabilizationonsolutionaccuracy.
AnefficientstrategyforaposterioriestimationfortheStokesflowproblem
canberealizedbybuildingonthePoissonerrorestimationstrategyandsolving
ananalogouslocalproblemforeachvelocitycomponent.Notethatthevelocity
approximationtypicallyhasadiscontinuousnormalderivativeacrossinterele-
mentboundaries,soitisnaturaltogeneralizethefluxjump(5)togiveastress
jumpacrossedge E adjoiningelementsT and S:
[[∇(cid:5)u − p (cid:5)I]] :=((∇(cid:5)u − p (cid:5)I)| −(∇(cid:5)u − p (cid:5)I)| )(cid:5)n .
h h h h T h h S E,T
ACMTransactionsonMathematicalSoftware,Vol.33,No.2,Article14,Publicationdate:June2007.
Algorithm866:IFISS,AMatlabToolboxforModellingIncompressibleFlow • 9
(cid:5)
IfaC0pressureapproximationisused,thenthejumpin p Iiszero.Exactlyas
h
for the characterization (6) for the Poisson problem, the mixed approximation
errorfunctions (cid:5)e:= (cid:5)u− (cid:5)u ∈Xand(cid:7) := p− p ∈ M canbeshowntosatisfya
h h
setoflocalizedStokesproblems:
(cid:6) (cid:7)
(cid:2) (cid:18) (cid:19) (cid:2) (cid:2)
(∇(cid:5)e,∇(cid:5)v) −((cid:7),∇· (cid:5)v) = (R(cid:5) , (cid:5)v) − (cid:9)R(cid:5) , (cid:5)v(cid:10)
T T T T E E
T∈Th (cid:2) T(cid:2)∈Th E∈E(T)
− (q,∇· (cid:5)e) = (R ,q)
T T T
T∈Th T∈Th
for all ((cid:5)v,q) ∈ X × M. Here, the equidistributed stress jump is R(cid:5) :=
E
1[[∇(cid:5)u − p (cid:5)I]], and the elementwise interior residuals are given by R(cid:5) :=
2 h h T
{∇2(cid:5)u −∇p }| and R :={∇· (cid:5)u }| ,respectively.
h h T T h T
The IFISS software computes an estimate of the approximation error when-
ever Q velocity approximation is used (that is, if either stabililized Q –P or
1 1 0
stabilizedQ –Q mixedapproximationisused).Usingthespace Q(cid:5) :=(Q )2,
1 1 T T
thevelocityerrorsaregivenbysolvingapairof5×5Poissonproblems.Specif-
ically, (cid:5)e ∈ Q(cid:5) iscomputedsatisfying
T T
(cid:2)
(∇(cid:5)e ,∇(cid:5)v) =(R(cid:5) , (cid:5)v) − (cid:9)R(cid:5) , (cid:5)v(cid:10) ∀(cid:5)v∈ Q(cid:5) , (16)
T T T T E E T
E∈E(T)
for every element in the grid. The local error estimator is given by the combi-
nationofthe‘energynorm’ofthevelocityerrorandtheL normoftheelement
2
divergenceerror,thatis,
η2 :=(cid:12)∇(cid:5)e (cid:12)2 +(cid:12)R (cid:12)2 . (17)
T T T T T
(cid:15) (cid:16)
(cid:3) 1/2
Theglobalerrorestimatorisη:= η2 .
T∈Th T
ThismethodoferrorestimationwasintroducedbyKayandSilvester[1999],
and is a modification of an approach devised by Ainsworth and Oden [1997].
Foranassessmentofitseffectiveness,seeElmanetal.[2005b,section5.4.2].
2.4 TheNavier-StokesEquations
Thesteady-stateNavier-Stokesequationscanbewrittenas
−ν∇2(cid:5)u+ (cid:5)u·∇(cid:5)u+∇p = (cid:5)f, (18)
∇· (cid:5)u = 0,
whereν >0isagivenconstantcalledthekinematicviscosity.Associatedbound-
aryconditionsaregivenby
(cid:5)u= w(cid:5) on∂(cid:2) , ν∂(cid:5)u − (cid:5)np= (cid:5)0on∂(cid:2) . (19)
D ∂n N
This system is the basis for computational modelling of the flow of an incom-
pressible Newtonian fluid such as air or water. The presence of the nonlinear
convection term (cid:5)u·∇(cid:5)u means that boundary value problems associated with
theNavier-Stokesequationscanhavemorethanonesolution.
ACMTransactionsonMathematicalSoftware,Vol.33,No.2,Article14,Publicationdate:June2007.