Lect-12-P4 - M.Sc. in Computational Science Fundamentals of...

Info iconThis preview shows page 1. Sign up to view the full content.

View Full Document Right Arrow Icon
This is the end of the preview. Sign up to access the rest of the document.

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, Belfield. January–April, 2004. Lecture 12 Semi-Lagrangian 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 semi-Lagrangian 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 satisfied automatically by the scheme. In a fully Lagrangian scheme, the trajectories of actual physical parcels of fluid would be followed throughout the motion. The problem with this aproach, is that the distribution of representative parcels rapidly becomes highly non-uniform. In the semi-Lagrangian scheme the individual parcels are followed only for a single time-step. 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 semi-Lagrangian 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 semi-Lagrangian method was pioneered by Andr´ Robert, e the renowned Canadian meteorologist. The first operational implementation of such a scheme was in 1982 at the Irish Meteorological Service. Semi-Lagrangian 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 fluid flow 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 semi-Lagrangian 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 | | | ******* | | | | | | ************* | | | +--------+--------*******************--------+--------+ n-1 | | ************************* | | | | ******************************* | | +--------*************************************--------+ n-2 | ******************************************* | | ************************************************* | ******************************************************* n-3 | | | | | | | m-3 m-2 m-1 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 influence 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 influence the value Ym is called the numerical domain of n dependence of Ym. Parcel coming from Outside Domain of Dependence +--------+--------+--------+--------+--------•--------+ n | | | | | • ******* | | | | | •| ************* | +--------+--------+--------+--•-----******************* n-1 | | | •| ********************** | | |• | ************************* +--------+-----•--+--------**************************** n-2 | |• | ******************************* | • | | ********************************** •--------+--------************************************* n-3 | | | | | | | m-5 m-4 m-3 m-2 m-1 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 fluid 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 fulfilled by the semi-Lagrangian scheme. 9 The central idea of the Lagrangian scheme is to represent the physical trajectory of the fluid 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 | | | | | ◦ | | | | | | ◦| | | +--------+--------+--------+++•++++++--------+--------+ n-1 | | | ◦| | | | | | |◦ | | | | +--------++++++◦+++--------+--------+--------+--------+ n-2 | | | | | | | m-5 m-4 m-3 m-2 m-1 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 fluid 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 m-5 m-4 m-3 m-2 m-1 m m+1 The distance travelled in time ∆t is s = c∆t. We define the integer and fractional parts of s as follows p = [s] = Integral part of s α = s − p = Fractional part of s Note that, by definition, 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: |A|2 = = = = = |exp(−ikp∆x)|2 · |(1 − α) + α exp(−ik ∆x)|2 |(1 − α) + α cos k ∆x − iα sin k ∆x|2 [(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 |A|2 = 1 − 4α(1 − α) = (1 − 2α)2 < 1 . Taking the smallest value of 1 − cos k ∆x gives |A|2 = 1 . In either case, |A|2 ≤ 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 sufficient accuracy, but can be much larger than for an Eulerian scheme. Typically, ∆t is about six times larger for a semi-Lagrangian scheme than for an Eulerian scheme. This is a substantial gain in computational efficiency. 17 18 ...
View Full Document

This note was uploaded on 01/31/2011 for the course MATH 21A taught by Professor Osserman during the Spring '07 term at UC Davis.

Ask a homework question - tutors are online