test_sdis_volumic_power.c (17132B)
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 19 #include <rsys/clock_time.h> 20 #include <rsys/double3.h> 21 #include <rsys/math.h> 22 #include <star/ssp.h> 23 24 /* 25 * The scene is composed of a solid cube/square whose temperature is unknown. 26 * The convection coefficient with the surrounding fluid is null everywhere The 27 * Temperature of the +/- X faces are fixed to T0, and the solid has a volumic 28 * power of P0. This test computes the temperature of a probe position pos into 29 * the solid and check that at t=inf it is is equal to: 30 * 31 * T(pos) = P0 / (2*LAMBDA) * (A^2/4 - (pos-0.5)^2) + T0 32 * 33 * with LAMBDA the conductivity of the solid and A the size of the cube/square, 34 * i.e. 1. 35 * 36 * 3D 2D 37 * 38 * ///// (1,1,1) ///// (1,1) 39 * +-------+ +-------+ 40 * /' /| | | 41 * +-------+ T0 T0 T0 42 * T0 +.....|.+ | | 43 * |, |/ +-------+ 44 * +-------+ (0,0) ///// 45 * (0,0,0) ///// 46 */ 47 48 #define N 10000 /* #realisations */ 49 50 #define T0 320 51 #define LAMBDA 0.1 52 #define P0 10 53 #define DELTA 1.0/60.0 54 55 /******************************************************************************* 56 * Helper functions 57 ******************************************************************************/ 58 static const char* 59 algo_cstr(const enum sdis_diffusion_algorithm diff_algo) 60 { 61 const char* cstr = "none"; 62 63 switch(diff_algo) { 64 case SDIS_DIFFUSION_DELTA_SPHERE: cstr = "delta sphere"; break; 65 case SDIS_DIFFUSION_WOS: cstr = "WoS"; break; 66 default: FATAL("Unreachable code.\n"); break; 67 } 68 return cstr; 69 } 70 71 /******************************************************************************* 72 * Media 73 ******************************************************************************/ 74 struct solid { 75 double lambda; 76 double rho; 77 double cp; 78 double delta; 79 double vpower; 80 double initial_temperature; 81 double t0; 82 }; 83 84 static double 85 fluid_get_temperature 86 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 87 { 88 (void)data; 89 CHK(vtx != NULL); 90 return SDIS_TEMPERATURE_NONE; 91 } 92 93 static double 94 solid_get_calorific_capacity 95 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 96 { 97 CHK(vtx != NULL); 98 return ((struct solid*)sdis_data_cget(data))->cp; 99 } 100 101 static double 102 solid_get_thermal_conductivity 103 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 104 { 105 CHK(vtx != NULL); 106 return ((struct solid*)sdis_data_cget(data))->lambda; 107 } 108 109 static double 110 solid_get_volumic_mass 111 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 112 { 113 CHK(vtx != NULL); 114 return ((struct solid*)sdis_data_cget(data))->rho; 115 } 116 117 static double 118 solid_get_delta 119 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 120 { 121 CHK(vtx != NULL); 122 return ((struct solid*)sdis_data_cget(data))->delta; 123 } 124 125 static double 126 solid_get_temperature 127 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 128 { 129 double t0; 130 CHK(vtx != NULL); 131 CHK(data != NULL); 132 t0 = ((const struct solid*)sdis_data_cget(data))->t0; 133 if(vtx->time > t0) { 134 return SDIS_TEMPERATURE_NONE; 135 } else { 136 return ((const struct solid*)sdis_data_cget(data))->initial_temperature; 137 } 138 } 139 140 static double 141 solid_get_volumic_power 142 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 143 { 144 (void)data; 145 CHK(vtx != NULL); 146 return ((struct solid*)sdis_data_cget(data))->vpower; 147 } 148 149 /******************************************************************************* 150 * Interfaces 151 ******************************************************************************/ 152 struct interf { 153 double temperature; 154 }; 155 156 static double 157 interface_get_temperature 158 (const struct sdis_interface_fragment* frag, struct sdis_data* data) 159 { 160 const struct interf* interf = sdis_data_cget(data); 161 CHK(frag && data); 162 return interf->temperature; 163 } 164 165 static double 166 interface_get_convection_coef 167 (const struct sdis_interface_fragment* frag, struct sdis_data* data) 168 { 169 CHK(frag && data); 170 return 0; 171 } 172 173 /******************************************************************************* 174 * Helper functions 175 ******************************************************************************/ 176 static void 177 solve 178 (struct sdis_scene* scn, 179 struct ssp_rng* rng, 180 struct solid* solid) 181 { 182 char dump[128]; 183 struct time t0, t1, t2; 184 struct sdis_estimator* estimator; 185 struct sdis_estimator* estimator2; 186 struct sdis_green_function* green; 187 struct sdis_solve_probe_args solve_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT; 188 struct sdis_mc T = SDIS_MC_NULL; 189 struct sdis_mc time = SDIS_MC_NULL; 190 size_t nreals; 191 size_t nfails; 192 enum sdis_scene_dimension dim; 193 const int nsimuls = 4; 194 int isimul; 195 ASSERT(scn && rng && solid); 196 197 OK(sdis_scene_get_dimension(scn, &dim)); 198 FOR_EACH(isimul, 0, nsimuls) { 199 const enum sdis_diffusion_algorithm algo = (isimul / 2) % 2 200 ? SDIS_DIFFUSION_WOS 201 : SDIS_DIFFUSION_DELTA_SPHERE; 202 const int steady = (isimul % 2) == 0; 203 double power = P0 == SDIS_VOLUMIC_POWER_NONE ? 0 : P0; 204 205 /* Restore power value */ 206 solid->vpower = P0; 207 208 solve_args.diff_algo = algo; 209 solve_args.position[0] = ssp_rng_uniform_double(rng, 0.1, 0.9); 210 solve_args.position[1] = ssp_rng_uniform_double(rng, 0.1, 0.9); 211 solve_args.position[2] = 212 dim == SDIS_SCENE_2D ? 0 : ssp_rng_uniform_double(rng, 0.1, 0.9); 213 214 solve_args.nrealisations = N; 215 if(steady) { 216 solve_args.time_range[0] = solve_args.time_range[1] = INF; 217 } else { 218 solve_args.time_range[0] = 100 * (double)isimul; 219 solve_args.time_range[1] = 4 * solve_args.time_range[0]; 220 } 221 222 time_current(&t0); 223 OK(sdis_solve_probe(scn, &solve_args, &estimator)); 224 time_sub(&t0, time_current(&t1), &t0); 225 time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); 226 227 OK(sdis_estimator_get_realisation_count(estimator, &nreals)); 228 OK(sdis_estimator_get_failure_count(estimator, &nfails)); 229 OK(sdis_estimator_get_temperature(estimator, &T)); 230 OK(sdis_estimator_get_realisation_time(estimator, &time)); 231 232 if(steady) { 233 const double x = solve_args.position[0] - 0.5; 234 const double ref = power / (2 * LAMBDA) * (1.0 / 4.0 - x * x) + T0; 235 printf 236 ("Steady temperature - pos: %g, %g, %g; Power: %g algo: %s " 237 "= %g ~ %g +/- %g\n", 238 SPLIT3(solve_args.position), power, algo_cstr(algo), ref, T.E, T.SE); 239 CHK(eq_eps(T.E, ref, T.SE * 3)); 240 } else { 241 printf( 242 "Mean temperature - pos: %g, %g, %g; t in [%g %g]; power: %g; algo: %s " 243 "~ %g +/- %g\n", 244 SPLIT3(solve_args.position), SPLIT2(solve_args.time_range), 245 power, algo_cstr(algo), T.E, T.SE); 246 } 247 248 printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); 249 printf("Elapsed time = %s\n", dump); 250 printf("Time per realisation (in usec) = %g +/- %g\n\n", time.E, time.SE); 251 252 CHK(nfails + nreals == N); 253 CHK(nfails <= N/1000); 254 255 /* Check green function */ 256 time_current(&t0); 257 OK(sdis_solve_probe_green_function(scn, &solve_args, &green)); 258 time_current(&t1); 259 OK(sdis_green_function_solve(green, &estimator2)); 260 time_current(&t2); 261 262 OK(sdis_estimator_get_realisation_count(estimator2, &nreals)); 263 OK(sdis_estimator_get_failure_count(estimator2, &nfails)); 264 OK(sdis_estimator_get_temperature(estimator2, &T)); 265 266 if(steady) { 267 const double x = solve_args.position[0] - 0.5; 268 const double ref = power / (2 * LAMBDA) * (1.0 / 4.0 - x * x) + T0; 269 printf 270 ("Green steady temperature - pos: %g, %g, %g; Power: %g algo: %s" 271 "= %g ~ %g +/- %g\n", 272 SPLIT3(solve_args.position), power, algo_cstr(algo), ref, T.E, T.SE); 273 } else { 274 printf( 275 "Green mean temperature - pos: %g, %g, %g; t: [%g %g]; power: %g; algo: %s " 276 "~ %g +/- %g\n", 277 SPLIT3(solve_args.position), SPLIT2(solve_args.time_range), 278 power, algo_cstr(algo), T.E, T.SE); 279 } 280 281 printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); 282 time_sub(&t0, &t1, &t0); 283 time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); 284 printf("Green estimation time = %s\n", dump); 285 time_sub(&t1, &t2, &t1); 286 time_dump(&t1, TIME_ALL, NULL, dump, sizeof(dump)); 287 printf("Green solve time = %s\n", dump); 288 289 check_green_function(green); 290 check_estimator_eq(estimator, estimator2); 291 check_green_serialization(green, scn); 292 293 OK(sdis_estimator_ref_put(estimator)); 294 OK(sdis_estimator_ref_put(estimator2)); 295 printf("\n"); 296 297 /* Check same green used at a different power level */ 298 solid->vpower = power = 3 * P0; 299 300 time_current(&t0); 301 OK(sdis_solve_probe(scn, &solve_args, &estimator)); 302 time_sub(&t0, time_current(&t1), &t0); 303 time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); 304 305 OK(sdis_estimator_get_realisation_count(estimator, &nreals)); 306 OK(sdis_estimator_get_failure_count(estimator, &nfails)); 307 OK(sdis_estimator_get_temperature(estimator, &T)); 308 OK(sdis_estimator_get_realisation_time(estimator, &time)); 309 310 if(steady) { 311 const double x = solve_args.position[0] - 0.5; 312 const double ref = power / (2 * LAMBDA) * (1.0 / 4.0 - x * x) + T0; 313 printf 314 ("Steady temperature - pos: %g, %g, %g; Power: %g algo: %s " 315 "= %g ~ %g +/- %g\n", 316 SPLIT3(solve_args.position), power, algo_cstr(algo), ref, T.E, T.SE); 317 CHK(eq_eps(T.E, ref, T.SE * 3)); 318 } else { 319 printf( 320 "Mean temperature - pos: %g, %g, %g; t in [%g %g]; power: %g; algo: %s " 321 "~ %g +/- %g\n", 322 SPLIT3(solve_args.position), SPLIT2(solve_args.time_range), 323 power, algo_cstr(algo), T.E, T.SE); 324 } 325 326 printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); 327 printf("Elapsed time = %s\n", dump); 328 printf("Time per realisation (in usec) = %g +/- %g\n", time.E, time.SE); 329 330 CHK(nfails + nreals == N); 331 CHK(nfails <= N/1000); 332 333 time_current(&t0); 334 OK(sdis_green_function_solve(green, &estimator2)); 335 time_sub(&t0, time_current(&t1), &t0); 336 time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); 337 printf("Green solve time = %s\n", dump); 338 339 check_green_function(green); 340 check_estimator_eq(estimator, estimator2); 341 342 OK(sdis_estimator_ref_put(estimator)); 343 OK(sdis_estimator_ref_put(estimator2)); 344 OK(sdis_green_function_ref_put(green)); 345 346 printf("\n\n"); 347 } 348 } 349 350 static void 351 check_null_power_term_with_green(struct sdis_scene* scn, struct solid* solid) 352 { 353 struct sdis_solve_probe_args solve_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT; 354 struct sdis_mc T = SDIS_MC_NULL; 355 struct sdis_estimator* estimator = NULL; 356 struct sdis_green_function* green = NULL; 357 double x = 0; 358 double ref = 0; 359 ASSERT(scn && solid); 360 361 solve_args.position[0] = 0.5; 362 solve_args.position[1] = 0.5; 363 solve_args.position[2] = 0; 364 solve_args.nrealisations = N; 365 solve_args.time_range[0] = 366 solve_args.time_range[1] = INF; 367 368 solid->vpower = 0; 369 OK(sdis_solve_probe_green_function(scn, &solve_args, &green)); 370 371 solid->vpower = P0; 372 OK(sdis_green_function_solve(green, &estimator)); 373 OK(sdis_estimator_get_temperature(estimator, &T)); 374 375 x = solve_args.position[0] - 0.5; 376 ref = solid->vpower / (2*LAMBDA) * (1.0/4.0 - x *x) + T0; 377 printf("Green steady temperature at (%g, %g, %g) with Power=%g = %g ~ %g +/- %g\n", 378 SPLIT3(solve_args.position), solid->vpower, ref, T.E, T.SE); 379 CHK(eq_eps(ref, T.E, 3*T.SE)); 380 381 OK(sdis_estimator_ref_put(estimator)); 382 OK(sdis_green_function_ref_put(green)); 383 } 384 385 /******************************************************************************* 386 * Test 387 ******************************************************************************/ 388 int 389 main(int argc, char** argv) 390 { 391 struct sdis_data* data = NULL; 392 struct sdis_device* dev = NULL; 393 struct sdis_medium* fluid = NULL; 394 struct sdis_medium* solid = NULL; 395 struct sdis_interface* interf_adiabatic = NULL; 396 struct sdis_interface* interf_T0 = NULL; 397 struct sdis_scene* box_scn = NULL; 398 struct sdis_scene* square_scn = NULL; 399 struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT; 400 struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER; 401 struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER; 402 struct sdis_interface_shader interf_shader = SDIS_INTERFACE_SHADER_NULL; 403 struct sdis_interface* box_interfaces[12 /*#triangles*/]; 404 struct sdis_interface* square_interfaces[4/*#segments*/]; 405 struct interf* interf_props = NULL; 406 struct solid* solid_props = NULL; 407 struct ssp_rng* rng = NULL; 408 (void)argc, (void)argv; 409 410 OK(sdis_device_create(&SDIS_DEVICE_CREATE_ARGS_DEFAULT, &dev)); 411 412 fluid_shader.temperature = fluid_get_temperature; 413 OK(sdis_fluid_create(dev, &fluid_shader, NULL, &fluid)); 414 415 /* Setup the solid shader */ 416 solid_shader.calorific_capacity = solid_get_calorific_capacity; 417 solid_shader.thermal_conductivity = solid_get_thermal_conductivity; 418 solid_shader.volumic_mass = solid_get_volumic_mass; 419 solid_shader.delta = solid_get_delta; 420 solid_shader.temperature = solid_get_temperature; 421 solid_shader.volumic_power = solid_get_volumic_power; 422 423 /* Create the solid medium */ 424 OK(sdis_data_create(dev, sizeof(struct solid), 16, NULL, &data)); 425 solid_props = sdis_data_get(data); 426 solid_props->lambda = LAMBDA; 427 solid_props->cp = 2; 428 solid_props->rho = 25; 429 solid_props->delta = DELTA; 430 solid_props->vpower = P0; 431 solid_props->t0 = 0; 432 solid_props->initial_temperature = T0; 433 OK(sdis_solid_create(dev, &solid_shader, data, &solid)); 434 OK(sdis_data_ref_put(data)); 435 436 /* Setup the interface shader */ 437 interf_shader.convection_coef = interface_get_convection_coef; 438 interf_shader.front.temperature = interface_get_temperature; 439 440 /* Create the adiabatic interface */ 441 OK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data)); 442 interf_props = sdis_data_get(data); 443 interf_props->temperature = SDIS_TEMPERATURE_NONE; 444 OK(sdis_interface_create 445 (dev, solid, fluid, &interf_shader, data, &interf_adiabatic)); 446 OK(sdis_data_ref_put(data)); 447 448 /* Create the T0 interface */ 449 OK(sdis_data_create(dev, sizeof(struct interf), 16, NULL, &data)); 450 interf_props = sdis_data_get(data); 451 interf_props->temperature = T0; 452 OK(sdis_interface_create 453 (dev, solid, fluid, &interf_shader, data, &interf_T0)); 454 OK(sdis_data_ref_put(data)); 455 456 /* Release the media */ 457 OK(sdis_medium_ref_put(solid)); 458 OK(sdis_medium_ref_put(fluid)); 459 460 /* Map the interfaces to their box triangles */ 461 box_interfaces[0] = box_interfaces[1] = interf_adiabatic; /* Front */ 462 box_interfaces[2] = box_interfaces[3] = interf_T0; /* Left */ 463 box_interfaces[4] = box_interfaces[5] = interf_adiabatic; /* Back */ 464 box_interfaces[6] = box_interfaces[7] = interf_T0; /* Right */ 465 box_interfaces[8] = box_interfaces[9] = interf_adiabatic; /* Top */ 466 box_interfaces[10]= box_interfaces[11]= interf_adiabatic; /* Bottom */ 467 468 /* Map the interfaces to their square segments */ 469 square_interfaces[0] = interf_adiabatic; /* Bottom */ 470 square_interfaces[1] = interf_T0; /* Left */ 471 square_interfaces[2] = interf_adiabatic; /* Top */ 472 square_interfaces[3] = interf_T0; /* Right */ 473 474 /* Create the box scene */ 475 scn_args.get_indices = box_get_indices; 476 scn_args.get_interface = box_get_interface; 477 scn_args.get_position = box_get_position; 478 scn_args.nprimitives = box_ntriangles; 479 scn_args.nvertices = box_nvertices; 480 scn_args.context = box_interfaces; 481 OK(sdis_scene_create(dev, &scn_args, &box_scn)); 482 483 /* Create the square scene */ 484 scn_args.get_indices = square_get_indices; 485 scn_args.get_interface = square_get_interface; 486 scn_args.get_position = square_get_position; 487 scn_args.nprimitives = square_nsegments; 488 scn_args.nvertices = square_nvertices; 489 scn_args.context = square_interfaces; 490 OK(sdis_scene_2d_create(dev, &scn_args, &square_scn)); 491 492 /* Release the interfaces */ 493 OK(sdis_interface_ref_put(interf_adiabatic)); 494 OK(sdis_interface_ref_put(interf_T0)); 495 496 /* Solve */ 497 OK(ssp_rng_create(NULL, SSP_RNG_MT19937_64, &rng)); 498 printf(">> Box scene\n"); 499 solve(box_scn, rng, solid_props); 500 printf(">> Square scene\n"); 501 solve(square_scn, rng, solid_props); 502 503 /* Check green registration with a null power term */ 504 check_null_power_term_with_green(box_scn, solid_props); 505 506 OK(sdis_scene_ref_put(box_scn)); 507 OK(sdis_scene_ref_put(square_scn)); 508 OK(sdis_device_ref_put(dev)); 509 OK(ssp_rng_ref_put(rng)); 510 511 CHK(mem_allocated_size() == 0); 512 return 0; 513 } 514