Merge pull request #35 from MennoVeerman/update_to_reference_integrat…

…ed_background

better TOD-TOA handling and 1D ray tracing

config/snellius.cmake

@@ -68,8 +68,9 @@ set(LIBS ${NETCDF_LIB_C} ${HDF5_LIB} ${SZIP_LIB})
 
 if(USECUDA)
     set(CUDA_PROPAGATE_HOST_FLAGS OFF)
+    set(CMAKE_CUDA_ARCHITECTURES 80)
     set(LIBS ${LIBS} -rdynamic $ENV{EBROOTCUDA}/lib64/libcufft.so)
-    set(USER_CUDA_FLAGS "-arch=sm_80 -std=c++14 --expt-relaxed-constexpr")
+    set(USER_CUDA_FLAGS "-arch=sm_80 -std=c++14 -O3 --expt-relaxed-constexpr")
     set(USER_CUDA_FLAGS_RELEASE "-DNDEBUG")
     add_definitions(-DRTE_RRTMGP_GPU_MEMPOOL_CUDA)
 endif()

include_rt/Raytracer.h

@@ -23,6 +23,7 @@ class Raytracer
 
         void trace_rays(
                 const int qrng_gpt_offset,
+                const bool switch_independent_column,
                 const Int photons_per_pixel,
                 const Raytracer_definitions::Vector<int> grid_cells,
                 const Raytracer_definitions::Vector<Float> grid_d,

include_rt/Raytracer_bw.h

@@ -96,6 +96,18 @@ class Raytracer_bw
                 const Float total_source,
                 Array_gpu<Float,3>& XYZ);
 
+        void accumulate_clouds(
+                const Vector<Float>& grid_d,
+                const Vector<int>& grid_cells,
+                const Array_gpu<Float,2>& lwp,
+                const Array_gpu<Float,2>& iwp,
+                const Array_gpu<Float,2>& tau_cloud,
+                const Camera& camera,
+                Array_gpu<Float,2>& liwp_cam,
+                Array_gpu<Float,2>& tauc_cam,
+                Array_gpu<Float,2>& dist_cam,
+                Array_gpu<Float,2>& zen_cam);
+
     private:
 
 };

include_rt/raytracer_functions.h

@@ -20,6 +20,8 @@ namespace Raytracer_functions
     template<typename T> static inline __host__ __device__
     Vector<T> operator-(const Vector<T> v1, const Vector<T> v2) { return Vector<T>{v1.x-v2.x, v1.y-v2.y, v1.z-v2.z}; }
     template<typename T> static inline __host__ __device__
+    Vector<T> operator-(const Vector<T> v, const Float s) { return Vector<T>{v.x-s, v.y-s, v.z-s}; }
+    template<typename T> static inline __host__ __device__
     Vector<T> operator+(const Vector<T> v, const Float s) { return Vector<T>{s+v.x, s+v.y, s+v.z}; }
     template<typename T> static inline __host__ __device__
     Vector<T> operator+(const Vector<T> v1, const Vector<T> v2) { return Vector<T>{v1.x+v2.x, v1.y+v2.y, v1.z+v2.z}; }

include_rt_kernels/raytracer_kernels.h

@@ -20,6 +20,7 @@ constexpr Float k_null_gas_min = Float(1.e-3);
 
 __global__
 void ray_tracer_kernel(
+        const Bool independent_column,
         const Int photons_to_shoot,
         const Int qrng_grid_x,
         const Int qrng_grid_y,

include_rt_kernels/raytracer_kernels_bw.h

@@ -25,6 +25,11 @@ 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);
+constexpr Float sun_solid_angle = Float(6.799910294339209e-05); // 2.*M_PI*(1-cos_half_angle);
+constexpr Float sun_solid_angle_reciprocal = Float(14706.07635563193); 
+
 
 struct Grid_knull
 {
@@ -116,4 +121,19 @@ void ray_tracer_kernel_bw(
         const Float* __restrict__ mie_phase,
         const Float* __restrict__ mie_phase_ang,
         const int mie_table_size);
+
+__global__
+void accumulate_clouds_kernel(
+    const Float* __restrict__ lwp, 
+    const Float* __restrict__ iwp, 
+    const Float* __restrict__ tau_cloud,
+    const Vector<Float> grid_d,
+    const Vector<Float> grid_size,
+    const Vector<int> grid_cells,
+    Float* __restrict__ liwp_cam, 
+    Float* __restrict__ tauc_cam, 
+    Float* __restrict__ dist_cam, 
+    Float* __restrict__ zen_cam, 
+    const Camera camera);
+
 #endif

include_test/Radiation_solver_bw.h

@@ -145,8 +145,10 @@ class Radiation_solver_shortwave
                 Array<Float,3>& sw_bnd_flux_dn_dir, Array<Float,3>& sw_bnd_flux_net) const;
 
         void load_mie_tables(
-                const std::string& file_name_mie,
-                const bool switch_broadband);
+                const std::string& file_name_mie_bb,
+                const std::string& file_name_mie_im,
+                const bool switch_broadband,
+                const bool switch_image);
 
         #ifdef __CUDACC__
         void solve_gpu(
@@ -157,6 +159,8 @@ class Radiation_solver_shortwave
                 const bool switch_lu_albedo,
                 const bool switch_delta_cloud,
                 const bool switch_delta_aerosol,
+                const bool switch_cloud_cam,
+                const bool switch_raytracing,
                 const Vector<int>& grid_cells,
                 const Vector<Float>& grid_d,
                 const Vector<int>& kn_grid,
@@ -175,7 +179,11 @@ class Radiation_solver_shortwave
                 const Array_gpu<Float,2>& rh,
                 const Aerosol_concs_gpu& aerosol_concs,
                 const Camera& camera,
-                Array_gpu<Float,3>& XYZ);
+                Array_gpu<Float,3>& XYZ,
+                Array_gpu<Float,2>& liwp_cam,
+                Array_gpu<Float,2>& tauc_cam,
+                Array_gpu<Float,2>& dist_cam,
+                Array_gpu<Float,2>& zen_cam);
         #endif
 
         #ifdef __CUDACC__
