The Legendre–Burgers equation: When artificial dissipation fails

Applied Mathematics and Computation 217 (2010) 1949–1964

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

The Legendre–Burgers equation: When artificial dissipation fails John P. Boyd Department of Atmospheric, Oceanic and Space Science, University of Michigan, 2455 Hayward Avenue, Ann Arbor, MI 48109, United States

a r t i c l e

i n f o

Keywords: Shock Burgers equation Legendre polynomials Legendre viscosity Spectral viscosity Hyperviscosity

a b s t r a c t Artificial viscosity is a common device for stabilizing flows with shocks and fronts. The computational diffusion smears the frontal zone over a small distance l where l is chosen so that the discretization has a couple of grid points in the front, and thus is able to resolve the shock. Spectral element methods use a Legendre spectral viscosity whose effect is to damp the coefficient of Pn(x) by some amount that depends only on the degree n of the Legendre polynomial. Legendre viscosity is better than ordinary diffusion because it does not require spurious boundary conditions, does not increase the temporal stiffness of the differential equations, and can be applied locally on an element-by-element basis. Unfortunately, Legendre diffusion is equivalent to a diffusion with a spatially-varying coefficient that goes to zero at the boundaries. Using the simplest example, one in which the second derivative of Burgers equation is replaced by the Legendre operator to give the ‘‘Legendre– Burgers” equation, ut + uux = m[(1  x2)ux]x, we show that the width of the computational front can similarly tend to zero at the endpoints, causing a numerical catastrophe. Ó 2010 Elsevier Inc. All rights reserved.

1. Introduction 1.1. Observed disasters in spectral element models with spectral vanishing viscosity Spectral element methods, employing Legendre polynomial expansions within each subdomain, have been quite successful in fluid mechanics. However, a common problem in all numerical methods for flows of high Reynolds number are (i) the spurious build-up of energy near the aliasing limit (variously called ‘‘spectral blocking” or ‘‘appearance of the [dreaded] 2-h waves”) and (ii) the spontaneous formation of shocks, which are called ‘‘fronts” in geophysics. The usual strategy for managing these difficulties is to add artificial, scale-selective damping, an idea that goes back to von Neuman and Richtmyer more than half a century ago [59,45]. With spectral elements, the most popular implementation is usually called ‘‘spectral vanishing viscosity” or SVV for short. In a Legendre–Galerkin discretization, the damping is a diagonal matrix which multiplies the nth Legendre coefficient an by a function q(n); in the absence of other physics, the damping would cause an to decay as an = an(0)exp(q(n)t). The ‘‘spectral vanishing” adjective means that the damping q(n) is chosen to vanish as n ? 0 so that the artificial damping does not affect the lowest Legendre coefficients at all (except indirectly). Such a ‘‘spectral vanishing” damping does not destroy the spectral accuracy of spectral elements for a smooth, shock-free solution, does not require spurious additional boundary conditions, and, with suitable postprocessing, allows recovery of spectral accuracy away from the shock even for flows with fronts. Although this artifice has many antecedents, ‘‘spectral vanishing viscosity” (SVV) was coined by Tadmor, who provided the first rigorous proof of its effectiveness for Fourier methods in 1989 [54]. Maday et al. [36] subsequently extended the idea to Legendre polynomial methods including spectral elements.

E-mail address: [email protected] 0096-3003/$ - see front matter Ó 2010 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2010.06.051

1950

J.P. Boyd / Applied Mathematics and Computation 217 (2010) 1949–1964

Articles on the theory and application of SVV with a Fourier basis include [37,55,51,15,56,57,29,21,24,38]. Legendre single-domain and spectral element applications and theory include [36,6,40,27,21,25,35,2,14,42,43,28]. Chebyshev polynomials or the Chebyshev eigenoperator have been combined with SVV in [1,33,34,18,3,47–49,53,42,43,50]. SVV has even been applied with hierarchal finite elements [12] and spherical harmonics [20]. The sheer length of the list of references shows that spectral vanishing viscosity has been very popular. Unfortunately, experience with atmospheric models has shown that spectral viscosity can sometimes fail in the sense that it is insufficient to damp spurious oscillations. Already, spectral viscosity has been purged from a version of the spectral element dynamics model at the National Center for Atmospheric Research because unphysical precipitation developed at the edges of each subdomain, and this could be cured only by retreating to a standard diffusive damping (Thomas, private communication). The spectral element model at the Naval Postgraduate School, developed quite independently, has similar spurious precipitation (Giraldo, private communication). The false, patchy rain is a symptom of unphysical oscillations in the vertical velocity near the edges of elements, but why is it difficult to find a Legendre spectral viscosity that can suppress such spurious wiggles without simultaneously damping all components of the motion so strongly as to turn the free-flowing shell of air into a spherical annulus of honey? We propose a two-part strategy to address this issue. The second part, which is still in progress, is a mostly numerical study comparing different choices of q(n) and analyzing their limitations. Because the most popular and widely used forms of spectral viscosity are exponential q(n) which are not analytic but merely C1, theoretical analysis will necessarily be limited. Therefore, in this article, we have pursued a simple strategy, which is to thoroughly analyze what we shall dub the ‘‘Legendre–Burgers” equation:

  ut þ c0 ux þ uux ¼ m ½1  x2 ux x :

ð1Þ

This is similar to the widely-studied and well-understood Burgers equation, differing only in that the second derivative damping of Burgers equation is replaced by the singular differential operator whose eigenfunctions are the Legendre polynomials. The close similarity to Burgers’ equation will be exploited throughout this work to understand its Legendre generalization. The dissipation in the Legendre–Burgers equation is scale-selective, affecting the nth coefficient in proportion to the eigenvalue of the dissipative operator, n(n + 1). It also is diagonal in a Legendre–Galerkin discretization, and does not require additional, unphysical boundary conditions. The only difference from a typical spectral vanishing viscosity is that the usual SVV’s are even more scale-selective with q(n) equal to zero or at least exponentially small as n ? 0. However, to paraphrase the old ‘‘it’s not a bug, it’s a feature” joke about software, the modest scale-selectivity of the Legendre viscosity in (1) is actually a virtue: if this damping-of-all-Legendre-coefficients sometimes fails to prevent discontinuity formation, then it seems rather unlikely that the weaker damping of a standard SVV, which dissipates only the Legendre coefficients of high degree, will somehow succeed. 1.2. The worrying physics of vanishing diffusion In an earlier note [7], we pointed out that a dissipation which acts independently on each Legendre coefficient can be regarded as being a function of the eigenoperator for the Legendre polynomials. To be precise, if a function u(x) has the Legendre expansion

uðxÞ ¼

1 X

an P n ðxÞ

ð2Þ

n¼0

and the action of the damping operator D has the Legendre expansion

DuðxÞ ¼

1 X

an qðnÞPn ðxÞ;

ð3Þ

n¼0

then using the property

LPn ðxÞ ¼ nðn þ 1ÞPn ðxÞ;

ð4Þ

  d d ð1  x2 Þ dx dx

ð5Þ

where

L

is the ‘‘Legendre operator”, it follows through the usual theory of linear operators that

D¼q

hpffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi i  1  4L  1 =2 :

ð6Þ

In computations, one uses only the ‘‘spectral viscosity function” q(n), and not the rather complicated expression for D in terms of L. Nevertheless, (6) shows explicitly that the properties of a Legendre spectral viscosity follow from those of the Legendre operator.

J.P. Boyd / Applied Mathematics and Computation 217 (2010) 1949–1964

1951

The first worry is that the Legendre operator, when added as a damping term as in the Legendre–Burgers equation, is a spatially variable diffusion l(x) with the diffusion going to zero at both endpoints, x = ±1 as noted in [7]. This is alarming because it suggests, in agreement with the disastrous unphysical oscillations near x = ±1 in atmospheric spectral element codes, that the Legendre diffusion is not damping at all at the boundaries. Before investigating this worry further, it is must be noted that the singularity of the Legendre operator and the vanishing of diffusion at the boundaries has its virtues, too. Unlike a conventional diffusion, Legendre viscosities and spectral viscosities do not require additional spurious boundary conditions. For example, the axisymmetric component of the diffusion equation on the surface of a sphere is

