htrdr_planets_solve_volrad_budget.c (15163B)
1 /* Copyright (C) 2018-2019, 2022-2025 Centre National de la Recherche Scientifique 2 * Copyright (C) 2020-2022 Institut Mines Télécom Albi-Carmaux 3 * Copyright (C) 2022-2025 Institut Pierre-Simon Laplace 4 * Copyright (C) 2022-2025 Institut de Physique du Globe de Paris 5 * Copyright (C) 2018-2025 |Méso|Star> (contact@meso-star.com) 6 * Copyright (C) 2022-2025 Observatoire de Paris 7 * Copyright (C) 2022-2025 Université de Reims Champagne-Ardenne 8 * Copyright (C) 2022-2025 Université de Versaille Saint-Quentin 9 * Copyright (C) 2018-2019, 2022-2025 Université Paul Sabatier 10 * 11 * This program is free software: you can redistribute it and/or modify 12 * it under the terms of the GNU General Public License as published by 13 * the Free Software Foundation, either version 3 of the License, or 14 * (at your option) any later version. 15 * 16 * This program is distributed in the hope that it will be useful, 17 * but WITHOUT ANY WARRANTY; without even the implied warranty of 18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 19 * GNU General Public License for more details. 20 * 21 * You should have received a copy of the GNU General Public License 22 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 23 24 #include "planets/htrdr_planets_c.h" 25 #include "planets/htrdr_planets_source.h" 26 27 #include "core/htrdr_accum.h" 28 #include "core/htrdr_log.h" 29 #include "core/htrdr_ran_wlen_discrete.h" 30 #include "core/htrdr_ran_wlen_planck.h" 31 #include "core/htrdr_solve_buffer.h" 32 33 #include <star/smsh.h> 34 #include <star/ssp.h> 35 36 #include <rsys/clock_time.h> 37 38 /******************************************************************************* 39 * Helper functions 40 ******************************************************************************/ 41 static void 42 spectral_sampling 43 (struct htrdr_planets* cmd, 44 const struct htrdr_solve_item_args* args, 45 double* out_wlen, /* Sampled wavelength [nm] */ 46 double* out_wlen_pdf, /* [nm^-1] */ 47 size_t* out_iband, /* Spectral band in which the sampled wavelength falls */ 48 size_t* out_iquad) /* Sampled quadrature point in the spectral band */ 49 { 50 size_t iband_range[2] = {0, 0}; 51 size_t iband = 0; 52 size_t iquad = 0; 53 double r0, r1, r2; /* Random Numbers */ 54 double wlen[2] = {0,0}; /* [nm] */ 55 double pdf = 0; /* [nm^1] */ 56 57 /* Preconditions */ 58 ASSERT(cmd && args && out_wlen && out_wlen_pdf && out_iband && out_iquad); 59 60 r0 = ssp_rng_canonical(args->rng); 61 r1 = ssp_rng_canonical(args->rng); 62 r2 = ssp_rng_canonical(args->rng); 63 64 /* Sample a wavelength with respect to the type of spectral integration */ 65 switch(cmd->spectral_domain.type) { 66 /* Longwave */ 67 case HTRDR_SPECTRAL_LW: 68 wlen[0] = htrdr_ran_wlen_planck_sample(cmd->planck, r0, r1, &pdf); 69 break; 70 /* Shortwave */ 71 case HTRDR_SPECTRAL_SW: 72 if(htrdr_planets_source_does_radiance_vary_spectrally(cmd->source)) { 73 wlen[0] = htrdr_ran_wlen_discrete_sample(cmd->discrete, r0, r1, &pdf); 74 } else { 75 wlen[0] = htrdr_ran_wlen_planck_sample(cmd->planck, r0, r1, &pdf); 76 } 77 break; 78 default: FATAL("Unreachable code\n"); break; 79 } 80 wlen[1] = wlen[0]; 81 82 /* Find the band the wavelength belongs to */ 83 RNATM(find_bands(cmd->atmosphere, wlen, iband_range)); 84 ASSERT(iband_range[0] == iband_range[1]); 85 iband = iband_range[0]; 86 87 /* Sample a quadrature point */ 88 RNATM(band_sample_quad_pt(cmd->atmosphere, r2, iband, &iquad)); 89 90 *out_wlen = wlen[0]; 91 *out_wlen_pdf = pdf; 92 *out_iband = iband; 93 *out_iquad = iquad; 94 } 95 96 static INLINE void 97 position_sampling 98 (const struct htrdr_solve_item_args* args, 99 const struct smsh_desc* desc, 100 double pos[3]) 101 { 102 const double* v0 = NULL; 103 const double* v1 = NULL; 104 const double* v2 = NULL; 105 const double* v3 = NULL; 106 107 /* Preconditions */ 108 ASSERT(args && desc && pos); 109 110 /* Retrieve the vertices of the tetrahedron */ 111 v0 = desc->nodes + desc->cells[args->item_id*4/*#vertices*/+0]*3/*#coords*/; 112 v1 = desc->nodes + desc->cells[args->item_id*4/*#vertices*/+1]*3/*#coords*/; 113 v2 = desc->nodes + desc->cells[args->item_id*4/*#vertices*/+2]*3/*#coords*/; 114 v3 = desc->nodes + desc->cells[args->item_id*4/*#vertices*/+3]*3/*#coords*/; 115 116 ssp_ran_tetrahedron_uniform(args->rng, v0, v1, v2, v3, pos, NULL/*pdf*/); 117 } 118 119 static double /* [W/m^2/sr/m] */ 120 get_source 121 (struct htrdr_planets* cmd, 122 const double pos[3], 123 const double wlen) /* [nm] */ 124 { 125 struct rnatm_cell_pos cell_pos = RNATM_CELL_POS_NULL; 126 double temperature = 0; /* [K] */ 127 double source = 0; /* [W/m^2/sr/m] */ 128 const double wlen_m = wlen * 1.e-9; /* Wavelength [m] */ 129 ASSERT(cmd && pos); 130 131 switch(cmd->spectral_domain.type) { 132 case HTRDR_SPECTRAL_SW: 133 /* In shortwave, the source is external to the system */ 134 source = 0; /* [W/m^2/sr/m] */ 135 break; 136 137 case HTRDR_SPECTRAL_LW: 138 RNATM(fetch_cell(cmd->atmosphere, pos, RNATM_GAS, &cell_pos)); 139 140 if(SUVM_PRIMITIVE_NONE(&cell_pos.prim)) { 141 /* The position is not in the gas */ 142 source = 0; /* [W/m^2/sr/m] */ 143 144 } else { 145 /* Fetch the source temperature */ 146 RNATM(cell_get_gas_temperature(cmd->atmosphere, &cell_pos, &temperature)); 147 source = htrdr_planck_monochromatic(wlen_m, temperature); /* [W/m^2/sr/m] */ 148 } 149 break; 150 151 default: FATAL("Unreachable code\n"); break; 152 } 153 154 return source; /* [W/m^2/sr/m] */ 155 } 156 157 /* Return the total absorption coefficient, 158 * i.e. the sum of the gas and aerosol ka */ 159 static double 160 get_ka 161 (struct htrdr_planets* cmd, 162 const double pos[3], 163 const size_t iband, /* Spectral band */ 164 const size_t iquad) /* Quadrature point */ 165 { 166 struct rnatm_cell_pos cells[RNATM_MAX_COMPONENTS_COUNT]; 167 struct rnatm_get_radcoef_args get_k_args = RNATM_GET_RADCOEF_ARGS_NULL; 168 double ka = 0; 169 170 ASSERT(cmd && pos); 171 172 get_k_args.cells = cells; 173 get_k_args.iband = iband; 174 get_k_args.iquad = iquad; 175 get_k_args.radcoef = RNATM_RADCOEF_Ka; 176 177 /* Retrieve the list of aerosol and gas cells in which pos lies */ 178 RNATM(fetch_cell_list(cmd->atmosphere, pos, get_k_args.cells, NULL)); 179 180 /* Retrive the total absorption coefficient */ 181 RNATM(get_radcoef(cmd->atmosphere, &get_k_args, &ka)); 182 183 return ka; 184 } 185 186 static void 187 realisation 188 (struct htrdr_planets* cmd, 189 const struct htrdr_solve_item_args* args, 190 const struct smsh_desc* volrad_mesh_desc, 191 double weights[PLANETS_VOLRAD_WEIGHTS_COUNT]) /* [W/m^3] */ 192 { 193 struct planets_compute_radiance_args rad_args = 194 PLANETS_COMPUTE_RADIANCE_ARGS_NULL; 195 196 /* Spectral integration */ 197 double wlen = 0; /* Wavelength [nm] */ 198 double wlen_pdf_nm = 0; /* Wavelength pdf [nm^-1] */ 199 double wlen_pdf_m = 0; /* Wavelength pdf [m^-1] */ 200 size_t iband = 0; /* Spectral band */ 201 size_t iquad = 0; /* Quadrature point */ 202 203 /* Spatial & angular integration */ 204 double dir_src[3] = {0,0,0}; /* Direction toward the source */ 205 double dir[3] = {0,0,0}; 206 double pos[3] = {0,0,0}; 207 double dir_src_pdf = 0; 208 double dir_pdf = 0; 209 210 double S = 0; /* Source [W/m^2/sr/m] */ 211 double L_direct = 0; /* Direct radiance [W/m^2/sr/m] */ 212 double L_diffuse = 0; /* Diffuse radiance [W/m^2/sr/m] */ 213 double ka = 0; /* Absorption coefficient */ 214 215 /* Preconditions */ 216 ASSERT(cmd && args && volrad_mesh_desc); 217 ASSERT(cmd->output_type == HTRDR_PLANETS_ARGS_OUTPUT_VOLUMIC_RADIATIVE_BUDGET); 218 219 /* Initialise the weights */ 220 memset(weights, 0, sizeof(double)*PLANETS_VOLRAD_WEIGHTS_COUNT); 221 222 spectral_sampling(cmd, args, &wlen, &wlen_pdf_nm, &iband, &iquad); 223 position_sampling(args, volrad_mesh_desc, pos); 224 ssp_ran_sphere_uniform(args->rng, dir, &dir_pdf); 225 226 S = get_source(cmd, pos, wlen); /* [W/m^2/sr/m] */ 227 228 ka = get_ka(cmd, pos, iband, iquad); 229 wlen_pdf_m = wlen_pdf_nm * 1.e9; /* Transform pdf from nm^-1 to m^-1 */ 230 231 /* Compute the radiance in W/m^2/sr/m */ 232 d3_set(rad_args.path_org, pos); 233 rad_args.rng = args->rng; 234 rad_args.ithread = args->ithread; 235 rad_args.wlen = wlen; /* [nm] */ 236 rad_args.iband = iband; 237 rad_args.iquad = iquad; 238 239 if(cmd->spectral_domain.type == HTRDR_SPECTRAL_LW) { 240 /* In the longwave (radiation due to the medium), simply sample a radiative 241 * path for the sampled direction and position: the radiance is considered 242 * as purely diffuse. */ 243 d3_set(rad_args.path_dir, dir); 244 L_diffuse = planets_compute_radiance(cmd, &rad_args); /* [W/m^2/sr/m] */ 245 246 /* Calculate the weights [W/m^3] */ 247 weights[PLANETS_VOLRAD_DIRECT] = 0.0; 248 weights[PLANETS_VOLRAD_DIFFUSE] = ka * (L_diffuse - S) / (wlen_pdf_m * dir_pdf); 249 250 } else { 251 /* In the so-called shortwave region (actually, the radiation due the 252 * external source) is decomposed in its direct and diffuse components */ 253 254 dir_src_pdf = htrdr_planets_source_sample_direction 255 (cmd->source, args->rng, pos, dir_src); 256 257 d3_set(rad_args.path_dir, dir_src); 258 rad_args.component = PLANETS_RADIANCE_CPNT_DIRECT; 259 L_direct = planets_compute_radiance(cmd, &rad_args); /* [W/m^2/sr/m] */ 260 261 d3_set(rad_args.path_dir, dir); 262 rad_args.component = PLANETS_RADIANCE_CPNT_DIFFUSE; 263 L_diffuse = planets_compute_radiance(cmd, &rad_args); /* [W/m^2/sr/m] */ 264 265 /* Calculate the weights [W/m^3] */ 266 weights[PLANETS_VOLRAD_DIRECT] = ka * (L_direct - S) / (wlen_pdf_m * dir_src_pdf); 267 weights[PLANETS_VOLRAD_DIFFUSE] = ka * (L_diffuse - S) / (wlen_pdf_m * dir_pdf); 268 } 269 270 /* Calculate the weights [W/m^3] */ 271 weights[PLANETS_VOLRAD_TOTAL] = 272 weights[PLANETS_VOLRAD_DIRECT] 273 + weights[PLANETS_VOLRAD_DIFFUSE]; 274 } 275 276 static void 277 solve_volumic_radiative_budget 278 (struct htrdr* htrdr, 279 const struct htrdr_solve_item_args* args, 280 void* data) 281 { 282 /* Volumic mesh on which volumic radiative budget is estimated */ 283 struct smsh_desc volrad_mesh_desc = SMSH_DESC_NULL; 284 285 /* Miscellaneous */ 286 struct htrdr_planets* cmd = NULL; 287 struct planets_voxel_radiative_budget* voxel = data; 288 size_t i = 0; 289 290 /* Preconditions */ 291 ASSERT(htrdr && args && data); 292 (void)htrdr, (void)cmd; 293 294 cmd = args->context; 295 ASSERT(cmd); 296 ASSERT(cmd->output_type == HTRDR_PLANETS_ARGS_OUTPUT_VOLUMIC_RADIATIVE_BUDGET); 297 298 SMSH(get_desc(cmd->volrad_mesh, &volrad_mesh_desc)); 299 300 /* Initialse voxel accumulators to 0 */ 301 *voxel = PLANETS_VOXEL_RADIATIVE_BUDGET_NULL; 302 303 FOR_EACH(i, 0, args->nrealisations) { 304 /* Time recording */ 305 struct time t0, t1; 306 307 /* Monte Carlo weights */ 308 double w[PLANETS_VOLRAD_WEIGHTS_COUNT] = {0}; /* [W/m^3] */ 309 double usec = 0; /* [us] */ 310 311 int iweight = 0; 312 313 /* Start of realisation time recording */ 314 time_current(&t0); 315 316 /* Run the realisation */ 317 realisation(cmd, args, &volrad_mesh_desc, w); 318 319 /* Stop time recording */ 320 time_sub(&t0, time_current(&t1), &t0); 321 usec = (double)time_val(&t0, TIME_NSEC) * 0.001; 322 323 FOR_EACH(iweight, 0, PLANETS_VOLRAD_WEIGHTS_COUNT){ 324 /* Update the volumic radiance budget accumulator */ 325 voxel->volrad_budget[iweight].sum_weights += w[iweight]; 326 voxel->volrad_budget[iweight].sum_weights_sqr += w[iweight]*w[iweight]; 327 voxel->volrad_budget[iweight].nweights += 1; 328 } 329 330 /* Update the per realisation time accumulator */ 331 voxel->time.sum_weights += usec; 332 voxel->time.sum_weights_sqr += usec*usec; 333 voxel->time.nweights += 1; 334 } 335 } 336 337 static res_T 338 write_buffer 339 (struct htrdr_planets* cmd, 340 struct htrdr_buffer* buf, 341 struct htrdr_accum* out_budget, /* W/m^3 */ 342 struct htrdr_accum* out_time, /* us */ 343 FILE* stream) 344 { 345 struct htrdr_buffer_layout layout = HTRDR_BUFFER_LAYOUT_NULL; 346 size_t x = 0; 347 348 /* Preconditions */ 349 ASSERT(cmd && buf && out_budget && out_time && stream); 350 351 htrdr_buffer_get_layout(buf, &layout); 352 ASSERT(layout.height == 1); 353 (void)cmd; 354 355 /* Reset global accumulators */ 356 *out_budget = HTRDR_ACCUM_NULL; 357 *out_time = HTRDR_ACCUM_NULL; 358 359 FOR_EACH(x, 0, layout.width) { 360 struct planets_voxel_radiative_budget* voxel = htrdr_buffer_at(buf, x, 0); 361 struct htrdr_estimate estim_time = HTRDR_ESTIMATE_NULL; 362 struct htrdr_accum* budget = NULL; 363 int iestim = 0; 364 365 budget = voxel->volrad_budget; 366 FOR_EACH(iestim, 0, PLANETS_VOLRAD_WEIGHTS_COUNT) { 367 struct htrdr_estimate estim_budget = HTRDR_ESTIMATE_NULL; 368 369 /* Write the estimate of the volumic radiative budget */ 370 htrdr_accum_get_estimation(&budget[iestim], &estim_budget); 371 fprintf(stream, "%g %g ", estim_budget.E, estim_budget.SE); 372 373 /* Write the accumulator of the volumic radiative budget */ 374 fprintf(stream, "%g %g ", 375 budget[iestim].sum_weights, 376 budget[iestim].sum_weights_sqr); 377 } 378 379 /* Write the number of realisations. 380 * It must be the same for all components */ 381 ASSERT(budget[PLANETS_VOLRAD_TOTAL].nweights 382 == budget[PLANETS_VOLRAD_DIRECT].nweights); 383 ASSERT(budget[PLANETS_VOLRAD_TOTAL].nweights 384 == budget[PLANETS_VOLRAD_DIFFUSE].nweights); 385 fprintf(stream, "%lu ", (unsigned long)budget[PLANETS_VOLRAD_TOTAL].nweights); 386 387 /* Write the estimate of the per realisation time */ 388 htrdr_accum_get_estimation(&voxel->time, &estim_time); 389 fprintf(stream, "%g %g\n", estim_time.E, estim_time.SE); 390 391 /* Update the overall (total) volumic radiative budget accumulator */ 392 out_budget->sum_weights += budget[PLANETS_VOLRAD_TOTAL].sum_weights; 393 out_budget->sum_weights_sqr += budget[PLANETS_VOLRAD_TOTAL].sum_weights_sqr; 394 out_budget->nweights += budget[PLANETS_VOLRAD_TOTAL].nweights; 395 396 /* Update the overall per realisation time accumulator */ 397 out_time->sum_weights += voxel->time.sum_weights; 398 out_time->sum_weights_sqr += voxel->time.sum_weights_sqr; 399 out_time->nweights += voxel->time.nweights; 400 } 401 return RES_OK; 402 } 403 404 /******************************************************************************* 405 * Local functions 406 ******************************************************************************/ 407 res_T 408 planets_solve_volrad_budget(struct htrdr_planets* cmd) 409 { 410 struct htrdr_solve_buffer_args args = HTRDR_SOLVE_BUFFER_ARGS_NULL; 411 struct htrdr_accum acc_budget = HTRDR_ACCUM_NULL; 412 struct htrdr_accum acc_time = HTRDR_ACCUM_NULL; 413 struct htrdr_estimate estim_budget = HTRDR_ESTIMATE_NULL; 414 struct htrdr_estimate estim_time = HTRDR_ESTIMATE_NULL; 415 res_T res = RES_OK; 416 417 /* Preconditions */ 418 ASSERT(cmd); 419 ASSERT(cmd->output_type == HTRDR_PLANETS_ARGS_OUTPUT_VOLUMIC_RADIATIVE_BUDGET); 420 421 args.solve_item = solve_volumic_radiative_budget; 422 args.buffer_layout = cmd->buf_layout; 423 args.nrealisations = cmd->spt; 424 args.context = cmd; 425 426 res = htrdr_solve_buffer(cmd->htrdr, &args, cmd->buf); 427 if(res != RES_OK) goto error; 428 429 /* Nothing more to do on non master processes */ 430 if(htrdr_get_mpi_rank(cmd->htrdr) != 0) goto exit; 431 432 /* Write outut data */ 433 res = write_buffer(cmd, cmd->buf, &acc_budget, &acc_time, cmd->output); 434 if(res != RES_OK) goto error; 435 436 htrdr_accum_get_estimation(&acc_time, &estim_time); 437 htrdr_accum_get_estimation(&acc_budget, &estim_budget); 438 439 /* Write overall results */ 440 htrdr_log(cmd->htrdr, "Time per radiative path (in µs): %g +/- %g\n", 441 estim_time.E, estim_time.SE); 442 htrdr_log(cmd->htrdr, "Volumic radiative budget (in W/m³): %g +/- %g\n", 443 estim_budget.E, estim_budget.SE); 444 445 exit: 446 return res; 447 error: 448 goto exit; 449 }