test_sdis_volumic_power4.c (13138B)
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/clock_time.h> 19 #include <rsys/math.h> 20 21 #define Tf 100.0 22 #define Power 10000.0 23 #define H 50.0 24 #define LAMBDA 100.0 25 #define DELTA 0.2/*(1.0/2.0)*/ 26 #define N 100000 27 #define LENGTH 10000.0 28 29 /* 30 * The 2D scene is a solid slabs stretched along the X dimension to simulate a 31 * 1D case. The slab has a volumic power and has a convective exchange with 32 * surrounding fluid whose temperature is fixed to Tf. 33 * 34 * 35 * _\ Tf 36 * / / 37 * \__/ 38 * 39 * ... -----Hboundary----- ... 40 * 41 * Lambda, Power 42 * 43 * ... -----Hboundary----- ... 44 * 45 * _\ Tf 46 * / / 47 * \__/ 48 * 49 */ 50 51 static const double vertices_2d[4/*#vertices*/*2/*#coords per vertex*/] = { 52 LENGTH,-0.5, 53 -LENGTH,-0.5, 54 -LENGTH, 0.5, 55 LENGTH, 0.5 56 }; 57 58 static const double vertices_3d[8/*#vertices*/*3/*#coords per vertex*/] = { 59 -LENGTH,-0.5,-LENGTH, 60 LENGTH,-0.5,-LENGTH, 61 -LENGTH, 0.5,-LENGTH, 62 LENGTH, 0.5,-LENGTH, 63 -LENGTH,-0.5, LENGTH, 64 LENGTH,-0.5, LENGTH, 65 -LENGTH, 0.5, LENGTH, 66 LENGTH, 0.5, LENGTH 67 }; 68 69 /******************************************************************************* 70 * Geometry 71 ******************************************************************************/ 72 static void 73 get_position_2d(const size_t ivert, double pos[2], void* context) 74 { 75 (void)context; 76 CHK(pos); 77 pos[0] = vertices_2d[ivert*2+0]; 78 pos[1] = vertices_2d[ivert*2+1]; 79 } 80 81 static void 82 get_position_3d(const size_t ivert, double pos[3], void* context) 83 { 84 (void)context; 85 CHK(pos); 86 pos[0] = vertices_3d[ivert*3+0]; 87 pos[1] = vertices_3d[ivert*3+1]; 88 pos[2] = vertices_3d[ivert*3+2]; 89 } 90 91 static void 92 get_interface(const size_t iprim, struct sdis_interface** bound, void* context) 93 { 94 struct sdis_interface** interfaces = context; 95 CHK(context && bound); 96 *bound = interfaces[iprim]; 97 } 98 99 /******************************************************************************* 100 * Solid medium 101 ******************************************************************************/ 102 struct solid { 103 double cp; 104 double lambda; 105 double rho; 106 double delta; 107 double volumic_power; 108 double temperature; 109 }; 110 111 static double 112 solid_get_calorific_capacity 113 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 114 { 115 CHK(data != NULL && vtx != NULL); 116 return ((const struct solid*)sdis_data_cget(data))->cp; 117 } 118 119 static double 120 solid_get_thermal_conductivity 121 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 122 { 123 CHK(data != NULL && vtx != NULL); 124 return ((const struct solid*)sdis_data_cget(data))->lambda; 125 } 126 127 static double 128 solid_get_volumic_mass 129 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 130 { 131 CHK(data != NULL && vtx != NULL); 132 return ((const struct solid*)sdis_data_cget(data))->rho; 133 } 134 135 static double 136 solid_get_delta 137 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 138 { 139 CHK(data != NULL && vtx != NULL); 140 return ((const struct solid*)sdis_data_cget(data))->delta; 141 } 142 143 static double 144 solid_get_temperature 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))->temperature; 149 } 150 151 static double 152 solid_get_volumic_power 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))->volumic_power; 157 } 158 159 /******************************************************************************* 160 * Fluid medium 161 ******************************************************************************/ 162 struct fluid { 163 double temperature; 164 }; 165 166 static double 167 fluid_get_temperature 168 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 169 { 170 const struct fluid* fluid; 171 CHK(data != NULL && vtx != NULL); 172 fluid = sdis_data_cget(data); 173 return fluid->temperature; 174 } 175 176 /******************************************************************************* 177 * Interfaces 178 ******************************************************************************/ 179 struct interf { 180 double h; 181 double temperature; 182 }; 183 184 static double 185 interface_get_convection_coef 186 (const struct sdis_interface_fragment* frag, struct sdis_data* data) 187 { 188 CHK(frag && data); 189 return ((const struct interf*)sdis_data_cget(data))->h; 190 } 191 192 static double 193 interface_get_temperature 194 (const struct sdis_interface_fragment* frag, struct sdis_data* data) 195 { 196 CHK(frag && data); 197 return ((const struct interf*)sdis_data_cget(data))->temperature; 198 } 199 200 /******************************************************************************* 201 * Test 202 ******************************************************************************/ 203 int 204 main(int argc, char** argv) 205 { 206 char dump[128]; 207 struct time t0, t1; 208 struct solid* solid_param = NULL; 209 struct fluid* fluid_param = NULL; 210 struct interf* interf_param = NULL; 211 struct sdis_device* dev = NULL; 212 struct sdis_data* data = NULL; 213 struct sdis_medium* fluid1 = NULL; 214 struct sdis_medium* fluid2 = NULL; 215 struct sdis_medium* solid = NULL; 216 struct sdis_scene* scn_2d = NULL; 217 struct sdis_scene* scn_3d = NULL; 218 struct sdis_estimator* estimator = NULL; 219 struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT; 220 struct sdis_fluid_shader fluid_shader = SDIS_FLUID_SHADER_NULL; 221 struct sdis_solid_shader solid_shader = SDIS_SOLID_SHADER_NULL; 222 struct sdis_interface_shader interf_shader = SDIS_INTERFACE_SHADER_NULL; 223 struct sdis_interface* interf_adiabatic = NULL; 224 struct sdis_interface* interf_solid_fluid1 = NULL; 225 struct sdis_interface* interf_solid_fluid2 = NULL; 226 struct sdis_interface* interfaces[12/*#max primitives*/]; 227 struct sdis_mc T = SDIS_MC_NULL; 228 struct sdis_solve_probe_args solve_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT; 229 size_t nreals, nfails; 230 double pos[3]; 231 double Tref; 232 double x; 233 (void)argc, (void)argv; 234 235 OK(sdis_device_create(&SDIS_DEVICE_CREATE_ARGS_DEFAULT, &dev)); 236 237 /* Create the fluid medium */ 238 fluid_shader.temperature = fluid_get_temperature; 239 fluid_shader.calorific_capacity = dummy_medium_getter; 240 fluid_shader.volumic_mass = dummy_medium_getter; 241 242 OK(sdis_data_create 243 (dev, sizeof(struct fluid), ALIGNOF(struct fluid), NULL, &data)); 244 fluid_param = sdis_data_get(data); 245 fluid_param->temperature = Tf; 246 OK(sdis_fluid_create(dev, &fluid_shader, data, &fluid1)); 247 OK(sdis_data_ref_put(data)); 248 249 OK(sdis_data_create 250 (dev, sizeof(struct fluid), ALIGNOF(struct fluid), NULL, &data)); 251 fluid_param = sdis_data_get(data); 252 fluid_param->temperature = Tf; 253 OK(sdis_fluid_create(dev, &fluid_shader, data, &fluid2)); 254 OK(sdis_data_ref_put(data)); 255 256 /* Setup the solid shader */ 257 solid_shader.calorific_capacity = solid_get_calorific_capacity; 258 solid_shader.thermal_conductivity = solid_get_thermal_conductivity; 259 solid_shader.volumic_mass = solid_get_volumic_mass; 260 solid_shader.delta = solid_get_delta; 261 solid_shader.temperature = solid_get_temperature; 262 solid_shader.volumic_power = solid_get_volumic_power; 263 264 /* Create the solid medium */ 265 OK(sdis_data_create 266 (dev, sizeof(struct solid), ALIGNOF(struct solid), NULL, &data)); 267 solid_param = sdis_data_get(data); 268 solid_param->cp = 500000; 269 solid_param->rho = 1000; 270 solid_param->lambda = LAMBDA; 271 solid_param->delta = DELTA; 272 solid_param->volumic_power = Power; 273 solid_param->temperature = SDIS_TEMPERATURE_NONE; 274 OK(sdis_solid_create(dev, &solid_shader, data, &solid)); 275 OK(sdis_data_ref_put(data)); 276 277 /* Setup the interface shader */ 278 interf_shader.convection_coef = interface_get_convection_coef; 279 interf_shader.front.temperature = interface_get_temperature; 280 281 /* Create the adiabatic interface */ 282 OK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf), 283 NULL, &data)); 284 interf_param = sdis_data_get(data); 285 interf_param->h = 0; 286 interf_param->temperature = SDIS_TEMPERATURE_NONE; 287 OK(sdis_interface_create(dev, solid, fluid1, &interf_shader, data, 288 &interf_adiabatic)); 289 OK(sdis_data_ref_put(data)); 290 291 /* Create the solid fluid1 interface */ 292 OK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf), 293 NULL, &data)); 294 interf_param = sdis_data_get(data); 295 interf_param->h = H; 296 interf_param->temperature = SDIS_TEMPERATURE_NONE; 297 OK(sdis_interface_create(dev, solid, fluid1, &interf_shader, data, 298 &interf_solid_fluid1)); 299 OK(sdis_data_ref_put(data)); 300 301 /* Create the solid fluid2 interface */ 302 OK(sdis_data_create (dev, sizeof(struct interf), ALIGNOF(struct interf), 303 NULL, &data)); 304 interf_param = sdis_data_get(data); 305 interf_param->h = H; 306 interf_param->temperature = SDIS_TEMPERATURE_NONE; 307 OK(sdis_interface_create(dev, solid, fluid2, &interf_shader, data, 308 &interf_solid_fluid2)); 309 OK(sdis_data_ref_put(data)); 310 311 /* Release the media */ 312 OK(sdis_medium_ref_put(fluid1)); 313 OK(sdis_medium_ref_put(fluid2)); 314 OK(sdis_medium_ref_put(solid)); 315 316 #if 0 317 dump_segments(stdout, vertices, nvertices, indices, nsegments); 318 exit(0); 319 #endif 320 321 /* Map the interfaces to their square segments */ 322 interfaces[0] = interf_solid_fluid2; /* Bottom */ 323 interfaces[1] = interf_adiabatic; /* Left */ 324 interfaces[2] = interf_solid_fluid1; /* Top */ 325 interfaces[3] = interf_adiabatic; /* Right */ 326 327 /* Create the 2D scene */ 328 scn_args.get_indices = square_get_indices; 329 scn_args.get_interface = get_interface; 330 scn_args.get_position = get_position_2d; 331 scn_args.nprimitives = square_nsegments; 332 scn_args.nvertices = square_nvertices; 333 scn_args.context = interfaces; 334 OK(sdis_scene_2d_create(dev, &scn_args, &scn_2d)); 335 336 /* Map the interfaces to their box triangles */ 337 interfaces[0] = interfaces[1] = interf_adiabatic; /* Front */ 338 interfaces[2] = interfaces[3] = interf_adiabatic; /* Left */ 339 interfaces[4] = interfaces[5] = interf_adiabatic; /* Back */ 340 interfaces[6] = interfaces[7] = interf_adiabatic; /* Right */ 341 interfaces[8] = interfaces[9] = interf_solid_fluid1; /* Top */ 342 interfaces[10]= interfaces[11]= interf_solid_fluid2; /* Bottom */ 343 344 /* Create the 3D scene */ 345 scn_args.get_indices = box_get_indices; 346 scn_args.get_interface = get_interface; 347 scn_args.get_position = get_position_3d; 348 scn_args.nprimitives = box_ntriangles; 349 scn_args.nvertices = box_nvertices; 350 scn_args.context = interfaces; 351 OK(sdis_scene_create(dev, &scn_args, &scn_3d)); 352 353 /* Release the interfaces */ 354 OK(sdis_interface_ref_put(interf_adiabatic)); 355 OK(sdis_interface_ref_put(interf_solid_fluid1)); 356 OK(sdis_interface_ref_put(interf_solid_fluid2)); 357 358 pos[0] = 0; 359 pos[1] = 0; 360 pos[2] = 0; 361 362 x = pos[1]; 363 Tref = -Power / (2*LAMBDA) * x*x + Tf + Power/(2*H) + Power/(8*LAMBDA); 364 365 solve_args.nrealisations = N; 366 solve_args.position[0] = pos[0]; 367 solve_args.position[1] = pos[1]; 368 solve_args.position[2] = pos[2]; 369 solve_args.time_range[0] = INF; 370 solve_args.time_range[1] = INF; 371 372 printf(">>> 2D\n"); 373 374 time_current(&t0); 375 OK(sdis_solve_probe(scn_2d, &solve_args, &estimator)); 376 time_sub(&t0, time_current(&t1), &t0); 377 time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); 378 printf("Elapsed time = %s\n", dump); 379 380 OK(sdis_estimator_get_temperature(estimator, &T)); 381 OK(sdis_estimator_get_realisation_count(estimator, &nreals)); 382 OK(sdis_estimator_get_failure_count(estimator, &nfails)); 383 printf("Temperature at (%g %g) = %g ~ %g +/- %g [%g %g]\n", 384 SPLIT2(pos), Tref, T.E, T.SE, T.E-3*T.SE, T.E+3*T.SE); 385 printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); 386 OK(sdis_estimator_ref_put(estimator)); 387 CHK(nfails + nreals == N); 388 CHK(nfails < N/1000); 389 CHK(eq_eps(T.E, Tref, T.SE*3)); 390 391 printf("\n>>> 3D\n"); 392 393 time_current(&t0); 394 OK(sdis_solve_probe(scn_3d, &solve_args, &estimator)); 395 time_sub(&t0, time_current(&t1), &t0); 396 time_dump(&t0, TIME_ALL, NULL, dump, sizeof(dump)); 397 printf("Elapsed time = %s\n", dump); 398 399 OK(sdis_estimator_get_temperature(estimator, &T)); 400 OK(sdis_estimator_get_realisation_count(estimator, &nreals)); 401 OK(sdis_estimator_get_failure_count(estimator, &nfails)); 402 printf("Temperature at (%g %g) = %g ~ %g +/- %g [%g %g]\n", 403 SPLIT2(pos), Tref, T.E, T.SE, T.E-3*T.SE, T.E+3*T.SE); 404 printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); 405 OK(sdis_estimator_ref_put(estimator)); 406 CHK(nfails + nreals == N); 407 CHK(nfails < N/1000); 408 CHK(eq_eps(T.E, Tref, T.SE*3)); 409 410 OK(sdis_scene_ref_put(scn_2d)); 411 OK(sdis_scene_ref_put(scn_3d)); 412 OK(sdis_device_ref_put(dev)); 413 414 CHK(mem_allocated_size() == 0); 415 return 0; 416 } 417