sdis_green.c (54499B)
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_device_c.h" 17 #include "sdis_estimator_c.h" 18 #include "sdis_green.h" 19 #include "sdis_interface_c.h" 20 #include "sdis_log.h" 21 #include "sdis_medium_c.h" 22 #include "sdis_misc.h" 23 #include "sdis_radiative_env_c.h" 24 #include "sdis_scene_c.h" 25 #include "sdis_source_c.h" 26 27 #include <star/ssp.h> 28 29 #include <rsys/cstr.h> 30 #include <rsys/dynamic_array.h> 31 #include <rsys/hash_table.h> 32 #include <rsys/mem_allocator.h> 33 #include <rsys/ref_count.h> 34 #include <rsys/dynamic_array.h> 35 36 #include <limits.h> 37 38 /* Path used to estimate the green function */ 39 struct sdis_green_path { 40 /* Internal data. Should not be accessed */ 41 void* green__; 42 size_t id__; 43 }; 44 45 struct power_term { 46 double term; /* Power term computed during green estimation */ 47 unsigned id; /* Identifier of the medium of the term */ 48 }; 49 50 #define POWER_TERM_NULL__ {DBL_MAX, UINT_MAX} 51 static const struct power_term POWER_TERM_NULL = POWER_TERM_NULL__; 52 53 static INLINE void 54 power_term_init(struct mem_allocator* allocator, struct power_term* term) 55 { 56 ASSERT(term); (void)allocator; 57 *term = POWER_TERM_NULL; 58 } 59 60 /* Generate the dynamic array of power terms */ 61 #define DARRAY_NAME power_term 62 #define DARRAY_DATA struct power_term 63 #define DARRAY_FUNCTOR_INIT power_term_init 64 #include <rsys/dynamic_array.h> 65 66 struct flux_term { 67 double term; 68 unsigned id; /* Id of the interface of the flux term */ 69 enum sdis_side side; 70 }; 71 #define FLUX_TERM_NULL__ {DBL_MAX, UINT_MAX, SDIS_SIDE_NULL__} 72 static const struct flux_term FLUX_TERM_NULL = FLUX_TERM_NULL__; 73 74 static INLINE void 75 flux_term_init(struct mem_allocator* allocator, struct flux_term* term) 76 { 77 ASSERT(term); (void)allocator; 78 *term = FLUX_TERM_NULL; 79 } 80 81 /* Generate the dynamic array of flux terms */ 82 #define DARRAY_NAME flux_term 83 #define DARRAY_DATA struct flux_term 84 #define DARRAY_FUNCTOR_INIT flux_term_init 85 #include <rsys/dynamic_array.h> 86 87 static INLINE void 88 extflux_terms_init 89 (struct mem_allocator* allocator, 90 struct sdis_green_external_flux_terms* terms) 91 { 92 ASSERT(terms); (void)allocator; 93 *terms = SDIS_GREEN_EXTERNAL_FLUX_TERMS_NULL; 94 } 95 96 /* Generate the dynamic array of the external flux terms */ 97 #define DARRAY_NAME extflux_terms 98 #define DARRAY_DATA struct sdis_green_external_flux_terms 99 #define DARRAY_FUNCTOR_INIT extflux_terms_init 100 #include <rsys/dynamic_array.h> 101 102 struct green_path { 103 double elapsed_time; 104 struct darray_extflux_terms extflux_terms; /* List of external flux terms */ 105 struct darray_flux_term flux_terms; /* List of flux terms */ 106 struct darray_power_term power_terms; /* List of volumic power terms */ 107 union { 108 struct sdis_rwalk_vertex vertex; 109 struct sdis_interface_fragment fragment; 110 struct sdis_radiative_ray ray; 111 } limit; 112 unsigned limit_id; /* Identifier of the limit medium/interface */ 113 enum sdis_green_path_end_type end_type; 114 115 /* Indices of the last accessed medium/interface. Used to speed up the access 116 * to the medium/interface. */ 117 uint16_t ilast_medium; 118 uint16_t ilast_interf; 119 }; 120 121 static INLINE void 122 green_path_init(struct mem_allocator* allocator, struct green_path* path) 123 { 124 ASSERT(path); 125 darray_extflux_terms_init(allocator, &path->extflux_terms); 126 darray_flux_term_init(allocator, &path->flux_terms); 127 darray_power_term_init(allocator, &path->power_terms); 128 path->elapsed_time = -INF; 129 path->limit.vertex = SDIS_RWALK_VERTEX_NULL; 130 path->limit.fragment = SDIS_INTERFACE_FRAGMENT_NULL; 131 path->limit.ray = SDIS_RADIATIVE_RAY_NULL; 132 path->limit_id = UINT_MAX; 133 path->end_type = SDIS_GREEN_PATH_END_TYPES_COUNT__; 134 path->ilast_medium = UINT16_MAX; 135 path->ilast_interf = UINT16_MAX; 136 } 137 138 static INLINE void 139 green_path_release(struct green_path* path) 140 { 141 ASSERT(path); 142 darray_flux_term_release(&path->flux_terms); 143 darray_power_term_release(&path->power_terms); 144 darray_extflux_terms_release(&path->extflux_terms); 145 } 146 147 static INLINE res_T 148 green_path_copy(struct green_path* dst, const struct green_path* src) 149 { 150 res_T res = RES_OK; 151 ASSERT(dst && src); 152 dst->elapsed_time = src->elapsed_time; 153 dst->limit = src->limit; 154 dst->limit_id = src->limit_id; 155 dst->end_type = src->end_type; 156 dst->ilast_medium = src->ilast_medium; 157 dst->ilast_interf = src->ilast_interf; 158 res = darray_extflux_terms_copy(&dst->extflux_terms, &src->extflux_terms); 159 if(res != RES_OK) return res; 160 res = darray_flux_term_copy(&dst->flux_terms, &src->flux_terms); 161 if(res != RES_OK) return res; 162 res = darray_power_term_copy(&dst->power_terms, &src->power_terms); 163 if(res != RES_OK) return res; 164 return RES_OK; 165 } 166 167 static INLINE res_T 168 green_path_copy_and_clear(struct green_path* dst, struct green_path* src) 169 { 170 res_T res = RES_OK; 171 ASSERT(dst && src); 172 dst->elapsed_time = src->elapsed_time; 173 dst->limit = src->limit; 174 dst->limit_id = src->limit_id; 175 dst->end_type = src->end_type; 176 dst->ilast_medium = src->ilast_medium; 177 dst->ilast_interf = src->ilast_interf; 178 res = darray_extflux_terms_copy_and_clear 179 (&dst->extflux_terms, &src->extflux_terms); 180 if(res != RES_OK) return res; 181 res = darray_flux_term_copy_and_clear(&dst->flux_terms, &src->flux_terms); 182 if(res != RES_OK) return res; 183 res = darray_power_term_copy_and_clear(&dst->power_terms, &src->power_terms); 184 if(res != RES_OK) return res; 185 return RES_OK; 186 187 } 188 189 static INLINE res_T 190 green_path_copy_and_release(struct green_path* dst, struct green_path* src) 191 { 192 res_T res = RES_OK; 193 ASSERT(dst && src); 194 dst->elapsed_time = src->elapsed_time; 195 dst->limit = src->limit; 196 dst->limit_id = src->limit_id; 197 dst->end_type = src->end_type; 198 dst->ilast_medium = src->ilast_medium; 199 dst->ilast_interf = src->ilast_interf; 200 res = darray_extflux_terms_copy_and_release 201 (&dst->extflux_terms, &src->extflux_terms); 202 if(res != RES_OK) return res; 203 res = darray_flux_term_copy_and_release(&dst->flux_terms, &src->flux_terms); 204 if(res != RES_OK) return res; 205 res = darray_power_term_copy_and_release(&dst->power_terms, &src->power_terms); 206 if(res != RES_OK) return res; 207 return RES_OK; 208 } 209 210 static res_T 211 green_path_write(const struct green_path* path, FILE* stream) 212 { 213 size_t sz = 0; 214 res_T res = RES_OK; 215 ASSERT(path && stream); 216 217 #define WRITE(Var, N) { \ 218 if(fwrite((Var), sizeof(*(Var)), (N), stream) != (N)) { \ 219 res = RES_IO_ERR; \ 220 goto error; \ 221 } \ 222 } (void)0 223 224 /* Write elapsed time */ 225 WRITE(&path->elapsed_time, 1); 226 227 /* Write the list of external flux terms */ 228 sz = darray_extflux_terms_size_get(&path->extflux_terms); 229 WRITE(&sz, 1); 230 WRITE(darray_extflux_terms_cdata_get(&path->extflux_terms), sz); 231 232 /* Write the list of flux terms */ 233 sz = darray_flux_term_size_get(&path->flux_terms); 234 WRITE(&sz, 1); 235 WRITE(darray_flux_term_cdata_get(&path->flux_terms), sz); 236 237 /* Write the list of power terms */ 238 sz = darray_power_term_size_get(&path->power_terms); 239 WRITE(&sz, 1); 240 WRITE(darray_power_term_cdata_get(&path->power_terms), sz); 241 242 /* Write the limit point */ 243 WRITE(&path->limit, 1); 244 WRITE(&path->limit_id, 1); 245 WRITE(&path->end_type, 1); 246 247 /* Write miscellaneous data */ 248 WRITE(&path->ilast_medium, 1); 249 WRITE(&path->ilast_interf, 1); 250 251 #undef WRITE 252 253 exit: 254 return res; 255 error: 256 goto exit; 257 } 258 259 static res_T 260 green_path_read(struct green_path* path, FILE* stream) 261 { 262 size_t sz = 0; 263 res_T res = RES_OK; 264 ASSERT(path && stream); 265 266 #define READ(Var, N) { \ 267 if(fread((Var), sizeof(*(Var)), (N), stream) != (N)) { \ 268 if(feof(stream)) { \ 269 res = RES_BAD_ARG; \ 270 } else if(ferror(stream)) { \ 271 res = RES_IO_ERR; \ 272 } else { \ 273 res = RES_UNKNOWN_ERR; \ 274 } \ 275 goto error; \ 276 } \ 277 } (void)0 278 279 /* Read elapsed time */ 280 READ(&path->elapsed_time, 1); 281 282 /* Read the list of external flux terms */ 283 READ(&sz, 1); 284 res = darray_extflux_terms_resize(&path->extflux_terms, sz); 285 if(res != RES_OK) goto error; 286 READ(darray_extflux_terms_data_get(&path->extflux_terms), sz); 287 288 /* Read the list of flux terms */ 289 READ(&sz, 1); 290 res = darray_flux_term_resize(&path->flux_terms, sz); 291 if(res != RES_OK) goto error; 292 READ(darray_flux_term_data_get(&path->flux_terms), sz); 293 294 /* Read the list of power tems */ 295 READ(&sz, 1); 296 res = darray_power_term_resize(&path->power_terms, sz); 297 if(res != RES_OK) goto error; 298 READ(darray_power_term_data_get(&path->power_terms), sz); 299 300 /* Read the limit point */ 301 READ(&path->limit, 1); 302 READ(&path->limit_id, 1); 303 READ(&path->end_type, 1); 304 305 /* Read the miscellaneous data */ 306 READ(&path->ilast_medium, 1); 307 READ(&path->ilast_interf, 1); 308 309 #undef READ 310 311 exit: 312 return res; 313 error: 314 goto exit; 315 } 316 317 /* Generate the dynamic array of green paths */ 318 #define DARRAY_NAME green_path 319 #define DARRAY_DATA struct green_path 320 #define DARRAY_FUNCTOR_INIT green_path_init 321 #define DARRAY_FUNCTOR_RELEASE green_path_release 322 #define DARRAY_FUNCTOR_COPY green_path_copy 323 #define DARRAY_FUNCTOR_COPY_AND_RELEASE green_path_copy_and_release 324 #include <rsys/dynamic_array.h> 325 326 /* Generate the hash table that maps and id to an interface */ 327 #define HTABLE_NAME interf 328 #define HTABLE_KEY unsigned 329 #define HTABLE_DATA struct sdis_interface* 330 #include <rsys/hash_table.h> 331 332 /* Generate the hash table that maps an id to a medium */ 333 #define HTABLE_NAME medium 334 #define HTABLE_KEY unsigned 335 #define HTABLE_DATA struct sdis_medium* 336 #include <rsys/hash_table.h> 337 338 struct sdis_green_function { 339 struct htable_medium media; 340 struct htable_interf interfaces; 341 struct darray_green_path paths; /* List of paths used to estimate the green */ 342 343 size_t npaths_valid; 344 size_t npaths_invalid; 345 hash256_T signature; 346 347 struct accum realisation_time; /* Time per realisation */ 348 349 enum ssp_rng_type rng_type; 350 FILE* rng_state; 351 352 ref_T ref; 353 struct sdis_scene* scn; 354 }; 355 356 /******************************************************************************* 357 * Helper functions 358 ******************************************************************************/ 359 static INLINE res_T 360 check_sdis_green_function_create_from_stream_args 361 (const struct sdis_green_function_create_from_stream_args* args) 362 { 363 if(!args || !args->scene || !args->stream) 364 return RES_BAD_ARG; 365 return RES_OK; 366 } 367 368 static res_T 369 ensure_medium_registration 370 (struct sdis_green_function* green, 371 struct sdis_medium* mdm) 372 { 373 unsigned id; 374 res_T res = RES_OK; 375 ASSERT(green && mdm); 376 377 id = medium_get_id(mdm); 378 if(htable_medium_find(&green->media, &id)) goto exit; 379 380 res = htable_medium_set(&green->media, &id, &mdm); 381 if(res != RES_OK) goto error; 382 383 SDIS(medium_ref_get(mdm)); 384 385 exit: 386 return res; 387 error: 388 goto exit; 389 } 390 391 static res_T 392 ensure_interface_registration 393 (struct sdis_green_function* green, 394 struct sdis_interface* interf) 395 { 396 unsigned id; 397 res_T res = RES_OK; 398 ASSERT(green && interf); 399 400 id = interface_get_id(interf); 401 if(htable_interf_find(&green->interfaces, &id)) goto exit; 402 403 res = htable_interf_set(&green->interfaces, &id, &interf); 404 if(res != RES_OK) goto error; 405 406 SDIS(interface_ref_get(interf)); 407 408 exit: 409 return res; 410 error: 411 goto exit; 412 } 413 414 static FINLINE struct sdis_medium* 415 green_function_fetch_medium 416 (struct sdis_green_function* green, const unsigned medium_id) 417 { 418 struct sdis_medium* const* pmdm = NULL; 419 ASSERT(green); 420 pmdm = htable_medium_find(&green->media, &medium_id); 421 ASSERT(pmdm); 422 return *pmdm; 423 } 424 425 static FINLINE struct sdis_interface* 426 green_function_fetch_interf 427 (struct sdis_green_function* green, const unsigned interf_id) 428 { 429 struct sdis_interface* const* pinterf = NULL; 430 ASSERT(green); 431 pinterf = htable_interf_find(&green->interfaces, &interf_id); 432 ASSERT(pinterf); 433 return *pinterf; 434 } 435 436 static double /* [K] */ 437 green_path_power_contribution 438 (struct sdis_green_function* green, 439 const struct green_path* path) 440 { 441 double temperature = 0; /* [K] */ 442 size_t i=0, n=0; 443 444 ASSERT(green && path); 445 446 n = darray_power_term_size_get(&path->power_terms); 447 FOR_EACH(i, 0, n) { 448 struct sdis_rwalk_vertex vtx = SDIS_RWALK_VERTEX_NULL; 449 const struct power_term* power_term = NULL; 450 const struct sdis_medium* medium = NULL; 451 double power = 0; /* [W] */ 452 453 power_term = darray_power_term_cdata_get(&path->power_terms) + i; 454 medium = green_function_fetch_medium(green, power_term->id); 455 456 /* Dummy argument used only to satisfy the function profile used to recover 457 * power. Its position is unused, since power is assumed to be constant in 458 * space, and its time is set to infinity, since the green function is 459 * assumed to be evaluated at steady state */ 460 vtx.time = INF; 461 power = solid_get_volumic_power(medium, &vtx); 462 463 temperature += power_term->term * power; /* [K] */ 464 } 465 return temperature; /* [K] */ 466 } 467 468 static double /* [K] */ 469 green_path_flux_contribution 470 (struct sdis_green_function* green, 471 const struct green_path* path) 472 { 473 double temperature = 0; 474 size_t i=0, n=0; 475 476 ASSERT(green && path); 477 478 n = darray_flux_term_size_get(&path->flux_terms); 479 FOR_EACH(i, 0, n) { 480 struct sdis_interface_fragment frag = SDIS_INTERFACE_FRAGMENT_NULL; 481 const struct flux_term* flux_term = NULL; 482 const struct sdis_interface* interf = NULL; 483 double flux = 0; /* [W/m^2] */ 484 485 flux_term = darray_flux_term_cdata_get(&path->flux_terms) + i; 486 interf = green_function_fetch_interf(green, flux_term->id); 487 488 /* Interface fragment. Its position is unused, since flux is assumed to be 489 * constant in space, and its time is set to infinity, since the green 490 * function is assumed to be evaluated at steady state */ 491 frag.time = INF; 492 frag.side = flux_term->side; 493 flux = interface_side_get_flux(interf, &frag); 494 495 temperature += flux_term->term * flux; /* [K] */ 496 } 497 return temperature; /* [K] */ 498 } 499 500 static double /* [K] */ 501 green_path_external_flux_contribution 502 (struct sdis_green_function* green, 503 const struct green_path* path) 504 { 505 const struct sdis_source* extsrc = NULL; 506 double value = 0; 507 size_t i=0, n=0; 508 509 ASSERT(green && path); 510 511 if((extsrc = green->scn->source) == NULL) return 0; 512 513 n = darray_extflux_terms_size_get(&path->extflux_terms); 514 FOR_EACH(i, 0, n) { 515 const struct sdis_green_external_flux_terms* extflux = NULL; 516 double power = 0; /* [W] */ 517 double diffrad = 0; /* [W/m^2/sr] */ 518 519 extflux = darray_extflux_terms_cdata_get(&path->extflux_terms)+i; 520 power = source_get_power(extsrc, extflux->time); 521 diffrad = source_get_diffuse_radiance(extsrc, extflux->time, extflux->dir); 522 523 value += extflux->term_wrt_power * power; /* [K] */ 524 value += extflux->term_wrt_diffuse_radiance * diffrad; /* [K] */ 525 } 526 return value; /* [K] */ 527 } 528 529 static res_T 530 green_function_solve_path 531 (struct sdis_green_function* green, 532 const size_t ipath, 533 double* weight) 534 { 535 const struct green_path* path = NULL; 536 const struct sdis_medium* medium = NULL; 537 const struct sdis_interface* interf = NULL; 538 struct sdis_scene* scn = NULL; 539 struct sdis_rwalk_vertex vtx = SDIS_RWALK_VERTEX_NULL; 540 struct sdis_interface_fragment frag = SDIS_INTERFACE_FRAGMENT_NULL; 541 struct sdis_radiative_ray ray = SDIS_RADIATIVE_RAY_NULL; 542 double power; 543 double flux; 544 double external_flux; 545 double end_temperature; 546 res_T res = RES_OK; 547 ASSERT(green && ipath < darray_green_path_size_get(&green->paths) && weight); 548 549 path = darray_green_path_cdata_get(&green->paths) + ipath; 550 if(path->end_type == SDIS_GREEN_PATH_END_ERROR) { /* Rejected path */ 551 res = RES_BAD_OP; 552 goto error; 553 } 554 555 power = green_path_power_contribution(green, path); 556 flux = green_path_flux_contribution(green, path); 557 external_flux = green_path_external_flux_contribution(green, path); 558 559 /* Compute path's end temperature */ 560 switch(path->end_type) { 561 case SDIS_GREEN_PATH_END_AT_INTERFACE: 562 interf = green_function_fetch_interf(green, path->limit_id); 563 frag = path->limit.fragment; 564 end_temperature = interface_side_get_temperature(interf, &frag); 565 break; 566 case SDIS_GREEN_PATH_END_AT_RADIATIVE_ENV: 567 SDIS(green_function_get_scene(green, &scn)); 568 ray = path->limit.ray; 569 end_temperature = radiative_env_get_temperature(green->scn->radenv, &ray); 570 break; 571 case SDIS_GREEN_PATH_END_IN_VOLUME: 572 medium = green_function_fetch_medium(green, path->limit_id); 573 vtx = path->limit.vertex; 574 end_temperature = medium_get_temperature(medium, &vtx); 575 break; 576 default: FATAL("Unreachable code.\n"); break; 577 } 578 579 if(SDIS_TEMPERATURE_IS_UNKNOWN(end_temperature)) { 580 log_err(green->scn->dev, 581 "%s: unknown boundary/initial condition.\n", FUNC_NAME); 582 res = RES_BAD_ARG; 583 goto error; 584 } 585 586 /* Compute the path weight */ 587 *weight = power + flux + external_flux + end_temperature; 588 589 exit: 590 return res; 591 error: 592 goto exit; 593 } 594 595 static res_T 596 write_media(struct sdis_green_function* green, FILE* stream) 597 { 598 struct htable_medium_iterator it, it_end; 599 size_t nmedia = 0; 600 res_T res = RES_OK; 601 ASSERT(green && stream); 602 603 #define WRITE(Var) { \ 604 if(fwrite((Var), sizeof(*(Var)), 1, stream) != 1) { \ 605 res = RES_IO_ERR; \ 606 goto error; \ 607 } \ 608 } (void)0 609 610 nmedia = htable_medium_size_get(&green->media); 611 WRITE(&nmedia); 612 613 htable_medium_begin(&green->media, &it); 614 htable_medium_end(&green->media, &it_end); 615 while(!htable_medium_iterator_eq(&it, &it_end)) { 616 const struct sdis_medium* mdm = *htable_medium_iterator_data_get(&it); 617 htable_medium_iterator_next(&it); 618 WRITE(&mdm->id); 619 WRITE(&mdm->type); 620 } 621 622 #undef WRITE 623 624 exit: 625 return res; 626 error: 627 goto exit; 628 } 629 630 static res_T 631 read_media(struct sdis_green_function* green, FILE* stream) 632 { 633 size_t nmedia = 0; 634 size_t imedium = 0; 635 res_T res = RES_OK; 636 ASSERT(green && stream); 637 638 #define READ(Var) { \ 639 if(fread((Var), sizeof(*(Var)), 1, stream) != 1) { \ 640 if(feof(stream)) { \ 641 res = RES_BAD_ARG; \ 642 } else if(ferror(stream)) { \ 643 res = RES_IO_ERR; \ 644 } else { \ 645 res = RES_UNKNOWN_ERR; \ 646 } \ 647 goto error; \ 648 } \ 649 } (void)0 650 651 READ(&nmedia); 652 FOR_EACH(imedium, 0, nmedia) { 653 struct name* name = NULL; 654 struct sdis_medium* mdm = NULL; 655 struct fid id; 656 enum sdis_medium_type mdm_type; 657 658 READ(&id); 659 READ(&mdm_type); 660 661 name = flist_name_get(&green->scn->dev->media_names, id); 662 if(!name) { 663 log_err(green->scn->dev, "%s: a Stardis medium is missing.\n", 664 FUNC_NAME); 665 res = RES_BAD_ARG; 666 goto error; 667 } 668 669 mdm = name->mem; 670 671 if(mdm_type != mdm->type) { 672 log_err(green->scn->dev, "%s: inconsistency between the a Stardis medium " 673 "and its serialised data.\n", FUNC_NAME); 674 res = RES_BAD_ARG; 675 goto error; 676 } 677 678 res = ensure_medium_registration(green, mdm); 679 if(res != RES_OK) goto error; 680 } 681 682 #undef READ 683 684 exit: 685 return res; 686 error: 687 goto exit; 688 } 689 690 static res_T 691 write_interfaces(struct sdis_green_function* green, FILE* stream) 692 { 693 struct htable_interf_iterator it, it_end; 694 size_t ninterfaces = 0; 695 res_T res = RES_OK; 696 ASSERT(green && stream); 697 698 #define WRITE(Var) { \ 699 if(fwrite((Var), sizeof(*(Var)), 1, stream) != 1) { \ 700 res = RES_IO_ERR; \ 701 goto error; \ 702 } \ 703 } (void)0 704 ninterfaces = htable_interf_size_get(&green->interfaces); 705 WRITE(&ninterfaces); 706 707 htable_interf_begin(&green->interfaces, &it); 708 htable_interf_end(&green->interfaces, &it_end); 709 while(!htable_interf_iterator_eq(&it, &it_end)) { 710 const struct sdis_interface* interf = *htable_interf_iterator_data_get(&it); 711 htable_interf_iterator_next(&it); 712 WRITE(&interf->id); 713 WRITE(&interf->medium_front->id); 714 WRITE(&interf->medium_front->type); 715 WRITE(&interf->medium_back->id); 716 WRITE(&interf->medium_back->type); 717 } 718 #undef WRITE 719 720 exit: 721 return res; 722 error: 723 goto exit; 724 } 725 726 static res_T 727 read_interfaces(struct sdis_green_function* green, FILE* stream) 728 { 729 size_t ninterfs = 0; 730 size_t iinterf = 0; 731 res_T res = RES_OK; 732 ASSERT(green && stream); 733 734 #define READ(Var) { \ 735 if(fread((Var), sizeof(*(Var)), 1, stream) != 1) { \ 736 if(feof(stream)) { \ 737 res = RES_BAD_ARG; \ 738 } else if(ferror(stream)) { \ 739 res = RES_IO_ERR; \ 740 } else { \ 741 res = RES_UNKNOWN_ERR; \ 742 } \ 743 goto error; \ 744 } \ 745 } (void)0 746 747 READ(&ninterfs); 748 FOR_EACH(iinterf, 0, ninterfs) { 749 struct name* name = NULL; 750 struct sdis_interface* interf = NULL; 751 struct sdis_medium* mdm_front = NULL; 752 struct sdis_medium* mdm_back = NULL; 753 struct fid id; 754 struct fid mdm_front_id; 755 struct fid mdm_back_id; 756 enum sdis_medium_type mdm_front_type; 757 enum sdis_medium_type mdm_back_type; 758 759 READ(&id); 760 READ(&mdm_front_id); 761 READ(&mdm_front_type); 762 READ(&mdm_back_id); 763 READ(&mdm_back_type); 764 765 name = flist_name_get(&green->scn->dev->interfaces_names, id); 766 if(!name) { 767 log_err(green->scn->dev, "%s: a Stardis interface is missing.\n", 768 FUNC_NAME); 769 res = RES_BAD_ARG; 770 goto error; 771 } 772 773 interf = name->mem; 774 mdm_front = flist_name_get(&green->scn->dev->media_names, mdm_front_id)->mem; 775 mdm_back = flist_name_get(&green->scn->dev->media_names, mdm_back_id)->mem; 776 777 if(mdm_front != interf->medium_front 778 || mdm_back != interf->medium_back 779 || mdm_front_type != interf->medium_front->type 780 || mdm_back_type != interf->medium_back->type) { 781 log_err(green->scn->dev, "%s: inconsistency between the a Stardis interface " 782 "and its serialised data.\n", FUNC_NAME); 783 res = RES_BAD_ARG; 784 goto error; 785 } 786 787 res = ensure_interface_registration(green, interf); 788 if(res != RES_OK) goto error; 789 } 790 791 #undef READ 792 793 exit: 794 return res; 795 error: 796 goto exit; 797 } 798 799 static res_T 800 write_paths_list(struct sdis_green_function* green, FILE* stream) 801 { 802 size_t npaths = 0; 803 size_t ipath = 0; 804 res_T res = RES_OK; 805 ASSERT(green && stream); 806 807 #define WRITE(Var) { \ 808 if(fwrite((Var), sizeof(*(Var)), 1, stream) != 1) { \ 809 res = RES_IO_ERR; \ 810 goto error; \ 811 } \ 812 } (void)0 813 npaths = darray_green_path_size_get(&green->paths); 814 WRITE(&npaths); 815 FOR_EACH(ipath, 0, npaths) { 816 const struct green_path* path = NULL; 817 path = darray_green_path_cdata_get(&green->paths) + ipath; 818 819 res = green_path_write(path, stream); 820 if(res != RES_OK) goto error; 821 } 822 #undef WRITE 823 824 exit: 825 return res; 826 error: 827 goto exit; 828 } 829 830 static res_T 831 read_paths_list(struct sdis_green_function* green, FILE* stream) 832 { 833 size_t npaths = 0; 834 size_t ipath = 0; 835 res_T res = RES_OK; 836 837 #define READ(Var) { \ 838 if(fread((Var), sizeof(*(Var)), 1, stream) != 1) { \ 839 if(feof(stream)) { \ 840 res = RES_BAD_ARG; \ 841 } else if(ferror(stream)) { \ 842 res = RES_IO_ERR; \ 843 } else { \ 844 res = RES_UNKNOWN_ERR; \ 845 } \ 846 goto error; \ 847 } \ 848 } (void)0 849 850 READ(&npaths); 851 res = darray_green_path_resize(&green->paths, npaths); 852 if(res != RES_OK) goto error; 853 854 FOR_EACH(ipath, 0, npaths) { 855 struct green_path* path = NULL; 856 path = darray_green_path_data_get(&green->paths) + ipath; 857 858 res = green_path_read(path, stream); 859 if(res != RES_OK) goto error; 860 } 861 #undef READ 862 863 exit: 864 return res; 865 error: 866 goto exit; 867 } 868 869 static void 870 green_function_clear(struct sdis_green_function* green) 871 { 872 struct htable_medium_iterator it_medium, end_medium; 873 struct htable_interf_iterator it_interf, end_interf; 874 ASSERT(green); 875 876 /* Clean up medium hash table */ 877 htable_medium_begin(&green->media, &it_medium); 878 htable_medium_end(&green->media, &end_medium); 879 while(!htable_medium_iterator_eq(&it_medium, &end_medium)) { 880 struct sdis_medium* medium; 881 medium = *htable_medium_iterator_data_get(&it_medium); 882 SDIS(medium_ref_put(medium)); 883 htable_medium_iterator_next(&it_medium); 884 } 885 htable_medium_clear(&green->media); 886 887 /* Clean up the interface hash table */ 888 htable_interf_begin(&green->interfaces, &it_interf); 889 htable_interf_end(&green->interfaces, &end_interf); 890 while(!htable_interf_iterator_eq(&it_interf, &end_interf)) { 891 struct sdis_interface* interf; 892 interf = *htable_interf_iterator_data_get(&it_interf); 893 SDIS(interface_ref_put(interf)); 894 htable_interf_iterator_next(&it_interf); 895 } 896 htable_interf_clear(&green->interfaces); 897 898 /* Clean up the registered paths */ 899 darray_green_path_clear(&green->paths); 900 } 901 902 static void 903 green_function_release(ref_T* ref) 904 { 905 struct sdis_scene* scn; 906 struct sdis_green_function* green; 907 ASSERT(ref); 908 green = CONTAINER_OF(ref, struct sdis_green_function, ref); 909 scn = green->scn; 910 green_function_clear(green); 911 htable_medium_release(&green->media); 912 htable_interf_release(&green->interfaces); 913 darray_green_path_release(&green->paths); 914 if(green->rng_state) fclose(green->rng_state); 915 MEM_RM(scn->dev->allocator, green); 916 SDIS(scene_ref_put(scn)); 917 } 918 919 /******************************************************************************* 920 * Exported functions 921 ******************************************************************************/ 922 res_T 923 sdis_green_function_ref_get(struct sdis_green_function* green) 924 { 925 if(!green) return RES_BAD_ARG; 926 ref_get(&green->ref); 927 return RES_OK; 928 } 929 930 res_T 931 sdis_green_function_ref_put(struct sdis_green_function* green) 932 { 933 if(!green) return RES_BAD_ARG; 934 ref_put(&green->ref, green_function_release); 935 return RES_OK; 936 } 937 938 res_T 939 sdis_green_function_solve 940 (struct sdis_green_function* green, 941 struct sdis_estimator** out_estimator) 942 { 943 struct sdis_estimator* estimator = NULL; 944 size_t npaths; 945 size_t ipath; 946 size_t N = 0; /* #realisations */ 947 double accum = 0; 948 double accum2 = 0; 949 res_T res = RES_OK; 950 951 if(!green || !out_estimator) { 952 res = RES_BAD_ARG; 953 goto error; 954 } 955 956 npaths = darray_green_path_size_get(&green->paths); 957 958 /* Create the estimator */ 959 res = estimator_create(green->scn->dev, SDIS_ESTIMATOR_TEMPERATURE, &estimator); 960 if(res != RES_OK) goto error; 961 962 /* Solve the green function */ 963 FOR_EACH(ipath, 0, npaths) { 964 double w; 965 966 res = green_function_solve_path(green, ipath, &w); 967 if(res == RES_BAD_OP) continue; 968 if(res != RES_OK) goto error; 969 970 accum += w; 971 accum2 += w*w; 972 ++N; 973 } 974 975 /* Setup the estimated temperature */ 976 estimator_setup_realisations_count(estimator, npaths, N); 977 estimator_setup_temperature(estimator, accum, accum2); 978 estimator_setup_realisation_time 979 (estimator, green->realisation_time.sum, green->realisation_time.sum2); 980 981 exit: 982 if(out_estimator) *out_estimator = estimator; 983 return res; 984 error: 985 if(estimator) { 986 SDIS(estimator_ref_put(estimator)); 987 estimator = NULL; 988 } 989 goto exit; 990 } 991 992 res_T 993 sdis_green_function_write(struct sdis_green_function* green, FILE* stream) 994 { 995 struct ssp_rng* rng = NULL; 996 hash256_T hash; 997 res_T res = RES_OK; 998 999 if(!green || !stream) { 1000 res = RES_BAD_ARG; 1001 goto error; 1002 } 1003 1004 #define WRITE(Var, Nb) { \ 1005 if(fwrite((Var), sizeof(*(Var)), (Nb), stream) != (Nb)) { \ 1006 res = RES_IO_ERR; \ 1007 goto error; \ 1008 } \ 1009 } (void)0 1010 1011 if(green->rng_type == SSP_RNG_TYPE_NULL) { 1012 log_err(green->scn->dev, 1013 "%s: could not write a green function with an unknown RNG type.\n", 1014 FUNC_NAME); 1015 res = RES_BAD_ARG; 1016 goto error; 1017 } 1018 1019 WRITE(&SDIS_GREEN_FUNCTION_VERSION, 1); 1020 1021 res = scene_compute_hash(green->scn, hash); 1022 if(res != RES_OK) goto error; 1023 1024 WRITE(hash, sizeof(hash256_T)); 1025 WRITE(green->signature, sizeof(hash256_T)); 1026 1027 res = write_media(green, stream); 1028 if(res != RES_OK) goto error; 1029 res = write_interfaces(green, stream); 1030 if(res != RES_OK) goto error; 1031 res = write_paths_list(green, stream); 1032 if(res != RES_OK) goto error; 1033 1034 WRITE(&green->npaths_valid, 1); 1035 WRITE(&green->npaths_invalid, 1); 1036 WRITE(&green->realisation_time, 1); 1037 WRITE(&green->rng_type, 1); 1038 #undef WRITE 1039 1040 /* Create a temporary RNG used to serialise the RNG state */ 1041 res = ssp_rng_create(green->scn->dev->allocator, green->rng_type, &rng); 1042 if(res != RES_OK) goto error; 1043 rewind(green->rng_state); 1044 res = ssp_rng_read(rng, green->rng_state); 1045 if(res != RES_OK) goto error; 1046 res = ssp_rng_write(rng, stream); 1047 if(res != RES_OK) goto error; 1048 1049 exit: 1050 if(rng) SSP(rng_ref_put(rng)); 1051 return res; 1052 error: 1053 goto exit; 1054 } 1055 1056 res_T 1057 sdis_green_function_create_from_stream 1058 (struct sdis_green_function_create_from_stream_args* args, 1059 struct sdis_green_function** out_green) 1060 { 1061 hash256_T hash0, hash1; 1062 struct sdis_green_function* green = NULL; 1063 struct ssp_rng* rng = NULL; 1064 int version = 0; 1065 res_T res = RES_OK; 1066 1067 if(!out_green) { res = RES_BAD_ARG; goto error; } 1068 res = check_sdis_green_function_create_from_stream_args(args); 1069 if(res != RES_OK) goto error; 1070 1071 res = green_function_create(args->scene, args->signature, &green); 1072 if(res != RES_OK) goto error; 1073 1074 #define READ(Var, Nb) { \ 1075 if(fread((Var), sizeof(*(Var)), (Nb), args->stream) != (Nb)) { \ 1076 if(feof(args->stream)) { \ 1077 res = RES_BAD_ARG; \ 1078 } else if(ferror(args->stream)) { \ 1079 res = RES_IO_ERR; \ 1080 } else { \ 1081 res = RES_UNKNOWN_ERR; \ 1082 } \ 1083 goto error; \ 1084 } \ 1085 } (void)0 1086 1087 READ(&version, 1); 1088 if(version != SDIS_GREEN_FUNCTION_VERSION) { 1089 log_err(green->scn->dev, 1090 "%s: unexpected green function version %d. Expecting a green function " 1091 "in version %d.\n", 1092 FUNC_NAME, version, SDIS_GREEN_FUNCTION_VERSION); 1093 res = RES_BAD_ARG; 1094 goto error; 1095 } 1096 1097 res = scene_compute_hash(green->scn, hash0); 1098 if(res != RES_OK) goto error; 1099 1100 READ(hash1, sizeof(hash256_T)); 1101 if(!hash256_eq(hash0, hash1)) { 1102 log_err(green->scn->dev, 1103 "%s: the submitted scene does not match the scene used to estimate the " 1104 "green function.\n", FUNC_NAME); 1105 res = RES_BAD_ARG; 1106 goto error; 1107 } 1108 1109 READ(hash1, sizeof(hash256_T)); 1110 if(!hash256_eq(hash1, green->signature)) { 1111 log_err(green->scn->dev, 1112 "%s: the input signature does not match the saved signature\n", 1113 FUNC_NAME); 1114 res = RES_BAD_ARG; 1115 goto error; 1116 } 1117 1118 res = read_media(green, args->stream); 1119 if(res != RES_OK) goto error; 1120 res = read_interfaces(green, args->stream); 1121 if(res != RES_OK) goto error; 1122 res = read_paths_list(green, args->stream); 1123 if(res != RES_OK) goto error; 1124 1125 READ(&green->npaths_valid, 1); 1126 READ(&green->npaths_invalid, 1); 1127 READ(&green->realisation_time, 1); 1128 READ(&green->rng_type, 1); 1129 #undef READ 1130 1131 /* Create a temporary RNG used to deserialise the RNG state */ 1132 res = ssp_rng_create(green->scn->dev->allocator, green->rng_type, &rng); 1133 if(res != RES_OK) goto error; 1134 res = ssp_rng_read(rng, args->stream); 1135 if(res != RES_OK) goto error; 1136 res = ssp_rng_write(rng, green->rng_state); 1137 if(res != RES_OK) goto error; 1138 1139 exit: 1140 if(rng) SSP(rng_ref_put(rng)); 1141 if(out_green) *out_green = green; 1142 return res; 1143 error: 1144 if(green) { 1145 SDIS(green_function_ref_put(green)); 1146 green = NULL; 1147 } 1148 goto exit; 1149 } 1150 1151 res_T 1152 sdis_green_function_get_scene 1153 (const struct sdis_green_function* green, 1154 struct sdis_scene** scn) 1155 { 1156 if(!green || !scn) return RES_BAD_ARG; 1157 ASSERT(green->npaths_valid != SIZE_MAX); 1158 *scn = green->scn; 1159 return RES_OK; 1160 } 1161 1162 res_T 1163 sdis_green_function_get_paths_count 1164 (const struct sdis_green_function* green, size_t* npaths) 1165 { 1166 if(!green || !npaths) return RES_BAD_ARG; 1167 ASSERT(green->npaths_valid != SIZE_MAX); 1168 *npaths = green->npaths_valid; 1169 return RES_OK; 1170 } 1171 1172 res_T 1173 sdis_green_function_get_invalid_paths_count 1174 (const struct sdis_green_function* green, size_t* nfails) 1175 { 1176 if(!green || !nfails) return RES_BAD_ARG; 1177 ASSERT(green->npaths_invalid != SIZE_MAX); 1178 *nfails = green->npaths_invalid; 1179 return RES_OK; 1180 } 1181 1182 res_T 1183 sdis_green_function_get_signature 1184 (const struct sdis_green_function* green, hash256_T signature) 1185 { 1186 if(!green || !signature) return RES_BAD_ARG; 1187 memcpy(signature, green->signature, sizeof(hash256_T)); 1188 return RES_OK; 1189 } 1190 1191 res_T 1192 sdis_green_function_for_each_path 1193 (struct sdis_green_function* green, 1194 sdis_process_green_path_T func, 1195 void* context) 1196 { 1197 size_t npaths; 1198 size_t ipath; 1199 res_T res = RES_OK; 1200 1201 if(!green || !func) { 1202 res = RES_BAD_ARG; 1203 goto error; 1204 } 1205 1206 npaths = darray_green_path_size_get(&green->paths); 1207 FOR_EACH(ipath, 0, npaths) { 1208 struct sdis_green_path path_handle; 1209 const struct green_path* path = darray_green_path_cdata_get(&green->paths)+ipath; 1210 1211 if(path->end_type == SDIS_GREEN_PATH_END_ERROR) continue; 1212 1213 path_handle.green__ = green; 1214 path_handle.id__ = ipath; 1215 1216 res = func(&path_handle, context); 1217 if(res != RES_OK) goto error; 1218 } 1219 1220 exit: 1221 return res; 1222 error: 1223 goto exit; 1224 } 1225 1226 res_T 1227 sdis_green_path_get_elapsed_time 1228 (struct sdis_green_path* path_handle, double* elapsed) 1229 { 1230 const struct green_path* path = NULL; 1231 struct sdis_green_function* green = NULL; 1232 res_T res = RES_OK; 1233 1234 if(!path_handle || !elapsed) { 1235 res = RES_BAD_ARG; 1236 goto error; 1237 } 1238 1239 green = path_handle->green__; 1240 ASSERT(path_handle->id__ < darray_green_path_size_get(&green->paths)); 1241 1242 path = darray_green_path_cdata_get(&green->paths) + path_handle->id__; 1243 *elapsed = path->elapsed_time; 1244 1245 exit: 1246 return res; 1247 error: 1248 goto exit; 1249 } 1250 1251 res_T 1252 sdis_green_path_get_end 1253 (struct sdis_green_path* path_handle, 1254 struct sdis_green_path_end* end) 1255 { 1256 const struct green_path* path = NULL; 1257 struct sdis_green_function* green = NULL; 1258 res_T res = RES_OK; 1259 1260 if(!path_handle || !end) { 1261 res = RES_BAD_ARG; 1262 goto error; 1263 } 1264 1265 green = path_handle->green__; 1266 ASSERT(path_handle->id__ < darray_green_path_size_get(&green->paths)); 1267 1268 path = darray_green_path_cdata_get(&green->paths) + path_handle->id__; 1269 end->type = path->end_type; 1270 1271 switch(path->end_type) { 1272 case SDIS_GREEN_PATH_END_AT_INTERFACE: 1273 end->data.itfrag.intface = green_function_fetch_interf(green, path->limit_id); 1274 end->data.itfrag.fragment = path->limit.fragment; 1275 break; 1276 case SDIS_GREEN_PATH_END_AT_RADIATIVE_ENV: 1277 end->data.radenvray.radenv = green->scn->radenv; 1278 end->data.radenvray.ray = path->limit.ray; 1279 break; 1280 case SDIS_GREEN_PATH_END_IN_VOLUME: 1281 end->data.mdmvert.medium = green_function_fetch_medium(green, path->limit_id); 1282 end->data.mdmvert.vertex = path->limit.vertex; 1283 break; 1284 case SDIS_GREEN_PATH_END_ERROR: 1285 res = RES_BAD_OP; 1286 goto error; 1287 break; 1288 default: FATAL("Unreachable code.\n"); break; 1289 } 1290 1291 exit: 1292 return res; 1293 error: 1294 goto exit; 1295 } 1296 1297 res_T 1298 sdis_green_path_get_green_function 1299 (struct sdis_green_path* path_handle, 1300 struct sdis_green_function** out_green) 1301 1302 { 1303 struct sdis_green_function* green = NULL; 1304 res_T res = RES_OK; 1305 1306 if(!path_handle || !out_green) { 1307 res = RES_BAD_ARG; 1308 goto error; 1309 } 1310 1311 green = path_handle->green__; 1312 ASSERT(path_handle->id__ < darray_green_path_size_get(&green->paths)); 1313 1314 *out_green = green; 1315 1316 exit: 1317 return res; 1318 error: 1319 goto exit; 1320 } 1321 1322 res_T 1323 sdis_green_function_get_power_terms_count 1324 (const struct sdis_green_path* path_handle, 1325 size_t* nterms) 1326 { 1327 const struct green_path* path = NULL; 1328 struct sdis_green_function* green = NULL; 1329 res_T res = RES_OK; 1330 1331 if(!path_handle || !nterms) { 1332 res = RES_BAD_ARG; 1333 goto error; 1334 } 1335 1336 green = path_handle->green__; (void)green; 1337 ASSERT(path_handle->id__ < darray_green_path_size_get(&green->paths)); 1338 1339 path = darray_green_path_cdata_get(&green->paths) + path_handle->id__; 1340 1341 *nterms = darray_power_term_size_get(&path->power_terms); 1342 1343 exit: 1344 return res; 1345 error: 1346 goto exit; 1347 } 1348 1349 res_T 1350 sdis_green_function_get_flux_terms_count 1351 (const struct sdis_green_path* path_handle, 1352 size_t* nterms) 1353 { 1354 const struct green_path* path = NULL; 1355 struct sdis_green_function* green = NULL; 1356 res_T res = RES_OK; 1357 1358 if(!path_handle || !nterms) { 1359 res = RES_BAD_ARG; 1360 goto error; 1361 } 1362 1363 green = path_handle->green__; (void)green; 1364 ASSERT(path_handle->id__ < darray_green_path_size_get(&green->paths)); 1365 1366 path = darray_green_path_cdata_get(&green->paths) + path_handle->id__; 1367 1368 *nterms = darray_flux_term_size_get(&path->flux_terms); 1369 1370 exit: 1371 return res; 1372 error: 1373 goto exit; 1374 } 1375 1376 res_T 1377 sdis_green_function_get_external_flux_terms_count 1378 (const struct sdis_green_path* path_handle, 1379 size_t* nterms) 1380 { 1381 const struct green_path* path = NULL; 1382 struct sdis_green_function* green = NULL; 1383 res_T res = RES_OK; 1384 1385 if(!path_handle || !nterms) { 1386 res = RES_BAD_ARG; 1387 goto error; 1388 } 1389 1390 green = path_handle->green__; (void)green; 1391 ASSERT(path_handle->id__ < darray_green_path_size_get(&green->paths)); 1392 1393 path = darray_green_path_cdata_get(&green->paths) + path_handle->id__; 1394 1395 *nterms = darray_extflux_terms_size_get(&path->extflux_terms); 1396 1397 exit: 1398 return res; 1399 error: 1400 goto exit; 1401 } 1402 1403 res_T 1404 sdis_green_path_for_each_power_term 1405 (struct sdis_green_path* path_handle, 1406 sdis_process_medium_power_term_T func, 1407 void* context) 1408 { 1409 const struct green_path* path = NULL; 1410 struct sdis_green_function* green = NULL; 1411 const struct power_term* terms = NULL; 1412 size_t i, n; 1413 res_T res = RES_OK; 1414 1415 if(!path_handle || !func) { 1416 res = RES_BAD_ARG; 1417 goto error; 1418 } 1419 1420 green = path_handle->green__; 1421 ASSERT(path_handle->id__ < darray_green_path_size_get(&green->paths)); 1422 1423 path = darray_green_path_cdata_get(&green->paths) + path_handle->id__; 1424 1425 n = darray_power_term_size_get(&path->power_terms); 1426 terms = darray_power_term_cdata_get(&path->power_terms); 1427 FOR_EACH(i, 0, n) { 1428 struct sdis_medium* mdm; 1429 mdm = green_function_fetch_medium(green, terms[i].id); 1430 res = func(mdm, terms[i].term, context); 1431 if(res != RES_OK) goto error; 1432 } 1433 1434 exit: 1435 return res; 1436 error: 1437 goto exit; 1438 } 1439 1440 res_T 1441 sdis_green_path_for_each_flux_term 1442 (struct sdis_green_path* path_handle, 1443 sdis_process_interface_flux_term_T func, 1444 void* context) 1445 { 1446 const struct green_path* path = NULL; 1447 struct sdis_green_function* green = NULL; 1448 const struct flux_term* terms = NULL; 1449 size_t i, n; 1450 res_T res = RES_OK; 1451 1452 if(!path_handle || !func) { 1453 res = RES_BAD_ARG; 1454 goto error; 1455 } 1456 1457 green = path_handle->green__; 1458 ASSERT(path_handle->id__ < darray_green_path_size_get(&green->paths)); 1459 1460 path = darray_green_path_cdata_get(&green->paths) + path_handle->id__; 1461 1462 n = darray_flux_term_size_get(&path->flux_terms); 1463 terms = darray_flux_term_cdata_get(&path->flux_terms); 1464 FOR_EACH(i, 0, n) { 1465 struct sdis_interface* interf; 1466 interf = green_function_fetch_interf(green, terms[i].id); 1467 1468 res = func(interf, terms[i].side, terms[i].term, context); 1469 if(res != RES_OK) goto error; 1470 } 1471 1472 exit: 1473 return res; 1474 error: 1475 goto exit; 1476 } 1477 1478 res_T 1479 sdis_green_path_for_each_external_flux_terms 1480 (struct sdis_green_path* path_handle, 1481 sdis_process_external_flux_terms_T func, 1482 void* context) 1483 { 1484 const struct green_path* path = NULL; 1485 struct sdis_green_function* green = NULL; 1486 size_t i, n; 1487 res_T res = RES_OK; 1488 1489 if(!path_handle || !func) { 1490 res = RES_BAD_ARG; 1491 goto error; 1492 } 1493 1494 green = path_handle->green__; 1495 ASSERT(path_handle->id__ < darray_green_path_size_get(&green->paths)); 1496 1497 path = darray_green_path_cdata_get(&green->paths) + path_handle->id__; 1498 1499 n = darray_extflux_terms_size_get(&path->extflux_terms); 1500 if(n && !green->scn->source) { 1501 /* In can't have external flux terms without an external source */ 1502 log_err(green->scn->dev, "%s: the external source is missing\n", FUNC_NAME); 1503 res = RES_BAD_ARG; 1504 goto error; 1505 } 1506 1507 FOR_EACH(i, 0, n) { 1508 const struct sdis_green_external_flux_terms* terms = NULL; 1509 terms = darray_extflux_terms_cdata_get(&path->extflux_terms) + i; 1510 1511 res = func(green->scn->source, terms, context); 1512 if(res != RES_OK) goto error; 1513 } 1514 1515 exit: 1516 return res; 1517 error: 1518 goto exit; 1519 } 1520 1521 1522 /******************************************************************************* 1523 * Local functions 1524 ******************************************************************************/ 1525 res_T 1526 green_function_create 1527 (struct sdis_scene* scn, 1528 const hash256_T signature, 1529 struct sdis_green_function** out_green) 1530 { 1531 struct sdis_green_function* green = NULL; 1532 res_T res = RES_OK; 1533 ASSERT(scn && out_green); 1534 1535 green = MEM_CALLOC(scn->dev->allocator, 1, sizeof(*green)); 1536 if(!green) { 1537 res = RES_MEM_ERR; 1538 goto error; 1539 } 1540 ref_init(&green->ref); 1541 SDIS(scene_ref_get(scn)); 1542 green->scn = scn; 1543 htable_medium_init(scn->dev->allocator, &green->media); 1544 htable_interf_init(scn->dev->allocator, &green->interfaces); 1545 darray_green_path_init(scn->dev->allocator, &green->paths); 1546 green->npaths_valid = SIZE_MAX; 1547 green->npaths_invalid = SIZE_MAX; 1548 memcpy(green->signature, signature, sizeof(hash256_T)); 1549 1550 /* TODO replace the tmpfile. tmpfile can only be called a limited number of 1551 * times while one could create a huge amount of green functions at the same 1552 * time (e.g. for image rendering) */ 1553 green->rng_state = tmpfile(); 1554 if(!green->rng_state) { 1555 res = RES_IO_ERR; 1556 goto error; 1557 } 1558 1559 exit: 1560 *out_green = green; 1561 return res; 1562 error: 1563 if(green) { 1564 SDIS(green_function_ref_put(green)); 1565 green = NULL; 1566 } 1567 goto exit; 1568 } 1569 1570 res_T 1571 green_function_merge_and_clear 1572 (struct sdis_green_function* dst, struct sdis_green_function* src) 1573 { 1574 struct htable_medium_iterator it_medium, end_medium; 1575 struct htable_interf_iterator it_interf, end_interf; 1576 struct green_path* paths_src; 1577 struct green_path* paths_dst; 1578 size_t npaths_src; 1579 size_t npaths_dst; 1580 size_t npaths; 1581 size_t i; 1582 res_T res = RES_OK; 1583 ASSERT(dst && src); 1584 1585 if(dst == src) goto exit; 1586 1587 npaths_src = darray_green_path_size_get(&src->paths); 1588 npaths_dst = darray_green_path_size_get(&dst->paths); 1589 npaths = npaths_src + npaths_dst; 1590 1591 res = darray_green_path_resize(&dst->paths, npaths); 1592 if(res != RES_OK) goto error; 1593 1594 paths_src = darray_green_path_data_get(&src->paths); 1595 paths_dst = darray_green_path_data_get(&dst->paths) + npaths_dst; 1596 1597 FOR_EACH(i, 0, darray_green_path_size_get(&src->paths)) { 1598 res = green_path_copy_and_clear(&paths_dst[i], &paths_src[i]); 1599 if(res != RES_OK) goto error; 1600 } 1601 1602 htable_medium_begin(&src->media, &it_medium); 1603 htable_medium_end(&src->media, &end_medium); 1604 while(!htable_medium_iterator_eq(&it_medium, &end_medium)) { 1605 struct sdis_medium* medium; 1606 medium = *htable_medium_iterator_data_get(&it_medium); 1607 res = ensure_medium_registration(dst, medium); 1608 if(res != RES_OK) goto error; 1609 htable_medium_iterator_next(&it_medium); 1610 } 1611 1612 htable_interf_begin(&src->interfaces, &it_interf); 1613 htable_interf_end(&src->interfaces, &end_interf); 1614 while(!htable_interf_iterator_eq(&it_interf, &end_interf)) { 1615 struct sdis_interface* interf; 1616 interf = *htable_interf_iterator_data_get(&it_interf); 1617 res = ensure_interface_registration(dst, interf); 1618 if(res != RES_OK) goto error; 1619 htable_interf_iterator_next(&it_interf); 1620 } 1621 1622 green_function_clear(src); 1623 exit: 1624 return res; 1625 error: 1626 goto exit; 1627 } 1628 1629 res_T 1630 green_function_redux_and_clear 1631 (struct sdis_green_function* dst, 1632 struct sdis_green_function* greens[], 1633 const size_t ngreens) 1634 { 1635 size_t i; 1636 res_T res = RES_OK; 1637 ASSERT(dst && greens); 1638 1639 FOR_EACH(i, 0, ngreens) { 1640 res = green_function_merge_and_clear(dst, greens[i]); 1641 if(res != RES_OK) goto error; 1642 } 1643 1644 exit: 1645 return res; 1646 error: 1647 goto exit; 1648 } 1649 1650 res_T 1651 green_function_finalize 1652 (struct sdis_green_function* green, 1653 struct ssp_rng_proxy* proxy, 1654 const struct accum* time) 1655 { 1656 size_t i, n; 1657 res_T res = RES_OK; 1658 1659 if(!green || !proxy || !time) { 1660 res = RES_BAD_ARG; 1661 goto error; 1662 } 1663 1664 /* Save the RNG state */ 1665 SSP(rng_proxy_get_type(proxy, &green->rng_type)); 1666 res = ssp_rng_proxy_write(proxy, green->rng_state); 1667 if(res != RES_OK) goto error; 1668 1669 /* Compute the number of valid/invalid green paths */ 1670 green->npaths_valid = 0; 1671 n = darray_green_path_size_get(&green->paths); 1672 FOR_EACH(i, 0, n) { 1673 const struct green_path* path = darray_green_path_cdata_get(&green->paths)+i; 1674 green->npaths_valid += (path->end_type != SDIS_GREEN_PATH_END_ERROR); 1675 } 1676 green->npaths_invalid = n - green->npaths_valid; 1677 1678 ASSERT(green->npaths_valid == time->count); 1679 green->realisation_time = *time; 1680 1681 exit: 1682 return res; 1683 error: 1684 goto exit; 1685 } 1686 1687 res_T 1688 green_function_create_path 1689 (struct sdis_green_function* green, 1690 struct green_path_handle* handle) 1691 { 1692 size_t n; 1693 res_T res = RES_OK; 1694 ASSERT(green && handle); 1695 1696 n = darray_green_path_size_get(&green->paths); 1697 res = darray_green_path_resize(&green->paths, n+1); 1698 if(res != RES_OK) return res; 1699 1700 handle->green = green; 1701 handle->path = darray_green_path_data_get(&green->paths) + n; 1702 return RES_OK; 1703 } 1704 1705 res_T 1706 green_path_set_limit_interface_fragment 1707 (struct green_path_handle* handle, 1708 struct sdis_interface* interf, 1709 const struct sdis_interface_fragment* frag, 1710 const double elapsed_time) 1711 { 1712 res_T res = RES_OK; 1713 ASSERT(handle && interf && frag); 1714 ASSERT(handle->path->end_type == SDIS_GREEN_PATH_END_TYPES_COUNT__); 1715 res = ensure_interface_registration(handle->green, interf); 1716 if(res != RES_OK) return res; 1717 handle->path->elapsed_time = elapsed_time; 1718 handle->path->limit.fragment = *frag; 1719 handle->path->limit_id = interface_get_id(interf); 1720 handle->path->end_type = SDIS_GREEN_PATH_END_AT_INTERFACE; 1721 return RES_OK; 1722 } 1723 1724 res_T 1725 green_path_set_limit_vertex 1726 (struct green_path_handle* handle, 1727 struct sdis_medium* mdm, 1728 const struct sdis_rwalk_vertex* vert, 1729 const double elapsed_time) 1730 { 1731 res_T res = RES_OK; 1732 ASSERT(handle && mdm && vert); 1733 ASSERT(handle->path->end_type == SDIS_GREEN_PATH_END_TYPES_COUNT__); 1734 res = ensure_medium_registration(handle->green, mdm); 1735 if(res != RES_OK) return res; 1736 handle->path->elapsed_time = elapsed_time; 1737 handle->path->limit.vertex = *vert; 1738 handle->path->limit_id = medium_get_id(mdm); 1739 handle->path->end_type = SDIS_GREEN_PATH_END_IN_VOLUME; 1740 return RES_OK; 1741 } 1742 1743 res_T 1744 green_path_set_limit_radiative_ray 1745 (struct green_path_handle* handle, 1746 const struct sdis_radiative_ray* ray, 1747 const double elapsed_time) 1748 { 1749 ASSERT(handle); 1750 ASSERT(handle->path->end_type == SDIS_GREEN_PATH_END_TYPES_COUNT__); 1751 handle->path->elapsed_time = elapsed_time; 1752 handle->path->limit.ray = *ray; 1753 handle->path->end_type = SDIS_GREEN_PATH_END_AT_RADIATIVE_ENV; 1754 return RES_OK; 1755 } 1756 1757 res_T 1758 green_path_reset_limit(struct green_path_handle* handle) 1759 { 1760 ASSERT(handle); 1761 handle->path->elapsed_time = -INF; 1762 handle->path->end_type = SDIS_GREEN_PATH_END_TYPES_COUNT__; 1763 return RES_OK; 1764 } 1765 1766 res_T 1767 green_path_add_power_term 1768 (struct green_path_handle* handle, 1769 struct sdis_medium* mdm, 1770 const struct sdis_rwalk_vertex* vtx, 1771 const double val) 1772 { 1773 struct green_path* path; 1774 struct power_term* terms; 1775 size_t nterms; 1776 size_t iterm; 1777 unsigned id; 1778 res_T res = RES_OK; 1779 ASSERT(handle && mdm && vtx); 1780 1781 /* Unused position and time: the current implementation of the green function 1782 * assumes that the power is constant in space and time per medium. */ 1783 (void)vtx; 1784 1785 res = ensure_medium_registration(handle->green, mdm); 1786 if(res != RES_OK) goto error; 1787 1788 path = handle->path; 1789 terms = darray_power_term_data_get(&path->power_terms); 1790 nterms = darray_power_term_size_get(&path->power_terms); 1791 id = medium_get_id(mdm); 1792 iterm = SIZE_MAX; 1793 1794 /* Early find */ 1795 if(path->ilast_medium < nterms && terms[path->ilast_medium].id == id) { 1796 iterm = path->ilast_medium; 1797 } else { 1798 /* Linear search of the submitted medium */ 1799 FOR_EACH(iterm, 0, nterms) if(terms[iterm].id == id) break; 1800 } 1801 1802 /* Add the power term to the path wrt the submitted medium */ 1803 if(iterm < nterms) { 1804 terms[iterm].term += val; 1805 } else { 1806 struct power_term term = POWER_TERM_NULL__; 1807 term.term = val; 1808 term.id = id; 1809 res = darray_power_term_push_back(&handle->path->power_terms, &term); 1810 if(res != RES_OK) goto error; 1811 } 1812 1813 /* Register the slot into which the last accessed medium lies */ 1814 CHK(iterm < UINT16_MAX); 1815 path->ilast_medium = (uint16_t)iterm; 1816 1817 exit: 1818 return res; 1819 error: 1820 goto exit; 1821 } 1822 1823 res_T 1824 green_path_add_flux_term 1825 (struct green_path_handle* handle, 1826 struct sdis_interface* interf, 1827 const struct sdis_interface_fragment* frag, 1828 const double val) 1829 { 1830 struct green_path* path; 1831 struct flux_term* terms; 1832 size_t nterms; 1833 size_t iterm; 1834 unsigned id; 1835 res_T res = RES_OK; 1836 ASSERT(handle && interf && frag && val >= 0); 1837 1838 res = ensure_interface_registration(handle->green, interf); 1839 if(res != RES_OK) goto error; 1840 1841 path = handle->path; 1842 terms = darray_flux_term_data_get(&path->flux_terms); 1843 nterms = darray_flux_term_size_get(&path->flux_terms); 1844 id = interface_get_id(interf); 1845 iterm = SIZE_MAX; 1846 1847 /* Early find */ 1848 if(path->ilast_interf < nterms 1849 && terms[path->ilast_interf].id == id 1850 && terms[path->ilast_interf].side == frag->side) { 1851 iterm = path->ilast_interf; 1852 } else { 1853 /* Linear search of the submitted interface */ 1854 FOR_EACH(iterm, 0, nterms) { 1855 if(terms[iterm].id == id && terms[iterm].side == frag->side) { 1856 break; 1857 } 1858 } 1859 } 1860 1861 /* Add the flux term to the path wrt the submitted interface */ 1862 if(iterm < nterms) { 1863 terms[iterm].term += val; 1864 } else { 1865 struct flux_term term = FLUX_TERM_NULL__; 1866 term.term = val; 1867 term.id = id; 1868 term.side = frag->side; 1869 res = darray_flux_term_push_back(&handle->path->flux_terms, &term); 1870 if(res != RES_OK) goto error; 1871 } 1872 1873 /* Register the slot into which the last accessed interface lies */ 1874 CHK(iterm < UINT16_MAX); 1875 path->ilast_interf = (uint16_t)iterm; 1876 1877 exit: 1878 return res; 1879 error: 1880 goto exit; 1881 } 1882 1883 res_T 1884 green_path_add_external_flux_terms 1885 (struct green_path_handle* handle, 1886 const struct sdis_green_external_flux_terms* terms) 1887 { 1888 res_T res = RES_OK; 1889 ASSERT(handle && terms); 1890 1891 res = darray_extflux_terms_push_back(&handle->path->extflux_terms, terms); 1892 if(res != RES_OK) { 1893 log_err(handle->green->scn->dev, 1894 "%s: cannot store external flux terms -- %s\n", 1895 FUNC_NAME, res_to_cstr(res)); 1896 goto error; 1897 } 1898 1899 exit: 1900 return res; 1901 error: 1902 goto exit; 1903 }