test_sdis_picard.c (26540B)
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 <string.h> 20 21 #define N 10000 22 23 /* This test consists in solving the stationary temperature profile in a solid 24 * slab surrounded by two different radiative temperatures (left / right). The 25 * conductivity of the solid material is known, as well as its thickness and 26 * the source term (volumic power density). 27 * 28 * The purpose is to test the Picard radiative transfer algorithm, that can be 29 * compared with analytic results. This algorithm can use a possibly 30 * non-uniform reference temperature field. When the reference temperature 31 * field is uniform and the picard order set to 1, results should be identical 32 * to the classical Monte-Carlo algorithm (using a linearized radiative 33 * transfer scheme). 34 * 35 * Y 36 * | (0.1,1) 37 * o--- X +----------+------+ (1.1,1) 38 * |##########| | 39 * |##########| | 40 * 280K E=1|##########| E=1 | 350K 41 * |##########| | 42 * |##########| | 43 * (-1,-1) +----------+------+ 44 * (0,-1) 45 * 46 * 47 * 48 * Y (0.1, 1, 1) 49 * | +----------+------+ (1.1,1,1) 50 * o--- X /##########/' /| 51 * / +----------+------+ | 52 * Z |##########|*' | | 350K 53 * |##########|*' | | 54 * 280K E=1|##########|*'E=1 | | 55 * |##########|*+....|.+ 56 * |##########|/ |/ 57 * (-1,-1,-1) +----------+------+ 58 * (0,-1,-1) 59 * 60 * lambda = 1.15 W/(m.K) 61 * rho = 1000 kg.m^-3 62 * cp = 800 J/(kg.K) 63 * emissivity = 1 64 * 65 * Basic Tref = 300 K 66 * probe = 0.05 0 0 m 67 * (power = 1000 W.m^-3) */ 68 69 enum interface_type { 70 ADIABATIC, 71 SOLID_FLUID_mX, 72 SOLID_FLUID_pX, 73 BOUNDARY_pX, 74 INTERFACES_COUNT__ 75 }; 76 77 /******************************************************************************* 78 * Geometry 79 ******************************************************************************/ 80 struct geometry { 81 const double* positions; 82 const size_t* indices; 83 struct sdis_interface** interfaces; 84 }; 85 86 static const double vertices_2d[6/*#vertices*/*2/*#coords par vertex*/] = { 87 0.1, -1.0, 88 0.0, -1.0, 89 0.0, 1.0, 90 0.1, 1.0, 91 1.1, -1.0, 92 1.1, 1.0 93 }; 94 static const size_t nvertices_2d = sizeof(vertices_2d) / (sizeof(double)*2); 95 96 static const size_t indices_2d[7/*#segments*/*2/*#indices per segment*/] = { 97 0, 1, /* Solid -Y */ 98 1, 2, /* Solid -X */ 99 2, 3, /* Solid +Y */ 100 3, 0, /* Solid +X */ 101 102 4, 0, /* Right fluid -Y */ 103 3, 5, /* Right fluid +Y */ 104 5, 4 /* Right fluid +X */ 105 }; 106 static const size_t nprimitives_2d = sizeof(indices_2d) / (sizeof(size_t)*2); 107 108 static const double vertices_3d[12/*#vertices*/*3/*#coords per vertex*/] = { 109 0.0,-1.0,-1.0, 110 0.1,-1.0,-1.0, 111 0.0, 1.0,-1.0, 112 0.1, 1.0,-1.0, 113 0.0,-1.0, 1.0, 114 0.1,-1.0, 1.0, 115 0.0, 1.0, 1.0, 116 0.1, 1.0, 1.0, 117 1.1,-1.0,-1.0, 118 1.1, 1.0,-1.0, 119 1.1,-1.0, 1.0, 120 1.1, 1.0, 1.0 121 }; 122 static const size_t nvertices_3d = sizeof(vertices_3d) / (sizeof(double)*3); 123 124 static const size_t indices_3d[22/*#triangles*/*3/*#indices per triangle*/] = { 125 0, 2, 1, 1, 2, 3, /* Solid -Z */ 126 0, 4, 2, 2, 4, 6, /* Solid -X */ 127 4, 5, 6, 6, 5, 7, /* Solid +Z */ 128 3, 7, 1, 1, 7, 5, /* Solid +X */ 129 2, 6, 3, 3, 6, 7, /* Solid +Y */ 130 0, 1, 4, 4, 1, 5, /* Solid -Y */ 131 132 1, 3, 8, 8, 3, 9, /* Right fluid -Z */ 133 5, 10, 7, 7, 10, 11, /* Right fluid +Z */ 134 9, 11, 8, 8, 11, 10, /* Right fluid +X */ 135 3, 7, 9, 9, 7, 11, /* Right fluid +Y */ 136 1, 8, 5, 5, 8, 10 /* Right fluid -Y */ 137 }; 138 static const size_t nprimitives_3d = sizeof(indices_3d) / (sizeof(size_t)*3); 139 140 static void 141 get_indices_2d(const size_t iseg, size_t ids[2], void* ctx) 142 { 143 struct geometry* geom = ctx; 144 CHK(ctx != NULL); 145 ids[0] = geom->indices[iseg*2+0]; 146 ids[1] = geom->indices[iseg*2+1]; 147 } 148 149 static void 150 get_indices_3d(const size_t itri, size_t ids[3], void* ctx) 151 { 152 struct geometry* geom = ctx; 153 CHK(ctx != NULL); 154 ids[0] = geom->indices[itri*3+0]; 155 ids[1] = geom->indices[itri*3+1]; 156 ids[2] = geom->indices[itri*3+2]; 157 } 158 159 static void 160 get_position_2d(const size_t ivert, double pos[2], void* ctx) 161 { 162 struct geometry* geom = ctx; 163 CHK(ctx != NULL); 164 pos[0] = geom->positions[ivert*2+0]; 165 pos[1] = geom->positions[ivert*2+1]; 166 } 167 168 static void 169 get_position_3d(const size_t ivert, double pos[3], void* ctx) 170 { 171 struct geometry* geom = ctx; 172 CHK(ctx != NULL); 173 pos[0] = geom->positions[ivert*3+0]; 174 pos[1] = geom->positions[ivert*3+1]; 175 pos[2] = geom->positions[ivert*3+2]; 176 } 177 178 static void 179 get_interface(const size_t iprim, struct sdis_interface** bound, void* ctx) 180 { 181 struct geometry* geom = ctx; 182 CHK(ctx != NULL); 183 *bound = geom->interfaces[iprim]; 184 } 185 186 /******************************************************************************* 187 * media 188 ******************************************************************************/ 189 struct solid { 190 double lambda; 191 double rho; 192 double cp; 193 double volumic_power; 194 }; 195 196 static double 197 solid_get_calorific_capacity 198 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 199 { 200 const struct solid* solid = sdis_data_cget(data); 201 CHK(vtx && solid); 202 return solid->cp; 203 } 204 205 static double 206 solid_get_thermal_conductivity 207 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 208 { 209 const struct solid* solid = sdis_data_cget(data); 210 CHK(vtx && solid); 211 return solid->lambda; 212 } 213 214 static double 215 solid_get_volumic_mass 216 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 217 { 218 const struct solid* solid = sdis_data_cget(data); 219 CHK(vtx && solid); 220 return solid->rho; 221 } 222 223 static double 224 solid_get_delta 225 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 226 { 227 CHK(vtx && data); 228 return 0.0025; 229 } 230 231 static double 232 solid_get_temperature 233 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 234 { 235 CHK(vtx && data); 236 return SDIS_TEMPERATURE_NONE; 237 } 238 239 static double 240 solid_get_volumic_power 241 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 242 { 243 const struct solid* solid = sdis_data_cget(data); 244 CHK(vtx && solid); 245 return solid->volumic_power; 246 } 247 248 static void 249 create_solid 250 (struct sdis_device* dev, 251 const struct solid* solid_props, 252 struct sdis_medium** solid) 253 { 254 struct sdis_data* data = NULL; 255 struct sdis_solid_shader shader = SDIS_SOLID_SHADER_NULL; 256 CHK(dev && solid_props && solid); 257 258 OK(sdis_data_create 259 (dev, sizeof(struct solid), ALIGNOF(struct solid), NULL, &data)); 260 memcpy(sdis_data_get(data), solid_props, sizeof(struct solid)); 261 shader.calorific_capacity = solid_get_calorific_capacity; 262 shader.thermal_conductivity = solid_get_thermal_conductivity; 263 shader.volumic_mass = solid_get_volumic_mass; 264 shader.delta = solid_get_delta; 265 shader.temperature = solid_get_temperature; 266 shader.volumic_power = solid_get_volumic_power; 267 OK(sdis_solid_create(dev, &shader, data, solid)); 268 OK(sdis_data_ref_put(data)); 269 } 270 271 static void 272 create_fluid(struct sdis_device* dev, struct sdis_medium** fluid) 273 { 274 struct sdis_fluid_shader shader = DUMMY_FLUID_SHADER; 275 OK(sdis_fluid_create(dev, &shader, NULL, fluid)); 276 } 277 278 /******************************************************************************* 279 * Interface 280 ******************************************************************************/ 281 struct interf { 282 double temperature; 283 double h; 284 double emissivity; 285 double specular_fraction; 286 double Tref; 287 }; 288 289 static double 290 interface_get_temperature 291 (const struct sdis_interface_fragment* frag, struct sdis_data* data) 292 { 293 const struct interf* interf = sdis_data_cget(data); 294 CHK(frag && interf); 295 return interf->temperature; 296 } 297 298 static double 299 interface_get_convection_coef 300 (const struct sdis_interface_fragment* frag, struct sdis_data* data) 301 { 302 const struct interf* interf = sdis_data_cget(data); 303 CHK(frag && interf); 304 return interf->h; 305 } 306 307 static double 308 interface_get_emissivity 309 (const struct sdis_interface_fragment* frag, 310 const unsigned source_id, 311 struct sdis_data* data) 312 { 313 const struct interf* interf = sdis_data_cget(data); 314 (void)source_id; 315 CHK(frag && interf); 316 return interf->emissivity; 317 } 318 319 static double 320 interface_get_specular_fraction 321 (const struct sdis_interface_fragment* frag, 322 const unsigned source_id, 323 struct sdis_data* data) 324 { 325 const struct interf* interf = sdis_data_cget(data); 326 (void)source_id; 327 CHK(frag && interf); 328 return interf->specular_fraction; 329 } 330 331 static double 332 interface_get_Tref 333 (const struct sdis_interface_fragment* frag, struct sdis_data* data) 334 { 335 const struct interf* interf = sdis_data_cget(data); 336 CHK(frag && interf); 337 return interf->Tref; 338 } 339 340 static void 341 create_interface 342 (struct sdis_device* dev, 343 struct sdis_medium* front, 344 struct sdis_medium* back, 345 const struct interf* interf, 346 struct sdis_interface** out_interf) 347 { 348 struct sdis_interface_shader shader = SDIS_INTERFACE_SHADER_NULL; 349 struct sdis_data* data = NULL; 350 351 CHK(interf != NULL); 352 353 shader.front.temperature = interface_get_temperature; 354 shader.back.temperature = interface_get_temperature; 355 if(sdis_medium_get_type(front) != sdis_medium_get_type(back)) { 356 shader.convection_coef = interface_get_convection_coef; 357 } 358 if(sdis_medium_get_type(front) == SDIS_FLUID) { 359 shader.front.emissivity = interface_get_emissivity; 360 shader.front.specular_fraction = interface_get_specular_fraction; 361 shader.front.reference_temperature = interface_get_Tref; 362 } 363 if(sdis_medium_get_type(back) == SDIS_FLUID) { 364 shader.back.emissivity = interface_get_emissivity; 365 shader.back.specular_fraction = interface_get_specular_fraction; 366 shader.back.reference_temperature = interface_get_Tref; 367 } 368 shader.convection_coef_upper_bound = MMAX(0, interf->h); 369 370 OK(sdis_data_create 371 (dev, sizeof(struct interf), ALIGNOF(struct interf), NULL, &data)); 372 memcpy(sdis_data_get(data), interf, sizeof(*interf)); 373 374 OK(sdis_interface_create(dev, front, back, &shader, data, out_interf)); 375 OK(sdis_data_ref_put(data)); 376 } 377 378 /******************************************************************************* 379 * Radiative environment 380 ******************************************************************************/ 381 struct radenv { 382 double temperature; /* [K] */ 383 double reference; /* [K] */ 384 }; 385 386 static double 387 radenv_get_temperature 388 (const struct sdis_radiative_ray* ray, 389 struct sdis_data* data) 390 { 391 (void)ray; 392 return ((const struct radenv*)sdis_data_cget(data))->temperature; 393 } 394 395 static double 396 radenv_get_reference_temperature 397 (const struct sdis_radiative_ray* ray, 398 struct sdis_data* data) 399 { 400 (void)ray; 401 return ((const struct radenv*)sdis_data_cget(data))->reference; 402 } 403 404 static struct sdis_radiative_env* 405 create_radenv(struct sdis_device* dev) 406 { 407 struct sdis_radiative_env_shader shader = SDIS_RADIATIVE_ENV_SHADER_NULL; 408 struct sdis_radiative_env* radenv = NULL; 409 struct sdis_data* data = NULL; 410 struct radenv* env = NULL; 411 412 OK(sdis_data_create(dev, sizeof(struct radenv), ALIGNOF(radenv), NULL, &data)); 413 env = sdis_data_get(data); 414 env->temperature = 300; 415 env->reference = 300; 416 417 shader.temperature = radenv_get_temperature; 418 shader.reference_temperature = radenv_get_reference_temperature; 419 OK(sdis_radiative_env_create(dev, &shader, data, &radenv)); 420 OK(sdis_data_ref_put(data)); 421 return radenv; 422 } 423 424 /******************************************************************************* 425 * Helper functions 426 ******************************************************************************/ 427 struct reference_result { 428 double T; /* Temperature at the center of the solid [K] */ 429 double T1; /* Temperature on the left boundary of the solid [K] */ 430 double T2; /* Temperature on the right boundary of the solid [K] */ 431 }; 432 static const struct reference_result REFERENCE_RESULT_NULL = {0,0,0}; 433 434 static void 435 test_picard 436 (struct sdis_scene* scn, 437 const size_t picard_order, 438 const enum sdis_diffusion_algorithm algo, 439 const struct reference_result* ref) 440 { 441 struct sdis_solve_probe_args probe_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT; 442 struct sdis_solve_boundary_args bound_args = SDIS_SOLVE_BOUNDARY_ARGS_DEFAULT; 443 struct sdis_mc mc = SDIS_MC_NULL; 444 struct sdis_estimator* estimator = NULL; 445 enum sdis_scene_dimension dim; 446 size_t prims[2]; 447 enum sdis_side sides[2]; 448 CHK(scn && ref && picard_order >= 1); 449 450 OK(sdis_scene_get_dimension(scn, &dim)); 451 switch(dim) { 452 case SDIS_SCENE_2D: printf(">>> 2D\n"); break; 453 case SDIS_SCENE_3D: printf(">>> 3D\n"); break; 454 default: FATAL("Unreachable code.\n"); break; 455 } 456 457 probe_args.nrealisations = N; 458 probe_args.position[0] = 0.05; 459 probe_args.position[1] = 0; 460 probe_args.position[2] = 0; 461 probe_args.picard_order = picard_order; 462 probe_args.diff_algo = algo; 463 OK(sdis_solve_probe(scn, &probe_args, &estimator)); 464 OK(sdis_estimator_get_temperature(estimator, &mc)); 465 printf("Temperature at `%g %g %g' = %g ~ %g +/- %g\n", 466 SPLIT3(probe_args.position), ref->T, mc.E, mc.SE); 467 CHK(eq_eps(ref->T, mc.E, mc.SE*3)); 468 OK(sdis_estimator_ref_put(estimator)); 469 470 switch(dim) { 471 case SDIS_SCENE_2D: 472 prims[0] = 1; sides[0] = SDIS_BACK; 473 bound_args.nprimitives = 1; 474 break; 475 case SDIS_SCENE_3D: 476 prims[0] = 2; sides[0] = SDIS_BACK; 477 prims[1] = 3; sides[1] = SDIS_BACK; 478 bound_args.nprimitives = 2; 479 break; 480 default: FATAL("Unreachable code.\n"); break; 481 } 482 bound_args.nrealisations = N; 483 bound_args.primitives = prims; 484 bound_args.sides = sides; 485 bound_args.picard_order = picard_order; 486 OK(sdis_solve_boundary(scn, &bound_args, &estimator)); 487 OK(sdis_estimator_get_temperature(estimator, &mc)); 488 printf("T1 = %g ~ %g +/- %g\n", ref->T1, mc.E, mc.SE); 489 CHK(eq_eps(ref->T1, mc.E, mc.SE*3)); 490 OK(sdis_estimator_ref_put(estimator)); 491 492 switch(dim) { 493 case SDIS_SCENE_2D: 494 prims[0] = 3; sides[0] = SDIS_BACK; 495 bound_args.nprimitives = 1; 496 break; 497 case SDIS_SCENE_3D: 498 prims[0] = 6; sides[0] = SDIS_BACK; 499 prims[1] = 7; sides[1] = SDIS_BACK; 500 bound_args.nprimitives = 2; 501 break; 502 default: FATAL("Unreachable code.\n"); break; 503 } 504 bound_args.nrealisations = N; 505 bound_args.primitives = prims; 506 bound_args.sides = sides; 507 bound_args.picard_order = picard_order; 508 OK(sdis_solve_boundary(scn, &bound_args, &estimator)); 509 OK(sdis_estimator_get_temperature(estimator, &mc)); 510 printf("T2 = %g ~ %g +/- %g\n", ref->T2, mc.E, mc.SE); 511 CHK(eq_eps(ref->T2, mc.E, mc.SE*3)); 512 OK(sdis_estimator_ref_put(estimator)); 513 } 514 515 static void 516 register_heat_paths(struct sdis_scene* scn, const size_t picard_order, FILE* stream) 517 { 518 struct sdis_solve_probe_args probe_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT; 519 struct sdis_estimator* estimator = NULL; 520 CHK(scn && picard_order >= 1 && stream); 521 522 523 probe_args.nrealisations = 10; 524 probe_args.position[0] = 0.05; 525 probe_args.position[1] = 0; 526 probe_args.position[2] = 0; 527 probe_args.picard_order = picard_order; 528 probe_args.register_paths = SDIS_HEAT_PATH_ALL; 529 printf("Register %lu heat paths.\n", probe_args.nrealisations); 530 OK(sdis_solve_probe(scn, &probe_args, &estimator)); 531 dump_heat_paths(stream, estimator); 532 OK(sdis_estimator_ref_put(estimator)); 533 } 534 535 static void 536 create_scene_3d 537 (struct sdis_device* dev, 538 struct sdis_interface* interfaces[INTERFACES_COUNT__], 539 struct sdis_radiative_env* radenv, 540 struct sdis_scene** scn) 541 { 542 struct geometry geom; 543 struct sdis_interface* prim_interfaces[32]; 544 struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT; 545 546 CHK(dev && interfaces && radenv && scn); 547 548 /* Setup the per primitive interface of the solid medium */ 549 prim_interfaces[0] = prim_interfaces[1] = interfaces[ADIABATIC]; 550 prim_interfaces[2] = prim_interfaces[3] = interfaces[SOLID_FLUID_mX]; 551 prim_interfaces[4] = prim_interfaces[5] = interfaces[ADIABATIC]; 552 prim_interfaces[6] = prim_interfaces[7] = interfaces[SOLID_FLUID_pX]; 553 prim_interfaces[8] = prim_interfaces[9] = interfaces[ADIABATIC]; 554 prim_interfaces[10] = prim_interfaces[11] = interfaces[ADIABATIC]; 555 556 /* Setup the per primitive interface for the right fluid */ 557 prim_interfaces[12] = prim_interfaces[13] = interfaces[BOUNDARY_pX]; 558 prim_interfaces[14] = prim_interfaces[15] = interfaces[BOUNDARY_pX]; 559 prim_interfaces[16] = prim_interfaces[17] = interfaces[BOUNDARY_pX]; 560 prim_interfaces[18] = prim_interfaces[19] = interfaces[BOUNDARY_pX]; 561 prim_interfaces[20] = prim_interfaces[21] = interfaces[BOUNDARY_pX]; 562 563 /* Create the scene */ 564 geom.positions = vertices_3d; 565 geom.indices = indices_3d; 566 geom.interfaces = prim_interfaces; 567 scn_args.get_indices = get_indices_3d; 568 scn_args.get_interface = get_interface; 569 scn_args.get_position = get_position_3d; 570 scn_args.nprimitives = nprimitives_3d; 571 scn_args.nvertices = nvertices_3d; 572 scn_args.t_range[0] = 280; 573 scn_args.t_range[1] = 350; 574 scn_args.radenv = radenv; 575 scn_args.context = &geom; 576 OK(sdis_scene_create(dev, &scn_args, scn)); 577 } 578 579 static void 580 create_scene_2d 581 (struct sdis_device* dev, 582 struct sdis_interface* interfaces[INTERFACES_COUNT__], 583 struct sdis_radiative_env* radenv, 584 struct sdis_scene** scn) 585 { 586 struct geometry geom; 587 struct sdis_interface* prim_interfaces[10/*#segment*/]; 588 struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT; 589 590 CHK(dev && interfaces && radenv && scn); 591 592 /* Setup the per primitive interface of the solid medium */ 593 prim_interfaces[0] = interfaces[ADIABATIC]; 594 prim_interfaces[1] = interfaces[SOLID_FLUID_mX]; 595 prim_interfaces[2] = interfaces[ADIABATIC]; 596 prim_interfaces[3] = interfaces[SOLID_FLUID_pX]; 597 598 /* Setup the per primitive interface of the fluid on the right of the medium */ 599 prim_interfaces[4] = interfaces[BOUNDARY_pX]; 600 prim_interfaces[5] = interfaces[BOUNDARY_pX]; 601 prim_interfaces[6] = interfaces[BOUNDARY_pX]; 602 603 /* Create the scene */ 604 geom.positions = vertices_2d; 605 geom.indices = indices_2d; 606 geom.interfaces = prim_interfaces; 607 scn_args.get_indices = get_indices_2d; 608 scn_args.get_interface = get_interface; 609 scn_args.get_position = get_position_2d; 610 scn_args.nprimitives = nprimitives_2d; 611 scn_args.nvertices = nvertices_2d; 612 scn_args.t_range[0] = 280; 613 scn_args.t_range[1] = 350; 614 scn_args.radenv = radenv; 615 scn_args.context = &geom; 616 OK(sdis_scene_2d_create(dev, &scn_args, scn)); 617 } 618 619 /******************************************************************************* 620 * Test 621 ******************************************************************************/ 622 int 623 main(int argc, char** argv) 624 { 625 FILE* stream = NULL; 626 627 struct sdis_device* dev = NULL; 628 struct sdis_radiative_env* radenv = NULL; 629 struct sdis_scene* scn_2d = NULL; 630 struct sdis_scene* scn_3d = NULL; 631 struct sdis_medium* solid = NULL; 632 struct sdis_medium* fluid = NULL; 633 struct sdis_medium* dummy = NULL; 634 struct sdis_interface* interfaces[INTERFACES_COUNT__]; 635 636 struct radenv* radenv_props = NULL; 637 struct solid solid_props; 638 struct solid* psolid_props; 639 struct reference_result ref = REFERENCE_RESULT_NULL; 640 struct interf interf_props; 641 struct interf* pinterf_props[INTERFACES_COUNT__]; 642 643 double t_range[2]; 644 645 size_t i; 646 (void)argc, (void)argv; 647 648 OK(sdis_device_create(&SDIS_DEVICE_CREATE_ARGS_DEFAULT, &dev)); 649 650 radenv = create_radenv(dev); 651 radenv_props = sdis_data_get(sdis_radiative_env_get_data(radenv)); 652 653 /* Solid medium */ 654 solid_props.lambda = 1.15; 655 solid_props.rho = 1000; 656 solid_props.cp = 800; 657 solid_props.volumic_power = SDIS_VOLUMIC_POWER_NONE; 658 create_solid(dev, &solid_props, &solid); 659 660 /* Dummy solid medium */ 661 solid_props.lambda = 0; 662 solid_props.rho = 1000; 663 solid_props.cp = 800; 664 solid_props.volumic_power = SDIS_VOLUMIC_POWER_NONE; 665 create_solid(dev, &solid_props, &dummy); 666 667 /* Fluid medium */ 668 create_fluid(dev, &fluid); 669 670 /* Create the adiabatic interface for the solid */ 671 interf_props.temperature = SDIS_TEMPERATURE_NONE; 672 interf_props.h = -1; 673 interf_props.emissivity = -1; 674 interf_props.specular_fraction = -1; 675 interf_props.Tref = SDIS_TEMPERATURE_NONE; 676 create_interface(dev, solid, dummy, &interf_props, interfaces+ADIABATIC); 677 678 /* Create the interface between the solid and the fluid */ 679 interf_props.temperature = SDIS_TEMPERATURE_NONE; 680 interf_props.h = 0; 681 interf_props.emissivity = 1; 682 interf_props.specular_fraction = 0; 683 interf_props.Tref = 280; 684 create_interface(dev, solid, fluid, &interf_props, interfaces+SOLID_FLUID_mX); 685 interf_props.Tref = 350; 686 create_interface(dev, solid, fluid, &interf_props, interfaces+SOLID_FLUID_pX); 687 688 /* Create the interface for the fluid on the right */ 689 interf_props.temperature = 350; 690 interf_props.h = -1; 691 interf_props.emissivity = 1; 692 interf_props.specular_fraction = 0; 693 interf_props.Tref = 350; 694 create_interface(dev, fluid, dummy, &interf_props, interfaces+BOUNDARY_pX); 695 696 /* Fetch pointers toward the solid and the interfaces */ 697 psolid_props = sdis_data_get(sdis_medium_get_data(solid)); 698 FOR_EACH(i, 0, INTERFACES_COUNT__) { 699 pinterf_props[i] = sdis_data_get(sdis_interface_get_data(interfaces[i])); 700 } 701 702 create_scene_2d(dev, interfaces, radenv, &scn_2d); 703 create_scene_3d(dev, interfaces, radenv, &scn_3d); 704 705 CHK((stream = tmpfile()) != NULL); 706 707 /* Test picard1 with a constant Tref <=> regular linearisation */ 708 printf("Test Picard1 with a constant Tref of 300 K\n"); 709 ref.T = 314.99999999999989; 710 ref.T1 = 307.64122364709766; 711 ref.T2 = 322.35877635290217; 712 pinterf_props[SOLID_FLUID_mX]->Tref = 300; 713 pinterf_props[SOLID_FLUID_pX]->Tref = 300; 714 pinterf_props[BOUNDARY_pX]->Tref = 300; 715 radenv_props->temperature = 280; 716 radenv_props->reference = 300; 717 test_picard(scn_2d, 1/*Picard order*/, SDIS_DIFFUSION_WOS, &ref); 718 test_picard(scn_3d, 1/*Picard order*/, SDIS_DIFFUSION_WOS, &ref); 719 printf("\n"); 720 721 /* Test picard1 using T4 as a reference */ 722 printf("Test Picard1 using T4 as a reference\n"); 723 ref.T = 320.37126474482994; 724 ref.T1 = 312.12650299072266; 725 ref.T2 = 328.61602649893723; 726 pinterf_props[SOLID_FLUID_mX]->Tref = ref.T1; 727 pinterf_props[SOLID_FLUID_pX]->Tref = ref.T2; 728 pinterf_props[BOUNDARY_pX]->Tref = 350; 729 radenv_props->temperature = 280; 730 radenv_props->reference = 280; 731 test_picard(scn_2d, 1/*Picard order*/, SDIS_DIFFUSION_DELTA_SPHERE, &ref); 732 test_picard(scn_3d, 1/*Picard order*/, SDIS_DIFFUSION_DELTA_SPHERE, &ref); 733 printf("\n"); 734 735 /* Test picard2 */ 736 printf("Test Picard2 with a constant Tref of 300K\n"); 737 ref.T = 320.37126474482994; 738 ref.T1 = 312.12650299072266; 739 ref.T2 = 328.61602649893723; 740 pinterf_props[SOLID_FLUID_mX]->Tref = 300; 741 pinterf_props[SOLID_FLUID_pX]->Tref = 300; 742 pinterf_props[BOUNDARY_pX]->Tref = 300; 743 radenv_props->temperature = 280; 744 radenv_props->reference = 300; 745 test_picard(scn_2d, 2/*Picard order*/, SDIS_DIFFUSION_WOS, &ref); 746 test_picard(scn_3d, 2/*Picard order*/, SDIS_DIFFUSION_WOS, &ref); 747 printf("\n"); 748 749 t_range[0] = 200; 750 t_range[1] = 500; 751 752 OK(sdis_scene_set_temperature_range(scn_2d, t_range)); 753 OK(sdis_scene_set_temperature_range(scn_3d, t_range)); 754 755 /* Test picard3 */ 756 printf("Test Picard3 with a delta T of 300K\n"); 757 ref.T = 416.4023; 758 ref.T1 = 372.7557; 759 ref.T2 = 460.0489; 760 pinterf_props[BOUNDARY_pX]->temperature = t_range[1]; 761 pinterf_props[SOLID_FLUID_mX]->Tref = 350; 762 pinterf_props[SOLID_FLUID_pX]->Tref = 450; 763 pinterf_props[BOUNDARY_pX]->Tref = pinterf_props[BOUNDARY_pX]->temperature; 764 radenv_props->temperature = t_range[0]; 765 radenv_props->reference = t_range[0]; 766 test_picard(scn_2d, 3/*Picard order*/, SDIS_DIFFUSION_WOS, &ref); 767 test_picard(scn_3d, 3/*Picard order*/, SDIS_DIFFUSION_WOS, &ref); 768 register_heat_paths(scn_2d, 3/*Picard order*/, stream); 769 register_heat_paths(scn_3d, 3/*Picard order*/, stream); 770 printf("\n"); 771 772 t_range[0] = 280; 773 t_range[1] = 350; 774 OK(sdis_scene_set_temperature_range(scn_2d, t_range)); 775 OK(sdis_scene_set_temperature_range(scn_3d, t_range)); 776 pinterf_props[BOUNDARY_pX]->temperature = t_range[1]; 777 778 /* Add volumic power */ 779 psolid_props->volumic_power = 1000; 780 781 /* Test picard1 with a volumic power and constant Tref */ 782 printf("Test Picard1 with a volumic power of 1000 W/m^3 and a constant Tref " 783 "of 300 K\n"); 784 ref.T = 324.25266420769509; 785 ref.T1 = 315.80693133305368; 786 ref.T2 = 330.52448403885825; 787 pinterf_props[SOLID_FLUID_mX]->Tref = 300; 788 pinterf_props[SOLID_FLUID_pX]->Tref = 300; 789 pinterf_props[BOUNDARY_pX]->Tref = 300; 790 radenv_props->temperature = t_range[0]; 791 radenv_props->reference = 300; 792 test_picard(scn_2d, 1/*Picard order*/, SDIS_DIFFUSION_DELTA_SPHERE, &ref); 793 test_picard(scn_3d, 1/*Picard order*/, SDIS_DIFFUSION_DELTA_SPHERE, &ref); 794 printf("\n"); 795 796 /* Test picard1 with a volumic power and T4 a the reference */ 797 printf("Test Picard1 with a volumic power of 1000 W/m^3 and T4 as a reference\n"); 798 ref.T = 327.95981050850446; 799 ref.T1 = 318.75148773193359; 800 ref.T2 = 334.99422024159708; 801 pinterf_props[SOLID_FLUID_mX]->Tref = ref.T1; 802 pinterf_props[SOLID_FLUID_pX]->Tref = ref.T2; 803 pinterf_props[BOUNDARY_pX]->Tref = 350; 804 radenv_props->temperature = 280; 805 radenv_props->reference = 280; 806 test_picard(scn_2d, 1/*Picard order*/, SDIS_DIFFUSION_WOS, &ref); 807 test_picard(scn_3d, 1/*Picard order*/, SDIS_DIFFUSION_WOS, &ref); 808 printf("\n"); 809 810 /* Release memory */ 811 OK(sdis_radiative_env_ref_put(radenv)); 812 OK(sdis_scene_ref_put(scn_2d)); 813 OK(sdis_scene_ref_put(scn_3d)); 814 OK(sdis_interface_ref_put(interfaces[ADIABATIC])); 815 OK(sdis_interface_ref_put(interfaces[SOLID_FLUID_mX])); 816 OK(sdis_interface_ref_put(interfaces[SOLID_FLUID_pX])); 817 OK(sdis_interface_ref_put(interfaces[BOUNDARY_pX])); 818 OK(sdis_medium_ref_put(fluid)); 819 OK(sdis_medium_ref_put(solid)); 820 OK(sdis_medium_ref_put(dummy)); 821 OK(sdis_device_ref_put(dev)); 822 CHK(fclose(stream) == 0); 823 824 CHK(mem_allocated_size() == 0); 825 return 0; 826 }