test_sdis_conducto_radiative.c (21627B)
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/math.h> 20 #include <star/ssp.h> 21 22 /* The scene is composed of a solid cube whose temperature is unknown. The cube 23 * faces on +/-X are in contact with a fluid and their convection coefficient 24 * is null while their emissivity is 1. The left and right fluids are enclosed 25 * by surfaces whose emissivity are null excepted for the faces orthogonal to 26 * the X axis that are fully emissive and whose temperature is known. The 27 * medium that surrounds the solid cube and the 2 fluids is a solid with a null 28 * conductivity. 29 * 30 * Y (1, 1, 1) 31 * | +------+----------+------+ (1.5,1,1) 32 * o--- X /' /##########/' /| 33 * / +------+----------+------+ | 34 * Z | ' |##########|*' | | 310K 35 * | ' |##########|*' | | 36 * 300K | ' E=1|##########|*'E=1 | | 37 * | +....|##########|*+....|.+ 38 * |/ |##########|/ |/ 39 * (-1.5,-1,-1) +------+----------+------+ 40 * (-1,-1,-1) 41 */ 42 43 /******************************************************************************* 44 * Geometry 45 ******************************************************************************/ 46 struct geometry { 47 const double* positions; 48 const size_t* indices; 49 struct sdis_interface** interfaces; 50 }; 51 52 static const double vertices[16/*#vertices*/*3/*#coords per vertex*/] = { 53 -1.0,-1.0,-1.0, 54 1.0,-1.0,-1.0, 55 -1.0, 1.0,-1.0, 56 1.0, 1.0,-1.0, 57 -1.0,-1.0, 1.0, 58 1.0,-1.0, 1.0, 59 -1.0, 1.0, 1.0, 60 1.0, 1.0, 1.0, 61 -1.5,-1.0,-1.0, 62 1.5,-1.0,-1.0, 63 -1.5, 1.0,-1.0, 64 1.5, 1.0,-1.0, 65 -1.5,-1.0, 1.0, 66 1.5,-1.0, 1.0, 67 -1.5, 1.0, 1.0, 68 1.5, 1.0, 1.0, 69 }; 70 static const size_t nvertices = sizeof(vertices) / (sizeof(double)*3); 71 72 static const size_t indices[32/*#triangles*/*3/*#indices per triangle*/] = { 73 0, 2, 1, 1, 2, 3, /* Solid back face */ 74 0, 4, 2, 2, 4, 6, /* Solid left face*/ 75 4, 5, 6, 6, 5, 7, /* Solid front face */ 76 3, 7, 1, 1, 7, 5, /* Solid right face */ 77 2, 6, 3, 3, 6, 7, /* Solid top face */ 78 0, 1, 4, 4, 1, 5, /* Solid bottom face */ 79 80 8, 10, 0, 0, 10, 2, /* Left fluid back face */ 81 8, 12, 10, 10, 12, 14, /* Left fluid left face */ 82 12, 4, 14, 14, 4, 6, /* Left fluid front face */ 83 10, 14, 2, 2, 14, 6, /* Left fluid top face */ 84 8, 0, 12, 12, 0, 4, /* Left fluid bottom face */ 85 86 1, 3, 9, 9, 3, 11, /* Right fluid back face */ 87 5, 13, 7, 7, 13, 15, /* Right fluid front face */ 88 11, 15, 9, 9, 15, 13, /* Right fluid right face */ 89 3, 7, 11, 11, 7, 15, /* Right fluid top face */ 90 1, 9, 5, 5, 9, 13 /* Right fluid bottom face */ 91 }; 92 static const size_t ntriangles = sizeof(indices) / (sizeof(size_t)*3); 93 94 static void 95 get_indices(const size_t itri, size_t ids[3], void* ctx) 96 { 97 struct geometry* geom = ctx; 98 CHK(ctx != NULL); 99 ids[0] = geom->indices[itri*3+0]; 100 ids[1] = geom->indices[itri*3+1]; 101 ids[2] = geom->indices[itri*3+2]; 102 } 103 104 static void 105 get_position(const size_t ivert, double pos[3], void* ctx) 106 { 107 struct geometry* geom = ctx; 108 CHK(ctx != NULL); 109 pos[0] = geom->positions[ivert*3+0]; 110 pos[1] = geom->positions[ivert*3+1]; 111 pos[2] = geom->positions[ivert*3+2]; 112 } 113 114 static void 115 get_interface(const size_t itri, struct sdis_interface** bound, void* ctx) 116 { 117 struct geometry* geom = ctx; 118 CHK(ctx != NULL); 119 *bound = geom->interfaces[itri]; 120 } 121 122 /******************************************************************************* 123 * Media 124 ******************************************************************************/ 125 struct solid { 126 double lambda; 127 double initial_temperature; 128 double t0; 129 }; 130 131 static double 132 temperature_unknown(const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 133 { 134 CHK(vtx != NULL); (void)data; 135 return SDIS_TEMPERATURE_NONE; 136 } 137 138 static double 139 solid_get_calorific_capacity 140 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 141 { 142 CHK(vtx != NULL); (void)data; 143 return 1; 144 } 145 146 static double 147 solid_get_thermal_conductivity 148 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 149 { 150 CHK(vtx != NULL); 151 CHK(data != NULL); 152 return ((const struct solid*)sdis_data_cget(data))->lambda; 153 } 154 155 static double 156 solid_get_volumic_mass 157 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 158 { 159 CHK(vtx != NULL); (void)data; 160 return 1; 161 } 162 163 static double 164 solid_get_delta 165 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 166 { 167 CHK(vtx != NULL); (void)data; 168 return 1.0/10.0; 169 } 170 171 static double 172 solid_get_temperature 173 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 174 { 175 double t0; 176 CHK(vtx != NULL); 177 CHK(data != NULL); 178 t0 = ((const struct solid*)sdis_data_cget(data))->t0; 179 if(vtx->time > t0) { 180 return SDIS_TEMPERATURE_NONE; 181 } else { 182 return ((const struct solid*)sdis_data_cget(data))->initial_temperature; 183 } 184 } 185 186 /******************************************************************************* 187 * Interface 188 ******************************************************************************/ 189 struct interfac { 190 double temperature; 191 double convection_coef; 192 double emissivity; 193 double specular_fraction; 194 double Tref; 195 }; 196 197 static double 198 interface_get_temperature 199 (const struct sdis_interface_fragment* frag, struct sdis_data* data) 200 { 201 CHK(data != NULL && frag != NULL); 202 return ((const struct interfac*)sdis_data_cget(data))->temperature; 203 } 204 205 static double 206 interface_get_convection_coef 207 (const struct sdis_interface_fragment* frag, struct sdis_data* data) 208 { 209 CHK(data != NULL && frag != NULL); 210 return ((const struct interfac*)sdis_data_cget(data))->convection_coef; 211 } 212 213 static double 214 interface_get_emissivity 215 (const struct sdis_interface_fragment* frag, 216 const unsigned source_id, 217 struct sdis_data* data) 218 { 219 (void)source_id; 220 CHK(data != NULL && frag != NULL); 221 return ((const struct interfac*)sdis_data_cget(data))->emissivity; 222 } 223 224 static double 225 interface_get_specular_fraction 226 (const struct sdis_interface_fragment* frag, 227 const unsigned source_id, 228 struct sdis_data* data) 229 { 230 (void)source_id; 231 CHK(data != NULL && frag != NULL); 232 return ((const struct interfac*)sdis_data_cget(data))->specular_fraction; 233 } 234 235 static double 236 interface_get_Tref 237 (const struct sdis_interface_fragment* frag, struct sdis_data* data) 238 { 239 CHK(data != NULL && frag != NULL); 240 return ((const struct interfac*)sdis_data_cget(data))->Tref; 241 } 242 243 /******************************************************************************* 244 * Helper functions 245 ******************************************************************************/ 246 static void 247 create_interface 248 (struct sdis_device* dev, 249 struct sdis_medium* front, 250 struct sdis_medium* back, 251 const struct interfac* interf, 252 struct sdis_interface** out_interf) 253 { 254 struct sdis_interface_shader shader = SDIS_INTERFACE_SHADER_NULL; 255 struct sdis_data* data = NULL; 256 257 CHK(interf != NULL); 258 259 shader.front.temperature = interface_get_temperature; 260 shader.back.temperature = interface_get_temperature; 261 if(sdis_medium_get_type(front) != sdis_medium_get_type(back)) { 262 shader.convection_coef = interface_get_convection_coef; 263 } 264 if(sdis_medium_get_type(front) == SDIS_FLUID) { 265 shader.front.emissivity = interface_get_emissivity; 266 shader.front.specular_fraction = interface_get_specular_fraction; 267 shader.front.reference_temperature = interface_get_Tref; 268 } 269 if(sdis_medium_get_type(back) == SDIS_FLUID) { 270 shader.back.emissivity = interface_get_emissivity; 271 shader.back.specular_fraction = interface_get_specular_fraction; 272 shader.back.reference_temperature = interface_get_Tref; 273 } 274 shader.convection_coef_upper_bound = MMAX(0, interf->convection_coef); 275 276 OK(sdis_data_create(dev, sizeof(struct interfac), ALIGNOF(struct interfac), 277 NULL, &data)); 278 *((struct interfac*)sdis_data_get(data)) = *interf; 279 280 OK(sdis_interface_create(dev, front, back, &shader, data, out_interf)); 281 OK(sdis_data_ref_put(data)); 282 } 283 284 /******************************************************************************* 285 * Test that the evaluation of the green function failed with a picard order 286 * greater than 1, i.e. when one want to handle the non-linearties of the 287 * system. 288 ******************************************************************************/ 289 static void 290 test_invalidity_picardN_green 291 (struct sdis_scene* scn, 292 struct sdis_medium* solid) 293 { 294 struct sdis_solve_probe_args probe = SDIS_SOLVE_PROBE_ARGS_DEFAULT; 295 struct sdis_solve_boundary_args bound = SDIS_SOLVE_BOUNDARY_ARGS_DEFAULT; 296 struct sdis_solve_medium_args mdm = SDIS_SOLVE_MEDIUM_ARGS_DEFAULT; 297 struct sdis_solve_probe_boundary_args probe_bound = 298 SDIS_SOLVE_PROBE_BOUNDARY_ARGS_DEFAULT; 299 300 struct sdis_green_function* green = NULL; 301 CHK(scn); 302 303 CHK(probe.picard_order == 1); 304 CHK(probe_bound.picard_order == 1); 305 CHK(bound.picard_order == 1); 306 CHK(mdm.picard_order == 1); 307 308 probe.position[0] = 0; 309 probe.position[1] = 0; 310 probe.position[2] = 0; 311 probe.picard_order = 2; 312 BA(sdis_solve_probe_green_function(scn, &probe, &green)); 313 314 probe_bound.iprim = 2; /* Solid left */ 315 probe_bound.uv[0] = 0.3; 316 probe_bound.uv[1] = 0.3; 317 probe_bound.side = SDIS_FRONT; 318 probe_bound.picard_order = 2; 319 BA(sdis_solve_probe_boundary_green_function(scn, &probe_bound, &green)); 320 321 bound.primitives = &probe_bound.iprim; 322 bound.sides = &probe_bound.side; 323 bound.nprimitives = 1; 324 bound.picard_order = 2; 325 BA(sdis_solve_boundary_green_function(scn, &bound, &green)); 326 327 mdm.medium = solid; 328 mdm.picard_order = 2; 329 BA(sdis_solve_medium_green_function(scn, &mdm, &green)); 330 } 331 332 /******************************************************************************* 333 * Test 334 ******************************************************************************/ 335 int 336 main(int argc, char** argv) 337 { 338 struct mem_allocator allocator; 339 struct interfac interf; 340 struct geometry geom; 341 struct sdis_data* data = NULL; 342 struct sdis_device* dev = NULL; 343 struct sdis_medium* fluid = NULL; 344 struct sdis_medium* solid = NULL; 345 struct sdis_medium* solid2 = NULL; 346 struct sdis_interface* interfaces[5] = {NULL}; 347 struct sdis_interface* prim_interfaces[32/*#triangles*/]; 348 struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT; 349 struct sdis_device_create_args dev_args = SDIS_DEVICE_CREATE_ARGS_DEFAULT; 350 struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER; 351 struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER; 352 struct sdis_scene* scn = NULL; 353 struct ssp_rng* rng = NULL; 354 const int nsimuls = 4; 355 int isimul; 356 const double emissivity = 1;/* Emissivity of the side +/-X of the solid */ 357 const double lambda = 0.1; /* Conductivity of the solid */ 358 const double Tref = 300; /* Reference temperature */ 359 const double T0 = 300; /* Fixed temperature on the left side of the system */ 360 const double T1 = 310; /* Fixed temperature on the right side of the system */ 361 const double thickness = 2.0; /* Thickness of the solid along X */ 362 double t_range[2]; 363 double Ts0, Ts1, hr, tmp; 364 struct interfac* p_intface; 365 (void)argc, (void)argv; 366 367 OK(mem_init_proxy_allocator(&allocator, &mem_default_allocator)); 368 dev_args.allocator = &allocator; 369 OK(sdis_device_create(&dev_args, &dev)); 370 371 /* Create the fluid medium */ 372 fluid_shader.temperature = temperature_unknown; 373 OK(sdis_fluid_create(dev, &fluid_shader, NULL, &fluid)); 374 375 /* Create the solid medium */ 376 OK(sdis_data_create(dev, sizeof(struct solid), ALIGNOF(struct solid), 377 NULL, &data)); 378 ((struct solid*)sdis_data_get(data))->lambda = lambda; 379 ((struct solid*)sdis_data_get(data))->t0 = 0; 380 ((struct solid*)sdis_data_get(data))->initial_temperature = (T0 + T1) / 2; 381 solid_shader.calorific_capacity = solid_get_calorific_capacity; 382 solid_shader.thermal_conductivity = solid_get_thermal_conductivity; 383 solid_shader.volumic_mass = solid_get_volumic_mass; 384 solid_shader.delta = solid_get_delta; 385 solid_shader.temperature = solid_get_temperature; 386 OK(sdis_solid_create(dev, &solid_shader, data, &solid)); 387 OK(sdis_data_ref_put(data)); 388 389 /* Create the surrounding solid medium */ 390 OK(sdis_data_create(dev, sizeof(struct solid), ALIGNOF(struct solid), 391 NULL, &data)); 392 ((struct solid*)sdis_data_get(data))->lambda = 0; 393 solid_shader.calorific_capacity = solid_get_calorific_capacity; 394 solid_shader.thermal_conductivity = solid_get_thermal_conductivity; 395 solid_shader.volumic_mass = solid_get_volumic_mass; 396 solid_shader.delta = solid_get_delta; 397 solid_shader.temperature = temperature_unknown; 398 OK(sdis_solid_create(dev, &solid_shader, data, &solid2)); 399 OK(sdis_data_ref_put(data)); 400 401 /* Create the interface that forces to keep in conduction */ 402 interf.temperature = SDIS_TEMPERATURE_NONE; 403 interf.convection_coef = -1; 404 interf.emissivity = -1; 405 interf.specular_fraction = -1; 406 interf.Tref = Tref; 407 create_interface(dev, solid, solid2, &interf, interfaces+0); 408 409 /* Create the interface that emits radiative heat from the solid */ 410 interf.temperature = SDIS_TEMPERATURE_NONE; 411 interf.convection_coef = 0; 412 interf.emissivity = emissivity; 413 interf.specular_fraction = 1; 414 interf.Tref = Tref; 415 create_interface(dev, solid, fluid, &interf, interfaces+1); 416 417 /* Create the interface that forces the radiative heat to bounce */ 418 interf.temperature = SDIS_TEMPERATURE_NONE; 419 interf.convection_coef = 0; 420 interf.emissivity = 0; 421 interf.specular_fraction = 1; 422 interf.Tref = Tref; 423 create_interface(dev, fluid, solid2, &interf, interfaces+2); 424 425 /* Create the interface with a limit condition of T0 Kelvin */ 426 interf.temperature = T0; 427 interf.convection_coef = 0; 428 interf.emissivity = 1; 429 interf.specular_fraction = 1; 430 interf.Tref = T0; 431 create_interface(dev, fluid, solid2, &interf, interfaces+3); 432 433 /* Create the interface with a limit condition of T1 Kelvin */ 434 interf.temperature = T1; 435 interf.convection_coef = 0; 436 interf.emissivity = 1; 437 interf.specular_fraction = 1; 438 interf.Tref = T1; 439 create_interface(dev, fluid, solid2, &interf, interfaces+4); 440 441 /* Setup the per primitive interface of the solid medium */ 442 prim_interfaces[0] = prim_interfaces[1] = interfaces[0]; 443 prim_interfaces[2] = prim_interfaces[3] = interfaces[1]; 444 prim_interfaces[4] = prim_interfaces[5] = interfaces[0]; 445 prim_interfaces[6] = prim_interfaces[7] = interfaces[1]; 446 prim_interfaces[8] = prim_interfaces[9] = interfaces[0]; 447 prim_interfaces[10] = prim_interfaces[11] = interfaces[0]; 448 449 /* Setup the per primitive interface of the fluid on the left of the medium */ 450 prim_interfaces[12] = prim_interfaces[13] = interfaces[2]; 451 prim_interfaces[14] = prim_interfaces[15] = interfaces[3]; 452 prim_interfaces[16] = prim_interfaces[17] = interfaces[2]; 453 prim_interfaces[18] = prim_interfaces[19] = interfaces[2]; 454 prim_interfaces[20] = prim_interfaces[21] = interfaces[2]; 455 456 /* Setup the per primitive interface of the fluid on the right of the medium */ 457 prim_interfaces[22] = prim_interfaces[23] = interfaces[2]; 458 prim_interfaces[24] = prim_interfaces[25] = interfaces[2]; 459 prim_interfaces[26] = prim_interfaces[27] = interfaces[4]; 460 prim_interfaces[28] = prim_interfaces[29] = interfaces[2]; 461 prim_interfaces[30] = prim_interfaces[31] = interfaces[2]; 462 463 /* Create the scene */ 464 geom.positions = vertices; 465 geom.indices = indices; 466 geom.interfaces = prim_interfaces; 467 scn_args.get_indices = get_indices; 468 scn_args.get_interface = get_interface; 469 scn_args.get_position = get_position; 470 scn_args.nprimitives = ntriangles; 471 scn_args.nvertices = nvertices; 472 scn_args.t_range[0] = MMIN(T0, T1); 473 scn_args.t_range[1] = MMAX(T0, T1); 474 scn_args.context = &geom; 475 OK(sdis_scene_create(dev, &scn_args, &scn)); 476 477 hr = 4.0 * BOLTZMANN_CONSTANT * Tref*Tref*Tref * emissivity; 478 479 /* Run the simulations */ 480 p_intface 481 = (struct interfac*)sdis_data_get(sdis_interface_get_data(interfaces[4])); 482 OK(ssp_rng_create(&allocator, SSP_RNG_KISS, &rng)); 483 FOR_EACH(isimul, 0, nsimuls) { 484 struct sdis_mc T = SDIS_MC_NULL; 485 struct sdis_mc time = SDIS_MC_NULL; 486 struct sdis_estimator* estimator; 487 struct sdis_estimator* estimator2; 488 struct sdis_green_function* green; 489 struct sdis_solve_probe_args solve_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT; 490 double ref = SDIS_TEMPERATURE_NONE; 491 size_t nreals = 0; 492 size_t nfails = 0; 493 const size_t N = 10000; 494 double T1b; 495 int steady = (isimul % 2) == 0; 496 497 /* Reset temperature */ 498 p_intface->temperature = T1; 499 500 solve_args.nrealisations = N; 501 if(steady) { 502 solve_args.time_range[0] = solve_args.time_range[1] = INF; 503 } else { 504 solve_args.time_range[0] = 100 * (double)isimul; 505 solve_args.time_range[1] = 4 * solve_args.time_range[0]; 506 } 507 solve_args.position[0] = ssp_rng_uniform_double(rng, -0.9, 0.9); 508 solve_args.position[1] = ssp_rng_uniform_double(rng, -0.9, 0.9); 509 solve_args.position[2] = ssp_rng_uniform_double(rng, -0.9, 0.9); 510 511 OK(sdis_solve_probe(scn, &solve_args, &estimator)); 512 OK(sdis_estimator_get_realisation_count(estimator, &nreals)); 513 OK(sdis_estimator_get_failure_count(estimator, &nfails)); 514 OK(sdis_estimator_get_temperature(estimator, &T)); 515 OK(sdis_estimator_get_realisation_time(estimator, &time)); 516 517 tmp = lambda / (2 * lambda + thickness * hr) * (T1 - T0); 518 Ts0 = T0 + tmp; 519 Ts1 = T1 - tmp; 520 if(steady) { 521 double u = (solve_args.position[0] + 1) / thickness; 522 ref = u * Ts1 + (1 - u) * Ts0; 523 printf("Steady temperature at (%g, %g, %g) with T1=%g = %g ~ %g +/- %g\n", 524 SPLIT3(solve_args.position), p_intface->temperature, ref, T.E, T.SE); 525 } else { 526 printf("Mean temperature at (%g, %g, %g) with t in [%g %g]" 527 " and T1=%g ~ %g +/- %g\n", 528 SPLIT3(solve_args.position), SPLIT2(solve_args.time_range), 529 p_intface->temperature, T.E, T.SE); 530 } 531 printf("Time per realisation (in usec) = %g +/- %g\n", time.E, time.SE); 532 printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); 533 534 CHK(nfails + nreals == N); 535 CHK(nfails <= N/1000); 536 if(steady) CHK(eq_eps(T.E, ref, 3*T.SE) == 1); 537 538 /* Check green function */ 539 OK(sdis_solve_probe_green_function(scn, &solve_args, &green)); 540 OK(sdis_green_function_solve(green, &estimator2)); 541 check_green_function(green); 542 check_estimator_eq(estimator, estimator2); 543 check_green_serialization(green, scn); 544 545 OK(sdis_estimator_ref_put(estimator)); 546 OK(sdis_estimator_ref_put(estimator2)); 547 printf("\n"); 548 549 /* Check same green used at a different temperature */ 550 p_intface->temperature = T1b = T1 + ((double)isimul + 1) * 10; 551 t_range[0] = MMIN(T0, T1b); 552 t_range[1] = MMAX(T0, T1b); 553 OK(sdis_scene_set_temperature_range(scn, t_range)); 554 555 OK(sdis_solve_probe(scn, &solve_args, &estimator)); 556 OK(sdis_estimator_get_realisation_count(estimator, &nreals)); 557 OK(sdis_estimator_get_failure_count(estimator, &nfails)); 558 OK(sdis_estimator_get_temperature(estimator, &T)); 559 OK(sdis_estimator_get_realisation_time(estimator, &time)); 560 561 tmp = lambda / (2 * lambda + thickness * hr) * (T1b - T0); 562 Ts0 = T0 + tmp; 563 Ts1 = T1b - tmp; 564 565 if(steady) { 566 double u = (solve_args.position[0] + 1) / thickness; 567 ref = u * Ts1 + (1 - u) * Ts0; 568 printf("Steady temperature at (%g, %g, %g) with T1=%g = %g ~ %g +/- %g\n", 569 SPLIT3(solve_args.position), p_intface->temperature, ref, T.E, T.SE); 570 } else { 571 printf("Mean temperature at (%g, %g, %g) with t in [%g %g]" 572 " and T1=%g ~ %g +/- %g\n", 573 SPLIT3(solve_args.position), SPLIT2(solve_args.time_range), 574 p_intface->temperature, T.E, T.SE); 575 } 576 printf("Time per realisation (in usec) = %g +/- %g\n", time.E, time.SE); 577 printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); 578 579 CHK(nfails + nreals == N); 580 CHK(nfails <= N/1000); 581 if(steady) CHK(eq_eps(T.E, ref, 3*T.SE) == 1); 582 583 OK(sdis_green_function_solve(green, &estimator2)); 584 check_green_function(green); 585 check_estimator_eq(estimator, estimator2); 586 587 OK(sdis_estimator_ref_put(estimator)); 588 OK(sdis_estimator_ref_put(estimator2)); 589 OK(sdis_green_function_ref_put(green)); 590 591 solve_args.nrealisations = 10; 592 solve_args.register_paths = SDIS_HEAT_PATH_ALL; 593 594 OK(sdis_solve_probe(scn, &solve_args, &estimator)); 595 OK(sdis_estimator_ref_put(estimator)); 596 597 printf("\n\n"); 598 } 599 600 test_invalidity_picardN_green(scn, solid); 601 602 /* Release memory */ 603 OK(sdis_scene_ref_put(scn)); 604 OK(sdis_interface_ref_put(interfaces[0])); 605 OK(sdis_interface_ref_put(interfaces[1])); 606 OK(sdis_interface_ref_put(interfaces[2])); 607 OK(sdis_interface_ref_put(interfaces[3])); 608 OK(sdis_interface_ref_put(interfaces[4])); 609 OK(sdis_medium_ref_put(fluid)); 610 OK(sdis_medium_ref_put(solid)); 611 OK(sdis_medium_ref_put(solid2)); 612 OK(sdis_device_ref_put(dev)); 613 OK(ssp_rng_ref_put(rng)); 614 615 check_memory_allocator(&allocator); 616 mem_shutdown_proxy_allocator(&allocator); 617 CHK(mem_allocated_size() == 0); 618 return 0; 619 }