updating backward raytracer

Author: Menno Veerman

include_rt/Raytracer_bw.h

@@ -20,7 +20,7 @@ class Raytracer_bw
 
         void trace_rays(
                 const Int photons_to_shoot,
-                const int n_col_x, const int n_col_y, const int n_lay,
+                const int n_col_x, const int n_col_y, const int nz, const int n_lay,
                 const Float dx_grid, const Float dy_grid, const Float dz_grid,
                 const Array_gpu<Float,1>& z_lev,
                 const Optical_props_2str_rt& optical_props,
@@ -29,7 +29,7 @@ class Raytracer_bw
                 const Float zenith_angle,
                 const Float azimuth_angle,
                 const Float tod_inc_direct,
-                const Float tod_inc_diffuse.
+                const Float tod_inc_diffuse,
                 const Float toa_factor,
                 const Float rayleigh,
                 const Array_gpu<Float,2>& col_dry,
@@ -41,6 +41,7 @@ class Raytracer_bw
                 const Int photons_to_shoot,
                 const int ncol_x, const int ncol_y, const int nz, const int nlay,
                 const Float dx_grid, const Float dy_grid, const Float dz_grid,
+                const Array_gpu<Float,1>& z_lev,
                 const Array_gpu<Float,2>& tau_gas,
                 const Array_gpu<Float,2>& ssa_gas,
                 const Array_gpu<Float,2>& asy_gas,
@@ -53,6 +54,11 @@ class Raytracer_bw
                 const Array_gpu<Float,1>& cam_data,
                 Array_gpu<Float,2>& flux_camera);
 
+        void add_camera(
+                const int cam_nx, const int cam_ny,
+                const Array_gpu<Float,2>& flux_camera,
+                Array_gpu<Float,2>& radiance);
+
         void add_xyz_camera(
                 const int cam_nx, const int cam_ny,
                 const Array_gpu<Float,1>& xyz_factor,

include_rt/gpt_combine_kernel_launcher_cuda_rt.h

@@ -31,33 +31,35 @@
 
 namespace gpt_combine_kernel_launcher_cuda_rt
 {
-    
+
     void get_from_gpoint(const int ncol, const int igpt,
                   Float* var_full, const Float* var_sub);
 
     void add_from_gpoint(const int ncol, const int nlay,
                   Float* var1_full, Float* var2_full, Float* var3_full,  Float* var4_full, Float* var5_full,
                   const Float* var1_sub, const Float* var2_sub, const Float* var3_sub, const Float* var4_sub, const Float* var5_sub);
-    
+
     void add_from_gpoint(const int ncol, const int nlay,
                   Float* var1_full, Float* var2_full, Float* var3_full,  Float* var4_full,
                   const Float* var1_sub, const Float* var2_sub, const Float* var3_sub, const Float* var4_sub);
-    
+
     void add_from_gpoint(const int ncol, const int nlay,
                   Float* var1_full, Float* var2_full,
                   const Float* var1_sub, const Float* var2_sub);
-    
-    
+
+    void add_from_gpoint(const int ncol, const int nlay,
+                  Float* var1_full, const Float* var1_sub);
+
     void add_from_gpoint(const int ncol, const int nlay,
                   Float* var1_full, Float* var2_full, Float* var3_full,
                   const Float* var1_sub, const Float* var2_sub, const Float* var3_sub);
 
-    
+
     void get_from_gpoint(const int ncol, const int nlay, const int igpt,
                   Float* var1_full, Float* var2_full, Float* var3_full,  Float* var4_full,
                   const Float* var1_sub, const Float* var2_sub, const Float* var3_sub, const Float* var4_sub);
 
-    
+
     void get_from_gpoint(const int ncol, const int nlay, const int igpt,
                   Float* var1_full, Float* var2_full, Float* var3_full,
                   const Float* var1_sub, const Float* var2_sub, const Float* var3_sub);

include_rt/raytracer_kernels_bw.h

@@ -7,17 +7,17 @@
 
 #ifdef RTE_RRTMGP_SINGLE_PRECISION
 //using Float = float;
-constexpr int block_size= 512;
-constexpr int grid_size = 1024;
+constexpr int block_size= 256;
+constexpr int grid_size = 512;
 #else
 //using Float = double;
-constexpr int block_size = 512;
+constexpr int block_size = 256;
 constexpr int grid_size = 256;
 #endif
 using Int = unsigned long long;
-constexpr int ngrid_x = 20;
-constexpr int ngrid_y = 20;
-constexpr int ngrid_z = 35;
+constexpr int ngrid_x = 48;
+constexpr int ngrid_y = 48;
+constexpr int ngrid_z = 36;
 constexpr Float k_null_gas_min = Float(1.e-3);
 
 
@@ -45,16 +45,17 @@ void ray_tracer_kernel_bw(
         const Int photons_to_shoot,
         const Grid_knull* __restrict__ k_null_grid,
         Float* __restrict__ camera_count,
+        Float* __restrict__ camera_diff,
         Float* __restrict__ camera_shot,
         int* __restrict__ counter,
         const int cam_nx, const int cam_ny, const Float* __restrict__ cam_data,
         const Optics_ext* __restrict__ k_ext, const Optics_scat* __restrict__ ssa_asy,
+        const Optics_ext* __restrict__ k_ext_bg, const Optics_scat* __restrict__ ssa_asy_bg,
+        const Float* __restrict__ z_lev_bg,
         const Float* __restrict__ surface_albedo,
         const Float mu,
-        const Float bg_trans,
-        const Float tod_frac_diffuse,
         const Float x_size, const Float y_size, const Float z_size,
         const Float dx_grid, const Float dy_grid, const Float dz_grid,
         const Float dir_x, const Float dir_y, const Float dir_z,
-        const int itot, const int jtot, const int ktot);
+        const int itot, const int jtot, const int ktot, const int nbg);
 #endif

include_test/Radiation_solver_bw.h

@@ -153,11 +153,13 @@ class Radiation_solver_shortwave
                 const Array_gpu<Float,1>& grid_dims,
                 Array_gpu<Float,2>& col_dry,
                 const Array_gpu<Float,2>& sfc_alb_dir, const Array_gpu<Float,2>& sfc_alb_dif,
-                const Array_gpu<Float,1>& tsi_scaling, const Array_gpu<Float,1>& mu0,
+                const Array_gpu<Float,1>& tsi_scaling,
+                const Array_gpu<Float,1>& mu0, const Array_gpu<Float,1>& azi,
                 const Array_gpu<Float,2>& lwp, const Array_gpu<Float,2>& iwp,
                 const Array_gpu<Float,2>& rel, const Array_gpu<Float,2>& rei,
                 const Array_gpu<Float,1>& cam_data,
                 Array_gpu<Float,3>& XYZ);
