test_sdis_conducto_radiative_2d.c (18302B)
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/ssp.h> 20 21 /* The scene is composed of a solid square whose temperature is unknown. The 22 * square segments on +/-X are in contact with a fluid and their convection 23 * coefficient is null while their emissivity is 1. The left and right fluids 24 * are enclosed by segments whose emissivity are null excepted for the segments 25 * orthogonal to the X axis that are fully emissive and whose temperature is 26 * known. The medium that surrounds the solid square and the 2 fluids is a 27 * solid with a null conductivity. 28 * 29 * (1, 1) 30 * +-----+----------+-----+ (1.5,1,1) 31 * | |##########| | 32 * | |##########| | 33 * 300K | E=1 |##########| E=1 | 310K 34 * | |##########| | 35 * | |##########| | 36 * (-1.5,-1) +-----+----------+-----+ 37 * (-1,-1) 38 */ 39 40 /******************************************************************************* 41 * Geometry 42 ******************************************************************************/ 43 struct geometry { 44 const double* positions; 45 const size_t* indices; 46 struct sdis_interface** interfaces; 47 }; 48 49 static const double vertices[8/*#vertices*/*2/*#coords par vertex*/] = { 50 1.0, -1.0, 51 -1.0, -1.0, 52 -1.0, 1.0, 53 1.0, 1.0, 54 1.5, -1.0, 55 -1.5, -1.0, 56 -1.5, 1.0, 57 1.5, 1.0 58 }; 59 static const size_t nvertices = sizeof(vertices) / (sizeof(double)*2); 60 61 static const size_t indices[10/*#segments*/*2/*#indices per segment*/] = { 62 0, 1, /* Solid bottom segment */ 63 1, 2, /* Solid left segment */ 64 2, 3, /* Solid top segment */ 65 3, 0, /* Solid right segment */ 66 67 1, 5, /* Left fluid bottom segment */ 68 5, 6, /* Left fluid left segment */ 69 6, 2, /* Left fluid top segment */ 70 71 4, 0, /* Right fluid bottom segment */ 72 3, 7, /* Right fluid top segment */ 73 7, 4 /* Right fluid right segment */ 74 }; 75 static const size_t nsegments = sizeof(indices) / (sizeof(size_t)*2); 76 77 static void 78 get_indices(const size_t iseg, size_t ids[2], void* ctx) 79 { 80 struct geometry* geom = ctx; 81 CHK(ctx != NULL); 82 ids[0] = geom->indices[iseg*2+0]; 83 ids[1] = geom->indices[iseg*2+1]; 84 } 85 86 static void 87 get_position(const size_t ivert, double pos[2], void* ctx) 88 { 89 struct geometry* geom = ctx; 90 CHK(ctx != NULL); 91 pos[0] = geom->positions[ivert*2+0]; 92 pos[1] = geom->positions[ivert*2+1]; 93 } 94 95 static void 96 get_interface(const size_t iseg, struct sdis_interface** bound, void* ctx) 97 { 98 struct geometry* geom = ctx; 99 CHK(ctx != NULL); 100 *bound = geom->interfaces[iseg]; 101 } 102 103 /******************************************************************************* 104 * Media 105 ******************************************************************************/ 106 struct solid { 107 double lambda; 108 }; 109 110 static double 111 temperature_unknown(const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 112 { 113 CHK(vtx != NULL); (void)data; 114 return SDIS_TEMPERATURE_NONE; 115 } 116 117 static double 118 solid_get_calorific_capacity 119 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 120 { 121 CHK(vtx != NULL); (void)data; 122 return 1; 123 } 124 125 static double 126 solid_get_thermal_conductivity 127 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 128 { 129 CHK(vtx != NULL); 130 CHK(data != NULL); 131 return ((const struct solid*)sdis_data_cget(data))->lambda; 132 } 133 134 static double 135 solid_get_volumic_mass 136 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 137 { 138 CHK(vtx != NULL); (void)data; 139 return 1; 140 } 141 142 static double 143 solid_get_delta 144 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 145 { 146 CHK(vtx != NULL); (void)data; 147 return 1.0/10.0; 148 } 149 150 /******************************************************************************* 151 * Interface 152 ******************************************************************************/ 153 struct interfac { 154 double convection_coef; 155 struct { 156 double temperature; 157 double emissivity; 158 double specular_fraction; 159 double reference_temperature; 160 } front, back; 161 }; 162 163 static const struct interfac INTERFACE_NULL = { 164 0, {SDIS_TEMPERATURE_NONE, -1, -1, -1}, {-1, -1, -1, -1} 165 }; 166 167 static double 168 interface_get_temperature 169 (const struct sdis_interface_fragment* frag, struct sdis_data* data) 170 { 171 const struct interfac* interf; 172 double T = SDIS_TEMPERATURE_NONE; 173 CHK(data != NULL && frag != NULL); 174 interf = sdis_data_cget(data); 175 switch(frag->side) { 176 case SDIS_FRONT: T = interf->front.temperature; break; 177 case SDIS_BACK: T = interf->back.temperature; break; 178 default: FATAL("Unreachable code.\n"); break; 179 } 180 return T; 181 } 182 183 static double 184 interface_get_convection_coef 185 (const struct sdis_interface_fragment* frag, struct sdis_data* data) 186 { 187 const struct interfac* interf; 188 CHK(data != NULL && frag != NULL); 189 interf = sdis_data_cget(data); 190 return interf->convection_coef; 191 } 192 193 static double 194 interface_get_emissivity 195 (const struct sdis_interface_fragment* frag, 196 const unsigned source_id, 197 struct sdis_data* data) 198 { 199 const struct interfac* interf; 200 double e = -1; 201 (void)source_id; 202 CHK(data != NULL && frag != NULL); 203 interf = sdis_data_cget(data); 204 switch(frag->side) { 205 case SDIS_FRONT: e = interf->front.emissivity; break; 206 case SDIS_BACK: e = interf->back.emissivity; break; 207 default: FATAL("Unreachable code.\n"); break; 208 } 209 return e; 210 } 211 212 static double 213 interface_get_specular_fraction 214 (const struct sdis_interface_fragment* frag, 215 const unsigned source_id, 216 struct sdis_data* data) 217 { 218 const struct interfac* interf; 219 double f = -1; 220 (void)source_id; 221 CHK(data != NULL && frag != NULL); 222 interf = sdis_data_cget(data); 223 switch(frag->side) { 224 case SDIS_FRONT: f = interf->front.specular_fraction; break; 225 case SDIS_BACK: f = interf->back.specular_fraction; break; 226 default: FATAL("Unreachable code.\n"); break; 227 } 228 return f; 229 } 230 231 static double 232 interface_get_reference_temperature 233 (const struct sdis_interface_fragment* frag, struct sdis_data* data) 234 { 235 const struct interfac* interf; 236 double T = SDIS_TEMPERATURE_NONE; 237 CHK(data != NULL && frag != NULL); 238 interf = sdis_data_cget(data); 239 switch(frag->side) { 240 case SDIS_FRONT: T = interf->front.reference_temperature; break; 241 case SDIS_BACK: T = interf->back.reference_temperature; break; 242 default: FATAL("Unreachable code.\n"); break; 243 } 244 return T; 245 } 246 247 /******************************************************************************* 248 * Helper functions 249 ******************************************************************************/ 250 static void 251 create_interface 252 (struct sdis_device* dev, 253 struct sdis_medium* front, 254 struct sdis_medium* back, 255 const struct interfac* interf, 256 struct sdis_interface** out_interf) 257 { 258 struct sdis_interface_shader shader = SDIS_INTERFACE_SHADER_NULL; 259 struct sdis_data* data = NULL; 260 const enum sdis_medium_type type_f = sdis_medium_get_type(front); 261 const enum sdis_medium_type type_b = sdis_medium_get_type(back); 262 263 CHK(interf != NULL); 264 265 shader.back.temperature = interface_get_temperature; 266 shader.front.temperature = interface_get_temperature; 267 268 if(type_f != type_b) { 269 shader.convection_coef = interface_get_convection_coef; 270 } 271 if(type_f == SDIS_FLUID) { 272 shader.front.emissivity = interface_get_emissivity; 273 shader.front.specular_fraction = interface_get_specular_fraction; 274 shader.front.reference_temperature = interface_get_reference_temperature; 275 } 276 if(type_b == SDIS_FLUID) { 277 shader.back.emissivity = interface_get_emissivity; 278 shader.back.specular_fraction = interface_get_specular_fraction; 279 shader.back.reference_temperature = interface_get_reference_temperature; 280 } 281 shader.convection_coef_upper_bound = MMAX(0, interf->convection_coef); 282 283 OK(sdis_data_create(dev, sizeof(struct interfac), ALIGNOF(struct interfac), 284 NULL, &data)); 285 *((struct interfac*)sdis_data_get(data)) = *interf; 286 287 OK(sdis_interface_create(dev, front, back, &shader, data, out_interf)); 288 OK(sdis_data_ref_put(data)); 289 } 290 291 /******************************************************************************* 292 * Test that the evaluation of the green function failed with a picard order 293 * greater than 1, i.e. when one want to handle the non-linearties of the 294 * system. 295 ******************************************************************************/ 296 static void 297 test_invalidity_picardN_green 298 (struct sdis_scene* scn, 299 struct sdis_medium* solid) 300 { 301 struct sdis_solve_probe_args probe = SDIS_SOLVE_PROBE_ARGS_DEFAULT; 302 struct sdis_solve_boundary_args bound = SDIS_SOLVE_BOUNDARY_ARGS_DEFAULT; 303 struct sdis_solve_medium_args mdm = SDIS_SOLVE_MEDIUM_ARGS_DEFAULT; 304 struct sdis_solve_probe_boundary_args probe_bound = 305 SDIS_SOLVE_PROBE_BOUNDARY_ARGS_DEFAULT; 306 307 struct sdis_green_function* green = NULL; 308 CHK(scn); 309 310 CHK(probe.picard_order == 1); 311 CHK(probe_bound.picard_order == 1); 312 CHK(bound.picard_order == 1); 313 CHK(mdm.picard_order == 1); 314 315 probe.position[0] = 0; 316 probe.position[1] = 0; 317 probe.picard_order = 2; 318 BA(sdis_solve_probe_green_function(scn, &probe, &green)); 319 320 probe_bound.iprim = 1; /* Solid left */ 321 probe_bound.uv[0] = 0.5; 322 probe_bound.side = SDIS_FRONT; 323 probe_bound.picard_order = 2; 324 BA(sdis_solve_probe_boundary_green_function(scn, &probe_bound, &green)); 325 326 bound.primitives = &probe_bound.iprim; 327 bound.sides = &probe_bound.side; 328 bound.nprimitives = 1; 329 bound.picard_order = 2; 330 BA(sdis_solve_boundary_green_function(scn, &bound, &green)); 331 332 mdm.medium = solid; 333 mdm.picard_order = 2; 334 BA(sdis_solve_medium_green_function(scn, &mdm, &green)); 335 } 336 337 /******************************************************************************* 338 * Test 339 ******************************************************************************/ 340 int 341 main(int argc, char** argv) 342 { 343 struct mem_allocator allocator; 344 struct interfac interf; 345 struct geometry geom; 346 struct ssp_rng* rng = NULL; 347 struct sdis_scene* scn = NULL; 348 struct sdis_data* data = NULL; 349 struct sdis_device* dev = NULL; 350 struct sdis_medium* fluid = NULL; 351 struct sdis_medium* solid = NULL; 352 struct sdis_medium* solid2 = NULL; 353 struct sdis_interface* interfaces[5] = {NULL}; 354 struct sdis_interface* prim_interfaces[10/*#segment*/]; 355 struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT; 356 struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER; 357 struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER; 358 const size_t nsimuls = 4; 359 size_t isimul; 360 const double emissivity = 1;/* Emissivity of the side +/-X of the solid */ 361 const double lambda = 0.1; /* Conductivity of the solid */ 362 const double Tref = 300; /* Reference temperature */ 363 const double T0 = 300; /* Fixed temperature on the left side of the system */ 364 const double T1 = 310; /* Fixed temperature on the right side of the system */ 365 const double thickness = 2.0; /* Thickness of the solid along X */ 366 double Ts0, Ts1, hr, tmp; 367 (void)argc, (void)argv; 368 369 OK(mem_init_proxy_allocator(&allocator, &mem_default_allocator)); 370 OK(sdis_device_create(&SDIS_DEVICE_CREATE_ARGS_DEFAULT, &dev)); 371 372 /* Create the fluid medium */ 373 fluid_shader.temperature = temperature_unknown; 374 OK(sdis_fluid_create(dev, &fluid_shader, NULL, &fluid)); 375 376 /* Create the solid medium */ 377 OK(sdis_data_create 378 (dev, sizeof(struct solid), ALIGNOF(struct solid), NULL, &data)); 379 ((struct solid*)sdis_data_get(data))->lambda = lambda; 380 solid_shader.calorific_capacity = solid_get_calorific_capacity; 381 solid_shader.thermal_conductivity = solid_get_thermal_conductivity; 382 solid_shader.volumic_mass = solid_get_volumic_mass; 383 solid_shader.delta = solid_get_delta; 384 solid_shader.temperature = temperature_unknown; 385 OK(sdis_solid_create(dev, &solid_shader, data, &solid)); 386 OK(sdis_data_ref_put(data)); 387 388 /* Create the surrounding solid medium */ 389 OK(sdis_data_create(dev, sizeof(struct solid), ALIGNOF(struct solid), 390 NULL, &data)); 391 ((struct solid*)sdis_data_get(data))->lambda = 0; 392 solid_shader.calorific_capacity = solid_get_thermal_conductivity; 393 solid_shader.thermal_conductivity = solid_get_thermal_conductivity; 394 solid_shader.volumic_mass = solid_get_volumic_mass; 395 solid_shader.delta = solid_get_delta; 396 OK(sdis_solid_create(dev, &solid_shader, data, &solid2)); 397 OK(sdis_data_ref_put(data)); 398 399 /* Create the interface that forces to keep in conduction */ 400 interf = INTERFACE_NULL; 401 create_interface(dev, solid, solid2, &interf, interfaces+0); 402 403 /* Create the interface that emits radiative heat from the solid */ 404 interf = INTERFACE_NULL; 405 interf.back.temperature = SDIS_TEMPERATURE_NONE; 406 interf.back.emissivity = emissivity; 407 interf.back.specular_fraction = -1; /* Should not be fetched */ 408 interf.back.reference_temperature = Tref; 409 create_interface(dev, solid, fluid, &interf, interfaces+1); 410 411 /* Create the interface that forces the radiative heat to bounce */ 412 interf = INTERFACE_NULL; 413 interf.front.temperature = SDIS_TEMPERATURE_NONE; 414 interf.front.emissivity = 0; 415 interf.front.specular_fraction = 1; 416 interf.front.reference_temperature = Tref; 417 create_interface(dev, fluid, solid2, &interf, interfaces+2); 418 419 /* Create the interface with a limit condition of T0 Kelvin */ 420 interf = INTERFACE_NULL; 421 interf.front.temperature = T0; 422 interf.front.emissivity = 1; 423 interf.front.specular_fraction = 1; 424 interf.front.reference_temperature = T0; 425 create_interface(dev, fluid, solid2, &interf, interfaces+3); 426 427 /* Create the interface with a limit condition of T1 Kelvin */ 428 interf = INTERFACE_NULL; 429 interf.front.temperature = T1; 430 interf.front.emissivity = 1; 431 interf.front.specular_fraction = 1; 432 interf.front.reference_temperature = T1; 433 create_interface(dev, fluid, solid2, &interf, interfaces+4); 434 435 /* Setup the per primitive interface of the solid medium */ 436 prim_interfaces[0] = interfaces[0]; 437 prim_interfaces[1] = interfaces[1]; 438 prim_interfaces[2] = interfaces[0]; 439 prim_interfaces[3] = interfaces[1]; 440 441 /* Setup the per primitive interface of the fluid on the left of the medium */ 442 prim_interfaces[4] = interfaces[2]; 443 prim_interfaces[5] = interfaces[3]; 444 prim_interfaces[6] = interfaces[2]; 445 446 /* Setup the per primitive interface of the fluid on the right of the medium */ 447 prim_interfaces[7] = interfaces[2]; 448 prim_interfaces[8] = interfaces[2]; 449 prim_interfaces[9] = interfaces[4]; 450 451 /* Create the scene */ 452 geom.positions = vertices; 453 geom.indices = indices; 454 geom.interfaces = prim_interfaces; 455 scn_args.get_indices = get_indices; 456 scn_args.get_interface = get_interface; 457 scn_args.get_position = get_position; 458 scn_args.nprimitives = nsegments; 459 scn_args.nvertices = nvertices; 460 scn_args.context = &geom; 461 scn_args.t_range[0] = MMIN(T0, T1); 462 scn_args.t_range[1] = MMAX(T0, T1); 463 OK(sdis_scene_2d_create(dev, &scn_args, &scn)); 464 465 hr = 4*BOLTZMANN_CONSTANT * Tref*Tref*Tref * emissivity; 466 tmp = lambda/(2*lambda + thickness*hr) * (T1 - T0); 467 Ts0 = T0 + tmp; 468 Ts1 = T1 - tmp; 469 470 /* Run the simulations */ 471 OK(ssp_rng_create(&allocator, SSP_RNG_KISS, &rng)); 472 FOR_EACH(isimul, 0, nsimuls) { 473 struct sdis_mc T = SDIS_MC_NULL; 474 struct sdis_mc time = SDIS_MC_NULL; 475 struct sdis_estimator* estimator; 476 struct sdis_estimator* estimator2; 477 struct sdis_green_function* green; 478 struct sdis_solve_probe_args solve_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT; 479 double ref, u; 480 size_t nreals = 0; 481 size_t nfails = 0; 482 const size_t N = 10000; 483 484 solve_args.nrealisations = N; 485 solve_args.position[0] = ssp_rng_uniform_double(rng, -0.9, 0.9); 486 solve_args.position[1] = ssp_rng_uniform_double(rng, -0.9, 0.9); 487 solve_args.time_range[0] = INF; 488 solve_args.time_range[1] = INF; 489 490 OK(sdis_solve_probe(scn, &solve_args, &estimator)); 491 OK(sdis_estimator_get_realisation_count(estimator, &nreals)); 492 OK(sdis_estimator_get_failure_count(estimator, &nfails)); 493 OK(sdis_estimator_get_temperature(estimator, &T)); 494 OK(sdis_estimator_get_realisation_time(estimator, &time)); 495 496 u = (solve_args.position[0] + 1) / thickness; 497 ref = u * Ts1 + (1-u) * Ts0; 498 printf("Temperature at (%g, %g) = %g ~ %g +/- %g\n", 499 SPLIT2(solve_args.position), ref, T.E, T.SE); 500 printf("Time per realisation (in usec) = %g +/- %g\n", time.E, time.SE); 501 printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); 502 503 CHK(nfails + nreals == N); 504 CHK(nfails < N/1000); 505 CHK(eq_eps(T.E, ref, 3*T.SE) == 1); 506 507 /* Check green function */ 508 OK(sdis_solve_probe_green_function(scn, &solve_args, &green)); 509 OK(sdis_green_function_solve(green, &estimator2)); 510 check_green_function(green); 511 check_estimator_eq(estimator, estimator2); 512 check_green_serialization(green, scn); 513 514 OK(sdis_estimator_ref_put(estimator)); 515 OK(sdis_estimator_ref_put(estimator2)); 516 OK(sdis_green_function_ref_put(green)); 517 518 solve_args.nrealisations = 10; 519 solve_args.register_paths = SDIS_HEAT_PATH_ALL; 520 521 OK(sdis_solve_probe(scn, &solve_args, &estimator)); 522 OK(sdis_estimator_ref_put(estimator)); 523 524 printf("\n\n"); 525 } 526 527 test_invalidity_picardN_green(scn, solid); 528 529 /* Release memory */ 530 OK(sdis_scene_ref_put(scn)); 531 OK(sdis_interface_ref_put(interfaces[0])); 532 OK(sdis_interface_ref_put(interfaces[1])); 533 OK(sdis_interface_ref_put(interfaces[2])); 534 OK(sdis_interface_ref_put(interfaces[3])); 535 OK(sdis_interface_ref_put(interfaces[4])); 536 OK(sdis_medium_ref_put(fluid)); 537 OK(sdis_medium_ref_put(solid)); 538 OK(sdis_medium_ref_put(solid2)); 539 OK(ssp_rng_ref_put(rng)); 540 OK(sdis_device_ref_put(dev)); 541 542 check_memory_allocator(&allocator); 543 mem_shutdown_proxy_allocator(&allocator); 544 CHK(mem_allocated_size() == 0); 545 546 return 0; 547 } 548