ut ðx; tÞ ¼ m Lu;

L

  d d ½1  x2  dx dx

ð7Þ

The solution is

uðx; tÞ ¼

1 X

an ð0Þ expðnðn þ 1ÞtÞ Pn ðxÞ

ð8Þ

n¼0

Each Legendre component damps out independently of all the others, and no explicit boundary conditions are necessary or even possible. In contrast, it is well-known that the discretized second derivative operator is very ‘‘stiff” with eigenvalues growing as N4 where N is the Legendre truncation [8]. Because the computational diffusion coefficient is typically O(1/N2), the stiffness of the second derivative operator requires that explicit time-marching methods use a timestep which is O(1/N2). In contrast, the Legendre diffusion has a largest eigenvalue of mN(N + 1), so the Legendre diffusion is not stiff. 1.3. The ever-shrinking Taylor shock The vanishing-diffusion worry is amplified by the known physics of Burgers equation,

ut þ c0 ux þ uux ¼ l uxx

ð9Þ

where l is a constant. Burgers equation has the exact travelling shock solution, the so-called Taylor shock,

uðx; tÞ ¼ k 

  J J tanh ½x  ðk þ c0 Þt 2 4l

ð10Þ

where the constant J > 0 is the jump across the shock and k is a second arbitrary constant, this one with no restriction on sign. The Taylor shock plays a role much larger role in the theory of Burgers equation than its status as an exact solution would indicate. The reason is that one can show through the method of matched asymptotic expansions that if the viscosity m is small, all fronts in Burgers equation can be approximated by matching a Taylor shock of appropriate jump, width and location (as the ‘‘inner” approximation) to a solution of the inviscid Burgers equation, which is the ‘‘outer” approximation [5,4]. The damping in the Legendre–Burgers equation is spatially variable rather than constant. However, similar singular perturbation arguments, developed further below with both theoretical arguments and numerical illustrations, show that a slowly-propagating shock will adjust its frontal width to the local value l(x) of the viscosity: the Taylor shock will still be an accurate approximation. This is very alarming because the Taylor shock formula shows that the width of the frontal zone goes to zero when the strength of the viscosity goes to zero. However, the viscosity vanishes at both endpoints for the Legendre–Burgers equation. This suggests that Legendre diffusion and Legendre spectral vanishing viscosities might allow a flow to be singular at the boundaries, wrecking the computation. We shall show below that this disaster sometimes happens at least for the Legendre–Burgers equation. Similarly, one expects that if the spectral viscosity is low, Legendre–Burgers solutions will also have fronts of universal form, approximated by the Legendre–Burgers equivalent of the Taylor shock. (Fast-moving shocks may exit the domain before evolving to the form predicted by matched asymptotic expansions, but these are given their own analysis below.) Thus, our study will not attempt to exhaustively explore all possible initial conditions, but will instead focus on what happens to a single mature shock as it propagates. 1.4. Scientific rhetoric and its limitations A pure mathematician would rightly object that our analysis of the Legendre–Burgers equation cannot prove anything about the properties of spectral vanishing viscosities. We are indulging in scientific rhetoric rather than scientific proof. For this reason, this present article is necessarily incomplete, and we are working hard on a follow-up that will study various flavors of the ‘‘Burgers-SVV” equation. An infamous example of such proof-by-extension is an article by Marvin Minsky and Seymour Papert that rigorously proved that so-called ‘‘single-layer perceptrons” were useless. They argued that multi-layer perceptrons would fail, too, even though they had not actually analyzed such things. The field of neural networks went into a kind of Dark Ages for a decade,

1952

J.P. Boyd / Applied Mathematics and Computation 217 (2010) 1949–1964

only to be resurrected when it was shown that multi-layer perceptrons are a powerful computational tool. Neural networks today are flourishing because even though the Papert–Minsky theorem is completely true, the inference drawn from it is completely wrong. Still, our alarming results for the Legendre–Burgers equation are a strong and worrying motivation for further studies of spectral vanishing viscosity. The only reason that exponential q(n) are used in SVV is to improve the scale-selectivity and reduce dissipation of smooth parts of the flow. In some sense, a Legendre diffusion, which damps all Legendre-coefficients and not just those of large degree, should be at least as good at ‘‘singularity-prevention” and shock-smoothing as spectral vanishing viscosity. So, it is fair to say that our Legendre–Burgers equation results are a sound basis for a conjecture about the limitations of spectral vanishing viscosity even though this present article is far short of a proof. 2. Theory of the Legendre–Burgers equation 2.1. Initial and boundary conditions The Legendre–Burgers equation,

  ut þ c0 ux þ uux ¼ m ½1  x2 ux x ;

uðx; 0Þ ¼ Q ðxÞ x 2 ½1; 1

ð11Þ

is well-posed with the same conditions for non-zero m as for zero m: one must impose a boundary condition on the inflow boundaries only, that is,

uð1; tÞ ¼ uL ðtÞ iff ðuð1; tÞ þ c0 Þ > 0;

uð1; tÞ ¼ uR ðtÞ iff ðuð1; tÞ þ c0 Þ < 0:

ð12Þ

2.2. Invariance theorems The following theorem implies that to discuss travelling shocks, it is sufficient to fix the amplitude of the jump J across the front at an arbitrary value since solutions for other values of J may be obtained by a simple dilation. Theorem 1. If v(x, t) solves the Legendre–Burgers equation for c0 ¼ c00 and m = m0 , then

uðx; tÞ  k v ðx; ktÞ is a solution for c0 ¼

kc00

ð13Þ 0

and m = km for arbitrary values of the positive constant k.

Proof. Substitution of the assumed solution into the differential equation followed by cancellation of common factors of k. h The next theorem shows that there is no loss of generality in discussing only shocks which are antisymmetric about the value of x where u = 0: unsymmetric shocks can be obtained by adding the same constant to u(x, t) and (with a minus sign) to the speed c0 as follows. Theorem 2. If v(x, t) solves the Legendre–Burgers equation for c0 ¼ c00 , then

uðx; tÞ  k þ v ðx; tÞ

ð14Þ

is a solution for c0 ¼ c00  k. Proof. Elementary substitution. h 2.3. Stationary shock By direct substitution, one can verify the Legendre–Burgers equation has bounded solutions which are independent of time

u ¼ c0 



J J arctanhðxÞ  p ; tanh 2 4m

ð15Þ

where J and p are arbitrary constants of integration. Theorem 3 (Singularities of the stationary shock). The function

u ¼ c0 



J J tanh arctanhðxÞ ; 2 4m

ð16Þ

J.P. Boyd / Applied Mathematics and Computation 217 (2010) 1949–1964

1953

which is the special case p = 0 of the general stationary Legendre–Burgers solution, near x = 1 is

uð1  zÞ  c0 

J 2

   

J J 1  2½z=2J=ð4mÞ exp z þ 2½z=2J=ð2mÞ exp z 8m 4m

ð17Þ

and therefore has a high order fractional power branch at x = 1. Proof. Substitute the usual expansions

tanhðqÞ  1  2 expð2qÞ þ 2 expð4qÞ þ    ; 1 1 2 z þ  arctanhð1  zÞ  ð1=2Þ logðz=2Þ  z  4 16

ð18Þ ð19Þ

into the analytical form of the function. h Although the singularities of this stationary shock are very weak in the sense that all derivatives up to order J/(4m) are bounded at the endpoints, it is still alarming that the shock is singular at all. 2.4. Nonsingularity of fast-moving shocks A simple argument shows that when the shock is moving sufficiently fast, it will not be singular at the endpoints as it exits. Suppose a travelling front has a width which is O(m) around x = 0 (where the viscosity is a maximum). The most extreme reduction in viscosity is an instantaneous jump to zero, which reduces the Legendre–Burgers equation to the onedimensional advection equation,

