test_sdis_unsteady_atm.c (30904B)
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.h" 17 #include "test_sdis_utils.h" 18 19 #include <rsys/clock_time.h> 20 #include <rsys/mem_allocator.h> 21 #include <rsys/double3.h> 22 #include <rsys/math.h> 23 #include <star/ssp.h> 24 25 /* 26 * The physical configuration is the following: a slab of fluid with known 27 * thermophysical properties but unknown temperature is located between a 28 * "ground" and a slab of solid, with also a unknown temperature profile. On 29 * the other side of the solid slab, is an "atmosphere" with known temperature, 30 * and known radiative temperature. 31 * 32 * Solving the system means: finding the temperature of the ground, of the 33 * fluid, of the boundaries, and also the temperature inside the solid, at 34 * various locations (the 1D slab is discretized in order to obtain the 35 * reference) 36 * 37 * The reference for this system comes from a numerical method and is not 38 * analytic. Thus the compliance test MC VS reference is not the usual |MC - 39 * ref| <= 3*sigma but is |MC -ref| <= (Tmax -Tmin) * 0.01. 40 * 41 * 3D 2D 42 * 43 * /////////////// /////////////// 44 * +-----------+-+ +-----------+-+ 45 * /' / /| | | | 46 * +-----------+-+ | HA _\ <---- TR | _\ | | HA _\ <---- TR 47 * | | _\ | | | / / <---- TR TG| HG / / HC | | / / <---- TR 48 * | |HG / / HC| | | TA \__/ <---- TR | \__/ | | TA \__/ <---- TR 49 * TG| | \__/ | | | | | | 50 * | +.........|.|.+ +-----------+-+ 51 * |, |/|/ /////H///////E/ 52 * +-----------+-+ 53 * //////H//////E/ 54 */ 55 56 #define XH 3 57 #define XHpE 3.2 58 #define XE (XHpE - XH) 59 60 #define T0_SOLID 300 61 #define T0_FLUID 300 62 63 #define N 10000 /* #realisations */ 64 65 #define TG 310 66 #define HG 400 67 68 #define HC 400 69 70 #define TA 290 71 #define HA 400 72 #define TR 260 73 74 #define TMAX (MMAX(MMAX(MMAX(T0_FLUID, T0_SOLID), MMAX(TG, TA)), TR)) 75 #define TMIN (MMIN(MMIN(MMIN(T0_FLUID, T0_SOLID), MMIN(TG, TA)), TR)) 76 #define EPS ((TMAX-TMIN)*0.01) 77 78 /* hr = 4.0 * BOLTZMANN_CONSTANT * Tref * Tref * Tref * epsilon 79 * Tref = (hr / (4 * 5.6696e-8 * epsilon)) ^ 1/3, hr = 6 */ 80 #define TREF 297.974852286 81 82 #define RHO_F 1.3 83 #define CP_F 1005 84 #define RHO_S 2400 85 #define CP_S 800 86 #define LAMBDA 0.6 87 88 #define X_PROBE (XH + 0.2 * XE) 89 90 #define DELTA (XE/40.0) 91 92 /******************************************************************************* 93 * Box geometry 94 ******************************************************************************/ 95 static const double model3d_vertices[12/*#vertices*/*3/*#coords per vertex*/] = { 96 0, 0, 0, 97 XH, 0, 0, 98 XHpE, 0, 0, 99 0, XHpE, 0, 100 XH, XHpE, 0, 101 XHpE, XHpE, 0, 102 0, 0, XHpE, 103 XH, 0, XHpE, 104 XHpE, 0, XHpE, 105 0, XHpE, XHpE, 106 XH, XHpE, XHpE, 107 XHpE, XHpE, XHpE 108 }; 109 static const size_t model3d_nvertices = sizeof(model3d_vertices)/(sizeof(double)*3); 110 111 /* The following array lists the indices toward the 3D vertices of each 112 * triangle. 113 * ,3---,4---,5 ,3----4----5 ,4 114 * ,' | ,' | ,'/| ,'/| \ | \ | ,'/| 115 * 9----10---11 / | 9' / | \ | \ | 10 / | Y 116 * |', |', | / ,2 | / ,0---,1---,2 | / ,1 | 117 * | ',| ',|/,' |/,' | ,' | ,' |/,' o--X 118 * 6----7----8' 6----7'---8' 7 / 119 * Front, right Back, left and Internal Z 120 * and Top faces bottom faces face */ 121 static const size_t model3d_indices[22/*#triangles*/*3/*#indices per triangle*/] = { 122 0, 3, 1, 1, 3, 4, 1, 4, 2, 2, 4, 5, /* -Z */ 123 0, 6, 3, 3, 6, 9, /* -X */ 124 6, 7, 9, 9, 7, 10, 7, 8, 10, 10, 8, 11, /* +Z */ 125 5, 11, 8, 8, 2, 5, /* +X */ 126 3, 9, 10, 10, 4, 3, 4, 10, 11, 11, 5, 4, /* +Y */ 127 0, 1, 7, 7, 6, 0, 1, 2, 8, 8, 7, 1, /* -Y */ 128 4, 10, 7, 7, 1, 4 /* Inside */ 129 }; 130 static const size_t model3d_ntriangles = sizeof(model3d_indices)/(sizeof(size_t)*3); 131 132 static INLINE void 133 model3d_get_indices(const size_t itri, size_t ids[3], void* context) 134 { 135 (void)context; 136 CHK(ids); 137 CHK(itri < model3d_ntriangles); 138 ids[0] = model3d_indices[itri * 3 + 0]; 139 ids[1] = model3d_indices[itri * 3 + 1]; 140 ids[2] = model3d_indices[itri * 3 + 2]; 141 } 142 143 static INLINE void 144 model3d_get_position(const size_t ivert, double pos[3], void* context) 145 { 146 (void)context; 147 CHK(pos); 148 CHK(ivert < model3d_nvertices); 149 pos[0] = model3d_vertices[ivert * 3 + 0]; 150 pos[1] = model3d_vertices[ivert * 3 + 1]; 151 pos[2] = model3d_vertices[ivert * 3 + 2]; 152 } 153 154 static INLINE void 155 model3d_get_interface(const size_t itri, struct sdis_interface** bound, void* context) 156 { 157 struct sdis_interface** interfaces = context; 158 CHK(context && bound); 159 CHK(itri < model3d_ntriangles); 160 *bound = interfaces[itri]; 161 } 162 163 /******************************************************************************* 164 * Square geometry 165 ******************************************************************************/ 166 static const double model2d_vertices[6/*#vertices*/ * 2/*#coords per vertex*/] = { 167 XHpE, 0, 168 XH, 0, 169 0, 0, 170 0, XHpE, 171 XH, XHpE, 172 XHpE, XHpE 173 }; 174 static const size_t model2d_nvertices = sizeof(model2d_vertices)/(sizeof(double)*2); 175 176 static const size_t model2d_indices[7/*#segments*/ * 2/*#indices per segment*/] = { 177 0, 1, 1, 2, /* Bottom */ 178 2, 3, /* Left */ 179 3, 4, 4, 5, /* Top */ 180 5, 0, /* Right */ 181 4, 1 /* Inside */ 182 }; 183 static const size_t model2d_nsegments = sizeof(model2d_indices)/(sizeof(size_t)*2); 184 185 static INLINE void 186 model2d_get_indices(const size_t iseg, size_t ids[2], void* context) 187 { 188 (void)context; 189 CHK(ids); 190 CHK(iseg < model2d_nsegments); 191 ids[0] = model2d_indices[iseg * 2 + 0]; 192 ids[1] = model2d_indices[iseg * 2 + 1]; 193 } 194 195 static INLINE void 196 model2d_get_position(const size_t ivert, double pos[2], void* context) 197 { 198 (void)context; 199 CHK(pos); 200 CHK(ivert < model2d_nvertices); 201 pos[0] = model2d_vertices[ivert * 2 + 0]; 202 pos[1] = model2d_vertices[ivert * 2 + 1]; 203 } 204 205 static INLINE void 206 model2d_get_interface 207 (const size_t iseg, struct sdis_interface** bound, void* context) 208 { 209 struct sdis_interface** interfaces = context; 210 CHK(context && bound); 211 CHK(iseg < model2d_nsegments); 212 *bound = interfaces[iseg]; 213 } 214 215 /******************************************************************************* 216 * Media 217 ******************************************************************************/ 218 struct solid { 219 double lambda; 220 double rho; 221 double cp; 222 double delta; 223 double temperature; 224 double t0; 225 }; 226 227 static double 228 solid_get_calorific_capacity 229 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 230 { 231 struct solid* solid; 232 CHK(vtx && data); 233 solid = ((struct solid*)sdis_data_cget(data)); 234 return solid->cp; 235 } 236 237 static double 238 solid_get_thermal_conductivity 239 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 240 { 241 struct solid* solid; 242 CHK(vtx && data); 243 solid = ((struct solid*)sdis_data_cget(data)); 244 return solid->lambda; 245 } 246 247 static double 248 solid_get_volumic_mass 249 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 250 { 251 struct solid* solid; 252 CHK(vtx && data); 253 solid = ((struct solid*)sdis_data_cget(data)); 254 return solid->rho; 255 } 256 257 static double 258 solid_get_delta 259 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 260 { 261 struct solid* solid; 262 CHK(vtx && data); 263 solid = ((struct solid*)sdis_data_cget(data)); 264 return solid->delta; 265 } 266 267 static double 268 solid_get_temperature 269 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 270 { 271 struct solid* solid; 272 CHK(vtx && data); 273 solid = ((struct solid*)sdis_data_cget(data)); 274 if(vtx->time <= solid->t0) 275 return solid->temperature; 276 return SDIS_TEMPERATURE_NONE; 277 } 278 279 struct fluid { 280 double rho; 281 double cp; 282 double t0; 283 double temperature; 284 }; 285 286 static double 287 fluid_get_temperature 288 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 289 { 290 struct fluid* fluid; 291 CHK(vtx && data); 292 fluid = ((struct fluid*)sdis_data_cget(data)); 293 if(vtx->time <= fluid->t0) 294 return fluid->temperature; 295 return SDIS_TEMPERATURE_NONE; 296 } 297 298 static double 299 fluid_get_volumic_mass 300 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 301 { 302 struct fluid* fluid; 303 CHK(vtx && data); 304 fluid = ((struct fluid*)sdis_data_cget(data)); 305 return fluid->rho; 306 } 307 308 static double 309 fluid_get_calorific_capacity 310 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 311 { 312 struct fluid* fluid; 313 CHK(vtx && data); 314 fluid = ((struct fluid*)sdis_data_cget(data)); 315 return fluid->cp; 316 } 317 318 /******************************************************************************* 319 * Interfaces 320 ******************************************************************************/ 321 struct interf { 322 double temperature; 323 double emissivity; 324 double h; 325 double Tref; 326 }; 327 328 static double 329 interface_get_temperature 330 (const struct sdis_interface_fragment* frag, struct sdis_data* data) 331 { 332 const struct interf* interf; 333 CHK(frag && data); 334 interf = sdis_data_cget(data); 335 return interf->temperature; 336 } 337 338 static double 339 interface_get_convection_coef 340 (const struct sdis_interface_fragment* frag, struct sdis_data* data) 341 { 342 const struct interf* interf; 343 CHK(frag && data); 344 interf = sdis_data_cget(data); 345 return interf->h; 346 } 347 348 static double 349 interface_get_emissivity 350 (const struct sdis_interface_fragment* frag, 351 const unsigned source_id, 352 struct sdis_data* data) 353 { 354 const struct interf* interf; 355 (void)source_id; 356 CHK(frag && data); 357 interf = sdis_data_cget(data); 358 return interf->emissivity; 359 } 360 361 static double 362 interface_get_Tref 363 (const struct sdis_interface_fragment* frag, 364 struct sdis_data* data) 365 { 366 const struct interf* interf; 367 CHK(frag && data); 368 interf = sdis_data_cget(data); 369 return interf->Tref; 370 } 371 372 /******************************************************************************* 373 * Radiative environment 374 ******************************************************************************/ 375 static double 376 radenv_get_temperature 377 (const struct sdis_radiative_ray* ray, 378 struct sdis_data* data) 379 { 380 (void)ray, (void)data; 381 return TR; /* [K] */ 382 } 383 384 static double 385 radenv_get_reference_temperature 386 (const struct sdis_radiative_ray* ray, 387 struct sdis_data* data) 388 { 389 (void)ray, (void)data; 390 return TR; /* [K] */ 391 } 392 393 static struct sdis_radiative_env* 394 create_radenv(struct sdis_device* sdis) 395 { 396 struct sdis_radiative_env_shader shader = SDIS_RADIATIVE_ENV_SHADER_NULL; 397 struct sdis_radiative_env* radenv = NULL; 398 399 shader.temperature = radenv_get_temperature; 400 shader.reference_temperature = radenv_get_reference_temperature; 401 OK(sdis_radiative_env_create(sdis, &shader, NULL, &radenv)); 402 return radenv; 403 } 404 405 /******************************************************************************* 406 * Helper functions 407 ******************************************************************************/ 408 static void 409 create_interface 410 (struct sdis_device* dev, 411 struct sdis_medium* front, 412 struct sdis_medium* back, 413 const struct interf* interf, 414 struct sdis_interface** out_interf) 415 { 416 struct sdis_interface_shader shader = SDIS_INTERFACE_SHADER_NULL; 417 struct sdis_data* data = NULL; 418 419 CHK(interf != NULL); 420 421 shader.front.temperature = interface_get_temperature; 422 shader.back.temperature = interface_get_temperature; 423 if(sdis_medium_get_type(front) != sdis_medium_get_type(back)) { 424 shader.convection_coef = interface_get_convection_coef; 425 shader.convection_coef_upper_bound = interf->h; 426 } 427 if(sdis_medium_get_type(front) == SDIS_FLUID) { 428 shader.front.emissivity = interface_get_emissivity; 429 shader.front.reference_temperature = interface_get_Tref; 430 } 431 if(sdis_medium_get_type(back) == SDIS_FLUID) { 432 shader.back.emissivity = interface_get_emissivity; 433 shader.back.reference_temperature = interface_get_Tref; 434 } 435 436 OK(sdis_data_create(dev, sizeof(struct interf), ALIGNOF(struct interf), 437 NULL, &data)); 438 *((struct interf*)sdis_data_get(data)) = *interf; 439 440 OK(sdis_interface_create(dev, front, back, &shader, data, out_interf)); 441 OK(sdis_data_ref_put(data)); 442 } 443 444 static void 445 solve_tbound1 446 (struct sdis_scene* scn, 447 struct ssp_rng* rng) 448 { 449 char dump[128]; 450 struct time t0, t1; 451 struct sdis_estimator* estimator; 452 struct sdis_solve_probe_boundary_args solve_args 453 = SDIS_SOLVE_PROBE_BOUNDARY_ARGS_DEFAULT; 454 struct sdis_mc T = SDIS_MC_NULL; 455 struct sdis_mc time = SDIS_MC_NULL; 456 size_t nreals; 457 size_t nfails; 458 enum sdis_scene_dimension dim; 459 const double t[] = { 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000 }; 460 const double ref[sizeof(t) / sizeof(*t)] = { 461 290.046375, 289.903935, 289.840490, 289.802690, 289.777215, 289.759034, 462 289.745710, 289.735826, 289.728448, 289.722921 463 }; 464 const int nsimuls = sizeof(t) / sizeof(*t); 465 int isimul; 466 ASSERT(scn && rng); 467 468 OK(sdis_scene_get_dimension(scn, &dim)); 469 470 solve_args.nrealisations = N; 471 solve_args.side = SDIS_FRONT; 472 FOR_EACH(isimul, 0, nsimuls) { 473 solve_args.time_range[0] = solve_args.time_range[1] = t[isimul]; 474 if(dim == SDIS_SCENE_2D) { 475 solve_args.iprim = model2d_nsegments - 1; 476 solve_args.uv[0] = ssp_rng_uniform_double(rng, 0, 1); 477 } else { 478 double u = ssp_rng_uniform_double(rng, 0, 1); 479 double v = ssp_rng_uniform_double(rng, 0, 1); 480 double w = ssp_rng_uniform_double(rng, 0, 1); 481 double x = 1 / (u + v + w); 482 solve_args.iprim = (isimul % 2) ? 10 : 11; /* +X face */ 483 solve_args.uv[0] = u * x; 484 solve_args.uv[1] = v * x; 485 } 486 487 time_current(&t0); 488 OK(sdis_solve_probe_boundary(scn, &solve_args, &estimator)); 489 time_sub(&t0, time_current(&t1), &t0); 490 time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); 491 492 OK(sdis_estimator_get_realisation_count(estimator, &nreals)); 493 OK(sdis_estimator_get_failure_count(estimator, &nfails)); 494 OK(sdis_estimator_get_temperature(estimator, &T)); 495 OK(sdis_estimator_get_realisation_time(estimator, &time)); 496 497 switch(dim) { 498 case SDIS_SCENE_2D: 499 printf("Unstationary temperature at (%lu/%g) at t=%g = %g ~ %g +/- %g\n", 500 (unsigned long)solve_args.iprim, solve_args.uv[0], t[isimul], 501 ref[isimul], T.E, T.SE); 502 break; 503 case SDIS_SCENE_3D: 504 printf("Unstationary temperature at (%lu/%g,%g) at t=%g = %g ~ %g +/- %g\n", 505 (unsigned long)solve_args.iprim, SPLIT2(solve_args.uv), t[isimul], 506 ref[isimul], T.E, T.SE); 507 break; 508 default: FATAL("Unreachable code.\n"); break; 509 } 510 printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); 511 printf("Elapsed time = %s\n", dump); 512 printf("Time per realisation (in usec) = %g +/- %g\n\n", time.E, time.SE); 513 514 CHK(eq_eps(T.E, ref[isimul], EPS)); 515 /*CHK(eq_eps(T.E, ref[isimul], T.SE*3));*/ 516 517 OK(sdis_estimator_ref_put(estimator)); 518 } 519 } 520 521 static void 522 solve_tbound2 523 (struct sdis_scene* scn, 524 struct ssp_rng* rng) 525 { 526 char dump[128]; 527 struct time t0, t1; 528 struct sdis_estimator* estimator; 529 struct sdis_solve_probe_boundary_args solve_args 530 = SDIS_SOLVE_PROBE_BOUNDARY_ARGS_DEFAULT; 531 struct sdis_mc T = SDIS_MC_NULL; 532 struct sdis_mc time = SDIS_MC_NULL; 533 size_t nreals; 534 size_t nfails; 535 enum sdis_scene_dimension dim; 536 const double t[] = { 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000 }; 537 const double ref[sizeof(t) / sizeof(*t)] = { 538 309.08032, 309.34626, 309.46525, 309.53625, 309.58408, 309.618121, 539 309.642928, 309.661167, 309.674614, 309.684524 540 }; 541 const int nsimuls = sizeof(t) / sizeof(*t); 542 int isimul; 543 ASSERT(scn && rng); 544 545 OK(sdis_scene_get_dimension(scn, &dim)); 546 547 solve_args.nrealisations = N; 548 solve_args.side = SDIS_FRONT; 549 FOR_EACH(isimul, 0, nsimuls) { 550 solve_args.time_range[0] = solve_args.time_range[1] = t[isimul]; 551 if(dim == SDIS_SCENE_2D) { 552 solve_args.iprim = model2d_nsegments - 1; 553 solve_args.uv[0] = ssp_rng_uniform_double(rng, 0, 1); 554 } else { 555 double u = ssp_rng_uniform_double(rng, 0, 1); 556 double v = ssp_rng_uniform_double(rng, 0, 1); 557 double w = ssp_rng_uniform_double(rng, 0, 1); 558 double x = 1 / (u + v + w); 559 solve_args.iprim = (isimul % 2) ? 20 : 21; /* Internal face */ 560 solve_args.uv[0] = u * x; 561 solve_args.uv[1] = v * x; 562 } 563 564 time_current(&t0); 565 OK(sdis_solve_probe_boundary(scn, &solve_args, &estimator)); 566 time_sub(&t0, time_current(&t1), &t0); 567 time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); 568 569 OK(sdis_estimator_get_realisation_count(estimator, &nreals)); 570 OK(sdis_estimator_get_failure_count(estimator, &nfails)); 571 OK(sdis_estimator_get_temperature(estimator, &T)); 572 OK(sdis_estimator_get_realisation_time(estimator, &time)); 573 574 switch(dim) { 575 case SDIS_SCENE_2D: 576 printf("Unstationary temperature at (%lu/%g) at t=%g = %g ~ %g +/- %g\n", 577 (unsigned long)solve_args.iprim, solve_args.uv[0], t[isimul], 578 ref[isimul], T.E, T.SE); 579 break; 580 case SDIS_SCENE_3D: 581 printf("Unstationary temperature at (%lu/%g,%g) at t=%g = %g ~ %g +/- %g\n", 582 (unsigned long)solve_args.iprim, SPLIT2(solve_args.uv), t[isimul], 583 ref[isimul], T.E, T.SE); 584 break; 585 default: FATAL("Unreachable code.\n"); break; 586 } 587 printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); 588 printf("Elapsed time = %s\n", dump); 589 printf("Time per realisation (in usec) = %g +/- %g\n\n", time.E, time.SE); 590 591 CHK(nfails + nreals == N); 592 CHK(nfails <= N/1000); 593 CHK(eq_eps(T.E, ref[isimul], EPS)); 594 /*CHK(eq_eps(T.E, ref[isimul], T.SE*3));*/ 595 596 OK(sdis_estimator_ref_put(estimator)); 597 } 598 } 599 600 static void 601 solve_tsolid 602 (struct sdis_scene* scn, 603 struct ssp_rng* rng) 604 { 605 char dump[128]; 606 struct time t0, t1; 607 struct sdis_estimator* estimator; 608 struct sdis_solve_probe_args solve_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT; 609 struct sdis_mc T = SDIS_MC_NULL; 610 struct sdis_mc time = SDIS_MC_NULL; 611 size_t nreals; 612 size_t nfails; 613 enum sdis_scene_dimension dim; 614 const double t[] = { 0, 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000 }; 615 const double ref[sizeof(t) / sizeof(*t)] = { 616 300, 300.87408, 302.25832, 303.22164, 303.89954, 304.39030, 304.75041, 617 305.01595, 305.21193, 305.35641, 305.46271 618 }; 619 const int nsimuls = sizeof(t) / sizeof(*t); 620 int isimul; 621 ASSERT(scn && rng); 622 623 OK(sdis_scene_get_dimension(scn, &dim)); 624 625 solve_args.nrealisations = N; 626 FOR_EACH(isimul, 0, nsimuls) { 627 solve_args.time_range[0] = solve_args.time_range[1] = t[isimul]; 628 solve_args.position[0] = X_PROBE; 629 solve_args.position[1] = ssp_rng_uniform_double(rng, 0.1*XHpE, 0.9*XHpE); 630 631 if(dim == SDIS_SCENE_3D) 632 solve_args.position[2] = ssp_rng_uniform_double(rng, 0.1*XHpE, 0.9*XHpE); 633 634 time_current(&t0); 635 OK(sdis_solve_probe(scn, &solve_args, &estimator)); 636 time_sub(&t0, time_current(&t1), &t0); 637 time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); 638 639 OK(sdis_estimator_get_realisation_count(estimator, &nreals)); 640 OK(sdis_estimator_get_failure_count(estimator, &nfails)); 641 OK(sdis_estimator_get_temperature(estimator, &T)); 642 OK(sdis_estimator_get_realisation_time(estimator, &time)); 643 644 switch (dim) { 645 case SDIS_SCENE_2D: 646 printf("Unstationary temperature at (%g,%g) at t=%g = %g ~ %g +/- %g\n", 647 SPLIT2(solve_args.position), t[isimul], ref[isimul], T.E, T.SE); 648 break; 649 case SDIS_SCENE_3D: 650 printf("Unstationary temperature at (%g,%g,%g) at t=%g = %g ~ %g +/- %g\n", 651 SPLIT3(solve_args.position), t[isimul], ref[isimul], T.E, T.SE); 652 break; 653 default: FATAL("Unreachable code.\n"); break; 654 } 655 printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); 656 printf("Elapsed time = %s\n", dump); 657 printf("Time per realisation (in usec) = %g +/- %g\n\n", time.E, time.SE); 658 659 CHK(nfails + nreals == N); 660 CHK(nfails <= N / 1000); 661 CHK(eq_eps(T.E, ref[isimul], EPS)); 662 /*CHK(eq_eps(T.E, ref[isimul], T.SE*3));*/ 663 664 OK(sdis_estimator_ref_put(estimator)); 665 } 666 } 667 668 static void 669 solve_tfluid 670 (struct sdis_scene* scn) 671 { 672 char dump[128]; 673 struct time t0, t1; 674 struct sdis_estimator* estimator; 675 struct sdis_solve_probe_args solve_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT; 676 struct sdis_mc T = SDIS_MC_NULL; 677 struct sdis_mc time = SDIS_MC_NULL; 678 size_t nreals; 679 size_t nfails; 680 enum sdis_scene_dimension dim; 681 double eps; 682 const double t[] = { 0, 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000 }; 683 const double ref[sizeof(t) / sizeof(*t)] = { 684 300, 309.53905, 309.67273, 309.73241, 309.76798, 309.79194, 309.80899, 685 309.82141, 309.83055, 309.83728, 309.84224 686 }; 687 const int nsimuls = sizeof(t) / sizeof(*t); 688 int isimul; 689 ASSERT(scn); 690 691 OK(sdis_scene_get_dimension(scn, &dim)); 692 693 solve_args.nrealisations = N; 694 solve_args.position[0] = XH * 0.5; 695 solve_args.position[1] = XH * 0.5; 696 solve_args.position[2] = XH * 0.5; 697 FOR_EACH(isimul, 0, nsimuls) { 698 solve_args.time_range[0] = solve_args.time_range[1] = t[isimul]; 699 700 time_current(&t0); 701 OK(sdis_solve_probe(scn, &solve_args, &estimator)); 702 time_sub(&t0, time_current(&t1), &t0); 703 time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); 704 705 OK(sdis_estimator_get_realisation_count(estimator, &nreals)); 706 OK(sdis_estimator_get_failure_count(estimator, &nfails)); 707 OK(sdis_estimator_get_temperature(estimator, &T)); 708 OK(sdis_estimator_get_realisation_time(estimator, &time)); 709 710 switch (dim) { 711 case SDIS_SCENE_2D: 712 printf("Unstationary fluid temperature at t=%g = %g ~ %g +/- %g\n", 713 t[isimul], ref[isimul], T.E, T.SE); 714 break; 715 case SDIS_SCENE_3D: 716 printf("Unstationary fluid temperature at t=%g = %g ~ %g +/- %g\n", 717 t[isimul], ref[isimul], T.E, T.SE); 718 break; 719 default: FATAL("Unreachable code.\n"); break; 720 } 721 printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); 722 printf("Elapsed time = %s\n", dump); 723 printf("Time per realisation (in usec) = %g +/- %g\n\n", time.E, time.SE); 724 725 CHK(nfails + nreals == N); 726 CHK(nfails <= N / 1000); 727 728 eps = EPS; 729 CHK(eq_eps(T.E, ref[isimul], eps)); 730 731 OK(sdis_estimator_ref_put(estimator)); 732 } 733 } 734 735 /******************************************************************************* 736 * Test 737 ******************************************************************************/ 738 int 739 main(int argc, char** argv) 740 { 741 struct sdis_data* data = NULL; 742 struct sdis_device* dev = NULL; 743 struct sdis_medium* fluid = NULL; 744 struct sdis_medium* fluid_A = NULL; 745 struct sdis_medium* solid = NULL; 746 struct sdis_medium* dummy_solid = NULL; 747 struct sdis_interface* interf_adiabatic_1 = NULL; 748 struct sdis_interface* interf_adiabatic_2 = NULL; 749 struct sdis_interface* interf_TG = NULL; 750 struct sdis_interface* interf_P = NULL; 751 struct sdis_interface* interf_TA = NULL; 752 struct sdis_radiative_env* radenv = NULL; 753 struct sdis_scene* box_scn = NULL; 754 struct sdis_scene* square_scn = NULL; 755 struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT; 756 struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER; 757 struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER; 758 struct sdis_interface* model3d_interfaces[22 /*#triangles*/]; 759 struct sdis_interface* model2d_interfaces[7/*#segments*/]; 760 struct interf interf_props; 761 struct solid* solid_props = NULL; 762 struct fluid* fluid_props = NULL; 763 struct ssp_rng* rng = NULL; 764 (void)argc, (void)argv; 765 766 OK(sdis_device_create(&SDIS_DEVICE_CREATE_ARGS_DEFAULT, &dev)); 767 768 radenv = create_radenv(dev); 769 770 /* Setup the solid shader */ 771 solid_shader.calorific_capacity = solid_get_calorific_capacity; 772 solid_shader.thermal_conductivity = solid_get_thermal_conductivity; 773 solid_shader.volumic_mass = solid_get_volumic_mass; 774 solid_shader.delta = solid_get_delta; 775 solid_shader.temperature = solid_get_temperature; 776 777 /* Create the solid media */ 778 OK(sdis_data_create(dev, sizeof(struct solid), 16, NULL, &data)); 779 solid_props = sdis_data_get(data); 780 solid_props->lambda = LAMBDA; 781 solid_props->cp = CP_S; 782 solid_props->rho = RHO_S; 783 solid_props->delta = DELTA; 784 solid_props->t0 = 0; 785 solid_props->temperature = T0_SOLID; 786 OK(sdis_solid_create(dev, &solid_shader, data, &solid)); 787 OK(sdis_data_ref_put(data)); 788 789 /* Create a dummy solid media to be used outside the model */ 790 OK(sdis_data_create(dev, sizeof(struct solid), 16, NULL, &data)); 791 solid_props = sdis_data_get(data); 792 solid_props->lambda = 0; 793 solid_props->cp = 1; 794 solid_props->rho = 1; 795 solid_props->delta = 1; 796 solid_props->t0 = INF; 797 solid_props->temperature = SDIS_TEMPERATURE_NONE; 798 OK(sdis_solid_create(dev, &solid_shader, data, &dummy_solid)); 799 OK(sdis_data_ref_put(data)); 800 801 /* Setup the fluid shader */ 802 fluid_shader.calorific_capacity = fluid_get_calorific_capacity; 803 fluid_shader.volumic_mass = fluid_get_volumic_mass; 804 fluid_shader.temperature = fluid_get_temperature; 805 806 /* Create the internal fluid media */ 807 OK(sdis_data_create(dev, sizeof(struct fluid), 16, NULL, &data)); 808 fluid_props = sdis_data_get(data); 809 fluid_props->cp = CP_F; 810 fluid_props->rho = RHO_F; 811 fluid_props->t0 = 0; 812 fluid_props->temperature = T0_FLUID; 813 OK(sdis_fluid_create(dev, &fluid_shader, data, &fluid)); 814 OK(sdis_data_ref_put(data)); 815 816 /* Create the 'A' fluid media */ 817 OK(sdis_data_create(dev, sizeof(struct fluid), 16, NULL, &data)); 818 fluid_props = sdis_data_get(data); 819 fluid_props->cp = 1; 820 fluid_props->rho = 1; 821 fluid_props->t0 = INF; 822 fluid_props->temperature = TA; 823 OK(sdis_fluid_create(dev, &fluid_shader, data, &fluid_A)); 824 OK(sdis_data_ref_put(data)); 825 826 /* Create the adiabatic interfaces */ 827 interf_props.temperature = SDIS_TEMPERATURE_NONE; 828 interf_props.h = 0; 829 interf_props.emissivity = 0; 830 interf_props.Tref = TREF; 831 create_interface(dev, fluid, dummy_solid, &interf_props, &interf_adiabatic_1); 832 create_interface(dev, solid, dummy_solid, &interf_props, &interf_adiabatic_2); 833 834 /* Create the P interface */ 835 interf_props.temperature = SDIS_TEMPERATURE_NONE; 836 interf_props.h = HC; 837 interf_props.emissivity = 1; 838 interf_props.Tref = TREF; 839 create_interface(dev, fluid, solid, &interf_props, &interf_P); 840 841 /* Create the TG interface */ 842 interf_props.temperature = TG; 843 interf_props.h = HG; 844 interf_props.emissivity = 1; 845 interf_props.Tref = TG; 846 create_interface(dev, fluid, dummy_solid, &interf_props, &interf_TG); 847 848 /* Create the TA interface */ 849 interf_props.temperature = SDIS_TEMPERATURE_NONE; 850 interf_props.h = HA; 851 interf_props.emissivity = 1; 852 interf_props.Tref = TREF; 853 create_interface(dev, solid, fluid_A, &interf_props, &interf_TA); 854 855 /* Release the media */ 856 OK(sdis_medium_ref_put(solid)); 857 OK(sdis_medium_ref_put(dummy_solid)); 858 OK(sdis_medium_ref_put(fluid)); 859 OK(sdis_medium_ref_put(fluid_A)); 860 861 /* Front */ 862 model3d_interfaces[0] = interf_adiabatic_1; 863 model3d_interfaces[1] = interf_adiabatic_1; 864 model3d_interfaces[2] = interf_adiabatic_2; 865 model3d_interfaces[3] = interf_adiabatic_2; 866 /* Left */ 867 model3d_interfaces[4] = interf_TG; 868 model3d_interfaces[5] = interf_TG; 869 /* Back */ 870 model3d_interfaces[6] = interf_adiabatic_1; 871 model3d_interfaces[7] = interf_adiabatic_1; 872 model3d_interfaces[8] = interf_adiabatic_2; 873 model3d_interfaces[9] = interf_adiabatic_2; 874 /* Right */ 875 model3d_interfaces[10] = interf_TA; 876 model3d_interfaces[11] = interf_TA; 877 /* Top */ 878 model3d_interfaces[12] = interf_adiabatic_1; 879 model3d_interfaces[13] = interf_adiabatic_1; 880 model3d_interfaces[14] = interf_adiabatic_2; 881 model3d_interfaces[15] = interf_adiabatic_2; 882 /* Bottom */ 883 model3d_interfaces[16] = interf_adiabatic_1; 884 model3d_interfaces[17] = interf_adiabatic_1; 885 model3d_interfaces[18] = interf_adiabatic_2; 886 model3d_interfaces[19] = interf_adiabatic_2; 887 /* Inside */ 888 model3d_interfaces[20] = interf_P; 889 model3d_interfaces[21] = interf_P; 890 891 /* Bottom */ 892 model2d_interfaces[0] = interf_adiabatic_2; 893 model2d_interfaces[1] = interf_adiabatic_1; 894 /* Left */ 895 model2d_interfaces[2] = interf_TG; 896 /* Top */ 897 model2d_interfaces[3] = interf_adiabatic_1; 898 model2d_interfaces[4] = interf_adiabatic_2; 899 /* Right */ 900 model2d_interfaces[5] = interf_TA; 901 /* Contact */ 902 model2d_interfaces[6] = interf_P; 903 904 /* Create the box scene */ 905 scn_args.get_indices = model3d_get_indices; 906 scn_args.get_interface = model3d_get_interface; 907 scn_args.get_position = model3d_get_position; 908 scn_args.nprimitives = model3d_ntriangles; 909 scn_args.nvertices = model3d_nvertices; 910 scn_args.context = model3d_interfaces; 911 scn_args.radenv = radenv; 912 scn_args.t_range[0] = MMIN(MMIN(MMIN(MMIN(T0_FLUID, T0_SOLID), TA), TG), TR); 913 scn_args.t_range[1] = MMAX(MMAX(MMAX(MMAX(T0_FLUID, T0_SOLID), TA), TG), TR); 914 OK(sdis_scene_create(dev, &scn_args, &box_scn)); 915 916 /* Create the square scene */ 917 scn_args.get_indices = model2d_get_indices; 918 scn_args.get_interface = model2d_get_interface; 919 scn_args.get_position = model2d_get_position; 920 scn_args.nprimitives = model2d_nsegments; 921 scn_args.nvertices = model2d_nvertices; 922 scn_args.context = model2d_interfaces; 923 scn_args.radenv = radenv; 924 scn_args.t_range[0] = MMIN(MMIN(MMIN(MMIN(T0_FLUID, T0_SOLID), TA), TG), TR); 925 scn_args.t_range[1] = MMAX(MMAX(MMAX(MMAX(T0_FLUID, T0_SOLID), TA), TG), TR); 926 OK(sdis_scene_2d_create(dev, &scn_args, &square_scn)); 927 928 /* Release the interfaces */ 929 OK(sdis_interface_ref_put(interf_adiabatic_1)); 930 OK(sdis_interface_ref_put(interf_adiabatic_2)); 931 OK(sdis_interface_ref_put(interf_TG)); 932 OK(sdis_interface_ref_put(interf_P)); 933 OK(sdis_interface_ref_put(interf_TA)); 934 935 /* Solve */ 936 OK(ssp_rng_create(NULL, SSP_RNG_KISS, &rng)); 937 printf(">> Box scene\n"); 938 solve_tfluid(box_scn); 939 solve_tbound1(box_scn, rng); 940 solve_tbound2(box_scn, rng); 941 solve_tsolid(box_scn, rng); 942 printf("\n>> Square scene\n"); 943 solve_tfluid(square_scn); 944 solve_tbound1(box_scn, rng); 945 solve_tbound2(box_scn, rng); 946 solve_tsolid(square_scn, rng); 947 948 OK(sdis_radiative_env_ref_put(radenv)); 949 OK(sdis_scene_ref_put(box_scn)); 950 OK(sdis_scene_ref_put(square_scn)); 951 OK(sdis_device_ref_put(dev)); 952 OK(ssp_rng_ref_put(rng)); 953 954 CHK(mem_allocated_size() == 0); 955 return 0; 956 }