Skip to content

7_Chebyshev_grids

Jump to bottom
BenMql edited this page Jul 2, 2024 · 6 revisions

For setting properly the output of the code, but also when postprocessing data, it is necessary to have some notions on the collocation grids used in the code. While the collocation grid along the horizontal is a classic, equispaced, Fourier collocation grid, the Chebyshev collocation grid used along the vertical has subtleties that we detail here. Since v1.3.0, the user can choose between two grids. These grids are completely equivalent from the viewpoint of computations. They do not play any role for enforcing boundary conditions, which are imposed in spectral space by basis recombination. The two kinds of grid differ by the presence of endpoints. Crudely: the Gauss--Chebyshev grid contains only inner-points, with a (ever so slightly) better resolution in the bulk, whereas the Gauss--Lobatto grid contains the end-points, which offers the possibility to export the value of the field along the top and bottom bounding surfaces.

The Gauss--Chebyshev grid, or "inner points" grid

This is Coral's default mode. When the user requests NC Chebyshev polynomials along z as parameter #03 in the coral.parameters.in file, (linear and quadratic) fields are computed on a grid composed of NG = 3 NC / 2 points. (The reason is to avoid modes aliasing around Nyquist frequency [see e.g. Boyd's book on spectral methods].) For a computational domain centered around zCenter and with a gap zGap (parameters #08 and #07, respectively), the i-th point of the grid (starting at i=1, in Fortran spirit) has altitude:

z_i = zCenter + zGap/2 * cos( [2*i-1]*pi /[2 NG] )    for i = 1..NG

For Python which uses 0-based indexing, a possible definition is:

zGaussChebyshev = zCenter + zGap/2. * np.cos( (2.*np.arange(NG)+1)*np.pi / (2.*NG) )

In particular, notice that:

  • the grid is not equispaced
  • the grid in inverted: i=1 is the highest point, i=NG is the lowest point
  • the grid consist of inner points only: i=1 is slightly below the top and i=NG is slightly above the bottom. In particular, there are no endpoints, where boundary conditions are applied. In case of Dirichlet boundary conditions, the last bullet provides the clear advantage of not computing the fields where we know in advance that they cancel. Hence we do not waste grid-points where the are not needed because boundary conditions are enforced. For Neumann boundary conditions, this grid does not allow to export the values of the field along the bounding surfaces

The Gauss--Lobatto grid, or "end points" grid

Coral can be run in this mode with the --grid-with-endpoints flag, or alternatively the --gauss-lobatto-grid flag. When the user requests NC Chebyshev polynomials along z as parameter #03 in the coral.parameters.in file, alias-free fields are obtained on a grid with NG = 3 NC / 2 + 1 points. For better performance, NC should be choosen with small prime factors decomposition (2,3,5, or 7, to stay clear from trouble). For a computational domain centered around zCenter and with a gap zGap (parameters #08 and #07, respectively), the i-th point of the grid has altitude, using 1-based indexing:

z_i = zCenter + zGap/2 * cos( [i-1]*pi /[NG-1] )    for i = 1..NG

A possible Python (0-based indexing) definition is:

zGaussLobatto = zCenter + zGap/2. * np.cos( (np.arange(NG))*np.pi / (NG-1.))

In particular, notice that:

  • the grid is not equispaced, and inverted
  • the grid has the end points: i=1 and i=NG correspond to the top and the bottom of the interval, respectively. In case of Dirichlet boundary conditions, the first and last points are wasteful, as the fields vanish. Because two points lie at the end of the interval, only NG-2 are in the bulk. The advantage of the Gauss-Lobatto grid is that is offers the possibility to export the field at both boundaries.

For instance, if benchmark 01 is run with NC=24 Chebyshev polynomials, the mean flows uFull and vFull at the (stress-free) top and bottom surfaces can be exported by adding the following lines to the file coral.timeseries:

>>Horizontal_avg :: uFull ;; Position= 1 <<
>>Horizontal_avg :: uFull ;; Position= 37 <<
>>Horizontal_avg :: vFull ;; Position= 1 <<
>>Horizontal_avg :: vFull ;; Position= 37 <<

Remarks

The CheckPoints obtained with one grid cannot be used for restarting a run with the other grid. However, they can be interpolated, either naively or by using DCT-based methods.

For reference, transforms between physical and spectral space use DCT-II and DCT-III for the inner, Gauss-Chebyshev grid.

For the Gauss-Lobatto (endpoints) grid, the DCT-I is used for both forward and backward transforms.

Clone this wiki locally