test_sdis_solve_boundary.c (23926B)
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 21 /* 22 * The scene is composed of a solid cube/square whose temperature is unknown. 23 * The convection coefficient with the surrounding fluid is null exepted for 24 * the +X face whose value is 'H'. The Temperature of the -X face is fixed to 25 * Tb. This test computes the temperature on the +X face and check that it is 26 * equal to: 27 * 28 * T = (H*Tf + LAMBDA/A * Tb) / (H+LAMBDA/A) 29 * 30 * with Tf the temperature of the surrounding fluid, lambda the conductivity of 31 * the cube and A the size of the cube/square, i.e. 1. 32 * 33 * 3D 2D 34 * 35 * ///// (1,1,1) ///// (1,1) 36 * +-------+ +-------+ 37 * /' /| _\ | | _\ 38 * +-------+ | / / Tf Tb | / / Tf 39 * Tb +.....|.+ \__/ | | \__/ 40 * |, |/ +-------+ 41 * +-------+ (0,0) ///// 42 * (0,0,0) ///// 43 */ 44 45 #define N 10000 /* #realisations */ 46 #define N_dump 10 /* #dumped paths */ 47 48 #define Tf 310.0 49 #define Tb 300.0 50 #define H 0.5 51 #define LAMBDA 0.1 52 53 /******************************************************************************* 54 * Media 55 ******************************************************************************/ 56 struct fluid { 57 double temperature; 58 }; 59 60 static double 61 fluid_get_temperature 62 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 63 { 64 CHK(data != NULL && vtx != NULL); 65 return ((const struct fluid*)sdis_data_cget(data))->temperature; 66 } 67 68 static double 69 solid_get_calorific_capacity 70 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 71 { 72 (void)data; 73 CHK(vtx != NULL); 74 return 2.0; 75 } 76 77 static double 78 solid_get_thermal_conductivity 79 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 80 { 81 (void)data; 82 CHK(vtx != NULL); 83 return LAMBDA; 84 } 85 86 static double 87 solid_get_volumic_mass 88 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 89 { 90 (void)data; 91 CHK(vtx != NULL); 92 return 25.0; 93 } 94 95 static double 96 solid_get_delta 97 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 98 { 99 (void)data; 100 CHK(vtx != NULL); 101 return 1.0/20.0; 102 } 103 104 static double 105 solid_get_temperature 106 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 107 { 108 (void)data; 109 CHK(vtx != NULL); 110 if(vtx->time > 0) return SDIS_TEMPERATURE_NONE; 111 return Tf; 112 } 113 114 /******************************************************************************* 115 * Interfaces 116 ******************************************************************************/ 117 struct interf { 118 double temperature; 119 double hc; 120 }; 121 122 static double 123 interface_get_temperature 124 (const struct sdis_interface_fragment* frag, struct sdis_data* data) 125 { 126 const struct interf* interf = sdis_data_cget(data); 127 CHK(frag && data); 128 return interf->temperature; 129 } 130 131 static double 132 interface_get_convection_coef 133 (const struct sdis_interface_fragment* frag, struct sdis_data* data) 134 { 135 const struct interf* interf = sdis_data_cget(data); 136 CHK(frag && data); 137 return interf->hc; 138 } 139 140 /******************************************************************************* 141 * Helper function 142 ******************************************************************************/ 143 static void 144 check_estimator 145 (const struct sdis_estimator* estimator, 146 const size_t nrealisations, /* #realisations */ 147 const double ref) 148 { 149 struct sdis_mc T = SDIS_MC_NULL; 150 struct sdis_mc time = SDIS_MC_NULL; 151 size_t nreals; 152 size_t nfails; 153 CHK(estimator && nrealisations); 154 155 OK(sdis_estimator_get_temperature(estimator, &T)); 156 OK(sdis_estimator_get_realisation_time(estimator, &time)); 157 OK(sdis_estimator_get_realisation_count(estimator, &nreals)); 158 OK(sdis_estimator_get_failure_count(estimator, &nfails)); 159 printf("%g ~ %g +/- %g\n", ref, T.E, T.SE); 160 printf("Time per realisation (in usec) = %g +/- %g\n", time.E, time.SE); 161 printf("#failures = %lu/%lu\n\n", 162 (unsigned long)nfails, (unsigned long)nrealisations); 163 CHK(nfails + nreals == nrealisations); 164 CHK(nfails < N/1000); 165 CHK(eq_eps(T.E, ref, 3*T.SE)); 166 } 167 168 /******************************************************************************* 169 * Test 170 ******************************************************************************/ 171 int 172 main(int argc, char** argv) 173 { 174 FILE* fp = NULL; 175 struct sdis_data* data = NULL; 176 struct sdis_device* dev = NULL; 177 struct sdis_medium* fluid = NULL; 178 struct sdis_medium* solid = NULL; 179 struct sdis_interface* interf_adiabatic = NULL; 180 struct sdis_interface* interf_Tb = NULL; 181 struct sdis_interface* interf_H = NULL; 182 struct sdis_scene* box_scn = NULL; 183 struct sdis_scene* square_scn = NULL; 184 struct sdis_estimator* estimator = NULL; 185 struct sdis_estimator* estimator2 = NULL; 186 struct sdis_green_function* green = NULL; 187 struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT; 188 struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER; 189 struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER; 190 struct sdis_interface_shader interf_shader = SDIS_INTERFACE_SHADER_NULL; 191 struct sdis_interface* box_interfaces[12 /*#triangles*/]; 192 struct sdis_interface* square_interfaces[4/*#segments*/]; 193 struct sdis_solve_probe_boundary_args probe_args = 194 SDIS_SOLVE_PROBE_BOUNDARY_ARGS_DEFAULT; 195 struct sdis_solve_boundary_args bound_args = SDIS_SOLVE_BOUNDARY_ARGS_DEFAULT; 196 struct ssp_rng* rng = NULL; 197 struct interf* interf_props = NULL; 198 struct fluid* fluid_param; 199 double pos[3]; 200 double ref; 201 size_t prims[4]; 202 enum sdis_side sides[4]; 203 int is_master_process = 0; 204 (void)argc, (void)argv; 205 206 create_default_device(&argc, &argv, &is_master_process, &dev); 207 208 /* Temporary file used to dump heat paths */ 209 CHK((fp = tmpfile()) != NULL); 210 211 /* Create the fluid medium */ 212 OK(sdis_data_create 213 (dev, sizeof(struct fluid), ALIGNOF(struct fluid), NULL, &data)); 214 fluid_param = sdis_data_get(data); 215 fluid_param->temperature = Tf; 216 fluid_shader.temperature = fluid_get_temperature; 217 OK(sdis_fluid_create(dev, &fluid_shader, data, &fluid)); 218 OK(sdis_data_ref_put(data)); 219 220 /* Create the solid medium */ 221 solid_shader.calorific_capacity = solid_get_calorific_capacity; 222 solid_shader.thermal_conductivity = solid_get_thermal_conductivity; 223 solid_shader.volumic_mass = solid_get_volumic_mass; 224 solid_shader.delta = solid_get_delta; 225 solid_shader.temperature = solid_get_temperature; 226 OK(sdis_solid_create(dev, &solid_shader, NULL, &solid)); 227 228 /* Setup the interface shader */ 229 interf_shader.convection_coef = interface_get_convection_coef; 230 interf_shader.front.temperature = interface_get_temperature; 231 interf_shader.front.emissivity = NULL; 232 interf_shader.front.specular_fraction = NULL; 233 interf_shader.back = SDIS_INTERFACE_SIDE_SHADER_NULL; 234 235 /* Create the adiabatic interface */ 236 OK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data)); 237 interf_props = sdis_data_get(data); 238 interf_props->hc = 0; 239 interf_props->temperature = SDIS_TEMPERATURE_NONE; 240 OK(sdis_interface_create 241 (dev, solid, fluid, &interf_shader, data, &interf_adiabatic)); 242 OK(sdis_data_ref_put(data)); 243 244 /* Create the Tb interface */ 245 OK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data)); 246 interf_props = sdis_data_get(data); 247 interf_props->hc = 0; 248 interf_props->temperature = Tb; 249 OK(sdis_interface_create 250 (dev, solid, fluid, &interf_shader, data, &interf_Tb)); 251 OK(sdis_data_ref_put(data)); 252 253 /* Create the H interface */ 254 OK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data)); 255 interf_props = sdis_data_get(data); 256 interf_props->hc = H; 257 interf_props->temperature = SDIS_TEMPERATURE_NONE; 258 OK(sdis_interface_create 259 (dev, solid, fluid, &interf_shader, data, &interf_H)); 260 OK(sdis_data_ref_put(data)); 261 262 /* Release the media */ 263 OK(sdis_medium_ref_put(solid)); 264 OK(sdis_medium_ref_put(fluid)); 265 266 /* Map the interfaces to their box triangles */ 267 box_interfaces[0] = box_interfaces[1] = interf_adiabatic; /* Front */ 268 box_interfaces[2] = box_interfaces[3] = interf_Tb; /* Left */ 269 box_interfaces[4] = box_interfaces[5] = interf_adiabatic; /* Back */ 270 box_interfaces[6] = box_interfaces[7] = interf_H; /* Right */ 271 box_interfaces[8] = box_interfaces[9] = interf_adiabatic; /* Top */ 272 box_interfaces[10]= box_interfaces[11]= interf_adiabatic; /* Bottom */ 273 274 /* Map the interfaces to their square segments */ 275 square_interfaces[0] = interf_adiabatic; /* Bottom */ 276 square_interfaces[1] = interf_Tb; /* Lef */ 277 square_interfaces[2] = interf_adiabatic; /* Top */ 278 square_interfaces[3] = interf_H; /* Right */ 279 280 /* Create the box scene */ 281 scn_args.get_indices = box_get_indices; 282 scn_args.get_interface = box_get_interface; 283 scn_args.get_position = box_get_position; 284 scn_args.nprimitives = box_ntriangles; 285 scn_args.nvertices = box_nvertices; 286 scn_args.context = box_interfaces; 287 OK(sdis_scene_create(dev, &scn_args, &box_scn)); 288 289 /* Create the square scene */ 290 scn_args.get_indices = square_get_indices; 291 scn_args.get_interface = square_get_interface; 292 scn_args.get_position = square_get_position; 293 scn_args.nprimitives = square_nsegments; 294 scn_args.nvertices = square_nvertices; 295 scn_args.context = square_interfaces; 296 OK(sdis_scene_2d_create(dev, &scn_args, &square_scn)); 297 298 /* Release the interfaces */ 299 OK(sdis_interface_ref_put(interf_adiabatic)); 300 OK(sdis_interface_ref_put(interf_Tb)); 301 OK(sdis_interface_ref_put(interf_H)); 302 303 ref = (H*Tf + LAMBDA * Tb) / (H + LAMBDA); 304 305 #define SOLVE sdis_solve_probe_boundary 306 #define GREEN sdis_solve_probe_boundary_green_function 307 308 probe_args.nrealisations = N; 309 probe_args.uv[0] = 0.3; 310 probe_args.uv[1] = 0.3; 311 probe_args.iprim = 6; 312 probe_args.time_range[0] = INF; 313 probe_args.time_range[1] = INF; 314 probe_args.side = SDIS_FRONT; 315 316 BA(SOLVE(NULL, &probe_args, &estimator)); 317 BA(SOLVE(box_scn, NULL, &estimator)); 318 BA(SOLVE(box_scn, &probe_args, NULL)); 319 probe_args.nrealisations = 0; 320 BA(SOLVE(box_scn, &probe_args, &estimator)); 321 probe_args.nrealisations = N; 322 probe_args.iprim = 12; 323 BA(SOLVE(box_scn, &probe_args, &estimator)); 324 probe_args.iprim = 6; 325 probe_args.side = SDIS_SIDE_NULL__; 326 BA(SOLVE(box_scn, &probe_args, &estimator)); 327 probe_args.side = SDIS_FRONT; 328 probe_args.time_range[0] = probe_args.time_range[1] = -1; 329 BA(SOLVE(box_scn, &probe_args, &estimator)); 330 probe_args.time_range[0] = 1; 331 BA(SOLVE(box_scn, &probe_args, &estimator)); 332 probe_args.time_range[1] = 0; 333 BA(SOLVE(box_scn, &probe_args, &estimator)); 334 probe_args.time_range[0] = probe_args.time_range[1] = INF; 335 probe_args.picard_order = 0; 336 BA(SOLVE(box_scn, &probe_args, &estimator)); 337 probe_args.picard_order = 1; 338 339 OK(SOLVE(box_scn, &probe_args, &estimator)); 340 OK(sdis_scene_get_boundary_position 341 (box_scn, probe_args.iprim, probe_args.uv, pos)); 342 if(is_master_process) { 343 printf("Boundary temperature of the box at (%g %g %g) = ", SPLIT3(pos)); 344 check_estimator(estimator, N, ref); 345 } 346 347 /* Check RNG type */ 348 probe_args.rng_state = NULL; 349 probe_args.rng_type = SSP_RNG_TYPE_NULL; 350 BA(SOLVE(box_scn, &probe_args, &estimator2)); 351 probe_args.rng_type = 352 SDIS_SOLVE_PROBE_BOUNDARY_ARGS_DEFAULT.rng_type == SSP_RNG_THREEFRY 353 ? SSP_RNG_MT19937_64 : SSP_RNG_THREEFRY; 354 OK(SOLVE(box_scn, &probe_args, &estimator2)); 355 if(is_master_process) { 356 struct sdis_mc T, T2; 357 check_estimator(estimator2, N, ref); 358 OK(sdis_estimator_get_temperature(estimator, &T)); 359 OK(sdis_estimator_get_temperature(estimator2, &T2)); 360 CHK(T2.E != T.E); 361 OK(sdis_estimator_ref_put(estimator2)); 362 } 363 364 /* Check RNG state */ 365 OK(ssp_rng_create(NULL, SSP_RNG_THREEFRY, &rng)); 366 OK(ssp_rng_discard(rng, 31415926535)); /* Move the RNG state */ 367 probe_args.rng_state = rng; 368 probe_args.rng_type = SSP_RNG_TYPE_NULL; 369 OK(SOLVE(box_scn, &probe_args, &estimator2)); 370 OK(ssp_rng_ref_put(rng)); 371 if(is_master_process) { 372 struct sdis_mc T, T2; 373 check_estimator(estimator2, N, ref); 374 OK(sdis_estimator_get_temperature(estimator, &T)); 375 OK(sdis_estimator_get_temperature(estimator2, &T2)); 376 CHK(T2.E != T.E); 377 OK(sdis_estimator_ref_put(estimator2)); 378 } 379 380 probe_args.rng_state = SDIS_SOLVE_PROBE_BOUNDARY_ARGS_DEFAULT.rng_state; 381 probe_args.rng_type = SDIS_SOLVE_PROBE_BOUNDARY_ARGS_DEFAULT.rng_type; 382 383 BA(GREEN(NULL, &probe_args, &green)); 384 BA(GREEN(box_scn, NULL, &green)); 385 BA(GREEN(box_scn, &probe_args, NULL)); 386 probe_args.nrealisations = 0; 387 BA(GREEN(box_scn, &probe_args, &green)); 388 probe_args.nrealisations = N; 389 probe_args.iprim = 12; 390 BA(GREEN(box_scn, &probe_args, &green)); 391 probe_args.iprim = 6; 392 probe_args.side = SDIS_SIDE_NULL__; 393 BA(GREEN(box_scn, &probe_args, &green)); 394 probe_args.side = SDIS_FRONT; 395 OK(GREEN(box_scn, &probe_args, &green)); 396 397 if(!is_master_process) { 398 CHK(estimator == NULL); 399 CHK(green == NULL); 400 } else { 401 check_green_function(green); 402 OK(sdis_green_function_solve(green, &estimator2)); 403 check_estimator(estimator2, N, ref); 404 check_green_serialization(green, box_scn); 405 406 OK(sdis_green_function_ref_put(green)); 407 OK(sdis_estimator_ref_put(estimator)); 408 OK(sdis_estimator_ref_put(estimator2)); 409 } 410 411 /* Dump paths */ 412 probe_args.nrealisations = N_dump; 413 probe_args.register_paths = SDIS_HEAT_PATH_ALL; 414 OK(SOLVE(box_scn, &probe_args, &estimator)); 415 if(is_master_process) { 416 dump_heat_paths(fp, estimator); 417 OK(sdis_estimator_ref_put(estimator)); 418 } 419 420 /* The external fluid cannot have an unknown temperature */ 421 fluid_param->temperature = SDIS_TEMPERATURE_NONE; 422 BA(SOLVE(box_scn, &probe_args, &estimator)); 423 fluid_param->temperature = Tf; 424 425 probe_args.nrealisations = N; 426 probe_args.register_paths = SDIS_HEAT_PATH_NONE; 427 probe_args.uv[0] = 0.5; 428 probe_args.iprim = 4; 429 430 BA(SOLVE(square_scn, &probe_args, &estimator)); 431 probe_args.iprim = 3; 432 OK(SOLVE(square_scn, &probe_args, &estimator)); 433 434 OK(GREEN(square_scn, &probe_args, &green)); 435 if(is_master_process) { 436 check_green_function(green); 437 OK(sdis_green_function_solve(green, &estimator2)); 438 check_estimator(estimator2, N, ref); 439 check_green_serialization(green, square_scn); 440 441 OK(sdis_estimator_ref_put(estimator)); 442 OK(sdis_estimator_ref_put(estimator2)); 443 OK(sdis_green_function_ref_put(green)); 444 } 445 446 /* The external fluid cannot have an unknown temperature */ 447 fluid_param->temperature = SDIS_TEMPERATURE_NONE; 448 BA(SOLVE(square_scn, &probe_args, &estimator)); 449 fluid_param->temperature = Tf; 450 451 /* Right-side temperature at initial time */ 452 probe_args.time_range[0] = 0; 453 probe_args.time_range[1] = 0; 454 455 probe_args.iprim = 6; 456 OK(SOLVE(box_scn, &probe_args, &estimator)); 457 if(is_master_process) { 458 check_estimator(estimator, N, Tf); 459 OK(sdis_estimator_ref_put(estimator)); 460 } 461 462 probe_args.iprim = 3; 463 OK(SOLVE(square_scn, &probe_args, &estimator)); 464 if(is_master_process) { 465 check_estimator(estimator, N, Tf); 466 OK(sdis_estimator_ref_put(estimator)); 467 } 468 469 #undef F 470 #undef SOLVE 471 #undef GREEN 472 473 sides[0] = SDIS_FRONT; 474 sides[1] = SDIS_FRONT; 475 sides[2] = SDIS_FRONT; 476 sides[3] = SDIS_FRONT; 477 478 bound_args.nrealisations = N; 479 bound_args.sides = sides; 480 bound_args.primitives = prims; 481 bound_args.nprimitives = 2; 482 bound_args.time_range[0] = INF; 483 bound_args.time_range[1] = INF; 484 485 #define SOLVE sdis_solve_boundary 486 #define GREEN sdis_solve_boundary_green_function 487 prims[0] = 6; 488 prims[1] = 7; 489 490 BA(SOLVE(NULL, &bound_args, &estimator)); 491 BA(SOLVE(box_scn, NULL, &estimator)); 492 BA(SOLVE(box_scn, &bound_args, NULL)); 493 bound_args.nrealisations = 0; 494 BA(SOLVE(box_scn, &bound_args, &estimator)); 495 bound_args.nrealisations = N; 496 bound_args.primitives = NULL; 497 BA(SOLVE(box_scn, &bound_args, &estimator)); 498 bound_args.primitives = prims; 499 bound_args.sides = NULL; 500 BA(SOLVE(box_scn, &bound_args, &estimator)); 501 bound_args.sides = sides; 502 bound_args.nprimitives = 0; 503 BA(SOLVE(box_scn, &bound_args, &estimator)); 504 bound_args.nprimitives = 2; 505 prims[0] = 12; 506 BA(SOLVE(box_scn, &bound_args, &estimator)); 507 prims[0] = 6; 508 sides[0] = SDIS_SIDE_NULL__; 509 BA(SOLVE(box_scn, &bound_args, &estimator)); 510 sides[0] = SDIS_FRONT; 511 bound_args.time_range[0] = bound_args.time_range[1] = -1; 512 BA(SOLVE(box_scn, &bound_args, &estimator)); 513 bound_args.time_range[0] = 1; 514 BA(SOLVE(box_scn, &bound_args, &estimator)); 515 bound_args.time_range[1] = 0; 516 BA(SOLVE(box_scn, &bound_args, &estimator)); 517 bound_args.time_range[0] = bound_args.time_range[1] = INF; 518 bound_args.picard_order = 0; 519 BA(SOLVE(box_scn, &bound_args, &estimator)); 520 bound_args.picard_order = 1; 521 522 /* Average temperature on the right side of the box */ 523 OK(SOLVE(box_scn, &bound_args, &estimator)); 524 if(is_master_process) { 525 printf("Average temperature of the right side of the box = "); 526 check_estimator(estimator, N, ref); 527 } 528 529 /* Check RNG type */ 530 bound_args.rng_state = NULL; 531 bound_args.rng_type = SSP_RNG_TYPE_NULL; 532 BA(SOLVE(box_scn, &bound_args, &estimator2)); 533 bound_args.rng_type = 534 SDIS_SOLVE_BOUNDARY_ARGS_DEFAULT.rng_type == SSP_RNG_THREEFRY 535 ? SSP_RNG_MT19937_64 : SSP_RNG_THREEFRY; 536 OK(SOLVE(box_scn, &bound_args, &estimator2)); 537 if(is_master_process) { 538 struct sdis_mc T, T2; 539 check_estimator(estimator2, N, ref); 540 OK(sdis_estimator_get_temperature(estimator, &T)); 541 OK(sdis_estimator_get_temperature(estimator2, &T2)); 542 CHK(T2.E != T.E); 543 OK(sdis_estimator_ref_put(estimator2)); 544 } 545 546 /* Check RNG state */ 547 OK(ssp_rng_create(NULL, SSP_RNG_THREEFRY, &rng)); 548 OK(ssp_rng_discard(rng, 31415926535)); /* Move the RNG state */ 549 bound_args.rng_state = rng; 550 bound_args.rng_type = SSP_RNG_TYPE_NULL; 551 OK(SOLVE(box_scn, &bound_args, &estimator2)); 552 OK(ssp_rng_ref_put(rng)); 553 if(is_master_process) { 554 struct sdis_mc T, T2; 555 check_estimator(estimator2, N, ref); 556 OK(sdis_estimator_get_temperature(estimator, &T)); 557 OK(sdis_estimator_get_temperature(estimator2, &T2)); 558 CHK(T2.E != T.E); 559 OK(sdis_estimator_ref_put(estimator2)); 560 } 561 562 /* Restore args */ 563 bound_args.rng_state = SDIS_SOLVE_BOUNDARY_ARGS_DEFAULT.rng_state; 564 bound_args.rng_type = SDIS_SOLVE_BOUNDARY_ARGS_DEFAULT.rng_type; 565 566 BA(GREEN(NULL, &bound_args, &green)); 567 BA(GREEN(box_scn, NULL, &green)); 568 BA(GREEN(box_scn, &bound_args, NULL)); 569 bound_args.nrealisations = 0; 570 BA(GREEN(box_scn, &bound_args, &green)); 571 bound_args.nrealisations = N; 572 bound_args.primitives = NULL; 573 BA(GREEN(box_scn, &bound_args, &green)); 574 bound_args.primitives = prims; 575 bound_args.sides = NULL; 576 BA(GREEN(box_scn, &bound_args, &green)); 577 bound_args.sides = sides; 578 bound_args.nprimitives = 0; 579 BA(GREEN(box_scn, &bound_args, &green)); 580 bound_args.nprimitives = 2; 581 prims[0] = 12; 582 BA(GREEN(box_scn, &bound_args, &green)); 583 prims[0] = 6; 584 sides[0] = SDIS_SIDE_NULL__; 585 BA(GREEN(box_scn, &bound_args, &green)); 586 sides[0] = SDIS_FRONT; 587 588 OK(GREEN(box_scn, &bound_args, &green)); 589 if(!is_master_process) { 590 CHK(estimator == NULL); 591 CHK(green == NULL); 592 } else { 593 check_green_function(green); 594 OK(sdis_green_function_solve(green, &estimator2)); 595 check_estimator(estimator2, N, ref); 596 check_green_serialization(green, box_scn); 597 598 OK(sdis_green_function_ref_put(green)); 599 OK(sdis_estimator_ref_put(estimator)); 600 OK(sdis_estimator_ref_put(estimator2)); 601 } 602 603 /* Dump path */ 604 bound_args.nrealisations = N_dump; 605 bound_args.register_paths = SDIS_HEAT_PATH_ALL; 606 607 /* Check simulation error handling when paths are registered */ 608 fluid_param->temperature = SDIS_TEMPERATURE_NONE; 609 BA(SOLVE(box_scn, &bound_args, &estimator)); 610 611 /* Dump path */ 612 fluid_param->temperature = Tf; 613 OK(SOLVE(box_scn, &bound_args, &estimator)); 614 if(is_master_process) { 615 dump_heat_paths(fp, estimator); 616 OK(sdis_estimator_ref_put(estimator)); 617 } 618 619 /* Switch in 2D */ 620 bound_args.nrealisations = N; 621 bound_args.register_paths = SDIS_HEAT_PATH_NONE; 622 bound_args.nprimitives = 1; 623 prims[0] = 4; 624 BA(SOLVE(square_scn, &bound_args, &estimator)); 625 626 /* Average temperature on the right side of the square */ 627 prims[0] = 3; 628 OK(SOLVE(square_scn, &bound_args, &estimator)); 629 if(is_master_process) { 630 printf("Average temperature of the right side of the square = "); 631 check_estimator(estimator, N, ref); 632 } 633 634 OK(GREEN(square_scn, &bound_args, &green)); 635 if(is_master_process) { 636 check_green_function(green); 637 OK(sdis_green_function_solve(green, &estimator2)); 638 check_estimator(estimator2, N, ref); 639 check_green_serialization(green, square_scn); 640 641 OK(sdis_green_function_ref_put(green)); 642 OK(sdis_estimator_ref_put(estimator)); 643 OK(sdis_estimator_ref_put(estimator2)); 644 } 645 646 /* Dump path */ 647 bound_args.nrealisations = N_dump; 648 bound_args.register_paths = SDIS_HEAT_PATH_ALL; 649 OK(SOLVE(square_scn, &bound_args, &estimator)); 650 if(is_master_process) { 651 dump_heat_paths(fp, estimator); 652 OK(sdis_estimator_ref_put(estimator)); 653 } 654 655 bound_args.register_paths = SDIS_HEAT_PATH_NONE; 656 bound_args.nrealisations = N; 657 658 /* Average temperature on the left+right sides of the box */ 659 prims[0] = 2; 660 prims[1] = 3; 661 prims[2] = 6; 662 prims[3] = 7; 663 664 ref = (ref + Tb) / 2; 665 666 bound_args.nprimitives = 4; 667 OK(SOLVE(box_scn, &bound_args, &estimator)); 668 if(is_master_process) { 669 printf("Average temperature of the left+right sides of the box = "); 670 check_estimator(estimator, N, ref); 671 } 672 673 OK(GREEN(box_scn, &bound_args, &green)); 674 if(is_master_process) { 675 check_green_function(green); 676 OK(sdis_green_function_solve(green, &estimator2)); 677 check_estimator(estimator2, N, ref); 678 check_green_serialization(green, box_scn); 679 680 OK(sdis_green_function_ref_put(green)); 681 OK(sdis_estimator_ref_put(estimator)); 682 OK(sdis_estimator_ref_put(estimator2)); 683 } 684 685 /* Average temperature on the left+right sides of the square */ 686 prims[0] = 1; 687 prims[1] = 3; 688 bound_args.nprimitives = 2; 689 OK(SOLVE(square_scn, &bound_args, &estimator)); 690 if(is_master_process) { 691 printf("Average temperature of the left+right sides of the square = "); 692 check_estimator(estimator, N, ref); 693 } 694 695 OK(GREEN(square_scn, &bound_args, &green)); 696 if(is_master_process) { 697 check_green_function(green); 698 OK(sdis_green_function_solve(green, &estimator2)); 699 check_estimator(estimator2, N, ref); 700 check_green_serialization(green, square_scn); 701 702 OK(sdis_green_function_ref_put(green)); 703 OK(sdis_estimator_ref_put(estimator)); 704 OK(sdis_estimator_ref_put(estimator2)); 705 } 706 707 /* Right-side temperature at initial time */ 708 bound_args.time_range[0] = 0; 709 bound_args.time_range[1] = 0; 710 711 prims[0] = 6; 712 prims[1] = 7; 713 bound_args.nprimitives = 2; 714 OK(SOLVE(box_scn, &bound_args, &estimator)); 715 if(is_master_process) { 716 check_estimator(estimator, N, Tf); 717 OK(sdis_estimator_ref_put(estimator)); 718 } 719 720 prims[0] = 3; 721 bound_args.nprimitives = 1; 722 OK(SOLVE(square_scn, &bound_args, &estimator)); 723 if(is_master_process) { 724 check_estimator(estimator, N, Tf); 725 OK(sdis_estimator_ref_put(estimator)); 726 } 727 728 #undef SOLVE 729 #undef GREEN 730 731 OK(sdis_scene_ref_put(box_scn)); 732 OK(sdis_scene_ref_put(square_scn)); 733 free_default_device(dev); 734 735 CHK(fclose(fp) == 0); 736 737 CHK(mem_allocated_size() == 0); 738 return 0; 739 } 740