sadist_flux_with_h.c (9378B)
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 #include <star/sstl.h> 22 23 #include <rsys/cstr.h> 24 #include <rsys/float3.h> 25 #include <rsys/math.h> 26 27 #include <errno.h> 28 #include <stdio.h> 29 #include <string.h> /* strerror */ 30 #include <wait.h> /* WIFEXITED, WEXITSTATUS */ 31 32 /* 33 * The configuration is a rectangle with conductivity lambda and unknown 34 * temperature. Its upper and lower boundaries are adiabatic, while its left and 35 * right boundaries have fixed fluxes phi1 and phi2, respectively. The left 36 * boundary also has convective exchange with the surrounding fluid, whose 37 * temperature is Text. In steady state, the temperature at a given position in 38 * the solid rectangle is: 39 * 40 * T(x) = phi2/lambda*x + (Text + (phi1 + phi2)/h) 41 * 42 * with h the convective coefficient on the left boundary 43 * 44 * Text /////// (1,1) 45 * +---------+ 46 * | | 47 * h _\ |--> <--| 48 * / / phi1-> <-phi2 49 * \__/ |--> <--| 50 * | | 51 * +---------+ 52 * Y (0,0) /////// 53 * |__ X 54 * / 55 * Z 56 * 57 * This test is performed with different system descriptions in order to verify 58 * how Stardis processes them. Interfaces with imposed flows are defined either 59 * as boundary conditions or as fluid-solid connections with a surrounding fluid 60 * whose temperature is known. The objective is to verify that Stardis correctly 61 * processes these different descriptions. 62 */ 63 64 #define FILENAME1 "scene.txt" 65 #define FILENAME2 "scene2.txt" 66 67 /* Geometry */ 68 #define CUBOID_LEN_X 0.2 /* [m] */ 69 #define CUBOID_LEN_Y 0.5 /* [m] */ 70 #define CUBOID_LEN_Z 1.0 /* [m] */ 71 72 /* Physical properties */ 73 #define LAMBDA 25.0 /* [W.m^-1.K^-1] */ 74 #define RHO 7500.0 /* [kg.m^-3] */ 75 #define CP 500.0 /* [J.K^-1.kg^-1] */ 76 #define Text 373.15 /* [K] */ 77 #define PHI1 1000.0 /* [W.m^-2] */ 78 #define PHI2 5000.0 /* [W.m^-2] */ 79 #define H 10.0 /* [W.K^-1.m^-2] */ 80 81 /* Probe position */ 82 #define PX 0.1 83 #define PY 0.25 84 #define PZ 0.0 85 86 /* The commands to be executed */ 87 #define COMMAND "stardis -a wos -V3 -p"STR(PX)","STR(PY)","STR(PZ)" -n1000 -M" 88 89 enum facet_flag { 90 LX = BIT(0), /* Lower X */ 91 UX = BIT(1), /* Upper X */ 92 LY = BIT(2), /* Lower Y */ 93 UY = BIT(3), /* Upper Y */ 94 LZ = BIT(4), /* Lower Z */ 95 UZ = BIT(5), /* Upper Z */ 96 X = LX|UX, /* Lower and upper X */ 97 Y = LY|UY, /* Lower and upper Y */ 98 Z = LZ|UZ /* Lower and upper Z */ 99 }; 100 101 /******************************************************************************* 102 * Helper functions 103 ******************************************************************************/ 104 static enum facet_flag 105 get_facet_flag(const struct sstl_facet* facet) 106 { 107 ASSERT(facet); 108 109 if(facet->vertices[0][0] == facet->vertices[1][0] 110 && facet->vertices[0][0] == facet->vertices[2][0]) { 111 return facet->vertices[0][0] > 0 ? UX : LX; 112 } 113 114 if(facet->vertices[0][1] == facet->vertices[1][1] 115 && facet->vertices[0][1] == facet->vertices[2][1]) { 116 return facet->vertices[0][1] > 0 ? UY : LY; 117 } 118 119 if(facet->vertices[0][2] == facet->vertices[1][2] 120 && facet->vertices[0][2] == facet->vertices[2][2]) { 121 return facet->vertices[0][2] > 0 ? UY : LY; 122 } 123 124 FATAL("Unreachable code\n"); 125 return 0; 126 } 127 128 static res_T 129 write_faces 130 (FILE* fp, 131 const struct s3dut_mesh* mesh, 132 const unsigned int facet_mask) 133 { 134 struct sstl_writer_create_args args = SSTL_WRITER_CREATE_ARGS_DEFAULT__; 135 struct sstl_writer* writer = NULL; 136 struct s3dut_mesh_data data; 137 size_t i = 0; 138 res_T res = RES_OK; 139 ASSERT(fp && mesh); 140 141 S3DUT(mesh_get_data(mesh, &data)); 142 143 args.filename = "Cuboid facets"; 144 args.stream = fp; 145 args.type = SSTL_ASCII; 146 args.verbose = 1; 147 if((res = sstl_writer_create(&args, &writer)) != RES_OK) goto error; 148 149 FOR_EACH(i, 0, data.nprimitives) { 150 struct sstl_facet facet = SSTL_FACET_NULL; 151 const float offset[3] = {(float)PX*0.5f, (float)PY*0.5f, (float)PZ*0.5f}; 152 enum facet_flag flag = 0; 153 154 f3_set_d3(facet.vertices[0], data.positions + data.indices[i*3+0] * 3); 155 f3_set_d3(facet.vertices[1], data.positions + data.indices[i*3+1] * 3); 156 f3_set_d3(facet.vertices[2], data.positions + data.indices[i*3+2] * 3); 157 f3_add(facet.vertices[0], facet.vertices[0], offset); 158 f3_add(facet.vertices[1], facet.vertices[1], offset); 159 f3_add(facet.vertices[2], facet.vertices[2], offset); 160 flag = get_facet_flag(&facet); 161 162 if(!(flag & facet_mask)) continue; /* Discard the facet */ 163 164 if((res = sstl_write_facet(writer, &facet)) != RES_OK) goto error; 165 } 166 167 exit: 168 if(writer) SSTL(writer_ref_put(writer)); 169 return res; 170 error: 171 goto exit; 172 } 173 174 static res_T 175 setup_scene(FILE* fp) 176 { 177 ASSERT(fp); 178 fprintf(fp, "SOLID cuboid %f %f %f AUTO 300 UNKNOWN 0 BACK cuboid.stl\n", 179 LAMBDA, RHO, CP); 180 fprintf(fp, "H_BOUNDARY_FOR_SOLID Adiabatic 0 0 0 0 0 cuboid_yYzZ.stl\n"); 181 fprintf(fp, "HF_BOUNDARY_FOR_SOLID HFlux 0 0 0 %f %f %f cuboid_x.stl\n", 182 H, Text, PHI1); 183 fprintf(fp, "F_BOUNDARY_FOR_SOLID Flux %f cuboid_X.stl\n", PHI2); 184 return RES_OK; 185 } 186 187 static res_T 188 setup_scene2(FILE* fp) 189 { 190 fprintf(fp, "SOLID cuboid %f %f %f AUTO 300 UNKNOWN 0 BACK cuboid.stl\n", 191 LAMBDA, RHO, CP); 192 fprintf(fp, "FLUID Extern 1 1 %f %f FRONT cuboid.stl\n", Text, Text); 193 fprintf(fp, "SOLID_FLUID_CONNECTION Adiabatic 0 0 0 0 cuboid_yYzZ.stl\n"); 194 fprintf(fp, "F_SOLID_FLUID_CONNECTION HFlux 0 0 0 %f %f cuboid_x.stl\n", H, PHI1); 195 fprintf(fp, "F_SOLID_FLUID_CONNECTION Flux 0 0 0 0 %f cuboid_X.stl\n", PHI2); 196 return RES_OK; 197 } 198 199 static int 200 init(void) 201 { 202 struct s3dut_mesh* mesh = NULL; 203 FILE* fp_cuboid = NULL; 204 FILE* fp_cuboid_x = NULL; 205 FILE* fp_cuboid_X = NULL; 206 FILE* fp_cuboid_yYzZ = NULL; 207 FILE* fp_scene1 = NULL; 208 FILE* fp_scene2 = NULL; 209 res_T res = RES_OK; 210 int err = 0; 211 212 #define FOPEN(Name) { \ 213 if(!(fp_##Name = fopen(STR(Name)".stl", "w"))) { \ 214 perror(STR(Name)".stl"); \ 215 goto error; \ 216 } \ 217 } (void)0 218 FOPEN(cuboid); 219 FOPEN(cuboid_x); 220 FOPEN(cuboid_X); 221 FOPEN(cuboid_yYzZ); 222 #undef FOPEN 223 224 if(!(fp_scene1 = fopen(FILENAME1, "w"))) { perror(FILENAME1); goto error; } 225 if(!(fp_scene2 = fopen(FILENAME2, "w"))) { perror(FILENAME2); goto error; } 226 227 res = s3dut_create_cuboid(NULL, CUBOID_LEN_X, CUBOID_LEN_Y, CUBOID_LEN_Z, &mesh); 228 if(res != RES_OK) { 229 fprintf(stderr, "Error creating the cuboid -- %s\n", res_to_cstr(res)); 230 goto error; 231 } 232 233 if((res = write_faces(fp_cuboid, mesh, X|Y|Z)) != RES_OK) goto error; 234 if((res = write_faces(fp_cuboid_x, mesh, LX)) != RES_OK) goto error; 235 if((res = write_faces(fp_cuboid_X, mesh, UX)) != RES_OK) goto error; 236 if((res = write_faces(fp_cuboid_yYzZ, mesh, Y|Z)) != RES_OK) goto error; 237 if((res = setup_scene(fp_scene1)) != RES_OK) goto error; 238 if((res = setup_scene2(fp_scene2)) != RES_OK) goto error; 239 240 exit: 241 if(mesh) S3DUT(mesh_ref_put(mesh)); 242 if(fp_cuboid) fclose(fp_cuboid); 243 if(fp_cuboid_x) fclose(fp_cuboid_x); 244 if(fp_cuboid_X) fclose(fp_cuboid_X); 245 if(fp_scene1) fclose(fp_scene1); 246 if(fp_scene2) fclose(fp_scene2); 247 return err; 248 error: 249 err = 1; 250 goto exit; 251 } 252 253 static double 254 analytical_solution(void) 255 { 256 return PHI2 / LAMBDA * PX + (Text + (PHI1 + PHI2)/H); 257 } 258 259 static int 260 run(const char* command) 261 { 262 FILE* output = NULL; 263 double E = 0; 264 double SE = 0; 265 double ref = 0; 266 int n = 0; 267 int err = 0; 268 int status = 0; 269 270 printf("%s\n", command); 271 272 if(!(output = popen(command, "r"))) { 273 fprintf(stderr, "Error executing stardis -- %s\n", strerror(errno)); 274 err = errno; 275 goto error; 276 } 277 278 if((n = fscanf(output, "%lf %lf %*d %*d", &E, &SE)), n != 2 && n != EOF) { 279 fprintf(stderr, "Error reading the output stream -- %s\n", strerror(errno)); 280 err = errno; 281 goto error; 282 } 283 284 /* Check command exit status */ 285 if((status=pclose(output)), output=NULL, status) { 286 if(status == -1) err = errno; 287 else if(WIFEXITED(status)) err = WEXITSTATUS(status); 288 else if(WIFSIGNALED(status)) err = WTERMSIG(status); 289 else if(WIFSTOPPED(status)) err = WSTOPSIG(status); 290 else FATAL("Unreachable code.\n"); 291 goto error; 292 } 293 294 ref = analytical_solution(); 295 printf("T = %g ~ %g +/- %g\n", ref, E, SE); 296 if(!eq_eps(ref, E, SE*3)) { 297 err = 1; 298 goto error; 299 } 300 301 exit: 302 if(output) pclose(output); 303 return err; 304 error: 305 goto exit; 306 } 307 308 /******************************************************************************* 309 * The test 310 ******************************************************************************/ 311 int 312 main(int argc, char** argv) 313 { 314 int err = 0; 315 (void)argc, (void)argv; 316 317 if((err = init())) goto error; 318 if((err = run(COMMAND STR(FILENAME1)))) goto error; 319 if((err = run(COMMAND STR(FILENAME2)))) goto error; 320 321 exit: 322 return err; 323 error: 324 goto exit; 325 }