Introduction
In this workshop we will introduce you to the least-squares spectral element method. As you can see from the lecture notes, this method is a combination of the weak formulation derived from minimizing the residual in the least-squares sense and the approximation by spectral elements.
In setting up this workshop we had to decide whether we wanted to let you become familiar with the method in a complex setting – non-linear, time-dependent calculations – or restrict ourselves to simpler problems, which allows you to get a better understanding of the method. We have chosen for the second line for several reasons: 1) The method is quite new, so it is hard to find references which address the more theoretical issues, and 2) The extension to more complicated problems merely requires time (and a fair amount of book-keeping).
The goal of this workshop is to understand the method by combining theoretical deriva- tions with actual computational results.
This workshop has also been given at the Von Karman Institute for Fluid Dynamics in the summer of 2006. In this workshop students spend 2 whole days on the subject. In order to retain the balance between theory, programming and analyzing the results as much as possible, for some problems the answers are already provided. This allows you to spend the two hours allocated in this Burgers course to play with the sample programs.
The main focus will be on the 1D problem, but the code for the hyperbolic and elliptic 2D problems are also provided for those interested.
If you have any questions during or after this workshop, do not hesitate to contact me (M.I.Gerritsma@TUDelft.nl). We hope you have a nice week and look forward to your appraisal.
The exercises for this part of the Burgers course are scheduled on Friday 25 January at 15.45 – 17.30.
1D LSQSEM problem
We will use the 1D model problem of the lecture notes to get some practical experi- ence with different forms of least-squares spectral element method. Before we start the exercises we will repeat some theory from the lecture serie.
Let us consider a first order scalar ordinary differential equation in the interval [−1, 1]
given by
u0(x) = g(x) , (1)
u(−1) = 0 , (2)
where u0 = du/dx and g is a given function.
Consider the norm for u ku0k2L2(−1,1)=
Z 1
−1
¡u0(x)¢2
dx ≤ ∞ , ∀ u ∈ U =©
ku0kL2(−1,1)| u(−1) = 0ª
. (3) Now let us consider all functions for which
kuk2H1(−1,1):= kuk2L2(−1,1)+ ku0k2L2(−1,1) < ∞ . (4) The set of all these functions obviously forms a linear function space denoted by H1(−1, 1) where the superscript 1 denotes the highest derivative that must be square integrable.
In order to satisfy the boundary conditions we choose the linear subspace H =©
u ∈ H1(−1, 1) | u(−1) = 0ª
. (5)
Now it is clear that
ku0k2L2(−1,1)≤ kuk2H1(−1,1), ∀u ∈ H . (6) So the differential operator is bounded – continuous – in H with constant C2 = 1. In order to show coercivity we use the Poincar´e inequality
u(x) = u(−1) + Z x
−1
u0(ξ) dξ = Z x
−1
u0(ξ) dξ , ∀u ∈ H . (7) Therefore
kuk2L2(−1,1) = Z 1
−1
u2dx = Z 1
−1
½Z x
−1
u0(ξ) dξ
¾2 dx
≤ Z 1
−1
½Z x
−1
u0(ξ)2dξ
¾ dx
≤ Z 1
−1
½Z 1
−1
u0(ξ)2dξ
¾ dx
= Z 1
−1
ku0k2L2(−1,1)dx
= ku0k2
Z 1 dx
So
kuk2H1(−1,1)= kuk2L2(−1,1)+ ku0k2L2(−1,1) ≤ 3ku0k2L2(−1,1). (9) Which gives C1 = 1/3. So we have the desired norm equivalence
1
3kuk2H1(−1,1)≤ ku0kL22(−1,1) ≤ kuk2H1(−1,1). (10) Now that we have established norm-equivalence for all u ∈ H we can insert the function u − uex in which case we obtain
1
3ku − uexk2H1(−1,1) ≤ ku0− (uex)0k2L2(−1,1) = ku0− gk2L2(−1,1)≤ ku − uexk2H1(−1,1). (11) The first inequality in (11) tells us that if we have sequence {un}n≥1 for which
n→∞lim k (un)0− gkL2(−1,1)= 0 (12)
in the L2-norm, then the error
n→∞lim kun− uexkH1(−1,1)= 0 (13)
converges in the H1-norm. The second inequality in (11) tells us that the norm of the residual is less or equal to the H1-norm of the best approximation to uex in the H1-norm.
Problem Minimize the L2-norm of the residual, given by
°°u0− g°
°2
L2 = 1 2
Z 1
−1
(u0− g)2dx (14)
by means of variational analysis
²→0lim d d²
°°
°°d
dx(u + ²v) − g
°°
°°
2 L2
= 0 , ∀v ∈ H (15)
Answer Z 1
−1
u0v0dx = Z 1
−1
v0g dx , ∀v ∈ H
The next step is to discretize this formulation. As in finite element methods the domain is divided in Ncellssubdomains. On each of these subdomains the solutions is approximated by a Lagrange polynomial (see the lecture notes for these polynomials) of order P :
uhe(x) =
P +1X
j=1
ujhj(ξ) , (16)
with uj the unknown coefficients and for ξ = 2x−x∆x1 − 1, locale coordinate system, with x1the coordinate of the left vortex-node of the element and ∆x denotes the length of the element.
If we substitute (16) into the derived minimization. By choosing v = hi(ξ) with 1 ≤ i ≤ (P + 1) the element matrix can be derived.
Both the element matrix and the element load vector have already been implemented in the Matlab file mainExercise1.m, which is located in the directort LSSEM_1D. When composing the entries of the element a number of integrals will have to be evaluated. In finite element methods this can be done exactly, however when a spectral approximation is used in general these integrals have to be evaluated numerically. In standard spectral element methods these integrals are evaluated using a Gauss integration:
Z 1
−1
f (ξ) dξ ≈
N +1X
i=1
wif (ξi) , (17)
with wi the N + 1 Gauss-Lobatto-Legendre weights and ξi the GLL-roots. Often N = P is chosen, however it is possible, and for some cases necessary, the choose N ≥ P (keep in mind that (17) is exact when f(x) is a polynomial up to order 2N − 1). The Matlab routines gl_roots(n) and gl_weights(n) are available to calculate the GLL-roots and weights. If you would need to calculate the value of the basis functions h(ξ) or the derivative of the first derivative in a certain node the functions derivativeMatrices1D.m and derivativeMatricesRandomNodes1D.m will be useful.
Problem Implement the parts of the code where you calculate (numerically) the L2- norm of the error and the L2-norm of the residual
Problem Implement also the H1-norm of the error and verify the a priori estimate (11).
Problem Test your program by inserting a polynomial exact solution. Both the error and the residual should be zero when the polynomial degree in the approximation is cho- sen sufficiently high
A numerical method is called optimal when the highest possible convergence rate is attained. The convergence rate cannot exceed the convergence rate obtained by direct interpolation of the exact solution, otherwise our numerical method would approximate the exact solution faster. Let us therefore establish the convergence rate of interpolation of the exact solution.
Problem Assume that the exact solution is square integrable over the interval [−1, 1].
Then the exact solution can be written as uex(x) =
X∞ i=0
αiLi(x) , x ∈ [−1, 1] ,
where Li(x) is the ith Legendre polynomial and αi is called the Legendre coefficient.
Assume that we want to approximate this function with a finite dimensional representa- tion of the form
uN(x) = XN i=0
˜
αiLi(x) , x ∈ [−1, 1] .
Show that the finite dimensional approximation which minimizes the L2-norm of the error kuex− uNk2L2[−1,1] is given by ˜αi= αi, for i = 0, . . . , N .
Hint: make use of the fact that the Legendre polynomials are orthogonal over the interval
[−1, 1], i.e. Z 1
−1
Li(x)Lj(x) dx = 2 2i + 1δij , where δij is the Kronecker delta defined by
δij =
1 if i = j 0 if i 6= j Answer The error is given by
uex− uN = X∞ i=0
βiLi = XN
i=0
(αi− ˜αi) Li+ X∞ i=N +1
αiLi.
So βi = αi− ˜αi for i = 0, . . . , N and βi = αi for i = N + 1, . . . , ∞ Now taking the L2-norm of this error gives
kuex− uNk2 = Z 1
−1
X∞ i=0
βiLi X∞ j=0
βjLjdx
= X∞
i=0
βi X∞ j=0
βj Z 1
−1
LiLjdx
= X∞
i=0
βi2 2 2i + 1
= XN
i=0
2
2i + 1(αi− ˜αi)2+ X∞ i=N +1
2 2i + 1α2i .
Now the αi’s are given, so the only degrees of freedom we have to minimize the error are the ˜αi’s and we readily see that by taking ˜αi= αithe L2-norm of the error is minimized.
So the error for th optimal L2 approximation is given by kuex− uNk2=
X∞ i=N +1
2 2i + 1α2i .
Now assume that the exact solution is a member of the Sobolev space Hk(−1, 1) (the function is square integrable and all its derivatives up to the k-th derivative are square integrable) then the error can bounded by
kuex− uNk2 = X∞ i=N +1
2 2i + 1
(i − s)!
(i + s)!
(i + s)!
(i − s)!α2i ≤ (N + 1 − s)!
(N + 1 + s)!
X∞ i=N +1
2 2i + 1
(i + s)!
(i − s)!α2i , where s = min(N + 1, k). Now we use the fact that
|uex|2Hs
w(−1,1)= Z 1
−1
¡1 − x2¢2µ dsuex
dxs
¶2 dx =
X∞ i=N +1
2 2i + 1
(i + s)!
(i − s)!α2i . This allows us to bound the interpolation error as
kuex− uNk2 ≤ (N + 1 − s)!
(N + 1 + s)!|uex|2Hs w(−1,1).
If the exact solution is very smooth, i.e. k very large, we have that s = N + 1, so kuex− uNk2 ≤ 1
(2N + 2)!|uex|2HN +1 w (−1,1). Problem Use the error estimate kuex−uNk2L2(−1,1)in terms of |uex|2
HwN +1(−1,1)to derive the error over an arbitrary interval [a, b] and show that for a fixed polynomial degree N , the error converges as
kuex− uNkL2(a,b)= O¡
(b − a)N +1¢ .
Answer Let
ξ = 2x − (a + b)
b − a , =⇒ dξ = 2
b − adx , d
dξ = b − a 2
d dx . So we have that
kuex− uNk2L2(a,b) = Z b
a
¡uex(x) − uN(x)¢2 dx
= Z 1
−1
¡uex(ξ) − uN(ξ)¢2 b − a 2 dξ
= b − a
2 kuex− uNkL2(−1,1). Furthermore we have,
dN +1uex dxN +1 =
µ 2 b − a
¶N +1
dN +1uex dξN +1 , and therefore
|uex|2HN +1 w (a,b)=
µ 2 b − a
¶2N +1
|uex|2HN +1 w (−1,1) . Combining these results gives
kuex− uNk2L2(a,b) = b − a
2 kuex− uNkL2(−1,1)
≤ b − a 2
1
(2N + 2)!|uex|2HN +1 w (−1,1)
= b − a 2
µb − a 2
¶2N +1 1
(2N + 2)!|uex|2HN +1 w (a,b)
Therefore
kuex− uNkL2(a,b)≤
µb − a 2
¶N +1s 1
(2N + 2)!|uex|HN +1
w (a,b)= O¡
(b − a)N +1¢ .
Problem Perform convergence studies to establish the convergence rate of the LSQSEM method that you have implemented. Does it make a difference whether you use polyno- mials of odd degree or even degree. Is LSQSEM optimal for this problem?
Problem Try different right hand side functions and associated exact solutions to see how the method responds to this change. Try for instance, the exact solution
eκx− e−κ
eκ− e−κ , κ > 0 .
What happens if κ is chosen sufficiently large? For instance κ = 100. Perform a conver- gence study in which you play with the number of elements and the polynomial degree.
Do you see wiggles? Do the wiggles pollute the whole domain?
Problem Perform a similar convergence study using the spectral Galerkin method.
How do you need to modify the code to convert it to a Galerkin Method (NB Also change the solver!). Does it make a difference whether you use polynomials of odd degree or even degree. Is the Galerkin method optimal for this problem?
Derivation of the LSQSEM for the Poisson equation
Having analyzed the 1D model problem we now turn to a two-dimensional problem. The Poisson equation we wish to consider in this workshop is given by
∆φ = f , for x ∈ Ω , (18)
with ∂φ
∂n = g , for x ∈ ∂Ω . (19)
Here Ω ∈ R2 is a bounded – not necessarily simply-connected – domain with boundary
∂Ω. n is the outward unit normal to ∂Ω. The functions f and g ∈ L2(Ω) are given. We need to solve this problem for the unknown function φ, which we will refer from now on to as the potential.
Although the Poisson equation plays an important role in many branches of physics, we will mainly interpret this equation as the equation describing an incompressible, irrota- tional, inviscid flow in the absence of body forces. The velocity field generated by such a potential flow is given by
~u = ~∇φ ⇐⇒ u = ∂φ
∂x , v = ∂φ
∂y . (20)
Problem Rewrite the Poisson problem as a first order system of equations and write down the norm of the residuals in the L2 norm. For which functions do these norms make sense?
Answer First order system
∂u
∂x +∂v
∂y = f , u = ∂φ
∂x , v = ∂φ
∂y .
Physically the first equation expresses conservation of mass for an incompressible flow and the remaining equation define the velocity field. The L2-norm of the residuals is given by
k∂u
∂x +∂v
∂y− f k2+ ku −∂φ
∂xk2+ kv − ∂φ
∂yk2 .
These equations only make sense if u ∈ L2, v ∈ L2, ~∇ · ~u ∈ L2 and ~∇φ ∈ L2. Note that we do not require that ∂u/∂x and ∂v/∂y separately belong to L2!!
Note that for an irrotational flow we have that the vorticity ξ = ∂v/∂x − ∂u/∂y , is zero.
Problem Why don’t we add the condition that the flow is irrotational?
Answer The fact that we derive the velocity field from a potential function guarantees that the resulting flow field is irrotational
∂v
∂x−∂u
∂y = ∂
∂x
∂φ
∂y − ∂
∂y
∂φ
∂x ≡ 0 .
Alternatively, we can also employ the conventional Galerkin method to solve the Poisson equation. An advantage of this approach is that we do not need to convert the second order differential equation to an equivalent first order system. Having solved (18) for φ allows one obtain the velocity components in a post-processing step. In general the accuracy of these derived quantities is one order lower than the order of accuracy of the primal variable φ. Additionally, the derived quantities are in general not continuous across element boundaries.
Problem Write down the weak formulation for the Poisson problem using the Galerkin method without introducing auxiliary variables. For what functions φ does this formula- tion makes sense? Argue that continuity of φ across element boundaries is sufficient to have a conforming approximation.
Answer Multiply the Poisson equation with an arbitrary test function ψ and apply integration by parts
0 = Z
Ω
(∆φ − f ) ψ dΩ = Z
∂Ω
∂φ
∂nψ dS − Z
Ω
³∇φ, ~~ ∇ψ
´ dΩ −
Z
Ω
f ψ dΩ .
These integrals make sense if φ, ψ ∈ H1(Ω) and f ∈ L2(Ω). Discontinuities are square in- tegrable, but the derivative of discontinuities – Dirac delta distributions – are not square integrable, therefore we do not allow discontinuities in this formulation. This implies, that for a conforming method, the solution between elements needs to be continuous.
The derived quantities, such as a velocity field, will then in general be discontinuous between elements.
If we first rewrite the Poisson equation as a first order system and then we apply the Galerkin method, we obtain a so-called mixed method. If we follow this approach, we start with the first order system
~u − ~∇φ = 0 , in Ω , (21)
and
∇ · ~u = f , in Ω .~ (22)
Now multiplying (21) by ~v ∈ W and integrating, and multiplying (22) by ψ ∈ H and integrating by parts, we obtain the mixed Galerkin formulation. Here W and H are admissible function spaces to be determined.
(~u, ~v) −
³∇φ, ~v~
´
= 0 , ∀~v ∈ W , (23)
−
³∇ψ, ~u~
´
= (f, ψ) . (24)
For convenience we assume homogeneous boundary conditions. The mixed Galerkin method produces a saddle-point variational problem and therefore existence and unique- ness require that the function spaces W and H satisfy the LBB condition given by
sup
06=φ∈H
R
Ω∇φ · ~u dΩ~ µR
Ω
¯¯
¯~∇φ
¯¯
¯2 dΩ
¶1/2 ≥ γk~ukL2(Ω) , ∀~u ∈ W ,
where γ > 0. This condition puts restrictions on the function spaces that can be used. If incompatible function spaces are used, the method may not produce sensible results. If appropriate function spaces are used, and we assume that (φ, ~u) ∈ Hr+1(Ω) × [Hr(Ω)]d then we have the following error estimate
|φ − φh|1+ k~u − ~uhk0≤ Chr(|φ|r+1+ k~ukr) .
So even for a compatible approximation, the approximation of the velocity field is always one order lower than the approximation of the potential.
Inspecting the mixed Galerkin formulation we also see that the method is non-positive
Problem We have derived the L2-norm of the residual of the first order system. In order to find the solution which minimizes the norm of the residuals, we need to apply variational analysis to obtain the weak least-squares formulation. Assume that the sum of the squared residuals is given by J (φ, u, v), then variational analysis gives the three equations
∂J (φ + ²ψ, u, v)
∂²
¯¯
¯¯
²=0
= 0 , ∀ψ ∈ W ,
∂J (φ, u + ²˜u, v)
∂²
¯¯
¯¯
²=0
= 0 , ∀˜u ∈ U ,
and ∂J (φ, u, v + ²˜v)
∂²
¯¯
¯¯
²=0
= 0 , ∀˜v ∈ V . Explicitly derive the weak least-squares formulation.
Answer Applying the above recipe to the φ gives 2
µ u −∂φ
∂x,∂ψ
∂x
¶ + 2
µ v −∂φ
∂y,∂ψ
∂y
¶
= 0 , ∀ψ ∈ W , and to u
2 µ∂u
∂x+∂v
∂y − f,∂ ˜u
∂x
¶ + 2
µ u −∂φ
∂x, ˜u
¶
= 0 , ∀˜u ∈ U , and finally to v
2 µ∂u
∂x +∂v
∂y− f,∂˜v
∂y
¶ + 2
µ v −∂φ
∂y, ˜v
¶
= 0 , ∀˜v ∈ V ,
Problem Show that the weak system you have derived by requiring that the total residual is minimized, is symmetric
Answer By taking a linear combination of the above equations we have µ
u − ∂φ
∂x, ˜u −∂ψ
∂x
¶ +
µ v −∂φ
∂y, ˜v −∂ψ
∂y
¶ +
µ∂u
∂x +∂v
∂y − f,∂ ˜u
∂x+∂˜v
∂y
¶
= 0 ,
∀ψ ∈ W , ∀˜u ∈ U and ˜v ∈ V .
Here we see that interchanging the trial functions φ, u and v with the test functions ψ, ˜u and ˜v leaves the formulation unaltered, thus showing that the formulation is symmetric.
Having derived a symmetric formulation, we still do not know whether minimization of the residuals actually means that the error with respect to the exact solution also decreases. In order to show this, we need to derive the a priori estimates which relate the L2-norm of the residuals to the error-norm.
Problem Define the following function spaces H1(Ω) =
n
φ ∈ L2(Ω) and ~∇φ ∈£
L2(Ω)¤2o
, (25)
and
H(div, Ω) = n
~u ∈£
L2(Ω)¤2
and ~∇ · ~u ∈ L2(Ω) o
. (26)
The associated norms are defined by kφk2H1(Ω)=
Z
Ω
"
φ2+ µ∂φ
∂x
¶2 +
µ∂φ
∂y
¶2#
dΩ , (27)
and
k~uk2H(div,Ω)= Z
Ω
"
u2+ v2+ µ∂u
∂x+∂v
∂y
¶2#
dΩ . (28)
Show that the sum of the residuals squared is bounded from above by the kφk2H1(Ω)+ k~uk2
H(div,Ω).
Answer The sum of the residuals squared is given by k∂u
∂x +∂v
∂yk2+ ku − ∂φ
∂xk2+ kv − ∂φ
∂yk2 = k~∇ · ~uk2+ kuk2+ kvk2+ k~∇φk2− 2
µ u,∂φ
∂x
¶
− 2 µ
v,∂φ
∂y
¶
Now use the fact 0 ≤ (a + b)2 which implies that −2ab ≤ a2+ b2. This inequality allows us to bound
−2 µ
u,∂φ
∂x
¶
≤ kuk2+ k∂φ
∂xk2 , and
−2 µ
v,∂φ
∂y
¶
≤ kvk2+ k∂φ
∂yk2. Combining gives
k∂u
∂x+ ∂v
∂yk2+ ku −∂φ
∂xk2+ kv − ∂φ
∂yk2 ≤ 2k~uk2+ k~∇ · ~uk2+ 2k~∇φk2. Now using that k~∇φk2≤ kφk2H1(Ω) and k~uk2+ k~∇ · ~uk2 = k~uk2
H(div,Ω)we obtain k∂u
∂x+∂v
∂yk2+ ku − ∂φ
∂xk2+ kv − ∂φ
∂yk2 ≤ 2 h
kφk2H1(Ω)+ k~uk2H(div,Ω)
i . This demonstrates that the residuals are bounded in the norm of unknowns.
This estimate can be interpreted in several ways. First, for all functions φ and ~u which are measurable in the norms given, the residuals are finite as well. If the right hand side
A second implication is that if the error in the associated norms goes to zero, the residuals must go to zero as well.
The last interpretation can be stated as follows: Let W be the space of all functions for which the sum of the residuals yields a finite value and let V be the space of all functions for which kφk2H1(Ω)+ k~uk2
H(div,Ω) is finite. Then the inequality you just derived stated that W ⊂ V.
We now want to derive the reverse, namely that V ⊂ W and that when the norm of the residuals goes to zero, the error goes to zero as well.
Problem The other estimate – coercivity – is usually harder to establish, so let me guide you through the steps involved:
Show that
k~∇ · ~vk2L2+ k~∇ψ − ~vk2L2 ≥ Ckψk2L2 + k~vk2L2 + 2
³∇ · ~v, ψ~
´
. (29)
Hints: Use the following estimates
kf k2 = (f, f ) , where the inner-product of functions (f, g) is defined as
(f, g) = Z
Ω
f (x)g(x) dΩ .
Use integration by parts to derive (29). Also use the Friedrichs inequality which states that on a bounded domain Ω for φ ∈ H1(Ω) with φ = 0 at the boundary ∂Ω we have
kφk2L2(Ω) ≤ Ck~∇φk2L2(Ω). (30) Restrict therefore your attention to homogeneous problems where ψ = 0 at the boundary.
(Is this a restriction?)
Answer Using the Hints we have
k~∇ · ~vk2+ k~v − ~∇ψk2 ≥ k~v − ~∇ψk2
= k~vk2+ k~∇ψk2− 2
³
~v, ~∇ψ
´
= k~vk2+ k~∇ψk2+ 2
³∇ · ~v, ψ~
´
integration by parts
≥ k~vk2+ Ckψk2+ 2
³∇ · ~v, ψ~
´
Friedrichs inequality
Furthermore we have that
k~∇ · ~vk2+ k~u − ~∇ψk2 ≥ k~∇ · ~vk2, (31) and therefore also for C > 0
1 C
h
k~∇ · ~vk2+ k~u − ~∇ψk2 i
≥ 1
Ck~∇ · ~vk2 . (32)
If we combine this estimate with (29) we have that µ
1 + 1 C
¶ h
k~∇ · ~vk2+ k~u − ~∇ψk2 i
≥ 1
Ck~∇ · ~vk2+ k~vk2+ Ckψk2+ 2
³∇ · ~v, ψ~
´
= k 1
√C
∇ · ~v +~ √
Cψk2+ k~vk2
≥ k~vk2 (33)
Problem Every norm should satisfy the triangle inequality ka + bk ≤ kak + kbk .
Use the triangle inequality, (33) and the fact that
k~∇ · ~vk2+ k~u − ~∇ψk2 ≥ k~u − ~∇ψk2, to show that there exists a constant C such that
C h
k~∇ · ~vk2+ k~u − ~∇ψk2 i
≥ k~∇ψk . (34)
Answer
k~∇ · ~vk2+ k~u − ~∇ψk2≥ k~u − ~∇ψk2 =⇒
q
k~∇ · ~vk2+ k~u − ~∇ψk2≥ k~u − ~∇ψk . Furthermore, we have from (33) that
C1 q
k~∇ · ~vk2+ k~u − ~∇ψk2 ≥ k~vk , where C1 =p
1 + 1/C. Therefore we have (1 + C1)
q
k~∇ · ~vk2+ k~u − ~∇ψk2≥ k~u − ~∇ψk + k~vk . Now use the fact that
k~∇ψk = k − ~∇ψk = k~v − ~∇ψ + ~vk ≤ k~v − ~∇ψk + k~vk . Here we use the triangle inequality with a = ~v − ~∇ψ and b = ~v.
Now using the Friedrichs inequality again we have that C2
h
k~∇ · ~vk2+ k~u − ~∇ψk2 i
≥ kψk . (35)
Problem Use (31), (33), (34) and (35) to show coercivity in the spaces ψ ∈ H1(Ω) and
Answer Combining the four aforementioned estimates gives k~∇ · ~vk2+ k~u − ~∇ψk2 ≥ C3
h
k~uk2+ k~∇ · ~vk2+ kψk2+ k~∇ψk2 i
= C3 h
k~vk2H(div,Ω)+ kψk2H1(Ω)
i
. (36)
This coercivity result implies that if the L2-norm of the residual goes to zero, the error in φ converges in the H1-norm and the error in the velocity converges in the H(div)-norm.
The H(div)-norm is a weaker norm than the H1-norm. When we are going to play with this method, you will see what the implications of this weaker norm are.
The coercivity result also states that the space of all functions W for which the residuals are defined contains the space of functions V = H1(Ω) × H(div, Ω). So, we have shown now that the both spaces, W and V contain the same functions, which essentially means that both spaces are equivalent.
Moreover, the two inequalities that we derived can be combined as C1
h
k~vk2H(div,Ω)+ kψk2H1(Ω)
i
≤ k~∇ · ~vk2+ k~u − ~∇ψk2≤ C2 h
k~vk2H(div,Ω)+ kψk2H1(Ω)
i . (37)
Problem Implement the least-squares formulation in the sample program that will be provided. Select a function φ(x, y) and compute the associated right hand side vector f . Also compute – manually – the velocity components u and v.
Check your implementation by choosing for φ a polynomial solution in x and y. This allows you to heck whether for a given polynomial degree of the approximating function space both the residual norm and the error norm become zero.
Then choose a more exotic function φ and perform a similar convergence analysis as was done for the 1D problem, i.e. h-convergence and p-convergence.
Plot kφh− φexkL2, kuh− uexkL2, kvh− vexkL2 and the norm of the residuals as a function of the mesh-size h and as a function of the polynomial degree p.
Answer φ converges optimally with h-refinement, i.e. the convergence rate will be hP +1, if P is the polynomial degree used for the spectral elements. u and v will converge sub-optimal, i.e. the convergence rate will be hα, with α < P + 1.
With p refinement exponential convergence should result in all error norms.
Optimal LSQSEM
In the previous section is what argued that adding the equation for an irrotational flow
∂v
∂x−∂u
∂y = 0 ,
is superfluous. However, numerically it may make a difference.
Problem Add the irrotationality condition to your discrete system in the LSQSEM program. Perform again the convergence studies. How does the addition of this extra equation affect the error norms for φ, u and v?
Answer All variables converge optimally with h-refinement, i.e. the convergence rate will be hP +1, if P is the polynomial degree used for the spectral elements. With p refinement exponential convergence should result in all error norms.
In many practical problems where we encounter Poisson problems we are not so much interested in the potential φ, but more in the ~∇φ. We have seen that the optimal least-squares method converges optimally for these variables of interest. This cannot be achieved with the Galerkin or finite difference/finite volume methods. Furthermore, the velocity field produced by the optimal least-squares method is continuous, so we do not have multi-valued solutions at element boundaries. And finally, we do not need to bother about compatible function space as expressed by the LBB condition. These features make the LSQSEM formulation an interesting alternative to existing discretization methods.
And then ...
You have had some experience with the least-squares spectral element method. The main building block is the norm equivalence between the residual norm and the error in suitable norms. For many equations of interest such a priori estimates can be found in the literature. An important tool in deriving such error estimates is the so-called ADN theory for elliptic systems. However, much needs to be done on the mathematical justification of the method.
If we do not have a priori estimates all we can do is hope for the best. Suppose there exists exotic norms for which norm-equivalence can be established. Then use you favorite least- squares discretization and remember that at the discrete level all norms are equivalent.
So even if you do not have a formulation in the appropriate norms, you will still get reasonable results with your method of choice. However, optimal convergence cannot be guaranteed in that case ...
We hope you enjoyed this workshop. Hopefully, the techniques presented inspire you to use LSQSEM, or at least borrow some of its tricks to improve your work.
Additional information can found in the references at the end of the lecture notes
