Analysis Of Finite Elements And Finite Differences For Shallow Water .

B. Neta / Finitedifferencesfor shallowwaterequations150Then the Fouriertransformsof the input and outputare relatedby9{J(k) J(k)Fi(k)’(4.2)This is the representationof the system in the transformis given by the convolution theoremdomain.The representationin thephysical domain(4.3)d( x is known as impulse response of the system to an input 6(x), since (x) lrn4(x- Y) (Y) dr.(4.4)-ccThe magnitudeand phase of J(k)can be found by taking a sinusoidalinputyi( x) ejkd.r,where j \r-1.y,,(x)(4.5)Then (k,)ejk? 1 (ki)1 ejkcb ejk,x 1 (k,)1 ejk#(x kcb/kt).(4.6)Thus, the magnitude141 of the transfer function is the factor by which the amplitudeof ais amplified or attenuated.The argumentk,/k,is the shift of the phase of theinput.In general, the differentfrequenciesare amplified/attenuatedand shifted in phase bydifferent amounts. The distortion (change in shape) of the output depends on both effects.In electrical engineering,any linear invariant system which can be described by a transferfunction is called a filter. When a continuous process is discretized, one obtains a discrete orsinusoidaldigital filter.For example,2 f(x),wheref(x)WIis input and y(x)y(x Ax) y(x)is output,is a filter. A digital filter for this could be Axf(x)(nothing more than Euler’s method). The transferby Fourier transform. To this end, one recalls(4.8)functiony (x Ax) /m y(x Ax) e-jkx dx ejkAxy (x).--CDThe transferfunction6(k) ; of digital filters can be determined(4.9)for the filter is(4.10)

151B. Neta / Finite differences for shallow water equationsand for the above digital filter is (easy exercise)e-jk(6”(k) Ax/2(4.11)sin ik Ax ’j;AxThe main difference between continuous and digital filters is the so-called aliasing. A samplingdevice at interval Ax is not capable of resolving waves of frequency greater than r/Axand ifenergy is present in a sampled continuous signal at such high frequencies,it will erroneouslyberesolved into a frequency lower than n/Ax.5. Transfer functionanalysisfor one-dimensionalshallow water equations5.1. Analytic caseThe one-dimensional[48]:linearizedshallow water equations; fu 0,with no mean flow were analyzeddh; H; O.in(5.1)This model is especially importantin the study of geostrophicadjustment.This process hasbeen studied in some detail from several approachesin [l-3,45,47,69].The process ofgeostrophicadjustmentis importantbecauseit is the primary mechanismby which theatmospherereacts to errors in the initial data. The process takes place by means of wavepropagationfrom local regions of initial imbalance, leaving behind a steady semigeostrophicbalance flow. Cahn [3] showed that a sudden perturbationcaused by an addition of momentumwill cause the disturbanceamplitude to first grow linearly with time followed by an asymptoticstatesh -Lcos(ft \&The speed of propagationgh-0f(.if7-r).of energy is ,/&.(54Obukhov[40] showed that for two-dimensionalflow(5.3)Schoenstadt[48] has shown that the amplitude distortion in this system is governed by one ofthe three factors l/v, k/a or k/v2, or the square of one of these were v/k, the phase speed,is given bylJ kfk/1 R;k2,(54

B. Neta / Finite differences for shallow water equations152u,L’,hu, L‘, hi-lu,i 1,hhi lu, L’hi-liA grid13,hi lB grid 1,h11hu, 1i-luiL‘, hu, hi li-lC gridL’u,hiI’u,hi lD gridFig. 6.wheredigHfR,, -is the Rossby radius of deformation.When drawing the amplitude coefficients,it is clear that the terms with the coefficientl/vhave their low frequenciesleast affected (low-pass filter), the terms with coefficientk/v arehigh-pass filters and those with k/v* are band-pass filters.5.2. Semi-discrete caseSchoenstadt[48] analyzed second- and fourth-orderfinite-differenceand finite-elementschemes for which the variables are either unstaggered(A grid) or staggered (B, C, D grids)(see Fig. 6). Finite-elementD grid is not given there. See also [25].Neta and Navon [34] analyzed the Turkel-Zwasexplicit large time step scheme [61]. Theby r/v, S/v, Q-S/V*. For thefilter coefficientsl/v,k/v, k/v2 are replaced, respectively,Turkel-Zwasscheme the filter coefficients change but not always in the same manner. In thisThe filter coefficientsfor eachcase there are extra coefficientsnamely r/v*and rS/v2.scheme can be found in [34,48]. The phase speed for the Turkel-Zwasscheme is(WwhereT (2 cos sk AX)(5.7)andsin sk Axs SAX’(54Schoenstadtconcludes that Scheme A produces the greatest distortion especially at highfrequencies.Other methods very accurately approximatethe transfer function for high-passfilter. Scheme B overstates the amount of energy distributed into the short wave, while Scheme

B. Neta / Finite differences for shallow water equationsu,v,hu,v,h153hh11UP’u,v,hI,hhLnFig. 7(b). Grid B.Fig. 7(a). Grid A.finiteC understatesthis. The finite-elementmethods appear more accurate than fourth-orderdifference.The group velocity for Schemes A, D is reversed for short waves. Based on theseone-dimensionalresults. Scheme C is best.6. Transfer functionanalysis for two-dimensionalshallow water equationsThe transfer function analysis initiated by Schoenstadt[48] for the one-dimensionalcase wasextendedto two dimensionsin [32,33,64]. The transfer coefficientsare functionsof bothwavenumbersk, 1. The methods analyzed are the same as in the one-dimensionalcase. What issurprisingis that the Turkel-Zwasscheme must be modifiedto get convergence.Thismodificationwas not enough to make the method competitivewith finite-elementor fourthorder C scheme. See Fig. 7 for two-dimensionalarrangementof variables. The filter coeffi-h”hh”””Fig. 7(c). Grid C.nh”Fig. 7(d). Grid D.h

B. Neta / Finite differences for shallow water equations154Table 5SchemeCY,ff A211B211c21cos ;kd cos ld,B(k, 1)sin kddsin ;kdIcospiD21cos ;kd cos ldA411B411c41cos ;kd cos ldD41cm ;kd cos ldsin kdidsin kdpcos;ldd8 sin kd -sin 2kd6d-sin :kd 27 sin tkd12d- sin kd 27 sin kd12d8 sin kd - sin 2kdTZ134 cos ksd cos lsd)FET (3 cos kd 2 cos ;kd cos Id);(3 cos kd 2 cos ;kd cos Id)FER (2 cos;(2 coskd)(2 cosId);ldkd)(2 cosId)sin ksd6da; RlQ2(k,1) p’(A k))7cos ;ldsd2( sin kd sin ikd cos ld) (2 coscients for all schemesare a,/(a,v ,P(k, I /(a,v),P(l, k)/(a,v),c ,p(l, )/((Y:v and P(k, OpU; k)/(azv21.The frequency1, is given bycfv fcos i/d3dsin kdIdIda,P(k,l)/(azv2),(6.1)where ayx, (Y,,, p for each scheme are listed in Table 5. The only exception to the above is in thefinite-elementscheme on isosceles triangles (FET), in which P(E, k) should be replaced by(cos ikd sin Zd)/d.7. Rossby wave frequenciesandThe hydrostaticprimitive equationnumericalmodels that are used for atmosphereoceanographicprediction permit inertial gravity waves, Rossby waves and advective effects. Theinfluence of a numerical method on each of these types of motion is most easily analyzed byseparating the linearizedpredictionequations into vertical modes with an equivalentdepthanalysis (see, e.g., [13]). In this case the equations for each vertical mode are just the linearizedshallow water equations (7.1). Arakawa and Lamb [l] analyzed inertial gravity wave motions for

B. Neta/ Finite differencesfar shallow water equations155grids A, B, C, D. They found that the geostrophicadjustment for grids A, D is poor and thatthe adjustment for grids B, C is good.In this section we discuss the treatmentof Rossby waves in vorticity-divergenceshallowwater formulationswith various finite-elementand finite-differenceschemes. For comparisonthe finite-differenceprimitive equations solution for grids A, B and C are also included (seeB',641).7.1.Primitiue formThe linearizedshallow water equationson a beta plane can be writtenas follows:&A0,(7.1)al!ah; “b g- 0,ay(7.2)z-fi; ggaudY-1The frequencycan be obtainedWF - o.(7.3)byPO*6 E (YK2’where (Y, I, 6, E depend on the method used and where PO is an averageThese parametersare listed in [38] and reproducedin Table 6.(7.4)value of df/dy.7.2. Vorticity-diuergence formThe vorticity-divergenceequation set, which is obtained by differentiatingrespect to x and y, respectively,and combining, can be writtenz fD pc 0,(7.1) and (7.2) with(7.5)dl-f Su y ) o.(7.6);(7.7)aD HD O.

B. Neta156Table 6The operatorsfor the various numerical/ Finite differencesschemesVorticity-divergenceSchemefor shallow water equationsfor the shallow water equationsformFinite larStaniforthand Mitchellt55,561a1i(3 cos2(sin9P6P2Ek*x*X sin 2costxcos XcosY)3Axsin*fXY) (2 cossinZLX(t A;)2(; Ax)*(3 cosx-4cos xcosY)2 Ay*xx2 cosY)sin Xz (2 cos Y)sin’ Y---- (2 cos(:AY) (2 cosX)(2 cossin X-&2 cosY)Y)sin’ X--- 5 cosY)Y)(t Ax)sin’ Y cosX)(DAY) 5 cos X)(5 cos Y)X)(Y’*X k;(2 cosAx,Y lAy.To isolate the Rossby mode more easily we apply the quasi-geostrophicset (7SH7.7)approximationto the(7.9)ahz HD O,where f,, and PO are evaluated at an approximatecentral latitude.before with an appropriatevalue for the parametersIY, 9, 6, E.(7.10)The frequencyis given as7.3. Semi-implicit fourth orderStaniforth and Mitchell [56] introduced a finite-elementscheme in which the first derivativesare approximatedto fourth order. The frequency in this case is given by(7.11)The parametersLX,CX’,I/I, 6, E are listed in Table 6.Neta and Williams [38] have shown that the finite-elementmethod based on isoscelestriangles (FET) and on rectangles (FER) are as good as the method due to StaniforthandMitchell GM). These three methods are better than all others.

B. Neta / Finite differences for shallow water equations157Table d orderFourth1sin X14 sin Xorder1 sin 2X3 Ax6 Axcos2X-16cosX 15AXsin’ XABC1sin X-cosAXsin2X1sin Xcos2;x cos2;Ysin Xcos2YAxsin’;X (4 Ax) 6 Ax’(4 Ax) AY)’Axsin’fX ;(l cos(f Ax)sin2LY (l cossin2Y6 Ay28. Semi-Lagrangian YAX2cos2Y-16cosY 15sin21Y*Y)sin” Y (5 AY)*x (t AY)Ay2methodsSemi-Lagrangianmethods were first proposedelsewhere. See recent review [53].Consider the advection-diffusionequationa# t u(x,atwhere II is the diffusion[42-4414(x,formFinite differencesFinite differences(tequationt) b(x,t) vcoefficient.t At) - 4(x -a,t)in [21,46], and improvedv2 , Iapproximationmay be written1(x,t Ar) v v24 1(r-a )],2in [59] andt),A semi-LagrangianAtrecentlyas(84wherea(x) At V(X - ;a,(84t ; At).Suppose x is a grid point of a given regular mesh. Suppose that at these mesh points weknow 4 at time t and V at time t i At. A semi-Lagrangianalgorithm is then the following.(1) Solve (8.3) iterativelyaci ‘)(x)for a(x)by At V(x - zci), t 3 At),and by an interpolationof I/ between(2) Evaluatev V* J I cX,tj at meshmethod. Use interpolationto evaluate(4 - ’ At V’4)1 (x,t Af) i 0, 1,. .,mesh points.points by finite-difference,the right-hand side of(4 ’Atv’4)1 (x-ay,t).(8.4)finite-elementor spectral(8.5)

158B. Neta / Finite differences for shallow water equations(3) Solve the above Helmholtz equation at mesh points (SOR, ADI, finite Fourierfinite elements).(4) Repeat the steps to obtain 4(x, t 2 At).transform,Remarks. (1) Eulerian schemes suffer from serious numerical dispersion and diffusion.(2) High-ord er E u 1erian schemes tend to produce artificial extrema because of the unrealisticphase speeds of the high wave number.(3) Lagrangianschemes suffer from distortionof the initial grid after a long integration(more than twelve hours).In two dimensions it is appropriateto consider the pseudo staggering suggested in [7]. “Thescheme works well with the primitive form of the equations, uses unstaggered grid but does notpropagate small scale energy in the wrong direction, works well with variable resolution and ascomputationallyefficient as staggered formulationusing the primitive form of the equations.”9. Domain decompositionOne of the major developing areas in numerical analysis is parallel computationwhich offersthe possibility of significantly faster computationalspeeds. Domain decompositionmethods arebased on subdivision of the domain into several (maybe overlapping)subdomains and solvingthe problem on several subdomainsin parallel. The methods can be regarded as divide-andconquer algorithms (see [65]). The interactionsbetween the solutions on the subdomains leadto an iterative technique.When the number of subdomainsis large, one can improve theconvergenceof this iterative procedureby using a coarse grid to obtain starting values of thesolution on the interfaces. In this respect, the methods are similar to multigrid schemes. Thecrucial point is how to pass informationfrom one domain to other processors. Two differentapproacheswere followed in the literature.The first approach is based on decompositionsofthe domain into contiguous regions (see, e.g., [12] and referencesthere). The second is basedon having overlappingregions (Schwartz alternatingmethod, see, e.g., [22-241 and referencesthere).The main difficulty of such parallel techniques is in the initial assignment of values to theinterfaces between subdomains. The more accurate such values are, the faster the convergence.We have already mentioned the idea based on multigrid. There are several other possibilities toaccelerateconvergence(see the proceedingsof four internationalconferenceson domaindecomposition[5,6,14,X]). According to Lions [24] the approach based on optima1 control asfollowed by Glowinski et al. [16,17] is a bit faster, at least for Laplace equations.Neta and Okamoto 1351 have suggested an accelerationbased on boundary elements. Thesetwo approachesfor shallow water equations are now under investigationby Neta and Navon.10. ConclusionsThis paper reviewed analyses of finite-elementand finite-differenceschemes for the solutionof the shallow water equations. The superiority of certain finite-elementmethods is indicated.Some recent results in parallel computationsare included.

B. Neta / Finite differences for shallow water equations159References[I] A. Arakawa and V.R. Lamb, Computationaldesign of the basic dynamical processes of the UCLA GeneralCirculation Model, in: Methods Comput. Phys. 17 (Academic Press, New York, 1977) 173-265.[2] W. Blumen, Geostrophicadjustment,Reu. Geophys. Space Phys. 10 (1972) 485-528.[3] A. Cahn, An investigationof the free oscillations of a simple current system, J. Meteorology 2 (1945) 113-119.[4] M.M. Cecchi, Domain decompositionmethod for shallow water flow problems, in: H. Niki and M. Kawahara,eds., Proc. Internat. Conf on Computational Methods in Flow Analysis (Okayama Univ. of Science, Japan, 1988)1322-1329.[5] T.F. Chan, R. Glowinski, J. Periaux and O.B. Widlund, eds., Proc. 2nd Internat. Syrnp. on Domain Decomposition Methods (SIAM, Philadelphia,PA, 1989).[6] T.F. Chan, R. Glowinski, J. Periaux and O.B. Widlund, eds., Proc. 3rd Internat. Symp. on Domain Decomposition Methods for PDEs (SIAM, Philadelphia,PA, 1989).[7] J. Cot& S. Gravel and A. Staniforth,Improving tegrationschemes by pseudo-staggering,Monthly Weather Reo. 118 (1990) 2718-2731.[8] J. C8tC and A. Staniforth, A two-time-levelsemi-Lagrangiansemi-implicit scheme for spectral models, MonthlyWeather Ret,. 116 (1988) 2003-2012.[9] M.P.J. Cullen and C.D. Hall, Forecastingand general circulation results from finite element models, Quart. .I.Roy. Meteorological Sot. 105 (1979) 571-592.[lo] C.A.J. Fletcher, Computational Galerkin Methods (Springer, New York, 1984).[ 1 l] M.G.G. Foreman, A two dimensional dispersion analysis of selected methods for solving the linearized shallowwater equations, J. Comput. Phys. 56 (1984) 287-323.[12] D. Funaro, Multidomainspectral approximationof elliptic equations,Numer. Methods Partial DifferentialEquations 2 (1986) 187-205.[13] A.E. Gill, Atmosphere-OceanDynamics (Academic Press, New York, 1982).[14] R. Glowinski, G.H. Golub, G.A. Meurant and J. Periaux, eds., Proc. 1st Internat. Symp. on Domain Decomposition Methods for Partial Differential Equations (SIAM, Philadelphia,PA, 1988).[15] R. Glowinski, Y.A. Kuznetsov,G. Meurant, J. Periaux and 0. Widlund, eds., Proc. 4th Internat. Symp. onDomain Decomposition Methods for Partial Differential Equations (SIAM, Philadelphia,PA, to appear).[16] R. Glowinski and P. Le Tallec, AugmentedLagrangian interpretationof the nonoverlappingSchwarz alternating scheme, in: T.F. Chan et al., eds., Proc. 3rd Internat. Symp. on Domain Decomposition Methods for PDEs(SIAM, Philadelphia,PA, 1989) 224-231.[17] R. Glowinski, Q.V. Dinh and J. Periaux, Domain decompositionmethods for nonlinearproblemsin fluiddynamics, Comput. Methods Appl. Mech. Engrg. 40 (1) (1983) 27-109.[18] G.J. Haltiner and R.T. Williams, Numerical Prediction and Dynamic Meteorology (Wiley, New York, 1980).[19] R.G. Kelley and R.T. Williams, A finite elementpredictionmodel with variable elementsizes, NavalPostgraduateSchool Report NPS-63 Wu76101, Monterey, CA, 1976.[20] I. Kinnmark,The Shallow- Water Wat,e Equations: Formulation, Analysis, and Application (Springer,Berlin,1985).,[21] T.N. Krishnamurti,Numerical integrationof primitive equations by a quasi-Lagrangianadvective scheme, J.Appl. Met. l(1962) 508-521.[22] P.L. Lions, On the Schwarz alternatingmethod, in: R. Glowinski et al., eds., Proc. 1st Internat. Symp. onDomain Decomposition Methods for Partial Diff erential Equations (SIAM, Philadelphia,PA, 1988) l-42.[23] P.L. Lions, On the Schwarz alternating method II: Stochastic interpretationand order properties,in: T.F. Chanet al., eds., Proc. 2nd Internat. Symp. on Domain Decomposition Methods (SIAM, Philadelphia,PA 1989) 47-70.[24] P.L. Lions, On the Schwarz alternating method III: A variant for nonoverlappingsubdomains,in: T.F. Chan etal., eds., Proc. 3rd Internat. Symp. on Domain Decomposition Methods for PDEs (SIAM, Philadelphia,PA, 1989)202-223.[25] F. Mcsinger, Dependenceof vorticity analogue and the Rossby wave phase speed on the choice of horizontalgrid, Sci. Math. 10 (1979) 5-1.5.[26] A.R. Mitchell and R. Wait, The Finite Element Method in Partial Diff erential Equations (Wiley, London, 1977).[27] K.B. Monk, Studies of baroclinicflow over topographyusing semi-Lagrangian,semi-implicitmodel, Dept.Meteorology,Naval PostgraduateSchool, Monterey, CA, 1989.

160B. Neta / Finite differences for shallow water equations[28] I.M. Navon, A Numerov-Galerkintechnique applied to a finite-elementshallow-waterequations model withexact conservationof integral constraints,in: T. Kawai, ed., Proc. 4th Internat. Symp. on Finite Element Methodsin Flow Problems (Chuo Univ., Tokyo, 1982) 75-85.[29] I.M. Navon, A Numerov-Galerkintechnique applied to a finite element shallow-waterequations model withenforced conservationof integral invariants and selective lumping, J. Comput. Phys. 52 (1983) 313-339.[30] I.M. Navon, A review of finite-elementmethods for solving the search Institute Report FSU-SCRI-88-49,Tallahassee,FL, 1988.[31] I.M. Navon and H.A. Riphagen,An implicit compact fourth-orderalgorithm for solving the shallow-waterequations in conservation-lawform, Monthly Weather Rec. 107 (1979) 1107-1127.[32] B. Neta, Analysis of the Turkel-Zwasscheme for the shallow-waterequations,in: W.F. Ames, ed., I&ICSTrans. on Scientific Computing ‘88, Vols. 1.1 and

In this review article we discuss analyses of finite-element and finite-difference approximations of the shallow water equations. An extensive bibliography is given. 0. Introduction In this article we review analyses of finite-element and finite-difference methods for the approximation of the shallow water equations.