test_sdis_volumic_power3_2d.c (15554B)
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 Pw 10000.0 /* Volumic power */ 21 #define LAMBDA 10.0 /* Lambda of the middle slab */ 22 #define LAMBDA1 1.0 /* Lambda of the upper slab */ 23 #define LAMBDA2 LAMBDA1 /* Lambda of the lower slab */ 24 #define T1 373.15 /* Temperature of the upper fluid */ 25 #define T2 273.15 /* Temperature of the lower fluid */ 26 #define H1 5.0 /* Convection coef between the upper slab and the fluid */ 27 #define H2 10.0 /* Convection coef between the lower slab and the fluid */ 28 #define DELTA 0.01 /* Delta of the middle slab */ 29 #define DELTA1 0.02 /* Delta of the upper slab */ 30 #define DELTA2 0.07 /* Delta of the lower slab */ 31 #define L 0.2 /* Size of the middle slab */ 32 #define L1 0.4 /* Size of the upper slab */ 33 #define L2 1.4 /* Size of the lower slab */ 34 35 #define N 10000 /* #realisations */ 36 37 /* Analitically computed temperatures wrt the previous parameters.*/ 38 #define Tp1 648.6217 39 #define Tp2 335.4141 40 #define Ta 1199.5651 41 #define Tb 1207.1122 42 43 /* Fixed temperatures */ 44 #define Tsolid1_fluid SDIS_TEMPERATURE_NONE /*Tp1*/ 45 #define Tsolid2_fluid SDIS_TEMPERATURE_NONE /*Tp2*/ 46 #define Tsolid_solid1 SDIS_TEMPERATURE_NONE /*Ta*/ 47 #define Tsolid_solid2 SDIS_TEMPERATURE_NONE /*Tb*/ 48 49 #define PROBE_POS 1.8 50 51 /* 52 * The 2D scene is composed of 3 stacked solid slabs whose middle slab has a 53 * volumic power. The +/-X sides of the slabs are stretched far away to 54 * simulate a 1D case. The upper and lower bounds of the "sandwich" has a 55 * convective exchange with the surrounding fluid whose temperature is known. 56 * 57 * _\ T1 58 * / / 59 * \__/ 60 * ... -----H1------ ... Tp1 61 * LAMBDA1 62 * 63 * ... ------------- ... Ta 64 * LAMBDA, Pw 65 * ... ------------- ... Tb 66 * 67 * LAMBDA2 68 * 69 * 70 * ... -----H2------ ... Tp2 71 * _\ T2 72 * / / 73 * \__/ 74 */ 75 76 static const double vertices[8/*#vertices*/*2/*#coords per vertex*/] = { 77 -100000.5, 0.0, /* 0 */ 78 -100000.5, 1.4, /* 1 */ 79 -100000.5, 1.6, /* 2 */ 80 -100000.5, 2.0, /* 3 */ 81 100000.5, 2.0, /* 4 */ 82 100000.5, 1.6, /* 5 */ 83 100000.5, 1.4, /* 6 */ 84 100000.5, 0.0 /* 7 */ 85 }; 86 static const size_t nvertices = sizeof(vertices)/(sizeof(double)*2); 87 88 static const size_t indices[10/*#segments*/*2/*#indices per segment*/]= { 89 0, 1, 90 1, 2, 91 2, 3, 92 3, 4, 93 4, 5, 94 5, 6, 95 6, 7, 96 7, 0, 97 6, 1, 98 2, 5 99 }; 100 static const size_t nsegments = sizeof(indices)/(sizeof(size_t)*2); 101 102 /******************************************************************************* 103 * Geometry 104 ******************************************************************************/ 105 static void 106 get_indices(const size_t iseg, size_t ids[2], void* context) 107 { 108 (void)context; 109 CHK(ids); 110 ids[0] = indices[iseg*2+0]; 111 ids[1] = indices[iseg*2+1]; 112 } 113 114 static void 115 get_position(const size_t ivert, double pos[2], void* context) 116 { 117 (void)context; 118 CHK(pos); 119 pos[0] = vertices[ivert*2+0]; 120 pos[1] = vertices[ivert*2+1]; 121 } 122 123 static void 124 get_interface(const size_t iseg, struct sdis_interface** bound, void* context) 125 { 126 struct sdis_interface** interfaces = context; 127 CHK(context && bound); 128 *bound = interfaces[iseg]; 129 } 130 131 /******************************************************************************* 132 * Solid medium 133 ******************************************************************************/ 134 struct solid { 135 double cp; 136 double lambda; 137 double rho; 138 double delta; 139 double volumic_power; 140 double temperature; 141 }; 142 143 static double 144 solid_get_calorific_capacity 145 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 146 { 147 CHK(data != NULL && vtx != NULL); 148 return ((const struct solid*)sdis_data_cget(data))->cp; 149 } 150 151 static double 152 solid_get_thermal_conductivity 153 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 154 { 155 CHK(data != NULL && vtx != NULL); 156 return ((const struct solid*)sdis_data_cget(data))->lambda; 157 } 158 159 static double 160 solid_get_volumic_mass 161 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 162 { 163 CHK(data != NULL && vtx != NULL); 164 return ((const struct solid*)sdis_data_cget(data))->rho; 165 } 166 167 static double 168 solid_get_delta 169 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 170 { 171 CHK(data != NULL && vtx != NULL); 172 return ((const struct solid*)sdis_data_cget(data))->delta; 173 } 174 175 static double 176 solid_get_temperature 177 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 178 { 179 CHK(data != NULL && vtx != NULL); 180 return ((const struct solid*)sdis_data_cget(data))->temperature; 181 } 182 183 static double 184 solid_get_volumic_power 185 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 186 { 187 CHK(data != NULL && vtx != NULL); 188 return ((const struct solid*)sdis_data_cget(data))->volumic_power; 189 } 190 191 /******************************************************************************* 192 * Fluid medium 193 ******************************************************************************/ 194 struct fluid { 195 double temperature_lower; 196 double temperature_upper; 197 }; 198 199 static double 200 fluid_get_temperature 201 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 202 { 203 const struct fluid* fluid; 204 CHK(data != NULL && vtx != NULL); 205 fluid = sdis_data_cget(data); 206 return vtx->P[1] < 1 ? fluid->temperature_lower : fluid->temperature_upper; 207 } 208 209 210 /******************************************************************************* 211 * Interfaces 212 ******************************************************************************/ 213 struct interf { 214 double h; 215 double temperature; 216 }; 217 218 static double 219 interface_get_convection_coef 220 (const struct sdis_interface_fragment* frag, struct sdis_data* data) 221 { 222 CHK(frag && data); 223 return ((const struct interf*)sdis_data_cget(data))->h; 224 } 225 226 static double 227 interface_get_temperature 228 (const struct sdis_interface_fragment* frag, struct sdis_data* data) 229 { 230 CHK(frag && data); 231 return ((const struct interf*)sdis_data_cget(data))->temperature; 232 } 233 234 /******************************************************************************* 235 * Test 236 ******************************************************************************/ 237 int 238 main(int argc, char** argv) 239 { 240 struct solid* solid_param = NULL; 241 struct fluid* fluid_param = NULL; 242 struct interf* interf_param = NULL; 243 struct sdis_device* dev = NULL; 244 struct sdis_data* data = NULL; 245 struct sdis_medium* fluid = NULL; 246 struct sdis_medium* solid = NULL; 247 struct sdis_medium* solid1 = NULL; 248 struct sdis_medium* solid2 = NULL; 249 struct sdis_scene* scn = NULL; 250 struct sdis_estimator* estimator = NULL; 251 struct sdis_fluid_shader fluid_shader = SDIS_FLUID_SHADER_NULL; 252 struct sdis_solid_shader solid_shader = SDIS_SOLID_SHADER_NULL; 253 struct sdis_interface_shader interf_shader = SDIS_INTERFACE_SHADER_NULL; 254 struct sdis_interface* interf_solid_adiabatic = NULL; 255 struct sdis_interface* interf_solid1_adiabatic = NULL; 256 struct sdis_interface* interf_solid2_adiabatic = NULL; 257 struct sdis_interface* interf_solid_solid1 = NULL; 258 struct sdis_interface* interf_solid_solid2 = NULL; 259 struct sdis_interface* interf_solid1_fluid = NULL; 260 struct sdis_interface* interf_solid2_fluid = NULL; 261 struct sdis_interface* interfaces[10/*#segment*/]; 262 struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT; 263 struct sdis_solve_probe_args solve_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT; 264 struct sdis_mc T = SDIS_MC_NULL; 265 double Tref; 266 double time_range[2]; 267 double pos[2]; 268 size_t nfails; 269 size_t nreals; 270 (void)argc, (void)argv; 271 272 OK(sdis_device_create(&SDIS_DEVICE_CREATE_ARGS_DEFAULT, &dev)); 273 274 time_range[0] = INF; 275 time_range[1] = INF; 276 277 /* Create the fluid medium */ 278 fluid_shader.temperature = fluid_get_temperature; 279 fluid_shader.calorific_capacity = dummy_medium_getter; 280 fluid_shader.volumic_mass = dummy_medium_getter; 281 OK(sdis_data_create 282 (dev, sizeof(struct fluid), ALIGNOF(struct fluid), NULL, &data)); 283 fluid_param = sdis_data_get(data); 284 fluid_param->temperature_upper = T1; 285 fluid_param->temperature_lower = T2; 286 OK(sdis_fluid_create(dev, &fluid_shader, data, &fluid)); 287 OK(sdis_data_ref_put(data)); 288 289 /* Setup the solid shader */ 290 solid_shader.calorific_capacity = solid_get_calorific_capacity; 291 solid_shader.thermal_conductivity = solid_get_thermal_conductivity; 292 solid_shader.volumic_mass = solid_get_volumic_mass; 293 solid_shader.delta = solid_get_delta; 294 solid_shader.temperature = solid_get_temperature; 295 solid_shader.volumic_power = solid_get_volumic_power; 296 297 /* Create the medium of the upper slab */ 298 OK(sdis_data_create 299 (dev, sizeof(struct solid), ALIGNOF(struct solid), NULL, &data)); 300 solid_param = sdis_data_get(data); 301 solid_param->cp = 500000; 302 solid_param->rho = 1000; 303 solid_param->lambda = LAMBDA1; 304 solid_param->delta = DELTA1; 305 solid_param->volumic_power = SDIS_VOLUMIC_POWER_NONE; 306 solid_param->temperature = SDIS_TEMPERATURE_NONE; 307 OK(sdis_solid_create(dev, &solid_shader, data, &solid1)); 308 OK(sdis_data_ref_put(data)); 309 310 /* Create the medium of the lower slab */ 311 OK(sdis_data_create 312 (dev, sizeof(struct solid), ALIGNOF(struct solid), NULL, &data)); 313 solid_param = sdis_data_get(data); 314 solid_param->cp = 500000; 315 solid_param->rho = 1000; 316 solid_param->lambda = LAMBDA2; 317 solid_param->delta = DELTA2; 318 solid_param->volumic_power = SDIS_VOLUMIC_POWER_NONE; 319 solid_param->temperature = SDIS_TEMPERATURE_NONE; 320 OK(sdis_solid_create(dev, &solid_shader, data, &solid2)); 321 OK(sdis_data_ref_put(data)); 322 323 /* Create the medium of the middle slab */ 324 OK(sdis_data_create 325 (dev, sizeof(struct solid), ALIGNOF(struct solid), NULL, &data)); 326 solid_param = sdis_data_get(data); 327 solid_param->cp = 500000; 328 solid_param->rho = 1000; 329 solid_param->lambda = LAMBDA; 330 solid_param->delta = DELTA; 331 solid_param->volumic_power = Pw; 332 solid_param->temperature = SDIS_TEMPERATURE_NONE; 333 OK(sdis_solid_create(dev, &solid_shader, data, &solid)); 334 OK(sdis_data_ref_put(data)); 335 336 interf_shader.front.temperature = interface_get_temperature; 337 338 /* Create the solid/solid1 interface */ 339 OK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf), 340 NULL, &data)); 341 interf_param = sdis_data_get(data); 342 interf_param->temperature = Tsolid_solid1; 343 OK(sdis_interface_create(dev, solid, solid1, &interf_shader, 344 data, &interf_solid_solid1)); 345 OK(sdis_data_ref_put(data)); 346 347 /* Create the solid/solid2 interface */ 348 OK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf), 349 NULL, &data)); 350 interf_param = sdis_data_get(data); 351 interf_param->temperature = Tsolid_solid2; 352 OK(sdis_interface_create(dev, solid, solid2, &interf_shader, 353 data, &interf_solid_solid2)); 354 OK(sdis_data_ref_put(data)); 355 356 /* Setup the interface shader */ 357 interf_shader.convection_coef = interface_get_convection_coef; 358 interf_shader.front.temperature = interface_get_temperature; 359 360 /* Create the adiabatic interfaces */ 361 OK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf), 362 NULL, &data)); 363 interf_param = sdis_data_get(data); 364 interf_param->h = 0; 365 interf_param->temperature = SDIS_TEMPERATURE_NONE; 366 OK(sdis_interface_create(dev, solid, fluid, &interf_shader, data, 367 &interf_solid_adiabatic)); 368 OK(sdis_interface_create(dev, solid1, fluid, &interf_shader, data, 369 &interf_solid1_adiabatic)); 370 OK(sdis_interface_create(dev, solid2, fluid, &interf_shader, data, 371 &interf_solid2_adiabatic) ); 372 OK(sdis_data_ref_put(data)); 373 374 /* Create the solid1 fluid interace */ 375 OK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf), 376 NULL, &data)); 377 interf_param = sdis_data_get(data); 378 interf_param->h = H1; 379 interf_param->temperature = Tsolid1_fluid; 380 OK(sdis_interface_create(dev, solid1, fluid, &interf_shader, data, 381 &interf_solid1_fluid)); 382 OK(sdis_data_ref_put(data)); 383 384 /* Create the solid2 fluid 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 = H2; 389 interf_param->temperature = Tsolid2_fluid; 390 OK(sdis_interface_create(dev, solid2, fluid, &interf_shader, data, 391 &interf_solid2_fluid)); 392 OK(sdis_data_ref_put(data)); 393 394 /* Release the media */ 395 OK(sdis_medium_ref_put(fluid)); 396 OK(sdis_medium_ref_put(solid)); 397 OK(sdis_medium_ref_put(solid1)); 398 OK(sdis_medium_ref_put(solid2)); 399 400 /* Map the interfaces to their square segments */ 401 interfaces[0] = interf_solid2_adiabatic; 402 interfaces[1] = interf_solid_adiabatic; 403 interfaces[2] = interf_solid1_adiabatic; 404 interfaces[3] = interf_solid1_fluid; 405 interfaces[4] = interf_solid1_adiabatic; 406 interfaces[5] = interf_solid_adiabatic; 407 interfaces[6] = interf_solid2_adiabatic; 408 interfaces[7] = interf_solid2_fluid; 409 interfaces[8] = interf_solid_solid2; 410 interfaces[9] = interf_solid_solid1; 411 412 #if 0 413 dump_segments(stdout, vertices, nvertices, indices, nsegments); 414 exit(0); 415 #endif 416 417 /* Create the scene */ 418 scn_args.get_indices = get_indices; 419 scn_args.get_interface = get_interface; 420 scn_args.get_position = get_position; 421 scn_args.nprimitives = nsegments; 422 scn_args.nvertices = nvertices; 423 scn_args.context = interfaces; 424 OK(sdis_scene_2d_create(dev, &scn_args, &scn)); 425 426 /* Release the interfaces */ 427 OK(sdis_interface_ref_put(interf_solid_adiabatic)); 428 OK(sdis_interface_ref_put(interf_solid1_adiabatic)); 429 OK(sdis_interface_ref_put(interf_solid2_adiabatic)); 430 OK(sdis_interface_ref_put(interf_solid_solid1)); 431 OK(sdis_interface_ref_put(interf_solid_solid2)); 432 OK(sdis_interface_ref_put(interf_solid1_fluid)); 433 OK(sdis_interface_ref_put(interf_solid2_fluid)); 434 435 pos[0] = 0; 436 pos[1] = PROBE_POS; 437 438 if(pos[1] > 0 && pos[1] < L2) { /* Lower slab */ 439 Tref = Tp2 + (Tb - Tp2) * pos[1] / L2; 440 } else if(pos[1] > L2 && pos[1] < L2 + L) { /* Middle slab */ 441 Tref = 442 (Ta + Tb) / 2 443 + (Ta - Tb)/L * (pos[1] - (L2+L/2)) 444 + Pw * (L*L/4.0 - pow((pos[1] - (L2+L/2)), 2)) / (2*LAMBDA); 445 } else if(pos[1] > L2 + L && pos[1] < L2 + L1 + L) { 446 Tref = Ta + (Tp1 - Ta) / L1 * (pos[1] - (L+L2)); 447 } else { 448 FATAL("Unreachable code.\n"); 449 } 450 451 solve_args.nrealisations = N; 452 solve_args.position[0] = pos[0]; 453 solve_args.position[1] = pos[1]; 454 solve_args.time_range[0] = time_range[0]; 455 solve_args.time_range[1] = time_range[1]; 456 OK(sdis_solve_probe(scn, &solve_args, &estimator)); 457 OK(sdis_estimator_get_temperature(estimator, &T)); 458 OK(sdis_estimator_get_failure_count(estimator, &nfails)); 459 OK(sdis_estimator_get_realisation_count(estimator, &nreals)); 460 printf("Temperature at (%g %g) = %g ~ %g +/- %g [%g, %g]\n", 461 SPLIT2(pos), Tref, T.E, T.SE, T.E-3*T.SE, T.E+3*T.SE); 462 printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); 463 OK(sdis_estimator_ref_put(estimator)); 464 465 CHK(nfails + nreals == N); 466 CHK(nfails < N/1000); 467 CHK(eq_eps(T.E, Tref, T.SE*3)); 468 469 OK(sdis_scene_ref_put(scn)); 470 OK(sdis_device_ref_put(dev)); 471 472 CHK(mem_allocated_size() == 0); 473 return 0; 474 } 475