Sampling NWP grids on the GPU

Every GPU-based visualization algorithm in Met.3D — sections, isosurfaces, direct volume rendering, normal curves — needs the same primitive: give me the value of this data field at an arbitrary position in longitude/latitude/pressure space. How expensive that primitive is depends entirely on the vertical level type of the grid, and it is the main reason why pressure level data renders considerably faster than model level data.

This page describes what is uploaded to the GPU, how a field is sampled, and how empty-space skipping accelerates the raycaster.

Source code: src/glsl/volume_sample_utils.glsl, src/glsl/volume_pressure_utils.glsl, src/glsl/volume_hybrid_utils.glsl, src/glsl/volume_auxiliarypressure_utils.glsl, src/glsl/volume_acceleration_utils.glsl, src/data/structuredgrid.cpp.

Sampling of data fields in GPU shaders.

Fig. 55 Sampling of data fields in GPU shaders. (a) Since the elevation of the model levels (green lines, level indices k in green) changes with surface pressure (black curve at the bottom), the data value for a given pressure value p can be located at different levels in the grid (the red line marks p = 600 hPa). (b) For each vertex of a horizontal section, the model levels k and k+1 are found by binary search; the scalar value is linearly interpolated in ln(p) between these two levels. (c) Pressure level grids are rectilinear (left), which allows trilinear hardware interpolation between the grid points surrounding a sample position (red dot). For model level grids (right), the sample position can lie between different model levels k for two adjacent grid columns, which prohibits hardware interpolation. (d) Example of how the surface orography affects the vertical displacement of the grid points in a vertical section. Figure reproduced from Rautenhaus et al. (2015, doi:10.5194/gmd-8-2329-2015), Figs. 9 and 11, © Author(s) 2015. CC Attribution 3.0 License.

What is uploaded to the GPU

The unit of upload is one data field: one variable, one ensemble member, one time step (see the smallest data entity assumption). For that field, the following GPU resources are created on demand by MStructuredGrid and registered with MGLResourcesManager:

Resource

Content

3-D data texture

The scalar values, nlon x nlat x nlev, single precision.

1-D axis texture

Longitude and latitude values, and for pressure level grids also the pressure of each level.

Surface pressure texture (2-D)

For hybrid sigma-pressure grids: the companion surface pressure field.

Hybrid coefficients

For hybrid sigma-pressure grids: the \(a_k\) and \(b_k\) arrays.

3-D pressure texture

For auxiliary pressure grids only: the full 3-D pressure field.

Note what is not uploaded for hybrid sigma-pressure grids: a 3-D field of pressure values. Since \(p(k,j,i) = a_k + b_k \cdot p_\text{sfc}(j,i)\), the pressure of any grid point can be computed in the shader from the 2-D surface pressure field and the coefficient arrays, which saves an entire 3-D texture per field. Auxiliary pressure grids cannot do this and pay the extra texture.

Sampling pressure level grids

A pressure level grid is rectilinear: the pressure of a grid point depends only on its level index, not on its horizontal position. The whole field can therefore be sampled with a single texture fetch, letting the GPU perform trilinear interpolation in hardware.

The only complication is that pressure levels are usually not equidistant, so the mapping from a pressure value to the vertical texture coordinate is not a linear scaling. Met.3D solves this with a lookup table:

  • MRegularLonLatStructuredPressureGrid::getPressureTexCoordTexture1D() builds a 1-D texture with 2048 entries, regularly spaced in ln(p) between the bottom and the top level of the grid.

  • Each entry holds the vertical texture coordinate belonging to that ln(p), obtained by linear interpolation in ln(p) between the known level pressures.

  • In the shader, one fetch into this LUT converts ln(p) into the texture coordinate, and one fetch into the 3-D texture returns the interpolated value.

Interpolating in ln(p) rather than in p matches the roughly exponential vertical structure of the atmosphere and is what the ECMWF FULLPOS interpolation routines do as well.

Sampling hybrid sigma-pressure grids

Model level grids are not rectilinear. Because the level elevation follows the surface pressure, the continuous level index that corresponds to a given pressure differs between neighbouring grid columns — panel (c) of the figure above. In the worst case, a sample position lies between different model levels in each of its four surrounding columns. Hardware interpolation requires a single common texture coordinate for all of them and is therefore unusable.

