Merge pull request #43 from MennoVeerman/main

fix bug in null collision grid creation and some addition to rt_lite

Author: Chiil

include_rt_kernels/rte_solver_kernels_cuda_rt.h

@@ -38,6 +38,8 @@ namespace Rte_solver_kernels_cuda_rt
     void apply_BC(const int ncol, const int nlay, const int ngpt, const Bool top_at_1, Float* gpt_flux_dn);
 
     void apply_BC(const int ncol, const int nlay, const int ngpt, const Bool top_at_1, const Float* inc_flux_dif, Float* gpt_flux_dn);
+    
+    void apply_BC(const int ncol, const int nlay, const int ngpt, const Bool top_at_1, const Float inc_flux, Float* gpt_flux_dn);
 
     void sw_solver_2stream(
             const int ncol, const int nlay, const int ngpt, const Bool top_at_1,

python/set_virtual_camera.py

@@ -91,7 +91,13 @@
         cam[var][:] = args[var]
 
 if not args['sza'] is None:
+    try:
+        ncf.createVariable('sza','f4',ncf['mu0'].dimensions) 
+    except:
+        pass
+    ncf['sza'][:] = np.deg2rad(args['sza'])
     ncf['mu0'][:] = np.cos(np.deg2rad(args['sza']))
+
 if not args['azi'] is None:
     ncf['azi'][:] = np.deg2rad(args['azi'])
 

src_cuda_rt/Raytracer.cu

@@ -53,18 +53,18 @@ namespace
             const Float fz = Float(grid_cells.z) / Float(kn_grid.z);
 
             const int x0 = grid_x*fx;
-            const int x1 = floor((grid_x+1)*fx);
+            const int x1 = min(grid_cells.x-1, int(floor((grid_x+1)*fx)));
             const int y0 = grid_y*fy;
-            const int y1 = floor((grid_y+1)*fy);
+            const int y1 = min(grid_cells.y-1, int(floor((grid_y+1)*fy)));
             const int z0 = grid_z*fz;
-            const int z1 = floor((grid_z+1)*fz);
+            const int z1 = min(grid_cells.z-1, int(floor((grid_z+1)*fz)));
 
             const int ijk_grid = grid_x + grid_y*kn_grid.x + grid_z*kn_grid.y*kn_grid.x;
             Float k_null = k_ext_null_min;
 
-            for (int k=z0; k<z1; ++k)
-                for (int j=y0; j<y1; ++j)
-                    for (int i=x0; i<x1; ++i)
+            for (int k=z0; k<=z1; ++k)
+                for (int j=y0; j<=y1; ++j)
+                    for (int i=x0; i<=x1; ++i)
                     {
                         const int ijk_in = i + j*grid_cells.x + k*grid_cells.x*grid_cells.y;
                         k_null = max(k_null, k_ext[ijk_in]);

src_kernels_cuda_rt/rte_solver_kernels_launchers_rt.cu

@@ -41,6 +41,16 @@ namespace Rte_solver_kernels_cuda_rt
     }
 
 
+    void apply_BC(const int ncol, const int nlay, const int ngpt, const Bool top_at_1, const Float inc_flux, Float* gpt_flux_dn)
+    {
+        const int block_col = 32;
+        const int grid_col = ncol/block_col + (ncol%block_col > 0);
+
+        dim3 grid_gpu(grid_col);
+        dim3 block_gpu(block_col);
+        apply_BC_kernel<<<grid_gpu, block_gpu>>>(ncol, nlay, ngpt, top_at_1, inc_flux, gpt_flux_dn);
+    }
+    
     void apply_BC(const int ncol, const int nlay, const int ngpt, const Bool top_at_1, const Float* inc_flux_dif, Float* gpt_flux_dn)
     {
         const int block_col = 32;

src_kernels_cuda_rt/rte_solver_kernels_rt.cu

@@ -305,6 +305,18 @@ void apply_BC_kernel_lw(const int isfc, int ncol, const int nlay, const int ngpt
     }
 }
 
+__global__
+void apply_BC_kernel(const int ncol, const int nlay, const int ngpt, const Bool top_at_1, const Float inc_flux, Float* __restrict__ flux_dn)
+{
+    const int icol = blockIdx.x*blockDim.x + threadIdx.x;
+    if ( (icol < ncol) )
+    {
+        const int idx_out = icol + ((top_at_1 ? 0 : (nlay * ncol))); 
+        const int idx_in = icol;
+        flux_dn[idx_out] = inc_flux;
+    }
+}
+
 __global__
 void apply_BC_kernel(const int ncol, const int nlay, const int ngpt, const Bool top_at_1, const Float* __restrict__ inc_flux, Float* __restrict__ flux_dn)
 {

src_test/test_rt_lite.cu

@@ -30,6 +30,8 @@
 #include "Raytracer_bw.h"
 #include "raytracer_kernels_bw.h"
 #include "types.h"
+#include "rte_solver_kernels_cuda_rt.h"
+#include "Rte_sw_rt.h"
 #include "tools_gpu.h"
 
 
@@ -129,20 +131,22 @@ void solve_radiation(int argc, char** argv)
     // Parse the command line options.
     std::map<std::string, std::pair<bool, std::string>> command_line_switches {
         {"raytracing"        , { true,  "Use forward raytracer for irradiances. '--raytracing 256': use 256 rays per pixel" }},
-        {"bw_raytracing"     , { true,  "Use backward raytracer radiances. '--raytracing 256': use 256 rays per pixel" }},
+        {"bw-raytracing"     , { true,  "Use backward raytracer radiances. '--raytracing 256': use 256 rays per pixel" }},
+        {"two-stream"        , { true,  "Perform two-stream computations"}},
         {"cloud-mie"         , { false, "Use Mie tables for cloud scattering in ray tracer"  }},
         {"independent-column", { false, "run raytracer in independent column mode"}},
         {"profiling"         , { false, "Perform additional profiling run."         }} };
 
     std::map<std::string, std::pair<int, std::string>> command_line_ints {
         {"raytracing", {32, "Number of rays initialised at TOD per pixel."}}, 
-        {"bw_raytracing", {32, "Number of rays initialised at per camera pixel."}}} ;
+        {"bw-raytracing", {32, "Number of rays initialised at per camera pixel."}}} ;
 
     if (parse_command_line_options(command_line_switches, command_line_ints, argc, argv))
         return;
 
     const bool switch_raytracing         = command_line_switches.at("raytracing"        ).first;
-    const bool switch_bw_raytracing      = command_line_switches.at("bw_raytracing"     ).first;
+    const bool switch_bw_raytracing      = command_line_switches.at("bw-raytracing"     ).first;
+    const bool switch_two_stream         = command_line_switches.at("two-stream"        ).first;
     const bool switch_cloud_mie          = command_line_switches.at("cloud-mie"         ).first;
     const bool switch_independent_column = command_line_switches.at("independent-column").first;
     const bool switch_profiling          = command_line_switches.at("profiling"         ).first;
@@ -207,6 +211,7 @@ void solve_radiation(int argc, char** argv)
     // Read the atmospheric fields.
     const Array<Float,2> tot_tau(input_nc.get_variable<Float>("tot_tau", {n_lay, ny, nx}), {ncol, n_lay});
     const Array<Float,2> tot_ssa(input_nc.get_variable<Float>("tot_ssa", {n_lay, ny, nx}), {ncol, n_lay});
+    const Array<Float,2> tot_asy(input_nc.get_variable<Float>("tot_asy", {n_lay, ny, nx}), {ncol, n_lay});
 
     Array<Float,2> cld_tau({ncol, n_lay});	
     Array<Float,2> cld_ssa({ncol, n_lay});	
@@ -225,13 +230,16 @@ void solve_radiation(int argc, char** argv)
     aer_asy = std::move(input_nc.get_variable<Float>("aer_asy", {n_lay, ny, nx}));
 
     // read albedo, solar angles, and top-of-domain fluxes
-    Array<Float,2> sfc_albedo({1,ncol});
-    sfc_albedo.fill(input_nc.get_variable<Float>("albedo"));
+    Array<Float,2> sfc_albedo({ncol, 1});
+    sfc_albedo = std::move(input_nc.get_variable<Float>("albedo", {ny, nx}));
+    
     const Float zenith_angle = input_nc.get_variable<Float>("sza");
     const Float azimuth_angle = input_nc.get_variable<Float>("azi");
     const Float tod_dir = input_nc.get_variable<Float>("tod_direct");
 
-        
+    Array<Float,1> mu0({ncol});	
+    mu0.fill(cos(zenith_angle));
+
     Camera camera;
     if (switch_bw_raytracing)
     {
@@ -267,6 +275,10 @@ void solve_radiation(int argc, char** argv)
     Array_gpu<Float,3> flux_abs_dif({nx, ny, nz});
     Array_gpu<Float,2> radiance({camera.nx, camera.ny});
 
+    Array_gpu<Float,2> flux_dn_2stream;
+    Array_gpu<Float,2> flux_up_2stream;
+    Array_gpu<Float,2> flux_dn_dir_2stream;
+    
     // empty arrays (mie scattering not (yet) supported in lite version)
     Array<Float,2> mie_cdfs_c;
     Array<Float,3> mie_angs_c;
@@ -317,34 +329,69 @@ void solve_radiation(int argc, char** argv)
     lum_c.fill(Float(0.));
     Array_gpu<Float,1> land_use_map(lum_c);
 
-
+    //// GPU arrays
+    Array_gpu<Float,2> tot_tau_g(tot_tau);
+    Array_gpu<Float,2> tot_ssa_g(tot_ssa);
+    Array_gpu<Float,2> tot_asy_g(tot_asy);
+    Array_gpu<Float,2> cld_tau_g(cld_tau);
+    Array_gpu<Float,2> cld_ssa_g(cld_ssa);
+    Array_gpu<Float,2> cld_asy_g(cld_asy);
+    Array_gpu<Float,2> aer_tau_g(aer_tau);
+    Array_gpu<Float,2> aer_ssa_g(aer_ssa);
+    Array_gpu<Float,2> aer_asy_g(aer_asy);
+    Array_gpu<Float,2> sfc_albedo_g(sfc_albedo);
+    Array_gpu<Float,1> mu0_g(mu0);
+   
     ////// CREATE THE OUTPUT FILE //////
     // Create the general dimensions and arrays.
     Status::print_message("Preparing NetCDF output file.");
 
     Netcdf_file output_nc("rt_lite_output.nc", Netcdf_mode::Create);
-    if (switch_raytracing)
+    if (switch_raytracing || switch_two_stream)
     {
         output_nc.add_dimension("x", nx);
         output_nc.add_dimension("y", ny);
-        output_nc.add_dimension("z", n_z_in);
+        output_nc.add_dimension("z", nz);
+        output_nc.add_dimension("lev", n_lay+1);
     }
     if (switch_bw_raytracing)
     {
         output_nc.add_dimension("nx", camera.nx);
         output_nc.add_dimension("ny", camera.ny);
     }
 
-    //// GPU arrays
-    Array_gpu<Float,2> tot_tau_g(tot_tau);
-    Array_gpu<Float,2> tot_ssa_g(tot_ssa);
-    Array_gpu<Float,2> cld_tau_g(cld_tau);
-    Array_gpu<Float,2> cld_ssa_g(cld_ssa);
-    Array_gpu<Float,2> cld_asy_g(cld_asy);
-    Array_gpu<Float,2> aer_tau_g(aer_tau);
-    Array_gpu<Float,2> aer_ssa_g(aer_ssa);
-    Array_gpu<Float,2> aer_asy_g(aer_asy);
-    Array_gpu<Float,2> sfc_albedo_g(sfc_albedo);
+
+    if (switch_two_stream)
+    {
+        flux_up_2stream.set_dims({ncol, n_lay+1});
+        flux_dn_2stream.set_dims({ncol, n_lay+1});
+        flux_dn_dir_2stream.set_dims({ncol, n_lay+1});
+        
+        Rte_sw_rt rte_sw;
+        Rte_solver_kernels_cuda_rt::apply_BC(ncol, n_lay, 1, 0, tod_dir * cos(zenith_angle), flux_dn_dir_2stream.ptr());
+        Rte_solver_kernels_cuda_rt::apply_BC(ncol, n_lay, 1, 0, flux_dn_2stream.ptr());
+
+        Rte_solver_kernels_cuda_rt::sw_solver_2stream(
+            ncol, n_lay, 1, 0,
+            tot_tau_g.ptr(),
+            tot_ssa_g.ptr(),
+            tot_asy_g.ptr(),
+            mu0_g.ptr(),
+            sfc_albedo_g.ptr(), sfc_albedo_g.ptr(),
+            flux_up_2stream.ptr(), flux_dn_2stream.ptr(), flux_dn_dir_2stream.ptr());
+
+        Array<Float,2> flux_up_2stream_c(flux_up_2stream);
+        Array<Float,2> flux_dn_2stream_c(flux_dn_2stream);
+        Array<Float,2> flux_dn_dir_2stream_c(flux_dn_dir_2stream);
+
+        auto nc_up_2stream = output_nc.add_variable<Float>("flux_up_2stream" , {"lev", "y", "x"});
+        auto nc_dn_2stream = output_nc.add_variable<Float>("flux_dn_2stream" , {"lev", "y", "x"});
+        auto nc_dn_dir_2stream = output_nc.add_variable<Float>("flux_dn_dir_2stream" , {"lev", "y", "x"});
+
+        nc_up_2stream.insert(flux_up_2stream_c  .v(), {0, 0, 0});
+        nc_dn_2stream.insert(flux_dn_2stream_c  .v(), {0, 0, 0});
+        nc_dn_dir_2stream.insert(flux_dn_dir_2stream_c  .v(), {0, 0, 0});
+    }
 
     if (switch_raytracing)
     {