@@ -186,6 +194,8 @@ class Radiation_solver_shortwave
                 const bool switch_lu_albedo,
                 const bool switch_delta_cloud,
                 const bool switch_delta_aerosol,
+                const bool switch_cloud_cam,
+                const bool switch_raytracing,
                 const Vector<int>& grid_cells,
                 const Vector<Float>& grid_d,
                 const Vector<int>& kn_grid,
@@ -204,7 +214,11 @@ class Radiation_solver_shortwave
                 const Array_gpu<Float,2>& rh,
                 const Aerosol_concs_gpu& aerosol_concs,
                 const Camera& camera,
-                Array_gpu<Float,2>& radiance);
+                Array_gpu<Float,2>& radiance,
+                Array_gpu<Float,2>& liwp_cam,
+                Array_gpu<Float,2>& tauc_cam,
+                Array_gpu<Float,2>& dist_cam,
+                Array_gpu<Float,2>& zen_cam);
 
         int get_n_gpt_gpu() const { return this->kdist_gpu->get_ngpt(); };
         int get_n_bnd_gpu() const { return this->kdist_gpu->get_nband(); };
@@ -232,11 +246,16 @@ class Radiation_solver_shortwave
         std::unique_ptr<Optical_props_2str_rt> cloud_optical_props;
         std::unique_ptr<Optical_props_2str_rt> aerosol_optical_props;
 
-        Array_gpu<Float,1> mie_phase_angs;
+        Array_gpu<Float,1> mie_phase_angs_bb;
+        Array_gpu<Float,3> mie_cdfs_bb;
+        Array_gpu<Float,4> mie_angs_bb;
+        Array_gpu<Float,4> mie_phase_bb;
 
+        Array_gpu<Float,1> mie_phase_angs;
         Array_gpu<Float,3> mie_cdfs;
         Array_gpu<Float,4> mie_angs;
         Array_gpu<Float,4> mie_phase;
+
         #endif
 };
 #endif

include_test/Radiation_solver_rt.h

