1.3.COMPARISON WITH FINITE ELEMENT METHODS or finite difference methods for many classes of problems.However,they are most useful when the geometry of the problem is fairly smooth and regular. So-called"spectral element"methods gain the best of both worlds by hybridizing spec- tral and finite element methods.The domain is subdivided into elements,just as in finite elements,to gain the flexibility and matrix sparsity of finite elements.At the same time, the degree of the polynomial p in each subdomain is sufficiently high to retain the high accuracy and low storage of spectral methods.(Typically,p =6 to 8,but spectral element codes are almost always written so that p is an arbitrary,user-choosable parameter.) It turns out that most of the theory for spectral elements is the same as for global spec- tral methods,that is,algorithms in which a single expansion is used everywhere.Conse- quently,we shall concentrate on spectral methods in the early going.The final chapter will describe how to match expansions in multiple subdomains. Low order finite elements can be derived,justified and implemented without knowl- edge of Fourier or Chebyshev convergence theory.However,as the order is increased,it turns out that ad hoc schemes become increasingly ill-conditioned and ill-behaved.The only practical way to implement "nice"high order finite elements,where "high order" generally means sixth or higher order,is to use the technology of spectral methods. Similarly,it turns out that the easiest way to match spectral expansions across subdo- main walls is to use the variational formalism of finite elements. Thus,it really doesn't make much sense to ask:Are finite elements or spectral methods better?For sixth or higher order,they are essentially the same.The big issue is:Does one need high order,or is second or fourth order sufficient? h-refinement Smaller h increase subdivide only where polynomial high resolution degree p needed p-refinement r-refinement Figure 1.2:Schematic of three types of finite elements
1.3. COMPARISON WITH FINITE ELEMENT METHODS 5 or finite difference methods for many classes of problems. However, they are most useful when the geometry of the problem is fairly smooth and regular. So-called “spectral element” methods gain the best of both worlds by hybridizing spectral and finite element methods. The domain is subdivided into elements, just as in finite elements, to gain the flexibility and matrix sparsity of finite elements. At the same time, the degree of the polynomial p in each subdomain is sufficiently high to retain the high accuracy and low storage of spectral methods. (Typically, p = 6 to 8, but spectral element codes are almost always written so that p is an arbitrary, user-choosable parameter.) It turns out that most of the theory for spectral elements is the same as for global spectral methods, that is, algorithms in which a single expansion is used everywhere. Consequently, we shall concentrate on spectral methods in the early going. The final chapter will describe how to match expansions in multiple subdomains. Low order finite elements can be derived, justified and implemented without knowledge of Fourier or Chebyshev convergence theory. However, as the order is increased, it turns out that ad hoc schemes become increasingly ill-conditioned and ill-behaved. The only practical way to implement “nice” high order finite elements, where “high order” generally means sixth or higher order, is to use the technology of spectral methods. Similarly, it turns out that the easiest way to match spectral expansions across subdomain walls is to use the variational formalism of finite elements. Thus, it really doesn’t make much sense to ask: Are finite elements or spectral methods better? For sixth or higher order, they are essentially the same. The big issue is: Does one need high order, or is second or fourth order sufficient? Smaller h h-refinement h p-refinement increase polynomial degree p r-refinement subdivide only where high resolution needed Figure 1.2: Schematic of three types of finite elements
6 CHAPTER 1.INTRODUCTION 1.4 Comparisons with Finite Difference Method:Why Spec- tral Methods are Accurate and Memory-Minimizing Finite difference methods approximate the unknown u(r)by a sequence of overlapping polynomials which interpolate u(x)at a set of grid points.The derivative of the local interpolant is used to approximate the derivative of u(x).The result takes the form of a weighted sum of the values of u(z)at the interpolation points. Spectral One high-order polynomial for WHOLE domain Finite Difference Multiple Overlapping Low-Order Polynomials Finite Element/Spectral Element Non-Overlapping Polynomials, One per Subdomain Figure 1.3:Three types of numerical algorithms.The thin,slanting lines illustrate all the grid points(black circles)that directly affect the estimates of derivatives at the points shown above the lines by open circles.The thick black vertical lines in the bottom grid are the subdomain walls. The most accurate scheme is to center the interpolating polynomial on the grid point where the derivative is needed.Quadratic,three-point interpolation and quartic,five-point interpolation give df/dz [f(+h)-f(x-h)]/(2h)+O(h2) (1.12) d/dx≈[-f(x+2h)+8f(x+h)-8f(x-h)+f(x-2h)1/(12h)+O(h4) (1.13)
6 CHAPTER 1. INTRODUCTION 1.4 Comparisons with Finite Difference Method: Why Spectral Methods are Accurate and Memory-Minimizing Finite difference methods approximate the unknown u(x) by a sequence of overlapping polynomials which interpolate u(x) at a set of grid points. The derivative of the local interpolant is used to approximate the derivative of u(x). The result takes the form of a weighted sum of the values of u(x) at the interpolation points. Spectral Finite Element/Spectral Element One high-order polynomial for WHOLE domain Finite Difference Multiple Overlapping Low-Order Polynomials Non-Overlapping Polynomials, One per Subdomain Figure 1.3: Three types of numerical algorithms. The thin, slanting lines illustrate all the grid points (black circles) that directly affect the estimates of derivatives at the points shown above the lines by open circles. The thick black vertical lines in the bottom grid are the subdomain walls. The most accurate scheme is to center the interpolating polynomial on the grid point where the derivative is needed. Quadratic, three-point interpolation and quartic, five-point interpolation give df /dx ≈ [f(x + h) − f(x − h)]/(2h) + O(h2) (1.12) df /dx ≈ [−f(x + 2h)+8f(x + h) − 8f(x − h) + f(x − 2h)]/(12h) + O(h4) (1.13)
1.4.COMPARISONS WITH FINITE DIFFERENCES 16 2引 0. 0.4 0 -04 08 1.2 -1.6 T 27 为 Figure 1.4:Weights w;in the approximation df/drlz=xo≈∑，w;f(ro+jh)where zo=-π and h=/5.In each group,the Fourier weights are the open,leftmost bars.Middle,cross- hatched bars (j=+1,+2 only):Fourth-order differences.Rightmost,solid bars (j =+1 only):weights for second order differences. The function O(),the "Landau gauge symbol",denotes that in order-of-magnitude,the errors are proportional to h2 and h4,respectively. Since the pseudospectral method is based on evaluating the residual function only at the"selected points".{i).we can take the grid point values of the approximate solution, the set fuN(zi)),as the unknowns instead of the series coefficients.Given the value of a function at (N+1)points,we can compute the (N+1)series coefficients {an}through polynomial or trigonometric interpolation.Indeed,this symbolic equation series coefficients{an)grid point valuesfuN(i)} (1.14) is one of the most important themes we will develop in this course,though the mechanics of interpolation will be deferred to Chapters 4 and 5. Similarly,the finite element and spectral element algorithms approximate derivatives as a weighted sum of grid point values.However,only those points which lie within a given subdomain contribute directly to the derivative approximations in that subdomain. (Because the solution in one subdomain is matched to that in the other subdomain,there is an indirect connection between derivatives at a point and the whole solution,as true of finite differences,too.)Figure 1.3 compares the regions of direct dependency in derivative formulas for the three families of algorithms. Figs.1.4 and 1.5 compare the weights of each point in the second and fourth-order fi- nite difference approximations with the N =10 Fourier pseudospectral weights.Since the basis functions can be differentiated analytically and since each spectral coefficient an is determined by all the grid point values of u(),it follows that the pseudospectral differen- tiation rules are not 3-point formulas,like second-order finite differences,or even 5-point formulas,like the fourth-order expressions;rather,the pseudospectral rules are N-point formulas.To equal the accuracy of the pseudospectral procedure for N =10,one would need a tenth-order finite difference or finite element method with an error of O(h10). As N is increased,the pseudospectral method benefits in two ways.First,the intervalh
1.4. COMPARISONS WITH FINITE DIFFERENCES 7 Figure 1.4: Weights wj in the approximation df /dx |x=x0≈ P j wjf(x0 + jh) where x0 = π and h = π/5. In each group, the Fourier weights are the open, leftmost bars. Middle, crosshatched bars (j = ±1, ±2 only): Fourth-order differences. Rightmost, solid bars (j = ±1 only): weights for second order differences. The function O( ), the “Landau gauge symbol”, denotes that in order-of-magnitude, the errors are proportional to h2 and h4, respectively. Since the pseudospectral method is based on evaluating the residual function only at the “selected points”, {xi}, we can take the grid point values of the approximate solution, the set {uN (xi)}, as the unknowns instead of the series coefficients. Given the value of a function at (N+1) points, we can compute the (N + 1) series coefficients {an} through polynomial or trigonometric interpolation. Indeed, this symbolic equation series coefficients{an} ⇐⇒ grid point values{uN (xi)} (1.14) is one of the most important themes we will develop in this course, though the mechanics of interpolation will be deferred to Chapters 4 and 5. Similarly, the finite element and spectral element algorithms approximate derivatives as a weighted sum of grid point values. However, only those points which lie within a given subdomain contribute directly to the derivative approximations in that subdomain. (Because the solution in one subdomain is matched to that in the other subdomain, there is an indirect connection between derivatives at a point and the whole solution, as true of finite differences, too.) Figure 1.3 compares the regions of direct dependency in derivative formulas for the three families of algorithms. Figs.1.4 and 1.5 compare the weights of each point in the second and fourth-order fi- nite difference approximations with the N = 10 Fourier pseudospectral weights. Since the basis functions can be differentiated analytically and since each spectral coefficient an is determined by all the grid point values of u(x), it follows that the pseudospectral differentiation rules are not 3-point formulas, like second-order finite differences, or even 5-point formulas, like the fourth-order expressions; rather, the pseudospectral rules are N-point formulas. To equal the accuracy of the pseudospectral procedure for N = 10, one would need a tenth-order finite difference or finite element method with an error of O(h10). As N is increased, the pseudospectral method benefits in two ways. First, the interval h
8 CHAPTER 1.INTRODUCTION 2 多 ● 2 6 8 不 2n Figure 1.5:Same as previous figure except for the second derivative.Hollow bars:pseu- dospectral.Cross-hatched bars:Fourth-order differences.Solid bars:Second-order differ- ences. between grid points becomes smaller-this would cause the error to rapidly decrease even if the order of the method were fixed.Unlike finite difference and finite element methods, however,the order is not fixed.When N increases from 10 to 20,the error becomes O(h20) in terms of the new,smaller h.Since h is O(1/N),we have Pseudospectral errorO[(1/N)N] (1.15) The error is decreasing faster than any finite power of N because the power in the error formula is always increasing,too.This is "infinite order"or "exponential"convergence.! This is the magic of pseudospectral methods.When many decimal places of accu- racy are needed,the contest between pseudospectral algorithms and finite difference and finite element methods is not an even battle but a rout:pseudospectral methods win hands-down.This is part of the reason that physicists and quantum chemists,who must judge their calculations against experiments accurate to as many as fourteen decimal places (atomic hydrogen maser),have always preferred spectral methods. However,even when only a crude accuracy of perhaps 5%is needed,the high order of pseudospectral methods makes it possible to obtain this modest error with about half as many degrees of freedom-in each dimension-as needed by a fourth order method.In other words,spectral methods,because of their high accuracy,are memory-minimizing.Problems that require high resolution can often be done satisfactorily by spectral methods when a three-dimensional second order finite difference code would fail because the need for eight or ten times as many grid points would exceed the core memory of the available computer. "Tis true that virtual memory gives almost limitless memory capacity in theory.In prac- tice,however,swapping multi-megabyte blocks of data to and from the hard disk is very slow.Thus,in a practical (as opposed to theoretical)sense,virtual storage is not an option when core memory is exhausted.The Nobel Laureate Ken Wilson has observed that be- cause of this,memory is a more severe constraint on computational problem-solving than Chapter 2 shows show that the convergence is always exponential for well-behaved functions.but(1.15)is usually too optimistic.The error in an N-point method is O(M[nh")where M(n)is a proportionality constant: we ignored the (slow)growth of this constant with n to derive (1.15)
8 CHAPTER 1. INTRODUCTION Figure 1.5: Same as previous figure except for the second derivative. Hollow bars: pseudospectral. Cross-hatched bars: Fourth-order differences. Solid bars: Second-order differences. between grid points becomes smaller – this would cause the error to rapidly decrease even if the order of the method were fixed. Unlike finite difference and finite element methods, however, the order is not fixed. When N increases from 10 to 20, the error becomes O(h20) in terms of the new, smaller h. Since h is O(1/N), we have Pseudospectral error ≈ O[(1/N) N ] (1.15) The error is decreasing faster than any finite power of N because the power in the error formula is always increasing, too. This is “infinite order” or “exponential” convergence.1 This is the magic of pseudospectral methods. When many decimal places of accuracy are needed, the contest between pseudospectral algorithms and finite difference and finite element methods is not an even battle but a rout: pseudospectral methods win hands-down. This is part of the reason that physicists and quantum chemists, who must judge their calculations against experiments accurate to as many as fourteen decimal places (atomic hydrogen maser), have always preferred spectral methods. However, even when only a crude accuracy of perhaps 5% is needed, the high order of pseudospectral methods makes it possible to obtain this modest error with about half as many degrees of freedom – in each dimension – as needed by a fourth order method. In other words, spectral methods, because of their high accuracy, are memory-minimizing. Problems that require high resolution can often be done satisfactorily by spectral methods when a three-dimensional second order finite difference code would fail because the need for eight or ten times as many grid points would exceed the core memory of the available computer. ’Tis true that virtual memory gives almost limitless memory capacity in theory. In practice, however, swapping multi-megabyte blocks of data to and from the hard disk is very slow. Thus, in a practical (as opposed to theoretical) sense, virtual storage is not an option when core memory is exhausted. The Nobel Laureate Ken Wilson has observed that because of this, memory is a more severe constraint on computational problem-solving than 1Chapter 2 shows show that the convergence is always exponential for well-behaved functions, but (1.15) is usually too optimistic. The error in an N-point method is O(M[n]hn) where M(n) is a proportionality constant; we ignored the (slow) growth of this constant with n to derive (1.15)
1.5.PARALLEL COMPUTERS 9 CPU time.It is easy to beg a little more time on a supercomputer,or to continue a job on your own workstation for another night,but if one runs out of memory,one is simply stuck-unless one switches to an algorithm that uses a lot less memory,such as a spectral method. For this reason,pseudospectral methods have triumphed in metereology,which is most emphatically an area where high precision is impossible! The drawbacks of spectral methods are three-fold.First,they are usually more difficult to program than finite difference algorithms.Second,they are more costly per degree of freedom than finite difference procedures.Third,irregular domains inflict heavier losses of accuracy and efficiency on spectral algorithms than on lower-order alternatives.Over the past fifteen years,however,numerical modellers have learned the "right"way to im- plement pseudospectral methods so as to minimize these drawbacks. 1.5 Parallel Computers The current generation of massively parallel machines is communications-limited.That is to say,each processor is a workstation-class chip capable of tens of megaflops or faster,but the rate of interprocessor transfers is considerably slower. Spectral elements function very well on massively parallel machines.One can assign a single large element with a high order polynomial approximation within it to a single processor.A three-dimensional element of degree p has roughly p3 internal degrees of freedom,but the number of grid points on its six walls is O(6p2).It is these wall values that must be shared with other elements-i.e.,other processors-so that the numerical solution is everywhere continuous.As p increases,the ratio of internal grid points to boundary grid points increases,implying that more and more of the computations are internal to the element,and the shared boundary values become smaller and smaller compared to the total number of unknowns.Spectral elements generally require more computation per unknown than low order methods,but this is irrelevant when the slowness of interprocessor data transfers,rather than CPU time,is the limiting factor. To do the same calculation with low order methods,one would need roughly eight times as many degrees of freedom in three dimensions.That would increase the interpro- cessor communication load by at least a factor of four.The processors would likely have a lot of idle time:After applying low order finite difference formulas quickly throughout its assigned block of unknowns,each processor is then idled while boundary values from neighboring elements are communicated to it. Successful applications of spectral elements to complicated fluid flows on massively parallel machines have been given by Fischer(1990,1994a,b,1997)Iskandarani,Haidvogel and Boyd (1994),Taylor,Tribbia and Iskandarani(1997)and Curchitser,Iskandarani and Haidvogel(1998),among others. 1.6 Choice of basis functions Now that we have compared spectral methods with other algorithms,we can return to some fundamental issues in understanding spectral methods themselves. An important question is:What sets of "basis functions"n(z)will work?It is obvi- ous that we would like our basis sets to have a number of properties:(i)easy to compute (ii)rapid convergence and (iii)completeness,which means that any solution can be repre- sented to arbitrarily high accuracy by taking the truncation N to be sufficiently large
1.5. PARALLEL COMPUTERS 9 CPU time. It is easy to beg a little more time on a supercomputer, or to continue a job on your own workstation for another night, but if one runs out of memory, one is simply stuck – unless one switches to an algorithm that uses a lot less memory, such as a spectral method. For this reason, pseudospectral methods have triumphed in metereology, which is most emphatically an area where high precision is impossible! The drawbacks of spectral methods are three-fold. First, they are usually more difficult to program than finite difference algorithms. Second, they are more costly per degree of freedom than finite difference procedures. Third, irregular domains inflict heavier losses of accuracy and efficiency on spectral algorithms than on lower-order alternatives. Over the past fifteen years, however, numerical modellers have learned the “right” way to implement pseudospectral methods so as to minimize these drawbacks. 1.5 Parallel Computers The current generation of massively parallel machines is communications-limited. That is to say, each processor is a workstation-class chip capable of tens of megaflops or faster, but the rate of interprocessor transfers is considerably slower. Spectral elements function very well on massively parallel machines. One can assign a single large element with a high order polynomial approximation within it to a single processor. A three-dimensional element of degree p has roughly p3 internal degrees of freedom, but the number of grid points on its six walls is O(6p2). It is these wall values that must be shared with other elements – i. e., other processors – so that the numerical solution is everywhere continuous. As p increases, the ratio of internal grid points to boundary grid points increases, implying that more and more of the computations are internal to the element, and the shared boundary values become smaller and smaller compared to the total number of unknowns. Spectral elements generally require more computation per unknown than low order methods, but this is irrelevant when the slowness of interprocessor data transfers, rather than CPU time, is the limiting factor. To do the same calculation with low order methods, one would need roughly eight times as many degrees of freedom in three dimensions. That would increase the interprocessor communication load by at least a factor of four. The processors would likely have a lot of idle time: After applying low order finite difference formulas quickly throughout its assigned block of unknowns, each processor is then idled while boundary values from neighboring elements are communicated to it. Successful applications of spectral elements to complicated fluid flows on massively parallel machines have been given by Fischer(1990, 1994a,b, 1997) Iskandarani, Haidvogel and Boyd (1994), Taylor, Tribbia and Iskandarani(1997) and Curchitser, Iskandarani and Haidvogel(1998), among others. 1.6 Choice of basis functions Now that we have compared spectral methods with other algorithms, we can return to some fundamental issues in understanding spectral methods themselves. An important question is: What sets of ”basis functions” φn(x) will work? It is obvious that we would like our basis sets to have a number of properties: (i) easy to compute (ii) rapid convergence and (iii) completeness, which means that any solution can be represented to arbitrarily high accuracy by taking the truncation N to be sufficiently large