Lect-11-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 11 Time Integration Methods 2 Introduction • We have studied various simple solutions of the shallow water equations by making approximations. • In particular, by means of the perturbation method the equations have been linearised, making them amenable to analytical investigation. • However, to obtain solutions in the general case, it is necessary to solve the full nonlinear system. • In numerical weather prediction (NWP) the fully nonlinear primitive equations are solved by numerical means. • In the atmosphere, the nonlinear advection process is a dominant factor. • To get some idea of the methods used, we look at the simple problem of formulating time-integration algorithms for the solution of the simple advection equation. • The development of accurate and efficient numerical schemes is a huge area of ongoing research, and we can do no more than to introduce the subject. ´ • Much of our research effort in Met Eireann over the past ten years has been in this area, and a brief outline of this work will be given. Discretization Methods There are several distinct approaches to the formulation of computer methods for solving differential equations. We will confine ourselves to the finite difference methods. Other approaches include finite element method and the spectral method. The central idea of the finite difference approach is to approximate the derivatives in the equation by differences between adjacent points in space or time, and thereby reduce the differential equation to a difference equation. • An analytical problem becomes an algebraic one. • A problem with an infinite degree of freedom is replaced by one with a finite degree of freedom. • A continuous problem goes over to a discrete one. 5 The Finite Difference Method We start by looking at the Taylor expansion of f (x): f (x + ∆x) = f (x) + f (x).∆x + 1 f (x)∆x2 + [O(∆x3)] 2 1 f (x − ∆x) = f (x) − f (x).∆x + 2 f (x)∆x2 + [O(∆x3)] (1) (2) The higher order terms, represented by O(∆x3), become less important as ∆x becomes smaller. We can neglect these and use (1) or (2) to get an approximation for the derivative of f (x) as follows: f (x + ∆x ) − f (x ) + O(∆x) = fF + O(∆x) ∆x f ( x ) − f ( x − ∆x ) f (x ) = + O(∆x) = fB + O(∆x) . ∆x These are called the forward and backward differences. Keeping only leading terms, we incur errors of order O(∆x). f (x ) = 6 We can do better than this: subtracting (2) from (1) yields the following: f ( x + ∆x ) − f ( x − ∆x ) + O(∆x2) = fC + O(∆x2) f (x ) = 2∆x which is seen to be of order O(∆x2), therefore more accurate for small ∆x. Adding (1) and (2) gives the corresponding expression for the second derivative: f (x + ∆x) − 2f (x) + f (x − ∆x) + O(∆x2) f (x ) = ∆x 2 These centered differences are of accuracy O(∆x2). We can continue taking more and more terms, but obviously there is a trade-off between accuracy and efficiency. Fourth-order accurate schemes are sometimes used in NWP, but second order accuracy is more popular. Exercise: Consider the function f (x) = A sin(2πx/L). We know that the derivative is (2π/L)A cos(2πx/L). • Show that a forward difference approximation gives sin(π ∆x/L) fF (x) = (2π/L)A cos[2π (x + ∆x/2)/L] · π ∆x/L whereas the centered difference yields: sin(2π ∆x/L) fC (x) = (2π/L)A cos[2πx/L] · 2π ∆x/L • Compare these to the true derivative f (x) and investigate their behaviour for small ∆x. • Demonstrate thus that the centered difference is of higher order accuracy. Grid Resolution and Accuracy The size of the gridstep ∆x determines the accuracy of the numerical scheme. For the simple sine function the error depended on the ratio (∆x/L). For synoptic scale waves in the atmosphere a typical value of L is 1000 km. To make the ratio equal to 0.1 we need to have a grid size of about 100 km. This is larger than the typical gridsizes used in operational NWP. The higher the resolution, that is, the smaller the grid-size, the heavier the computational burden. There is a trade-off between resolution and accuracy. Linear Computational Instability We consider the equation describing the conservation of a quantity Y (x, t) following the motion of a fluid flow in one space dimension: dY ∂Y ∂Y = 0. ≡ +u dt ∂t ∂x If the velocity is taken to be constant, u = c, or if we linearise ¯ about a mean flow u = c, the equation becomes ∂Y ∂Y +c = 0. ∂t ∂x This simple equation we will call the linear advection equation. It is analogous to a factor of the usual wave equation ∂2 ∂2 − c2 2 ∂t2 ∂x Y= ∂ ∂ +c ∂t ∂x ∂ ∂ −c Y = 0, ∂t ∂x and the general solution is Y = Y (x − ct). 9 10 Since the advection equation is linear, we can construct a general solution from Fourier components Y = a exp[ik (x − ct)] ; k = 2π/L . We take the following initial condition for Y : Y (x, 0) = a exp[ikx] Next, we approximate the differential equation by a finite difference equation using centered differences for the space and time derivatives. Let the variables x and t be represented by the horizontal and vertical axes. Positive time corresponds to the upper half plane. The initial data occur on the x-axis. The continuous variables are replaced by discrete gridpoints at their integral values and the problem is solved on a finite difference grid. Space-Time Grid: n=5 n=4 n=3 n=2 n=1 n=0 Space axis horizontal Time axis vertical +——–+——–+——–+——–+——–+——–+ I I I I I I I I I I I I I I +——–+——–+——–+——–+——–+——–+ I I I I I I I I I I I I I I +——–+——–+——–+——–+——–+——–+ I I I I I I I I I I I I I I +——–+——–+——–+——–+——–+——–+ I I I I I I I I I I I I I I +——–+——–+——–+——–+——–+——–+ I I I I I I I I I I I I I I +——–+——–+——–+——–+——–+——–+ m=-3 m=-2 m=-1 m=0 m=1 m=2 m=3 We denote the value of Y at a grid point by: n Y (m∆x, n∆t) = Ym . Inter-dependency of Points n+1 n n-1 + X + m-1 X O X m + X + m+1 = 0. Then the finite difference approximation to the differential equation may be written as follows: n n Ym+1 − Ym−1 2∆t n n Ym+1 − Ym−1 +c 2∆x Solving for the value at time (n + 1)∆t gives n n Ym+1 = Ym−1 − c ∆t n n (Ym+1 − Ym−1) ∆x The value at the time (n + 1)∆t is obtained by adding a term to the value at (n − 1)∆t; the method is known as the leapfrog method because of this jump over the time n∆t. Note the ratio c ∆t λ≡ . ∆x The value of this will be found to be crucial. 13 The evaluation of the equation at point O involves values of n the variable at points X. Solving for Ym+1 thus requires n Ym−1 , n Ym−1 and n Ym+1 . Note how the leapfrog scheme splits the grid into two independent subgrids. 14 Solving for the value at time (n + 1)∆t gives c ∆t n n n n Ym+1 = Ym−1 − (Ym+1 − Ym−1) ∆x Assuming that we know the solution up to time n∆t, the values at time (n + 1)∆t can be calculated, and the solution advanced one timestep in this way. Then the whole procedure can be repeated to advance the solution to (n + 2)∆t, and so on. Under what conditions does the solution of the finite difference equation approximate that of the original differential equation? Intuitively, we would expect that a good approximation would be obtained provided the grid steps ∆x and ∆t are small enough. However, it turns out that this is not enough, and that the value of the ratio λ = c∆t/∆x is critical. This surprising result has important practical implications. The CFL Stability Criterion Let us assume a solution of the finite difference equation in the form n Ym = a An exp[ikm∆x] Substituting this in the equation and dividing by a common factor gives A2 + (2iλ sin k ∆x)A − 1 = 0 where λ = c∆t/∆x. This is a quadratic equation for the amplitude A, with solutions A± = −iλ sin k ∆x ± 1 − λ2 sin2 k ∆x . Question: The roots of a quadratic equation may be either real or complex, depending on the value of the coefficients. We consider in turn the two possible cases. Case I: |λ| ≤ 1 The quantity under the square-root sign is positive, so the modulus of A is given by |A|2 = 1 − λ2 sin2 k ∆x + (λ sin k ∆x)2 = 1 . The modulus of A is seen to be unity. Thus, we may write A = exp(iψ ) . The two values of the phase are ψ1 = − arcsin(λ sin k ∆x) and ψ2 = π − ψ1 . For small λ, we have ψ1 ≈ −λk ∆x = −kc∆t and ψ2 ≈ π + λk ∆x = π + kc∆t The solution of the equation may now be written n Ym = D exp(iψ1n) + E exp[i(−ψ1 + π )n] exp(ikm∆x) Taking account of the initial conditions, we get the solution n Ym = (a − E ) exp[ik (m∆x + ψ1n/k )] + (−1)nE exp[ik (m∆x − ψ1n/k )] Physical Mode Computational Mode Exercise: Check the algebra leading to this solution. The first term of this solution is called the physical mode and corresponds to the solution of the differential equation. The second term is called the computational mode. It arises through the use of centered differences resulting in the approximation of a first order differential equation by a second order difference equation (with an extra solution). If the ratio λ is small, the finite difference solution is given approximately by Y ≈ a exp[ik (m∆x − cn∆t)] which is just the analytical solution. Thus, we would expect good results from the finite difference approximation in this case. 17 18 where D and E are arbitrary constants. The centered difference approximation cannot be used for the first timestep because it would involve values at time t = −∆t which are not known. Instead, we must use another approximation for the first step, typically an uncentered forward time step. Thereafter, the leapfrog scheme can be used. In essence, this amounts to specifing another “initial contition”, the computational initial condition, at t = ∆t. The value of this determines the amplitude of the computational mode. It should be chosen to minimize this. The above solution implies 0 Ym = a exp(ikm∆x) 1 Ym = [a exp(iψ1) − 2E cos ψ1] exp(ikm∆x) 1 0 Requiring E = 0, we find that Ym = exp(iψ1)Ym. In this simple Case II: |λ| > 1 Recall that the roots of the quadratic are A± = −iλ sin k ∆x ± 1 − λ2 sin2 k ∆x . If |λ| > 1, there will be some wavelengths for which λ2 sin2 k ∆x > 1 . Then the two roots of the quadratic are pure imaginary A = i −λ sin k ∆x ± λ2 sin2 k ∆x − 1 and therefore either |A+| > 1 or |A−| > 1, i.e., the modulus of one of the roots will exceed unity. In that case the amplitude of the corresponding solution will increase without limit as time increases. That is, the amplitude of the solution of the finite difference equation will grow without bound for large time. This phenomenon is called computational instability. case, we can eliminate the computational mode. In general, it is much more difficult. In case of computational instability, the solution of the finite difference equation cannot possibly resemble the physical solution. The physical solution remains of constant amplitude for all time. The numerical solution grows without limit with time. We thus require that |λ| ≤ 1. The condition for stability is known as the CFL Criterion: Unconditionally Stable Schemes. ´ A large part of the research effort in Met Eireann recently has been devoted to the development of integration schemes which are free of the CFL constraint. The semi-Lagrangian scheme for advection is based on the idea of approximating the Lagrangian form of the 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. 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. We discuss semi-Lagrangian schemes in the next lecture. 22 c ∆t ≤1 ∆x after Courant, Friedichs and Lewy (1928), who first published the result. It implies that, if we refine the space grid, that is, decrease ∆x, we must also shorten the time step ∆t. Thus, halving the grid size in a two dimensional domain results in an eightfold increase in computation time. 21 ...
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