test_sdis_solve_boundary_flux.c (22026B)
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 <star/ssp.h> 20 #include <rsys/math.h> 21 22 /* 23 * The scene is composed of a solid cube/square whose temperature is unknown. 24 * The convection coefficient with the surrounding fluid is null excepted for 25 * the X faces whose value is 'H'. The Temperature T of the -X face is fixed 26 * to Tb. The ambient radiative temperature is 0 excepted for the X faces 27 * whose value is 'Trad'. 28 * This test computes temperature and fluxes on the X faces and check that 29 * they are equal to: 30 * 31 * T(+X) = (H*Tf + Hrad*Trad + LAMBDA/A * Tb) / (H+Hrad+LAMBDA/A) 32 * with Hrad = 4 * BOLTZMANN_CONSTANT * Tref^3 * epsilon 33 * T(-X) = Tb 34 * 35 * CF = H * (Tf - T) 36 * RF = Hrad * (Trad - T) 37 * TF = CF + RF 38 * 39 * with Tf the temperature of the surrounding fluid, lambda the conductivity of 40 * the cube and A the size of the cube/square, i.e. 1. 41 * 42 * 3D 43 * 44 * ///////(1,1,1) 45 * +-------+ 46 * /' /| _\ <----- 47 * -----> _\ +-------+ | / / H,Tf <----- Trad 48 * Trad -----> H,Tf / / Tb +.....|.+ \__/ <----- 49 * -----> \__/ |, |/ 50 * +-------+ 51 * (0,0,0)/////// 52 * 53 * 54 * 2D 55 * 56 * ///////(1,1) 57 * +-------+ 58 * -----> _\ | | _\ <----- 59 * Trad -----> H,Tf / / Tb | / / H,Tf <----- Trad 60 * -----> \__/ | | \__/ <----- 61 * +-------+ 62 * (0,0)/////// 63 */ 64 65 #define N 100000 /* #realisations */ 66 67 #define Tf 300.0 68 #define Tb 0.0 69 #define H 0.5 70 #define Trad 300.0 71 #define LAMBDA 0.1 72 #define EPSILON 1.0 73 74 #define Tref 300.0 75 #define Hrad (4 * BOLTZMANN_CONSTANT * Tref * Tref * Tref * EPSILON) 76 77 /******************************************************************************* 78 * Media 79 ******************************************************************************/ 80 struct fluid { 81 double temperature; 82 }; 83 84 static double 85 fluid_get_temperature 86 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 87 { 88 CHK(data != NULL && vtx != NULL); 89 return ((const struct fluid*)sdis_data_cget(data))->temperature; 90 } 91 92 93 static double 94 solid_get_calorific_capacity 95 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 96 { 97 (void) data; 98 CHK(vtx != NULL); 99 return 2.0; 100 } 101 102 static double 103 solid_get_thermal_conductivity 104 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 105 { 106 (void) data; 107 CHK(vtx != NULL); 108 return LAMBDA; 109 } 110 111 static double 112 solid_get_volumic_mass 113 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 114 { 115 (void) data; 116 CHK(vtx != NULL); 117 return 25.0; 118 } 119 120 static double 121 solid_get_delta 122 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 123 { 124 (void) data; 125 CHK(vtx != NULL); 126 return 1.0 / 40.0; 127 } 128 129 static double 130 solid_get_temperature 131 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 132 { 133 (void) data; 134 CHK(vtx != NULL); 135 return SDIS_TEMPERATURE_NONE; 136 } 137 138 /******************************************************************************* 139 * Interfaces 140 ******************************************************************************/ 141 struct interf { 142 double temperature; 143 double emissivity; 144 double hc; 145 double reference_temperature; 146 }; 147 148 static double 149 interface_get_temperature 150 (const struct sdis_interface_fragment* frag, struct sdis_data* data) 151 { 152 const struct interf* interf = sdis_data_cget(data); 153 CHK(frag && data); 154 return interf->temperature; 155 } 156 157 static double 158 interface_get_emissivity 159 (const struct sdis_interface_fragment* frag, 160 const unsigned source_id, 161 struct sdis_data* data) 162 { 163 const struct interf* interf = sdis_data_cget(data); 164 (void)source_id; 165 CHK(frag && data); 166 return interf->emissivity; 167 } 168 169 static double 170 interface_get_convection_coef 171 (const struct sdis_interface_fragment* frag, struct sdis_data* data) 172 { 173 const struct interf* interf = sdis_data_cget(data); 174 CHK(frag && data); 175 return interf->hc; 176 } 177 178 static double 179 interface_get_reference_temperature 180 (const struct sdis_interface_fragment* frag, struct sdis_data* data) 181 { 182 const struct interf* interf = sdis_data_cget(data); 183 CHK(frag && data); 184 return interf->reference_temperature; 185 } 186 187 /******************************************************************************* 188 * Radiative environment 189 ******************************************************************************/ 190 static double 191 radenv_get_temperature 192 (const struct sdis_radiative_ray* ray, 193 struct sdis_data* data) 194 { 195 (void)ray, (void)data; 196 return Trad; 197 } 198 199 static double 200 radenv_get_reference_temperature 201 (const struct sdis_radiative_ray* ray, 202 struct sdis_data* data) 203 { 204 (void)ray, (void)data; 205 return Trad; 206 } 207 208 static struct sdis_radiative_env* 209 create_radenv(struct sdis_device* dev) 210 { 211 struct sdis_radiative_env_shader shader = SDIS_RADIATIVE_ENV_SHADER_NULL; 212 struct sdis_radiative_env* radenv = NULL; 213 214 shader.temperature = radenv_get_temperature; 215 shader.reference_temperature = radenv_get_reference_temperature; 216 OK(sdis_radiative_env_create(dev, &shader, NULL, &radenv)); 217 return radenv; 218 } 219 220 /******************************************************************************* 221 * Helper function 222 ******************************************************************************/ 223 static void 224 check_estimator 225 (const struct sdis_estimator* estimator, 226 const size_t nrealisations, /* #realisations */ 227 const double T, 228 const double CF, 229 const double RF, 230 const double TF) 231 { 232 struct sdis_mc V = SDIS_MC_NULL; 233 enum sdis_estimator_type type; 234 size_t nreals; 235 size_t nfails; 236 CHK(estimator && nrealisations); 237 238 OK(sdis_estimator_get_temperature(estimator, &V)); 239 OK(sdis_estimator_get_realisation_count(estimator, &nreals)); 240 OK(sdis_estimator_get_failure_count(estimator, &nfails)); 241 printf("T = %g ~ %g +/- %g\n", T, V.E, V.SE); 242 CHK(eq_eps(V.E, T, 3 * (V.SE ? V.SE : FLT_EPSILON))); 243 OK(sdis_estimator_get_type(estimator, &type)); 244 if(type == SDIS_ESTIMATOR_FLUX) { 245 OK(sdis_estimator_get_convective_flux(estimator, &V)); 246 printf("Convective flux = %g ~ %g +/- %g\n", CF, V.E, V.SE); 247 CHK(eq_eps(V.E, CF, 3 * (V.SE ? V.SE : FLT_EPSILON))); 248 OK(sdis_estimator_get_radiative_flux(estimator, &V)); 249 printf("Radiative flux = %g ~ %g +/- %g\n", RF, V.E, V.SE); 250 CHK(eq_eps(V.E, RF, 3 * (V.SE ? V.SE : FLT_EPSILON))); 251 OK(sdis_estimator_get_total_flux(estimator, &V)); 252 printf("Total flux = %g ~ %g +/- %g\n", TF, V.E, V.SE); 253 CHK(eq_eps(V.E, TF, 3 * (V.SE ? V.SE : FLT_EPSILON))); 254 } 255 printf("#failures = %lu/%lu\n", 256 (unsigned long) nfails, (unsigned long) nrealisations); 257 CHK(nfails + nreals == nrealisations); 258 CHK(nfails < N / 1000); 259 } 260 261 /******************************************************************************* 262 * Test 263 ******************************************************************************/ 264 int 265 main(int argc, char** argv) 266 { 267 struct sdis_data* data = NULL; 268 struct sdis_device* dev = NULL; 269 struct sdis_medium* fluid1 = NULL; 270 struct sdis_medium* fluid2 = NULL; 271 struct sdis_medium* solid = NULL; 272 struct sdis_interface* interf_adiabatic = NULL; 273 struct sdis_interface* interf_Tb = NULL; 274 struct sdis_interface* interf_H = NULL; 275 struct sdis_radiative_env* radenv = NULL; 276 struct sdis_scene* box_scn = NULL; 277 struct sdis_scene* square_scn = NULL; 278 struct sdis_estimator* estimator = NULL; 279 struct sdis_estimator* estimator2 = NULL; 280 struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT; 281 struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER; 282 struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER; 283 struct sdis_interface_shader interf_shader = SDIS_INTERFACE_SHADER_NULL; 284 struct sdis_interface* box_interfaces[12 /*#triangles*/]; 285 struct sdis_interface* square_interfaces[4/*#segments*/]; 286 struct sdis_solve_probe_boundary_flux_args probe_args = 287 SDIS_SOLVE_PROBE_BOUNDARY_FLUX_ARGS_DEFAULT; 288 struct sdis_solve_boundary_flux_args bound_args = 289 SDIS_SOLVE_BOUNDARY_FLUX_ARGS_DEFAULT; 290 struct interf* interf_props = NULL; 291 struct fluid* fluid_args = NULL; 292 struct ssp_rng* rng = NULL; 293 enum sdis_estimator_type type; 294 double pos[3]; 295 double analyticT, analyticCF, analyticRF, analyticTF; 296 size_t prims[2]; 297 int is_master_process; 298 (void)argc, (void)argv; 299 300 create_default_device(&argc, &argv, &is_master_process, &dev); 301 radenv = create_radenv(dev); 302 303 /* Create the fluid medium. In fact, create two fluid media, even if they are 304 * identical, to check that Robin's boundary conditions can be defined with 305 * several media, without compromising path sampling. Convective paths cannot 306 * be sampled in enclosures with several media, but radiative paths can be, 307 * and it should be possible to calculate the boundary flux without any 308 * problem */ 309 OK(sdis_data_create 310 (dev, sizeof(struct fluid), ALIGNOF(struct fluid), NULL, &data)); 311 fluid_args = sdis_data_get(data); 312 fluid_args->temperature = Tf; 313 fluid_shader.temperature = fluid_get_temperature; 314 OK(sdis_fluid_create(dev, &fluid_shader, data, &fluid1)); 315 OK(sdis_fluid_create(dev, &fluid_shader, data, &fluid2)); 316 OK(sdis_data_ref_put(data)); 317 318 /* Create the solid_medium */ 319 solid_shader.calorific_capacity = solid_get_calorific_capacity; 320 solid_shader.thermal_conductivity = solid_get_thermal_conductivity; 321 solid_shader.volumic_mass = solid_get_volumic_mass; 322 solid_shader.delta = solid_get_delta; 323 solid_shader.temperature = solid_get_temperature; 324 OK(sdis_solid_create(dev, &solid_shader, NULL, &solid)); 325 326 /* Setup the interface shader */ 327 interf_shader.convection_coef = interface_get_convection_coef; 328 interf_shader.front.temperature = interface_get_temperature; 329 interf_shader.front.specular_fraction = NULL; 330 interf_shader.back = SDIS_INTERFACE_SIDE_SHADER_NULL; 331 332 /* Create the adiabatic interface */ 333 OK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data)); 334 interf_props = sdis_data_get(data); 335 interf_props->hc = 0; 336 interf_props->temperature = SDIS_TEMPERATURE_NONE; 337 interf_props->emissivity = 0; 338 OK(sdis_interface_create 339 (dev, solid, fluid1, &interf_shader, data, &interf_adiabatic)); 340 OK(sdis_data_ref_put(data)); 341 342 /* Create the Tb interface */ 343 OK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data)); 344 interf_props = sdis_data_get(data); 345 interf_props->hc = H; 346 interf_props->temperature = Tb; 347 interf_props->emissivity = EPSILON; 348 interf_props->reference_temperature = Tb; 349 interf_shader.back.emissivity = interface_get_emissivity; 350 interf_shader.back.reference_temperature = interface_get_reference_temperature; 351 OK(sdis_interface_create 352 (dev, solid, fluid1, &interf_shader, data, &interf_Tb)); 353 interf_shader.back.emissivity = NULL; 354 OK(sdis_data_ref_put(data)); 355 356 /* Create the H interface */ 357 OK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data)); 358 interf_props = sdis_data_get(data); 359 interf_props->hc = H; 360 interf_props->temperature = SDIS_TEMPERATURE_NONE; 361 interf_props->emissivity = EPSILON; 362 interf_props->reference_temperature = Tref; 363 interf_shader.back.emissivity = interface_get_emissivity; 364 interf_shader.back.reference_temperature = interface_get_reference_temperature; 365 OK(sdis_interface_create 366 (dev, solid, fluid2, &interf_shader, data, &interf_H)); 367 interf_shader.back.emissivity = NULL; 368 OK(sdis_data_ref_put(data)); 369 370 /* Release the media */ 371 OK(sdis_medium_ref_put(solid)); 372 OK(sdis_medium_ref_put(fluid1)); 373 OK(sdis_medium_ref_put(fluid2)); 374 375 /* Map the interfaces to their box triangles */ 376 box_interfaces[0] = box_interfaces[1] = interf_adiabatic; /* Front */ 377 box_interfaces[2] = box_interfaces[3] = interf_Tb; /* Left */ 378 box_interfaces[4] = box_interfaces[5] = interf_adiabatic; /* Back */ 379 box_interfaces[6] = box_interfaces[7] = interf_H; /* Right */ 380 box_interfaces[8] = box_interfaces[9] = interf_adiabatic; /* Top */ 381 box_interfaces[10] = box_interfaces[11] = interf_adiabatic; /* Bottom */ 382 383 /* Map the interfaces to their square segments */ 384 square_interfaces[0] = interf_adiabatic; /* Bottom */ 385 square_interfaces[1] = interf_Tb; /* Lef */ 386 square_interfaces[2] = interf_adiabatic; /* Top */ 387 square_interfaces[3] = interf_H; /* Right */ 388 389 /* Create the box scene */ 390 scn_args.get_indices = box_get_indices; 391 scn_args.get_interface = box_get_interface; 392 scn_args.get_position = box_get_position; 393 scn_args.nprimitives = box_ntriangles; 394 scn_args.nvertices = box_nvertices; 395 scn_args.t_range[0] = MMIN(MMIN(Tf, Trad), Tb); 396 scn_args.t_range[1] = MMAX(MMAX(Tf, Trad), Tb); 397 scn_args.radenv = radenv; 398 scn_args.context = box_interfaces; 399 OK(sdis_scene_create(dev, &scn_args, &box_scn)); 400 401 /* Create the square scene */ 402 scn_args.get_indices = square_get_indices; 403 scn_args.get_interface = square_get_interface; 404 scn_args.get_position = square_get_position; 405 scn_args.nprimitives = square_nsegments; 406 scn_args.nvertices = square_nvertices; 407 scn_args.t_range[0] = MMIN(MMIN(Tf, Trad), Tb); 408 scn_args.t_range[1] = MMAX(MMAX(Tf, Trad), Tb); 409 scn_args.radenv = radenv; 410 scn_args.context = square_interfaces; 411 OK(sdis_scene_2d_create(dev, &scn_args, &square_scn)); 412 413 /* Release the interfaces */ 414 OK(sdis_interface_ref_put(interf_adiabatic)); 415 OK(sdis_interface_ref_put(interf_Tb)); 416 OK(sdis_interface_ref_put(interf_H)); 417 418 analyticT = (H*Tf + Hrad*Trad + LAMBDA * Tb) / (H + Hrad + LAMBDA); 419 analyticCF = H * (Tf - analyticT); 420 analyticRF = Hrad * (Trad - analyticT); 421 analyticTF = analyticCF + analyticRF; 422 423 #define SOLVE sdis_solve_probe_boundary_flux 424 probe_args.nrealisations = N; 425 probe_args.iprim = 6; 426 probe_args.uv[0] = 0.3; 427 probe_args.uv[1] = 0.3; 428 probe_args.time_range[0] = INF; 429 probe_args.time_range[1] = INF; 430 BA(SOLVE(NULL, &probe_args, &estimator)); 431 BA(SOLVE(box_scn, NULL, &estimator)); 432 BA(SOLVE(box_scn, &probe_args, NULL)); 433 probe_args.nrealisations = 0; 434 BA(SOLVE(box_scn, &probe_args, &estimator)); 435 probe_args.nrealisations = N; 436 probe_args.iprim = 12; 437 BA(SOLVE(box_scn, &probe_args, &estimator)); 438 probe_args.iprim = 6; 439 probe_args.uv[0] = probe_args.uv[1] = 1; 440 BA(SOLVE(box_scn, &probe_args, &estimator)); 441 probe_args.uv[0] = probe_args.uv[1] = 0.3; 442 probe_args.time_range[0] = probe_args.time_range[1] = -1; 443 BA(SOLVE(box_scn, &probe_args, &estimator)); 444 probe_args.time_range[0] = 1; 445 BA(SOLVE(box_scn, &probe_args, &estimator)); 446 probe_args.time_range[1] = 0; 447 BA(SOLVE(box_scn, &probe_args, &estimator)); 448 probe_args.time_range[1] = INF; 449 BA(SOLVE(box_scn, &probe_args, &estimator)); 450 probe_args.time_range[0] = INF; 451 OK(SOLVE(box_scn, &probe_args, &estimator)); 452 453 if(!is_master_process) { 454 CHK(estimator == NULL); 455 } else { 456 OK(sdis_estimator_get_type(estimator, &type)); 457 CHK(type == SDIS_ESTIMATOR_FLUX); 458 459 OK(sdis_scene_get_boundary_position 460 (box_scn, probe_args.iprim, probe_args.uv, pos)); 461 printf("Boundary values of the box at (%g %g %g) = ", SPLIT3(pos)); 462 check_estimator(estimator, N, analyticT, analyticCF, analyticRF, analyticTF); 463 } 464 465 /* Check the RNG type */ 466 probe_args.rng_state = NULL; 467 probe_args.rng_type = SSP_RNG_TYPE_NULL; 468 BA(SOLVE(box_scn, &probe_args, &estimator2)); 469 probe_args.rng_type = 470 SDIS_SOLVE_PROBE_BOUNDARY_FLUX_ARGS_DEFAULT.rng_type == SSP_RNG_THREEFRY 471 ? SSP_RNG_MT19937_64 : SSP_RNG_THREEFRY; 472 OK(SOLVE(box_scn, &probe_args, &estimator2)); 473 if(is_master_process) { 474 struct sdis_mc T, T2; 475 check_estimator(estimator2, N, analyticT, analyticCF, analyticRF, analyticTF); 476 OK(sdis_estimator_get_temperature(estimator, &T)); 477 OK(sdis_estimator_get_temperature(estimator2, &T2)); 478 CHK(T2.E != T.E); 479 OK(sdis_estimator_ref_put(estimator2)); 480 } 481 482 /* Check RNG state */ 483 OK(ssp_rng_create(NULL, SSP_RNG_THREEFRY, &rng)); 484 OK(ssp_rng_discard(rng, 31415926535)); /* Move the RNG state */ 485 probe_args.rng_state = rng; 486 probe_args.rng_type = SSP_RNG_TYPE_NULL; 487 OK(SOLVE(box_scn, &probe_args, &estimator2)); 488 OK(ssp_rng_ref_put(rng)); 489 if(is_master_process) { 490 struct sdis_mc T, T2; 491 check_estimator(estimator2, N, analyticT, analyticCF, analyticRF, analyticTF); 492 OK(sdis_estimator_get_temperature(estimator, &T)); 493 OK(sdis_estimator_get_temperature(estimator2, &T2)); 494 CHK(T2.E != T.E); 495 OK(sdis_estimator_ref_put(estimator2)); 496 } 497 498 if(estimator) OK(sdis_estimator_ref_put(estimator)); 499 500 /* Restore arguments */ 501 probe_args.rng_state = SDIS_SOLVE_PROBE_BOUNDARY_FLUX_ARGS_DEFAULT.rng_state; 502 probe_args.rng_type = SDIS_SOLVE_PROBE_BOUNDARY_FLUX_ARGS_DEFAULT.rng_type; 503 504 probe_args.uv[0] = 0.5; 505 probe_args.iprim = 4; 506 BA(SOLVE(square_scn, &probe_args, &estimator)); 507 probe_args.iprim = 3; 508 OK(SOLVE(square_scn, &probe_args, &estimator)); 509 if(is_master_process) { 510 OK(sdis_scene_get_boundary_position 511 (square_scn, probe_args.iprim, probe_args.uv, pos)); 512 printf("Boundary values of the square at (%g %g) = ", SPLIT2(pos)); 513 check_estimator(estimator, N, analyticT, analyticCF, analyticRF, analyticTF); 514 OK(sdis_estimator_ref_put(estimator)); 515 } 516 517 #undef F 518 #undef SOLVE 519 520 #define SOLVE sdis_solve_boundary_flux 521 prims[0] = 6; 522 prims[1] = 7; 523 bound_args.nrealisations = N; 524 bound_args.primitives = prims; 525 bound_args.nprimitives = 2; 526 bound_args.time_range[0] = INF; 527 bound_args.time_range[1] = INF; 528 529 BA(SOLVE(NULL, &bound_args, &estimator)); 530 BA(SOLVE(box_scn, NULL, &estimator)); 531 BA(SOLVE(box_scn, &bound_args, NULL)); 532 bound_args.primitives = NULL; 533 BA(SOLVE(box_scn, &bound_args, &estimator)); 534 bound_args.primitives = prims; 535 bound_args.nprimitives = 0; 536 BA(SOLVE(box_scn, &bound_args, &estimator)); 537 bound_args.nprimitives = 2; 538 bound_args.time_range[0] = bound_args.time_range[1] = -1; 539 BA(SOLVE(box_scn, &bound_args, &estimator)); 540 bound_args.time_range[0] = 1; 541 BA(SOLVE(box_scn, &bound_args, &estimator)); 542 bound_args.time_range[1] = 0; 543 BA(SOLVE(box_scn, &bound_args, &estimator)); 544 bound_args.time_range[1] = INF; 545 BA(SOLVE(box_scn, &bound_args, &estimator)); 546 bound_args.time_range[0] = INF; 547 prims[0] = 12; 548 BA(SOLVE(box_scn, &bound_args, &estimator)); 549 prims[0] = 6; 550 OK(SOLVE(box_scn, &bound_args, &estimator)); 551 552 if(!is_master_process) { 553 CHK(estimator == NULL); 554 } else { 555 /* Average temperature on the right side of the box */ 556 printf("Average values of the right side of the box = "); 557 check_estimator(estimator, N, analyticT, analyticCF, analyticRF, analyticTF); 558 } 559 560 /* Check the RNG type */ 561 bound_args.rng_state = NULL; 562 bound_args.rng_type = SSP_RNG_TYPE_NULL; 563 BA(SOLVE(box_scn, &bound_args, &estimator2)); 564 bound_args.rng_type = 565 SDIS_SOLVE_PROBE_BOUNDARY_FLUX_ARGS_DEFAULT.rng_type == SSP_RNG_THREEFRY 566 ? SSP_RNG_MT19937_64 : SSP_RNG_THREEFRY; 567 OK(SOLVE(box_scn, &bound_args, &estimator2)); 568 if(is_master_process) { 569 struct sdis_mc T, T2; 570 check_estimator(estimator2, N, analyticT, analyticCF, analyticRF, analyticTF); 571 OK(sdis_estimator_get_temperature(estimator, &T)); 572 OK(sdis_estimator_get_temperature(estimator2, &T2)); 573 CHK(T2.E != T.E); 574 OK(sdis_estimator_ref_put(estimator2)); 575 } 576 577 /* Check RNG state */ 578 OK(ssp_rng_create(NULL, SSP_RNG_THREEFRY, &rng)); 579 OK(ssp_rng_discard(rng, 31415926535)); /* Move the RNG state */ 580 bound_args.rng_state = rng; 581 bound_args.rng_type = SSP_RNG_TYPE_NULL; 582 OK(SOLVE(box_scn, &bound_args, &estimator2)); 583 OK(ssp_rng_ref_put(rng)); 584 if(is_master_process) { 585 struct sdis_mc T, T2; 586 check_estimator(estimator2, N, analyticT, analyticCF, analyticRF, analyticTF); 587 OK(sdis_estimator_get_temperature(estimator, &T)); 588 OK(sdis_estimator_get_temperature(estimator2, &T2)); 589 CHK(T2.E != T.E); 590 OK(sdis_estimator_ref_put(estimator2)); 591 } 592 593 if(estimator) OK(sdis_estimator_ref_put(estimator)); 594 595 /* Restore arguments */ 596 bound_args.rng_state = SDIS_SOLVE_BOUNDARY_FLUX_ARGS_DEFAULT.rng_state; 597 bound_args.rng_type = SDIS_SOLVE_BOUNDARY_FLUX_ARGS_DEFAULT.rng_type; 598 599 /* Average temperature on the right side of the square */ 600 prims[0] = 4; 601 bound_args.nprimitives = 1; 602 BA(SOLVE(square_scn, &bound_args, &estimator)); 603 prims[0] = 3; 604 OK(SOLVE(square_scn, &bound_args, &estimator)); 605 if(is_master_process) { 606 printf("Average values of the right side of the square = "); 607 check_estimator(estimator, N, analyticT, analyticCF, analyticRF, analyticTF); 608 OK(sdis_estimator_ref_put(estimator)); 609 } 610 611 /* Flux computation on Dirichlet boundaries is not available yet. 612 * Once available, the expected total flux is the same we expect on the right 613 * side (as the other sides are adiabatic). */ 614 prims[0] = 2; 615 prims[1] = 3; 616 bound_args.nprimitives = 2; 617 BA(SOLVE(box_scn, &bound_args, &estimator)); 618 prims[0] = 1; 619 bound_args.nprimitives = 1; 620 BA(SOLVE(square_scn, &bound_args, &estimator)); 621 #undef SOLVE 622 623 OK(sdis_radiative_env_ref_put(radenv)); 624 OK(sdis_scene_ref_put(box_scn)); 625 OK(sdis_scene_ref_put(square_scn)); 626 free_default_device(dev); 627 628 CHK(mem_allocated_size() == 0); 629 return 0; 630 }