Merge pull request #39 from MennoVeerman/main

some changes to rt_lite and fixing cloud optics bug

Author: Chiil

include_rt_kernels/raytracer_kernels_bw.h

@@ -23,7 +23,6 @@ constexpr int bw_kernel_grid = 1024;
 constexpr int bw_kernel_block = 256;
 constexpr int bw_kernel_grid = 256;
 #endif
-constexpr Float k_null_gas_min = Float(1.e-3);
 
 // sun has a half angle of .266 degrees
 constexpr Float cos_half_angle = Float(0.9999891776066407); // cos(half_angle);
@@ -60,24 +59,21 @@ struct Camera
         const Float yaw = yaw_deg / Float(180.) * M_PI;
         const Float pitch = pitch_deg / Float(180.) * M_PI;
         const Float roll = roll_deg / Float(180.) * M_PI;
-        mx = {cos(yaw)*sin(pitch), (cos(yaw)*cos(pitch)*sin(roll)-sin(yaw)*cos(roll)), (cos(yaw)*cos(pitch)*cos(roll)+sin(yaw)*sin(roll))};
-        my = {sin(yaw)*sin(pitch), (sin(yaw)*cos(pitch)*sin(roll)+cos(yaw)*cos(roll)), (sin(yaw)*cos(pitch)*cos(roll)-cos(yaw)*sin(roll))};
-        mz = {-cos(pitch), sin(pitch)*sin(roll), sin(pitch)*cos(roll)};
+        mx = {cos(yaw)*cos(pitch), (cos(yaw)*sin(pitch)*sin(roll)-sin(yaw)*cos(roll)), (cos(yaw)*sin(pitch)*cos(roll)+sin(yaw)*sin(roll))};
+        my = {sin(yaw)*cos(pitch), (sin(yaw)*sin(pitch)*sin(roll)+cos(yaw)*cos(roll)), (sin(yaw)*sin(pitch)*cos(roll)-cos(yaw)*sin(roll))};
+        mz = {-sin(pitch), cos(pitch)*sin(roll), cos(pitch)*cos(roll)};
     }
 
     void setup_normal_camera(const Camera camera)
     {
         if (!fisheye)
         {
-            const Vector<Float> dir_tmp = {0, 0, 1};
-            const Vector<Float> dir_cam = normalize(Vector<Float>({dot(camera.mx,dir_tmp), dot(camera.my,dir_tmp), dot(camera.mz,dir_tmp)*Float(-1)}));
-            Vector<Float> dir_up;
-            if ( (int(dir_cam.z)==1) || (int(dir_cam.z)==-1) )
-                dir_up = {1, 0, 0};
-            else
-                dir_up = {0, 0, 1};
-
-            cam_width = normalize(cross(dir_cam, dir_up));
+            const Vector<Float> dir_tmp = {1, 0, 0};
+            const Vector<Float> dir_cam = normalize(Vector<Float>({dot(camera.mx,dir_tmp), dot(camera.my,dir_tmp), dot(camera.mz,dir_tmp)}));
+            const Vector<Float> dir_up_tmp = {0, 0, 1};
+            const Vector<Float> dir_up = normalize(Vector<Float>({dot(camera.mx,dir_up_tmp), dot(camera.my,dir_up_tmp), dot(camera.mz,dir_up_tmp)}));
+
+            cam_width = Float(-1) * normalize(cross(dir_cam, dir_up));
             cam_height = normalize(cross(dir_cam, cam_width));
             cam_depth = dir_cam / tan(fov/Float(180)*M_PI/Float(2.));
 
@@ -91,17 +87,18 @@ struct Camera
     // size of output arrays, either number of horizontal and vertical pixels, or number of zenith/azimuth angles of fisheye lens. Default to 1024
     int ny = 1024;
     int nx = 1024;
+    Int npix;
 };
 
 
 __global__
 void ray_tracer_kernel_bw(
         const int igpt,
-        const int photons_per_pixel,
+        const Int photons_per_pixel,
         const Grid_knull* __restrict__ k_null_grid,
         Float* __restrict__ camera_count,
         Float* __restrict__ camera_shot,
-        int* __restrict__ counter,
+        Int* __restrict__ counter,
         const Float* __restrict__ k_ext, const Optics_scat* __restrict__ scat_asy,
         const Float* __restrict__ k_ext_bg, const Optics_scat* __restrict__ scat_asy_bg,
         const Float* __restrict__ z_lev_bg,

python/set_virtual_camera.py

@@ -11,6 +11,7 @@
 import netCDF4 as nc
 import numpy as np
 import argparse
+import os
 
 camera_variables = {
     "yaw": "Horizontal direction of camera, 0: east, 90: north, 180/-180: weast, -90/270: south",
@@ -39,7 +40,10 @@
 args = vars(parser.parse_args())
 
 # open netcdf file
-ncf = nc.Dataset(args['name'],"r+")
+if os.path.isfile(args['name']):
+    ncf = nc.Dataset(args['name'],"r+")
+else:
+    print("file does not exist")
 
 # add group if it does not exsist yet
 if not "camera-settings" in ncf.groups:
@@ -94,7 +98,7 @@
 print("Camera settings:")
 for v in camera_variables.keys():
     print("{:6}{:>8}".format(v, str(cam[v][:])))
-print("{:6}{:>8}".format("sza", str(np.round(np.rad2deg(np.arccos(ncf['mu0'][0,0])),1))))
-print("{:6}{:>8}".format("azi", str(np.round(np.rad2deg(ncf['azi'][0,0]),1))))
+print("{:6}{:>8}".format("sza", str(np.round(np.rad2deg(np.arccos(ncf['mu0'][:].flatten()[0])),1))))
+print("{:6}{:>8}".format("azi", str(np.round(np.rad2deg(ncf['azi'][:].flatten()[0]),1))))
 
 ncf.close()

src_cuda_rt/Cloud_optics_rt.cu

@@ -43,11 +43,10 @@ namespace
 
             if (mask[idx])
             {
-                const Float reff= min(Float(21.5), re[idx]);
-                const int index = min(int((reff - offset) / step_size) + 1, nsteps-1) - 1;
+                const int index = min(int((re[idx] - offset) / step_size) + 1, nsteps-1) - 1;
 
                 const int idx_ib = index + ibnd*nsteps;
-                const Float fint = (reff - offset) /step_size - (index);
+                const Float fint = (re[idx] - offset) /step_size - (index);
 
                 const Float tau_local = cwp[idx] *
                                      (tau_table[idx_ib] + fint * (tau_table[idx_ib+1] - tau_table[idx_ib]));

src_cuda_rt/Raytracer_bw.cu

@@ -456,11 +456,12 @@ void Raytracer_bw::trace_rays(
 
     Array_gpu<Float,2> camera_count({camera.nx, camera.ny});
     Array_gpu<Float,2> shot_count({camera.nx, camera.ny});
-    Array_gpu<int,1> counter({1});
+    Array_gpu<Int,1> counter({1});
 
     Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(camera.nx, camera.ny, camera_count.ptr());
     Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(camera.nx, camera.ny, shot_count.ptr());
-    Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(1, counter.ptr());
+    
+    counter.fill(Int(0));
 
     // domain sizes
     const Vector<Float> grid_size = grid_d * grid_cells;
@@ -604,11 +605,11 @@ void Raytracer_bw::trace_rays_bb(
 
     Array_gpu<Float,2> camera_count({camera.nx, camera.ny});
     Array_gpu<Float,2> shot_count({camera.nx, camera.ny});
-    Array_gpu<int,1> counter({1});
+    Array_gpu<Int,1> counter({1});
 
     Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(camera.nx, camera.ny, camera_count.ptr());
     Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(camera.nx, camera.ny, shot_count.ptr());
-    Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(1, counter.ptr());
+    counter.fill(Int(0));
 
     // domain sizes
     const Vector<Float> grid_size = grid_d * grid_cells;

src_kernels_cuda_rt/raytracer_kernels_bw.cu

@@ -175,7 +175,7 @@ namespace
             Photon& photon,
             Float* __restrict__ camera_count,
             Float* __restrict__ camera_shot,
-            const int ij_cam, const int n,
+            const Int ij_cam, const int n,
             Random_number_generator<Float>& rng,
             const Vector<Float>& sun_direction,
             const Grid_knull* __restrict__ k_null_grid,
@@ -278,11 +278,11 @@ namespace
 __global__
 void ray_tracer_kernel_bw(
         const int igpt,
-        const int photons_per_pixel,
+        const Int photons_per_pixel,
         const Grid_knull* __restrict__ k_null_grid,
         Float* __restrict__ camera_count,
         Float* __restrict__ camera_shot,
-        int* __restrict__ counter,
+        Int* __restrict__ counter,
         const Float* __restrict__ k_ext, const Optics_scat* __restrict__ scat_asy,
         const Float* __restrict__ k_ext_bg, const Optics_scat* __restrict__ scat_asy_bg,
         const Float* __restrict__ z_lev_bg,
@@ -338,13 +338,13 @@ void ray_tracer_kernel_bw(
 
     const Float s_min = max(max(grid_size.z, grid_size.x), grid_size.y) * Float_epsilon;
     const Float s_min_bg = max(max(grid_size.x, grid_size.y), z_top) * Float_epsilon;
-
-    while (counter[0] < camera.nx*camera.ny*photons_per_pixel)
+    
+    while (counter[0] < camera.npix*photons_per_pixel)
     {
-        const int count = atomicAdd(&counter[0], 1);
-        const int ij_cam = count / photons_per_pixel;
+        const Int count = atomicAdd(&counter[0], 1);
+        const Int ij_cam = count / photons_per_pixel;
 
-        if (ij_cam >= camera.nx*camera.ny)
+        if (ij_cam >= camera.npix)
             return;
 
         Float weight;