test_sdis_transcient.c (23327B)
1 /* Copyright (C) 2016-2019 |Meso|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 <string.h> 20 21 /* 22 * The scene is composed of a solid cuboid whose temperature is fixed on its 6 23 * faces. The test consist in checking that the estimated temperature at a 24 * given temperature and time is compatible with the reference temperature 25 * computed by analytically evaluating the green function. 26 * 27 * The test is performed on 3 scenes that actually represent the same system. 28 * The first scene is simply the cuboid, as it. The second scene is the same 29 * cuboid but this time formed by 2 sub cuboid with strictly the same physical 30 * properties. Finally, the last scene is the same cube with successive 31 * cubes into it used only to add fictive boundaries. 32 */ 33 34 static const double vertices[12/*#vertices*/*3/*#coords per vertex*/] = { 35 0.0, 0.0, 0.0, 36 0.5, 0.0, 0.0, 37 0.0, 1.0, 0.0, 38 0.5, 1.0, 0.0, 39 0.0, 0.0, 1.0, 40 0.5, 0.0, 1.0, 41 0.0, 1.0, 1.0, 42 0.5, 1.0, 1.0, 43 1.0, 0.0, 0.0, 44 1.0, 1.0, 0.0, 45 1.0, 0.0, 1.0, 46 1.0, 1.0, 1.0 47 }; 48 static const size_t nvertices = sizeof(vertices) / (sizeof(double)*3); 49 50 /* The following array lists the indices toward the 3D vertices of each 51 * triangle. 52 * ,2---,3 ,3---,9 | ,2----3 ,3----9 53 * ,' | ,'/| ,' | ,'/| | ,'/| \ | ,'/| \ | Y 54 * 6----7' / | 7---11' / | | 6' / | \ | 7' / | \ | | 55 * |', | / ,1 |', | / ,8 | | / ,0---,1 | / ,1---,8 o--X 56 * | ',|/,' | ',|/,' | |/,' | ,' |/,' | ,' / 57 * 4----5 5---10 | 4----5' 5---10' Z 58 */ 59 static const size_t indices[22/*#triangles*/*3/*#indices per triangle*/] = { 60 0, 4, 2, 2, 4, 6, /* X min */ 61 3, 7, 5, 5, 1, 3, /* X mid */ 62 9,11,10,10, 8, 9, /* X max */ 63 0, 5, 4, 0, 1, 5, 1,10, 5, 1, 8,10, /* Y min */ 64 2, 6, 7, 2, 7, 3, 3, 7,11, 3,11, 9, /* Y max */ 65 0, 2, 1, 1, 2, 3, 1, 3, 8, 8, 3, 9, /* Z min */ 66 4, 5, 6, 6, 5, 7, 5,10, 7, 7,10,11 /* Z max */ 67 }; 68 static const size_t ntriangles = sizeof(indices) / (sizeof(size_t)*3); 69 70 /******************************************************************************* 71 * Box geometry functions 72 ******************************************************************************/ 73 struct context { 74 const double* vertices; 75 const size_t* indices; 76 struct sdis_interface* interfs[22]; 77 const double* scale; 78 }; 79 80 static void 81 get_position(const size_t ivert, double pos[3], void* context) 82 { 83 struct context* ctx = context; 84 CHK(ctx); 85 pos[0] = ctx->vertices[ivert*3+0] * ctx->scale[0]; 86 pos[1] = ctx->vertices[ivert*3+1] * ctx->scale[1]; 87 pos[2] = ctx->vertices[ivert*3+2] * ctx->scale[2]; 88 } 89 90 static void 91 get_indices(const size_t itri, size_t ids[3], void* context) 92 { 93 struct context* ctx = context; 94 CHK(ctx); 95 ids[0] = ctx->indices[itri*3+0]; 96 ids[1] = ctx->indices[itri*3+1]; 97 ids[2] = ctx->indices[itri*3+2]; 98 } 99 100 static void 101 get_interface(const size_t itri, struct sdis_interface** bound, void* context) 102 { 103 struct context* ctx = context; 104 CHK(ctx); 105 *bound = ctx->interfs[itri]; 106 } 107 108 /******************************************************************************* 109 * Setup a scene composed of a box that successively contains smaller boxes 110 ******************************************************************************/ 111 struct matriochka_context { 112 struct sdis_interface* interfs[13]; 113 const double* scale; 114 size_t nboxes; 115 }; 116 117 static void 118 matriochka_position(const size_t ivert, double pos[3], void* context) 119 { 120 struct matriochka_context* ctx = context; 121 const size_t ibox = ctx->nboxes - ivert / box_nvertices - 1; 122 const double* verts = box_vertices; 123 const double box_szmin = 1.0 / (double)ctx->nboxes; 124 const size_t i = ivert % box_nvertices; 125 CHK(ibox <= ctx->nboxes); 126 pos[0] = ((verts[i*3+0]-0.5)*(double)(ibox+1)*box_szmin + 0.5)*ctx->scale[0]; 127 pos[1] = ((verts[i*3+1]-0.5)*(double)(ibox+1)*box_szmin + 0.5)*ctx->scale[1]; 128 pos[2] = ((verts[i*3+2]-0.5)*(double)(ibox+1)*box_szmin + 0.5)*ctx->scale[2]; 129 } 130 131 static void 132 matriochka_indices(const size_t itri, size_t ids[3], void* context) 133 { 134 struct matriochka_context* ctx = context; 135 const size_t i = itri % box_ntriangles; 136 const size_t ibox = ctx->nboxes - itri / box_ntriangles - 1; 137 CHK(ibox <= ctx->nboxes); 138 (void)context; 139 ids[0] = box_indices[i*3+0] + ibox*box_nvertices; 140 ids[1] = box_indices[i*3+1] + ibox*box_nvertices; 141 ids[2] = box_indices[i*3+2] + ibox*box_nvertices; 142 } 143 144 static void 145 matriochka_interface 146 (const size_t itri, struct sdis_interface** bound, void* context) 147 { 148 struct matriochka_context* ctx = context; 149 const size_t ibox = ctx->nboxes - itri / box_ntriangles - 1; 150 const size_t i = itri % box_ntriangles; 151 CHK(ibox < ctx->nboxes); 152 *bound = ibox != 0 ? ctx->interfs[12] : ctx->interfs[i]; 153 } 154 155 static INLINE void 156 dump_matriochkas(FILE* stream, struct matriochka_context* ctx) 157 { 158 size_t i; 159 ASSERT(ctx && stream); 160 FOR_EACH(i, 0, ctx->nboxes*box_nvertices) { 161 double pos[3]; 162 matriochka_position(i, pos, ctx); 163 fprintf(stream, "v %g %g %g\n", SPLIT3(pos)); 164 } 165 FOR_EACH(i, 0, ctx->nboxes*box_ntriangles) { 166 size_t ids[3]; 167 matriochka_indices(i, ids, ctx); 168 fprintf(stream, "f %lu %lu %lu\n", 169 (unsigned long)ids[0]+1, 170 (unsigned long)ids[1]+1, 171 (unsigned long)ids[2]+1); 172 } 173 } 174 175 /******************************************************************************* 176 * Solid medium 177 ******************************************************************************/ 178 struct solid { 179 double cp; 180 double lambda; 181 double rho; 182 double delta; 183 double init_temperature; 184 }; 185 186 static double 187 solid_get_calorific_capacity 188 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 189 { 190 CHK(data != NULL && vtx != NULL); 191 return ((const struct solid*)sdis_data_cget(data))->cp; 192 } 193 194 static double 195 solid_get_thermal_conductivity 196 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 197 { 198 CHK(data != NULL && vtx != NULL); 199 return ((const struct solid*)sdis_data_cget(data))->lambda; 200 } 201 202 static double 203 solid_get_volumic_mass 204 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 205 { 206 CHK(data != NULL && vtx != NULL); 207 return ((const struct solid*)sdis_data_cget(data))->rho; 208 } 209 210 static double 211 solid_get_delta 212 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 213 { 214 CHK(data != NULL && vtx != NULL); 215 return ((const struct solid*)sdis_data_cget(data))->delta; 216 } 217 218 static double 219 solid_get_temperature 220 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 221 { 222 CHK(data != NULL && vtx != NULL); 223 if(vtx->time <= 0) { 224 return ((const struct solid*)sdis_data_cget(data))->init_temperature; 225 } else { 226 return SDIS_TEMPERATURE_NONE; 227 } 228 } 229 230 /******************************************************************************* 231 * Interface 232 ******************************************************************************/ 233 struct interf { 234 double temperature; 235 }; 236 237 static double 238 interface_get_temperature 239 (const struct sdis_interface_fragment* frag, struct sdis_data* data) 240 { 241 const struct interf* interf = sdis_data_cget(data); 242 CHK(frag && data); 243 return interf->temperature; 244 } 245 246 static struct sdis_interface* 247 create_interface 248 (struct sdis_device* dev, 249 struct sdis_medium* front, 250 struct sdis_medium* back, 251 const struct sdis_interface_shader* interf_shader, 252 const double temperature) 253 { 254 struct sdis_data* data = NULL; 255 struct sdis_interface* interf = NULL; 256 struct interf* interf_props = NULL; 257 258 OK(sdis_data_create 259 (dev, sizeof(struct interf), ALIGNOF(struct interf), NULL, &data)); 260 interf_props = sdis_data_get(data); 261 interf_props->temperature = temperature; 262 OK(sdis_interface_create 263 (dev, front, back, interf_shader, data, &interf)); 264 OK(sdis_data_ref_put(data)); 265 return interf; 266 } 267 268 /******************************************************************************* 269 * Analytical solution 270 ******************************************************************************/ 271 static void 272 fourier_pq 273 (const size_t nterms_pq, 274 const double pos[3], 275 const double sz[3], 276 const int i0, 277 const int i1, 278 const int i2, 279 double green[2]) 280 { 281 size_t p, q; 282 CHK(green); 283 284 green[0] = 0; 285 green[1] = 0; 286 287 FOR_EACH(p, 0, nterms_pq+1) { 288 FOR_EACH(q, 0, nterms_pq+1) { 289 const double p2 = (double)(2*p + 1); 290 const double q2 = (double)(2*q + 1); 291 double L_sqr, L, tmp; 292 L_sqr = PI * PI * ((p2*p2)/(sz[i1]*sz[i1]) + (q2*q2)/(sz[i2]*sz[i2])); 293 L = sqrt(L_sqr); 294 tmp = sin(PI*p2*pos[i1]/sz[i1]) 295 * sin(PI*q2*pos[i2]/sz[i2]) 296 / (sinh(sz[i0]*L)*(p2*q2)); 297 if(tmp != 0) { 298 green[0] += sinh(L*(sz[i0]-pos[i0]))*tmp; 299 green[1] += sinh(L*pos[i0])*tmp; 300 } 301 } 302 } 303 } 304 305 /* This function computes the Green function between a given probe 306 * position/time in a parallelepipedic box and each face of this box (within 307 * the model of a homogeneous boundary condition on each face). */ 308 static void 309 green_analytical 310 (const double box_size[3], 311 const double probe[3], 312 const double time, 313 const double rho, 314 const double cp, 315 const double lambda, 316 double green[7]) 317 { 318 const size_t nterms_fs = 20; /* #terms in the Fourier expansion series */ 319 const size_t nterms_pq = 100; /* #terms in double p/q sums */ 320 const size_t nt_pq = (nterms_fs - 1)/2; 321 double Gs[7], Gi[7], Gtmp[7]; 322 size_t i, m, n, o, p, q; 323 double a, b, c; 324 double alpha; 325 326 CHK(box_size && probe && time >= 0 && green); 327 328 if(time == 0) { 329 memset(green, 0, sizeof(double[7])); 330 green[6] = 1; 331 return; 332 } 333 334 memset(Gs, 0, sizeof(double[7])); 335 memset(Gi, 0, sizeof(double[7])); 336 memset(Gtmp, 0, sizeof(double[7])); 337 338 /* Steady state solution */ 339 fourier_pq(nterms_pq, probe, box_size, 0, 1, 2, Gtmp+0); /* Faces 0 and 1 */ 340 fourier_pq(nterms_pq, probe, box_size, 1, 0, 2, Gtmp+2); /* Faces 2 and 3 */ 341 fourier_pq(nterms_pq, probe, box_size, 2, 1, 0, Gtmp+4); /* Faces 4 and 5 */ 342 FOR_EACH(i, 0, 6) Gs[i] += 16 * Gtmp[i] / (PI * PI); 343 344 alpha=lambda/(rho*cp); 345 a=box_size[0]; 346 b=box_size[1]; 347 c=box_size[2]; 348 349 /* Transient solution */ 350 FOR_EACH(m, 0, nterms_fs+1) { 351 const double beta = PI*(double)m/a; 352 const double beta_sqr = beta*beta; 353 const int m_is_even = (m%2 == 0); 354 355 FOR_EACH(n, 0, nterms_fs+1) { 356 const double gamma = PI*(double)n/b; 357 const double gamma_sqr = gamma * gamma; 358 const int n_is_even = (n%2==0); 359 360 FOR_EACH(o, 0, nterms_fs+1) { 361 const double eta = PI*(double)o/c; 362 const double eta_sqr = eta*eta; 363 const double zeta = alpha*(beta_sqr+gamma_sqr+eta_sqr); 364 const int o_is_even = (o%2==0); 365 double Fxyzt; 366 367 memset(Gtmp, 0, sizeof(double[7])); 368 FOR_EACH(p, 0, nt_pq+1) { 369 FOR_EACH(q, 0, nt_pq+1) { 370 const double p2 = (double)(2*p + 1); 371 const double q2 = (double)(2*q + 1); 372 const double Lx_sqr = PI*PI*((p2*p2)/(b*b)+(q2*q2)/(c*c)); 373 const double Ly_sqr = PI*PI*((p2*p2)/(a*a)+(q2*q2)/(c*c)); 374 const double Lz_sqr = PI*PI*((p2*p2)/(b*b)+(q2*q2)/(a*a)); 375 const double pq = p2*q2; 376 double itg[11] = {0}; 377 378 itg[1] = 2*q-o+1 == 0 ? -c/2.0 : 0; 379 itg[2] = eta / (eta_sqr+Lz_sqr); 380 itg[4] = 2*p-n+1 == 0 ? -b/2.0 : 0; 381 itg[5] = gamma / (gamma_sqr+Ly_sqr); 382 itg[7] = beta / (beta_sqr+Lx_sqr); 383 itg[9] = 2*p-m+1 == 0 ? -a/2.0 : 0; 384 itg[10] = 2*q-m+1 == 0 ? -a/2.0 : 0; 385 386 if((2*q-o+1==0) && (2*p-n+1==0)) { 387 const double z1 = itg[1]*itg[4]*itg[7]; 388 Gtmp[0] += z1/pq; 389 Gtmp[1] += (z1*pow(-1.0,(double)(m+1)))/pq; 390 } 391 if((2*q-o+1==0) && (2*p-m+1==0)) { 392 const double z2 = itg[1]*itg[9]*itg[5]; 393 Gtmp[2] += z2/pq; 394 Gtmp[3] += (z2*pow(-1.0,(double)(n+1)))/pq; 395 } 396 if((2*p-n+1==0) && (2*q-m+1==0)) { 397 const double z3 = itg[4]*itg[10]*itg[2]; 398 Gtmp[4] += z3/pq; 399 Gtmp[5] += (z3*pow(-1.0,(double)(o+1)))/pq; 400 } 401 } 402 } 403 Fxyzt = 404 sin(probe[0]*beta) 405 * sin(probe[1]*gamma) 406 * sin(probe[2]*eta) 407 * exp(-zeta*time); 408 409 FOR_EACH(i, 0, 6) { 410 Gi[i] += Gtmp[i] * Fxyzt; 411 } 412 if((!m_is_even) && (!n_is_even) && (!o_is_even)) { 413 Gi[6] += 8.0 * Fxyzt/(beta*gamma*eta); 414 } 415 } 416 } 417 } 418 419 /* Gi[i], i=0,5: Green of boundary index i */ 420 FOR_EACH(i, 0, 6) { 421 Gi[i] = -128.0 * Gi[i]/(a*b*c*PI*PI); 422 } 423 424 /* Gi[6]: Green of the initial condition */ 425 Gi[6] = 8.0*Gi[6]/(a*b*c); 426 427 /* Computing total Green function */ 428 FOR_EACH(i, 0, 7) { 429 green[i] = Gs[i] + Gi[i]; 430 } 431 } 432 433 static double 434 temperature_analytical 435 (const double temperature_bounds[6], 436 const double temperature_init, 437 const double box_size[3], 438 const double probe[3], 439 const double time, 440 const double rho, 441 const double cp, 442 const double lambda) 443 { 444 double green[7]; 445 double temperature = 0; 446 size_t i; 447 CHK(temperature_bounds && temperature_init && box_size && probe); 448 green_analytical(box_size, probe, time, rho, cp, lambda, green); 449 450 FOR_EACH(i, 0, 6) { 451 printf("Green for face %lu: %g\n", (unsigned long)i, green[i]); 452 temperature += green[i] * temperature_bounds[i]; 453 } 454 temperature += green[6] * temperature_init; 455 return temperature; 456 } 457 458 /******************************************************************************* 459 * Main function 460 ******************************************************************************/ 461 int 462 main(int argc, char** argv) 463 { 464 struct sdis_device* dev = NULL; 465 struct sdis_scene* box_scn = NULL; 466 struct sdis_scene* box2_scn = NULL; 467 struct sdis_scene* box_matriochka_scn = NULL; 468 struct sdis_medium* fluid = NULL; 469 struct sdis_medium* solid = NULL; 470 struct sdis_data* data = NULL; 471 struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT; 472 struct sdis_fluid_shader fluid_shader = SDIS_FLUID_SHADER_NULL; 473 struct sdis_solid_shader solid_shader = SDIS_SOLID_SHADER_NULL; 474 struct sdis_interface_shader interf_shader = SDIS_INTERFACE_SHADER_NULL; 475 struct sdis_interface* interfs[7] = {NULL}; 476 struct sdis_estimator* estimator = NULL; 477 struct sdis_mc temperature = SDIS_MC_NULL; 478 struct sdis_solve_probe_args solve_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT; 479 struct solid* solid_param = NULL; 480 struct context ctx; 481 struct matriochka_context matriochka_ctx; 482 const size_t nrealisations = 10000; 483 const size_t nmatriochkas = 5; 484 size_t nfails = 0; 485 double probe[3]; 486 double time[2]; 487 double Tbounds[6]; 488 double Tinit; 489 double Tref; 490 double boxsz[3]; 491 double rho, cp, lambda; 492 (void)argc, (void)argv; 493 494 /* System description */ 495 rho = 1700; 496 cp = 800; 497 lambda = 1.15; 498 Tbounds[0] = 280; /* Xmin */ 499 Tbounds[1] = 290; /* Xmax */ 500 Tbounds[2] = 310; /* Ymin */ 501 Tbounds[3] = 270; /* Ymax */ 502 Tbounds[4] = 300; /* Zmin */ 503 Tbounds[5] = 320; /* Zmax */ 504 Tinit = 100; 505 boxsz[0] = 0.3; 506 boxsz[1] = 0.1; 507 boxsz[2] = 0.2; 508 509 OK(sdis_device_create(&SDIS_DEVICE_CREATE_ARGS_DEFAULT, &dev)); 510 511 /* Create the fluid medium */ 512 fluid_shader = DUMMY_FLUID_SHADER; 513 OK(sdis_fluid_create(dev, &fluid_shader, NULL, &fluid)); 514 515 /* Setup the solid shader */ 516 solid_shader.calorific_capacity = solid_get_calorific_capacity; 517 solid_shader.thermal_conductivity = solid_get_thermal_conductivity; 518 solid_shader.volumic_mass = solid_get_volumic_mass; 519 solid_shader.delta = solid_get_delta; 520 solid_shader.temperature = solid_get_temperature; 521 522 /* Create the solid medium */ 523 OK(sdis_data_create 524 (dev, sizeof(struct solid), ALIGNOF(struct solid), NULL, &data)); 525 solid_param = sdis_data_get(data); 526 solid_param->rho = rho; 527 solid_param->cp = cp; 528 solid_param->lambda = lambda; 529 solid_param->delta = 1.0/30.0 * MMIN(MMIN(boxsz[0], boxsz[1]), boxsz[2]); 530 solid_param->init_temperature = Tinit; 531 OK(sdis_solid_create(dev, &solid_shader, data, &solid)); 532 533 /* Setup the interface shader */ 534 interf_shader.front.temperature = interface_get_temperature; 535 536 /* Create the interfaces */ 537 interfs[0] = create_interface(dev, solid, fluid, &interf_shader, Tbounds[0]); 538 interfs[1] = create_interface(dev, solid, fluid, &interf_shader, Tbounds[1]); 539 interfs[2] = create_interface(dev, solid, fluid, &interf_shader, Tbounds[2]); 540 interfs[3] = create_interface(dev, solid, fluid, &interf_shader, Tbounds[3]); 541 interfs[4] = create_interface(dev, solid, fluid, &interf_shader, Tbounds[4]); 542 interfs[5] = create_interface(dev, solid, fluid, &interf_shader, Tbounds[5]); 543 interfs[6] = create_interface(dev, solid, solid, &interf_shader, SDIS_TEMPERATURE_NONE); 544 545 /* Setup the box scene context */ 546 ctx.indices = box_indices; 547 ctx.vertices = box_vertices; 548 ctx.interfs[0] = ctx.interfs[1] = interfs[4]; /* Zmin */ 549 ctx.interfs[2] = ctx.interfs[3] = interfs[0]; /* Xmin */ 550 ctx.interfs[4] = ctx.interfs[5] = interfs[5]; /* Zmax */ 551 ctx.interfs[6] = ctx.interfs[7] = interfs[1]; /* Xmax */ 552 ctx.interfs[8] = ctx.interfs[9] = interfs[3]; /* Ymax */ 553 ctx.interfs[10] = ctx.interfs[11] = interfs[2]; /* Ymin */ 554 ctx.scale = boxsz; 555 556 /* Create the box scene */ 557 scn_args.get_indices = get_indices; 558 scn_args.get_interface = get_interface; 559 scn_args.get_position = get_position; 560 scn_args.nprimitives = box_ntriangles; 561 scn_args.nvertices = box_nvertices; 562 scn_args.context = &ctx; 563 OK(sdis_scene_create(dev, &scn_args, &box_scn)); 564 565 /* Setup the box2 scene context */ 566 ctx.indices = indices; 567 ctx.vertices = vertices; 568 ctx.interfs[0] = ctx.interfs[1] = interfs[0]; /* Xmin */ 569 ctx.interfs[2] = ctx.interfs[3] = interfs[6]; /* Xmid */ 570 ctx.interfs[4] = ctx.interfs[5] = interfs[1]; /* Xmax */ 571 ctx.interfs[6] = ctx.interfs[7] = interfs[2]; /* Ymin */ 572 ctx.interfs[8] = ctx.interfs[9] = interfs[2]; /* Ymin */ 573 ctx.interfs[10] = ctx.interfs[11] = interfs[3]; /* Ymax */ 574 ctx.interfs[12] = ctx.interfs[13] = interfs[3]; /* Ymax */ 575 ctx.interfs[14] = ctx.interfs[15] = interfs[4]; /* Zmin */ 576 ctx.interfs[16] = ctx.interfs[17] = interfs[4]; /* Zmin */ 577 ctx.interfs[18] = ctx.interfs[19] = interfs[5]; /* Zmax */ 578 ctx.interfs[20] = ctx.interfs[21] = interfs[5]; /* Zmax */ 579 ctx.scale = boxsz; 580 581 /* Create the box2 scene */ 582 scn_args.nprimitives = ntriangles; 583 scn_args.nvertices = nvertices; 584 OK(sdis_scene_create(dev, &scn_args, &box2_scn)); 585 586 /* Setup the matriochka context */ 587 matriochka_ctx.interfs[0] = matriochka_ctx.interfs[1] = interfs[4]; /* Zmin */ 588 matriochka_ctx.interfs[2] = matriochka_ctx.interfs[3] = interfs[0]; /* Xmin */ 589 matriochka_ctx.interfs[4] = matriochka_ctx.interfs[5] = interfs[5]; /* Zmax */ 590 matriochka_ctx.interfs[6] = matriochka_ctx.interfs[7] = interfs[1]; /* Xmax */ 591 matriochka_ctx.interfs[8] = matriochka_ctx.interfs[9] = interfs[3]; /* Ymax */ 592 matriochka_ctx.interfs[10] = matriochka_ctx.interfs[11] = interfs[2]; /* Ymin */ 593 matriochka_ctx.interfs[12] = interfs[6]; /* The remaining internal triangles */ 594 matriochka_ctx.scale = boxsz; 595 matriochka_ctx.nboxes = nmatriochkas; 596 597 /* Create the matriochka scene */ 598 scn_args.get_indices = matriochka_indices; 599 scn_args.get_interface = matriochka_interface; 600 scn_args.get_position = matriochka_position; 601 scn_args.nprimitives = box_ntriangles*nmatriochkas; 602 scn_args.nvertices = box_nvertices*nmatriochkas; 603 scn_args.context = &matriochka_ctx; 604 OK(sdis_scene_create(dev, &scn_args, &box_matriochka_scn)); 605 606 /* Setup and run the simulation */ 607 probe[0] = 0.1; 608 probe[1] = 0.06; 609 probe[2] = 0.130; 610 time[0] = time[1] = 500; /* Observation time range */ 611 612 /* Compute the solution */ 613 Tref = temperature_analytical 614 (Tbounds, Tinit, boxsz, probe, time[0], rho, cp, lambda); 615 616 /* Run simulation on regular scene */ 617 solid_param->delta = 1.0/30.0 * MMIN(MMIN(boxsz[0], boxsz[1]), boxsz[2]); 618 solve_args.nrealisations = nrealisations; 619 solve_args.position[0] = probe[0]; 620 solve_args.position[1] = probe[1]; 621 solve_args.position[2] = probe[2]; 622 solve_args.time_range[0] = time[0]; 623 solve_args.time_range[1] = time[1]; 624 OK(sdis_solve_probe(box_scn, &solve_args, &estimator)); 625 626 OK(sdis_estimator_get_failure_count(estimator, &nfails)); 627 OK(sdis_estimator_get_temperature(estimator, &temperature)); 628 printf("Temperature at (%g, %g, %g) m at %g s (delta = %g) = %g ~ %g +/- %g\n", 629 SPLIT3(probe), time[0], solid_param->delta, Tref, temperature.E, 630 temperature.SE); 631 printf("#failures = %lu/%lu\n", 632 (unsigned long)nfails, (unsigned long)nrealisations); 633 CHK(eq_eps(Tref, temperature.E, temperature.SE*3)); 634 OK(sdis_estimator_ref_put(estimator)); 635 636 /* Run simulation on split scene */ 637 solid_param->delta = 1.0/30.0 * MMIN(MMIN(boxsz[0]/2.0, boxsz[1]), boxsz[2]); 638 OK(sdis_solve_probe(box2_scn, &solve_args, &estimator)); 639 OK(sdis_estimator_get_failure_count(estimator, &nfails)); 640 OK(sdis_estimator_get_temperature(estimator, &temperature)); 641 printf("Temperature at (%g, %g, %g) m at %g s (delta = %g) = %g ~ %g +/- %g\n", 642 SPLIT3(probe), time[0], solid_param->delta, Tref, temperature.E, 643 temperature.SE); 644 printf("#failures = %lu/%lu\n", 645 (unsigned long)nfails, (unsigned long)nrealisations); 646 CHK(eq_eps(Tref, temperature.E, temperature.SE*3)); 647 OK(sdis_estimator_ref_put(estimator)); 648 649 /* Run simulation on matriochkas */ 650 solid_param->delta = MMIN(MMIN(boxsz[0], boxsz[1]), boxsz[2]); 651 solid_param->delta /= (double)nmatriochkas; 652 solid_param->delta *= 1.0/30.0; 653 OK(sdis_solve_probe(box_matriochka_scn, &solve_args, &estimator)); 654 OK(sdis_estimator_get_failure_count(estimator, &nfails)); 655 OK(sdis_estimator_get_temperature(estimator, &temperature)); 656 printf("Temperature at (%g, %g, %g) m at %g s (delta = %g) = %g ~ %g +/- %g\n", 657 SPLIT3(probe), time[0], solid_param->delta, Tref, temperature.E, 658 temperature.SE); 659 printf("#failures = %lu/%lu\n", 660 (unsigned long)nfails, (unsigned long)nrealisations); 661 CHK(eq_eps(Tref, temperature.E, temperature.SE*3)); 662 OK(sdis_estimator_ref_put(estimator)); 663 664 OK(sdis_interface_ref_put(interfs[0])); 665 OK(sdis_interface_ref_put(interfs[1])); 666 OK(sdis_interface_ref_put(interfs[2])); 667 OK(sdis_interface_ref_put(interfs[3])); 668 OK(sdis_interface_ref_put(interfs[4])); 669 OK(sdis_interface_ref_put(interfs[5])); 670 OK(sdis_interface_ref_put(interfs[6])); 671 OK(sdis_data_ref_put(data)); 672 OK(sdis_medium_ref_put(solid)); 673 OK(sdis_medium_ref_put(fluid)); 674 OK(sdis_device_ref_put(dev)); 675 OK(sdis_scene_ref_put(box_scn)); 676 OK(sdis_scene_ref_put(box2_scn)); 677 OK(sdis_scene_ref_put(box_matriochka_scn)); 678 679 CHK(mem_allocated_size() == 0); 680 return 0; 681 }