sadist_external_flux.c (8435B)
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 <rsys/math.h> 21 22 #include <math.h> 23 #include <string.h> /* strerror */ 24 #include <wait.h> /* WIFEXITED, WEXITSTATUS */ 25 26 /* 27 * The system consists of 2 parallelepipeds: a vertical one called the wall, and 28 * a horizontal one representing the floor. The wall is a black body, while the 29 * floor is a perfectly reflective surface. The surrounding fluid has a fixed 30 * temperature and, finally, an external spherical source represents the sun. 31 * This test calculates the steady temperature at a position in the wall and 32 * compares it with the analytical solution given for a perfectly diffuse or 33 * specular ground. 34 * 35 * (-0.1,1500) 36 * +---+ External source 37 * | E=1 T_FLUID ## 38 * Probe x | _\ #### 39 * | | / / ## 40 * +---+ \__/ 41 * (0,500) 42 * 43 * (0,0) 44 * Y *--------E=0------------- - - - 45 * | | 46 * o--X *------------------------ - - - 47 * / (0,-1) 48 * Z 49 * 50 */ 51 #define FILENAME_GROUND "ground.stl" 52 #define FILENAME_WALL "wall.stl" 53 #define FILENAME_SCENE "scene.txt" 54 #define FILENAME_SCENE_PROG "scene2.txt" 55 56 #define SOURCE_ELEVATION 30.0 /* [degree] */ 57 #define SOURCE_DISTANCE 1.5e11 /* [m] */ 58 #define SOURCE_RADIUS 6.5991756e8 59 #define SOURCE_POWER 3.845e26 /* [W] */ 60 61 /* Probe position */ 62 #define PX -0.05 63 #define PY 1000.0 64 #define PZ 0.0 65 66 /* The Commands and the expected temperatures */ 67 #define COMMAND1 "stardis -i -V3 -p "STR(PX)","STR(PY)","STR(PZ)" -M "FILENAME_SCENE 68 #define COMMAND2 "stardis -i -V3 -p "STR(PX)","STR(PY)","STR(PZ)\ 69 " -M "FILENAME_SCENE_PROG" -n 100000" 70 #define T_REF1 375.88 /* [K] */ 71 #define T_REF2 417.77 /* [K] */ 72 73 static const double ground_vertices[] = { 74 0.0, -1.0, -1.0e6, 75 1.0e12, -1.0, -1.0e6, 76 0.0, 1.0, -1.0e6, 77 1.0e12, 1.0, -1.0e6, 78 0.0, -1.0, 1.0e6, 79 1.0e12, -1.0, 1.0e6, 80 0.0, 1.0, 1.0e6, 81 1.0e12, 1.0, 1.0e6 82 }; 83 static const size_t ground_nvertices = sizeof(ground_vertices)/(sizeof(double)*3); 84 85 static const double wall_vertices[] = { 86 -0.1, 500.0, -500.0, 87 0.0, 500.0, -500.0, 88 -0.1, 1500.0, -500.0, 89 0.0, 1500.0, -500.0, 90 -0.1, 500.0, 500.0, 91 0.0, 500.0, 500.0, 92 -0.1, 1500.0, 500.0, 93 0.0, 1500.0, 500.0 94 }; 95 static const size_t wall_nvertices = sizeof(wall_vertices)/(sizeof(double)*3); 96 97 /* Ground and the wall indices */ 98 static const size_t indices[] = { 99 0, 2, 1, 1, 2, 3, /* -z */ 100 0, 4, 2, 2, 4, 6, /* -x */ 101 4, 5, 6, 6, 5, 7, /* +z */ 102 3, 7, 5, 5, 1, 3, /* +x */ 103 2, 6, 7, 7, 3, 2, /* +y */ 104 0, 1, 5, 5, 4, 0, /* -y */ 105 }; 106 107 static const size_t ntriangles = sizeof(indices) / (sizeof(size_t)*3); 108 109 /******************************************************************************* 110 * Helper functions 111 ******************************************************************************/ 112 static void 113 setup_scene(FILE* fp) 114 { 115 const double elevation = MDEG2RAD(SOURCE_ELEVATION); 116 double pos[3]; 117 ASSERT(fp); 118 119 fprintf(fp, "SOLID ground 1.15 1700 800 5.0e-3 0 UNKNOWN 0 FRONT " 120 FILENAME_GROUND "\n"); 121 fprintf(fp, "SOLID wall 1.15 1700 800 5.0e-3 0 UNKNOWN 0 FRONT " 122 FILENAME_WALL "\n"); 123 fprintf(fp, "FLUID dummy 1 1 300 300 BACK "FILENAME_GROUND" BACK "FILENAME_WALL "\n"); 124 fprintf(fp, "SOLID_FLUID_CONNECTION fluid_ground 0 0 0 0 "FILENAME_GROUND "\n"); 125 fprintf(fp, "SOLID_FLUID_CONNECTION fluid_wall 0 1 0 10 "FILENAME_WALL "\n"); 126 fprintf(fp, "TRAD 0 0\n"); 127 128 pos[0] = cos(elevation) * SOURCE_DISTANCE; 129 pos[1] = sin(elevation) * SOURCE_DISTANCE; 130 pos[2] = 0; 131 fprintf(fp, "SPHERICAL_SOURCE 6.5991756e8 %f %f %f 3.845e26 0\n", 132 pos[0], pos[1], pos[2]); 133 } 134 135 static void 136 setup_scene_prog(FILE* fp) 137 { 138 const double elevation = MDEG2RAD(SOURCE_ELEVATION); 139 double pos[3]; 140 ASSERT(fp); 141 142 fprintf(fp, "PROGRAM source libsadist_spherical_source.so\n"); 143 fprintf(fp, "SOLID ground 1.15 1700 800 5.0e-3 0 UNKNOWN 0 FRONT " 144 FILENAME_GROUND "\n"); 145 fprintf(fp, "SOLID wall 1.15 1700 800 5.0e-3 0 UNKNOWN 0 FRONT " 146 FILENAME_WALL "\n"); 147 fprintf(fp, "FLUID dummy 1 1 300 300 BACK "FILENAME_GROUND" BACK "FILENAME_WALL "\n"); 148 fprintf(fp, "SOLID_FLUID_CONNECTION fluid_ground 0 0 1 0 "FILENAME_GROUND "\n"); 149 fprintf(fp, "SOLID_FLUID_CONNECTION fluid_wall 0 1 0 10 "FILENAME_WALL "\n"); 150 fprintf(fp, "TRAD 0 0\n"); 151 152 pos[0] = cos(elevation) * SOURCE_DISTANCE; 153 pos[1] = sin(elevation) * SOURCE_DISTANCE; 154 pos[2] = 0; 155 156 fprintf(fp, "SPHERICAL_SOURCE_PROG "STR(SOURCE_RADIUS)" source " 157 "PROG_PARAMS -d 0 -p %f,%f,%f -w "STR(SOURCE_POWER)"\n", 158 pos[0], pos[1], pos[2]); 159 } 160 161 static int 162 init(void) 163 { 164 FILE* fp_ground = NULL; 165 FILE* fp_wall = NULL; 166 FILE* fp_scene = NULL; 167 FILE* fp_scene_prog = NULL; 168 int err = 0; 169 170 if((fp_ground = fopen(FILENAME_GROUND, "w")) == NULL) { 171 fprintf(stderr, "Error opening the `"FILENAME_GROUND"' file -- %s\n", 172 strerror(errno)); 173 err = errno; 174 goto error; 175 } 176 if((fp_wall = fopen(FILENAME_WALL, "w")) == NULL) { 177 fprintf(stderr, "Error opening the `"FILENAME_WALL "' file -- %s\n", 178 strerror(errno)); 179 err = errno; 180 goto error; 181 } 182 if((fp_scene = fopen(FILENAME_SCENE, "w")) == NULL) { 183 fprintf(stderr, "Error opening the `"FILENAME_SCENE"' file -- %s\n", 184 strerror(errno)); 185 err = errno; 186 goto error; 187 } 188 if((fp_scene_prog = fopen(FILENAME_SCENE_PROG, "w")) == NULL) { 189 fprintf(stderr, "Error opening the `"FILENAME_SCENE_PROG"' file -- %s\n", 190 strerror(errno)); 191 err = errno; 192 goto error; 193 } 194 195 sadist_write_stl(fp_ground, ground_vertices, ground_nvertices, indices, ntriangles); 196 sadist_write_stl(fp_wall, wall_vertices, wall_nvertices, indices, ntriangles); 197 setup_scene(fp_scene); 198 setup_scene_prog(fp_scene_prog); 199 200 exit: 201 if(fp_ground && fclose(fp_ground)) { perror("fclose"); if(!err) err = errno; } 202 if(fp_wall && fclose(fp_wall)) { perror("fclose"); if(!err) err = errno; } 203 if(fp_scene && fclose(fp_scene)) { perror("fclose"); if(!err) err = errno; } 204 if(fp_scene_prog && fclose(fp_scene_prog)) { perror("fclose"); if(!err) err = errno; } 205 return err; 206 error: 207 goto exit; 208 } 209 210 static int 211 run(const char* command, const double ref) 212 { 213 FILE* output = NULL; 214 double E = 0; 215 double SE = 0; 216 int n = 0; 217 int err = 0; 218 int status = 0; 219 220 printf("%s\n", command); 221 222 if(!(output = popen(command, "r"))) { 223 fprintf(stderr, "Error executing stardis -- %s\n", strerror(errno)); 224 err = errno; 225 goto error; 226 } 227 228 if((n = fscanf(output, "%lf %lf %*d %*d", &E, &SE)), n != 2 && n != EOF) { 229 fprintf(stderr, "Error reading the output stream -- %s\n", strerror(errno)); 230 err = errno; 231 goto error; 232 } 233 234 /* Check command exit status */ 235 if((status=pclose(output)), output=NULL, status) { 236 if(status == -1) err = errno; 237 else if(WIFEXITED(status)) err = WEXITSTATUS(status); 238 else if(WIFSIGNALED(status)) err = WTERMSIG(status); 239 else if(WIFSTOPPED(status)) err = WSTOPSIG(status); 240 else FATAL("Unreachable code.\n"); 241 goto error; 242 } 243 244 printf("T = %g ~ %g +/- %g\n", ref, E, SE); 245 if(!eq_eps(ref, E, SE*3)) { 246 err = 1; 247 goto error; 248 } 249 250 exit: 251 if(output) pclose(output); 252 return err; 253 error: 254 goto exit; 255 } 256 257 /******************************************************************************* 258 * The test 259 ******************************************************************************/ 260 int 261 main(int argc, char** argv) 262 { 263 int err = 0; 264 (void)argc, (void)argv; 265 266 if((err = init())) goto error; 267 if((err = run(COMMAND1, T_REF1))) goto error; 268 if((err = run(COMMAND2, T_REF2))) goto error; 269 270 exit: 271 return err; 272 error: 273 goto exit; 274 }