Unformatted text preview: M.Sc. in Computational Science Fundamentals of Atmospheric Modelling
´ Peter Lynch, Met Eireann
Mathematical Computation Laboratory (Opp. Room 30) Dept. of Maths. Physics, UCD, Belﬁeld. January–April, 2004. Lecture 12 SemiLagrangian Advection 2 Introduction
• We have studied the Eulerian leapfrog scheme and found it to be conditionally stable. • The criterion for stability was the CFL condition The Basic Idea
The semiLagrangian scheme for advection is based on the idea of approximating the Lagrangian time derivative. It is so formulated that the numerical domain of dependence always includes the physical domain of dependence. This necessary condition for stability is satisﬁed automatically by the scheme. In a fully Lagrangian scheme, the trajectories of actual physical parcels of ﬂuid would be followed throughout the motion. The problem with this aproach, is that the distribution of representative parcels rapidly becomes highly nonuniform. In the semiLagrangian scheme the individual parcels are followed only for a single timestep. After each step, we revert to a uniform grid. c ∆t ≤ 1. ∆x
• For high spatial resolution (small ∆x) this severly limits the maximum time step ∆t that is allowed. • In numerical weather prediction (NWP), timeliness of the forecast is of the essence. • In this lecture, we study an alternative approach to time integration, which is unconditionally stable and so, free from restrictions of the CFL condition. The semiLagrangian algorithm has enabled us to integrate the primitive equations using a time step of 15 minutes. This can be compared to a typical timestep of 2.5 minutes for conventional schemes. The consequential saving of computation time means that the operational numerical guidance is available to the forecasters much earlier than would otherwise be the case. The semiLagrangian method was pioneered by Andr´ Robert, e the renowned Canadian meteorologist. The ﬁrst operational implementation of such a scheme was in 1982 at the Irish Meteorological Service. SemiLagrangian advection schemes are now in widespread use in all the main Numerical Weather Prediction centres. Eulerian and Lagrangian Approach
We consider the linear advection equation which describes the conservation of a quantity Y (x, t) following the motion of a ﬂuid ﬂow in one space dimension with constant advecting velocity c. This may be written in either of two alternative forms: ∂Y ∂Y +c =0 ⇐ Eulerian Form ∂t ∂x dY =0 ⇐ Lagrangian Form dt The general solution is Y = Y (x − ct). To develop numerical solution methods, we may start from either the Eulerian or the Lagrangian form of the equation. For the semiLagrangian scheme, we choose the latter. 5 6 Since the advection equation is linear, we can construct a general solution from Fourier components Y = a exp[ik (x − ct)] ; k = 2π/L . This expression may be separated into the product of a function of space and a function of time: Y = a × exp(−iωt) × exp(ikx) ; ω = kc . Therefore, in analysing the properties of numerical schemes, we seek a solution of the form
n Ym = a × exp(−iωn∆t) × exp(ikm∆x) = aAnexp(ikm∆x) Numerical Domain of Dependence. where A = exp(−iω ∆t). The character of the solution depends on the modulus of A: If A < 1, the solution decays with time. If A = 1, the solution is neutral with time. If A > 1, the solution grows with time. In the third case (growing solution), the scheme is unstable. Space axis horizontal Time axis vertical +++O+++ n    *******       *************    ++*******************++ n1   *************************     *******************************   +*************************************+ n2  *******************************************   *************************************************  ******************************************************* n3        m3 m2 m1 m m+1 m+2 m+3
n For the Eulerian Leapfrom Scheme, the value Ym at time n∆t and position m∆x depends on values within the area depicted by asterisks. Values outside this region have no n inﬂuence on Ym. Numerical Domain of Dependence
n Each computed value Ym depends on previously computed values and on the initial conditions. The set of points which n inﬂuence the value Ym is called the numerical domain of n dependence of Ym. Parcel coming from Outside Domain of Dependence
+++++•+ n      • *******      • *************  ++++•******************* n1    • **********************   •  ************************* ++•+**************************** n2  •  *******************************  •   ********************************** •+************************************* n3        m5 m4 m3 m2 m1 m m+1 The line of bullets (•) represents a parcel trajectory.
n The value everywhere on the trajectory is Ym. Since the parcel originates outside the numerical domain of dependence, the Eulerian scheme cannot model it correctly.
10 It is clear on physical grounds that if the parcel of ﬂuid arriving at point m∆x at time n∆t originates outside the numerical domain of dependence, the numerical scheme cannot yield an accurate result: the necessary information is not available to the scheme. Worse again, the numerical solution may bear absolutely no relationship to the physical solution and may grow exponentially with time even when the true solution is bounded. A necessary condition for avoidance of this phenomenon is that the numerical domain of dependence should include the physical trajectory. This condition is fulﬁlled by the semiLagrangian scheme.
9 The central idea of the Lagrangian scheme is to represent the physical trajectory of the ﬂuid parcel. We consider a parcel arriving at gridpoint m∆x at the new time (n + 1)∆t and ask whence it has come. The departure point will not normally be a grid point. Therefore, the value at the departure point must be calculated by interpolation from surrounding points. But this interpolation ensures that the trajectory falls within the numerical domain of dependence. We will show that this leads to a numerically stable scheme. Interpolation using Surrounding Points
+++++◦+ n      ◦       ◦   ++++++•++++++++ n1    ◦      ◦     +++++++◦+++++++ n2        m5 m4 m3 m2 m1 m m+1 The line of circles (◦) represents a parcel trajectory. At time (n − 1)∆t the parcel is at (•), which is not a gridpoint. The value at the departure point is obtained by interpolation from surrounding points. Thus we ensure that the physical trajectory is within the domain of numerical dependence. The advection equation in Lagrangian form may be written dY = 0. dt In physical terms, this equation says that the value of Y is constant for a ﬂuid parcel. Applying the equation over the time interval [n∆t, (n + 1)∆t], we get Value of Y at Value of Y at point m∆x at = departure point time (n + 1)∆t at time n∆t or, in symbolic form,
n Ym+1 = Y•n Interpolation using Surrounding Points
+++++◦+ n+1      ◦       ◦   ++++++•++++++++ n m5 m4 m3 m2 m1 m m+1 The distance travelled in time ∆t is s = c∆t. We deﬁne the integer and fractional parts of s as follows p = [s] = Integral part of s α = s − p = Fractional part of s Note that, by deﬁnition, 0 ≤ α < 1. So, the departure point falls between the grid points m − p − 1 and m − p. In the picture, p = 1 and α ≈ 2/3. A linear interpolation gives
n n Y•n = αYm−p−1 + (1 − α)Ym−p . where Y•n represents the value at the departure point, normally not a grid point. Check: Verify that this gives sensible values in the limits α = 0 and α = 1.
13 14 Numerical Stability of the Scheme
The discrete equation may be written
n n n Ym+1 = αYm−p−1 + (1 − α)Ym−p . Again, A = α exp[ik (−p − 1)∆x] + (1 − α) exp[ik (−p)∆x] = exp(−ikp∆x) · [(1 − α) + α exp(−ik ∆x)] Now consider the squared modulus of A: A2 = = = = = exp(−ikp∆x)2 · (1 − α) + α exp(−ik ∆x)2 (1 − α) + α cos k ∆x − iα sin k ∆x2 [(1 − α) + α cos k ∆x]2 + [sin k ∆x]2 (1 − α)2 + 2(1 − α)α cos k ∆x + α2 cos2 k ∆x + α2 sin2 k ∆x 1 − 2α(1 − α)[1 − cos k ∆x] . Let us look for a solution of the form
n Ym = a An exp(ikm∆x) . Substituting into the equation we get aAn+1 exp(ikm∆x) = αaAn exp[ik (m − p − 1)∆x] + (1 − α)aAn exp[ik (m − p)∆x] or, removing common terms, A = α exp[ik (−p − 1)∆x] + (1 − α) exp[ik (−p)∆x] We note that 0 ≤ (1 − cos θ) ≤ 2. Taking the largest value of 1 − cos k ∆x gives A2 = 1 − 4α(1 − α) = (1 − 2α)2 < 1 . Taking the smallest value of 1 − cos k ∆x gives A2 = 1 . In either case, A2 ≤ 1, so there is numerical stability. Discussion and Conclusion
We have determined the departure point by interpolation. This ensures that 0 ≤ α < 1. This in turn ensures that A ≤ 1. In other words, we have unconditional numerical stability. The implication is that the time step is unlimited. In contradistinction to the Eulerian scheme there is no CFL criterion. Of course, we must consider accuracy as well as stability The time step ∆t is chosen to ensure suﬃcient accuracy, but can be much larger than for an Eulerian scheme. Typically, ∆t is about six times larger for a semiLagrangian scheme than for an Eulerian scheme. This is a substantial gain in computational eﬃciency.
17 18 ...
View
Full Document
 Spring '07
 Osserman
 Math, Numerical Analysis, Eulerian Leapfrom Scheme, m3 m2 m1, numerical domain

Click to edit the document details