Nonuniform dz by magpowell · Pull Request #79 · microhh/rte-rrtmgp-cpp

magpowell · 2026-07-23T21:28:50Z

Overview

The Monte Carlo ray tracer required an evenly spaced vertical grid. This PR removes that restriction (enables arbitrary, non-uniform vertical grids) for the forward shortwave and longwave tracers, the backward tracer, and the test drivers.

Changes

The simplest change is we replace the constant dz multiplies within the tau to k_ext and count to flux conversion calculations with a local dz.

The null collision grid requires changes too:

It is now also non-uniform, with null-cell boundaries on every n-th model interface. Placing boundaries on model interfaces keeps k_null >= k_ext by construction; boundaries at fixed heights would cut through layers and could locally underestimate cloud extinction.
At a potential collision we need both k_null and k_ext to decide probabilistically whether a true collision occurred. We know the z of the potential collision, but we need the physical layer that contains it (for k_ext) and the null cell that contains it (for k_null). On an equidistant grid this was a cheap int(z/dz). Replaced with a lookup table binned at the finest layer thickness: one read plus at most one correction, a small fixed cost per photon instead of a search.

Longwave emission positions are now sampled from the interface heights; the emission energy weighting is tau based and needs no change.

Tests

Cabauw, uniform grid, old code vs this PR

When the layer count divides the null cell count the old and new code are bit identical. The old and new code agree to roundoff (heating rates within 1e-7 K/day, per pixel fluxes within 1e-2 W/m2).
When the layer count does not divide the null cell count, the new code places null-cell boundaries at layer edges while the old code put them at fixed heights inside layers. Both bounds are valid, but the two codes take different random walks through the same atmosphere (HR and SW flux agree at Monte Carlo noise).

Same null grid boundaries:

~~Different null grid boundaries:~~
When the layer count does not divide the null cell count, the old code placed boundaries at fixed heights inside layers. With the constant-dz fast path the new code defaults back to those same boundaries on uniform grids, so this case is also bit identical.

RCEMIP on a uniform and stretched grid

We hold the inputs fixed (same cloud scenes) and only vary the grids (mass conserving remapping). Using the offline ray tracer on RCEMIP scenes, we ran 3D RT on a 144 level stretched grid and on a uniform 40 m grid with 806 levels, remapping the scene conservatively between them.
Shortwave and longwave fluxes and heating agree at Monte Carlo noise and the energy budget closes. Mean differences over twelve scenes are below 0.1 W/m2.
How do we know it's MC noise and not a grid effect? A grid effect would not depend on the samples per pixel (spp), but the MC noise does: going from 256 to 4096 spp (16x) shrinks the per pixel RMSE about 4x (direct 2.0 to 0.6 W/m2, diffuse 1.3 to 0.3 W/m2), as expected for MC noise.

SW surface fluxes (one timeslice):

LW heating rates and surface flux (one timeslice):

The spike in the heating rate difference sits at the cloud boundary: the mass conserving remap smears the cloud edge within one 40 m cell and slightly shifts the location of peak heating rate. The difference cancels in the column.

Pixel-wise comparison (12 timeslices):

Performance

~~On uniform grids the radiation step costs 3.8% more.~~ On uniform grids performance is unchanged (uses constant-dz path). Stretched grids need far fewer levels which improves overall performance. The RCEMIP scenes here use 144 levels where the equivalent uniform grid needs 806, a 5.6x reduction in cells, and the ray tracer runs ~1.2x faster on the same scene.

magpowell · 2026-07-23T21:32:43Z

