test_sdis_solve_probe_list.c (23085B)
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/float3.h> 20 21 #include <star/s3d.h> 22 #include <star/s3dut.h> 23 #include <star/ssp.h> 24 25 #ifdef SDIS_ENABLE_MPI 26 #include <mpi.h> 27 #endif 28 29 /* 30 * The system is a trilinear profile of steady state temperature, i.e. at each 31 * point of the system we can analytically calculate the temperature. We immerse 32 * a super shape in this temperature field which represents a solid in which we 33 * want to evaluate by Monte Carlo the temperature at several positions. On the 34 * Monte Carlo side, the temperature of the super shape is unknown. Only the 35 * boundary temperature is set to the temperature of the aforementioned 36 * trilinear profile. Thus, we should find by Monte Carlo the temperature 37 * defined by the trilinear profile. 38 * 39 * T(z) /\ <-- T(x,y,z) 40 * | T(y) ___/ \___ 41 * |/ \ . T=? / 42 * o--- T(x) /_ __ _\ 43 * \/ \/ 44 * 45 */ 46 47 /******************************************************************************* 48 * Helper functions 49 ******************************************************************************/ 50 static double 51 trilinear_profile(const double pos[3]) 52 { 53 /* Range in X, Y and Z in which the trilinear profile is defined */ 54 const double lower = -3; 55 const double upper = +3; 56 57 /* Upper temperature limit in X, Y and Z [K]. Lower temperature limit is 58 * implicitly 0 */ 59 const double a = 333; /* Upper temperature limit in X [K] */ 60 const double b = 432; /* Upper temperature limit in Y [K] */ 61 const double c = 579; /* Upper temperature limit in Z [K] */ 62 63 double x, y, z; 64 65 /* Check pre-conditions */ 66 CHK(pos); 67 CHK(lower <= pos[0] && pos[0] <= upper); 68 CHK(lower <= pos[1] && pos[1] <= upper); 69 CHK(lower <= pos[2] && pos[2] <= upper); 70 71 x = (pos[0] - lower) / (upper - lower); 72 y = (pos[1] - lower) / (upper - lower); 73 z = (pos[2] - lower) / (upper - lower); 74 return a*x + b*y + c*z; 75 } 76 77 static INLINE float 78 rand_canonic(void) 79 { 80 return (float)rand() / (float)((unsigned)RAND_MAX+1); 81 } 82 83 static INLINE void 84 check_intersection 85 (const double val0, 86 const double eps0, 87 const double val1, 88 const double eps1) 89 { 90 double interval0[2], interval1[2]; 91 double intersection[2]; 92 interval0[0] = val0 - eps0; 93 interval0[1] = val0 + eps0; 94 interval1[0] = val1 - eps1; 95 interval1[1] = val1 + eps1; 96 intersection[0] = MMAX(interval0[0], interval1[0]); 97 intersection[1] = MMIN(interval0[1], interval1[1]); 98 CHK(intersection[0] <= intersection[1]); 99 } 100 101 /******************************************************************************* 102 * Super shape 103 ******************************************************************************/ 104 static struct s3dut_mesh* 105 create_super_shape(void) 106 { 107 struct s3dut_mesh* mesh = NULL; 108 struct s3dut_super_formula f0 = S3DUT_SUPER_FORMULA_NULL; 109 struct s3dut_super_formula f1 = S3DUT_SUPER_FORMULA_NULL; 110 const double radius = 1; 111 const unsigned nslices = 256; 112 113 f0.A = 1.5; f0.B = 1; f0.M = 11.0; f0.N0 = 1; f0.N1 = 1; f0.N2 = 2.0; 114 f1.A = 1.0; f1.B = 2; f1.M = 3.6; f1.N0 = 1; f1.N1 = 2; f1.N2 = 0.7; 115 OK(s3dut_create_super_shape(NULL, &f0, &f1, radius, nslices, nslices/2, &mesh)); 116 117 return mesh; 118 } 119 120 /******************************************************************************* 121 * View, i.e. acceleration structure used to query geometry. In this test it is 122 * used to calculate the delta parameter and to sample the probe positions in 123 * the supershape. 124 ******************************************************************************/ 125 static void 126 view_get_indices(const unsigned itri, unsigned ids[3], void* ctx) 127 { 128 struct s3dut_mesh_data* mesh_data = ctx; 129 CHK(ids && mesh_data && itri < mesh_data->nprimitives); 130 /* Flip the indices to ensure that the normal points into the super shape */ 131 ids[0] = (unsigned)mesh_data->indices[itri*3+0]; 132 ids[1] = (unsigned)mesh_data->indices[itri*3+2]; 133 ids[2] = (unsigned)mesh_data->indices[itri*3+1]; 134 } 135 136 static void 137 view_get_position(const unsigned ivert, float pos[3], void* ctx) 138 { 139 struct s3dut_mesh_data* mesh_data = ctx; 140 CHK(pos && mesh_data && ivert < mesh_data->nvertices); 141 pos[0] = (float)mesh_data->positions[ivert*3+0]; 142 pos[1] = (float)mesh_data->positions[ivert*3+1]; 143 pos[2] = (float)mesh_data->positions[ivert*3+2]; 144 } 145 146 static struct s3d_scene_view* 147 create_view(struct s3dut_mesh* mesh) 148 { 149 struct s3d_vertex_data vdata = S3D_VERTEX_DATA_NULL; 150 struct s3d_device* s3d = NULL; 151 struct s3d_scene* scn = NULL; 152 struct s3d_shape* shape = NULL; 153 struct s3d_scene_view* view = NULL; 154 155 struct s3dut_mesh_data mesh_data; 156 157 OK(s3dut_mesh_get_data(mesh, &mesh_data)); 158 159 vdata.usage = S3D_POSITION; 160 vdata.type = S3D_FLOAT3; 161 vdata.get = view_get_position; 162 OK(s3d_device_create(NULL, NULL, 0, &s3d)); 163 OK(s3d_shape_create_mesh(s3d, &shape)); 164 OK(s3d_mesh_setup_indexed_vertices(shape, (unsigned)mesh_data.nprimitives, 165 view_get_indices, (unsigned)mesh_data.nvertices, &vdata, 1, &mesh_data)); 166 OK(s3d_scene_create(s3d, &scn)); 167 OK(s3d_scene_attach_shape(scn, shape)); 168 OK(s3d_scene_view_create(scn, S3D_TRACE, &view)); 169 170 OK(s3d_device_ref_put(s3d)); 171 OK(s3d_shape_ref_put(shape)); 172 OK(s3d_scene_ref_put(scn)); 173 174 return view; 175 } 176 177 static double 178 view_compute_delta(struct s3d_scene_view* view) 179 { 180 float S = 0; /* Surface */ 181 float V = 0; /* Volume */ 182 183 OK(s3d_scene_view_compute_area(view, &S)); 184 OK(s3d_scene_view_compute_volume(view, &V)); 185 CHK(S > 0 && V > 0); 186 187 return (4.0*V/S)/30.0; 188 } 189 190 static void 191 view_sample_position 192 (struct s3d_scene_view* view, 193 const double delta, 194 double pos[3]) 195 { 196 /* Ray variables */ 197 const float dir[3] = {0, 1, 0}; 198 const float range[2] = {0, FLT_MAX}; 199 float org[3]; 200 201 /* View variables */ 202 float low[3]; 203 float upp[3]; 204 205 OK(s3d_scene_view_get_aabb(view, low, upp)); 206 207 /* Sample a position in the supershape by uniformly sampling a position within 208 * its bounding box and rejecting it if it is not in the supershape or if it 209 * is too close to its boundaries. */ 210 for(;;) { 211 struct s3d_hit hit = S3D_HIT_NULL; 212 float N[3] = {0, 0, 0}; /* Normal */ 213 214 pos[0] = low[0] + rand_canonic()*(upp[0] - low[0]); 215 pos[1] = low[1] + rand_canonic()*(upp[1] - low[1]); 216 pos[2] = low[2] + rand_canonic()*(upp[2] - low[2]); 217 218 org[0] = (float)pos[0]; 219 org[1] = (float)pos[1]; 220 org[2] = (float)pos[2]; 221 OK(s3d_scene_view_trace_ray(view, org, dir, range, NULL, &hit)); 222 223 if(S3D_HIT_NONE(&hit)) continue; 224 225 f3_normalize(N, hit.normal); 226 if(f3_dot(N, dir) > 0) continue; 227 228 OK(s3d_scene_view_closest_point(view, org, (float)INF, NULL, &hit)); 229 CHK(!S3D_HIT_NONE(&hit)); 230 231 /* Sample position is in the super shape and is not too close to its 232 * boundaries */ 233 if(hit.distance > delta) break; 234 } 235 } 236 237 /******************************************************************************* 238 * Scene, i.e. the system to simulate 239 ******************************************************************************/ 240 struct scene_context { 241 struct s3dut_mesh_data mesh_data; 242 struct sdis_interface* interf; 243 }; 244 static const struct scene_context SCENE_CONTEXT_NULL = {{0}, NULL}; 245 246 static void 247 scene_get_indices(const size_t itri, size_t ids[3], void* ctx) 248 { 249 struct scene_context* context = ctx; 250 CHK(ids && context && itri < context->mesh_data.nprimitives); 251 /* Flip the indices to ensure that the normal points into the super shape */ 252 ids[0] = (unsigned)context->mesh_data.indices[itri*3+0]; 253 ids[1] = (unsigned)context->mesh_data.indices[itri*3+2]; 254 ids[2] = (unsigned)context->mesh_data.indices[itri*3+1]; 255 } 256 257 static void 258 scene_get_interface(const size_t itri, struct sdis_interface** interf, void* ctx) 259 { 260 struct scene_context* context = ctx; 261 CHK(interf && context && itri < context->mesh_data.nprimitives); 262 *interf = context->interf; 263 } 264 265 static void 266 scene_get_position(const size_t ivert, double pos[3], void* ctx) 267 { 268 struct scene_context* context = ctx; 269 CHK(pos && context && ivert < context->mesh_data.nvertices); 270 pos[0] = context->mesh_data.positions[ivert*3+0]; 271 pos[1] = context->mesh_data.positions[ivert*3+1]; 272 pos[2] = context->mesh_data.positions[ivert*3+2]; 273 } 274 275 static struct sdis_scene* 276 create_scene 277 (struct sdis_device* sdis, 278 const struct s3dut_mesh* mesh, 279 struct sdis_interface* interf) 280 { 281 struct sdis_scene* scn = NULL; 282 struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT; 283 struct scene_context context = SCENE_CONTEXT_NULL; 284 285 OK(s3dut_mesh_get_data(mesh, &context.mesh_data)); 286 context.interf = interf; 287 288 scn_args.get_indices = scene_get_indices; 289 scn_args.get_interface = scene_get_interface; 290 scn_args.get_position = scene_get_position; 291 scn_args.nprimitives = context.mesh_data.nprimitives; 292 scn_args.nvertices = context.mesh_data.nvertices; 293 scn_args.context = &context; 294 OK(sdis_scene_create(sdis, &scn_args, &scn)); 295 return scn; 296 } 297 298 /******************************************************************************* 299 * Solid, i.e. medium of the super shape 300 ******************************************************************************/ 301 #define SOLID_PROP(Prop, Val) \ 302 static double \ 303 solid_get_##Prop \ 304 (const struct sdis_rwalk_vertex* vtx, \ 305 struct sdis_data* data) \ 306 { \ 307 (void)vtx, (void)data; /* Avoid the "unused variable" warning */ \ 308 return Val; \ 309 } 310 SOLID_PROP(calorific_capacity, 500.0) /* [J/K/Kg] */ 311 SOLID_PROP(thermal_conductivity, 25.0) /* [W/m/K] */ 312 SOLID_PROP(volumic_mass, 7500.0) /* [kg/m^3] */ 313 SOLID_PROP(temperature, SDIS_TEMPERATURE_NONE) /* [K] */ 314 #undef SOLID_PROP 315 316 static double 317 solid_get_delta 318 (const struct sdis_rwalk_vertex* vtx, 319 struct sdis_data* data) 320 { 321 const double* delta = sdis_data_get(data); 322 (void)vtx; /* Avoid the "unused variable" warning */ 323 return *delta; 324 } 325 326 static struct sdis_medium* 327 create_solid(struct sdis_device* sdis, struct s3d_scene_view* view) 328 { 329 struct sdis_solid_shader shader = SDIS_SOLID_SHADER_NULL; 330 struct sdis_medium* solid = NULL; 331 struct sdis_data* data = NULL; 332 double* pdelta = NULL; 333 334 OK(sdis_data_create(sdis, sizeof(double), ALIGNOF(double), NULL, &data)); 335 pdelta = sdis_data_get(data); 336 *pdelta = view_compute_delta(view); 337 338 shader.calorific_capacity = solid_get_calorific_capacity; 339 shader.thermal_conductivity = solid_get_thermal_conductivity; 340 shader.volumic_mass = solid_get_volumic_mass; 341 shader.delta = solid_get_delta; 342 shader.temperature = solid_get_temperature; 343 OK(sdis_solid_create(sdis, &shader, data, &solid)); 344 345 OK(sdis_data_ref_put(data)); 346 return solid; 347 } 348 349 /******************************************************************************* 350 * Dummy environment, i.e. environment surrounding the super shape. It is 351 * defined only for Stardis compliance: in Stardis, an interface must divide 2 352 * media. 353 ******************************************************************************/ 354 static struct sdis_medium* 355 create_dummy(struct sdis_device* sdis) 356 { 357 struct sdis_fluid_shader shader = SDIS_FLUID_SHADER_NULL; 358 struct sdis_medium* dummy = NULL; 359 360 shader.calorific_capacity = dummy_medium_getter; 361 shader.volumic_mass = dummy_medium_getter; 362 shader.temperature = dummy_medium_getter; 363 OK(sdis_fluid_create(sdis, &shader, NULL, &dummy)); 364 return dummy; 365 } 366 367 /******************************************************************************* 368 * Interface: its temperature is fixed to the trilinear profile 369 ******************************************************************************/ 370 static double 371 interface_get_temperature 372 (const struct sdis_interface_fragment* frag, 373 struct sdis_data* data) 374 { 375 (void)data; /* Avoid the "unused variable" warning */ 376 return trilinear_profile(frag->P); 377 } 378 379 static struct sdis_interface* 380 create_interface 381 (struct sdis_device* sdis, 382 struct sdis_medium* front, 383 struct sdis_medium* back) 384 { 385 struct sdis_interface* interf = NULL; 386 struct sdis_interface_shader shader = SDIS_INTERFACE_SHADER_NULL; 387 388 shader.front.temperature = interface_get_temperature; 389 shader.back.temperature = interface_get_temperature; 390 OK(sdis_interface_create(sdis, front, back, &shader, NULL, &interf)); 391 return interf; 392 } 393 394 /******************************************************************************* 395 * Validations 396 ******************************************************************************/ 397 static void 398 check_probe_list_api 399 (struct sdis_scene* scn, 400 const struct sdis_solve_probe_list_args* in_args) 401 { 402 struct sdis_solve_probe_list_args args = *in_args; 403 struct sdis_estimator_buffer* estim_buf = NULL; 404 405 /* Check API */ 406 BA(sdis_solve_probe_list(NULL, &args, &estim_buf)); 407 BA(sdis_solve_probe_list(scn, NULL, &estim_buf)); 408 BA(sdis_solve_probe_list(scn, &args, NULL)); 409 args.nprobes = 0; 410 BA(sdis_solve_probe_list(scn, &args, &estim_buf)); 411 args.nprobes = in_args->nprobes; 412 args.probes = NULL; 413 BA(sdis_solve_probe_list(scn, &args, &estim_buf)); 414 } 415 416 /* Check the estimators against the analytical solution */ 417 static void 418 check_estimator_buffer 419 (const struct sdis_estimator_buffer* estim_buf, 420 const struct sdis_solve_probe_list_args* args) 421 { 422 struct sdis_mc T = SDIS_MC_NULL; 423 size_t iprobe = 0; 424 425 /* Variables used to check the global estimation results */ 426 size_t total_nrealisations = 0; 427 size_t total_nrealisations_ref = 0; 428 size_t total_nfailures = 0; 429 size_t total_nfailures_ref = 0; 430 double sum = 0; /* Global sum of weights */ 431 double sum2 = 0; /* Global sum of squared weights */ 432 double N = 0; /* Number of (successful) realisations */ 433 double E = 0; /* Expected value */ 434 double V = 0; /* Variance */ 435 double SE = 0; /* Standard Error */ 436 437 /* Check the results */ 438 FOR_EACH(iprobe, 0, args->nprobes) { 439 const struct sdis_estimator* estimator = NULL; 440 size_t probe_nrealisations = 0; 441 size_t probe_nfailures = 0; 442 double probe_sum = 0; 443 double probe_sum2 = 0; 444 double ref = 0; 445 446 /* Fetch result */ 447 OK(sdis_estimator_buffer_at(estim_buf, iprobe, 0, &estimator)); 448 OK(sdis_estimator_get_temperature(estimator, &T)); 449 OK(sdis_estimator_get_realisation_count(estimator, &probe_nrealisations)); 450 OK(sdis_estimator_get_failure_count(estimator, &probe_nfailures)); 451 452 /* Check probe estimation */ 453 ref = trilinear_profile(args->probes[iprobe].position); 454 printf("T(%g, %g, %g) = %g ~ %g +/- %g\n", 455 SPLIT3(args->probes[iprobe].position), ref, T.E, T.SE); 456 CHK(eq_eps(ref, T.E, 3*T.SE)); 457 458 /* Check miscellaneous results */ 459 CHK(probe_nrealisations == args->probes[iprobe].nrealisations); 460 461 /* Compute the sum of weights and sum of squared weights of the probe */ 462 probe_sum = T.E * (double)probe_nrealisations; 463 probe_sum2 = (T.V + T.E*T.E) * (double)probe_nrealisations; 464 465 /* Update the global estimation */ 466 total_nrealisations_ref += probe_nrealisations; 467 total_nfailures_ref += probe_nfailures; 468 sum += probe_sum; 469 sum2 += probe_sum2; 470 } 471 472 /* Check the overall estimate of the estimate buffer. This estimate is not 473 * really a result expected by the caller since the probes are all 474 * independent. But to be consistent with the estimate buffer API, we need to 475 * provide it. And so we check its validity */ 476 OK(sdis_estimator_buffer_get_temperature(estim_buf, &T)); 477 OK(sdis_estimator_buffer_get_realisation_count(estim_buf, &total_nrealisations)); 478 OK(sdis_estimator_buffer_get_failure_count(estim_buf, &total_nfailures)); 479 480 CHK(total_nrealisations == total_nrealisations_ref); 481 CHK(total_nfailures == total_nfailures_ref); 482 483 N = (double)total_nrealisations_ref; 484 E = sum / N; 485 V = sum2 / N - E*E; 486 SE = sqrt(V/N); 487 check_intersection(E, SE, T.E, T.SE); 488 } 489 490 /* Check that the buffers store statistically independent estimates */ 491 static void 492 check_estimator_buffer_combination 493 (const struct sdis_estimator_buffer* estim_buf1, 494 const struct sdis_estimator_buffer* estim_buf2, 495 const struct sdis_solve_probe_list_args* args) 496 { 497 size_t iprobe; 498 499 /* Check that the 2 series of estimates are compatible but not identical */ 500 FOR_EACH(iprobe, 0, args->nprobes) { 501 const struct sdis_estimator* estimator1 = NULL; 502 const struct sdis_estimator* estimator2 = NULL; 503 size_t nrealisations1 = 0; 504 size_t nrealisations2 = 0; 505 struct sdis_mc T1 = SDIS_MC_NULL; 506 struct sdis_mc T2 = SDIS_MC_NULL; 507 double sum = 0; /* Sum of weights */ 508 double sum2 = 0; /* Sum of squared weights */ 509 double E = 0; /* Expected value */ 510 double V = 0; /* Variance */ 511 double SE = 0; /* Standard Error */ 512 double N = 0; /* Number of realisations */ 513 double ref = 0; /* Analytical solution */ 514 515 /* Retrieve the estimation results */ 516 OK(sdis_estimator_buffer_at(estim_buf1, iprobe, 0, &estimator1)); 517 OK(sdis_estimator_buffer_at(estim_buf2, iprobe, 0, &estimator2)); 518 OK(sdis_estimator_get_temperature(estimator1, &T1)); 519 OK(sdis_estimator_get_temperature(estimator2, &T2)); 520 521 /* Estimates are not identical... */ 522 CHK(T1.E != T2.E); 523 CHK(T1.SE != T2.SE); 524 525 /* ... but are compatible ... */ 526 check_intersection(T1.E, 3*T1.SE, T2.E, 3*T2.SE); 527 528 /* We can combine the 2 estimates since they must be numerically 529 * independent. We can therefore verify that their combination is valid with 530 * respect to the analytical solution */ 531 OK(sdis_estimator_get_realisation_count(estimator1, &nrealisations1)); 532 OK(sdis_estimator_get_realisation_count(estimator2, &nrealisations2)); 533 534 sum = 535 T1.E*(double)nrealisations1 536 + T2.E*(double)nrealisations2; 537 sum2 = 538 (T1.V + T1.E*T1.E)*(double)nrealisations1 539 + (T2.V + T2.E*T2.E)*(double)nrealisations2; 540 N = (double)(nrealisations1 + nrealisations2); 541 E = sum / N; 542 V = sum2 / N - E*E; 543 SE = sqrt(V/N); 544 545 ref = trilinear_profile(args->probes[iprobe].position); 546 printf("T(%g, %g, %g) = %g ~ %g +/- %g\n", 547 SPLIT3(args->probes[iprobe].position), ref, E, SE); 548 CHK(eq_eps(ref, E, 3*SE)); 549 } 550 } 551 552 static void 553 check_probe_list 554 (struct sdis_scene* scn, 555 struct s3d_scene_view* view, 556 const int is_master_process) 557 { 558 #define NPROBES 10 559 560 /* RNG state */ 561 const char* rng_state_filename = "rng_state"; 562 struct ssp_rng* rng = NULL; 563 enum ssp_rng_type rng_type = SSP_RNG_TYPE_NULL; 564 565 /* Probe variables */ 566 struct sdis_solve_probe_args probes[NPROBES]; 567 struct sdis_solve_probe_list_args args = SDIS_SOLVE_PROBE_LIST_ARGS_DEFAULT; 568 size_t iprobe = 0; 569 570 /* Estimations */ 571 struct sdis_estimator_buffer* estim_buf = NULL; 572 struct sdis_estimator_buffer* estim_buf2 = NULL; 573 574 /* Misc */ 575 double delta = 0; 576 FILE* fp = NULL; 577 578 delta = view_compute_delta(view); 579 580 /* Setup the list of probes to calculate */ 581 args.probes = probes; 582 args.nprobes = NPROBES; 583 FOR_EACH(iprobe, 0, NPROBES) { 584 probes[iprobe] = SDIS_SOLVE_PROBE_ARGS_DEFAULT; 585 probes[iprobe].nrealisations = 10000; 586 view_sample_position(view, delta, probes[iprobe].position); 587 } 588 589 check_probe_list_api(scn, &args); 590 591 /* Solve the probes */ 592 OK(sdis_solve_probe_list(scn, &args, &estim_buf)); 593 594 if(!is_master_process) { 595 CHK(estim_buf == NULL); 596 } else { 597 check_estimator_buffer(estim_buf, &args); 598 599 /* Check the use of a non-default rng_state. Run the calculations using the 600 * final RNG state from the previous estimate. Thus, the estimates are 601 * statically independent and can be combined. To do this, write the RNG 602 * state to a file so that it is read by all processes for the next test. */ 603 OK(sdis_estimator_buffer_get_rng_state(estim_buf, &rng)); 604 OK(ssp_rng_get_type(rng, &rng_type)); 605 CHK(fp = fopen(rng_state_filename, "w")); 606 CHK(fwrite(&rng_type, sizeof(rng_type), 1, fp) == 1); 607 OK(ssp_rng_write(rng, fp)); 608 CHK(fclose(fp) == 0); 609 } 610 611 #ifdef SDIS_ENABLE_MPI 612 /* Synchronize processes. Wait for the master process to write RNG state to 613 * be used */ 614 { 615 struct sdis_device* sdis = NULL; 616 int is_mpi_used = 0; 617 OK(sdis_scene_get_device(scn, &sdis)); 618 OK(sdis_device_is_mpi_used(sdis, &is_mpi_used)); 619 if(is_mpi_used) { 620 CHK(MPI_Barrier(MPI_COMM_WORLD) == MPI_SUCCESS); 621 } 622 } 623 #endif 624 625 /* Read the RNG state */ 626 CHK(fp = fopen(rng_state_filename, "r")); 627 CHK(fread(&rng_type, sizeof(rng_type), 1, fp) == 1); 628 OK(ssp_rng_create(NULL, rng_type, &rng)); 629 OK(ssp_rng_read(rng, fp)); 630 CHK(fclose(fp) == 0); 631 632 /* Run a new calculation using the read RNG state */ 633 args.rng_state = rng; 634 OK(sdis_solve_probe_list(scn, &args, &estim_buf2)); 635 OK(ssp_rng_ref_put(rng)); 636 637 if(!is_master_process) { 638 CHK(estim_buf2 == NULL); 639 } else { 640 check_estimator_buffer_combination(estim_buf, estim_buf2, &args); 641 OK(sdis_estimator_buffer_ref_put(estim_buf)); 642 OK(sdis_estimator_buffer_ref_put(estim_buf2)); 643 } 644 645 #undef NPROBES 646 } 647 648 /******************************************************************************* 649 * The test 650 ******************************************************************************/ 651 int 652 main(int argc, char** argv) 653 { 654 /* Stardis */ 655 struct sdis_device* sdis = NULL; 656 struct sdis_interface* interf = NULL; 657 struct sdis_medium* solid = NULL; 658 struct sdis_medium* dummy = NULL; /* Medium surrounding the solid */ 659 struct sdis_scene* scn = NULL; 660 661 /* Miscellaneous */ 662 struct s3d_scene_view* view = NULL; 663 struct s3dut_mesh* super_shape = NULL; 664 int is_master_process = 0; 665 666 (void)argc, (void)argv; /* Avoid the "unused variable" warning */ 667 668 create_default_device(&argc, &argv, &is_master_process, &sdis); 669 670 super_shape = create_super_shape(); 671 view = create_view(super_shape); 672 673 solid = create_solid(sdis, view); 674 dummy = create_dummy(sdis); 675 interf = create_interface(sdis, solid, dummy); 676 scn = create_scene(sdis, super_shape, interf); 677 678 check_probe_list(scn, view, is_master_process); 679 680 OK(s3dut_mesh_ref_put(super_shape)); 681 OK(s3d_scene_view_ref_put(view)); 682 OK(sdis_medium_ref_put(solid)); 683 OK(sdis_medium_ref_put(dummy)); 684 OK(sdis_interface_ref_put(interf)); 685 OK(sdis_scene_ref_put(scn)); 686 687 free_default_device(sdis); 688 689 CHK(mem_allocated_size() == 0); 690 691 return 0; 692 }