test_sdis_flux2.c (17706B)
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 /* This test consists in solving the temperature profile in a solid slab 20 * surrounded by two different convective and radiative temperatures. The 21 * conductivity of the solid material is known, as well as its thickness. A net 22 * flux is fixed on the left side of the slab. 23 * 24 * 25 * Y ////(0.1,1,1) 280K 26 * | +----------+-------+ (1.1,1,1) 27 * o--- X /##########/'E=1 /| 28 * / +----------+-------+ | 29 * Z /_ |##########|*' _\ | | 30 * ---> \ \ E=1|##########|*'/ / |280K 31 * 320K-> \__/ |##########|*'\__/ | | 32 * ---> 330K |##########|*'290K | | 33 * --\|##########|*+.... |.+ 34 * 10000W.m^-2 --/|##########|/ |/ 35 * (-1,-1,-1) +----------+-------+ 36 * (0,-1,-1)/////// 280K 37 * 38 */ 39 40 enum interface_type { 41 ADIABATIC, 42 FIXED_TEMPERATURE, 43 SOLID_FLUID_WITH_FLUX, 44 SOLID_FLUID, 45 INTERFACES_COUNT__ 46 }; 47 48 struct probe { 49 double x; 50 double time; 51 double ref; /* Analytically computed temperature */ 52 }; 53 54 /******************************************************************************* 55 * Geometry 56 ******************************************************************************/ 57 struct geometry { 58 const double* positions; 59 const size_t* indices; 60 struct sdis_interface** interfaces; 61 }; 62 63 static const double vertices_3d[12/*#vertices*/*3/*#coords per vertex*/] = { 64 0.0,-1.0,-1.0, 65 0.1,-1.0,-1.0, 66 0.0, 1.0,-1.0, 67 0.1, 1.0,-1.0, 68 0.0,-1.0, 1.0, 69 0.1,-1.0, 1.0, 70 0.0, 1.0, 1.0, 71 0.1, 1.0, 1.0, 72 1.1,-1.0,-1.0, 73 1.1, 1.0,-1.0, 74 1.1,-1.0, 1.0, 75 1.1, 1.0, 1.0 76 }; 77 static const size_t nvertices_3d = sizeof(vertices_3d) / (sizeof(double)*3); 78 79 static const size_t indices_3d[22/*#triangles*/*3/*#indices per triangle*/] = { 80 0, 2, 1, 1, 2, 3, /* Solid -Z */ 81 0, 4, 2, 2, 4, 6, /* Solid -X */ 82 4, 5, 6, 6, 5, 7, /* Solid +Z */ 83 3, 7, 1, 1, 7, 5, /* Solid +X */ 84 2, 6, 3, 3, 6, 7, /* Solid +Y */ 85 0, 1, 4, 4, 1, 5, /* Solid -Y */ 86 87 1, 3, 8, 8, 3, 9, /* Right fluid -Z */ 88 5, 10, 7, 7, 10, 11, /* Right fluid +Z */ 89 9, 11, 8, 8, 11, 10, /* Right fluid +X */ 90 3, 7, 9, 9, 7, 11, /* Right fluid +Y */ 91 1, 8, 5, 5, 8, 10 /* Right fluid -Y */ 92 }; 93 static const size_t nprimitives_3d = sizeof(indices_3d) / (sizeof(size_t)*3); 94 95 static void 96 get_indices_3d(const size_t itri, size_t ids[3], void* ctx) 97 { 98 struct geometry* geom = ctx; 99 CHK(ctx != NULL); 100 ids[0] = geom->indices[itri*3+0]; 101 ids[1] = geom->indices[itri*3+1]; 102 ids[2] = geom->indices[itri*3+2]; 103 } 104 105 static void 106 get_position_3d(const size_t ivert, double pos[3], void* ctx) 107 { 108 struct geometry* geom = ctx; 109 CHK(ctx != NULL); 110 pos[0] = geom->positions[ivert*3+0]; 111 pos[1] = geom->positions[ivert*3+1]; 112 pos[2] = geom->positions[ivert*3+2]; 113 } 114 115 static void 116 get_interface(const size_t iprim, struct sdis_interface** bound, void* ctx) 117 { 118 struct geometry* geom = ctx; 119 CHK(ctx != NULL); 120 *bound = geom->interfaces[iprim]; 121 } 122 123 /******************************************************************************* 124 * Solid media 125 ******************************************************************************/ 126 struct solid { 127 double lambda; 128 double rho; 129 double cp; 130 double volumic_power; 131 }; 132 133 static double 134 solid_get_calorific_capacity 135 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 136 { 137 const struct solid* solid = sdis_data_cget(data); 138 CHK(vtx && solid); 139 return solid->cp; 140 } 141 142 static double 143 solid_get_thermal_conductivity 144 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 145 { 146 const struct solid* solid = sdis_data_cget(data); 147 CHK(vtx && solid); 148 return solid->lambda; 149 } 150 151 static double 152 solid_get_volumic_mass 153 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 154 { 155 const struct solid* solid = sdis_data_cget(data); 156 CHK(vtx && solid); 157 return solid->rho; 158 } 159 160 static double 161 solid_get_delta 162 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 163 { 164 CHK(vtx && data); 165 return 0.005; 166 } 167 168 static double 169 solid_get_temperature 170 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 171 { 172 const struct solid* solid = sdis_data_cget(data); 173 CHK(vtx && solid); 174 175 if(vtx->time > 0) { 176 return SDIS_TEMPERATURE_NONE; 177 } else { 178 /* The initial temperature is a linear profile between T1 and T2, where T1 179 * and T2 are the temperature on the left and right slab boundary, 180 * respectively. */ 181 const double T1 = 306.334; 182 const double T2 = 294.941; 183 const double u = CLAMP(vtx->P[0] / 0.1, 0.0, 1.0); 184 return u*(T2 - T1) + T1; 185 } 186 } 187 188 static void 189 create_solid 190 (struct sdis_device* dev, 191 const struct solid* solid_props, 192 struct sdis_medium** solid) 193 { 194 struct sdis_data* data = NULL; 195 struct sdis_solid_shader shader = SDIS_SOLID_SHADER_NULL; 196 CHK(dev && solid_props && solid); 197 198 OK(sdis_data_create 199 (dev, sizeof(struct solid), ALIGNOF(struct solid), NULL, &data)); 200 memcpy(sdis_data_get(data), solid_props, sizeof(struct solid)); 201 shader.calorific_capacity = solid_get_calorific_capacity; 202 shader.thermal_conductivity = solid_get_thermal_conductivity; 203 shader.volumic_mass = solid_get_volumic_mass; 204 shader.delta = solid_get_delta; 205 shader.temperature = solid_get_temperature; 206 OK(sdis_solid_create(dev, &shader, data, solid)); 207 OK(sdis_data_ref_put(data)); 208 } 209 210 /******************************************************************************* 211 * Fluid media 212 ******************************************************************************/ 213 struct fluid { 214 double temperature; 215 }; 216 217 static double 218 fluid_get_temperature 219 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 220 { 221 const struct fluid* fluid = sdis_data_cget(data); 222 CHK(vtx && fluid); 223 return fluid->temperature; 224 } 225 226 static void 227 create_fluid 228 (struct sdis_device* dev, 229 const struct fluid* fluid_props, 230 struct sdis_medium** fluid) 231 { 232 struct sdis_data* data = NULL; 233 struct sdis_fluid_shader shader = DUMMY_FLUID_SHADER; 234 CHK(dev && fluid_props && fluid); 235 236 OK(sdis_data_create 237 (dev, sizeof(struct fluid), ALIGNOF(struct fluid), NULL, &data)); 238 memcpy(sdis_data_get(data), fluid_props, sizeof(struct fluid)); 239 shader.temperature = fluid_get_temperature; 240 OK(sdis_fluid_create(dev, &shader, data, fluid)); 241 OK(sdis_data_ref_put(data)); 242 } 243 244 /******************************************************************************* 245 * Interface 246 ******************************************************************************/ 247 struct interf { 248 double h; 249 double emissivity; 250 double phi; 251 double temperature; 252 double Tref; 253 }; 254 255 static double 256 interface_get_temperature 257 (const struct sdis_interface_fragment* frag, struct sdis_data* data) 258 { 259 const struct interf* interf = sdis_data_cget(data); 260 CHK(frag && data); 261 return interf->temperature; 262 } 263 264 static double 265 interface_get_reference_temperature 266 (const struct sdis_interface_fragment* frag, struct sdis_data* data) 267 { 268 const struct interf* interf = sdis_data_cget(data); 269 CHK(frag && interf); 270 return interf->Tref; 271 } 272 273 static double 274 interface_get_convection_coef 275 (const struct sdis_interface_fragment* frag, struct sdis_data* data) 276 { 277 const struct interf* interf = sdis_data_cget(data); 278 CHK(frag && interf); 279 return interf->h; 280 } 281 282 static double 283 interface_get_emissivity 284 (const struct sdis_interface_fragment* frag, 285 const unsigned source_id, 286 struct sdis_data* data) 287 { 288 const struct interf* interf = sdis_data_cget(data); 289 (void)source_id; 290 CHK(frag && interf); 291 return interf->emissivity; 292 } 293 294 static double 295 interface_get_specular_fraction 296 (const struct sdis_interface_fragment* frag, 297 const unsigned source_id, 298 struct sdis_data* data) 299 { 300 (void)source_id; 301 CHK(frag && data); 302 return 0; /* Unused */ 303 } 304 305 static double 306 interface_get_flux 307 (const struct sdis_interface_fragment* frag, struct sdis_data* data) 308 { 309 const struct interf* interf = sdis_data_cget(data); 310 CHK(frag && interf); 311 return interf->phi; 312 } 313 314 static void 315 create_interface 316 (struct sdis_device* dev, 317 struct sdis_medium* front, 318 struct sdis_medium* back, 319 const struct interf* interf, 320 struct sdis_interface** out_interf) 321 { 322 struct sdis_interface_shader shader = SDIS_INTERFACE_SHADER_NULL; 323 struct sdis_data* data = NULL; 324 325 CHK(interf != NULL); 326 327 shader.front.temperature = interface_get_temperature; 328 shader.back.temperature = interface_get_temperature; 329 shader.front.flux = interface_get_flux; 330 shader.back.flux = interface_get_flux; 331 332 if(sdis_medium_get_type(front) != sdis_medium_get_type(back)) { 333 shader.convection_coef = interface_get_convection_coef; 334 } 335 if(sdis_medium_get_type(front) == SDIS_FLUID) { 336 shader.front.emissivity = interface_get_emissivity; 337 shader.front.specular_fraction = interface_get_specular_fraction; 338 shader.front.reference_temperature = interface_get_reference_temperature; 339 } 340 if(sdis_medium_get_type(back) == SDIS_FLUID) { 341 shader.back.emissivity = interface_get_emissivity; 342 shader.back.specular_fraction = interface_get_specular_fraction; 343 shader.back.reference_temperature = interface_get_reference_temperature; 344 } 345 shader.convection_coef_upper_bound = MMAX(0, interf->h); 346 347 OK(sdis_data_create 348 (dev, sizeof(struct interf), ALIGNOF(struct interf), NULL, &data)); 349 memcpy(sdis_data_get(data), interf, sizeof(*interf)); 350 351 OK(sdis_interface_create(dev, front, back, &shader, data, out_interf)); 352 OK(sdis_data_ref_put(data)); 353 } 354 355 /******************************************************************************* 356 * Create the radiative environment 357 ******************************************************************************/ 358 static double 359 radenv_get_temperature 360 (const struct sdis_radiative_ray* ray, 361 struct sdis_data* data) 362 { 363 (void)ray, (void)data; 364 return 320; /* [K] */ 365 } 366 367 static double 368 radenv_get_reference_temperature 369 (const struct sdis_radiative_ray* ray, 370 struct sdis_data* data) 371 { 372 (void)ray, (void)data; 373 return 300; /* [K] */ 374 } 375 376 static struct sdis_radiative_env* 377 create_radenv(struct sdis_device* sdis) 378 { 379 struct sdis_radiative_env_shader shader = SDIS_RADIATIVE_ENV_SHADER_NULL; 380 struct sdis_radiative_env* radenv = NULL; 381 382 shader.temperature = radenv_get_temperature; 383 shader.reference_temperature = radenv_get_reference_temperature; 384 OK(sdis_radiative_env_create(sdis, &shader, NULL, &radenv)); 385 return radenv; 386 } 387 388 /******************************************************************************* 389 * Create scene 390 ******************************************************************************/ 391 static void 392 create_scene_3d 393 (struct sdis_device* dev, 394 struct sdis_interface* interfaces[INTERFACES_COUNT__], 395 struct sdis_radiative_env* radenv, 396 struct sdis_scene** scn) 397 { 398 struct geometry geom; 399 struct sdis_interface* prim_interfaces[32]; 400 struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT; 401 402 CHK(dev && interfaces && radenv && scn); 403 404 /* Setup the per primitive interface of the solid medium */ 405 prim_interfaces[0] = prim_interfaces[1] = interfaces[ADIABATIC]; 406 prim_interfaces[2] = prim_interfaces[3] = interfaces[SOLID_FLUID_WITH_FLUX]; 407 prim_interfaces[4] = prim_interfaces[5] = interfaces[ADIABATIC]; 408 prim_interfaces[6] = prim_interfaces[7] = interfaces[SOLID_FLUID]; 409 prim_interfaces[8] = prim_interfaces[9] = interfaces[ADIABATIC]; 410 prim_interfaces[10] = prim_interfaces[11] = interfaces[ADIABATIC]; 411 412 /* Setup the per primitive interface for the right fluid */ 413 prim_interfaces[12] = prim_interfaces[13] = interfaces[FIXED_TEMPERATURE]; 414 prim_interfaces[14] = prim_interfaces[15] = interfaces[FIXED_TEMPERATURE]; 415 prim_interfaces[16] = prim_interfaces[17] = interfaces[FIXED_TEMPERATURE]; 416 prim_interfaces[18] = prim_interfaces[19] = interfaces[FIXED_TEMPERATURE]; 417 prim_interfaces[20] = prim_interfaces[21] = interfaces[FIXED_TEMPERATURE]; 418 419 /* Create the scene */ 420 geom.positions = vertices_3d; 421 geom.indices = indices_3d; 422 geom.interfaces = prim_interfaces; 423 scn_args.get_indices = get_indices_3d; 424 scn_args.get_interface = get_interface; 425 scn_args.get_position = get_position_3d; 426 scn_args.nprimitives = nprimitives_3d; 427 scn_args.nvertices = nvertices_3d; 428 scn_args.t_range[0] = 300; 429 scn_args.t_range[1] = 300; 430 scn_args.context = &geom; 431 scn_args.radenv = radenv; 432 OK(sdis_scene_create(dev, &scn_args, scn)); 433 } 434 435 /******************************************************************************* 436 * Helper functions 437 ******************************************************************************/ 438 static void 439 check(struct sdis_scene* scn, const struct probe* probe) 440 { 441 struct sdis_solve_probe_args args = SDIS_SOLVE_PROBE_ARGS_DEFAULT; 442 struct sdis_estimator* estimator = NULL; 443 struct sdis_mc T = SDIS_MC_NULL; 444 445 CHK(scn && probe); 446 447 args.nrealisations = 10000; 448 args.picard_order = 1; 449 args.position[0] = probe->x; 450 args.position[1] = 0; 451 args.position[2] = 0; 452 args.time_range[0] = probe->time; 453 args.time_range[1] = probe->time; 454 455 OK(sdis_solve_probe(scn, &args, &estimator)); 456 OK(sdis_estimator_get_temperature(estimator, &T)); 457 OK(sdis_estimator_ref_put(estimator)); 458 459 printf("Temperature at x=%g, t=%g: %g ~ %g +/- %g\n", 460 probe->x, probe->time, probe->ref, T.E, T.SE); 461 462 CHK(eq_eps(probe->ref, T.E, T.SE*3)); 463 } 464 465 /******************************************************************************* 466 * Test 467 ******************************************************************************/ 468 int 469 main(int argc, char** argv) 470 { 471 struct sdis_device* dev = NULL; 472 struct sdis_radiative_env* radenv = NULL; 473 struct sdis_scene* scn_3d = NULL; 474 struct sdis_medium* solid = NULL; 475 struct sdis_medium* dummy = NULL; 476 struct sdis_medium* fluid1 = NULL; 477 struct sdis_medium* fluid2 = NULL; 478 struct sdis_interface* interfaces[INTERFACES_COUNT__]; 479 480 struct solid solid_props; 481 struct fluid fluid_props; 482 struct interf interf_props; 483 484 const struct probe probes[] = { 485 {0.01, 1000.0, 481.72005748728628}, 486 {0.05, 1000.0, 335.19469995601020}, 487 {0.09, 1000.0, 299.94436943411478}, 488 {0.01, 2000.0, 563.21759568607558}, 489 {0.05, 2000.0, 392.79827670626440}, 490 {0.09, 2000.0, 324.89742556243448}, 491 {0.01, 3000.0, 620.25242712533577}, 492 {0.05, 3000.0, 444.73414407361213}, 493 {0.09, 3000.0, 359.44045704073852}, 494 {0.01, 4000.0, 665.65935222224493}, 495 {0.05, 4000.0, 490.32470982110840}, 496 {0.09, 4000.0, 393.89924931902408}, 497 {0.01, 10000.0, 830.44439052891505}, 498 {0.05, 10000.0, 664.82771620162805}, 499 {0.09, 10000.0, 533.92442748613928} 500 }; 501 const size_t nprobes = sizeof(probes)/sizeof(*probes); 502 size_t iprobe; 503 504 (void)argc, (void)argv; 505 506 OK(sdis_device_create(&SDIS_DEVICE_CREATE_ARGS_DEFAULT, &dev)); 507 508 radenv = create_radenv(dev); 509 510 /* Solid medium */ 511 solid_props.lambda = 1.15; 512 solid_props.rho = 1700; 513 solid_props.cp = 800; 514 create_solid(dev, &solid_props, &solid); 515 516 /* Dummy solid medium */ 517 solid_props.lambda = 0; 518 solid_props.rho = 1700; 519 solid_props.cp = 800; 520 create_solid(dev, &solid_props, &dummy); 521 522 /* Fluid media */ 523 fluid_props.temperature = 330; 524 create_fluid(dev, &fluid_props, &fluid1); 525 fluid_props.temperature = 290; 526 create_fluid(dev, &fluid_props, &fluid2); 527 528 /* Adiabatic interfaces */ 529 interf_props.h = 0; 530 interf_props.emissivity = 0; 531 interf_props.phi = SDIS_FLUX_NONE; 532 interf_props.temperature = SDIS_TEMPERATURE_NONE; 533 interf_props.Tref = SDIS_TEMPERATURE_NONE; 534 create_interface(dev, solid, dummy, &interf_props, &interfaces[ADIABATIC]); 535 536 /* Interfaces with a fixed temperature */ 537 interf_props.h = 1; 538 interf_props.emissivity = 1; 539 interf_props.phi = SDIS_FLUX_NONE; 540 interf_props.temperature = 280; 541 interf_props.Tref = 300; 542 create_interface 543 (dev, fluid2, dummy, &interf_props, &interfaces[FIXED_TEMPERATURE]); 544 545 /* Interfaces with a fixed flux */ 546 interf_props.h = 2; 547 interf_props.emissivity = 1; 548 interf_props.phi = 10000; 549 interf_props.temperature = SDIS_TEMPERATURE_NONE; 550 interf_props.Tref = 300; 551 create_interface 552 (dev, solid, fluid1, &interf_props, &interfaces[SOLID_FLUID_WITH_FLUX]); 553 554 interf_props.h = 8; 555 interf_props.emissivity = 1; 556 interf_props.phi = SDIS_FLUX_NONE; 557 interf_props.temperature = SDIS_TEMPERATURE_NONE; 558 interf_props.Tref = 300; 559 create_interface 560 (dev, solid, fluid2, &interf_props, &interfaces[SOLID_FLUID]); 561 562 create_scene_3d(dev, interfaces, radenv, &scn_3d); 563 564 FOR_EACH(iprobe, 0, nprobes) { 565 check(scn_3d, &probes[iprobe]); 566 } 567 568 /* Release memory */ 569 OK(sdis_radiative_env_ref_put(radenv)); 570 OK(sdis_scene_ref_put(scn_3d)); 571 OK(sdis_medium_ref_put(solid)); 572 OK(sdis_medium_ref_put(dummy)); 573 OK(sdis_medium_ref_put(fluid1)); 574 OK(sdis_medium_ref_put(fluid2)); 575 OK(sdis_interface_ref_put(interfaces[ADIABATIC])); 576 OK(sdis_interface_ref_put(interfaces[FIXED_TEMPERATURE])); 577 OK(sdis_interface_ref_put(interfaces[SOLID_FLUID_WITH_FLUX])); 578 OK(sdis_interface_ref_put(interfaces[SOLID_FLUID])); 579 OK(sdis_device_ref_put(dev)); 580 581 CHK(mem_allocated_size() == 0); 582 return 0; 583 }