test_sdis_solve_probe_boundary_list.c (16868B)
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 18 #include "test_sdis_utils.h" 19 #include "test_sdis_mesh.h" 20 21 #include <star/s3dut.h> 22 #include <rsys/math.h> 23 24 #include <math.h> 25 26 /* 27 * The system is a trilinear profile of the temperature at steady state, i.e. at 28 * each point of the system we can calculate the temperature analytically. Two 29 * forms are immersed in this temperature field: a super shape and a sphere 30 * included in the super shape. On the Monte Carlo side, the temperature is 31 * unknown everywhere except on the surface of the super shape whose 32 * temperature is defined from the aformentionned trilinear profile. We will 33 * estimate the temperature at the sphere boundary at several probe points. We 34 * should find by Monte Carlo the temperature of the trilinear profile at the 35 * position of the probe. It's the test. 36 * 37 * /\ <-- T(x,y,z) 38 * ___/ \___ 39 * T(z) \ __ / 40 * | T(y) \ / \ / 41 * |/ T=? *__/ \ 42 * o--- T(x) /_ __ _\ 43 * \/ \/ 44 */ 45 46 /******************************************************************************* 47 * Helper functions 48 ******************************************************************************/ 49 static double 50 trilinear_profile(const double pos[3]) 51 { 52 /* Range in X, Y and Z in which the trilinear profile is defined */ 53 const double lower = -4; 54 const double upper = +4; 55 56 /* Upper temperature limit in X, Y and Z [K]. Lower temperature limit is 57 * implicitly 0 */ 58 const double a = 333; /* Upper temperature limit in X [K] */ 59 const double b = 432; /* Upper temperature limit in Y [K] */ 60 const double c = 579; /* Upper temperature limit in Z [K] */ 61 62 double x, y, z; 63 64 /* Check pre-conditions */ 65 CHK(pos); 66 CHK(lower <= pos[0] && pos[0] <= upper); 67 CHK(lower <= pos[1] && pos[1] <= upper); 68 CHK(lower <= pos[2] && pos[2] <= upper); 69 70 x = (pos[0] - lower) / (upper - lower); 71 y = (pos[1] - lower) / (upper - lower); 72 z = (pos[2] - lower) / (upper - lower); 73 return a*x + b*y + c*z; 74 } 75 76 static INLINE float 77 rand_canonic(void) 78 { 79 return (float)rand() / (float)((unsigned)RAND_MAX+1); 80 } 81 82 static void 83 sample_sphere(double pos[3]) 84 { 85 const double phi = rand_canonic() * 2 * PI; 86 const double v = rand_canonic(); 87 const double cos_theta = 1 - 2 * v; 88 const double sin_theta = 2 * sqrt(v * (1 - v)); 89 pos[0] = cos(phi) * sin_theta; 90 pos[1] = sin(phi) * sin_theta; 91 pos[2] = cos_theta; 92 } 93 94 static INLINE void 95 check_intersection 96 (const double val0, 97 const double eps0, 98 const double val1, 99 const double eps1) 100 { 101 double interval0[2], interval1[2]; 102 double intersection[2]; 103 interval0[0] = val0 - eps0; 104 interval0[1] = val0 + eps0; 105 interval1[0] = val1 - eps1; 106 interval1[1] = val1 + eps1; 107 intersection[0] = MMAX(interval0[0], interval1[0]); 108 intersection[1] = MMIN(interval0[1], interval1[1]); 109 CHK(intersection[0] <= intersection[1]); 110 } 111 112 static void 113 mesh_add_super_shape(struct mesh* mesh) 114 { 115 struct s3dut_mesh* sshape = NULL; 116 struct s3dut_mesh_data sshape_data; 117 struct s3dut_super_formula f0 = S3DUT_SUPER_FORMULA_NULL; 118 struct s3dut_super_formula f1 = S3DUT_SUPER_FORMULA_NULL; 119 const double radius = 2; 120 const unsigned nslices = 256; 121 122 f0.A = 1.5; f0.B = 1; f0.M = 11.0; f0.N0 = 1; f0.N1 = 1; f0.N2 = 2.0; 123 f1.A = 1.0; f1.B = 2; f1.M = 3.6; f1.N0 = 1; f1.N1 = 2; f1.N2 = 0.7; 124 OK(s3dut_create_super_shape(NULL, &f0, &f1, radius, nslices, nslices/2, &sshape)); 125 OK(s3dut_mesh_get_data(sshape, &sshape_data)); 126 mesh_append(mesh, sshape_data.positions, sshape_data.nvertices, 127 sshape_data.indices, sshape_data.nprimitives, NULL); 128 OK(s3dut_mesh_ref_put(sshape)); 129 } 130 131 static void 132 mesh_add_sphere(struct mesh* mesh) 133 { 134 struct s3dut_mesh* sphere = NULL; 135 struct s3dut_mesh_data sphere_data; 136 const double radius = 1; 137 const unsigned nslices = 128; 138 139 OK(s3dut_create_sphere(NULL, radius, nslices, nslices/2, &sphere)); 140 OK(s3dut_mesh_get_data(sphere, &sphere_data)); 141 mesh_append(mesh, sphere_data.positions, sphere_data.nvertices, 142 sphere_data.indices, sphere_data.nprimitives, NULL); 143 OK(s3dut_mesh_ref_put(sphere)); 144 } 145 146 /******************************************************************************* 147 * Solid, i.e. medium of super shape and sphere 148 ******************************************************************************/ 149 #define SOLID_PROP(Prop, Val) \ 150 static double \ 151 solid_get_##Prop \ 152 (const struct sdis_rwalk_vertex* vtx, \ 153 struct sdis_data* data) \ 154 { \ 155 (void)vtx, (void)data; /* Avoid the "unused variable" warning */ \ 156 return Val; \ 157 } 158 SOLID_PROP(calorific_capacity, 500.0) /* [J/K/Kg] */ 159 SOLID_PROP(thermal_conductivity, 25.0) /* [W/m/K] */ 160 SOLID_PROP(volumic_mass, 7500.0) /* [kg/m^3] */ 161 SOLID_PROP(temperature, SDIS_TEMPERATURE_NONE/*<=> unknown*/) /* [K] */ 162 SOLID_PROP(delta, 1.0/20.0) /* [m] */ 163 164 static struct sdis_medium* 165 create_solid(struct sdis_device* sdis) 166 { 167 struct sdis_solid_shader shader = SDIS_SOLID_SHADER_NULL; 168 struct sdis_medium* solid = NULL; 169 170 shader.calorific_capacity = solid_get_calorific_capacity; 171 shader.thermal_conductivity = solid_get_thermal_conductivity; 172 shader.volumic_mass = solid_get_volumic_mass; 173 shader.delta = solid_get_delta; 174 shader.temperature = solid_get_temperature; 175 OK(sdis_solid_create(sdis, &shader, NULL, &solid)); 176 return solid; 177 } 178 179 /******************************************************************************* 180 * Dummy environment, i.e. environment surrounding the super shape. It is 181 * defined only for Stardis compliance: in Stardis, an interface must divide 2 182 * media. 183 ******************************************************************************/ 184 static struct sdis_medium* 185 create_dummy(struct sdis_device* sdis) 186 { 187 struct sdis_fluid_shader shader = SDIS_FLUID_SHADER_NULL; 188 struct sdis_medium* dummy = NULL; 189 190 shader.calorific_capacity = dummy_medium_getter; 191 shader.volumic_mass = dummy_medium_getter; 192 shader.temperature = dummy_medium_getter; 193 OK(sdis_fluid_create(sdis, &shader, NULL, &dummy)); 194 return dummy; 195 } 196 197 /******************************************************************************* 198 * Interfaces 199 ******************************************************************************/ 200 static double 201 interface_get_temperature 202 (const struct sdis_interface_fragment* frag, 203 struct sdis_data* data) 204 { 205 (void)data; /* Avoid the "unused variable" warning */ 206 return trilinear_profile(frag->P); 207 } 208 209 static struct sdis_interface* 210 create_interface 211 (struct sdis_device* sdis, 212 struct sdis_medium* front, 213 struct sdis_medium* back) 214 { 215 struct sdis_interface* interf = NULL; 216 struct sdis_interface_shader shader = SDIS_INTERFACE_SHADER_NULL; 217 218 /* Fluid/solid interface: fix the temperature to the temperature profile */ 219 if(sdis_medium_get_type(front) != sdis_medium_get_type(back)) { 220 shader.front.temperature = interface_get_temperature; 221 shader.back.temperature = interface_get_temperature; 222 } 223 224 OK(sdis_interface_create(sdis, front, back, &shader, NULL, &interf)); 225 return interf; 226 } 227 228 /******************************************************************************* 229 * Scene, i.e. the system to simulate 230 ******************************************************************************/ 231 struct scene_context { 232 const struct mesh* mesh; 233 size_t sshape_end_id; /* Last triangle index of the super shape */ 234 struct sdis_interface* sshape; 235 struct sdis_interface* sphere; 236 }; 237 #define SCENE_CONTEXT_NULL__ {NULL, 0, 0, 0} 238 static const struct scene_context SCENE_CONTEXT_NULL = SCENE_CONTEXT_NULL__; 239 240 static void 241 scene_get_indices(const size_t itri, size_t ids[3], void* ctx) 242 { 243 struct scene_context* context = ctx; 244 CHK(ids && context && itri < mesh_ntriangles(context->mesh)); 245 /* Flip the indices to ensure that the normal points into the super shape */ 246 ids[0] = (unsigned)context->mesh->indices[itri*3+0]; 247 ids[1] = (unsigned)context->mesh->indices[itri*3+2]; 248 ids[2] = (unsigned)context->mesh->indices[itri*3+1]; 249 } 250 251 static void 252 scene_get_interface(const size_t itri, struct sdis_interface** interf, void* ctx) 253 { 254 struct scene_context* context = ctx; 255 CHK(interf && context && itri < mesh_ntriangles(context->mesh)); 256 if(itri < context->sshape_end_id) { 257 *interf = context->sshape; 258 } else { 259 *interf = context->sphere; 260 } 261 } 262 263 static void 264 scene_get_position(const size_t ivert, double pos[3], void* ctx) 265 { 266 struct scene_context* context = ctx; 267 CHK(pos && context && ivert < mesh_nvertices(context->mesh)); 268 pos[0] = context->mesh->positions[ivert*3+0]; 269 pos[1] = context->mesh->positions[ivert*3+1]; 270 pos[2] = context->mesh->positions[ivert*3+2]; 271 } 272 273 static struct sdis_scene* 274 create_scene(struct sdis_device* sdis, struct scene_context* ctx) 275 276 { 277 struct sdis_scene* scn = NULL; 278 struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT; 279 280 scn_args.get_indices = scene_get_indices; 281 scn_args.get_interface = scene_get_interface; 282 scn_args.get_position = scene_get_position; 283 scn_args.nprimitives = mesh_ntriangles(ctx->mesh); 284 scn_args.nvertices = mesh_nvertices(ctx->mesh); 285 scn_args.context = ctx; 286 OK(sdis_scene_create(sdis, &scn_args, &scn)); 287 return scn; 288 } 289 290 /******************************************************************************* 291 * Validations 292 ******************************************************************************/ 293 static void 294 check_probe_boundary_list_api 295 (struct sdis_scene* scn, 296 const struct sdis_solve_probe_boundary_list_args* in_args) 297 { 298 struct sdis_solve_probe_boundary_list_args args = *in_args; 299 struct sdis_estimator_buffer* estim_buf = NULL; 300 301 /* Check API */ 302 BA(sdis_solve_probe_boundary_list(NULL, &args, &estim_buf)); 303 BA(sdis_solve_probe_boundary_list(scn, NULL, &estim_buf)); 304 BA(sdis_solve_probe_boundary_list(scn, &args, NULL)); 305 args.nprobes = 0; 306 BA(sdis_solve_probe_boundary_list(scn, &args, &estim_buf)); 307 args.nprobes = in_args->nprobes; 308 args.probes = NULL; 309 BA(sdis_solve_probe_boundary_list(scn, &args, &estim_buf)); 310 } 311 312 /* Check the estimators against the analytical solution */ 313 static void 314 check_estimator_buffer 315 (const struct sdis_estimator_buffer* estim_buf, 316 struct sdis_scene* scn, 317 const struct sdis_solve_probe_boundary_list_args* args) 318 { 319 struct sdis_mc T = SDIS_MC_NULL; 320 size_t iprobe = 0; 321 322 /* Variables used to check the global estimation results */ 323 size_t total_nrealisations = 0; 324 size_t total_nrealisations_ref = 0; 325 size_t total_nfailures = 0; 326 size_t total_nfailures_ref = 0; 327 double sum = 0; /* Global sum of weights */ 328 double sum2 = 0; /* Global sum of squared weights */ 329 double N = 0; /* Number of (successful) realisations */ 330 double E = 0; /* Expected value */ 331 double V = 0; /* Variance */ 332 double SE = 0; /* Standard Error */ 333 334 /* Check the results */ 335 FOR_EACH(iprobe, 0, args->nprobes) { 336 const struct sdis_estimator* estimator = NULL; 337 size_t probe_nrealisations = 0; 338 size_t probe_nfailures = 0; 339 double pos[3] = {0,0,0}; 340 double probe_sum = 0; 341 double probe_sum2 = 0; 342 double ref = 0; 343 344 OK(sdis_scene_get_boundary_position(scn, args->probes[iprobe].iprim, 345 args->probes[iprobe].uv, pos)); 346 347 /* Fetch result */ 348 OK(sdis_estimator_buffer_at(estim_buf, iprobe, 0, &estimator)); 349 OK(sdis_estimator_get_temperature(estimator, &T)); 350 OK(sdis_estimator_get_realisation_count(estimator, &probe_nrealisations)); 351 OK(sdis_estimator_get_failure_count(estimator, &probe_nfailures)); 352 353 /* Check probe estimation */ 354 ref = trilinear_profile(pos); 355 printf("T(%g, %g, %g) = %g ~ %g +/- %g\n", SPLIT3(pos), ref, T.E, T.SE); 356 CHK(eq_eps(ref, T.E, 3*T.SE)); 357 358 /* Check miscellaneous results */ 359 CHK(probe_nrealisations == args->probes[iprobe].nrealisations); 360 361 /* Compute the sum of weights and sum of squared weights of the probe */ 362 probe_sum = T.E * (double)probe_nrealisations; 363 probe_sum2 = (T.V + T.E*T.E) * (double)probe_nrealisations; 364 365 /* Update the global estimation */ 366 total_nrealisations_ref += probe_nrealisations; 367 total_nfailures_ref += probe_nfailures; 368 sum += probe_sum; 369 sum2 += probe_sum2; 370 } 371 372 /* Check the overall estimate of the estimate buffer. This estimate is not 373 * really a result expected by the caller since the probes are all 374 * independent. But to be consistent with the estimate buffer API, we need to 375 * provide it. And so we check its validity */ 376 OK(sdis_estimator_buffer_get_temperature(estim_buf, &T)); 377 OK(sdis_estimator_buffer_get_realisation_count(estim_buf, &total_nrealisations)); 378 OK(sdis_estimator_buffer_get_failure_count(estim_buf, &total_nfailures)); 379 380 CHK(total_nrealisations == total_nrealisations_ref); 381 CHK(total_nfailures == total_nfailures_ref); 382 383 N = (double)total_nrealisations_ref; 384 E = sum / N; 385 V = sum2 / N - E*E; 386 SE = sqrt(V/N); 387 check_intersection(E, SE, T.E, T.SE); 388 } 389 390 static void 391 check_probe_boundary_list(struct sdis_scene* scn, const int is_master_process) 392 { 393 #define NPROBES 10 394 395 /* Estimations */ 396 struct sdis_estimator_buffer* estim_buf = NULL; 397 398 /* Probe variables */ 399 struct sdis_solve_probe_boundary_args probes[NPROBES]; 400 struct sdis_solve_probe_boundary_list_args args = 401 SDIS_SOLVE_PROBE_BOUNDARY_LIST_ARGS_DEFAULT; 402 size_t iprobe; 403 404 /* Miscellaneous */ 405 struct sdis_scene_find_closest_point_args closest_pt_args = 406 SDIS_SCENE_FIND_CLOSEST_POINT_ARGS_NULL; 407 408 (void)is_master_process; 409 410 /* Setup the list of probes to calculate */ 411 args.probes = probes; 412 args.nprobes = NPROBES; 413 FOR_EACH(iprobe, 0, NPROBES) { 414 sample_sphere(closest_pt_args.position); 415 closest_pt_args.radius = INF; 416 417 probes[iprobe] = SDIS_SOLVE_PROBE_BOUNDARY_ARGS_DEFAULT; 418 probes[iprobe].nrealisations = 10000; 419 probes[iprobe].side = SDIS_FRONT; 420 421 OK(sdis_scene_find_closest_point 422 (scn, &closest_pt_args, &probes[iprobe].iprim, probes[iprobe].uv)); 423 } 424 425 check_probe_boundary_list_api(scn, &args); 426 427 /* Solve the probes */ 428 OK(sdis_solve_probe_boundary_list(scn, &args, &estim_buf)); 429 if(!is_master_process) { 430 CHK(estim_buf == NULL); 431 } else { 432 check_estimator_buffer(estim_buf, scn, &args); 433 OK(sdis_estimator_buffer_ref_put(estim_buf)); 434 } 435 436 #undef NPROBES 437 } 438 439 /******************************************************************************* 440 * The test 441 ******************************************************************************/ 442 int 443 main(int argc, char** argv) 444 { 445 /* Stardis */ 446 struct sdis_device* dev = NULL; 447 struct sdis_interface* solid_fluid = NULL; 448 struct sdis_interface* solid_solid = NULL; 449 struct sdis_medium* fluid = NULL; 450 struct sdis_medium* solid = NULL; 451 struct sdis_scene* scn = NULL; 452 453 /* Miscellaneous */ 454 struct scene_context ctx = SCENE_CONTEXT_NULL; 455 struct mesh mesh = MESH_NULL; 456 size_t sshape_end_id = 0; /* Last index of the super shape */ 457 int is_master_process = 1; 458 (void)argc, (void)argv; 459 460 create_default_device(&argc, &argv, &is_master_process, &dev); 461 462 /* Setup the mesh */ 463 mesh_init(&mesh); 464 mesh_add_super_shape(&mesh); 465 sshape_end_id = mesh_ntriangles(&mesh); 466 mesh_add_sphere(&mesh); 467 468 /* Setup physical properties */ 469 fluid = create_dummy(dev); 470 solid = create_solid(dev); 471 solid_fluid = create_interface(dev, solid, fluid); 472 solid_solid = create_interface(dev, solid, solid); 473 474 /* Create the scene */ 475 ctx.mesh = &mesh; 476 ctx.sshape_end_id = sshape_end_id; 477 ctx.sshape = solid_fluid; 478 ctx.sphere = solid_solid; 479 scn = create_scene(dev, &ctx); 480 481 check_probe_boundary_list(scn, is_master_process); 482 483 mesh_release(&mesh); 484 OK(sdis_interface_ref_put(solid_fluid)); 485 OK(sdis_interface_ref_put(solid_solid)); 486 OK(sdis_medium_ref_put(fluid)); 487 OK(sdis_medium_ref_put(solid)); 488 OK(sdis_scene_ref_put(scn)); 489 490 free_default_device(dev); 491 492 CHK(mem_allocated_size() == 0); 493 return 0; 494 }