ut þ c0 ux þ uux ¼ 0;

uðx; 0Þ ¼ Q ðxÞ:

ð20Þ

This has the exact solution by the method of characteristics [60]

uðx; tÞ ¼ Q ðn½x; tÞ;

x ¼ n þ ½c0 þ Q ðnÞt:

ð21Þ

Unfortunately, this is only an implicit solution. However, the wave ‘‘breaks” at a time tB when characteristics first cross and ux is unbounded for some x. An explicit u(x, t) is not necessary to determine the time of breaking, which is

tB ¼ 

1 minx Q x ðxÞ

ð22Þ

provided that the minimum slope is negative. Unfortunately, if the flow has already evolved to a shock with a frontal zone of width O(m), then in the complete absence of viscosity, u(x, t) will become singular in a very short O(m) time. In particular, if u(x, 0) is a Taylor shock centered on the middle of the interval with a width controlled by the viscosity at x = 0,

Q ðxÞ ¼ u0 



J J tanh x ; 2 4m

ð23Þ

then

tB ¼

8m J2

ð24Þ

:

The good news is that if

c0 þ u0 >

J2 ; 8m

ð25Þ

then the shock will exit the domain before breaking. The bad news is that if m is very small, as true in most applications, then the shock must be moving very rapidly to satisfy this criterion. Spectral element domain decomposition helps because if the domain is divided into M elements, then the time for the shock to traverse a single element is reduced by 1/M. For fixed viscosity m, a shock travelling at any fixed speed c0 can be desingularized, guaranteed, by employing a sufficiently large number of elements M. The shock will zip across an element too fast to break. Instead, as it propagates through rapidly-varying viscosity, the shock will assume a frontal width determined by the average viscosity in an element. One complication is that when M is increased with the polynomial degree in each element fixed, the viscosity coefficient is usually reduced to reduce the width of the shock zone to exploit the increase in spatial resolution. For the Burgers and Legendre–Burgers equations, keeping a fixed number of grid points across the shock requires choosing the viscosity coefficient to be the inverse of the average grid spacing h, which implies m  O(1/M). Thus, increasing the number of elements does not help, for these equations, unless the viscosity is fixed.

1954

J.P. Boyd / Applied Mathematics and Computation 217 (2010) 1949–1964

2.5. Moving coordinate system and the narrow shock approximation The speed of a Taylor shock is independent of the viscosity. If

u0 ¼ ð1=2Þðuþ þ u Þ;

ð26Þ

where u+ and u are the values of u(x, t) on either side of the frontal zone (more precisely, the asymptotic limits of u as we move away from the front), then the front will travel at the speed c0 + u0. It is convenient to render the front stationary by shifting to the moving coordinate

X  x  ct  d;

T  t;

ð27Þ

where c and d are arbitrary. The convenient choices are (i) c = c0 + u0 to match the speed of the shock, and (ii) d is the location of the front at t = 0 so that the front will be centered on X = 0 for all time. It is convenient to redefine the unknown so that the shock is (roughly) antisymmetric about the middle of the frontal zone by setting

wðx; tÞ  uðx; tÞ  u0 :

ð28Þ

Using the identities

@ @ @ ¼ c ; @t @T @X

@ @ ¼ ; @x @X

ð29Þ

the Legendre–Burgers equation becomes

  wT þ wwX ¼ m ½1  ðX þ cTÞ2 wX : X

ð30Þ

In this reference frame, written in terms of the new unknown, the front is stationary and roughly antisymmetric in the sense that w(X, T)  w(X, T) where the only reason for the approximate equality is that the Taylor shock, which is strictly antisymmetric about X = 0, is only an inner approximation (in the language of matched asymptotic expansions) to a general w(X, T) which need not be antisymmetric. If the viscosity coefficient m is very small, which is the usual case since the Legendre viscosity is added only for computational stability, then the frontal zone will be very narrow, of width O(m). The first derivative in the viscosity operator is then O(m) smaller than the second derivative. Similarly dropping O(m) corrections in the coefficient of the second derivative terms gives, with relative error O(m),

wT þ wwX ¼ mð1  c2 T 2 ÞwXX ;

ð31Þ

which takes the form of a Burgers equations with a viscosity that depends upon time only. Burgers-with-time-varying-viscosity has been studied extensively in the acoustics literature by Crighton and co-workers [16,17,52,46,44]. Unfortunately, no useful analytic solutions are known. Burgers equation can be recast as a linear diffusion through the Hopf-Cole transformation, but no such linearizing transformation is known here. In earlier applications to acoustics, the damping strengthens with time, so the theory of Scott [52] cannot be directly applied to the present case where the damping is weakening. Thus, this reduction to a Burgers equation with time-dependent coefficient is mostly useful in a negative sense to show that previous studies of generalized Burgers equations are no help here. 2.6. Analytic approximations to travelling shocks If a shock is moving slowly and the width of the frontal zone is very narrow, which will be true everywhere if m  1, then it is at least plausible that the shock will adjust itself to the local value of the viscosity. This gives the approximation

u

  J J ðx  c0 tÞ ½Spatially-Varying Taylor Shock: tanh 2 4m 1  x2

ð32Þ

Substitution into the differential equation shows that the residual is O(m) smaller than the other terms in the Legendre–Burgers equation except where all those terms are exponentially small. More important, this remains true even as jxj ? 1. A second approximate of slightly different form can be derived by recalling that the Legendre–Burgers equation has a stationary (independent of time) solution. If we make the ansatz of allowing the phase shift constant in the stationary front to vary with time, then

uðx; tÞ  

  J J tanh farctanhðxÞ  arctanhðs½tÞg ; 2 4m

ð33Þ

where s(t) is the location of the center of the shock where u = 0. With this ansatz, the nonlinear and Legendre viscosity terms cancel exactly, and the residual is

R ¼ ut þ c0 ux ¼ 

   J2 J ds 2 2 2 ð1  x arctanhðxÞ  arctanhðs½tÞ Þ þ c ð1  s ðtÞÞ : sech f g 0 4m dt 8mð1  x2 Þð1  s2 ðtÞÞ

ð34Þ

J.P. Boyd / Applied Mathematics and Computation 217 (2010) 1949–1964

1955

Because the last factor is spatially-varying, it is not possible to make the residual zero everywhere, and thus generate an exact solution. However, the hyperbolic secant function ensures that the residual is exponentially small outside a narrow frontal zone centered at x = s(t). It is possible to make the residual zero at x = s(t), and therefore O(m) smaller in the frontal zone than the individual terms of the Legendre–Burgers equation, by choosing s(t) = c0t, which reduces the residual to

R ¼ c0

   2  J2 J 2 f arctanhðxÞ  arctanhðs½tÞ g x  c20 t2 ; sech 4m 8mð1  x2 Þð1  s2 ðtÞÞ

ð35Þ

where for simplicity the phase of the shock is chosen such that the shock is centered on x = 0 at t = 0. The residual of this more complicated approximation is simpler than for the Spatially-Varying Taylor Shock approximation, and is proportional to c0. This implies that for slowly-moving shocks, what we shall dub the ‘‘Arctanh Approximation” is better. However, when c0 is O(1), there is little to choose between them — both are O(m) approximations, uniformly accurate in time and space, though the Arctanh Approximation seems to have a little smaller maximum pointwise error. By using an elementary trigonometric identity, the Arctanh Approximation can be put in the form

uðx; tÞ  



 J J x  c0 t tanh arctanh ½Arctanh Approximation: 2 4m 1  xc0 t

ð36Þ

Numerical simulations in the next section show that these approximations are very accurate except perhaps very close to the exit boundary. 3. Numerical solutions of the Legendre–Burgers equation 3.1. Numerical algorithms The Legendre–Burgers equation was discretized by a Legendre polynomial collocation method. This reduces the partial differential equation to a system of ordinary differential equations in time of the form !

