Table Of ContentShishkin mesh simulation: A new stabilization
3 technique for convection-diffusion problems and
1
0 complementary experiments
2
n
a Bosco Garc´ıa-Archilla∗
J
7 January 18, 2013
1
]
A
Abstract
N
. We present together in this document the paper B. Garc´ıa-Archilla,
h
Shiskin mesh simulation: A new stabilization technique for convection-
t
a diffusionproblems,Comput. MethodsAppl. Mech. Engrg.,256(2013),1–
m
16,andthemanuscriptB.Garc´ıa-Archilla,Shiskinmeshsimulation: Com-
[ plementaryexperiments,whichcontainssomeoftheexperimentsthat,for
the sake of brevity were not included in the first one. This second paper
2
with complementary experiments is neither self-contained nor intended to
v
bepublishedinanymajorjournal,butisintendedtobeaccessibletoany-
1
onewishingtolearnmoreontheperformanceoftheSMSmethod. Follow-
2
3 ingaprincipleofreproduciblecomputationalresearch,thesourcecodesof
4 theexperimentsinthepresentpapersareavailablefromtheauthoronre-
. questoron”http://personal.us.es/bga/bga files/software bga.html”.
2
Excluded are the codes corresponding to Example 7 which is subject of
1
current resarch
2
1
:
v
i
X
r
a
∗Departamento de Matema´tica Aplicada II, Universidad de Sevilla, Sevilla, Spain. Re-
searchsupportedbySpanishMECundergrantMTM2009-07849([email protected])
1
Shishkin mesh simulation: A new stabilization
technique for convection-diffusion problems
Bosco Garc´ıa-Archilla1
Abstract
Anewstabilizationprocedureispresented. Itisbasedonasimulation
of the interaction between the coarse and fine parts of a Shishkin mesh,
but can be applied on coarse and irregular meshes and on domains with
nontrivial geometries. The technique, which does not require adjusting
any parameter, can be applied to different stabilized and non stabilized
methods. Numerical experiments show it to obtain oscillation-free ap-
proximations onproblems withboundaryand internallayers,on uniform
and nonuniform meshes and on domains with curvedboundaries.
Key words. Convection-dominatedproblems. Stabilized methods.
Finite-element methods. Galerkin method. SUPG method.
1 Introduction
The numerical solution of convection-diffusionproblems when convection dom-
inates is, despite more than 30 years of research, a challenging problem nowa-
days. Standard finite-element or finite-difference methods typically suffer from
unphysicalor spurious oscillations unless meshes are takenso fine that are use-
lessforallpracticalpurposes. Thereasonisthepresenceoflayersorthinregions
where solutions change fast. Modification of standard methods, known as sta-
bilized methods have been proposed in the literature, from upwind methods
35 years ago [45], to strongly-consistent stabilized methods like the stream-
line upwind/Petrov-Galerkin (SUPG) method [7], also known as the stream-
line diffusion finite element method (SDFEM), or the Galerkin least-squares
(GALS) method [24]. More recently, local projection stabilization(LPS) meth-
ods,[4],[6],[18],continuousinteriorpenalty(CIP)methods[8],ordiscontinuous
Galerkin(DG)methods[23],[35]havebeenintroduced,tociteafewofthemany
techniquesproposed(see[36],[38]forasurveyofmethods). Itmustbenoticed,
however,thatcomputationalstudies(seee.g.,[2],[26])findithardtoputapar-
ticular method above the others. It must be also mentioned that most of these
methods depend on at least one parameter about which there is no unanimous
agreement on its optimal choice in practical problems [27].
1Departamento de Matema´tica Aplicada II, Universidad de Sevilla, 41092 Sevilla, Spain.
ResearchsupportedbySpanishMECundergrantMTM2009-07849([email protected])
1
A different approach is to use layer-adapted meshes. Among these we cite
Shishkin meshes (described below) [31], [37], which have received considerable
attention in recent years [14], [15], [16], [28], [32], [33] [39], [43] and [47]. How-
ever, it is generally acknowledged that the main drawback of Shishkin meshes
isthedifficultytodesignthemondomainswithnontrivialgeometries,although
some works overcoming this difficulty can be found in the literature [46], [28].
The method we propose, however, does not suffer from the above indicated
drawbacks: It does not depend on parameters and, although it is based on the
ideaofsimulatingaShishkinmesh,theexperimentswepresentshowitproduces
excellent results on domains with nontrivial geometries.
We consider the problem
ε∆u+b u+cu=f, in Ω, (1)
− ·∇
∂u
u=g , in ∂Ω , =g , in ∂Ω . (2)
1 D 2 N
∂n
Here, Ω is a bounded domain in Rd, d = 1,2,3, its boundary ∂Ω being the
disjointunionofΓ andΓ , bandcaregivenfunctionsandε>0isaconstant
D N
diffusioncoefficient. We assumethatΓ− ∂Ω ,Γ− being the inflowboundary
D
of Ω Rd, i.e., the set of points x ∂Ω s⊂uch that b(x) n(x)<0.
⊂ ∈ ·
It is well-known if ε sup b(x) x Ω ( being the euclidean norm)
≪ {| | | ∈ } |·|
boundarylayersarelikelytodevelopalong∂Ω Γ−,althoughtheyhavedifferent
\
structureonΓ0 = x ∂Ω b(x) n(x)=0 andΓ+ = x ∂Ω b(x) n(x)>0 .
{ ∈ | · } { ∈ | · }
As already mentioned, these boundary layers, when present, are responsible of
the spurious oscillations that pollute the numerical approximations obtained
with standard methods unless extremely fine meshes are used. For uniform
meshes, oscillations typically disappear when the mesh P´eclet number
b h
Pe= k kL∞(Ω)2
2ε
(h being the mesh size) is of the order of 1.
Letusbrieflydescribenowtheideaofthemethodweproposeinthefollowing
simple problem:
L(u) εu′′(x)+b(x)u′(x)+c(x)=f(x), 0<x<1, (3)
≡−
u(0)=u(1)=0. (4)
In (3) we assume that b, c and f are sufficiently smooth functions, and that
0<β < min b(x), 0 min c(x). (5)
x∈[0,1] ≤x∈[0,1]
The standard Galerkin linear finite-element method for (3-4) on a partition or
mesh 0=x <x <...<x =1 of [0,1] obtains a continuous piecewise linear
0 1 J
approximation u (x) to u. As it is customary, h denotes the mesh diameter,
h
h = max h , where h = x x , for j = 1,...,J. The approximation
1≤j≤J j j j j−1
−
can be expressed as U(x)=u ϕ (x)+ +u ϕ (x), where the ϕ (x) are
1 1 J−1 J−1 j
···
2
the basis or hat (piecewise linear) functions taking value 1 at the node x and
j
0 in the rest of the nodes of the partition (thus, U(x ) = u ). The values u ,
j j j
j =1,...,J 1, are obtained by solving the linear system of equations
−
a(u ,ϕ )=(f,ϕ ), i=1,...,J 1, (6)
h i i
−
where, a is the bilinear form associated with (3), which is given by
a(v,w)=ε(v′,w′)+(bv′+c,w),
(, ) being the standard inner product in L2(0,1),
· ·
1
(f,g)= f(x)g(x)dx.
Z
0
The Shishkin mesh with J =2N nodes is composed of two uniform meshes
with N subintervals on each side of the transition point x =1 σ, where
N
−
1 2
σ =min( , εlogN),
2 β
for an adequate constant β, that is, x = j(1 σ)/N, for j = 0,...,N, and
j
−
x = x +jσ/N, for j = 1,...,N. Let us consider the coarse and fine grid
N+j N
parts of the Galerkin approximation given by
U (x)=u ϕ (x)+ +u ϕ (x),
c 1 1 N−1 N−1
···
U (x)=u ϕ (x)+ +u ϕ (x),
f N+1 N+1 2N−1 2N−1
···
so that U +u ϕ +U is the Galerkin approximation on the Shishkin mesh.
c N N f
Since for i = 1,...,J 1, the support of the basis function ϕ is [x ,x ],
i i−1 i+1
−
wehavea(U ,ϕ )=0anda(U ,ϕ )=0,forj =1,...,N 1. Consequently
c N+j f j
−
the system (6) on the Shishkin mesh can be rewritten as
a(U ,ϕ ) =(f,ϕ ), i=1,...,N 2, (7)
c i i
−
a(U ,ϕ ) +u a(ϕ ,ϕ ) =(f,ϕ ), i=N 1, (8)
c i N N i i
−
a(U ,ϕ )+u a(ϕ ,ϕ )+a(U ,ϕ )=(f,ϕ ), (9)
c N N N N f N N
u a(ϕ ,ϕ ) + a(U ,ϕ ) =(f,ϕ ), i=N +1,...,2N 1.
N N i f i i
−
(10)
We notice that were it not for the presence of the u a(ϕ ,ϕ ) in (8), the
N N i
system (7–8) would be the equations
a(U,ϕ )=(f,ϕ ), i=1,...,N 1, (11)
i i
−
of the Galerkin approximation U =U ϕ + +U ϕ for the problem
1 1 N−1 N−1
···
εu′′(x)+b(x)u′(x)+c(x)=f(x), 0<x<1 σ, (12)
− −
u(0)=u(1 σ)=0. (13)
−
3
1
0.8
0.6
0.4
0.2
0
−0.2
−0.4
−0.6
−0.8
−1
0 0.2 0.4 0.6 0.8 1
Figure 1: Galerkin approximation on a uniform mesh with N = 9 (continuous
line) to the solution of (12–13) with ε = 10−8, σ = 4εlog(2 N), b = f = 1,
∗
andc=0. The U partofthe GalerkinapproximationonaShishkinmeshwith
c
J =18(brokenline)forsameεandf. Circlesarethevaluesofthetruesolution
on the nodes.
The Galerkin approximation U for this problem, unless εN > 1/2, is likely to
have spurious oscillations of large amplitude as we show in Fig. 1 for ε=10−8,
σ = 4εlog(J), b(x) = f(x) = 1, c = 0 and N = 9. It is however the presence
of u a(ϕ ,ϕ ) in equation (8) that suppresses the oscillations, as we can see
N N i
inFig.1,wherethecomponentU ofthe GalerkinapproximationonaShishkin
c
gridwithJ =2N =18isalsoshown(discontinuousline)togetherwiththetrue
solution at the nodes of the coarse part of the mesh.
It is remarkable that just by adding the value
α∗ =u a(ϕ ,ϕ ) (14)
N N−1 N
to the last equation of the Galerkin method for (12–13) we get the oscillation-
free approximation U . Obviously, in order to have the value of α∗ we have
c
to solve the whole system (7–10). In the present paper, we introduce a tech-
niquetoapproximateα∗ withouttheneedtocomputethewholeapproximation
on the Shishkin grid. In Fig. 1, the approximation computed with the esti-
mated α∗ is indistinguishable from U . Numerical experiments in the present
c
paper show that, in two-dimensional problems, the oscillation-free approxima-
tion on a coarse mesh can be obtained by this technique at half the computa-
tional cost of a Shishkin grid, and a more substantial gain can be expected in
three-dimensional problems.
Furthermore, this technique can be extended when the grid is no part of
any Shishkin grid, while, at the same time, managing to get rid of the spurious
oscillations. This allows to obtain accurate approximations on domains with
nontrivialboundaries, where Shishkin meshes may be difficult to construct. In
spite of this, we call the new technique Shishkin mesh simulation (SMS), since
it was derived, as described above, in an attempt to simulate Shishkin grids.
We must mention, however, that in the present paper we only consider the
case of dominant convection, both in the analysis and in the numerical experi-
4
ments. Thequestionofhowtomodifythemethod(ifnecessary)whenthemesh
P´eclet number Pe tends to one will be addressed elsewhere.
It is well-known that the Galerkin method is a far from ideal method in
convection-diffusionproblems. Let us also notice that despite the good proper-
ties of stabilized methods developed in recent years, the SUPG method is still
consideredthe standardapproach[27]. For this reason,in the numericalexper-
iments presented below, we compare the new method with the SUPG method.
Therestofthepaperisasfollows. InSection2wedescribetheSMSmethod
in detail. In Section 3 we present a limited analysis of the new technique.
Section 4 contains the numerical experiments and Section 5 the conclusions.
2 The new technique: Shishkin mesh simulation
2.1 The one-dimensional case
Weconsider(3–4)satisfying(5). Givenapartition0=x <x <...<x =1,
0 1 J
of[0,1],wedenotebyX thespaceofcontinuouspiecewiselinearfunctions,and
h
by V the subspace of X of functions taking zero values at x = 0 and x = 1,
h h
so that we can express
X =span ϕ V span ϕ .
h 0 h J
{ }⊕ ⊕ { }
We consider the operator L given by
h
L (v )=bv′ +cv , v V .
h h h h h ∈ h
(See alsoRemark1below). We denote byu V the standardGalerkinlinear
h h
∈
finite-elementapproximation,andbyu˜ the SMSapproximation. This isfound
h
as the solution of the least-squares problem
min L (u˜ ) f , (15)
u˜h∈Vh,α∈Rk h h − kL2(0,xJ−1)
subject to the restriction
a(u˜ ,ϕ )+αϕ (x )=(f,ϕ), ϕ V . (16)
h h h J−1 h h
∈
Observethatsinceϕ (x )=0forϕ notproportionaltoϕ ,andrecalling
h J−1 h J−1
howwedefinedα∗ in(14),equation(16)issimilarto(7–8),andthen, theleast-
squares problem (15) is the way of finding the value α hopefully close to α∗.
Notice alsothatthe restriction(16) is,infact, asetofasmanyindependent
restrictions as interior nodes or nodal basis functions in V . Consequently, in
h
theoptimalityconditions,theremustbeaLagrangemultiplierforeveryinterior
node. We gather all this multipliers in a function z V whose value at every
h h
∈
node is that of the corresponding Lagrange multiplier. Thus, the optimality
conditions of the SMS approximation can be written as the following linear
5
problem: find u˜ ,z V and α R such that
h h h
∈ ∈
L (u˜ ),L (ϕ ) a (ϕ ,z )=(f,L (ϕ )) , ϕ V ,
h h h h L2(0,xJ−1)− h h h h h L2(0,xJ−1) h ∈ h
(cid:0) (cid:1) (17)
z (x )=0, (18)
h J−1
a(u˜ ,ϕ )+αϕ (x )=(f,ϕ), ϕ V , (19)
h h h J−1 h h
∈
where here and in the sequel (, ) denotes the standard inner product
L2(I)
· ·
in L2(I).
Remark 1 Observe that for linear elements the operator L coincides with
h
J
L (v )= ( εv′′+bv′ +cv ) , (20)
h h − h h h |(xj−1,xj)
Xj=1
that is, L applied element by element. This expression of L is better suited
h
to the SMS method for higher-orderelements, a topic that will be studied else-
where.
Remark 2 The fact that the approximation u˜ is found by solving the least-
h
squares problem (15) may suggest a possible relation with the GALS method
(compare for example least-squares problem (15) with that in [36, p. 327]).
However, they are very different methods, since for the examples in Section 4
in the present paper, SUPG and GALS methods are identical (see e.g. [24]),
and, as shown in Section 4, the SMS method and the SUPG method produce
markedly different results.
2.2 The multidimensional case
We will assume that b(x) = 0 for x Ω, and that every characteristic (i.e.,
6 ∈
solution of dx/dt =b(x)) in Ω enters and leaves Ω in finite time. Also for sim-
plicity,wewillassumethatthedomainΩin(1–2)hasapolygonalorpolyhedral
boundary.
Let a triangulation of it, that is, a partition of the closure Ω of Ω in
h
T
n-simpliceswith only faces andverticesin common. Foreveryτ , let (τ)
h
∈T N
be the set of its vertices, and let = (τ) be the set of vertices of .
Nh ∪τ∈ThN Th
Similarly to the one-dimensionalcase, let X the space of continuous piece-
h
wise linear polynomials. We express
X =X− V X+,
h h ⊕ h⊕ h
where ϕ = 0 on Γ if ϕ V , and for ϕ in X− (resp. X+), if ϕ (x) = 0
h D h ∈ h h h h h 6
for x , then x Γ− (resp. Γ Γ−). In the standard Galerkin method,
h D
first an∈eNlement uD ∈X− X+ is se\lected such that the restriction uD to
h ∈ h ⊕ h h|ΓD
Γ is a good approximation to the Dirichlet data g in (2). This restriction is
D 1
typically the interpolant or the L2(Γ )-projection onto the restrictions to Γ
D D
of functions in X . Then, the Galerkin approximation u uD+V satisfies
h h ∈ h h
a(u ,ϕ )=(f,ϕ )+ε g ,ϕ , ϕ V , (21)
h h h h 2 hiΓN ∀ h ∈ h
6
where here and in the sequel, , denotes the L2 inner product on Γ ∂Ω,
Γ
h· ·i ⊆
and
a(v,w)=ε( v, w)+(b v+cv,w), v,w H1(Ω).
∇ ∇ ·∇ ∈
For the new method, similarly to the one-dimensional case, we consider
L (v )=b v +cv .
h h h h
·∇
In order to describe the new approximation, we must set up the multidimen-
sional version of the last interval (x ,x ) in the one-dimensional case. For
J−1 J
this purpose we denote
Γ0+ =(Γ+ Γ0) Γ .
D ∪ ∩ D
WeconsiderasuitableΩ+withΓ0+ ∂Ω+(tobespecifiedbelow)andletusde-
noteby andN thesehtofvertDices⊂in∂Ωh+ Γ andtheirindices,respectively,
Nδ δ h\ D
that is
= ∂Ω+ Γ ), N = j N x ,
Nδ N h\ D δ { ∈ | j ∈Nδ}
andbyRNδ wewillrefer(cid:0)tothesetofrealvectorsoftheform(tj)j∈Nδ. Thepoints
in and the set Ω+ will play a role similar to that of x and (x ,x ) in
Nδ h J−1 J−1 J
the one-dimensional case. This can be seen in Fig. 2, where, shadowed in grey,
we show the set Ω+ for a triangulation of the unit square for b = [1,1]T, with
h
Γ = ∂Ω, so that Γ0 Γ+ is consists of the sides y = 1 and x = 1. We show
D
∪
the points in marked with circles.
δ
N
Figure2: AtriangulationoftheunitsquarewithΓ0 Γ+ markedwithathicker
line, the set Ω+ in grey and the the points of wi∪th circles.
h Nδ
The approximationu˜ uD+V , is then foundby solvingthe least-squares
h ∈ h h
problem
min L (u˜ ) f , (22)
u˜h∈uDh+Vh,t∈RNδ(cid:13) h h − (cid:13)L2(Ω\Ω+h)
(cid:13) (cid:13)
subject to the restriction
a(u˜ ,ϕ )+ t ϕ (x )=(f,ϕ )+ε g ,ϕ , ϕ V . (23)
h h j h j h h 2 hiΓN h ∈ h
jX∈Nδ
7
Theoptimalityconditionsofthis problemcanbe writtenasthe followinglinear
problem: find u˜h ∈uDh +Vh, zh ∈Vh and t∈RNδ such that
L (u˜ ),L (ϕ ) a (ϕ ,z )=(L (ϕ ),f) , ϕ V , (24)
h h h h L2(Ω\Ω+h)− h h h h h L2(Ω\Ω+h) h ∈ h
(cid:0) (cid:1)
z (x )=0, j N , (25)
h j δ
∈
a(u˜ ,ϕ )+ t ϕ (x )=(f,ϕ )+ε g ,ϕ , ϕ V ,
h h j h j h h 2 hiΓN h ∈ h
jX∈Nδ
(26)
where, as in the one-dimensional case, z is the Lagrange multiplier of restric-
h
tion (23).
LetusspecifynowthesetΩ+. TheobviouschoiceofsettingΩ+ =B where
h h h
B = τ,
h
τ∩(Γ[0+)6=∅
D
as in Fig. 2, may lead to an unstable method if there are nodes x interior
h
∈N
toB ,asitcanbeeasilycheckedinanumericalexperiment. Tounderstandwhy
h
Figure3: AtriangulationoftheunitsquarewithΓ0 Γ+ markedwithathicker
line and the points of with circles. In grey, the∪sets B (left plot) and Ω+
Nδ h h
(right plot) for b=[1,1]T.
insuchacasetheSMSmethodmaybeunstable,considerthelimitcaseε=0,b
constantandc=0(sothatthebilinearformaisskew-symmetric)andconsider
a mesh as that depicted in Fig. 3, where there is one node x interior to B
i h
who has only one neighbour node x . Then, on the one hand, the basis
k δ
∈ N
function ϕ of the node x satisfies that a (ϕ ,ϕ ) = 0 for all ϕ V , except
i i h i h h h
∈
when ϕ is a multiple of the basis function ϕ of the node x . Furthermore,
h k k
since the support of ϕ is contained in B , we have that L (ϕ ) =0.
i h k h i kL2(Ω\Bh)
Consequently, taking v˜ = ϕ , z = 0, t = 0 for l = k and t = a(ϕ ,ϕ ) we
h i h l k i k
6 −
have a nontrivial solution of
L (v˜ ),L (ϕ ) a (ϕ ,z )=0, ϕ V , (27)
(cid:0) h h h h (cid:1)L2(Ω\Ω+h)− h h h h ∈ h
z (x )=0, j N , (28)
h j δ
∈
a(v˜ ,ϕ )+ t ϕ (x )=0, ϕ V , (29)
h h j h j h h
∈
jX∈Nδ
8
if we set Ω+ = B . It can be easily checked that this is also the case if the
h h
nodex interiortoB isconnectedtomorethatonenodein . Forε>0,one
i h δ
N
canexpectauniquesolutionbut,aswehavecheckedinpractice,itisasolution
where u˜(x )=O(1/ε).
i
ToavoidthissituationweremovefromB theupwindtriangleofanyinterior
h
point. Moreprecisely,foreverynodex ,letusdenotebyτ (x )itsupwind
i h − i
∈N
triangle, that is
τ (x )= τ x λb τ = , λ 0 .
− i h i
{ ∈T |{ − }∩ 6 ∅ → }
Then we define
Ω+ =B τ (x ) .
h h\ − i
(cid:0) [ ◦ (cid:1)
b
xi∈Nh∩Bh
Ifthere aretwo upwindtrianglesbecause x λb isandedgefor λsmallerthan
j
−
some λ , we may select the first one of them after ordering all the elements.
0
With this choice of Ω+, as we will show in Section 3.2.1, the few cases were
h
therearenontrivialsolutionsof(27–29)haveaneasyfixfromthecomputational
point of view. We also note that a slightly different choice of Ω+ is taken
h
in Example 7 in Section 4, where the vector field b has closed integral curves.
2.3 Further Extensions
In the same way that the Shishkin meshes are not restricted to the standard
Galerkin method, the new technique, originally motivated by Shishkin meshes,
is not restricted to the Galerkin method either. In this section we comment on
how to use it with some other methods of widespread use, such as the SUPG,
GALS, LPS and CIP methods.
Inthesemethodstheapproximationu uD+V ,insteadofsatisfying(21),
h ∈ h h
satisfies
a (u ,ϕ )=(f,ϕ ) +ε g ,ϕ , ϕ V , (30)
h h h h h h 2 hiΓN ∀ h ∈ h
where a and (, ) are mesh-dependent bilinear forms. For the SUPG method
h h
· ·
they are given by
a (v ,ϕ )=a(v ,ϕ )+ δ (L (v ),b ϕ ) , (31)
h h h h h τ |τ h · h τ
τX∈Th
(f,ϕ ) =(f,ϕ )+ δ (f,b ϕ ) , (32)
h h h τ h τ
·
τX∈Th
whereδ isandadequatelychosenparameterandL denotesLrestrictedtoτ.
τ |τ
In the GALS method the term b ϕ is replaced by L (ϕ ) (see e. g., [36]
·∇ h |τ h
for a full description of these methods). The extension of the SMS technique
to these methods is as follows: we solve the least-squares problem (22) but we
replace the restriction (23) by
a (u˜ ,ϕ )+ t ϕ (x )=(f,ϕ ) +ε g ,ϕ . ϕ V . (33)
h h h j h j h h h 2 hiΓN h ∈ h
jX∈Nδ
9
