sadist_conducto_radiative.c (6457B)
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/rsys.h> 21 #include <rsys/math.h> 22 23 #include <errno.h> 24 #include <stdio.h> 25 #include <string.h> /* strerror */ 26 #include <wait.h> /* WIFEXITED, WEXITSTATUS */ 27 28 /* 29 * The scene is a solid cube whose +/-X faces exchange with the radiative 30 * environment. All other faces are adiabatic. The radiative environment is set 31 * at T0 along the -X direction and T1 along the +X direction. 32 * 33 * This test calculates the steady temperature at a given position of the probe 34 * in the solid and validates it against its analytical solution. 35 * 36 * //////(1,1,1) 37 * +-------+ 38 * ---> /' /| <--- 39 * T0 K ---> +-------+ | E=1 <--- T1 K 40 * Y ---> E=1 | +.....|.+ <--- 41 * | |, |/ 42 * o--X +-------+ 43 * / (0,0,0) //// 44 * Z 45 */ 46 47 #define FILENAME_CUBE "cube.stl" 48 #define FILENAME_ADIABATIC "adiabatic_boundary.stl" 49 #define FILENAME_RADIATIVE "radiative_boundary.stl" 50 #define FILENAME_SCENE "scene.txt" 51 52 #define EMISSIVITY 1.0 53 #define TREF 315.0 /* [K] */ 54 #define T0 280.0 /* [K] */ 55 #define T1 350.0 /* [K] */ 56 #define LAMBDA 0.1 57 58 #define BOLTZMANN_CONSTANT 5.6696e-8 /* [W/m^2/K^4] */ 59 #define HR (4.0*BOLTZMANN_CONSTANT * TREF*TREF*TREF * EMISSIVITY) 60 61 /* Probe position */ 62 #define PX 0.2 63 #define PY 0.4 64 #define PZ 0.6 65 66 /* The commands to be executed */ 67 #define COMMAND1 "stardis -V3 -p "STR(PX)","STR(PY)","STR(PZ)" -M "FILENAME_SCENE 68 #define COMMAND2 COMMAND1 " -a wos" 69 70 static const double vertices[] = { 71 0.0, 0.0, 0.0, 72 1.0, 0.0, 0.0, 73 0.0, 1.0, 0.0, 74 1.0, 1.0, 0.0, 75 0.0, 0.0, 1.0, 76 1.0, 0.0, 1.0, 77 0.0, 1.0, 1.0, 78 1.0, 1.0, 1.0 79 }; 80 static const size_t nvertices = sizeof(vertices) / (sizeof(double)*3); 81 82 static const size_t indices[] = { 83 0, 4, 2, 2, 4, 6, /* -x */ 84 3, 7, 5, 5, 1, 3, /* +x */ 85 0, 1, 5, 5, 4, 0, /* -y */ 86 2, 6, 7, 7, 3, 2, /* +y */ 87 0, 2, 1, 1, 2, 3, /* -z */ 88 4, 5, 6, 6, 5, 7 /* +z */ 89 }; 90 static const size_t ntriangles = sizeof(indices) / (sizeof(size_t)*3); 91 92 /******************************************************************************* 93 * Helper functions 94 ******************************************************************************/ 95 static void 96 setup_scene(FILE* fp) 97 { 98 ASSERT(fp); 99 fprintf(fp, "PROGRAM radenv_1d libsadist_radenv_1d.so\n"); 100 fprintf(fp, "PROGRAM solid_fluid libsadist_solid_fluid.so\n"); 101 fprintf(fp, "FLUID environment 1 1 0 UNKNOWN BACK %s\n", FILENAME_CUBE); 102 fprintf(fp, "SOLID cube %g 1 1 0.05 0 UNKNOWN 0 FRONT %s\n", 103 LAMBDA, FILENAME_CUBE); 104 fprintf(fp, "SOLID_FLUID_CONNECTION_PROG radiative solid_fluid %s " 105 "PROG_PARAMS -e %g -r 300\n", FILENAME_RADIATIVE, EMISSIVITY); 106 fprintf(fp, "SOLID_FLUID_CONNECTION adiabatic 0 0 0 0 %s\n", 107 FILENAME_ADIABATIC); 108 fprintf(fp, "TRAD_PROG radenv_1d PROG_PARAMS -r %g,%g -t %g,%g\n", 109 TREF, TREF, T0, T1); 110 } 111 112 static int 113 init(void) 114 { 115 FILE* fp_cube = NULL; 116 FILE* fp_adiabatic = NULL; 117 FILE* fp_radiative = NULL; 118 FILE* fp_scene = NULL; 119 int err = 0; 120 121 #define FOPEN(Fp, Filename) \ 122 if(((Fp) = fopen((Filename), "w")) == NULL) { \ 123 fprintf(stderr, "Error opening `%s' -- file %s\n", \ 124 (Filename), strerror(errno)); \ 125 err = errno; \ 126 goto error; \ 127 } (void)0 128 FOPEN(fp_cube, FILENAME_CUBE); 129 FOPEN(fp_adiabatic, FILENAME_ADIABATIC); 130 FOPEN(fp_radiative, FILENAME_RADIATIVE); 131 FOPEN(fp_scene, FILENAME_SCENE); 132 #undef FOPEN 133 134 sadist_write_stl(fp_cube, vertices, nvertices, indices, ntriangles); 135 sadist_write_stl(fp_adiabatic, vertices, nvertices, indices+12, ntriangles-4); 136 sadist_write_stl(fp_radiative, vertices, nvertices, indices, 4); 137 setup_scene(fp_scene); 138 139 exit: 140 #define FCLOSE(Fp) \ 141 if((Fp) && fclose(Fp)) { perror("fclose"); if(!err) err = errno; } (void)0 142 FCLOSE(fp_cube); 143 FCLOSE(fp_adiabatic); 144 FCLOSE(fp_radiative); 145 FCLOSE(fp_scene); 146 #undef FCLOSE 147 return err; 148 error: 149 goto exit; 150 } 151 152 static double 153 analytical_solution(void) 154 { 155 const double tmp = LAMBDA / (2*LAMBDA + HR) * (T1 - T0); 156 const double Ts0 = T0 + tmp; 157 const double Ts1 = T1 - tmp; 158 const double ref = PX*Ts1 + (1 - PX)*Ts0; 159 return ref; 160 } 161 162 static int 163 run(const char* command) 164 { 165 FILE* output = NULL; 166 double E = 0; 167 double SE = 0; 168 double ref = 0; 169 int n = 0; 170 int err = 0; 171 int status = 0; 172 173 printf("%s\n", command); 174 175 if(!(output = popen(command, "r"))) { 176 fprintf(stderr, "Error executing stardis -- %s\n", strerror(errno)); 177 err = errno; 178 goto error; 179 } 180 181 if((n = fscanf(output, "%lf %lf %*d %*d", &E, &SE)), n != 2 && n != EOF) { 182 fprintf(stderr, "Error reading the output stream -- %s\n", strerror(errno)); 183 err = errno; 184 goto error; 185 } 186 187 /* Check command exit status */ 188 if((status=pclose(output)), output=NULL, status) { 189 if(status == -1) err = errno; 190 else if(WIFEXITED(status)) err = WEXITSTATUS(status); 191 else if(WIFSIGNALED(status)) err = WTERMSIG(status); 192 else if(WIFSTOPPED(status)) err = WSTOPSIG(status); 193 else FATAL("Unreachable code.\n"); 194 goto error; 195 } 196 197 ref = analytical_solution(); 198 printf("T = %g ~ %g +/- %g\n", ref, E, SE); 199 if(!eq_eps(ref, E, SE*3)) { 200 err = 1; 201 goto error; 202 } 203 204 exit: 205 if(output) pclose(output); 206 return err; 207 error: 208 goto exit; 209 } 210 211 /******************************************************************************* 212 * The test 213 ******************************************************************************/ 214 int 215 main(int argc, char** argv) 216 { 217 int err = 0; 218 (void)argc, (void)argv; 219 220 if((err = init())) goto error; 221 if((err = run(COMMAND1))) goto error; 222 if((err = run(COMMAND2))) goto error; 223 224 exit: 225 return err; 226 error: 227 goto exit; 228 }