The sampling is consequently split:

  1. Four vertical interpolations, one per surrounding grid column, each yielding the value at the requested pressure in that column.

  2. One bilinear horizontal interpolation between those four values, implemented in software in the shader.

A naive implementation would run a binary search per column, as the section shaders do. The raycaster instead extends the LUT idea into two dimensions, which is roughly twice as fast:

  • MLonLatHybridSigmaPressureGrid::getPressureTexCoordTexture2D() builds a 2-D texture whose rows correspond to discrete surface pressure values and whose columns are regularly spaced in ln(p).

  • The current table covers 1200 surface pressure values from 1050 hPa down to 450 hPa in steps of 0.5 hPa, with 2048 vertical entries each (the first two entries of each row store the table boundaries). At 4 bytes per entry this is roughly 10 MB of GPU memory.

  • Bilinear hardware interpolation of this table interpolates in both surface pressure and ln(p) at once, so a single fetch turns (psfc, ln p) into a vertical texture coordinate, and the actual value lookup can again use hardware interpolation in the vertical.

  • The table depends only on the grid, not on the variable, and is therefore shared by all fields on the same grid.

That comes to five texture accesses per column. Each row of the table is regularly spaced in ln(p) between that column’s own top and bottom pressure, both of which move with the surface pressure, so ln_p has to be normalised before it can be used as a row coordinate. The two constants needed for that (the top of the column and its extent in ln(p)) are stored in the first two entries of the row.

Sampling auxiliary pressure grids

Auxiliary pressure grids store the pressure of every grid point in a separate 3-D field, so no formula and no lookup table can shortcut the vertical search: the pressure profile of a column is only known by reading it. Sampling (volume_auxiliarypressure_utils.glsl) therefore falls back to the approach the hybrid path avoids — a binary search per grid column, with unfiltered point fetches throughout:

  1. sampleAuxiliaryPressureColumnAtP() fetches the bottom entry of the column’s pressure profile to check that p is inside the data volume at all.

  2. It then binary-searches the pressure texture for the two levels enclosing p, one texelFetch per iteration.

  3. It fetches the two enclosing pressures, interpolates in ln(p), and fetches the two enclosing scalar values.

  4. sampleAuxiliaryPressureVolumeAtPosPressure() repeats that for all four surrounding columns and blends the results bilinearly in software.

None of this can use hardware filtering, because the vertical neighbours are found by search rather than addressed by a texture coordinate.

Table 1 Texture accesses for one sample

Level type

Fetches per sample

Composition

Pressure levels

~2

One 1-D LUT fetch, one filtered 3-D fetch; the GPU performs the trilinear interpolation.

Hybrid sigma-pressure

~20

Per column: one surface pressure fetch, three fetches into the 2-D LUT (two of them for the row’s normalisation constants), one filtered 3-D fetch. Four columns, then software bilinear.

Auxiliary pressure

~48

Per column: log2(nlev) binary search steps plus three further pressure fetches and two scalar fetches (~12 for 137 levels). Four columns, then software bilinear.

Note

Sampling a hybrid sigma-pressure grid costs several times as much as sampling a pressure level grid, and auxiliary pressure grids are the most expensive of the three: their cost grows with log2(nlev), and the binary search diverges between neighbouring fragments, which GPUs handle badly. If raycasting performance matters more than staying on native levels, regridding to pressure levels is the most effective single change.

Empty-space skipping

Most of a forecast volume is irrelevant for a given isosurface. Met.3D therefore maintains a min/max acceleration structure per data field:

  • MStructuredGrid::getMinMaxAccelTexture3D() divides the longitude/latitude/pressure extent of the field into a uniform grid of 32 x 32 x 32 bricks and stores the minimum and maximum data value of each brick in an RG 3-D texture. It is computed once per field, on demand, and cached like any other GPU resource.

  • While traversing a ray, the shader determines the brick it is currently in (findCurrentAccelerationStructureBrick()) and tests whether the requested isovalue can occur inside it (canSkipBrick()). If not, the ray jumps straight to the brick exit.

  • In direct volume rendering the same structure is used against the transfer function: if the transfer function is fully transparent below or above a certain value, bricks whose entire value range falls into the transparent part contribute nothing and are skipped as well.

Effectiveness depends strongly on the data and on the chosen isovalue. A field whose values are spread evenly over the domain gives the structure little to skip, and rendering times vary accordingly with camera position and isovalue.