sadist_custom_conductive_path.c (8472B)
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/double2.h> 25 #include <rsys/double3.h> 26 #include <rsys/math.h> 27 28 #include <errno.h> 29 #include <float.h> 30 #include <getopt.h> 31 #include <string.h> /* strerror */ 32 #include <wait.h> /* WIFEXITED, WEXITSTATUS */ 33 34 /* 35 * The system is a trilinear profile of the temperature at steady state, i.e. at 36 * each point of the system we can calculate the temperature analytically. Two 37 * forms are immersed in this temperature field: a super shape and a sphere 38 * included in the super shape. On the Monte Carlo side, the temperature is 39 * unknown everywhere except on the surface of the super shape whose 40 * temperature is defined from the aformentionned trilinear profile. 41 * 42 * We will estimate the temperature at the position of a probe in solids by 43 * providing a user-side function to sample the conductive path in the sphere. 44 * We should find the temperature of the trilinear profile at the probe position 45 * by Monte Carlo, independently of this coupling with an external path sampling 46 * routine. 47 * 48 * 49 * /\ <-- T(x,y,z) 50 * ___/ \___ 51 * T(z) \ __ / 52 * | T(y) T=?/. \ / 53 * |/ / \__/ \ 54 * o--- T(x) /_ __ _\ 55 * \/ \/ 56 */ 57 58 #define FILENAME_SSHAPE "sshape.stl" 59 #define FILENAME_SPHERE "sphere.stl" 60 #define FILENAME_SCENE "scene.txt" 61 62 /* Temperature at the upper bound of the X, Y and Z axis. The temperature at the 63 * lower bound is implicitly 0 K */ 64 #define TX 333.0 /* [K] */ 65 #define TY 432.0 /* [K] */ 66 #define TZ 579.0 /* [K] */ 67 68 /* Probe position */ 69 #define PX 0.125 70 #define PY 0.250 71 #define PZ 0.375 72 73 /* Commands to calculate probe temperature */ 74 #define COMMAND "stardis -V3 -p "STR(PX)","STR(PY)","STR(PZ)" -M "FILENAME_SCENE 75 76 /* Axis Aligned Bounding Box */ 77 struct aabb { 78 double lower[3]; 79 double upper[3]; 80 }; 81 #define AABB_NULL__ {{DBL_MAX,DBL_MAX,DBL_MAX}, {-DBL_MAX,-DBL_MAX,-DBL_MAX}} 82 static const struct aabb AABB_NULL = AABB_NULL__; 83 84 85 /******************************************************************************* 86 * Helper functions 87 ******************************************************************************/ 88 static void 89 aabb_update(struct aabb* aabb, const struct s3dut_mesh_data* mesh) 90 { 91 size_t ivertex = 0; 92 ASSERT(aabb); 93 94 FOR_EACH(ivertex, 0, mesh->nvertices) { 95 const double* vertex = mesh->positions + ivertex*3; 96 aabb->lower[0] = MMIN(aabb->lower[0], vertex[0]); 97 aabb->lower[1] = MMIN(aabb->lower[1], vertex[1]); 98 aabb->lower[2] = MMIN(aabb->lower[2], vertex[2]); 99 aabb->upper[0] = MMAX(aabb->upper[0], vertex[0]); 100 aabb->upper[1] = MMAX(aabb->upper[1], vertex[1]); 101 aabb->upper[2] = MMAX(aabb->upper[2], vertex[2]); 102 } 103 } 104 105 static void 106 setup_sshape(FILE* stream, struct aabb* aabb) 107 { 108 struct s3dut_mesh* sshape = NULL; 109 struct s3dut_mesh_data sshape_data; 110 struct s3dut_super_formula f0 = S3DUT_SUPER_FORMULA_NULL; 111 struct s3dut_super_formula f1 = S3DUT_SUPER_FORMULA_NULL; 112 const double radius = 2; 113 const unsigned nslices = 256; 114 115 f0.A = 1.5; f0.B = 1; f0.M = 11.0; f0.N0 = 1; f0.N1 = 1; f0.N2 = 2.0; 116 f1.A = 1.0; f1.B = 2; f1.M = 3.6; f1.N0 = 1; f1.N1 = 2; f1.N2 = 0.7; 117 S3DUT(create_super_shape(NULL, &f0, &f1, radius, nslices, nslices/2, &sshape)); 118 S3DUT(mesh_get_data(sshape, &sshape_data)); 119 120 aabb_update(aabb, &sshape_data); 121 sadist_write_stl(stream, sshape_data.positions, sshape_data.nvertices, 122 sshape_data.indices, sshape_data.nprimitives); 123 124 S3DUT(mesh_ref_put(sshape)); 125 } 126 127 static void 128 setup_sphere(FILE* stream, struct aabb* aabb) 129 { 130 struct s3dut_mesh* sphere = NULL; 131 struct s3dut_mesh_data sphere_data; 132 const double radius = 1; 133 const unsigned nslices = 128; 134 135 S3DUT(create_sphere(NULL, radius, nslices, nslices/2, &sphere)); 136 S3DUT(mesh_get_data(sphere, &sphere_data)); 137 138 aabb_update(aabb, &sphere_data); 139 sadist_write_stl(stream, sphere_data.positions, sphere_data.nvertices, 140 sphere_data.indices, sphere_data.nprimitives); 141 142 S3DUT(mesh_ref_put(sphere)); 143 } 144 145 static void 146 setup_scene 147 (FILE* fp, 148 const char* sshape, 149 const char* sphere, 150 const struct aabb* aabb, 151 struct sadist_trilinear_profile* profile) 152 { 153 double low, upp; 154 ASSERT(sshape && sphere && aabb && profile); 155 156 fprintf(fp, "PROGRAM solid_sphere libsadist_solid_sphere.so\n"); 157 fprintf(fp, "PROGRAM trilinear_profile libsadist_trilinear_profile.so\n"); 158 fprintf(fp, "SOLID SuperShape 25 7500 500 0.05 0 UNKNOWN 0 BACK %s FRONT %s\n", 159 sshape, sphere); 160 fprintf(fp, "SOLID_PROG Sphere solid_sphere BACK %s ", sphere); 161 162 /* StL lambda rho cp radius */ 163 fprintf(fp, "PROG_PARAMS -m %s -l25 -r7500 -c500 -R1\n", FILENAME_SPHERE); 164 165 low = MMIN(MMIN(aabb->lower[0], aabb->lower[1]), aabb->lower[2]); 166 upp = MMAX(MMAX(aabb->upper[0], aabb->upper[1]), aabb->upper[2]); 167 fprintf(fp, "T_BOUNDARY_FOR_SOLID_PROG Dirichlet trilinear_profile %s", sshape); 168 fprintf(fp, " PROG_PARAMS -b %g,%g -t %g,%g,%g\n", low, upp, TX, TY, TZ); 169 170 d3_splat(profile->lower, low); 171 d3_splat(profile->upper, upp); 172 d2(profile->a, 0, TX); 173 d2(profile->b, 0, TY); 174 d2(profile->c, 0, TZ); 175 } 176 177 static int 178 init(struct sadist_trilinear_profile* profile) 179 { 180 struct aabb aabb = AABB_NULL; 181 FILE* fp_sshape = NULL; 182 FILE* fp_sphere = NULL; 183 FILE* fp_scene = NULL; 184 int err = 0; 185 186 #define FOPEN(Fp, Filename) \ 187 if(((Fp) = fopen((Filename), "w")) == NULL) { \ 188 fprintf(stderr, "Error opening `%s' -- file %s\n", \ 189 (Filename), strerror(errno)); \ 190 err = errno; \ 191 goto error; \ 192 } (void)0 193 FOPEN(fp_sshape, FILENAME_SSHAPE); 194 FOPEN(fp_sphere, FILENAME_SPHERE); 195 FOPEN(fp_scene, FILENAME_SCENE); 196 #undef FOPEN 197 198 setup_sshape(fp_sshape, &aabb); 199 setup_sphere(fp_sphere, &aabb); 200 setup_scene(fp_scene, FILENAME_SSHAPE, FILENAME_SPHERE, &aabb, profile); 201 202 exit: 203 #define FCLOSE(Fp) \ 204 if((Fp) && fclose(Fp)) { perror("fclose"); if(!err) err = errno; } (void)0 205 FCLOSE(fp_sshape); 206 FCLOSE(fp_sphere); 207 FCLOSE(fp_scene); 208 #undef FCLOSE 209 return err; 210 error: 211 goto exit; 212 } 213 214 static int 215 run(const struct sadist_trilinear_profile* profile) 216 { 217 const double P[3] = {PX, PY, PZ}; 218 219 FILE* output = NULL; 220 double E = 0; 221 double SE = 0; 222 double ref = 0; 223 224 int n = 0; 225 int err = 0; 226 int status = 0; 227 228 printf("%s\n", COMMAND); 229 230 if(!(output = popen(COMMAND, "r"))) { 231 fprintf(stderr, "Error executing stardis -- %s\n", strerror(errno)); 232 err = errno; 233 goto error; 234 } 235 236 if((n = fscanf(output, "%lf %lf %*d %*d", &E, &SE)), n != 2 && n != EOF) { 237 fprintf(stderr, "Error reading the output stream -- %s\n", strerror(errno)); 238 err = errno; 239 goto error; 240 } 241 242 /* Check command exit status */ 243 if((status=pclose(output)), output=NULL, status) { 244 if(status == -1) err = errno; 245 else if(WIFEXITED(status)) err = WEXITSTATUS(status); 246 else if(WIFSIGNALED(status)) err = WTERMSIG(status); 247 else if(WIFSTOPPED(status)) err = WSTOPSIG(status); 248 else FATAL("Unreachable code.\n"); 249 goto error; 250 } 251 252 ref = sadist_trilinear_profile_temperature(profile, P); 253 printf("T = %g ~ %g +/- %g\n", ref, E, SE); 254 if(!eq_eps(ref, E, SE*3)) { 255 err = 1; 256 goto error; 257 } 258 259 exit: 260 if(output) pclose(output); 261 return err; 262 error: 263 goto exit; 264 } 265 266 /******************************************************************************* 267 * The test 268 ******************************************************************************/ 269 int 270 main(int argc, char** argv) 271 { 272 struct sadist_trilinear_profile profile = SADIST_TRILINEAR_PROFILE_NULL; 273 int err; 274 (void)argc, (void)argv; 275 276 if((err = init(&profile))) goto error; 277 if((err = run(&profile))) goto error; 278 279 exit: 280 return err; 281 error: 282 goto exit; 283 }