test_sdis_custom_solid_path_sampling.c (19534B)
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_mesh.h" 18 #include "test_sdis_utils.h" 19 20 #include <star/s3d.h> 21 #include <star/s3dut.h> 22 #include <star/ssp.h> 23 24 #include <rsys/double3.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. 33 * 34 * We will estimate the temperature at the position of a probe in solids by 35 * providing a user-side function to sample the conductive path in the sphere. 36 * We should find the temperature of the trilinear profile at the probe position 37 * by Monte Carlo, independently of this coupling with an external path sampling 38 * routine. 39 * 40 * 41 * /\ <-- T(x,y,z) 42 * ___/ \___ 43 * T(z) \ __ / 44 * | T(y) T=?/. \ / 45 * |/ / \__/ \ 46 * o--- T(x) /_ __ _\ 47 * \/ \/ 48 */ 49 50 #define NREALISATIONS 10000 51 #define SPHERE_RADIUS 1 52 53 /******************************************************************************* 54 * Helper functions 55 ******************************************************************************/ 56 static double 57 trilinear_profile(const double pos[3]) 58 { 59 /* Range in X, Y and Z in which the trilinear profile is defined */ 60 const double lower = -4; 61 const double upper = +4; 62 63 /* Upper temperature limit in X, Y and Z [K]. Lower temperature limit is 64 * implicitly 0 */ 65 const double a = 333; /* Upper temperature limit in X [K] */ 66 const double b = 432; /* Upper temperature limit in Y [K] */ 67 const double c = 579; /* Upper temperature limit in Z [K] */ 68 69 double x, y, z; 70 71 /* Check pre-conditions */ 72 CHK(pos); 73 CHK(lower <= pos[0] && pos[0] <= upper); 74 CHK(lower <= pos[1] && pos[1] <= upper); 75 CHK(lower <= pos[2] && pos[2] <= upper); 76 77 x = (pos[0] - lower) / (upper - lower); 78 y = (pos[1] - lower) / (upper - lower); 79 z = (pos[2] - lower) / (upper - lower); 80 return a*x + b*y + c*z; 81 } 82 83 /******************************************************************************* 84 * Scene view 85 ******************************************************************************/ 86 /* Parameter for get_inidices/get_vertices functions. Although its layout is the 87 * same as that of the mesh data structure, we have deliberately defined a new 88 * data type since mesh functions cannot be invoked. This is because the form's 89 * member variables are not necessarily extensible arrays, as is the case with 90 * mesh */ 91 struct shape { 92 double* pos; 93 size_t* ids; 94 size_t npos; 95 size_t ntri; 96 }; 97 #define SHAPE_NULL__ {NULL, NULL, 0, 0} 98 static const struct shape SHAPE_NULL = SHAPE_NULL__; 99 100 static void 101 get_position(const unsigned ivert, float pos[3], void* ctx) 102 { 103 const struct shape* shape = ctx; 104 CHK(shape && pos && ivert < shape->npos); 105 f3_set_d3(pos, shape->pos + ivert*3); 106 } 107 108 static void 109 get_indices(const unsigned itri, unsigned ids[3], void* ctx) 110 { 111 const struct shape* shape = ctx; 112 CHK(shape && ids && itri < shape->ntri); 113 ids[0] = (unsigned)shape->ids[itri*3+0]; 114 ids[1] = (unsigned)shape->ids[itri*3+1]; 115 ids[2] = (unsigned)shape->ids[itri*3+2]; 116 } 117 118 static struct s3d_scene_view* 119 create_view(struct shape* shape_data) 120 { 121 /* Star3D */ 122 struct s3d_vertex_data vdata = S3D_VERTEX_DATA_NULL; 123 struct s3d_device* dev = NULL; 124 struct s3d_scene* scn = NULL; 125 struct s3d_shape* shape = NULL; 126 struct s3d_scene_view* view = NULL; 127 128 OK(s3d_device_create(NULL, NULL, 0, &dev)); 129 130 /* Shape */ 131 vdata.usage = S3D_POSITION; 132 vdata.type = S3D_FLOAT3; 133 vdata.get = get_position; 134 OK(s3d_shape_create_mesh(dev, &shape)); 135 OK(s3d_mesh_setup_indexed_vertices(shape, (unsigned)shape_data->ntri, 136 get_indices, (unsigned)shape_data->npos, &vdata, 1, shape_data)); 137 138 /* Scene view */ 139 OK(s3d_scene_create(dev, &scn)); 140 OK(s3d_scene_attach_shape(scn, shape)); 141 OK(s3d_scene_view_create(scn, S3D_TRACE|S3D_GET_PRIMITIVE, &view)); 142 143 /* Clean up */ 144 OK(s3d_device_ref_put(dev)); 145 OK(s3d_scene_ref_put(scn)); 146 OK(s3d_shape_ref_put(shape)); 147 148 return view; 149 } 150 151 /******************************************************************************* 152 * Mesh, i.e. supershape and sphere 153 ******************************************************************************/ 154 static void 155 mesh_add_super_shape(struct mesh* mesh) 156 { 157 struct s3dut_mesh* sshape = NULL; 158 struct s3dut_mesh_data sshape_data; 159 struct s3dut_super_formula f0 = S3DUT_SUPER_FORMULA_NULL; 160 struct s3dut_super_formula f1 = S3DUT_SUPER_FORMULA_NULL; 161 const double radius = 2; 162 const unsigned nslices = 256; 163 164 f0.A = 1.5; f0.B = 1; f0.M = 11.0; f0.N0 = 1; f0.N1 = 1; f0.N2 = 2.0; 165 f1.A = 1.0; f1.B = 2; f1.M = 3.6; f1.N0 = 1; f1.N1 = 2; f1.N2 = 0.7; 166 OK(s3dut_create_super_shape(NULL, &f0, &f1, radius, nslices, nslices/2, &sshape)); 167 OK(s3dut_mesh_get_data(sshape, &sshape_data)); 168 mesh_append(mesh, sshape_data.positions, sshape_data.nvertices, 169 sshape_data.indices, sshape_data.nprimitives, NULL); 170 OK(s3dut_mesh_ref_put(sshape)); 171 } 172 173 static void 174 mesh_add_sphere(struct mesh* mesh) 175 { 176 struct s3dut_mesh* sphere = NULL; 177 struct s3dut_mesh_data sphere_data; 178 const double radius = SPHERE_RADIUS; 179 const unsigned nslices = 128; 180 181 OK(s3dut_create_sphere(NULL, radius, nslices, nslices/2, &sphere)); 182 OK(s3dut_mesh_get_data(sphere, &sphere_data)); 183 mesh_append(mesh, sphere_data.positions, sphere_data.nvertices, 184 sphere_data.indices, sphere_data.nprimitives, NULL); 185 OK(s3dut_mesh_ref_put(sphere)); 186 } 187 188 /******************************************************************************* 189 * Custom conductive path 190 ******************************************************************************/ 191 struct custom_solid { 192 struct s3d_scene_view* view; /* Star-3D view of the shape */ 193 const struct shape* shape; /* Raw data */ 194 }; 195 196 static void 197 setup_solver_primitive 198 (struct sdis_scene* scn, 199 const struct shape* shape, 200 const struct s3d_hit* user_hit, 201 struct s3d_primitive* prim) 202 { 203 struct sdis_primkey key = SDIS_PRIMKEY_NULL; 204 const double *v0, *v1, *v2; 205 float v0f[3], v1f[3], v2f[3]; 206 struct s3d_attrib attr0, attr1, attr2; 207 208 v0 = shape->pos + shape->ids[user_hit->prim.prim_id*3+0]*3; 209 v1 = shape->pos + shape->ids[user_hit->prim.prim_id*3+1]*3; 210 v2 = shape->pos + shape->ids[user_hit->prim.prim_id*3+2]*3; 211 sdis_primkey_setup(&key, v0, v1, v2); 212 OK(sdis_scene_get_s3d_primitive(scn, &key, prim)); 213 214 /* Check that the primitive on the solver side is the same as that on the 215 * user side. On the solver side, vertices are stored in simple precision in 216 * Star-3D view. We therefore need to take care of this conversion to check 217 * that the vertices are the same */ 218 OK(s3d_triangle_get_vertex_attrib(prim, 0, S3D_POSITION, &attr0)); 219 OK(s3d_triangle_get_vertex_attrib(prim, 1, S3D_POSITION, &attr1)); 220 OK(s3d_triangle_get_vertex_attrib(prim, 2, S3D_POSITION, &attr2)); 221 f3_set_d3(v0f, v0); 222 f3_set_d3(v1f, v1); 223 f3_set_d3(v2f, v2); 224 225 /* The vertices have been inverted on the user's side to reverse the normal 226 * orientation. Below it is taken into account */ 227 CHK(f3_eq(v0f, attr0.value)); 228 CHK(f3_eq(v1f, attr2.value)); 229 CHK(f3_eq(v2f, attr1.value)); 230 } 231 232 /* Implementation of a Walk on Sphere algorithm for an origin-centered spherical 233 * geometry of radius SPHERE_RADIUS */ 234 static res_T 235 sample_steady_diffusive_path 236 (struct sdis_scene* scn, 237 struct ssp_rng* rng, 238 struct sdis_path* path, 239 struct sdis_data* data) 240 { 241 struct custom_solid* solid = NULL; 242 struct s3d_hit hit = S3D_HIT_NULL; 243 244 double pos[3]; 245 const double epsilon = SPHERE_RADIUS * 1.e-6; /* Epsilon shell */ 246 247 CHK(scn && rng && path && data); 248 249 solid = sdis_data_get(data); 250 251 d3_set(pos, path->vtx.P); 252 253 do { 254 /* Distance from the geometry center to the current position */ 255 const double dst = d3_len(pos); 256 257 double dir[3] = {0,0,0}; 258 double r = 0; /* Radius */ 259 260 r = SPHERE_RADIUS - dst; 261 CHK(dst > 0); 262 263 if(r > epsilon) { 264 /* Uniformly sample a new position on the surrounding sphere */ 265 ssp_ran_sphere_uniform(rng, dir, NULL); 266 267 /* Move to the new position */ 268 d3_muld(dir, dir, r); 269 d3_add(pos, pos, dir); 270 271 /* The current position is in the epsilon shell: 272 * move it to the nearest interface position */ 273 } else { 274 float posf[3]; 275 276 d3_set(dir, pos); 277 d3_normalize(dir, dir); 278 d3_muld(pos, dir, SPHERE_RADIUS); 279 280 /* Map the position to the sphere geometry */ 281 f3_set_d3(posf, pos); 282 OK(s3d_scene_view_closest_point(solid->view, posf, (float)INF, NULL, &hit)); 283 } 284 285 /* The calculation is performed in steady state, so the path necessarily stops 286 * at a boundary */ 287 } while(S3D_HIT_NONE(&hit)); 288 289 /* Setup the path state */ 290 d3_set(path->vtx.P, pos); 291 path->weight = 0; 292 path->at_limit = 0; 293 path->prim_2d = S2D_PRIMITIVE_NULL; 294 setup_solver_primitive(scn, solid->shape, &hit, &path->prim_3d); 295 296 return RES_OK; 297 } 298 299 /******************************************************************************* 300 * The solids, i.e. media of the super shape and the sphere 301 ******************************************************************************/ 302 #define SOLID_PROP(Prop, Val) \ 303 static double \ 304 solid_get_##Prop \ 305 (const struct sdis_rwalk_vertex* vtx, \ 306 struct sdis_data* data) \ 307 { \ 308 (void)vtx, (void)data; /* Avoid the "unused variable" warning */ \ 309 return Val; \ 310 } 311 SOLID_PROP(calorific_capacity, 500.0) /* [J/K/Kg] */ 312 SOLID_PROP(thermal_conductivity, 25.0) /* [W/m/K] */ 313 SOLID_PROP(volumic_mass, 7500.0) /* [kg/m^3] */ 314 SOLID_PROP(temperature, SDIS_TEMPERATURE_NONE/*<=> unknown*/) /* [K] */ 315 SOLID_PROP(delta, 1.0/40.0) /* [m] */ 316 317 static struct sdis_medium* 318 create_solid(struct sdis_device* sdis) 319 { 320 struct sdis_solid_shader shader = SDIS_SOLID_SHADER_NULL; 321 struct sdis_medium* solid = NULL; 322 323 shader.calorific_capacity = solid_get_calorific_capacity; 324 shader.thermal_conductivity = solid_get_thermal_conductivity; 325 shader.volumic_mass = solid_get_volumic_mass; 326 shader.delta = solid_get_delta; 327 shader.temperature = solid_get_temperature; 328 329 OK(sdis_solid_create(sdis, &shader, NULL, &solid)); 330 return solid; 331 } 332 333 static struct sdis_medium* 334 create_custom 335 (struct sdis_device* sdis, 336 struct s3d_scene_view* view, 337 const struct shape* shape) 338 { 339 /* Stardis variables */ 340 struct sdis_solid_shader shader = SDIS_SOLID_SHADER_NULL; 341 struct sdis_medium* solid = NULL; 342 struct sdis_data* data = NULL; 343 344 /* Mesh variables */ 345 struct custom_solid* custom_solid = NULL; 346 const size_t sz = sizeof(struct custom_solid); 347 const size_t al = ALIGNOF(struct custom_solid); 348 349 OK(sdis_data_create(sdis, sz, al, NULL, &data)); 350 custom_solid = sdis_data_get(data); 351 custom_solid->view = view; 352 custom_solid->shape = shape; 353 354 shader.calorific_capacity = solid_get_calorific_capacity; 355 shader.thermal_conductivity = solid_get_thermal_conductivity; 356 shader.volumic_mass = solid_get_volumic_mass; 357 shader.delta = solid_get_delta; 358 shader.temperature = solid_get_temperature; 359 shader.sample_path = sample_steady_diffusive_path; 360 361 OK(sdis_solid_create(sdis, &shader, data, &solid)); 362 OK(sdis_data_ref_put(data)); 363 364 return solid; 365 } 366 367 /******************************************************************************* 368 * Dummy environment, i.e. environment surrounding the super shape. It is 369 * defined only for Stardis compliance: in Stardis, an interface must divide 2 370 * media. 371 ******************************************************************************/ 372 static struct sdis_medium* 373 create_dummy(struct sdis_device* sdis) 374 { 375 struct sdis_fluid_shader shader = SDIS_FLUID_SHADER_NULL; 376 struct sdis_medium* dummy = NULL; 377 378 shader.calorific_capacity = dummy_medium_getter; 379 shader.volumic_mass = dummy_medium_getter; 380 shader.temperature = dummy_medium_getter; 381 OK(sdis_fluid_create(sdis, &shader, NULL, &dummy)); 382 return dummy; 383 } 384 385 /******************************************************************************* 386 * Interface: its temperature is fixed to the trilinear profile 387 ******************************************************************************/ 388 static double 389 interface_get_temperature 390 (const struct sdis_interface_fragment* frag, 391 struct sdis_data* data) 392 { 393 (void)data; /* Avoid the "unused variable" warning */ 394 return trilinear_profile(frag->P); 395 } 396 397 static struct sdis_interface* 398 create_interface 399 (struct sdis_device* sdis, 400 struct sdis_medium* front, 401 struct sdis_medium* back) 402 { 403 struct sdis_interface* interf = NULL; 404 struct sdis_interface_shader shader = SDIS_INTERFACE_SHADER_NULL; 405 406 /* Fluid/solid interface: fix the temperature to the temperature profile */ 407 if(sdis_medium_get_type(front) != sdis_medium_get_type(back)) { 408 shader.front.temperature = interface_get_temperature; 409 shader.back.temperature = interface_get_temperature; 410 } 411 412 OK(sdis_interface_create(sdis, front, back, &shader, NULL, &interf)); 413 return interf; 414 } 415 416 /******************************************************************************* 417 * Scene, i.e. the system to simulate 418 ******************************************************************************/ 419 struct scene_context { 420 const struct mesh* mesh; 421 size_t sshape_end_id; /* Last triangle index of the super shape */ 422 struct sdis_interface* sshape; 423 struct sdis_interface* sphere; 424 }; 425 #define SCENE_CONTEXT_NULL__ {NULL, 0, 0, 0} 426 static const struct scene_context SCENE_CONTEXT_NULL = SCENE_CONTEXT_NULL__; 427 428 static void 429 scene_get_indices(const size_t itri, size_t ids[3], void* ctx) 430 { 431 struct scene_context* context = ctx; 432 CHK(ids && context && itri < mesh_ntriangles(context->mesh)); 433 /* Flip the indices to ensure that the normal points into the super shape */ 434 ids[0] = (unsigned)context->mesh->indices[itri*3+0]; 435 ids[1] = (unsigned)context->mesh->indices[itri*3+2]; 436 ids[2] = (unsigned)context->mesh->indices[itri*3+1]; 437 } 438 439 static void 440 scene_get_interface(const size_t itri, struct sdis_interface** interf, void* ctx) 441 { 442 struct scene_context* context = ctx; 443 CHK(interf && context && itri < mesh_ntriangles(context->mesh)); 444 if(itri < context->sshape_end_id) { 445 *interf = context->sshape; 446 } else { 447 *interf = context->sphere; 448 } 449 } 450 451 static void 452 scene_get_position(const size_t ivert, double pos[3], void* ctx) 453 { 454 struct scene_context* context = ctx; 455 CHK(pos && context && ivert < mesh_nvertices(context->mesh)); 456 pos[0] = context->mesh->positions[ivert*3+0]; 457 pos[1] = context->mesh->positions[ivert*3+1]; 458 pos[2] = context->mesh->positions[ivert*3+2]; 459 } 460 461 static struct sdis_scene* 462 create_scene(struct sdis_device* sdis, struct scene_context* ctx) 463 { 464 struct sdis_scene* scn = NULL; 465 struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT; 466 467 scn_args.get_indices = scene_get_indices; 468 scn_args.get_interface = scene_get_interface; 469 scn_args.get_position = scene_get_position; 470 scn_args.nprimitives = mesh_ntriangles(ctx->mesh); 471 scn_args.nvertices = mesh_nvertices(ctx->mesh); 472 scn_args.context = ctx; 473 OK(sdis_scene_create(sdis, &scn_args, &scn)); 474 return scn; 475 } 476 477 /******************************************************************************* 478 * Validations 479 ******************************************************************************/ 480 static void 481 check_probe(struct sdis_scene* scn, const int is_master_process) 482 { 483 struct sdis_solve_probe_args args = SDIS_SOLVE_PROBE_ARGS_DEFAULT; 484 struct sdis_mc T = SDIS_MC_NULL; 485 struct sdis_estimator* estimator = NULL; 486 double ref = 0; 487 488 args.position[0] = SPHERE_RADIUS*0.125; 489 args.position[1] = SPHERE_RADIUS*0.250; 490 args.position[2] = SPHERE_RADIUS*0.375; 491 args.nrealisations = NREALISATIONS; 492 493 OK(sdis_solve_probe(scn, &args, &estimator)); 494 495 if(!is_master_process) return; 496 497 OK(sdis_estimator_get_temperature(estimator, &T)); 498 499 ref = trilinear_profile(args.position); 500 501 printf("T(%g, %g, %g) = %g ~ %g +/- %g\n", 502 SPLIT3(args.position), ref, T.E, T.SE); 503 504 CHK(eq_eps(ref, T.E, T.SE*3)); 505 OK(sdis_estimator_ref_put(estimator)); 506 } 507 508 /******************************************************************************* 509 * The test 510 ******************************************************************************/ 511 int 512 main(int argc, char** argv) 513 { 514 /* Stardis */ 515 struct sdis_device* dev = NULL; 516 struct sdis_interface* solid_dummy = NULL; 517 struct sdis_interface* custom_solid = NULL; 518 struct sdis_medium* solid = NULL; 519 struct sdis_medium* custom = NULL; 520 struct sdis_medium* dummy = NULL; /* Medium surrounding the solid */ 521 struct sdis_scene* scn = NULL; 522 523 /* Star3D */ 524 struct shape shape = SHAPE_NULL; 525 struct s3d_scene_view* sphere_view = NULL; 526 527 /* Miscellaneous */ 528 struct scene_context ctx = SCENE_CONTEXT_NULL; 529 struct mesh mesh = MESH_NULL; 530 size_t sshape_end_id = 0; /* Last index of the super shape */ 531 int is_master_process = 1; 532 (void)argc, (void)argv; 533 534 create_default_device(&argc, &argv, &is_master_process, &dev); 535 536 /* Mesh */ 537 mesh_init(&mesh); 538 mesh_add_super_shape(&mesh); 539 sshape_end_id = mesh_ntriangles(&mesh); 540 mesh_add_sphere(&mesh); 541 542 /* Create a view of the sphere's geometry. This will be used to couple custom 543 * solid path sampling to the solver */ 544 shape.pos = mesh.positions; 545 shape.ids = mesh.indices + sshape_end_id*3; 546 shape.npos = mesh_nvertices(&mesh); 547 shape.ntri = mesh_ntriangles(&mesh) - sshape_end_id/* #sshape triangles*/; 548 sphere_view = create_view(&shape); 549 550 /* Physical properties */ 551 dummy = create_dummy(dev); 552 solid = create_solid(dev); 553 custom = create_custom(dev, sphere_view, &shape); 554 solid_dummy = create_interface(dev, solid, dummy); 555 custom_solid = create_interface(dev, custom, solid); 556 557 /* Scene */ 558 ctx.mesh = &mesh; 559 ctx.sshape_end_id = sshape_end_id; 560 ctx.sshape = solid_dummy; 561 ctx.sphere = custom_solid; 562 scn = create_scene(dev, &ctx); 563 564 check_probe(scn, is_master_process); 565 566 mesh_release(&mesh); 567 568 OK(s3d_scene_view_ref_put(sphere_view)); 569 OK(sdis_interface_ref_put(solid_dummy)); 570 OK(sdis_interface_ref_put(custom_solid)); 571 OK(sdis_medium_ref_put(custom)); 572 OK(sdis_medium_ref_put(dummy)); 573 OK(sdis_medium_ref_put(solid)); 574 OK(sdis_scene_ref_put(scn)); 575 576 free_default_device(dev); 577 578 CHK(mem_allocated_size() == 0); 579 return 0; 580 }