Atmospheric and Oceanic Modeling - Spatial Discretization and Numerical Dispersion

Lecture 2 Spatial discretization and dispersion We have out-lined the basic method of forming finite difference approxima­ tions to equations (at least for one-dimensional). Now we’ll examine how the choice of discretization can impact the behavior of solutions and in par­ ticular we will examined wave motions and propagation on discrete grids. For the purposes of this section we will treat the time-derivative as perfect; this allows us to analyze the dynamical impact of spatial difference without worrying about the stability and accuracy of the time differencing. 2.1 Advection equation The simplest equation to consider is the scalar advection by a constant flow field: �t ω + c�x ω = 0 For the solution at any time t + �t is ω(x, t + �t) = ω(x − �tc, t) which means that shape is preserved and the whole solution simply translates with time. If we represent the spatial structure of the solution using Fourier series, each mode can be seen to be a traveling wave: ω(x, t) = ωo ei(kx−�t) and the frequency of each mode satisfies a dispersion relation obtained by substituting the form of the solution into the governing equation: � = ck The phase speed �/k = c and group speed �k � = c are equal and constant the waves are said to be non-dispersive. Approximating the spatial derivative with the centered second order deriva­ tive yields: c i dt ω + ∂i ω = 0 �x or c (ωi+1 − ωi−1 ) = 0 d t ωi + 2�x Substituting in our solution of the form ei(kx−�t) gives: � c � ik�x e − e−ik�x 2�x ci = − sin k�x �x −i� = − In the limit of small �x, the frequency of long waves (k�x << 1) is: �= c �x�0 sin k�x = ck �x so the approximation appears to converge. However, the numerical waves are dispersive since the group and phase speeds are not constant with wave number. Also, the frequency of the grid-scale waves (wavelength of 2�x � k = η/�x) is always zero regardless of the resolution. Now consider the fourth order accurate representation of the spatial derivative: � �i c 1 ∂i ω − ∂ii ω = 0 dt ω + �x 6 or c ∂i (ωi−2 − 8ωi−1 + 8ωi+1 − ωi+2 ) = 0 dt ω + 12�x Substituting in a wave solution yields the dispersion relation: � c � 2ik�x −e + 8eik�x − 8e−ik�x + e−2ik�x 12�x ci = (− sin 2k�x + 8 sin k�x) 6�x −i� = These dispersion relations are plotted in Fig. 2.1 and they show how the frequency tends to zero for the grid-scale modes. We plot the non-dimensional 3.5 3 ω’ =πk’ ω’2=2skck 2 ω’4=2skck(1+2/3sk ) 2 2 ω’6=2skck(1+2/3sk (1+4/5sk )) ω’=ωΔx/c 2.5 2 1.5 1 0.5 0 0 0.2 0.4 k’=kΔx/π 0.6 0.8 1 Figure 2.1: Dispersion relations for constant flow advection in the contin­ uum, and using second and fourth order spatial differences. frequency, � � = ��x/c, as a function of the non-dimensional wave number, k � = k�x/η. The fourth order method is more accurate or closer to the continuum for a larger range of wave numbers but because the frequency then has to drop to zero over a shorted span of wave numbers the group speed error is larger. Finally, the sixth order method: c 1 1 ∂i ω − ∂ii ω − ∂ii ω dt ω + 5 �x 6 has the dispersion relation: � ⎞ � ⎞ ��i ⎟⎟ c 4 k�x 4 k�x ) 1 + sin ( ) �= sin (k�x) 1 + sin ( 2 5 2 �x 6 Again, when plotted we see that the frequencies are more accurate but the grid-scale mode is stationary and so the group speed of the shortest waves is even worse than with the fourth order scheme. a) u,η u,η u,η i−1 i i+1 x b) η η u u η x i−1 i−1/2 i i+1/2 i+1 Figure 2.2: A regular grid with the gravity wave variables u and � co-located on an non-staggered grid (a) and on a staggered grid (b). The pattern that emerges is that the grid-scale modes are always station­ ary and that increasing the order of the scheme does not necessarily improve the behavior of the shortest scales. 2.2 Gravity waves Consider the adjustment of a non-rotating, constant depth shallow water layer. Dropping non-linear terms, the governing equations in one dimension are: �t u = −g�x � �t � = −H�x u where g is gravity, H is the nominal depth of the layer, u is the flow and � the free-surface displacement. These first order equations can be re-arranged into a second order wave equation: �tt � = gH�xx � If we substitute in a wave-like solution, � = �o ei(kx−�t) we obtain an equation for the frequency: � � 2 = gHk 2 or � = ± gHk which is quadratic and has two roots, meaning that there is a pair of modes (waves can travel in either direction). � The waves are non-dispersive and have group and phase speeds of cg = ± gH. Now consider the second order spatial discretization of the gravity wave equations on a regular grid corresponding to Fig. 2.2a: g ∂i � i �x H ∂i ui �t � = − �x (2.1) �t u = − (2.2) (2.3) We can follow the same procedure of substituting the first equation into the second to obtain: gH �tt � = ∂ii � ii �x2 The discrete operator on the right hand side in full form is: ∂ii � ii = 1 (�i−2 − 2�i + �i+2 ) 4 and so substituting in the solution � = �o ei(kx−�t) into the discrete wave equation gives: � gH � −i2k�x i2k�x e − 2 + e 4�x2 gH (2 cos 2k�x − 2) = 4�x2 k�x k�x 4gH cos2 = − sin2 2 2 2 �x −� 2 = Here, we’ve expressed the frequency as a function of k�x/2 since 0 � sin k�x � 1 and is monotonic. Note that the largest discrete wave num­ 2 ber permitted is the Nyquist mode, k = η/�x. In order to plot these relations, we will use a non-dimensional wave num�k ��x ber k � = �x and plot � . The curves for the continuum and non-staggered gH grid are plotted in Fig. 2.3. From these curves we note the following: • the numerical gravity waves are dispersive while the continuum is nondispersive, • a false extrema in the frequency occurs at k = � 2�x or k � = 12 , • the effective group speed has the wrong sign for modes higher than � k = 2�x (or wavelengths shorter than 4�x), 3.5 3 ω’=πk’ ω’=2sin(πk’/2)cos(πk’/2) ω’=2sin(πk’/2) ω’=ωΔx/cg 2.5 2 1.5 1 0.5 0 0 0.2 0.4 k’=kΔx/π 0.6 0.8 1 Figure 2.3: Dispersion of numerical gravity waves for the non-staggered grid (green) and the staggered grid (red). The continuum (� = k) is plotted for comparison (blue). • the Nyquist mode (k = frequency). � �x or k � = 1) is stationary (i.e. has zero The last two points are symptomatic of the presence of a computational or non-physical mode and these often occur when averaging operators (e.g. � i ) are present. In this case, the two-grid length wave (a pattern of plus, minus, plus, minus, ...) is not “seen” by the discrete operators and so is allowed to be stationary. When the group speed has the wrong sign, energy is not radiated properly in the adjustment process and is often seen as “noise” near sources of energy. 2.2.1 The staggered grid We saw earlier that a staggered difference (e.g. ∂i �) is four times more accurate than the non-staggered difference (e.g. ∂i � i ). Staggered differences become natural choices if we stagger the model variables in space as shown in Fig. 2.2b. The discrete model then becomes: g ∂i � �x H ∂i u �t � = − �x �t u = − The same analysis can be performed as before but here we’ll introduce the matrix method which performs the same process of elimination but more symbolically. The system can be written: ⎞ g ∂ �x i �t H ∂ �x i �t ⎟⎞ u � ⎟ =0 On substituting a wave solution, the operators �t and ∂i have “response functions” of −i� and 2i sin k�x/2 respectively. ⎞ −i� 2iH sin k�x �x 2 2ig �x sin k�x 2 −i� ⎟⎞ u � ⎟ =0 The determinant of the operator matrix is: �2 − 4gH k�x =0 sin2 2 2 �x The dispersion relation is plotted in red in Fig. 2.3 and we note the following: • when compared to the continuum we see that the numerical modes are still dispersive, • there is no false extrema, unlike the non-staggered grid, • the group speed is therefore of the correct sign everywhere, even if reduced. The remarkable improvement in behavior of the staggered grid over the nonstaggered grid is due to the complete absence of operators with null-spaces (such as the averaging operator). The formal truncation error is smaller but the improved accuracy is not responsible for the qualitative improvements. 2.2.2 Inertia-gravity waves in 1-D We’ll now consider the geostrophic adjustment problem which introduces rotation. The relevant linear shallow water equations in one dimension are: �t u − f v + g�x � = 0 �t v + f u = 0 �t � + H�x u = 0 (2.4) (2.5) (2.6) where gradients in the y direction have been dropped. This system adds one equation and one degree of freedom over the previous gravity wave problem. The dispersion relation can be obtained using the matrix method, represent­ ing the system by: � ⎛� ⎛ �t −f g�x u � � ⎜ �t 0 ⎜ � f ⎝� v ⎝ = 0 H�x 0 �t � and taking the determinant of the operator matrix �� � � −i� �� �� f � � Hik ⎛� ⎠ −f gik �� �=0 ⎜� −i� 0 ⎝� = 0 � � � 2 = f 2 + gHk 2 0 −i� � For discretizing the model, we consider the four possible configurations of variables on a regular grid, illustrated in Fig. 2.4 and labeled A, B, C and D. The four discrete models based on each grid are then: • A-grid model g ∂i � i = 0 �x �t v + f u = 0 H �t � + ∂i ui = 0 �x �t u − f v + • B-grid model g ∂i � = 0 �x �t v + f u = 0 H �t � + ∂i u = 0 �x �t u − f v + A) u,v,η u,v,η u,v,η i i+1 x i−1 B) η η u,v u,v η x i−1 C) v,η i−1/2 i i+1/2 i+1 u v,η u v,η x i−1 D) u,η i−1/2 i i+1/2 i+1 v u,η v u,η x i−1 i−1/2 i i+1/2 i+1 Figure 2.4: The inertia-gravity wave variables u, v and � in the four possible configurations (A, B, C and D) on a finite difference grid. • C-grid model g ∂i � = 0 �x �t v + f ui = 0 H �t � + ∂i u = 0 �x �t u − f v i + • D-grid model g ∂i � i = 0 �x �t v + f ui = 0 H ∂i ui = 0 �t � + �x �t u − f v i + The corresponding dispersion relations are: A: �2 4L2d 2 2 + s c = 1 f2 �x2 k k r=(2L /Δ x)=2 d 12 10 2 2 2 ω=1+π r k /4 2 2 ω =1+r2s c A k k 2 2 ωB=1+r sk 2 2 2 ω =c +r s C k k 2 2 2 ω =c +r2s c D k Continuum k k ω’=ω/f 8 6 B C 4 2 A a) D 0 0 0.2 0.4 k’=kΔx/π 0.6 0.8 1 r=(2L /Δ x)=0.5 d 1.8 1.6 1.4 2 2 2 ω=1+π r k /4 2 2 ω =1+r2s c A k k 2 2 ω =1+r s B k 2 2 2 ωC=ck +r sk 2 2 2 2 ω =c +r s c D k Continuum B k k ω’=ω/f 1.2 A 1 0.8 0.6 0.4 C b) 0.2 0 0 D 0.2 0.4 k’=kΔx/π 0.6 0.8 1 Figure 2.5: The one dimensional dispersion relations for inertia-gravity waves in the continuum and for the four grids, (A, B, C and D) for (a) r = 2 and (b) r = 1/2. B: C: D: �2 4L2d 2 = 1 + s f2 �x2 k �2 4L2d 2 2 = c + s k �x2 k f2 �2 4L2d 2 2 2 = c + s c k �x2 k k f2 � and where Ld = gH/f is the radius of deformation and sk = sin k�x 2 k�x ck = cos 2 and short-hand of the common factors. The coefficient in front of the gravity-wave contribution (sk2 ) is a resolution parameter: r2 = � 2Ld �x �2 and is controls the form of the some of the dispersion relations. The dis­ persion relation of the continuous and four numerical inertia-gravity waves are plotted in Fig. 2.5. In panel (a) the resolution parameter is larger than one meaning the deformation radius is resolved. In panel (b) the resolution parameter is less than one meaning the deformation radius is not resolved. In the continuum, the frequency increases monotonically with wave number, the waves are dispersive and the group and phase speeds are pointed in the direction of the wave vector. None of the numerical schemes behave precisely in this manner: • Scheme A, regardless of resolution has a false maximum at k � = 1/2 so that grid-scale waves (k � = 1) oscillate with the inertial frequency rather than propagate as gravity waves, • Scheme D has false maxima but is worse in that the grid-scale waves become stationary, • Scheme C has the wrong slope (group speed in wrong direction) if r < 1 but is qualitatively correct if r > 1, • Scheme B is qualitatively correct for both values of r and so appears to be the best scheme. This analysis suggests that staggering variables in the form of the B grid is most likely to avoid computational modes when solving these one-dimensional shallow water equations. However, the results of an analysis in two dimen­ sions are different.