updates and bug fixes to bw implementation

include_rt/raytracer_functions.h

@@ -87,104 +87,62 @@ namespace Raytracer_functions
     }
 
     __device__
-    inline Float mie_sample_angle(const Float* mie_cdf, const Float* mie_lut, const Float random_number, const Float r_eff, const int n_mie)
+    inline int find_index(const float* mie_cdf, const int size, const float random_number)
     {
-        // interpolation over effective radius. Currently, r_eff should range between 2.5 and 21.5 (similar to RRTMGP) OR
-        // be exactly 100 micrometer for optical effects such as rainbows
-        const int r_idx = (r_eff == Float(100.)) ? 20 : min(max(int(r_eff-2.5), 0), 18);
-        const Float r_rest = fmod(r_eff-Float(2.5),Float(1.));
+        int left = 0;
+        int right = size - 1;
 
-        int i = 0;
-        while (random_number < mie_cdf[i])
-        {
-            ++i;
-        }
+        while (left < right) {
+            int mid = left + (right - left) / 2;
 
-        // sampled scattering angle
-        Float ang;
-        if (r_idx < 20)
-        {
-            if (i==0)
-            {
-                const Float ang_lwr = mie_lut[r_idx*n_mie]*(1-r_rest);
-                const Float ang_upr = mie_lut[(r_idx+1)*n_mie]*r_rest;
-                ang = ang_lwr + ang_upr;
-            }
-            else
-            {
-                const int midx_lwr = r_idx*n_mie;
-                const int midx_upr = (r_idx+1)*n_mie;
-                const Float dr = abs(mie_cdf[i] - mie_cdf[i-1]);
-
-                const Float ang_lwr = (abs(random_number - mie_cdf[i])*mie_lut[(i-1)+midx_lwr] + abs(mie_cdf[i-1]-random_number)*mie_lut[i+midx_lwr]) / dr;
-                const Float ang_upr = (abs(random_number - mie_cdf[i])*mie_lut[(i-1)+midx_upr] + abs(mie_cdf[i-1]-random_number)*mie_lut[i+midx_upr]) / dr;
-                ang = ang_lwr * (1-r_rest) + ang_upr * r_rest;
+            if (random_number >= mie_cdf[mid]) {
+                right = mid;
+            } else {
+                left = mid + 1;
             }
         }
-        else
-        {
-            if (i==0)
-            {
-                ang = mie_lut[r_idx*n_mie];
-            }
-            else
-            {
-                const int midx = r_idx*n_mie;
-                const Float dr = abs(mie_cdf[i] - mie_cdf[i-1]);
 
-                ang = (abs(random_number - mie_cdf[i])*mie_lut[(i-1)+midx] + abs(mie_cdf[i-1]-random_number)*mie_lut[i+midx]) / dr;
-            }
-        }
+        return left - 1;
+    }
+
+    __device__
+    inline Float mie_sample_angle(const Float* mie_cdf, const Float* mie_lut, const Float random_number, const Float r_eff, const int n_mie)
+    {
+        // interpolation over effective radius. Currently, r_eff should range between 2.5 and 21.5 (similar to RRTMGP)
+        const int r_idx = min(max(int(r_eff-2.5), 0), 18);
+        const Float r_rest = fmod(r_eff-Float(2.5),Float(1.));
+
+        const int i = min(max(0, find_index(mie_cdf, n_mie, random_number)), n_mie - 2);
+
+        const int midx_lwr = r_idx*n_mie;
+        const int midx_upr = (r_idx+1)*n_mie;
+        const Float dr = abs(mie_cdf[i+1] - mie_cdf[i]);
+
+        const Float ang_lwr = (abs(random_number - mie_cdf[i+1])*mie_lut[(i)+midx_lwr] + abs(mie_cdf[i]-random_number)*mie_lut[i+midx_lwr+1]) / dr;
+        const Float ang_upr = (abs(random_number - mie_cdf[i+1])*mie_lut[(i)+midx_upr] + abs(mie_cdf[i]-random_number)*mie_lut[i+midx_upr+1]) / dr;
+        const Float ang = ang_lwr * (1-r_rest) + ang_upr * r_rest;
         return ang;
     }
 
     __device__
     inline Float mie_interpolate_phase_table(const Float* mie_phase, const Float* mie_lut, const Float scat_ang, const Float r_eff, const int n_mie)
     {
-        // interpolation over effective radius. Currently, r_eff should range between 2.5 and 21.5 (similar to RRTMGP) OR
-        // be exactly 100 micrometer for optical effects such as rainbows
-        const int r_idx = (r_eff == Float(100.)) ? 20 : min(max(int(r_eff-2.5), 0), 18);
+        // interpolation over effective radius. Currently, r_eff should range between 2.5 and 21.5 (similar to RRTMGP)
+        const int r_idx = min(max(int(r_eff-2.5), 0), 18);
         const Float r_rest = fmod(r_eff-Float(2.5),Float(1.));
 
         // interpolation between 1800 equally spaced scattering angles between 0 and PI (both inclusive).
-        const Float d_pi = Float(1.74629942e-03);
-        const int i = min(max(0, int(1800-(scat_ang/d_pi+1))), 1798);
+        constexpr Float d_pi = Float(1.74629942e-03);
+        const int i = min(max(0, int(scat_ang/d_pi)), 1798);
 
-        // probability (of scattering at angle scat_ang)
-        Float prob;
-        if (r_idx < 20)
-        {
-            if (i==0)
-            {
-                const Float prob_lwr = mie_lut[r_idx*n_mie]*(1-r_rest);
-                const Float prob_upr = mie_lut[(r_idx+1)*n_mie]*r_rest;
-                prob = prob_lwr + prob_upr;
-            }
-            else
-            {
-                const int midx_lwr = r_idx*n_mie;
-                const int midx_upr = (r_idx+1)*n_mie;
-                const Float dr = abs(mie_phase[i] - mie_phase[i-1]);
-
-                const Float prob_lwr = (abs(scat_ang - mie_phase[i])*mie_lut[(i-1)+midx_lwr] + abs(mie_phase[i-1]-scat_ang)*mie_lut[i+midx_lwr]) / dr;
-                const Float prob_upr = (abs(scat_ang - mie_phase[i])*mie_lut[(i-1)+midx_upr] + abs(mie_phase[i-1]-scat_ang)*mie_lut[i+midx_upr]) / dr;
-                prob = prob_lwr * (1-r_rest) + prob_upr * r_rest;
-            }
-        }
-        else
-        {
-            if (i==0)
-            {
-                prob = mie_lut[r_idx*n_mie];
-            }
-            else
-            {
-                const int midx = r_idx*n_mie;
-                const Float dr = abs(mie_phase[i] - mie_phase[i-1]);
+        const int midx_lwr = r_idx*n_mie;
+        const int midx_upr = (r_idx+1)*n_mie;
+        const Float dr = abs(mie_phase[i+1] - mie_phase[i]);
+
+        const Float prob_lwr = (abs(scat_ang - mie_phase[i+1])*mie_lut[(i)+midx_lwr] + abs(mie_phase[i]-scat_ang)*mie_lut[i+1+midx_lwr]) / dr;
+        const Float prob_upr = (abs(scat_ang - mie_phase[i+1])*mie_lut[(i)+midx_upr] + abs(mie_phase[i]-scat_ang)*mie_lut[i+1+midx_upr]) / dr;
+        const Float prob = prob_lwr * (1-r_rest) + prob_upr * r_rest;
 
-                prob = (abs(scat_ang - mie_phase[i])*mie_lut[(i-1)+midx] + abs(mie_phase[i-1]-scat_ang)*mie_lut[i+midx]) / dr;
-            }
-        }
         return prob;
     }
 

include_rt_kernels/raytracer_kernels_bw.h

@@ -17,8 +17,8 @@ using namespace Raytracer_functions;
 
 
 #ifdef RTE_USE_SP
-constexpr int bw_kernel_block= 512;
-constexpr int bw_kernel_grid = 1024;
+constexpr int bw_kernel_block= 128;
+constexpr int bw_kernel_grid = 2048;
 #else
 constexpr int bw_kernel_block = 256;
 constexpr int bw_kernel_grid = 256;

python/set_virtual_camera.py

@@ -63,7 +63,7 @@
     cam["yaw"][:]   = 0
     cam["pitch"][:] = -90
     cam["roll"][:]  = 0
-    cam["cam_type"][:]= 1
+    cam["cam_type"][:]= 0
     cam["fov"][:] = 80
     cam["px"][:] = 0
     cam["py"][:] = 0
@@ -76,7 +76,7 @@
     cam["yaw"][:] = 0
     cam["pitch"][:] = 0
     cam["roll"][:]  = 0
-    cam["cam_type"][:]= 0
+    cam["cam_type"][:]= 1
     cam["fov"][:] = 80
     cam["px"][:] = 0.
     cam["py"][:] = 0.

src_kernels_cuda_rt/raytracer_kernels_bw.cu

@@ -280,7 +280,8 @@ namespace
             const Vector<Float>& normal,
             const Phase_kind kind)
     {
-        const Float cos_angle = dot(photon.direction, sun_direction);
+        const Float cos_angle = min(max(Float(-1.), dot(photon.direction, sun_direction)), Float(1.));
+
         if (kind == Phase_kind::HG)
         {
             return henyey_phase(g, cos_angle) * sun_solid_angle;

src_test/test_rte_rrtmgp_bw.cu

@@ -401,7 +401,7 @@ void solve_radiation(int argc, char** argv)
             rel = std::move(input_nc.get_variable<Float>("rel", {n_lay, n_col_y, n_col_x}));
         }
 
-        if(switch_ice_cloud_optics)
+        if (switch_ice_cloud_optics)
         {
             iwp.set_dims({n_col, n_lay});
             iwp = std::move(input_nc.get_variable<Float>("iwp", {n_lay, n_col_y, n_col_x}));

src_test/test_rte_rrtmgp_rt.cu

@@ -483,7 +483,6 @@ void solve_radiation(int argc, char** argv)
             const Array<Float,2>& gas = aerosol_concs.get_vmr(aerosol_name);
             if (gas.size() > 1) {
                 if (gas.get_dims()[0] == 1){
-                    printf("----");
                     aerosol_concs.set_vmr(aerosol_name, aerosol_concs.get_vmr(aerosol_name).subset({ {{1,n_col}, {1, n_lay}}} ));
                 }
             }