test_sdis_volumic_power2.c (15363B)
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 #define DELTA 0.01 23 #define DELTA_PSQUARE 0.01 24 25 struct reference { 26 double pos[3]; 27 double temperature_2d; /* In celcius */ 28 double temperature_3d; /* In celcius */ 29 }; 30 31 /* Temperature in Celcius. The reference is computed by EDF with Syrthes 32 * #realisations: 100000 33 * 34 * >>> Check 1 35 * 0.85 0 = 190.29 ~ 190.198 +/- 0.572596; #failures: 46 36 * 0.65 0 = 259.95 ~ 259.730 +/- 0.678251; #failures: 73 37 * 0.45 0 = 286.33 ~ 285.287 +/- 0.693572; #failures: 74 38 * 0.25 0 = 235.44 ~ 235.672 +/- 0.710927; #failures: 61 39 * 0.05 0 = 192.33 ~ 192.464 +/- 0.693148; #failures: 70 40 *-0.15 0 = 156.82 ~ 157.526 +/- 0.668902; #failures: 43 41 *-0.35 0 = 123.26 ~ 124.234 +/- 0.634061; #failures: 31 42 *-0.55 0 = 90.250 ~ 91.0285 +/- 0.566423; #failures: 32 43 * 44 * >>> Check 2 45 * 0.85 0 = 678.170 ~ 671.302 +/- 4.03424; #failures: 186 46 * 0.65 0 = 1520.84 ~ 1523.42 +/- 5.38182; #failures: 442 47 * 0.45 0 = 1794.57 ~ 1790.60 +/- 5.44808; #failures: 528 48 * 0.25 0 = 1429.74 ~ 1419.80 +/- 5.33467; #failures: 406 */ 49 50 static const double vertices[16/*#vertices*/*3/*#coords per vertex*/] = { 51 -0.5,-1.0,-0.5, 52 -0.5, 1.0,-0.5, 53 0.5, 1.0,-0.5, 54 0.5,-1.0,-0.5, 55 -0.5,-1.0, 0.5, 56 -0.5, 1.0, 0.5, 57 0.5, 1.0, 0.5, 58 0.5,-1.0, 0.5, 59 -0.1, 0.4,-0.5, 60 -0.1, 0.6,-0.5, 61 0.1, 0.6,-0.5, 62 0.1, 0.4,-0.5, 63 -0.1, 0.4, 0.5, 64 -0.1, 0.6, 0.5, 65 0.1, 0.6, 0.5, 66 0.1, 0.4, 0.5 67 }; 68 static const size_t nvertices = sizeof(vertices)/(sizeof(double)*3); 69 70 static const size_t indices[36/*#triangles*/*3/*#indices per triangle*/]= { 71 0, 4, 5, 5, 1, 0, /* Cuboid left */ 72 1, 5, 6, 6, 2, 1, /* Cuboid top */ 73 6, 7, 3, 3, 2, 6, /* Cuboid right */ 74 0, 3, 7, 7, 4, 0, /* Cuboid bottom */ 75 /* Cuboid back */ 76 0, 1, 9, 9, 8, 0, 77 1, 2, 10, 10, 9, 1, 78 2, 3, 11, 11, 10, 2, 79 3, 0, 8, 8, 11, 3, 80 /* Cuboid front */ 81 5, 4, 12, 12, 13, 5, 82 5, 13, 14, 14, 6, 5, 83 6, 14, 15, 15, 7, 6, 84 7, 15, 12, 12, 4, 7, 85 8, 12, 13, 13, 9, 8, /* Cube left */ 86 9, 13, 14, 14, 10, 9, /* Cube top */ 87 14, 15, 11, 11, 10, 14, /* Cube right */ 88 8, 11, 15, 15, 12, 8, /* Cube bottom */ 89 8, 9, 10, 10, 11, 8, /* Cube back */ 90 12, 15, 14, 14, 13, 12 /* Cube front */ 91 }; 92 static const size_t ntriangles = sizeof(indices)/(sizeof(size_t)*3); 93 94 /******************************************************************************* 95 * Geometry 96 ******************************************************************************/ 97 static void 98 get_indices(const size_t itri, size_t ids[3], void* context) 99 { 100 (void)context; 101 CHK(ids); 102 ids[0] = indices[itri*3+0]; 103 ids[1] = indices[itri*3+1]; 104 ids[2] = indices[itri*3+2]; 105 } 106 107 static void 108 get_position(const size_t ivert, double pos[3], void* context) 109 { 110 (void)context; 111 CHK(pos); 112 pos[0] = vertices[ivert*3+0]; 113 pos[1] = vertices[ivert*3+1]; 114 pos[2] = vertices[ivert*3+2]; 115 } 116 117 static void 118 get_interface(const size_t itri, struct sdis_interface** bound, void* context) 119 { 120 struct sdis_interface** interfaces = context; 121 CHK(context && bound); 122 *bound = interfaces[itri/2]; 123 } 124 125 /******************************************************************************* 126 * Solid medium 127 ******************************************************************************/ 128 struct solid { 129 double cp; 130 double lambda; 131 double rho; 132 double delta; 133 double P; 134 double T; 135 }; 136 137 static double 138 solid_get_calorific_capacity 139 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 140 { 141 CHK(data != NULL && vtx != NULL); 142 return ((const struct solid*)sdis_data_cget(data))->cp; 143 } 144 145 static double 146 solid_get_thermal_conductivity 147 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 148 { 149 CHK(data != NULL && vtx != NULL); 150 return ((const struct solid*)sdis_data_cget(data))->lambda; 151 } 152 153 static double 154 solid_get_volumic_mass 155 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 156 { 157 CHK(data != NULL && vtx != NULL); 158 return ((const struct solid*)sdis_data_cget(data))->rho; 159 } 160 161 static double 162 solid_get_delta 163 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 164 { 165 CHK(data != NULL && vtx != NULL); 166 return ((const struct solid*)sdis_data_cget(data))->delta; 167 } 168 169 static double 170 solid_get_temperature 171 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 172 { 173 CHK(data != NULL && vtx != NULL); 174 return ((const struct solid*)sdis_data_cget(data))->T; 175 } 176 177 static double 178 solid_get_volumic_power 179 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 180 { 181 CHK(data != NULL && vtx != NULL); 182 return ((const struct solid*)sdis_data_cget(data))->P; 183 } 184 185 /******************************************************************************* 186 * Fluid medium 187 ******************************************************************************/ 188 struct fluid { 189 double temperature; 190 }; 191 192 static double 193 fluid_get_temperature 194 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 195 { 196 CHK(data != NULL && vtx != NULL); 197 return ((const struct fluid*)sdis_data_cget(data))->temperature; 198 } 199 200 /******************************************************************************* 201 * Interfaces 202 ******************************************************************************/ 203 struct interf { 204 double h; 205 double temperature; 206 }; 207 208 static double 209 interface_get_convection_coef 210 (const struct sdis_interface_fragment* frag, struct sdis_data* data) 211 { 212 CHK(frag && data); 213 return ((const struct interf*)sdis_data_cget(data))->h; 214 } 215 216 static double 217 interface_get_temperature 218 (const struct sdis_interface_fragment* frag, struct sdis_data* data) 219 { 220 CHK(frag && data); 221 return ((const struct interf*)sdis_data_cget(data))->temperature; 222 } 223 224 /******************************************************************************* 225 * Helper functions 226 ******************************************************************************/ 227 static void 228 check(struct sdis_scene* scn, const struct reference refs[], const size_t nrefs) 229 { 230 struct sdis_estimator* estimator = NULL; 231 struct sdis_solve_probe_args solve_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT; 232 struct sdis_mc T = SDIS_MC_NULL; 233 size_t nreals; 234 size_t nfails; 235 size_t i; 236 237 solve_args.time_range[0] = INF; 238 solve_args.time_range[1] = INF; 239 solve_args.nrealisations = N; 240 241 FOR_EACH(i, 0, nrefs) { 242 double Tc; 243 solve_args.position[0] = refs[i].pos[0]; 244 solve_args.position[1] = refs[i].pos[1]; 245 solve_args.position[2] = refs[i].pos[2]; 246 247 OK(sdis_solve_probe(scn, &solve_args, &estimator)); 248 OK(sdis_estimator_get_temperature(estimator, &T)); 249 OK(sdis_estimator_get_realisation_count(estimator, &nreals)); 250 OK(sdis_estimator_get_failure_count(estimator, &nfails)); 251 Tc = T.E - 273.15; /* Convert in Celcius */ 252 printf("Temperature at (%g %g %g) = %g ~ %g +/- %g [%g, %g]\n", 253 SPLIT3(refs[i].pos), refs[i].temperature_2d, Tc, T.SE, Tc-3*T.SE, Tc+3*T.SE); 254 printf("#realisations: %lu; #failures: %lu\n", 255 (unsigned long)nreals, (unsigned long)nfails); 256 /*CHK(eq_eps(Tc, refs[i].temperature, T.SE*3));*/ 257 OK(sdis_estimator_ref_put(estimator)); 258 } 259 } 260 261 int 262 main(int argc, char** argv) 263 { 264 struct solid* solid_param = NULL; 265 struct fluid* fluid_param = NULL; 266 struct interf* interf_param = NULL; 267 struct sdis_device* dev = NULL; 268 struct sdis_data* data = NULL; 269 struct sdis_medium* fluid1 = NULL; 270 struct sdis_medium* fluid2 = NULL; 271 struct sdis_medium* solid1 = NULL; 272 struct sdis_medium* solid2 = NULL; 273 struct sdis_scene* scn = NULL; 274 struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT; 275 struct sdis_fluid_shader fluid_shader = SDIS_FLUID_SHADER_NULL; 276 struct sdis_solid_shader solid_shader = SDIS_SOLID_SHADER_NULL; 277 struct sdis_interface_shader interf_shader = SDIS_INTERFACE_SHADER_NULL; 278 struct sdis_interface* interf_solid1_adiabatic = NULL; 279 struct sdis_interface* interf_solid2_adiabatic = NULL; 280 struct sdis_interface* interf_solid1_solid2 = NULL; 281 struct sdis_interface* interf_solid1_fluid1 = NULL; 282 struct sdis_interface* interf_solid1_fluid2 = NULL; 283 struct sdis_interface* interfaces[18 /*#rectangles*/]; 284 /* In celcius. Computed by EDF with Syrthes */ 285 const struct reference refs1[] = { /* Lambda1=1, Lambda2=10, Pw = 10000 */ 286 {{0, 0.85, 0}, 190.29, 189.13}, 287 {{0, 0.65, 0}, 259.95, 247.09}, 288 {{0, 0.45, 0}, 286.33, 308.42}, 289 {{0, 0.25, 0}, 235.44, 233.55}, 290 {{0, 0.05, 0}, 192.33, 192.30}, 291 {{0,-0.15, 0}, 156.82, 156.98}, 292 {{0,-0.35, 0}, 123.26, 123.43}, 293 {{0,-0.55, 0}, 90.250, 90.040} 294 }; 295 const struct reference refs2[] = { /* Lambda1=0.1, Lambda2=10, Pw=10000 */ 296 {{0, 0.85}, 678.170, 0}, 297 {{0, 0.65}, 1520.84, 0}, 298 {{0, 0.45}, 1794.57, 0}, 299 {{0, 0.25}, 1429.74, 0} 300 }; 301 (void)argc, (void)argv; 302 303 OK(sdis_device_create(&SDIS_DEVICE_CREATE_ARGS_DEFAULT, &dev)); 304 305 /* Setup the fluid shader */ 306 fluid_shader.temperature = fluid_get_temperature; 307 fluid_shader.calorific_capacity = dummy_medium_getter; 308 fluid_shader.volumic_mass = dummy_medium_getter; 309 310 /* Create the fluid1 medium */ 311 OK(sdis_data_create 312 (dev, sizeof(struct fluid), ALIGNOF(struct fluid), NULL, &data)); 313 fluid_param = sdis_data_get(data); 314 fluid_param->temperature = 373.15; 315 OK(sdis_fluid_create(dev, &fluid_shader, data, &fluid1) ); 316 OK(sdis_data_ref_put(data)); 317 318 /* Create the fluid2 medium */ 319 OK(sdis_data_create 320 (dev, sizeof(struct fluid), ALIGNOF(struct fluid), NULL, &data)); 321 fluid_param = sdis_data_get(data); 322 fluid_param->temperature = 273.15; 323 OK(sdis_fluid_create(dev, &fluid_shader, data, &fluid2)); 324 OK(sdis_data_ref_put(data)); 325 326 /* Setup the solid shader */ 327 solid_shader.calorific_capacity = solid_get_calorific_capacity; 328 solid_shader.thermal_conductivity = solid_get_thermal_conductivity; 329 solid_shader.volumic_mass = solid_get_volumic_mass; 330 solid_shader.delta = solid_get_delta; 331 solid_shader.temperature = solid_get_temperature; 332 solid_shader.volumic_power = solid_get_volumic_power; 333 334 /* Create the solid1 medium */ 335 OK(sdis_data_create 336 (dev, sizeof(struct solid), ALIGNOF(struct solid), NULL, &data) ); 337 solid_param = sdis_data_get(data); 338 solid_param->cp = 500000; 339 solid_param->rho = 1000; 340 solid_param->lambda = 1; 341 solid_param->delta = DELTA; 342 solid_param->P = SDIS_VOLUMIC_POWER_NONE; 343 solid_param->T = SDIS_TEMPERATURE_NONE; 344 OK(sdis_solid_create(dev, &solid_shader, data, &solid1)); 345 OK(sdis_data_ref_put(data)); 346 347 /* Create the solid2 medium */ 348 OK(sdis_data_create 349 (dev, sizeof(struct solid), ALIGNOF(struct solid), NULL, &data)); 350 solid_param = sdis_data_get(data); 351 solid_param->cp = 500000; 352 solid_param->rho = 1000; 353 solid_param->lambda = 10; 354 solid_param->delta = DELTA_PSQUARE; 355 solid_param->P = Pw; 356 solid_param->T = SDIS_TEMPERATURE_NONE; 357 OK(sdis_solid_create(dev, &solid_shader, data, &solid2)); 358 OK(sdis_data_ref_put(data)); 359 360 /* Create the solid1/solid2 interface */ 361 OK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf), 362 NULL, &data)); 363 OK(sdis_interface_create(dev, solid2, solid1, &SDIS_INTERFACE_SHADER_NULL, 364 NULL, &interf_solid1_solid2)); 365 OK(sdis_data_ref_put(data)); 366 367 /* Setup the interface shader */ 368 interf_shader.convection_coef = interface_get_convection_coef; 369 370 /* Create the adiabatic interfaces */ 371 OK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf), 372 NULL, &data)); 373 interf_param = sdis_data_get(data); 374 interf_param->h = 0; 375 OK(sdis_interface_create(dev, solid1, fluid1, &interf_shader, data, 376 &interf_solid1_adiabatic)); 377 OK(sdis_interface_create(dev, solid2, fluid1, &interf_shader, data, 378 &interf_solid2_adiabatic)); 379 OK(sdis_data_ref_put(data)); 380 381 /* Setup the interface shader */ 382 interf_shader.front.temperature = interface_get_temperature; 383 384 /* Create the solid1/fluid1 interface */ 385 OK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf), 386 NULL, &data)); 387 interf_param = sdis_data_get(data); 388 interf_param->h = 5; 389 interf_param->temperature = SDIS_TEMPERATURE_NONE; 390 OK(sdis_interface_create(dev, solid1, fluid1, &interf_shader, data, 391 &interf_solid1_fluid1)); 392 OK(sdis_data_ref_put(data)); 393 394 /* Create the solid1/fluid2 interace */ 395 OK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf), 396 NULL, &data)); 397 interf_param = sdis_data_get(data); 398 interf_param->h = 10; 399 interf_param->temperature = SDIS_TEMPERATURE_NONE; 400 OK(sdis_interface_create(dev, solid1, fluid2, &interf_shader, data, 401 &interf_solid1_fluid2)); 402 OK(sdis_data_ref_put(data)); 403 404 /* Map the interfaces to their faces */ 405 interfaces[0] = interf_solid1_adiabatic; 406 interfaces[1] = interf_solid1_fluid1; 407 interfaces[2] = interf_solid1_adiabatic; 408 interfaces[3] = interf_solid1_fluid2; 409 interfaces[4] = interf_solid1_adiabatic; 410 interfaces[5] = interf_solid1_adiabatic; 411 interfaces[6] = interf_solid1_adiabatic; 412 interfaces[7] = interf_solid1_adiabatic; 413 interfaces[8] = interf_solid1_adiabatic; 414 interfaces[9] = interf_solid1_adiabatic; 415 interfaces[10] = interf_solid1_adiabatic; 416 interfaces[11] = interf_solid1_adiabatic; 417 interfaces[12] = interf_solid1_solid2; 418 interfaces[13] = interf_solid1_solid2; 419 interfaces[14] = interf_solid1_solid2; 420 interfaces[15] = interf_solid1_solid2; 421 interfaces[16] = interf_solid2_adiabatic; 422 interfaces[17] = interf_solid2_adiabatic; 423 424 /* Create the scene */ 425 scn_args.get_indices = get_indices; 426 scn_args.get_interface = get_interface; 427 scn_args.get_position = get_position; 428 scn_args.nprimitives = ntriangles; 429 scn_args.nvertices = nvertices; 430 scn_args.context = interfaces; 431 OK(sdis_scene_create(dev, &scn_args, &scn)); 432 433 #if 0 434 dump_mesh(stdout, vertices, nvertices, indices, ntriangles); 435 exit(0); 436 #endif 437 438 printf(">>> Check 1\n"); 439 check(scn, refs1, sizeof(refs1)/sizeof(struct reference)); 440 441 /* Update the scene */ 442 OK(sdis_scene_ref_put(scn)); 443 data = sdis_medium_get_data(solid1); 444 solid_param = sdis_data_get(data); 445 solid_param->lambda = 0.1; 446 OK(sdis_scene_create(dev, &scn_args, &scn)); 447 448 printf("\n>>> Check 2\n"); 449 check(scn, refs2, sizeof(refs2)/sizeof(struct reference)); 450 451 /* Release the interfaces */ 452 OK(sdis_interface_ref_put(interf_solid1_adiabatic)); 453 OK(sdis_interface_ref_put(interf_solid2_adiabatic)); 454 OK(sdis_interface_ref_put(interf_solid1_fluid1)); 455 OK(sdis_interface_ref_put(interf_solid1_fluid2)); 456 OK(sdis_interface_ref_put(interf_solid1_solid2)); 457 458 /* Release the media */ 459 OK(sdis_medium_ref_put(fluid1)); 460 OK(sdis_medium_ref_put(fluid2)); 461 OK(sdis_medium_ref_put(solid1)); 462 OK(sdis_medium_ref_put(solid2)); 463 464 OK(sdis_scene_ref_put(scn)); 465 OK(sdis_device_ref_put(dev)); 466 467 CHK(mem_allocated_size() == 0); 468 return 0; 469 } 470