test_sdis_solid_random_walk_robustness.c (15153B)
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 <star/s3dut.h> 20 #include <rsys/math.h> 21 22 #define Tfluid 350.0 /* Temperature of the fluid in Kelvin */ 23 #define LAMBDA 10.0 /* Thermal conductivity */ 24 #define Hcoef 1.0 /* Convection coefficient */ 25 #define Pw 10000.0 /* Volumetric power */ 26 #define Nreals 10000 /* #realisations */ 27 28 /******************************************************************************* 29 * Helper functions 30 ******************************************************************************/ 31 static double 32 trilinear_temperature(const double pos[3]) 33 { 34 const double a = 333; 35 const double b = 432; 36 const double c = 579; 37 double x, y, z; 38 CHK(pos[0] >= -10.0 && pos[0] <= 10.0); 39 CHK(pos[1] >= -10.0 && pos[1] <= 10.0); 40 CHK(pos[2] >= -10.0 && pos[2] <= 10.0); 41 42 x = (pos[0] + 10.0) / 20.0; 43 y = (pos[1] + 10.0) / 20.0; 44 z = (pos[2] + 10.0) / 20.0; 45 46 return a*x + b*y + c*z; 47 } 48 49 static double 50 volumetric_temperature(const double pos[3], const double upper[3]) 51 { 52 const double beta = - 1.0 / 3.0 * Pw / (2*LAMBDA); 53 const double temp = 54 beta * (pos[0]*pos[0] - upper[0]*upper[0]) 55 + beta * (pos[1]*pos[1] - upper[1]*upper[1]) 56 + beta * (pos[2]*pos[2] - upper[2]*upper[2]); 57 CHK(temp > 0); 58 return temp; 59 } 60 61 static const char* 62 algo_cstr(const enum sdis_diffusion_algorithm diff_algo) 63 { 64 const char* cstr = "none"; 65 66 switch(diff_algo) { 67 case SDIS_DIFFUSION_DELTA_SPHERE: cstr = "delta sphere"; break; 68 case SDIS_DIFFUSION_WOS: cstr = "WoS"; break; 69 default: FATAL("Unreachable code.\n"); break; 70 } 71 return cstr; 72 } 73 74 /******************************************************************************* 75 * Geometry 76 ******************************************************************************/ 77 struct context { 78 struct s3dut_mesh_data msh; 79 struct sdis_interface* interf; 80 }; 81 82 static void 83 get_indices(const size_t itri, size_t ids[3], void* context) 84 { 85 const struct context* ctx = context; 86 CHK(ids && ctx && itri < ctx->msh.nprimitives); 87 ids[0] = ctx->msh.indices[itri*3+0]; 88 ids[2] = ctx->msh.indices[itri*3+1]; 89 ids[1] = ctx->msh.indices[itri*3+2]; 90 } 91 92 static void 93 get_position(const size_t ivert, double pos[3], void* context) 94 { 95 const struct context* ctx = context; 96 CHK(pos && ctx && ivert < ctx->msh.nvertices); 97 pos[0] = ctx->msh.positions[ivert*3+0]; 98 pos[1] = ctx->msh.positions[ivert*3+1]; 99 pos[2] = ctx->msh.positions[ivert*3+2]; 100 } 101 102 static void 103 get_interface(const size_t itri, struct sdis_interface** bound, void* context) 104 { 105 const struct context* ctx = context; 106 CHK(bound && ctx && itri < ctx->msh.nprimitives); 107 *bound = ctx->interf; 108 } 109 110 static struct sdis_scene* 111 create_scene(struct sdis_device* dev, struct sdis_interface* interf) 112 { 113 /* Star-3DUT */ 114 struct s3dut_super_formula f0 = S3DUT_SUPER_FORMULA_NULL; 115 struct s3dut_super_formula f1 = S3DUT_SUPER_FORMULA_NULL; 116 struct s3dut_mesh* msh = NULL; 117 118 /* Stardis */ 119 struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT; 120 struct sdis_scene* scn = NULL; 121 122 /* Miscellaneous */ 123 struct context ctx; 124 125 ASSERT(dev && interf); 126 127 /* Create the solid super shape */ 128 f0.A = 1; f0.B = 1; f0.M = 20; f0.N0 = 1; f0.N1 = 1; f0.N2 = 5; 129 f1.A = 1; f1.B = 1; f1.M = 7; f1.N0 = 1; f1.N1 = 2; f1.N2 = 5; 130 OK(s3dut_create_super_shape(NULL, &f0, &f1, 1, 128, 64, &msh)); 131 OK(s3dut_mesh_get_data(msh, &ctx.msh)); 132 133 /*dump_mesh(stdout, ctx.msh.positions, 134 ctx.msh.nvertices, ctx.msh.indices, ctx.msh.nprimitives);*/ 135 136 /* Create the scene */ 137 ctx.interf = interf; 138 scn_args.get_indices = get_indices; 139 scn_args.get_interface = get_interface; 140 scn_args.get_position = get_position; 141 scn_args.nprimitives = ctx.msh.nprimitives; 142 scn_args.nvertices = ctx.msh.nvertices; 143 scn_args.context = &ctx; 144 OK(sdis_scene_create(dev, &scn_args, &scn)); 145 146 OK(s3dut_mesh_ref_put(msh)); 147 148 return scn; 149 } 150 151 /******************************************************************************* 152 * Interface 153 ******************************************************************************/ 154 enum profile { 155 PROFILE_UNKNOWN, 156 PROFILE_TRILINEAR, 157 PROFILE_VOLUMETRIC_POWER, 158 PROFILE_COUNT__ 159 }; 160 struct interf { 161 enum profile profile; 162 double upper[3]; /* Upper bound of the scene */ 163 double h; 164 }; 165 166 static double 167 interface_get_temperature 168 (const struct sdis_interface_fragment* frag, struct sdis_data* data) 169 { 170 const struct interf* interf; 171 double temperature; 172 CHK(data != NULL && frag != NULL); 173 interf = sdis_data_cget(data); 174 switch(interf->profile) { 175 case PROFILE_UNKNOWN: 176 temperature = SDIS_TEMPERATURE_NONE; 177 break; 178 case PROFILE_VOLUMETRIC_POWER: 179 temperature = volumetric_temperature(frag->P, interf->upper); 180 break; 181 case PROFILE_TRILINEAR: 182 temperature = trilinear_temperature(frag->P); 183 break; 184 default: FATAL("Unreachable code.\n"); break; 185 } 186 return temperature; 187 } 188 189 static double 190 interface_get_convection_coef 191 (const struct sdis_interface_fragment* frag, struct sdis_data* data) 192 { 193 CHK(frag != NULL); 194 return ((const struct interf*)sdis_data_cget(data))->h;; 195 } 196 197 static struct sdis_interface* 198 create_interface 199 (struct sdis_device* dev, 200 struct sdis_medium* front, 201 struct sdis_medium* back) 202 { 203 struct sdis_interface_shader interf_shader = SDIS_INTERFACE_SHADER_NULL; 204 struct sdis_interface* interf = NULL; 205 struct sdis_data* data = NULL; 206 struct interf* interf_param = NULL; 207 ASSERT(dev && front && back); 208 209 OK(sdis_data_create 210 (dev, sizeof(struct interf), ALIGNOF(struct interf), NULL, &data)); 211 interf_param = sdis_data_get(data); 212 interf_param->h = 1; 213 214 interf_shader.convection_coef = interface_get_convection_coef; 215 interf_shader.front.temperature = interface_get_temperature; 216 OK(sdis_interface_create(dev, front, back, &interf_shader, data, &interf)); 217 OK(sdis_data_ref_put(data)); 218 219 return interf; 220 } 221 222 /******************************************************************************* 223 * Fluid 224 ******************************************************************************/ 225 static double 226 fluid_get_temperature 227 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 228 { 229 CHK(vtx != NULL); 230 (void)data; 231 return Tfluid; 232 } 233 static struct sdis_medium* 234 create_fluid(struct sdis_device* dev) 235 { 236 struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER; 237 struct sdis_medium* fluid = NULL; 238 ASSERT(dev); 239 240 /* Create the fluid medium */ 241 fluid_shader.temperature = fluid_get_temperature; 242 OK(sdis_fluid_create(dev, &fluid_shader, NULL, &fluid)); 243 244 return fluid; 245 } 246 247 /******************************************************************************* 248 * Solid API 249 ******************************************************************************/ 250 struct solid { 251 double delta; 252 double cp; 253 double lambda; 254 double rho; 255 double temperature; 256 double power; 257 }; 258 259 static double 260 solid_get_calorific_capacity 261 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 262 { 263 CHK(data != NULL && vtx != NULL); 264 return ((const struct solid*)sdis_data_cget(data))->cp; 265 } 266 267 static double 268 solid_get_thermal_conductivity 269 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 270 { 271 CHK(data != NULL && vtx != NULL); 272 return ((const struct solid*)sdis_data_cget(data))->lambda; 273 } 274 275 static double 276 solid_get_volumic_mass 277 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 278 { 279 CHK(data != NULL && vtx != NULL); 280 return ((const struct solid*)sdis_data_cget(data))->rho; 281 } 282 283 static double 284 solid_get_delta 285 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 286 { 287 CHK(data != NULL && vtx != NULL); 288 return ((const struct solid*)sdis_data_cget(data))->delta; 289 } 290 291 static double 292 solid_get_temperature 293 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 294 { 295 CHK(data != NULL && vtx != NULL); 296 return ((const struct solid*)sdis_data_cget(data))->temperature; 297 } 298 299 static double 300 solid_get_volumetric_power 301 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 302 { 303 CHK(data != NULL && vtx != NULL); 304 return ((const struct solid*)sdis_data_cget(data))->power; 305 } 306 307 static struct sdis_medium* 308 create_solid(struct sdis_device* dev) 309 { 310 struct sdis_solid_shader solid_shader = SDIS_SOLID_SHADER_NULL; 311 struct sdis_medium* solid = NULL; 312 struct sdis_data* data = NULL; 313 struct solid* solid_param; 314 ASSERT(dev); 315 316 OK(sdis_data_create 317 (dev, sizeof(struct solid), ALIGNOF(struct solid), NULL, &data)); 318 solid_param = sdis_data_get(data); 319 solid_param->delta = 0.01; 320 solid_param->lambda = 10; 321 solid_param->cp = 1.0; 322 solid_param->rho = 1.0; 323 solid_param->temperature = SDIS_TEMPERATURE_NONE; 324 solid_param->power = SDIS_VOLUMIC_POWER_NONE; 325 326 solid_shader.calorific_capacity = solid_get_calorific_capacity; 327 solid_shader.thermal_conductivity = solid_get_thermal_conductivity; 328 solid_shader.volumic_mass = solid_get_volumic_mass; 329 solid_shader.delta = solid_get_delta; 330 solid_shader.temperature = solid_get_temperature; 331 solid_shader.volumic_power = solid_get_volumetric_power; 332 OK(sdis_solid_create(dev, &solid_shader, data, &solid)); 333 OK(sdis_data_ref_put(data)); 334 335 return solid; 336 } 337 338 /******************************************************************************* 339 * Solve functions 340 ******************************************************************************/ 341 static void 342 check_estimation 343 (const struct sdis_estimator* estimator, const double Tref) 344 { 345 struct sdis_mc T = SDIS_MC_NULL; 346 size_t nreals; 347 size_t nfails; 348 349 OK(sdis_estimator_get_temperature(estimator, &T)); 350 OK(sdis_estimator_get_realisation_count(estimator, &nreals)); 351 OK(sdis_estimator_get_failure_count(estimator, &nfails)); 352 CHK(nfails <= Nreals * 0.0005); 353 printf("T = %g ~ %g +/- %g [%g, %g]; #failures = %lu / %lu\n", 354 Tref, T.E, T.SE, T.E - 3*T.SE, T.E + 3*T.SE, 355 (unsigned long)nfails, (unsigned long)Nreals); 356 CHK(eq_eps(T.E, Tref, T.SE*3)); 357 } 358 359 static void 360 check_probe 361 (struct sdis_scene* scn, 362 const enum profile profile, 363 const enum sdis_diffusion_algorithm diff_algo) 364 { 365 struct sdis_solve_probe_args solve_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT; 366 struct sdis_estimator* estimator = NULL; 367 double lower[3], upper[3]; 368 double Tref; 369 ASSERT(scn); 370 371 printf("algo: %s\n", algo_cstr(diff_algo)); 372 373 solve_args.nrealisations = Nreals; 374 solve_args.position[0] = 0; 375 solve_args.position[1] = 0; 376 solve_args.position[2] = 0; 377 solve_args.time_range[0] = INF; 378 solve_args.time_range[1] = INF; 379 solve_args.diff_algo = diff_algo; 380 solve_args.register_paths = SDIS_HEAT_PATH_FAILURE; 381 OK(sdis_solve_probe(scn, &solve_args, &estimator)); 382 383 switch(profile) { 384 case PROFILE_TRILINEAR: 385 Tref = trilinear_temperature(solve_args.position); 386 break; 387 case PROFILE_VOLUMETRIC_POWER: 388 OK(sdis_scene_get_aabb(scn, lower, upper)); 389 Tref = volumetric_temperature(solve_args.position, upper); 390 break; 391 default: FATAL("Unreachable code.\n"); break; 392 } 393 394 check_estimation(estimator, Tref); 395 OK(sdis_estimator_ref_put(estimator)); 396 } 397 398 static void 399 check_medium 400 (struct sdis_scene* scn, 401 struct sdis_medium* medium, 402 const enum sdis_diffusion_algorithm diff_algo) 403 { 404 struct sdis_solve_medium_args solve_mdm_args = SDIS_SOLVE_MEDIUM_ARGS_DEFAULT; 405 struct sdis_estimator* estimator = NULL; 406 ASSERT(scn); 407 408 printf("algo: %s\n", algo_cstr(diff_algo)); 409 410 solve_mdm_args.nrealisations = Nreals; 411 solve_mdm_args.medium = medium; 412 solve_mdm_args.time_range[0] = INF; 413 solve_mdm_args.time_range[1] = INF; 414 solve_mdm_args.diff_algo = diff_algo; 415 OK(sdis_solve_medium(scn, &solve_mdm_args, &estimator)); 416 417 check_estimation(estimator, Tfluid); 418 /*dump_heat_paths(stdout, estimator);*/ 419 OK(sdis_estimator_ref_put(estimator)); 420 } 421 422 /******************************************************************************* 423 * Main function 424 ******************************************************************************/ 425 int 426 main(int argc, char** argv) 427 { 428 struct sdis_device* dev = NULL; 429 struct sdis_medium* fluid = NULL; 430 struct sdis_medium* solid = NULL; 431 struct sdis_interface* interf = NULL; 432 struct sdis_scene* scn = NULL; 433 struct interf* interf_param = NULL; 434 struct solid* solid_param = NULL; 435 double lower[3]; 436 double upper[3]; 437 double spread; 438 (void)argc, (void)argv; 439 440 OK(sdis_device_create(&SDIS_DEVICE_CREATE_ARGS_DEFAULT, &dev)); 441 442 fluid = create_fluid(dev); 443 solid = create_solid(dev); 444 interf = create_interface(dev, solid, fluid); 445 scn = create_scene(dev, interf); 446 447 solid_param = sdis_data_get(sdis_medium_get_data(solid)); 448 interf_param = sdis_data_get(sdis_interface_get_data(interf)); 449 450 OK(sdis_scene_get_medium_spread(scn, solid, &spread)); 451 OK(sdis_scene_get_aabb(scn, lower, upper)); 452 453 /* Compute the delta of the solid random walk */ 454 solid_param->delta = 0.4 / spread; 455 456 /* Setup the upper boundary required to solve the trilinear profile */ 457 interf_param->upper[0] = upper[0]; 458 interf_param->upper[1] = upper[1]; 459 interf_param->upper[2] = upper[2]; 460 461 /* Launch probe estimation with trilinear profile set at interfaces */ 462 solid_param->delta = 0.4 / spread; 463 interf_param->profile = PROFILE_TRILINEAR; 464 check_probe(scn, PROFILE_TRILINEAR, SDIS_DIFFUSION_DELTA_SPHERE); 465 solid_param->delta = 1e-5 / spread; /* Make life difficult for WoS */ 466 check_probe(scn, PROFILE_TRILINEAR, SDIS_DIFFUSION_WOS); 467 468 /* Launch probe estimation with volumetric power profile set at interfaces */ 469 solid_param->delta = 0.4 / spread; 470 solid_param->power = Pw; 471 interf_param->profile = PROFILE_VOLUMETRIC_POWER; 472 check_probe(scn, PROFILE_VOLUMETRIC_POWER, SDIS_DIFFUSION_DELTA_SPHERE); 473 solid_param->delta = 1e-5 / spread; /* Make life difficult for WoS */ 474 check_probe(scn, PROFILE_VOLUMETRIC_POWER, SDIS_DIFFUSION_WOS); 475 solid_param->power = SDIS_VOLUMIC_POWER_NONE; 476 477 /* Launch medium integration. Do not use an analytic profile as a boundary 478 * condition but a Robin boundary condition */ 479 interf_param->profile = PROFILE_UNKNOWN; 480 solid_param->delta = 0.4 / spread; 481 check_medium(scn, solid, SDIS_DIFFUSION_DELTA_SPHERE); 482 483 /* Contrary to previous WoS tests, don't reduce the delta parameter to avoid 484 * prohibitive increases in calculation time: too small a delta would trap the 485 * path in the solid */ 486 check_medium(scn, solid, SDIS_DIFFUSION_WOS); 487 488 /* Release */ 489 OK(sdis_device_ref_put(dev)); 490 OK(sdis_medium_ref_put(fluid)); 491 OK(sdis_medium_ref_put(solid)); 492 OK(sdis_interface_ref_put(interf)); 493 OK(sdis_scene_ref_put(scn)); 494 495 CHK(mem_allocated_size() == 0); 496 return 0; 497 }