!

!

! !

! ! ! ! ! d~ u ¼ c0 D ~ u ~ u  D~ u  m M1 E M ~ u; dt

ð37Þ ! !

where ~ uðtÞ contains the values of u(x, t) on the Legendre–Lobatto grid, D is the usual Legendre differentiation matrix, the dot ! !

denotes element-wise multiplication in the nonlinear term, M is the matrix that converts grid point values to Legendre-coef! !

ficients, and E is a diagonal matrix with Ejj = (j + 1)(j + 2) with j = 1, . . ., N. The differential equations are advanced in time by fourth order Runge–Kutta time integration. To impose the boundary conditions, either or both of the components of the differential equation, those corresponding to the endpoints, u(±1, t), are dropped and replaced by a Dirichlet boundary condition. The inflow conditions are imposed only where c0 + u(±1) is directed into the domain. Thus, an inflow condition u(1, t) = uleft is imposed only when c0 + uleft > 0, and a right inflow condition u(1, t) = uright only when c0 + uright < 0. The Legendre–Lobatto grid points, differentiation matrix and grid-to-spectral conversion matrix are given in [8] and numerous other texts and so are omitted here. The Legendre discretization was chosen because any form of Legendre viscosity can be trivially imposed by changing the ! !

elements of the diagonal matrix E . However, the method is far from optimal because no Fast Fourier Transform is available to reduce the cost of evaluating the right-hand side of the ODE system. Being a firm believer in the adage ‘‘a good workstation never sleeps”, we valued ease of programming over efficiency, and let the production computations run while we slept the sleep of the just. 3.2. Canonical problem The full parameter space of the Legendre–Burgers equation is vast. However, for present purposes, the interesting solutions are fronts. When the viscosity coefficient m is small, the solution will steepen due to nonlinearity until it reaches an equilibrium between nonlinearity and the local viscosity which can be approximated, in the language of matched asymptotics, by the Taylor shock. Of course, if the wave moves through the domain very rapidly, it may exit before this steepening has ever evolved to a Taylor shock. However, we have already proved that sufficiently fast-moving waves never break. To focus closely on solutions that may become singular as they exit, the numerical results presented here are the solutions to a simple canonical problem which consists of an antisymmetric Taylor shock initially centered in the middle of the domain where the viscosity is strongest and varying most slowly. All velocity comes from the explicit propagation term c0ux; there is no loss of generality in this because an unsymmetric shock can always be converted into an equivalent antisymmetric Taylor shock by merely adding or subtracting the same constant from both u(x,t) and the propagation velocity as shown in Theorem 2.

1956

J.P. Boyd / Applied Mathematics and Computation 217 (2010) 1949–1964

