sdis_solve_medium_Xd.h (26757B)
1 /* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) 2 * 3 * This program is free software: you can redistribute it and/or modify 4 * it under the terms of the GNU General Public License as published by 5 * the Free Software Foundation, either version 3 of the License, or 6 * (at your option) any later version. 7 * 8 * This program is distributed in the hope that it will be useful, 9 * but WITHOUT ANY WARRANTY; without even the implied warranty of 10 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 11 * GNU General Public License for more details. 12 * 13 * You should have received a copy of the GNU General Public License 14 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 15 16 #include "sdis_c.h" 17 #include "sdis_device_c.h" 18 #include "sdis_estimator_c.h" 19 #include "sdis_interface_c.h" 20 #include "sdis_log.h" 21 #include "sdis_green.h" 22 #include "sdis_realisation.h" 23 #include "sdis_scene_c.h" 24 25 #include <rsys/algorithm.h> 26 #include <rsys/clock_time.h> 27 #include <rsys/dynamic_array.h> 28 29 #include <omp.h> 30 31 #include "sdis_Xd_begin.h" 32 33 #ifndef SDIS_SOLVE_MEDIUM_XD_H 34 #define SDIS_SOLVE_MEDIUM_XD_H 35 36 /* 37 * Define the data structures and functions that are not generic to the 38 * SDIS_XD_DIMENSION parameter 39 */ 40 41 struct enclosure_cumul { 42 unsigned enc_id; 43 double cumul; 44 }; 45 46 /* Define the darray_enclosure_cumul dynamic array */ 47 #define DARRAY_NAME enclosure_cumul 48 #define DARRAY_DATA struct enclosure_cumul 49 #include <rsys/dynamic_array.h> 50 51 /******************************************************************************* 52 * Helper functions 53 ******************************************************************************/ 54 static INLINE int 55 cmp_double_to_enc_cumuls(const void* a, const void* b) 56 { 57 const double key = *(const double*)a; 58 const struct enclosure_cumul* enc_cumul = (const struct enclosure_cumul*)b; 59 if(key < enc_cumul->cumul) return -1; 60 if(key > enc_cumul->cumul) return +1; 61 return 0; 62 } 63 64 static res_T 65 compute_medium_enclosure_cumulative 66 (struct sdis_scene* scn, 67 const struct sdis_medium* mdm, 68 struct darray_enclosure_cumul* cumul) 69 { 70 struct htable_enclosure_iterator it, end; 71 double accum = 0; 72 res_T res = RES_OK; 73 ASSERT(scn && mdm && cumul); 74 75 darray_enclosure_cumul_clear(cumul); 76 77 htable_enclosure_begin(&scn->enclosures, &it); 78 htable_enclosure_end(&scn->enclosures, &end); 79 while(!htable_enclosure_iterator_eq(&it, &end)) { 80 struct enclosure_cumul enc_cumul; 81 const struct enclosure* enc = htable_enclosure_iterator_data_get(&it); 82 const unsigned* enc_id = htable_enclosure_iterator_key_get(&it); 83 htable_enclosure_iterator_next(&it); 84 85 if(sdis_medium_get_id(mdm) != enc->medium_id) continue; 86 87 accum += enc->V; 88 enc_cumul.enc_id = *enc_id; 89 enc_cumul.cumul = accum; 90 res = darray_enclosure_cumul_push_back(cumul, &enc_cumul); 91 if(res != RES_OK) goto error; 92 } 93 94 if(darray_enclosure_cumul_size_get(cumul) == 0) { 95 log_err(scn->dev, 96 "%s: there is no enclosure that encompasses the submitted medium.\n", 97 FUNC_NAME); 98 res = RES_BAD_ARG; 99 goto error; 100 } 101 102 exit: 103 return res; 104 error: 105 darray_enclosure_cumul_clear(cumul); 106 goto exit; 107 } 108 109 static unsigned 110 sample_medium_enclosure 111 (const struct darray_enclosure_cumul* cumul, struct ssp_rng* rng) 112 { 113 const struct enclosure_cumul* enc_cumuls = NULL; 114 const struct enclosure_cumul* enc_cumul_found = NULL; 115 double r; 116 size_t i; 117 size_t sz; 118 ASSERT(cumul && rng && darray_enclosure_cumul_size_get(cumul)); 119 120 sz = darray_enclosure_cumul_size_get(cumul); 121 enc_cumuls = darray_enclosure_cumul_cdata_get(cumul); 122 if(sz == 1) { 123 enc_cumul_found = enc_cumuls; 124 } else { 125 /* Generate an uniform random number in [0, cumul[ */ 126 r = ssp_rng_canonical(rng); 127 r = r * enc_cumuls[sz-1].cumul; 128 129 enc_cumul_found = search_lower_bound 130 (&r, enc_cumuls, sz, sizeof(*enc_cumuls), cmp_double_to_enc_cumuls); 131 ASSERT(enc_cumul_found); 132 133 /* search_lower_bound returns the first entry that is not less than `r'. 134 * The following code discards entries that are also equal to `r'. */ 135 i = (size_t)(enc_cumul_found - enc_cumuls); 136 while(enc_cumuls[i].cumul == r && i < sz) ++i; 137 ASSERT(i < sz); 138 139 enc_cumul_found = enc_cumuls + i; 140 } 141 return enc_cumul_found->enc_id; 142 } 143 144 static INLINE res_T 145 check_solve_medium_args(const struct sdis_solve_medium_args* args) 146 { 147 if(!args) return RES_BAD_ARG; 148 149 /* Check the medium */ 150 if(!args->medium) return RES_BAD_ARG; 151 152 /* Check #realisations */ 153 if(!args->nrealisations || args->nrealisations > INT64_MAX) { 154 return RES_BAD_ARG; 155 } 156 157 /* Check time range */ 158 if(args->time_range[0] < 0 || args->time_range[1] < args->time_range[0]) { 159 return RES_BAD_ARG; 160 } 161 if(args->time_range[1] > DBL_MAX 162 && args->time_range[0] != args->time_range[1]) { 163 return RES_BAD_ARG; 164 } 165 166 /* Check picard order */ 167 if(args->picard_order < 1) { 168 return RES_BAD_ARG; 169 } 170 171 /* Check the RNG type */ 172 if(!args->rng_state && args->rng_type >= SSP_RNG_TYPES_COUNT__) { 173 return RES_BAD_ARG; 174 } 175 176 /* Check the diffusion algorithm */ 177 if((unsigned)args->diff_algo >= SDIS_DIFFUSION_ALGORITHMS_COUNT__) { 178 return RES_BAD_ARG; 179 } 180 181 return RES_OK; 182 } 183 184 static INLINE res_T 185 check_compute_power_args(const struct sdis_compute_power_args* args) 186 { 187 if(!args) return RES_BAD_ARG; 188 189 /* Check the medium */ 190 if(!args->medium) return RES_BAD_ARG; 191 192 /* Check #realisations */ 193 if(!args->nrealisations || args->nrealisations > INT64_MAX) { 194 return RES_BAD_ARG; 195 } 196 197 /* Check time range */ 198 if(args->time_range[0] < 0 || args->time_range[1] < args->time_range[0]) { 199 return RES_BAD_ARG; 200 } 201 if(args->time_range[1] > DBL_MAX 202 && args->time_range[0] != args->time_range[1]) { 203 return RES_BAD_ARG; 204 } 205 206 /* Check the RNG type */ 207 if(!args->rng_state && args->rng_type >= SSP_RNG_TYPES_COUNT__) { 208 return RES_BAD_ARG; 209 } 210 211 return RES_OK; 212 } 213 214 #endif /* !SDIS_SOLVE_MEDIUM_XD_H */ 215 216 /******************************************************************************* 217 * Helper functions generic to the SDIS_XD_DIMENSION parameter 218 ******************************************************************************/ 219 static res_T 220 XD(sample_enclosure_position) 221 (struct sdis_scene* scn, 222 const unsigned enc_id, 223 struct ssp_rng* rng, 224 double pos[DIM]) 225 { 226 const size_t MAX_NCHALLENGES = 1000; 227 228 const struct enclosure* enc = NULL; 229 float lower[DIM], upper[DIM]; 230 size_t ichallenge; 231 size_t i; 232 res_T res = RES_OK; 233 ASSERT(scn && rng && pos); 234 ASSERT(enc_id != ENCLOSURE_ID_NULL); 235 236 enc = scene_get_enclosure(scn, enc_id); 237 SXD(scene_view_get_aabb(enc->sXd(view), lower, upper)); 238 239 FOR_EACH(i, 0, DIM) { 240 if(lower[i] > upper[i] || eq_epsf(lower[i], upper[i], 1.e-6f)) { 241 res = RES_BAD_ARG; /* Degenerated enclosure */ 242 goto error; 243 } 244 } 245 246 FOR_EACH(ichallenge, 0, MAX_NCHALLENGES) { 247 struct sXd(hit) hit = SXD_HIT_NULL; 248 const float dir[3] = {1,0,0}; 249 const float range[2] = {FLT_MIN, FLT_MAX}; 250 float org[DIM]; 251 252 /* Generate an uniform position into the enclosure AABB */ 253 FOR_EACH(i, 0, DIM) { 254 org[i] = ssp_rng_uniform_float(rng, lower[i], upper[i]); 255 } 256 257 /* Check that pos lies into the enclosure; trace a ray and check that it 258 * hits something and that the normal points towards the traced ray 259 * direction (enclosure normals point inword the enclosure) */ 260 SXD(scene_view_trace_ray(enc->sXd(view), org, dir, range, NULL, &hit)); 261 if(!SXD_HIT_NONE(&hit) && fX(dot)(dir, hit.normal) < 0) { 262 dX_set_fX(pos, org); 263 break; 264 } 265 } 266 267 if(ichallenge >= MAX_NCHALLENGES) { 268 res = RES_BAD_ARG; 269 goto error; 270 } 271 272 exit: 273 return res; 274 error: 275 goto exit; 276 } 277 278 /******************************************************************************* 279 * Local function 280 ******************************************************************************/ 281 static res_T 282 XD(solve_medium) 283 (struct sdis_scene* scn, 284 const struct sdis_solve_medium_args* args, 285 struct sdis_green_function** out_green, /* May be NULL <=> No green func */ 286 struct sdis_estimator** out_estimator) /* May be NULL <=> No estimator */ 287 { 288 /* Time registration */ 289 struct time time0, time1; 290 char buf[128]; /* Temporary buffer used to store formated time */ 291 292 /* Device variables */ 293 struct mem_allocator* allocator = NULL; 294 size_t nthreads = 0; 295 296 /* Stardis variables */ 297 struct sdis_estimator* estimator = NULL; 298 struct sdis_green_function* green = NULL; 299 struct sdis_green_function** per_thread_green = NULL; 300 301 /* Random number generator */ 302 struct ssp_rng_proxy* rng_proxy = NULL; 303 struct ssp_rng** per_thread_rng = NULL; 304 305 /* Miscellaneous */ 306 struct darray_enclosure_cumul cumul; 307 struct accum* per_thread_acc_temp = NULL; 308 struct accum* per_thread_acc_time = NULL; 309 size_t nrealisations = 0; 310 int64_t irealisation; 311 int32_t* progress = NULL; /* Per process progress bar */ 312 int pcent_progress = 1; /* Percentage requiring progress update */ 313 int is_master_process = 1; 314 int cumul_is_init = 0; 315 int register_paths = SDIS_HEAT_PATH_NONE; 316 ATOMIC nsolved_realisations = 0; 317 ATOMIC res = RES_OK; 318 319 if(!scn) { res = RES_BAD_ARG; goto error; } 320 if(!out_estimator && !out_green) { res = RES_BAD_ARG; goto error; } 321 res = check_solve_medium_args(args); 322 if(res != RES_OK) goto error; 323 res = XD(scene_check_dimensionality)(scn); 324 if(res != RES_OK) goto error; 325 326 if(out_green && args->picard_order != 1) { 327 log_err(scn->dev, "%s: the evaluation of the green function does not make " 328 "sense when dealing with the non-linearities of the system; i.e. picard " 329 "order must be set to 1 while it is currently set to %lu.\n", 330 FUNC_NAME, (unsigned long)args->picard_order); 331 res = RES_BAD_ARG; 332 goto error; 333 } 334 335 #ifdef SDIS_ENABLE_MPI 336 is_master_process = !scn->dev->use_mpi || scn->dev->mpi_rank == 0; 337 #endif 338 339 nthreads = scn->dev->nthreads; 340 allocator = scn->dev->allocator; 341 342 /* Update the progress bar every percent if escape sequences are allowed in 343 * log messages or only every 10 percent when only plain text is allowed. 344 * This reduces the number of lines of plain text printed */ 345 pcent_progress = scn->dev->no_escape_sequence ? 10 : 1; 346 347 /* Create the per thread RNGs */ 348 res = create_per_thread_rng 349 (scn->dev, args->rng_state, args->rng_type, &rng_proxy, &per_thread_rng); 350 if(res != RES_OK) goto error; 351 352 /* Allocate the per process progress status */ 353 res = alloc_process_progress(scn->dev, &progress); 354 if(res != RES_OK) goto error; 355 356 /* Create the per thread accumulators */ 357 per_thread_acc_temp = MEM_CALLOC(allocator, nthreads, sizeof(struct accum)); 358 per_thread_acc_time = MEM_CALLOC(allocator, nthreads, sizeof(struct accum)); 359 if(!per_thread_acc_temp) { res = RES_MEM_ERR; goto error; } 360 if(!per_thread_acc_time) { res = RES_MEM_ERR; goto error; } 361 362 /* Compute the enclosure cumulative */ 363 darray_enclosure_cumul_init(scn->dev->allocator, &cumul); 364 cumul_is_init = 1; 365 res = compute_medium_enclosure_cumulative(scn, args->medium, &cumul); 366 if(res != RES_OK) goto error; 367 368 if(out_green) { 369 res = create_per_thread_green_function 370 (scn, args->signature, &per_thread_green); 371 if(res != RES_OK) goto error; 372 } 373 374 /* Create the estimator on the master process only. No estimator is needed 375 * for non master process */ 376 if(out_estimator && is_master_process) { 377 res = estimator_create(scn->dev, SDIS_ESTIMATOR_TEMPERATURE, &estimator); 378 if(res != RES_OK) goto error; 379 } 380 381 /* Synchronise the processes */ 382 process_barrier(scn->dev); 383 384 #define PROGRESS_MSG "Solving medium temperature: " 385 print_progress(scn->dev, progress, PROGRESS_MSG); 386 387 /* Begin time registration of the computation */ 388 time_current(&time0); 389 390 /* Here we go! Launch the Monte Carlo estimation */ 391 nrealisations = compute_process_realisations_count(scn->dev, args->nrealisations); 392 register_paths = out_estimator && is_master_process 393 ? args->register_paths : SDIS_HEAT_PATH_NONE; 394 omp_set_num_threads((int)scn->dev->nthreads); 395 #pragma omp parallel for schedule(static) 396 for(irealisation = 0; irealisation < (int64_t)nrealisations; ++irealisation) { 397 struct probe_realisation_args realis_args = PROBE_REALISATION_ARGS_NULL; 398 struct time t0, t1; 399 const int ithread = omp_get_thread_num(); 400 struct ssp_rng* rng = per_thread_rng[ithread]; 401 struct accum* acc_temp = &per_thread_acc_temp[ithread]; 402 struct accum* acc_time = &per_thread_acc_time[ithread]; 403 struct green_path_handle* pgreen_path = NULL; 404 struct green_path_handle green_path = GREEN_PATH_HANDLE_NULL; 405 struct sdis_heat_path* pheat_path = NULL; 406 struct sdis_heat_path heat_path; 407 double weight; 408 double time; 409 double pos[DIM]; 410 size_t n; 411 unsigned enc_id = ENCLOSURE_ID_NULL; 412 int pcent; 413 res_T res_local = RES_OK; 414 res_T res_simul = RES_OK; 415 416 if(ATOMIC_GET(&res) != RES_OK) continue; /* An error occurred */ 417 418 time_current(&t0); 419 420 time = sample_time(rng, args->time_range); 421 if(out_green) { 422 res_local = green_function_create_path(per_thread_green[ithread], &green_path); 423 if(res_local != RES_OK) { ATOMIC_SET(&res, res_local); goto error_it; } 424 pgreen_path = &green_path; 425 } 426 427 if(register_paths) { 428 heat_path_init(scn->dev->allocator, &heat_path); 429 pheat_path = &heat_path; 430 } 431 432 /* Uniformly Sample an enclosure that surround the submitted medium and 433 * uniformly sample a position into it */ 434 enc_id = sample_medium_enclosure(&cumul, rng); 435 res_local = XD(sample_enclosure_position)(scn, enc_id, rng, pos); 436 if(res_local != RES_OK) { 437 log_err(scn->dev, "%s: could not sample a medium position.\n", FUNC_NAME); 438 ATOMIC_SET(&res, res_local); 439 goto error_it; 440 } 441 442 /* Run a probe realisation */ 443 realis_args.rng = rng; 444 realis_args.enc_id = enc_id; 445 realis_args.time = time; 446 realis_args.picard_order = args->picard_order; 447 realis_args.green_path = pgreen_path; 448 realis_args.heat_path = pheat_path; 449 realis_args.irealisation = (size_t)irealisation; 450 realis_args.diff_algo = args->diff_algo; 451 dX(set)(realis_args.position, pos); 452 res_simul = XD(probe_realisation)(scn, &realis_args, &weight); 453 if(res_simul != RES_OK && res_simul != RES_BAD_OP) { 454 ATOMIC_SET(&res, res_simul); 455 goto error_it; 456 } 457 458 /* Finalize the registered path */ 459 if(pheat_path) { 460 pheat_path->status = res_simul == RES_OK 461 ? SDIS_HEAT_PATH_SUCCESS 462 : SDIS_HEAT_PATH_FAILURE; 463 464 /* Check if the path must be saved regarding the register_paths mask */ 465 if(!(register_paths & (int)pheat_path->status)) { 466 heat_path_release(pheat_path); 467 pheat_path = NULL; 468 } else { /* Register the sampled path */ 469 res_local = estimator_add_and_release_heat_path(estimator, pheat_path); 470 if(res_local != RES_OK) { 471 ATOMIC_SET(&res, res_local); 472 goto error_it; 473 } 474 } 475 pheat_path = NULL; 476 } 477 478 /* Stop time registration */ 479 time_sub(&t0, time_current(&t1), &t0); 480 481 /* Update accumulators */ 482 if(res_simul == RES_OK) { 483 const double usec = (double)time_val(&t0, TIME_NSEC) * 0.001; 484 acc_temp->sum += weight; acc_temp->sum2 += weight*weight; ++acc_temp->count; 485 acc_time->sum += usec; acc_time->sum2 += usec*usec; ++acc_time->count; 486 } 487 488 /* Update progress */ 489 n = (size_t)ATOMIC_INCR(&nsolved_realisations); 490 pcent = (int)((double)n * 100.0 / (double)nrealisations + 0.5/*round*/); 491 #pragma omp critical 492 if(pcent/pcent_progress > progress[0]/pcent_progress) { 493 progress[0] = pcent; 494 print_progress_update(scn->dev, progress, PROGRESS_MSG); 495 } 496 exit_it: 497 if(pheat_path) heat_path_release(pheat_path); 498 continue; 499 error_it: 500 goto exit_it; 501 } 502 /* Synchronise processes */ 503 process_barrier(scn->dev); 504 505 res = gather_res_T(scn->dev, (res_T)res); 506 if(res != RES_OK) goto error; 507 508 print_progress_completion(scn->dev, progress, PROGRESS_MSG); 509 #undef PROGRESS_MSG 510 511 /* Report computation time */ 512 time_sub(&time0, time_current(&time1), &time0); 513 time_dump(&time0, TIME_ALL, NULL, buf, sizeof(buf)); 514 log_info(scn->dev, "Medium temperature solved in %s.\n", buf); 515 516 /* Gather the RNG proxy sequence IDs and ensure that the RNG proxy state of 517 * the master process is greater than the RNG proxy state of all other 518 * processes */ 519 res = gather_rng_proxy_sequence_id(scn->dev, rng_proxy); 520 if(res != RES_OK) goto error; 521 522 /* Setup the estimated temperature */ 523 if(out_estimator) { 524 struct accum acc_temp; 525 struct accum acc_time; 526 527 time_current(&time0); 528 529 #define GATHER_ACCUMS(Msg, Acc) { \ 530 res = gather_accumulators(scn->dev, Msg, per_thread_##Acc, &Acc); \ 531 if(res != RES_OK) goto error; \ 532 } (void)0 533 GATHER_ACCUMS(MPI_SDIS_MSG_ACCUM_TEMP, acc_temp); 534 GATHER_ACCUMS(MPI_SDIS_MSG_ACCUM_TIME, acc_time); 535 #undef GATHER_ACCUMS 536 537 time_sub(&time0, time_current(&time1), &time0); 538 time_dump(&time0, TIME_ALL, NULL, buf, sizeof(buf)); 539 log_info(scn->dev, "Accumulators gathered in %s.\n", buf); 540 541 /* Return an estimator only on master process */ 542 if(is_master_process) { 543 ASSERT(acc_temp.count == acc_time.count); 544 estimator_setup_realisations_count(estimator, args->nrealisations, acc_temp.count); 545 estimator_setup_temperature(estimator, acc_temp.sum, acc_temp.sum2); 546 estimator_setup_realisation_time(estimator, acc_time.sum, acc_time.sum2); 547 res = estimator_save_rng_state(estimator, rng_proxy); 548 if(res != RES_OK) goto error; 549 } 550 } 551 552 if(out_green) { 553 time_current(&time0); 554 555 res = gather_green_functions 556 (scn, rng_proxy, per_thread_green, per_thread_acc_time, &green); 557 if(res != RES_OK) goto error; 558 559 time_sub(&time0, time_current(&time1), &time0); 560 time_dump(&time0, TIME_ALL, NULL, buf, sizeof(buf)); 561 log_info(scn->dev, "Green functions gathered in %s.\n", buf); 562 563 /* Return a green function only on master process */ 564 if(!is_master_process) { 565 SDIS(green_function_ref_put(green)); 566 green = NULL; 567 } 568 } 569 570 exit: 571 if(per_thread_rng) release_per_thread_rng(scn->dev, per_thread_rng); 572 if(per_thread_green) release_per_thread_green_function(scn, per_thread_green); 573 if(progress) free_process_progress(scn->dev, progress); 574 if(per_thread_acc_temp) MEM_RM(scn->dev->allocator, per_thread_acc_temp); 575 if(per_thread_acc_time) MEM_RM(scn->dev->allocator, per_thread_acc_time); 576 if(cumul_is_init) darray_enclosure_cumul_release(&cumul); 577 if(rng_proxy) SSP(rng_proxy_ref_put(rng_proxy)); 578 if(out_estimator) *out_estimator = estimator; 579 if(out_green) *out_green = green; 580 return (res_T)res; 581 error: 582 if(green) { SDIS(green_function_ref_put(green)); green = NULL; } 583 if(estimator) { SDIS(estimator_ref_put(estimator)); estimator = NULL; } 584 goto exit; 585 } 586 587 static res_T 588 XD(compute_power) 589 (struct sdis_scene* scn, 590 const struct sdis_compute_power_args* args, 591 struct sdis_estimator** out_estimator) 592 { 593 /* Time registration */ 594 struct time time0, time1; 595 char buf[128]; /* Temporary buffer used to store formated time */ 596 597 /* Device variables */ 598 struct mem_allocator* allocator = NULL; 599 size_t nthreads = 0; 600 601 /* Stardis variables */ 602 struct sdis_estimator* estimator = NULL; 603 604 /* Random number generator */ 605 struct ssp_rng_proxy* rng_proxy = NULL; 606 struct ssp_rng** per_thread_rng = NULL; 607 608 /* Miscellaneous */ 609 struct darray_enclosure_cumul cumul; 610 struct accum* per_thread_acc_mpow = NULL; 611 struct accum* per_thread_acc_time = NULL; 612 double spread = 0; 613 size_t nrealisations = 0; 614 int64_t irealisation = 0; 615 int32_t* progress = NULL; /* Per process progress bar */ 616 int cumul_is_init = 0; 617 int is_master_process = 1; 618 ATOMIC nsolved_realisations = 0; 619 ATOMIC res = RES_OK; 620 621 if(!scn) { res = RES_BAD_ARG; goto error; } 622 if(!out_estimator) { res = RES_BAD_ARG; goto error; } 623 res = check_compute_power_args(args); 624 if(res != RES_OK) goto error; 625 res = XD(scene_check_dimensionality)(scn); 626 if(res != RES_OK) goto error; 627 628 if(sdis_medium_get_type(args->medium) != SDIS_SOLID) { 629 log_err(scn->dev, "Could not compute mean power on a non solid medium.\n"); 630 res = RES_BAD_ARG; 631 goto error; 632 } 633 634 #ifdef SDIS_ENABLE_MPI 635 is_master_process = !scn->dev->use_mpi || scn->dev->mpi_rank == 0; 636 #endif 637 638 nthreads = scn->dev->nthreads; 639 allocator = scn->dev->allocator; 640 641 /* Create the per thread RNGs */ 642 res = create_per_thread_rng 643 (scn->dev, args->rng_state, args->rng_type, &rng_proxy, &per_thread_rng); 644 if(res != RES_OK) goto error; 645 646 /* Allocate the per process progress status */ 647 res = alloc_process_progress(scn->dev, &progress); 648 if(res != RES_OK) goto error; 649 650 /* Create the per thread accumulators */ 651 per_thread_acc_mpow = MEM_CALLOC(allocator, nthreads, sizeof(struct accum)); 652 per_thread_acc_time = MEM_CALLOC(allocator, nthreads, sizeof(struct accum)); 653 if(!per_thread_acc_mpow) { res = RES_MEM_ERR; goto error; } 654 if(!per_thread_acc_time) { res = RES_MEM_ERR; goto error; } 655 656 /* Compute the cumulative of the spreads of the enclosures surrounding the 657 * submitted medium */ 658 darray_enclosure_cumul_init(scn->dev->allocator, &cumul); 659 cumul_is_init = 1; 660 res = compute_medium_enclosure_cumulative(scn, args->medium, &cumul); 661 if(res != RES_OK) goto error; 662 663 /* Fetch the overall medium spread from the unormalized enclosure cumulative */ 664 spread = darray_enclosure_cumul_cdata_get(&cumul) 665 [darray_enclosure_cumul_size_get(&cumul)-1].cumul; 666 667 /* Create the estimator on the master process only. No estimator is needed 668 * for non master process */ 669 if(is_master_process) { 670 res = estimator_create(scn->dev, SDIS_ESTIMATOR_POWER, &estimator); 671 if(res != RES_OK) goto error; 672 } 673 674 /* Synchronise the processes */ 675 process_barrier(scn->dev); 676 677 #define PROGRESS_MSG "Computing mean power: " 678 print_progress(scn->dev, progress, PROGRESS_MSG); 679 680 /* Begin time registration of the computation */ 681 time_current(&time0); 682 683 /* Here we go! Launch the Monte Carlo estimation */ 684 nrealisations = compute_process_realisations_count(scn->dev, args->nrealisations); 685 omp_set_num_threads((int)scn->dev->nthreads); 686 #pragma omp parallel for schedule(static) 687 for(irealisation = 0; irealisation < (int64_t)nrealisations; ++irealisation) { 688 struct time t0, t1; 689 struct sdis_rwalk_vertex vtx = SDIS_RWALK_VERTEX_NULL; 690 const int ithread = omp_get_thread_num(); 691 struct ssp_rng* rng = per_thread_rng[ithread]; 692 struct accum* acc_mpow = &per_thread_acc_mpow[ithread]; 693 struct accum* acc_time = &per_thread_acc_time[ithread]; 694 double power = 0; 695 double usec = 0; 696 size_t n = 0; 697 unsigned enc_id = ENCLOSURE_ID_NULL; 698 int pcent = 0; 699 res_T res_local = RES_OK; 700 701 if(ATOMIC_GET(&res) != RES_OK) continue; /* An error occurred */ 702 703 /* Begin time registration */ 704 time_current(&t0); 705 706 /* Sample the time */ 707 vtx.time = sample_time(rng, args->time_range); 708 709 /* Uniformly Sample an enclosure that surround the submitted medium and 710 * uniformly sample a position into it */ 711 enc_id = sample_medium_enclosure(&cumul, rng); 712 res_local = XD(sample_enclosure_position)(scn, enc_id, rng, vtx.P); 713 if(res_local != RES_OK) { 714 log_err(scn->dev, "%s: could not sample a medium position.\n", FUNC_NAME); 715 ATOMIC_SET(&res, res_local); 716 goto error_it; 717 } 718 719 /* Fetch the volumic power */ 720 power = solid_get_volumic_power(args->medium, &vtx); 721 722 /* Stop time registration */ 723 time_sub(&t0, time_current(&t1), &t0); 724 usec = (double)time_val(&t0, TIME_NSEC) * 0.001; 725 726 /* Update accumulators */ 727 acc_mpow->sum += power; acc_mpow->sum2 += power*power; ++acc_mpow->count; 728 acc_time->sum += usec; acc_time->sum2 += usec*usec; ++acc_time->count; 729 730 /* Update progress */ 731 n = (size_t)ATOMIC_INCR(&nsolved_realisations); 732 pcent = (int)((double)n * 100.0 / (double)nrealisations + 0.5/*round*/); 733 #pragma omp critical 734 if(pcent > progress[0]) { 735 progress[0] = pcent; 736 print_progress_update(scn->dev, progress, PROGRESS_MSG); 737 } 738 exit_it: 739 continue; 740 error_it: 741 goto exit_it; 742 } 743 /* Synchronise the processes */ 744 process_barrier(scn->dev); 745 746 res = gather_res_T(scn->dev, (res_T)res); 747 if(res != RES_OK) goto error; 748 749 print_progress_update(scn->dev, progress, PROGRESS_MSG); 750 log_info(scn->dev, "\n"); 751 #undef PROGRESS_MSG 752 753 /* Report computation time */ 754 time_sub(&time0, time_current(&time1), &time0); 755 time_dump(&time0, TIME_ALL, NULL, buf, sizeof(buf)); 756 log_info(scn->dev, "Mean power computed in in %s.\n", buf); 757 758 /* Gather the RNG proxy sequence IDs and ensure that the RNG proxy state of 759 * the master process is greater than the RNG proxy state of all other 760 * processes */ 761 res = gather_rng_proxy_sequence_id(scn->dev, rng_proxy); 762 if(res != RES_OK) goto error; 763 764 /* Setup the estimated mean power */ 765 { 766 struct accum acc_mpow; 767 struct accum acc_time; 768 769 time_current(&time0); 770 771 #define GATHER_ACCUMS(Msg, Acc) { \ 772 res = gather_accumulators(scn->dev, Msg, per_thread_##Acc, &Acc); \ 773 if(res != RES_OK) goto error; \ 774 } (void)0 775 GATHER_ACCUMS(MPI_SDIS_MSG_ACCUM_MEAN_POWER, acc_mpow); 776 GATHER_ACCUMS(MPI_SDIS_MSG_ACCUM_TIME, acc_time); 777 #undef GATHER_ACCUMS 778 779 time_sub(&time0, time_current(&time1), &time0); 780 time_dump(&time0, TIME_ALL, NULL, buf, sizeof(buf)); 781 log_info(scn->dev, "Accumulators gathered in %s.\n", buf); 782 783 /* Return an estimator only on master process */ 784 if(is_master_process) { 785 ASSERT(acc_mpow.count == acc_time.count); 786 estimator_setup_realisations_count(estimator, args->nrealisations, acc_mpow.count); 787 estimator_setup_realisation_time(estimator, acc_time.sum, acc_time.sum2); 788 estimator_setup_power 789 (estimator, acc_mpow.sum, acc_mpow.sum2, spread, args->time_range); 790 res = estimator_save_rng_state(estimator, rng_proxy); 791 if(res != RES_OK) goto error; 792 } 793 } 794 795 exit: 796 if(per_thread_rng) release_per_thread_rng(scn->dev, per_thread_rng); 797 if(progress) free_process_progress(scn->dev, progress); 798 if(per_thread_acc_mpow) MEM_RM(scn->dev->allocator, per_thread_acc_mpow); 799 if(per_thread_acc_time) MEM_RM(scn->dev->allocator, per_thread_acc_time); 800 if(cumul_is_init) darray_enclosure_cumul_release(&cumul); 801 if(rng_proxy) SSP(rng_proxy_ref_put(rng_proxy)); 802 if(out_estimator) *out_estimator = estimator; 803 return (res_T)res; 804 error: 805 if(estimator) { SDIS(estimator_ref_put(estimator)); estimator = NULL; } 806 goto exit; 807 } 808 809 #include "sdis_Xd_end.h"