3.3.2: Spectral space differences
The use of spatial finite differences, as we saw in 3.3.1, introduces errors in the space derivatives, resulting in a computational phase speed slower than the true phase speed, especially for short waves. One approach to solve this problem is to use a spectral expansion for the space discretization, rather than finite differences. The basis functions chosen for the discretization should be orthogonal and satisfy the boundary conditions. For example, in one dimension, periodic boundary conditions suggest the use of complex Fourier series.
Consider a periodic domain of length L, and scale x by =(2/L). If we use discrete complex Fourier series truncated to include wavenumbers up to K, the spectral representation is:
. (3.3.2)
where A-k=Ak*, and the star represents the complex conjugate.
Alternatively, (3.3.2) can be written with real Fourier series as
where
There are 2M+1 distinct real coefficients that are determined by
.(3.3.3)
where we have used the orthogonality property
(1 if k=l and 0 otherwise). If JM=2M+1, the grid representation (3.3.2) and the spectral representation (3.3.3) contain the same number of degrees of freedom, and the same information.
Then, in the wave equation we can discretize U in space as in (3.3.2), and then compute the space derivative analytically:
If we neglect the time discretization errors, as before, and assume solutions of the form , we find that c'=c, i.e., the computational phase speed is equal to the true speed. The space discretization based on a spectral representation is extremely accurate (the space truncation errors are of "infinite" order). This is because the space derivatives are computed analytically, not numerically, as is the case in finite differences.
If the PDE is nonlinear, for example , then both the grid-point ("physical space") representation and the spectral representation are very useful: nonlinear products are computed efficiently in physical space, and derivatives are computed efficiently and accurately in spectral space. This is the so-called transform method used for spectral models: the space derivative is computed in spectral space, then U is transformed back into grid space, and the product is computed in grid space. We will see later that in order to avoid nonlinear instability introduced by aliasing of wavenumbers beyond K that appear in quadratic terms, the grid representation requires about 3/2 as many points as the minimum number of points required for a linear transform (JM=2K). For this reason the new values of U at time (n+1)t are stored in their spectral representation, which is more compact.
We can use von Neumann’s criterion to determine the maximum time step allowed for stability using for example the leap-frog time scheme. The FDE is
Assuming solutions for the wave equation of the form , we obtain that the amplification factor is , and in order to have we need to satisfy the stability condition .
Since the highest wavenumber present corresponds to L=2x, the stability criterion for spectral models is therefore . So, the stability criterion is more restrictive for spectral models than for finite difference models, but this is compensated by the fact that the accuracy, especially for shorter waves is much higher, and therefore fewer short waves need to be included.
Global models use as basis functions spherical harmonics, which are the eigenfunctions of Laplace's equation on the sphere:
The spherical harmonics are products of Fourier series in longitude and associated Legendre polynomials in latitude:
,
where =sin, m is the zonal wavenumber and n is the “total” wavenumber in spherical coordinates (as suggested by the Laplace equation). Pnm are the associated Legendre polynomials in x=cos , where =/2- is the colatitude. For example, the P00=1; P10=cos; P11=sin; P20=1/2 (3cos2-1); P21=3sincos; P21=3sin2;…
Using "triangular" truncation
the spatial resolution is uniform throughout the sphere. This has a major advantage over finite differences based on a latitude-longitude grid, where the convergence of the meridians at the poles requires very small time steps. Although there are solutions for this “pole problem” for finite differences, the natural approach to solve the pole problem for global models is the use of spherical harmonics. Williamson and Laprise (1998) provide a comprehensive description of numerical methods for global models.
Fig. 3.7 (reproduced from Williamson and Laprise, 1998) shows the shape of three spherical harmonics with total wavenumber n=6, and zonal wavenumber m=0, 3 and 6. Note that the distance between neighboring maxima and minima is similar for the 3 harmonics, and is associated with the “total” (2-dimensional) wavenumber n.
(Fig. 3.7 here)
3.3.3: Semilagrangian schemes
Another numerical method that has become very popular for NWP is the semi-lagrangian scheme. The equations of motion, as we have seen, can in general be written as conservation equations
du/dt=S(u),
where the left hand side of the equation represents a total time derivative (following and individual parcel) of the vector of dependent variables u. The total time derivative (also known as individual, substantial or lagrangian time derivative) is conserved for a parcel, except for the changes introduced by the source or sink S.
In a lagrangian scheme, one follows individual parcels (by transporting them with the 3-dimensional fluid velocity), and then add the source term. However, this is not convenient because one has to keep track of many individual parcels, and with time they may “bunch up” in certain areas of the fluid, and leave others without parcels to track. In the semi-lagrangian scheme, we use a regular grid, like in the previous schemes discussed so far (which are denoted eulerian). However, at every new time step we find out where the parcel arriving at every grid point (called arrival point AP) came from in the previous time step (called departure point DP). The conserved value of u along this trajectory is obtained interpolating the values of the grid points surrounding the departure point. It is evident from the schematic Fig. 3.7 that, because there is no extrapolation, the semilagrangian scheme is absolutely stable with respect to advection.
Fig. 3.8: Schematic of the semilagrangian scheme.
The semilagrangian scheme can be then written, for example, as
(Ujn+1)AP=(Un-1)DP+2t S(Un)MP
in a 3-level time scheme, or in a two time level as
(Ujn+1)AP=(Un)DP+t S(Un)DP.
For a nonlinear equation ,
QAP=QDP+S(Q)
In this case, the departure point has to be determined from the trajectory integrated between the departure and the arrival points, for example as
.
Since u evolves with time, UAP and UDP are not known until the departure point has been determined, so that this is an implicit equation that needs to be solved iteratively. For 3-level semi-lagrangian schemes, the approximation
also has to be solved iteratively for UMP, but it is simpler than for the two-level time scheme.
The accuracy of the SL scheme depends on the accuracy of the determination of the DP, and on the determination of the value of UDP and the other conserved quantities Q by interpolation from the neighboring points. A linear interpolation between neighboring points results in excessive smoothing, especially for the shortest waves. For this reason cubic interpolation is preferred (Williamson and Laprise, 1998). This is a costly overhead of semilagrangian schemes. Despite the additional costs, in practice this scheme has been found to be accurate and efficient (see the general review of semilagrangian methods by Staniforth and Cote, 1991). Purser and Leslie, 1991, Leslie and Purser, 1995 proposed a “cascade” method that allows efficient high order interpolation between the distorted Lagrangian grid and the regular Eulerian grid. This has allowed them to propose a forward trajectory semilagrangian approach instead of the conventional backward trajectory that we have so far described, which has additional advantages. (See Staniforth and Cote, 1991, Bates et al, 1995, Purser and Leslie, 1996, Williamson and Laprise, 1998 for further details). Combining it with a semi-implicit treatment (next sub-section) of gravity waves, as first suggested by Robert (1982) and Robert et al (1985), increases the efficiency of the semilagrangian scheme. Laprise et al (1997) have documented a "mesoscale compressible community" model, which is non hydrostatic, 3-level semilagrangian, and uses the semi-implicit approach for both the elastic terms (3-dimensional divergence) and the gravity wave terms. As such, it is a flexible model that can be used for many scales.