Skip to content

Commit

Permalink
Compiling code and excessive allocations for rt absorbance fields.
Browse files Browse the repository at this point in the history
  • Loading branch information
@Chiil
Chiil committed Oct 24, 2024
1 parent dbdf1cd commit 02c1caf
Show file tree
Hide file tree
Showing 5 changed files with 93 additions and 135 deletions.
4 changes: 2 additions & 2 deletions include_rt/Fluxes_rt.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -58,7 +58,7 @@ class Fluxes_rt
class Fluxes_broadband_rt : public Fluxes_rt
{
public:
Fluxes_broadband_rt(const int ncol_x, const int ncol_y, const int nlev);
Fluxes_broadband_rt(const int ncol_x, const int ncol_y, const int n_z, const int nlev);
virtual ~Fluxes_broadband_rt() {};

virtual void net_flux();
Expand Down Expand Up @@ -113,7 +113,7 @@ class Fluxes_broadband_rt : public Fluxes_rt
class Fluxes_byband_rt : public Fluxes_broadband_rt
{
public:
Fluxes_byband_rt(const int ncol_x, const int ncol_y, const int nlev, const int nbnd);
Fluxes_byband_rt(const int ncol_x, const int ncol_y, const int n_z, const int nlev, const int nbnd);
virtual ~Fluxes_byband_rt() {};

virtual void reduce(
Expand Down
65 changes: 6 additions & 59 deletions src_cuda_rt/Fluxes_rt.cu
View file
Original file line number Diff line number Diff line change
Expand Up @@ -116,61 +116,8 @@ namespace
}
}

//namespace rrtmgp_kernel_launcher
//{
//
// void sum_broadband(
// int ncol, int nlev, int ngpt,
// const Array<Float,3>& spectral_flux, Array<Float,2>& broadband_flux)
// {
// rrtmgp_kernels::sum_broadband(
// &ncol, &nlev, &ngpt,
// const_cast<Float*>(spectral_flux.ptr()),
// broadband_flux.ptr());
// }
//
//
// void net_broadband(
// int ncol, int nlev,
// const Array<Float,2>& broadband_flux_dn, const Array<Float,2>& broadband_flux_up,
// Array<Float,2>& broadband_flux_net)
// {
// rrtmgp_kernels::net_broadband_precalc(
// &ncol, &nlev,
// const_cast<Float*>(broadband_flux_dn.ptr()),
// const_cast<Float*>(broadband_flux_up.ptr()),
// broadband_flux_net.ptr());
// }
//
//
// void sum_byband(
// int ncol, int nlev, int ngpt, int nbnd,
// const Array<int,2>& band_lims,
// const Array<Float,3>& spectral_flux,
// Array<Float,3>& byband_flux)
// {
// rrtmgp_kernels::sum_byband(
// &ncol, &nlev, &ngpt, &nbnd,
// const_cast<int*>(band_lims.ptr()),
// const_cast<Float*>(spectral_flux.ptr()),
// byband_flux.ptr());
// }
//
//
// void net_byband(
// int ncol, int nlev, int nband,
// const Array<Float,3>& byband_flux_dn, const Array<Float,3>& byband_flux_up,
// Array<Float,3>& byband_flux_net)
// {
// rrtmgp_kernels::net_byband_precalc(
// &ncol, &nlev, &nband,
// const_cast<Float*>(byband_flux_dn.ptr()),
// const_cast<Float*>(byband_flux_up.ptr()),
// byband_flux_net.ptr());
// }


Fluxes_broadband_rt::Fluxes_broadband_rt(const int ncol_x, const int ncol_y, const int nlev) :

Fluxes_broadband_rt::Fluxes_broadband_rt(const int ncol_x, const int ncol_y, const int n_z, const int nlev) :
flux_up ({ncol_x*ncol_y, nlev}),
flux_dn ({ncol_x*ncol_y, nlev}),
flux_dn_dir ({ncol_x*ncol_y, nlev}),
Expand All @@ -180,8 +127,8 @@ Fluxes_broadband_rt::Fluxes_broadband_rt(const int ncol_x, const int ncol_y, con
flux_sfc_dir({ncol_x, ncol_y}),
flux_sfc_dif({ncol_x, ncol_y}),
flux_sfc_up ({ncol_x, ncol_y}),
flux_abs_dir({ncol_x, ncol_y, nlev-1}),
flux_abs_dif({ncol_x, ncol_y, nlev-1})
flux_abs_dir({ncol_x, ncol_y, n_z}),
flux_abs_dif({ncol_x, ncol_y, n_z})
{}


Expand Down Expand Up @@ -259,8 +206,8 @@ void Fluxes_broadband_rt::reduce(
}


Fluxes_byband_rt::Fluxes_byband_rt(const int ncol_x, const int ncol_y, const int nlev, const int nbnd) :
Fluxes_broadband_rt(ncol_x, ncol_y, nlev),
Fluxes_byband_rt::Fluxes_byband_rt(const int ncol_x, const int ncol_y, const int n_z, const int nlev, const int nbnd) :
Fluxes_broadband_rt(ncol_x, ncol_y, n_z, nlev),
bnd_flux_up ({ncol_x * ncol_y, nlev, nbnd}),
bnd_flux_dn ({ncol_x * ncol_y, nlev, nbnd}),
bnd_flux_dn_dir({ncol_x * ncol_y, nlev, nbnd}),
Expand Down
89 changes: 46 additions & 43 deletions src_test/Radiation_solver_bw.cu
View file
Original file line number Diff line number Diff line change
Expand Up @@ -508,6 +508,9 @@ void Radiation_solver_longwave::solve_gpu(
Array_gpu<Float,2>& lw_flux_up, Array_gpu<Float,2>& lw_flux_dn, Array_gpu<Float,2>& lw_flux_net,
Array_gpu<Float,3>& lw_bnd_flux_up, Array_gpu<Float,3>& lw_bnd_flux_dn, Array_gpu<Float,3>& lw_bnd_flux_net)
{
throw std::runtime_error("Longwave raytracing is not implemented");

/*
const int n_col = p_lay.dim(1);
const int n_lay = p_lay.dim(2);
const int n_lev = p_lev.dim(2);
Expand Down Expand Up @@ -624,10 +627,10 @@ void Radiation_solver_longwave::solve_gpu(
}
}
}
*/
}



Float get_x(const Float wv)
{
const Float a = (wv - Float(442.0)) * ((wv < Float(442.0)) ? Float(0.0624) : Float(0.0374));
Expand Down Expand Up @@ -744,10 +747,10 @@ void Radiation_solver_shortwave::load_mie_tables(
const int n_bnd_sw = this->get_n_bnd_gpu();
const int n_re = mie_nc.get_dimension_size("r_eff");
const int n_mie = mie_nc.get_dimension_size("n_ang");

Array<Float,3> mie_cdf(mie_nc.get_variable<Float>("phase_cdf", {n_bnd_sw, n_mie}), {n_mie, 1, n_bnd_sw});
Array<Float,4> mie_ang(mie_nc.get_variable<Float>("phase_cdf_angle", {n_bnd_sw, n_re, n_mie}), {n_mie, n_re, 1, n_bnd_sw});

Array<Float,4> mie_phase(mie_nc.get_variable<Float>("phase", {n_bnd_sw, n_re, n_mie}), {n_mie, n_re, 1, n_bnd_sw});
Array<Float,1> mie_phase_ang(mie_nc.get_variable<Float>("phase_angle", {n_mie}), {n_mie});

Expand All @@ -757,15 +760,15 @@ void Radiation_solver_shortwave::load_mie_tables(
this->mie_phase_bb = mie_phase;
this->mie_phase_angs_bb = mie_phase_ang;
}

if (switch_image)
{
Netcdf_file mie_nc(file_name_mie_im, Netcdf_mode::Read);
const int n_bnd_sw = this->get_n_bnd_gpu();
const int n_re = mie_nc.get_dimension_size("r_eff");
const int n_mie = mie_nc.get_dimension_size("n_ang");
const int n_sub = mie_nc.get_dimension_size("sub_band");

Array<Float,3> mie_cdf(mie_nc.get_variable<Float>("phase_cdf", {n_bnd_sw, n_sub, n_mie}), {n_mie, n_sub, n_bnd_sw});
Array<Float,4> mie_ang(mie_nc.get_variable<Float>("phase_cdf_angle", {n_bnd_sw, n_sub, n_re, n_mie}), {n_mie, n_re, n_sub, n_bnd_sw});

Expand All @@ -777,7 +780,7 @@ void Radiation_solver_shortwave::load_mie_tables(

this->mie_phase = mie_phase;
this->mie_phase_angs = mie_phase_ang;

}

}
Expand Down Expand Up @@ -874,13 +877,13 @@ void Radiation_solver_shortwave::solve_gpu(
}
Array_gpu<Float,2> albedo;
if (!switch_lu_albedo) albedo = sfc_alb.subset({{ {band, band}, {1, n_col}}});

if (!tune_step && (! (band == 10 || band == 11 || band ==12))) continue;

const Float solar_source_band = kdist_gpu->band_source(band_limits_gpt({1,band}), band_limits_gpt({2,band}));

constexpr int n_col_block = 1<<14; // 2^14

Array_gpu<Float,1> toa_src_temp({n_col_block});
auto gas_optics_subset = [&](
const int col_s, const int n_col_subset)
Expand All @@ -899,27 +902,27 @@ void Radiation_solver_shortwave::solve_gpu(
toa_src_temp,
col_dry);
};

const int n_blocks = n_col / n_col_block;
const int n_col_residual = n_col % n_col_block;

std::unique_ptr<Optical_props_arry_rt> optical_props_block =
std::make_unique<Optical_props_2str_rt>(n_col_block, n_lay, *kdist_gpu);

for (int n=0; n<n_blocks; ++n)
{
const int col_s = n*n_col_block;
gas_optics_subset(col_s, n_col_block);
}

if (n_col_residual > 0)
{
const int col_s = n_blocks*n_col_block;
gas_optics_subset(col_s, n_col_residual);
}

if (tune_step) return;

toa_src.fill(toa_src_temp({1}) * tsi_scaling({1}));
if (switch_cloud_optics)
{
Expand All @@ -932,11 +935,11 @@ void Radiation_solver_shortwave::solve_gpu(
rel,
rei,
*cloud_optical_props);

if (switch_delta_cloud)
cloud_optical_props->delta_scale();
}

// Add the cloud optical props to the gas optical properties.
add_to(
dynamic_cast<Optical_props_2str_rt&>(*optical_props),
Expand All @@ -948,7 +951,7 @@ void Radiation_solver_shortwave::solve_gpu(
Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_col, n_lay, cloud_optical_props->get_ssa().ptr());
Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_col, n_lay, cloud_optical_props->get_g().ptr());
}

if (switch_aerosol_optics)
{
if (band > previous_band)
Expand All @@ -959,11 +962,11 @@ void Radiation_solver_shortwave::solve_gpu(
aerosol_concs_subset,
rh, p_lev,
*aerosol_optical_props);

if (switch_delta_aerosol)
aerosol_optical_props->delta_scale();
}

// Add the aerosol optical props to the gas optical properties.
add_to(
dynamic_cast<Optical_props_2str_rt&>(*optical_props),
Expand All @@ -975,61 +978,61 @@ void Radiation_solver_shortwave::solve_gpu(
Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_col, n_lay, aerosol_optical_props->get_ssa().ptr());
Gas_optics_rrtmgp_kernels_cuda_rt::zero_array(n_col, n_lay, aerosol_optical_props->get_g().ptr());
}


/* rrtmgp's bands are quite broad, we divide each spectral band in three equally broad spectral intervals
and run each g-point for each spectral interval, using the mean rayleigh scattering coefficient of each spectral interval
in stead of RRTMGP's rayleigh scattering coefficients.
The contribution of each spectral interval to the spectral band is based on the integrated (<>) Planck source function:
<Planck(spectral interval)> / <Planck(spectral band)>, with a sun temperature of 5778 K. This is not entirely accurate because
the sun is not a black body radiatior, but the approximations comes close enough.
*/

// number of intervals
const Array<Float, 2>& band_limits_wn(this->kdist_gpu->get_band_lims_wavenumber());
const int nwv = n_sub;

const Float wv1 = 1. / band_limits_wn({2,band}) * Float(1.e7);
const Float wv2 = 1. / band_limits_wn({1,band}) * Float(1.e7);
const Float dwv = (wv2-wv1)/Float(nwv);

//
const Float total_planck = Planck_integrator(wv1,wv2);

for (int iwv=0; iwv<nwv; ++iwv)
{
const Float wv1_sub = wv1 + iwv*dwv;
const Float wv2_sub = wv1 + (iwv+1)*dwv;
const Float local_planck = Planck_integrator(wv1_sub,wv2_sub);

// use RRTMGPs scattering coefficients if solving per band instead of subbands
const Float rayleigh = (nwv==1) ? 0 : rayleigh_mean(wv1_sub, wv2_sub);
const Float toa_factor = local_planck / total_planck * Float(1.)/solar_source_band;

// XYZ factors
Array<Float,1> xyz_factor({3});
xyz_factor({1}) = xyz_irradiance(wv1_sub,wv2_sub,&get_x);
xyz_factor({2}) = xyz_irradiance(wv1_sub,wv2_sub,&get_y);
xyz_factor({3}) = xyz_irradiance(wv1_sub,wv2_sub,&get_z);
Array_gpu<Float,1> xyz_factor_gpu(xyz_factor);

if (switch_lu_albedo)
{
if (albedo.size() == 0) albedo.set_dims({1, n_col});
spectral_albedo(n_col, wv1, wv2, land_use_map, albedo);
}

const Float zenith_angle = std::acos(mu0({1}));
const Float azimuth_angle = azi({1});

if (switch_cloud_mie)
{
mie_cdfs_sub = mie_cdfs.subset({{ {1, n_mie}, {iwv+1,iwv+1}, {band, band} }});
mie_angs_sub = mie_angs.subset({{ {1, n_mie}, {1, n_re}, {iwv+1,iwv+1}, {band, band} }});
mie_phase_sub = mie_phase.subset({{ {1, n_mie}, {1, n_re}, {iwv+1,iwv+1}, {band, band} }});
}

raytracer.trace_rays(
igpt,
photons_per_pixel, n_lay,
Expand Down Expand Up @@ -1059,18 +1062,18 @@ void Radiation_solver_shortwave::solve_gpu(
gas_concs.get_vmr("h2o"),
camera,
flux_camera);

raytracer.add_xyz_camera(
camera,
xyz_factor_gpu,
flux_camera,
XYZ);

}
previous_band = band;
}
}
}