+        #endif
 
         #ifdef __CUDACC__
         void solve_gpu_bb(
@@ -171,7 +173,8 @@ class Radiation_solver_shortwave
                 const Array_gpu<Float,1>& grid_dims,
                 Array_gpu<Float,2>& col_dry,
                 const Array_gpu<Float,2>& sfc_alb_dir, const Array_gpu<Float,2>& sfc_alb_dif,
-                const Array_gpu<Float,1>& tsi_scaling, const Array_gpu<Float,1>& mu0,
+                const Array_gpu<Float,1>& tsi_scaling,
+                const Array_gpu<Float,1>& mu0, const Array_gpu<Float,1>& azi,
                 const Array_gpu<Float,2>& lwp, const Array_gpu<Float,2>& iwp,
                 const Array_gpu<Float,2>& rel, const Array_gpu<Float,2>& rei,
                 const Array_gpu<Float,1>& cam_data,

src_cuda_rt/Raytracer.cu

@@ -123,7 +123,7 @@ namespace
 
     __global__
     void count_to_flux_3d(
-            const int ncol_x, const int ncol_y, const int nlay, const Float photons_per_col, 
+            const int ncol_x, const int ncol_y, const int nlay, const Float photons_per_col,
             const Float dz_grid, const Float toa_src,
             const Float* __restrict__ count_1, const Float* __restrict__ count_2,
             Float* __restrict__ flux_1, Float* __restrict__ flux_2)
@@ -266,7 +266,7 @@ void Raytracer::trace_rays(
     // number of photons per thread, this should a power of 2 and nonzero
     Float photons_per_thread_tmp = std::max(Float(1), static_cast<Float>(photons_total) / (grid_size * block_size));
     Int photons_per_thread = pow(Float(2.), std::floor(std::log2(Float(photons_per_thread_tmp))));
-
+    printf(" raytracing \n");
     ray_tracer_kernel<<<grid, block>>>(
             photons_per_thread,
             qrng_grid_x,