2008
Discontinuous RMMC Galerkin methods Lecture 7
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 nite8.4 Scattering
y
-0.25 -0.5 -0.75
8.4 Scattering about a vertical cylinder in-0.75nite-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 aect 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 ow. We seek to model the scat-0.5 tering of an incident wave eld which is propagating toward a bottom-mounted rigid cylin-0.75 der positioned in the middle of a nite-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 nite-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 whats 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 7
Lets summarize the smooth case Gibbs Phenomenon Filtering and post-processing Limiting TVD Runge-Kutta Methods A few theoretical results Solving the Euler equations
Lets summarize
We have achieved a good understanding
The theoretical support for DG for conservation
laws is very solid. The requirements for exact integration is expensive. It seems advantageous to consider a nodal approach in combination with dissipation. Dissipation can be implemented using a lter There is a complete error-theory for smooth problems. ... but we have forgotten the unpleasant issue What about discontinuous solutions?
associated with the appearance of a shock. In Fig. 5.4 we show the polynomial regarding accuracy andsimple discontinuity situations. representation of the stability in such
Before considering the full problem, let a Gibbs Phenomenon shock.us illustratewefundamental problem u(x) = associated with the appearancesign(x), x [1,Fig. 5.4 of a In 1], show the polynomial
and plots of the simple discontinuity representation of the pointwise error for increasing
Let us rst consider a simple approximation
u(x) = sign(x), x [1, 1],
0 h
number of terms, N , in the expansion. The results in Fig. 5.4 illustrate three unfortunate eects of the shock:
0 0.5 1 1.5 1 uh(x): N=32
uh(x): N=64
1.5 1 0.5 u(x)
|uh u|
1.5 and plots of the pointwise error for increasing numberN=8 terms, N , in the of 10 N=16 1 expansion. The results in(x): N=16 5.4 illustrate three unfortunate eects of the u Fig. N=128 N=32 N=64 u(x)=sign(x) shock: 0.5 10
1
u(x)
102
100 N=128 N=64 0.5 10
1
uh(x): N=16
0 x 0.5 1
N=8 N=16
0.5 1
0.5 u(x)=sign(x)
10
3
N=32
1
0 x
0 Fig. 5.4. On the left is shown the approximation of a sign function using a Legendre uh(x): N=64 expansion with increasing number of terms, illustrating the Gibbs phenomenon. On uh(x): N=32 0.5 the right is shown the pointwise error of the expansion. 102 1 1.5 1
Overshoot does not go away with N First order point wise accuracy Oscillations are global
0.5 0 x 0.5 1 Gibbs Phenomenon0.5 |uh u| 103 1
0 x
0.5
1
understand whether the combination of the Gibbs oscillation with propagation nderstand whether combination of thethe Gibbs oscillation with propagation tand whether is the use of the high-order basis on the elements that the the combination of Gibbs oscillation with propagation s. destroys the After all, it estroys solutionsolution globally rather reduces the the accuracyrstorderorder the solution globally rather rather reduces accuracy to to rst or or reduces the accuracy to rst order ys the globally scheme for smooth high-order accuracy of theora high-order accurate solution.discussed with hope of recovering a high-order accurate solutions, as no hope of recovering ith no solution. o hope of recovering a high-order accurate solution. (term N in Eq. (5.10)) ter 4. We write the problem on the skew-symmetric form Eq. We write the problem on the skew-symmetric form (term N1 in 1 (5.10)) write the doto the oscillationsof s rstBut problem on the skew-symmetric form (term N1 in Eq. (5.10)) return the simpler case destroy the nice behavior? 1 1 u u u u1 u u1 auau1 1 + 1 au + Lu a u x 2 x a 2 = u u 1 +u 2 a+ + u 1 x u ax u = + Lu = 0,= 0, t t x 2 a 2 t Lu = 0, 0, + a 2 + + a(x)x = x u = + t t 2 xt 2 x x 2 t ndand assume that a(x)smooth. Both a(x) and and u(x, t) are considered periassume that a(x) is is smooth. Both a(x) u(x, t) are considered perid odic that a(x) However, we assume that the the initial condition, u(x,is in some simplicity. However, we a(x) and u(x, t) condition, stasume simplicity. previously. For assume we initialestablished u(x, 0), for is smooth - but u(x,0) is not dic for a(x)detail is smooth. Boththis casethat have are considered peri-0), is ossibly with the introducelter, depending oninitial condition, u(x, 0), is r nonsmooth to use of a Gibbs phenomenon. simplicity. However, we assume phenomenon. onsmooth to introduce the the Gibbsthat the the details of the ux tion. to introduce the Gibbs phenomenon. ooth Introduce theadjoint problem Introduce the the adjoint problem Dene adjoint problem the lter shortly, let us attempt to e we returnadjoint problem of roduce the to the impact v v nd whether the combination of the Gibbsoscillation with propagation L L v = 0, v= v t reduces 0, accuracy to rst order the solution globally or rather t the L v = 0, hope of recovering = high-order. accurate solution. where (Lu, = a smooth v(x,0) the adjoint problem we we assume (u, Lt solving v) here (Lu, v) v) (u, L v) . For For solving the adjoint problem assume solved withthe skew-symmetric form (term N in Eq. (5.10)) rite the problemconditions. 1 smooth initial on mooth initial conditions.
Gibbs Phenomenon
(Lu, v) = (u, L v) . For solving the adjoint problem we assume Clearly, u 1 au immediately immediately have that h We We conditions. have that1 initialu 1 we have u + a immediately2have + 2 x 2 ax u = t + Lu = 0, t d x that d (u, = 0 (u(t), v(t)) = = (u(0), v(0)) (u, v) v) 0 = (u(t), v(t)) (u(0), v(0)) . . me that d dtis dt Both a(x) and u(x, are considered peria(x) smooth. t) simplicity. (u, v) = we (u(t), v(t)) initial condition, u(x, 0), is However, 0 assume that the = (u(0), v(0)) . Assuming integration of of all terms central uxes (i.e., no no aliasing ssuming dt exact integrationall terms and and central uxes (i.e.,aliasing or or exact th to introduce the Gibbs phenomenon.
to Galerkin orthogonality. of lastdual can be controlled and v(0) smooth, term inally, since the approximationThethe of theproblem is stableby the Cauchy-smooth, Finally, since (0), v (0)) smoothnessv(0)) + C(u)hNis stable and v(0) Schwarz inequality approximation v(0) problem (uh the and the (u(0), of dual to obtain +1 N q |v(0)|,q . h ,h we have we have hN +1hN +1 N q |v(0)| . (uh (0),the v (t) (u(0), v(0)) + C(u)hN +1 is; stable and v(0) smooth, ,h v(t) approximation ,h C(t) |v(t)||v(t)| Finally, since vh (0))v(t) ,h of the dual problem ,q ;,q h vh (t) C(t) ,q Nq Nq we have
shock 5and 138 Now, dt the solution cannot be better than rst order. We therefore assume Nonlinear problems consider that uh (0) is understood to be the polynomial interpolation (uh exact (0)),h = of all v(0)) and central uxes (i.e., noof the projection Assuming (0), vhintegration (u(0),terms + (uh (0) u(0), vh (0)),h aliasing or at the grid points, case, there is no quantization error caused by the grid. of u(0). In that we have forever lost the information of the location of the dissipation), centralvuxes, = (u(0), v(0))rst(uh (0) Weu(0), scheme wethe(0), similar be also for + h solution cannot statementhave semidiscrete vh (0)),h Usingand(uobtain ha(0)),h webetter than the order. v(0)),hassume shock Now, consider + (u(0), vh (0) therefore . that uh (0) is understood to be the polynomial interpolation of the projection (u (t), v (t)),h = (uh (0), v+(0)) vh (0) v(0)) . h (u(0), . of u(0). In (u (0), h (0))h is = (u(0), v(0)) error caused by thevgrid.,h that case, there no quantization + (u ,h u(0), (0)) vh irst, one realizes that the second term on the right-hand side hvanishes due h ,h h (0) ,h Now, consider proceeding, exactly what the Cauchyo BeforeFirst, one realizes that theasecond term on the right-hand uh (0)vanishes due Galerkin orthogonality. The last term can be controlled byside means, Consider we need to be bit careful about(u(0), vh (0) v(0)),h . given that u(0) is assumed to be discontinuous. Ifcan simply read this function to Galerkin and (0)) smoothness of term +we u(0), v (0)) by the Cauchyorthogonality. The last + (u to be controlled chwarz inequality(0), vh the,h = (u(0), v(0))v(0) h (0)obtain h ,h (uh Schwarz inequality and the smoothness of v(0) to obtain First, one realizes that the second term on the right-hand side vanishes due + (u(0), v (0) v(0)),h . tohGalerkin orthogonality. The last+ C(u)hNh+1 N q |v(0)|,q .the Cauchyterm can be N +1 q by controlled (u (0), vh (0)),h (u(0), v(0)) N h (0)) Schwarz(uh (0), vthat the second term on +right-hand side |v(0)|,q . inequality ,h (u(0), v(0)) of C(u)h First, also have and the smoothness the v(0) to obtain vanishes due one realizes We
d Now,of u(0). grid points,= 0 haveis no quantization error causedthe location of the consider (u, v)case, there forever lost the information of by the grid. at the In that we (u(t), v(t)) = (u(0), v(0)) .
h
Gibbs Phenomenon
Finally, since the approximation of the dual problem is stable and v(0) smooth, N +1 hat is, we have safely exchange v (t)v . This. results in in ; can that is, we can safelyv(t) v for vh for hC(t) h results exchange v This we |v(t)|,q h ,h q hN +1 N v(t) vh (5.18) (usafelyv(t))(t) ,hv C(t) v(t))|v(t)|,q ;+ , (5.18) that is, we can h (t), (uh (t), v(t)),h =h(u(t), results in exchange (u(t), .N q v(t)), ,h = for v This +
Combining it all, we obtain
(uh depends only on the smoothness where where small and depends v(t)),h = (u(t), v(t)) + , v(x,v(x, t). (5.18) is very is very small and(t), only on the smoothness of of t).
that is, we can safely exchange v for vh . This results in
5.6 Gibbs Phenomenon Problems with discontinuous solutions
13
order accurate in a pointwise sense, in agreement with the simpler analysis o The solution This has been demonstrated the linear problem.is spectrally accurate ! also for much more comple ... [98, 99, hidden problemsbut it is100]. 5.6.1 Filtering This also shows
that the high-order accuracy is maintained -- the oscillations are not noise to In light of the above, it is reasonable to consider ways ! recover some of th Recall
accuracy hidden in the oscillatory solutions. As it turns out, we already hav How do we disposal. such a tool at ourrecover the accurate solution? Consider
Np
uh (x) =
n=1
un Pn1 (x), un =
1 1
u(x)Pn1 (x) dx.
If u(x) has a discontinuity in [1, 1], the analysis in Section 4.5 shows that
One easily shows that
1 q u u(x) H nun . n n Thus, a manifestation of the Gibbs phenomenon, or rather lack of regularit
q
un . Thus, a manifestation of the Gibbs phenomenon, or rather lack of regularity, n is a slow decay of the expansion coecients. This suggests that we could 1 manifestation of the Gibbs phenomenon, or rather lack of regularity, . un attempt to modify the expansion coecients to decay faster in the hope of n w decay of the expansion coecients. This suggests that we could So there is a close connection between smoothness recovering a more rapidly convergent expansion and, thus, a more accurate to modify the expansionacoecients to decay faster Gibbs hope of and rapidlyThus, the expansion and,of thea in theaccurate decayconvergent expansion coefcients.phenomenon, or rat for manifestation thus, more approximation. ng a more is a slow decay of the expansion coecients. This su mation. In Fig. 5.5 we illustrate the impact of using the exponential lter
Filtering
n If u(x) has a discontinuity in [1, 1], the analysis in Sec 1
un
.
() = exp( s to obtain the ltered expansion ) approximation. n the ltered expansion Fig. 5.5 we illustrate the impact of using the expo Np Consider In n1 F uh (x) = un Pn1 (x). Np N n 1 n=1 F () = exp( s ) uh (x) = un Pn1 (x). N n=1 In Section 5.3 this was used to recover stability by controlling the impact of to obtain the ltered expansion on 5.3the aliasing errors. this was used to recover stability by controlling the impact of Example We sing errors. consider the sequence of functions Np consider the sequence of functions F x n1 uh 0 (i1) un Pn1 (x) cos(x), 1 x (x) = (0) (i) u = u N (s) ds, u = x cos(x), 1 cos(x), 0(i) x 1, (i1) n=1 1 x0 < u = u (s) ds, u(0) = cos(x), 0 < x 1, In Section 5.3 this 1 used to H which is constructed such that u(p)was p [1, 1]. recover stability by cont s constructed such that u(p) theppointwise error associated with three dierent lter In Fig. 5.5 we show H [1, 1].
we impact modify exponential lter ig. 5.5Perhaps attempt convince the expansion coecients to we illustrate thecan to of using thethe expansion do () = exp( s ) recovering a more rapidly convergent expansion decay faster ?
decay and, t
Filtering
140
0
5 Nonlinear problems
100 102 104 106 8 10 1010 12 10 1014 1016 1 10 102 104 106 8 10 1010 12 10 1014 1016 1 100 102 104 106 108 1010 1012 1014 1016 1 100 102 104 106 8 10 1010 12 10 1014 16 10 1
0
10 2 10 104 6 10 108 1010 1012 14 10 1016 1 100 2 10 104 106 108 10 10 12 10 1014 1016 1 100 102 104 106 108 1010 1012 1014 1016 1 100 102 104 106 8 10 1010 12 10 1014 16 10 1
0.5
0 x
0.5
1
0.5
0 x
0.5
1
100 102 4 10 6 10 8 10 1010 1012 14 10 1016 1 100 2 10 104 106 108 10 10 12 10 1014 1016 1 100 102 104 106 108 1010 1012 1014 1016 1 100 102 104 106 8 10 1010 12 10 1014 16 10 1
u(0)
0.5 0 x 0.5 1
u
0.5 0 x 0.5 1
(1)
0.5
0 x
0.5
1
0.5
0 x
0.5
1
u(2)
0.5 0 x 0.5 1
0.5
0 x
0.5
1
0.5
0 x
0.5
1
u(3)
0.5 0 x 0.5 1
0.5
0 x
0.5
1
0.5
0 x
0.5
1
Filtering
This achieves exactly what we hoped for
Improves the accuracy away from the problem spot Does not destroy matter at the problem spot
... but does not help there. This suggests a strategy:
Use a lter to stabilize the scheme but do not
remove the oscillations. Postprocess the data after the end of the computation.
the is, compute with the oscillations that clearly are not little) that computation with the lter as much (or rather as noise, at is, compute with the oscillations that clearly are not noise, tprocess the solution when needed. To illustrate the prospects of u solution when needed. To illustrate the prospects of u2 ocess the + h, let us consider an = 0, x [1, 1], example.
et us consider an example.
Filtering t x
We solve Burgers equation
WeConsider equation uous solve Burgers Burgers initial condition
equation
142
2.5
5 Nonlinear problems
100 102
ontinuous initial condition
u u2 + 2 = 0, 2, x 0.5 x [1, 1], u0 (x) ut u0) = 0, x [1, 1], = u(x, x = +
2
104 106 |uuh| 0.5 0 x 0.5 1
1.5
108 1010
t
x
1, x > 0.5.
1
1012 1014
tinuous initial condition u0 (x) = u(x, 0) =
-Hugoniot conditions, we easily0.5 that the shock propa2, x nd u0 (x) = stant speed of u(x, 0) = is,1, xx 0.5 solution is 3; that 2, the exact > 0.5.
2.5 2
0.5 1
1016 1 0.5 0 x 0.5 1
100 102 104 106 108 1010 1012 1014 1016 1 100 102 104 106 108 1010 1012 1014 1016 1
ers equation is solved using a discontinuous Galerkin method on to dene in appropriate boundary conditions. with aliasing thethe computation with ux and 8 order polyno0seequidistant elements, each of the an N =we add a lter equation is solved to connect the elements. To stabilize bility. ax-Friedrichs ux using a discontinuous Galerkin method on the K = 20 equidistantcomputation ofwith an Nand8we add a lter h aliasing in the elements, with = ux and N = polyno- Eq. se an exponential lter each the 36 = order 0 (see c da ity. Lax-Friedrichs ux to connect the elements. To stabilize the
1.5
ine-Hugoniot appropriate easily nd conditions. to dene theconditions, we boundary that the shock propaconstant speed u(x,using is, the exact solution is of quation is solved3;t)that 0a discontinuous Galerkin method on = u (x 3t), Overltering leads to aliasing dene the appropriate boundaryux and we add a lter o use to in the computation of the conditions.
0.5 1 2.5 0.5 0 x 0.5 1
ankine-Hugoniot conditions, we 3t), nd that the shock propau(x, t) = u0 (x easily he constant speed of 3; that is, the exact solution is
1.5 1
1, x > 0.5.
|uuh|
0.5
0 x
0.5
1
every stage of the RK method.
Limited ltering looks much better
1
0.5 1
0.5
0 x
0.5
1
|uuh|
severe smearing.3t), u(x, t) = u0 (x
2
0.5
0 x
0.5
1
similar spirit, one can lter in physical space, often known d achieve comparable results [138, 298]. Fig. 5.7. On the left is shown the purely polynomial solution of Burgers equa with N = 256 and the right shows the Pad-Legendre reconstructed solution e can be obtained (i.e., exponential convergence and L = 8.from M = 20 away ntinuity), by reexpanding the computed solutions using An forms). In this Pade seeks two 10 e.g., Pad alternative -case, oneltering local polye 10 QL of order M and L, respectively, dened on Dk , such
Filtering
0.2 1 0.8 0.6 0.4 0.2
0
0.2 0.4 0.6 0.8
1
0.2 1 0.8 0.6 0.4 0.2
0
0.2 0.4 0.6 0.8
0
2
104
xD :
k
uk (x) h
RM (x) , = QL (x)
106 108
10
10 + 1. There are several ways of dening in which sense the 144 5 Nonlinear problems the same, with the most obvious one being 10
12
1.2
1.2
[0, M + L] :
0.8 0.6 0.4 0.2
1
a much more rapid convergence even close to the point of
To fully recover, the ns, one can use a quadrature to evaluate the inner product. shock location is of the Pad reconstruction approach over straightforward e required (see text). tic reduction of the Gibbs oscillations, a reduced smearing
0.6 0.4 0.2
D
k
uk QL h
RM
0.8
1
0.8 0.6 0.4 0.2 0 0.2 0.4 0.6 0.8 1 m dx = 0. P Fig. 5.8. Pointwise error for the reconstructed solution to Burgers equation.
1014 1
use N = 256 and M = 20 for the numerator and the three curves represent, f top to bottom, L = 0, L = 4, and L = 8, respectively.
Eliminates oscillations and improves accuracy .. but no improvement at the point olve Burgers equation
102
As evidence of the possibility of dramatically increasing the pointwise a racy away from the shock, we show in Fig. 5.8 the pointwise error for M = 0 0 and for increasing values of L (i.e., by increasing the order of denominato 0.2 0.2 the 0.8 1 1 0.8 0.6 0.4 0.2 0 0.2 0.4 0.6 0.8 1 1 0.8 0.6 0.4 0.2 0 0.2 0.4 0.6rational approximation). This clearly conrms the enhanced converge in the piecewise smooth regions away from the point of discontinuity. Fig. 5.7. On the left is shown the purely polynomial solution of Burgers equation with N = 256 and the right shows the Pad-Legendre reconstructed solution with e These Pad-based expansion techniques are discussed at length in [ e M = 20 and L = 8. 160, 161] for Legendre-Pad-reconstruction, relevant to the present discuss e Whereas the above techniques oer signicant improvements, they unable to recover spectral accuracy all the way up to and including 100
Limiting
So for some/many problems, we could simply leave the oscillations -- and then postprocess. However, for some applications (.. and advisors) this is not acceptable
Unphysical values (negative densities) Articial events (think combustion) Visually displeasing (.. for the advisor).
So we are looking for a way to completely remove the oscillations: Limiting
dt solution is a temdt s present problems for many others; forxexample, if the 2 |u | elds ure or a density, the oscillations may cause these = to u . unphysical take u + f(u ) (5.19 2 t x integrate value of . Now, multiply Eq. (5.19) examples(ux ))x a slightly higher over the domain to Furthermore, there could be by ( where and x obtain would change uniqueness, we must also ensure uniform boundedness to Apart variablefrom solution the solution dramatically (e.g., in a combustion ss guarantee existence of the solution; that is, where things could ignite articially). We d interested in guaranteeing uniform boundedness are he ultimate solution1 would ( (ux ))x fxif onedx = develop a process dx. naturally be (u ) could ( (ux ))x u ux L + xx dt the oscillations entirely, without adversely impacting the would eliminate |u| dx. u L1 C, u L1 = acy of the computation. To understand what is required to achieve this, We realize that low the discussion in [60] loosely and return to the continuous viscosity Consider ion (Nonlinearproblems = 2 (u )u f (u )u dx = 0, (u))x fx (u )dx x xx x 146 5 x (u) u + f (u ) = 2 u . and dene (5.19) = |u| t x x from solution uniqueness, we must alsohave to Eq. (5.19), we dene (u) = |u|. since To see that this holds forwe ensure uniform boundedness to (u)u 0. Furthermore, solutions We have solution; that boundaries, we note that Assuming ntee existence of thesimple periodic is, d ( (ux ))xx u = = ux uxt = )(u|u)2|dx = d ux L1 . xx x dx 0, dx = dx (ux u ( (ux )) utu xx1dx |u| dx. L dt dt L1 C, |ux |
Limiting
establishing uniform boundedness ( (u ))x and integrate over the domain to Now, multiply Eq. (5.19) and one easily provesby x
obtain d ux dt We realize that
L1
+
d u x L1 0. ( (ux ))x fx (u ) dx = dt
To translate this into the discrete case, let us rst consider the semidiscrete
( (u ))x u dx. x xx
Limiting
We would like to repeat this for the discrete scheme. Consider rst the N=0 FV scheme
establishing uniform boundedness d establishing uniform boundedness L1 0. ux dt d d ux L1 0. To translate this into the discrete dt ux let us 0. consider the semidiscrete case, L1 rst dt caseTo translate this into the discrete case, let nite volume scheme). In this with a constant basis (i.e., a rst-order us rst consider the semidiscrete case, we with a constant basis (i.e., acase, case have the into the discrete rst-order rst consider the semidiscrete To translate this one-dimensional method let usnite volume scheme). In this casecase, we have the one-dimensional method with a constantk basis (i.e., a rst-order nite volume scheme). In this duh + k k+1 hone-dimensional method k , uk1 ) = 0, (5.20) case, we have the dt duk f (uh , uh ) f (uh h h h + f (uk , uk+1 ) f (uk , uk1 ) = 0, (5.20) h h dt K cells h hlocated equidistantly with a spacing of duk where we assume that the are h h + f (uk , uk+1 ) f (uk , uk1 ) = 0, h h h h and f (a, b) is a monotone ux. h located equidistantly with a spacing(5.20) where we assume that the K cells are of dt Toand f (a, b) continuous case, we dene h mimic the is a monotone ux. where we assume that the K cells are dene equidistantly with a spacing of To mimic the continuous case, we located ux. h h and f (a, b) is a monotone k+1 uk uh uk uk1 1 k h h vh = k+1 k k uh dene uh uk1 To mimic the continuous case, we uh h h h k h 1 vh = h h h k+1 u uk uk k1 1 as a rst-order approximation toh( (ux ))x . MultiplyinguEq. (5.20) by this k h h h . Multiplying Eq. (5.20) by this = vhterms yields ( (ux ))x andas a rst-order approximation to h summing all h h and summing all terms yields K d approximation k k+1 )) . Multiplying Eq. (5.20) by this k as a rst-order h |T V + K |ud vh f to h , uhwith f (uk , uk1 ) = 0, 147 (u ( (ux discontinuous ) 5.6 Problems k+1 x h kh solutions k k dt all |uh |T V + v and summing terms yields h f (uh , uh ) f (uh , uk1 ) = 0, k=1 h dt K k=1 K where the discrete total |uh |T V = norm is uk |. as variation |uk+1 dened h h ddiscrete total variation norm k+1 where the |u | f + v k k=1(uk , u is dened (uk , uk1 ) = 0, ) f as
Multiply with
and sum over all elements to get
dt
h TV
h
h
h
h
h
1 k |uh f (b) f (a) = (f (a) + Ca),|T V = = |uh (b) hCb). (f u |. 2 2 and, thus, uniform boundedness. Similar results can be obtained with other k=1 monotone nondecreasing due to monotonicity, one easily shows e f + and f are uxes. Using To go term this is monotone, consider the two key questions Consider thatbeyond uxsimple case, we must one easily proves the the
e
Limiting
+
d |uh |KV 0, T 1k+1 dt
wherequickly becomes very dux |into 0, we have split thetechnical. increasing and decreasing components; th |uh T V To appreciate thedt challenges is, for the Lax-Friedrichs ux, in this, let us consider a high-order scheme but require that the local means, or cell averages, uk are uniformly bounded. h thus, So In thisboundedness. Similar results-- but what about N>0 uniform N=0we consider for case, everything f (a) + f (b) beCobtained with other is ne can + n (a b), f (a, b) = otone uxes. 2 2 dk uh we k k+1 consider the two key questions must k k1 To go beyond this simpleh case, f (ur , ul ) f (ul , ur ) = 0, + dt then hat happens when a higher-order basis is used and how does the time where we have introduced the f + (a) + f k and uk = 1 left and right u ration come a Forward (a, b) = method of +ltime,rwethe into play.k f Euler notation in (b), n as get using5 Nonlinear uh, respectively. Iffwe (a) +rst-order forward Euler method to use a 148 limit value of problems f (b), n = that Postponingintegrate in time, we havefor a while, we rst recall 1, a higherthe latter question r basis may h uk,n+1 oscillationsk,n , uk+1,n ) fviolate k1,n )bound on the where introduce uk,n + f (ur and, thus, (uk,n , urthe = 0, l l t severely complicates the analysis and1 complete analysis 1 + variation. This f (a) = (f (a) + Ca), f (b) = a(f (b) Cb). 2 and with a monotone ux, 2 kly becomes very technical. one can show that [60] To appreciate the challengesn+1 | this,|n | us+due to monotonicity, one easily sho let Since f + and f are nondecreasing consider a high-order scheme | in u u = 0,
TV TV
k vh (f + (uk ) come+intok1 ) + (uk+1uk1 ) (uk )) 0, integration f (uk ,huk+1 ) f f (uk , ) f h k f (u play. h hh vh h h h Postponing the latter question for while, a we rst recall that a higherk onsideringorder few possible(f + (uk ) f + (uHence, f violate )the f (uk )), the the basis may introduce oscillations k1 ) thus, (uk+1 bound on = vh combinations. h and, + we obtain h h h
of what happens when a higher-order basis is used and how does the time
0
and therefore severely complicates the analysis and a complete analysis total variation. This
k
h uk,n+1 uk,n + f (uk,n , uk+1,n ) f (uk,n , uk1,n ) = 0, r r l l t
Limiting
with a monotone ux, one can show that [60]
Resulting in
|n+1 |T V |n |T V + = 0, u u
re
=
However, the monotone ux is not enough to K guarantee uniform k+1/2,n boundedness through 0 k+1/2,n k+1,n k,n
( u
k=1
t + h
That is the job of the limiter -- which must satisfy K
(k1/2,n ) (k+1/2,n ) u u f + (uk,n ) f + (uk1,n ) r r
) ( u
)
p(u
) p(u
)
k=1 Ensures uniform boundedness/control oscillations K Does notk+1/2,n conservation k+1,n violate t ( u ) (k1/2,n ) f (ul u ) f (uk,n ) . l hDoes not change the formal/high-order accuracy k=1
e (u) = |u|,
(k+1/2,n ) = u
This turns out to be hard ! uk+1,n uk,n
h
,
cannot be avoided. Consider an approximation, uh , to a smooth solution, u(x). If u(x) locally is monotone, then uh will behave locally like a straight line and it is easy to see that Eqs. (5.21)-(5.23) are trivially obeyed. However, if u(x) has a local smooth extrema, then the conditions are not necessarily fullled. Two tasks solution is In this case, theat hand wrongfully identied as having an oscillation and the limiter will alter the local solution and reduce the formal accuracy to rst order. This is a general property of slope limiters if the solution is required Detect Naturally, cells to be TVDM.troubledfor the simple example above, there is no need to use Limit the slope to a similar eect will occur in any computation of a slope limiter. However, eliminate oscillations problems with shocks and regions of smooth behavior. To nd a practical way of modifying the local mean of the solution to Dene the minmod function guarantee the TVDM property, we dene the minmod function
Limiting
m(a1 , . . . , am ) =
s min1im |ai |, |s| = 1 0, otherwise,
1 s= m
m
sign(ai ).
i=1
(5.24) To appreciate what this function does, consider a situation where we have If arguments, (a1 , a the Then m(a1 a2 , a3 ) will threea are slopes, 2 , a3 ).minmod , function return a zero unless the
Returns the minimum slope is all have the same sign Returns slope zero if the slopes are different
piecewise linear solution; that return; Inserting these into Eqs. (5.21)-(5.23) is, reveal that they suce to guarantee the k h xk )(uk x , uh (x) = uk + (x small )timestep, t. TVDM property of the solution for suciently 0 h This canwhere xnow dene dene families of slopek .limiters with slope limited be explored to the center coordinate of We can k represents the interface uxes as D We dene the slightly dif0 solution ferent properties. As a v k = step,m(k uk , uk uk1 , uk+1 uk ), rst uk assume that the solution isrepresented(5.25) by a hin which lcase k+1 solution is uh h the h h l h Let linear solution; that is, k piecewiseus assume N=1 k + m(uk uk ,kuk hk1 ,ukk+1k ukk1 u h uh uh 1 k u u ). (5.26) v k =ku
Limiting
where We solution
which can be shown to satisfy Eqs. (5.21)-(5.23). A slightly more dissipak TVDM property the classic MUSCL (Monotone Upstream-centered Scheme for center coordinate of Dk . small timestep, t. x0 represents theof the solution for sucientlyWe dene the slope limited tive limiter is
= uk + (x reveal that , h uk (x) (5.21)-(5.23) xk )(uh )xthey suce to guarantee the 0 Inserting these into h Eqs.
uh (x) r uh + (x x0r = h )m (uh )h , x h
h
h/2 k
h,
h/2
h
,
(5.27)
where k represents the center coordinate of which can sligthly less satisfy Eqs. (5.21)-(5.23).. A slightly more dissipa0 or a be xshown to dissipative limiter D We dene the slope limited solution tive limiter is the classic MUSCL (Monotone Upstream-centered Scheme for Conservation Laws) limiter [218, 301]. k uk+1 uk uk uk1 h h h h 1 k k k
k
, , uh (x) = hu+ (x= uh 0 )m x0 )m x ,(uh )x , u h (x) x+ (x (uh ) h h h x0 h uk (x) = uk + (x h/2k )(uk )x , h/2 h
ferent properties. As a rst step, assume that the solution is represented by a uk uk+1 k uuk uk k1 k+1 h uh h h uh k1 k uh piecewise linear solution; that k k is, 1 k k 1 k k k h h
have can beclassiclimiter [218, families of slope limiters with slightly difThis the exploredMUSCL limiter Conservation Laws) to dene 301].
,,
(5.28) (5.27)
tive limiter is the classic MUSCL (Monotone Upstream-centered Scheme for There are many other [218, 301]. they are similar types but Conservation Laws) limiter uk+1 uk uk uk1 h h h h 1 k k k k , uh (x) = uh + (x x0 )m (uh )x , , (5.28)
uk1 h 1 k k k k , uh (x) = uh +shown xto)m (uh )x , (5.21)-(5.23). A slightly more dissipa, (5.28) which can be (x 0 satisfy Eqs. h h
uh (x) = uh + (x x0 )m (uh )x ,
uk+1 h
h/2
,
uk h
uk h
h/2
,
(5.27)
t change the formal accuracy of the method.
t two are easy to satisfy, the third property causes problems. To s consider a simple example.
Limiting
Consider We solve the simple problem
u u + = 0, x [1, 1], t x
1.5 1 0.5 u(x,t) 0
5.6 Problems with discontinuous solutions
149 Smooth initial condition
0.5 1 1.5 1
0.5
0 x
0.5
1
Reduction to 1st order at local smooth extrema
Fig. 5.9. Solution to linear wave equation using N = 1 and K = 50 and a simple slope limiter. The dashed line is the exact solution at T = 10 and the solid line is the computed solution.
Limiting
Introduce the TVB minmod
to address this is to relax the condition on decay of the total variation require that the total variation of the mean is just bounded, called the TV condition. Following [285], this can be achieved by slightly modifying denition of the minmod function, m(), as
m(a1 , . . . , am ) = m a1 , a2 + M h2 sign(a2 ), . . . , am + M h2 sign(am ) , (5
where M is a constant that should be an upper bound on the second deriva 154 5 M estimates maximum curvature Nonlinear problems at the local extrema. This is naturally not easy to estimate a priori. Too s 1.5 a value of M implies higher local dissipation and order reduction, whereas high a value of M 1 reintroduces the oscillations.
0.5 u(x,t)
minmodB.m
function mfunc = minmodB(v,M,h) 0
0.5 % function mfunc = minmodB(v,M,h) % Purpose: Implement the TVB modified midmod function. v is a vector 1
mfunc = v(1,:); 1.5 1 ids = find(abs(mfunc) > 0.5 M*h.^2);0 x
0.5
1
Fig. 5.10. Solution to the linear wave equation using N = 1 and K = 50 and a if(size(ids,2)>0) with M = 20. The dashed line is the exact solution at t = 10 TVBM slope limiter
the is, compute with the oscillations that clearly are not little) that computation with the lter as much (or rather as noise, at is, compute with the oscillations that clearly are not noise, tprocess the solution when needed. To illustrate the prospects of u solution when needed. To illustrate the prospects of u2 ocess the + h, let us consider an = 0, x [1, 1], example.
et us consider an example.
Limiting t x
5.6 Problems with discontinuous solutions
100 102 104 106 |uuh| 108 1010 1012 1014 1016 1 100 102 104 106 108 1010 1012 1014 1016 1 100 102 104 106
8 10 12 14 16
155
We solve Burgers equation
2
WeConsider equation uous solve Burgers Burgers initial condition
equation
2.5
2
ontinuous initial condition
u u + 2 = 0, 2, x 0.5 x [1, 1], u0 (x) ut u0) = 0, x [1, 1], = u(x, x = +
1.5
t
x
1, x > 0.5.
1
K=20
0.5 0 x 0.5 1
tinuous initial condition u0 (x) = u(x, 0) =
0.5 1 2.5
0.5
0 x
0.5
1
ers equation is solved using a discontinuous Galerkin method on 1.5 10 to dene in appropriate boundary conditions. with aliasing thethe computation with ux and 8 order polyno- 10 0seequidistant elements, each of the an N =we add a lter 1 equation is solved using a discontinuous Galerkin method on the 10 bility. ax-Friedrichs ux to connect the elements. To stabilize 10 0.5 K = 20 equidistantcomputation ofwith an Nand8we add 0 a lter 1 10 h aliasing in the elements, with = ux and N0.5 = polyno- Eq. 0.5 se an exponential lter each the 36 = order 0 (see c x da ity. Lax-Friedrichs ux to connect the elements. To stabilize the
ine-Hugoniot appropriate easily nd conditions. to dene theconditions, we boundary that the shock propaconstant speed u(x,using 0 (x the exact of quation is dissipativealimiting solution is 3t), =u Too solved3;t)that is,discontinuous Galerkin method on aliasing dene the appropriate boundaryux and we add a lter o use to in the computation of the conditions.
0.5 1 2.5 0.5 0 x 0.5 1
ankine-Hugoniot conditions, we 3t), nd that 1the shock propau(x, t) = u0 (x easily he constant speed of 3; that is, the exact solution is
1, x > 0.5.
1.5
K=100
|uuh|
-Hugoniot conditions, we easily0.5 that the shock propa2, x nd u0 (x) = stant speed of u(x, 0) = is,1, xx 0.5 solution is 3; that 2, the exact > 0.5.
2
0.5
0 x
0.5
1
.. but no oscillations!
K=100
|uuh|
leads to severe smearing. u(x, t) = u0 (x 3t),
2
1
0.5
every stage of the RK method.
0 x
0.5
1
Fig. 5.11. Solution of Burgers equation with a shock using limiting. In all examples,
Limiting
But what about N>1?
Compare limited and nonlimited interface values If equal, no limiting is needed. If different, reduce to N=1 and apply slope limiting
156 5 Nonlinear problems
2.5 10
0
102 2 10 10
|uuh|
4 6 8 10 12
1.5
10 10
1
10
1014 0.5 1 0.5 0 x 0.5 1 10
16
1
0.5
0 x
0.5
1
Fig. 5.12. Solution of Burgers equation with a shock using limiting. We use K = 20 and N = 8 order elements and the generalized slope limiter, N . The results are shown at T = 0.4, with the dashed curve being the exact solution and the right gure showing the pointwise error.
Limiting
General remarks on limiting
The development of a limiting technique that avoid
local reduction to 1st order accuracy is likely the most important outstanding problem in DG all have some limitations -- restricted to simple/ equidistant grids, not TVD/TVB etc challenging
There are a number of techniques around but they The extensions to 2D/3D and general grids are very
es that maintain the TVDM or TVBM property, provided thi estion becomes whether one can design high-order temporal integratio for TVD Runge-Kutta Euler method. hemesthe rst-orderthe TVDMmethods that maintain forward or TVBM property, provided this can b nsider the semidiscrete scheme own for the rst-order forward Euler method.
Consider the semidiscrete scheme Consider again the semi-discrete dt
ssume that we can establish the TV property using a forwar d assume which we just discussed TV propertyschemes forward Eul that Foris, we can establish the TVD/TVB using a as d; that is, thod; that
n n+1 (uh , tn h n+1 |uh |uh un+1 = uun+ tLh (unntn ), ),|un+1 |T V|T V |un |T Vn |T V . = hh + tL h h , . u h h h
scheme d d uh = Lh (uh , t), dtuh = Lh (uh , t),
t consider an explicit RK method with s stages of the form form us consider an explicit RK method with s stages of the .. but this is just 1st order in time -- we want (0) =nun (0) v = uhigh-order accuracy h v h i1 i = 1, . . . , s : (i) = i1 ij v (j) + ij tLh (v (j) , tn + j t) . (5.3 v (i) j=0 i 1, . . . , s(s) v = j=0 ij v (j) + ij tLh (v (j) , tn + j t) . = n+1 : uh v we n+1 Do =(s) have to redo it all ?
uh
=v
early, we must nd (ij , ij , j ) such that the order conditions (see, e.g , we must nd ( , and if ) such that the order conditions (s y,143, 144]) are satised, , additional degrees of freedom are availab
TVD Runge-Kutta+methods |u =u tL (u , t ), u
n+1 h n h h n h n
and assume that we can establish the TV property using a forward Eule method; that is,
n+1 |T V h
Assume we canan explicitERKmethod with on theofform Let us consider nd a RK method s stages the form
(0) n v = uh i = 1, . . . , s : v (i) = n+1 uh = v (s)
i1 j=0
|un |T V . h
ij v (j) + ij tLh (v (j) , tn + j t) .
(5.30
such that the order conditions Clearly, we Coefcients must ndto ij , ij , j )additionalconditions found ( satisfy order degrees of freedom are (see, e.g [40, 143, 144]) are satised, and if available
i1
we can attempt to optimize the scheme in some way. For consistency, we mus 158 5 as Write this Nonlinear problems have
i1
v (i) =
j=0
ij v (j) +
ij il = 1. tLh (v (j) , tn + j t) . ij l=0
A Ifcloser look is TVD/TVB for t RK method in Eq. (5.30) reveals that at the form of the using an Euler method, we can directly the scheme E positive, 0 )that allij , high order also, provided simply maximum timestep (ij , onif are result atij >the RK method is we use a a convex combination o Clearlyij rely forward Euler steps since we can write the stages as
tRK min
ij
ij tE . ij
When optimizing the scheme, the objective should be to maximize the fraction in front of tE to minimize the cost of the time integration.
TVD Runge-Kutta+methods |u =u tL (u , t ), u
n+1 h n h h n h n
and assume that we can establish the TV property using a forward Eule method; that is,
n+1 |T V h
Assume we canan explicitERKmethod with on theofform Let us consider nd a RK method s stages the form
(0) n v = uh i = 1, . . . , s : v (i) = n+1 uh = v (s)
i1 j=0
|un |T V . h
ij v (j) + ij tLh (v (j) , tn + j t) .
(5.30
such that the order conditions Clearly, we Coefcients must ndto ij , ij , j )additionalconditions found ( satisfy order degrees of freedom are (see, e.g [40, 143, 144]) are satised, and if available
i1
we can attempt to optimize the scheme in some way. For consistency, we mus 158 5 as Write this Nonlinear problems have
i1
v (i) =
j=0
ij v (j) +
ij il = 1. tLh (v (j) , tn + j t) . ij l=0
The scheme is a convex combination of Euler steps ij and When stabilitythe scheme, the objective should be to maximize the fraction the optimizing of the high-order methods follows
tRK min in front of tE to minimize the cost of the time integration.
A Ifcloser look is TVD/TVB for t RK method in Eq. (5.30) reveals that at the form of the using an Euler method, we can directly the scheme E positive, 0 )that allij , high order also, provided simply maximum timestep (ij , onif are result atij >the RK method is we use a a convex combination o Clearlyij rely forward Euler steps since we can write the stages as
ij tE . ij
TVD Runge-Kutta methods
... but do such schemes exits ?
h h
Methods of this kind are known as strong stability-preserving Runge-Kutta in(SSP-RK) or TVD-RK methods and can be used with advantage for problems front of tE to minimize the cost of the time integration. Methods of this kind are known as strong they guarantee that no additional with strong shocks and discontinuities, asstability-preserving Runge-Kutta (SSP-RK) or TVD-RK methods partcanthe used with advantage for problems oscillations are introduced as and of be time-integration process. with strongsecond-order discontinuities, as they guarantee that no additional For a shocks and two-stage SSP-RK scheme, the optimal scheme is [139, oscillations are introduced as part of the time-integration process. 140] For a second-order two-stage SSP-RK scheme, the optimal scheme is [139, 140] (5.31) v (1) = un + tLh (un , tn ),
2nd order
1 n t v (1)n+1u= +(2) = h (un , nn+ v (1) + tLh (v (1) , tn + t) , (5.31) uh= h v tL huh ), 5.7 Strong stability-preserving Runge-Kutta metho 1 2 n n+1 (2) uh + v (1) + tLh (v (1) , tn + t) , uh = v = and the optimal third-order three-stage SSP-RK scheme is given as 2
and the optimal third-order three-stage SSP-RK scheme is given as
(1)
un+1 = v (5) = 0.00683325884039un + 0.51723167208978v (2) h h
(3)
3rd order
No 4th order, 4 stage scheme is possible - but there are other options (not implicit) With lter/limiting
(2)
2 Both schemesmethods 3are known (e.g., low-storage forms and is themultistep sche are optimal in the sense that the maximum timestep SSP same as schemes are optimalEuler method. the maximum forward in the sense that Boththat of therefer to [140] for an overview of these timestep is and same references methods the more Unfortunately, one can show possible to construct as that of the forwardaEuler method.[139] that it is not be applied at ...