if (switch_cloud_cam)
{
cloud_optics_gpu->cloud_optics(
Expand All @@ -1080,13 +1083,13 @@ void Radiation_solver_shortwave::solve_gpu(
rel,
rei,
*cloud_optical_props);

raytracer.accumulate_clouds(
grid_d,
grid_cells,
lwp,
iwp,
dynamic_cast<Optical_props_2str_rt&>(*cloud_optical_props).get_tau(),
dynamic_cast<Optical_props_2str_rt&>(*cloud_optical_props).get_tau(),
camera,
liwp_cam,
tauc_cam,
Expand Down Expand Up @@ -1136,7 +1139,7 @@ void Radiation_solver_shortwave::solve_gpu_bb(

const int n_mie = (switch_cloud_mie) ? this->mie_angs_bb.dim(1) : 0;
const int n_re = (switch_cloud_mie) ? this->mie_angs_bb.dim(2) : 0;

const int cam_nx = radiance.dim(1);
const int cam_ny = radiance.dim(2);
const Bool top_at_1 = p_lay({1, 1}) < p_lay({1, n_lay});
Expand Down Expand Up @@ -1342,13 +1345,13 @@ void Radiation_solver_shortwave::solve_gpu_bb(
rel,
rei,
*cloud_optical_props);

raytracer.accumulate_clouds(
grid_d,
grid_cells,
lwp,
iwp,
dynamic_cast<Optical_props_2str_rt&>(*cloud_optical_props).get_tau(),
dynamic_cast<Optical_props_2str_rt&>(*cloud_optical_props).get_tau(),
camera,
liwp_cam,
tauc_cam,
Expand Down
Loading

0 comments on commit 02c1caf

Please sign in to comment.