test_sdis_utils.c (16652B)
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 "test_sdis_utils.h" 17 #include <rsys/math.h> 18 19 enum heat_vertex_attrib { 20 HEAT_VERTEX_BRANCH_ID, 21 HEAT_VERTEX_WEIGHT, 22 HEAT_VERTEX_TIME, 23 HEAT_VERTEX_TYPE 24 }; 25 26 /******************************************************************************* 27 * Helper functions 28 ******************************************************************************/ 29 struct green_accum { double sum, sum2; }; 30 31 static res_T 32 count_green_paths(struct sdis_green_path* path, void* ctx) 33 { 34 CHK(path && ctx); 35 *((size_t*)ctx) += 1; 36 return RES_OK; 37 } 38 39 static res_T 40 accum_power_terms(struct sdis_medium* mdm, const double power_term, void* ctx) 41 { 42 struct sdis_solid_shader shader = SDIS_SOLID_SHADER_NULL; 43 struct sdis_rwalk_vertex vtx = SDIS_RWALK_VERTEX_NULL; 44 struct sdis_data* data = NULL; 45 double* power = ctx; /* Power contribution [K] */ 46 47 CHK(mdm && ctx); 48 CHK(sdis_medium_get_type(mdm) == SDIS_SOLID); 49 50 OK(sdis_solid_get_shader(mdm, &shader)); 51 data = sdis_medium_get_data(mdm); 52 vtx.time = INF; 53 54 *power += power_term * shader.volumic_power(&vtx, data); 55 return RES_OK; 56 } 57 58 static res_T 59 accum_flux_terms 60 (struct sdis_interface* interf, 61 const enum sdis_side side, 62 const double flux_term, 63 void* ctx) 64 { 65 struct sdis_interface_shader shader = SDIS_INTERFACE_SHADER_NULL; 66 struct sdis_interface_fragment frag = SDIS_INTERFACE_FRAGMENT_NULL; 67 struct sdis_data* data = NULL; 68 double* flux = ctx; /* Flux contribution [K] */ 69 double phi; 70 71 CHK(interf && ctx); 72 73 OK(sdis_interface_get_shader(interf, &shader)); 74 data = sdis_interface_get_data(interf); 75 frag.time = INF; 76 frag.side = side; 77 78 phi = side == SDIS_FRONT 79 ? shader.front.flux(&frag, data) 80 : shader.back.flux(&frag, data); 81 82 *flux += flux_term * phi; 83 return RES_OK; 84 } 85 86 static res_T 87 accum_extflux 88 (struct sdis_source* source, 89 const struct sdis_green_external_flux_terms* terms, 90 void* ctx) 91 { 92 struct sdis_spherical_source_shader shader = SDIS_SPHERICAL_SOURCE_SHADER_NULL; 93 struct sdis_data* data = NULL; 94 double* extflux = ctx; /* External flux contribution [K] */ 95 double power = 0; /* [W] */ 96 double diffuse_radiance = 0; /* [W/m^2/sr] */ 97 98 CHK(source && terms && ctx); 99 100 data = sdis_source_get_data(source); 101 OK(sdis_spherical_source_get_shader(source, &shader)); 102 power = shader.power(terms->time, data); 103 if(shader.diffuse_radiance) { 104 diffuse_radiance = shader.diffuse_radiance(terms->time, terms->dir, data); 105 } 106 107 *extflux += terms->term_wrt_power * power; 108 *extflux += terms->term_wrt_diffuse_radiance * diffuse_radiance; 109 return RES_OK; 110 } 111 112 static res_T 113 solve_green_path(struct sdis_green_path* path, void* ctx) 114 { 115 struct sdis_green_path_end end = SDIS_GREEN_PATH_END_NULL; 116 117 struct sdis_rwalk_vertex vtx = SDIS_RWALK_VERTEX_NULL; 118 struct sdis_interface_fragment frag = SDIS_INTERFACE_FRAGMENT_NULL; 119 struct sdis_radiative_ray ray = SDIS_RADIATIVE_RAY_NULL; 120 121 struct sdis_solid_shader solid = SDIS_SOLID_SHADER_NULL; 122 struct sdis_fluid_shader fluid = SDIS_FLUID_SHADER_NULL; 123 struct sdis_interface_shader interf = SDIS_INTERFACE_SHADER_NULL; 124 struct sdis_radiative_env_shader radenv = SDIS_RADIATIVE_ENV_SHADER_NULL; 125 126 struct sdis_green_function* green = NULL; 127 struct sdis_scene* scn = NULL; 128 struct green_accum* acc = NULL; 129 struct sdis_data* data = NULL; 130 enum sdis_medium_type type; 131 double power = 0; 132 double flux = 0; 133 double extflux = 0; 134 double time, temp = 0; 135 double weight = 0; 136 CHK(path && ctx); 137 138 acc = ctx; 139 140 BA(sdis_green_path_get_green_function(NULL, NULL)); 141 BA(sdis_green_path_get_green_function(path, NULL)); 142 BA(sdis_green_path_get_green_function(NULL, &green)); 143 OK(sdis_green_path_get_green_function(path, &green)); 144 145 BA(sdis_green_function_get_scene(NULL, NULL)); 146 BA(sdis_green_function_get_scene(NULL, &scn)); 147 BA(sdis_green_function_get_scene(green, NULL)); 148 OK(sdis_green_function_get_scene(green, &scn)); 149 150 BA(sdis_green_path_for_each_power_term(NULL, accum_power_terms, &power)); 151 BA(sdis_green_path_for_each_power_term(path, NULL, &acc)); 152 OK(sdis_green_path_for_each_power_term(path, accum_power_terms, &power)); 153 154 BA(sdis_green_path_for_each_flux_term(NULL, accum_flux_terms, &flux)); 155 BA(sdis_green_path_for_each_flux_term(path, NULL, &acc)); 156 OK(sdis_green_path_for_each_flux_term(path, accum_flux_terms, &flux)); 157 158 BA(sdis_green_path_for_each_external_flux_terms(NULL, &accum_extflux, &extflux)); 159 BA(sdis_green_path_for_each_external_flux_terms(path, NULL, &extflux)); 160 OK(sdis_green_path_for_each_external_flux_terms(path, &accum_extflux, &extflux)); 161 162 BA(sdis_green_path_get_elapsed_time(NULL, NULL)); 163 BA(sdis_green_path_get_elapsed_time(path, NULL)); 164 BA(sdis_green_path_get_elapsed_time(NULL, &time)); 165 OK(sdis_green_path_get_elapsed_time(path, &time)); 166 167 BA(sdis_green_path_get_end(NULL, NULL)); 168 BA(sdis_green_path_get_end(NULL, &end)); 169 BA(sdis_green_path_get_end(path, NULL)); 170 OK(sdis_green_path_get_end(path, &end)); 171 switch(end.type) { 172 case SDIS_GREEN_PATH_END_AT_INTERFACE: 173 frag = end.data.itfrag.fragment; 174 OK(sdis_interface_get_shader(end.data.itfrag.intface, &interf)); 175 data = sdis_interface_get_data(end.data.itfrag.intface); 176 temp = frag.side == SDIS_FRONT 177 ? interf.front.temperature(&frag, data) 178 : interf.back.temperature(&frag, data); 179 break; 180 181 case SDIS_GREEN_PATH_END_AT_RADIATIVE_ENV: 182 ray = end.data.radenvray.ray; 183 OK(sdis_radiative_env_get_shader(end.data.radenvray.radenv, &radenv)); 184 data = sdis_radiative_env_get_data(end.data.radenvray.radenv); 185 temp = radenv.temperature(&ray, data); 186 break; 187 188 case SDIS_GREEN_PATH_END_IN_VOLUME: 189 vtx = end.data.mdmvert.vertex; 190 type = sdis_medium_get_type(end.data.mdmvert.medium); 191 data = sdis_medium_get_data(end.data.mdmvert.medium); 192 if(type == SDIS_FLUID) { 193 OK(sdis_fluid_get_shader(end.data.mdmvert.medium, &fluid)); 194 temp = fluid.temperature(&vtx, data); 195 } else { 196 OK(sdis_solid_get_shader(end.data.mdmvert.medium, &solid)); 197 temp = solid.temperature(&vtx, data); 198 } 199 break; 200 201 default: FATAL("Unreachable code.\n"); break; 202 } 203 204 weight = temp + power + extflux + flux; 205 acc->sum += weight; 206 acc->sum2 += weight*weight; 207 208 return RES_OK; 209 } 210 211 static size_t 212 heat_path_get_vertices_count(const struct sdis_heat_path* path) 213 { 214 size_t istrip = 0; 215 size_t nstrips = 0; 216 size_t nvertices = 0; 217 CHK(path); 218 219 OK(sdis_heat_path_get_line_strips_count(path, &nstrips)); 220 FOR_EACH(istrip, 0, nstrips) { 221 size_t n; 222 OK(sdis_heat_path_line_strip_get_vertices_count(path, istrip, &n)); 223 nvertices += n; 224 } 225 return nvertices; 226 } 227 228 static void 229 dump_heat_path_positions(FILE* stream, const struct sdis_heat_path* path) 230 { 231 size_t istrip, nstrips; 232 size_t ivert, nverts; 233 234 CHK(stream && path); 235 236 OK(sdis_heat_path_get_line_strips_count(path, &nstrips)); 237 FOR_EACH(istrip, 0, nstrips) { 238 OK(sdis_heat_path_line_strip_get_vertices_count(path, istrip, &nverts)); 239 FOR_EACH(ivert, 0, nverts) { 240 struct sdis_heat_vertex vtx; 241 OK(sdis_heat_path_line_strip_get_vertex(path, istrip, ivert, &vtx)); 242 fprintf(stream, "%g %g %g\n", SPLIT3(vtx.P)); 243 } 244 } 245 } 246 247 static void 248 dump_heat_path_line_strip 249 (FILE* stream, 250 const struct sdis_heat_path* path, 251 const size_t istrip, 252 const size_t offset) 253 { 254 size_t ivert, nverts; 255 256 CHK(stream); 257 258 OK(sdis_heat_path_line_strip_get_vertices_count(path, istrip, &nverts)); 259 fprintf(stream, "%lu", (unsigned long)nverts); 260 FOR_EACH(ivert, 0, nverts) { 261 fprintf(stream, " %lu", (unsigned long)(ivert + offset)); 262 } 263 fprintf(stream, "\n"); 264 } 265 266 static void 267 dump_heat_path_vertex_attribs 268 (FILE* stream, 269 const struct sdis_heat_path* path, 270 const enum heat_vertex_attrib attr) 271 { 272 size_t ivert, nverts; 273 size_t istrip, nstrips; 274 CHK(stream && path); 275 276 OK(sdis_heat_path_get_line_strips_count(path, &nstrips)); 277 FOR_EACH(istrip, 0, nstrips) { 278 OK(sdis_heat_path_line_strip_get_vertices_count(path, istrip, &nverts)); 279 FOR_EACH(ivert, 0, nverts) { 280 struct sdis_heat_vertex vtx; 281 OK(sdis_heat_path_line_strip_get_vertex(path, istrip, ivert, &vtx)); 282 switch(attr) { 283 case HEAT_VERTEX_BRANCH_ID: 284 fprintf(stream, "%i\n", vtx.branch_id); 285 break; 286 case HEAT_VERTEX_WEIGHT: 287 fprintf(stream, "%g\n", vtx.weight); 288 break; 289 case HEAT_VERTEX_TIME: 290 fprintf(stream, "%g\n", 291 IS_INF(vtx.time) || vtx.time > FLT_MAX ? -1 : vtx.time); 292 break; 293 case HEAT_VERTEX_TYPE: 294 switch(vtx.type) { 295 case SDIS_HEAT_VERTEX_CONDUCTION: fprintf(stream, "0.0\n"); break; 296 case SDIS_HEAT_VERTEX_CONVECTION: fprintf(stream, "0.5\n"); break; 297 case SDIS_HEAT_VERTEX_RADIATIVE: fprintf(stream, "1.0\n"); break; 298 default: FATAL("Unreachable code.\n"); break; 299 } 300 break; 301 default: FATAL("Unreachable code.\n"); break; 302 } 303 } 304 } 305 } 306 307 /******************************************************************************* 308 * Local function 309 ******************************************************************************/ 310 void 311 check_green_function(struct sdis_green_function* green) 312 { 313 double time_range[2]; 314 struct sdis_estimator* estimator; 315 struct sdis_mc mc; 316 struct green_accum accum = {0, 0}; 317 double E, V, SE; 318 size_t nreals; 319 size_t nfails; 320 size_t n; 321 322 time_range[0] = time_range[1] = INF; 323 324 OK(sdis_green_function_solve(green, &estimator)); 325 326 BA(sdis_green_function_get_paths_count(NULL, &n)); 327 BA(sdis_green_function_get_paths_count(green, NULL)); 328 OK(sdis_green_function_get_paths_count(green, &n)); 329 OK(sdis_estimator_get_realisation_count(estimator, &nreals)); 330 CHK(n == nreals); 331 332 BA(sdis_green_function_get_invalid_paths_count(NULL, &n)); 333 BA(sdis_green_function_get_invalid_paths_count(green, NULL)); 334 OK(sdis_green_function_get_invalid_paths_count(green, &n)); 335 OK(sdis_estimator_get_failure_count(estimator, &nfails)); 336 CHK(n == nfails); 337 338 n = 0; 339 BA(sdis_green_function_for_each_path(NULL, count_green_paths, &n)); 340 BA(sdis_green_function_for_each_path(green, NULL, &n)); 341 OK(sdis_green_function_for_each_path(green, count_green_paths, &n)); 342 CHK(n == nreals); 343 344 OK(sdis_green_function_for_each_path(green, solve_green_path, &accum)); 345 346 E = accum.sum / (double)n; 347 V = MMAX(0, accum.sum2 / (double)n - E*E); 348 SE = sqrt(V/(double)n); 349 OK(sdis_estimator_get_temperature(estimator, &mc)); 350 351 printf("Green: rebuild = %g +/- %g; solved = %g +/- %g\n", 352 E, SE, mc.E, mc.SE); 353 354 CHK(E + SE >= mc.E - mc.SE); 355 CHK(E - SE <= mc.E + mc.SE); 356 357 OK(sdis_estimator_get_realisation_time(estimator, &mc)); 358 printf("Green per realisation time (in usec) = %g +/- %g\n", 359 mc.E, mc.SE); 360 361 OK(sdis_estimator_ref_put(estimator)); 362 } 363 364 void 365 dump_heat_paths(FILE* stream, const struct sdis_estimator* estimator) 366 { 367 const struct sdis_heat_path* path; 368 size_t ipath; 369 size_t npaths; 370 size_t nstrips_overall; 371 size_t nvertices; 372 size_t offset; 373 CHK(stream && estimator); 374 375 OK(sdis_estimator_get_paths_count(estimator, &npaths)); 376 CHK(npaths); 377 378 /* Header */ 379 fprintf(stream, "# vtk DataFile Version 2.0\n"); 380 fprintf(stream, "Heat paths\n"); 381 fprintf(stream, "ASCII\n"); 382 fprintf(stream, "DATASET POLYDATA\n"); 383 384 /* Compute the overall number of vertices and the overall number of strips */ 385 nvertices = 0; 386 nstrips_overall = 0; 387 FOR_EACH(ipath, 0, npaths) { 388 size_t n; 389 OK(sdis_estimator_get_path(estimator, ipath, &path)); 390 nvertices += heat_path_get_vertices_count(path); 391 OK(sdis_heat_path_get_line_strips_count(path, &n)); 392 nstrips_overall += n; 393 } 394 395 /* Write path positions */ 396 fprintf(stream, "POINTS %lu double\n", (unsigned long)nvertices); 397 FOR_EACH(ipath, 0, npaths) { 398 OK(sdis_estimator_get_path(estimator, ipath, &path)); 399 dump_heat_path_positions(stream, path); 400 } 401 402 /* Write the segment of the paths */ 403 fprintf(stream, "LINES %lu %lu\n", 404 (unsigned long)(nstrips_overall), 405 (unsigned long)(nstrips_overall + nvertices)); 406 offset = 0; 407 FOR_EACH(ipath, 0, npaths) { 408 size_t path_istrip; 409 size_t path_nstrips; 410 OK(sdis_estimator_get_path(estimator, ipath, &path)); 411 OK(sdis_heat_path_get_line_strips_count(path, &path_nstrips)); 412 FOR_EACH(path_istrip, 0, path_nstrips) { 413 size_t n; 414 dump_heat_path_line_strip(stream, path, path_istrip, offset); 415 OK(sdis_heat_path_line_strip_get_vertices_count(path, path_istrip, &n)); 416 offset += n; 417 } 418 } 419 420 fprintf(stream, "POINT_DATA %lu\n", (unsigned long)nvertices); 421 422 /* Write the type of the random walk vertices */ 423 fprintf(stream, "SCALARS Vertex_Type float 1\n"); 424 fprintf(stream, "LOOKUP_TABLE vertex_type\n"); 425 FOR_EACH(ipath, 0, npaths) { 426 OK(sdis_estimator_get_path(estimator, ipath, &path)); 427 dump_heat_path_vertex_attribs(stream, path, HEAT_VERTEX_TYPE); 428 } 429 fprintf(stream, "LOOKUP_TABLE vertex_type 3\n"); 430 fprintf(stream, "0.0 1.0 1.0 1.0\n"); /* 0.0 = Magenta: conduction */ 431 fprintf(stream, "1.0 1.0 0.0 1.0\n"); /* 0.5 = Yellow: convection */ 432 fprintf(stream, "1.0 0.0 1.0 1.0\n"); /* 1.0 = Purple: radiative */ 433 434 /* Write the weights of the random walk vertices */ 435 fprintf(stream, "SCALARS Weight double 1\n"); 436 fprintf(stream, "LOOKUP_TABLE default\n"); 437 FOR_EACH(ipath, 0, npaths) { 438 OK(sdis_estimator_get_path(estimator, ipath, &path)); 439 dump_heat_path_vertex_attribs(stream, path, HEAT_VERTEX_WEIGHT); 440 } 441 442 /* Write the time of the random walk vertices */ 443 fprintf(stream, "SCALARS Time double 1\n"); 444 fprintf(stream, "LOOKUP_TABLE default\n"); 445 FOR_EACH(ipath, 0, npaths) { 446 OK(sdis_estimator_get_path(estimator, ipath, &path)); 447 dump_heat_path_vertex_attribs(stream, path, HEAT_VERTEX_TIME); 448 } 449 450 /* Write the branch id of the random walk vertices */ 451 fprintf(stream, "SCALARS BranchID int 1\n"); 452 fprintf(stream, "LOOKUP_TABLE default\n"); 453 FOR_EACH(ipath, 0, npaths) { 454 OK(sdis_estimator_get_path(estimator, ipath, &path)); 455 dump_heat_path_vertex_attribs(stream, path, HEAT_VERTEX_BRANCH_ID); 456 } 457 458 /* Write path type */ 459 fprintf(stream, "CELL_DATA %lu\n", (unsigned long)nstrips_overall); 460 fprintf(stream, "SCALARS Path_Type float 1\n"); 461 fprintf(stream, "LOOKUP_TABLE path_type\n"); 462 FOR_EACH(ipath, 0, npaths) { 463 size_t path_istrip; 464 size_t path_nstrips; 465 enum sdis_heat_path_flag status = SDIS_HEAT_PATH_NONE; 466 OK(sdis_estimator_get_path(estimator, ipath, &path)); 467 OK(sdis_heat_path_get_status(path, &status)); 468 OK(sdis_heat_path_get_line_strips_count(path, &path_nstrips)); 469 FOR_EACH(path_istrip, 0, path_nstrips) { 470 switch(status) { 471 case SDIS_HEAT_PATH_SUCCESS: fprintf(stream, "0.0\n"); break; 472 case SDIS_HEAT_PATH_FAILURE: fprintf(stream, "1.0\n"); break; 473 default: FATAL("Unreachable code.\n"); break; 474 } 475 } 476 } 477 fprintf(stream, "LOOKUP_TABLE path_type 2\n"); 478 fprintf(stream, "0.0 0.0 1.0 1.0\n"); /* 0.0 = Blue: success */ 479 fprintf(stream, "1.0 0.0 0.0 1.0\n"); /* 1.0 = Red: failure */ 480 } 481 482 void 483 check_green_serialization 484 (struct sdis_green_function* green, 485 struct sdis_scene* scn) 486 { 487 struct sdis_green_function_create_from_stream_args args = 488 SDIS_GREEN_FUNCTION_CREATE_FROM_STREAM_ARGS_DEFAULT; 489 FILE* stream = NULL; 490 struct sdis_estimator *e1 = NULL; 491 struct sdis_estimator *e2 = NULL; 492 struct sdis_green_function* green2 = NULL; 493 494 CHK(green && scn); 495 stream = tmpfile(); 496 CHK(stream); 497 498 OK(sdis_green_function_write(green, stream)); 499 500 rewind(stream); 501 args.scene = scn; 502 args.stream = stream; 503 OK(sdis_green_function_create_from_stream(&args, &green2)); 504 CHK(!fclose(stream)); 505 check_green_function(green2); 506 507 OK(sdis_green_function_solve(green, &e1)); 508 OK(sdis_green_function_solve(green2, &e2)); 509 check_estimator_eq_strict(e1, e2); 510 511 OK(sdis_estimator_ref_put(e1)); 512 OK(sdis_estimator_ref_put(e2)); 513 OK(sdis_green_function_ref_put(green2)); 514 }