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