@@ -100,13 +100,13 @@ class Radiation_solver_shortwave
         void solve_gpu(
                 const bool switch_fluxes,
                 const bool switch_raytracing,
+                const bool switch_independent_column,
                 const bool switch_cloud_optics,
                 const bool switch_cloud_mie,
                 const bool switch_aerosol_optics,
                 const bool switch_single_gpt,
                 const bool switch_delta_cloud,
                 const bool switch_delta_aerosol,
-                const bool switch_clear_sky_tod,
                 const int single_gpt,
                 const Int ray_count,
                 const Vector<int> grid_cells,

python/set_virtual_camera.py

@@ -61,24 +61,26 @@
     cam["roll"][:]  = 0
     cam["fisheye"][:]= 1
     cam["f_zoom"][:]= 1
+    cam["fov"][:] = 80
     cam["px"][:] = 0
     cam["py"][:] = 0
     cam["pz"][:] = 0
-    cam["nx"][:] = 1024
-    cam["ny"][:] = 1024
+    cam["nx"][:] = 256
+    cam["ny"][:] = 256
 
 # example imagery settings
 if args['image']:
     cam["yaw"][:] = 0
     cam["pitch"][:] = 0
     cam["roll"][:]  = 0
     cam["fisheye"][:]= 0
+    cam["f_zoom"][:]= 1
     cam["fov"][:] = 80
     cam["px"][:] = 0.
     cam["py"][:] = 0.
     cam["pz"][:] = 500.
-    cam["nx"][:] = 1024
-    cam["ny"][:] = 1024
+    cam["nx"][:] = 256
+    cam["ny"][:] = 256
 
 for var in camera_variables:
     if not args[var] is None:

src_cuda_rt/Gas_optics_rrtmgp_rt.cu

@@ -831,33 +831,38 @@ void Gas_optics_rrtmgp_rt::find_relevant_gases_gpt(
     const int minor_start_lower = first_last_minor_lower({1, igpt});
     const int minor_end_lower   = first_last_minor_lower({2, igpt});
 
-    for (int imnr=minor_start_lower+1; imnr<=minor_end_lower+1; ++imnr)
+    if (minor_start_lower >= 0)
     {
-        const int minor_gas = idx_minor_lower({imnr});
-        add_if_not_in_vector(gases, minor_gas);
-
-        if ( (minor_scales_with_density_lower({imnr})) && (idx_minor_scaling_lower({imnr}) > 0))
+        for (int imnr=minor_start_lower+1; imnr<=minor_end_lower+1; ++imnr)
         {
-            const int scale_gas = idx_minor_scaling_lower({imnr});
-            add_if_not_in_vector(gases, scale_gas);
+            const int minor_gas = idx_minor_lower({imnr});
+            add_if_not_in_vector(gases, minor_gas);
+
+            if ( (minor_scales_with_density_lower({imnr})) && (idx_minor_scaling_lower({imnr}) > 0))
+            {
+                const int scale_gas = idx_minor_scaling_lower({imnr});
+                add_if_not_in_vector(gases, scale_gas);
+            }
         }
     }
 
     const int minor_start_upper = first_last_minor_upper({1, igpt});
     const int minor_end_upper   = first_last_minor_upper({2, igpt});
 
-    for (int imnr=minor_start_upper+1; imnr<=minor_end_upper+1; ++imnr)
+    if (minor_start_upper >= 0)
     {
-        const int minor_gas = idx_minor_upper({imnr});
-        add_if_not_in_vector(gases, minor_gas);
-
-        if ( (minor_scales_with_density_upper({imnr})) && (idx_minor_scaling_upper({imnr}) > 0))
+        for (int imnr=minor_start_upper+1; imnr<=minor_end_upper+1; ++imnr)
         {
-            const int scale_gas = idx_minor_scaling_upper({imnr});
-            add_if_not_in_vector(gases, scale_gas);
+            const int minor_gas = idx_minor_upper({imnr});
+            add_if_not_in_vector(gases, minor_gas);
+
+            if ( (minor_scales_with_density_upper({imnr})) && (idx_minor_scaling_upper({imnr}) > 0))
+            {
+                const int scale_gas = idx_minor_scaling_upper({imnr});
+                add_if_not_in_vector(gases, scale_gas);
+            }
         }
     }
-
 }
 
 __global__

src_cuda_rt/Raytracer.cu

@@ -86,8 +86,8 @@ namespace
         const int icol_x = blockIdx.x*blockDim.x + threadIdx.x;
         const int icol_y = blockIdx.y*blockDim.y + threadIdx.y;
         const int iz = blockIdx.z*blockDim.z + threadIdx.z;
-
-        if ( (icol_x < grid_cells.x) && (icol_y < grid_cells.y) && (iz < grid_cells.z) )
+        
+        if ( (icol_x < grid_cells.x) && (icol_y < grid_cells.y) && (iz < (grid_cells.z - 1)) )
         {
             const int idx = icol_x + icol_y*grid_cells.x + iz*grid_cells.y*grid_cells.x;
             const Float kext_tot = tau_tot[idx] / grid_d.z;
@@ -106,6 +106,61 @@ namespace
         }
     }
 
