sdis_solve_boundary_Xd.h (33204B)
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_medium_c.h" 23 #include "sdis_misc.h" 24 #include "sdis_realisation.h" 25 #include "sdis_scene_c.h" 26 #include "sdis_heat_path_boundary_c.h" /* check_Tref_<2d|3d> */ 27 28 #include <rsys/clock_time.h> 29 #include <star/ssp.h> 30 #include <omp.h> 31 32 #include "sdis_Xd_begin.h" 33 34 struct XD(boundary_context) { 35 struct sXd(scene_view)* view; 36 const size_t* primitives; 37 }; 38 static const struct XD(boundary_context) XD(BOUNDARY_CONTEXT_NULL) = { 39 NULL, NULL 40 }; 41 42 /******************************************************************************* 43 * Non generic helper functions 44 ******************************************************************************/ 45 #ifndef SDIS_SOLVE_BOUNDARY_XD 46 #define SDIS_SOLVE_BOUNDARY_XD 47 48 static INLINE res_T 49 check_solve_boundary_args(const struct sdis_solve_boundary_args* args) 50 { 51 if(!args) return RES_BAD_ARG; 52 53 /* Check #realisations */ 54 if(!args->nrealisations || args->nrealisations > INT64_MAX) { 55 return RES_BAD_ARG; 56 } 57 58 /* Check the list of primitives */ 59 if(!args->primitives || !args->sides || !args->nprimitives) { 60 return RES_BAD_ARG; 61 } 62 63 /* Check time range */ 64 if(args->time_range[0] < 0 || args->time_range[1] < args->time_range[0]) { 65 return RES_BAD_ARG; 66 } 67 if(args->time_range[1] > DBL_MAX 68 && args->time_range[0] != args->time_range[1]) { 69 return RES_BAD_ARG; 70 } 71 72 /* Check picard order */ 73 if(args->picard_order < 1) { 74 return RES_BAD_ARG; 75 } 76 77 /* Check RNG type */ 78 if(!args->rng_state && args->rng_type >= SSP_RNG_TYPES_COUNT__) { 79 return RES_BAD_ARG; 80 } 81 82 /* Check the diffusion algorithm */ 83 if((unsigned)args->diff_algo >= SDIS_DIFFUSION_ALGORITHMS_COUNT__) { 84 return RES_BAD_ARG; 85 } 86 87 return RES_OK; 88 } 89 90 static INLINE res_T 91 check_solve_boundary_flux_args(const struct sdis_solve_boundary_flux_args* args) 92 { 93 if(!args) return RES_BAD_ARG; 94 95 /* Check #realisations */ 96 if(!args->nrealisations || args->nrealisations > INT64_MAX) { 97 return RES_BAD_ARG; 98 } 99 100 /* Check the list of primitives */ 101 if(!args->primitives || !args->nprimitives) { 102 return RES_BAD_ARG; 103 } 104 105 /* Check time range */ 106 if(args->time_range[0] < 0 || args->time_range[1] < args->time_range[0]) { 107 return RES_BAD_ARG; 108 } 109 if(args->time_range[1] > DBL_MAX 110 && args->time_range[0] != args->time_range[1]) { 111 return RES_BAD_ARG; 112 } 113 114 /* Check picard order */ 115 if(args->picard_order < 1) { 116 return RES_BAD_ARG; 117 } 118 119 /* Check RNG type */ 120 if(!args->rng_state && args->rng_type >= SSP_RNG_TYPES_COUNT__) { 121 return RES_BAD_ARG; 122 } 123 124 return RES_OK; 125 } 126 127 #endif /* SDIS_SOLVE_BOUNDARY_XD */ 128 129 /******************************************************************************* 130 * Helper functions 131 ******************************************************************************/ 132 static INLINE void 133 XD(boundary_get_indices)(const unsigned iprim, unsigned ids[DIM], void* context) 134 { 135 unsigned i; 136 (void)context; 137 ASSERT(ids); 138 FOR_EACH(i, 0, DIM) ids[i] = iprim*DIM + i; 139 } 140 141 static INLINE void 142 XD(boundary_get_position)(const unsigned ivert, float pos[DIM], void* context) 143 { 144 struct XD(boundary_context)* ctx = context; 145 struct sXd(primitive) prim; 146 struct sXd(attrib) attr; 147 const unsigned iprim_id = ivert / DIM; 148 const unsigned iprim_vert = ivert % DIM; 149 unsigned iprim; 150 ASSERT(pos && context); 151 152 iprim = (unsigned)ctx->primitives[iprim_id]; 153 SXD(scene_view_get_primitive(ctx->view, iprim, &prim)); 154 #if DIM == 2 155 s2d_segment_get_vertex_attrib(&prim, iprim_vert, SXD_POSITION, &attr); 156 #else 157 s3d_triangle_get_vertex_attrib(&prim, iprim_vert, SXD_POSITION, &attr); 158 #endif 159 ASSERT(attr.type == SXD_FLOATX); 160 fX(set)(pos, attr.value); 161 } 162 163 /******************************************************************************* 164 * Local functions 165 ******************************************************************************/ 166 static res_T 167 XD(solve_boundary) 168 (struct sdis_scene* scn, 169 const struct sdis_solve_boundary_args* args, 170 struct sdis_green_function** out_green, 171 struct sdis_estimator** out_estimator) 172 { 173 /* Time registration */ 174 struct time time0, time1; 175 char buf[128]; /* Temporary buffer used to store formated time */ 176 177 /* Device variables */ 178 struct mem_allocator* allocator = NULL; 179 size_t nthreads = 0; 180 181 /* Stardis variables */ 182 struct sdis_estimator* estimator = NULL; 183 struct sdis_green_function* green = NULL; 184 struct sdis_green_function** per_thread_green = NULL; 185 186 /* Random number generator */ 187 struct ssp_rng_proxy* rng_proxy = NULL; 188 struct ssp_rng** per_thread_rng = NULL; 189 190 /* XD variables */ 191 struct XD(boundary_context) ctx = XD(BOUNDARY_CONTEXT_NULL); 192 struct sXd(vertex_data) vdata = SXD_VERTEX_DATA_NULL; 193 struct sXd(scene)* scene = NULL; 194 struct sXd(shape)* shape = NULL; 195 struct sXd(scene_view)* view = NULL; 196 197 /* Miscellaneous */ 198 struct accum* per_thread_acc_temp = NULL; 199 struct accum* per_thread_acc_time = NULL; 200 size_t nrealisations = 0; 201 int64_t irealisation = 0; 202 size_t i; 203 int32_t* progress = NULL; /* Per process progress bar */ 204 int pcent_progress = 1; /* Percentage requiring progress update */ 205 int register_paths = SDIS_HEAT_PATH_NONE; 206 int is_master_process = 1; 207 ATOMIC nsolved_realisations = 0; 208 ATOMIC res = RES_OK; 209 210 if(!scn) { res = RES_BAD_ARG; goto error; } 211 if(!out_estimator && !out_green) { res = RES_BAD_ARG; goto error; } 212 res = check_solve_boundary_args(args); 213 if(res != RES_OK) goto error; 214 res = XD(scene_check_dimensionality)(scn); 215 if(res != RES_OK) goto error; 216 217 /* Check the submitted primitive indices */ 218 FOR_EACH(i, 0, args->nprimitives) { 219 res = scene_check_primitive_index(scn, args->primitives[i]); 220 if(res != RES_OK) goto error; 221 } 222 223 /* Check the submitted primitive sides */ 224 FOR_EACH(i, 0, args->nprimitives) { 225 if((unsigned)args->sides[i] >= SDIS_SIDE_NULL__) { 226 log_err(scn->dev, 227 "%s: invalid side for the primitive `%lu'.\n", 228 FUNC_NAME, (unsigned long)args->primitives[i]); 229 res = RES_BAD_ARG; 230 goto error; 231 } 232 } 233 234 if(out_green && args->picard_order != 1) { 235 log_err(scn->dev, "%s: the evaluation of the green function does not make " 236 "sense when dealing with the non-linearities of the system; i.e. picard " 237 "order must be set to 1 while it is currently set to %lu.\n", 238 FUNC_NAME, (unsigned long)args->picard_order); 239 res = RES_BAD_ARG; 240 goto error; 241 } 242 243 #ifdef SDIS_ENABLE_MPI 244 is_master_process = !scn->dev->use_mpi || scn->dev->mpi_rank == 0; 245 #endif 246 247 /* Update the progress bar every percent if escape sequences are allowed in 248 * log messages or only every 10 percent when only plain text is allowed. 249 * This reduces the number of lines of plain text printed */ 250 pcent_progress = scn->dev->no_escape_sequence ? 10 : 1; 251 252 /* Create the Star-XD shape of the boundary */ 253 #if SDIS_XD_DIMENSION == 2 254 res = s2d_shape_create_line_segments(scn->dev->sXd(dev), &shape); 255 #else 256 res = s3d_shape_create_mesh(scn->dev->sXd(dev), &shape); 257 #endif 258 if(res != RES_OK) goto error; 259 260 /* Initialise the boundary shape with the triangles/segments of the 261 * submitted primitives */ 262 ctx.primitives = args->primitives; 263 ctx.view = scn->sXd(view); 264 vdata.usage = SXD_POSITION; 265 vdata.get = XD(boundary_get_position); 266 #if SDIS_XD_DIMENSION == 2 267 vdata.type = S2D_FLOAT2; 268 res = s2d_line_segments_setup_indexed_vertices(shape, 269 (unsigned)args->nprimitives, boundary_get_indices_2d, 270 (unsigned)(args->nprimitives*2), &vdata, 1, &ctx); 271 #else 272 vdata.type = S3D_FLOAT3; 273 res = s3d_mesh_setup_indexed_vertices(shape, 274 (unsigned)args->nprimitives, boundary_get_indices_3d, 275 (unsigned)(args->nprimitives*3), &vdata, 1, &ctx); 276 #endif 277 if(res != RES_OK) goto error; 278 279 /* Create and setup the boundary Star-XD scene */ 280 res = sXd(scene_create)(scn->dev->sXd(dev), &scene); 281 if(res != RES_OK) goto error; 282 res = sXd(scene_attach_shape)(scene, shape); 283 if(res != RES_OK) goto error; 284 res = sXd(scene_view_create)(scene, SXD_SAMPLE, &view); 285 if(res != RES_OK) goto error; 286 287 nthreads = scn->dev->nthreads; 288 allocator = scn->dev->allocator; 289 290 /* Create the per thread RNGs */ 291 res = create_per_thread_rng 292 (scn->dev, args->rng_state, args->rng_type, &rng_proxy, &per_thread_rng); 293 if(res != RES_OK) goto error; 294 295 /* Allocate the per process progress status */ 296 res = alloc_process_progress(scn->dev, &progress); 297 if(res != RES_OK) goto error; 298 299 /* Create the per thread accumulators */ 300 per_thread_acc_temp = MEM_CALLOC(allocator, nthreads, sizeof(struct accum)); 301 per_thread_acc_time = MEM_CALLOC(allocator, nthreads, sizeof(struct accum)); 302 if(!per_thread_acc_temp) { res = RES_MEM_ERR; goto error; } 303 if(!per_thread_acc_time) { res = RES_MEM_ERR; goto error; } 304 305 /* Create the per thread green function */ 306 if(out_green) { 307 res = create_per_thread_green_function 308 (scn, args->signature, &per_thread_green); 309 if(res != RES_OK) goto error; 310 } 311 312 /* Create the estimator on the master process only. No estimator is needed 313 * for non master process */ 314 if(out_estimator && is_master_process) { 315 res = estimator_create(scn->dev, SDIS_ESTIMATOR_TEMPERATURE, &estimator); 316 if(res != RES_OK) goto error; 317 } 318 319 /* Synchronise the processes */ 320 process_barrier(scn->dev); 321 322 #define PROGRESS_MSG "Solving surface temperature: " 323 print_progress(scn->dev, progress, PROGRESS_MSG); 324 325 /* Begin time registration of the computation */ 326 time_current(&time0); 327 328 /* Here we go! Launch the Monte Carlo estimation */ 329 nrealisations = compute_process_realisations_count(scn->dev, args->nrealisations); 330 register_paths = out_estimator && is_master_process 331 ? args->register_paths : SDIS_HEAT_PATH_NONE; 332 omp_set_num_threads((int)scn->dev->nthreads); 333 #pragma omp parallel for schedule(static) 334 for(irealisation=0; irealisation<(int64_t)nrealisations; ++irealisation) { 335 struct boundary_realisation_args realis_args = BOUNDARY_REALISATION_ARGS_NULL; 336 struct time t0, t1; 337 const int ithread = omp_get_thread_num(); 338 struct ssp_rng* rng = per_thread_rng[ithread]; 339 struct accum* acc_temp = &per_thread_acc_temp[ithread]; 340 struct accum* acc_time = &per_thread_acc_time[ithread]; 341 struct green_path_handle* pgreen_path = NULL; 342 struct green_path_handle green_path = GREEN_PATH_HANDLE_NULL; 343 struct sdis_heat_path* pheat_path = NULL; 344 struct sdis_heat_path heat_path; 345 struct sXd(primitive) prim; 346 enum sdis_side side; 347 size_t iprim; 348 double w = NaN; 349 double uv[DIM-1]; 350 float st[DIM-1]; 351 double time; 352 size_t n; 353 int pcent; 354 res_T res_local = RES_OK; 355 res_T res_simul = RES_OK; 356 357 if(ATOMIC_GET(&res) != RES_OK) continue; /* An error occurred */ 358 359 /* Begin time registration */ 360 time_current(&t0); 361 362 time = sample_time(rng, args->time_range); 363 if(out_green) { 364 res_local = green_function_create_path 365 (per_thread_green[ithread], &green_path); 366 if(res_local != RES_OK) { 367 ATOMIC_SET(&res, res_local); 368 goto error_it; 369 } 370 pgreen_path = &green_path; 371 } 372 373 if(register_paths) { 374 heat_path_init(scn->dev->allocator, &heat_path); 375 pheat_path = &heat_path; 376 } 377 378 /* Sample a position onto the boundary */ 379 #if SDIS_XD_DIMENSION == 2 380 res_local = s2d_scene_view_sample 381 (view, 382 ssp_rng_canonical_float(rng), 383 ssp_rng_canonical_float(rng), 384 &prim, st); 385 uv[0] = (double)st[0]; 386 #else 387 res_local = s3d_scene_view_sample 388 (view, 389 ssp_rng_canonical_float(rng), 390 ssp_rng_canonical_float(rng), 391 ssp_rng_canonical_float(rng), 392 &prim, st); 393 d2_set_f2(uv, st); 394 #endif 395 if(res_local != RES_OK) { 396 ATOMIC_SET(&res, res_local); 397 goto error_it; 398 } 399 400 /* Map from boundary scene to sdis scene */ 401 ASSERT(prim.prim_id < args->nprimitives); 402 iprim = args->primitives[prim.prim_id]; 403 side = args->sides[prim.prim_id]; 404 405 /* Invoke the boundary realisation */ 406 realis_args.rng = rng; 407 realis_args.iprim = iprim; 408 realis_args.time = time; 409 realis_args.picard_order = args->picard_order; 410 realis_args.side = side; 411 realis_args.green_path = pgreen_path; 412 realis_args.heat_path = pheat_path; 413 realis_args.irealisation = (size_t)irealisation; 414 realis_args.diff_algo = args->diff_algo; 415 realis_args.uv[0] = uv[0]; 416 #if SDIS_XD_DIMENSION == 3 417 realis_args.uv[1] = uv[1]; 418 #endif 419 res_simul = XD(boundary_realisation)(scn, &realis_args, &w); 420 421 /* Fatal error */ 422 if(res_simul != RES_OK && res_simul != RES_BAD_OP) { 423 ATOMIC_SET(&res, res_simul); 424 goto error_it; 425 } 426 427 /* Register heat path */ 428 if(pheat_path) { 429 pheat_path->status = res_simul == RES_OK 430 ? SDIS_HEAT_PATH_SUCCESS 431 : SDIS_HEAT_PATH_FAILURE; 432 433 /* Check if the path must be saved regarding the register_paths mask */ 434 if(!(register_paths & (int)pheat_path->status)) { 435 heat_path_release(pheat_path); 436 pheat_path = NULL; 437 } else { /* Register the sampled path */ 438 res_local = estimator_add_and_release_heat_path(estimator, pheat_path); 439 if(res_local != RES_OK) { 440 ATOMIC_SET(&res, res_local); 441 goto error_it; 442 } 443 pheat_path = NULL; 444 } 445 } 446 447 /* Stop time registration */ 448 time_sub(&t0, time_current(&t1), &t0); 449 450 /* Update accumulators */ 451 if(res_simul == RES_OK) { 452 const double usec = (double)time_val(&t0, TIME_NSEC) * 0.001; 453 acc_temp->sum += w; acc_temp->sum2 += w*w; ++acc_temp->count; 454 acc_time->sum += usec; acc_time->sum2 += usec*usec; ++acc_time->count; 455 } 456 457 /* Update progress */ 458 n = (size_t)ATOMIC_INCR(&nsolved_realisations); 459 pcent = (int)((double)n * 100.0 / (double)nrealisations + 0.5/*round*/); 460 #pragma omp critical 461 if(pcent/pcent_progress > progress[0]/pcent_progress) { 462 progress[0] = pcent; 463 print_progress_update(scn->dev, progress, PROGRESS_MSG); 464 } 465 exit_it: 466 if(pheat_path) heat_path_release(pheat_path); 467 continue; 468 error_it: 469 goto exit_it; 470 } 471 /* Synchronise processes */ 472 process_barrier(scn->dev); 473 474 res = gather_res_T(scn->dev, (res_T)res); 475 if(res != RES_OK) goto error; 476 477 print_progress_completion(scn->dev, progress, PROGRESS_MSG); 478 #undef PROGRESS_MSG 479 480 /* Report computation time */ 481 time_sub(&time0, time_current(&time1), &time0); 482 time_dump(&time0, TIME_ALL, NULL, buf, sizeof(buf)); 483 log_info(scn->dev, "Surface probe temperature solved in %s.\n", buf); 484 485 /* Gather the RNG proxy sequence IDs and ensure that the RNG proxy state of 486 * the master process is greater than the RNG proxy state of all other 487 * processes */ 488 res = gather_rng_proxy_sequence_id(scn->dev, rng_proxy); 489 if(res != RES_OK) goto error; 490 491 /* Setup the estimated temperature */ 492 if(out_estimator) { 493 struct accum acc_temp; 494 struct accum acc_time; 495 496 time_current(&time0); 497 498 #define GATHER_ACCUMS(Msg, Acc) { \ 499 res = gather_accumulators(scn->dev, Msg, per_thread_##Acc, &Acc); \ 500 if(res != RES_OK) goto error; \ 501 } (void)0 502 GATHER_ACCUMS(MPI_SDIS_MSG_ACCUM_TEMP, acc_temp); 503 GATHER_ACCUMS(MPI_SDIS_MSG_ACCUM_TIME, acc_time); 504 #undef GATHER_ACCUMS 505 506 time_sub(&time0, time_current(&time1), &time0); 507 time_dump(&time0, TIME_ALL, NULL, buf, sizeof(buf)); 508 log_info(scn->dev, "Accumulators gathered in %s.\n", buf); 509 510 /* Return an estimator only on master process */ 511 if(is_master_process) { 512 ASSERT(acc_temp.count == acc_time.count); 513 estimator_setup_realisations_count(estimator, args->nrealisations, acc_temp.count); 514 estimator_setup_temperature(estimator, acc_temp.sum, acc_temp.sum2); 515 estimator_setup_realisation_time(estimator, acc_time.sum, acc_time.sum2); 516 res = estimator_save_rng_state(estimator, rng_proxy); 517 if(res != RES_OK) goto error; 518 } 519 } 520 521 /* Setup the green function */ 522 if(out_green) { 523 time_current(&time0); 524 525 res = gather_green_functions 526 (scn, rng_proxy, per_thread_green, per_thread_acc_time, &green); 527 if(res != RES_OK) goto error; 528 529 time_sub(&time0, time_current(&time1), &time0); 530 time_dump(&time0, TIME_ALL, NULL, buf, sizeof(buf)); 531 log_info(scn->dev, "Green functions gathered in %s.\n", buf); 532 533 /* Return a green function only on master process */ 534 if(!is_master_process) { 535 SDIS(green_function_ref_put(green)); 536 green = NULL; 537 } 538 } 539 540 exit: 541 if(per_thread_rng) release_per_thread_rng(scn->dev, per_thread_rng); 542 if(per_thread_green) release_per_thread_green_function(scn, per_thread_green); 543 if(progress) free_process_progress(scn->dev, progress); 544 if(per_thread_acc_temp) MEM_RM(scn->dev->allocator, per_thread_acc_temp); 545 if(per_thread_acc_time) MEM_RM(scn->dev->allocator, per_thread_acc_time); 546 if(scene) SXD(scene_ref_put(scene)); 547 if(shape) SXD(shape_ref_put(shape)); 548 if(view) SXD(scene_view_ref_put(view)); 549 if(rng_proxy) SSP(rng_proxy_ref_put(rng_proxy)); 550 if(out_green) *out_green = green; 551 if(out_estimator) *out_estimator = estimator; 552 return (res_T)res; 553 error: 554 if(estimator) { SDIS(estimator_ref_put(estimator)); estimator = NULL; } 555 if(green) { SDIS(green_function_ref_put(green)); green = NULL; } 556 goto exit; 557 } 558 559 static res_T 560 XD(solve_boundary_flux) 561 (struct sdis_scene* scn, 562 const struct sdis_solve_boundary_flux_args* args, 563 struct sdis_estimator** out_estimator) 564 { 565 /* Time registration */ 566 struct time time0, time1; 567 char buf[128]; /* Temporary buffer used to store formated time */ 568 569 /* Stardis variables */ 570 struct sdis_estimator* estimator = NULL; 571 572 /* XD variables */ 573 struct XD(boundary_context) ctx = XD(BOUNDARY_CONTEXT_NULL); 574 struct sXd(vertex_data) vdata = SXD_VERTEX_DATA_NULL; 575 struct sXd(scene)* scene = NULL; 576 struct sXd(shape)* shape = NULL; 577 struct sXd(scene_view)* view = NULL; 578 579 /* Random number generator */ 580 struct ssp_rng_proxy* rng_proxy = NULL; 581 struct ssp_rng** per_thread_rng = NULL; 582 583 /* Per thread accumulators */ 584 struct accum* per_thread_acc_tp = NULL; /* Temperature accumulator */ 585 struct accum* per_thread_acc_ti = NULL; /* Realisation time accumulator */ 586 struct accum* per_thread_acc_fl = NULL; /* Flux accumulator */ 587 struct accum* per_thread_acc_fc = NULL; /* Convective flux accumulator */ 588 struct accum* per_thread_acc_fr = NULL; /* Radiative flux accumulator */ 589 struct accum* per_thread_acc_fi = NULL; /* Imposed flux accumulator */ 590 591 /* Gathered accumulators */ 592 struct accum acc_tp = ACCUM_NULL; 593 struct accum acc_ti = ACCUM_NULL; 594 struct accum acc_fl = ACCUM_NULL; 595 struct accum acc_fc = ACCUM_NULL; 596 struct accum acc_fr = ACCUM_NULL; 597 struct accum acc_fi = ACCUM_NULL; 598 599 /* Miscellaneous */ 600 size_t nrealisations = 0; 601 int64_t irealisation; 602 size_t i; 603 int32_t* progress = NULL; /* Per process progress bar */ 604 int pcent_progress = 1; /* Percentage requiring progress update */ 605 int is_master_process = 1; 606 ATOMIC nsolved_realisations = 0; 607 ATOMIC res = RES_OK; 608 609 if(!scn || !out_estimator) { res = RES_BAD_ARG; goto error; } 610 611 res = check_solve_boundary_flux_args(args); 612 if(res != RES_OK) goto error; 613 res = XD(scene_check_dimensionality)(scn); 614 if(res != RES_OK) goto error; 615 616 FOR_EACH(i, 0, args->nprimitives) { 617 res = scene_check_primitive_index(scn, args->primitives[i]); 618 if(res != RES_OK) goto error; 619 } 620 621 /* Create the Star-XD shape of the boundary */ 622 #if SDIS_XD_DIMENSION == 2 623 res = s2d_shape_create_line_segments(scn->dev->sXd(dev), &shape); 624 #else 625 res = s3d_shape_create_mesh(scn->dev->sXd(dev), &shape); 626 #endif 627 if(res != RES_OK) goto error; 628 629 /* Initialise the boundary shape with the triangles/segments of the 630 * submitted primitives */ 631 ctx.primitives = args->primitives; 632 ctx.view = scn->sXd(view); 633 vdata.get = XD(boundary_get_position); 634 #if SDIS_XD_DIMENSION == 2 635 vdata.usage = S2D_POSITION; 636 vdata.type = S2D_FLOAT2; 637 res = s2d_line_segments_setup_indexed_vertices(shape, 638 (unsigned)args->nprimitives, XD(boundary_get_indices), 639 (unsigned)(args->nprimitives*2), &vdata, 1, &ctx); 640 #else /* DIM == 3 */ 641 vdata.usage = S3D_POSITION; 642 vdata.type = S3D_FLOAT3; 643 res = s3d_mesh_setup_indexed_vertices(shape, (unsigned)args->nprimitives, 644 XD(boundary_get_indices), (unsigned)(args->nprimitives*3), &vdata, 1, &ctx); 645 #endif 646 if(res != RES_OK) goto error; 647 648 /* Create and setup the boundary Star-XD scene */ 649 res = sXd(scene_create)(scn->dev->sXd(dev), &scene); 650 if(res != RES_OK) goto error; 651 res = sXd(scene_attach_shape)(scene, shape); 652 if(res != RES_OK) goto error; 653 res = sXd(scene_view_create)(scene, SXD_SAMPLE, &view); 654 if(res != RES_OK) goto error; 655 656 #ifdef SDIS_ENABLE_MPI 657 is_master_process = !scn->dev->use_mpi || scn->dev->mpi_rank == 0; 658 #endif 659 660 /* Update the progress bar every percent if escape sequences are allowed in 661 * log messages or only every 10 percent when only plain text is allowed. 662 * This reduces the number of lines of plain text printed */ 663 pcent_progress = scn->dev->no_escape_sequence ? 10 : 1; 664 665 /* Create the per thread RNGs */ 666 res = create_per_thread_rng 667 (scn->dev, args->rng_state, args->rng_type, &rng_proxy, &per_thread_rng); 668 if(res != RES_OK) goto error; 669 670 /* Allocate the per process progress status */ 671 res = alloc_process_progress(scn->dev, &progress); 672 if(res != RES_OK) goto error; 673 674 /* Create the per thread accumulator */ 675 #define ALLOC_ACCUMS(Dst) { \ 676 Dst = MEM_CALLOC(scn->dev->allocator, scn->dev->nthreads, sizeof(*Dst)); \ 677 if(!Dst) { res = RES_MEM_ERR; goto error; } \ 678 } (void)0 679 ALLOC_ACCUMS(per_thread_acc_tp); 680 ALLOC_ACCUMS(per_thread_acc_ti); 681 ALLOC_ACCUMS(per_thread_acc_fc); 682 ALLOC_ACCUMS(per_thread_acc_fl); 683 ALLOC_ACCUMS(per_thread_acc_fr); 684 ALLOC_ACCUMS(per_thread_acc_fi); 685 #undef ALLOC_ACCUMS 686 687 if(is_master_process) { 688 /* Create the estimator */ 689 res = estimator_create(scn->dev, SDIS_ESTIMATOR_FLUX, &estimator); 690 if(res != RES_OK) goto error; 691 } 692 693 /* Synchronise the processes */ 694 process_barrier(scn->dev); 695 696 #define PROGRESS_MSG "Solving surface flux: " 697 print_progress(scn->dev, progress, PROGRESS_MSG); 698 699 /* Begin time registration of the computation */ 700 time_current(&time0); 701 702 /* Here we go! Launch the Monte Carlo estimation */ 703 nrealisations = compute_process_realisations_count(scn->dev, args->nrealisations); 704 omp_set_num_threads((int)scn->dev->nthreads); 705 #pragma omp parallel for schedule(static) 706 for(irealisation = 0; irealisation < (int64_t)nrealisations; ++irealisation) { 707 struct boundary_flux_realisation_args realis_args = 708 BOUNDARY_FLUX_REALISATION_ARGS_NULL; 709 struct time t0, t1; 710 const int ithread = omp_get_thread_num(); 711 struct ssp_rng* rng = per_thread_rng[ithread]; 712 struct accum* acc_temp = &per_thread_acc_tp[ithread]; 713 struct accum* acc_time = &per_thread_acc_ti[ithread]; 714 struct accum* acc_flux = &per_thread_acc_fl[ithread]; 715 struct accum* acc_fcon = &per_thread_acc_fc[ithread]; 716 struct accum* acc_frad = &per_thread_acc_fr[ithread]; 717 struct accum* acc_fimp = &per_thread_acc_fi[ithread]; 718 struct sXd(primitive) prim; 719 struct sdis_interface_fragment frag = SDIS_INTERFACE_FRAGMENT_NULL; 720 const struct sdis_interface* interf; 721 const struct sdis_medium *fmd, *bmd; 722 enum sdis_side solid_side, fluid_side; 723 struct bound_flux_result result = BOUND_FLUX_RESULT_NULL__; 724 double epsilon, hc, hr, imposed_flux, imposed_temp, Tref; 725 size_t iprim; 726 double uv[DIM - 1]; 727 float st[DIM - 1]; 728 double time; 729 int flux_mask = 0; 730 size_t n; 731 int pcent; 732 res_T res_local = RES_OK; 733 res_T res_simul = RES_OK; 734 735 if(ATOMIC_GET(&res) != RES_OK) continue; /* An error occurred */ 736 737 /* Begin time registration */ 738 time_current(&t0); 739 740 time = sample_time(rng, args->time_range); 741 742 /* Sample a position onto the boundary */ 743 #if SDIS_XD_DIMENSION == 2 744 res_local = s2d_scene_view_sample 745 (view, 746 ssp_rng_canonical_float(rng), 747 ssp_rng_canonical_float(rng), 748 &prim, st); 749 uv[0] = (double)st[0]; 750 #else 751 res_local = s3d_scene_view_sample 752 (view, 753 ssp_rng_canonical_float(rng), 754 ssp_rng_canonical_float(rng), 755 ssp_rng_canonical_float(rng), 756 &prim, st); 757 d2_set_f2(uv, st); 758 #endif 759 if(res_local != RES_OK) { ATOMIC_SET(&res, res_local); continue; } 760 761 /* Map from boundary scene to sdis scene */ 762 ASSERT(prim.prim_id < args->nprimitives); 763 iprim = args->primitives[prim.prim_id]; 764 765 interf = scene_get_interface(scn, (unsigned)iprim); 766 fmd = interface_get_medium(interf, SDIS_FRONT); 767 bmd = interface_get_medium(interf, SDIS_BACK); 768 if(!fmd || !bmd 769 || ( !(fmd->type == SDIS_FLUID && bmd->type == SDIS_SOLID) 770 && !(fmd->type == SDIS_SOLID && bmd->type == SDIS_FLUID))) { 771 log_err(scn->dev, "%s: Attempt to compute a flux at a %s-%s interface.\n", 772 FUNC_NAME, 773 (fmd->type == SDIS_FLUID ? "fluid" : "solid"), 774 (bmd->type == SDIS_FLUID ? "fluid" : "solid")); 775 ATOMIC_SET(&res, RES_BAD_ARG); 776 continue; 777 } 778 solid_side = (fmd->type == SDIS_SOLID) ? SDIS_FRONT : SDIS_BACK; 779 fluid_side = (fmd->type == SDIS_FLUID) ? SDIS_FRONT : SDIS_BACK; 780 781 /* Build interface fragment on the fluid side of the primitive */ 782 res_local = XD(build_interface_fragment) 783 (&frag, scn, (unsigned)iprim, uv, fluid_side); 784 if(res_local!= RES_OK) { ATOMIC_SET(&res, res_local); continue; } 785 786 /* Fetch interface parameters */ 787 epsilon = interface_side_get_emissivity(interf, SDIS_INTERN_SOURCE_ID, &frag); 788 hc = interface_get_convection_coef(interf, &frag); 789 Tref = interface_side_get_reference_temperature(interf, &frag); 790 if(epsilon <= 0) { 791 hr = 0; 792 } else { 793 res_local = XD(check_Tref)(scn, frag.P, Tref, FUNC_NAME); 794 if(res_local != RES_OK) { ATOMIC_SET(&res, &res_local); continue; } 795 hr = 4.0 * BOLTZMANN_CONSTANT * Tref * Tref * Tref * epsilon; 796 } 797 798 frag.side = solid_side; 799 imposed_flux = interface_side_get_flux(interf, &frag); 800 imposed_temp = interface_side_get_temperature(interf, &frag); 801 if(SDIS_TEMPERATURE_IS_KNOWN(imposed_temp)) { 802 /* Flux computation on T boundaries is not supported yet */ 803 log_err(scn->dev, "%s: Attempt to compute a flux at a Dirichlet boundary " 804 "(not available yet).\n", FUNC_NAME); 805 ATOMIC_SET(&res, RES_BAD_ARG); 806 continue; 807 } 808 809 /* Fluid, Radiative and Solid temperatures */ 810 flux_mask = 0; 811 if(hr > 0) flux_mask |= FLUX_FLAG_RADIATIVE; 812 if(hc > 0) flux_mask |= FLUX_FLAG_CONVECTIVE; 813 814 /* Invoke the boundary flux realisation */ 815 realis_args.rng = rng; 816 realis_args.iprim = iprim; 817 realis_args.time = time; 818 realis_args.picard_order = args->picard_order; 819 realis_args.solid_side = solid_side; 820 realis_args.flux_mask = flux_mask; 821 realis_args.irealisation = (size_t)irealisation; 822 realis_args.diff_algo = args->diff_algo; 823 realis_args.uv[0] = uv[0]; 824 #if SDIS_XD_DIMENSION == 3 825 realis_args.uv[1] = uv[1]; 826 #endif 827 res_simul = XD(boundary_flux_realisation)(scn, &realis_args, &result); 828 829 /* Stop time registration */ 830 time_sub(&t0, time_current(&t1), &t0); 831 832 if(res_simul != RES_OK && res_simul != RES_BAD_OP) { 833 ATOMIC_SET(&res, res_simul); 834 continue; 835 } else if(res_simul == RES_OK) { /* Update accumulators */ 836 const double usec = (double)time_val(&t0, TIME_NSEC) * 0.001; 837 /* Convective flux from fluid to solid */ 838 const double w_conv = hc > 0 ? hc * (result.Tfluid - result.Tboundary) : 0; 839 /* Radiative flux from ambient to solid */ 840 const double w_rad = hr > 0 ? hr * (result.Tradiative - result.Tboundary) : 0; 841 /* Imposed flux that goes _into_ the solid */ 842 const double w_imp = (imposed_flux != SDIS_FLUX_NONE) ? imposed_flux : 0; 843 const double w_total = w_conv + w_rad + w_imp; 844 /* Temperature */ 845 acc_temp->sum += result.Tboundary; 846 acc_temp->sum2 += result.Tboundary*result.Tboundary; 847 ++acc_temp->count; 848 /* Time */ 849 acc_time->sum += usec; 850 acc_time->sum2 += usec*usec; 851 ++acc_time->count; 852 /* Overwall flux */ 853 acc_flux->sum += w_total; 854 acc_flux->sum2 += w_total*w_total; 855 ++acc_flux->count; 856 /* Convective flux */ 857 acc_fcon->sum += w_conv; 858 acc_fcon->sum2 += w_conv*w_conv; 859 ++acc_fcon->count; 860 /* Radiative flux */ 861 acc_frad->sum += w_rad; 862 acc_frad->sum2 += w_rad*w_rad; 863 ++acc_frad->count; 864 /* Imposed flux */ 865 acc_fimp->sum += w_imp; 866 acc_fimp->sum2 += w_imp*w_imp; 867 ++acc_fimp->count; 868 } 869 870 /* Update progress */ 871 n = (size_t)ATOMIC_INCR(&nsolved_realisations); 872 pcent = (int)((double)n * 100.0 / (double)nrealisations + 0.5/*round*/); 873 #pragma omp critical 874 if(pcent/pcent_progress > progress[0]/pcent_progress) { 875 progress[0] = pcent; 876 print_progress_update(scn->dev, progress, PROGRESS_MSG); 877 } 878 } 879 /* Synchronise processes */ 880 process_barrier(scn->dev); 881 882 res = gather_res_T(scn->dev, (res_T)res); 883 if(res != RES_OK) goto error; 884 885 print_progress_completion(scn->dev, progress, PROGRESS_MSG); 886 #undef PROGRESS_MSG 887 888 /* Report computation time */ 889 time_sub(&time0, time_current(&time1), &time0); 890 time_dump(&time0, TIME_ALL, NULL, buf, sizeof(buf)); 891 log_info(scn->dev, "Surface flux solved in %s.\n", buf); 892 893 /* Gather the RNG proxy sequence IDs and ensure that the RNG proxy state of 894 * the master process is greater than the RNG proxy state of all other 895 * processes */ 896 res = gather_rng_proxy_sequence_id(scn->dev, rng_proxy); 897 if(res != RES_OK) goto error; 898 899 time_current(&time0); 900 #define GATHER_ACCUMS(Msg, Acc) { \ 901 res = gather_accumulators(scn->dev, Msg, per_thread_##Acc, &Acc); \ 902 if(res != RES_OK) goto error; \ 903 } (void)0 904 GATHER_ACCUMS(MPI_SDIS_MSG_ACCUM_TEMP, acc_tp); 905 GATHER_ACCUMS(MPI_SDIS_MSG_ACCUM_TIME, acc_ti); 906 GATHER_ACCUMS(MPI_SDIS_MSG_ACCUM_FLUX_CONVECTIVE, acc_fc); 907 GATHER_ACCUMS(MPI_SDIS_MSG_ACCUM_FLUX_IMPOSED, acc_fi); 908 GATHER_ACCUMS(MPI_SDIS_MSG_ACCUM_FLUX_RADIATIVE, acc_fr); 909 GATHER_ACCUMS(MPI_SDIS_MSG_ACCUM_FLUX_TOTAL, acc_fl); 910 #undef GATHER_ACCUMS 911 912 time_sub(&time0, time_current(&time1), &time0); 913 time_dump(&time0, TIME_ALL, NULL, buf, sizeof(buf)); 914 log_info(scn->dev, "Accumulators gathered in %s.\n", buf); 915 916 /* Setup the estimated values */ 917 if(is_master_process) { 918 ASSERT(acc_tp.count == acc_fl.count); 919 ASSERT(acc_tp.count == acc_ti.count); 920 ASSERT(acc_tp.count == acc_fr.count); 921 ASSERT(acc_tp.count == acc_fc.count); 922 ASSERT(acc_tp.count == acc_fi.count); 923 estimator_setup_realisations_count(estimator, args->nrealisations, acc_tp.count); 924 estimator_setup_temperature(estimator, acc_tp.sum, acc_tp.sum2); 925 estimator_setup_realisation_time(estimator, acc_ti.sum, acc_ti.sum2); 926 estimator_setup_flux(estimator, FLUX_CONVECTIVE, acc_fc.sum, acc_fc.sum2); 927 estimator_setup_flux(estimator, FLUX_RADIATIVE, acc_fr.sum, acc_fr.sum2); 928 estimator_setup_flux(estimator, FLUX_IMPOSED, acc_fi.sum, acc_fi.sum2); 929 estimator_setup_flux(estimator, FLUX_TOTAL, acc_fl.sum, acc_fl.sum2); 930 931 res = estimator_save_rng_state(estimator, rng_proxy); 932 if(res != RES_OK) goto error; 933 } 934 935 exit: 936 if(per_thread_rng) release_per_thread_rng(scn->dev, per_thread_rng); 937 if(per_thread_acc_tp) MEM_RM(scn->dev->allocator, per_thread_acc_tp); 938 if(per_thread_acc_ti) MEM_RM(scn->dev->allocator, per_thread_acc_ti); 939 if(per_thread_acc_fc) MEM_RM(scn->dev->allocator, per_thread_acc_fc); 940 if(per_thread_acc_fr) MEM_RM(scn->dev->allocator, per_thread_acc_fr); 941 if(per_thread_acc_fl) MEM_RM(scn->dev->allocator, per_thread_acc_fl); 942 if(per_thread_acc_fi) MEM_RM(scn->dev->allocator, per_thread_acc_fi); 943 if(progress) free_process_progress(scn->dev, progress); 944 if(scene) SXD(scene_ref_put(scene)); 945 if(shape) SXD(shape_ref_put(shape)); 946 if(view) SXD(scene_view_ref_put(view)); 947 if(rng_proxy) SSP(rng_proxy_ref_put(rng_proxy)); 948 if(out_estimator) *out_estimator = estimator; 949 return (res_T)res; 950 error: 951 if(estimator) { SDIS(estimator_ref_put(estimator)); estimator = NULL; } 952 goto exit; 953 } 954 955 #include "sdis_Xd_end.h"