sadist_unsteady.c (6918B)
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 <star/s3dut.h> 21 22 #include <rsys/math.h> 23 #include <rsys/mem_allocator.h> 24 25 #include <errno.h> 26 #include <stdio.h> /* popen */ 27 #include <string.h> /* strerror */ 28 #include <wait.h> /* WIFEXITED, WEXITSTATUS */ 29 30 /* 31 * The system is an unsteady-state temperature profile, meaning that at any 32 * point, at any time, we can analytically calculate the temperature. We 33 * immerse in this temperature field a supershape representing a solid in which 34 * we want to evaluate the temperature by Monte Carlo at a given position and 35 * observation time. On the Monte Carlo side, the temperature of the supershape 36 * is unknown. Only the boundary temperature is fixed at the temperature of the 37 * unsteady trilinear profile mentioned above. Monte Carlo would then have to 38 * find the temperature defined by the unsteady profile. 39 * 40 * T(z) /\ <-- T(x,y,z,t) 41 * | T(y) ___/ \___ 42 * |/ \ . T=? / 43 * o--- T(x) /_ __ _\ 44 * \/ \/ 45 */ 46 47 #define FILENAME_SSHAPE "sshape.stl" 48 #define FILENAME_SCENE "scene.txt" 49 50 #define LAMBDA 0.1 /* [W/(m.K)] */ 51 #define RHO 25.0 /* [kg/m^3] */ 52 #define CP 2.0 /* [J/K/kg)] */ 53 #define DELTA 1.0/20.0 /* [m/fp_to_meter] */ 54 55 /* Parameters of the unsteady profile */ 56 #define UPROF_A 0.0 57 #define UPROF_B1 10.0 58 #define UPROF_B2 1000.0 59 #define UPROF_K (PI/4.0) 60 61 /* Probe position */ 62 #define PX 0.2 /* [m/fp_to_meter] */ 63 #define PY 0.3 /* [m/fp_to_meter] */ 64 #define PZ 0.4 /* [m/fp_to_meter] */ 65 66 /* Observation time */ 67 #define TIME 5 /* [s] */ 68 69 #define NREALISATIONS 100000 70 71 #define COMMAND "stardis -V3 -p "STR(PX)","STR(PY)","STR(PZ)","STR(TIME) \ 72 " -n "STR(NREALISATIONS)" -M "FILENAME_SCENE " -I -INF" 73 #define COMMAND_DSPHERE COMMAND" -a dsphere" 74 #define COMMAND_WOS COMMAND" -a wos" 75 76 /******************************************************************************* 77 * Helper functions 78 ******************************************************************************/ 79 static void 80 setup_sshape(FILE* stream) 81 { 82 struct s3dut_mesh* sshape = NULL; 83 struct s3dut_mesh_data sshape_data; 84 struct s3dut_super_formula f0 = S3DUT_SUPER_FORMULA_NULL; 85 struct s3dut_super_formula f1 = S3DUT_SUPER_FORMULA_NULL; 86 const double radius = 1; 87 const unsigned nslices = 256; 88 89 f0.A = 1.5; f0.B = 1; f0.M = 11.0; f0.N0 = 1; f0.N1 = 1; f0.N2 = 2.0; 90 f1.A = 1.0; f1.B = 2; f1.M = 3.6; f1.N0 = 1; f1.N1 = 2; f1.N2 = 0.7; 91 S3DUT(create_super_shape(NULL, &f0, &f1, radius, nslices, nslices/2, &sshape)); 92 S3DUT(mesh_get_data(sshape, &sshape_data)); 93 94 sadist_write_stl(stream, sshape_data.positions, sshape_data.nvertices, 95 sshape_data.indices, sshape_data.nprimitives); 96 97 S3DUT(mesh_ref_put(sshape)); 98 } 99 100 static void 101 setup_scene 102 (FILE* fp, 103 const char* sshape, 104 struct sadist_unsteady_profile* profile) 105 { 106 ASSERT(sshape && profile); 107 108 fprintf(fp, "PROGRAM unsteady_profile libsadist_unsteady_profile.so\n"); 109 fprintf(fp, "SOLID SuperShape %g %g %g %g 0 UNKNOWN 0 BACK %s\n", 110 LAMBDA, RHO, CP, DELTA, sshape); 111 112 fprintf(fp, "T_BOUNDARY_FOR_SOLID_PROG Dirichlet unsteady_profile %s", sshape); 113 fprintf(fp, " PROG_PARAMS -p %g,%g,%g -k %g,%g,%g -m %g,%g,%g\n", 114 UPROF_A, UPROF_B1, UPROF_B2, UPROF_K, UPROF_K, UPROF_K, LAMBDA, RHO, CP); 115 116 profile->A = UPROF_A; 117 profile->B1 = UPROF_B1; 118 profile->B2 = UPROF_B2; 119 profile->kx = UPROF_K; 120 profile->ky = UPROF_K; 121 profile->kz = UPROF_K; 122 profile->lambda = LAMBDA; 123 profile->rho = RHO; 124 profile->cp = CP; 125 } 126 127 static int 128 init(struct sadist_unsteady_profile* profile) 129 { 130 FILE* fp_sshape = NULL; 131 FILE* fp_scene = NULL; 132 int err = 0; 133 134 if((fp_sshape = fopen(FILENAME_SSHAPE, "w")) == NULL) { 135 fprintf(stderr, "Error opening the `"FILENAME_SSHAPE"' file -- %s\n", 136 strerror(errno)); 137 err = errno; 138 goto error; 139 } 140 141 if((fp_scene = fopen(FILENAME_SCENE, "w")) == NULL) { 142 fprintf(stderr, "Error opening the `"FILENAME_SCENE"' file -- %s\n", 143 strerror(errno)); 144 err = errno; 145 goto error; 146 } 147 148 setup_sshape(fp_sshape); 149 setup_scene(fp_scene, FILENAME_SSHAPE, profile); 150 151 exit: 152 if(fp_sshape && fclose(fp_sshape)) { perror("fclose"); if(!err) err = errno; } 153 if(fp_scene && fclose(fp_scene)) { perror("fclose"); if(!err) err = errno; } 154 return err; 155 error: 156 goto exit; 157 } 158 159 static int 160 run(struct sadist_unsteady_profile* profile, const char* command) 161 { 162 const double P[3] = {PX, PY, PZ}; 163 FILE* output = NULL; 164 double ref = 0; 165 double E = 0; 166 double SE = 0; 167 int n = 0; 168 int err = 0; 169 int status = 0; 170 171 printf("%s\n", command); 172 173 if(!(output = popen(command, "r"))) { 174 fprintf(stderr, "Error executing stardis -- %s\n", strerror(errno)); 175 fprintf(stderr, "\t"COMMAND"\n"); 176 err = errno; 177 goto error; 178 } 179 180 if((n = fscanf(output, "%lf %lf %*d %*d", &E, &SE)), n != 2 && n != EOF) { 181 fprintf(stderr, "Error reading the output stream -- %s\n", strerror(errno)); 182 err = errno; 183 goto error; 184 } 185 186 /* Check command exit status */ 187 if((status=pclose(output), output=NULL, status)) { 188 if(status == -1) err = errno; 189 else if(WIFEXITED(status)) err = WEXITSTATUS(status); 190 else if(WIFSIGNALED(status)) err = WTERMSIG(status); 191 else if(WIFSTOPPED(status)) err = WSTOPSIG(status); 192 else FATAL("Unreachable code.\n"); 193 goto error; 194 } 195 196 ref = sadist_unsteady_profile_temperature(profile, P, TIME); 197 printf("T = %g ~ %g +/- %g\n", ref, E, SE); 198 if(!eq_eps(ref, E, SE*3)) { 199 err = 1; 200 goto error; 201 } 202 203 exit: 204 if(output) pclose(output); 205 return err; 206 error: 207 goto exit; 208 } 209 210 /******************************************************************************* 211 * The test 212 ******************************************************************************/ 213 int 214 main(int argc, char** argv) 215 { 216 struct sadist_unsteady_profile profile = SADIST_UNSTEADY_PROFILE_NULL; 217 int err = 0; 218 (void)argc, (void)argv; 219 220 if((err = init(&profile))) goto error; 221 if((err = run(&profile, COMMAND_DSPHERE))) goto error; 222 if((err = run(&profile, COMMAND_WOS))) goto error; 223 224 exit: 225 if(mem_allocated_size() != 0) { 226 fprintf(stderr, "Memory leaks: %lu bytes\n", mem_allocated_size()); 227 if(!err) err = -1; 228 } 229 return err; 230 error: 231 goto exit; 232 }