A N A L Y S I S OF A G A L E R K I N - C H A R A C T E R I S T I C ALGORITHM FOR THE POTENTIAL V O R T I C I T Y - S T R E A M FUNCTION EQUATIONS By RODOLFO BERMEJO Dipl. Universidad M.Sc. Naval Architecture P o l i t e c n i c a de Madrid ,1977 The U n i v e r s i t y o f B r i t i s h Columbia,1986 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE STUDIES (Department o f Mathematics and I n s t . o f A p p l i e d We accept t h i s Mathematics) t h e s i s conforming to the r e q u i r e d standard THE UNIVERSITY OF BRITISH COLUMBIA January 1990 © R o d o l f o Bermejo, 1990 In presenting degree this thesis in partial fulfilment of requirements for of this thesis for scholarly department or by his or her I further agree that permission for purposes representatives. may be granted advanced It is permission. Department of The University of British Columbia Vancouver, Canada extensive by the head of understood that publication of this thesis for financial gain shall not be allowed without DE-6 (2788) an at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. copying the copying my or my written ABSTRACT In this thesis we Galerkin-Characteristic vorticity equations is stage a two material by yield stage. the a grid the of the analyze integrate algorithm. In the the efficient consists points by of first and c o n s e r v a t i v e The stage curves methods. for the this dependent interpolation of algorithm the approximated algorithm spline f o r At = 0(h). proposed Particle updating cubic characteristic equations. and a potential The m e t h o d the p o t e n t i a l v o r t i c i t y i s algorithm the of component stable inductive computationally an at foot to Galerkin-Characteristic Such variable method and a b a r o c l i n i c ocean. d e r i v a t i v e of combining This of develop the is at advective unconditionally The e r r o r analysis with converges with 2 respect to L -norm order 0(h); shows however, sufficiently smooth that the algorithm i n t h e maximum n o r m functions the foot it is of proved the that for characteristic 4 curves are The the superconvergent second Lagrangian space—time with algorithm of the is flow representation the boundary condition, These other global error and is 0(h) estimates estimates 0(h) for represent for ii the /At). a onto projection the analysis the for the an of Cartesian coordinated E l e m e n t s . The e r r o r 2 t o L —norm s h o w s t h a t 2 of to of order 0(h Finite respect condition. the representation component respect of Eularian Crank—Nicholson stage stage points with for this approximation the free-slip no—slip boundary improvement vorticity with previously reported in the literature. The evolutionary component + K of the 2 global that along error is depends the on equal the to K(At derivatives Characteristic. Since quasi-conservative quantitiy, general Numerical small. theoretical results for h), of the one the can the a constant quantity vorticity conclude experiments is advective potential both stages of iii where that K illustrate method. is is a in our TABLE OF CONTENTS Page ABSTRACT i i TABLE OF CONTENTS iv LIST OF TABLES vi v i i LIST OF FIGURES ACKNOWLEDGEMENTS viii 1 CHAPTER I. INTRODUCTION CHAPTER II. THE GOVERNING §1. Governing §2. Notation EQUATIONS AND THEIR NUMERICAL DISCRETIZATION 8 Equations 7 14 §3.1. T h e Weak S o l u t i o n F o r m u l a t i o n 18 §3.2. The F i n i t e 29 Element Approximation CHAPTER III. A GALERKIN-CHARACTERISTIC ALGORITHM FOR THE HYPERBOLIC STAGE 38 §1. Preliminaries 39 §2. Description of the Algorithm 44 2.1. The c o n t i n u o u s 2.2. The D i s c r e t e Problem 44 Problem 47 §3. Properties 58 §4. Error Analysis 60 §5. Numerical Experiments 70 iv CHAPTER IV. ERROR ANALYSIS OF THE POTENTIAL VORTICITY AND BAROTROPIC-BAROCLINIC MODE EQUATIONS 81 §1. Preliminaries 82 §2. Error Analysis of the Potential Vorticity Equation 86 §3. Error Analysis of the Barotropic and Baroclinic Modes..102 §4. Numerical Experiments 108 Appendix to Chapter IV 109 CHAPTER V. CONCLUSIONS 126 REFERENCES 130 APPENDIX 135 v LIST OF TABLES Table Page 1. Time E v o l u t i o n o f t h e Cone 73 2. Time E v o l u t i o n o f the S l o t t e d C y l i n d e r 73 3. Logarithm of r.m.s. Errors 112 4. Logarithm of r.m.s. Errors 112 5. logarithm of r.m.s. E r r o r s vi 112 LIST OF FIGURES Figure Page 1. The Two Layer Model 9 2. Regularization of the Particle Point (X^.t) that w i l l occupied the grid Point x 3a. at t + At Rotating Cone after 0 Revolutions 3b. Rotating Cone after 6 Revolutions 74 75 ..76 4a. Slotted Cylinder after 0 Revolutions 77 4b. Slotted Cylinder after 1/8 of a Revolution 78 4c. Slotted Cylinder after 1 Revolution 79 4d. Slotted Cylinder after 6 Revolutions 80 5a. Exact Stream Function Solution at Re = 1000 113 5b. Calculated Stream Function Solution at Re = 1000 114 5c. Pointwise Error distribution of the Stream—Function....115 6a. Exact Vorticity Solution at Re = 1000 116 6b. Calculated Vorticity Solution at Re = 1000 117 6c. Poitnwise error distribution of the Vorticity 118 vii ACKNOWLEDGEMENTS I would l i k e t o acknowledge, William f o r allowing me t h e freedom t o roam f a r a n d w i d e i n s e a r c h o f a n u l t i m a t e scheme; members of my s u p e r v i s o r , the following people; the Hsieh, with gratitude, my s u p e r v i s o r y Foreman a n d B r i a n Seymour of the thesis; were who crucial taught Andrew On Staniforth, me t h e b e a u t y a more U r i Ascher, Mike f o r t h e i r a d v i c e t o improve t h e t e x t i n the process framework o f f i n i t e committee, whose support o f my r e s e a r c h ; of doing and guidance and John numerical Heywood, analysis i n the elements. personal note, whose companionship during the writing of this made i t a l l possible, this thesis my t h a n k s to Martin and c o m p e t i t i v e n e s s thesis. Finally, in particular i s dedicated. viii were and Alex, very helpful t o my f a m i l y , who to Maria Jose, t o whom CHAPTER I INTRODUCTION During in the past t h e r e has been a steady t h e u s e o f n u m e r i c a l methods govern the ocean observational scientific study evidences, reasoning from models Navier-Stokes equations latter only space-time Though scales simple, a success is that formulation large of among computer a some along with the physics have very long periods one the of and a allow s m a l l enough features time able in of to to our decade. lead such to a the reasonably to minimally the m i d - l a t i t u d e to perform order of dynamics. of the use of The range the past which state (QGM). reasons assumptions, one i s of dimensional significantly during important Thus, to to hierarchy t o m i d - l a t i t u d e ocean t h e QG e q u a t i o n s , circulation. This heat related and models models dynamics the physical the mass, contributed viewpoint, of three quasi-geostrophic the ocean of technology hierarchy the time steps w i t h a g r i d spacing resolve ocean the QGM numerical interaction use significant of Active of the ocean dynamics. which retain understanding From to which produced equations down the equations numerics, has increase to solve dynamics. d i f f e r e n t aspects spans for 25 y e a r s calculations study the time e v o l u t i o n o f t h e phenomena. Most o f t h e p r e s e n t models 1 i n the ocean modeling community use standard The time finite differences, discretization leap-frog scheme the possible by one for is time carried most o f e x c e p t i o n of both the the discretization for the viscous explicitness the scheme without O(At), however, instead of The s p a t i a l e r r o r of order of main by of and the where of using this is to significantly stands scheme a with lagged time keep the altering the time e r r o r for called are the l e a p - f r o g h so which is of would i n QGM w i t h f r e e - s l i p b o u n d a r y 0(h), space. the equations, terms so in terms purpose by d o i n g (At ) as out viscous The CFL c o n d i t i o n ; time terms step. of in order suggest. conditions generic the grid is spacing interval. As of we w i l l the QG see in equations Chapter consists equation -transport-diffusion vorticity q coupled II, w i t h an the of mathematical structure semi-linear parabolic a equation- elliptic for the equation potential for the stream f u n c t i o n </». It is well known transport-diffusion the numerical equation is operator - Therefore, that equation analyst. This a combination and the numerical poses is of dissipative whether the a due global - o r p a r a b o l i c depends of against term circulation also problems a geographical to the the other. the g l o b a l 2 the - that the Laplacian of the the character. the such to an advective operator equation Furthermore, of of problem on t h e r e l a t i v e structure d i s t r i b u t i o n or fact terms structure more h y p e r b o l i c one challenging hyperbolic terms solution is magnitude in ocean equation There . are, has for instance, narrow dominant ones, significant domain, away from the efficient the the be difficulties understandable area appear under several [22], [28], [31])) marching small. why in such found with the schemes main free of o f ' these convective Reynolds terms algorithms of numbers, dynamics. The the seem is in equations to be the computational nonsymmetric l i n e a r approximate operators for of problems which i s as t h e new a l g o r i t h m s is [19], of generating dissipative The main of the flows at mechanism produced unpleasant. [21], For of the by the instance, result with a an e x c e s s i v e l y high poorly small time step the In the so called designed for dealing with such Galerkin-Characteristic Method, b a s e d on c o m b i n i n g t h e method o f c h a r a c t e r i s t i c s w i t h standard [31]). as well They the approximate i n t e g r a t i o n . One a as dynamics. for dominant is very computation which, terms spectrum, it an a [13], consequences very so, techniques. convective resolved are designing excessive the a spatial form objective the the And fluid i n t r o d u c e d by s t a n d a r d upwinding emphasis in [8], playing considerations algorithms (cf. [5], the parameterized These computational versions of the algorithm. become still areas boundaries, very terms terms i n broad solid should convective diffusive contrast, numerical research effects the transport-diffusion perfectly stable where with In terms out active but role. dissipative point areas finite element these papers, convergence and procedure the s t a b i l i t y of conservation 3 (cf. [3], [13], t h e method properties , as have [27], well been proved under evaluated In assumption this thesis, the we of the the h y p e r b o l i c it has been The n o v e l t y purely of equation operator. used our for assuming that of the inner which The in is the grid at time by t , one dimensional hyperbolic viewed first as a Consequently, first stage generated At. On by the progression space. stage particle given yields the the In this seeks in cell a datum v : the the algorithm curves second via over stage the 4 the the n this a purely can be method. hyperbolic time involves stage. t into a mapping diffusion s t a r t i n g w i t h the output w from the p r e v i o u s and element purely each level algorithm w = v o X , X being the time procedure finite m to location constructed the -> R , characteristic hand, at its respect, by of Explicitly, particle Galerkin generated an o u t p u t other of by is the been covered w i t h a approximation projected researchers. methods. Then is many any combination - approximation new b y a and n is means, c h a r a c t e r i s t i c curve a particle we not assigns level Specifically, curves by Particle algorithm to t h e way we a p p r o x i m a t e stage and algorithm characteristic years 20 is first the idea over an the the corresponding finite are approximate along At. products This form of operating the c o m p u t a t i o n a l domain has grid, the analyze equation. along algorithm Galerkin-Characteristic node and methods hyperbolic rectangular propose inductive algorithm. transport-diffusion for all transport-diffusion a two s t a g e typical of that exactly. integrate use the of Q step a time mechanism, The output v of this next time time projection into second level. the parabolic of with Euler the any the the input f o r second standard While we, of difficult nor conditions algorithm of the expensive. to be the for for discussed in of are two the Chapter as the the flow scheme for employ the reasons. to get, in First, in stage, scheme IV, the interested hyperbolic we w i s h at representation, authors Crank-Nicholson Second, of discretization contrast, scheme our thought Eulerian of in be stage representation time most the f i r s t may space-time Crank-Nicholson context stage Lagrangian scheme, implementation in The equations. backward the is Cartesian coordinated using stage the is neither under certain a higher accuracy time. The solved spatial is d i s c r e t i z a t i o n of made by functions on a g r i d of this using purely C°-finite composed kind hyperbolic of constructed. We show the upstream l o c a t i o n of spline elements, polynomials clear how integrals as for on the With example, can of the get elements as the a linear bilinear The m a i n is is on 1 be basis at the solution is solution at to cubic of C°-finite triangles, quadratic elements, relationship products and evaluation of algorithms. the G a l e r k i n - C h a r a c t e r i s t i c method 5 apparent classes to advantage equivalent interpolation classical The with particle C -finite simple inner equations projected p a r t i c l e other or of rectangles. the p a r t i c l e s rectangles, one of stage that interpolation. elements finite first the set it is not between the simple such linear integrals proposed in in [13] and [31] is eventually have a very lead demonstrated. approximation t h o u g h no show in [27] ideas exists, Chapter compute III. related the still In curves meteorological literature use of The briefly thesis first III performance Chapter 2 L -norm and of IV of integration as follows. In the on of the in for numerical a the formulation. analysis of order illustrate to estimates the with the numerical the equations. respect algorithm vorticity we two-layer in potential the II Some of the heuristic algorithm. stage that of the the e r r o r second to the f o o t the a l g o r i t h m w i t h l i n e a r h y p e r b o l i c we g i v e used Chapter the f o r m u l a t i o n and presented in [34]. their of that theoretically construct to we mention based substantiate as classical utilized QG e q u a t i o n s stage are previously but integrals are also at our interpolation the should [27] apply mention of continuum devoted hyperbolic experiments In the ocean is organized can may al algorithm, we i n t e r p o l a t i o n made i n is introduce stable [34]), we et linear the v a r i a b l e s (cf. Morton methodology We been thesis cubic spline baroclinic Chapter this has which we connection, method. i n t e r p o l a t i o n of In and a [19]), inner product products characteristic arguments. get particle inner Galerkin-Characteristic cubic spline the them this with as Nevertheless, between and (cf. algorithm to evaluate relationship algorithms process t o an u n s t a b l e recently particle costly to in equations. the the Also, 2 L -norm error estimates and v e l o c i t i e s a r e of the studied. 6 approximate stream functions Each chapter the s e c t i o n s of divided into sections and, when needed, The subject are f u r t h e r divided into subsections. each chapter technical is is i n t r o d u c e d by a s h o r t details, of the main description, aspects and free results of of the chapter. The criterion chapter consists f o l l o w e d by denotes refer to of have are also proof appendix, writing equation first attached the an number first of been a Chapter of Lemma placed at IV 2.3 is of the If chapter V. Two we deals of t h e end o f with this the theorem. 7 The the do it Appendices the we h a v e throughout the by Conclusions first In to number one technicalities chapter. thesis, Thus, we h a v e technical completeness. each section such section. closed w i t h the i n e q u a l i t i e s which a r e employed the Lax-Milgram number in the c h a p t e r f o l l o w e d by the for and equations s e c t i o n 3. different Chapter included the 12 o f The t h e s i s in the the e q u a t i o n i n t h e number o f presented to to e q u a t i o n number the equation. which the of t h e number o f (3.12) writing adopted is of second collected thesis and CHAPTER II THE GOVERNING EQUATIONS AND T H E I R NUMERICAL D I S C R E T I Z A T I O N This descriptive equations them. and In the proposed addition, throughout the chapter we thesis, introduces numerical also present and the the algorithm the spaces basic where continuum to integrate notation the used solutions are sought. 1. Governing We w i s h the Equations to apply the proposed n u m e r i c a l a l g o r i t h m c i r c u l a t i o n of to f a c i l i t a t e example two of layers a mid-latitude baroclinic the mathematical a n a l y s i s such of an ocean. Thus, densities mid-latitude, our and p 6-plane, ocean. we c h o o s e ocean rectangular, study In order the model respectively 2 > to simplest consists confined flat bottom with solid of to domain a D A rotating Figure with 1 angular illustrates coordinate system, the The model. each layer subindices upper and velocity stream are 1 and lower as the well as the and for whereas 8 the of geometrical and by stand layers, arrangement functions denoted 2 Cj a n d relative u , Cartesian parameters of vorticities of respectively. magnitudes 1- s t a n d s our boundaries. for related the to The the magnitudes Fig. 1 The Two L a y e r Model 9 evaluated at The acting flow on through is driven free in the v e r s i o n by variable The of wind coupling the stress between interface. We [20]. The z-dependence central Thus, defining finite the of T(t,x) layers use potential vorticity-stream employing direction. each l a y e r a surface. c a n be f o u n d i n vertical by displacement modelled we interface. is the the equations which the the the QG function variables differences potential is in the vorticity at by have i 2 Dq -fif- - A V'q. H where Q D Dt f = T df/dy - 5.^ 2 qAx,t) = qXx.O) = f = Qx[0, = 0, on t h e i n <? r , i boundary =1,2, (1.2.1) T, - , (1.2.2) (1.2.3) T], 8 dt 2\Q\sin9Q approximated its + S. cq. tangent and = f + (3y is the Coriolis by p r o j e c t i n g the e a r t h ' s plane at the reference are parameters parameter, spherical latitude 6 = whose m a g n i t u d e s 10 which geometry 6^, and depend on is onto /3 = | f i | , Qg and t h e r a d i u s g'= of ( the e a r t h )g, p here g (cf. [30]). denotes the acceleration due to 2 gravity, _ . the f o r c i n g F curlT = H A„ i s ti is a over two Define layer *(x,t) ocean, = $ H V(x, barotropic we may = t) and modes, as - ifi (x,t) and \jt lead 2*2 „ + , 2 (1.3.1) H $(x,t) , are (1.3.2) the respectively. decompose so called According the baroclinic to McWilliams baroclinic mode * in and [26], the fashion a considerations the g l o b a l of follows. ¥ ( x , t ) = * (x, t) + C(t)V By combinations e q u a t i o n s w h i c h a r e now f o r m u l a t e d . (x,t) l*l ' further following linear and $ ( x , t ) t) *(x,t) assumed symbol. of e l l i p t i c *(x, c o e f f i c i e n t which i s l i n e a r bottom f r i c t i o n c o e f f i c i e n t , the Kronecker to a set where eddy v i s c o s i t y t h e w h o l e d o m a i n Q, a constant 8. . i s For ' the l a t e r a l constant c l of s (x). mass c o n s e r v a t i o n (1.4) [26], ^(x,t) satisfies condition J f fi Kx, t)dQ = 0 , V t € [0,T] 11 . (1.5) In view o f (1.4) and (1.5) the time dependent determined a t any i n s t a n t (x, of (1.1), CV 2 , V t e [0,7]. 2 (1.6) (x)dn ( 1 . 3 ) a n d ( 1 . 4 ) we o b t a i n - A ; * (x) s = 0 i n Q_ , 1 2 * CV t)da = - — f * By v i r t u e is t by * C(t) function C(t) s (1.7.1) |_ = 1 , V t . r - A ;* Cx,t) = q 2 a * a |_ r 2 (1.7.2) - q in Q 2 T , = 0 , Vt , (1.8.1) (1.8.2) where f r// +ff; 2 V * = * 3 ffq = V #| Vt . The s y s t e m s ( 1 . 7 ) - ( 1 . 9 ) * or a (x,0) ( i - 8 3 ) + ff g vVx,0 = 0 • 2 A ii - / i n Q_ , (1.9.1) (1.9.2) are closed with the i n i t i a l = $(x,0; = 0 in Q , equivalently 12 conditions (1.10.1) ip^x.O) Finally, = \p (x,0) the velocity d e t e r m i n e d a t any u 1.2. The = curty i stream along On , and the boundary other T; at that hand, , the is for As vorticity formulation well vorticity physical however, the of n the boundary, t instant meaning transport in the boundaries, r state that t, the constant ( 1 . 12) the N-S is an solid denoting the the is that across a the solid value of the boundaries of the of the boundary assuming condition that , normal the function-relative ( 1 . 13) and tangent r e s p e c t i v e l y the f o l l o w i n g r e l a t i o n s 13 boundary from the values additional Specifically, = g.(x,t) zero equations, the of mode stream one can e a s i l y deduce, the and any barotropic unknown"on the stream f u n c t i o n . with (1.8.2) i = 1,2. known is u.l and is is = C (t) r boundary. for layer (1.11) ip^ a r e , i n t e g r a t e d normal on each . and vertically velocity at relation (1.7.2) (1.9.2) domain; (u^,v.) i = 1,2 condition relative u^= Conditions. functions the (1.10.2) t by the conditions <p.\ in vector instant On B o u n d a r y boundary = 0 2 vectors hold on 30. — |_ = -g.ot Sn T , "i (1.14.1) dip at If in the * other terms, In on the free-slip if wall, so boundary * stress on the s o l i d w a l l , relative vorticity is and 0 is = 0, known conditions for then there exists on the Now, the case g^ boundary is on generated. = 0 is = as denoted the value i A , or shear literature assuming that the p o t e n t i a l v o r t i c i t y 0 - -j— C(t) g n . is i n f l u i d dynamics condition. vorticity i d e n t i c a l l y zero relative vorticity condition while boundary relative i s no s h e a r contrast, boundary by n o - s l i p the the condition (1. 14.2) 6 then there wall. stress The 0, — |_ = g.on = 0 . T i 1,2 of the are f q.. ul 2. = X + f 1 , V t e 10,TL (1.15) Notation We now throughout introduce the some thesis. Let basic fi be notation an open which subset is of R used 2 with 2 Lipschitz space. Let which have boundary C^Cfi; T, be continuous w h i c h a r e bounded R is the set the of two real dimensional valued partial derivatives o n fi. On C^(Q) 14 of we i n t r o d u c e Euclidean functions order at the norm on least fi k, \v\. = ' k Boldface symbols sequences; sup xeQ \a\<k ' ° M \D v\ denote either the meaning w i l l multi-index notation. (2.1) . a If matrices , vectors be c l e a r f r o m t h e c o n t e x t , IN denotes the set of a or is a nonnegative integers a We h a v e := ( a , , a „ , . . . . a ) , a . e IN, 1 2 n l 1 = 1,2 . . ,n the f o l l o w i n g d e f i n i t i o n s |a| = a . + a + . . . + a , 1 2 n a £ (3 i f f (a.+ ( a = D a a > Vi i = a i~ i 0 + i ~ 1 8 i (a'.)(a>) 1 2 Ca - &) x i |3J a/ = - P a d o n lP(Q), i = 1,2, . . ,n = 3x , iiuii !) , ; 0) , Vi = 1,2, . ..,n , an n 1 s p < co , the space measurable of all equivalence functions. The norm d e f i n e d by = p ' 2 r e a l - v a l u e d Lebesgue p < co i s , ' ' 3 n a ai We d e n o t e b y LP(Q) of V , , an. (x ) , n 3x, 2 classes ' - p. , a i , , a.2.. = (x, ;(x_ ; 1 2 « , Vi ...(a = max(a. i = 1,2, . . ,n n ( { l"l <*n J p 15 P , (2.2) with Hull When p = = ess 2 we d r o p sup\u\, Vx € Q . the subindices p and Q if there is no confusion. For e a c h i n t e g e r m 2: 0 a n d r e a l Sobolev spaces W ' (Q) = m P (Q), p, 1 < p < 00 , we d e f i n e t h e as v e L I T Q ; : D v € L (Q), j P l a l s m . l s p s . o J . P a (2.3) T h e s p a c e M '(Q) i s a Banach s p a c e w i t h t h e norm m,P Ml n m ' P ' = ( Q £ IID°Vll \a\sm _ ) P ' , 1 s p < 00 , 1 / p (2.4.1) Q and t h e seminorm | V » O = < L I l^ln ' O (2.4.2) When p =00 llvll Also, _ = the space represents max (ess I a I sin W^^Cfl.) sup\D%\), is the completion of separable t h e space Vx e Q . for 1 (fen) (2.4.3) < in p < 00 a n d t h e norm II II m,p,Q Let X be a Banach space with norm r e p r e s e n t s a f u n c t i o n d e f i n e d o n Q^, we s e t 16 II II. If v(x,t) L (0,T;X) = P | v: | \\vl\ dt < oo | P (2.5.1) , 1/p = (T \\v\\ dt VO Ivll P L (0,T;X) P (2.5.2) T 1^(0, T;X) = | v: J di IID*vll' m < oo, a ^ m _ -\ 2/2 llvll {/"(O.TjX) a h e r e D.denotes ' a 3 t (2.5.4) a=0 a a ess supWD v\\ < oo , a ^ m j - , f 0, 17 Hull = max ess supWD v\\ , m Ottm [0,1] ^'"[O.TjX] In (2.5.3) 0 (2.6.1) (2.6.2) particular L (0,T;X) P = H°(0,T;X), L (0,T;X) = W°' °°(0, T;X). m S o m e t i m e s we u s e t h e s h o r t h a n d n o t a t i o n , i f c o n f u s i o n does not arise, L (0,T;X) P For p = 2, s L (X), P IT^'^n; iT(0,T;X) = ^(X). = [/"(n), w h e r e lF(Q) o f o r d e r m. The s u b s p a c e H^(Q) 17 of lf(a) is the H i l b e r t i s d e f i n e d by space = 0, 0 £ 3. T h e Weak S o l u t i o n We formulate problems with, (1.2), in this section (1.7)-(1.10) we d e f i n e t h e c l a s s e s H^(£2); subsets however, of solutions it H*(n) (2.7) m-2 Formulation f o r m u l a t i o n o f t h e weak is I<x| ^ t h e weak and t h e v e l o c i t i e s u^. convenient which are we a r e s e e k i n g . of To of f u n c t i o n s which a r e used form of the problems. is solution specifically begin i n the Our b a s i c to introduce some related the space closed to the Such subsets a r e : (3. 1.1) (3. 1.2) Thus, the stream f u n c t i o n s i/». € S , w h e r e a s t h e p o t e n t i a l a n d l c . v o r t i c i t i e s b e l o n g t o S^. T h e s t e a d y c o m p o n e n t * (x) relative of the b a r o c l i n i c other hand, the mode time V(x,t) is a n d t h e b a r o t r o p i c mode $(x,t) In (1.2), order we developed the to motivate introduce i n Chap. l i n e a r -advection w i t h C = 1. On t h e c b a r o c l i n i c c o m p o n e n t * (x, t) 3. dependent III. a r e i n HQ(Q) f o r a n y t € the here in S algorithm some Recall equation 18 ideas to solve which equation are the c h a r a c t e r i s t i c [0,T]. further curves of at (3.2.1) q(x,0) satisfy = q (x), Q t h e system dX(x,S;T) = U(X(X,S;T),X) dx Xix, , (3.2.2) =x , s;s) w h e r e u i s t h e v e l o c i t y v e c t o r w i t h d i v u = 0. The transformation conditions to be transformation equal (cf. of t o 1 almost To solve -» into defines, X(x,s;t) specified fi below, itself with under a certain quasi-isometric Jacobian determinant everywhere. the potential vorticity equation ( 1 . 2 ) we s e t [31]) ^ = ^ (x,t:T),r)\ (X where denotes The hand s i d e q x right along - the consideration the total derivative. of (3.3) represents characteristic suggests (3.3) , T=t the curves following algorithm. Suppose that a p a r t i t i o n o f 19 [0,T] the time e v o l u t i o n of of the two-stage flow. This inductive ? •• = \ 0 = t is <t Q s p e c i f i e d . We s h a l l by T upon is define Q -» R , i m =T1 the semi-discretization induced 0 s n £ m, we a s s u m e = 1,2, 2 that a ' r e f l e c t i o n ' map E : H (Q) n L (Q) S such < t i n t h e f o l l o w i n g two s t e p s : (1.2) 1) G i v e n q\: there < t that > H (Q' ) r\ L (n' ) , a S functions C^(Q) a a r e mapped onto functions, C^(Q') and £ i s l i n e a r and continuous, a n d Q' i s d e f i n e d b y n+1 n 1 ' * Thus, 2 a n y f u n c t i o n q : Q' q. Consequently, > R c a n be e x p r e s s e d a s =Eq. 3.4 t h e o u t p u t from s t e p one i s then given by 5? = S ^ V * ' W V ' V . The method this thesis, chapter. to approximate f u n c t i o n q^.(X,.(x,t ^hi q } i s one o f t h e c o n t r i b u t i o n s 1 and i t w i l l F o r now, i t hi approximate values, space where be d e s c r i b e d is sufficient , ; t ), t n+1 does (3.5) n n ), the solutions 20 to where not belong in detail i n the next indicate that the subindex to a finite a r e sought; of the h denotes dimensional however, the f u n c t i o n q\ . m u s t belong make t h e f u n c t i o n a }, be i n 1 hi from other notation t o H,. It i s the fashion t h a t makes o u r m e t h o d d i f f e r e n t h Galerkin-Characteristic i s now i n o r d e r : methods. is 1 X.(x, l t 2) computed t n+1 from 1 such n+1 l q (x) 1 i2~^~ A q ( i + i 1 t To (3.5), that q - 2 = A n solve g" n+1 ^i + J bi , e 9 / ) + f = - ) + 5 i - 4 - f v<rgf / f n - r ^ - c " i the boundary + ' Z / 2 ' Q) in T values i points 2 such n e T > , = HQ(Q) ® M . that Vu e W 21 ve ; + Ztf)), 1 V i G = H oC ^ n (3.6.1) (3.6.2) 1,2 ,; t n+1 e ). n A." o f t h e r e l a t i v e v o r t i c i t i e s t h e method o f G l o w i n s k i ff rn; and Vr 0 and c o - w o r k e r s L e t M be a complementary subspace + the the following: e HICQ) 7 Find ^ " ^ e at that u Q t/l we a d o p t [16]). takes 1 = X_uvdQ , and X ? = X.(x,t (u,v) update g? on A 1 " V * ^ ' + g"t where i ^1 values - g.CX") n+1 d the Given the output from i word ). n Find q , ( - A We s e t g ? = g.CX ?,) t o r e m i n d u s ^1 g? i n w h i c h we o f H^(n); (cf. [11], i.e., (3.7.n lfj - a " 1* , n) = ( V 01" , 1 Next, V +1 n+1 ; - J M u ds . (3.7.2) we calculate the barotropic and baroclinic modes, the stream functions and the velocities. i) Find 9 (x) e S ,c = 1, such that V<£ e HI(Q) s c 0 CV* , 70 ; + A ( * . <f>) = 0. (3.8.1) 2 s ii) Find <b (x), $ n+1 a n + 1 (x) e n\(Q) such that V0 e fl^ffi; C/ c /n.n+l „. . . 2, n+,l ,. , n+1 n+1 , (7* , 70; + X (9 , <(>) = -(q - q_ ; , 3 3. 1 c. ,_ „ (3.8.2) T (V9 , V<f>) = -(b ,4>) n+1 ,n+l b where , n+1 H = „ n+1 ., n+1 l l 22 q + H H +H 2 (3.8.3) , - f q 2 iii) f ^ From (1.4), * (x)dn n+1 = - - ^ 1 * (x)dQ (3.8.1), (3.8.2) . and (3.9) (3.9) we are able to determine the baroclinic mode ¥ * ( x ) by n+ V n + 1 (x) = <H (x) a n+1 + c"" " * (x) . s 1 (3.10) 1 Once the baroclinic and barotropic modes have been calculated, at instant t ,, n+1 the stream 22 functions 1 and 2 are obtained i n view o f (1.3.1) .n+1 * 2 ,n+l " 0 by s o l v i n g n+1 = * 2 n+1 and (1.3.2) n+1 n+1 H+H Finally, u (3.11) 2 n+1Cx) i s obtained from n+1 , n+1 n+1, i Equation Remark. beginning Our (3.12) (3.8.1) once next admit shall assume objective ourselves of t h e systems i n H*(Q), f o r any t that t h e boundary values are sufficiently smooth. prove the existence ( 3 . 6 ) f o r t h e upper proof f o r the lower main tool used in layer spaces V ' (Q) m P (T), Towards is T. We restrict and uniqueness of the layer Also, (i = 1), since the t h e same. T h e the Lax-Milgram t h i s end l e t us f i r s t P —m W theorem introduce and t h e trace W^' (Q), where T i s t h e boundary space e and we o f t h e spaces 1 £ p < co, t h e d u a l (3.6) and t h e f o r c i n g ( i = 2) i s e s s e n t i a l l y the proofs ( c f . [2] and A p p e n d i x ) . spaces that solution to solutions i s t o show a uniqe function For and f o r a l l a t t h e o f t h e computations. (3.8) the d u a l i s solved o f fi. D' is (Q) t h e space of 0 continuous l i n e a r f u n c t i o n a l s d e f i n e d o n w"l' (Q) w i t h t h e n o r m P 23 IfII where p ' It is , = 0 satisfies usual to p if 1 J 1 / 2 -1 -1 + p' [ < ' " „ „ Hull _ in, p, Q f > = i , and space (3.13.1: , u * 0 , l r <f, u> = <f,u> a s d u a l i t y denote belongs to the dual The sup ,m,p uelt^j -m, p ,ii fudfi . pairing, where f o f u. t r a c e spaces a r e d e f i n e d by = j <7 e L (Q) (D ; 3 u € ^(Q) 2 fa yyx = g on V J , f o gj where y . = dn J — i s the trace operator of j - t h order which can J be e x t e n d e d to continuous (D. linear operators mapping I n t h e f o l l o w i n g we u s e t h e s p a c e H ( Q) m (D H onto with norm d e f i n e d b y '•7J ,,o 1 o r b m-1/2, 2, T = i n f im n • l l u l 1 (3.13.2) uelf (Q) >2>& 1 Remark. F o r f u r t h e r p r o p e r t i e s o n t h e t r a c e o p e r a t o r s a n d t h e s p a c e s H" ^ ^^ (T) 1 Proposition see Appendix 2 For q ^. € 3.1. 1 1 of Chapter H 1 / 2 1/2 Proof. Since exists a q^ e H (Q) s u c h Next, H (D X consider is the range that the b i l i n e a r form 24 and ( D problem (3.6) has a unique solution IV. F X I + 1 / Z € H~ltl), in H*(n). of 1 there o n T, f o r i = 1, 2. i n H (Q), t h e n 1. = c j ^ u, v e HQ(&) > Ku.v) = (u,v) + v(Vu.Vv), LtA where It is well known coercive; for H v = that (cf. is, [29]) that there exist \I(u,v)\ constants < C Hull \I(u,u)\ Also, is I(u,v) the problem subindex i) n+1 that £ C n+1 continuous we h a v e to functional; n+1/2 to examine for i if = 1 - q (x) nfl Q to , e) + (F , n+1/2 u,v e H*(Q). (we drop the e H^(Q) a n d e) . ve e H^Q) apply, the Lax-Milgram theorem (3.14) defines i.e., > (q(x"), (F a n d H*- c o e r c i v e Q the r i g h t hand belongs . 2 1 0 8 0 H e H (SI) s u c h t h a t q (x) - Kq , order Hull 0 q , B) = (qCx"), B) - vCVqCx"), V8) - I(q - so > 0 such HQ- is: F i n d q (x) that , and Hull H In continuous v i n Hg(Q) u, Now, is I(u,v) side the of it . (3.14) remains to a continuous show linear mapping 8) - v((Vq(X ), n VO) - I(q , 8) + ( Q \ ) T1+1/2 F Q H *(&). We f i r s t n o t e t h a t Kqg> 0) i s c o n t i n u o u s , 8) - I(q , Q B) e H~ (Q). 1 inequality 25 Next, . by the Schwarz vCVqOf"), VQ) £ C\\Vq(X )\\ IIV6II , n t where C i s the a constant and X quasi-isometric Theorem 1.1.7 of is in , which homeomorphism [25] x —> is the By X. n image of Q virtue of and t h e f a c t \\q(X )\\ = \\q \\ , n since the Jacobian homeomorphism HVqll i s is the r i g h t With satisfied, the follows • P r o p o s i t i o n 3.2. in in HQ(Q), we have that (3.14) is and , i n H ^(0.). of the uniqueness Lax-Milgram of q ^(x) n+ theorem H^(Q) e (3.8.2) and (3.8.3) have a unique whereas the unique solution of (3.8.1) lies ^(Q). Proof. before, the 2 assumptions Problems quasi-isometric everywhere, v e ; e H~ (Q) 1 existence the Hence e) - vdqof ), the of 1 almost n hand s i d e of all solution to t o IIVq ll. (qCx"), so determinant equal equivalent n We f i r s t p r e s e n t let Problem e H*(Q) (3.8.1) Find * s the proof such for that i n the f o l l o w i n g - * 0 0 € HI(Q) 0 26 the Problem 1 o n T. (3.8.1). Next, fashion: s u c h t h a t V</> e H\(Q) 0 we As recast w* - v ), v^; + Q A C* - * , 2 S 0) = -cv*, 0 Vc5; - A 0) ^ , 2 (3.15) Define now t h e b i l i n e a r f o r m I'(u,v) = CVu, Vv) + \ (u, I'(u,v) I' (u, (3. 15) b e l o n g s guarantees Next, a n d H*- c o e r c i v e , t o H *(£l). Therefore, we p r o c e e d of essentially t o prove problems t h e same i n H*(Q). the existence (3.8.2) hand side of theorem of the solution V (x). and uniqueness and (3.8.3). f o r both problems. Secondly, so the right t h eLax-Milgram t h e e x i s t e n c e and uniqueness solutions are v). 2 i s continuous v) as of The p r o o f First, note t h e b i l i n e a r f o r m I"(u,v) the is that cr? i + defined by = CVu, Vv) + u C u , v) , u, v € H (a), I"(u,v) Q 2 where p. = (3.8.3), A is assumptions f o r Problem continuous (3.8.2) and u a n d HQ(Q)- o f the Lax-Milgram = 0 coercive. f o r Problem Since theorem a r e s a t i s f i e d , a l l the then i t 1. follows * t h e e x i s t e n c e and uniqueness i n HQ(Q). 3.2 and (3.11) and u n i q u e n e s s o f t h e stream s u b s e t s S„ a n d S „ . o f H (Q), l H 1+ H 2 ' C 2 the existence n 1 = ~0-nr 2 27 one deduces a °2 are defined by = and f u n c t i o n s "A"*^ d ip^^in t h e r e s p e c t i v e l y , w h e r e C , a n d C_ 1 °1 ¥ • From P r o p o s i t i o n C of the solutions a n d ^ ^ given 2 by (3.9): Remark On Regularity With reference Lax-Milgram II Hq of the to the Solutions. regularity theorem s t a t e s of 1,2,0 11 i < c ii » ll* |l a i,2,n ~ 4 n w h e r e Cy that g n 1,2,0. n 2-1,2,0 ' s C lib" -1,2,0 6 C^, and are constants w h i c h d e p e n d o n Q. 2 However, one c a n assume s c a l e wind condition qAx.O) regularity relative by i s a smooth the regularity In f a c t , vanishes at suppose theory that . for elliptic tha and t h e i n i t i a l polygon to expect i f one assumes t h e boundary, g because Q i s a bounded i t i s reasonable of the solutions. we m i g h t function, Since m corners, 2 F e L (0,T;L (0)), = f e C (0). vorticity t that stress no r e e n t r a n t instant the : " II* 1i s 1,2,0 ll* |I with solutions n n large the then = f + c" € H problems higher that the a t any (T) a n d (cf. [17]) we have " ^ " 2 , 2 ^ ^ 7 - ! where the constant " " F Ai"3/2,2,0 + depends \> (3 - 16 o n 0. When t h e r e l a t i v e v o r t i c i t y d o e s n o t v a n i s h a t the boundary is 3/2 not c l e a r w h e t h e r q^ e H f r o m o u r f o r m u l a t i o n t h e most 2*2 2 /2 i s t h a t g^e H (T); t h e r e f o r e , (D; we c a n s a y , i n v i e w o f ( 3 . 7 ) , 28 n by v i r t u e o f t h e t r a c e theorems q Next, we examine sufficiently that from ll+2 e H ^), the then w i t h no r e e n t r a n t there our problem corners, number is a If k equations fi £ 0 "•"'W^n x 3 two d i m e n s i o n a l polygon ( 3 . 16. 3 ) k 9 c i s g i v e n by ( k > * r*2 k.2,n' a n n n X.2,n a ) n b • t h e n ¥ , $ . e H^(&); n 1 imposed by t h e geometry n however, o f t h e d o m a i n fi i m p l y the that g H (Q). n e (cf.[17]) 8< ' q } a n d q " a r e i n H^(Q) $ such s n c is ^ ( ) , $ , w i t h Q s u f f i c i e n t l y smooth, n "^W.n* constraints a real modes. 2 s Since the e H (Q) \(i (x) The r e g u l a r i t y o f $ so of of e l l i p t i c fi that (3. 16.2) is theory we h a v e . n regularity the r e g u l a r i t y in Vr 1 n l smooth, but fT (Q); (cf. [29]), (3.16.4) 2 a 3.2 The F i n i t e Element We now proceed approximation first Approximation to of equations introduce the parameter, partition fi vertices into = C*^' x 2j*' the finite ( 3 . 6 ) - ( 3 . 12). Towards finite discretization of describe element 0 < h small (J> 29 £ spaces. < 1, rectangular e 1 2 > 1 - Let and this h - 1, e n d , we denote let elements i element 1 - a be fi with g J a - J. and s i d e s p a r a l l e l t o the c o o r d i n a t e axes. (i, i j), and directions, are local respectively. in the h = diaCQ^), e j sense that each v e r t e x k = in the identifiers The p a r t i t i o n 25^ = U there exists p = sup For a constant c dia(Sph and fi is e x^ regular <r > 0 s u c h t h a t for K the Sl^) h p. Given positive e integers numbers o f vertices t h e number of K and (or g r i d elements Q h = { h H In v H equal The is (i 1 •• (Q) n "l™- oh = h s ° C e (3.17.1), degree points) s v h\ Xg = constant, the interfaces its so that is a This fact of with <t>Xx ); *k V (x X 2 } v n e = x 1 < to is equal g dimensional x ) 2 to spaces h of basis are h) > = h s c n» • is o. h polynomials in x ^, The d i m e n s i o n the lines set ( linear w r i t i n g <p,(x) ^.1: x^ of \ P^\ y of r e s t r i c t i o n of piecewise by 2^i x ' X = constant 4>^(x) a l o n g 1 ~ k ~ K 1 ~ of is and any of polynomial in as product a tensor one i.e., t ( > t>/ ' > i the permits 2 < equal D € H space canonical variable. x 1 r n> \ s i n t e r f a c e between elements 1 p e = 2 i n our computations). t o K and <t> ( ) and N are defined ch = c the in is h w i t h the p a r t i t i o n H where e > i n D, , t h e f i n i t e e associated N 1 ~ J ' 1 ~ J ~ 30 I J <t>.(x. ; 1 lr where the = 5. , lr 4>^( ^ a x d n finite *j (3.7.1), piecewise [16] we analogue : u. | _ wi 1 H. h dimension points. of h The j * s * J linear the polynomials in according to Next, define the space space W^, the introduced H of = 0 Vfi e D, , (2 n T = 0 fi e h e " in is ( 3 . 17. 3) equal to the number is defined boundary by = i1' of the l f P i 6 r [ _ 0, o t h e r w i s e v consists of h = \ w , 1 s i s K., w. . € M , v. .(P.) I hi b hi h hi l . of (3.17.2) n The c a n o n i c a l b a s i s B, o f M, support }• , ' = H, ® M, . Oh h h B i s r ± I , by M. = •{ u, e H, h h h The , js respectively. x^, Pirbnneau dimensional = 8. are and and ) 2s (pXx^) variables Glowinski <f> .(x. elements adjacent 1. I ' to boundary. The spaces ^Qh^ a r e m e m D e r the s f a m ily of" finite 1 00 dimensional spaces (Q)(HQ(CI)). properties approximate Al) If H.c h Q^c V ' We now f o r m u l a t e of ^(Q), inf „ * "h (cf. [2], 1 + I>ntz0, llu - all i n the e r r o r complete and analysis in inverse of the [12]). t h e n Vu <= H (Q), £ Ch \\u\\ p s, 2, is some a p p r o x i m a t i o n which are useful solution which (Q) fi € 31 r*0, r " , r , 2, fi 0<s^ min(r.m) the min(1+1-s, r-s) . where p ^ 1 00 A2) If i n a d d i t i o n , we a s s u m e t h a t u e W ' (Q.) t h e n llu - Z « inf „ - Chllull, * "h e A3) For 2 s r * 1+1 inf I llu - *ll + hllu - ^11 + hf Let k^, llu - a» such n be a r b i t r a r y that W , ,. II) J + hllu oo, f2 then + _ ; I £ Ch llull 2, oo, fi numbers there exists r r,_ 2, a n d p, q € [l,a>] s u c h a constant that //^ c C = C(<r, k^.k^, p, q) e H, , h II IIv. II. _ h k2, ^ „,N/q-N/p+kl-k2.. _ ^ Ch q, .. IIv. It, , h kl,p, fi fi w i t h W t h e d i m e n s i o n o f t h e d o m a i n fi. 12) IIv II. _ * C h l l v , II , n l,co,i! n 2 , 2 , fi _ 1 ., „ ^ 2-JV/2. . n oo, Let in H, h compute Q .(x), * n l to the approximations ,. n 2, 2 , fi , Cx,), ¥ . (x.) a n d <J>fYx,) d e n o t e n sh cj.Tx;, i -1A-1/N,. fi * ah finite as (x), s h Y (x) a element linear and solutions combinations o f H.. T h u s h 32 $ Cx^, of the approximations respectively. we expand the basis To such functions ^ = pik*k k *sn =l*sk*k k (x) * an n (x) Likewise, in (3 * 0h * 0n H ( 3 . 13) the ^" V * (x) = E k *2h =l* 2n*k k ih (3 • stream (x) n XJ) h functions ' l e S„, , Ch , a r e known, sh can 3) 18 4) also ' (x) = f, i =1,2 , (3 (3.6)-(3.8) is defined ah = Q° h = 0 , with C = 1, s u c h t h a t T h e n f o r n > 0, a s s u m i n g * 18 be (3.18.5) element approximation of G i v e n Q° hs 2) - - (3 follows: f i n d #, 18 ' " H (x) 1 8 1 ) ' - S l*k+k k of - as n The f i n i t e = ( 3 * Cn (x) (x) view approximated \n € =l*ak*k k (x) *h ( x ) solve (?", for 33 , , i//?^ a n d ' 18 as 6) 1) < + At i 2 - f < 8 Q (x)l 3 i ^ ' + Q* b n+1 . . n) ^ J > +1 6 (W ( Q =* T ^ n U ^ 1 2 > ^ * V ^ V ' % (3.20) n uh h . that , ^ ah + "oh . e tf , . , V<0 ) + A T * ah ^ V -,n+l jn+1 „ , * , , $> e such ah n On x + = - ( Q - Q l 2 , <t> ), V<p e H^^ h h (3.21.1) ' V CT V , K = tf* V ' b 1 h ,n+l * i h B Finally, H 2 T - * r the values derivatives ah .n+1 \ we h a v e V *h + e ( " 3 cf * h 2 - 1 2 sh . (3.21.4) H .n+1 ,n+l 1 n+1 .n+1 h • ^ h = - n r r r * h % • T $ + 2 t o compute o f V^^h W e w o u the v e l o c i t i e s ^ ^ < and obtain i n each (3.12). velocity t o H , but they would 34 ) (3.21.3) o f the stream f u n c t i o n s would not belong Oh ' H . J = + 2 ' 1^ a h d f i = 1 1 layer ( 3 - 2 2 ) from I f we t o o k t h e vectors which 2 b e i n L (Q); i n o r d e r to force them derivatives to of 0, be „, in o n t o H,. n l,Zh we H^, For project i = 1, n+1 ^ 2 the define n+1 V' lh dx^ (3.23.2) • n+ lh orthogonally dXg 2 The o r t h o g o n a l p r o j e c t i o n f r o m L (SI) - (V J in n+ It remains vorticity is used. V*}* , ih to on describe the boundary This is in the (3.7.2), we h a v e f o r is Be rectangles exactly This l the h the computation when the n o - s l i p = -\ support integrals to of by by of the relative boundary condition the d i s c r e t e analogue taking — = 0 of in that J] J-u,ds, + he h adjacent Thus, M^. i = 1,2, V^^Vu.ds he where space given (3.23.3) the done by c o n s t r u c t i n g (3.7.2) f h is V0. e H , . h oh = 0, 1 onto l of the h h M^. Since boundary, (3.23.1) yields 35 Vu, h e N, h ' is Be then element (3.24.1) by we composed can element of compute of Be. the ^n+1 _ = ih where wall is is |n| point , .n+1 .n+1 ^ - ( 0 . , .,- i/»., ) . 2 , ih,W ^ih,I ( |n| ; to locally a context as r along the normal the first interior point second order formula Wood's formula. triangles or yields an algebraic easily solved boundary T the d i s t a n c e V . 3 A., n+1 ih,I - —=, 2 curved by points In the is, a more elements linear in known as system „„ d i r e c t i o n from Formula I. in finite general set support of of method general, not since very difference up; such the as (3.24.1) which large the (3.24.2) Be, equations Cholesky ^ (3.24.2) can number (see be of [16] for details). By computing (3.19)-(3.21) systems of one i W obtains [ s Q S[* ] a n+1 T[* ] sparse, ] the denotes banded, entries given 7 = [RJ that following appear algebraic in linear T 1 } = [ R = [R i2 3 4 a ( 3 . 2 5 . 1) , 1 = [R J n+1 [ integrals equations: S//* here the column symmetric, (3.25.2) ' ] (3.25.3) ] , (3.25.4) , matrix. The positive matrices definite S, V. l matrices and T are with by t . . = f V<t> .V<t> .dQ , u u 36 ( 3 . 2 6 . 1) s. . = f CV0 .V©>, . + \ <p, .<t>, .)dQ , (3.26.2) 2 V l i J = \ *hi*hj ( ^hi^hj + Q V 2ij = J / ^ h A j where v = (—^-AtA ) u The the (2,2) is unique. Conjugate the quadrature positive definite, v "*hi +hJ V a n d u = (1 + c a l c u l a t i o n of Gauss + then the f i n i t e The systems Gradient (JCG) (3.25) method. 37 (3.26.4) ' )da Ate) . integrals rule. (3.26.3) ' )da is Since performed are element s o l u t i o n e x i s t s and solved the using matrices are all by by the Jacobi CHAPTER III A G A L E R K I N - C H A R A C T E R I S T I C ALGORITHM FOR THE H Y P E R B O L I C STAGE This Chapter is d e v o t e d t o t h e f o r m u l a t i o n and a n a l y s i s the f i r s t stage of the algorithm. prototype equation the the solution curves is of that of the under solution coincides solution. Our weak main s o l u t i o n by methods. is An very t h i s purpose, linear advection which flow. For constant weak (transport) along interesting here is of with the e n d we c o v e r equation, this the equation the weak classical approximation a c o m b i n a t i o n of G a l e r k i n and Towards t h i s as characteristic assumptions everywhere concern the property regularity almost we t a k e of of the Characteristic t h e c o m p u t a t i o n a l d o m a i n Q, 2 n £ fij £ R [Q^ is supposed t o be a r e c t a n g u l a r domain), with a r e c t a n g u l a r g r i d and a s s i g n one p a r t i c l e t o e a c h node o f grid. the p a r t i c l e s the that T h e n we t r a c e b a c k the p o s i t i o n s c h a r a c t e r i s t i c curves the magnitude the p a r t i c l e s , grid points particles at X(x^, t x; particle The + of so x^ at t) the values instant denotes of would be interval x t+r; the values t). position t + T is linear 38 assume transported variable by at the p o s s e s s e d by the Here and at and along in the instant sequel, t of the a t x. . k "representation a is the dependent t +T a r e the instant mathematical a time the dependent v a r i a b l e the p o i n t s X(x^, that at idealization during of the of combination this of physical Dirac delta functions is c e n t e r e d a t X ( x ^ , t + T; t ) . B u t t h i s not very convenient regularize it substitutes smoother functions. which as piecewise element space t ) . Therefore, given with their view, this of limited centroid centered achieved the p a r t i c l e one to extension We choose functions the departure delta amounts velocities. the basis at if f o r the Dirac parcels functions s) of points solution at instant = Yp<- (t)<t> (x k k where a^t) denote s;t). Since ^ . C x ~ ^^ k' ' ^ S the p a r t i c l e equivalent to computing interpolation of ^k^ ^' = X the X(x^, s = t + x the values point that t) ^ s; s)), at the points n e n To d e t e r m i n e s o l u t i o n onto the grid We p r o v e k o f w(x, e l e m e n t s p a c e H^. we p r o j e c t t). S - X(* > k the values x to the f i n i t e u(x, by of s o we must by u(x, of is functions point particles computations, This piecewise physical are advected finite is a point smoother t+x; mathematically. From replacing f o r numerical representation belongs the weights a n d show t h a t of values this w f x , s) aAt) of by cubic a given algorithm XYx^, is u^t) this is spline functional conservative, 2 unconditionally Moreover, for stable sufficiently superconvergent 1. in the smooth L -norm functions and convergent. the s o l u t i o n at the foot of the c h a r a c t e r i s t i c is curves. Preliminaries Let the computational w i t h L i p s c h i t z boundary f, domain R 2 39 Q be a n open subset i s the two-dimensional of R 2 Euclidean space. We now B-splines as I = and 2 [c,d] positive introduce well as of that and / I a = x definitions that A < c = x For positive r,s spline functions I - 11 which < x A = 1 {x 2 and [a,b], be I, J .} J 2J 1 are satisfy = b II < x 21 Let 1 2 ix 1 1 1 2 splines Q = I x I . Let = 1 respectively, < x of their properties. the domain such & partitions some o f assume integers ^ the < < x 22 = d . 2J integers, we of order r over define the linear space of as I l S . (I ) = J S(x ) e H (I ):D S(x ) = 0 for x in [x .,x . ]\. r,Ai l l p l l l ii n +i ' r for i = 1, , I. Here = -{ f e (f (^) H^d^) r _ 1 f r e LPaj continuous, 2 \ , (I ) = L (I ) V P i i Specifically, degree r - derivatives The is '• D f is absolutely 2 D and r 1 a which 2 is up t o o r d e r r linear space defined similarly. product spline splines S of S(x^) € continuous - and a with P°ly n o m i l a continuous 2. spline We a r e 1 S functions of order i n t e r e s t e d i n the c l a s s . © S . d e f i n e d as r,Ai s,A2 40 follows. s of over 1^ tensor of S . ® S . = •{ S(x ,x ) = S (x )S (x ) : S e S . , r.Ai s,A2 ' 1 2 1 1 2 2 1 r,Ai We n o t e . 9 s r, A i w h e r e k = min(r, the of order S a result functions the notation, to Q of the tensor To this spline e n d l e t S(Q) be Furthermore, of ^ case r a Typically (&) h is a special [36] f o r t h e s p l i n e F o r a n y f e S(Q) o n e f . e S. , (Q) s u c h t h a t Sf h k,h ii f, ii s e n f ii , f e L hp p i i i ) Quasi-linearity (f = There Stability. II C f +f ) . 1 2 h product b y S. approximation linear c S(Q) L (£i) assume space P that S of for 1 satisfies properties. Uniqueness. ii) we d e n o t e a lemma w h i c h i n which S i s defined. the f o l l o w i n g cn; , s). due t o Widlund S. p 6 s , A2 p s n, m = 0 , l , . . . k . i) c i/ ' . . We now f o r m u l a t e s, A2 procedure i . to simplify restriction with . V s,A2' 2 that s In S e S +f ) 1 2 P there exists one and o n l y f.. h is a constant C, such that en; n sen; . in LP II s p || f (Q). II l h i p + II f -f 2 II h , 2 f p e ,f i 2 sen; iv) functions Optimal For accuracy. a l l f II f . - f II h p where k denotes * Ch |f | the order k p,k t , of the spline. 41 sufficiently smooth Lemma 1.1.[36]. i)-iv), Let the approximation S satisfy then II f.-f II h p ^ Ch \f\ p,r (1.1) . r 0 < r < k. for An ^ {x procedure alternate a n d {x . } {x ^.> i i l,...,I+r = constructive definition as f o l l o w s [7]. Enlarge involves B-splines 2 . }* 2 and J to non-decreasing of S . r, A i and the sequences sequences {x , } I + r and , the a d d i t i o n a l points being otherwise a r b i t r a r y . For B-spline) and j = of order r 1, . . . ,J+s , t h e i-th B - s p l i n e (s) f o r t h e s e s e q u e n c e s is B. (x) = (x . - x .)[x ., . . . ,x . ](x-x)/\ i i +r i,r i i i i ii+r (j-th l + (1.2) f o r a l l x i n I . S i m i l a r l y f o r B . (y), f o r a l l y i n I . i J>s 2 H e r e [x ^, .. .x^ ^ ] i s t h e r - t h d i v i d e d d i f f e r e n c e o n t h e f u n c t i o n r~ l (x^- x) + = max(0, x^- x) keeping x fixed. From (1.2). the following recurrence r e l a t i o n i s obtained. (x) = B. 1 , X B. Cx; = i,r i i ~ " 0 , otherwise. 1 1 — X x . n+r-i ' X X . - - - x . i i X ii+i (1.3.1) X B. (x) + i , r - i — X — B . (x). x . - x . 12+1,r-i 12+r n+i (1.3.2) We r e c a l l thesis some p r o p e r t i e s o f t h e B - s p l i n e s (cf.[7]). 42 which a r e needed i n t h e Bl) B. i,r BZ) F o r any = 0 f o r x not (x) Cx ,x 1 2 i n [x 11 .,x . ]. 11+r in Q ) E^/V^'V " • j=i I 1 - (1 J 1=1 4) B3) S. ,(Q) = span k,h B4) Let iB. > i,r * the (x )B . l j,s non-decreasing (x)\. '.\ l 'i,j=i 2 sequence {x .} I + r x {x 111 the interior •{B. 1 knots (x )B . 1 V X = ix • y 1 • For J the s t r i c t l y of form J + S 1 B-splines increasing sequence •) of data points, J the spline 2 J 1 S(x l , x ) = V c . .B. 2 fx ;B. Cx ) 1 j,s 2 i,r IJ , (1.5) a g r e e w i t h the g i v e n f u n c t i o n f a t x i f and o n l y Schoenberg-Whitney s o l u t i o n of theorem (1.5) (cf. [7]) i f and o n l y B ° # e q u i v a l e n t l y , i f and o n l y ' 2J X 1 1 < *2J 1 < 43 the unique ' V ( 1 - 6 1 if x . < y . < x . 1 1 guaranties if if l.r<*M> J.s<*zJ> B or, sequence J .} 1, J—1 2 x ix 111 will the corresponding (x )\-\' ._ JrS 1 of 2 X J , V i + 2J S' + , (1.6.2) r V J ' ' ) 2. Description of the Algorithm The C o n t i n u o u s 2.1. Consider the Problem Cauchy equation f o r u(x,i) » at + . u w(x,0) ^ assume t h a t denotes v w the scalar = o Dt in fi , (2.1.1) , (2.1.2) the v e l o c i t y f i e l d i s incompressible, = 0 in Q t o t h e boundary = u. i.e., , (2.2) T o f fi, i . e . , 0 . (2.3) require u ( x , t ) € L Yo,r;l/ <'fi,>,> . compute characteristic the solution curves system o f d i f f e r e n t i a l of X{x,s;t) equations 44 (2.4) ,co a To advection the material d e r i v a t i v e of w i n the flow u-n| We a l s o = = U (x) o V-u and t a n g e n t for i n t h e c y l i n d e r Q^. , d where problem (2.1) of (2.1) we introduce which satisfy the the X(x,s;s) = x We state solution without of proof (2.5) and the (2.5.2) existence summarize some P r o p o s i t i o n 2 . 1 . Under condition solution is in integer s C° (0,T;W °(Q)). regularity (2.4) there Furthermore, i,C of m k a unique (2.5) X(x,s;t) w h e r e X(x,s;t) - denotes f l u i d which i s n as u(X(x,s;x),x)dx s the p o s i t i o n at point x at We n e e d f u r t h e r then for a J J a e H , 1 c a n be e x p r e s s e d = S for some m t x that X(x,s;t) e C°(0,T; L (Q)). a of CD the results. exists assume k £ 1 u is in L (0,T;W ' (Q)), |a| s Jt, t -> D X(x,s;t) now uniqueness t -» X(x, s;t) of the system (2.5) such that The s o l u t i o n of and results of at instant (2.6) , time t of the p a r t i c l e s. the s o l u t i o n of (2.5) which we formulate. Lemma 2.1. \s-t\ sufficiently Under assumptions small, (2. 2),(2. 3) and (2.4) and for x -» X(x,s;t) homeomorphism of Q into itself is a with Jacobian quasi-isometric determinant equal to 1 a . e. Moreover, L~ \x - y | £ \X(x,s;t) - X(y,s;t)\ 1 w h e r e L = exp(\s-t|•|u| 1 n ) . i , «>, Q 45 £ L\x - y\ , (2.7) Proof. F o r that and a l l s,t there is is a unique. mapping in 2. 1 continuous, as guaranties On t h e o t h e r h a n d , follows the and its since X(x,s;t) that inequality (2.5) as (2.6). The transformation is by again continuous. Liouville's J(x,s; t) ~ 1 a . e . X(x,s;t) ^ 1, can being is itself. of one It is is to J(x,s;t) guaranties that and II-II 2,fi proved (see a quasi-isometric is a determinant = 11-11 then to of the above = det(DX(x, s; t)). t -> J(x,s;t) condition (2.2) is yield [25], Th. mapping 1.1.7) of II • II class _ and that simple of the matter to C '*(Q), ml II • II * -> m are m,p,Q check transformation x is that 1 if a.e. the then « condition well to quasi-isometric classical s o l u t i o n of (2.1) is given by (2.8) o is it obtain w(x, t) = u (X(x, t;0)) . It = x. 2,Q . The u n i q u e This and inequality m.p.n Jacobian a • be It one easy a determinant w h i c h maps Q i n t o Q , t h e n o r m s equivalent. as X(x,s;t). (2.3), Gronwall's x -> X(x,s;t) formula exists = X(X(x, s; t), t; s) applying as (cf.[18]) considered neighborhood inverse given 2.1 be mapping Jacobian Proposition Remark. I t this can the v e l o c i t y s a t i s f i e s Hence, homeomorphism. a maps fi i n t o (2.7) 2.1 U(x) f o r w h i c h X(x,s;t) into U(x) well from P r o p o s i t i o n s, X(x,s;t) fixing from Proposition we h a v e neighborhood By ^CxJ (0,7) holds known under [29] very that 46 for weak regularity some integer assumptions. m i 1, u e Z. CO r;l/ ' fn;; e o n J problem (2.8). p and u € \T' (^) o P ( 2 . 1 ) belongs t o L°°CO, T;]/ ' (Q)) 1 Now, l e t <£(•) € 2)(R ), J P solution of the and i s g i v e n where D ( R ) = \ <p e C°°(R ) 2 has compact s u p p o r t t h e weak 2 2 by : <p }• , then f u(x,t)<p(x)dx = f u (y)^(X(y.O;t))dy R R , (2.9) J s i n c e t h e J a c o b i a n d e t e r m i n a n t o f t h e t r a n s f o r m a t i o n x -» X ( x , t i s e q u a l t o 1. a . e . I n ( 2 . 9 ) y = X(x,t;0) 2.2. . The D i s c r e t e Problem L e t W be a p o s i t i v e i n t e g e r , At = T/M and t = mAt f o r 0 s m s W - l . L e t u^(x, t.) be an a p p r o x i m a t i o n t o u(x,t). that the u,(x,t) h satisfies approximate interval conditions(2.2)-(2.4). trajectories [t , s] a r e given, of the points For s x T For T The i n t h e time the departure point w h i c h r e a c h e s t h e p o i n t x a t i n s t a n t s=t error e(x) (2.10) , s At x = A t , X,(x,s;t ) denotes n m trajectory = t , m+i a c c o r d i n g t o ( 2 . 6 ) , by X, ( x , s ; s - x , ) = x - T u, CX, f x , s ; s-c), s-c)dc n o h h 0 s We assume = X(x, s;S-T) X,(X S;S-T) - L n 47 m+i o f the committed in satisfies Lemma approximating the f o l l o w i n g 2.2. conditions e(x) that (2.2)-(2.4), Tn Subtracting |e(x)| ^ r ' ' T o ^ T and i ( e x u(x,t) fulfill m in the interval the error T 1 any two set finite of step from 1 1 that we a p p r o x i m a t e First, this a) d>, (x) i s the bound t h e weak a particle is Towards t h i s • s o l u t i o n u(x, t) b y of projected as (2.8) onto e n d we i n t r o d u c e defined <p (x) (2.11). approximation approximation functions k b) supp$ (x) . ' inequality yields t 1,2, .. .K and 1 £ i with (2.10) |e(e)|de 03, Q < e l e m e n t s p a c e H^. k i t follows (2.11) u(X,(x,s;s-c),s-c)\dc h then cut-off 1) . Iu, CX, (x, s;s-c), s-e) h h process. up, - u(X(x,s;s-c),s-c)\dc Gronwall's instant _) |u, (X, (x, s;s-c), s-c) h h |Vu| + Si 0 Applying (2.6) P Q 1 o H (2.10) is bounded by Proof. is u,(x,t) h then for any integer 03, a by for any real x s u c h t h a t 0 s x ± A t |e(x)| * At trajectories lemma. Assume [0,M-1] and the follows. the the For set k = l basis of £ I, 1 s J < J bilinear in x and x 1 = ^ - h . 11-1 _ n l i _ ' 1 s x n ^ ^x l i £ h,. l 2 n 2 j_ ' , - h l i c) <f>^(x - x^) = <t> (x .), Q^x) k 48 h 1 2 J . -'' V i z x 2j - i e 'J' £ h 2 "i^k^ ^k-i' X .. 2J c a n o n i c a Thus t i ~ = <*k H Notice that intervals — Qj^x) shifted support of 4> (x) t o 1 , 2 k CR ; is clear , 2 C4/u ;j ^ C x j d x k the k J <f> (x) e l / ii) that ( 2 . 12. 1) = 1 , (2.12.2.) = meas(supp<p (x)). is. IS. X, Cx, , t ; t ) be h k m+i m which the . ,h .]x[-h . ,h .] , V i, j . I t 11-1 11 2J-1 2j 2 compact s u p p o r t i n R and s a t i s f y i) Let is suppfy^Cx) [-h have where u (x) reaches consider the 0. Cx k <f> (x), w h e r e is d>, (x - X, C x , , t k h k m+i make s e n s e grid knot x, k - X , C x , , t ^ ;t )) h k m+i m S( •) k the departure p o i n t * is the ; t )) m instant 5(x trajectory t . Let us X, C x , , t ^ ;t )) h k m+i m - measure. is the m+i = Dirac which at of required ^ Another for our * property of algorithm to the f o l l o w i n g : If iii) then supp^Cx - for some k and k m+I supp<p.(x) K n supp^Or X (* .t ;t )) h j m r\ supp<j> .(x) * 0 , J - V X j ' W ' V ; * (2.13) Let us check Assume regularity that At of = i i i )is OCh;, the f i n i t e satisfied. by virtue of (2.7), e l e m e n t p a r t i t i o n one 49 (2.11) has and the 0 IX. Or ,.t h j m+i it ) - X , ( x , , t : t ) \ m h £ Since supp<t>,(x) from (2.14) grids, that in not for k (3.14) to iii) in (I, 0(hAt + m+i + \X(x.,t j + * es i s uniform, m h ~ X,(x h ) - X(x |u - u. | . - A t h oo, Q o f o r d e r 0(h), For this there ;t m+\ - X (x m \ X ( x . , t : t J k m+i m t k k ,t an in compute the ;t Before doing simplification. solution Cartesian components B. Ax, ) 1,4 lp and ^(x^) a B. .(x„ j,4 2q t b y A^and A^ Let " { i ^ - j grid let points at is the main the proposed points introduce we s e t result index j are X, (x, , h k t m+1 t , respectively. be t h e v a l u e s instant t^^, 50 \x^pW 2qt' x 1 m ~ P Such m a t r i c e s a r e of the approximate and A a this to ,; t ). m+1 m notational hk . J points of algorithm X,Cx, , h k some ~ *' with matrices B. 1,4 1 .(x.) 1 ~ g invertible. solution symmetric m )\ )\ clear ensures points (xf, xf, _,), a n d d e n o t e t h e hkl hkz ), g e n e r a t e d b y t h e B - s p l i n e s the g r i d at the us Hereafter, r to s t a t e a d e s c r i p t i o n of so whose m )\ m . (2. 14) i n suppcp. (x). k conditions weak element ;t ;t m+\ structured also least m+i m+i then i t condition at t j ; reasonably exists |u-u, | - A t H h e h oo, Q s e c t i o n and p r o v i d e the ,t - , t ) + satisfied. . . . ,K) j such We a r e now £ J, is ^ ;t ) f o r a n y y m+i m a n d Bj \X(x t h a t x , e suppd> .(x - X,(x.,t ; t ) ) . Next, k J h j m+i m t h a t y e supp<p (x), t h e n i f one t a k e s y f o r x . i n k J i t f o l l o w s t h a t suppd>, (x-X, (x, , t , ; t )) approximates k h k m+1 m assume X,(y,t h £ m \x . - x , | + OChAt J k necessarily any zn+i r\ supp<t> .(x) that (1,...,K) up k at positive definite matrix with entries a = rs below the proof of t h e Theorem (<t> , d> ). The r s product (see 2.1) A[u> ] = [B] (2.15) m+1 satisfies of [B] the theorem represent the approximate that follows. the values at 2.1. instant s o l u t i o n w, (x, t ) a t h Theorem In (2.15) t of the entries b, k the p r o j e c t i o n of the departure p o i n t s xf\ . m nk -JB. Jx,)B. Jx„)\ .' . ,of ' i , 4 1 j,4 2 i , j = l for the space S „ ,(Q). Then at any instant t , Consider the set I J { B-spline basis 0 s m £ W-2, \ 4, n the entries m of the matrix [B] are given by j ^ J / ^ A ^ ^ X ^ ; . t h e coefficients Moreover, Vk i n c"! . satisfy the 2,2 K. (2.16.1) relation A ^ A ^ c " ; = A/"w 7 . 3 For the d e f i n i t i o n of theorem . In t h e m a t r i x A see below general, a d v a n t a g e o u s t o t a k e A= Theorem 2. 1 (2.16.2) m for well the proof structured states that at any the it is grids A f o r r e a s o n s t o be e x p l a i n e d simply of below. instant t the m approximate by solution performing departure Thus, consists cubic w ^Cx, t^ ^) + spline i s updated at i n t e r p o l a t i o n of the g r i d wCx, points t ) at the points the a l g o r i t h m of t o compute the f o l l o w i n g steps: 51 the f i r s t hyperbolic stage Given i) [u J : Vm e efficient of the compute [0, M-l], way o f tensor m performing product by s o l v i n g [c ] this step of matrices. is (2.16.2). using ~ L e t A[w m the The most properties * ] = R . The column * vector R can be a r r a n g e d I and j < j. as I s an I x J as Likewise, m a t r i x C such an I x J m a t r i x R = (r ^J, the column v e c t o r [c ] m 1 s i < is- arranged that CA* = Y . Then, for Finally, I s for J s J, 1 s i solve < j , A * Since l [ c solve ij ] * = I y ij ] ' 1 ~ J ~ J • t h e m a t r i c e s A^ a n d A ^ a r e d i a g o n a l l y dominant, they can be e f f i c i e n t l y i n v e r t e d by Gauss e l i m i n a t i o n w i t h o u t pivoting (cf. [7]). number step The of long operations at this is O(IxJ). At ii) any instant t^, determine by sequential or binary s e a r c h t h e i n t e r v a l s w h e r e t h e p o i n t s (X™, ) l i e . i i i ) F o r e a c h p a i r (xf „ x? x, < x? ,s x , , x. _ , hk hp! hg2 lp-1 hpl lp 2q-l X, _ s x „ , e v a l u a t e the B - s p l i n e basis ng2 2q v. = B using iv) 4 i CX the r e l a t i o n s For j = q, q+1, ), m p l f o r i = p, ( 1 . 3 . 1 ) and q+2, q+3 p+2, (1.3.2). form p+3 , ~ m d . = y v.c . . . J vO i• =p Evaluate 52 i ij p+2, p+3 J Y, d.B . JX? .) . j t J J' J 4 hc 2 2 vi) Finally, . , . . to o b t a i n The u, k solve the system (2.15) in operations O(IxJ) m+1 number of long operations (multiplications and 2 divisions) Here, r number denotes of departure than 49K of t a k e n by the operations points that of + 0(K). the is is iii) of to carry 64K + 0(K). is — v) order standard Thus, B-splines spline steps the This the expensive i n t e r p o l a t i o n procedure; Proof or the s o l u t i o n is of Theorem 2. 1 Mas-Gallic and approximation approximate of the the [24] weak initial combination of Dirac solution condition view local o(x, t) 0 is use bicubic nature is of not smooth. Raviart define w K larger the standard to may at which when t h e d o m a i n not g l o b a l l y we total slightly of the According Raviart the procedure the t h e B - s p l i n e s may o f f e r some a d v a n t a g e s rectangular is point however, point. computations spline than per Hence, number from a computational more + o(r) spline. out bicubic 4r ^y^ a as [32] particle follows. by and a We linear measures k=i for some s e t substituting of. p o i n t s iy^P^^ (2.17) into (2.9) 53 y^ e Q and P « k n a direct calculation k e R yields . By K = I P u ° 6 C x - X (x u (x,t) k Here ^ C x ^ , t ; 0 ; = y In k [32, Th.4.1.] solution i s proved t h e convergence to regularize product with function. degree a have cut-off been proposed we a r e i n t e r e s t e d o f smoothness conforming (2.18) tothe ( 2 . 9 ) f o r a n y <f> e C CR .). F r o m a p r a c t i c a l i t i s convenient however, of 2 of view functions (2.18) . 0 weak . ,t;t)) h k=\ P finite approximation types of i n the literature and which cut-off Next, at instant cut-off (cf. [32]); functions are suitable elements. t o u(x,t) (2.18) by a c o n v o l u t i o n Several i n using point t o work we o f low with define C°- particle as t m+i K u (x,t ^ ) = Tp.u..(t p where m+i w,,Ct ) hk m . ^ k h k r n k=i ;t ) , t ), m+i m m function at k process move vertices their h k {x^}, b y p a r c e l s centroids. the value i s a grid point. to replacing the characteristic of limited We d e n o t e b y w ,(x, t o f w (x, t ) which p m+i (2.19) m+i of value to a weight • « The r e g u l a r i z a t i o n point curves ph :t ^ )) , m+i as an approximate . amounts k u denotes p, and x o f (2.19) along - X (x,,t i s t o be understood w(X, Cx, , t n )8(x particles, which passing through the extension moving with ) theregularized form m+i i s constructed as follows, K u> .(x,t ph r ^ ) = Y.p,u,.(t m+i .^^khkrn k=i )8{x 54 - Xjx.,t h k ;t )) m+i m+1 A *(pl 1 k 4>.(X)) k K k=i We note right Young's theorem 2 is away ° hk m k that each h Y w ,(x, t ph m+i convolution k m+i 2 is ) p o i n t s X^(x,s;t) point a r e i n fi a s o n e d e d u c e s 2 1 t h a t w ,(x,t m+i xf, w h i c h a r r i v e s of . Fig.2 a t v e r t e x x, a t i n s t a n t k we l o o k f o r approximating Since H.. h r d>. (x-X,(x, , t k h t h e n we c a n s e t u ,(x, t k ' ; coincides )) m+i m+i m+i h i t is depicts a ^ process at t thesolution ;t that the m+i ) = w, (x, t ph & the regularization hie space by ' > ' f r o m Lemma 2 . 1 , m+i J representation Now, * and given ) e U ' (Q) f o r a n y t ph graphical P^^^ since — i n £ ( L (R ) ) . B y t h e d e f i n i t i o n o f <p clear 2 i n L (R ) operator 2 m+i m+i (2.20) i n the w i t h <b.(x), k ) where K % ' m+i (x and s u c h t ) = Z k k=i \ Q) ( x ) ' X 6 Q ' ( 2 ' 2 1 ) that (w,(x,t ) , <p.) = Cw, (y,t ), 4>.(y - X, f x , , t ;t)) h m+i k h m k h k m+i m m+ , V k, m, m (2.22) where (u, v) From (2.22) algebraic = uvdQ . one obtains linear K, step i n the proof the entries the ' weights satisfy the system A[u Our n e x t that b of m + 1 ] = [B] . i s t o show t h a t (2.15) f o r a n y k, [B] a r e t h e i n t e r p o l a n t 55 1 values k £ of a bicubic of the tensor prove of spline this at product assertion generality. a = the p o i n t s x < 1 Let < x By of f u n c t i o n s for virtue of a , b be the on r e c t a n g u l a r t h e one d i m e n s i o n a l t h e end p o i n t s = b and i X^. properties g r i d s we case without of the consider w k r J loss interval I K(X ) = I " f <(>Ax - X )<f> (x)dx , V k,r', k can (2.23; r a where (x h ; + k - x k • v the symbol Z, is either then by v i r t u e of derivatives in -{x^}-. generate polynomials k+i h ] ' +'XLk as - x (— h. • )l/> T k+i cubic polynomial uniform each one has possesses the in variable, (2.14) Such the W [ X b first interval i.e. a well that for and second [x collection ,x linear space with continuous > of ^ (X)dx if direction or structured order and ], is grid, an r continuous with breaking r+i x 3 derivatives 56 . Now, any k t h e r e 6 piecewise 4 in coordinate r points e ) c a n now be w r i t t e n piecewise a smoothly KCX^J x ' ] r is that x <t> (x)dx ^(^^) progressively such K(X )4> (X)dx a grid h otherwise. - x J the k' k [x • K(X, Clearly, k € 0, Dropping x > cubic °f up polynomials piecewise to second cubic order at the breaking theorem points ( c f . [7] ) -{x^. , P i,x with the s p l i n e S(x) in S characterize is obvious analog space of (2.23), .Cn). 4, n have one to the In linear of two there x hk r that dimensional exists a cubic case Since combination given of grids a n y S(x) of by e S tensor order hk the in ; = AlJ"] r A are the d i r e c t i o n A = A. a Xf,= then i t f o l l o w s entries -X, ), d> (x)). k r Hence, ascertain takes SCx where the in hk we if Curry-Schoenberg s u c h t h a t S(X?, ) = K()CT, ). I n h S(x) that S JQ) 4, coincides to 3 r> case According values the two to ). I t r dimensional S(x , (2.24) the inner which are products uniform (<P ( in each , (Q) c a n b e e x p r e s s e d 4, n products of cubic x k as B-splines, t h e n f o r a n y k, k Taking hk (^^^ , L . ii L I=IJ=I ^"zhk^ ^ ik' = X X lhk i,4 J 2k^ a D O V 4 e a n d u s 2hk ing (2.23) yields A*® A* [c ] = A[o> ] m m 2 1 Remark. I f a n o n u n i f o r m g r i d i n e a c h d i r e c t i o n i s u s e d , * A. T h i s time the implies step; however, algorithm limitation) several that from a p r a c t i c a l requires one the m a t r i x A has might well be reasons f o r doing t o be point structured tempted so: 57 to then recomputed A every of view, and since grids (this may be take A = A. There a are i) S a v i n g s is i n c o m p u t a t i o n a l t i m e a n d c o r e memory. calculated once and for all at the The m a t r i x A beginning of the the next computations. ii) The conservation section, would is be conserved stability iii) the 3. If better a>(x,t), a s we s h a l l s e e achieved. up In to is condition uniform fact, if A were O(At). - O t h e r and t h e convergence t h e CFL grid of is do n o t less in used, properties u>(x, t) as the suffer. than 1, i n each d i r e c t i o n or then A = A whether not. Properties In this algorithm. section we study Specifically, some we show important that it is features of the conservative and 2 unconditionally stable i n the L -norm. Conservation. Theorem 3 . 1 . algorithm Proof. L e t Let = meas(suppQ^). For 1 s k £ K and m € (2. 15)-(2. 16) conserves £ ^-J^t • k k -m+i1 = A/wm+i Then 58 [0,M-l] Now, recalling follows (1.4), (2.15) and (2.16) and with A = A it immediately r- m+i m k r o i l k Stability. Theorem 3.2. Algorithm (2.15)-(2. 16) is unconditionally stable L -norm. Proof. B y v i r t u e o f Multiplying by ( 2 . 2 2 ) we c a n w r i t e and summing o v e r k indices gives II u.(x, tm+i ^ ; i l s II u.(y,t f , )" h h m )\\ II Zv.wkf 0 ,k Cy - x hk 2 Now,it of remains the ) is to estimate inequality. the the supports may write + 1 To shifted of the second begin support <p (x) and with, of 0 Cy term on t h e we recall - xf that Let <f> (x). right £ hand side suppip^Cy - and denote SE. ,) r e s p e c t i v e l y . Then we SE J From ) ( 2 . 1 4 ) we h a v e that i f supp<j> . r\ supp<p i s J supp<pj n supp<f> = supp<j>Xy - X^j) n k where we h a v e assumed \u - u, | that 59 not empty, then k supp<p (y - X™^) + k _ s 0(h). W i t h O(hbt), this in the information and taking overlapping is finite into account one e a s i l y that the amount of the gets Hence ujx.t h m+i )\l s a + CM) II u.(x,t h ), which m )\\ , (3.2) and so s t a b i l i t y . 4. E r r o r Analysis Let u (x) = m+1 u(x,t m+i according to (2.8) satisfies o> (x) = u (X (x n m+1 where in xf is h E^, h a shorthand where macroelement t m E^ v e r t e x x, , c o n s i d e r k the of the this in the proof J transformation of ) ) = C^OO m f o r X.(x,t h support those n m+i of C4. 1) • ;t <p^, elements m ). For that which any x is, meet the at the transformation x -» x - x By ;t notation is composed m+i k + X* k = x + « the element stability, E, k we d e n o t e m k is m+1 60 (4.2] shifted the s h i f t e d SE, . L e t u s i n t r o d u c e t h e f u n c t i o n o) (x) k . which by a" , . hk As element by 3 i s d e f i n e d by w We also w (x) consider ( x ; = w (x + a , , ; . the approximation i n thef i n i t e ml (4.3) hk o> * (x) t o t h e s o l u t i o n m element space 1 . J^^Cx) i s defined by (cf.[28]) ( u < x J > , <p.) = ( J"(X? ), <fi. ) , V k . h k h h k m+1 (4.4) 2 Let IT^ b e t h e o r t h o g o n a l respect p r o j e c t i o n from t o t h e inner product Cn u, ) = Cu, h Then (^^(x) (u, v), , V* € # X with i.e., h . (4.5) may b e e x p r e s s e d a s «f Cx; = +1 hh h Next, L (Q) o n t o (4.6) cxf; . h f o r a n y x i n £ , we c o n s i d e r theparticle approximation t o (J^*^(x) i n tf, , w h i c h we d e f i n e a s ZT (x) h X where t h e v a l u e s The inner =u l u At HAx hk m k h - Xll ) hk 1 . (4.7) w. , (t ) a r e o b t a i n e d b y t h e r e l a t i o n products hk m on the right by i n t e g r a l s o f t h e form 61 hand s i d e of (4.8) a r e given h ~< ~^ - *^ (y) (y • dy E k By v i r t u e o f (4.2) (4 - we may w r i t e m = y and c o n s e q u e n t l y X the i n t e g r a l s L'E,whJ V x formulation (4.10) formulation according write as (4.8) ( l e t us hand the side recently inner render to , V k . By as an area-weighting virtue of (4.10) hk we , V k . « ),4> ) + (4.10) k may (4.11) m+i, . m+i, . -m+i, . fx) = w f x > - w, fx) . , . .„« (4. 12) n (4.4) is shown t h a t of the analysed of described products become viewed [27]. % computations is be k formulation perform (4.9) define v The can ,<p ) = ( (x % Finally, of a*)*.(x)dx hk'^k + k The %k + the in in inner [19]. by several and a products Morton approximation (4.4) [31] of et the on al right quadrature method the right [27] have hand side rules the f o r m u l a t i o n c o n d i t i o n a l l y unstable. to might We w i s h to 2 ascertain, in the w> (x) b y u^* (x) h m+1 X J L -norm, at a nJ y the error instant t 62 m+i incurred . Towards by replacing this end we 9) set m+1, . m+i, . m+1, . -m+i, . -m+i, , m+i, . w ( x ) - w, (x) = w f x ) - w, f x ) + w. f x ) - w, ( x ; . (4. 13) h h h h By the triangle inequality „ m+i, , m+i, ... .. m+i, , - m + i , II w ( x ; - u, ( x ; l i s II w ( x ; - w, ( x ; i n h + II ZT (x) The e s t i m a t e given by the f o l l o w i n g Lemma in of the f i r s t 4.1. For L (0,T;W ' (n)) m term on t h e r i g h t side (4.14) o f (4.14) is lemma. u(x,t) in and using 1 co . - uT\x)\\ X L°°CO, T;!^* ' (Q)), u(x,t) 1 2 C ° finite elements of degree k, k £ 1 m+i, , - m + i , .„ „ m, . - m , .„ 0) (x) - u, (x)\\ < II u (x) ujx)\\ h h + |wl . |u - u | _ A t + C h l.oo,C? h oo,Q Proof. S e e [31] . To e s t i m a t e observe Let by v i r t u e term on t h e r i g h t o f (4.10) that by t h e t r i a n g l e -m+i w, h L (4.15: t h e method side used o f (4. 1 4 ) i n [27]. We inequality L the f i r s t From _ . k+i,p,Q hand m+i„ _ ,, -m+i ~m+i„ „ ~m+i - u II s II w - w II + II u h h h h us study inequality. l l uV • t h e second we c a n a d o p t k + 1 L term on t h e r i g h t the relationship 63 u m+i„ II . h side ,„ (4. 16) o f t h e above it follows that sTUx) - Sf*Vx;n * ii n so that Using we h a v e (4.12) II iM"; - »l(x n h n to estimate the right , * h hand and the t r i a n g l e i n e q u a l i t y side The first term h h u hk h side + hk - v " Y x + a™ ; i l . ( 4 . 1 8 ) hk + II Ax?; h on the r i g h t (4.17).' «?;» m h of we o b t a i n ii uJVxf; - J?rx«™ ; » s II cx?; - < A x + (4.i7) hk o f t h e above inequality is bounded a s II u (xf) h m From (2.14) - c / Y x + « J ;|| - C|xf - f x + a " ; | _|vw| _ . ( 4 . 1 9 ) hk h h k oo, Q co, Q and (4.2) IX? - fx + a" )l = OihLt + |u h hk H e r e a f t e r we a s s u m e t h a t II (Ax?; h The s e c o n d - (Ax u. I A t ) h |u - u , | £ 0 ( h ) , n + a ,;il hk m term on t h e r i g h t £ C|vu| _ oo, Q hand s i d e then (4.20) (4.19) becomes hAt.(4.21) o f (4.18) yields = J (x? - (x + a™ ; ; - D i A F f x ; - ; d e o n hk 0 Q 64 . , (4.22) where F (x) = x + af, + Gfxf - (x + a*)) 8 hk h hk . a (4.23) Thus ii AxJ; - v (x m a ;n m k 2 s E f ixjj- k E. x - a i [jiVV Cx;;i dedfie m 2 2 k e o k Since FQ(X) image o f £, b y t h e t r a n s f o r m a t i o n x -> X, Cx, t ; t ,), t h e n we k ' h m+i m 6 may change employing is + a quasi-isometry the v a r i a b l e s in of £^ onto the second £ k , £ being k integral to the obtain, (4.20), II v OC?) - v (x + «™ ; i l £ C h A t l l W l l , m m m h hk By v i r t u e o f t h e a p p r o x i m a t i o n p r o p e r t i e s Sect. (cf. Chap.II, 3.2) II v (X?) m Finally, h - v (x + a" )\\ s CAth \\ u> \\ m 2 hk (~u>l (x) +1 - co^Cx), Collecting n , term on t h e r i g h t (4.24) hand (4.16). <p ) '= (Z> (y) - *(y). m k t h e same a r g u m e n t ~m+i w, n m 2,2,Q we come t o e s t i m a t e t h e s e c o n d side of the inequality Using of h a from the s t a b i l i t y proof m+i _ -m m - w, II £ (1 + C A t ) II w, - w, II . n h h n the estimates \(y M (4.21), 65 - X ^ , V k yields ,. (4.25) ( 4 . 2 4 ) a n d ( 4 . 2 5 ) we o b t a i n -m+i u, m+i„ ^ , , „ . . , „ -m m„ - w, II s (1 + C At)II w , - w , II h f h i + C hAtlVwl 2 Now,from H (4.14), m + . + C Ath H w H 2 oo, Q + (1 + CAt)-{ll u " - u w°= w°= <j°, t h e n h (4.27) and G r o n w a l l ' s II u - w inequality K(h,At) following Theorem = exp(CT). 2,2,Q . ]• oo,C? 2 oo, Q (4.27) substitutions into This _ m,Q , 1 (4.28) result brings us to the theorem 4.1. For w(x, t ) in L°°(0, T ; l / ' ( Q ) 2 2 L ° ° ( 0 , 2 \ V ( f 2 ) ) a n d A t = 0(h), algorithm n y 1 , c °(Q)),u in (2. 15)-(2. 16) with C° ,CO finite . give II w II 1 + hAt|Vw| where successive h II s K -7-r— -{h At h + K hAt|Vcj| n II + II w ° - u™ II}- . m that (4.26) that 2.2.Q oo,Q 1 . 2,2,Q J 2 K _ m 3 f A t ; h | II u II + 1 u h Assuming n ( 4 . 1 5 ) a n d ( 4 . 2 6 ) we f i n d _ ™ || < 1 w n elements of degree k = 1 c o n v e r g e s i n t h e L°°(0, T; L (Q)) Z norm with error 0 ( h ) . Clearly, the computations estimate mentioned (4.28) in [27] is suboptimal. show that for Numerical sufficiently 3 smooth functions o t h e r hand, the error i n the previous committed analysis p a i d t o t h e f a c t t h a t we a r e u s i n g the values at the departure points, 66 may b e 0 ( h ). On t h e no c o n s i d e r a t i o n h a s b e e n cubic spline to interpolate and i t i s w e l l known that this might points. In spline £ order in norm. r 1 a> At any i n s t a n t = +1 where S solution is I h spline h operator = j hj that and u(x,t) i s r s oo, solution PS »(X?.) U h interpolant at the points D, (x) the f i n i t e 1 £ s is u> (x) m+X of and p i s hj the approximate i d e n t i f i e d as the 1 - j £ K } 1 space (4.29) , hj f r o m t h e s p a c e l/ = { dimensional € /Y^.Then m+i, . m+i, . m+i, , m+i, . (x) - w, (x) = w (x) - Io) (x) n r where Iu d e n o t e s right hand whereas side the of term important the truncation error A further decomposition - is P S U W . ) , o f w. T h e f i r s t represents the concern us because and c o n t r o l s instance, (4.30) pS(Ax? J h hj - the interpolant second estimation w i l l 1 e n d , we a s s u m e this satisfying, m+1 7 2 the ^'"(Q)) + Iu> (x) iV * by t h e c u b i c h (S^CX^J^Xx) , those we now e s t i m a t e the approximate the cubic prolongation step played algorithm integers t 4 as 0 ( h ) a t m+i o>* (x) h u the r o l e n r and s T; W ' (£})), as h i g h Towards Z m simply into our L (0,T;U ' (Q)) € 4. error to ascertain i n t h e maximum u(x,t) r a pointwise interpolation error L (0, yield (4.30) term on the the approximation error evolutionary error. Its i t accumulates a t each time aspects of t h e method and t h e n u m e r i c a l as, for dissipation. yields = p o A x " ! ; - c A x J . ; ; + (u> (x".) m P 67 - scAx"\;; p(Sw (xl.) - ScAx?.» hj h hj m + Hence (dropping t h e s u b i n d e x Q) . m+i, . u) m+i, ^ . m+i, . m+i, >. (x)\ £ u (x) - Iu (x)\ n oo oo (x) - w, r + |cAx°!; J - a> (X? .) I m hj co + i ( A x ? J - s<Axf.) i + IStAxf J hj The t e r m B l g i v e s hj - SfcAx? J h hj with bilinear f i n i t e m+1 co mX , (4.31) elements £ Ch \a> (x) - Iw (x)\ h j oo . (B1..B4) | 2 , (4.32) 00 by i n t e r p o l a t i o n theory. find (see a l s o \d (X .) n - m J In regard <Ax? J hj IcAxfh j J of i t remains the splines | ca £ ClVul - oo, Q M u . - u| h StAxfh .) | < j oo t o bound 2 . 2 we B4. n . At) oo, Q (4.33) Lemma 1 . 1 , i . e . , Ch |w| 's, From S oo . (4.34) thes t a b i l i t y property one g e t s iscAxf j - S C A X J J I hj From Lemma [2],[28]) A bound f o r B3 i s o b t a i n e d b y u s i n g Finally, t o B2, from (4.30)-(4.35) h hj £ ci<Ax; - <Ax;i h co and assuming application of Gronwall's that inequality 68 w°= . (4.35) co co°, i t f o l l o w s o n l w m i_ m+ij + o- = min(2,s) where (4.36) K-^fh " At 0 n s oo,Q and K a + |vw| hAt) , n oo,Q positive (4.36) constant. The result i s f o r m u l a t e d as a theorem. T h e o r e m 4 . 2 . For u(x,t) e L (0, T; W ' °°(Q) n W (n)), a 2 L °(0,T;V ' '(Q)), s an integer C 1 a min(2,s); algorithm a(x,t) e S,m such that 1 £ s £ 4 a n d cr = (2.15)-(2. 16) with C° bilinear finite grid and A t = 0(h) converges elements on a rectangular in the maximum norm as given by (4.36). Note gives a that super framework, estimate regions will at of if u is s u f f i c i e n t l y smooth, convergent result, in the departure points Lemma 1.1 is local the s a y s=4, finite . Moreover, , we e x p e c t element since the that in order. Remark. T h e o r e m 4 . 2 may e x p l a i n some o f t h e n u m e r i c a l m e n t i o n e d b y M o r t o n e t a l . [27] u s i n g L a g r a n g e - G a l e r k i n C° f i n i t e linear to elements a d v e c t i o n problem. to integrate In t h i s simple the case one i s equivalent they at the foot of the c h a r a c t e r i s t i c curves. u^= u , then yields £ ;-norm o f and e r r o r 0(h*/At). 69 methods were able b u t as to i n t e r p o l a t i o n by cubic splines (4.34) results dimensional compute e x a c t l y t h e i n t e g r a l s o f t h e i n n e r p r o d u c t s ; T h e o r e m .2.1 s h o w s t h i s 2 those w h e r e u(x, t) i s s m o o t h e n o u g h t h e i n t e r p o l a t i o n a t X?^. be o f h i g h with (4.34) estimate I f one t a k e s i n t h e L^CO,T: 5. Numerical Experiments To i l l u s t r a t e t h e c o n s e r v a t i o n a n d s t a b i l i t y p r o p e r t i e s our scheme analysis as well as we c o n s i d e r advection of parameters a cone of this to verify is in a fixed rotating velocity field Q. The are r ) f a p o s i t i v e parameter, X = (X - X ) 1 use 2.9at steps the uniform h = (1/63), with to an the to complete initial shows of a The (5.2) 1 grid on maximum the domain velocity. revolution height activity is the 2 is [-1,1]x[-1,1], at the observed at time e v o l u t i o n of The equal c o n d i t i o n whereas F i g . 3 b revolutions. wiggle point and H = 100, -15h, r = 8h a n d t h e CFL c o n d i t i o n e q u a l = Q (5.1) ) , fix ) . 2 square X r + X 0 u = (-Six We error the / first the is ( 1 The of one " * ' * " L 0 , otherwise , where r validity two m o d e l p r o b l e m s . example r.o; = {-I o(x.O) the of to the 96. displays vertex is base of F i g . 3a time shows the cone a f t e r 6 now of the cone e v e r y number the 48 87 and cone. some Table time steps 1 (1/2 2 of and a revolution) maximum values have maximum and divided by for and the beginning been in terms minimum divided minimum to decrease mass Tw values. by values H. We o b s e r v e dissipative of their at that effects, as The square mass and initial different there is they are time passes. 70 , mass square values, instants mass For of of while have conservation. strong at example, the after , mass the been As very 1/4 of a r e v o l u t i o n the r a t e of time step 6-th revolution. the is 6*10 v e r t e x of to our decreasing 5 We a l s o the cone error dissipation observe is analysis to OflO 3*10 at the the error in of mass end of per the committed time step. 0(h/At) is square 5 that ) per which of This this at conforms case since u(x, t) € V ' . 1 Our second 'slotted* to test idea our 0 0 test high order this algorithm [28] of the complete fact a the high H We is initial to Lip ', where observe in 4. consider the by Zalesak used to is the figures control This is the in the i n i t i a l the wiggles a remarkable The a b i l i t y o f under due Fig.4a and the steps to emphasize 1 , p , but the does space. Figs.4b,4c,4d respectively. ability of the c o n d i t i o n and keep generated achievement at of a designed to t h e scheme to keep the 2 i t s _L - s t a b i l i t y . 71 The =4.2 l/ of schemes. [28]. (See[33]). outstanding The and time We 1979 [37] (m,p) L i p s c h i t z of specifically to of 96. not in CFL number in performance 1 and 6 r e v o l u t i o n s discontinuities. is The is m,P a control those the algorithms. difference equal Lip - under not are condition t o m a i n t a i n the shape is reported now scheme discontinuities. compare c a n a p p l y Lemma 1.1 these strict transport finite = show t h e c y l i n d e r a f t e r 1/8, which one. c o n d i t i o n , where h = 0.01, revolution we s t i l l to results cylinder 1 is experiment this However, We severe corrected order that belong more experiment this shows t h e i n i t i a l of flux with using parameters height a c y l i n d e r p r o b l e m w h i c h was p r o p o s e d b y Z a l e s a k behind Munz is handle Table 2 the scheme strong wiggles shows the time evolution of the r e v o l u t i o n experiment. are ,in some experiment. square of cylinder mass progressively exists is similar time time at the passes. to those initial For steps of drawn from t h i s mass c o n s e r v a t i o n . strong as 24 The c o n c l u s i o n s respects, There every of a results the cone The d i s s i p a t i o n stages instance, to 6 of decrease the rate of -4 dissipation is 10 the 5 after further upper but few after Large and time that asymptotically a f t e r the f i r s t second decrease. external, indicate 10 to revolution. wiggles lower, steps their revolution, appear r i n g s of they are maximum + 10% . 72 The trend to decrease is towards concentrated around the c y l i n d e r and damped and out minimum and the the values to a the slot; results tend N°of time steps 0 48 96 144 192 240 288 336 384 432 480 528 576 Conservation mass Conservation s q u a r e mass .66972175E+04 .99999999E+00 .10000000E+01 .33569349E+06 . 99706674E+00 .99512210E+00 .99338941E+00 .99178720E+00 .99027655E+00 .98883330E+00 .98744122E+00 .98608897E+00 .98476850E+00 .98347391E+00 .98220088E+00 .98094615E+00 99999999E+OO .99999999E+00 .99999961E+00 . 10000007E+01 . 10000021E+01 . 10000032E+01 . 10000037E+01 .10000035E+01 .10000026E+01 .10000010E+01 Maximum .10000000E+03 .92601255E+00 .91213264E+00 .90342497E+00 .89622887E+00 .88989364E+00 . 88427672E+00 .87930044E+00 .87488915E+00 .87096729E+00 .86746379E+00 .86431485E+00 .86144648E+00 Minimum - - -. -. -. -. -. -. -. -. 0 956667E- 2 101286E- 1 109947E- 1 113813E- 1 115438E- 1 118749E- 1 120397E- 1 120812E- 1 120317E- 1 119158E- 1 117518E- 1 118635E- 1 T A B L E 1. T i m e e v o l u t i o n o f t h e c o n e N°of time steps 0 24 48 96 120 144 168 192 216 240 264 288 312 336 360 384 408 432 456 480 504 528 552 576 Conservation mass .24120000E+04 .99999843E+00 .10000013E+01 .10000034E+01 .10000019E+01 .99999732E+00 .99999179E+00 .99999869E+00 .99999837E+01 .99999827E+01 .99999836E+01 .99999872E+01 .99999224E+01 .99999870E+01 .10000062E+01 .10000145E+01 .10000233E+01 .10000322E+01 .10000410E+01 .10000497E+01 . 10000579E+01 .10000655E+01 .10000725E+01 .10000788E+01- Conservation s q u a r e mass Maximum .96480000E+06 .40000000E+01 .995827685+00 .11428758E+01 .94983986E+00 .11180027E+01 .94457594E+00 .11307321E+01 .11386835E+01 .94061257E+00 .93735256E+00 .11387771E+01 .93453740E+00 .11348827E+01 .93203288E+00 .11289023E+01 .92976022E+00 .11239973E+01 .92766902E+00 . .11269107E+01 .11285815E+01 .92572481E+00 .92390286E+00 .11293093E+01 .92218474E+00 .11293036E+01 .92055625E+00 .11287153E+01 .11276561E+01 .91900618E+00 .11262104E+01 .91752550E+00 .11244438E+01 .91610680E+00 .11224079E+01 .91474390E+00 .11201438E+01 .91343160E+00 .91216546E+00 .11176853E+01 .11150600E+01 .91094165E+00 .11167015E+01 .90975685E+00 .11185765E+01 .90860816E+00 .11291797E+01 .90749300E+00 TABLE 2. Time e v o l u t i o n o f t h e s l o t t e d c y l i n d e r . 73 Minimum 0 - 150388E+0 -- 131663E+0 171569E+0 - 189159E+0 -. 1 9 5 7 6 2 E + 0 -. 1 9 7 6 2 6 E + 0 -. 1 9 6 1 7 1 E + 0 -. 1 9 3 0 1 0 E + 0 -. 1 8 8 8 2 8 E + 0 -. 1 8 4 0 1 7 E + 0 -. 1 7 8 8 1 5 E + 0 -. 1 7 3 3 7 3 E + 0 -. 1 6 7 7 9 3 E + 0 -. 1 6 3 3 2 9 E + 0 -. 1 6 0 6 6 9 E + 0 -. 1 5 7 7 5 7 E + ) -. 1 5 4 6 5 1 E + 0 -. 1 5 1 3 9 6 E + 0 -. 148027E+0 -. 1 4 4 5 7 3 E + 0 -. 141044E+0 -. 137493E+0 -. 133902E+0 F i g . 2 R e g u l a r i z a t i o n of the P a r t i c l e 74 Points 3a In °nd c 7S on fo, the C One F i g , 3b The Cone a f t e r 6 R e v o l u t i o n s 76 F i g . 4a I n i t i a l C o n d i t i o n f o r t h e C y l i n d e r 77 Fig. 4b The C y l i n d e r a f t e r 1/8 78 of a Revolution F i g . 4c The C y l i n d e r a f t e r 1 79 Revolution F i g . 4d The C y l i n d e r a f t e r 6 R e v o l u t i o n s 80 CHAPTER IV ERROR ANALYSIS OF THE POTENTIAL VORTICITY AND BAROCLINIC-BAROTROPIC MODE EQUATIONS The o b j e c t o f of the second potential t h i s chapter stage of the vorticities. approximations of is to develop an e r r o r algorithm We also used for estimate the stream f u n c t i o n s the and analysis computing the error the in the v e l o c i t i e s . In 2 both cases of the e r r o r a n a l y s i s the e r r o r namely, error, then employed elements in boundary the the we approximation prove that vorticities, to our case), approximation error the estimates element methods of [6] condition. vorticity are used. of The ), with no-slip the order given in [6] and error by for of of (linear finite free-slip conditions These results linear relative t h e o r d e r 0(h) ) when improvement 81 an no-slip in w i t h the estimate the for the method finite vorticity proved that the a p p r o x i m a t i o n e r r o r of two the when 0(h). - thinks evolutionary boundary [14] function [14] g i v e s 1/2 o r d e r . 0(h 0(h of the r e l a t i v e v o r t i c i t y i s boundary optimal, stream the space 2 is the is in one composed contributed variables is and If approximation is are used; improve In the error the which discretize conditions formulation. i n the L -norm. c o m m i t t e d i n the a p p r o x i m a t i o n as parts, potential is the boundary estimates for free-slip relative conditions of the approximation error is Characteristic method, equations the in circumvents respect the to component the suspected that products. Chap.8]. caused of by error, our a second Galerkin terms the of the variables, those terms. With method leads to a (by r e a s o n s to be time f i n i t e - d i f f e r e n c e component which is of variable 0(h) methods error term i s the first compute analysis the on the The error analysis of of methodology of stage. It is i n h e r i t e d b y many G a l e r k i n w h i c h do n o t based analysis in to potential given the e l l i p t i c the f i n i t e exactly element in the - inner vorticity [10] and equations [35, follows literature. Preliminaries We which introduce are vorticity us the is the standard the n o n - l i n e a r conventional is the c o n t r i b u t e d by the p a r t i c l e a p p r o x i m a t i o n this Our equations of there of Characteristic 1. that is evolution that o p t i m a l a n d more a c c u r a t e than and fact derivative evolutionary However, o r d e r 0(h) the lumping material which i s methods. by to difficulties the seen below) due needed in equations. consider u some the In error order the p a r a b o l i c to definitions analysis of simplify our and the results potential exposition let equation 2 - V u = F i n QKIO.T] , t ul preliminary r u(x, = 0) 0 Vt , = u (x) n (1.1) . 82 We d e f i n e the R i t z projection with respect CVPjU, In fact, solution as t h e o r t h o g o n a l to the inner product Vx) = (Vu, VX) PjU i s HQh p r o j e c t i o n P^ o n t o , ^X e H the f i n i t e so that . Qh element of- t h e corresponding (Vu, Vv), (1.2) approximation elliptic problem of the whose exact 2 solution i s u . T h e L - p r o j e c t i o n IT^ o n t o projection with respect ClT u, h Some results x) = (u, t o the inner product concerning established i n the appendix w h i c h we t h i n k so that w), (1.3) the approximation P , pertinent we i n t r o d u c e (u, X) , Vx e Hh . projection Next, i s the orthogonal to our of this properties error of the analysis, are chapter. the "discrete Laplacian" o f a s a n o p e r a t o r f r o m H, o n t o o p e r a t o r A^, itself, by n C A 0 , x) = ~(V<JJ, Vx), tf>, X e H h . QH (1.4.1) Note- t h a t A,P, h For, with x € f A ^ u , The u, h x) = ILV h 2 . Green's formula = -(VPf, semidiscrete : [o,T] 1 (1.4.2) yields, Vx) = - ( V u , Vx) = (V u, ) equation > H ,, takes 2 X o f (1.1) which t h e form oh 83 = CTT^u, consists of x) • finding ( u or, hf x ) since W r - x = ) ( U h ' F x = F h = " X h T h e c o m p l e t e d i s c r e t e s o l u t i o n if +7 : * e • - 5 that m + k V s.(kL.)F.(t l h h i=l 1 n functions uniformly numbers which + x .k) , n= 0,1,2,.., I 1 . 1 . The order and only r(X) ii) (1.6) and - l s . C X H a r e i 1 m 1 1 on t h e e i g e n v a l u e s -{ ^y T m a r e distinct of real 1]. Definition p if a r e bounded i n k a n d h, a n d where i n [0, = e discretization (1.6) is accurate of if + 0(X ) A and for time as A —> 0 , P+1 0 s 1 < p m . 1 t j E T s a ; = -4iyCe - E Aj- ) * 0(X ) as X -> 0 . A m i=l To X P l j=0 i l l u s t a t e t h e meaning J - o f r(X), m = 1, x 1 and t h e r e l a t i o n s = i) 1/2 , r(X) ' ' v a n d ii) = yield 84 and sAX) an example w i t h t h e C r a n k - N i c h o l s o n n (1.5.2) of (1.1) i s defined by 1 nu = rCkhJlf h rational i) ' n w h e r e w i t h k t h e t i m e s t e p a n d t = n k , r(X) ^ n kA^, 0h H are a l l i n H, , Oh 0 T —> H such On 1 v > ( 1 ) Q if x) h- = (u (x), u (x,0) 1 (F t h e f a c t o r s on t h e l e f t V " V h (f = h' ) scheme. For this \ Y , s(X) 1 - A/2 ' + we s h a l l 2 f o r A —> 0 = give scheme 1 1 - A/2 ' 1 + A/2 A 1 - A/2 = e 1 = A (e = A 1 2(2 - A/2; - X 2 i - A ; + OCA; (e - 2 1) + O C A ) x = 2A C e - 2 - A - A / 2 ; + 0(1) 3 1 4(1 - A / 2 ; observe latter n 1 1 - A/2 Also, ,*3. + 0(X ) that inequalities A 2 |rCA;i < imply the |sCA;i 1, < 1, unconditional for any A. The stability of the scheme. Let us now introduce the f u n c t i o n s v ^, 0 £ 1 £ p, defined by 1! 1 T r v cA; = i A ,( r „a 4 ) > J p> It follows is accurate vAX) Definition accurate from of order P 1.2. p if - X p-1 J 1.1 that the discretization (1.6) if ) as X —> 0, 1= 0,1,.., p n - J 1 1 of order p^, p Remark. C r a n k " > - I. r.s (X) , 1= 0 1=1 and o n l y The time v^X) ^ p Definition = 0(X o X\ J L -Z-Ar discretization (1.6) is strictly £ p if = 0 . 1=0,.. -,P -1 0 Nicholson and 85 backward Euler schemes are accurate and s t r i c t l y accurate of order p = 2 and p = 1, respectvely. 2. Error Analysis In the error potential the of the Potential Vorticity vorticities upper layer, solutions subindex analysis is of the approximate ( C h a p . I I . 3.20) since essentially t h e same. i from t h e equation, solutions of the we r e s t r i c t o u r s e l v e s the analysis are r e f e r r e d t o t h e upper Equations of the lower Hereafter understanding that to layer we d r o p the the variables layer. 2 1/2 Let vector us d e f i n e £ = (1 + | u | ) , where u i s a n d |o| i s t h e E u c l i d e a n norm. n u m b e r s f r o m t h e i n t e r v a l [0, 2n] c o s a . = u ./£ , j = 1, 2. W i t h J of q c a n be w r i t t e n d i r e c t i o n T = (cosa^, , cosa.^, dq ^ „ c Recall L e t oc^, <x^, such this be r e a l t h a t c o s a ^ = 1/^ , notation the total J derivative Dq the v e l o c i t y from Chapter the d e r i v a t i v e i n the Thus cosa^). dq . II as ,. dq(X(x,t;t), an equivalent t) ,„ .. d e f i n i t i o n of the total derivative, Dq Dt In the terms upper of = _dq_ ds v ( x ( X t t ; s ) t v ' s ;| ' 's=t (2.2) (2.1) the equation of the p o t e n t i a l v o r t i c i t y i n layer is: 86 g dq(X( g *' t : s ) q l ' r q b q(x, 0) The = s ) A V q(X(x,t;s),s) , Vt (2.3.1) g=t , (2.3.2) = f . (2.3.3) corresponding (suppressing + F(X(x, t; s) \ , 2 H semidiscrete version of from the n o t a t i o n the dependence (2.3) is o n X^(x, t;s) f o r brevity): W VSTQ*l *n > + -- F (2 = 0 , r Q*(x, 0) Q Q (2.4.3) = f , 1 denotes ''r " |u^| 15 (2.4.2) * here 4 %n 2 1/2 ) a d with From C h a p t e r canonical n t h e f u n c t i o n Q - Q', Q' e H (Q) s u c h II F h Vh = u^ being ((3.18.1)) basis of H, n ? + h ^ f- W ' V that ° Vt + the s e m i d i s c r e t e the expresion a p p r o x i m a t i o n o f u. *n o f Q i n terms o f the is uh Q* (x) = I Q* \(*> k n k Also, for from Q is Chapter (dropping II (3.20) the the superindex 87 • Vt complete n e T . discrete * for brevity): equation Q (x) - Q(X") n+1 ( R V Z> (F K 0 ; , T where y^, y ^ H r e parameters Q(X^) such (QCx, 2 h V f r o m /'O, 1] s u c h t h a t t) V = (2.5) from Chapter = Q(x, , Q(^ ) ), + ' &H We f u r t h e r r e c a l l is H ' ^h Oh a Ux) A CJ(vlQ n + III, S e c t i o n 2, £ y^, and that = I Q 4> (x> n k k that t ; , <(> ) = (<?"(>;, 0 C y n k - k t n + 2 ; r ^ J ) , Vk Hence A[Q] = [B] , and from Theorem 2.1 of Chapter are g i v e n by c u b i c s p l i n e III , the e n t r i e s i n t e r p o l a t i o n o f A[Q ]. n Lemma 2 . 1 . Assuming that bounded at any instant scheme (2.5) is unconditionally stable 2 in the L -norm. Proof. T a k i n g i n ( 2 . 5 ) jn+1 <P = y Q (x) + V Q(X^) n + h 2 2 yields + A )\V(? Q + n+1 H 1 o f [B] 7 Q(X^ )\\ £ \\F \\ \\<t> 2 2 n+K 2 h 88 t, the Now, (Q , since Q(X?)) < l/2(\\Q n+1 II n+1 + IIQII ; a n d t h e 2 second 2 h term on t h e r i g h t IIQ n + I H 2 hand s i d e - HQH is positive, < 2AtCIIF 2 n + K llCllQ we f i n d n + I H that + IIQIi; . Hence IIQ H n + 1 s s where ||Q|| + 2 A t l l F ' | l n + (l + C A t ; i l Q H the s t a b i l i t y We n e x t which scheme. theorem of , C (2.6) Chapter III has the accuracy we r e s t r i c t ourselves To t h e a u t h o r ' s scheme to the knowledge of used to t h e scheme the case classical this i s analysed C h a r a c t e r i s t i c methods. been to • t o prove corresponds Crank-Nicholson - n + inequality. wish simplicity, 1/2, + 2Atllf ' ll n obtain the l a s t For C is that y^ = n^ = Crank-Nicholson the f i r s t i n the context We p r o v e (2.5). f o r At time of that Galerkin = 0(h), our 2 method g i v e s an e v o l u t i o n a r y e r r o r of order OCAt + h) instead ) a s one would expect from t h e C r a n k - N i c h o l s o n scheme. 2 o f OCAt On the other easy to Euler scheme while our hand, following establish scheme is that of the steps the time order of accuracy 0(M+h). In this our proof of the sense, it is backward the latter i s o p t i m a l when a p p l i e d i n c o n j u n c t i o n w i t h o u r m e t h o d , the Crank-Nicholson main accurate theorem than the shows, scheme is suboptimal. the Crank-Nicholson backward Euler steps. 89 scheme However, scheme for as is more larger time In our analysis we need L a p l a c i a n operator which i s Lemma 2.2 The a property discrete Laplacian positive inner product space with respect definite ) x £ ^Oh ^ Let e = -(Vij>, Vx) X s u c in discrete lemma. -A^ operator operator is a considered 2 to the L -inner Proof. R e c a l l t h e d e f i n i t i o n o f A, , h h the formulated i n the next self-adjoint, (A iP, of as product. i.e., , <P, X e HQh that n then x) (A ip, h so = i f A^i/( = 0 , By Vx) -CV0, virtue = (*p, Lhx) t h e n ip = 0 . of Lemma = (ip, yp) exists a = H0II 2 , m there 2.2 set of eigenfuctions of functions KO "{XJ^J w n i c are n the orthonormal -A^, with the a s s o c i a t e d s e t i^j} of e i g e n v a l u e s . The f u n c t i o n s •{Xj} c o n s t i t u t e a n o r t h o n o r m a l b a s i s o f /7_, . T h u s , a n y V, e Hn, i s expressed Oh as h = V Assume that G(X) I h' k V*k (V is an h * arbitrary function 2 spectrum of -A^ . F o r v e L (Q) we G(Lh) = I G(-X j Oh )(v, J 90 set )x, , Xj J J defined on the then IGCA. ; v l l s n and by t h e P a r s e v a l ' s IIGCA, ; i l llvll , h relation \G(h,)\\ = max\G(-\ .) I h j j For example, I I I I = max\-X .\ . The s c h e m e Q where the (2.5), with y = r boundary 1 / 2 c a n b written e + ksCkA, ;F. (X? ), h h h = r(kL.)Q(X?) h h n+1 = 2 (2.7) +1/2 conditions have already as been imposed, so n+1 is Q i n HQ^* I + -^-A — - r(kA.) = r 1 S '*V = Lemma 2 . 3 . A t d e n o t e d b y k, ~ — * 2 + = A A and - J h . £—+ , 2 s J s K_ , 4"* • 2 J " - -k- > J r defined l / i t h P^ 2+2 s 2, the following l s J * which ' 9 ) 1 s s s , (2.10) .Sf & f \L X where, for our problem, the constant C = KRe , constant ( 2 holds. S X o . S q - q\\ + h\\V(P q - q)\\ £ Ch \\q\\ IP K by (1.2) and q € H (Q), inequality (2.8) is independent number. 91 with K of h, a n d Re is the another Reynolds S e e A p p e n d i x C h a p t e r IV Proof. We now with proceed respect solution the rate Q(x,t). If the smooth of by in the of of main in time that to error, of is the type vorticity. which no-slip boundary t h e n " the the is is the error free-slip approximation error 2 ) error fact that the regularity respect conditions with in with is 0(h). no—slip this case conditions introduces regardless the solve of type This of error 92 term boundary finite while elements in to shows linear finite loss and potential 2 L —norm, Chapter boundary t h e e q u a t i o n o f q(x, t). by conditions belongs of an the linear boundary q(x,t) approximation to and the produced time of space,depends that with a perpetrated prove we boundary convergence in rate error for 0(h approximation in the g l o b a l is conditions conditions Thus, conditions elements of q(x,t) quadratic of chapter approximate (see t h e method u s e d t o d i s c r e t i z e the s o l u t i o n on this the space follows) q(x,t) of vorticity and A further analysis approximation result convergence potential 2.1 Q(x,t) space. the function Theorem convergence that state to reasonably linear to • the is rate due to (Q.) a c c o r d i n g II. The of conditions with the of the to particle order employed 0(h) to Suppose Theorem 2 . 1 . i) € q(x,t) ii) q e t .... m ) 2 S 3 d q .2...1. ^ - e L (H ) , ,3 dx r denotes the derivative along the characteristic Then II? - 0 1 1 , C^h*{ \%n ^ 2 where } S + C and C^^ are directly free—slip bounday + ^ Si2tQ "^Th (H ) for 1 £ s £ 2, Vt e 10, T] a> 2 2,„s. „2 d q .2,.2. e L (H ), V ^ - e L (L ), . 2 dx T curves. 1 S L (H ) dq —r— " where p| W ' (n), ff (T2; boundary conditions 02 k 2 t + 2 h C ( 0 3 H H S according and to Chp.II L ^ L , Q proportional conditions + llgll» ) = 1 (3.16.1) S • T to Re, s (H ) ( + - 2 and s for U ) = 2 no—slip and (3.16.2), respectively. Proof. Q Set n - q(x, t ) = (Q n - P q(x, t )) 2 = e" + P " + ( P gC"x, t ) - q(x, . t^) (2.12) B y Lemma 2 . 3 l l p H £ Ch \\q(x, n It remains S t n n = C h { llg ll S s 2 Q to estimate 9 n 0 e H„, + llg " 2 . We f i r s t t L note r / / s ; } that . (2.13) the R i t z Oh p r o j e c t i o n o p e r a t o r P ^ g (x, satisfies t> = ^(x, w t h e f o l l o w i n g Lemma. 93 t) € , for g e , Lemma 2 . 4 . Let P^ be the Ritz projection operator defined in (1.2). Then M* a an p eJS= e L_ 2 dr dx q Let Proof. . q be Sx a differential + \^\ ) . of arc along . • the characteristics. |5x| = LU1 2 1/2 Then « / ,. q 1 liaxi^o » L P q(X) virtue interval of the Lemma, n = = \h ( F p ^-af + \ ^ " Wi* ( P q* = q - q', q' 2 S each rds * p (\(-.t ),s)\ ;S s=t <I<X <*.t:s).s)\ 1)} e L (H ) in ^ = ~ 4s 1 = F. (X, (x, t;s),s)\ h h s=t where satisfies t) equation n+1 * W / w, (x, ,/ t h e s e m i d i s c r e t e It , t ST) dP q* = ?—-w-— h ' Wh + - 5x) above 3w ' - A t 2 £ (X 1 1 5 x 1 - P q(X rjr-, 2 By T C T M * P q (X) - P q (X - 8r) 1 im i' |Sx|-»(A At->0 lim I - q (X) e h sst t * s,t £ t , , n n+1 such that q' \ = q b (2.14) Vt € 10, T] Remark. N o t e t h a t A, w, makes s e n s e b e c a u s e w, (X, ) i s d e f i n e d b y h h h h - 94 wAX?) h h (VP^, The is , w h e r e b y P,q*(X?) I n = P,q*(X?) 1 h V) solution given, 1 / 1 + = (Vq*(fy, X V) X we d e n o t e , X e H Qh l / o f the corresponding , x£ € Q . totally 1 the s o l u t i o n o f discrete scheme according t o (2.7), by = 2 - + ( 2 - 1 5 ) w h e r e £ , , = r ( k A , ,) a n d kh h S..F.Cxf; kh h h = sfkA, h We f u r t h e r d e c o m p o s e 9 = (Q n n - l/ ) h CX. (x, t h n+1 +k/2),t+k/2) n n 9 as n + (W" - wjx, 1 t h Here = z" + t )) n R n . (2. 16) II - (2.17) x z " + = ^u2<x?) kh h 3 + u s . ( i khh h - P ;|^ , 2 3T 7 hence i+l„ ^ IIZ" II +I * ,r- h,;/«ni. . „ / T h t h e scheme i s stable, also, we a s s u m e that theorem o f Chapter IZCx";il n pieces I hk Since All _ o »fl<7* ii £nzcxJJ;ii + kis^i ie ma - ry|f n~ together then |£ | i s b o u n d e d . I I I we h a v e * | £ . . l < 1 a n d |S | < 2 ; kh kh (1 + CJOIIZ !! . 71 give 95 From the stability \\7 \\ - n+1 where C l i z " l l £ CkWT^W + kC'WCI and inequality C are - Pjp- 1 constants. II , OT By virtue of Gronwall's a n d Lemmae 2 . 3 a n d 2 . 4 IIZ""*" !! £ C h | l l q l l 1 + S 0 w h e r e C = KReexp(K'T), 1 | -4 K a n d K* _ L L T L Cff ; 2 } S ' are constants (2.18) independent of h and A t . The m o s t d i f f i c u l t estimate of R. n n + = 1 R n + 1 V Note part of the proof that R ^ t - w,(x, h t .) n+1 X + (E.,w,(X?) kh h h satisfies n+ - w,( , h .) n+1 i s concerned w i t h the = E , , C W C X " ; - w, e x " ; ; + kh h h h + kS.jP.Z, pkh 1 h ST - A A,P.q*)(X, ) H h 1 h n u (Al, To e s t i m a t e A l we n o t e that \\Ai\\ £ \\w(x?) - wjx?)\\ h h h \E,,\ kh < 1 , £ iii/cx"; - W\/X";II h h h term on the second i i i / r x " ; - w.(x?)\i h h h according the last proof to £ (i the s t a b i l i t y term on the r i g h t given in Chapter A2, A3) (2.19) . (2.20) so + llw.fx"; - w.fx";il h h h h The f i r s t i n e q u a l i t y i n (2.20) + CJOM/ 1 theorem hand III 96 - wjx, h of side to . t ;n , n Chapter of i s bounded by (2.20) evaluate III. (2.21) To bound we m i m i c (4.16). the Thus, identifying functions Chapter III u, (4.16) Present (x) J^Ax) w, (x) > w, => w > h^^h^ from ' v h we h a v e Analysis h h (4.16) and (4.17) o f Chapter » V C X " ; - wjX?)\\ h h h h U III that £ UwjX?) - w*(X?)\\ + \\w*(X?) - v.(X?)\\ . h h h h h h h h (2.22) From (4.7)-(4.10) of Chapter V Next, we n o t e Chapter III that takes III it > (w = + follows V > * • v the function v(x) defined the following expression t ) - a(x, t^) by in (4.11) the in present situation: v Using n + 1 (x) = q(x, the procedure n + 1 Pl described in Chapter III a n d Lemma 2.3 yields llw. (X?) - v.(X?)\\ £ CJik\Vq\ h h h h 2 Thus, from (2.20) and (2.22) \\W(X?) - wjX?)\\ h h n + C-Wc|v*g| 2 oo, Qj. + _ + C.h k\\q(x, t ; i l _ „ Q 3 n s, 2,u S oo, T t h e term A l i s bounded by £ (1 + C.kH]/ 1 1 - wjx, t )\\ + n n t )\\ , V[t , t ]. n s, 2, Q n n+1 CjfkWqCx, 3 97 (2.23) The t e r m A2 may b e w r i t t e n a s Taylor expansions Y(X(x,k;k),k) 2 k along the c h a r a c t e r i s t i c curves = Y(X(x,k;0),0) + k-^Y(X(x,k; 2 Jc Y(x(x,k;s),s)\_ d ds +J ^0 =Q s), s) | s = Q + 2 2 / K ycxrx,ic;0;,e;de ds give * 2 = " l{ " P "V + tn+2 et — ^ , tn ' T ; T Likewise, - d T t o e s t i m a t e A 3 we e x p a n d Characteristic curves tn U. n 's-0 A .3 ^ q CXJx.t dx r J k ;T),T) i n Taylor 2^ .2 's=0 . (2.24) ds >> dr\ J series along the and o b t a i n 2 1 tn+2 dx h 3 2 ( t + n 4" - ^V/172 + 98 i N } • (2 - 25) Then, A2 + A3 = 2 k j - l-jr j=0 J d * <%) + * J + tfj + J , q , (2.26) ds J where and , i/kA^P* , stand f o r the i n t e g r a l remainders TJ^, i n (2.24) (2.25). Let wjj d e n o t e Pf^—. <J(^> • The f u n c t i o n s KkA ds J A)w , J J admit a s p e c t r a l r e p r e s e n t a t i o n i n terms o f t h e e i g e n f u n c t i o n s of kA^A^. I n d o i n g so, the factors terms o f t h e e i g e n v a l u e s o f k A^A^ b y lj(kA^A^) are given in 1 (X) = 1 - r(X) + s(X) = 0 Q = 1 - s(X) + (X/2)s(X) = 0 1JX) 1AX) 2 so W ' that = 1 - s(X) = v the f i r s t ' X 1- X ' term on the r i g h t hand side of (2.26) becomes j=0 where J ds J i s a p o s i t i v e constant We now t u r n less than our a t t e n t i o n towards 1. t h e terms 3? n „ „. Note that tn+2 + IIK^H s e y e 3* \\p - U L ||ds . ^ tn Also 99 ds (2.27.2) llft"ll £ In view that of c 6 \ k (2.23) l l i i in and n+1 A 2 ds P " (2,27.1.2.3) |£,, | < 1 because \\R \\ = l l l / H h 3 t h e scheme is into account we h a v e 5 n 2 .2 * V^l-V + C k sap II t£s£t , n+1 4 taking stable, ^ _ ,s,„ , 1 - (2.27.3) - wjx, t JII - III/ - wjx, t )\\ £ h n+1 h n n + 2 h ' S and C.kW]/ - wjx, t )\\ + C.hk\Vq\ _ 2 d ^ oo,Q + C„h kllqCx, 3 An+1 ^ r + C 5 k ds ^tn < l^-iL* 2 1 . ... t ;il n s, 2, n ^ „ + , 3 ds II V H*,^"' ^ Hh 1,2 ds 1 (2.28) d a* %-\\ , 2 ds 2 It remains (1.4.2) t o bound the terms ,2 * h To e s t i m a t e is d a* a n d IIP, 5-11. 2 , 3 ds 3 From we h a v e IIA P P^o IIA.P, h 2 in (cf.[38]) IIP \\P J -±4_|| 1 , 2 ds q ds then 3 3 * ii * ds By G r o n w a l l ' s d 3 II we n o t e t h a t there exists q 2d 3 # H D _ ,2 * £ IIV -4-ll . ,2 ds by 3 * 3 * cnypI- -^n s C I I V - ^ I I . ds i n e q u a l i t y (2.28) 100 3 Friedrich's C independent d ds d * IIVP^-^-^H £ I I V — - ^ l l ds ds the d i s c r e t e a constant * becomes . Since inequality of h such that s C,h t llgll.co. <>. + C_h llq ll , + 6 n L (H ) T 0 s,2,Q. - w.(x, t h n+1 1 t n S S n + IIV ds ds' n c m \v \ 8 n q (2.29) Hence, putting yields together the estimate (2.13), (2.16), (2.18) and (2.29) (2.11). Remarks. Notice i) the that derivatives flow. For less rapid of along so The along last curves the the evolutionary ). F r o m in evolution of error of this then point the and time e r r o r remark ii) one is the curves of than view, time of of in of the q is the t time more a c c u r a t e w i t h l a r g e r variables this elements, if that time 101 the steps one along term At gets the belongs than which i s takes scheme particle rather = to given ui 0(h ), closer to steps. the backward E u l e r now Oik), deduces from thereby, Crank-Nicholson iii) A similar analysis the norms larger the approximation of "optimality" with larger that of the f i n i t e its 1, of the variation curves comes the f l o w ; s £ the the use (2.29) 1/2 m flows permits the approximation e r r o r of s are characteristic characteristic term characteristic 0(h the 2 accuracy. to by multiplying k dominated such approximation to terms t h e scheme l o s s of ii) q convection direction, without the which is scheme optimal. the C r a n k - N i c h o l s o n steps. shows From scheme this is Our iv) methods inner analysis for the products proposed shows that Galerkin-Characteristic convection-diffusion are not by Morton et exactly al equation calculated, [27] and Benque as et in in al which the [4] the methods , lead to 2 an error estimate appearing in regardless of might the the norm w i t h a evolutionary the degree of term of component the f i n i t e of order the elements used. 0(h) error, Thus one t h e a p p r o x i m a t e G a l e r k i n - C h a r a c t e r i s t i c methods 2 methods i n t h e L - norm. Error Analysis In this second of t h e B a r o t r o p i c and B a r o c l i n i c Modes section barotropic, $, order we elliptic II). of the approximation of Our a d e t a i l e d account their respectively. CV 2 analysis of completeness and the equations Chapter For study approximation a n d b a r o c l i n i c , * , modes. of 4> L - classify a s 0(h) 3. in is elliptic based on problems we w r i t e h e r e element of B o t h modes s a t i s f y (recall equations i t c a n be f o u n d finite error the by f i n i t e the equations to (1.7)-(1.9) standard in Ciarlet the theory elements; [12]. satisfied by approximations Thus - A ;* 2 s = 0 , in Q_ , T 102 (3. 1) (V - A j)* 2 = 2 a Q - q l 2 , in Q T (3.2) * |_ = 0, a I V $ = b , 2 Vt in , Q r (3.3) * l = 0 , Vt r , where H l l q H\ - B 22 H + q ~ +H 2 ' F ( Cx, t ; + c a ; * *fx, t; = * fx; , a a n d C(t) is V where ¥ , sh For rv^, C V instant = (b , V where * h ah Note element A C* £ 2 + %n> V x2( with c - ch ; % + € S any II a h 4 ) (3.5) (1.6). element approximations ™sh> V ( " s d e t e r m i n e d b y Chap. The f i n i t e 3 = 0 satisfy ' *h V *0h ' € <f> ) = h 4> ), h - 6 ) 1. e !P ( h e r e a f t e r we d r o p , ( 3 % -< Q - l e H 4> ), Q , h 2 0 h the superindex V0 h £ // , 0 h , n), (3.7) (3.8) a n d $ € /7_, . Oh that 0^ - Q approximations £ in (3.6) of - 103 and in b ft in (3.2) (3.8) and b are finite in (3.3), respectively. Due to this fact some preliminary required before proceeding with the standard L e t <p(h) a n d <p b e t h e s o l u t i o n s 2 V <p(h) + c<p(h) = & , work is analysis. of the e l l i p t i c problems i n fi h (3. 9 . i : <p(h)\ = 0 . T 2 \7 <p + c<p = 6 , i n fi (3.9.2) <p\ = 0 , r where c i s a p o s i t i v e c o n s t a n t , such that element 9^ —» 9 weakly approximation to in 1 1 9. e H, c H (Q), 9 e H (Q) a n d h h 2 L (Q). L e t be t h e f i n i t e <p(h) i n Since fi is a convex 2 polygonal domain, t h e n b y r e g u l a r i t y t h e o r y (p(h) , cp € H (Q) f] and Hg(Q) * 2.2.a 9n where C are constants Lemma 3 . 2 . Proof. From £ c 2 2 -V (<p - <p(h)) " - independent <p —> cp strongly (3.9.1) l l e ( 3 - 1 0 - o f h. in H^(Q). and (3.9.2) one o b t a i n s + c(<p - <p(h)) = 0 - 9 h i n fi , (3.11) <p - <p(h)\ = 0 104 2 ) Since the 9 ^ —> 6 w e a k l y i n L (Q), <p(h) —> <p w e a k l y i n H (£2) a n d b y usual compactness argument lim \l<p - <p(h)\\ Lemma 3.2. (3.9.2), Let <p(h) and &p be the solutions Let <p^ e respectively. approximation " m = 0. be the finite to <p(h). Then the following * ~ ^i,2,n - V c , l e V"J where C = C(Q) is a positive of (3.9.1) and ^ ( h ) inequality element holds • ~ ^i,2,n <3 12) constant. Proof. L e t u. , </>, € H . , t h e n h (v( h Uh n - <p ), v u ; 9h h h c(<p - <t> , u ) = + h h h ( 3 . 13) (V(<p(h) - <t> ), Vp ) + c(<p(h) - <t> , p ) . h Taking h h = <p^ - <f>^ a n d u s i n g Friedrich's h inequality (3.13) yields ^h-^i,2,Q^ i^ K We n e x t ( 3 write " Vl,2,fi From -^i,2,n- (h) * I* - *< ti,2,Q h) (3.14) 105 + ** (h) ~ Vl,2,fi - 1 4 ) " - *hh,2.a * ^ ~ * ( h ) l i.2,a + 2 K i n f ^ ( h ) *h 1.2.a - l *n On eH (3.15) It remains (3.15). on t o bound Note (3.11) that the f i r s t term on t h e r i g h t by v i r t u e o f t h e Lax-Milgram hand side of theorem aplied one o b t a i n s " • ^ l 2 , 0 l 1 4 j f 3 , f l - f l h - U D l From t h e d e f i n i t i o n o f t h e d u a l norm by a n a p p l i c a t i o n o f Schwarz ( s e e Chap. and F r i e d r i c h ' s II) following i n e q u a l i t i e s one gets '* " * By s u b s t i t u t i n g We ( h ) i *V 1.2.a (3.16) into finite element "V • (3.15) y i e l d s a r e now i n c o n d i t i o n s standard 9 to apply (3 (3.12). theory 16) • t h e arguments approximation - for of the elliptic problems. Theorem 3 . 1 . Suppose i) i) Conditions ii) Regularity - iii) of Theorem conditions 2.1 (3.15.3) hold. and (3.15.4) of Chapter hold. Then for any instant "*a " *ah" J * C / { *1 Q t e T n ~ 2Q ti + 106 W*l ~ V'-.Q. } + *" °" ' T II II* - * II £ Ch i Z ™ s - * s h * * where are dependent higher lib, II + | v f g , - q.)\ t\ 1 Z 2 n y 03, _ i + H.O.T. J / W 3 / 2 . 2 . T ' C positive ( constants independent on Re, g = 1 on the boundary order of h but T and H.O.T. 3 - 1 ? ) directly stands for terms. 2 and * i n t h e L -norm Proof. We f i r s t e s t i m a t e t h e e r r o r f o r * a using d u a l i t y ideas. method (see C i a r l e t doing, I* S p e c i f i c a l l y , we e m p l o y t h e A u b i n - N i t s c h e we h a v e - * a [ 1 2 ] , Theorems t h a t a t any i n s t a n t .11 £ ah 3 . 2 . 4 and 3 . 2 . 5 ) . In t tf,hll* - * JI, . _ , 1 a ah 1,2, fi (3.18.1) II* - * J I £ W _ h l l * - $ j | , h 2 h 1,2,0. where are constants From "°"l Lemma 2 fi * Chap. II, '*a ' n 3.2, HQ(^)> Sect. *ah" so (3.18.2) i n d e p e n d e n t o f h. the equivalence the approximation between I°Ij 2 properties of fi a n 0 " H^ ( S e e 3 . 2 ) a n d T h e o r e m 2.1 we g e t * l C h 2 { "l Q ~ <V + { " h" | V ^1 b II ++ IvCg, - q )\ II* - #|| £ CJi \ lib. | V h 2 I h 2 ^ 1 1" 2 0 -*2 o3,Q } )[ T oo,Q n T °- + } + + H.O.T ' 2 The estimate of in t h e L -norm is given a p p l i c a t i o n o f T h e o r e m 8 . 7 i n Oden a n d R e d d y 107 by a [ 2 9 ] . Thus, direct 3/2, Now, and f r o m Theorem (3.22)) 3.1 i t follows 2, Q and ' the inmediatly relations (Chap. II (3.11) that (3.19) To estimate (3.12), the velocities (3.23.2)) and Theorem 3. 1 barotropic the this relations chapter to (Chap. that II get (3.20) 0(h). l states use (3.16) of llu. - U.ll = l we the approximate a n d b a r o c l i n i c modes c o n v e r g e solutions of quadratically the in the 2 L —norm is to the obtained are for linear velocity exact combinations components the reasonable to expect fact, 4. N u m e r i a c l Experiments In order first to v e r i f y the they shows is stream that will the of functions modes. functions, that the v e l o c i t y same d e g r e e expressed stream (3.20) of of are of In The the approximate derivative rate. solution. in because However, terms then it converge rate convergence of is they since the the first numerically with a lower of convergence error estimates of order. the v a l i d i t y of our 108 for the vorticity results the of and the stream some n u m e r i c a l purpose situation of of such a illustrate function, examples. experiments baroclinic our Hopefully, with computations oceanographical We have boundary to be calculated II). case higher In The is as The our experiments, of the is the domain r o t a t i o n a n d number numerical homogeneous c a s e . an exact solution Under of feature these (OJ,I/») solution is no-slip is because test the is, . *'r *n'r = 0 109 • estimate (1000 square the it do Note not is easy 2 2 sin —^— , for order that km from introduce nonrotating numerically UX (see a 2 nx have of with V 0 = w , 2 l —^— they to a space layers to The procedure. density. conditions to ip(x,t) = A(t)sin = thesis. the one error solution. The e x a c t the realistic unknown, our fi our additional with are because no r o t a t i o n a n d c o n s t a n t any of to i n the f u t u r e . numerical xlOOOkm),with analysis, wish more first O(Reh) no m a t t e r i f u b e l o n g s 1. 'real' this in method a only in be c a r r i e d o u t the boundary one than we the that c o m p u t a t i o n a l l y more d i f f i c u l t . par second simulate method the here emphasize analysis two r e a s o n s . condition is to instead this test the v o r t i c i t y at Chapter this to condition for valuesof not proposed situations will decided t h i s boundary is theoretical method present We s h o u l d ocean, Galerkin-Characteristic we to and choose calculated _ 2 27rx. 7ix_ 0 2nx_ o)(x,t) = — — A ( t ) [ c o s — ^ — sin —j— + cos—j— _ nx. sin —^—J , u(x, t ) \ = u , r where A(t) = A(1 - b exp(-t/T)). B y d i r e c t s u b s t i t u t i o n o f w a n d iji i n t o t h e e q u a t i o n = D W Dt A.A + H F o n e c o m p u t e t h e f o r c i n g f u n c t i o n F(x, t). We h a v e h, run several Reynolds number At C F L = max|u|—r—. for Re, a n d C o u r a n t — F r i e d r i c h — L e w y Tables . h experiments w i t h d i f f e r e n t values 3,4,5 give ~, the experiments, decreased always h 2 i n order constant. takes t h e v a l u e s o f ln— —~ ~ " f ~ h " f j l f | | A(t) = 0.64A, w i t h A t c o n s t a n t f o r except those t o make with h = 2 the s o l u t i o n Consequently, the values A linear in condition t h e v o r t i c i t y a n d s t r e a m f u n c t i o n i n e a c h e x p e r i m e n t . The e x p e r i m e n t s were r u n u n t i l all of 50, w h e r e converge. T a b l e s 3 , 4, 5 s h o w s m a x | u | was t h e CFL c o n d i t i o n i n c r e a s e s 20, 25, .. . ,40 a n d d e c r e a s e s regression i t was of the values as h 2 as = 50. of the v o r t i c i t y error that II w - o> II L IT—!T N = A. + B .hRe , 1 = 3, 4, 5 , Hull where 4 3 0 -3.96, B c = -2.75, B 3 0 1 = 0.95, A,, 4 = -3.19, B-. = 0.87, A 4 5 c = 0.80 . According equal 1 t o Theorem t o 1; h o w e v e r , B^ a n d B^. T h i s 2. 1 t h e c o e f f i c i e n t s B should be there exists discrepancy i n the values of might be due t o t h e s m a l l 110 number of points = used i n f i t t i n g the s t r a i g h t l i n e i n b o t h c a s e s . The v a l u e w h i c h was a d j u s t e d w i t h a number o f By than the one used a linear Similarly, stream function in B^ and regression error shows , B^ for that points is very the the these twice closed values of larger to 1. of the are proportional the vorticity to 2 h Re, with Figures 5a and 5b show the in exact and the case. approximate r e s p e c t i v e l y of the s t r e a m f u n c t i o n i n the r u n Re = solutions, 1000, s i m i l a r d i s c r e p a n c i e s as h * = 50 At = 2.30h. F i g u r e 5c i s the p o i n t w i s e e r r o r of the s t r e a m f u n c t i o n . F i g u r e s 6a and 6b r e p r e s e n t the e x a c t and approximate solutions of the the F i g u r e 6c i s the p o i n t w i s e e r r o r in this near figure the of present the boundaries, boundaries. type that This end conditions experiments conditions vorticity deserves we recomended i n for have [4], inthe of t h e v o r t i c i t y . particularly fact vorticity errors the more the used are spline the procedure. type we have a l s o o b s e r v e d n a t u r a l end c o n d i t i o n s gave v e r y poor r e s u l t s . Ill and with second run. We o b s e r v e larger eastern testing same at and northern different In the of end that the h" 1 CFL Vort. Error Str.Func. Error 20 1.27 -1.25 - 4 . 26 25 1.60 -1.42 -4.71 30 1. 90 - 1 . 55 - 5 . 01 35 2.22 - 1 . 70 - 5 . 51 40 2.37 - 1 . 90 -5.73 50 2.30 - 2 . 10 -6.20 T a b l e 3. h" Logarithm of Re. 100 r.m.s.errors CFL Vort. 35 2. 22 -.63 - 3 . 50 40 2. 37 -.78 -3.77 II 50 2. 30 - . 90 - 4 . 20 II 1 T a b l e 4. h" 1 Error Str.Func. Logarithm of CFL Vort. Error r.m.s. Error Re. errors Str.Func. Error Re. 2. 22 -.36 - 2 . 50 40 2. 37 - . 48 - 2 . 81 II 50 2.30 -3.25 II - . 65 Logarithm of r.m.s. 112 Number 1000 35 T a b l e 5. Number errors 3000 Number F i g . 5a E x a c t Stream F u n c t i o n 113 Solution Fig. 5b Approximated Stream 114 Function 5c P o i n t w i s e E r r o r f o r Stream 115 Function Fl 9. 6a F y , „ + 116 117 118 APPENDIX The main purpose order t o do so their approximation of we n e e d some The T r a c e appendix results properties are enunciated without 1. this TO CHAPTER by IV is to prove Lemma 2.3; on t h e t r a c e o p e r a t o r s finite elements. These in y . and results proof. Theorem 2 Let n boundary £ 1), be T. a bounded We d e f i n e Lipschitz the continuous trace operators subset y ^ and y of (j R an with integer f o r a f u n c t i o n u d e f i n e d i n Q by y u = u| Q , r (A. 1) 2 _ (A.2) 1=1 j where are Theorem Al. number. Then can H S be j. • be a trace to positive defined continuous the y^., linear mapping and Magenes consequence constants c of of [20 p. this and c above j < s operators y 4. H (Q) as 0 s operator rrS/n\ linear a the u n i t normal v e c t o r on operators Moreover, Lions of domain i- continuous As SI the extended 2 Proof. the components Let ^^ (D. 3 i j onto := S-l/2 and s 1/2, j mapping \t Q> f|j n = - T. a integer, H (Q) VJ' ' ' '*p^ ..S-j-1/2 ^ such 119 that onto i s a . (F). 42] theorem real • we have that " there exist c . Hull _ _ £ lly .ull . _ _ £ c . llu II 1 s,2,Q j s-j-l/2,2,Q 2 The n o r m in H I f II r-l/2 (D r-l/2,2,Y is defined _ _ = inf , • s,2,Q (A.3) n by \lv\l _ _ . r , 2, fi (A. 4 ) r v=f 0 In we h a v e Kerj^ 2. Boundary Kerj = H (Q) the f a m i l y of be a f i n i t e finite dimensional H. c P,Cn; dimensional space such C H, , n 1+1 the approximation p r o p e r t i e s inf „ X*H llu - *|| spaces {Q^}- that . s, 2, of 2: 02 ff^ £ 0 (Chap.II! < Ch "lluII „ _ fi r , 2, fi (A. 5) 0 h a- = min( 1 + 1-s, (A.5) holds for Let //Yn; 2 Recall . 2 2 characterization Families. Consider e the f o l l o w i n g a r b i t r a r y r, r-s) s e v e n r, 120 s £ 0 (cf. [1]). Theorem A2. Let Q. be functions U e fund on T such ions a domain with property For a proof The of following boundary \yj}J Q spaces (A.5) > m - this there be given following i) a r on the a finite ( £ h w h e r e C is H^(D of Aziz [1]. properties and Schatz of [7]. o f o r d e r q^, space • H^(Q), the Let and let then the such that 0 £ 1 - q . J 1/2, 0 £ J J lly .u - y .(/II ; < J J a. be h V1/ J W ' a real J number 1/2. of such h and y . that <7j_^ 2/2 + < a < J T h e n f o r u € /Y°Y£2; lly .u - y l/ll . _ £ Ch *f ' J (A. 6 ) 2 independent 0 < u . £ 1 - q . J J , J J a constant Let V j=0 Bramble dimensional numbers j=o h to J i _ i of q .+u . J inf UeH approximation trace operators 1-1 q .+u .+A . C £ h lly ull and families traces - q . - u . J J 1-1 ii) generate the > 0 . due e u ., A . r e a l J J A . £ J + 1/2 J h Then i n e q u a l i t i e s hold: For inf above. theorem see Babuska-and results H^CD 1/2 be a s y s t e m o f n o r m a l = as that 1 + 1/2 Proof. defined 3 2 r 121 oc-1/2 Hull a ' _ _ . > 2 n (A. 7 ) + I 4 . P r o o f o f Lemma 2 . 3 To p r o v e Lemma 2 . 3 we r e c a l l s a t i s f i e s f o r q € H (Q) f| H (Q) S (Vw, V ) v | We f i r s t r - q b the R i t z , V* 6 H , Qh P^g = w (A.8.1) (A.8.2) . h e s t i m a t e the term llvYw - q)\\. two a u x i l i a r y projection , 1 £ s £ 1+1 2 = CVg, V*; X that Towards t h i s e n d , we d e f i n e problems. L e t w, e H,(Q) 1 be the f i n i t e element a p p r o x i m a t i o n h -A V q = f , 2 H u q\ = q T h of ( A . 9 . 1) • (A.9.2) Then (Vw., where w V) X = CVg, V*; (Vw , Vx) - , *X e H , n , (A.10) e H (Q) i s s u c h t h a t w^lp = q . f o L e t w^ be the f i n i t e element a p p r o x i m a t i o n -V g' 2 g'l r = 0 , = q (A.10.1) - q b h of h (A.10.2] . Then (Vw y V) X = -fVw , 4 Vx) , Vx e H 122 Qh , (A.11) w h e r e w . € H (Q) i s s u c h t h a t w. |_ = q, , - q , . 4 4Y oh o By the theory Theor. of finite elements f o re l l i p t i c 0 C depends the t r a c e on A . By v i r t u e n theorem , w Since ( c f . [26], 8.5) we h a v e liq - wJI. _ _ £ Ch "llqll . _ , cr = mind, 1 1,2,0. s, 2,Q where problems 3" i C ] > 2 i n using A. 6 with q o f the Lax-Milgram "b"Vi/2 q ^ i s thef i n i t e s-1), ,r l 2 • (A. 12) theorem and (A - 13) e l e m e n t a p p r o x i m a t i o n o f q ^ e H^Cfi), = 0 and X^. = r llvJI. , _ £ Chr~1/2\\q.\\ , _ i C_h .3 2,2, fi 2 ^b r+2/2,2, T 2 S_1 llqll , _ . (A. 14) s, 2, fi Now, IIVCW -q )\\2 = cvcw - q;, vcw - q;; = = CVCw - q;, ~V(V - w ) + V(w - q)) 1 £ IIVCw - q)\\ IIVCw - v ;il + HVfw - q)\\ IIVCw^- qjll . Hence Wv - q;il * II fv - V"2,2,fi = " 3 i,2,n w From then (A.12) a n d (A.14), U + i ~ + nw llw^ - qH *h.2.a thef o l l o w i n g estimate IIVCw - q;il £ Ch S 2 llqll s 2 fl 123 • lf2>Q holds (A. 15) It remains to estimate the defining an auxiliary e l l i p t i c term llv - qll. We proceed by problem with homogeneous boundary 2 2 conditions. Let \p e H (U) and tp e L (Q). Furthemore, -A V \I) n = <p in f2 , 2 U (A. 16. 1) 0 1 = 0 . (A.16.2) By regularity theory \\if>\l 22 n s Cllf>ll . (A. 18) Then Cw - q, <p) = -(w - q,V ip) = CVCw - q), Vip) - <-|jj-p 2 For any x e 7 (^ ~ Q HQh(Q) (w - q, <p) = CVCw - q;, V C 0 - - < -| L n ' - q;> p (A.19) By virtue of the following items: i) Approximation properties of H ^ - IVC0 ii) s K hll^ll 2 2 n = /C hll<pll 2 The trace theorem l l I/2,2 r V" 2 2,Q V" l S ( iii) 2 (A. 5) l S ) Inequality A.6 and the trace theorem 124 where K_ are constants, the inequality (A.19) becomes (w - q, <p) £ K h\\V(w 6 - q)\\ llipll +K h \\q\\ S 7 s 2 Q \\<p\\ . Taking <p = w - q and using (A. 15) the Lemma is proven. 125 • CHAPTER V CONCLUSIONS We have formulated and analysed an algorithm of Galerkin-Characteristic convection-dominated algorithm method, has formulation vorticity diffusion been similar combined to of methods the the one to problems. with a given in function integrate The mixed to equations equations) proposed finite [11], quasi-geostrophic — stream the type for element produce a (potential a mid—latitude baroclinic ocean. The integration our algorithm is of the transport-diffusion essentially equations a two stage process. by The f i r s t stage corresponds to the integration of the advection operator along with the characteristic curves Galerkin method. reinterpretation of - in terms yields of For this the usual particle of the stage, a computationally efficient we combination show with scheme, functional of the dependent variable at the departure The scheme the grid maximum norm for smooth functions at the foot of this point thus devised unconditionally stable and convergent. the grids, is values of of a points conservative, Our error analysis in stage proves that for sufficiently the approximate solution is the a which consists by the particles. that rectangular interpolating of splines in Galerkin-Characteristic method methodology, cubic flow characteristic curves. super convergent To assess the performance of our algorithm for the hyperbolic stage we have 126 carried out two types of advective experiments. The f i r s t one is a f a i r l y hard problem which consists of advecting a cone in a rotating flow f i e l d . Munz [28] this by problem obtained high has reported the results of resolution finite schemes of the type total variation diminishing difference (TVD) such as MUSCL, UNO, etc. , which are considered to be the best finite difference visual schemes for hyperbolic problems. A comparison of our results with those portrayed in figures 6 to 16 of [28] shows that: i ) Our scheme is able to keep the shape of the cone better than any high resolution scheme. ii) The 'numerical diffusion' of our scheme is lower than that of the high resolution f i n i t e difference schemes, with the possible exception of the superbee scheme. iii) Our scheme exhibits small wiggle activity at the base of the cone, whereas the high resolution f i n i t e difference schemes are wiggle free because they possess the TVD property. iv) Our scheme is able to use much larger time steps than any high resolution f i n i t e ours is unconditionally difference scheme because stable with respect to the 2 L —norm. This property is responsible for keeping the wiggles under s t r i c t control. Our second 'slotted' hard numerical cylinder and the consists in a rotating flow f i e l d . problem because surface experiment of the discontinuities 'slotted' region, 127 of advecting This on is the respectively. A a a very lateral visual comparison of our results with those reported in SHASTA and figures 6 modified SHASTA schemes, to 16, using the and [28], aforementioned [37], using specifically high resolution schemes, yields the same conclusion as in the cone experiment. The second stage corresponds the algorithm via to the time progression of the diffusion mechanism starting with the output of the previous stage. out in our analysis The time progression by the Crank-Nicholson scheme. is carried The error 2 analysis with respect to L -norm, based on developed in [10] and [35], reveals for At = 0(h) of a term of order 0(h) techniques the presence in the evolutionary component of the error. This term is due to the particle approximation of the first stage and is suspected Galerkin-Characteristic methods to be which inherited by those approximate the inner products by alternative techniques, such as the ones proposed in [27] and [5]. In this sense, the Crank—Nicholson scheme is suboptimal when used in conjunction with our method; however, for larger time steps it leads to more accurate results than the backward Euler scheme , which appears to be optimal for At = 0(h). Our error analysis ones used in [13], technique is more general than the [31] and [33], which cannot be used for the Crank—Nicholson scheme. The constants C Q1 and C in Q2 Chap.IV (2.11), which multiply the terms of the evolutionary error, are the product of the exponential constant of the Gronwall's inequality with the norm of characteristic the derivatives curves. For of the convection 128 variable along the dominated flows the speed of variation along the characteristics is less rapid than in the t direction, so the algorithm permits larger time steps without loss of accuracy. On the other hand, since the algorithm in the second stage is also unconditionally stable, then the exponential constant of the Gronwall's inequality is less than one; so that as time progresses the influence of the evolutionary decreases, error and approximate on the eventually solution is global for T accuracy of the — > oo the' error the approximation error. method of The the latter 2 type of error is 0(h ) with free-slip boundary conditions, and 0(h ) with no—slip boundary conditions. one order than higher literature. with respect Finally, 2 previous we prove to L -norm for These estimates estimates that the given approximation the stream functions is 2 which is the optimal, while the L —norm approximation velocities is 0(h). Numerical in experiments the error 2 0(h ), error with are for no-slip boundary conditions shows the validity of our error estimates. As a f i n a l remark, we should mention that our analysis can be extended to approximate the solution equations by the proposed algorithm. 129 of Navier—Stokes REFERENCES 1. Adams, R.A.: Sobolev Spaces. Academic Press, New York (1975) 2. Babuska, I. and Aziz, A.K.: "Survey Lectures on the Mathematical Foundations of the Finite Element Method," Edi., The Mathematical Foundations Aziz, A.K., Finite Element Method with Applications Differential Equations. of to Partial Academeic Press, N.Y. (1972) in the pp 5-359 3. Bardos, C. , Bercovier, M. , and Pironneau, 0. : The vortex method with f i n i t e elements. Math. Comp. 36, 119-136 (1981) 4. Behforooz, cubic G. H. and Papamichael, splineinterpolation. J. N. : Inst. End Maths. conditions for Applies., 23, 355-366. (1979) 5. Benque, A J. P. , Finite Element -Proceedings elements Ibler, of in B. , Keramsi, Method for A. , and Labadie, Navier-Stokes G. : Equations. the third international conference on f i n i t e flow problems. Banff, Alberta, Canada, 10-13 June(1980) 6. Bermejo, R. : "An Element Method," analysis a Quasi-Geostrophic in Graham, G.A.C. Continuum Mechanics Co. New York of and Malik, and its Applications. S. L. Finite Eds., Hemisphere Publ. (1989) pp 825-835 7. de Boor, C. : A Practical Guide to New York, Berlin, Heidelberg (1978) 130 Splines.Springer-Verlag. 8. Boris, J.P., and Book, D.L.: Flux corrected transport. SHASTA. J. Comp. Phys. 1 1 , 36-39 (1973) 9. Bramble, J.H. and Schatz, A.: Least squares method for 2mth 10. order elliptic 1-32. (1972) Brenner, boundary—value P. , Crouzeix, M. problems. Math. and Thomee, V. : Comp., Single 25, step methods for inhomogeneous linear differential equations Banach spaces. RAIRO, Anal. Numer., 11. Bristeau, M.O., methods for Glowinski, R. and Periaux, the Navier—Stokes simulation of compressible 1 6 , 5-26 equation. and in (1982) J. : Numerical Applications incompressible to the viscous flow. Comp. Phys. Rep., 6 , 73-183 (1987) 12. 13. Ciarlet, The P.: Finite Element Method for Problems. Amsterdam: North-Holland (1978) Douglas, J. , and Russell, T. F. : Elliptic Numerical methods for convection dominated diffusion problems based on combining the method of the characteristics with finite elements or f i n i t e differences. SIAM J. Numer. Anal. 1 9 . (1982) 14. Fix, G. , Gunzburger, "Mixed Finite Equations. , M., Nicolaides, R. Element Proc. 5th Approximation International and Peterson, J . : for Symposium Elements and Flow Problems, Austin, Texas, 15. Girault, V. and Raviart, P.-A.: Finite of the Navier—Stokes vol. Equations. Biharmonic on finite (1984) Element Approximation Lecture Notes in Mathematics, 749. Springer Verlag. Berlin (1979) 131 the 16. Glowinski, R. and Pironneau, 0. : Numerical methods for the first biharmonic equation and Stokes problem. SIAM. Rev., 17. Grisvard, P.: Boundary for the two—dimensional 17, 167-212 (1979) Value Problems in Non-Smooth Univ. of Maryland, Dept. of Math. Lecture Notes, Domains. 19 (1980) 18. Hale, J. K. : Ordinary Differential Equations. Robert E. Krieger Publish Co. Malabar 19. Hasbani, Y., Livne, E., and characteristics 2 n d Edition. Florida(1980) and Bercovier, M.: Finite elements applied to advection-diffusion equations. Comput. Fluids. 11, 71-83 (1983) 20. Holland, W.R. and Rhines, P.B.: An example of eddy-induced ocean circulation. J. Phys. Oceang. , 10, 21. Hughes, T. "Streamline J. Tezduyar, Upwind Navier-Stokes Fourth R., and Internat. Formulation First Conf. T. E. , for Order on Finite 1010-1031 (1980) and Brooks, A.: Advection-Diffusion, Hyperbolic Equations". Element Methods in Fluids, Tokyo (1982) 22. Johnson, element C. , Navaert, methods for U. , and linear Pitkaranta, hyperbolic J. : problems. Finite Comp. Meth. Appl. Mech. Engrg. 45, 285-312 (1984) 23. Lions, J.L. Non-Homogeneous and at Magenes, Applications, (1968) J 132 E. : Problems vols aux 1 and 2, Limites Dunod. Paris 24. Mas-Gallic, first S. , order and Raviart, symmetric P-A.: systems. A particle method for Numer. Sobolev Spaces. (1987)25. Maz'ja, V. G. : Math. 51,323-352 Springer Verlag. New York, Berlin, Heidelberg (1979) 26. McWilliams, models J.C.: A note on consistent in multipled connected domains. quasi-geostrophic Dyn. Atmos. Ocean. 1, 427-441 (1977) 27. Morton, K. the W., Priestley, A., and Suli, E.: S t a b i l i t y of Lagrange-Galerkin method with non-exact integration. RAIRO. MAN. 4, 225-250 (1988) 2 28. Munz, C-D.: On the numerical resolution schemes for J. 29. dissipation of high the hyperbolic conservation laws. Comp. Phys, 77. 18-39 (1988) Oden, J.T. and Mathematical Reddy, Theory An Introduction J.N.: of Finite Elements. to the Wiley-Interscience, New York (1976) 30. Pedlosky, J . : Geophysical Fluid Dynamics. Springer Verlag, New York (1979) 31. Pironneau, its 0. : On the transport applications to diffusion algorithm and the Navier-Stokes equations. Numer. Math. 38, 309-332 (1982) 32. Raviart, Numerical P-A.: An Methods analysis in Fluid of Dynamics. Lecture Notes in Mathematics, vol. Berlin (1985) 133 particle 1127. F. methods. Brezzi. Springer In Edi. Verlag, 33. Suli, E. : Convergence and nonlinear s t a b i l i t y of Lagrange-Galerkin method for the Navier-Stokes the equations. Numer. Math. 53, 459-483 (1988) 34. Temperton, C. , and Staniforth, A. : level semi-Lagrangian R. Meteorl. Soc. 35. Thomee, V.: An efficient two-time semi-implicit integration scheme. Q. J. 113, 1025-1039 (1987) Galerkin Finite Element Methods for Parabolic Problems. Lecture Notes in Mathematics, vol. 1054. Springer Verlag, Berlin (1984) 36. Widlund, 0. : piecewise On best polynomial error bounds for functions. approximation by Numer. Math. 27, 327-338 (1977) 37. Zalesak, S. T. : Fully multidimensional flux-corrected transport. J. Comp. Phys. 31, 335-362 (1979) 38. Zeniseck, in the A.: finite Discrete forms of Friedrichs'inequa;lities element method. 255-286 (1982) 134 RAIRO. Anal. Numer. 15, APPENDIX ON INEQUALITIES We collect here some inequalities which have been used throughout the thesis,as wellas the Lax-Milgram theorem. 1. Gronwall Inequalities If <f>, a are real valued and continuous ^ 0 is integrable for a s t £ b, fi(t) on [a, b] and (1) a then <p(t) £ oc(t) + f !B(s)a(s)(exp(\ p(u)du))ds , a £ t £ b . (2) For a proof of (2) see [18]. The discrete version of the above inequality goes as follows. Let ( )' ^ ^' ^ r? a C n n b e ^ r e e sequences numbers such that (c ) is monotonically n a +b £ n n n-1 n n * of positive increasing 1, A > 0 , real and (3) then a n + b £ c exp(Xn) , n ^ 0 . n n For a proof ao (4) see reference [15] Chap.V. 135 (4) 2. Schwarz Inequality 2 2 Let n c R and f.geL (Q), then f fgdfi s f f d£> f g d n . 2 J r> (5) 2 Jn Jn 3. Young Inequality Let e > 0, and a, b e R, then ab £ -:—a 4c + eb . (6) 4. Friedrichs Inequality If Q is connected and bounded at least then for each integer m £ 0, there exists in one direct ion, a constant K = K(m, Q) > 0, such that Ivll _ s K\v\ _ , Vv e fljjfn; . m (7) The proof of (7) can be found in [1]. A discrete version of this inequality is given in [38]. 5. Lax-Milgram Theorem Let a(u,v) be a continuous and elliptic V, i.e., \a(u,v)\ £ tfllull llvll , Vu, v e V \a(u,u)\ £ allull ^ , Vu e V. 2 136 form on the space Then , the problem a(u,v) has a uniqe Proof. See = <l,v> solution , Vv in e V and V. [2]. 137 1 e V
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Analysis of a Galerkin-Characteristic algorithm for...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Analysis of a Galerkin-Characteristic algorithm for the potential vorticity-stream function equations Bermejo, Rodolfo 1990
pdf
Page Metadata
Item Metadata
Title | Analysis of a Galerkin-Characteristic algorithm for the potential vorticity-stream function equations |
Creator |
Bermejo, Rodolfo |
Publisher | University of British Columbia |
Date Issued | 1990 |
Description | In this thesis we develop and analyze a Galerkin-Characteristic method to integrate the potential vorticity equations of a baroclinic ocean. The method proposed is a two stage inductive algorithm. In the first stage the material derivative of the potential vorticity is approximated by combining Galerkin-Characteristic and Particle methods. This yield a computationally efficient algorithm for this stage. Such an algorithm consists of updating the dependent variable at the grid points by cubic spline interpolation at the foot of the characteristic curves of the advective component of the equations. The algorithm is unconditionally stable and conservative for Δt = O(h). The error analysis with respect to L² -norm shows that the algorithm converges with order O(h); however, in the maximum norm it is proved that for sufficiently smooth functions the foot of the characteristic curves are superconvergent points of order O(h⁴ /Δt). The second stage of the algorithm is a projection of the Lagrangian representation of the flow onto the Cartesian space-time Eularian representation coordinated with Crank-Nicholson Finite Elements. The error analysis for this stage with respect to L²-norm shows that the approximation component of the global error is O(h²) for the free-slip boundary condition, and O(h) for the no-slip boundary condition. These estimates represent an improvement with respect to other estimates for the vorticity previously reported in the literature. The evolutionary component of the global error is equal to K(Δt² + h), where K is a constant that depends on the derivatives of the advective quantity along the Characteristic. Since the potential vorticity is a quasi-conservative quantitiy, one can conclude that K is in general small. Numerical experiments illustrate our theoretical results for both stages of the method. |
Subject |
Galerkin methods Characteristic functions Lagrange equations |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2011-01-10 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0080415 |
URI | http://hdl.handle.net/2429/30561 |
Degree |
Doctor of Philosophy - PhD |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-UBC_1990_A1 B47.pdf [ 5.49MB ]
- Metadata
- JSON: 831-1.0080415.json
- JSON-LD: 831-1.0080415-ld.json
- RDF/XML (Pretty): 831-1.0080415-rdf.xml
- RDF/JSON: 831-1.0080415-rdf.json
- Turtle: 831-1.0080415-turtle.txt
- N-Triples: 831-1.0080415-rdf-ntriples.txt
- Original Record: 831-1.0080415-source.json
- Full Text
- 831-1.0080415-fulltext.txt
- Citation
- 831-1.0080415.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0080415/manifest