test_sdis_volumic_power2_2d.c (15957B)
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 #include <rsys/math.h> 19 20 #define N 10000 /* #realisations */ 21 #define Pw 10000 /* Volumic power */ 22 23 /* H delta T */ 24 #define Tboundary1 SDIS_TEMPERATURE_NONE 25 #define Tboundary2 SDIS_TEMPERATURE_NONE 26 #define DELTA 0.01 27 #define Tref 286.83 /* In Celsius. Computed with Syrthes at the position 0.5 */ 28 29 /* Dirichlets */ 30 /*#define Tboundary1 373.15*/ 31 /*#define Tboundary2 273.15*/ 32 /*#define DELTA 0.01*/ 33 /*#define Tref 246.93*/ /* In Celsius. Computed with Syrthes at the position 0.5 */ 34 35 /* Temperature in Celcius. The reference is computed by EDF with Syrthes 36 * #realisations: 100000 37 * 38 * >>> Check1 39 * 0.85 = 190.29 ~ 189.322 +/- 0.566717; #failures: 51 40 * 0.65 = 259.95 ~ 259.995 +/- 0.674453; #failures: 82 41 * 0.45 = 286.33 ~ 285.928 +/- 0.691044; #failures: 76 42 * 0.25 = 235.44 ~ 234.672 +/- 0.700354; #failures: 80 43 * 0.05 = 192.33 ~ 191.977 +/- 0.690793; #failures: 64 44 *-0.15 = 156.82 ~ 155.765 +/- 0.660722; #failures: 40 45 *-0.35 = 123.26 ~ 122.973 +/- 0.621093; #failures: 29 46 *-0.55 = 90.250 ~ 90.3501 +/- 0.561255; #failures: 27 47 * 48 * >>> Check 2 49 * 0.85 = 678.170 ~ 662.616 +/- 3.97997; #failures: 221 50 * 0.65 = 1520.84 ~ 1486.35 +/- 5.25785; #failures: 474 51 * 0.45 = 1794.57 ~ 1767.21 +/- 5.36318; #failures: 584 52 * 0.25 = 1429.74 ~ 1401.39 +/- 5.25579; #failures: 465 53 * 54 * >>> Check 3 55 * 0.85 = 83.99 ~ 84.0098 +/- 0.115932; #failures: 51 56 * 0.65 = 73.90 ~ 73.9596 +/- 0.138835; #failures: 82 57 * 0.45 = 68.43 ~ 70.0292 +/- 0.144928; #failures: 76 58 * 0.25 = 60.61 ~ 61.4412 +/- 0.153980; #failures: 80 59 * 0.05 = 52.09 ~ 51.9452 +/- 0.158045; #failures: 64 60 *-0.15 = 42.75 ~ 42.9072 +/- 0.156546; #failures: 40 61 *-0.35 = 33.04 ~ 33.9338 +/- 0.149751; #failures: 29 62 *-0.55 = 24.58 ~ 24.7237 +/- 0.136441; #failures: 27 */ 63 64 /* 65 * _\ T1 66 * / / 67 * \__/ 68 * ///+-----H1-------+/// 69 * ///| |/// 70 * ///| +------+ |/// 71 * ///| |LAMBDA| |/// 72 * ///| | Pw | |/// 73 * ///| +------+ |/// 74 * ///| |/// 75 * ///| |/// 76 * ///| LAMBDA1 |/// 77 * ///| |/// 78 * ///| |/// 79 * ///| |/// 80 * ///+-----H2-------+/// 81 * _\ T2 82 * / / 83 * \__/ 84 */ 85 86 struct reference { 87 double pos[2]; 88 double temperature; /* In celcius */ 89 }; 90 91 static const double vertices[8/*#vertices*/*2/*#coords per vertex*/] = { 92 -0.5,-1.0, 93 -0.5, 1.0, 94 0.5, 1.0, 95 0.5,-1.0, 96 -0.1, 0.4, 97 -0.1, 0.6, 98 0.1, 0.6, 99 0.1, 0.4 100 }; 101 static const size_t nvertices = sizeof(vertices)/(sizeof(double)*2); 102 103 static const size_t indices[8/*#segments*/*2/*#indices per segment*/]= { 104 0, 1, /* Rectangle left */ 105 1, 2, /* Rectangle top */ 106 2, 3, /* Rectangle right */ 107 3, 0, /* Rectangle bottom */ 108 4, 5, /* Square left */ 109 5, 6, /* Square top */ 110 6, 7, /* Square right */ 111 7, 4 /* Square bottom */ 112 }; 113 static const size_t nsegments = sizeof(indices)/(sizeof(size_t)*2); 114 115 /******************************************************************************* 116 * Geometry 117 ******************************************************************************/ 118 static void 119 get_indices(const size_t iseg, size_t ids[2], void* context) 120 { 121 (void)context; 122 CHK(ids); 123 ids[0] = indices[iseg*2+0]; 124 ids[1] = indices[iseg*2+1]; 125 } 126 127 static void 128 get_position(const size_t ivert, double pos[2], void* context) 129 { 130 (void)context; 131 CHK(pos); 132 pos[0] = vertices[ivert*2+0]; 133 pos[1] = vertices[ivert*2+1]; 134 } 135 136 static void 137 get_interface(const size_t iseg, struct sdis_interface** bound, void* context) 138 { 139 struct sdis_interface** interfaces = context; 140 CHK(context && bound); 141 *bound = interfaces[iseg]; 142 } 143 144 /******************************************************************************* 145 * Solid medium 146 ******************************************************************************/ 147 struct solid { 148 double cp; 149 double lambda; 150 double rho; 151 double delta; 152 double P; 153 double T; 154 }; 155 156 static double 157 solid_get_calorific_capacity 158 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 159 { 160 CHK(data != NULL && vtx != NULL); 161 return ((const struct solid*)sdis_data_cget(data))->cp; 162 } 163 164 static double 165 solid_get_thermal_conductivity 166 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 167 { 168 CHK(data != NULL && vtx != NULL); 169 return ((const struct solid*)sdis_data_cget(data))->lambda; 170 } 171 172 static double 173 solid_get_volumic_mass 174 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 175 { 176 CHK(data != NULL && vtx != NULL); 177 return ((const struct solid*)sdis_data_cget(data))->rho; 178 } 179 180 static double 181 solid_get_delta 182 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 183 { 184 CHK(data != NULL && vtx != NULL); 185 return ((const struct solid*)sdis_data_cget(data))->delta; 186 } 187 188 static double 189 solid_get_temperature 190 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 191 { 192 CHK(data != NULL && vtx != NULL); 193 return ((const struct solid*)sdis_data_cget(data))->T; 194 } 195 196 static double 197 solid_get_volumic_power 198 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 199 { 200 CHK(data != NULL && vtx != NULL); 201 return ((const struct solid*)sdis_data_cget(data))->P; 202 } 203 204 /******************************************************************************* 205 * Fluid medium 206 ******************************************************************************/ 207 struct fluid { 208 double temperature; 209 }; 210 211 static double 212 fluid_get_temperature 213 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 214 { 215 CHK(data != NULL && vtx != NULL); 216 return ((const struct fluid*)sdis_data_cget(data))->temperature; 217 } 218 219 220 /******************************************************************************* 221 * Interfaces 222 ******************************************************************************/ 223 struct interf { 224 double h; 225 double temperature; 226 }; 227 228 static double 229 interface_get_convection_coef 230 (const struct sdis_interface_fragment* frag, struct sdis_data* data) 231 { 232 CHK(frag && data); 233 return ((const struct interf*)sdis_data_cget(data))->h; 234 } 235 236 static double 237 interface_get_temperature 238 (const struct sdis_interface_fragment* frag, struct sdis_data* data) 239 { 240 CHK(frag && data); 241 return ((const struct interf*)sdis_data_cget(data))->temperature; 242 } 243 244 /******************************************************************************* 245 * Helper functions 246 ******************************************************************************/ 247 static void 248 check(struct sdis_scene* scn, const struct reference refs[], const size_t nrefs) 249 { 250 struct sdis_estimator* estimator = NULL; 251 struct sdis_mc T = SDIS_MC_NULL; 252 struct sdis_solve_probe_args solve_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT; 253 size_t nreals; 254 size_t nfails; 255 size_t i; 256 257 solve_args.nrealisations = N; 258 solve_args.time_range[0] = INF; 259 solve_args.time_range[1] = INF; 260 261 FOR_EACH(i, 0, nrefs) { 262 double Tc; 263 solve_args.position[0] = refs[i].pos[0]; 264 solve_args.position[1] = refs[i].pos[1]; 265 266 OK(sdis_solve_probe(scn, &solve_args, &estimator)); 267 OK(sdis_estimator_get_temperature(estimator, &T)); 268 OK(sdis_estimator_get_realisation_count(estimator, &nreals)); 269 OK(sdis_estimator_get_failure_count(estimator, &nfails)); 270 Tc = T.E - 273.15; /* Convert in Celcius */ 271 printf("Temperature at (%g %g) = %g ~ %g +/- %g [%g, %g]\n", 272 SPLIT2(refs[i].pos), refs[i].temperature, Tc, T.SE, Tc-3*T.SE, Tc+3*T.SE); 273 printf("#realisations: %lu; #failures: %lu\n", 274 (unsigned long)nreals, (unsigned long)nfails); 275 /*CHK(eq_eps(Tc, refs[i].temperature, T.SE*3));*/ 276 OK(sdis_estimator_ref_put(estimator)); 277 } 278 } 279 280 /******************************************************************************* 281 * Test 282 ******************************************************************************/ 283 int 284 main(int argc, char** argv) 285 { 286 struct solid* solid_param = NULL; 287 struct fluid* fluid_param = NULL; 288 struct interf* interf_param = NULL; 289 struct sdis_device* dev = NULL; 290 struct sdis_data* data = NULL; 291 struct sdis_medium* fluid1 = NULL; 292 struct sdis_medium* fluid2 = NULL; 293 struct sdis_medium* solid1 = NULL; 294 struct sdis_medium* solid2 = NULL; 295 struct sdis_scene* scn = NULL; 296 struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT; 297 struct sdis_fluid_shader fluid_shader = SDIS_FLUID_SHADER_NULL; 298 struct sdis_solid_shader solid_shader = SDIS_SOLID_SHADER_NULL; 299 struct sdis_interface_shader interf_shader = SDIS_INTERFACE_SHADER_NULL; 300 struct sdis_interface* interf_adiabatic = NULL; 301 struct sdis_interface* interf_solid1_solid2 = NULL; 302 struct sdis_interface* interf_solid1_fluid1 = NULL; 303 struct sdis_interface* interf_solid1_fluid2 = NULL; 304 struct sdis_interface* interfaces[8 /*#segment*/]; 305 306 /* In celcius. Computed by EDF with Syrthes */ 307 const struct reference refs1[] = { /* Lambda1=1, Lambda2=10, Pw = 10000 */ 308 {{0, 0.85}, 190.29}, 309 {{0, 0.65}, 259.95}, 310 {{0, 0.45}, 286.33}, 311 {{0, 0.25}, 235.44}, 312 {{0, 0.05}, 192.33}, 313 {{0,-0.15}, 156.82}, 314 {{0,-0.35}, 123.26}, 315 {{0,-0.55}, 90.250} 316 }; 317 const struct reference refs2[] = { /* Lambda1=0.1, Lambda2=10, Pw=10000 */ 318 {{0, 0.85}, 678.17}, 319 {{0, 0.65}, 1520.84}, 320 {{0, 0.45}, 1794.57}, 321 {{0, 0.25}, 1429.74} 322 }; 323 const struct reference refs3[] = { /* Lambda1=1, Lambda2=10, Pw=NONE */ 324 {{0, 0.85}, 83.99}, 325 {{0, 0.65}, 73.90}, 326 {{0, 0.45}, 68.43}, 327 {{0, 0.25}, 60.61}, 328 {{0, 0.05}, 52.09}, 329 {{0,-0.15}, 42.75}, 330 {{0,-0.35}, 33.04}, 331 {{0,-0.55}, 24.58} 332 }; 333 (void)argc, (void)argv; 334 335 OK(sdis_device_create(&SDIS_DEVICE_CREATE_ARGS_DEFAULT, &dev)); 336 337 /* Setup the fluid shader */ 338 fluid_shader.temperature = fluid_get_temperature; 339 fluid_shader.calorific_capacity = dummy_medium_getter; 340 fluid_shader.volumic_mass = dummy_medium_getter; 341 342 /* Create the fluid1 medium */ 343 OK(sdis_data_create 344 (dev, sizeof(struct fluid), ALIGNOF(struct fluid), NULL, &data)); 345 fluid_param = sdis_data_get(data); 346 fluid_param->temperature = 373.15; 347 OK(sdis_fluid_create(dev, &fluid_shader, data, &fluid1)); 348 OK(sdis_data_ref_put(data)); 349 350 /* Create the fluid2 medium */ 351 OK(sdis_data_create 352 (dev, sizeof(struct fluid), ALIGNOF(struct fluid), NULL, &data)); 353 fluid_param = sdis_data_get(data); 354 fluid_param->temperature = 273.15; 355 OK(sdis_fluid_create(dev, &fluid_shader, data, &fluid2)); 356 OK(sdis_data_ref_put(data)); 357 358 /* Setup the solid shader */ 359 solid_shader.calorific_capacity = solid_get_calorific_capacity; 360 solid_shader.thermal_conductivity = solid_get_thermal_conductivity; 361 solid_shader.volumic_mass = solid_get_volumic_mass; 362 solid_shader.delta = solid_get_delta; 363 solid_shader.temperature = solid_get_temperature; 364 solid_shader.volumic_power = solid_get_volumic_power; 365 366 /* Create the solid1 medium */ 367 OK(sdis_data_create 368 (dev, sizeof(struct solid), ALIGNOF(struct solid), NULL, &data)); 369 solid_param = sdis_data_get(data); 370 solid_param->cp = 500000; 371 solid_param->rho = 1000; 372 solid_param->lambda = 1; 373 solid_param->delta = DELTA; 374 solid_param->P = SDIS_VOLUMIC_POWER_NONE; 375 solid_param->T = SDIS_TEMPERATURE_NONE; 376 OK(sdis_solid_create(dev, &solid_shader, data, &solid1)); 377 OK(sdis_data_ref_put(data)); 378 379 /* Create the solid2 medium */ 380 OK(sdis_data_create 381 (dev, sizeof(struct solid), ALIGNOF(struct solid), NULL, &data)); 382 solid_param = sdis_data_get(data); 383 solid_param->cp = 500000; 384 solid_param->rho = 1000; 385 solid_param->lambda = 10; 386 solid_param->delta = DELTA; 387 solid_param->P = Pw; 388 solid_param->T = SDIS_TEMPERATURE_NONE; 389 OK(sdis_solid_create(dev, &solid_shader, data, &solid2)); 390 OK(sdis_data_ref_put(data)); 391 392 /* Create the solid1/solid2 interface */ 393 OK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf), 394 NULL, &data)); 395 OK(sdis_interface_create(dev, solid2, solid1, &SDIS_INTERFACE_SHADER_NULL, 396 NULL, &interf_solid1_solid2)); 397 OK(sdis_data_ref_put(data)); 398 399 /* Setup the interface shader */ 400 interf_shader.convection_coef = interface_get_convection_coef; 401 402 /* Create the adiabatic interface */ 403 OK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf), 404 NULL, &data)); 405 interf_param = sdis_data_get(data); 406 interf_param->h = 0; 407 OK(sdis_interface_create(dev, solid1, fluid1, &interf_shader, data, 408 &interf_adiabatic)); 409 OK(sdis_data_ref_put(data)); 410 411 /* Setup the interface shader */ 412 interf_shader.front.temperature = interface_get_temperature; 413 414 /* Create the solid1/fluid1 interface */ 415 OK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf), 416 NULL, &data)); 417 interf_param = sdis_data_get(data); 418 interf_param->h = 5; 419 interf_param->temperature = Tboundary1; 420 OK(sdis_interface_create(dev, solid1, fluid1, &interf_shader, data, 421 &interf_solid1_fluid1)); 422 OK(sdis_data_ref_put(data)); 423 424 /* Create the solid1/fluid2 interace */ 425 OK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf), 426 NULL, &data)); 427 interf_param = sdis_data_get(data); 428 interf_param->h = 10; 429 interf_param->temperature = Tboundary2; 430 OK(sdis_interface_create(dev, solid1, fluid2, &interf_shader, data, 431 &interf_solid1_fluid2)); 432 OK(sdis_data_ref_put(data)); 433 434 435 /* Map the interfaces to their square segments */ 436 interfaces[0] = interf_adiabatic; 437 interfaces[1] = interf_solid1_fluid1; 438 interfaces[2] = interf_adiabatic; 439 interfaces[3] = interf_solid1_fluid2; 440 interfaces[4] = interf_solid1_solid2; 441 interfaces[5] = interf_solid1_solid2; 442 interfaces[6] = interf_solid1_solid2; 443 interfaces[7] = interf_solid1_solid2; 444 445 /* Create the scene */ 446 scn_args.get_indices = get_indices; 447 scn_args.get_interface = get_interface; 448 scn_args.get_position = get_position; 449 scn_args.nprimitives = nsegments; 450 scn_args.nvertices = nvertices; 451 scn_args.context = interfaces; 452 OK(sdis_scene_2d_create(dev, &scn_args, &scn)); 453 454 printf(">>> Check 1\n"); 455 check(scn, refs1, sizeof(refs1)/sizeof(struct reference)); 456 457 /* Update the scene */ 458 OK(sdis_scene_ref_put(scn)); 459 data = sdis_medium_get_data(solid1); 460 solid_param = sdis_data_get(data); 461 solid_param->lambda = 0.1; 462 OK(sdis_scene_2d_create(dev, &scn_args, &scn)); 463 464 printf("\n>>> Check 2\n"); 465 check(scn, refs2, sizeof(refs2)/sizeof(struct reference)); 466 467 /* Update the scene */ 468 OK(sdis_scene_ref_put(scn)); 469 data = sdis_medium_get_data(solid1); 470 solid_param = sdis_data_get(data); 471 solid_param->lambda = 1; 472 data = sdis_medium_get_data(solid2); 473 solid_param = sdis_data_get(data); 474 solid_param->lambda = 10; 475 solid_param->P = SDIS_VOLUMIC_POWER_NONE; 476 OK(sdis_scene_2d_create(dev, &scn_args, &scn)); 477 478 printf("\n>>> Check 3\n"); 479 check(scn, refs3, sizeof(refs3)/sizeof(struct reference)); 480 481 #if 0 482 dump_segments(stdout, vertices, nvertices, indices, nsegments); 483 exit(0); 484 #endif 485 486 /* Release the interfaces */ 487 OK(sdis_interface_ref_put(interf_adiabatic)); 488 OK(sdis_interface_ref_put(interf_solid1_fluid1)); 489 OK(sdis_interface_ref_put(interf_solid1_fluid2)); 490 OK(sdis_interface_ref_put(interf_solid1_solid2)); 491 492 /* Release the media */ 493 OK(sdis_medium_ref_put(fluid1)); 494 OK(sdis_medium_ref_put(fluid2)); 495 OK(sdis_medium_ref_put(solid1)); 496 OK(sdis_medium_ref_put(solid2)); 497 498 OK(sdis_scene_ref_put(scn)); 499 OK(sdis_device_ref_put(dev)); 500 501 CHK(mem_allocated_size() == 0); 502 return 0; 503 } 504