A pseudospectral quadrature method for Navier-Stokes equations on rotating spheres
AA PSEUDOSPECTRAL QUADRATURE METHOD FORNAVIER-STOKES EQUATIONS ON ROTATING SPHERES
M. GANESH, Q. T. LE GIA, AND I. H. SLOAN
Abstract.
In this work, we describe, analyze, and implement a pseudospec-tral quadrature method for a global computer modeling of the incompressiblesurface Navier-Stokes equations on the rotating unit sphere. Our spectrallyaccurate numerical error analysis is based on the Gevrey regularity of the so-lutions of the Navier-Stokes equations on the sphere. The scheme is designedfor convenient application of fast evaluation techniques such as the fast Fouriertransform (FFT), and the implementation is based on a stable adaptive timediscretization. Introduction
In this paper we develop a pseudospectral quadrature method for the surfaceNavier-Stokes partial differential equations (PDEs) on the rotating unit sphere.Whereas the finite element method is best suited for handling non-smooth pro-cesses, the spectral global basis computer models are very efficient and performextremely well for processes with smooth regularity. For example, exponential con-vergence properties of the global Fourier basis spectral Galerkin methods (withoutquadrature) for the Ginzburg-Landau and Navier-Stokes PDEs on two dimensionalperiodic cells are based on the Gevrey regularity of solutions of the PDEs [7, 20].The complex three dimensional flows in the atmosphere and oceans are con-sidered to be accurately modeled by the Navier-Stokes PDEs of fluid mechanicstogether with classical thermodynamics [22]. Difficulties in computer modeling inthese PDEs resulted in several simplified models for which spectral approximationsare well known [2, 15, 22]. A famous open problem is to prove the global regularityfor the three dimensional incompressible Navier-Stokes PDEs [26]. However, theprecise Gevrey regularity of the unique solution of the (practically relevant) sur-face Navier-Stokes PDEs on the rotating sphere was proved in [5]. (Because theEarth’s surface is an approximate sphere, a standard surface model, to study globalatmospheric circulation on large planets, is the sphere.)Consequently, a natural next step is to describe, analyze, and implement anexponentially converging pseudospectral method for the Navier-Stokes PDEs on therotating sphere. In addition to the continuous model regularity results in [5], thispaper is also motivated by the recent work [12], where discrete computer modeling ofthe Navier-Stokes PDEs on one-dimensional and toroidal domains [8] was extendedto the unit sphere.
Date : October 29, 2018.2000
Mathematics Subject Classification.
Primary 65M12; Secondary 76D05.
Key words and phrases.
Navier-Stokes equations, unit sphere, vector spherical harmonics. a r X i v : . [ m a t h . NA ] S e p M. GANESH, Q. T. LE GIA, AND I. H. SLOAN
This paper is concerned with both implementation of our algorithm and itsnumerical analysis. The main numerical analysis contributions compared to resultsin [5], for the continuous problem, and in [12], for a discrete problem, are as follows.For the spatially discrete pseudospectral quadrature Galerkin solutions of theNavier-Stokes equations, we prove (i) the stability (that is, uniform boundednessof approximate solutions, independent of their truncation parameter N ), see The-orem 4.1; and (ii) a spectrally accurate rate of convergence [that is, O ( N − s ) ac-curacy, with s depending on the smoothness of input data], see Theorem 4.2. Weachieve these results by first generalizing the main regularity result in [5] to com-plex valued times (see Theorem 2.2 and 2.3), and then using this to prove thespectral rate of convergence of the time derivative of a Stokes projection compar-ison function (see Theorem 3.2). This time-derivative error result plays a crucialrole in the proof of spectral convergence of the approximate solutions (see the proofof Theorem 4.2). We note that the main result in [12, page 978] establishes onlyconvergence of the semi-discrete Galerkin method without quadrature in L p norms,but does not establish either stability or rate of convergence of the scheme.The rate of convergence results, supported by numerical experiments, formed thecore part of research on the Navier-Stokes equations on two dimensional domainsover the last few decades, see [8, 9, 27] and references therein. There is a vast lit-erature on numerical methods and analysis for the Navier-Stokes PDE on boundedEuclidean domains (see [8, 9, 27] and references therein), but their counterparts onclosed manifolds are rarer (see [12] and references therein). The implementationof the scheme in [12, page 978] is based on a fixed time-step explicit Runge-Kuttamethod that has a small stability region for the systems of ordinary differentialequations arising from the spatial discretization.The outline of this paper is as follows. In the next section, we recall variousknown preliminary results associated with the Navier-Stokes PDEs on the unitsphere, in strong and weak form. In Section 3, we introduce essential computationaland numerical analysis tools required for the discretization and analysis of theNavier-Stokes equations. In Section 4, we describe and prove spectral accuracy ofa pseudospectral quadrature method and give implementation details required toapply the FFT and adaptive-in-time simulation of the Navier-Stokes equations. InSection 5, we demonstrate computationally the accuracy and applicability of thealgorithm for well known benchmark examples.2. Navier-Stokes equations on the rotating unit sphere
The surface Navier-Stokes equations (NSE) describing a tangential, incompress-ible atmospheric stream on the rotating two-dimensional unit sphere S ⊂ R canbe written as [5, 16, 17, 19, 28](2.1) ∂∂t u + ∇ u u − ν ∆ u + ω × u + 1 ρ Grad p = f , Div u = 0 , u | t =0 = u on S. Here u = u ( (cid:98) x , t ) = ( u ( (cid:98) x , t ) , u ( (cid:98) x , t ) , u ( (cid:98) x , t )) T is the unknown tangential divergence-free velocity field at (cid:98) x ∈ S and t ∈ [0 , T ], p = p ( (cid:98) x , t ) is the unknown pressure. Theknown components in (2.1) are the constant viscosity and density of the fluid, re-spectively denoted by ν, ρ , the normal vector field ω = ω ( (cid:98) x ) = ω ( (cid:98) x ) (cid:98) x for theCoriolis acceleration term, and the external flow driving vector field f = f ( (cid:98) x , t ). PSEUDOSPECTRAL METHOD FOR NSEs ON ROTATING SPHERES 3
The Coriolis function ω is given by ω ( (cid:98) x ) = 2Ω cos θ , where Ω is the angular veloc-ity of the rotating sphere, and θ is the angle between (cid:98) x and the north pole. Thevorticity of the flow associated with the NSE (2.1), in the curvilinear coordinatesystem, is a normal vector field, defined, for a fixed t ≥
0, by(2.2)
Vort u ( (cid:98) x , t ) = Curl (cid:98) x u ( (cid:98) x , t ) = (cid:98) x ∆Ψ( (cid:98) x , t ) , (cid:98) x ∈ S, for some scalar-valued vorticity stream function Ψ.All spatial derivative operators in (2.1)-(2.2) are surface differential operators,obtained by restricting the corresponding domain operators (defined in a neighbor-hood of S ) to the unit sphere, using standard differential geometry concepts onclosed manifolds in R [16, 17].Using the fact that the outward unit normal at (cid:98) x ∈ S is (cid:98) x , the Curl of a scalarfunction v , of a normal vector field w = w (cid:98) x , and of a tangential vector field v on S are respectively defined by(2.3) Curl v = − (cid:98) x × Grad v, Curl w = − (cid:98) x × Grad w, Curl (cid:98) x v = − (cid:98) x Div ( (cid:98) x × v ) . The surface diffusion operator acting on tangential vector fields on S is denoted by ∆ (known as the Laplace-Beltrami or Laplace-de Rham operator) and is definedas(2.4) ∆ v = Grad
Div v − Curl Curl (cid:98) x v . The following relations connecting the above operators will be used throughoutthe paper:(2.5) Div
Curl v = 0 , Curl (cid:98) x Curl v = − (cid:98) x ∆ v, ∆ Curl v = Curl ∆ v, (2.6)2 ∇ w v = − Curl ( w × v )+ Grad ( w · v ) − v Div w + w Div v − v × Curl (cid:98) x w − w × Curl (cid:98) x v . In particular, for tangential divergence-free vector fields, such as the solution u ofthe NSE, using (2.6), the nonlinear term in (2.1) can be written as(2.7) ∇ u u = Grad | u | − u × Curl (cid:98) x u . A weak formulation.
A standard technique for removing the scalar pressurefield from the Navier-Stokes equations is to multiply the first equation in (2.1) bytest functions v from a space with elements having properties of the unknownvelocity field u (in particular, Div v = 0) and then integrate to obtain a weakformulation. (The unknown p , can then be computed by solving a pressure Poissonequation, obtained by applying the surface divergence operator in (2.1).)To this end, we introduce the standard inner products on the space of all squareintegrable (i) scalar functions on S , denoted by L ( S ); and (ii) tangential vectorfields on S , denoted by L ( T S ):( v , v ) = ( v , v ) L ( S ) = (cid:90) S v v dS, v , v ∈ L ( S ) , (2.8) ( v , v ) = ( v , v ) L ( T S ) = (cid:90) S v · v dS, u , v ∈ L ( T S ) , (2.9)where dS = sin θdθdφ . Throughout the paper, the induced norm on L ( T S ) isdenoted by (cid:107) · (cid:107) and for other inner product spaces, say X with inner product( · , · ) X , the associated norm is denoted by (cid:107) · (cid:107) X . For example, for s >
0, standard
M. GANESH, Q. T. LE GIA, AND I. H. SLOAN norms in the scalar and vector valued functions Sobolev spaces H s ( S ) and H s ( T S )are denoted by (cid:107) · (cid:107) H s ( S ) and (cid:107) · (cid:107) H s ( T S ) , respectively. Since H ( T S ) = L ( T S ), (cid:107) · (cid:107) H ( T S ) = (cid:107) · (cid:107) .We have the following identities for appropriate scalar and vector fields [16,(2.4)-(2.6)]:( Grad ψ, v ) = − ( ψ, Div v ) , ( Curl ψ, v ) = ( ψ, Curl (cid:98) x v ) , (2.10) ( Curl Curl (cid:98) x w , z ) = ( Curl (cid:98) x w , Curl (cid:98) x z ) . (2.11)In (2.10), the L ( T S ) inner product is used on the left hand side and the L ( S )inner product is used on the right hand side. Throughout the paper, we identify anormal vector field w with a scalar field w and hence(2.12) ( ψ, w ) := ( ψ, w ) L ( S ) , w = (cid:98) x w, ψ, w ∈ L ( S ) . Using (2.5), smooth ( C ∞ ) tangential fields on S can be decomposed into twocomponents, one in the space of all divergence-free fields and the other through theHodge decomposition theorem [1]:(2.13) C ∞ ( T S ) = C ∞ ( T S ; Grad ) ⊕ C ∞ ( T S ; Curl ) , where(2.14) C ∞ ( T S ; Grad ) = { Grad ψ : ψ ∈ C ∞ ( S ) } , C ∞ ( T S ; Curl ) = { Curl ψ : ψ ∈ C ∞ ( S ) } . For s ≥
0, let C ∞ ,s ( T S ; Curl ) denote the closure of C ∞ ( T S ; Curl ) in the H s ( T S )norm. In particular, following [5] we introduce a simpler notation H = closure of C ∞ ( T S ; Curl ) in L ( T S ) = C ∞ , ( T S ; Curl ) ,V = closure of C ∞ ( T S ; Curl ) in H ( T S ) = C ∞ , ( T S ; Curl ) . Using the Gauss surface divergence theorem, for any scalar valued function v on S with Grad v ∈ L ( T S ), using (2.10), we have(2.15) (
Grad v, w ) = (cid:90) S Grad v · w dS = − (cid:90) S v · Div w dS = 0 , w ∈ V, and hence the unknown pressure can be eliminated from the first equation in (2.1)through the weak formulation.Following [16, Page 567], for the diffusion part of the NSE, we consider the Stokesoperator(2.16) A = Curl Curl (cid:98) x . Using (2.4) and (2.5), it is easy to see that the Stokes operator is the restric-tion of the vector Laplace-de Rham operator − ∆ on V ; A = − P Curl ∆ , where P Curl : L ( T S ) → H is the orthogonal projection onto the divergence-freetangent space.For each positive integer L = 1 , , . . . , the eigenvalue λ L and the correspondingeigenvectors of the Stokes operator A are given by(2.17) λ L = L ( L + 1) , Z L,m ( θ, ϕ ) = λ − / L Curl Y L,m ( θ, ϕ ) , m = − L, . . . , L, where Y L,m are the scalar orthonormal spherical harmonics of degree L , defined by(2.18) Y L,m ( θ, ϕ ) = (cid:20) (2 L + 1)4 π ( L − | m | )!( L + | m | )! (cid:21) / P mL (cos θ ) e imϕ , m = − L, . . . , L,
PSEUDOSPECTRAL METHOD FOR NSEs ON ROTATING SPHERES 5 with P mL being the associated Legendre polynomials so that Y L,m = ( − m Y L, − m .The spectral property AZ L,m = λ L Z L,m follows from the fact that Y L,m areeigenfunctions of the scalar Laplace-Beltrami operator − ∆ with eigenvalues λ L ,the definition of the Stokes operator A in (2.16), Z L,m in (2.17), (2.5), and (2.3).Since { Y L,m : L = 0 , , . . . ; m = − L, . . . , L } is an orthonormal basis for L ( S ), it iseasy to see that { Z L,m : L = 1 , . . . ; m = − L, . . . , L } is an orthonormal basis for H .Thus an arbitrary v ∈ H can be written as(2.19) v = ∞ (cid:88) L =1 L (cid:88) m = − L (cid:98) v L,m Z L,m , (cid:98) v L,m = (cid:90) S v · Z L,m dS = ( v , Z L,m ) . We consider a subset of H ,(2.20) D ( A s/ ) = (cid:40) v ∈ H : v = ∞ (cid:88) L =1 L (cid:88) m = − L (cid:98) v L,m Z L,m , ∞ (cid:88) L =1 L (cid:88) m = − L λ sL | (cid:98) v L,m | < ∞ (cid:41) , which is the divergence-free subset of the Sobolev space H s ( T S ). For every v ∈D ( A s/ ), we set(2.21) (cid:107) v (cid:107) H s ( T S ) = (cid:34) ∞ (cid:88) L =1 L (cid:88) m = − L λ sL | (cid:98) v L,m | (cid:35) / , and for v ∈ D ( A s/ ), we define(2.22) A s/ v := ∞ (cid:88) L =1 L (cid:88) m = − L λ s/ L (cid:98) v L,m Z L,m ∈ H. For a tangential vector field v on S , we define the Coriolis operator C ,(2.23) ( Cv )( (cid:98) x ) = ω ( (cid:98) x ) × v ( (cid:98) x ) = ω ( (cid:98) x )( (cid:98) x × v ) , ω ( (cid:98) x ) = 2Ω cos θ. To treat the nonlinear term in (2.1), we consider the trilinear form b on V × V × V ,defined as(2.24) b ( v , w , z ) = ( ∇ v w , z ) = (cid:90) S ∇ v w · z dS, v , w , z ∈ V. Using (2.6) and (2.10), for divergence free fields v , w , z , the trilinear form can bewritten as(2.25) b ( v , w , z ) = 12 (cid:90) S [ − v × w · Curl (cid:98) x z + Curl (cid:98) x v × w · z − v × Curl (cid:98) x w · z ] dS. Moreover [16, Lemma 2.1](2.26) b ( v , w , w ) = 0 , b ( v , z , w ) = − b ( v , w , z ) v , w , z ∈ V. Throughout the paper, the space V is equipped with the norm (cid:107) · (cid:107) V = ( A · , · ).Thus, using (2.4), (2.10), (2.16), and (2.25), a weak solution of the Navier-Stokesequations (2.1) is a vector field u ∈ L ([0 , T ]; V ) with u (0) = u that satisfies theweak form(2.27) ( u t , v ) + b ( u , u , v ) + ν ( Curl (cid:98) x u , Curl (cid:98) x v ) + ( Cu , v ) = ( f , v ) , v ∈ V. M. GANESH, Q. T. LE GIA, AND I. H. SLOAN
This weak formulation can be written in operator equation form on V ∗ , the adjointof V : Let f ∈ L ([0 , T ]; V ∗ ) and u ∈ H . Find a vector field u ∈ L ([0 , T ]; V ), with u t ∈ L ([0 , T ]; V ∗ ) such that(2.28) u t + ν Au + B ( u , u ) + Cu = f , u (0) = u , where the bilinear form B ( u , v ) ∈ V ∗ is defined by(2.29) ( B ( u , v ) , w ) = b ( u , v , w ) w ∈ V. In the subsequent error analysis, we need the following estimate for the nonlinearterm (see Lemma 6.1 in Appendix):(2.30) (cid:107) A − δ B ( u , v ) (cid:107) ≤ (cid:40) C (cid:107) A − δ u (cid:107)(cid:107) v (cid:107) ≤ C (cid:107) A / u (cid:107)(cid:107) v (cid:107) ,C (cid:107) u (cid:107)(cid:107) A − δ v (cid:107) ≤ C (cid:107) u (cid:107)(cid:107) A / v (cid:107) , δ ∈ (1 / , u , v ∈ V. In (2.30), as throughout the paper, C is a generic constant independent of u and v , (and the discretization parameter N introduced in Section 3). From (2.30) wededuce the weak Lipschitz continuity property(2.31) (cid:107) A − δ ( B ( v , v ) − B ( w , w )) (cid:107) ≤ C (cid:107) v − w (cid:107) , δ ∈ (1 / , , if (cid:107) A / v (cid:107) , (cid:107) A / w (cid:107) < C. The existence and uniqueness of the solution u ∈ L ([0 , T ]; V ) of the weak for-mulation (2.27) are discussed in [16, 17, 19]. A regular solution of the Navier-Stokesequations (2.1) on [0 , T ] is a tangential divergence-free velocity field u that satisfiesthe equation obtained by integrating in time the weak form (2.27), from t to t , foralmost every t , t ∈ [0 , T ]. In order to recall the existence, uniqueness, and Gevreyregularity of the regular solution, we need a few more additional details from [5].These are also needed as tools for analyzing our pseudospectral quadrature method.2.2. Gevrey regularity of regular solution.
The Gevrey class of functions oforder s > σ >
0, associated with the Stokes operator defined in (2.16),is denoted by G s/ σ and is defined as(2.32) G s/ σ := D ( A s/ e σA / ) ⊂ D ( A s/ ) . Using (2.20), the Gevrey space(2.33) G s/ σ = (cid:40) v ∈ D ( A s/ ) : v = ∞ (cid:88) L =1 L (cid:88) m = − L (cid:98) v L,m Z L,m , ∞ (cid:88) L =1 L (cid:88) m = − L λ sL e σλ / L | (cid:98) v L,m | < ∞ (cid:41) is a Hilbert space with respect to the inner product(2.34) (cid:104) v , w (cid:105) G s/ σ = ∞ (cid:88) L =1 L (cid:88) m = − L λ sL e σλ / L (cid:98) v L,m (cid:98) w L,m , v , w ∈ G s/ σ . First we recall the following result from [5, 29].
Theorem 2.1. If u ∈ D ( A s +1 / ) and f ∈ L ∞ ((0 , ∞ ); D ( A s e σ A / )) , for some s, σ > , then for all t > there exists a T ∗ > , depending only on ν, f , and (cid:107) A s +1 / u (cid:107) L ( T S ) , such that the NSE (2.1) on S have a unique regular solution u ( · , t ) and u ( · , t ) ∈ G s +1 / σ ( t ) , where σ ( t ) = min { t, T ∗ , σ } . PSEUDOSPECTRAL METHOD FOR NSEs ON ROTATING SPHERES 7
In addition, from the assumption and proof of Theorem 2.1 [5, page 355],(2.35) (cid:107) A s +1 / e σ ( t ) A / u ( t ) (cid:107) ≤ M , t > , where M depends on (cid:107) A s +1 / u (0) (cid:107) , sup t ≥ (cid:107) A s f ( · , t ) (cid:107) and ν but not on t . Thebound in (2.35) is useful for establishing the quality of approximation of the Stokesprojection of u in the next section (see Theorem 3.1). It is also convenient to havea similar bound for the time derivative of the Stokes projections with t in (2.35)replaced with certain complex times ζ ∈ C , to prove the power of approximationof the time derivative of the Stokes projection (see Theorem 3.2). To this end, weconsider the NSE extended to complex times ζ ,(2.36) d u dζ + ν Au + B ( u , u ) + Cu = f , Div u = 0 , u (0) = u , ζ ∈ C , on S, with standard complexification (see [9]) of all the spaces and operators introducedearlier. In the next theorem, we extend arguments used in [9] for the solution of theNSE on the plane to the case of the sphere. The arguments differ in an essentialway only for the nonlinear term. Theorem 2.2.
Let u ∈ D ( A s +1 / ) and f ∈ C ([0 , T ]; D ( A s e σ A / )) , with s ≥ / .Let the domain T be defined by T := { ζ = re iθ : 0 ≤ r ≤ T ; | θ | ≤ π/ } . We assume further that f ( · , ζ ) is analytic for ζ ∈ T and that (2.37) K := sup {(cid:107) A s e ψ ( r cos θ ) A / f ( · , ζ ) (cid:107) : ζ = re iθ ∈ T } < ∞ . Then there exists T ∗∗ > such that (2.38) (cid:107) A s +1 / e ψ ( r cos θ ) A / u ( ζ ) (cid:107) ≤ M , ζ ∈ T and | ζ | ≤ T ∗∗ , where M depends on (cid:107) A s +1 / u (0) (cid:107) and hence u ( · , ζ ) ∈ G s +1 / ψ ( r cos θ ) , where ψ ( x ) := min { x, T ∗∗ , σ } .Proof. Let ζ = re iθ with r > | θ | ≤ π/
4, and let u A ( ζ ) := A s +1 / e ψ ( r cos θ ) A / u ( ζ ) . The definition of ψ gives ddx ψ := ψ (cid:48) ( x ) ≤
1, and ψ ≤ σ , thus for fixed θ we find ddr u A ( ζ ) = ψ (cid:48) ( r cos θ ) cos θ A s +1 e ψ ( r cos θ ) A / u ( ζ )+ e iθ A s +1 / e ψ ( r cos θ ) A / d u dζ ( ζ ) . (2.39)Using (2.39) in 12 ddr (cid:107) u A ( ζ ) (cid:107) = (cid:60) (cid:18) ddr u A ( ζ ) , u A ( ζ ) (cid:19) , where (cid:60) ( · ) denotes the real-part function, we get12 ddr (cid:107) u A ( ζ ) (cid:107) = ψ (cid:48) ( r cos θ ) cos θ (cid:60) ( A / u A ( ζ ) , u A ( ζ ))+ (cid:60) e iθ ( A s +1 / e ψ ( r cos θ ) A / d u dζ ( ζ ) , u A ( ζ )) . (2.40) M. GANESH, Q. T. LE GIA, AND I. H. SLOAN
Using (2.36) for the last term in (2.40) together with the fact that( A s e ψ ( r cos θ ) A / Cu , A s +1 e ψ ( r cos θ ) A / u ) = 0(see [5, Lemma 1]), we find12 ddr (cid:107) u A ( ζ ) (cid:107) + ν cos θ (cid:107) A / u A ( ζ ) (cid:107) = ψ (cid:48) ( r cos θ ) cos θ ( A / u A , u A ) − (cid:60) e iθ ( A s +1 / e ψ A / B ( u , u ) , u A )+ (cid:60) e iθ ( A s e ψ A / f , A / u A ) , (2.41)where in (2.41) and below we write u A = u A ( ζ ) , u = u ( ζ ) , f = f ( ζ ) , ψ = ψ ( r cos θ ).From [5, Lemma 2], we have (with p = max { − s, / } = 7 /
4, since s ≥ / | ( A s +1 / e ψ A / B ( u , u ) , u A ) | ≤ C (cid:107) u A (cid:107) − p (cid:107) A / u A (cid:107) p . Applying ψ (cid:48) ≤
1, the Cauchy-Schwarz inequality, (2.42) and Young’s inequality( ab ≤ a q /q + b q (cid:48) /q (cid:48) with 1 /q + 1 /q (cid:48) = 1) with q = 2 and q = 2 / (2 − p ) in (2.41), weget 12 ddr (cid:107) u A (cid:107) + ν cos θ (cid:107) A / u A (cid:107) ≤ cos θ (cid:107) A / u A (cid:107)(cid:107) u A (cid:107) + C (cid:107) u A (cid:107) − p (cid:107) A / u A (cid:107) p + (cid:107) A s e ψ A / f (cid:107)(cid:107) A / u A (cid:107)≤ ν cos θ (cid:107) A / u A (cid:107) + cos θν (cid:107) u A (cid:107) + C − p (2 − p ) p (cid:16) pν cos θ (cid:17) p − p (cid:107) u A (cid:107) − p )2 − p + ν cos θ (cid:107) A / u A (cid:107) + 1 ν cos θ (cid:107) A s e ψ A / f (cid:107) + ν cos θ (cid:107) A / u A (cid:107) . Therefore, ddr (cid:107) u A (cid:107) ≤ θν (cid:107) u A (cid:107) + C (cid:18) ν cos θ (cid:19) p − p (cid:107) u A (cid:107) − p )2 − p + 2 ν cos θ (cid:107) A s e ψ A / f (cid:107) ≤ θν (1 + (cid:107) u A (cid:107) ) + C (cid:18) ν cos θ (cid:19) p − p (1 + (cid:107) u A (cid:107) ) − p − p +2 ν cos θ (cid:107) A s e ψ A / f (cid:107) . (2.43)Using | θ | ≤ π/
4, and − p − p = 5, in (2.43), we get(2.44) ddr (cid:107) u A (cid:107) ≤ C (1 + (cid:107) u A (cid:107) ) + 2 √ ν K. With | θ | ≤ π/ y ( r ) = 1 + (cid:107) u A (cid:107) . Then on using (2.44) and y ≥ ddr y ≤ Cy , and hence on integrating the inequality we find that y ( r ) ≤ y (0) PSEUDOSPECTRAL METHOD FOR NSEs ON ROTATING SPHERES 9 provided that 0 ≤ r ≤ C y (0)) = 1564 C (cid:107) A s +1 / u (0) (cid:107) ) . By setting T ∗∗ := 1564 C (cid:107) A s +1 / u (0) (cid:107) ) , M = 1 + 2 (cid:107) A s +1 / u (0) (cid:107) , we deduce that (2.38) holds for 0 ≤ r ≤ T ∗∗ . (cid:3) We can extend the bound (2.38) to a larger domain which contains the interval[0 , T ] using the following property of the NSE solution on the sphere [29]:(2.45) (cid:107) A s +1 / u ( · , t ) (cid:107) ≤ M for all t ∈ [0 , T ] , where the constant M depends only on (cid:107) A s +1 / u (cid:107) , sup ≤ t ≤ T (cid:107) A s f ( · , t ) (cid:107) and ν but not on T . Theorem 2.3.
Suppose u and f satisfy all the conditions in the domain T as inTheorem 2.2. Then (2.46) (cid:107) A s +1 / e ψ ( r cos θ ) A / u ( ζ ) (cid:107) ≤ M , ζ ∈ T and | Im ζ | ≤ T ∗∗ / √ , where M := 1 + 2 M , T ∗∗ depends on M , and hence u ( · , ζ ) ∈ G s +1 / ψ ( r cos θ ) for all ζ ∈ T with | Im ζ | ≤ T ∗∗ / √ .Proof. We proceed as in the proof of Theorem 2.2 to obtain the ordinary differentialequation ddr y ≤ Cy , where y ( r ) = 1 + (cid:107) A s +1 / e ψ ( r cos θ ) A / u ( re iθ ) (cid:107) . On integrating the inequality wefind that y ( r ) ≤ y (0)provided that 0 ≤ r ≤ C (cid:107) A s +1 / u (0) (cid:107) ) . We define T ( ρ ) := 1564 C ρ ) , ρ ≥ . If T ( (cid:107) A s +1 / u (0) (cid:107) ) ≥ T we have finished the proof. Otherwise, we let T ∗∗ = T ( M ) , where M is given in (2.45). For ζ = re iθ ∈ T , 0 ≤ r ≤ T ∗∗ ,(2.47) (cid:107) A s +1 / e ψ ( r cos θ ) A / u ( re iθ ) (cid:107) ≤ (cid:107) A s +1 / u (0) (cid:107) . Consequently, (2.46) holds for 0 ≤ r ≤ T ∗∗ with M = 1 + 2 (cid:107) A s +1 / u (0) (cid:107) .Next we consider the case ζ = T ∗∗ + re iθ , with r ∈ [0 , T ∗∗ ]. We define, for | θ | ≤ π/ v ( re iθ ) := u ( T ∗∗ + re iθ ) , r ∈ [0 , T ∗∗ ] . Using (cid:107) A s +1 / v (0) (cid:107) ≤ M we can apply the previous arguments to obtain (2.46)(with u replaced with v ) for 0 ≤ r ≤ T ∗∗ . We complete the proof and obtain thebound (2.46) by repeating the last argument n times, where n = (cid:100) T /T ∗∗ (cid:101) . (cid:3) Finite dimensional spaces and Stokes projections
Throughout the remainder of the paper, with s and σ as in Theorem 2.1 andTheorem 2.2, we assume that(3.1) u ∈ D ( A s +1 / ) , f ∈ C ([0 , T ]; D ( A s +1 / e σ A / )) , s ≥ / . Natural finite dimensional spaces (depending on a parameter
N >
0) in whichto seek approximations to u ( t ) are(3.2) V N := span { Z L,m : L = 1 , . . . , N ; m = − L, . . . , L } . The dimension of V N is N + 2 N . Let Π N : H → V N be the orthogonal projectionwith respect to the L ( T S ) inner product defined by(3.3) Π N ( v ) = N (cid:88) L =1 L (cid:88) m = − L (cid:98) v L,m Z L,m . Lemma 3.1.
Let α > be given. If v ∈ D ( A α/ ) then (3.4) (cid:107) v − Π N v (cid:107) ≤ N − α (cid:107) v (cid:107) H α ( T S ) . Proof.
Using (2.17), (2.19), and (2.21) we get (cid:107) v − Π N v (cid:107) = ∞ (cid:88) L = N +1 L (cid:88) m = − L | (cid:98) v L,m | ≤ N − α ∞ (cid:88) L = N +1 L (cid:88) m = − L λ αL | (cid:98) v L,m | ≤ N − α (cid:107) v (cid:107) H α ( T S ) . (cid:3) In particular, using (3.1), we get(3.5) (cid:107) f − Π N f (cid:107) ≤ N − (2 s +1) (cid:107) f (cid:107) H s +1 ( T S ) , t ∈ [0 , T ] . For computer implementation, the Fourier coefficients in (3.3) and all Galerkintype integrals in computational schemes for the NSE need to be approximated bycubature rules on the sphere, leading to a pseudospectral method. To this end, fora continuous scalar field ψ on S , we consider a Gauss-rectangle quadrature sum Q M ( ψ ) with quadrature points { (cid:98) ξ p,q = p ( θ p , φ q ) } and positive weights w p of theform(3.6) Q M ( ψ ) := 2 πM M (cid:88) q =1 M/ (cid:88) p =1 w p ψ ( (cid:98) ξ p,q ) = 2 πM M (cid:88) q =1 M/ (cid:88) p =1 w p ψ ( θ p , φ q ) , where M ≥ w p and cos θ p for p = 1 , . . . , M/ − ,
1] and φ q = 2 qπ/M , q = 1 , . . . , M . The rule(3.6) is exact when ψ is a polynomial of degree M − S , that is, Q M ψ = (cid:90) S ψ dS, ψ ∈ P M − . Hence corresponding to (2.8) and (2.9), we define discrete inner products for scalarand vector fields on the unit sphere as(3.7) ( v , v ) M = Q M ( v v ) , ( v , v ) M = Q M ( v · v ) . The choice of M is very important; we choose M such that all Galerkin integralswith polynomial terms in our scheme are evaluated exactly. In particular, with the PSEUDOSPECTRAL METHOD FOR NSEs ON ROTATING SPHERES 11 unknown tangential divergence-free velocity field sought in the polynomial space V N , and knowing that the NSE nonlinearity is quadratic, we choose M such that(3.8) ( B ( v , w ) , z ) = ( B ( v , w ) , z ) M , v , w , z ∈ V N . This holds, for example, if M = 3 N + 2. We define a computable counterpart of(3.3), using L N : H ∩ C ( T S ) → V N , a discrete orthogonal projection with respectto the M / L N ( v ) = N (cid:88) L =1 L (cid:88) m = − L (cid:98) v L,m,M Z L,m , (cid:98) v L,m,M = Q M ( v · Z L,m ) = ( v , Z L,m ) M . With M chosen to satisfy (3.8), it is easy to see that(3.10) L N ( v ) = Π N ( v ) = v , v ∈ V N , and for v ∈ H ∩ C ( T S ) and w ∈ H ∩ C k ( T S ),(3.11) (cid:107) L N ( v ) (cid:107) ≤ C (cid:107) v (cid:107) ∞ , (cid:107) w − L N ( w ) (cid:107) ≤ CN − k (cid:107) w (cid:107) C k ( T S ) , where the last two inequalities follow from simple arguments used in Theorem 13and Lemma 14 of [24]. In particular, since D ( A s +1 / ) ⊂ H ∩ C s ( T S ), for aninteger 2 s , using (3.1)(3.12) (cid:107) f − L N ( f ) (cid:107) ≤ CN − s (cid:107) f (cid:107) C s ( T S ) , t ∈ [0 , T ] . Next we consider the Stokes projection in V N of the exact unique regular solution u ( t ) := u ( ., t ) of (2.1). For each fixed t , the Stokes projection (cid:101) u N ∈ V N of u isdefined by(3.13) ( A (cid:101) u N , v ) = ( Au , v ) , v ∈ V N . Since Π N A = AΠ N , it follows that (cid:101) u N = Π N ( u ). Following standard techniquesin finite element analysis, the Stokes projection of u plays an important role as acomparison function in the main analysis in the next section. Theorem 3.1.
Let u and f satisfy (3.1) . Then (3.14) (cid:107) u − (cid:101) u N (cid:107) ≤ Cλ − s − / N +1 e − σ ( t ) λ / N +1 ≤ CN − s − e − σ ( t ) N , t ∈ [0 , T ] , where σ ( t ) is as in Theorem 2.1.Proof. Using the fact that (cid:101) u N = Π N u , we have (cid:107) u − (cid:101) u N (cid:107) = (cid:88) L>N (cid:88) | m |≤ L | (cid:98) u L,m | ≤ λ − s − N +1 e − σ ( t ) λ / N +1 (cid:88) L>N (cid:88) | m |≤ L λ s +1 L e σ ( t ) λ / L | (cid:98) u L,m | ≤ λ − s − N +1 e − σ ( t ) λ / N +1 (cid:107) A s +1 / e σ ( t ) A / u ( t ) (cid:107) ≤ Cλ − s − N +1 e − σ ( t ) λ / N +1 where in the last step we used (2.35). The last inequality in (3.14) follows from thefact that N ≤ λ N +1 = ( N + 1)( N + 2). (cid:3) Theorem 3.2.
Let u , f satisfy (3.1) . We assume further that f is analytic in T and (2.37) holds. Then for t ∈ (0 , T ) , (3.15) (cid:13)(cid:13)(cid:13)(cid:13) ddt ( u − (cid:101) u N ) (cid:13)(cid:13)(cid:13)(cid:13) ≤ Cλ − s − / N +1 e − ψ ( t ) λ / N +1 ≤ CN − s − e − ψ ( t ) N , where ψ ( t ) = min { (1 − / √ t, T ∗∗ , σ } , and T ∗∗ is as in Theorem 2.3.Proof. Let t ∈ (0 , T ) be fixed. Let p N ( ζ ) = ( u − (cid:101) u N )( ζ ) be the standard com-plexification of u − (cid:101) u N at ζ = re iθ . Using Theorem 2.3 and the Cauchy integralformula, d p N ( t ) dt = 12 πi (cid:90) Γ p N ( ζ )( t − ζ ) dζ, where for t >
0, Γ is a circle in the ζ plane with center ( t,
0) and radiusmin { t/ √ , T ∗∗ / √ , T − t } , a condition that ensures that ζ = re iθ ∈ Γ lies in theregion T with | Im ζ | ≤ T ∗∗ / √
2. Using the fact that (cid:101) u N = Π N u , for ζ = re iθ ∈ Γwe have, from Theorem 2.3, (cid:107) p N ( ζ ) (cid:107) = (cid:107) u ( ζ ) − (cid:101) u N ( ζ ) (cid:107) = (cid:88) L>N (cid:88) | m |≤ L | (cid:98) u L,m | ≤ λ − s − N +1 e − ψ ( r cos θ ) λ / N +1 (cid:88) L>N (cid:88) | m |≤ L λ s +1 L e ψ ( r cos θ ) λ / L | (cid:98) u L,m | ≤ λ − s − N +1 e − ψ ( r cos θ ) λ / N +1 (cid:107) A s +1 / e ψ ( r cos θ ) A / u ( ζ ) (cid:107) ≤ Cλ − s − N +1 e − ψ ( r cos θ ) λ / N +1 ≤ CN − s +1) e − ψ ( r cos θ ) N . (3.16)For ζ = re iθ ∈ Γ it is easily seen that r cos θ ≥ (1 − / √ t , and hence that(3.17) ψ ( r cos θ ) ≥ min { (1 − / √ t, T ∗∗ , σ } =: ψ ( t ) . On using (3.16) in (3.14) we get (cid:13)(cid:13)(cid:13)(cid:13) ddt ( u − (cid:101) u N ) ( t ) (cid:13)(cid:13)(cid:13)(cid:13) ≤ π (cid:90) Γ (cid:107) p N ( ζ ) (cid:107)| t − ζ | dζ ≤ CN − s − e − ψ ( t ) N . (cid:3) A pseudospectral quadrature method
We are now ready to describe, analyze, and implement a spectrally accuratescheme to compute approximate solutions of the NSE (2.1) in V N , through its weakformulation (2.27). The task is then to compute u N ( · , t ) ∈ V N for t ∈ [0 , T ] with u N (0 , (cid:98) x ) = L N u ( (cid:98) x ), (cid:98) x ∈ S , satisfying the (spatially) discrete system of ordinarydifferential equations(4.1) ddt ( u N , v ) M + b ( u N , u N , v ) M + ν ( Curl (cid:98) x u N , Curl (cid:98) x v ) M + ( Cu N , v ) M = ( f , v ) M , for all v ∈ V N and prove that the scheme is spectrally accurate with respect to theparameter N (that is, converges with rate determined by the smoothness of thegiven data), and demonstrate the theory with numerical experiments. PSEUDOSPECTRAL METHOD FOR NSEs ON ROTATING SPHERES 13
Using (3.2), the exactness properties of the discrete inner product, (3.8), (2.10),(2.16), (2.23), and (2.29), the system (4.1) can be written as (cid:18) ddt u N + ν Au N + B ( u N , u N ) + Cu N , Z L,m (cid:19) = ( f , Z L,m ) M ,L = 1 , · · · , N ; m = − L, · · · , L. (4.2)4.1. Stability and convergence analysis.
First we establish the stability of theapproximate solution u N of (4.24). That is, similar to [6, Proposition 9.1], we provethat max t ∈ [0 ,T ] (cid:107) u N (cid:107) V is uniformly bounded with the bound depending only on theinitial data, forcing function, and the viscosity term in (4.24). Theorem 4.1.
Let u and f satisfy (3.1). Let N ≥ be an integer. Let u N be thesolution of the pseudospectral quadrature Galerkin system (4.1) . Then there existsa constant C depending on ν, (cid:107) u (cid:107) V and (cid:107) f (cid:107) ∞ := max t ∈ [0 ,T ] (cid:107) f ( t ) (cid:107) C ( T S ) so that max t ∈ [0 ,T ] (cid:107) u N (cid:107) V ≤ C. Proof.
The proof follows by repeating the arguments described in the first fourpages of [6, Section 9] (proving [6, Proposition 9.1]), provided that we establish [6,Inequalities (9.3) and (9.13)] for our system (4.1) on the spherical surface withadditional Coriolis term and quadrature approximations.Using (2.26), (2.29) and the exactness of the quadrature rule (given by (3.7)-(3.8)), we have ( B ( u N , u N ) , u N ) M = 0. Using (4.27), the symmetry of the coeffi-cients of u N in (4.24) and the exactness of the quadrature, we get(4.3) ( Cu N , u N ) M = ( Cu N , u N ) = ( − i ) N (cid:88) L =1 λ − L (cid:88) | m |≤ L m | α L,m | = 0 . By taking v to be u N in (4.1) and using (cid:107) u N (cid:107) = ( u N , u N ) M , (cid:107) u N (cid:107) V = ( Au N , u N ) M ,(4.3), Young’s inequality and the fact that all the eigenvalues λ J of A (correspond-ing to eigenvectors in u N ) satisfy λ J ≥ λ = 2 , J = 1 , · · · , N , we obtain12 ddt (cid:107) u N (cid:107) + ν (cid:107) u N (cid:107) V = ( f , u N ) M ≤ (cid:107) f (cid:107) ∞ (cid:107) u N (cid:107) ≤ (cid:107) f (cid:107) ∞ ν + ν (cid:107) u N (cid:107) ≤ (cid:107) f (cid:107) ∞ ν + ν (cid:107) u N (cid:107) V , Hence, for our discrete system (4.1), we obtain [6, Inequality (9.3)]:(4.4) ddt (cid:107) u N (cid:107) + ν (cid:107) u N (cid:107) V ≤ (cid:107) f (cid:107) ∞ νλ . Again using (4.27), the exactness of the quadrature rule, eigenfunction propertiesof A , and the symmetry of the coefficients of u N in (4.24), we get(4.5) ( Cu N , Au N ) = ( − i ) N (cid:88) L =1 (cid:88) | m |≤ L m | α L,m | = 0 , By taking v to be Au N in equation (4.1), and using (cid:107) Au N (cid:107) = ( Au N , Au N ) M ,(3.8), and (4.3), we obtain(4.6) 12 ddt (cid:107) u N (cid:107) V + ν (cid:107) Au N (cid:107) + b ( u N , u N , Au N ) = ( f , Au N ) M . Using [16, Lemma 3.1], b ( u N , u N , Au N ) = 0. The term ( f , Au N ) M can be esti-mated by using the exactness of the quadrature and Young’s inequality:( f , Au N ) M ≤ (cid:107) f (cid:107) ∞ (cid:107) Au N (cid:107) ≤ (cid:107) f (cid:107) ∞ ν + ν (cid:107) Au N (cid:107) Hence we obtain a stronger version of [6, Inequality (9.13)] for our quadraturediscrete scheme (4.1): ddt (cid:107) u N (cid:107) V + ν (cid:107) Au N (cid:107) ≤ (cid:107) f (cid:107) ∞ νλ . Thus, the result follows from arguments in [6, Page 74-77]. (cid:3)
Next, using Theorem 3.1, 3.2, and 4.1, we prove the spectral convergence of thesolution u N of (4.2) to the solution of u of (2.27). Theorem 4.2.
Let u , f satisfy (2.37) , (3.1) . Then there exists a T > , de-pending only on ν, f , u and the uniform bound in (2.45) (and hence there exists < µ ( t ) < min { t, T , σ } ) such that for all t ∈ (0 , T ) , (4.7) (cid:107) u − u N (cid:107) ≤ C (cid:104) N − s − e − µ ( t ) N + (cid:107) Π N f − L N f (cid:107) (cid:105) . In particular, with s being an integer (4.8) (cid:107) u − u N (cid:107) ≤ CN − s . Proof.
Let w N = (cid:101) u N − u N , where the comparison function (cid:101) u N is the solution of(3.13). Since u − u N = p N + w N , where p N = u − (cid:101) u N , in view of Theorem 3.1and 3.2, existence of T and µ ( t ) follows and it is sufficient to show that (cid:107) w N (cid:107) ≤ C (cid:2) N − s − e − µ ( t ) N + (cid:107) Π N f − L N f (cid:107) (cid:3) . For any v ∈ V N , using (2.28), (3.13), and(4.2),(( w N ) t , v ) + ν ( Aw N , v ) + ( Cw N , v )= (( (cid:101) u N ) t , v ) + ν ( A (cid:101) u N , v ) + ( C (cid:101) u N , v ) − (( u N ) t , v ) − ν ( Au N , v ) − ( Cu N , v )= (( (cid:101) u N ) t , v ) + ν ( A (cid:101) u N , v ) + ( C (cid:101) u N , v ) − ( f , v ) M + ( B ( u N , u N ) , v )= (( (cid:101) u N ) t , v ) + ν ( Au , v ) + ( C (cid:101) u N , v ) − ( f , v ) M + ( B ( u N , u N ) , v )= (( (cid:101) u N − u ) t , v ) + ( f , v ) − ( f , v ) M + ( C (cid:101) u N − Cu , v ) + ( B ( u N , u N ) − B ( u , u ) , v ) . Using the orthogonal projection Π N in (3.2), we can write this relation in functionalform as d w N dt = ( − ν A − C ) w N − Π N [( p N ) t + Cp N ]+ Π N f − L N f + Π N [ B ( u N , u N ) − B ( u , u )] . Integrating with respect to t and using w N (0) = 0, we have w N ( t ) = (cid:90) t e − ( t − s )( ν A + C ) (cid:20) − Π N ( dds p N + Cp N ) + Π N f − L N f (cid:21) ( s ) ds + (cid:90) t e − ( t − s )( ν A + C ) Π N [ B ( u N , u N ) − B ( u , u )] ( s ) ds. (4.9)Let(4.10) R N ( (cid:15), t − s ) = (cid:107) ν (cid:15) Π N A (cid:15) e − ( t − s )( ν A + C ) (cid:107) . PSEUDOSPECTRAL METHOD FOR NSEs ON ROTATING SPHERES 15
On using (2.31), with δ ∈ (1 / , Π N , we get (cid:107) e − ( t − s )( ν A + C ) Π N [ B ( u N , u N ) − B ( u , u )] (cid:107) (4.11) ≤ R N ( δ, t − s ) ν δ (cid:107) A − δ Π N ( B ( u N , u N ) − B ( u , u )) (cid:107) ≤ CR N ( δ, t − s ) (cid:107) u N − u (cid:107) . Taking norms and using (cid:107) u N − u (cid:107) ≤ (cid:107) w N (cid:107) + (cid:107) p N (cid:107) together with (4.10), and (4.12)in (4.9), we obtain (cid:107) w N ( t ) (cid:107) ≤ (cid:20) (cid:107) ddt p N (cid:107) + (cid:107) Cp N (cid:107) + (cid:107) Π N f − L N f (cid:107) (cid:21) (cid:90) t R N (0 , t − s ) ds + C (cid:90) t R N ( δ, t − s )( (cid:107) p N ( s ) (cid:107) + (cid:107) w N ( s ) (cid:107) ) ds. Using Gronwall’s inequality, we obtain for each t ∈ [0 , T ], (cid:107) w N ( t ) (cid:107) ≤ C (cid:26)(cid:20) (cid:107) ddt p N (cid:107) + (cid:107) Cp N (cid:107) + (cid:107) Π N f − L N f (cid:107) (cid:21) (cid:90) t R N (0 , t − s ) ds + (cid:107) p N (cid:107) (cid:90) t R N ( δ, t − s ) ds (cid:27) . (4.12)For (cid:15) ∈ [0 , (cid:90) t R N ( (cid:15), t − s ) ds = (cid:90) t R N ( (cid:15), r ) dr. Using also (4.10) and (4.27), R N ( (cid:15), r ) ≤ max ≤ L ≤ N ; | m |≤ L | ( νλ L ) (cid:15) e − νλ L r e − irmλ − / L | ≤ max z ∈ [ νλ ,νλ N ] z (cid:15) e − rz . Thus(4.14) R N ( r ) ≤ ( νλ N ) (cid:15) e − νλ N r if r ≤ (cid:15)/ ( νλ N ) ,(cid:15) (cid:15) e − (cid:15) r − (cid:15) if (cid:15)/ ( νλ N ) ≤ r ≤ (cid:15)/ ( νλ ) , ( νλ ) (cid:15) e − νλ r if r ≥ (cid:15)/ ( νλ ) . With(4.15) I = [0 , (cid:15)/ ( νλ N )] ∩ [0 , t ] , I = [ (cid:15)/ ( νλ N ) , (cid:15)/ ( νλ )] ∩ [0 , t ] , I = [ (cid:15)/ ( νλ ) , t ] ∩ [0 , t ] , the interval of integration in (4.13) can be subdivided into these three intervals. Inparticular, using (4.14) and (4.15), we get(4.16) (cid:90) I R N ( r ) dr ≤ (cid:90) I ( νλ N ) (cid:15) e − νλ N r dr ≤ (1 − e − (cid:15) )( νλ N ) − (cid:15) ≤ C, (4.17) (cid:90) I R N ( r ) dr ≤ (cid:90) I (cid:15) (cid:15) e − (cid:15) r − (cid:15) dr ≤ (cid:15)e − (cid:15) (1 − (cid:15) ) (cid:34)(cid:18) νλ (cid:19) (1 − (cid:15) ) − (cid:18) νλ N (cid:19) (1 − (cid:15) ) (cid:35) ≤ C, and(4.18) (cid:90) I R N ( r ) dr ≤ (cid:90) I ( νλ ) (cid:15) e νλ r dr ≤ ( e − (cid:15) − e − νλ t )( νλ ) − (cid:15) ≤ C. Using (4.16)- (4.18) in (4.13), we get(4.19) (cid:90) t R N ( (cid:15), t − s ) ds ≤ C, (cid:15) ∈ [0 , . Substituting this in (4.12), we get(4.20) (cid:107) w N ( t ) (cid:107) ≤ C (cid:20) (cid:107) ddt p N (cid:107) + (cid:107) Cp N (cid:107) + (cid:107) p N (cid:107) + (cid:107) Π N f − L N f (cid:107) (cid:21) . The term (cid:107) Cp N (cid:107) in (4.20) can be simplified using (2.23) and the fact that p N istangential (and hence (cid:98) x · p N ( (cid:98) x ) = 0),[ (cid:98) x × p N ( (cid:98) x )] · (cid:104) ( (cid:98) x × p N ( (cid:98) x ) (cid:105) = [ (cid:98) x · (cid:98) x ] (cid:104) p N ( (cid:98) x ) · p N ( (cid:98) x ) (cid:105) = (cid:107) p N (cid:107) , and hence(4.21) (cid:107) w N ( t ) (cid:107) ≤ C (cid:20) (cid:107) ddt p N (cid:107) + (cid:107) p N (cid:107) + (cid:107) Π N f − L N f (cid:107) (cid:21) . Hence from Theorem 3.1 and 3.2, for all t ∈ (0 , T ) we have(4.22) (cid:107) u − u N (cid:107) ≤ C (cid:104) N − s − e − µ ( t ) N + (cid:107) Π N f − L N f (cid:107) (cid:105) . In particular, using (3.5) and (3.12) with 2 s being an integer, we get(4.23) (cid:107) Π N f − L N f (cid:107) ≤ (cid:107) Π N f − f (cid:107) + (cid:107) f − L N f (cid:107) ≤ CN − s . Now the result (4.8) follows from (4.22) and (4.23). (cid:3)
Adaptive and fast implementation of the pseudospectral method.
Having established spectrally accurate convergence of the spatially discrete scheme(4.2), for implementation of the scheme (4.2) to simulate stable and accurate solu-tions of (2.1) and compare with benchmark random flow simulations in the litera-ture, we need to discretize the time derivative operator ddt in (4.2). Further, at each discrete time step we develop a (FFT-based) fast evaluation technique to set upthe resulting fully discrete nonlinear system with spatial O ( N ) complexity. Firstwe consider discretization of ddt in (4.2).In order to the make (4.2) fully discrete in space and time, and hence computethe N + 2 N unknown time-dependent coefficients in the representation of thetangential divergence-free approximate real velocity vector field(4.24) u N ( (cid:98) x , t ) := N (cid:88) L =1 (cid:88) | m |≤ L α L,m ( t ) Z L,m ( (cid:98) x ) , α L,m = α L, − m , α L,m (0) = ( u , Z L,m ) M , for (cid:98) x ∈ S and t ≥
0, the standard fixed-time-step backward-Euler (or Crank-Nicolson) Galerkin approach could be used in (4.1), leading to a first-order (or asecond-order, respectively) in time non-adaptive scheme [13]. However, due to thecomplicated unknown flow behavior of the NSE solutions, when the initial statesare random, it is more efficient instead to integrate (4.1) using a combination ofmulti-order integration formulas that allow adaptive choice of time step, leading tocomputation of solutions with a specified accuracy in time. In this paper we followthe latter approach.For implementation purposes, we first need to substitute (4.24) in (4.2). Wewrite the resulting N + 2 N -dimensional system of ordinary differential equations PSEUDOSPECTRAL METHOD FOR NSEs ON ROTATING SPHERES 17 (ODE), for the unknown N + 2 N time dependent coefficients of u N ( (cid:98) x , t ) in (4.24),as(4.25) ddt α ( t ) = F ( t, α ( t )) . It is well known that such nonlinear ODE systems are stiff, and hence it isimportant to use time discretization techniques with a large stability region [14].Further, it is important to use high-order implicit formulas whenever possible, butthe high-order discretization formulas are appropriate only at those time stepswhere the unknown exact solution is smooth.For practical problems such as the Navier-Stokes equations (with initial randomstate) where the necessary spatial discretization (at each time step) is expensive,it is important to optimize computing time by simulating only up to a requiredaccuracy by choosing adaptive discrete time steps. Unlike the adaptive spatialmesh for elliptic PDEs based on the a posteriori estimates, adaptive time steps canbe computed by comparing numerical solutions obtained using two distinct orderformulas [14] and hence simulation using multi-order formulas is appropriate.In particular, for practical realization of large stiff nonlinear ODE systems, multi-order implicit backward differentiation formulas (BDF) and their generalizationssuch as the numerical differential formulas (NDF) are most appropriate. The im-plicit NDF formula of order p (NDFp) with a parameter κ p (so that κ p = 0 cor-responds to BDFp) for the system (4.25), with α n ≈ α ( t n ) and ∇ m denoting the m -th Newton backward difference operator, is(4.26) p (cid:88) m =1 m ∇ m α n +1 = h F ( t n +1 , α n +1 ) + κ p ∇ p +1 α n +1 p (cid:88) j =1 j . It is well known [14] that BDFp (and hence NDFp) are unstable for p >
6, andfor p = 6 the stability region is small and hence not practically useful in our case.Further the celebrated Dahlquist barrier [14] implies that BDFp (and hence NDFp)cannot be absolutely stable [that is, A ( α )-stable with α = 90 ◦ ] for p > p = 1 , , , , κ p = − . , − / , − . , − . ,
0) and theseare A ( α p )-stable, with respective α p = 90 ◦ , ◦ , ◦ , ◦ , ◦ . For p = 1 , , , , p = 3 , α p = 86 ◦ , ◦ ) and thesame stability angle for p = 1 , , F ( t n +1 , α n +1 ) in (4.26) and it is important to have an efficientmethod to set up the spatial part of the nonlinear system (4.2). Using the spectralproperties of the Stokes operator A given by (2.17), the linear second term in (4.2)is trivial to evaluate using the diagonal matrix consisting of the eigenvalues of A .The Coriolis term can be evaluated similarly using the identity [5, Equation (24)](4.27) ( C Curl Y J,k , Z L,m ) = (2Ω cos θ (cid:98) x × Curl Y J,k , Z L,m ) = − i mλ / L δ L,J δ k,m . For the first component in the nonlinear third term in (4.2), we use (2.7) to write B ( Z R,s , Z J,k ) = − √ λ R λ J P Curl (∆ Y R,s
Grad Y J,k ) , (4.28) ( B ( Z R,s , Z J,k ) , Z L,m ) = (cid:114) λ R λ J λ L ( Y R,s
Grad Y J,k , (cid:98) x × Grad Y L,m ) . (4.29)It is convenient to write Grad Y J,k and (cid:98) x × Grad Y L,m in terms of expressionssimilar to those in (2.18). Such explicit representations are also useful for theefficient evaluation of the N Fourier coefficients ( f , Z L,m ) M of the source term in(4.2), and eventually for the computation of the vorticity field.In order to express the tangential (and normal, needed for computing the ap-proximate vorticity from u N ) vector harmonics as a linear combination of the scalarharmonics (2.18), we first recall, from the classical quantum mechanics literature(see, for example, [30]), the covariant spherical basis vectors(4.30) e +1 = − √ , , T + i [0 , , T ) , e = [0 , , T , e − = 1 √ , , T − i [0 , , T ) , and the Clebsch-Gordan coefficients(4.31) C j,mj ,m ,j ,m := ( − ( m + j − j ) (cid:112) j + 1 (cid:18) j j jm m − m (cid:19) , where (cid:18) a b cα β γ (cid:19) are the Wigner 3j-symbols given, for example, by the Racahformula, (cid:18) a b cα β γ (cid:19) = ( − ( a − b − γ ) (cid:112) T ( abc ) (cid:112) ( a + α )!( a − α )!( b + β )!( b − β )!( c + γ )!( c − γ )! × (cid:88) t ( − t t !( c − b + t + α )!( c − a + t − β )!( a + b − c − t )!( a − t − α )!( b − t + β )! , where the sum is over all integers t for which the factorials in the denominatorall have nonnegative arguments. In particular, the number of terms in the sum is1 + min { a ± α, b ± β, c ± γ, a + b − c, b + c − a, c + a − b } . The triangle coefficient T ( abc ) is defined by T ( abc ) = (cid:20) ( a + b − c )!( a − b + c )!( − a + b + c )!( a + b + c + 1)! (cid:21) . Below, we require C j,mj ,m ,j ,m only for some j , m ∈ {− , , } , and using varioussymmetry and other known properties (such as C j,mj ,m ,j ,m = 0 unless the condi-tions | j − j | ≤ j ≤ j + j and m + m = m hold) of Wigner 3j-symbols, thesecoefficients can be efficiently pre-computed and stored.In our computation, we used the following basis functions for the tangential vec-tor fields: ( i ) Grad Y L,m , ( ii ) (cid:98) x × Grad Y J,k . For the vorticity components of u N , inaddition we used ( iii ) Vort Z
J,m = Curl (cid:98) x × Z J,m = λ / J (cid:98) x Y J,m = − λ − / J (cid:98) x ∆ Y J,m .In particular, using (4.24) our approximation to the vorticity in (2.2), for a fixed t ≥ (cid:98) x ∈ S , is(4.32) Vort u N ( (cid:98) x , t ) = Curl (cid:98) x u N ( (cid:98) x ) = (cid:98) x ∆Ψ N ( (cid:98) x ) , PSEUDOSPECTRAL METHOD FOR NSEs ON ROTATING SPHERES 19 where Ψ N ( (cid:98) x , t ) = − N (cid:88) L =1 (cid:88) | m |≤ L λ − / L α L,m ( t ) Y L,m ( (cid:98) x ) . To facilitate easy application of fast transforms to evaluate these functions atthe M = O ( N ) quadrature points { (cid:98) ξ p,q = p ( θ p , φ q ) } , we represent these threetypes of fields first as a linear combination of the covariant vectors in (4.30): Grad Y L,m = B +1 ,L,m e +1 + B ,L,m e + B − ,L,m e − , (4.33) ( (cid:98) x × Grad Y J,k ) = D +1 ,J,k e +1 + D ,J,k e + D − ,J,k e − , (4.34)With c L = ( L + 1) (cid:113) L L +1 , d L = L (cid:113) L +12 L +1 , these coefficients are explicitly given by B +1 ,L,m = (cid:110) c L C L,mL − ,m − , , P m − L − (cos θ ) + d L C L,mL +1 ,m − , , P m − L +1 (cos θ ) (cid:111) e i ( m − ϕ B ,L,m = (cid:110) c L C L,mL − ,m, , P mL − (cos θ ) + d L C L,mL +1 ,m, , P mL +1 (cos θ ) (cid:111) e imϕ B − ,L,m = (cid:110) c L C L,mL − ,m +1 , , − P m +1 L − (cos θ ) + d L C L,mL +1 ,m +1 , , − P m +1 L +1 (cos θ ) (cid:111) e i ( m +1) ϕ .D +1 ,J,k = i (cid:112) λ J C J,kJ,k − , , P k − J (cos θ ) e i ( k − ϕ ,D ,J,k = i (cid:112) λ J C J,kJ,k, , P kJ (cos θ ) e ikϕ ,D − ,J,k = i (cid:112) λ J C J,kJ,k +1 , , − P k +1 J (cos θ ) e i ( k +1) ϕ . Noting (i) the complex azimuthal exponential terms e ikϕ , e imϕ in (2.18) and(4.33)-(4.34) (via the above representations for B and D ) for | k | ≤ J, | m | ≤ L ; 1 ≤ L, J ≤ N , and (ii) the need to evaluate several O ( N ) sums, of the form in (4.24)and (4.28)-(4.28), at the equally spaced O ( N ) azimuthal quadrature points (see(3.6)), we reduce the complexity by O ( N ) in each of these sums, at each adaptive-time step (described below), by using the FFT for setting up the nonlinear system(4.2), similar to the approach in [4]. In our numerical experiments (see Section 5)for adaptive-time simulation of a flow induced by random initial states, we ob-served that such an efficient FFT based implementation reduced the (non-FFTcode) computing time substantially for the case N = 100, to simulate from t = 0to t = 60.In addition, by using the fast Legendre/spherical transforms along the latitudinaldirection (obtained, for example, by modifying the NFFT algorithm in [21] forevaluation of the Legendre functions in the above terms at O ( N ) non-uniformlatitudinal quadrature points), we could reduce the complexity by O ( N ). We didnot use the fast Legendre/spherical transforms in our implementation due to thespectral convergence of the scheme and the fact that N ≤
100 in our simulations.(In these complexity counts, we ignored O (log N ) and O (log N ) terms.)5. Numerical Experiments
We demonstrate the fully discrete pseudospectral quadrature algorithm by sim-ulating (i) a known solution test case with low to high frequency modes and (ii)a benchmark example [8, page 305] in which the unknown velocity and vorticityfields are generated by a random initial state.The first test example is useful to demonstrate that the pseudospectral quad-rature algorithm reproduces any number of high frequency modes in the solution (provided V N contains all these modes), with computational error dominated onlyby the chosen accuracy for the adaptive time evolution for the ordinary differentialsystem (4.25).5.1. Example 1.
Our test case first example is (2.1) with(5.1) u | t =0 ( (cid:98) x ) = u ( (cid:98) x ) = g (0)[ W ( (cid:98) x ) − W ( (cid:98) x )] , (5.2) g ( t ) = νe − t [sin(5 t )+cos(10 t )] , W = Z , +2 (cid:60) ( Z , ) , W = Z , +2 (cid:60) ( Z , + Z , ) , where Z L,m is given by (2.17), (cid:60) ( · ) denotes the real-part function, and the externalforce f ( (cid:98) x , t ) in (2.1) is chosen so that(5.3) u ( (cid:98) x , t ) = tg ( t ) N (cid:88) L =1 (cid:34) Z L, + 2 L (cid:88) m =1 (cid:60) ( Z L,m ) (cid:35) ( (cid:98) x ) + g ( t ) W ( (cid:98) x ) + ( t − g ( t ) W ( (cid:98) x ) , is the exact tangential divergence-free velocity field, solving the NSE (2.1). Theexact test field (5.2)-(5.3) has high oscillations both in space and time, and expo-nentially decays in time. Note the dependence on a parameter N , the maximumorder of the spherical harmonics in the exact solution.In our calculation of the approximate solution u N , we chose N = N , so that allfrequencies of the exact solution can be recovered. The solution (5.3) is then usedto validate our algorithm and code by numerical adaptive time-integration of (4.1),for various values of N = N . In particular, for a fixed integration tolerance error,we demonstrate in Figure 1 that all N modes in (5.3) can indeed be recovered bythe approximate solution u N , within the chosen error tolerance, for all N = N =70 , , , Example 2.
Having established the validity of our algorithm for a simpleknown solution, we use the same code to simulate unknown velocity and vorticityfields generated by (2.1) with random initial velocity field as in [8, 12] with angularvelocity of the rotation Ω = 1 and ν − = 10 ,
000 in (2.1).The initial flow and external force in this benchmark example satisfy the mainassumption (3.1) for any s, σ > u N is spectrally accurate and converges super-algebraically withorder given by (4.8) for any s >
0. As mentioned in Section 1, this is the mainadvantage of the present paper over the recent paper [12], where such spectrallyaccurate convergence results are neither discussed nor proved. On the other hand,convergence results for two-dimensional problems on a Euclidean plane, supportedby numerical experiments, have formed a core part of research on the NSE over thelast few decades, see [8, 9, 27] and references therein.The random initial tangential divergence-free velocity field, having propertiessimilar to those considered in [8, page 305] and [12, page 988] (but not exactlysame as in [8, 12], due to randomness), is a smooth function u = v ∈ G s/ σ , aGevrey class of order s and index σ , see (2.33), for any s, σ >
0, with Fouriercoefficients (cid:98) v L,m , L = 1 , , · · · , | m | ≤ L , defined by(5.4) (cid:98) v L,m = a L exp( iφ m ) L = 1 , . . . , m = 0 , . . . , L,a L ( − m exp( − iφ m ) L = 1 , . . . , m = − L, · · · , − , , L > m = − L, . . . , L,
PSEUDOSPECTRAL METHOD FOR NSEs ON ROTATING SPHERES 21
Figure 1. (cid:107) u − u N (cid:107) for Example 1 with a fixed time-integrationerror, N = N = 70 , , , φ = 0, and φ m ∈ (0 , π ) are random numbers for m >
0, and a L = b L / || b || ,with b ∈ R having components b L = 2 / (cid:2) L + ( νL ) . (cid:3) , L = 1 , . . . ,
20. Thevorticity stream function Ψ (see (2.2)) of
Vort u ( · ,
0) in Figure 2 demonstrates therandomness of the field at time t = 0.The external force field f = f ( (cid:98) x , t ) in (2.1) for our simulation is motivated bythat considered in [8, page 305] and is exactly same as that in [12, page 988].The source term f is a decaying tangential divergence-free field which, for any s, σ >
0, belongs to C ([0 , T ]; D ( A s +1 / e σ A / )) (for any T >
0) with the only non-zero Fourier coefficient being (cid:100) f ( t ) , . The Fourier coefficient (cid:100) f ( t ) , is a continuousfunction in time and is defined by (cid:100) f ( t ) , = (cid:40) , ≤ t ≤ , cos( πt/
5) exp( − ( t − / t > . The impact of the external force on the numerical velocity and its complement(generating the approximate inertial manifold) is demonstrated in Figure 3 withthe time evolution of the velocity field matching that of the external force, leadingto little change in evolution process of the velocity as the external force gets smallerand smaller.We chose a fixed relative error tolerance to be of accuracy at least O (10 − ), foradaptive time-integration solving the N + 2 N -dimensional system ordinary differ-ential equations (4.1) in time, using backward differential formulas with variable order (one to five) and variable adaptive time step sizes that meet the fixed errortolerance.As discussed in the next subsection, the high-frequency components of the solu-tion turn out to be unreliable as t increases, [8], presumably because of the timediscretization error. We therefore retained only the frequency components up tosome order N ≤ N , correspondingly, we define an additional approximation of u N ,defined in (4.1)-(4.24), by(5.5) u N ; N ( · , t ) := Π N u N ( · , t ) , N ≤ N. As in Theorem 4.2, for the spatially discrete pseudospectral quadrature method,assuming (3.1) and exact time integration of (4.1), and hence using (2.35), (3.4),and (4.8), we get spectral convergence: (cid:107) u − u N ; N (cid:107) ≤ (cid:107) u − Π N u (cid:107) + (cid:107) Π N ( u − u N ) (cid:107) ≤ C (cid:104) N − s +1)1 + N − s (cid:105) ≤ CN − s . Our simulated approximate velocity fields U ( t ) := u ( t ) , for t = 10 , , , , , , are in Figure 4–9 and the associated vorticity stream function Ψ ( t ) of Vort U ( t )(computed using (4.32)) are in Figure 10–15. These figures demonstrate that theinitial random flow with several smaller structures evolve into regular flow withlarger structures, similar to those observed in [8, page 307]. The choice of N willbe discussed in the next subsection.5.3. Energy spectrum of the solution. If u (and hence the number of modesin u ) is unknown, as it is in the case in the benchmark test Example 2, and if u N does not contain all of the modes in u , the higher modes of u N are usuallyless accurate than the chosen practical error tolerance and hence can even violateimportant physical properties of u (because of the error time-integration tolerancebeing not very small). In such cases, it is important to choose N < N , dependingon certain known physical properties of u .For a fixed time t , the L -th mode energy spectrum of a tangential divergence-freeflow u on the sphere is defined by(5.6) E ( L ) = E ( u , L ) = (cid:88) | m |≤ L | β L,m ( t ) | , u ( (cid:98) x , t ) := ∞ (cid:88) L =1 (cid:88) | m |≤ L β L,m ( t ) Z L,m ( (cid:98) x ) . Although the analytical form of the flow in Example 2 is not known, several in-vestigations have been carried out for such fields with initial spectrum of the L -thmode decaying with order L − or L − . In particular, it is well known (see [8]), forthis benchmark test case (on periodic two dimensional geometries), that the energyspectrum of the velocity has a power-law inertial range and an exponential decay(dissipation range) for wave numbers larger than the Kraichnan’s dissipation wavenumber. Further, several random smaller structures built in the initial randomvorticity evolve into regular flow with larger structures.With u being the unique solution of the NSE (2.1), let us decompose u = (cid:101) u N + w N , where (cid:101) u N = Π N ( u ) contains all modes lower or equal N and w N := u − Π N ( u ) contains all higher modes. The existence of a relation between w N and (cid:101) u N of a form w N = Φ( (cid:101) u N ) was established in [28]. The graph of Φ is known PSEUDOSPECTRAL METHOD FOR NSEs ON ROTATING SPHERES 23 as the inertial manifold of (2.1). For computational purposes, the higher modescan be computed efficiently using an approximate inertial manifold:(5.7) (cid:101) Φ (cid:102) N ( (cid:101) u N ) = ( ν A + C ) − (cid:16) Π (cid:102) N − Π N (cid:17) [ f − B ( (cid:101) u N , (cid:101) u N )] , (cid:102) N > N . This well known approximation (without the Coriolis term) was introduced in [10,11] for nonlinear dissipative systems, including the NSE on domains and we choose (cid:102) N = 2 N .For the benchmark test case, in the dissipation range, even L E ( (cid:101) Φ N ( (cid:101) u N ) , L )decays exponentially to zero, further justifying restriction of the infinite dimensionalrange after certain values of N . As discussed in [8, page 280, 307] (and repeatedin [12, page 991]), such a faster decay (one order higher than the L − decay known isKraichnan theory of turbulence) is expected due to the Reynolds number consideredin [8, 12] being much small than 25 , , ν in(2.1), the Reynolds number is O ( ν − ) with the order constant r given by the productof the mean fluid velocity and the characteristic length-scale. With ν − = 10 , E ( (cid:101) Φ N ( (cid:101) u N ) , L ), observed in Figures 16-17. Finally, the exponentialdecay of the energy spectrum show that N = 100 , N = 3 N/ O (10 − ). Figure 2.
Initial random vorticity stream function Ψ of
Vort u at t = 0.Time evolution of U ( t ) and f ( t ). Time evolution of (cid:101) Φ ( U ( t )) and f ( t ). Figure 3.
Impact of the external force on a numerical velocityand its complement.
PSEUDOSPECTRAL METHOD FOR NSEs ON ROTATING SPHERES 25
Figure 4.
Numerical velocity U ( t ), at t = 10. Figure 5.
Numerical velocity U ( t ), at t = 20. Figure 6.
Numerical velocity U ( t ), at t = 30. Figure 7.
Numerical velocity U ( t ), at t = 40. Figure 8.
Numerical velocity U ( t ), at t = 50. Figure 9.
Numerical velocity U ( t ), at t = 60. PSEUDOSPECTRAL METHOD FOR NSEs ON ROTATING SPHERES 27
Figure 10.
Vorticity stream function Ψ of Vort U at t = 10. Figure 11.
Vorticity stream function Ψ of Vort U at t = 20. Figure 12.
Vorticity stream function Ψ of Vort U at t = 30. Figure 13.
Vorticity stream function Ψ of Vort U at t = 40. Figure 14.
Vorticity stream function Ψ of Vort U at t = 50. Figure 15.
Vorticity stream function Ψ of Vort U at t = 60. PSEUDOSPECTRAL METHOD FOR NSEs ON ROTATING SPHERES 29
Energy spectrum E ( L ) of U ( t ). Energy spectrum E ( L ) of (cid:101) Φ ( U ( t )). Figure 16.
Energy spectra of velocity U ( t ) and Φ ( U ( t )).Energy spectrum L ∗ E ( L ) of U ( t ) Energy spectrum L ∗ E ( L ) of (cid:101) Φ ( U ( t )). Figure 17.
Scaled energy spectra L ∗ E ( L ) of U ( t ) and (cid:101) Φ ( U ( t )). Appendix
In this section, we generalize certain known domain case estimates, that arefundamental for the NSE analysis, to the spherical surface case and hence prove(2.30). Following [16, Page 574], for u ∈ C ∞ ( T S ), we extend u to the sphericallayer S × I , I = ( r , r ) , < r < < r < ∞ by the formula(6.1) (cid:101) u ( (cid:98) x ) = ϕ ( | (cid:98) x | ) u ( (cid:98) x / | (cid:98) x | ) , where ϕ ( t ) ∈ C ∞ ( I ), ϕ ( t ) ≥ , t ∈ I , and ϕ (1) = 1. We have (cid:90) S × I | (cid:101) u | p = (cid:90) r r ϕ p ( r ) (cid:90) S | u | p dS In other words(6.2) (cid:107) (cid:101) u (cid:107) L p ( S × I ) = c (cid:107) u (cid:107) L p ( T S ) , where c = c ( ϕ, r , r ). Suppose u , v , w are extended from S to the spherical layer S × I by (6.1). Then from [16, Lemma 4.3] we have(6.3) b ( u , v , w ) = cb ( (cid:101) u , (cid:101) v , (cid:101) w ) , On the sphere S , we have the following version of Sobolev embedding inequality[1](6.4) (cid:107) u (cid:107) L q ( T S ) ≤ C (cid:107) u (cid:107) H s ( T S ) , s < , q = 12 − s . The following nonlinearity estimate is an adaptation of [6, Proposition 6.1] for S . Proposition 6.1.
Let s , s , s ≥ be real numbers, and we assume that s + s + s ≥ and ( s , s , s ) (cid:54) = (0 , , , (0 , , , (1 , , . Then there exists a constantdepending on s , s , s such that | b ( u , v , w ) | ≤ C (cid:107) u (cid:107) H s ( T S ) (cid:107) v (cid:107) H s ( T S ) (cid:107) w (cid:107) H s ( T S ) or in the extrapolated form, writing H s ( T S ) = H sT S , for all u , v , w ∈ C ∞ ( T S ) | b ( u , v , w ) | ≤ C (cid:107) u (cid:107) s ] − s H [ s TS (cid:107) u (cid:107) s − [ s ] H [ s TS (cid:107) v (cid:107) s ] − s H [ s TS (cid:107) v (cid:107) s − [ s ] H [ s TS (cid:107) w (cid:107) s ] − s H [ s TS (cid:107) w (cid:107) s − [ s ] H [ s TS . Proof.
Let u , v , w ∈ C ∞ ( T S ) and (cid:101) u , (cid:101) v , (cid:101) w be their corresponding extension to thespherical layer (cid:101) Ω := S × I . Let us consider first the case s i < i = 1 , , q , q , q , q so that (cid:80) i =1 1 q i = 1 and q i = − s i for i = 1 , ,
3. Then, by H¨older’s inequality | b ( (cid:101) u , (cid:101) v , (cid:101) w ) | = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:88) i,j =1 (cid:90) (cid:101) Ω (cid:101) u j ∂ (cid:101) v i ∂x j (cid:101) w i (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≤ (cid:107) (cid:101) u (cid:107) L q (Ω) (cid:107)∇ (cid:101) v (cid:107) L q (Ω) (cid:107) (cid:101) w (cid:107) L q (Ω) (cid:107) (cid:107) L q (Ω) Restricting to the sphere by (6.3) and (6.2) then using the Sobolev embeddingtheorem on the sphere (6.4) we have | b ( u , v , w ) | ≤ C (cid:107) u (cid:107) L q TS (cid:107) Grad v (cid:107) L q TS (cid:107) w (cid:107) L q TS ≤ C (cid:107) u (cid:107) H s TS (cid:107) Grad v (cid:107) H s TS (cid:107) w (cid:107) H s TS . (cid:3) Lemma 6.1.
Let δ ∈ (1 / , be given and u , v ∈ V . Then there exists C ,independent of u and v , such that (cid:107) A − δ B ( u , v ) (cid:107) ≤ C (cid:26) (cid:107) A − δ u (cid:107)(cid:107) v (cid:107)(cid:107) u (cid:107)(cid:107) A − δ v (cid:107) , u , v ∈ V. PSEUDOSPECTRAL METHOD FOR NSEs ON ROTATING SPHERES 31
Proof.
Let u , v , w ∈ V . Hence from (2.26), and by using Proposition 6.1 with s = 0, s = 2 δ − > s = 2 − δ > | b ( u , v , A − δ w ) | = | b ( u , A − δ w , v ) | ≤ C (cid:107) u (cid:107)(cid:107) A − δ w (cid:107) H δ ( T S ) (cid:107) v (cid:107) H − δ ( T S ) . Since (cid:107) A − δ w (cid:107) H δ ( T S ) = (cid:107) A δ A − δ w (cid:107) = (cid:107) w (cid:107) , we get | b ( u , v , A − δ w ) | ≤ C (cid:107) u (cid:107)(cid:107) w (cid:107)(cid:107) v (cid:107) H − δ ( T S ) , w ∈ V. Since the inequality is true for all w ∈ V , we obtain the first bound (cid:107) A − δ B ( u , v ) (cid:107) ≤ C (cid:107) u (cid:107)(cid:107) v (cid:107) H − δ ( T S ) = C (cid:107) u (cid:107)(cid:107) A − δ v (cid:107) . To obtain the second bound, we again use (2.26) and Proposition 6.1 but with s = 2 − δ , s = 2 δ − s = 0, | b ( u , v , A − δ w ) | = | b ( u , A − δ w , v ) | ≤ C (cid:107) u (cid:107) H − δ ( T S ) (cid:107) A − δ w (cid:107) H δ ( T S ) (cid:107) v (cid:107) = C (cid:107) A − δ u (cid:107)(cid:107) w (cid:107)(cid:107) v (cid:107) , w ∈ V. Since the inequality is true for all w ∈ H , we obtain (cid:107) A − δ B ( u , v ) (cid:107) ≤ C (cid:107) A − δ u (cid:107)(cid:107) v (cid:107) . (cid:3) Acknowledgments: The support of the Australian Research Council under its Dis-covery and Centre of Excellence programs is gratefully acknowledged. The authorsthank Professors M. Farge and E. S. Titi [29] for valuable discussions.
References [1] T. Aubin,
Nonlinear analysis on manifolds, Monge-Amp ` e re Equations . Springer-Verlag, 1982.[2] J.P. Boyd, Chebyshev and Fourier Spectral Methods , Dover, 2001.[3] M.E. Brachet, M. Meneguzzi, H. Politano, P.L. Sulem,
The dynamics of freely decaying two-dimensional turbulence.
J. Fluid Mech. (1988), 333-349.[4] M.A.J. Chaplain, M. Ganesh, I.G. Graham,
Spatio-temporal pattern formation on sphericalsurfaces: numerical simulation and application to solid tumour growth . J. Math. Biology (2001), 387-423.[5] C. Cao, M.A. Rammaha, E.S. Titi, The Navier-Stokes equations on the rotating 2-D sphere:Gevrey regularity and asymptotic degrees of freedom . Zeit. Ang. Math. Phys. (1999), 341–360.[6] P. Constantin, C. Foias, Navier-Stokes equations , Chicago Lectures in Mathematics, The Uni-versity of Chicago Press, 1988.[7] A. Doelman, E.S. Titi,
Regularity of solutions and the convergence of the Galerkin method inthe Ginzburg-Landau equation.
Numer. Funct. Anal. Optim. (1993), 299–321.[8] A. Debussche, T. Dubois, R. Temam, The nonlinear Galerkin method: a multiscale methodapplied to the simulation of homogeneous turbulent flows.
Theor. Comp. Fluid Dyn. (1995),279–315.[9] C. Foias, O. Manley, R. Rosa, R. Temam, Navier-Stokes Equations and Turbulence . CambridgeUniversity Press, 2001.[10] C. Foias, M.S. Jolly, I.G. Kevrekidis, G.R. Sell, E.S. Titi,
On the computation of inertialmanifolds.
Phys. Lett. A (1988), 433-436.[11] C. Foias, O. Manley, R. Temam,
Modelling of the interaction of small and large eddies intwo dimensional turbulence flows,
RAIRO Mod´el. Math. Anal. Num´er. (1988), 93–118.[12] M.J. Fengler, W. Freeden, A nonlinear Galerkin scheme involving vector and tensor sphericalharmonics for solving the incompressible Navier-Stokes equation on the sphere.
SIAM J. Sci.Comp. (2005), 967–994.[13] M. Ganesh, K. Mustapha A fully discrete H -Galerkin method with quadrature for nonlinearadvection-diffusion-reaction equations . Numer. Algorithms, (2006), 355–383. [14] E. Hairer, G. Wanner Solving Ordinary Differential Equations II: Stiff and differential-algebraic problems . Springer-Verlag, 1982.[15] J. Hesthaven, S. Gottlieb, D. Gottlieb,
Spectral Methods for Time-Dependent Problems . Cam-bridge University Press, 2007.[16] A.A. Il’in,
The Navier-Stokes and Euler equations on two dimensional closed manifolds .Math. USSR Sbornik (1991), 559–579.[17] A.A. Il’in, Partially dissipative semigroups generated by the Navier-Stokes system on two-dimensional manifolds, and their attractors.
Russian Acad. Sci. Sbornik Mathematics (1994), 47–76.[18] A.A. Il’in, Navier-Stokes equations on the rotating sphere. A simple proof of the attractordimension estimate.
Nonlinearity (1994), 31–39.[19] A.A. Il’in, A.N. Filatov, On unique solvability of the Navier-Stokes equations on the twodimensional sphere.
Soviet Math. Dokl. (1989), 9–13.[20] D.A. Jones, E.S. Titi, A remark on quasi-stationary approximate inertial manifold for theNavier-Stokes equations.
SIAM J. Math. Anal. (1994),894–914.[21] S. Kunis, D. Potts, Fast spherical Fourier algorithms.
J. Comp. Appl. Math. (2003),75–98.[22] J. Norbury, I. Roulstone (Eds.),
Large-Scale Atmosphere-Ocean Dynamics, Vol. 1& 2 , Cam-bridge University Press, 2002.[23] S. A. Orszag,
Turbulence and transition: A progress report . Proc. 5th Int. Conf. Numer.Meth. Fluids Dyn., Lecture Notes in Phys., Vol. 59, Springer, 32–51, 1977.[24] M. Pieper,
Vector hyperinterpolation on the sphere , J. Approx. Theory. (2009), 173–186.[25] L. F. Shampine, M. W. Reichlet
The Matlab ODE suite . SIAM J. Sci. Comput., (1997),1–22.[26] T. Tao, Why global regularity for Navier-Stokes is hard . http://terrytao.wordpress.com/2007/03/18/why-global-regularity-for-navier-stokes-is-hard/ .[27] R. Temam, Navier-Stokes Equations, Theory and Numerical Analysis , Amer. Math. Soc.,2001.[28] R. Temam, S. Wang,
Inertial forms of Navier-Stokes equations on the sphere . J. Funct. Anal. (1993), 215–242.[29] E. S. Titi,
Private communication , (2009).[30] D. A. Varshalovich, A. N. Moskalev, V. K. Khersonskii,
Quantum Theory of Angular Mo-mentum , World Scientific, 1988.
Department of Mathematical and Computer Sciences, Colorado School of Mines,Golden, CO 80401
E-mail address : [email protected] School of Mathematics and Statistics, University of New South Wales, Sydney, NSW2052, Australia
E-mail address : [email protected] School of Mathematics and Statistics, University of New South Wales, Sydney, NSW2052, Australia
E-mail address ::