test_sdis_solve_medium_2d.c (15399B)
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/stretchy_array.h> 20 #include <rsys/math.h> 21 22 #include <string.h> 23 24 #define Tf0 300.0 25 #define Tf1 330.0 26 #define N 1000ul /* #realisations */ 27 #define Np 10000ul /* #realisations precise */ 28 29 /* 30 * The scene is composed of a square and a disk whose temperature is unknown. 31 * The square is surrounded by a fluid whose temperature is Tf0 while the disk 32 * is in a fluid whose temperature is Tf1. The temperature of the square 33 * and the disk are thus uniform and equal to Tf0 and Tf1, respectively. 34 * 35 * # # Tf1 +---------+ 36 * # # _\ | | _\ Tf0 37 * # # / / | | / / 38 * # # \__/ | | \__/ 39 * # # | | 40 * # # +---------+ 41 * 42 * This program performs 2 tests. In the first one, the square and the disk 43 * have different media; the medium solver estimates the temperature of 44 * the square or the one of the disk. In the second test, the scene is updated 45 * to use the same medium for the 2 shapes. When invoked on 46 * the right medium, the estimated temperature T is equal to : 47 * 48 * T = Tf0 * A0/(A0 + A1) + Tf1 * A1/(A0 + A1) 49 * 50 * with A0 and A1 the area of the shape and the area of the disk, respectively. 51 */ 52 53 /******************************************************************************* 54 * Geometry 55 ******************************************************************************/ 56 struct context { 57 const double* positions; 58 const size_t* indices; 59 size_t nsegments_interf0; /* #segments of the interface 0 */ 60 struct sdis_interface* interf0; 61 struct sdis_interface* interf1; 62 }; 63 64 static void 65 get_indices(const size_t iseg, size_t ids[2], void* context) 66 { 67 const struct context* ctx = context; 68 ids[0] = ctx->indices[iseg*2+0]; 69 ids[1] = ctx->indices[iseg*2+1]; 70 } 71 72 static void 73 get_position(const size_t ivert, double pos[2], void* context) 74 { 75 const struct context* ctx = context; 76 pos[0] = ctx->positions[ivert*2+0]; 77 pos[1] = ctx->positions[ivert*2+1]; 78 } 79 80 static void 81 get_interface(const size_t iseg, struct sdis_interface** bound, void* context) 82 { 83 const struct context* ctx = context; 84 *bound = iseg < ctx->nsegments_interf0 ? ctx->interf0 : ctx->interf1; 85 } 86 87 /******************************************************************************* 88 * Fluid medium 89 ******************************************************************************/ 90 struct fluid { 91 double temperature; 92 }; 93 94 static double 95 fluid_get_temperature 96 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 97 { 98 CHK(data != NULL && vtx != NULL); 99 return ((const struct fluid*)sdis_data_cget(data))->temperature; 100 } 101 102 /******************************************************************************* 103 * Solid medium 104 ******************************************************************************/ 105 struct solid { 106 double cp; 107 double lambda; 108 double rho; 109 double delta; 110 double temperature; 111 }; 112 113 static double 114 solid_get_calorific_capacity 115 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 116 { 117 CHK(data != NULL && vtx != NULL); 118 return ((const struct solid*)sdis_data_cget(data))->cp; 119 } 120 121 static double 122 solid_get_thermal_conductivity 123 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 124 { 125 CHK(data != NULL && vtx != NULL); 126 return ((const struct solid*)sdis_data_cget(data))->lambda; 127 } 128 129 static double 130 solid_get_volumic_mass 131 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 132 { 133 CHK(data != NULL && vtx != NULL); 134 return ((const struct solid*)sdis_data_cget(data))->rho; 135 } 136 137 static double 138 solid_get_delta 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))->delta; 143 } 144 145 static double 146 solid_get_temperature 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))->temperature; 151 } 152 153 struct interf { 154 double hc; 155 double epsilon; 156 double specular_fraction; 157 }; 158 159 static double 160 interface_get_convection_coef 161 (const struct sdis_interface_fragment* frag, struct sdis_data* data) 162 { 163 CHK(data != NULL && frag != NULL); 164 return ((const struct interf*)sdis_data_cget(data))->hc; 165 } 166 167 static double 168 interface_get_emissivity 169 (const struct sdis_interface_fragment* frag, 170 const unsigned source_id, 171 struct sdis_data* data) 172 { 173 (void)source_id; 174 CHK(data != NULL && frag != NULL); 175 return ((const struct interf*)sdis_data_cget(data))->epsilon; 176 } 177 178 static double 179 interface_get_specular_fraction 180 (const struct sdis_interface_fragment* frag, 181 const unsigned source_id, 182 struct sdis_data* data) 183 { 184 (void)source_id; 185 CHK(data != NULL && frag != NULL); 186 return ((const struct interf*)sdis_data_cget(data))->specular_fraction; 187 } 188 189 /******************************************************************************* 190 * Test 191 ******************************************************************************/ 192 int 193 main(int argc, char** argv) 194 { 195 struct sdis_mc T = SDIS_MC_NULL; 196 struct sdis_mc time = SDIS_MC_NULL; 197 struct sdis_device* dev = NULL; 198 struct sdis_medium* solid0 = NULL; 199 struct sdis_medium* solid1 = NULL; 200 struct sdis_medium* fluid0 = NULL; 201 struct sdis_medium* fluid1 = NULL; 202 struct sdis_interface* solid0_fluid0 = NULL; 203 struct sdis_interface* solid0_fluid1 = NULL; 204 struct sdis_interface* solid1_fluid1 = NULL; 205 struct sdis_scene* scn = NULL; 206 struct sdis_data* data = NULL; 207 struct sdis_estimator* estimator = NULL; 208 struct sdis_estimator* estimator2 = NULL; 209 struct sdis_green_function* green = NULL; 210 struct fluid* fluid_param = NULL; 211 struct solid* solid_param = NULL; 212 struct interf* interface_param = NULL; 213 struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT; 214 struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER; 215 struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER; 216 struct sdis_interface_shader interface_shader = SDIS_INTERFACE_SHADER_NULL; 217 struct sdis_solve_medium_args solve_args = SDIS_SOLVE_MEDIUM_ARGS_DEFAULT; 218 struct context ctx; 219 double a, a0, a1; 220 double ref; 221 double* positions = NULL; 222 size_t* indices = NULL; 223 size_t nverts; 224 size_t nreals; 225 size_t nfails; 226 size_t i; 227 int is_master_process; 228 (void)argc, (void)argv; 229 230 create_default_device(&argc, &argv, &is_master_process, &dev); 231 232 fluid_shader.temperature = fluid_get_temperature; 233 234 /* Create the fluid0 medium */ 235 OK(sdis_data_create 236 (dev, sizeof(struct fluid), ALIGNOF(struct fluid), NULL, &data)); 237 fluid_param = sdis_data_get(data); 238 fluid_param->temperature = Tf0; 239 OK(sdis_fluid_create(dev, &fluid_shader, data, &fluid0)); 240 OK(sdis_data_ref_put(data)); 241 242 /* Create the fluid1 medium */ 243 OK(sdis_data_create 244 (dev, sizeof(struct fluid), ALIGNOF(struct fluid), NULL, &data)); 245 fluid_param = sdis_data_get(data); 246 fluid_param->temperature = Tf1; 247 OK(sdis_fluid_create(dev, &fluid_shader, data, &fluid1)); 248 OK(sdis_data_ref_put(data)); 249 250 /* Setup the solid shader */ 251 solid_shader.calorific_capacity = solid_get_calorific_capacity; 252 solid_shader.thermal_conductivity = solid_get_thermal_conductivity; 253 solid_shader.volumic_mass = solid_get_volumic_mass; 254 solid_shader.delta = solid_get_delta; 255 solid_shader.temperature = solid_get_temperature; 256 257 /* Create the solid0 medium */ 258 OK(sdis_data_create 259 (dev, sizeof(struct solid), ALIGNOF(struct solid), NULL, &data)); 260 solid_param = sdis_data_get(data); 261 solid_param->cp = 1.0; 262 solid_param->lambda = 0.1; 263 solid_param->rho = 1.0; 264 solid_param->delta = 1.0/20.0; 265 solid_param->temperature = SDIS_TEMPERATURE_NONE; 266 OK(sdis_solid_create(dev, &solid_shader, data, &solid0)); 267 OK(sdis_data_ref_put(data)); 268 269 /* Create the solid1 medium */ 270 OK(sdis_data_create 271 (dev, sizeof(struct solid), ALIGNOF(struct solid), NULL, &data)); 272 solid_param = sdis_data_get(data); 273 solid_param->cp = 1.0; 274 solid_param->lambda = 1.0; 275 solid_param->rho = 1.0; 276 solid_param->delta = 1.0/20.0; 277 solid_param->temperature = SDIS_TEMPERATURE_NONE; 278 OK(sdis_solid_create(dev, &solid_shader, data, &solid1)); 279 OK(sdis_data_ref_put(data)); 280 281 /* Create the interfaces */ 282 OK(sdis_data_create(dev, sizeof(struct interf), 283 ALIGNOF(struct interf), NULL, &data)); 284 interface_param = sdis_data_get(data); 285 interface_param->hc = 0.5; 286 interface_param->epsilon = 0; 287 interface_param->specular_fraction = 0; 288 interface_shader.convection_coef = interface_get_convection_coef; 289 interface_shader.front = SDIS_INTERFACE_SIDE_SHADER_NULL; 290 interface_shader.back.temperature = NULL; 291 interface_shader.back.emissivity = interface_get_emissivity; 292 interface_shader.back.specular_fraction = interface_get_specular_fraction; 293 OK(sdis_interface_create 294 (dev, solid0, fluid0, &interface_shader, data, &solid0_fluid0)); 295 OK(sdis_interface_create 296 (dev, solid0, fluid1, &interface_shader, data, &solid0_fluid1)); 297 OK(sdis_interface_create 298 (dev, solid1, fluid1, &interface_shader, data, &solid1_fluid1)); 299 OK(sdis_data_ref_put(data)); 300 301 /* Setup the square geometry */ 302 (void)sa_add(positions, square_nvertices*2); 303 (void)sa_add(indices, square_nsegments*2); 304 memcpy(positions, square_vertices, square_nvertices*sizeof(double[2])); 305 memcpy(indices, square_indices, square_nsegments*sizeof(size_t[2])); 306 307 /* Transate the square in X */ 308 FOR_EACH(i, 0, square_nvertices) positions[i*2] += 2; 309 310 /* Setup a disk */ 311 nverts = 64; 312 FOR_EACH(i, 0, nverts) { 313 const double theta = (double)i * (2*PI)/(double)nverts; 314 const double r = 1; /* Radius */ 315 const double x = cos(theta) * r - 2/* X translation */; 316 const double y = sin(theta) * r + 0.5/* Y translation */; 317 sa_push(positions, x); 318 sa_push(positions, y); 319 } 320 FOR_EACH(i, 0, nverts) { 321 const size_t i0 = i + square_nvertices; 322 const size_t i1 = (i+1) % nverts + square_nvertices; 323 /* Flip the ids to ensure that the normals point inward the disk */ 324 sa_push(indices, i1); 325 sa_push(indices, i0); 326 } 327 328 /* Create the scene */ 329 ctx.positions = positions; 330 ctx.indices = indices; 331 ctx.nsegments_interf0 = square_nsegments; 332 ctx.interf0 = solid0_fluid0; 333 ctx.interf1 = solid1_fluid1; 334 scn_args.get_indices = get_indices; 335 scn_args.get_interface = get_interface; 336 scn_args.get_position = get_position; 337 scn_args.nprimitives = sa_size(indices)/2; 338 scn_args.nvertices = sa_size(positions)/2; 339 scn_args.context = &ctx; 340 OK(sdis_scene_2d_create(dev, &scn_args, &scn)); 341 342 OK(sdis_scene_get_medium_spread(scn, solid0, &a0)); 343 CHK(eq_eps(a0, 1.0, 1.e-6)); 344 OK(sdis_scene_get_medium_spread(scn, solid1, &a1)); 345 /* Rough estimation since the disk is coarsely discretized */ 346 CHK(eq_eps(a1, PI, 1.e-1)); 347 348 solve_args.nrealisations = N; 349 solve_args.time_range[0] = INF; 350 solve_args.time_range[1] = INF; 351 352 /* Estimate the temperature of the square */ 353 solve_args.medium = solid0; 354 OK(sdis_solve_medium(scn, &solve_args, &estimator)); 355 if(!is_master_process) { 356 CHK(estimator == NULL); 357 } else { 358 OK(sdis_estimator_get_temperature(estimator, &T)); 359 OK(sdis_estimator_get_realisation_time(estimator, &time)); 360 OK(sdis_estimator_get_realisation_count(estimator, &nreals)); 361 OK(sdis_estimator_get_failure_count(estimator, &nfails)); 362 printf("Square temperature = "STR(Tf0)" ~ %g +/- %g\n", T.E, T.SE); 363 printf("Time per realisation (in usec) = %g +/- %g\n", time.E, time.SE); 364 printf("#failures = %lu / %lu\n\n", (unsigned long)nfails, N); 365 CHK(eq_eps(T.E, Tf0, T.SE)); 366 CHK(nreals + nfails == N); 367 OK(sdis_estimator_ref_put(estimator)); 368 } 369 370 /* Estimate the temperature of the disk */ 371 solve_args.medium = solid1; 372 OK(sdis_solve_medium(scn, &solve_args, &estimator)); 373 if(is_master_process) { 374 OK(sdis_estimator_get_temperature(estimator, &T)); 375 OK(sdis_estimator_get_realisation_time(estimator, &time)); 376 OK(sdis_estimator_get_realisation_count(estimator, &nreals)); 377 OK(sdis_estimator_get_failure_count(estimator, &nfails)); 378 printf("Disk temperature = "STR(Tf1)" ~ %g +/- %g\n", T.E, T.SE); 379 printf("Time per realisation (in usec) = %g +/- %g\n", time.E, time.SE); 380 printf("#failures = %lu / %lu\n\n", (unsigned long)nfails, N); 381 CHK(eq_eps(T.E, Tf1, T.SE)); 382 CHK(nreals + nfails == N); 383 OK(sdis_estimator_ref_put(estimator)); 384 } 385 386 /* Create a new scene with the same medium for the disk and the square */ 387 OK(sdis_scene_ref_put(scn)); 388 ctx.interf0 = solid0_fluid0; 389 ctx.interf1 = solid0_fluid1; 390 OK(sdis_scene_2d_create(dev, &scn_args, &scn)); 391 392 OK(sdis_scene_get_medium_spread(scn, solid0, &a)); 393 CHK(eq_eps(a, a0+a1, 1.e-6)); 394 395 /* Estimate the temperature of the square and disk shapes */ 396 solve_args.medium = solid1; 397 solve_args.nrealisations = Np; 398 BA(sdis_solve_medium(scn, &solve_args, &estimator)); 399 solve_args.medium = solid0; 400 OK(sdis_solve_medium(scn, &solve_args, &estimator)); 401 if(is_master_process) { 402 OK(sdis_estimator_get_temperature(estimator, &T)); 403 OK(sdis_estimator_get_realisation_time(estimator, &time)); 404 OK(sdis_estimator_get_realisation_count(estimator, &nreals)); 405 OK(sdis_estimator_get_failure_count(estimator, &nfails)); 406 ref = Tf0 * a0/a + Tf1 * a1/a; 407 printf("Square + Disk temperature = %g ~ %g +/- %g\n", ref, T.E, T.SE); 408 printf("Time per realisation (in usec) = %g +/- %g\n", time.E, time.SE); 409 printf("#failures = %lu / %lu\n", (unsigned long)nfails, Np); 410 CHK(eq_eps(T.E, ref, 3*T.SE)); 411 CHK(nreals + nfails == Np); 412 } 413 414 /* Solve green */ 415 BA(sdis_solve_medium_green_function(NULL, &solve_args, &green)); 416 BA(sdis_solve_medium_green_function(scn, NULL, &green)); 417 BA(sdis_solve_medium_green_function(scn, &solve_args, NULL)); 418 OK(sdis_solve_medium_green_function(scn, &solve_args, &green)); 419 420 if(!is_master_process) { 421 CHK(green == NULL); 422 } else { 423 OK(sdis_green_function_solve(green, &estimator2)); 424 check_green_function(green); 425 check_estimator_eq(estimator, estimator2); 426 check_green_serialization(green, scn); 427 428 OK(sdis_green_function_ref_put(green)); 429 430 OK(sdis_estimator_ref_put(estimator)); 431 OK(sdis_estimator_ref_put(estimator2)); 432 } 433 434 /* Release */ 435 OK(sdis_medium_ref_put(solid0)); 436 OK(sdis_medium_ref_put(solid1)); 437 OK(sdis_medium_ref_put(fluid0)); 438 OK(sdis_medium_ref_put(fluid1)); 439 OK(sdis_interface_ref_put(solid0_fluid0)); 440 OK(sdis_interface_ref_put(solid0_fluid1)); 441 OK(sdis_interface_ref_put(solid1_fluid1)); 442 OK(sdis_scene_ref_put(scn)); 443 444 free_default_device(dev); 445 446 sa_release(positions); 447 sa_release(indices); 448 449 CHK(mem_allocated_size() == 0); 450 return 0; 451 } 452