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