From e15607bbc201d1b35e16672f49ef53458da3797f Mon Sep 17 00:00:00 2001 From: Menno Date: Sat, 2 Nov 2024 18:40:18 +0100 Subject: [PATCH 1/5] ensure rt_ite version also works with a background atmosphere --- src_test/test_rt_lite.cu | 45 +++++++++++++++++++++------------------- 1 file changed, 24 insertions(+), 21 deletions(-) diff --git a/src_test/test_rt_lite.cu b/src_test/test_rt_lite.cu index 66691a4..1f14741 100644 --- a/src_test/test_rt_lite.cu +++ b/src_test/test_rt_lite.cu @@ -137,7 +137,7 @@ void solve_radiation(int argc, char** argv) return; const bool switch_raytracing = command_line_switches.at("raytracing" ).first; - const bool switch_independent_column = command_line_switches.at("independent-column").first; + const bool switch_independent_column = command_line_switches.at("independent-column").first; const bool switch_profiling = command_line_switches.at("profiling" ).first; // Print the options to the screen. @@ -159,14 +159,18 @@ void solve_radiation(int argc, char** argv) Netcdf_file input_nc("rt_lite_input.nc", Netcdf_mode::Read); const int nx = input_nc.get_dimension_size("x"); const int ny = input_nc.get_dimension_size("y"); - const int nz = input_nc.get_dimension_size("z"); + const int n_z_in = input_nc.get_dimension_size("z"); + const int n_lay = input_nc.get_dimension_size("lay"); + + const int nz = (n_z_in < n_lay) ? n_z_in+1 : n_z_in; + const int ncol = nx*ny; const Vector grid_cells = {nx, ny, nz}; // Read the x,y,z dimensions if raytracing is enabled const Array grid_x(input_nc.get_variable("x", {nx}), {nx}); const Array grid_y(input_nc.get_variable("y", {ny}), {ny}); - const Array grid_z(input_nc.get_variable("z", {nz}), {nz}); + const Array grid_z(input_nc.get_variable("z", {n_z_in}), {n_z_in}); const Float dx = grid_x({2}) - grid_x({1}); const Float dy = grid_y({2}) - grid_y({1}); @@ -179,24 +183,24 @@ void solve_radiation(int argc, char** argv) const Vector kn_grid = {ngrid_x, ngrid_y, ngrid_z}; // Read the atmospheric fields. - const Array tot_tau(input_nc.get_variable("tot_tau", {nz, ny, nx}), {ncol, nz}); - const Array tot_ssa(input_nc.get_variable("tot_ssa", {nz, ny, nx}), {ncol, nz}); + const Array tot_tau(input_nc.get_variable("tot_tau", {n_lay, ny, nx}), {ncol, n_lay}); + const Array tot_ssa(input_nc.get_variable("tot_ssa", {n_lay, ny, nx}), {ncol, n_lay}); - Array cld_tau({ncol, nz}); - Array cld_ssa({ncol, nz}); - Array cld_asy({ncol, nz}); + Array cld_tau({ncol, n_lay}); + Array cld_ssa({ncol, n_lay}); + Array cld_asy({ncol, n_lay}); - cld_tau = std::move(input_nc.get_variable("cld_tau", {nz, ny, nx})); - cld_ssa = std::move(input_nc.get_variable("cld_ssa", {nz, ny, nx})); - cld_asy = std::move(input_nc.get_variable("cld_asy", {nz, ny, nx})); + cld_tau = std::move(input_nc.get_variable("cld_tau", {n_lay, ny, nx})); + cld_ssa = std::move(input_nc.get_variable("cld_ssa", {n_lay, ny, nx})); + cld_asy = std::move(input_nc.get_variable("cld_asy", {n_lay, ny, nx})); - Array aer_tau({ncol, nz}); - Array aer_ssa({ncol, nz}); - Array aer_asy({ncol, nz}); + Array aer_tau({ncol, n_lay}); + Array aer_ssa({ncol, n_lay}); + Array aer_asy({ncol, n_lay}); - aer_tau = std::move(input_nc.get_variable("aer_tau", {nz, ny, nx})); - aer_ssa = std::move(input_nc.get_variable("aer_ssa", {nz, ny, nx})); - aer_asy = std::move(input_nc.get_variable("aer_asy", {nz, ny, nx})); + aer_tau = std::move(input_nc.get_variable("aer_tau", {n_lay, ny, nx})); + aer_ssa = std::move(input_nc.get_variable("aer_ssa", {n_lay, ny, nx})); + aer_asy = std::move(input_nc.get_variable("aer_asy", {n_lay, ny, nx})); // read albedo, solar angles, and top-of-domain fluxes Array sfc_albedo({1,ncol}); @@ -204,7 +208,6 @@ void solve_radiation(int argc, char** argv) const Float zenith_angle = input_nc.get_variable("sza"); const Float azimuth_angle = input_nc.get_variable("azi"); const Float tod_dir = input_nc.get_variable("tod_direct"); - const Float tod_dif = input_nc.get_variable("tod_diffuse"); // output arrays Array_gpu flux_tod_dn({nx, ny}); @@ -228,7 +231,7 @@ void solve_radiation(int argc, char** argv) Netcdf_file output_nc("rt_lite_output.nc", Netcdf_mode::Create); output_nc.add_dimension("x", nx); output_nc.add_dimension("y", ny); - output_nc.add_dimension("z", nz); + output_nc.add_dimension("z", n_z_in); //// GPU arrays Array_gpu tot_tau_g(tot_tau); @@ -278,8 +281,8 @@ void solve_radiation(int argc, char** argv) sfc_albedo_g, zenith_angle, azimuth_angle, - tod_dir, - tod_dif, + tod_dir * std::cos(zenith_angle), + Float(0.), flux_tod_dn, flux_tod_up, flux_sfc_dir, From 4695823bd15968c3a64b35b56f2984335e5e0630 Mon Sep 17 00:00:00 2001 From: Menno Date: Tue, 5 Nov 2024 09:26:40 +0100 Subject: [PATCH 2/5] bw cam pixel counter Int instead of int, also precompute cam nx*ny --- include_rt_kernels/raytracer_kernels_bw.h | 5 +++-- src_cuda_rt/Raytracer_bw.cu | 9 +++++---- src_kernels_cuda_rt/raytracer_kernels_bw.cu | 16 ++++++++-------- src_test/test_rte_rrtmgp_bw.cu | 1 + 4 files changed, 17 insertions(+), 14 deletions(-) diff --git a/include_rt_kernels/raytracer_kernels_bw.h b/include_rt_kernels/raytracer_kernels_bw.h index d69471c..133d361 100644 --- a/include_rt_kernels/raytracer_kernels_bw.h +++ b/include_rt_kernels/raytracer_kernels_bw.h @@ -91,17 +91,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, diff --git a/src_cuda_rt/Raytracer_bw.cu b/src_cuda_rt/Raytracer_bw.cu index 6b8813f..85eff1b 100644 --- a/src_cuda_rt/Raytracer_bw.cu +++ b/src_cuda_rt/Raytracer_bw.cu @@ -456,11 +456,12 @@ void Raytracer_bw::trace_rays( Array_gpu camera_count({camera.nx, camera.ny}); Array_gpu shot_count({camera.nx, camera.ny}); - Array_gpu counter({1}); + Array_gpu 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 grid_size = grid_d * grid_cells; @@ -604,11 +605,11 @@ void Raytracer_bw::trace_rays_bb( Array_gpu camera_count({camera.nx, camera.ny}); Array_gpu shot_count({camera.nx, camera.ny}); - Array_gpu counter({1}); + Array_gpu 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 grid_size = grid_d * grid_cells; diff --git a/src_kernels_cuda_rt/raytracer_kernels_bw.cu b/src_kernels_cuda_rt/raytracer_kernels_bw.cu index abc2747..609b607 100644 --- a/src_kernels_cuda_rt/raytracer_kernels_bw.cu +++ b/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& rng, const Vector& 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; diff --git a/src_test/test_rte_rrtmgp_bw.cu b/src_test/test_rte_rrtmgp_bw.cu index ef57947..852ce34 100644 --- a/src_test/test_rte_rrtmgp_bw.cu +++ b/src_test/test_rte_rrtmgp_bw.cu @@ -321,6 +321,7 @@ void solve_radiation(int argc, char** argv) camera.nx = int(cam_in.get_variable("nx")); camera.ny = int(cam_in.get_variable("ny")); + camera.npix = Int(camera.nx * camera.ny); camera.setup_rotation_matrix(cam_in.get_variable("yaw"), cam_in.get_variable("pitch"), From a787779b5b403d93cb7d66a42db91a9d41f4d99d Mon Sep 17 00:00:00 2001 From: Menno Date: Tue, 5 Nov 2024 09:27:37 +0100 Subject: [PATCH 3/5] add bw raytracer to rt_lite --- include_rt_kernels/raytracer_kernels_bw.h | 1 - src_test/test_rt_lite.cu | 406 ++++++++++++++++------ 2 files changed, 296 insertions(+), 111 deletions(-) diff --git a/include_rt_kernels/raytracer_kernels_bw.h b/include_rt_kernels/raytracer_kernels_bw.h index 133d361..3663237 100644 --- a/include_rt_kernels/raytracer_kernels_bw.h +++ b/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); diff --git a/src_test/test_rt_lite.cu b/src_test/test_rt_lite.cu index 1f14741..aa49a86 100644 --- a/src_test/test_rt_lite.cu +++ b/src_test/test_rt_lite.cu @@ -27,6 +27,8 @@ #include "Array.h" #include "Raytracer.h" #include "raytracer_kernels.h" +#include "Raytracer_bw.h" +#include "raytracer_kernels_bw.h" #include "types.h" #include "tools_gpu.h" @@ -126,31 +128,51 @@ void solve_radiation(int argc, char** argv) ////// FLOW CONTROL SWITCHES ////// // Parse the command line options. std::map> command_line_switches { - {"raytracing" , { true, "Use raytracing for flux computation. '--raytracing 256': use 256 rays per pixel" }}, + {"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" }}, + {"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> command_line_ints { - {"raytracing", {32, "Number of rays initialised at TOD per pixel per quadraute."}}} ; + {"raytracing", {32, "Number of rays initialised at TOD per 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_raytracing = command_line_switches.at("raytracing" ).first; + const bool switch_bw_raytracing = command_line_switches.at("bw_raytracing" ).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; + const bool switch_profiling = command_line_switches.at("profiling" ).first; // Print the options to the screen. print_command_line_options(command_line_switches, command_line_ints); - Int photons_per_pixel = Int(command_line_ints.at("raytracing").first); - if (Float(int(std::log2(Float(photons_per_pixel)))) != std::log2(Float(photons_per_pixel))) + Int photons_per_pixel; + Int photons_per_pixel_bw; + if (switch_raytracing) { - std::string error = "number of photons per pixel should be a power of 2 "; - throw std::runtime_error(error); + photons_per_pixel = Int(command_line_ints.at("raytracing").first); + if (Float(int(std::log2(Float(photons_per_pixel)))) != std::log2(Float(photons_per_pixel))) + { + std::string error = "number of photons per pixel should be a power of 2 "; + throw std::runtime_error(error); + } + Status::print_message("Using "+ std::to_string(photons_per_pixel) + " rays per pixel"); } - Status::print_message("Using "+ std::to_string(photons_per_pixel) + " rays per pixel"); + if (switch_bw_raytracing) + { + photons_per_pixel_bw = Int(command_line_ints.at("bw_raytracing").first); + if (Float(int(std::log2(Float(photons_per_pixel_bw)))) != std::log2(Float(photons_per_pixel_bw))) + { + std::string error = "number of bw photons per pixel should be a power of 2 "; + throw std::runtime_error(error); + } + Status::print_message("Using "+ std::to_string(photons_per_pixel_bw) + " bw rays per pixel"); + } ////// READ THE ATMOSPHERIC DATA ////// @@ -165,13 +187,13 @@ void solve_radiation(int argc, char** argv) const int nz = (n_z_in < n_lay) ? n_z_in+1 : n_z_in; const int ncol = nx*ny; - const Vector grid_cells = {nx, ny, nz}; // Read the x,y,z dimensions if raytracing is enabled const Array grid_x(input_nc.get_variable("x", {nx}), {nx}); const Array grid_y(input_nc.get_variable("y", {ny}), {ny}); const Array grid_z(input_nc.get_variable("z", {n_z_in}), {n_z_in}); - + Array z_lev; + const Float dx = grid_x({2}) - grid_x({1}); const Float dy = grid_y({2}) - grid_y({1}); const Float dz = grid_z({2}) - grid_z({1}); @@ -209,7 +231,33 @@ void solve_radiation(int argc, char** argv) const Float azimuth_angle = input_nc.get_variable("azi"); const Float tod_dir = input_nc.get_variable("tod_direct"); - // output arrays + + Camera camera; + if (switch_bw_raytracing) + { + Netcdf_group cam_in = input_nc.get_group("camera-settings"); + camera.f_zoom = cam_in.get_variable("f_zoom"); + camera.fov = cam_in.get_variable("fov"); + camera.fisheye= int(cam_in.get_variable("fisheye")); + camera.position = {cam_in.get_variable("px"), + cam_in.get_variable("py"), + cam_in.get_variable("pz")}; + + camera.nx = int(cam_in.get_variable("nx")); + camera.ny = int(cam_in.get_variable("ny")); + camera.npix = Int(camera.nx * camera.ny); + + camera.setup_rotation_matrix(cam_in.get_variable("yaw"), + cam_in.get_variable("pitch"), + cam_in.get_variable("roll")); + camera.setup_normal_camera(camera); + + z_lev.set_dims({n_lay+1}); + z_lev = std::move(input_nc.get_variable("z_lev", {n_lay+1})); + + } + + // output arrays (setting all dimensions even if we only run FW or BW is note very memory efficiency, so deserves conditional dimension assigment at some point Array_gpu flux_tod_dn({nx, ny}); Array_gpu flux_tod_up({nx, ny}); Array_gpu flux_sfc_dir({nx, ny}); @@ -217,22 +265,76 @@ void solve_radiation(int argc, char** argv) Array_gpu flux_sfc_up({nx, ny}); Array_gpu flux_abs_dir({nx, ny, nz}); Array_gpu flux_abs_dif({nx, ny, nz}); - + Array_gpu radiance({camera.nx, camera.ny}); + + // empty arrays (mie scattering not (yet) supported in lite version) + Array mie_cdfs_c; + Array mie_angs_c; + Array mie_cdfs_bw_c; + Array mie_angs_bw_c; + Array mie_phase_bw_c; + Array mie_phase_angs_bw_c; + Array rel_c({ncol, n_lay}); + + if (switch_cloud_mie) + { + // lwp.set_dims({n_col, n_lay}); + // lwp = std::move(input_nc.get_variable("lwp", {n_lay, n_col_y, n_col_x})); + + const int n_re = input_nc.get_dimension_size("r_eff"); + const int n_mie = input_nc.get_dimension_size("n_ang"); + + mie_cdfs_c.set_dims({n_mie, 1}); + mie_angs_c.set_dims({n_mie, n_re, 1}); + mie_cdfs_bw_c.set_dims({n_mie, 1, 1}); + mie_angs_bw_c.set_dims({n_mie, n_re, 1, 1}); + mie_phase_bw_c.set_dims({n_mie, n_re, 1, 1}); + mie_phase_angs_bw_c.set_dims({n_mie}); + + mie_cdfs_c = std::move(input_nc.get_variable("phase_cdf", {1, 1, n_mie})); + mie_angs_c = std::move(input_nc.get_variable("phase_cdf_angle", {1, 1, n_re, n_mie})); + mie_cdfs_bw_c = std::move(input_nc.get_variable("phase_cdf", {1, 1, n_mie})); + mie_angs_bw_c = std::move(input_nc.get_variable("phase_cdf_angle", {1, 1, n_re, n_mie})); + mie_phase_bw_c = std::move(input_nc.get_variable("phase", {1, 1, n_re, n_mie})); + mie_phase_angs_bw_c = std::move(input_nc.get_variable("phase_angle", {n_mie})); + + rel_c = std::move(input_nc.get_variable("rel", {n_lay, ny, nx})); + } + else + { + rel_c.fill(Float(0.)); + } + + Array_gpu rel(rel_c); + Array_gpu mie_cdfs(mie_cdfs_c); + Array_gpu mie_angs(mie_angs_c); + Array_gpu mie_cdfs_bw(mie_cdfs_bw_c); + Array_gpu mie_angs_bw(mie_angs_bw_c); + Array_gpu mie_phase_bw(mie_phase_bw_c); + Array_gpu mie_phase_angs_bw(mie_phase_angs_bw_c); + + Array lum_c({ncol}); + lum_c.fill(Float(0.)); + Array_gpu land_use_map(lum_c); - // empty arrays (mie scattering not supported in lite version) - Array_gpu mie_cdfs_sub; - Array_gpu mie_angs_sub; - Array_gpu rel; ////// 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); - output_nc.add_dimension("x", nx); - output_nc.add_dimension("y", ny); - output_nc.add_dimension("z", n_z_in); - + if (switch_raytracing) + { + output_nc.add_dimension("x", nx); + output_nc.add_dimension("y", ny); + output_nc.add_dimension("z", n_z_in); + } + if (switch_bw_raytracing) + { + output_nc.add_dimension("nx", camera.nx); + output_nc.add_dimension("ny", camera.ny); + } + //// GPU arrays Array_gpu tot_tau_g(tot_tau); Array_gpu tot_ssa_g(tot_ssa); @@ -244,101 +346,185 @@ void solve_radiation(int argc, char** argv) Array_gpu aer_asy_g(aer_asy); Array_gpu sfc_albedo_g(sfc_albedo); - Raytracer raytracer; + if (switch_raytracing) + { + Raytracer raytracer; + const Vector grid_cells = {nx, ny, nz}; - // Solve the radiation. - Status::print_message("Starting the raytracer!!"); + // Solve the radiation. + Status::print_message("Starting the raytracer!!"); - auto run_solver = [&]() - { - cudaDeviceSynchronize(); - cudaEvent_t start; - cudaEvent_t stop; - cudaEventCreate(&start); - cudaEventCreate(&stop); - - cudaEventRecord(start, 0); - // do something. + auto run_solver = [&]() + { + cudaDeviceSynchronize(); + cudaEvent_t start; + cudaEvent_t stop; + cudaEventCreate(&start); + cudaEventCreate(&stop); + + cudaEventRecord(start, 0); + // do something. + + raytracer.trace_rays( + 0, + switch_independent_column, + photons_per_pixel, + grid_cells, + grid_d, + kn_grid, + mie_cdfs, + mie_angs, + tot_tau_g, + tot_ssa_g, + cld_tau_g, + cld_ssa_g, + cld_asy_g, + aer_tau_g, + aer_ssa_g, + aer_asy_g, + rel, + sfc_albedo_g, + zenith_angle, + azimuth_angle, + tod_dir * std::cos(zenith_angle), + Float(0.), + flux_tod_dn, + flux_tod_up, + flux_sfc_dir, + flux_sfc_dif, + flux_sfc_up, + flux_abs_dir, + flux_abs_dif); + + cudaEventRecord(stop, 0); + cudaEventSynchronize(stop); + float duration = 0.f; + cudaEventElapsedTime(&duration, start, stop); + + cudaEventDestroy(start); + cudaEventDestroy(stop); + + Status::print_message("Duration raytracer: " + std::to_string(duration) + " (ms)"); + }; + + // Tuning step; + run_solver(); + + //// Profiling step; + if (switch_profiling) + { + cudaProfilerStart(); + run_solver(); + cudaProfilerStop(); + } - raytracer.trace_rays( - 0, - switch_independent_column, - photons_per_pixel, - grid_cells, - grid_d, - kn_grid, - mie_cdfs_sub, - mie_angs_sub, - tot_tau_g, - tot_ssa_g, - cld_tau_g, - cld_ssa_g, - cld_asy_g, - aer_tau_g, - aer_ssa_g, - aer_asy_g, - rel, - sfc_albedo_g, - zenith_angle, - azimuth_angle, - tod_dir * std::cos(zenith_angle), - Float(0.), - flux_tod_dn, - flux_tod_up, - flux_sfc_dir, - flux_sfc_dif, - flux_sfc_up, - flux_abs_dir, - flux_abs_dif); - - cudaEventRecord(stop, 0); - cudaEventSynchronize(stop); - float duration = 0.f; - cudaEventElapsedTime(&duration, start, stop); - - cudaEventDestroy(start); - cudaEventDestroy(stop); - - Status::print_message("Duration raytracer: " + std::to_string(duration) + " (ms)"); - }; - - // Tuning step; - run_solver(); - - //// Profiling step; - if (switch_profiling) + // output arrays to cpu + Array flux_tod_dn_c(flux_tod_dn); + Array flux_tod_up_c(flux_tod_up); + Array flux_sfc_dir_c(flux_sfc_dir); + Array flux_sfc_dif_c(flux_sfc_dif); + Array flux_sfc_up_c(flux_sfc_up); + Array flux_abs_dir_c(flux_abs_dir); + Array flux_abs_dif_c(flux_abs_dif); + // Store the output. + Status::print_message("Storing the raytracer output."); + + auto nc_flux_tod_dn = output_nc.add_variable("flux_tod_dn" , {"y", "x"}); + auto nc_flux_tod_up = output_nc.add_variable("flux_tod_up" , {"y", "x"}); + auto nc_flux_sfc_dir = output_nc.add_variable("flux_sfc_dir", {"y", "x"}); + auto nc_flux_sfc_dif = output_nc.add_variable("flux_sfc_dif", {"y", "x"}); + auto nc_flux_sfc_up = output_nc.add_variable("flux_sfc_up" , {"y", "x"}); + auto nc_flux_abs_dir = output_nc.add_variable("abs_dir" , {"z", "y", "x"}); + auto nc_flux_abs_dif = output_nc.add_variable("abs_dif" , {"z", "y", "x"}); + + nc_flux_tod_dn .insert(flux_tod_dn_c .v(), {0, 0}); + nc_flux_tod_up .insert(flux_tod_up_c .v(), {0, 0}); + nc_flux_sfc_dir .insert(flux_sfc_dir_c .v(), {0, 0}); + nc_flux_sfc_dif .insert(flux_sfc_dif_c .v(), {0, 0}); + nc_flux_sfc_up .insert(flux_sfc_up_c .v(), {0, 0}); + nc_flux_abs_dir .insert(flux_abs_dir_c .v(), {0, 0, 0}); + nc_flux_abs_dif .insert(flux_abs_dif_c .v(), {0, 0, 0}); + + } + + if (switch_bw_raytracing) { - cudaProfilerStart(); + Raytracer_bw raytracer_bw; + const Vector grid_cells = {nx, ny, n_z_in}; + + // Solve the radiation. + Status::print_message("Starting the backward raytracer!!"); + + auto run_solver = [&]() + { + cudaDeviceSynchronize(); + cudaEvent_t start; + cudaEvent_t stop; + cudaEventCreate(&start); + cudaEventCreate(&stop); + + cudaEventRecord(start, 0); + + raytracer_bw.trace_rays_bb( + 0, + photons_per_pixel_bw, n_lay, + grid_cells, grid_d, kn_grid, + z_lev, + mie_cdfs_bw, + mie_angs_bw, + mie_phase_bw, + mie_phase_angs_bw, + rel, + tot_tau_g, + tot_ssa_g, + cld_tau_g, + cld_ssa_g, + cld_asy_g, + aer_tau_g, + aer_ssa_g, + aer_asy_g, + sfc_albedo_g, + land_use_map, + zenith_angle, + azimuth_angle, + tod_dir, + camera, + radiance); + + cudaEventRecord(stop, 0); + cudaEventSynchronize(stop); + float duration = 0.f; + cudaEventElapsedTime(&duration, start, stop); + + cudaEventDestroy(start); + cudaEventDestroy(stop); + + Status::print_message("Duration bw raytracer: " + std::to_string(duration) + " (ms)"); + }; + + // Tuning step; run_solver(); - cudaProfilerStop(); - } - // output arrays to cpu - Array flux_tod_dn_c(flux_tod_dn); - Array flux_tod_up_c(flux_tod_up); - Array flux_sfc_dir_c(flux_sfc_dir); - Array flux_sfc_dif_c(flux_sfc_dif); - Array flux_sfc_up_c(flux_sfc_up); - Array flux_abs_dir_c(flux_abs_dir); - Array flux_abs_dif_c(flux_abs_dif); - // Store the output. - Status::print_message("Storing the raytracer output."); - - auto nc_flux_tod_dn = output_nc.add_variable("flux_tod_dn" , {"y", "x"}); - auto nc_flux_tod_up = output_nc.add_variable("flux_tod_up" , {"y", "x"}); - auto nc_flux_sfc_dir = output_nc.add_variable("flux_sfc_dir", {"y", "x"}); - auto nc_flux_sfc_dif = output_nc.add_variable("flux_sfc_dif", {"y", "x"}); - auto nc_flux_sfc_up = output_nc.add_variable("flux_sfc_up" , {"y", "x"}); - auto nc_flux_abs_dir = output_nc.add_variable("abs_dir" , {"z", "y", "x"}); - auto nc_flux_abs_dif = output_nc.add_variable("abs_dif" , {"z", "y", "x"}); - - nc_flux_tod_dn .insert(flux_tod_dn_c .v(), {0, 0}); - nc_flux_tod_up .insert(flux_tod_up_c .v(), {0, 0}); - nc_flux_sfc_dir .insert(flux_sfc_dir_c .v(), {0, 0}); - nc_flux_sfc_dif .insert(flux_sfc_dif_c .v(), {0, 0}); - nc_flux_sfc_up .insert(flux_sfc_up_c .v(), {0, 0}); - nc_flux_abs_dir .insert(flux_abs_dir_c .v(), {0, 0, 0}); - nc_flux_abs_dif .insert(flux_abs_dif_c .v(), {0, 0, 0}); + //// Profiling step; + if (switch_profiling) + { + cudaProfilerStart(); + run_solver(); + cudaProfilerStop(); + } + + // output arrays to cpu + Array radiance_c(radiance); + // Store the output. + Status::print_message("Storing the bw raytracer output."); + + auto nc_radiance = output_nc.add_variable("radiance" , {"ny", "nx"}); + + nc_radiance.insert(radiance_c .v(), {0, 0}); + + + + } Status::print_message("###### Finished RAYTRACING #####"); } From 526ece539c954fd9241c8d2b70c257b9430bb868 Mon Sep 17 00:00:00 2001 From: Menno Date: Tue, 5 Nov 2024 09:29:20 +0100 Subject: [PATCH 4/5] improve bw camera rotation, now it fully allows all 3D rotation. --- include_rt_kernels/raytracer_kernels_bw.h | 21 +++++++++------------ python/set_virtual_camera.py | 10 +++++++--- 2 files changed, 16 insertions(+), 15 deletions(-) diff --git a/include_rt_kernels/raytracer_kernels_bw.h b/include_rt_kernels/raytracer_kernels_bw.h index 3663237..bc9d63e 100644 --- a/include_rt_kernels/raytracer_kernels_bw.h +++ b/include_rt_kernels/raytracer_kernels_bw.h @@ -59,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 dir_tmp = {0, 0, 1}; - const Vector dir_cam = normalize(Vector({dot(camera.mx,dir_tmp), dot(camera.my,dir_tmp), dot(camera.mz,dir_tmp)*Float(-1)})); - Vector 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 dir_tmp = {1, 0, 0}; + const Vector dir_cam = normalize(Vector({dot(camera.mx,dir_tmp), dot(camera.my,dir_tmp), dot(camera.mz,dir_tmp)})); + const Vector dir_up_tmp = {0, 0, 1}; + const Vector dir_up = normalize(Vector({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.)); diff --git a/python/set_virtual_camera.py b/python/set_virtual_camera.py index 849eddc..79cf166 100644 --- a/python/set_virtual_camera.py +++ b/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() From e018ced04ae2fb5cbc77fe80a973c205fc64b4ea Mon Sep 17 00:00:00 2001 From: Menno Veerman Date: Wed, 20 Nov 2024 18:38:34 +0100 Subject: [PATCH 5/5] don't limit effective radius to 21.5, that severely overestimates ice optical depth! --- src_cuda_rt/Cloud_optics_rt.cu | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src_cuda_rt/Cloud_optics_rt.cu b/src_cuda_rt/Cloud_optics_rt.cu index c5a3f90..f42e0aa 100644 --- a/src_cuda_rt/Cloud_optics_rt.cu +++ b/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]));