2008
Discontinuous RMMC Galerkin methods Lecture 8
1 0.75 0.5 0.25 0 1 0.75 0.5 0.25 0
Jan S Hesthaven Brown University [e-mail] about a vertical cylinder in a finite8.4 Scattering
y
-0.25 -0.5 -0.75
8.4 Scattering about a vertical cylinder in-0.75finite-width channel a
-0.5 0
1.203 1.142 1.081 1.019 0.958 0.896 0.835 0.774 0.712 0.651 0.590 0.528 0.467 0.405 0.344 0.283 0.221 0.160 0.099 0.037
-0.25 -0.5
-0.0028 -0.0072 -0.0117 -0.0162 -0.0207 -0.0252 -0.0297 -0.0342 -0.0387 -0.0432 -0.0477 -0.0522 -0.0567 -0.0612 -0.0657 -0.0702 -0.0747 -0.0791 -0.0836 -0.0881
y
147
-1 -1
x
0.5
1
-1 -1
-0.5
0
x
0.5
1
width channel
1 0.75 0.5 0.25 0 1
-0.25 -0.5 -0.75 -1 -1
A numerical study is carried 10out to both 0.75 9 test the consistency of the imposed bound8 0.5 ary conditions and investigate 7how the geo6 metric representation of the the domain may 0.25 affect the computed solution. In a numeri0 cal model based on the linear Pad (2,2) roe tational velocity version we set up a test -0.25 for open-channel flow. We seek to model the scat-0.5 tering of an incident wave field which is propagating toward a bottom-mounted rigid cylin-0.75 der positioned in the middle of a finite-width -1 -0.5 0.5 1 -1 channel. 0
0.02
0.02
0.01
0.01
1.95E-06 1.70E-06 1.45E-06 1.20E-06 9.45E-07 6.94E-07 4.43E-07 1.92E-07 -5.90E-08 -3.10E-07 -5.61E-07 -8.12E-07 -1.06E-06 -1.31E-06 -1.57E-06
y
0
y
0
y
-0.01
-0.02
-0.03
-0.03
-0.02
-0.01
0
x
0.01
0.02
0.03
y
-0.01 -0.02 -0.03 -0.03 -0.02 -0.01 0
Figure 8.12: Scattering of waves about a x Darmstadt International Workshop, October 2004 p.79 cylinder in a finite-width channel. Due to the symmetry, the solution is mathex
-0.5
0
0.5
x
0.01
0.02
0.03
1
A brief overview of what's to come
Lecture 1: Introduction, Motivation, History Lecture 2: Basic elements of DG-FEM Lecture 3: Linear systems and some theory Lecture 4: A bit more theory and discrete stability Lecture 5: Attention to implementations Lecture 6: Nonlinear problems and properties Lecture 7: Problems with discontinuities and shocks Lecture 8: Higher order/Global problems
Lecture 8
Let's summarize what we know so far The heat equation Higher order and mixed problems Elliptic problems Basic theoretical results What we know and where to go
Lets summarize
We have a thorough understanding of 1D 1st order problems
For the linear problem, the error analysis and
convergence theory is essentially complete. The theoretical support for DG for conservation laws is very solid. It seems advantageous to consider a nodal approach in combination with dissipation. Limiting is perhaps the most pressing open problem .... and we have a nice and flexible way to implement it all Time to move beyond the 1st order problem
Example 7.1. Consider the see shortly, the generalization to his chapter and, as we will linear heat equation Example 7.1. Consider the linear heat equation -order spatial operators is less direct than onewould expect. 2 2 withu periodic boundary conditions and u taken u 2 [288]. u , x [0, 2], the following simple example, u from = 2], , ux [0, = 2 t found 2 easily = [0, as t t x x2 , x x2],
The heat equation
u(x, t) = with periodic boundary conditions= sin(x). The exact solution is with periodic boundary conditions and u(x, 0) 0) = sin(x). 0) =exactsolution is e nsider theperiodicheat equation linear boundary conditions and u(x, and u(x, The sin(x). The exa with easilyeasily foundfound as found as as easily 2 u u on the u(x, Basedu(x, t) = previous discussions of t = , x [0, 2],t) = e-t -t sin(x). e-t sin(x). u(x, t) = e sin(x). t x2 is tempting to simply write the heat e Based on Based on thediscussions discussions of the discontinuous Galerkin the the Based on previous previous of the discontinuous Galerkin methods, it undary conditions previous discussions of The discontinuous Galerkin methods, and u(x, 0) = sin(x). the exact solution is We tempting write the heat equation as equation as is tempting be temptedthe write this as is tempting can to simply simplytoheat equation as u is to simply to write write the heat -
Let us consider the heat equation
t x u u(x, t) = e-t sin(x). u - u = 0, x - ux = 0, t x x u - ux = 0, t and then identify u as the flux in t evious discussions of the discontinuous Galerkin methods,xit t x u as and and then identify asx the the flux in the first order equation. The resulting u equation flux in the first order mplythen identify just use our scheme becomesequation. The resulting write the heat x as standard and then identify ux as the flux approach order equation. T scheme becomes in the first and then scheme becomes u scheme becomes kduk duk k = v k k ^ ukk , v kk h D v Sv h v kt - rx , kMkh h - k k = - h = n rk hh- M k (x) dx,- Sv h = = D uk ux du 0, h h dt Sv h = - - v k = D r uk , M (x) dx, k dt D^ vh - v h h duh Dk n dt k vnatureuof thekequation. The central flux vh - v k (x = Dr k ,order problem,ka=resulting n k as the the h y ux in each element, k. A first M choice - Sv h flux could k ^a simple central h Givenflux in the reasonable dt for the - D be in each the flux could be a simple central in each element, k. A reasonablethe heatfor element, no preferred direction of flux (i.e., reasonable v = {{vh }}) since choice equation has k. A reasonable choic seems= {{v }}) since the heat equationhas {{v }}) since the heat flux (i.e., in each element, k. A reasonable choice no preferred direction of si vk flux (i.e., v = for h flux could be a propagation. h the duh propagation. - Sv k = - ^ n v k - v k (x) dx, uk , Mk
244 7 Higher-order equations The heat equation
Lets see what happens when we run it
244 7 Table 7.1. Global L2 -errors for solving the heat equation using K elements, each N \K 20 with a local order of approximation, N , using10 scheme in Example 7.1. A 40 the `' marks that the algorithm is unstable. N \K 1 2 4 8 10 4.27E-1 5.00E-1 1.68E-1 7.46E-3
Table 7.1. Global L2 -errors for solving the heat equation using K elements, each with a local order of approximation, N , using the scheme in Example 7.1. A `' Higher-order equations marks that the algorithm is unstable. 80 4.38E-1 4.43E-1 1.26E-1 160 4.39E-1 4.42E-1
1 2 20 4.34E-1 4 4.58E-1 1.37E-1 8
8.60E-3
4.27E-1 4.34E-1 40 80 160 5.00E-1 4.58E-1 4.37E-1 4.38E-1 4.39E-1 1.68E-1 1.37E-1 4.46E-1 4.43E-1 4.42E-1 1.28E-1 1.26E-1 7.46E-3 8.60E-3
4.37E-1 4.46E-1 1.28E-1
0.8 0.6 0.4 0.2 u(x,0.8) 0 N=4 Exact N=2
0.8 0.6 0.4 0.2
N=1
N=1 N=4 Exact N=2
It does not work! It is weakly unstable
6.28
u(x,0.8)
1.57
-0.2 -0.4 -0.6 -0.8 0
0 -0.2 -0.4
3.14 x 4.71
addressing this. 7.1.1 The heat equation
guarantee well-posedness.can rewrite this problem as Since a(x)>0 we We write this as u u = aq, q = a , t x x u u = aq, q = a , t x x to recover a system of first order equations. This we can discretize using the techniques system of firstthe conservation laws; that canwe assume using(u, q) recover a developed for order equations. This we is, discretize that the Now follow our standard approach 246 7 Higher-order can be approximated asthe equations chniques developed for conservation laws; that is, we assume that (u, q)
q(x, t) uh (x, t) = qh (x, t)
K k=1
e idea introduced above. We consider u u = a(x) , t x x u u = a(x) , t x x with uWe know outer DG is good is easily shown that a(x) > 0 suffices = 0 at the that boundaries. It for 1st order systems. to guarantee the outer boundaries. It is easily shown that a(x) > 0 suffices th u = 0 at well-posedness. We rewrite this problem as
1.1 The heat equation Let us first return to the heat equation discussed in Example 7.1 and follow the idea introducednewheat equation discussed in Example 7.1 and follow t us We return a above. We -- consider first need to the idea consider
The heat equation
n be approximated as u(x, t)
uk (x, t) h = k qh (x, t)
K
Np
k=1 i=1
uk (xi , t) h k qh (xi , t)
k i (x),
where, as usual, we represent (u, q) by N -th-order piecewise polynomials on K elements. We recover the strong form
The heat equation
246
Treating this as a 1st order system we have
or the corresponding weak form
Here
u(x, t) uh (x, t) uk (x, t) uk (xi , t) k h h (x), = =N k k qh (x, t) K q(x, t) Higher-order equations qh (x, t) K p qh (xi , t) i 246 u(x, t) 7 uh (x, t) uk uk k=1 (x, t) k=1 i=1 (xi , t) k h h = = k k i (x), q (x, represent (u,hK by N -th-order K qh (xi , t) polynomials on q(x, t) usual,hwe t) q q) t) (x, Np where, as k=1 k=1 i=1 piecewise k u(x, t) uh (x, t) uh (x, t) uk (xi , t) k h (x), = = K elements. We recover the strong form k k where, as usual, we represent (u, q) by N -th-order piecewise polynomials ion qh (x, t) q(x, t) qh (x, t) qh (xi , t) k=1 k=1 i=1 k k k K elements. We duk recover the strong form k h ~ aqk - ^ n ( N h ) - ( aqh ) M usual, we h (x) dx, where, as dt = S represent (u, q) by aq-th-order piecewise polynomials on k k KMk duh = We a q k - theD ( form - ( aq k ) k (x) dx, elements. S recover strong aq k ) ~ ^ n h h h k k k dt k k k a k Dk ^ n M qk duh ua-k au - ( au) k (x) dx, h =S h ~h ^ h aq k M a = S q h - k n k ( h )k-( k h ) k (x) dx, aq D k k k dt uk - ^ n D auh - ( auh ) M qh = S (x) dx, k and the weak form h D k k k k = S a uk - ( kauk ) k (x) dx, and the weak formk duh ^ n M q auh - k ^ n ( aqh )h k (x) dx, (7.1) M k h = -(S ha )T q k + D h k dt du D k and the k h form ^ n ( aqh ) k (x) dx, (7.1) M weak = -(S a )T q k + h k dt k k k k D ^ ( q h = -( ~ a ) a k q MMk duh S-(ST u)T+k + n n auh ) k )(x) (x) dx, k dx. h = (7.1) h Dk ^ k ( aqh k k k dt~ a )T uk + ^ ( n D auh ) (x) dx. M q h = -(S h k As in Section 5.2, we have introduced the two special operators D k k k k ~ a )Tthe + special(operators dx. k ^ As in Section 5.2, weM q hintroduced uh two k n auh ) (x) have = -(S d a(x) k (x) d j (x) j a D k k ~ a= dx, S = dx. a(x) i (x) k i (x) Sij d k (x)dx k k As in Section d a(x)have introduced ij two special operators 5.2, we j (x) dx j a a the k D (x) Da(x) k (x) ~ Sij = dx, Sij = dx. i i k k dx dx D D We note that these operators arek (x) d a(x) j closely connected as d k (x) j a a k k ~ = Sij these operators areclosely connected= dx, Sij as dx. (x) a(x) i (x) i We note that k k xr dx dx a a D S D ~ a(x) i (x) j (x) . ij + Sji = xr a operators are closely connectedxas a l ~ We note that Sij + Sji = these a(x) i (x) j (x) . xl x Before defining the numerical flux, it is worth making a rfew observations. In a ~ +S a = Sij a(x) i (x) (x) . ji Before defining the numerical flux, it is worth making ja few observations. In
7 Higher-order equations
K
K
Np
~ S +S = The heat equation(x) d a(x)
246 a ~
ij
As in Section 5.2, we a introduced the two special operators have a xr d k (x) j Sij defining the numerical flux, it is worth= dx, Sij making a few k (x) dx. = (x) a(x) i i Before observations. In k k dx K dx N How do D choose the fluxes? K D p we
j 7 Higher-order equations k a
xl ji k
We note that these operators are closely Dk connected as a(x) i (x) j (x)
.
general, it is t) u(x, reasonable that the numerical fluxes can have idependencies as uh (x, t) uk (x, t) uk (x , t) k h h = = We note q(x, t) these qoperators are closely connected ask (x , t) i (x), that k i + h (x, t) - q (x, t)+ - q i=1 k=1, ( haqh ) , ( k=1 h ) , ( h auh ) ), ( aqh ) = (( aq) f au h xr a a ~ S S where, as(usual, )we=ij + aqji)-= aqa(x) (iau j )- , (au )+ ). represent (u, q) by N)+ , (x) piecewise. polynomials on -th-order (x) auh g(( ,( h h h h xl
K elements. We recover the strong form
The problem with numerical form is that the two first-order equations are In generic flux, it is worth making a few observations. Before definingduk this the k k k k h ~ athat the numerical )fluxes can k (x) dependencies as tightly coupled through the numerical flux and, hence,) have dx, ^ n ( aqh - ( aqh must be solved simul=S M general, it is reasonable q h - k dt taneously as a globally coupled D If, however, we restrict the generality system. - + - ( aqh ) as a (( aq , aqh k , au ( auh + of the numericalqflux = f uk - h ) n( au) - ( aukh ) , k (x) dx, ) ), ^ Mk k = S ( ) h h h ++ h- - + + k D- ( auh ) = f (( aqhh ) , , ( aqh ) ) , (, ( auh )h ), ( , (auauh ), ). au ( aq h ) = g(( aq )- ( aqh h) )
Problem: ^ n = the ) couples ( loss (x) locality (7.1) M tightly that q k canEverything q h + aflux and, hence, Hence,be solved simulcoupled through -(S numerical local-- aqh ) of dx, the auxiliary must dtrecovered through Dk we see be operation. h taneouslyq(x,a globally coupled system. If, however, we element to compute function, as t), is k truly local variable, used only on each restrict the generality a k ~ a u + However, if we qrestrict)Tit kasconditions.auk ) (x) dx. ^ M h= of the numerical flux as -(S h h the derivatives and impose boundary Dk n ( - Keeping in aq )we have aq - , (aqtwoof (the operators ) ), As in Sectionmind the finherent)properties +specialheat ,equation+(i.e., there 5.2, = (( introduced the h ) , auh ) ( auh ( h h k - is no preferred direction of propagation), it is natural to consider the simple + d a(x) j (x) d k (x) j a ( auh ) = g(( auh ) , ( auh ) ), a k k central ~ flux
k h a T k k k
- ( auh )this generic form auh )+ ), = The problem withduk g(( auh ) , ( is that the two first-order equations are
and the weak form
we cannote that these operators are closely connected as eliminate q-variable locally h We
k
function, q(x,stabilitytruly local variable, used onlyin the following. to compute Semidiscrete t), is a of this scheme is established on each element
(x) dx dx, Sij = k a(x) i (x) dx dx. i k D ( aqh ) = {{ aqh }}, ( auh ) D= {{ auh }}. we see that q can be recovered through a local operation. Hence, the auxiliary
Sij =
+(uMq h h , h ) = h M auh ) auh S - ( h ) h ) ( aq ) - ( aqh ) , h h , q ; - S uh + u( - ( ) note au q + Bh xl where Vh is h h space of)N -th-order polynomials with ), hif xl hhqh ) the element,numerical schem the (h auh = g(( dt hWe ,first auh that (u , on satisfies the support D. If out now k superscript for the localas T S a uh + simplicity. - ( auh ) xr , we the choose the test+ T Mq h - element for h ( auh ) functions h have left see that q k can be recovered through a local operation. Hence,0, ( , ) V , h x we Bh (uh , qh ; h , h ) = l the auxiliary h h h h ) satisfies a truly local scheme, then ote that if (uh , qhq(x, t), is the numerical variable, used only on each element to compute function, we have left out the = superscript qfor the local element for simplicity. h k u h , h = h , where
The heat equation
h 1 d h 2h ~ ~ 2 uT S a q h= 0,T S a uh = uT S a q h + uT (S a )T q h = uof D + h h Dh+ h = qhl + q h in the h q = u , is established h this scheme h - , r Semidiscrete stability following. h 2 dt at we recover where and use that 1 d 2 xr a T a T = a au qT - ( a aq ) u - ( au ) q . T u 2 ~ q h + h (S ) q h =h auh q ,h2 dt h D + qh D + r - l = 0, h h h q h + q h S uh = uh S u hh xT xr T a T a T a T a h
h the h (uh ,first h , and impose),boundary conditions.the test functions as derivatives h ) = if (u(h h ) VIf ,wenumerical scheme, then We q ; note that 0, h , q satisfies the B D. h now choose h andGiven hthe mind the of the heat-equation, a natural flux there use that in nature inherent properties of the heat equation (i.e., Keeping Bh (uh , qh ; h , h ) = 0, ( , h ) Vh , s the space T preferred T polynomials propagation), h Tiselement, toh consider = qh , simple with the x iscoulda q + central=fluxes support on it q natural =ruh , h the noof~N -th-orderS a u direction ofT S a q + uT (S a ) auh qh x , uh S be q h uh ~ = h h h ow choose the testhfunctions ashof N -th-order polynomials with support on thel element, where Vh is the space central flux and use that D. If we now choose the test{{ aq }}, ( au ) = {{ au }}. functions as we recover ( == h h h = u , aqh ) q ,
where V is the space of N -th-order polynomials with supp
But is it stable ?
~ uh S qh + qh S
~ auh qh x , uh = uh S q h + uh (S ) l q h = l where the central fluxes First, assume that a(x) is continuous and use = auh qh - ( aqh ) uh - ( auh ) qh d 1we recover uh 2 ( +aqh )2 =1 da{{qh2 l =(0, auh ) = a{{uh }}. qh D + r - }}, First, assume that= 0, is continuous and use the central 2 D uh D + qh D + r - l a(x) 2 dt 2 dt Consider r , where we have the term x where ( aqh ) = a{{qh }}, ( auh ) = a{{uh = auh qh - ( aq) uau( qhauh )aqh ) uh - (auh ) qh . h = h - - ( qh . h a - + Consider+ r , x uh qh fluxes - where we fluxes me that a(x) is continuous and = is continuous anduusehthe central have the term First, assume that r use the central + h q . a(x) - 2 a - + - r u h qh + u + qh . a{{q elements, we immediately recover au ) a{{qh }}, }}. ( aqh ) =over allh }},( (aqh ) h= = a{{uh( auh ) = a{{uh }}. = - h Summering 2 Consider x , where we have the term Summering over all elements, we immediately recover r , where we have the rterm 1 d uh 2 qh 2 = 0, + ,h ,h 2 dt 1 d a - + r+ - a u- q + + u+ q - . =- uh 2 + qh 2 = 0, h h h h
Computing the local energy in a single element yields
Stability
The heat equation
% compute right hand sides of the semi-discrete PDE rhsu = rx.*(Dr*q) - LIFT*(Fscale.*(nx.*dq)); return
So this show the routine to compute the right-hand sid In HeatCRHS1D.m, we is stable! for the semidiscrete discretization of the simplest problem with a(x) = 1. W Howimpose homogenous Dirichlet boundary conditions by defining the exterio about boundary conditions + ghost states (u+ , qh ), h
u+ h
In a similar fashion, Neumann conditions are imposed as
7.1 Higher-order = 0, [[uh ]] = 2n- u- 2 {{uh }} time-dependent^problems 249 Higher-order time-dependent problems - h = qh - Dirichlet {{qh }} = qh , [[qh ]] = 0. In a similar fashion, Neumann conditions are imposed as
-u- , h
7.1 + qh =
Neumann +
If the boundary conditions are inhomogeneous, this is straightforwardly mo If Inhomogeneous BC inhomogeneous, this is straightforwardly modthe boundary conditions are ified; that is, + - ified; that is, u+ = -u- + 2f (t), qh = qh , h h for a Dirichlet boundary condition u(x, t) = f (t).
uh = uh , qh = -qh
+ - u+ = +- , qh - -qh uh = h -
{{uh }} - u- , [[uh ]] = 0 = h }} = uh , [[uh ]] = 0 {{u {{q }} = 0,h [[q ]] = 2n- q - . h h - ^- h
{{qh }} = 0,
^ [[qh ]] = 2n qh .
for a Dirichlet boundary condition u(x, t) = f (t).
+ - u+ = -u- + 2f (t), qh = qh , h h
... and likewise for Neumann Heat1D.m
function [u,time] = Heat1D(u,FinalTime) Heat1D.m
The heat equation
Back to the example Looks good || u - uh ||
100 N=1 10-2 10-4 10-6 N=4 10-8 10-10 10-12 101 N=2 N=3
7.1 Higher-order time-dependent problems
2
h1
h3 h3 h5
.. but an even/odd pattern
K 1/h
102
Fig. 7.2. Convergence of the scheme for solving the heat equation using a centr flux.
Section 6.5 when solving Maxwell's equations using a central flux. For th heat problem, this behavior is confirmed in the following theorem [78]:
Theorem 7.3. Let u = uh - u and q = qh - q signify the pointwise erro for the heat equation with periodic boundaries and a constant coefficient a(x computed with Eq. (7.1) and central fluxes. Then (T )
2
+
T
(s)
2
ds Ch2N ,
|| u - uh ||
The heat equation
10-6
-8
10-4
N=3
h3
7.1 Higher-order time-dependent problems
100 N=1 10-2 h1 N=2 N=3
2
Back to the example 10 Looks good 10-10 10-12 101
|| u - uh ||
N=4
h3 h5
h3
10-4 10-6
-8
h K N=4 1/h .. but an even/odd 10 Fig. 7.2. Convergence of the scheme for solving the heat equation using a central pattern h
3
102
flux.
5
10-10
-12 Section 6.5 when solving Maxwell's10equations using a central flux. For the 101 102 heat problem, this behavior is confirmed in the following K 1/h theorem [78]:
Fig. - Theorem 7.3. Let u = uh 7.2.uConvergence ofhthe scheme for solving the heat equation using a centr and q = q - q signify the pointwise errors flux. for the heat equation with periodic boundaries and a constant coefficient a(x), computed with Eq. (7.1)Section 6.5 when solving Maxwell's equations using a central flux. For th and central fluxes. Then
for the heat equation with periodic boundaries and a constant coefficient a(x where C depends on thecomputed with Eq. (7.1)and N . Forfluxes. ThenC is O(h2 ). regularity of u, T , and central N even,
u (T Theorem 7.3. Let(s) = uhds u and q = qh - q signify the pointwise erro ) 2 + q u 2 - Ch2N , ,h ,h
0
T
heat problem, this behavior is confirmed in the following theorem [78]: T
The proof is technical and can be found (T ) [78].+The computational 2N , in (s) 2 ds Chresults
2
here Vh is the space of N -th-order polynomials with support on the element, in mind that we have considerable freedom in choosing the . If we now choose the test and we can use the stability considerations to guide this cho functions as h = u h , h = q h ,
the schemes Can we do anythingthewith optimal convergence obtain at each interface toproof of Theorem this? rates. improve on 7.2 we nd use that From
T ~ a uh S q h
The heat equation be gained by recalling that upwind fluxes mos inspiration can
+ + h = Recall theustability condition qT S h
a T ~ a uh S q h
uT (S h
e recover252
here
uh D + qh D + r - l = 0, dt 252 at each 7 2Higher-order equationsshows that interface. One easily
or, alternatively, through + - (One easily {{ , aq at, each interface. auh ) = a}}uh ( is)guaranteed. with the flu shows that this h = - a- qh onsider xr where we have the term ( auh ) = {{ a}}u+ , ( aqh ) = a- qh . h
fluxes rst, assume that a(x) is continuous and use the- - u centralaq ) = {{ a}}q + , ( auh ) = a h , ( h h or, alternatively, through 252 7 Higher-order equations ( aqh ) = a{{qh }}, ( auh ) = a{{uh }}.
this - - + 0 with the flux c is guaranteed r l - aq ) at each interface.hOneau h) u=- (a- u ) this isaq ) = {{a}}q + , flux choice auh = auh q - ( easily shows that ,qh . h ( ( guaranteed with the h h h h
7 Higher-order equations To guarantee stability, we must ensure 1 d 2 2
x r , )= h = q u q auh qh x aq ) u - ( au ) q . h h a-( l h h h h
a T
Stable choices
^ [[ auh ]], ( h ) {{ aqh }} [[ aqh ^ ( all elements, immediately aq = ummering overauh ) = {{ through recover we ^ ^ - or, alternatively, auh }} +
( auh ) = {{ auh }} + [[ auh ]], ( aqh ) = {{ aqh }} - [[ aqh ]],
- + solution - a A slightly more - ( auhq) + u+ q- - u is, ( aqh ) = {{ a}}q + , symmetric = a . r = uh h h h A slightly more symmetrichsolution is h 2
2 ^ can 1 d taken as nor -n,0, = ^ a}}u+ , Here, can be taken2 as + ^ h-n,{{where essential property - q -is the diffe be uh(au or = LDG flux (aq ) = ^ = n. ,h n q ,h where the the essential propertydifference ^ ^h ) Here, is theh Upwind/downwind -^ 2 dt a ^ h h in sign betweenthe two fluxes. in sign between the two fluxes.
The heat equation
Back to the example Looks good .. full order restored
it is often computationally advantageous t an aliased solution and stabilization if need
100 10
-2
10 N=1 N=2 N=3 N=4 h
2
10-
|| u - uh ||
10-4 10-6 10-8 10-10
|| u - uh ||
h3 h4 h5
10-
10-
10-
10-1 102
10-12 101
10-1
K 1/h
Fig. 7.3. On the left, we show convergence of s equation using an LDG flux. On the right is sh used for solving the problem with homogeneous
The heat equation
Back to the example Looks good -
it is often computationally advantageous t an aliased solution and stabilization if need
10 7.1 Higher-order time-dependent problems 10
-2 0
N=1
253
10
|| u - uh ||
h .. full order restored As elegant as this construction appears, the motivation for considering it 10
4 -8
10-6
N=4
|| u - uh ||
N=2 % compute right hand sides of the semi-discrete PDE rhsu = rx.*(Dr*q) - LIFT*(Fscale.*(nx.*dq)); 10-4 N=3 return
h
2
10-
h3
10-
10-
is whether one can improve the convergence rate for N being odd. In Fig. 7.3 h5 we show results for solving the heat 10-10 equation, and these appear to show that optimal order of accuracy is restored. 10-12 This observation is confirmed by the 101 following theorem, the proof of which 102 can be found in [78]. K 1/h
10-
10-1
10-1
and q = - h signify the pointwise errors Theorem 7.4. Let u = u - uh Fig. 7.3. qOnqthe left, we show convergence of s for the heat equation with periodic boundaries and a constant coefficient a(x), is sh equation using an LDG flux. On the right computed with Eq. (7.1) and LDG fluxes. solving the problem with homogeneous used for Then u (T )
2 ,h
+
0
T
q (s)
2 ,h
ds Ch2N +2 ,
where C depends on the regularity of u, T , and N .
as convection-diffusion problems, is of the type 7.1 Higher-order time-dependent problems With the approach for the heat equation in place, it is natural to consider 255 7.1.2 Extensions to mixedinstance, a very important class of problems, known and higher-order problems u further extensions. For u 7.1.2 Extensions to t + x f (u) = x a(x) x , mixed is of the type as convection-diffusion problems,and higher-order problems
With the approach for the heat equation in place, it is natural to consider subject to appropriate for u With the approach boundary and initial conditions.it is natural known the very equation in class further extensions. For instance, +heat (u) = a(x)place, of problems, to consider a f important u , WeTo develop a suitable and match whatx know it in firstcan now mix t of a forx problems, we write such we further extensions. Fordiscretizationtype instance, as convection-diffusion problems, is x thevery important class of problems, known
Higher order and mixed problems
Consider
orderconvection-diffusion problems, is initial conditions. as form appropriate boundary and of the type subject toas u To develop a suitableu f (u)= for uproblems, we write it in firstdiscretization a(x) + u f (u) - such ,= 0, ++ aqa(x) u , (7.2) t t x f (u) = x x x order form as t x x x 7 Higher-order equations 252 initial conditions. subject to appropriate boundary andu ,and initial conditions. qu +a = subject to appropriate boundary f (u) - aq = 0, (7.2) at each interface. One easily x t x for such problems, we write it in first- shows tha To develop a suitablesuitable discretization for such problems, we write it in firstTo develop a discretization have byorder combining the order form as form as knowledge we a u , from solving both conservation laws anda- u- , q= ( auh ) = h x heat equations. We can now discretize this exactly as above. The only problem u u f 0, = 0, that combining the knowledge we have- (u) solving both conservation laws and is aq) by requires attention + the choice from -numerical fluxes, (f (u) - (7.2) +(u) of the =aq or, alternatively, through (7.2) f aq x heat au) . can x discretize this exactly as above. The only problem and ( equations. We t nowt + u u au ) that Higher-order equations choice , the it is natural to (f a( - haq)= 252 requires attention isathe , discussions, numerical fluxes, use(u) monotone{{ a}}uh 7 However, based q =our q = a of on past x x and ( ; . flux for fau)for example, a Lax-Friedrichs flux like A slightly use at each interface.the knowledge wethat this is guaranteed with the fluxmonotonesolution is choice our past have from it is natural to moreasymmetric by However, based on we have from solving both conservation lawslaws and combining One easily shows discussions, solving both conservation and by combining the knowledge C flux for f ; for example,{{f (u)}}- thisC max(|fabove. = {{au }} + [[au ]] + ( [[u]], = {{ above. f (u) = now - u this exactly like as (u)|. The problem^ heat equations. aucan = adiscretizeaq flux as a}}q +auh ) only only problem We ) a Lax-Friedrichs ) exactly ,The heat equations. We can now discretize h , 2 h h h ( h h - aq) that C that requires requires attention choice choice of the numerical fluxes, (f- attention is the is the of the numerical fluxes, (f (u) (u) aq) for the parts corresponding or, alternatively, through {{f (u)}} + [[u]],toC max |f can be taken as n or -n, whe f (u) = ^ the dissipative operator, ( ^ aq) ^ and. Here, (u)|. 2 and ( Furthermore, . au)( au) and ( However, based on our au) , we rely on the results of the previous natural tochoose two fluxes. section and use monotone in it is sign between the either a However, based on(au parts {{ past discussions,natural - . use a monotone our past discussions, it is =dissipative choice ofathese fluxes, leadi + - q to operator, ( aq) Furthermore, the the h ) flux, Lax-Friedrichs ) to the a}}uh at central f ; for example, acorrespondingaqhflux theapurely diffusive problem, flux or flux for for LDG = which, , (least for like These h example, a Lax-Friedrichs flux like section and choose either a flux for f ; forau) , we rely on therates. of the previous discontinuous Galerkin (LDG) metho and (optimal convergence results yields A slightly more the LDG flux, which, atC central flux or symmetric solution is mixed approach, let us consider the even if it is purely diffusive problem, To illustrate the (u) = {{f (u)}} + least for themax |f (u)|.upwinding, folthis [[u]], C tively utilize f behavior of C
and rewrite as
Now choose fluxes as we know how
u(x,t)
boundary conditions. 1 t=0.5
0.75 0.5 0.25 0 t=0.0
1.25
t=1.0
||u-uh||
central flux or the LDG flux, which, at least for the purely diffusive problem, yields optimal convergence rates. yields illustrateconvergence rates. mixed approach, let us consider the folTo optimal the behavior of this lowing example: the behavior of this mixed approach, let us consider the folTo illustrate lowing example: Example 7.5. We consider the solution of the viscous Burgers' equation Example 7.5. We consider the solution of the viscous Burgers' equation u u2 2u + = 2 , x [-1, 1], t x u u x 2 2 u + = 2 , x [-1, 1], t x 2 x which has an exact traveling wave solution which has an exact traveling wave solution - t x + 0.5 u(x, t) = - tanh + 1. 2 - t x + 0.5 7.1 Higher-order time-dependent problems 257 u(x, t) = - tanh + 1. 2 To solve this, we use a standard DG method with a Lax-Friedrichs flux for the 2.25 nonlinearthis, we use au2 /2, and central fluxeswiththe dissipative operator. In To solve flux, f (u) = standard DG method for a Lax-Friedrichs flux for the 10-2 2 BurgersRHS1D.m,(u) = u2 /2, and central fluxes for the this approach N=1 the In for nonlinear flux, f we illustrate the implementation of dissipative operator. 1.75 evaluation of the right-hand side, the implementation of being approach for the BurgersRHS1D.m, we illustrate with the exact solution this used to impose 1.5 10-4 boundary conditions. evaluation of the right-hand side, with the exact solution being used to impose
Higher order and mixed problems
Consider viscous Burgers equation
t=1.5
N=2
10-6
N=3 N=4 102
10-8
-0.5 0 x 0.5 1
-0.25 -1
101
K
agreement with what we found what we found in the The possible loss The possible loss flux, in agreement with in the above example. above example. racy from the linear operatorlinearnot appeardoes not appear to impact this, although of accuracy from the does operator to impact this, although based on experimental evidence only. evidence only. this is based on experimental us briefly Let us brieflygeneral higher-order spatial operatorsspatial operators and how also consider also consider general higher-order and how etize these. As an example, consider the problem
Higher order and mixed problems
3u u = . t x3
to discretize these. As anorder dispersive equation Consider the 3rd example, consider the problem 3u u = . t x3
we need to understand what kind of boundary conditions are required First, we need to understand We cankind of re well-posedness boundary conditions doboundary conditions are required in a finite domain. what understand need? Whichwell-posedness in a finite domain. We this through this through wecan understand to ensure rgy method
the energy method 258 7 Higher-order2 equations
1 d u 2 dt
2
where = [xl , xr ]. To u = u x2 - ensure well-posedness we , need to impose boun l 2 dt x 2 x xl conditions like 2u u x = xl : On u or , and 2 x x
= u
u 1 u - , 1 d2 2 x 2 u x 2
2 xr
x must 2be bounded 1 u
r
2u x = xr : On u or . 2 x As expected, three conditions are needed, although it is essential that tw imposed on the left-hand side of the domain.
can are needed, To form is this problem, as for exactly cted, threefirst-order equations in the discretize essential that we proceed as heat 2equation; that ystem of conditions be done inalthough it the same way two be the for the u heat equation and u on expected, threethe elementwise operator, althoughl, :in ; essential as x sum side of the As the left-handforms q domain. a system of first-orderhequations ithOnhu, h ,xh ), that ,two b conditions are u needed, B (ux =qx p the form 2 and the is or h, h it as u scretize this problem, we proceed p forpthe Since Bh (uh ,and, write h , h , h ) = 0 for any N -th as , three first order terms.the. domain. qh ph ; p = , mposed on theequations q = side of x equation q left-handthe form = heat u u u 2 t x in x ystem of first-order x , x = ), u = or .2 . Vh when for , qqr, hOnis a solution, we can polynomial (h , h , h ) proceed t =(uh=the:pheatp equation and writ x h x x x Tobe discretizedthisa problem, we manner. Proving stability discretize in straightforward as n now u h , h , h ) =This , -ph , qh ). discretized in a straightforward manner. Proving sta (u As expected, three conditions ( = q , q = p , hpcan u . be A bit of manipulation yields a local eleme = now t as a exactlyt of x first-order equations equation; that are needed, although it is essent done in systemthe same way as for the heat in the form is, one x can done in exactly the same way as for the heat balance as x beimposed on the left-hand side of the domain. equation; that is he elementwise operator, Bh (uforms, ph ; h , h ,1 h ), as thestability the Proving sum (u of h , qh the elementwise as the equa n now be discretized in a straightforwardTo discretize operator, Buh , qh , ph ; h , h , for),the heatsum o manner. d p 2 h this problem, - u h ; hqh , h ) = 0 forhanykN -th-order l , u , DB = r we ; , , as =h0 for any N -th= , q ,.p proceed ) p firstheat equation; that is, one h h st order terms. Since B way as, for= , q = done in exactly the sameh (uh , qhthree the,asorder terms. Sincep h (uh equations inhthe form h 2 dt it a system of first-order h mial (h , h , hoperator,when h ,t,, qh , pxis, ahh ,ash ) sum canx h , qh , ph ) is a solution, we can c (u polynomial (h ), x we of the e elementwise ) Vh Bh (u (uh ph ; h )2 h , solution, Vh when choose qh the h, ph = (uh , -ph , qh ). Aelementwise p yields alocal elemen u ) = (uh , -p Since Bh (uh , qh , manipulation = 0 q +local-th-orderq(u q)= ). bit of p(h ,h ,h , )h -yields a u N bit of manipulation , p = . , h sthorder terms. h , qhbeAdiscretized=in a) straightforward=manner.- ph (ph )u .stabilit h;h uh for any h (qh ) + qh , h h This can now Proving t x x as as (h , h , h ) Vh when (uh , balance is2 a solution, we can choosex qh , ph ) ial 1 d can=beh ,doneqWeAcan uh manipulationcan thelocal discretized2 in= r - l , to stable is, on in 1 bit of 2 k = to guide now choice dofu fluxesa equation; that schem This straightforward manner. Pr h ) (u -ph , h ). exactly this r - yields a as for theh heat leading useD the samel ,way be elementwise Dk 2 dt 2 dt can Bh (u ,2 h ph all variables asfor as orms the elementwise operator,be donehinqexactly h , samewayyields thesum equatio , example, the use of central fluxes ,on ; the h h ), as the heat of th ph 2 1 d 2 , ) = - h h , = 2 h); u q + u (qh ) B h for p l, p . three first ph - 2 hterms.k(qh )r+ qB(u(utheqelementwise, h ,hh ) =hq(uh ,hqh ,-hp;hhN -th-ordet orderdt huh uh Since forms ) - ph,(ph-hh h operator, + 0(u ) any(phh., h ), as u q + D = h h h 1 terms. - B three first order + - Since + h (uh , qh , ph ; h , h , h ) = 0 for a 2 u,h q, + is when- + q , p ) isstable scheme We can use = to guide hph choicehofsolution, this fluxes leading to 2 (h , h , h ) Vh when (uh , qh ,the ) uh qa - ph ph , we can choos polynomial polynomial ( h use this = ph - uh qh +choice) use example,)the of )central fluxes Vh all for (uh , yields a solution, to guide the uh (qh of fluxes leading 2tohstableh ) on schemes; variables h h + q (u - ph (p . = ( bit h , use of central h , -ph , all variables ,ofh )manipulation yields a local elementwis (theh , h2) = (ufluxes on qh ).hAhh , hyields (uh , -ph , qh ). A bit of manipulation yields a loc e, h which, after fluxes balance as stable schemes; vanishes. This establishes energ 1 use this guide the choice of summation over alledges, + - + leading to uh 1 d h + balancetoas servation in the periodic case, as = 2 theqhfor u- qh2- p- p+ , h for continuous hequation. , the use of central= 1 u+ q - +variables yields , on uh Dk = r - l , fluxes h hall u-1 hd p- p+2 q+ - h h 0 h 2 uequation,r - lalso vanishes. LDG-type upwind = over can 2 dt use an This establishes energy , the heat summation we all edges, 1Similar towhich, after h Dk + - - + dt + 2 - p2 choosingedges,servation in the periodic case, as for the continuous after summation= 2 uh qh + uh qh - ph ph This establishes energy + uh (qh ) + qequation. ph (ph ) . over all vanishes. , = h - uh qh conh (uh ) - 2 Similar to the heat LDG-type upwindi n in the periodic case, asph the continuous= u- , equation, we + , (ph )use an - , for equation. h2 = q can (uh ) establishes energy con- also = ph (q ) h fter summation over all edges, - uh qh + uhuse this+ qguide ) h - ph (ph ) . leading to stab vanishes. This choosing can = we can also use an(qh ) to h (uh the choice of fluxes We lar to the heat equation, 2 the continuous equation. upwinding by LDG-type - + - n in the periodic or case, as for example, the(uh ) of central h ) = qon all hvariables yields use = uh , (q fluxes h , (p ) = ph , g to the ar an - We can heat h ) = utowe can or the (uhh)) = of ,,fluxes leading(ph ) = +- . schemes; fo useequation, , guide= q+ , choice p- upwindingqby, + to stable - + this - (qh ) also use (p LDG-type (qh ) = 1h u+ p (u = h h - h h + = u q- + u- q h - p -p , . h h h example, the use of central fluxes on all h ) = uh (qh ) = qh , (ph ) h= ph h h , yields - + - (u variables 2
x
Higher order and mixed problems
Write it as a 1st order system
To choose the fluxes, we consider the energy
Central fluxes yields
Alternative LDG-flux
tleus conclude with an example. of the dispersive equation 7.6. We consider the solution
he exact solution is given as 260 7 Higher-order equations roblem is solved using both the central fluxes and the upwind-style u(x, t) = cos( 3 t + in Fig. Centralresults shown x). 0 7.5. LDG flux fluxes, with0 the convergence flux
10 10
exactly as for the 2nd periodic boundary conditions and3 t + x). = cos(x). One easily shows order problem u(x, t) = cos( u(x, 0)
3u u ple 7.6. We consider the solution of [-1, 1], = , x the dispersive equation 3 Consider t x 3 u u eriodic boundary conditionsand u(x, 0) = 1], Convergence behavior = , x [-1, cos(x). One easily shows 3 t he exact solution is given as x
Higher order and mixed problems
problem 10-1 solved using both the central 10-1 is fluxes and the upwind-style N=1 fluxes, with the convergence results shown10-2 Fig. 7.5. in 10-2
N=1
||u-uh||
||u-uh||
10-3 10-4 10-5 10-6 10-7 101 N=2 N=3 102
10-3 10-4 10-5 10-6 10-7 101 102 N=3 N=2
K
K
Higher order and mixed problems
Few comments
The reformulation to a system of 1st order problems
is entirely general for any order operator
When combined with other operators, one chooses
fluxes for each operator according to the analysis. require two derivates rather than one.
The biggest problem is cost -- a 2nd order operator There are alternative `direct' ways but they tend to
be problem specific
nsider problems of the type With the ability to solve problems with higher spatial derivatives we ca u Now we of the type = 2 u - f (x), onsider problemscould consider solving a problem like t x2 2u u = - f (x), th appropriate boundary and initial conditions. For this problem, any of 2 t x e methods discussed in the previous section can be applied. with however, we arewe areandin a steady-state solution (i.e., utwe0), any If, appropriate boundary interested in the steady state = However, if interested initial conditions. For this problem, it is he natural to ore methodsbeconsiderin the previous section can be applied. may discussed the elliptic problem better off considering If, however, we are interested in a steady-state solution (i.e., ut = 0), it 2 u problem more natural to consider the elliptic f (x), = Higher-order equations x2 2u = f the d invert the discrete operator to recover(x), approximate solution, uh . We 2 We can use = 0 Poisson methods wewell known that under mild to u(0) = as the classic of thex the exactjust derived to = sin(x u(2) any and with cognize this problem. It is solution u(x) obtain to linear system CRHS1D.m thebuild the discrete approximation nd invert the discrete operator to a unique solution [59]. nditions on f (x), this problem has recover the approximate solution, uh . W ecognize this as the classic Poisson problem. Itproblem As an example, let us consider the following is well known that under mi Auh a unique solution [59]. onditions on f (x), this problem has= f h , 2u As an example, let us consider the following problem cond order problem by = - sin(x), x [0, 2], x2 simply calling HeatCRHS1D.m with the N
7.2 the ability problems ith Elliptic to solve problems with higher spatial derivatives we can
Elliptic problems
Elliptic problems
Assume we use a central flux.
When we try to solve we discover that A is singular!
0.01 0.005 imag() 0 0.4 0.3 0.2 0.1 0 v -0.1 -0.2 -0.3 -0.01 -15 0.4 0.3 -12 -9 -6 real() -3 0 3 -0.4 0 0.4 0.3 0.2 0.1 0 v -0.1 -0.2 -0.3 1.57 3.14 x 4.71 6.28 -0.4 0 1.57 3.14 x 4.71 6.28 1.57 3.14 x 4.71 6.28
trix is singular. This is in contrast to the continuous problem, which is uniquely solvable. The question arises of what is causing this sudden breakdown and why we did not experience problems when solving the heat equation. The singularity of the discrete operator indicates that it has at least one zero eigenvalue in contrast to the continuous case in which there are none once the boundary conditions are enforced. To understand the severity of the problem, we show in Fig. 7.6 the computed eigenvalues of A close to the real axis, confirming that all eigenvalues are negative and real but also that there is at least one zero eigenvalue. A closer inspection reveals that there is exactly one such zero eigenvalue and in Fig. 7.6 we also show the corresponding eigenvector, v, for N = 1, 2, 4 and K = 6.
-0.005
N=1 K=6
N=4 K=6
v
0.2 0.1 0 -0.1 -0.2 -0.3 -0.4 0
N=2 K=6
This also explains why this particular spurious mode caused no pro in the time-dependent problem as it is characterized by being constant in Elliptic problemsspace. As long as the solutions are well res and very rapidly varying in this mode will not cause any problems. However, for marginally resolved What is happening? lems, the situation may be a little more delicate and the spurious mod impact the accuracy. In such cases, the stabilization methods develop The case can also be useful to improve -- accuracy one the elliptic discontinuous basis is too richthe it allows for time-depe extra problems. null vector: With this added understanding, we can ask ourselves what to do it. An obvious solution would be, at least asymptotically in N and/ A local null vector with {{u}}=0 to disallow the eigenmode with the properties in Eq. (7.3) without i vanishing inside the element; that is, to remove the artificial null-space What can we Galerkin (DG) method now shows one of its stre The discontinuous do ? by allowing us to modify the choice of the numerical flux. In particul us consider a set of flux by penalizing this mode Change the numerical fluxes as q = {{q}} - [[u]], u = {{u}}.
For = 0 this reducesof DG central flux. The role of the added term The flexibility to the shows its strength! penalize the solution to disallow large jumps in u; clearly if we introduc
us the one-dimensional Poisson equation. % Set up operator
% Set RHS f = a simple example to illustrate consider -J.*((invV'*invV)*sin(x));
Elliptic problems
[A] = PoissonCstab1D(); solvec = A\f(:); d2 u u = reshape(solvec,Np,K);
Example 7.7. Consider the problem
the performance when
e 7.7. Does itthe problem Consider work? % Solve Problem dx2 = - sin(x), x [0, 2],with u(0) = u(2) = 0.
d2 u = - sin 2 dx
The routine for setting up the d 10 0 sonCstabRHS1D.m, and in PoissonC 0) = u(2) = 0. N=1 10 routine for setting up the discrete Poisson operator is shown to solve the problem routine needed in Pois10 N=2
0 -2
|| u - uh ||
20 bRHS1D.m, and in PoissonCstabDriver1D.m, we illustrate the driver -4 N=3 10 30 needed to solve the problem. Poisso h3 N=4 40 10-6 function [A] = PoissonCstab1D(); h4 10-8
-10
h2
PoissonCstab1D.m on [A] = PoissonCstab1D(); h5
10
50
tion [A] = PoissonCstab1D(); 10-12 101 102 ose: Set up symmetric Poisson 1/h matrix with estabilized K central fluxes Globals1D;
% 60 function [A] = PoissonCstab1D() % 70 Purpose: Set up symmetric Poiss % 80 0 central fluxes 20 40 60 80
nz = 622
What about the s1D; ros(K*Np); g = zeros(K*Np,1);
Fig. 7.7. On the left, we show the convergence rate for the solving the homogeneous Poisson problem with a stabilized central flux. On the right is shown the = zeros(K*Np,1 A = zeros(K*Np); g sparsity pattern for K = 20 and N = 3 of the discrete Poisson operator based on this flux.
other flux - the LDG flux?
The problem is solved with different number of elements and orders of
% Build matrix -- one column at a
One can naturally question whether the LDG flux, introduced to reco optimal convergence in Section 7.1 for the heat equation, can also be used discretize the Poisson problem. This is naturally possible and we show in 7.8 the results of solving the problem in Example 7.5 using a local stabili Consider the stabilized LDG flux flux of the kind
Elliptic problems
268
^ ^ with = 1 for all cases and = n. As0 in Fig. 7.7, we recover an opti 10 N=1 convergence rate. The code used to solve 10 problem is a direct combinat the 10 N=2 h 20 of PoissonCstabRHS1D.m and HeatLDGRHS1D.m. N=3 10 30 If we now compare the twoh plots of the nonzero elements in Fig. 7.7 N=4 10 Fig. 7.8, we note that the latter, based 40 the LDG fluxes, is considera on h 50 sparser that the former, based on the central fluxes. 10
0 -2 2 -4
^ ^ = {{qh }} equations qhHigher-order + [[qh ]] - [[uh ]], u = {{uh }} - [[uh ]], h 7
|| u - uh ||
3
-6
4
-8
h5
60 70 102 80 0 20 40 nz = 472 60 80
10
-10
10-12 101
K 1/h
Fig. 7.8. On the left, we show the convergence rate for the solving the homogeneous Poisson problem with a stabilized LDG flux. On the right is shown the sparsity pattern of the discrete Poisson operator based on this flux.
Works fine as expected - but we also note that A is much more sparse!
To appreciate the source this, Fig. On On the left, we show consider rate simplest schemes homogeneous Fig. 7.8.7.8.the left, we showofthis, the convergence the for thethe with N =homogeneou for To appreciate schemes with N Ellipticproblem the source aofthe convergence the the ratesolving homogeneous=0.0. sparsi problemselementconsider the simplestthethe solving the the In this case problemshow the convergenceLDGOn thesolvingis right isthe sparsity the scheme stabilized LDGbecomes for stabilized flux. for k rate flux. On Poisson Fig. 7.8. On the left,with with we a Poisson right shown shown
10-12 101
101
K 1/h K 1/h
K 1/h10
102
2
80
0
0 20
20
40 60 nz = 472 80 40 nz = 472 60 80 nz = 472
Why is one(qmore usparse-than, qthe, uother?= hfh , qh h ,k h k+1 h , uh )k qh (qh h k h , uh ) q , k-1 Consider the N=0 case
the a k+1 this Poissonpattern ofwithk discretekPoisson operator on this onk-1 flux.the sparsity problem stabilized k+1 flux. k the k flux. LDG pattern of the discrete Poisson operator On k-1 right is shown k based based q k qh , h h pattern of the discrete , k+1 , uoperator - qh (qon qk-1 flux. uh ) = hfhk k+1 h (qh Poissonkh , uh ) based k , this , uh , k-1 k
In this case the scheme for element k becomes
To N= To In case the the scheme for element appreciate the source of element k becomes In we first consider the central flux with the simplest schemes with N = 0. If thisthis case scheme for this, consider k becomes In this first consider the central flux with If we case the scheme for element k becomes k k - + k+1 k+1 k- qk (qk+1 - ,,,quh uh k h }} - (q-hq k-1 , , q k-1 , uk ={{u k , = , + qh (qh, qh ,qh kku+k ,=,k+1, uk+1 [[u, qh (q kuk-1 k-1 )=uhfh }}, hf k , h + q , uh + ) h - u+ (qu , hhk+1 ){{q ) h qk ,) k-1 ,uk hu hhu) =)hf k ,k-1 ) uh q - h q h ]], u, (uh , h h h , h (q , q h , h h (q - , h , uh ) = - h - h h k h u qhqhh ,hqh hku-k+1 uh {{qhk}}(qh [[uh ]], h h (uh , uh ) = {{uh }}, h h k-1 u k (uuuh k ) - uh (uk-1 h k) = hg ). , k+1 k-1 k k+1 k ) ,u we recover (uh , uhh (uh , u (uk , -huh) = , uh. h = hgh . h k h uk-2 (uh hgh uh h hk+2 - uh ) h uh - 2uh + uh uk+1 - uk-1 we recover k h If we first consider the central uk-2 with uk+1 - uh = fh . k+2 k-1 (7.4) k 2 flux + uh central + with flux h h h - (2h) central If we first consider 2uhflux h the with If we first consider the =+ h . fk + (7.4) - (2h)2 + - + hu (u- , u ) = {{uh }}, is, the , stencil {{q qh (q +, q the +, ) + h }} - [[uh would h - One notices thath- uh+ uh - =is wider than ]], - be expected;+that h h - h - u h qh u , uh ) = {{qh }} u (uh , ]], = {{u , uh qh (qh , qqh,(qh ,,uh ), =h{{qh }} - [[uh ]], - [[uhu+ )uh (uh h }}, ) = {{uh }}, h h h standard second-order finite is wider than would be the nearest neighbors One notices that the stencil difference method uses onlyexpected; that is, the we recover (7.4), standardrecover u central difference method -which the nearest levels of we in contrast to the k+2 - 2uk + uk-2 Eq. uk+1 usesk-1 engages two neighbors recover second-order finite stencil in we uhonly k (7.4) uk+2h 2uk+2 - 2uk + +k+1 h uk-1 uh - h which =uk-1 - u k stencil uh in elements. to the central +h2k-2 hin Eq. (7.4), uk+1 -engages two levels of contrast k fh . h h uk-2 k (7.4) h (2h) LDG fluxh + h = fh . h h h + If we = fh . (7.4 elements. instead consider 2the (2h) h 2 Onewe instead consider - + (2h) - than would - h+ notices that the stencil is wider be expected; the wider flux would be expected; thatthat is, the - the + One If notices thath (qh , stencil , ishLDG h than[[uh ]], u (uh , uh ) = u+ , is, the qh the ) = q - q h standard notices that, finitestencil is methodthan only thebeh expected; that is, th One second-order uh u difference wider uses would nearest neighbors standard second-order + - difference method uses only the+ finite nearest neighbors - in contrastqh (qh , qh , uh , u+ ) = qhin Eq. h ]], method, uses = u+ ,two nearestof to the central stencil- - [[u (7.4), which engages uh (u- standard second-orderhfinite difference which h uh ) only the levels neighbo h we recover the central in contrast to the stencil stencil in Eq. (7.4), engages two levels of elements. in the stencil in Eq. (7.4), which engages two levels elements. contrast tok+1 central k-1 we recover instead consider2uk + u flux uk+1 - uk-1 the stencil - the LDG If we uh elements. k h h h If we instead consider the LDG flux + h = fh . (7.5) k+1 consider the- k + h2 + uk-1 LDG k+1 - uk-1 + h - If we qinstead -u2uh ) = q - uh ]], u h , u )k= u+ , flux - + u (qhh, qh , - , uh - h + [[uh - (uh =hf + h h - h + + - + h h . (7.5)
u (u , uh ) - uh (uh , uh ) = hgk . h h k h k+1 uh (uh uh ) - u this, k-1 ) consider appreciate ,the source of(uk ,ofconsider the simplest schemesschemes= 0. To appreciate the source uhthis,= hgh . the simplest with N with h h
Using the central flux yields
Wide
Using the LDG flux yields
If the internal penalty flux,of sparsity as LDG no longer guarantees nearest neighwe compute the defined as Trying to one observes that dramatic in terms condition number, ,balance sparsity and the conditionin bor connectionsqhon{{(uh )x }}the internal penaltyCflux, remain sparser than the grids,h = {{u 2(A will defined as = general (ALDG ) h }}. - [[uh ]], u although it ); multidimensional 0, the flux qofdepends solely on the penaltyterm.283]. We note that for N = version h the central fluxes [252, This that is, moreconditionshortly, that operator/stencilqdoes, however, - [[uatis about the compact discrete the value of plays a larger discrete xoperators ]], price. number of the LDG-based= {{(uh ) }} come h a u = indicates, as we will discuss h role for h The the internal sparsity is for the twiceThe penalty fluxcondition previous two cases. but it comes at a of this is a that of the than a good on centralobserves that -- = 0, the If weNevertheless, to firstoperator basedthingapproach fluxes. The impact price solely compute the convince Wenumber, , for N works, we re- flux q depends ourselves that this one note that h slightly problem in Example 7.5of direct solvers and increased iteration counts when peat the worse accuracy using the internal penalty flux, and we show indicates, as we will rates discuss shortly, that the value the results in solver is used to solve observe 2(AC ); of hope, we the optimal an iterative Fig. 7.9. As one could (ALDG ) linear system.convergence and, when comparing with Figs. 7.7 and 7.8, a sparsity pattern in for the previous two the and flux Trying operator in between those obtainedthe penalty flux discrete operators is about to balance sparsityinternalLDG-based anda possible compromise is conditioning,than the discrete condition number of the the that We seek a flux balancingwith the central first convince ourselves that t is, the Nevertheless, to and conditioning? sparsity LDG flux. Furthermore, the condition number behaves as the internal the operator based on defined twice that ofpenalty flux, peat theas central fluxes. The impact of this is a problem in Example 7.5 using the interna (AC ) (AIP ); slightly worse accuracy= {{(uhresults [[uFig. u = As one could counts we ob direct solvers h ]], 7.9. {{uh iteration hope, when qh ofthe )x }} - in and increased }}. h that is, the internal penalty method appears to offer a suitable compromise an iterative solver is used to solve linear system. the vergence depends solely on the with Figs. 7.7 and between central fluxes N = 0, fluxes.flux q and, the implementation of and LDG the An example of when comparing penalty term. This WeTrying to balance sparsity and h conditioning, a possible compromise is note that for the internal penalty flux is shown in PoissonIPstabRHS1D.m. Internal penalty flux the discrete operator in between those obtainedfor indicates, aspenalty flux, defined as that the value of plays a larger role w the internal we will discuss shortly, LDG flux. previous two cases. the internal penalty flux than for 0theFurthermore, the condition number beh 10 N=1 qfirst {{(uh )x }} - [[uh ]], that this h }}. = convince ourselves u = {{u approach works, we re10 Nevertheless, to h h 10 (AC ) (AIP ); N=2 h 20 peat the problem in Example 7.5 q using the internal penalty flux, and we This We note that for N = 0, the flux 30h depends solely on the penalty term. show 10 N=3 the results we will discuss shortly, that the value of method rates of conh that is, the internal penalty optimal appears to o indicates, asin Fig. 7.9. As one could hope, we observe plays a larger role for 40 10 between central 7.7 and 7.8, a Mission vergence and, when comparing with previous two cases. sparsity pattern in the internal penalty flux than for50the Figs. fluxes and LDG fluxes. An examp N=4 10 internal penalty that shown central flux we retheNevertheless, to first between 60ourselvesflux isthis approach works,and the discrete operatorh in convincethose obtained with the in PoissonIPstabRHS accomplished LDGthe problemh in Example condition number behaves as flux, and we show flux. Furthermore, the 7.5 70 peat 10 using the internal penalty
Elliptic problems
0
-2
2
-4
|| u - uh ||
3
-6
-8
5
5
-10
0 20 40 60 80 the results in Fig.K 7.9. As one could hope, =IP );observe optimal rates of con10 C 1/h 0 nz 586 (A0 ) (A we 101 102
80
Elliptic problems
Remaining question: How do you choose ?
270
The analysis shows that:
7 Higher-order equations
Let us now return to the question of how to choose the penalty param ter, . For the central flux, in Section 7.2.2, for schemes based on cent As we will show rigorously > 0 suffices or LDG fluxes, it suffices to > 0 > 0 to ensure invertibility and a pur take For the LDG flux,48, 49]. suffices for the internal penalty-bas negative real spectrum [13, However, For mustIP 7 Higher-order require scheme, we the 270 flux, that must equations require one (N + 1)2 , question of how to choose th Let us C now return to theC 1, h ter, . As we will show rigorously is derived in [281]. to guarantee similar properties. The bound on Cin Section 7.2.2, for schem or LDG ensure invertibility take These suffices to guarantee stability, > 0 to ensure invertib Whereas these resultsfluxes, it suffices to and definiteness, the actual val but they negative the spectrum [13, accuracy may not give the best 48, 49]. However, the scheme. of may well affect both real conditioning and the accuracy of for the inter fact, as we will scheme, we must require that quite a different impact on t see, the actual value of has (N + 1)2 schemes. Generally, a good choice is C , C 1, h
PoissonIPstabRHS1D.m
ctions, let us develop some of the basic numerical flux choices approximations of Table 7.3. Overview of theory for DG 288 7 iptic problems.Higher-order equations u q h h To kee...