test_sdis_external_flux.c (19357B)
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 "test_sdis_utils.h" 17 #include "sdis.h" 18 19 #include <rsys/rsys.h> 20 21 /* 22 * The system consists of 2 parallelepipeds: a vertical one called the wall, and 23 * a horizontal one representing the floor. The wall is a black body, while the 24 * floor is a perfectly reflective surface. The surrounding fluid has a fixed 25 * temperature and, finally, an external spherical source represents the sun. 26 * This test calculates the steady temperature at a position in the wall and 27 * compares it with the analytical solution given for a perfectly diffuse or 28 * specular ground. 29 * 30 * (-0.1,1500) 31 * +---+ External source 32 * | E=1 T_FLUID ## 33 * Probe x | _\ #### 34 * | | / / ## 35 * +---+ \__/ 36 * (0,500) 37 * 38 * (0,0) 39 * Y *--------E=0------------- - - - 40 * | | 41 * o--X *------------------------ - - - 42 * / (0,-1) 43 * Z 44 * 45 */ 46 47 #define T_FLUID 300.0 /* [K] */ 48 #define T_REF 300.0 /* [K] */ 49 50 /******************************************************************************* 51 * Geometries 52 ******************************************************************************/ 53 static const double positions_2d[] = { 54 /* Ground */ 55 1.0e12, -1.0, 56 0.0, -1.0, 57 0.0, 0.0, 58 1.0e12, 0.0, 59 60 /* Wall */ 61 0.0, 500.0, 62 -0.1, 500.0, 63 -0.1, 1500.0, 64 0.0, 1500.0 65 }; 66 static const size_t nvertices_2d = sizeof(positions_2d) / (sizeof(double)*2); 67 68 static const double positions_3d[] = { 69 /* Ground */ 70 0.0, -1.0, -1.0e6, 71 1.0e12, -1.0, -1.0e6, 72 0.0, 1.0, -1.0e6, 73 1.0e12, 1.0, -1.0e6, 74 0.0, -1.0, 1.0e6, 75 1.0e12, -1.0, 1.0e6, 76 0.0, 1.0, 1.0e6, 77 1.0e12, 1.0, 1.0e6, 78 79 /* Wall */ 80 -0.1, 500.0, -500.0, 81 0.0, 500.0, -500.0, 82 -0.1, 1500.0, -500.0, 83 0.0, 1500.0, -500.0, 84 -0.1, 500.0, 500.0, 85 0.0, 500.0, 500.0, 86 -0.1, 1500.0, 500.0, 87 0.0, 1500.0, 500.0 88 }; 89 static const size_t nvertices_3d = sizeof(positions_3d) / (sizeof(double)*3); 90 91 static const size_t indices_2d[] = { 92 /* Ground */ 93 0, 1, /* -y */ 94 1, 2, /* -x */ 95 2, 3, /* +y */ 96 3, 0, /* +x */ 97 98 /* Wall */ 99 4, 5, /* -y */ 100 5, 6, /* -x */ 101 6, 7, /* +y */ 102 7, 4 /* +x */ 103 }; 104 static const size_t nsegments = sizeof(indices_2d) / (sizeof(size_t)*2); 105 106 static const size_t indices_3d[] = { 107 /* Ground */ 108 0, 2, 1, 1, 2, 3, /* -z */ 109 0, 4, 2, 2, 4, 6, /* -x */ 110 4, 5, 6, 6, 5, 7, /* +z */ 111 3, 7, 5, 5, 1, 3, /* +x */ 112 2, 6, 7, 7, 3, 2, /* +y */ 113 0, 1, 5, 5, 4, 0, /* -y */ 114 115 /* Wall */ 116 8, 10, 9, 9, 10, 11, /* -z */ 117 8, 12, 10, 10, 12, 14, /* -x */ 118 12, 13, 14, 14, 13, 15, /* +z */ 119 11, 15, 13, 13, 9, 11, /* +x */ 120 10, 14, 15, 15, 11, 10, /* +y */ 121 8, 9, 13, 13, 12, 8 /* -y */ 122 }; 123 static const size_t ntriangles = sizeof(indices_3d) / (sizeof(size_t)*3); 124 125 /******************************************************************************* 126 * Media 127 ******************************************************************************/ 128 #define MEDIUM_PROP(Type, Prop, Val) \ 129 static double \ 130 Type##_get_##Prop \ 131 (const struct sdis_rwalk_vertex* vtx, \ 132 struct sdis_data* data) \ 133 { \ 134 (void)vtx, (void)data; /* Avoid the "unused variable" warning */ \ 135 return Val; \ 136 } 137 MEDIUM_PROP(medium, volumic_mass, 1700) /* [kj/m^3] */ 138 MEDIUM_PROP(medium, calorific_capacity, 800) /* [J/K/Kg] */ 139 MEDIUM_PROP(solid, thermal_conductivity, 1.15) /* [W/m/K] */ 140 MEDIUM_PROP(solid, delta, 0.1/20.0) /* [m] */ 141 MEDIUM_PROP(solid, temperature, SDIS_TEMPERATURE_NONE/*<=> unknown*/) /* [K] */ 142 MEDIUM_PROP(fluid, temperature, T_FLUID) /* [K] */ 143 #undef MEDIUM_PROP 144 145 static struct sdis_medium* 146 create_solid(struct sdis_device* sdis) 147 { 148 struct sdis_solid_shader shader = SDIS_SOLID_SHADER_NULL; 149 struct sdis_medium* solid = NULL; 150 151 shader.calorific_capacity = medium_get_calorific_capacity; 152 shader.thermal_conductivity = solid_get_thermal_conductivity; 153 shader.volumic_mass = medium_get_volumic_mass; 154 shader.delta = solid_get_delta; 155 shader.temperature = solid_get_temperature; 156 OK(sdis_solid_create(sdis, &shader, NULL, &solid)); 157 return solid; 158 } 159 160 static struct sdis_medium* 161 create_fluid(struct sdis_device* sdis) 162 { 163 struct sdis_fluid_shader shader = SDIS_FLUID_SHADER_NULL; 164 struct sdis_medium* fluid = NULL; 165 166 shader.calorific_capacity = medium_get_calorific_capacity; 167 shader.volumic_mass = medium_get_volumic_mass; 168 shader.temperature = fluid_get_temperature; 169 OK(sdis_fluid_create(sdis, &shader, NULL, &fluid)); 170 return fluid; 171 } 172 173 /******************************************************************************* 174 * Interfaces 175 ******************************************************************************/ 176 struct interface { 177 double emissivity; 178 double specular_fraction; 179 double convection_coef; /* [W/m^2/K] */ 180 }; 181 182 #define INTERF_PROP(Prop, Val) \ 183 static double \ 184 interface_get_##Prop \ 185 (const struct sdis_interface_fragment* frag, \ 186 struct sdis_data* data) \ 187 { \ 188 (void)frag, (void)data; /* Avoid the "unused variable" warning */ \ 189 return Val; \ 190 } 191 INTERF_PROP(temperature, SDIS_TEMPERATURE_NONE/*<=> unknown*/) /* [K] */ 192 INTERF_PROP(reference_temperature, T_REF) /* [K] */ 193 #undef INTERF_PROP 194 195 #define INTERF_PROP(Prop) \ 196 static double \ 197 interface_get_##Prop \ 198 (const struct sdis_interface_fragment* frag, \ 199 const unsigned source_id, \ 200 struct sdis_data* data) \ 201 { \ 202 struct interface* interf_data = NULL; \ 203 (void)frag, (void)source_id; /* Avoid the "unused variable" warning */ \ 204 interf_data = sdis_data_get(data); \ 205 return source_id == SDIS_INTERN_SOURCE_ID ? 0 : interf_data->Prop; \ 206 } 207 INTERF_PROP(emissivity) 208 INTERF_PROP(specular_fraction) 209 #undef INTERF_PROP 210 211 static double /* [W/m^2/K] */ 212 interface_get_convection_coef 213 (const struct sdis_interface_fragment* frag, 214 struct sdis_data* data) 215 { 216 struct interface* interf_data = NULL; 217 (void)frag; /* Avoid the "unused variable" warning */ 218 interf_data = sdis_data_get(data); 219 return interf_data->convection_coef; 220 } 221 222 static struct sdis_interface* 223 create_interface 224 (struct sdis_device* sdis, 225 struct sdis_medium* front, 226 struct sdis_medium* back, 227 const double emissivity, 228 const double convection_coef, 229 struct interface** out_interf_data) /* May be NULL */ 230 { 231 struct sdis_data* data = NULL; 232 struct sdis_interface* interf = NULL; 233 struct sdis_interface_shader shader = SDIS_INTERFACE_SHADER_NULL; 234 struct interface* interf_data = NULL; 235 236 OK(sdis_data_create 237 (sdis, sizeof(struct interface), ALIGNOF(struct interface), NULL, &data)); 238 interf_data = sdis_data_get(data); 239 interf_data->emissivity = emissivity; 240 interf_data->convection_coef = convection_coef; /* [W/m^2/K] */ 241 interf_data->specular_fraction = 0; /* Diffuse */ 242 if(out_interf_data) *out_interf_data = interf_data; 243 244 shader.front.temperature = interface_get_temperature; 245 shader.back.temperature = interface_get_temperature; 246 shader.back.reference_temperature = interface_get_reference_temperature; 247 shader.back.emissivity = interface_get_emissivity; 248 shader.back.specular_fraction = interface_get_specular_fraction; 249 shader.convection_coef = interface_get_convection_coef; 250 shader.convection_coef_upper_bound = convection_coef; 251 OK(sdis_interface_create(sdis, front, back, &shader, data, &interf)); 252 OK(sdis_data_ref_put(data)); 253 254 return interf; 255 } 256 257 /******************************************************************************* 258 * External source 259 ******************************************************************************/ 260 static void 261 source_get_position 262 (const double time, 263 double pos[3], 264 struct sdis_data* data) 265 { 266 const double elevation = MDEG2RAD(30); /* [radian] */ 267 const double distance = 1.5e11; /* [m] */ 268 (void)time, (void)data; /* Avoid the "unusued variable" warning */ 269 270 pos[0] = cos(elevation) * distance; 271 pos[1] = sin(elevation) * distance; 272 pos[2] = 0; 273 } 274 275 static double 276 source_get_power(const double time, struct sdis_data* data) 277 { 278 (void)time, (void)data; /* Avoid the "unusued variable" warning */ 279 return 3.845e26; /* [W] */ 280 } 281 282 static struct sdis_source* 283 create_source(struct sdis_device* sdis) 284 { 285 struct sdis_spherical_source_shader shader = SDIS_SPHERICAL_SOURCE_SHADER_NULL; 286 struct sdis_source* src = NULL; 287 288 shader.position = source_get_position; 289 shader.power = source_get_power; 290 shader.radius = 6.5991756e8; /* [m] */ 291 OK(sdis_spherical_source_create(sdis, &shader, NULL, &src)); 292 return src; 293 } 294 295 /******************************************************************************* 296 * Radiative environment 297 ******************************************************************************/ 298 static double 299 radenv_get_temperature 300 (const struct sdis_radiative_ray* ray, 301 struct sdis_data* data) 302 { 303 (void)ray, (void)data; 304 return 0; /* [K] */ 305 } 306 307 static double 308 radenv_get_reference_temperature 309 (const struct sdis_radiative_ray* ray, 310 struct sdis_data* data) 311 { 312 (void)ray, (void)data; 313 return T_REF; /* [K] */ 314 } 315 316 static struct sdis_radiative_env* 317 create_radenv(struct sdis_device* sdis) 318 { 319 struct sdis_radiative_env_shader shader = SDIS_RADIATIVE_ENV_SHADER_NULL; 320 struct sdis_radiative_env* radenv = NULL; 321 322 shader.temperature = radenv_get_temperature; 323 shader.reference_temperature = radenv_get_reference_temperature; 324 OK(sdis_radiative_env_create(sdis, &shader, NULL, &radenv)); 325 return radenv; 326 } 327 328 /******************************************************************************* 329 * Scene 330 ******************************************************************************/ 331 struct scene_context { 332 struct sdis_interface* interf_ground; 333 struct sdis_interface* interf_wall; 334 }; 335 static const struct scene_context SCENE_CONTEXT_NULL = {NULL, NULL}; 336 337 static void 338 scene_get_indices_2d(const size_t iseg, size_t ids[2], void* ctx) 339 { 340 struct scene_context* context = ctx; 341 CHK(ids && context && iseg < nsegments); 342 ids[0] = (unsigned)indices_2d[iseg*2+0]; 343 ids[1] = (unsigned)indices_2d[iseg*2+1]; 344 } 345 346 static void 347 scene_get_indices_3d(const size_t itri, size_t ids[3], void* ctx) 348 { 349 struct scene_context* context = ctx; 350 CHK(ids && context && itri < ntriangles); 351 ids[0] = (unsigned)indices_3d[itri*3+0]; 352 ids[1] = (unsigned)indices_3d[itri*3+1]; 353 ids[2] = (unsigned)indices_3d[itri*3+2]; 354 } 355 356 static void 357 scene_get_interface_2d(const size_t iseg, struct sdis_interface** interf, void* ctx) 358 { 359 struct scene_context* context = ctx; 360 CHK(interf && context && iseg < nsegments); 361 *interf = iseg < 4 ? context->interf_ground : context->interf_wall; 362 } 363 364 static void 365 scene_get_interface_3d(const size_t itri, struct sdis_interface** interf, void* ctx) 366 { 367 struct scene_context* context = ctx; 368 CHK(interf && context && itri < ntriangles); 369 *interf = itri < 12 ? context->interf_ground : context->interf_wall; 370 } 371 372 static void 373 scene_get_position_2d(const size_t ivert, double pos[2], void* ctx) 374 { 375 struct scene_context* context = ctx; 376 CHK(pos && context && ivert < nvertices_2d); 377 pos[0] = positions_2d[ivert*2+0]; 378 pos[1] = positions_2d[ivert*2+1]; 379 } 380 381 static void 382 scene_get_position_3d(const size_t ivert, double pos[3], void* ctx) 383 { 384 struct scene_context* context = ctx; 385 CHK(pos && context && ivert < nvertices_3d); 386 pos[0] = positions_3d[ivert*3+0]; 387 pos[1] = positions_3d[ivert*3+1]; 388 pos[2] = positions_3d[ivert*3+2]; 389 } 390 391 static struct sdis_scene* 392 create_scene_2d 393 (struct sdis_device* sdis, 394 struct sdis_interface* interf_ground, 395 struct sdis_interface* interf_wall, 396 struct sdis_source* source, 397 struct sdis_radiative_env* radenv) 398 { 399 struct sdis_scene* scn = NULL; 400 struct sdis_source* src = NULL; 401 struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT; 402 struct scene_context context = SCENE_CONTEXT_NULL; 403 404 context.interf_ground = interf_ground; 405 context.interf_wall = interf_wall; 406 407 scn_args.get_indices = scene_get_indices_2d; 408 scn_args.get_interface = scene_get_interface_2d; 409 scn_args.get_position = scene_get_position_2d; 410 scn_args.nprimitives = nsegments; 411 scn_args.nvertices = nvertices_2d; 412 scn_args.t_range[0] = T_REF; /* [K] */ 413 scn_args.t_range[1] = T_REF; /* [K] */ 414 scn_args.source = source; 415 scn_args.radenv = radenv; 416 scn_args.context = &context; 417 OK(sdis_scene_2d_create(sdis, &scn_args, &scn)); 418 419 BA(sdis_scene_get_source(NULL, &src)); 420 BA(sdis_scene_get_source(scn, NULL)); 421 OK(sdis_scene_get_source(scn, &src)); 422 CHK(src == source); 423 424 return scn; 425 } 426 427 static struct sdis_scene* 428 create_scene_3d 429 (struct sdis_device* sdis, 430 struct sdis_interface* interf_ground, 431 struct sdis_interface* interf_wall, 432 struct sdis_source* source, 433 struct sdis_radiative_env* radenv) 434 { 435 struct sdis_scene* scn = NULL; 436 struct sdis_source* src = NULL; 437 struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT; 438 struct scene_context context = SCENE_CONTEXT_NULL; 439 440 context.interf_ground = interf_ground; 441 context.interf_wall = interf_wall; 442 443 scn_args.get_indices = scene_get_indices_3d; 444 scn_args.get_interface = scene_get_interface_3d; 445 scn_args.get_position = scene_get_position_3d; 446 scn_args.nprimitives = ntriangles; 447 scn_args.nvertices = nvertices_3d; 448 scn_args.t_range[0] = 0; /* [K] */ 449 scn_args.t_range[1] = 0; /* [K] */ 450 scn_args.source = source; 451 scn_args.radenv = radenv; 452 scn_args.context = &context; 453 OK(sdis_scene_create(sdis, &scn_args, &scn)); 454 455 BA(sdis_scene_get_source(NULL, &src)); 456 BA(sdis_scene_get_source(scn, NULL)); 457 OK(sdis_scene_get_source(scn, &src)); 458 CHK(src == source); 459 460 return scn; 461 } 462 463 /******************************************************************************* 464 * Validations 465 ******************************************************************************/ 466 static void 467 check 468 (struct sdis_scene* scn, 469 const size_t nrealisations, 470 const double analytical_ref, 471 const int is_master_process) 472 { 473 struct sdis_solve_probe_args probe_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT; 474 struct sdis_mc T = SDIS_MC_NULL; 475 struct sdis_estimator* estimator = NULL; 476 477 probe_args.position[0] = -0.05; 478 probe_args.position[1] = 1000; 479 probe_args.position[2] = 0; 480 probe_args.nrealisations = nrealisations; 481 OK(sdis_solve_probe(scn, &probe_args, &estimator)); 482 483 if(!is_master_process) return; 484 485 OK(sdis_estimator_get_temperature(estimator, &T)); 486 printf("T(%g, %g, %g) = %g ~ %g +/- %g\n", 487 SPLIT3(probe_args.position), analytical_ref, T.E, T.SE); 488 OK(sdis_estimator_ref_put(estimator)); 489 490 CHK(eq_eps(analytical_ref, T.E, 3*T.SE)); 491 } 492 493 static void 494 check_green 495 (struct sdis_scene* scn, 496 const size_t nrealisations, 497 const double analytical_ref, 498 const int is_master_process) 499 { 500 struct sdis_solve_probe_args probe_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT; 501 struct sdis_mc T = SDIS_MC_NULL; 502 struct sdis_green_function* green = NULL; 503 struct sdis_estimator* estimator = NULL; 504 505 probe_args.position[0] = -0.05; 506 probe_args.position[1] = 1000; 507 probe_args.position[2] = 0; 508 probe_args.nrealisations = nrealisations; 509 OK(sdis_solve_probe_green_function(scn, &probe_args, &green)); 510 511 if(!is_master_process) return; 512 513 OK(sdis_green_function_solve(green, &estimator)); 514 check_green_function(green); 515 516 OK(sdis_estimator_get_temperature(estimator, &T)); 517 518 printf("T(%g, %g, %g) = %g ~ %g +/- %g\n", 519 SPLIT3(probe_args.position), analytical_ref, T.E, T.SE); 520 OK(sdis_green_function_ref_put(green)); 521 OK(sdis_estimator_ref_put(estimator)); 522 523 CHK(eq_eps(analytical_ref, T.E, 3*T.SE)); 524 } 525 526 /******************************************************************************* 527 * The test 528 ******************************************************************************/ 529 int 530 main(int argc, char** argv) 531 { 532 struct sdis_device* dev = NULL; 533 struct sdis_medium* fluid = NULL; 534 struct sdis_medium* solid = NULL; 535 struct sdis_interface* interf_ground = NULL; 536 struct sdis_interface* interf_wall = NULL; 537 struct sdis_radiative_env* radenv = NULL; 538 struct sdis_scene* scn_2d = NULL; 539 struct sdis_scene* scn_3d = NULL; 540 struct sdis_source* src = NULL; 541 542 struct interface* ground_interf_data = NULL; 543 int is_master_process = 0; 544 (void) argc, (void)argv; 545 546 create_default_device(&argc, &argv, &is_master_process, &dev); 547 548 fluid = create_fluid(dev); 549 solid = create_solid(dev); 550 interf_ground = create_interface 551 (dev, solid, fluid, 0/*emissivity*/, 0/*h*/, &ground_interf_data); 552 interf_wall = create_interface 553 (dev, solid, fluid, 1/*emissivity*/, 10/*h*/, NULL); 554 src = create_source(dev); 555 radenv = create_radenv(dev); 556 scn_2d = create_scene_2d(dev, interf_ground, interf_wall, src, radenv); 557 scn_3d = create_scene_3d(dev, interf_ground, interf_wall, src, radenv); 558 559 ground_interf_data->specular_fraction = 0; /* Lambertian */ 560 check(scn_2d, 10000, 375.88, is_master_process); 561 check_green(scn_3d, 10000, 375.88, is_master_process); 562 563 ground_interf_data->specular_fraction = 1; /* Specular */ 564 check(scn_2d, 100000, 417.77, is_master_process); 565 check_green(scn_3d, 100000, 417.77, is_master_process); 566 567 OK(sdis_medium_ref_put(fluid)); 568 OK(sdis_medium_ref_put(solid)); 569 OK(sdis_interface_ref_put(interf_ground)); 570 OK(sdis_interface_ref_put(interf_wall)); 571 OK(sdis_radiative_env_ref_put(radenv)); 572 OK(sdis_source_ref_put(src)); 573 OK(sdis_scene_ref_put(scn_2d)); 574 OK(sdis_scene_ref_put(scn_3d)); 575 576 free_default_device(dev); 577 578 CHK(mem_allocated_size() == 0); 579 return 0; 580 }