Atmospheric and Oceanic Modeling - Space-Time

Lecture 8 Space-time discretization We have so far analyzed time and space discretizations while respectively treating the complimentary dimension as continuous. Now we will consider how to discretize both space and time at once. In general, the methods we have outlined so far will work in conjunction. For example, in section 2.1, the advection equation was discretized using second, fourth and sixth order differences. Since the motions are oscillatory we know that we need to use a time-stepping scheme that is [conditionally] stable for wave-like motions, such as leap-frog, Huen or Runge-Kutta. Here, we will analyze particular combinations of space-time discretiza­ tions that yield new properties. 5.1 Forward in Time, Upwind in Space The upwind scheme uses a side-difference in space biased in the upwind direction so that for positive flow (c > 0) the scheme is: � ωin+1 − ωin c � n n + ωi − ωi−1 =0 (5.1) �t �x To find the dispersion relation for the numerical solution we substitute in a wave solution of the form e−(�+i�)t+ikx : � e−(�+i�)�t = 1 − C 1 − e−ik�x � where C = c�t is the Courant number. Separating the imaginary and real �x components gives: e−��t sin ��t = C sin k�x e−��t cos ��t = 1 − C + C cos k�x Solving for e��t and tan ��t gives: C sin k�x 1 − C(1 − cos k�x) 2C sin k�x cos k�x 2 2 = 2 k�x 1 − 2C sin 2 tan ��t = and e−2��t = (1 − C(1 − cos k�x))2 + (C sin k�x)2 k�x = 1 − 4C(1 − C) sin2 2 varies between 0 (for long waves) and 1 (for the grid-scale Since sin2 k�x 2 waves) stability depends on the sign of the quantity 4C(1 − C); if either C < 0 or C > 1 then 4C(1 − C) < 0 and the solution grows with time. Therefore stability is conditional on: 0 � 4C(1 − C) � 1 or simply 0 � C � 1. The strongest damping occurs at C = 1/2 which max­ imizes 4C(1 − C). Since C < 0 is unstable, we can infer that the downwind difference scheme is unconditionally unstable. The waves are dispersive since tan ��t depends on k. If C < 1/2 then the frequency is less than k so that the scheme is decelerating and if C > 1/2 then the frequency is either larger than k or changes sign. C = 1/2 is a special point because the denominator becomes 1 − cos 2 k�x so that the 2 k�x whole expression becomes 2C tan 2 and the frequency is exactly correct. The upwind scheme also exhibits a special property of preserving extrema. This can be seen by re-arranging the difference equation for the future un­ known value: n + (1 − C)ωin ωin+1 = Cωi−1 n This is simply a linear interpolation between ωi−1 and ωin with the C being n+1 n the sliding parameter. Hence, ωi must fall on or between ωi−1 and ωin so long as 0 � C � 1. We know from section 1.3 that both the forward in time and side difference are both of first order accuracy. We also learnt in section 4.2 that the forward 1 C=1/2 0.8 0.6 C=1/3 0.4 0.2 ωΔt/Cπ C<<1 0 C=1 −0.2 −0.4 C=2/3 −0.6 −0.8 −1 0 0.1 0.2 0.3 0.4 0.5 kΔx/π 0.6 0.7 0.8 0.9 1 Figure 5.1: The frequency of the FTUS (or “upwind”) scheme for various Courant number C = �tc . C = 1/2 falls on w = k which is the frequency of �x the continuum. scheme is generally unstable when used for oscillatory motions. It is therefore some surprise the scheme is stable at all. One way of looking at why the scheme works is to express it as a FTCS (forward in time, centered in space) scheme with a specific amount of diffusion: ω n − ωin−1 |C| ωin+1 − 2ωin + ωin−1 ω − ω = −C i+1 − (5.2) 2 2 C |C| i ∂ii ω n (5.3) = −C∂i ω n + 2 In this form, the upwind scheme appears as a space centered derivative but �x2 that is required and with a diffusion term with diffusion coefficient |C| 2 �t sufficient enough to make the scheme stable. As a final note, advection schemes are often most useful when written in flux form (i.e. as a divergence of a flux). The last form allows us to write the scheme: ω n+1 − ω n 1 ∂i F =− �t �x where F is defined as |c| i F = cω n − ∂i ω n 2 n+1 n ⎤ � which a general way to write the upwind flux. 5.2 The Lax-Wendroff Method In a similar vain, the Lax-Wendroff method adds diffusion to the FTCS scheme: ⎡ ⎣ ωn+1 − ω n cC n 1 i n ∂i ω ∂i cω − (5.4) =− 2 �t �x where the diffusion term has an implied diffusion coefficient equal to c2 �t/2. That is, the last term is an approximation to: c2 �t �xx ω 2 Although the scheme is written as a forward difference in time, it is in fact second order accurate in time and space; the truncation error of the forward difference on the LHS is: �t �t c2 �t �t (−c�x ω) = �xx ω �tt ω = 2 2 2 The diffusion term therefore cancels the leading truncation error from the forward time difference. This is an example of how treating time and space together leads to a substantially different scheme than would be obtained by discretizing the dimensions independently. An alternative way of accessing the second order nature of the LaxWendroff scheme is to break it down into a two stage scheme: 1 �t c ∂i ω n 2 �x �tc �n+ 1 n = ω − ∂i ω 2 �x i ω �n+ 2 = ω n − ω n+1 (5.5) (5.6) Here, the time marching looks likes the mid-point second order Runge-Kutta method and the mid-point values are staggered in space. The dispersion relation for the Lax-Wendroff scheme is: 2C sin k�x cos k�x 2 2 tan ��t = 1 − 2C 2 sin2 k�x 2 and amplification: e−2��t = 1 − 4C 2 (1 − C 2 ) sin2 k�x 2 These expressions both take a similar form to those of the FTUS scheme except for the second order dependence on C. However, the Lax-Wendroff scheme does not conserve extrema like the FTUS method does. 5.3 Flux limiters The FTUS and Lax-Wendroff methods each have their advantages and dis­ advantages; the FTUS is only first order accurate and very diffusive but does conserve extrema while the Lax-Wendroff scheme doesn’t conserve extrema but is second order accurate. Our objective in this section is to try to blend these two schemes, capturing the desired features of each. First, we write the system in flux form: � 1 � n+1 1 (F 1 − Fi− 1 ) ωi − ωin = − 2 �t �x i+ 2 And now we cast the advective flux, F , as some unknown combination of the upwind flux, F U S , and Lax-Wendroff flux, F LW : Fi+ 1 = �i+ 1 F LW + (1 − �i+ 1 )F U S 2 2 2 US Fi+ = cωi 1 2 LW = cωi + Fi+ 1 2 c(1 − C) (ωi+1 − ωi ) 2 so that the advective flux is: Fi+ 1 = cωi + 2 c(1 − C) �i+1/2 (ωi+1 − ωi ) 2 The factor �i+ 1 is some function, yet to be determined. In some texts, 2 this form of the flux is justified by casting the Lax-Wendroff as above and arguing that the second term is a correction to the upwind flux and hence adjustable. Substituting into the prognostic equation gives: ωin+1 = ωin ⎤ � C(1 − C) C(1 − C) n n − C− �i− 1 (ωin − ωi−1 )− �i+ 1 (ωi+1 − ωin ) 2 2 2 2 ψ = 2r ψ(r ) ψ= 2 2 1 r 1 2 3 Figure 5.2: The two lines indicate the bounds on the limiter function, �(r) given by the constraints for the scheme to be TVD (total variance diminishing). The last term can be re-written as: C(1 − C) C(1 − C) �i+ 12 n n n �i+ 1 (ωi+1 (ω − ωi−1 ) − ωin ) = 2 2 2 ri+ 1 i 2 where ri+ 1 is the slope ratio which is defined as: 2 ri+ 1 = 2 n ) (ωin − ωi−1 n n (ωi+1 − ωi ) Now we will re-write the prognostic equation using the last expression: ωin+1 � ⎦ (1 − C) (1 − C) �i+ 12 � n n = ωin − C �1 − �i− 1 + (ωi − ωi−1 ) 2 2 2 ri+ 1 2 which looks likes an FTUS (upwind) scheme but with a modified Courant number. The upwind scheme is both monotone and stable if the “effective” Courant number is both non-negative and less than one: � ⎦ (1 − C) (1 − C) �i+ 12 � 0 � C �1 − �i− 1 + �1 2 2 2 ri+ 1 2 ψ= r ψ(r ) 2 1 ψ= 1 r 1 2 3 Figure 5.3: The two lines indicate the special limiter functions, �(r) = 1, which yields the Lax-Wendroff flux, and �(r) = r which yields the Warming and Beam flux. Since both these schemes are of second order, any linear combination of these schemes is also of second order. The shaded region between them is thus the space where �(r) will yield a second order scheme. or �i+ 1 −2 2 2 − �i− 1 � � 2 ri+ 1 C 1−C 2 Since C � 0 then the above is satisfied if the following stronger constraint is satisfied: � � � �� 1 � � i+ 2 � − �i− 1 �� � 2 � 2 � � ri+ 1 2 Now we will allow the “limiter function” � be a function of the slope ratio: � = �(r) If r < 0, the slope must have changed sign and indicates a local extrema. In this instance we should limit the advective flux to take the form of the upwind flux since it is the only (linear) scheme capable of conserving extrema. Therefore, for r < 0 we set �(r) = 0. If r > 0, the above inequality is satisfied when 0� �(r) � 2 and 0 � �(r) � 2 r ψ = 2r ψ(r ) ψ= r ψ= 2 2 Superbee Van Leer 1 ψ= 1 minmod r 1 2 3 Figure 5.4: The shaded region is the intersection between the TVD region and region of second order accuracy. The short dashed line is the Superbee limiter, the long dash line is the minmod limiter and the solid curve is the Van Leer limiter. and we must find functions that meet these criteria if the new scheme is to conserve extrema. The limits on �(r) are indicated by the shaded region in Fig. 5.2. This region is said to be TVD (total variance diminishing) which means that the norm of gradients of a field can not be increased by the scheme. It happens that the FTUS scheme is TVD. Although we haven’t directly use the TVD definition or cast the constraints as associated with TVD behaviour we have nevertheless derived constraints on a non-linear flux limiter that will make it TVD. A further constraint on the limiter function is the preferred order of accu­ racy of the resulting scheme. We don’t show the proof here but if the scheme can be shown to be an interpolation of the Lax-Wendroff method (�(r) = 1) and Warming and Beam method (which corresponds to �(r) = r) then the scheme will be of second order accuracy. �(r) must therefore fall in the region indicated in Fig. 5.3 to be second order accurate. Finally, if �(r) is chosen to fall in the intersection of these two regions, second order and TVD, then the resulting scheme will be both TVD and second order accurate. There are many possible limiters but the most widely used are: • Superbee: �(r) = max(0, min(1, 2r), min(2, r)) due to Roe, 1985, • minmod: �(r) = max(0, min(1, r)) due to who?, • Van Leer: �(r) = r+|r| 1+|r| due to Van Leer, 1974. each of which is shown in Fig. 5.4. A comparison of many different advection schemes is illustrated in Figs. 5.5 and 5.6. We have grouped the schemes in the following way: a) linear schemes with upwind bias (odd order), b) linear space centered schemes, c) second-order flux limited schemes (as discussed above) and d) higher order flux limited schemes. The reason for the grouping is the common behaviour within each group. The upwind biased schemes are smoother than the spacecentered schemes. All the linear schemes, except for the upwind scheme, have false extrema. The flux limited schemes are monotone. The second order flux limited schemes tend to have reduced amplitude of extrema due to diffusion. The higher order flux limited methods can correct this. 5.4 Modified equations In section 5.2 we implicitly made use of a “modified equation”. Modified equations provide a way of anticipating the behaviour of a numerical method by examing terms implied by a particular approximation. For a given differ­ ence equation, the corresponding modified equation is a continuous equation for which the same difference equation is a higher order approximation! For example, the difference equation 1 n+1 c n (ωi − ωin ) + (ωi − ωin−1 ) = 0 �x �t (5.7) is the F.T.U.S. scheme (same as equation 5.1), which is a first order approx­ imation in time and space to �t ω + c�x ω = 0. However, it is also a second order in time and space to the continuous equa­ tion �t c�x �tt ω + c�x ω − �xx ω = 0. (5.8) �t ω + 2 2 This is because the O(�t) truncation term arising from 1 n+1 �t �t2 �tt ω + (ωi − ωin ) = �t ω + �ttt ω + . . . 2! �t 3! exactly matches the second term in (5.8) leaving the O(�t2 ) as the timetruncation error. Similarly, the O(�x) truncation term arising from the side difference �x �x2 c n �xx ω + c (ωi − ωin−1 ) = c�x ω − c �xxx ω − . . . 2! �x 3! exactly matches the last term in (5.8). We can eliminate the second order time derivative from (5.8) by differ­ entiating the governing equation:�tt ω = c2 �xx . Thus, the first order approxi­ mation to the advection equation (5.7) is also a second order approximation to the modified equation �t ω + c�x ω − c�x c�t (1 − )�xx ω = 0. 2 �x This allows us to interpret the F.T.U.S. scheme; the modified equation dif­ feres from the governing equation by a diffusive term with diffusivity c�x c�t ). (1 − 2 �x This effective diffusivity is propotional to the flow c, it decreases with de­ creasing �x (i.e. higher resolution is less diffusive) and also decreases with the Courant number c�t/�x. This is consistent with the special case of c�t/�x = 1 for which the F.T.U.S. scheme gives the exact answer. Now we’ll take the approach a step further and derive a third order mod­ ified equation. To get there we will use the following relations which are simply repeated differentiations of the second-order modified equation (5.8): c�x (1 − C)�xxx ω 2 c�x −c�xxx ω + (1 − C)�xxxx ω = −c�xxx ω + O(�x) 2 c�x −c�xt ω + (1 − C)�xxt ω 2 ⎡ ⎣ c�x c�x (1 − C)�xxx ω + (1 − C) (−c�xxx ω) + O(�x2 ) −c −c�xx ω + 2 2 c2 �xx ω − c2 �x(1 − C)�xxx ω + O(�x2 ) −c�xxt ω + O(�x) = c2 �xxx + O(�x) −c�xtt ω + O(�x) = −c3 �xxx + O(�x) �xt ω = −c�xx ω + �xxt ω = �tt ω = = = �xtt ω = �ttt ω = where C = c�t . �x Now we substitute in for the first and second order Taylor series terms in the F.T.U.S scheme: 1 n+1 c n (ω − ωin−1 ) (ωi − ωin ) + �t �x i �t �x �t2 �x2 = �t ω + c�x ω + �tt ω − c �xx ω + �ttt ω + c �xxx ω + . . . 2! 2! 3! 3! c�x2 c�x = �t ω + c�x ω − (1 − C)�xx ω + (1 − C)(1 − 2C)�xxx ω (5.9) 2 6 keeping only O(�x2 ) terms. This is the modified equation to which (5.7) is an O(�x3 ) approximation. The second term is diffusive as we saw before with the second order modified equation. The last term in (5.9) causes waves to be dispersive (it leads to a −ik 3 term in the dispersion relation). Time = 1 CFL = 0.05 Num. pnts.=100 Analytic solution upwind−1 upwind−3 DST−3 Θ 1 0.5 0 0 0.1 0.2 0.3 0.4 0.5 X 0.6 0.5 X 0.6 0.5 X 0.6 0.5 X 0.6 0.7 0.8 0.9 1 Analytic solution centered−2 centered−4 Lax−Wendroff Θ 1 0.5 0 0 0.1 0.2 0.3 0.4 0.7 0.8 0.9 1 Analytic solution minmod Superbee van Leer (MC) van Leer (alb) Θ 1 0.5 0 0 0.1 0.2 0.3 0.4 0.7 0.8 0.9 1 Analytic solution DST Sweby µ=1 DST Sweby µ(c) Θ 1 0.5 0 0 0.1 0.2 0.3 0.4 0.7 0.8 0.9 1 Figure 5.5: Solution obtained with a small Courant number. The analytic solution is the thick solid line in each panel. The upwind scheme does not exhibit false extrema but is clearly very diffusive. The third order upwind is much better at preserving the shapes. The even order methods (second panel) have multiple false extrema but the Lax-Wendroff method is at least smoother. The third panel shows the second order limited solutions, all of which conserve extrema. The fourth panel shows some third order flux limited solutions which preserve amplitude better. Time = 1 CFL = 0.24 Num. pnts.=100 Analytic solution upwind−1 upwind−3 DST−3 Θ 1 0.5 0 0 0.1 0.2 0.3 0.4 0.5 X 0.6 0.5 X 0.6 0.5 X 0.6 0.5 X 0.6 0.7 0.8 0.9 1 Analytic solution centered−2 centered−4 Lax−Wendroff Θ 1 0.5 0 0 0.1 0.2 0.3 0.4 0.7 0.8 0.9 1 Analytic solution minmod Superbee van Leer (MC) van Leer (alb) Θ 1 0.5 0 0 0.1 0.2 0.3 0.4 0.7 0.8 0.9 1 Analytic solution DST Sweby µ=1 DST Sweby µ(c) Θ 1 0.5 0 0 0.1 0.2 0.3 0.4 0.7 0.8 0.9 1 Figure 5.6: Solution obtained with a large Courant number. As for Fig. 5.5 but notice that the noise levels in the unlimited schemes is worse.