This canonical initial condition defines a three-dimensional parameter space where the parameters are J, the strength of the jump across the shock, c0, the velocity of the shock, and m, the Legendre viscosity coefficient. However, it is easy to show that if we multiply all three of these parameters by an arbitrary constant k, then k can be divided out of the Legendre–Burgers equation if we dilate the time coordinate k by the same factor as expressed formally by Theorem 1 above. The parameter space is thus the projective plane in this three-dimensional space: it is sufficient to pick J = 2, as we shall do, and the dynamics in this plane will be the same (modulo a dilation in time) as for any other J = constant plane. Furthermore, it is sufficient to consider positive c0 only since flows with negative c0 are just the mirror image of flows with positive c0. The parameter plane falls into three regions. Specializing (25) to u0 = 0 [shocks with a zero mean value] and J = 2 [shocks asymptoting to ±1 at the ends of the frontal zone, we have theoretically proved that u(x, t) is nonsingular if

c0 >

1 2m

½Theoretical Nonsingularity Criterion:

ð38Þ

However, this argument was based on reducing the viscosity instantly to zero instead of the gradual reduction of the Legendre–Burgers equation, and therefore only furnishes a conservative lower bound. Our numerical calculations delineate where the non-singularity is not guaranteed by this theoretical criterion, but can be confidently asserted on the basis of self-consistent, convergent computations. The third region is the ‘‘zone of uncertainty” where the front steepens so much that numerical computations, with a reasonable number of grid points, are no longer reliable and we are left to wonder at each point in this region: Is the front truly singular, or merely very, very narrow as it exits the domain? 3.3. Verification of the analytic approximations Fig. 1 shows that both approximations are accurate for both small c0 and large c0. There is nothing special about our initial condition of a shock whose initial width matches that of the Taylor shock. Fig. 2 shows the error in the Arctanh Approximation when the initial condition was deliberately chosen to be double that which balances nonlinearity and dispersion. As expected, nonlinearity steepens the shock so that the flow spontaneously evolves to the form predicted by the Arctanh Approximation. However, these encouraging results cannot be the whole story because the analytic approximations always predict that u(x, t) will be discontinuous at x = 1 as the shock exits the domain, and the numerical solutions show that for fixed c0 and a sufficiently strong Legendre viscosity, the solutions are nonsingular. Fig. 3 compares the errors in the Arctanh Approximation for three different viscosities. The errors for small viscosities are flat until very near the exit time when the front becomes so narrow that it is underresolved, even with N = 800 Legendre polynomials. In contrast, the error for m = 1/50 (top curve, thin dashes) resembles a parabola rather than a straight line, and rises sharply to 2, the magnitude of the jump across the frontal zone. The reason for this parabola-and-rise-to-the jump error is that the flow for this (relatively) large viscosity is smooth whereas the Arctanh Approximation is discontinuous. For large viscosities — and m = 1/50 is not very large, but evidently large enough — errors can accumulate until higher order terms in m are no longer negligible. A better analytical approximation is highly desirable, but our efforts have so far been unsuccessful. However, the nearly-linear error curves for small m show no signs that the approximation is failing as the shock approaches the boundaries, narrowing and narrowing as it moves. If the Arctanh Approximation fails for m = 1/200, it must do so when the front has already become so narrow that even 1600 Legendre polynomials are insufficient to resolve it.

Loo error in analytic approximations

ν =0.0125 c0=0.125

ν =0.01 c0=1.5

0.05

0.02

0.04 0.015 0.03 0.01 0.02 0.005

0.01

0 0

0.2

0.4

0.6

time/exit time

0.8

1

0

0

0.2

0.4

0.6

0.8

1

time/exit time

Fig. 1. Maximum pointwise errors in the two analytic approximations for m = 1/80 and c0 = 1/8 [left] and m = 1/100 and c0 = 1.5 [right]. The errors at the initial time are zero because the approximations were used to generate the initial conditions, which centered the shocks at x = 0. The thick lower curve is the Arctanh Approximation.

1957

J.P. Boyd / Applied Mathematics and Computation 217 (2010) 1949–1964

Arctanh error, ν=0.005 N=1600 c0=1.5

0.5

Loo error

0.4 0.3 0.2 0.1 0 0.7

0.8

0.9

1

time/exit time Fig. 2. Maximum pointwise errors in the Arctanh Approximation from an initial shock centered at x = 0.7 with twice the width predicted by the approximation for that location. The long plateau in error is at roughly 0.0235. Although 1600 Legendre polynomials were used, the error jumps very close to the exit time because the front has become too narrow to be resolved even with this huge number of degrees of freedom.

Error in Arctanh Approx. vs Time, N=800

Error in Arctanh Approx. vs Time, N=800 0.2

0.2

ν=0.01 ν=0.01 0.15

ν=0.02

ν=0.02

Loo error

Loo error

ν=0.015

ν=0.015

0.15

0.1

0.1

0.05

0.05

0 0

0.2

0.4

0.6

0.8

time/exit time

1

0 0.9

0.92

0.94

0.96

0.98

1

time/exit time

Fig. 3. Maximum pointwise errors in the Arctanh Approximation for c0 = 1.5 and three different Legendre viscosities. The right panel is a zoom plot of the left.

3.4. Numerical results with one inflow boundary: c0 P 1 Fig. 4 shows a typical experiment in which a shock travelling at unit speed, c0 = 1, is computed as it moves from the center to the right boundary for eleven different viscosities. The initial condition of each case was the Taylor shock with its width adjusted to match the strength of the viscosity at the center of the interval. As m decreases, the spectrum of Legendre-coefficients versus degree becomes flatter and flatter (left panel). The right panel shows that even for provably nonsingular solutions, well-resolved by the basis of four hundred Legendre polynomials employed, the width of the front (inversely proportional to the maximum of ux in the front) narrows by a factor of thousands! Fig. 5 shows that the log-of-the-log of the slope amplification is so smooth it is well fit by a parabola, allowing extrapolation to m = 0. The extrapolated amplification rises sharply with the degree of the fit. When converted for the log-of-the-log to the slope amplification itself, the extrapolated values are 1053, 10250 and overflow. Polynomial extrapolation is fraught with peril. Nevertheless, the extrapolation clearly suggests that the Legendre–Burgers solution is singular with an infinitely narrow slope at exit when m = 0. There is no sign of singularity for finite, non-zero m. 3.5. Numerical results with two inflow boundaries: c0 < 1 When 0 < c0 < 1 and the initial condition asymptotes to u = 1 at x = 1 and u = 1 at x = 1, there is inflow from both the right and the left boundaries. If the right boundary condition is fixed at x = 1, then singularity is guaranteed. The initial

1958

J.P. Boyd / Applied Mathematics and Computation 217 (2010) 1949–1964

Legendre coefficients

slope amplification 104

100

ν=0.055 103 10-5

0 ν=0.1

102

0 .07

ν=0

5

10-10

0

200

101 0.05

400

0.1

ν

degree j

Fig. 4. Left panel: Legendre-coefficients at the time of maximum slope for each of eleven different Legendre viscosity coefficients, m from 0.055 to 0.10. c0 = 1, u(x, 0) = tanh (x/(2m)) with N = 400 polynomials used in the spectral computation. Right panel: the ratio of the maximum slope [very close to the time when the shock exits the domain] divided by the maximum of juxj at t = 0.

log10 (log10(slope amplification)) 5 4 3 2 1 0 -1

0

0.02

0.04

0.06

0.08

0.1

ν Fig. 5. The circles are the same data as the right panel of the previous graph. The dashed curves are the polynomial fits to the logarithm (base 10) of the logarithm of the slope amplification factor: quadratic polynomial fit [bottom, red], cubic [green] and quartic [top curve, magenta]. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

shock propagates propagates to the right at c0. When it reaches the right boundary, the ‘‘plateau” of u = 1 to the left of the front now fills the entire domain, so that the solution should be u(x, t)  1 for all t > 1/c0. However, the right boundary condition is fixed at u = 1, a legitimate choice. The shock is trying to exit to the right while the right inflow boundary condition is bringing in negative u from the right. The collision produces a jump discontinuity at x = 1 where the local Legendre viscosity, being zero at that point, can do nothing to intervene as illustrated in Fig. 6. Fig. 7 illustrates a typical case. The contour plot shows only one-fortieth of the full spatial domain, and a very short time interval to illustrate that, despite Legendre viscosity with the not-very-small coefficient of m = 1/10, the flow is well-resolved until the shocks collide in the upper right corner of the space–time plane.

J.P. Boyd / Applied Mathematics and Computation 217 (2010) 1949–1964

1959

u(-1,t)=1

u(1,t)= -1 x=1

x=-1

Fig. 6. When the shock is moving at a speed c0 < 1, there is inflow at the left boundary (where u(1, t) + c0 > 0) and at the right boundary (where u(1, t) + c0 < 0). The hollow arrow shows the rightward movement of the frontal zone; the thick black arrows depict the movements of the flat portions of the wave away from each of the inflow boundaries. When the shock reaches the right boundary, u(1, t) must jump discontinuously from 1 to the specified inflow boundary value, 1.

4. Regime diagram Fig. 8 summarizes the theoretical and numerical results above by displaying four different regimes in the m  c0 plane, that is, the plane in parameter space which is spanned by the Legendre viscosity coefficient m and the speed c0 of the shock. In the upper middle and right regimes, all shocks travelling sufficiently fast so that there is only one inflow boundary, the shock narrows greatly as it propagates towards a boundary, but never narrows to a width of zero. The numerical results described above suggest that for c0 > 1, there is singularity only when m = 0. However, the uncertainties of polynomial extrapolation are only part of the reason that the extrapolated small m is shaded by question marks. Another reason is that the maximum slope of the front rises exponentially fast — an exponential-of-an-exponential growth — that the front becomes so narrow as to be numerically unresolvable on an any reasonable grid for small m. When m < 1/50, exp(log(10)exp(log(10))) > 1 even with the most conservative extrapolation (by the quadratic polynomial in m, red dashes), which is equivalent to a slope amplification larger than 1010 and a frontal width of O(1010) of the domain width. When jc0j < 1, it is necessary to impose inflow boundary conditions on both boundaries to be consistent with our initial conditions, which are u(±1, t) = ±1. This leads to a shock collision at the right boundary (left boundary if c0 < 0), which is always singular. 5. The rate of convergence of Legendre series It was proved many decades ago that if a function u(x) is free of singularities everywhere within a region in the complexplane bounded by an ellipse with foci at x = ±1, but is singular on this ellipse, then the Legendre-coefficients an will decrease in proportion to exp (nl) where l is the quasi-radial coordinate in elliptic coordinates in the complex plane,

c0=1/2 N=200 1.0

1.95

1.9

t 1.85

0.8

0.6 0.4 0.2 0 -0.2 0 - .4 -0.6 -1.0

8

-0.

1.8

1.75 0.95

0.96

0.97

0.98

0.99

1

x Fig. 7. Zoom plot of the upper corner of the space–time domain, showing the convergence of the contour lines at the point x = 1, t = 2 where the shock tries to exit the domain at an amplitude of u = 1 while the boundary condition at x = 1 is fixed at u = 1. N = 200. m = 1/10, c0 = 1/2, u(x, 0) = tanh(5x). The boundary conditions are u(x = 1, t) = 1 and u(x = 1, t) = 1.

1960

J.P. Boyd / Applied Mathematics and Computation 217 (2010) 1949–1964

4

? ? ? ? ? ? ? ?? ?? ? ?

3.5 3

c0

2.5 2 1.5 1

Nonsingular numerically

Colliding shocks: inflow from both boundaries Singular

0.5 0

Nonsingular by theory

0

0.1

0.2

0.3

0.4

ν Fig. 8. Schematic showing the regimes of the problem of a shock propagating from the center of the domain to a boundary. In the upper right, the solution is provably nonsingular by a theoretical argument, which is that a fast-moving shock cannot break before exiting even if viscosity is ignored (c0 > (1/2m), Eq. (38)). The middle region is nonsingular because numerical solutions show that the solution is smooth even as it exits, as demonstrated by exponentiallydecaying Legendre-coefficients for u(x, t) for all t. The region shaded by question marks is uncertain: solutions may be singular, but the numerical solutions are underresolved even with 800 Legendre polynomials. The bottom region is where shocks collide at one boundary; the solution is always singular when this happens.

x ¼ coshðlÞ cosðgÞ;

l 2 ½0; 1; y ¼  sinhðlÞ sinðgÞ; g 2 ½0; 2p;

ð39Þ

If the convergence-limiting singularity is located at the point x0 + iy0, then

l ¼ Imfarccos½x0 þ iy0 g

ð40Þ

¼ log jz0 ðz20  1Þ1=2 j  pffiffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ log a þ a2  1 ;

ð41Þ ð42Þ

where the sign in the second line is chosen to make the argument of the logarithm > 1 so that l > 0 and where

a

1 2

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 ðx0 þ 1Þ2 þ y20 þ ðx0  1Þ2 þ y20 : 2

ð43Þ

The subtlety is that if the singularity is a distance  from the interval x 2 [  1, 1], then l  pif singularity is not too ffiffiffiffiffithe ffi close to the endpoints. However, if the singularity is very near the endpoints x = ±1, then l  2, which is a much faster rate of convergence [8]. In other words, Legendre series, like the Chebyshev series discussed in Chapter 2 of [8] are much more tolerant of singularities near the boundaries than near the center.

Legendre convergence for ν=1/10 0.07 0.06 0.05

μ

0.04 0.03 0.02 0.01 0

0

0.2

0.4

0.6

0.8

1

shock location s Fig. 9. The asymptotic rate of convergence l versus the shock location for m = 1/10. (The meaning of the ‘‘asymptotic rate of convergence” is that the Legendre-coefficients an fall proportional to exp(ln); When Legendre-coefficients are plotted on a log-linear scale versus degree, l is the slope of the envelope of the curve).

J.P. Boyd / Applied Mathematics and Computation 217 (2010) 1949–1964

1961

Heuristically, the same conclusion follows from examining the Legendre collocation grid, which has points separated by O(1/N) near the center of the interval but by only O(1/N2) near the boundary. Consequently, we cannot instantly conclude that merely because the shock is narrowing that the convergence of the Legendre series will necessarily decrease as the shock propagates toward the endpoints and simultaneously narrows. It is well-known that the function tanh(z) has poles at z = ±ip. The hyperbolic tangent is only an approximation to the Legendre–Burgers shock; nevertheless, it is instructive to examine the asymptotic rate of convergence l for the Spatially-Varying Taylor Shock function,

f ðxÞ  tanh



 xs : mð1  x2 Þ

ð44Þ

The convergence rate l is inversely proportional to m, so it suffices to plot l(s) for a singular arbitrary value of m as done in Fig. 9. Unfortunately, the increasing grid density near the endpoints is not sufficient to insulate the Legendre series from the pernicious effects of the narrowing of the shock. It is clear that as the front narrows, l falls rapidly. In terms of a log-linear plot of janj versus n where the an are the Legendre-coefficients, the slopes, which are equal to l, will becomes shallower and shallower as the front propagates closer and closer to the boundary at x = 1. Both Legendre convergence theory as well as the numerical results in the left panel of Fig. 4 show that the narrowing of the front greatly degrades the rate of convergence of the Legendre series, even for cases where the front is moving so rapidly that its width does not shrink all the way to zero. 6. Two philosophies of shock-capturing In meteorology, the strategy for capturing fronts via artificial damping is to to make the dissipation sufficiently strong so that the width of the frontal zone is three times the grid spacing. No further approximation is necessary because three grid points within the front are sufficient to resolve it. Of course, a viscosity sufficiently strong to smear the front this wide will suck energy out of even the largest length scales, converting a model of the troposphere into a model of the ‘‘honeysphere”, that is, a layer of gas as viscous as honey. To mitigate this, meteorologists and oceanographers long ago shifted to hyperviscosity: an artificial damping proportional to the pth power of the Laplace operator. A Fourier component exp(ikx) will be dissipated by a hyperviscosity mp@ 2p/@x2p as exp(mpk2pt). In order to damp the highest wavenumber resolvable on a grid of uniform spacing h, which is k = kalias = p/h, mp is chosen to be O(h2p) so that the waves of wavelength 2h are damped on an O(1/h) timescale where the timestep is typically O(1/h). The damping of the lowest wavenumber k is then O(exp (h2pt)), which is negligibly weak if p is large and h  1. Spectral viscosity, which replaces a damping that is a power of k by a function that is zero for small jkj, and grows exponentially for larger jkj, is similar in philosophy but even more scale-selective, preserving spectral accuracy in the smooth parts of the flow by damping the long waves not at all. The disadvantage of both hyperviscosity and spectral viscosity is that the fronts are no longer monotonic but oscillatory [6,26]. The philosophy of Maday et al. [36,40] differs dramatically from the meteorological perspective in that the spectral viscosity is used not to resolve the front, but only to stabilize the flow to force the discretized u(x, t) to converge to the correct solution, namely, one that respects the ‘‘entropy condition” as explained in these references. The issue of whether the front is resolved is not important; instead, a smooth, monotonic (or nearly monotonic) front is retrieved only by postprocessing using some kind of filter, [39,41,58], Padé reconstruction [19] and Gegenbauer or other polynomial reprojection [23,11,22]. To a meteorologist, this ‘‘stabilize-and-postprocess” strategy is completely impractical. The reason is that rain falls where the vertical velocity is upward. As air rises, the mean pressure falls, the parcel of air expands and it also cools adiabatically. The cooling reduces the amount of moisture that can remain as vapor. Condensation begins when the relative humidity reaches 100%. Thus, Gibbs oscillations in the vertical velocity can easily trigger spurious patches of precipitation that we may dub ‘‘underresolution rain”. In atmospheric general circulation models (GCMs), the code is divided into the ‘‘dynamical core” and the ‘‘physics modules”. The latter include not only parameterizations of clouds and precipitation, but also chemistry, radiative transfer, air-sea interaction and a host of other processes. The crucial point is that the output of the dynamics is, at every timestep, the input for the physics modules. Gibbs-corrupted dynamics will cause many of these physics modules, including precipitation, to fail catastrophically. These difficulties are not peculiar to geophysical fluid dynamics, but will arise in any discipline where fluid mechanics is combined with additional sensitive physical processes. For such ‘‘physics-rich” flows, postprocessing is not an operation that can be performed once at the end. Rather, the flow must be accurate at the end of each and every timestep. Either the postprocessing must be performed every timestep (expensive) or the front must be so wide that postprocessing is unnecessary. 7. Corner waves and all that: when dispersion and dissipation fail In the same way that Burgers equation is the most familiar dissipative regularization of the one-dimensional advection equation, the most familiar dispersive regularization is the KdV equation,

ut þ uux þ uxxx ¼ 0:

ð45Þ

1962

J.P. Boyd / Applied Mathematics and Computation 217 (2010) 1949–1964

Unfortunately, both of these regularizations are rather misleading. For both, the magnitude of the regularization term, either uxx [Burgers] or uxxx [KdV], increases without bound as the length scale goes to zero. A frontal zone can never shrink to zero and thus become singular because, no matter how strong the nonlinearity, the increasing narrowness of the frontal zone will amplify the regularizing term until it balances the nonlinearity, and then the front shrinks no further. Other choices of dispersion and dissipation allow waves to break if the nonlinearity is sufficiently strong. In this respect, the possibility of singular solutions of the Legendre–Burgers equation, and of similar equations obtained by a Legendre spectral vanishing viscosity, is not an isolated, surprising event. Rather, it is quite common. Examples of wave breaking and singularity formation due to insufficient dispersion are discussed in [32,9,10,13]. The existence of a ‘‘critical threshold” [in amplitude] for breaking in dissipative equations is analyzed in [30,31]. The model which is closest to the Legendre–Burgers equation is what we shall dub the ‘‘RCE–Burgers” equation, a Roseneau–Chapman–Enskog model, which is, defining F as the Fourier transform and UðkÞ  F ðuÞ as the Fourier transform of u,

!

2

ut þ uux ¼ m F

1

k

1þm

2 k2

Uðt; kÞ

 

1 ¼ m F 1 Uðt; kÞ : 2 1 þ m2 k xx

ð46Þ

The dissipation differs from that of the classical Burgers equation because of the denominator of 1 + m2k2 in wavenumber space; the rational function of wavenumber k ensures that the dissipation is bounded instead of unbounded, as true of a second derivative damping, as k ? 1. Liu and Tadmor [30] rigorously prove that some solutions of the RCE–Burgers equation develop infinite slope in finite time. The breaking-in-spite-of-dissipation behavior of the RCE–Burgers equation does not prove anything for the Legendre–Burgers equation. Liu and Tadmor’s theorem does falsify the proposition that damping always prohibits breaking in a generalized Burgers equation. For spatially periodic problems with dissipation operators of simple form in Fourier space, it is easy to see, as for the RCE– Burgers equation, that the dissipation operator is bounded, and therefore can be overwhelmed by the nonlinearity of a sufficiently large initial disturbance. But why is the Legendre diffusion insufficient when its effect on the nth Legendre polynomial grows as n2, without limit, as n ? 1? A partial answer is: truncation is not the same as dissipation. The most drastic form of a spectral vanishing viscosity, the most scale-selective, is a filter that is a step function in Legendre degree n — but this merely truncates the Legendre spectrum at a lower degree. Truncation by itself is hardly an adequate dissipation, or there would be no reason to bother with an explicit computational dissipation in the first place. Clearly, it is not enough to simply taper the tail of the Legendre series; one must examine exactly how such filtering alters the width of high gradient features of the function which is being filtered. 8. Summary Artificial diffusion that damps the coefficient of the nth Legendre polynomial by a damping function that depends only on n is very appealing for reasons noted in the abstract. By setting r(n/N) = 0 for sufficiently small n, one can preserve spectral accuracy in the smooth regions of the flow. Legendre damping can be applied on an element-by-element basis, strictly locally, and thus is ‘‘trivially parallel” in the language of multiprocessor computing. Legendre viscosities never require spurious, unphysical boundary conditions. Furthermore, if the strength of the spectral viscosity is chosen so that the highest retained Legendre-coefficients are damped on the same time scale as that of the inviscid dynamics, which is the obvious and conventional choice, then the addition of spectral viscosity does not require shortening the timestep below what is needed for stability and accuracy in timestepping the inviscid flow. In contrast, a conventional diffusion raises the order of the differential equation and, because the operator is nonsingular at the endpoints, requires additional boundary conditions. Insofar as the inviscid dynamics is concerned, these additional boundary conditions are spurious, an unfortunate complication. Spurious outflow boundary layers also result. Furthermore, the diffusion operator d2/dx2 must be discretized on a global basis that couples all the elements together. To pile bad news upon bad news, a diffusion with a coefficient of O(1/N2) where N is the truncation, chosen thus so that the shortest waves allowed on the grid (‘‘2-h waves”) are damped on an O(1) timescale, will require shortening the timestep to O(1/N2). The reason is that the diffusion operator has large spurious eigenvalues of size O(N4) [8]. Overall, it is far, far cheaper and easier to damp by means of a Legendre spectral viscosity, that is, an ? r(j/N)an for some filter function r with an as the coefficient of the nth Legendre polynomial, than it is add the usual Laplace operator, d2/dx2. It is therefore very disappointing that a Legendre viscosity sometimes fails to prevent shocks from breaking. Our answer is frustratingly ambiguous. Spectral viscosity has been successfully used in many applications; nothing here argues otherwise. Stationary fronts that remain always in the center of the domain or an element will not be singular. Fast-moving shocks will not break either, and instead have a width determined by the average viscosity in an element. It seems impossible to avoid the conclusion that the dissipative effectiveness of Legendre spectral viscosity really does drop to zero at the boundaries of an element or domain, and this is disastrous when slow-moving shocks exit. The ‘‘Gibbs rain” that occurs in spectral element climate models with hydrology, delineating the element walls in the precipitation graphs as spurious patches of rain, are different from under-resolved shock waves, but are likely due to the same cause as the underresolved fronts discussed here: at the boundaries, Legendre damping is a non-damping. The simplest remedy is to use a true diffusion instead of a Legendre diffusion, and this is what has been adopted for the NCAR spectral element models. Unfortunately, this is a very expensive option.

J.P. Boyd / Applied Mathematics and Computation 217 (2010) 1949–1964

1963

An alternative is to use what one might call ‘‘element-dithering”. When elements cover a spherical earth in latitude and longitude, one may, every few time steps, rotate the elements in the longitudinal direction by half the width of an element, then similarly rotate about an axis through the Greenwich meridian at the equator. Such half-an-element rotation implies that the weakly damped regions at the edges of the old elements are now strongly damped because they are at the center of the new elements. This is expensive for climate models because of the heavy bookkeeping required by all the physics and chemistry. However, it is only necessary to use the rotated elements for Legendre spectral viscosity. One could simply interpolate across element boundaries to apply the damping only while leaving the elements unchanged for all the other physics. To give an honest assessment of the cost and effectiveness of ‘‘element-dithering” requires an article in its own right and will not be attempted here. At present, we have no proven cheap and efficient remedy for spurious rain at element boundaries, only suggestions. But the first step in removing an obstacle is surely to identify it. Acknowledgements This work was supported by the National Science Foundation through Grant OCE0451951 and ATM 0723440. I thank Steve Thomas for reinvigorating my interest in this problem, Frank Giraldo and Mark Taylor for helpful conversations and the reviewers. References [1] Ø. Andreassen, I. Lie, C.E. Wasberg, The spectral viscosity method applied to simulation of waves in a stratified atmosphere, J. Comput. Phys. 110 (1994) 257–273. [2] H.M. Blackburn, S. Schmidt, Spectral element filtering techniques for large eddy simulation with dynamic estimation, J. Comput. Phys. 186 (2003) 610–629. [3] S. Bonazzola, E. Gourgoulhon, J.A. Marck, Spectral methods in general relativistic astrophysics, J. Comput. Appl. Math. 109 (1999) 433–473. [4] J.P. Boyd, The energy spectrum of fronts: the time evolution of shocks in Burgers’ equation, J. Atmos. Sci. 49 (1992) 128–139. [5] J.P. Boyd, Physical acoustics and the method of matched asymptotic expansions, Phys. Acoust. 11 (1992) 69–149. [6] J.P. Boyd, Hyperviscous shock layers and diffusion zones: monotonicity, spectral viscosity, and pseudospectral methods for high order differential equations, J. Sci. Comput. 9 (1994) 81–106. [7] J.P. Boyd, Two comments on filtering, J. Comput. Phys. 143 (1998) 283–288. [8] J.P. Boyd, Chebyshev and Fourier Spectral Methods, second ed., Dover, Mineola, New York, 2001. 665 pp.. [9] J.P. Boyd, A Legendre pseudospectral method for computing travelling waves with corners (slope discontinuities) in one space dimension with application to Whitham’s equation family, J. Comput. Phys. 189 (2003) 98–110. [10] J.P. Boyd, The cnoidal wave/corner wave/breaking wave scenario: a one-sided infinite-dimension bifurcation, Math. Comput. Simulat. 69 (2005) 235– 242. [11] J.P. Boyd, Trouble with Gegenbauer reconstruction for defeating Gibbs’ phenomenon: Runge phenomenon in the diagonal limit of Gegenbauer polynomial approximations, J. Comput. Phys. 204 (2005) 253–264. [12] M. Calhoun-Lopez, M.D. Gunzburger, A finite element multiresolution viscosity method for hyperbolic conservation laws, SIAM J. Numer. Anal. 43 (2005) 1988–2001. [13] R. Camassa, D.D. Holm, J.M. Hyman, A new integrable shallow water equation, in: T. Wu, J.W. Hutchinson (Eds.), Advances in Applied Mechanics, vol. 31, no. 31, Academic Press, New York, 1994, pp. 1–33. [14] M.H. Carpenter, On the conservation and convergence to weak solutions of global schemes, J. Sci. Comput. 18 (2003) 111–132. [15] G.Q. Chen, Q. Du, E. Tadmor, Spectral viscosity approximations to multidimensional scalar conservation laws, Math. Comput. 61 (1993) 629–643. [16] D.G. Crighton, Model equations of nonlinear acoustics, Ann. Rev. Fluid Mech. 11 (1979) 11–33. [17] D.G. Crighton, J.F. Scott, Asymptotic solutions of model equations in nonlinear acoustics, Philos. Trans. Roy. Soc. Lond. A 292 (1979) 101–134. [18] W.S. Don, D. Gottlieb, Spectral simulation of supersonic reactive flows, SIAM J. Numer. Anal. 35 (1998) 2370–2384. [19] T.A. Driscoll, B. Fornberg, A Padé-based algorithm for overcoming the Gibbs phenomenon, Numer. Algorithms 26 (2001) 77–92. [20] A. Gelb, J.P. Gleeson, Spectral viscosity for shallow water equations in spherical geometry, Mon. Weath. Rev. 129 (2001) 2346–2360. [21] A. Gelb, E. Tadmor, Enhanced spectral viscosity approximations for conservation laws, Appl. Numer. Math. 33 (2000) 3–21. [22] A. Gelb, J. Tanner, Robust reprojection methods for the resolution of the Gibbs phenomenon, Appl. Comput. Harm. Anal. 20 (2006) 3–25. [23] D. Gottlieb, C.-W. Shu, On the Gibbs phenomenon and its resolution, SIAM Rev. 39 (1997) 644–668. [24] J. Guermond, S. Prudhomme, Mathematical analysis of a spectral hyperviscosity LES model for the simulation of turbulent flows, ESAIM – Math. Model. Numer. Anal. 37 (2003) 893–908. [25] B. Guo, H. Ma, E. Tadmor, Spectral vanishing viscosity method for nonlinear conservation laws, SIAM J. Numer. Anal. 39 (2001) 1254–1268. [26] J. Jiménez, Hyperviscous vortices, J. Fluid Mech. 279 (1994) 169–176. [27] G. Karamanos, G.E. Karniadakis, A spectral vanishing viscosity method for large-eddy simulations, J. Comput. Phys. 163 (2000) 22–50. [28] R.M. Kirby, S.J. Sherwin, Stabilisation of spectral/hp element methods through spectral vanishing viscosity: application to fluid mechanics modelling, Comput. Methods Appl. Mech. Eng. 195 (2006) 3128–3144. [29] O. Lepsky, Spectral viscosity approximations to Hamilton–Jacobi solutions, SIAM J. Numer. Anal. 38 (2000) 1439–1453. [30] H.L. Liu, E. Tadmor, Critical thresholds in a convolution model for nonlinear conservation laws, SIAM J. Math. Anal. 33 (2001) 930–945. [31] H.L. Liu, E. Tadmor, Spectral dynamics of the velocity gradient field in restricted flows, Commun. Math. Phys. 228 (2002) 435–466. [32] H.L. Liu, E. Tadmor, Rotation prevents finite-time breakdown, Physica D 188 (2004) 262–276. [33] H.-P. Ma, Chebyshev–Legendre spectral viscosity method for nonlinear conservation laws, SIAM J. Numer. Anal. 35 (1998) 869–892. [34] H.-P. Ma, Chebyshev–Legendre super spectral viscosity method for nonlinear conservation laws, SIAM J. Numer. Anal. 35 (1998) 893–908. [35] X. Ma, V. Symeonidis, G.E. Karniadakis, A spectral vanishing viscosity method for stabilizing viscoelastic flows, J. Non-Newtonian Fluid Mech. 115 (2003) 125–155. [36] Y. Maday, S.M.O. Kaber, E. Tadmor, Legendre pseudospectral viscosity method for nonlinear conservation laws, SIAM J. Numer. Anal. 30 (1993) 321– 342. [37] Y. Maday, E. Tadmor, Analysis of the spectral viscosity method for periodic conservation laws, SIAM J. Numer. Anal. 26 (1989) 854–870. [38] A.A. Oberai, J. Wanderer, A dynamic multiscale viscosity method for the spectral approximation of conservation laws, Comput. Methods Appl. Mech. Eng. 195 (2006) 1778–1792. [39] S.M. Ould Kaber, Filtering non-periodic functions, in: C. Bernardi, Y. Maday (Eds.), Analysis, Algorithms and Applications of Spectral and High Order Methods for Partial Differential Equations, Selected Papers from the International Conference on Spectral and High Order Methods (ICOSAHOM ’92), Le Corum, Montpellier, France, 22–26 June 1992, North-Holland, Amsterdam, 1994, pp. 123–130.

1964

J.P. Boyd / Applied Mathematics and Computation 217 (2010) 1949–1964

[40] S.M. Ould Kaber, A Legendre pseudospectral viscosity method, J. Comput. Phys. 128 (1996) 165–180. [41] S.M. Ould Kaber, H. Vandeven, Reconstruction of a discontinuous function from its Legendre coefficients, C.R. Acad. Sci. – Math. 317 (1993) 527–532 (French with English abstract and summary). [42] R. Pasquetti, Spectral vanishing viscosity method for LES: sensitivity to the SVV control parameters, J. Turbulence 6 (2005), doi:10.1080/ 14685240500125476. [43] R. Pasquetti, Spectral vanishing viscosity method for large eddy simulation of turbulent flows, J. Sci. Comput. (2006). [44] C.S. Rao, P.L. Sachdev, M. Ramaswamy, Self-similar solutions of a generalized burgers equation with nonlinear damping, Nonlinear Anal. – Real World Appl. 4 (2003) 723–741. [45] R.D. Richtmyer, Difference Methods for Initial-Value Problems, Interscience Tracts in Pure and Applied Mathematics, vol. 4, Interscience/John Wiley, New York, 1957. [46] O.V. Rudenko, S.I. Soluyan, Theoretical Foundations of Nonlinear Acoustics (R.T. Beyer, Trans.), Consultants Bureau (Plenum), New York, 1977. [47] S.A. Sarra, Chebyshev super spectral viscosity method for a fluidized bed model, J. Comput. Phys. 186 (2003) 630–651. [48] S.A. Sarra, Spectral methods with postprocessing for numerical hyperbolic heat transfer, Numer. Heat Transfer A – Appl. 43 (2003) 717–730. [49] S.A. Sarra, The spectral signal processing suite, ACM Trans. Math. Software 29 (2003) 195–217. [50] S.A. Sarra, Digital total variation filtering as postprocessing for Chebyshev pseudospectral methods for conservation laws, Numer. Algorithms 41 (2006) 17–33. [51] S. Schochet, The rate of convergence of spectral-viscosity methods for periodic scalar conservation laws, SIAM J. Numer. Anal. 27 (1990) 1142–1159. [52] J.F. Scott, The long time asymptotics of solutions to the generalized Burgers equation, Proc. Roy. Soc. London A 373 (1981) 443–456. [53] S. Sirisup, G.E. Karniadakis, A spectral viscosity method for correcting the long-term behavior of pod models, J. Comput. Phys. 194 (2004) 92–116. [54] E. Tadmor, Convergence of spectral methods for nonlinear conservation laws, SIAM J. Numer. Anal. 26 (1989) 30–44. [55] E. Tadmor, Shock capturing by the spectral viscosity method, Comput. Methods Appl. Mech. Eng. 80 (1990) 197–208. [56] E. Tadmor, Superviscosity and spectral approximations of nonlinear conservation laws, in: M.J. Baines, K.W. Morton (Eds.), Numerical Methods for Fluid Dynamics IV, Oxford University Press, London, 1993, pp. 69–82. [57] E. Tadmor, Total variation and error estimates for spectral viscosity approximations, Math. Comput. 60 (1993) 245–256. [58] E. Tadmor, J. Tanner, Adaptive mollifiers for high resolution recovery of piecewise smooth data from its spectral information, Found. Comput. Math. 2 (2002) 155–189. [59] J. von Neumann, R.D. Richtmyer, A method for the numerical calculation of hydrodynamic shocks, J. Appl. Phys. 21 (1950) 231–237. [60] G.B. Whitham, Linear and Nonlinear Waves, John Wiley & Sons, New York, 1974.