stardis-compute-probe-boundary.c (35321B)
1 /* Copyright (C) 2018-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 "stardis-app.h" 17 #include "stardis-compute.h" 18 #include "stardis-output.h" 19 #include "stardis-prog-properties.h" 20 #include "stardis-solid.h" 21 #include "stardis-solid-prog.h" 22 #include "stardis-ssconnect.h" 23 24 #include <sdis.h> 25 26 #include <star/senc3d.h> 27 #include <star/ssp.h> 28 29 #include <rsys/logger.h> 30 #include <rsys/str.h> 31 #include <rsys/clock_time.h> 32 33 #include <strings.h> 34 35 struct filter_ctx { 36 float distance; 37 enum sg3d_property_type side; 38 }; 39 #define FILTER_CTX_DEFAULT__ {0.f, SG3D_INTFACE} 40 static const struct filter_ctx FILTER_CTX_DEFAULT = FILTER_CTX_DEFAULT__; 41 42 /******************************************************************************* 43 * Helper functions 44 ******************************************************************************/ 45 static INLINE const char* 46 sdis_side_to_cstr(const enum sdis_side side) 47 { 48 const char* cstr = NULL; 49 switch(side) { 50 case SDIS_FRONT: cstr = "FRONT"; break; 51 case SDIS_BACK: cstr = "BACK"; break; 52 case SDIS_SIDE_NULL__: cstr = "UNDEFINED"; break; 53 default: FATAL("Unreachable code.\n"); 54 } 55 return cstr; 56 } 57 58 static res_T 59 read_rng_state 60 (struct stardis* stardis, 61 const char* filename, 62 struct ssp_rng* rng) 63 { 64 FILE* fp = NULL; 65 res_T res = RES_OK; 66 ASSERT(stardis && filename && rng); 67 68 if(stardis->mpi_initialized && stardis->mpi_rank != 0) { 69 goto exit; /* Non master process. Nothing to do */ 70 } 71 72 if((fp = fopen(filename, "r")) == NULL) { 73 logger_print(stardis->logger, LOG_ERROR, 74 "Cannot open generator's state file ('%s').\n", filename); 75 res = RES_IO_ERR; 76 goto error; 77 } 78 79 res = ssp_rng_read(rng, fp); 80 if(res != RES_OK) { 81 logger_print(stardis->logger, LOG_ERROR, 82 "Cannot read random generator's state ('%s').\n", filename); 83 goto error; 84 } 85 86 exit: 87 if(fp) CHK(fclose(fp) == 0); 88 return res; 89 error: 90 goto exit; 91 } 92 93 static res_T 94 write_rng_state 95 (struct stardis* stardis, 96 const char* filename, 97 struct ssp_rng* rng_state) 98 { 99 FILE* fp = NULL; 100 res_T res = RES_OK; 101 ASSERT(stardis && filename && rng_state); 102 103 if(stardis->mpi_initialized && stardis->mpi_rank != 0) { 104 goto exit; /* Non master process. Nothing to do */ 105 } 106 107 if((fp = fopen(filename, "wb")) == NULL) { 108 logger_print(stardis->logger, LOG_ERROR, 109 "Cannot open generator's state file ('%s').\n", filename); 110 res = RES_IO_ERR; 111 goto error; 112 } 113 114 res = ssp_rng_write(rng_state, fp); 115 if(res != RES_OK) { 116 logger_print(stardis->logger, LOG_ERROR, 117 "Cannot write random generator's state ('%s').\n", filename); 118 res = RES_IO_ERR; 119 goto error; 120 } 121 122 exit: 123 if(fp) CHK(fclose(fp) == 0); 124 return res; 125 error: 126 goto exit; 127 } 128 129 /* Filter used from a point query to determine not only one of the closest 130 * point, but the better one if there are more than one. In some circumstances 131 * it is not possible to determine the medium we are in using a given hit, but 132 * it is possible using another equidistant hit : 133 * 134 * 135 * P C 136 * +.............+---trg 1--- 137 * | 138 * medium 1 trg 2 medium 2 139 * | 140 * 141 * C is the closest point from P, and 2 different hits at C are possible (one 142 * on each triangle). However, only hit on trg 2 allows to find out that P is 143 * in medium 1 using sign(PC.Ntrg1) as PC.Ntrg2 = 0. 144 * The following filter function aims at selecting the hit on trg2 regardless 145 * of the order in which the 2 triangles are checked. 146 * One unexpected case cannot be decided though, but it implies that the 147 * closest triangle has 2 different media on its sides and that P lies on the 148 * triangle's plane : 149 * 150 * P C medium 1 151 * + +---trg--- 152 * medium 2 */ 153 static int 154 hit_filter 155 (const struct s3d_hit* hit, 156 const float ray_org[3], 157 const float ray_dir[3], 158 const float ray_range[2], 159 void* ray_data, 160 void* filter_data) 161 { 162 struct filter_ctx* filter_ctx = ray_data; 163 164 (void)ray_org, (void)ray_range, (void)filter_data; 165 ASSERT(hit && filter_ctx); 166 ASSERT(hit->uv[0] == CLAMP(hit->uv[0], 0, 1)); 167 ASSERT(hit->uv[1] == CLAMP(hit->uv[1], 0, 1)); 168 169 /* That's not the closest point. Keep the previous one if it can be used to 170 * detect the medium (i.e. side != SG3D_INTFACE) */ 171 if(filter_ctx->distance == hit->distance && filter_ctx->side != SG3D_INTFACE) { 172 return 1; /* Skip */ 173 } 174 175 filter_ctx->distance = hit->distance; 176 177 if(filter_ctx->distance == 0) { 178 filter_ctx->side = SG3D_INTFACE; 179 } else { 180 float sign = 0; 181 float N[3] = {0,0,0}; /* Normalized normal */ 182 183 /* Calculate the dot product with normalized vectors limits the numerical 184 * inaccuracies on its sign */ 185 f3_normalize(N, hit->normal); 186 sign = f3_dot(ray_dir, N); 187 188 /* Star3D hit normals are left-handed */ 189 if(sign < 0) filter_ctx->side = SG3D_FRONT; 190 else if(sign > 0) filter_ctx->side = SG3D_BACK; 191 else/*sign == 0*/ filter_ctx->side = SG3D_INTFACE; 192 } 193 194 return 0; /* Keep */ 195 } 196 197 static INLINE res_T 198 find_closest_point 199 (struct stardis* stardis, 200 const double pos[3], 201 struct filter_ctx* filter_ctx, 202 size_t* iprim, 203 double uv[2]) 204 { 205 struct sdis_scene_find_closest_point_args closest_pt_args = 206 SDIS_SCENE_FIND_CLOSEST_POINT_ARGS_NULL; 207 res_T res = RES_OK; 208 ASSERT(stardis && pos && filter_ctx && iprim && uv); 209 210 /* Find the surface point closest to the input position */ 211 closest_pt_args.position[0] = pos[0]; 212 closest_pt_args.position[1] = pos[1]; 213 closest_pt_args.position[2] = pos[2]; 214 closest_pt_args.radius = INF; 215 closest_pt_args.filter_3d = hit_filter; 216 closest_pt_args.filter_data = filter_ctx; 217 ERR(sdis_scene_find_closest_point 218 (stardis->sdis_scn , &closest_pt_args, iprim, uv)); 219 if(*iprim == SDIS_PRIMITIVE_NONE) { 220 res = RES_BAD_ARG; 221 goto error; 222 } 223 224 exit: 225 return res; 226 error: 227 goto exit; 228 } 229 230 static res_T 231 check_move_to_solid_boundary 232 (const struct stardis* stardis, 233 const double pos[3], /* Original position */ 234 const double time, /* [s] */ 235 const struct description* desc, /* Solid medium in which pos lies */ 236 const size_t iprim, /* Triangle index to which to move */ 237 const double uv[2], /* Triangle coordinates to which to move */ 238 const double distance, /* Move distance */ 239 const int advice) 240 { 241 struct stardis_vertex vtx = STARDIS_VERTEX_NULL; 242 const char* prefix = ""; 243 const char* solid_name = ""; 244 double delta = 0; 245 res_T res = RES_OK; 246 247 /* Check pre-conditions */ 248 ASSERT(stardis && pos && time > 0 && desc && uv && distance >= 0); 249 250 /* Retrieve the delta and define the prefix of the solid for log messages */ 251 switch(desc->type) { 252 /* Regular solid, i.e. solid with constant properties */ 253 case DESC_MAT_SOLID: 254 delta = desc->d.solid->delta; 255 prefix = ""; 256 break; 257 258 /* Solid with programmed properties */ 259 case DESC_MAT_SOLID_PROG: 260 vtx.P[0] = pos[0]; 261 vtx.P[1] = pos[1]; 262 vtx.P[2] = pos[2]; 263 vtx.time = time; 264 delta = desc->d.solid_prog->delta(&vtx, desc->d.solid_prog->prog_data); 265 prefix = "programmed "; 266 break; 267 268 default: FATAL("Unreachable code.\n"); 269 } 270 271 solid_name = str_cget(get_description_name(desc)); 272 logger_print(stardis->logger, LOG_OUTPUT, "Probe was in %ssolid '%s'.\n", 273 prefix, solid_name); 274 275 /* The position is close from the triangle */ 276 if(distance < 0.5*delta) { 277 logger_print(stardis->logger, LOG_OUTPUT, 278 "Probe was %g delta from closest boundary.\n", 279 distance/delta); 280 281 /* Notice that the position is a little far from the triangle */ 282 } else if(distance < 2*delta) { 283 logger_print(stardis->logger, LOG_WARNING, 284 "Probe was %g delta from closest boundary.%s\n", 285 distance/delta, 286 (advice ? " Consider using -p instead of -P.\n" : "")); 287 288 /* The position is too far from the triangle */ 289 } else { 290 logger_print(stardis->logger, LOG_ERROR, 291 "Probe moved to (%g, %g, %g), primitive %lu, uv = (%g, %g). " 292 "Move is %g delta long. Use -p instead of -P.\n", 293 SPLIT3(pos), (unsigned long)iprim, SPLIT2(uv), distance/delta); 294 res = RES_BAD_ARG; 295 goto error; 296 } 297 298 exit: 299 return res; 300 error: 301 goto exit; 302 } 303 304 /* This function checks nothing. It only records the status. It is named as the 305 * one used to check the projection on the solid limit to make it symmetrical, 306 * and thus simplify the reading of sources */ 307 static res_T 308 check_move_to_fluid_boundary 309 (struct stardis* stardis, 310 const struct description* desc, /* Fluid medium in which pos lies */ 311 const double distance) /* Move distance */ 312 { 313 const char* prefix = ""; 314 const char* fluid_name = ""; 315 316 ASSERT(stardis && desc && distance >= 0); 317 318 switch(desc->type) { 319 case DESC_MAT_FLUID: prefix = ""; break; 320 case DESC_MAT_FLUID_PROG: prefix = "programmed "; break; 321 default: FATAL("Unreachable code.\n"); 322 } 323 324 fluid_name = str_cget(get_description_name(desc)); 325 logger_print(stardis->logger, LOG_OUTPUT, 326 "Probe was in %sfluid '%s'.\n", prefix, fluid_name); 327 logger_print(stardis->logger, LOG_OUTPUT, 328 "Probe distance from closest boundary was %g.\n", distance); 329 330 return RES_OK; 331 } 332 333 static res_T 334 move_to_boundary 335 (struct stardis* stardis, 336 const double pos[3], 337 const double time, /* [s] */ 338 const int is_probe_temp_computation, 339 size_t* out_iprim, 340 double uv[2]) 341 { 342 /* Position on boundary */ 343 struct filter_ctx filter_ctx = FILTER_CTX_DEFAULT; 344 double proj_pos[3] = {0,0,0}; 345 size_t iprim = 0; 346 347 /* Miscellaneous */ 348 size_t nvertices_close = 0; 349 res_T res = RES_OK; 350 351 /* Check pre-conditions */ 352 ASSERT(stardis && pos && time >= 0 && out_iprim && uv); 353 354 ERR(find_closest_point(stardis, pos, &filter_ctx, &iprim, uv)); 355 356 if(filter_ctx.side != SG3D_INTFACE) { 357 /* Properties */ 358 const struct description* desc_list = NULL; 359 const struct description* desc = NULL; 360 unsigned desc_ids[SG3D_PROP_TYPES_COUNT__]; 361 362 SG3D(geometry_get_unique_triangle_properties 363 (stardis->geometry.sg3d, (unsigned)iprim, desc_ids)); 364 365 desc_list = darray_descriptions_cdata_get(&stardis->descriptions); 366 367 /* Probe is outside the system */ 368 if(desc_ids[filter_ctx.side] == SG3D_UNSPECIFIED_PROPERTY) { 369 logger_print(stardis->logger, LOG_WARNING, 370 "Probe was outside the system.\n"); 371 372 /* Probe is in a medium */ 373 } else { 374 desc = desc_list + desc_ids[filter_ctx.side]; 375 376 switch(desc->type) { 377 case DESC_MAT_SOLID: 378 case DESC_MAT_SOLID_PROG: 379 ERR(check_move_to_solid_boundary 380 (stardis, pos, time, desc, iprim, uv, filter_ctx.distance, 381 is_probe_temp_computation)); 382 break; 383 case DESC_MAT_FLUID: 384 case DESC_MAT_FLUID_PROG: 385 ERR(check_move_to_fluid_boundary 386 (stardis, desc, filter_ctx.distance)); 387 break; 388 default: FATAL("Unreachable code.\n"); 389 } 390 } 391 } 392 393 SDIS(scene_get_boundary_position(stardis->sdis_scn, iprim, uv, proj_pos)); 394 395 /* Count the number of vertices that are close to the boundary position 396 * and issue a warning if necessary */ 397 nvertices_close += CLAMP(uv[0], 0.0005, 0.9995) != uv[0]; 398 nvertices_close += CLAMP(uv[1], 0.0005, 0.9995) != uv[1]; 399 if(nvertices_close) { 400 logger_print(stardis->logger, LOG_WARNING, 401 "Probe %s close to / on %s. " 402 "If computation fails, try moving it slightly.\n", 403 filter_ctx.distance == 0 ? "is" : "moved", 404 nvertices_close == 1 ? "an edge" : "a vertex"); 405 } 406 407 /* Probe is on a boundary */ 408 if(filter_ctx.distance == 0) { 409 logger_print(stardis->logger, LOG_OUTPUT, 410 "Probe is on primitive %lu, uv = (%g, %g), not moved.\n", 411 (unsigned long)iprim, SPLIT2(uv)); 412 413 /* Probe was projected on a boundary */ 414 } else { 415 logger_print(stardis->logger, LOG_OUTPUT, 416 "Probe moved to (%g, %g, %g), primitive %lu, uv = (%g, %g).\n", 417 SPLIT3(proj_pos), (unsigned long)iprim, SPLIT2(uv)); 418 } 419 420 exit: 421 *out_iprim = iprim; 422 return res; 423 error: 424 goto exit; 425 } 426 427 static res_T 428 setup_probe_side 429 (struct stardis* stardis, 430 const unsigned desc_ids[SG3D_PROP_TYPES_COUNT__], 431 const char* side_str, 432 const size_t iprim, 433 enum sdis_side *out_side) 434 { 435 const struct description* desc_list = NULL; 436 const struct description* desc_front = NULL; 437 const struct description* desc_back = NULL; 438 size_t ndescs = 0; 439 enum sdis_side side = SDIS_SIDE_NULL__; 440 res_T res = RES_OK; 441 (void)ndescs; /* Avoid "Unused variable" warnings in release */ 442 443 /* Check pre-conditions */ 444 ASSERT(stardis && side_str && desc_ids && out_side); 445 446 /* Fetch the properties */ 447 desc_list = darray_descriptions_cdata_get(&stardis->descriptions); 448 ndescs = darray_descriptions_size_get(&stardis->descriptions); 449 desc_front = desc_list + desc_ids[SG3D_FRONT]; 450 desc_back = desc_list + desc_ids[SG3D_BACK]; 451 452 /* No side specified */ 453 if(!side_str || !strlen(side_str)) { 454 side = SDIS_SIDE_NULL__; 455 456 /* Set probe to front side */ 457 } else if(!strcasecmp(side_str, "FRONT")) { 458 ASSERT(desc_ids[SG3D_FRONT] < ndescs && DESC_IS_MEDIUM(desc_front)); 459 side = SDIS_FRONT; 460 461 /* Set probe to back side */ 462 } else if(!strcasecmp(side_str, "BACK")) { 463 ASSERT(desc_ids[SG3D_BACK] < ndescs && DESC_IS_MEDIUM(desc_back)); 464 side = SDIS_BACK; 465 466 /* Set the probe to the side that points to the submitted medium name */ 467 } else { 468 unsigned med_id_probe = 0; /* Medium defined on the probe */ 469 unsigned med_id_front = 0; /* Medium on front side */ 470 unsigned med_id_back = 0; /* Medium on back side */ 471 ASSERT(DESC_IS_MEDIUM(desc_front) && DESC_IS_MEDIUM(desc_back)); 472 473 if(!find_medium_by_name(stardis, side_str, &med_id_probe)) { 474 logger_print(stardis->logger, LOG_ERROR, 475 "Cannot locate side from medium name '%s' (unknown medium)\n", 476 side_str); 477 res = RES_BAD_ARG; 478 goto error; 479 } 480 481 description_get_medium_id(desc_front, &med_id_front); 482 description_get_medium_id(desc_back, &med_id_back); 483 484 /* Invalid probe medium wrt the boundary on which it is located */ 485 if(med_id_probe != med_id_front 486 && med_id_probe != med_id_back) { 487 logger_print(stardis->logger, LOG_ERROR, 488 "Medium '%s' is not used at this interface (prim id=%lu)\n", 489 side_str, (unsigned long)iprim); 490 res = RES_BAD_ARG; 491 goto error; 492 } 493 494 /* The same medium is used on both sides: cannot differentiate */ 495 if(med_id_front == med_id_back) { 496 unsigned encs[2]; /* Identifier of the enclosures */ 497 498 ERR(senc3d_scene_get_triangle_enclosures 499 (stardis->senc3d_scn, (unsigned)iprim, encs)); 500 logger_print(stardis->logger, LOG_ERROR, 501 "Medium '%s' is used on both sides of this interface (prim id=%lu).\n", 502 side_str, (unsigned long)iprim); 503 logger_print(stardis->logger, LOG_ERROR, 504 "Side must be defined using either FRONT or BACK.\n"); 505 logger_print(stardis->logger, LOG_ERROR, 506 "FRONT side is related to enclosure %u, BACK side to enclosure %u.\n", 507 encs[SENC3D_FRONT], encs[SENC3D_BACK]); 508 509 res = RES_BAD_ARG; 510 goto error; 511 } 512 513 side = med_id_probe == med_id_front ? SDIS_FRONT : SDIS_BACK; 514 } 515 516 exit: 517 *out_side = side; 518 return res; 519 error: 520 side = SDIS_SIDE_NULL__; 521 goto exit; 522 } 523 524 /* This function checks the conformity between the potential thermal contact 525 * resistance at the probe location and the specified probe side. */ 526 static res_T 527 setup_thermal_contact_resistance 528 (struct stardis* stardis, 529 const unsigned desc_ids[SG3D_PROP_TYPES_COUNT__], 530 const enum sdis_side probe_side) 531 { 532 struct str log_msg; 533 const struct description* desc_list = NULL; 534 const struct description* desc_front = NULL; 535 const struct description* desc_back = NULL; 536 const struct description* desc_intface = NULL; 537 size_t ndescs = 0; 538 double tcr = 0; 539 res_T res = RES_OK; 540 (void)ndescs; /* Avoid "Unused variable" warnings in release */ 541 542 /* Check pre-conditions */ 543 ASSERT(stardis && desc_ids); 544 545 str_init(stardis->allocator, &log_msg); 546 547 /* Fetch the properties */ 548 desc_list = darray_descriptions_cdata_get(&stardis->descriptions); 549 ndescs = darray_descriptions_size_get(&stardis->descriptions); 550 desc_front = desc_list + desc_ids[SG3D_FRONT]; 551 desc_back = desc_list + desc_ids[SG3D_BACK]; 552 desc_intface = desc_list + desc_ids[SG3D_INTFACE]; 553 554 /* Get the thermal contact resistance between solid/solid connection if any */ 555 if(desc_ids[SG3D_INTFACE] != SG3D_UNSPECIFIED_PROPERTY 556 && desc_intface->type == DESC_SOLID_SOLID_CONNECT) { 557 ASSERT(desc_ids[SG3D_INTFACE] < ndescs); 558 tcr = desc_intface->d.ss_connect->tcr; 559 } 560 561 /* Warn if side defined and no resistance defined */ 562 if(tcr == 0 && probe_side != SDIS_SIDE_NULL__) { 563 logger_print(stardis->logger, LOG_WARNING, 564 "Specifying a compute side at an interface with no contact resistance " 565 "is meaningless.\n"); 566 } 567 568 #define GET_DESC_NAME(Desc) str_cget(get_description_name(Desc)) 569 570 /* A thermal contact resistance cannot be defined if probe side is NULL */ 571 if(tcr != 0 && probe_side == SDIS_SIDE_NULL__) { 572 logger_print(stardis->logger, LOG_ERROR, 573 "Cannot let probe computation side unspecified on an interface with a " 574 "non-nul thermal resistance.\n"); 575 576 /* Format the log string */ 577 if(desc_ids[SG3D_FRONT] != SG3D_UNSPECIFIED_PROPERTY) { 578 ASSERT(desc_ids[SG3D_FRONT] < ndescs); 579 ERR(str_append_printf(&log_msg, " FRONT: '%s'", GET_DESC_NAME(desc_front))); 580 } 581 if(desc_ids[SG3D_FRONT] != SG3D_UNSPECIFIED_PROPERTY 582 && desc_ids[SG3D_BACK] != SG3D_UNSPECIFIED_PROPERTY) { 583 ERR(str_append_char(&log_msg, ',')); 584 } 585 if(desc_ids[SG3D_BACK] != SG3D_UNSPECIFIED_PROPERTY) { 586 ASSERT(desc_ids[SG3D_BACK] < ndescs); 587 ERR(str_append_printf(&log_msg, " BACK: '%s'", GET_DESC_NAME(desc_back))); 588 } 589 590 /* Print error message */ 591 logger_print(stardis->logger, LOG_ERROR, 592 "Interface '%s',%s, resistance=%g K.m^2/W.\n", 593 GET_DESC_NAME(desc_intface), str_cget(&log_msg), tcr); 594 595 res = RES_BAD_ARG; 596 goto error; 597 } 598 599 /* Log that a calculation is done on a boundary with tcr */ 600 if(tcr > 0) { 601 const char* medium_name = probe_side == SDIS_FRONT 602 ? GET_DESC_NAME(desc_front) 603 : GET_DESC_NAME(desc_back); 604 605 logger_print(stardis->logger, LOG_OUTPUT, 606 "Probe computation on an interface with a thermal resistance = %g K.m^2/W " 607 "on %s side (medium is '%s').\n", 608 tcr, sdis_side_to_cstr(probe_side), medium_name); 609 } 610 611 #undef GET_DESC_NAME 612 613 exit: 614 str_release(&log_msg); 615 return res; 616 error: 617 goto exit; 618 } 619 620 static res_T 621 solve 622 (struct stardis* stardis, 623 struct time* start, 624 struct sdis_solve_probe_boundary_args* args, 625 res_T output[2]) 626 { 627 struct time t0, t1; 628 struct sdis_mc time = SDIS_MC_NULL; 629 struct dump_path_context ctx = DUMP_PATH_CONTEXT_NULL; 630 struct sdis_estimator* estimator = NULL; 631 const struct str* rng_in = NULL; 632 const struct str* rng_out = NULL; 633 struct ssp_rng* rng = NULL; 634 int is_master_process = 0; 635 res_T res = RES_OK; 636 ASSERT(stardis && args && output); 637 638 is_master_process = !stardis->mpi_initialized || stardis->mpi_rank == 0; 639 640 rng_in = &stardis->rndgen_state_in_filename; 641 rng_out = &stardis->rndgen_state_out_filename; 642 643 /* Read RNG state from file */ 644 if(!str_is_empty(rng_in)) { 645 ERR(ssp_rng_create(stardis->allocator, SSP_RNG_THREEFRY, &rng)); 646 ERR(read_rng_state(stardis, str_cget(rng_in), rng)); 647 args->rng_state = rng; 648 } 649 650 /* Run the calculation */ 651 time_current(&t0); 652 ERR(sdis_solve_probe_boundary(stardis->sdis_scn, args, &estimator)); 653 time_current(&t1); 654 655 /* No more to do for non master processes */ 656 if(!is_master_process) goto exit; 657 658 /* Per per realisation time */ 659 ERR(sdis_estimator_get_realisation_time(estimator, &time)); 660 ERR(print_computation_time(&time, stardis, start, &t0, &t1, NULL)); 661 662 /* Write outputs and save their status */ 663 ctx.stardis = stardis; 664 ctx.rank = 0; 665 output[0] = print_single_MC_result(estimator, stardis, stdout); 666 output[1] = sdis_estimator_for_each_path(estimator, dump_path, &ctx); 667 668 /* Write the resulting RNG state to a file */ 669 if(!str_is_empty(rng_out)) { 670 struct ssp_rng* rng_state = NULL; 671 ERR(sdis_estimator_get_rng_state(estimator, &rng_state)); 672 ERR(write_rng_state(stardis, str_cget(rng_out), rng_state)); 673 } 674 675 exit: 676 if(estimator) SDIS(estimator_ref_put(estimator)); 677 if(rng) SSP(rng_ref_put(rng)); 678 return res; 679 error: 680 goto exit; 681 } 682 683 static res_T 684 solve_list 685 (struct stardis* stardis, 686 struct time* start, 687 struct sdis_solve_probe_boundary_list_args* args, 688 res_T* output) 689 { 690 struct time t0, t1; 691 struct sdis_mc time = SDIS_MC_NULL; /* Time per realisation */ 692 struct sdis_estimator_buffer* buffer = NULL; 693 size_t i = 0; 694 size_t def[2] = {0, 0}; 695 int is_master_process = 0; 696 res_T res = RES_OK; 697 ASSERT(stardis && start && args); 698 699 is_master_process = !stardis->mpi_initialized || stardis->mpi_rank == 0; 700 701 /* Run the calculation */ 702 time_current(&t0); 703 ERR(sdis_solve_probe_boundary_list(stardis->sdis_scn, args, &buffer)); 704 time_current(&t1); 705 706 /* No more to do for non master processes */ 707 if(!is_master_process) goto exit; 708 709 /* Retrieve the number of solved probes */ 710 ERR(sdis_estimator_buffer_get_definition(buffer, def)); 711 ASSERT(def[0] == darray_probe_boundary_size_get(&stardis->probe_boundary_list)); 712 ASSERT(def[1] == 1); 713 714 ERR(sdis_estimator_buffer_get_realisation_time(buffer, &time)); 715 ERR(print_computation_time(&time, stardis, start, &t0, &t1, NULL)); 716 717 /* Print the estimated temperature of each probe */ 718 FOR_EACH(i, 0, def[0]) { 719 const struct stardis_probe_boundary* probe = NULL; 720 const struct sdis_estimator* estimator = NULL; 721 res_T res2 = RES_OK; 722 723 probe = darray_probe_boundary_cdata_get(&stardis->probe_boundary_list) + i; 724 ERR(sdis_estimator_buffer_at(buffer, i, 0, &estimator)); 725 726 res2 = print_single_MC_result_probe_boundary 727 (stardis, probe, estimator, stdout); 728 if(res2 != RES_OK && *output == RES_OK) *output = res2; 729 } 730 731 exit: 732 if(buffer) SDIS(estimator_buffer_ref_put(buffer)); 733 return res; 734 error: 735 goto exit; 736 } 737 738 static res_T 739 setup_solve_probe_boundary_flux_args 740 (struct stardis* stardis, 741 const struct stardis_probe_boundary* probe, 742 struct sdis_solve_probe_boundary_flux_args* args) 743 { 744 double uv[2] = {0, 0}; 745 size_t iprim = SIZE_MAX; 746 unsigned desc_ids[SG3D_PROP_TYPES_COUNT__]; 747 res_T res = RES_OK; 748 ASSERT(stardis && probe && args); 749 750 /* Calculate the probe position on the boundary */ 751 ERR(move_to_boundary(stardis, probe->position, probe->time[0], 0, &iprim, uv)); 752 753 ERR(sg3d_geometry_get_unique_triangle_properties(stardis->geometry.sg3d, 754 (unsigned)iprim, desc_ids)); 755 756 d3_set(stardis->probe, probe->position); 757 /* Setup of solve input parameters */ 758 759 args->nrealisations = stardis->samples; 760 args->iprim = iprim; 761 d2_set(args->uv, uv); 762 args->picard_order = stardis->picard_order; 763 d2_set(args->time_range, stardis->time_range); 764 args->diff_algo = stardis->diff_algo; 765 766 exit: 767 return res; 768 error: 769 goto exit; 770 } 771 772 static res_T 773 solve_flux_list 774 (struct stardis* stardis, 775 struct time* start, 776 const struct stardis_probe_boundary* probes, 777 size_t nprobes) 778 { 779 struct time t0, t1; 780 struct sdis_mc time = SDIS_MC_NULL; /* Time per realisation */ 781 struct sdis_estimator_buffer* buffer = NULL; 782 size_t i = 0; 783 struct sdis_estimator* estimator = NULL; 784 FILE* stream_r = NULL; 785 struct sdis_solve_probe_boundary_flux_args args; 786 res_T res = RES_OK; 787 ASSERT(stardis && start); 788 789 /* Input random state? */ 790 READ_RANDOM_STATE(&stardis->rndgen_state_in_filename); 791 792 time_current(&t0); 793 for(i = 0; i < nprobes; i++) { 794 args = SDIS_SOLVE_PROBE_BOUNDARY_FLUX_ARGS_DEFAULT; 795 ERR(setup_solve_probe_boundary_flux_args(stardis, &probes[i], &args)); 796 /* Run the calculation */ 797 ERR(sdis_solve_probe_boundary_flux(stardis->sdis_scn, &args, &estimator)); 798 /* Print the estimated temperature */ 799 ERR(print_single_MC_result(estimator, stardis, stdout)); 800 ERR(sdis_estimator_ref_put(estimator)); 801 } 802 time_current(&t1); 803 804 /* Output random state? */ 805 WRITE_RANDOM_STATE(&stardis->rndgen_state_out_filename); 806 807 ERR(print_computation_time(&time, stardis, start, &t0, &t1, NULL)); 808 809 exit: 810 if(buffer) SDIS(estimator_buffer_ref_put(buffer)); 811 return res; 812 error: 813 goto exit; 814 } 815 816 static res_T 817 solve_green 818 (struct stardis* stardis, 819 struct time* start, 820 struct sdis_solve_probe_boundary_args* args) 821 { 822 struct time t0/*calculation start*/, t1/*calculation end*/, t2/*output end*/; 823 struct sdis_green_function* green = NULL; 824 FILE* fp_green = NULL; 825 FILE* fp_path = NULL; 826 const struct str* rng_in = NULL; 827 struct ssp_rng* rng = NULL; 828 int is_master_process = 0; 829 res_T res = RES_OK; 830 831 ASSERT(stardis && args); 832 833 is_master_process = !stardis->mpi_initialized || stardis->mpi_rank == 0; 834 835 rng_in = &stardis->rndgen_state_in_filename; 836 837 /* Read RNG state from file */ 838 if(!str_is_empty(rng_in)) { 839 ERR(ssp_rng_create(stardis->allocator, SSP_RNG_THREEFRY, &rng)); 840 ERR(read_rng_state(stardis, str_cget(rng_in), rng)); 841 args->rng_state = rng; 842 } 843 844 /* Try to open output files to detect errors early */ 845 if(is_master_process && (stardis->mode & MODE_GREEN_BIN)) { 846 const char* green_filename = str_cget(&stardis->bin_green_filename); 847 const char* path_filename = str_cget(&stardis->end_paths_filename); 848 849 if((fp_green = fopen(green_filename, "wb")) == NULL) { 850 logger_print(stardis->logger, LOG_ERROR, 851 "Cannot open file '%s' for binary writing.\n", green_filename); 852 res = RES_IO_ERR; 853 goto error; 854 } 855 856 if(strlen(path_filename) != 0) { 857 if((fp_path = fopen(path_filename, "w")) == NULL) { 858 logger_print(stardis->logger, LOG_ERROR, 859 "Cannot open file '%s' for writing.\n", path_filename); 860 res = RES_IO_ERR; 861 goto error; 862 } 863 } 864 } 865 866 /* Run the Green estimation */ 867 time_current(&t0); /* Calculation starts */ 868 ERR(sdis_solve_probe_boundary_green_function(stardis->sdis_scn, args, &green)); 869 time_current(&t1); /* Calculation ends */ 870 871 /* No more to do for non master processes */ 872 if(is_master_process) goto exit; 873 874 /* Write ASCII Green */ 875 if(stardis->mode & MODE_GREEN_ASCII) { 876 ERR(dump_green_ascii(green, stardis, stdout)); 877 } 878 879 /* Write binary Green */ 880 if(stardis->mode & MODE_GREEN_BIN) { 881 ERR(dump_green_bin(green, stardis, fp_green)); 882 if(fp_path) { 883 ERR(dump_paths_end(green, stardis, fp_path)); 884 } 885 } 886 887 time_current(&t2); /* Output ends */ 888 889 ERR(print_computation_time(NULL, stardis, start, &t0, &t1, &t2)); 890 891 /* Note that the resulting RNG state is not written in an output file because 892 * the solver API does not provide a function to recover it. But in fact, the 893 * green function saves the RNG state after its estimation. Therefore, the API 894 * can be expected to provide such functionality soon. 895 * 896 * TODO write the RNG status of the Green function when it is available */ 897 898 exit: 899 if(fp_green) CHK(fclose(fp_green) == 0); 900 if(fp_path) CHK(fclose(fp_path) == 0); 901 if(green) SDIS(green_function_ref_put(green)); 902 if(rng) SSP(rng_ref_put(rng)); 903 return res; 904 error: 905 goto exit; 906 } 907 908 static res_T 909 compute_single_flux_probe_on_interface 910 (struct stardis* stardis, 911 struct time* start, 912 const struct stardis_probe_boundary* probe) 913 { 914 res_T res = RES_OK; 915 struct sdis_green_function* green = NULL; 916 struct sdis_estimator* estimator = NULL; 917 struct sdis_solve_probe_boundary_flux_args args 918 = SDIS_SOLVE_PROBE_BOUNDARY_FLUX_ARGS_DEFAULT; 919 FILE* stream_r = NULL; 920 struct time compute_start, compute_end; 921 922 ASSERT(stardis && start && probe); 923 924 ERR(setup_solve_probe_boundary_flux_args(stardis, probe, &args)); 925 926 /* Input random state? */ 927 READ_RANDOM_STATE(&stardis->rndgen_state_in_filename); 928 929 if(stardis->mpi_initialized && stardis->mpi_rank != 0) { 930 ERR(sdis_solve_probe_boundary_flux(stardis->sdis_scn, &args, &estimator)); 931 } else { 932 struct sdis_mc time = SDIS_MC_NULL; 933 time_current(&compute_start); 934 ERR(sdis_solve_probe_boundary_flux(stardis->sdis_scn, &args, &estimator)); 935 time_current(&compute_end); 936 ERR(sdis_estimator_get_realisation_time(estimator, &time)); 937 ERR(print_computation_time 938 (&time, stardis, start, &compute_start, &compute_end, NULL)); 939 940 ERR(print_single_MC_result(estimator, stardis, stdout)); 941 } 942 943 /* Output random state? */ 944 WRITE_RANDOM_STATE(&stardis->rndgen_state_out_filename); 945 946 end: 947 if(estimator) SDIS(estimator_ref_put(estimator)); 948 if(green) SDIS(green_function_ref_put(green)); 949 if(args.rng_state) SSP(rng_ref_put(args.rng_state)); 950 return res; 951 error: 952 goto end; 953 } 954 955 static res_T 956 setup_solve_probe_boundary_args 957 (struct stardis* stardis, 958 const struct stardis_probe_boundary* probe, 959 struct sdis_solve_probe_boundary_args* args) 960 { 961 enum sdis_side probe_side = SDIS_SIDE_NULL__; 962 double uv[2] = {0, 0}; 963 size_t iprim = SIZE_MAX; 964 unsigned desc_ids[SG3D_PROP_TYPES_COUNT__]; 965 res_T res = RES_OK; 966 ASSERT(stardis && probe && args && (stardis->mode && MODE_COMPUTE_TEMP_MAP_ON_SURF)); 967 968 /* Calculate the probe position on the boundary */ 969 ERR(move_to_boundary(stardis, probe->position, probe->time[0], 1, &iprim, uv)); 970 971 ERR(sg3d_geometry_get_unique_triangle_properties(stardis->geometry.sg3d, 972 (unsigned)iprim, desc_ids)); 973 974 ERR(setup_probe_side(stardis, desc_ids, probe->side, iprim, &probe_side)); 975 ERR(setup_thermal_contact_resistance(stardis, desc_ids, probe_side)); 976 977 /* Setup of solve input parameters */ 978 args->nrealisations = stardis->samples; 979 args->iprim = iprim; 980 args->uv[0] = uv[0]; 981 args->uv[1] = uv[1]; 982 args->time_range[0] = probe->time[0]; 983 args->time_range[1] = probe->time[1]; 984 args->picard_order = stardis->picard_order; 985 args->side = probe_side; 986 args->register_paths = stardis->dump_paths; 987 args->diff_algo = stardis->diff_algo; 988 989 /* The solver does not accept that the side of the interface on which the 990 * probe is placed is invalid. Below, the side is arbitrarily defined because 991 * at this point, Stardis has already arbitrated that this side does not 992 * matter (i.e. there is no thermal contact resistance) */ 993 if(args->side == SDIS_SIDE_NULL__) { 994 args->side = SDIS_FRONT; 995 } 996 997 exit: 998 return res; 999 error: 1000 goto exit; 1001 } 1002 1003 static res_T 1004 compute_single_probe_on_interface 1005 (struct stardis* stardis, 1006 struct time* start, 1007 const struct stardis_probe_boundary* probe) 1008 { 1009 struct sdis_solve_probe_boundary_args args 1010 = SDIS_SOLVE_PROBE_BOUNDARY_ARGS_DEFAULT; 1011 res_T output_status[2] = {RES_OK, RES_OK}; 1012 res_T res = RES_OK; 1013 ASSERT(stardis && start && probe); 1014 1015 ERR(setup_solve_probe_boundary_args(stardis, probe, &args)); 1016 1017 /* Run the calculation */ 1018 if(stardis->mode & (MODE_GREEN_ASCII | MODE_GREEN_BIN)) { 1019 ERR(solve_green(stardis, start, &args)); 1020 } else { 1021 ERR(solve(stardis, start, &args, output_status)); 1022 } 1023 1024 exit: 1025 res = (res != RES_OK ? res : output_status[0]); 1026 res = (res != RES_OK ? res : output_status[1]); 1027 return res; 1028 error: 1029 goto exit; 1030 } 1031 1032 static res_T 1033 compute_multiple_probes_on_interface 1034 (struct stardis* stardis, 1035 struct time* start) 1036 { 1037 /* Probes */ 1038 const struct stardis_probe_boundary* probes = NULL; 1039 struct sdis_solve_probe_boundary_args* solve_args = NULL; 1040 struct sdis_solve_probe_boundary_list_args solve_list_args = 1041 SDIS_SOLVE_PROBE_BOUNDARY_LIST_ARGS_DEFAULT; 1042 size_t nprobes = 0; 1043 1044 /* Miscellaneous */ 1045 res_T output_status = RES_OK; 1046 res_T res = RES_OK; 1047 size_t i= 0; 1048 ASSERT(stardis && start); 1049 1050 /* Fetch the list of probes arguments */ 1051 probes = darray_probe_boundary_cdata_get(&stardis->probe_boundary_list); 1052 nprobes = darray_probe_boundary_size_get(&stardis->probe_boundary_list); 1053 ASSERT(nprobes > 1); 1054 1055 solve_args = MEM_CALLOC(stardis->allocator, nprobes, sizeof(*solve_args)); 1056 if(!probes) { 1057 logger_print(stardis->logger, LOG_ERROR, 1058 "Argument list allocation error for resolving multiple probes " 1059 "on the boundary.\n"); 1060 res = RES_MEM_ERR; 1061 goto error; 1062 } 1063 1064 /* Setup the solve arguments */ 1065 FOR_EACH(i, 0, nprobes) { 1066 solve_args[i] = SDIS_SOLVE_PROBE_BOUNDARY_ARGS_DEFAULT; 1067 ERR(setup_solve_probe_boundary_args(stardis, &probes[i], &solve_args[i])); 1068 } 1069 solve_list_args.probes = solve_args; 1070 solve_list_args.nprobes = nprobes; 1071 1072 /* Run calculations */ 1073 ERR(solve_list(stardis, start, &solve_list_args, &output_status)); 1074 1075 exit: 1076 if(probes) MEM_RM(stardis->allocator, solve_args); 1077 res = (res != RES_OK ? res : output_status); 1078 return res; 1079 error: 1080 goto exit; 1081 } 1082 1083 static res_T 1084 compute_multiple_flux_probes_on_interface 1085 (struct stardis* stardis, 1086 struct time* start) 1087 { 1088 /* Probes */ 1089 const struct stardis_probe_boundary* probes = NULL; 1090 size_t nprobes = 0; 1091 res_T res = RES_OK; 1092 ASSERT(stardis && start); 1093 1094 /* Fetch the list of probes arguments */ 1095 probes = darray_probe_boundary_cdata_get(&stardis->probe_boundary_list); 1096 nprobes = darray_probe_boundary_size_get(&stardis->probe_boundary_list); 1097 ASSERT(nprobes > 1); 1098 1099 /* Run calculations */ 1100 ERR(solve_flux_list(stardis, start, probes, nprobes)); 1101 1102 exit: 1103 return res; 1104 error: 1105 goto exit; 1106 } 1107 1108 /******************************************************************************* 1109 * Local functions 1110 ******************************************************************************/ 1111 res_T 1112 compute_probe_on_interface(struct stardis* stardis, struct time* start) 1113 { 1114 int temp_flags = MODE_COMPUTE_PROBE_TEMP_ON_SURF 1115 | MODE_COMPUTE_LIST_PROBE_TEMP_ON_SURF; 1116 int flux_flags = MODE_COMPUTE_PROBE_FLUX_DNSTY_ON_SURF 1117 | MODE_COMPUTE_LIST_PROBE_FLUX_DNSTY_ON_SURF; 1118 res_T res = RES_OK; 1119 ASSERT(stardis && start); 1120 ASSERT((stardis->mode & temp_flags) != (stardis->mode & flux_flags)); /* xor */ 1121 1122 if(stardis->mode & temp_flags) { 1123 if(darray_probe_boundary_size_get(&stardis->probe_boundary_list) > 1) { 1124 res = compute_multiple_probes_on_interface(stardis, start); 1125 if(res != RES_OK) goto error; 1126 } else { 1127 const struct stardis_probe_boundary* probe 1128 = darray_probe_boundary_cdata_get(&stardis->probe_boundary_list); 1129 res = compute_single_probe_on_interface(stardis, start, probe); 1130 if(res != RES_OK) goto error; 1131 } 1132 } else if(stardis->mode & flux_flags) { 1133 if(darray_probe_boundary_size_get(&stardis->probe_boundary_list) > 1) { 1134 res = compute_multiple_flux_probes_on_interface(stardis, start); 1135 if(res != RES_OK) goto error; 1136 } else { 1137 const struct stardis_probe_boundary* probe 1138 = darray_probe_boundary_cdata_get(&stardis->probe_boundary_list); 1139 res = compute_single_flux_probe_on_interface(stardis, start, probe); 1140 if(res != RES_OK) goto error; 1141 } 1142 } 1143 1144 exit: 1145 return res; 1146 error: 1147 goto exit; 1148 }