sadist_probe_boundary.c (13771B)
1 /* Copyright (C) 2024 |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 #define _POSIX_C_SOURCE 200112L /* popen */ 17 18 #include "sadist.h" 19 20 #include <stardis/stardis-prog-properties.h> 21 22 #include <star/s3dut.h> 23 24 #include <rsys/cstr.h> 25 #include <rsys/double2.h> 26 #include <rsys/double3.h> 27 #include <rsys/math.h> 28 #include <rsys/mem_allocator.h> 29 30 #include <errno.h> 31 #include <float.h> 32 #include <getopt.h> 33 #include <string.h> /* strerror */ 34 #include <wait.h> /* WIFEXITED, WEXITSTATUS */ 35 36 /* Axis Aligned Bounding Box */ 37 struct aabb { 38 double lower[3]; 39 double upper[3]; 40 }; 41 #define AABB_NULL__ {{DBL_MAX,DBL_MAX,DBL_MAX}, {-DBL_MAX,-DBL_MAX,-DBL_MAX}} 42 static const struct aabb AABB_NULL = AABB_NULL__; 43 44 struct args { 45 unsigned nprobes; /* Number of probes to solve */ 46 int quit; 47 }; 48 #define ARGS_DEFAULT__ {1, 0} 49 static const struct args ARGS_DEFAULT = ARGS_DEFAULT__; 50 51 /* 52 * The system is a trilinear profile of the temperature at steady state, i.e. at 53 * each point of the system we can calculate the temperature analytically. Two 54 * forms are immersed in this temperature field: a super shape and a sphere 55 * included in the super shape. On the Monte Carlo side, the temperature is 56 * unknown everywhere except on the surface of the super shape whose 57 * temperature is defined from the aformentionned trilinear profile. We will 58 * estimate the temperature at the sphere boundary at a probe points. We 59 * should find by Monte Carlo the temperature of the trilinear profile at the 60 * position of the probe. It's the test. 61 * 62 * /\ <-- T(x,y,z) 63 * ___/ \___ 64 * T(z) \ __ / 65 * | T(y) \ / \ / 66 * |/ T=? *__/ \ 67 * o--- T(x) /_ __ _\ 68 * \/ \/ 69 */ 70 71 #define FILENAME_SSHAPE "sshape.stl" 72 #define FILENAME_SPHERE "sphere.stl" 73 #define FILENAME_PROBES "probes.txt" 74 #define FILENAME_SCENE "scene.txt" 75 76 /* Temperature at the upper bound of the X, Y and Z axis. The temperature at the 77 * lower bound is implicitly 0 K */ 78 #define TX 333.0 /* [K] */ 79 #define TY 432.0 /* [K] */ 80 #define TZ 579.0 /* [K] */ 81 82 /* Probe position */ 83 #define PX 1.0 84 #define PY 0.0 85 #define PZ 0.0 86 87 /* Commands to calculate 1 or N probes */ 88 #define COMMAND1 "stardis -V3 -P "STR(PX)","STR(PY)","STR(PZ)" -M "FILENAME_SCENE 89 #define COMMANDN "stardis -V3 -L "FILENAME_PROBES" -M "FILENAME_SCENE 90 #define COMMANDN_MPI "mpirun -n 2 "COMMANDN 91 92 /******************************************************************************* 93 * Helper functions 94 ******************************************************************************/ 95 static double 96 rand_canonic(void) 97 { 98 return (double)rand() / (double)((long)RAND_MAX - 1); 99 } 100 101 static int 102 is_mpi_enabled(void) 103 { 104 char buf[128] = {0}; 105 FILE* sdis = NULL; 106 FILE* grep = NULL; 107 108 if(!(sdis = popen("stardis -v", "r"))) 109 perror("stardis"), exit(errno); 110 111 if(!(grep = popen("grep -qe \" MPI is enabled\"", "w"))) 112 perror("grep"), exit(errno); 113 114 while(fgets(buf, sizeof(buf), sdis) && !feof(sdis)) { 115 if(fputs(buf, grep) == EOF) abort(); 116 } 117 118 if(ferror(sdis)) errno=EIO, perror("stardis"), exit(errno); 119 if(pclose(sdis) == -1) perror("stardis"), exit(errno); 120 return !pclose(grep); 121 } 122 123 static void 124 sample_unit_sphere(double p[3]) 125 { 126 const double phi = rand_canonic() * 2*PI; 127 const double v = rand_canonic(); 128 const double cos_theta = 1 - 2*v; 129 const double sin_theta = 2 * sqrt(v*(1-v)); 130 p[0] = cos(phi) * sin_theta; 131 p[1] = sin(phi) * sin_theta; 132 p[2] = cos_theta; 133 } 134 135 static void 136 aabb_update(struct aabb* aabb, const struct s3dut_mesh_data* mesh) 137 { 138 size_t ivertex = 0; 139 ASSERT(aabb); 140 141 FOR_EACH(ivertex, 0, mesh->nvertices) { 142 const double* vertex = mesh->positions + ivertex*3; 143 aabb->lower[0] = MMIN(aabb->lower[0], vertex[0]); 144 aabb->lower[1] = MMIN(aabb->lower[1], vertex[1]); 145 aabb->lower[2] = MMIN(aabb->lower[2], vertex[2]); 146 aabb->upper[0] = MMAX(aabb->upper[0], vertex[0]); 147 aabb->upper[1] = MMAX(aabb->upper[1], vertex[1]); 148 aabb->upper[2] = MMAX(aabb->upper[2], vertex[2]); 149 } 150 } 151 152 static void 153 setup_sshape(FILE* stream, struct aabb* aabb) 154 { 155 struct s3dut_mesh* sshape = NULL; 156 struct s3dut_mesh_data sshape_data; 157 struct s3dut_super_formula f0 = S3DUT_SUPER_FORMULA_NULL; 158 struct s3dut_super_formula f1 = S3DUT_SUPER_FORMULA_NULL; 159 const double radius = 2; 160 const unsigned nslices = 256; 161 162 f0.A = 1.5; f0.B = 1; f0.M = 11.0; f0.N0 = 1; f0.N1 = 1; f0.N2 = 2.0; 163 f1.A = 1.0; f1.B = 2; f1.M = 3.6; f1.N0 = 1; f1.N1 = 2; f1.N2 = 0.7; 164 S3DUT(create_super_shape(NULL, &f0, &f1, radius, nslices, nslices/2, &sshape)); 165 S3DUT(mesh_get_data(sshape, &sshape_data)); 166 167 aabb_update(aabb, &sshape_data); 168 sadist_write_stl(stream, sshape_data.positions, sshape_data.nvertices, 169 sshape_data.indices, sshape_data.nprimitives); 170 171 S3DUT(mesh_ref_put(sshape)); 172 } 173 174 static void 175 setup_sphere(FILE* stream, struct aabb* aabb) 176 { 177 struct s3dut_mesh* sphere = NULL; 178 struct s3dut_mesh_data sphere_data; 179 const double radius = 1; 180 const unsigned nslices = 128; 181 182 S3DUT(create_sphere(NULL, radius, nslices, nslices/2, &sphere)); 183 S3DUT(mesh_get_data(sphere, &sphere_data)); 184 185 aabb_update(aabb, &sphere_data); 186 sadist_write_stl(stream, sphere_data.positions, sphere_data.nvertices, 187 sphere_data.indices, sphere_data.nprimitives); 188 189 S3DUT(mesh_ref_put(sphere)); 190 } 191 192 static void 193 setup_scene 194 (FILE* fp, 195 const char* sshape, 196 const char* sphere, 197 const struct aabb* aabb, 198 struct sadist_trilinear_profile* profile) 199 { 200 double low, upp; 201 ASSERT(sshape && sphere && aabb && profile); 202 203 fprintf(fp, "PROGRAM trilinear_profile libsadist_trilinear_profile.so\n"); 204 fprintf(fp, "SOLID SuperShape 25 7500 500 0.05 0 UNKNOWN 0 BACK %s FRONT %s\n", 205 sshape, sphere); 206 fprintf(fp, "SOLID Sphere 25 7500 500 0.05 0 UNKNOWN 0 BACK %s\n", sphere); 207 208 low = MMIN(MMIN(aabb->lower[0], aabb->lower[1]), aabb->lower[2]); 209 upp = MMAX(MMAX(aabb->upper[0], aabb->upper[1]), aabb->upper[2]); 210 fprintf(fp, "T_BOUNDARY_FOR_SOLID_PROG Dirichlet trilinear_profile %s", sshape); 211 fprintf(fp, " PROG_PARAMS -b %g,%g -t %g,%g,%g\n", low, upp, TX, TY, TZ); 212 213 d3_splat(profile->lower, low); 214 d3_splat(profile->upper, upp); 215 d2(profile->a, 0, TX); 216 d2(profile->b, 0, TY); 217 d2(profile->c, 0, TZ); 218 } 219 220 static void 221 usage(const char* name, FILE* stream) 222 { 223 ASSERT(name && stream); 224 fprintf(stream, "usage: %s [-h] [-p probes_count]\n", name); 225 } 226 227 static int 228 args_init(struct args* args, int argc, char** argv) 229 { 230 int opt = 0; 231 int err = 0; 232 233 ASSERT(args); 234 while((opt = getopt(argc, argv, "hp:")) != -1) { 235 switch(opt) { 236 case 'h': 237 usage(argv[0], stdout); 238 args->quit = 1; 239 break; 240 case 'p': 241 err = (cstr_to_uint(optarg, &args->nprobes) != RES_OK); 242 if(!err && args->nprobes == 0) err = 1; 243 break; 244 default: err = 1; break; 245 } 246 if(err) { 247 if(optarg) { 248 fprintf(stderr, "%s: invalid option argument `%s' -- `%c'\n", 249 argv[0], optarg, opt); 250 } 251 goto error; 252 } 253 } 254 255 exit: 256 return err; 257 error: 258 usage(argv[0], stderr); 259 goto exit; 260 } 261 262 static int 263 init(struct sadist_trilinear_profile* profile) 264 { 265 struct aabb aabb = AABB_NULL; 266 FILE* fp_sshape = NULL; 267 FILE* fp_sphere = NULL; 268 FILE* fp_scene = NULL; 269 int err = 0; 270 271 if((fp_sshape = fopen(FILENAME_SSHAPE, "w")) == NULL) { 272 fprintf(stderr, "Error opening the `"FILENAME_SSHAPE"' file -- %s\n", 273 strerror(errno)); 274 err = errno; 275 goto error; 276 } 277 if((fp_sphere = fopen(FILENAME_SPHERE, "w")) == NULL) { 278 fprintf(stderr, "Error opening the `"FILENAME_SPHERE"' file -- %s\n", 279 strerror(errno)); 280 err = errno; 281 goto error; 282 } 283 if((fp_scene = fopen(FILENAME_SCENE, "w")) == NULL) { 284 fprintf(stderr, "Error opening the `"FILENAME_SCENE"' file -- %s\n", 285 strerror(errno)); 286 err = errno; 287 goto error; 288 } 289 290 setup_sshape(fp_sshape, &aabb); 291 setup_sphere(fp_sphere, &aabb); 292 setup_scene(fp_scene, FILENAME_SSHAPE, FILENAME_SPHERE, &aabb, profile); 293 294 exit: 295 if(fp_sshape && fclose(fp_sshape)) { perror("fclose"); if(!err) err = errno; } 296 if(fp_sphere && fclose(fp_sphere)) { perror("fclose"); if(!err) err = errno; } 297 if(fp_scene && fclose(fp_scene)) { perror("fclose"); if(!err) err = errno; } 298 return err; 299 error: 300 goto exit; 301 } 302 303 static int 304 init_probe_list(double* probes, const size_t nprobes) 305 { 306 FILE* fp_probes = NULL; 307 size_t i = 0; 308 int err = 0; 309 ASSERT(probes && nprobes); 310 311 if((fp_probes = fopen(FILENAME_PROBES, "w")) == NULL) { 312 fprintf(stderr, "Error opening `"FILENAME_PROBES"' file -- %s\n", 313 strerror(errno)); 314 err = errno; 315 goto error; 316 } 317 318 FOR_EACH(i, 0, nprobes) { 319 sample_unit_sphere(probes+i*3); 320 if(fprintf(fp_probes, "%g,%g,%g\n", SPLIT3(probes+i*3)) < 0) { 321 fprintf(stderr, "Error writing probes -- %s\n", strerror(errno)); 322 err = errno; 323 goto error; 324 } 325 } 326 327 exit: 328 if(fp_probes && fclose(fp_probes)) { perror("fclose"); if(!err) err = errno; } 329 return err; 330 error: 331 goto exit; 332 } 333 334 static int 335 run1(const struct sadist_trilinear_profile* profile) 336 { 337 const double P[3] = {PX, PY, PZ}; 338 FILE* output = NULL; 339 double ref = 0; 340 double E = 0; 341 double SE = 0; 342 int n = 0; 343 int err = 0; 344 int status = 0; 345 346 printf(COMMAND1"\n"); 347 348 if(!(output = popen(COMMAND1, "r"))) { 349 fprintf(stderr, "Error executing stardis -- %s\n", strerror(errno)); 350 fprintf(stderr, "\t"COMMAND1"\n"); 351 err = errno; 352 goto error; 353 } 354 355 if((n = fscanf(output, "%lf %lf %*d %*d", &E, &SE)), n != 2 && n != EOF) { 356 fprintf(stderr, "Error reading the output stream -- %s\n", strerror(errno)); 357 err = errno; 358 goto error; 359 } 360 361 /* Check command exit status */ 362 if((status=pclose(output), output=NULL, status)) { 363 if(status == -1) err = errno; 364 else if(WIFEXITED(status)) err = WEXITSTATUS(status); 365 else if(WIFSIGNALED(status)) err = WTERMSIG(status); 366 else if(WIFSTOPPED(status)) err = WSTOPSIG(status); 367 else FATAL("Unreachable code.\n"); 368 goto error; 369 } 370 371 ref = sadist_trilinear_profile_temperature(profile, P); 372 printf("T = %g ~ %g +/- %g\n", ref, E, SE); 373 if(!eq_eps(ref, E, SE*3)) { 374 err = 1; 375 goto error; 376 } 377 378 exit: 379 if(output) pclose(output); 380 return err; 381 error: 382 goto exit; 383 } 384 385 static int 386 runN(const struct sadist_trilinear_profile* profile, const size_t nprobes) 387 { 388 const char* command = NULL; 389 double* probes = NULL; /* Positions */ 390 double* results = NULL; /* Expected values and standard deviations */ 391 FILE* output = NULL; 392 size_t i = 0; 393 int err = 0; 394 int status = 0; 395 ASSERT(profile && nprobes); 396 397 if(!(probes = mem_calloc(nprobes, sizeof(double)*3))) { 398 err = errno = ENOMEM; 399 perror("mem_calloc"); 400 goto error; 401 } 402 403 if(!(results = mem_calloc(nprobes, sizeof(double)*2))) { 404 err = errno = ENOMEM; 405 perror("mem_calloc"); 406 goto error; 407 } 408 409 if((err = init_probe_list(probes, nprobes))) goto error; 410 411 if(is_mpi_enabled()) { 412 command = COMMANDN_MPI; 413 } else { 414 command = COMMANDN; 415 } 416 printf("%s\n", command); 417 418 if(!(output = popen(command, "r"))) { 419 fprintf(stderr, "Error executing stardis -- %s\n", strerror(errno)); 420 fprintf(stderr, "\t%s\n", command); 421 err = errno; 422 goto error; 423 } 424 425 /* Fetch the result */ 426 FOR_EACH(i, 0, nprobes) { 427 double* E = results + i*2 + 0; 428 double* SE = results + i*2 + 1; 429 int n = 0; 430 431 if((n = fscanf(output, "%lf %lf %*d %*d", E, SE)), n != 2 && n != EOF) { 432 perror("fscanf"); 433 err = errno; 434 goto error; 435 } 436 } 437 438 /* Check command exit status */ 439 if((status=pclose(output)), output = NULL, status) { 440 if(status == -1) err = errno; 441 else if(WIFEXITED(status)) err = WEXITSTATUS(status); 442 else if(WIFSIGNALED(status)) err = WTERMSIG(status); 443 else if(WIFSTOPPED(status)) err = WSTOPSIG(status); 444 else FATAL("Unreachable code.\n"); 445 goto error; 446 } 447 output = NULL; 448 449 /* Validate the calculations */ 450 FOR_EACH(i, 0, nprobes) { 451 const double* P = probes + i*3; 452 const double E = results[i*2 + 0]; 453 const double SE = results[i*2 + 1]; 454 const double ref = sadist_trilinear_profile_temperature(profile, P); 455 456 printf("T = %g ~ %g +/- %g\n", ref, E, SE); 457 if(!eq_eps(ref, E, SE*3)) { 458 err = 1; 459 goto error; 460 } 461 } 462 463 exit: 464 if(probes) mem_rm(probes); 465 if(results) mem_rm(results); 466 if(output) pclose(output); 467 return err; 468 error: 469 goto exit; 470 } 471 472 /******************************************************************************* 473 * The test 474 ******************************************************************************/ 475 int 476 main(int argc, char** argv) 477 { 478 struct args args = ARGS_DEFAULT; 479 struct sadist_trilinear_profile profile = SADIST_TRILINEAR_PROFILE_NULL; 480 int err = 0; 481 482 if((err = args_init(&args, argc, argv))) goto error; 483 if(args.quit) goto exit; 484 485 if((err = init(&profile))) goto error; 486 487 if(args.nprobes == 1) { 488 if((err = run1(&profile))) goto error; 489 } else { 490 if((err = runN(&profile, args.nprobes))) goto error; 491 } 492 493 exit: 494 if(mem_allocated_size() != 0) { 495 fprintf(stderr, "Memory leaks: %lu bytes\n", mem_allocated_size()); 496 if(!err) err = -1; 497 } 498 return err; 499 error: 500 goto exit; 501 }