test_sdis_contact_resistance_2.c (17074B)
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 #include "test_sdis_contact_resistance.h" 19 20 #include <rsys/clock_time.h> 21 #include <rsys/mem_allocator.h> 22 #include <rsys/double3.h> 23 #include <rsys/math.h> 24 #include <star/ssp.h> 25 26 /* 27 * The scene is composed of a solid cube/square of size L whose temperature is 28 * unknown. The cube is made of 2 solids that meet at x=e in ]0 L[ and the 29 * thermal contact resistance is R, thus T(X0-) differs from T(X0+). 30 * The faces are adiabatic exept at x=0 where T(0)=T0 and at x=L where T(L)=TL. 31 * At steady state: 32 * 33 * Flux(x0) = (T(x0+) - T(x0-)) / R 34 * 35 * 3D 2D 36 * 37 * /////////(L,L,L) /////////(L,L) 38 * +-------+ +-------+ 39 * /' / /| | ! | 40 * +-------+ TL T0 r TL 41 * | | ! | | | ! | 42 * T0 +.r...|.+ +-------+ 43 * |, ! |/ (0,0)///x=X0/// 44 * +-------+ 45 * (0,0,0)///x=X0/// 46 */ 47 48 #define N 10000 /* #realisations */ 49 50 #define T0 0.0 51 #define LAMBDA1 0.1 52 53 #define TL 100.0 54 #define LAMBDA2 0.2 55 56 #define DELTA1 X0/30.0 57 #define DELTA2 (L-X0)/30.0 58 59 /******************************************************************************* 60 * Media 61 ******************************************************************************/ 62 struct solid { 63 double lambda; 64 double rho; 65 double cp; 66 double delta; 67 }; 68 69 static double 70 fluid_get_temperature 71 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 72 { 73 (void)data; 74 CHK(vtx); 75 return SDIS_TEMPERATURE_NONE; 76 } 77 78 static double 79 solid_get_calorific_capacity 80 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 81 { 82 CHK(vtx && data); 83 return ((struct solid*)sdis_data_cget(data))->cp; 84 } 85 86 static double 87 solid_get_thermal_conductivity 88 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 89 { 90 CHK(vtx && data); 91 return ((struct solid*)sdis_data_cget(data))->lambda; 92 } 93 94 static double 95 solid_get_volumic_mass 96 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 97 { 98 CHK(vtx && data); 99 return ((struct solid*)sdis_data_cget(data))->rho; 100 } 101 102 static double 103 solid_get_delta 104 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 105 { 106 CHK(vtx && data); 107 return ((struct solid*)sdis_data_cget(data))->delta; 108 } 109 110 static double 111 solid_get_temperature 112 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 113 { 114 CHK(vtx && data); 115 return SDIS_TEMPERATURE_NONE; 116 } 117 118 /******************************************************************************* 119 * Interfaces 120 ******************************************************************************/ 121 struct interf { 122 double temperature; 123 double resistance; 124 }; 125 126 static double 127 interface_get_temperature 128 (const struct sdis_interface_fragment* frag, struct sdis_data* data) 129 { 130 const struct interf* interf = sdis_data_cget(data); 131 CHK(frag && data); 132 return interf->temperature; 133 } 134 135 static double 136 interface_get_convection_coef 137 (const struct sdis_interface_fragment* frag, struct sdis_data* data) 138 { 139 CHK(frag && data); 140 return 0; 141 } 142 143 static double 144 interface_get_contact_resistance 145 (const struct sdis_interface_fragment* frag, struct sdis_data* data) 146 { 147 const struct interf* interf = sdis_data_cget(data); 148 CHK(frag && data); 149 return interf->resistance; 150 } 151 152 /******************************************************************************* 153 * Helper functions 154 ******************************************************************************/ 155 static void 156 solve_probe 157 (struct sdis_scene* scn, 158 struct interf* interf_props, 159 struct ssp_rng* rng) 160 { 161 char dump[128]; 162 struct time t0, t1; 163 struct sdis_estimator* estimator; 164 struct sdis_solve_probe_boundary_args solve_args 165 = SDIS_SOLVE_PROBE_BOUNDARY_ARGS_DEFAULT; 166 struct sdis_mc T = SDIS_MC_NULL; 167 struct sdis_mc time = SDIS_MC_NULL; 168 size_t nreals; 169 size_t nfails; 170 double ref_L, ref_R; 171 enum sdis_scene_dimension dim; 172 int nsimuls; 173 int isimul; 174 ASSERT(scn && interf_props && rng); 175 176 OK(sdis_scene_get_dimension(scn, &dim)); 177 178 nsimuls = (dim == SDIS_SCENE_2D) ? 2 : 4; 179 FOR_EACH(isimul, 0, nsimuls) { 180 double ref; 181 double r = pow(10, ssp_rng_uniform_double(rng, -2, 2)); 182 183 interf_props->resistance = r; 184 185 ref_L = ( 186 T0 * (r * LAMBDA1 / X0) * (1 + r * LAMBDA2 / (L - X0)) 187 + TL * (r * LAMBDA2 / (L - X0)) 188 ) 189 / ((1 + r * LAMBDA1 / X0) * (1 + r * LAMBDA2 / (L - X0)) - 1); 190 191 ref_R = ref_L * (1 + r * LAMBDA1 / X0) - T0 * r * LAMBDA1 / X0; 192 193 if(dim == SDIS_SCENE_2D) 194 /* last segment */ 195 solve_args.iprim = model2d_nsegments - 1; 196 else 197 /* last 2 triangles */ 198 solve_args.iprim = model3d_ntriangles - ((isimul % 2) ? 2 : 1); 199 200 solve_args.uv[0] = ssp_rng_canonical(rng); 201 solve_args.uv[1] = (dim == SDIS_SCENE_2D) 202 ? 0 : ssp_rng_uniform_double(rng, 0, 1 - solve_args.uv[0]); 203 204 if(isimul < nsimuls / 2) { 205 solve_args.side = SDIS_FRONT; 206 ref = ref_L; 207 } else { 208 solve_args.side = SDIS_BACK; 209 ref = ref_R; 210 } 211 212 solve_args.nrealisations = N; 213 solve_args.time_range[0] = solve_args.time_range[1] = INF; 214 215 time_current(&t0); 216 OK(sdis_solve_probe_boundary(scn, &solve_args, &estimator)); 217 time_sub(&t0, time_current(&t1), &t0); 218 time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); 219 220 OK(sdis_estimator_get_realisation_count(estimator, &nreals)); 221 OK(sdis_estimator_get_failure_count(estimator, &nfails)); 222 OK(sdis_estimator_get_temperature(estimator, &T)); 223 OK(sdis_estimator_get_realisation_time(estimator, &time)); 224 225 switch(dim) { 226 case SDIS_SCENE_2D: 227 printf("Steady temperature at (%lu/%s/%g) with R=%g = %g ~ %g +/- %g\n", 228 (unsigned long)solve_args.iprim, 229 (solve_args.side == SDIS_FRONT ? "front" : "back"), 230 solve_args.uv[0], 231 r, ref, T.E, T.SE); 232 break; 233 case SDIS_SCENE_3D: 234 printf("Steady temperature at (%lu/%s/%g,%g) with R=%g = %g ~ %g +/- %g\n", 235 (unsigned long)solve_args.iprim, 236 (solve_args.side == SDIS_FRONT ? "front" : "back"), 237 SPLIT2(solve_args.uv), 238 r, ref, T.E, T.SE); 239 break; 240 default: FATAL("Unreachable code.\n"); break; 241 } 242 printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); 243 printf("Elapsed time = %s\n", dump); 244 printf("Time per realisation (in usec) = %g +/- %g\n\n", time.E, time.SE); 245 246 CHK(nfails + nreals == N); 247 CHK(nfails <= N/1000); 248 CHK(eq_eps(T.E, ref, T.SE * 3)); 249 250 OK(sdis_estimator_ref_put(estimator)); 251 } 252 } 253 static void 254 solve 255 (struct sdis_scene* scn, 256 struct interf* interf_props, 257 struct ssp_rng* rng) 258 { 259 char dump[128]; 260 struct time t0, t1; 261 struct sdis_estimator* estimator; 262 struct sdis_solve_boundary_args solve_args = SDIS_SOLVE_BOUNDARY_ARGS_DEFAULT; 263 struct sdis_mc T = SDIS_MC_NULL; 264 struct sdis_mc time = SDIS_MC_NULL; 265 const enum sdis_side all_front[] = { SDIS_FRONT, SDIS_FRONT }; 266 const enum sdis_side all_back[] = { SDIS_BACK, SDIS_BACK }; 267 size_t plist[2]; 268 size_t nreals; 269 size_t nfails; 270 double ref_L, ref_R; 271 enum sdis_scene_dimension dim; 272 int nsimuls; 273 int isimul; 274 ASSERT(scn && interf_props && rng); 275 276 OK(sdis_scene_get_dimension(scn, &dim)); 277 278 nsimuls = (dim == SDIS_SCENE_2D) ? 2 : 4; 279 FOR_EACH(isimul, 0, nsimuls) { 280 double ref; 281 double r = pow(10, ssp_rng_uniform_double(rng, -2, 2)); 282 283 interf_props->resistance = r; 284 285 ref_L = ( 286 T0 * (r * LAMBDA1 / X0) * (1 + r * LAMBDA2 / (L - X0)) 287 + TL * (r * LAMBDA2 / (L - X0)) 288 ) 289 / ((1 + r * LAMBDA1 / X0) * (1 + r * LAMBDA2 / (L - X0)) - 1); 290 291 ref_R = ref_L * (1 + r * LAMBDA1 / X0) - T0 * r * LAMBDA1 / X0; 292 293 if(dim == SDIS_SCENE_2D) { 294 /* last segment */ 295 solve_args.nprimitives = 1; 296 plist[0] = model2d_nsegments - 1; 297 } else { 298 /* last 2 triangles */ 299 solve_args.nprimitives = 2; 300 plist[0] = model3d_ntriangles - 2; 301 plist[1] = model3d_ntriangles - 1; 302 } 303 solve_args.primitives = plist; 304 305 if(isimul < nsimuls / 2) { 306 solve_args.sides = all_front; 307 ref = ref_L; 308 } else { 309 solve_args.sides = all_back; 310 ref = ref_R; 311 } 312 313 solve_args.nrealisations = N; 314 solve_args.time_range[0] = solve_args.time_range[1] = INF; 315 316 time_current(&t0); 317 OK(sdis_solve_boundary(scn, &solve_args, &estimator)); 318 time_sub(&t0, time_current(&t1), &t0); 319 time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); 320 321 OK(sdis_estimator_get_realisation_count(estimator, &nreals)); 322 OK(sdis_estimator_get_failure_count(estimator, &nfails)); 323 OK(sdis_estimator_get_temperature(estimator, &T)); 324 OK(sdis_estimator_get_realisation_time(estimator, &time)); 325 326 printf("Steady temperature at the %s side with R=%g = %g ~ %g +/- %g\n", 327 (solve_args.sides[0] == SDIS_FRONT ? "front" : "back"), 328 r, ref, T.E, T.SE); 329 printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); 330 printf("Elapsed time = %s\n", dump); 331 printf("Time per realisation (in usec) = %g +/- %g\n\n", time.E, time.SE); 332 333 CHK(nfails + nreals == N); 334 CHK(nfails <= N / 1000); 335 CHK(eq_eps(T.E, ref, T.SE * 3)); 336 337 OK(sdis_estimator_ref_put(estimator)); 338 } 339 } 340 341 /******************************************************************************* 342 * Test 343 ******************************************************************************/ 344 int 345 main(int argc, char** argv) 346 { 347 struct sdis_data* data = NULL; 348 struct sdis_device* dev = NULL; 349 struct sdis_medium* fluid = NULL; 350 struct sdis_medium* solid1 = NULL; 351 struct sdis_medium* solid2 = NULL; 352 struct sdis_interface* interf_adiabatic1 = NULL; 353 struct sdis_interface* interf_adiabatic2 = NULL; 354 struct sdis_interface* interf_T0 = NULL; 355 struct sdis_interface* interf_TL = NULL; 356 struct sdis_interface* interf_R = NULL; 357 struct sdis_scene* box_scn = NULL; 358 struct sdis_scene* square_scn = NULL; 359 struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT; 360 struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER; 361 struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER; 362 struct sdis_interface_shader interf_shader = SDIS_INTERFACE_SHADER_NULL; 363 struct sdis_interface* model3d_interfaces[22 /*#triangles*/]; 364 struct sdis_interface* model2d_interfaces[7/*#segments*/]; 365 struct interf* interf_props = NULL; 366 struct solid* solid_props = NULL; 367 struct ssp_rng* rng = NULL; 368 (void)argc, (void)argv; 369 370 OK(sdis_device_create(&SDIS_DEVICE_CREATE_ARGS_DEFAULT, &dev)); 371 372 fluid_shader.temperature = fluid_get_temperature; 373 OK(sdis_fluid_create(dev, &fluid_shader, NULL, &fluid)); 374 375 /* Setup the solid shader */ 376 solid_shader.calorific_capacity = solid_get_calorific_capacity; 377 solid_shader.thermal_conductivity = solid_get_thermal_conductivity; 378 solid_shader.volumic_mass = solid_get_volumic_mass; 379 solid_shader.delta = solid_get_delta; 380 solid_shader.temperature = solid_get_temperature; 381 382 /* Create the solid medium #1 */ 383 OK(sdis_data_create(dev, sizeof(struct solid), 16, NULL, &data)); 384 solid_props = sdis_data_get(data); 385 solid_props->lambda = LAMBDA1; 386 solid_props->cp = 2; 387 solid_props->rho = 25; 388 solid_props->delta = DELTA1; 389 OK(sdis_solid_create(dev, &solid_shader, data, &solid1)); 390 OK(sdis_data_ref_put(data)); 391 392 /* Create the solid medium #2 */ 393 OK(sdis_data_create(dev, sizeof(struct solid), 16, NULL, &data)); 394 solid_props = sdis_data_get(data); 395 solid_props->lambda = LAMBDA2; 396 solid_props->cp = 2; 397 solid_props->rho = 25; 398 solid_props->delta = DELTA2; 399 OK(sdis_solid_create(dev, &solid_shader, data, &solid2)); 400 OK(sdis_data_ref_put(data)); 401 402 /* Setup the interface shader */ 403 interf_shader.front.temperature = interface_get_temperature; 404 interf_shader.back.temperature = interface_get_temperature; 405 interf_shader.convection_coef = interface_get_convection_coef; 406 407 /* Create the adiabatic interfaces */ 408 OK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data)); 409 interf_props = sdis_data_get(data); 410 interf_props->temperature = SDIS_TEMPERATURE_NONE; 411 OK(sdis_interface_create 412 (dev, solid1, fluid, &interf_shader, data, &interf_adiabatic1)); 413 OK(sdis_interface_create 414 (dev, solid2, fluid, &interf_shader, data, &interf_adiabatic2)); 415 OK(sdis_data_ref_put(data)); 416 417 /* Create the T0 interface */ 418 OK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data)); 419 interf_props = sdis_data_get(data); 420 interf_props->temperature = T0; 421 OK(sdis_interface_create 422 (dev, solid1, fluid, &interf_shader, data, &interf_T0)); 423 OK(sdis_data_ref_put(data)); 424 425 /* Create the TL interface */ 426 OK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data)); 427 interf_props = sdis_data_get(data); 428 interf_props->temperature = TL; 429 OK(sdis_interface_create 430 (dev, solid2, fluid, &interf_shader, data, &interf_TL)); 431 OK(sdis_data_ref_put(data)); 432 433 /* Create the solid1-solid2 interface */ 434 interf_shader.convection_coef = NULL; 435 interf_shader.thermal_contact_resistance = interface_get_contact_resistance; 436 OK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data)); 437 interf_props = sdis_data_get(data); 438 interf_props->temperature = SDIS_TEMPERATURE_NONE; 439 OK(sdis_interface_create 440 (dev, solid1, solid2, &interf_shader, data, &interf_R)); 441 OK(sdis_data_ref_put(data)); 442 443 /* Release the media */ 444 OK(sdis_medium_ref_put(solid1)); 445 OK(sdis_medium_ref_put(solid2)); 446 OK(sdis_medium_ref_put(fluid)); 447 448 /* Map the interfaces to their box triangles */ 449 /* Front */ 450 model3d_interfaces[0] = interf_adiabatic1; 451 model3d_interfaces[1] = interf_adiabatic1; 452 model3d_interfaces[2] = interf_adiabatic2; 453 model3d_interfaces[3] = interf_adiabatic2; 454 /* Left */ 455 model3d_interfaces[4] = interf_T0; 456 model3d_interfaces[5] = interf_T0; 457 /* Back */ 458 model3d_interfaces[6] = interf_adiabatic1; 459 model3d_interfaces[7] = interf_adiabatic1; 460 model3d_interfaces[8] = interf_adiabatic2; 461 model3d_interfaces[9] = interf_adiabatic2; 462 /* Right */ 463 model3d_interfaces[10] = interf_TL; 464 model3d_interfaces[11] = interf_TL; 465 /* Top */ 466 model3d_interfaces[12] = interf_adiabatic1; 467 model3d_interfaces[13] = interf_adiabatic1; 468 model3d_interfaces[14] = interf_adiabatic2; 469 model3d_interfaces[15] = interf_adiabatic2; 470 /* Bottom */ 471 model3d_interfaces[16] = interf_adiabatic1; 472 model3d_interfaces[17] = interf_adiabatic1; 473 model3d_interfaces[18] = interf_adiabatic2; 474 model3d_interfaces[19] = interf_adiabatic2; 475 /* Inside */ 476 model3d_interfaces[20] = interf_R; 477 model3d_interfaces[21] = interf_R; 478 479 /* Map the interfaces to their square segments */ 480 /* Bottom */ 481 model2d_interfaces[0] = interf_adiabatic2; 482 model2d_interfaces[1] = interf_adiabatic1; 483 /* Left */ 484 model2d_interfaces[2] = interf_T0; 485 /* Top */ 486 model2d_interfaces[3] = interf_adiabatic1; 487 model2d_interfaces[4] = interf_adiabatic2; 488 /* Right */ 489 model2d_interfaces[5] = interf_TL; 490 /* Contact */ 491 model2d_interfaces[6] = interf_R; 492 493 /* Create the box scene */ 494 scn_args.get_indices = model3d_get_indices; 495 scn_args.get_interface = model3d_get_interface; 496 scn_args.get_position = model3d_get_position; 497 scn_args.nprimitives = model3d_ntriangles; 498 scn_args.nvertices = model3d_nvertices; 499 scn_args.context = model3d_interfaces; 500 OK(sdis_scene_create(dev, &scn_args, &box_scn)); 501 502 /* Create the square scene */ 503 scn_args.get_indices = model2d_get_indices; 504 scn_args.get_interface = model2d_get_interface; 505 scn_args.get_position = model2d_get_position; 506 scn_args.nprimitives = model2d_nsegments; 507 scn_args.nvertices = model2d_nvertices; 508 scn_args.context = model2d_interfaces; 509 OK(sdis_scene_2d_create(dev, &scn_args, &square_scn)); 510 511 /* Release the interfaces */ 512 OK(sdis_interface_ref_put(interf_adiabatic1)); 513 OK(sdis_interface_ref_put(interf_adiabatic2)); 514 OK(sdis_interface_ref_put(interf_T0)); 515 OK(sdis_interface_ref_put(interf_TL)); 516 OK(sdis_interface_ref_put(interf_R)); 517 518 /* Solve */ 519 OK(ssp_rng_create(NULL, SSP_RNG_KISS, &rng)); 520 printf(">> Box scene\n"); 521 solve_probe(box_scn, interf_props, rng); 522 solve(box_scn, interf_props, rng); 523 printf("\n>> Square scene\n"); 524 solve_probe(square_scn, interf_props, rng); 525 solve(square_scn, interf_props, rng); 526 527 OK(sdis_scene_ref_put(box_scn)); 528 OK(sdis_scene_ref_put(square_scn)); 529 OK(sdis_device_ref_put(dev)); 530 OK(ssp_rng_ref_put(rng)); 531 532 CHK(mem_allocated_size() == 0); 533 return 0; 534 } 535