Atmospheric and Oceanic Modeling - Two-Dimensions

Lecture 10 Advection in two dimensions 6.1 Stability of multiple terms (in multiple dimensions) When we analyzed the stability of time-stepping methods we tended to con­ sider either a single damping term or a single oscillatory term. Since the forward method is easy to analyze and stable for damping terms, consider the two-dimensional equation: �t ∂ = −κ∂ + ��xx ∂ + ��yy ∂ The second order space, forward in time difference equation is: � � �n ∂ n+1/2 = −κ + �jj ∂ �ii ∂ + 2 �x �y 2 As before, to do the stability analysis, we substitute in a solution of the form ei(kx−�t)−�t . The response function for the spatial operators are: k�x 2 ) = −4s2k 2 l�y 2 R(�j �j ∂) = (2i sin ) = −4s2l 2 R(�i �i ∂) = (2i sin and sl = sin l�y are convenient short-hand. The amplifi­ where sk = sin k�x 2 2 cation equation then is: e−��t = 1 − �tκ − 4��t 2 4��t 2 s − s �x2 k �y 2 l so that for stability we must satisfy: �tκ + 4��t 4��t + �1 �x2 �y 2 Each term contributes to the numerical stability additively. The consequence is that two-dimensional problems typically have a smaller stability thresh­ old than one-dimensional problems and similarly three-dimensional problems more so. This effect is generally true when ever the terms are treated uni­ formly by the time-stepping scheme. Schemes that treat some terms differ­ ently (say implicitly) can avoid this additive effect of stability criteria. In the following sections we’ll see how to avoid this effect for the forward advection schemes. 6.2 First order upwind in two dimensions Let’s apply the FTUS scheme to advection in two dimensions by constant flow (u, v). Without any other considerations the resulting scheme is: u � n v � n 1 � n+1 n n + ∂ − ∂i,j−1 =0 ∂i,j − ∂in−1,j + ∂i,j − ∂i,j �x �y i,j �t (6.1) where we have assumed the flow, u and v, to be in the positive x and y directions respectively. A stability analysis shows the scheme to be conditionally stable: u�t v�t u�t v�t � 0 and + �1 � 0 and �x �y �x �y but notice that the condition is on sum of the Courant number in each direction. The modified equation corresponding to 6.1 can be found by first reexpressing the scheme in terms of centered differences: 1 1 1 u 1 v n i j (�i ∂ − �ii ∂) + (�j ∂ − �jj ∂) = 0 (�n ∂ + �nn ∂) + �x 2 �y 2 �t 2 which is a second order approximation of �t ∂ + u�x ∂ + v�y ∂ = u�x u�t v�y v�t (1 − )�xx ∂ + (1 − )�yy ∂ − uv�t�xy ∂ �y 2 �x 2 The last term is an an-isotropic term; it looks different along different direc­ tions. In fact, it leads to diffusion in the direction tangential to the flow and anti-diffusion in the direction of flow. The effect of this term leads to the shape distortion which is apparent in the solution shown in Fig. 6.1a. Fig. 6.2a visualizes the fluxes implied by the two-dimensional FTUS scheme. Each flux is normal to a coordinate line (cell boundary) and sweeps out a length �tu or �tv so that the area (shaded regions) transported through the cell wall is �tu�y and �tv�x for the x and y directions re­ spectively. The western and southern fluxes are self explanatory. However, the northern and eastern flux appear to overlap in the area that is swept out (indicated by shaded region 3). The scheme is conservative because the sum of regions 1 and 2 minus the double counted region 3 is equal to �x�y. How­ ever, the lack of dependence on the cell to the south-west mean the scheme can not represent the characteristic properties of diagonal flow. 6.3 Corner transport upstream The CTU algorithm can be derived by various means but the most intu­ itive is based on the same principles as used in the finite volume method. Consider the Galilean translation of the fluid parcels resulting from constant flows (u, v). Following the location of the grid at time level n + 1 back in time along the characteristics we arrive at the picture indicated in Fig. 6.2b where the translated cell overlaps with the three upstream cells. The area of intersections represent the weights to attribute to the tracer on the original grid but at the previous time-level, n. The resulting scheme is then: n+1 n n ∂i,j = (1 − Cu )(1 − Cv )∂i,j + Cu (1 − Cv )∂in−1,j + (1 − Cu )Cv ∂i,j−1 + Cu Cv ∂in−1,j−1 where Cu = u�t/�x and Cu = u�t/�x are the Courant numbers in the x and y directions. Stability is conditional: 0 � Cu � 1 and 0 � Cv � 1 which is the same as satisfying the stability of FTUS in each direction allows a much longer time-step than using the two-dimensional FTUS. The CTU method can be written flux form with a contribution to the x and y fluxes from the corners. This is awkward to visualize and introduces 1st order upwind 1st order corner transport 0.8 0.8 0.6 0.6 0.4 0.4 0.2 0.2 0.2 0.4 0.6 0.8 0.2 2nd order Van Leer flux limiter 0.8 0.6 0.6 0.4 0.4 0.2 0.2 0.4 0.6 0.8 0.6 0.8 True solution 0.8 0.2 0.4 0.2 0.4 0.6 0.8 Figure 6.1: Advection by a uniform diagonal flow (u = v) using a) the FTUS applied in each direction, b) the two dimensional upstream corner transport method, c) the multi-dimensional second order Van Leer flux limiter and d) the true solution. a) b) C u Δx C u Δx 3 Δy 1 C v Δy 2 1 Δy 1 2 3 1 Δx C v Δy Δx Figure 6.2: A schematic of the volumes swept out by the fluxes in a) scheme 6.1, the simple two-dimensional upwind method, and b) the corner transport method. wider stencils for each flux than the simple FTUS flux. Another way to write the CTU method is via an intermediate step: � n n ∂i,j = (1 − Cu )∂i,j + Cu ∂i−1,j n+1 � � ∂i,j = (1 − Cv )∂i,j + Cv ∂i,j−1 which is exactly equivalent to previous form. This last form can be inter­ preted as applying FTUS successively in each direction, using the latest val­ ues at each stage. We can see here why the stability depends simply on the stability in each direction since both stages must be independently stable. 6.4 Multi-dimensional algorithm As a rule, the forward advection methods (Lax-Wendroff, flux limiters, ...) will all need a similar treatment to that that rendered the FTUS into the CTU method. A general approach is as follows: �t �t �i Fx (u, ∂ n ) + ∂ n �i u �x �x �t �t = ∂ n+1/m − �j Fy (v, ∂ n+1/m ) + ∂ n �j u �y �y ∂ n+1/m = ∂ n − ∂n+2/m .. . ∂ n+1 = ∂ n+m/m +∂ n � �t �t �i v + . . . �i u + �y �x � where Fx (u, ∂) is a general form of the advective flux. This scheme reproduces the CTU method for constant flow if Fx is the upwind flux. The last term in each equation account for the “apparent” divergence of the flow when treating each direction separately. This multi-dimensional method is a form of operator splitting. We should emphasize that it is not as necessary to consider multi-dimensional truncation errors when the time-difference is centered in time (e.g. using leap­ frog) since the erroneous cross terms that appeared in the two-dimensional FTUS scheme appeared from the O(�t) truncation term of the forward time difference. Such terms do appear in the modified equations of centered time difference but don’t generally take a diffusive form and can only appear at order O(�t2 ).