+    __device__
+    inline int height_to_int(const Float z, const Float* z_lev, const int* z_lut, const Float lut_dz, const int lut_size, const int ntot_max)
+    {
+        int k = z_lut[min(max(static_cast<int>(z / lut_dz), 0), lut_size-1)];


here we are finding which physical layer the possible collision occurs. k is an index for z_lev such that z_lev[k] <= z < z_lev[k+1]

magpowell · 2026-07-23T21:37:39Z

+        while (k > 0 && z < z_lev[k])
+            --k;
+        while (k < ntot_max-1 && z >= z_lev[k+1])
+            ++k;


this is a correction, a bin in the lookup table can overlap with at most two layers. the LUT stores the layer of the bin's lower edge so it's a possible upward correction.

magpowell · 2026-07-23T21:37:51Z

+    inline int height_to_int(const Float z, const Float* z_lev, const int* z_lut, const Float lut_dz, const int lut_size, const int ntot_max)
+    {
+        int k = z_lut[min(max(static_cast<int>(z / lut_dz), 0), lut_size-1)];
+        while (k > 0 && z < z_lev[k])


for floating point safety

magpowell · 2026-07-23T21:49:16Z

+            dz_min = std::min(dz_min, z_lev({k+2}) - z_lev({k+1}));
+
+        const Float zsize = z_lev({nz+1});
+        const int lut_size = static_cast<int>(std::ceil(zsize/dz_min));


we use the thinnest physical grid spacing (dz_min) to set the table vertical resolution

magpowell · 2026-07-23T22:15:03Z

                const Raytracer_definitions::Vector<int> grid_cells,
                const Raytracer_definitions::Vector<Float> grid_d,
                const Raytracer_definitions::Vector<int> kn_grid,
+                const Array_gpu<Float,1>& z_lev,


z_lev: physical grid interface heights
kn_z_lev: null-cell boundary heights (kn_ktot+1 vals), a subset of z_lev, so each null cell holds a whole number of layers (~ktot/kn_ktot each).
z_lut: height bin (k) to physical layer index table
kn_z_lut: height bin (k) to null cell index table

most of the code changes are just carrying these around. the ray tracer kernels can no longer derive k_ext and k_null from a single dz, each lookup needs the two tables. Once built, it's fairly cheap: per photon cost is one read, one compare.

magpowell · 2026-07-23T22:37:45Z

-            const Float kext_tot = tau_tot[idx] / grid_d.z;
-            const Float kext_cld = tau_cld[idx] / grid_d.z;
-            const Float kext_aer = tau_aer[idx] / grid_d.z;
+            const Float dz = z_lev[iz+1] - z_lev[iz];


switching to a local dz

magpowell · 2026-07-23T22:39:05Z

-            const Float kext_tot = tau_tot[idx] / grid_d.z;
-            const Float kext_cld = tau_cld[idx] / grid_d.z;
-            const Float kext_aer = tau_aer[idx] / grid_d.z;
+            const Float dz = z_lev[iz+1] - z_lev[iz];


magpowell · 2026-07-23T22:40:58Z


                const int km = k - 1;
-                photon.position.z = (km + rng()) * grid_d.z;
+                photon.position.z = z_lev[km] + rng() * (z_lev[km+1] - z_lev[km]);


again using local not fixed dz

magpowell · 2026-07-23T23:00:37Z

        Status::print_message("Storing the bw raytracer output.");

-        auto nc_radiance = output_nc.add_variable<Float>("radiance" , {"ny", "nx"});
+        auto nc_radiance = output_nc.add_variable<Float>("radiance" , {"py", "px"});


unrelated to this PR, but I think this is a bug

radiance is camera output, so it is sized by the camera pixel counts (py, px), not the domain (ny, nx). The old dims only work when the two happen to match. It breaks for a camera whose pixel count differs from the domain size. test_rte_rrtmgp_bw already writes radiance with py/px, so this just makes the two drivers consistent.

magpowell · 2026-07-23T23:02:02Z

+
+        Array<int,1> z_lut({lut_size});
+        Array<int,1> kn_z_lut({lut_size});
+        for (int i=0; i<lut_size; ++i)


Example: layers of 20, 25, 60, 200 m. LUT spacing 20 m.

Bin 0 (0 - 20 m) corresponds to layer 0 (0 - 20 m), LUT stores k = 0.
Bin 1 (20–40 m) lies inside layer 1 (20- 45 m), LUT stores k = 1
Bin 2 (40–60 m) includes the layers 1 and 2. LUT stores the layer of its lower edge: k = 1. A possible collision that falls within bin 2 would be eligible for a correction (++k) in the height_to_int func if z >= z_lev[k+1]
Bin 3 covers 60–80 m, which lies inside layer 2 (45–105 m), LUT stores k =2.
and so on

The lw/sw solve_gpu built a uniform vertical grid from grid_d.z instead of the input interface heights, and the lw plane-parallel fallback divided by a scalar dz; both broke the 3D absorption output on non-uniform grids.

magpowell · 2026-07-25T00:35:55Z

+        cudaMemcpy(&max_kext_gas, max_kext_gas_g, sizeof(Float), cudaMemcpyDeviceToHost);

-        const Float lowest_gas_mean_free_path = grid_d.z / max_tau_gas;
+        const Float lowest_gas_mean_free_path = Float(1.) / max_kext_gas;


Before the mean-free-path criterion divided a fixed spacing by the maximum layer tau, which was fine for uniform grid but for stretched grid biases the estimate low by dz/dz_min. Now we use the maximum of tau/dz, the largest extinction coefficient in the domain, which is grid-independent.

magpowell · 2026-07-28T22:42:40Z

Update: I added the constant-dz fast path Menno suggested. The kernels take a dz_constant template argument like independent_column. The path is chosen at runtime. The solver checks whether the interface heights are equally spaced, and if so it runs the original equidistant arithmetic, including the original null-cell boundaries at uniform spacing for any kn_ktot.

This simplifies what I wrote in the Cabauw test section. On uniform grids the new code now always uses the old null-cell boundaries, so the non-divisible case is also bit identical and the second test figure no longer applies.

magpowell · 2026-07-28T22:44:23Z

    const int mie_phase_table_size = mie_phase_ang.size();

+    // Constant dz implies the old uniform null-cell walls, valid for any kn_grid.
+    const bool dz_constant = vertical_spacing_is_constant(z_lev, grid_cells.z);


check if there is a uniform grid to set the argument passed to the ray tracer kernels

MennoVeerman

Code works like a charm, thanks for adding support for nonuniform dz grids!

Margaret Powell added 2 commits July 22, 2026 22:52

Support non-uniform vertical grids in the sw, lw and bw ray tracers

0afe164

Fix radiance output dimensions in test_rt_lite

348e3a4

magpowell commented Jul 23, 2026

View reviewed changes

Comment thread src_test/test_rt_lite.cu

magpowell commented Jul 23, 2026

View reviewed changes

magpowell force-pushed the nonuniform-dz branch from d82e854 to 1e66ab2 Compare July 25, 2026 00:10

Use input z_lev in the forward solvers of test_rte_rrtmgp_rt

1e71055

The lw/sw solve_gpu built a uniform vertical grid from grid_d.z instead of the input interface heights, and the lw plane-parallel fallback divided by a scalar dz; both broke the 3D absorption output on non-uniform grids.

magpowell force-pushed the nonuniform-dz branch 2 times, most recently from 2a5a777 to 3c06c15 Compare July 25, 2026 00:24

Base the LW ray-tracer selection on extinction, not layer optical depth

04913ac

magpowell force-pushed the nonuniform-dz branch from 3c06c15 to 04913ac Compare July 25, 2026 00:26

magpowell commented Jul 25, 2026

View reviewed changes

magpowell marked this pull request as ready for review July 27, 2026 19:58

Add a constant-dz fast path to the ray tracer kernels

a19c1cd

magpowell force-pushed the nonuniform-dz branch from 566c1be to a19c1cd Compare July 28, 2026 22:36

magpowell commented Jul 28, 2026

View reviewed changes

MennoVeerman approved these changes Jul 29, 2026

View reviewed changes

Uh oh!

Conversation

magpowell commented Jul 23, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Overview

Changes

Tests

Cabauw, uniform grid, old code vs this PR

RCEMIP on a uniform and stretched grid

Performance

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

magpowell Jul 23, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

magpowell Jul 23, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

magpowell Jul 25, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

magpowell commented Jul 28, 2026

Uh oh!

Choose a reason for hiding this comment

Uh oh!

MennoVeerman left a comment

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

magpowell commented Jul 23, 2026 •

edited

Loading

magpowell Jul 23, 2026 •

edited

Loading

magpowell Jul 23, 2026 •

edited

Loading

magpowell Jul 25, 2026 •

edited

Loading