test_sdis_contact_resistance.c (14425B)
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 * T(X0-) = (T0 * (R * LAMBDA1 / X0) * (1 + R * LAMBDA2 / (L - X0)) 34 * + TL * (R * LAMBDA2 / (L - X0))) 35 * / ((1 + R * LAMBDA1 / X0) * (1 + R * LAMBDA2 / (L - X0)) - 1) 36 * 37 * T(X0+) = T(X0-) * (1 + r * LAMBDA1 / X0) - T0 * r * LAMBDA1 / X0 38 * 39 * T(x) is linear between T(0) and T(X0-) if x in [0 X0[ 40 * T(x) is linear between T(X0+) and T(L) if x in ]X0 L] 41 * 42 * 3D 2D 43 * 44 * /////////(L,L,L) /////////(L,L) 45 * +-------+ +-------+ 46 * /' / /| | ! | 47 * +-------+ TL T0 r TL 48 * | | ! | | | ! | 49 * T0 +.r...|.+ +-------+ 50 * |, ! |/ (0,0)///x=X0/// 51 * +-------+ 52 * (0,0,0)///x=X0/// 53 */ 54 55 #define N 10000 /* #realisations */ 56 57 #define T0 0.0 58 #define LAMBDA1 0.1 59 60 #define TL 100.0 61 #define LAMBDA2 0.2 62 63 #define DELTA1 X0/25.0 64 #define DELTA2 (L-X0)/25.0 65 66 /******************************************************************************* 67 * Media 68 ******************************************************************************/ 69 struct solid { 70 double lambda; 71 double rho; 72 double cp; 73 double delta; 74 }; 75 76 static double 77 fluid_get_temperature 78 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 79 { 80 (void)data; 81 CHK(vtx); 82 return SDIS_TEMPERATURE_NONE; 83 } 84 85 static double 86 solid_get_calorific_capacity 87 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 88 { 89 CHK(vtx && data); 90 return ((struct solid*)sdis_data_cget(data))->cp; 91 } 92 93 static double 94 solid_get_thermal_conductivity 95 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 96 { 97 CHK(vtx && data); 98 return ((struct solid*)sdis_data_cget(data))->lambda; 99 } 100 101 static double 102 solid_get_volumic_mass 103 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 104 { 105 CHK(vtx && data); 106 return ((struct solid*)sdis_data_cget(data))->rho; 107 } 108 109 static double 110 solid_get_delta 111 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 112 { 113 CHK(vtx && data); 114 return ((struct solid*)sdis_data_cget(data))->delta; 115 } 116 117 static double 118 solid_get_temperature 119 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 120 { 121 CHK(vtx && data); 122 return SDIS_TEMPERATURE_NONE; 123 } 124 125 /******************************************************************************* 126 * Interfaces 127 ******************************************************************************/ 128 struct interf { 129 double temperature; 130 double resistance; 131 }; 132 133 static double 134 interface_get_temperature 135 (const struct sdis_interface_fragment* frag, struct sdis_data* data) 136 { 137 const struct interf* interf = sdis_data_cget(data); 138 CHK(frag && data); 139 return interf->temperature; 140 } 141 142 static double 143 interface_get_convection_coef 144 (const struct sdis_interface_fragment* frag, struct sdis_data* data) 145 { 146 CHK(frag && data); 147 return 0; 148 } 149 150 static double 151 interface_get_contact_resistance 152 (const struct sdis_interface_fragment* frag, struct sdis_data* data) 153 { 154 const struct interf* interf = sdis_data_cget(data); 155 CHK(frag && data); 156 return interf->resistance; 157 } 158 159 /******************************************************************************* 160 * Helper functions 161 ******************************************************************************/ 162 static void 163 solve 164 (struct sdis_scene* scn, 165 struct interf* interf_props, 166 struct ssp_rng* rng) 167 { 168 char dump[128]; 169 struct time t0, t1; 170 struct sdis_estimator* estimator; 171 struct sdis_solve_probe_args solve_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT; 172 struct sdis_mc T = SDIS_MC_NULL; 173 struct sdis_mc time = SDIS_MC_NULL; 174 size_t nreals; 175 size_t nfails; 176 double ref_L, ref_R; 177 enum sdis_scene_dimension dim; 178 const int nsimuls = 8; 179 int isimul; 180 ASSERT(scn && interf_props && rng); 181 182 OK(sdis_scene_get_dimension(scn, &dim)); 183 184 FOR_EACH(isimul, 0, nsimuls) { 185 double x, ref; 186 double r = pow(10, ssp_rng_uniform_double(rng, -2, 2)); 187 188 interf_props->resistance = r; 189 190 ref_L = ( 191 T0 * (r * LAMBDA1 / X0) * (1 + r * LAMBDA2 / (L - X0)) 192 + TL * (r * LAMBDA2 / (L - X0)) 193 ) 194 / ((1 + r * LAMBDA1 / X0) * (1 + r * LAMBDA2 / (L - X0)) - 1); 195 196 ref_R = ref_L * (1 + r * LAMBDA1 / X0) - T0 * r * LAMBDA1 / X0; 197 198 if(isimul % 2) { /* In solid 1 */ 199 x = ssp_rng_uniform_double(rng, 0.05 * X0, 0.95 * X0); 200 ref = T0 * (1 - x / X0) + ref_L * x / X0; 201 } else { /* In solid 2 */ 202 x = X0 + ssp_rng_uniform_double(rng, 0.05 * (L - X0), 0.95 * (L - X0)); 203 ref = ref_R * (1 - (x - X0) / (L - X0)) + TL * (x - X0) / (L - X0); 204 } 205 206 solve_args.position[0] = x; 207 solve_args.position[1] = ssp_rng_uniform_double(rng, 0.05 * L, 0.95 * L); 208 solve_args.position[2] = (dim == SDIS_SCENE_2D) 209 ? 0 : ssp_rng_uniform_double(rng, 0.05 * L, 0.95 * L); 210 211 solve_args.nrealisations = N; 212 solve_args.time_range[0] = solve_args.time_range[1] = INF; 213 214 time_current(&t0); 215 OK(sdis_solve_probe(scn, &solve_args, &estimator)); 216 time_sub(&t0, time_current(&t1), &t0); 217 time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); 218 219 OK(sdis_estimator_get_realisation_count(estimator, &nreals)); 220 OK(sdis_estimator_get_failure_count(estimator, &nfails)); 221 OK(sdis_estimator_get_temperature(estimator, &T)); 222 OK(sdis_estimator_get_realisation_time(estimator, &time)); 223 224 switch(dim) { 225 case SDIS_SCENE_2D: 226 printf("Steady temperature at (%g, %g) with R=%g = %g ~ %g +/- %g\n", 227 SPLIT2(solve_args.position), r, ref, T.E, T.SE); 228 break; 229 case SDIS_SCENE_3D: 230 printf("Steady temperature at (%g, %g, %g) with R=%g = %g ~ %g +/- %g\n", 231 SPLIT3(solve_args.position), r, ref, T.E, T.SE); 232 break; 233 default: FATAL("Unreachable code.\n"); break; 234 } 235 printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); 236 printf("Elapsed time = %s\n", dump); 237 printf("Time per realisation (in usec) = %g +/- %g\n\n", time.E, time.SE); 238 239 CHK(nfails + nreals == N); 240 CHK(nfails <= N/1000); 241 CHK(eq_eps(T.E, ref, T.SE * 3)); 242 243 OK(sdis_estimator_ref_put(estimator)); 244 } 245 } 246 247 /******************************************************************************* 248 * Test 249 ******************************************************************************/ 250 int 251 main(int argc, char** argv) 252 { 253 struct sdis_data* data = NULL; 254 struct sdis_device* dev = NULL; 255 struct sdis_medium* fluid = NULL; 256 struct sdis_medium* solid1 = NULL; 257 struct sdis_medium* solid2 = NULL; 258 struct sdis_interface* interf_adiabatic1 = NULL; 259 struct sdis_interface* interf_adiabatic2 = NULL; 260 struct sdis_interface* interf_T0 = NULL; 261 struct sdis_interface* interf_TL = NULL; 262 struct sdis_interface* interf_R = NULL; 263 struct sdis_scene* box_scn = NULL; 264 struct sdis_scene* square_scn = NULL; 265 struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT; 266 struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER; 267 struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER; 268 struct sdis_interface_shader interf_shader = SDIS_INTERFACE_SHADER_NULL; 269 struct sdis_interface* model3d_interfaces[22 /*#triangles*/]; 270 struct sdis_interface* model2d_interfaces[7/*#segments*/]; 271 struct interf* interf_props = NULL; 272 struct solid* solid_props = NULL; 273 struct ssp_rng* rng = NULL; 274 (void)argc, (void)argv; 275 276 OK(sdis_device_create(&SDIS_DEVICE_CREATE_ARGS_DEFAULT, &dev)); 277 278 fluid_shader.temperature = fluid_get_temperature; 279 OK(sdis_fluid_create(dev, &fluid_shader, NULL, &fluid)); 280 281 /* Setup the solid shader */ 282 solid_shader.calorific_capacity = solid_get_calorific_capacity; 283 solid_shader.thermal_conductivity = solid_get_thermal_conductivity; 284 solid_shader.volumic_mass = solid_get_volumic_mass; 285 solid_shader.delta = solid_get_delta; 286 solid_shader.temperature = solid_get_temperature; 287 288 /* Create the solid medium #1 */ 289 OK(sdis_data_create(dev, sizeof(struct solid), 16, NULL, &data)); 290 solid_props = sdis_data_get(data); 291 solid_props->lambda = LAMBDA1; 292 solid_props->cp = 2; 293 solid_props->rho = 25; 294 solid_props->delta = DELTA1; 295 OK(sdis_solid_create(dev, &solid_shader, data, &solid1)); 296 OK(sdis_data_ref_put(data)); 297 298 /* Create the solid medium #2 */ 299 OK(sdis_data_create(dev, sizeof(struct solid), 16, NULL, &data)); 300 solid_props = sdis_data_get(data); 301 solid_props->lambda = LAMBDA2; 302 solid_props->cp = 2; 303 solid_props->rho = 25; 304 solid_props->delta = DELTA2; 305 OK(sdis_solid_create(dev, &solid_shader, data, &solid2)); 306 OK(sdis_data_ref_put(data)); 307 308 /* Setup the interface shader */ 309 interf_shader.front.temperature = interface_get_temperature; 310 interf_shader.back.temperature = interface_get_temperature; 311 interf_shader.convection_coef = interface_get_convection_coef; 312 313 /* Create the adiabatic interfaces */ 314 OK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data)); 315 interf_props = sdis_data_get(data); 316 interf_props->temperature = SDIS_TEMPERATURE_NONE; 317 OK(sdis_interface_create 318 (dev, solid1, fluid, &interf_shader, data, &interf_adiabatic1)); 319 OK(sdis_interface_create 320 (dev, solid2, fluid, &interf_shader, data, &interf_adiabatic2)); 321 OK(sdis_data_ref_put(data)); 322 323 /* Create the T0 interface */ 324 OK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data)); 325 interf_props = sdis_data_get(data); 326 interf_props->temperature = T0; 327 OK(sdis_interface_create 328 (dev, solid1, fluid, &interf_shader, data, &interf_T0)); 329 OK(sdis_data_ref_put(data)); 330 331 /* Create the TL interface */ 332 OK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data)); 333 interf_props = sdis_data_get(data); 334 interf_props->temperature = TL; 335 OK(sdis_interface_create 336 (dev, solid2, fluid, &interf_shader, data, &interf_TL)); 337 OK(sdis_data_ref_put(data)); 338 339 /* Create the solid1-solid2 interface */ 340 interf_shader.convection_coef = NULL; 341 interf_shader.thermal_contact_resistance = interface_get_contact_resistance; 342 OK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data)); 343 interf_props = sdis_data_get(data); 344 interf_props->temperature = SDIS_TEMPERATURE_NONE; 345 OK(sdis_interface_create 346 (dev, solid1, solid2, &interf_shader, data, &interf_R)); 347 OK(sdis_data_ref_put(data)); 348 349 /* Release the media */ 350 OK(sdis_medium_ref_put(solid1)); 351 OK(sdis_medium_ref_put(solid2)); 352 OK(sdis_medium_ref_put(fluid)); 353 354 /* Map the interfaces to their box triangles */ 355 /* Front */ 356 model3d_interfaces[0] = interf_adiabatic1; 357 model3d_interfaces[1] = interf_adiabatic1; 358 model3d_interfaces[2] = interf_adiabatic2; 359 model3d_interfaces[3] = interf_adiabatic2; 360 /* Left */ 361 model3d_interfaces[4] = interf_T0; 362 model3d_interfaces[5] = interf_T0; 363 /* Back */ 364 model3d_interfaces[6] = interf_adiabatic1; 365 model3d_interfaces[7] = interf_adiabatic1; 366 model3d_interfaces[8] = interf_adiabatic2; 367 model3d_interfaces[9] = interf_adiabatic2; 368 /* Right */ 369 model3d_interfaces[10] = interf_TL; 370 model3d_interfaces[11] = interf_TL; 371 /* Top */ 372 model3d_interfaces[12] = interf_adiabatic1; 373 model3d_interfaces[13] = interf_adiabatic1; 374 model3d_interfaces[14] = interf_adiabatic2; 375 model3d_interfaces[15] = interf_adiabatic2; 376 /* Bottom */ 377 model3d_interfaces[16] = interf_adiabatic1; 378 model3d_interfaces[17] = interf_adiabatic1; 379 model3d_interfaces[18] = interf_adiabatic2; 380 model3d_interfaces[19] = interf_adiabatic2; 381 /* Inside */ 382 model3d_interfaces[20] = interf_R; 383 model3d_interfaces[21] = interf_R; 384 385 /* Map the interfaces to their square segments */ 386 /* Bottom */ 387 model2d_interfaces[0] = interf_adiabatic2; 388 model2d_interfaces[1] = interf_adiabatic1; 389 /* Left */ 390 model2d_interfaces[2] = interf_T0; 391 /* Top */ 392 model2d_interfaces[3] = interf_adiabatic1; 393 model2d_interfaces[4] = interf_adiabatic2; 394 /* Right */ 395 model2d_interfaces[5] = interf_TL; 396 /* Contact */ 397 model2d_interfaces[6] = interf_R; 398 399 /* Create the box scene */ 400 scn_args.get_indices = model3d_get_indices; 401 scn_args.get_interface = model3d_get_interface; 402 scn_args.get_position = model3d_get_position; 403 scn_args.nprimitives = model3d_ntriangles; 404 scn_args.nvertices = model3d_nvertices; 405 scn_args.context = model3d_interfaces; 406 OK(sdis_scene_create(dev, &scn_args, &box_scn)); 407 408 /* Create the square scene */ 409 scn_args.get_indices = model2d_get_indices; 410 scn_args.get_interface = model2d_get_interface; 411 scn_args.get_position = model2d_get_position; 412 scn_args.nprimitives = model2d_nsegments; 413 scn_args.nvertices = model2d_nvertices; 414 scn_args.context = model2d_interfaces; 415 OK(sdis_scene_2d_create(dev, &scn_args, &square_scn)); 416 417 /* Release the interfaces */ 418 OK(sdis_interface_ref_put(interf_adiabatic1)); 419 OK(sdis_interface_ref_put(interf_adiabatic2)); 420 OK(sdis_interface_ref_put(interf_T0)); 421 OK(sdis_interface_ref_put(interf_TL)); 422 OK(sdis_interface_ref_put(interf_R)); 423 424 /* Solve */ 425 OK(ssp_rng_create(NULL, SSP_RNG_KISS, &rng)); 426 printf(">> Box scene\n"); 427 solve(box_scn, interf_props, rng); 428 printf("\n>> Square scene\n"); 429 solve(square_scn, interf_props, rng); 430 431 OK(sdis_scene_ref_put(box_scn)); 432 OK(sdis_scene_ref_put(square_scn)); 433 OK(sdis_device_ref_put(dev)); 434 OK(ssp_rng_ref_put(rng)); 435 436 CHK(mem_allocated_size() == 0); 437 return 0; 438 } 439