+    __global__
+    void bundles_optical_props_tod(
+            const Vector<int> grid_cells, const Vector<Float> grid_d, const int n_lev,
+            const Float* __restrict__ tau_tot, const Float* __restrict__ ssa_tot,
+            const Float* __restrict__ tau_cld, const Float* __restrict__ ssa_cld, const Float* __restrict__ asy_cld,
+            const Float* __restrict__ tau_aer, const Float* __restrict__ ssa_aer, const Float* __restrict__ asy_aer,
+            Float* __restrict__ k_ext, Optics_scat* __restrict__ scat_asy)
+    {
+        const int icol_x = blockIdx.x*blockDim.x + threadIdx.x;
+        const int icol_y = blockIdx.y*blockDim.y + threadIdx.y;
+
+        const int z_tod = grid_cells.z - 1;
+
+        if ( (icol_x < grid_cells.x) && (icol_y < grid_cells.y))
+        {
+            Float tau_tot_sum = Float(0.);
+            Float tausca_tot_sum = Float(0.);
+
+            Float tausca_cld_sum = Float(0.);
+            Float tauscag_cld_sum = Float(0.);
+
+            Float tausca_aer_sum = Float(0.);
+            Float tauscag_aer_sum = Float(0.);
+
+            for (int iz=z_tod; iz<n_lev; ++iz)
+            {
+                const int idx = icol_x + icol_y*grid_cells.x + iz*grid_cells.y*grid_cells.x;
+                tau_tot_sum += tau_tot[idx];
+                tausca_tot_sum += tau_tot[idx] * ssa_tot[idx];
+
+                tausca_cld_sum += tau_cld[idx] * ssa_cld[idx];
+                tauscag_cld_sum += tau_cld[idx] * ssa_cld[idx] * asy_cld[idx];
+
+                tausca_aer_sum += tau_aer[idx] * ssa_aer[idx];
+                tauscag_aer_sum += tau_aer[idx] * ssa_aer[idx] * asy_aer[idx];
+            }
+
+            const int idx = icol_x + icol_y*grid_cells.x + z_tod*grid_cells.y*grid_cells.x;
+
+            const Float kext_tot = tau_tot_sum / grid_d.z;
+
+            const Float ksca_cld = tausca_cld_sum / grid_d.z;
+            const Float ksca_aer = tausca_aer_sum / grid_d.z;
+            const Float ksca_gas = tausca_tot_sum / grid_d.z - ksca_cld - ksca_aer;
+            k_ext[idx] = kext_tot;
+
+            scat_asy[idx].k_sca_gas = ksca_gas;
+            scat_asy[idx].k_sca_cld = ksca_cld;
+            scat_asy[idx].k_sca_aer = ksca_aer;
+
+            scat_asy[idx].asy_cld = tauscag_cld_sum / tausca_cld_sum;
+            scat_asy[idx].asy_aer = tauscag_aer_sum / tausca_aer_sum;
+        }
+    }
+
 
     __global__
     void count_to_flux_2d(
@@ -169,6 +224,7 @@ Raytracer::Raytracer()
 
 void Raytracer::trace_rays(
         const int igpt,
+        const bool switch_independent_column,
         const Int photons_per_pixel,
         const Vector<int> grid_cells,
         const Vector<Float> grid_d,
@@ -197,6 +253,8 @@ void Raytracer::trace_rays(
         Array_gpu<Float,3>& flux_abs_dir,
         Array_gpu<Float,3>& flux_abs_dif)
 {
+    const int n_lay = tau_total.dim(2);
+
     // set of block and grid dimensions used in data processing kernels - requires some proper tuning later
     const int block_col_x = 8;
     const int block_col_y = 8;
@@ -215,12 +273,21 @@ void Raytracer::trace_rays(
     Array_gpu<Float,3> k_ext({grid_cells.x, grid_cells.y, grid_cells.z});
     Array_gpu<Optics_scat,3> scat_asy({grid_cells.x, grid_cells.y, grid_cells.z});
 
+    // first on the whole grid expect the extra layer
     bundles_optical_props<<<grid_3d, block_3d>>>(
             grid_cells, grid_d,
             tau_total.ptr(), ssa_total.ptr(),
             tau_cloud.ptr(), ssa_cloud.ptr(), asy_cloud.ptr(),
             tau_aeros.ptr(), ssa_aeros.ptr(), asy_aeros.ptr(),
             k_ext.ptr(), scat_asy.ptr());
+
+    // second, integrate from TOD to TOA 
+    bundles_optical_props_tod<<<grid_2d, block_2d>>>(
+            grid_cells, grid_d, n_lay,
+            tau_total.ptr(), ssa_total.ptr(),
+            tau_cloud.ptr(), ssa_cloud.ptr(), asy_cloud.ptr(),
+            tau_aeros.ptr(), ssa_aeros.ptr(), asy_aeros.ptr(),
+            k_ext.ptr(), scat_asy.ptr());
 
     // create k_null_grid
     const int block_kn_x = 8;
@@ -303,6 +370,7 @@ void Raytracer::trace_rays(
 
     const int qrng_gpt_offset = (igpt-1) * rt_kernel_grid_size * rt_kernel_block_size * photons_per_thread;
     ray_tracer_kernel<<<grid, block,sizeof(Float)*mie_table_size>>>(
+            switch_independent_column,
             photons_per_thread,
             qrng_grid_x,
             qrng_grid_y,