Q&A
コンパイラの警告レベルを上げたら何かしらの警告は出ませんか？
前提・実現したいこと
C言語で数値計算を行なっているのですが，実行のタイミングによってnanになることがあります．
わかる方どうかお助けください．
発生している問題
出力値が実行のタイミングによってnanになったり数値になったいるする（数値は基本的に変わらない）
ソースコード
c
1#include <stdio.h> 2#include <stdlib.h> 3#define _USE_MATH_DEFINES 4#include <math.h> 5#include <sys/stat.h> 6#include<unistd.h> 7#include<sys/types.h> 8 9#define NX 10 10#define NY 10 11#define NZ 10 //odd number 12#define Tmax 1.0 13 14double fun_p(double a) 15{ 16 return pow(10, (5.006 + 6.806 - 30295 / a)); 17} 18 19double jac(double b, double c) 20{ 21 //return 0.83*pow((pow(0.83,2)-(pow((b-1.0e-3),2)+pow((c-1.0e-3),2))),(-0.5)); 22 return 1.0; 23} 24 25double ver(d, e){ 26 return sqrt(pow(0.83, 2) - (pow(d, 2) + pow(e, 2))) / 0.83; 27 return 1.0; 28} 29 30int lmin[NY + 2][NZ + 2], lmax[NY + 2][NZ + 2]; 31double temp[NX + 2][NY + 2][NZ + 2][2]; 32double pl[NX + 2][NY + 2], ph[NX + 2][NY + 2]; 33 34int main(void) 35{ 36 const char* file = "/Users/akitatomoya/Public/NewFolder"; 37 38 if (mkdir(file,S_IRWXU) == 0) 39 printf("Made a directory\n"); 40 //if (mkdir(folder) == 0) 41 //printf("Made a directry\n", folder); 42 43 //const char* file = "C:\NewFolder"; 44 45 double dt, dx, time, pp, lp, cp, lam, rho, 46 tinf, ss, phi, ep, sig, eta, dd, ff, p1, p2, 47 rs, x, y, z, x1, y1, z1; 48 double xmax, ymax, zmax, tem, dy, dz, eps, resd, tp, 49 ckx, cky, ckz, alpha, rs_max, rs_min, 50 pres,sss, ps, cfl; 51 int i, j, k, imax, jmax, kmax, n, size, count, iter; 52 int *int_p; 53 double *double_p; 54 char filename[20]; 55 FILE* fp; 56 FILE* fl; 57 //fopen_s(&fp, "C:\thrust.csv", "w"); 58 dd = 1.66e-3; 59 xmax = 2.0e-3; 60 ymax = 2.0e-3; 61 zmax = 2.0e-3; 62 rs = dd*0.5; 63 //dt = Tmax / (double)NT; 64 dx = xmax / (double)NX; 65 dy = ymax / (double)NY; 66 dz = zmax / (double)NZ; 67 imax = NX; 68 jmax = NY; 69 kmax = NZ; 70 time = 0; 71 tm = 1850.0 + 273.0; 72 tinf = 1850.0 + 273.0; 73 tem = 30.0 + 273.0; 74 phi = 0.7e-3; 75 ss = 0.25 * M_PI * pow(phi, 2); 76 ep = 0.32; 77 sig = 5.67e-8; 78 eta = 0.121; 79 x1 = xmax/2; 80 y1 = ymax / 2; 81 z1 = zmax / 2; 82 rho = 5800.0; 83 cp = 436.0; 84 lam = 94.61; 85 alpha = lam / (rho * cp); 86 cfl=0.1; 87 dt=cfl*pow(dx,2)/2/alpha; 88 ckx = alpha * dt / (dx * dx); 89 cky = alpha * dt / (dy * dy); 90 ckz = alpha * dt / (dz * dz); 91 printf("nt=%5.3e\n", Tmax/dt); 92 fl=fopen("thrust.csv","w"); 93 94 //initial condition 95 x = 0.0; 96 for (i = 0; i <= imax; i++){ 97 y = 0.0; 98 for (j = 0; j <= jmax; j++) { 99 z = 0.0; 100 count = 0; 101 for (k = 0; k <= kmax; k++) { 102 if (z >= z1 - sqrt(pow(rs, 2) - pow((y - y1), 2) - pow((x - x1), 2)) && z <= z1 + sqrt(pow(rs, 2) - pow((y - y1), 2) - pow((x - x1), 2))) { 103 lmax[i][j] = k; 104 count = count + 1; //calculate the element count 105 } 106 107 if (z >= z1 - sqrt(pow(rs, 2) - pow((y - y1), 2) - pow((x - x1), 2)) && z <= z1 + sqrt(pow(rs, 2) - pow((y - y1), 2) - pow((x - x1), 2))) temp[i][j][k][0] = tinf; 108 else temp[i][j][k][1] = tem; 109 z=z+dz; 110 } 111 lmin[i][j] = lmax[i][j] - count+1; 112 y=y+dy; 113 } 114 x=x+dx; 115 } 116 117 //time integration 118 for (n = 0; time <= Tmax; n++) 119 { 120 chdir("/Users/akitatomoya/Public/NewFolder"); 121 if(n%10000==0){ 122 sprintf(filename,"tec%03d.txt",n); 123 fp=fopen(filename, "w"); 124 fprintf(fp, "x,y,z,Temperature(K)\n"); 125 126 for (i = 0; i <= imax; i++) { 127 for (j = 0; j <= jmax; j++) { 128 for (k = 0; k <= kmax; k++) { 129 fprintf(fp, "%8.3e,", i * dx); 130 fprintf(fp, "%8.3e,", j * dy); 131 fprintf(fp, "%8.3e,", k * dz); 132 fprintf(fp, "%8.3e\n", temp[i][j][k][1]); 133 } 134 } 135 } 136 } 137 138 /**************************************************** 139 laser power 140 ****************************************************/ 141 if(time>0.1)lp=100; 142 if(time>0.9)lp=0; 143 ps = lp / ss; 144 145 146 /**************************************************** 147 explicit method 148 ****************************************************/ 149 x = dx; 150 for (i = 1; i <= imax; i++) { 151 y = dy; 152 for (j = 1; j <= jmax; j++) { 153 z = dz; 154 for (k = 1; k <= kmax; k++) { 155 //rho = 6782 - 0.27 *temp[i][j][k][0]; 156 cp = 436 - 0.0813 * (temp[i][j][k][0] - tm); 157 lam = 94.61 + 4.41e-2 * (temp[i][j][k][0] - tm); 158 alpha = lam / (rho * cp); 159 ckx = alpha * dt / (dx * dx); 160 cky = alpha * dt / (dy * dy); 161 ckz = alpha * dt / (dz * dz); 162 if (z >= z1 - sqrt(pow(rs, 2) - pow((y - y1), 2) - pow((x - x1), 2)) && z <= z1 + sqrt(pow(rs, 2) - pow((y - y1), 2) - pow((x - x1), 2))) { 163 if (i == lmin[j][k] || i == lmax[j][k] || j == lmin[i][k] || j == lmax[i][k] || k == lmin[i][j] || k == lmax[i][j]) { 164 if (i == lmin[j][k]) temp[i - 1][j][k][0] = temp[i][j][k][0] - dx*ep*sig*pow(temp[i][j][k][0], 4)/lam; 165 if (i == lmax[j][k]) temp[i + 1][j][k][0] = temp[i][j][k][0] - dx*ep*sig*pow(temp[i][j][k][0], 4)/lam; 166 if (j == lmin[i][k]) temp[i][j - 1][k][0] = temp[i][j][k][0] - dy*ep*sig*pow(temp[i][j][k][0], 4)/lam; 167 if (j == lmax[i][k]) temp[i][j + 1][k][0] = temp[i][j][k][0] - dy*ep*sig*pow(temp[i][j][k][0], 4)/lam; 168 if (k == lmin[i][j]) temp[i][j][k - 1][0] = temp[i][j][k][0] - dz*ep*sig*pow(temp[i][j][k][0], 4)/lam; 169 if (k == lmax[i][j]) temp[i][j][k + 1][0] = temp[i][j][k][0] - dz*ep*sig*pow(temp[i][j][k][0], 4)/lam + eta*ps*dz/lam; 170 } 171 temp[i][j][k][1] = temp[i][j][k][0]+ckx*(temp[i+1][j][k][0]-2*temp[i][j][k][0]+temp[i-1][j][k][0])+ 172 cky*(temp[i][j+1][k][0]-2*temp[i][j][k][0]+temp[i][j-1][k][0])+ckz*(temp[i][j][k+1][0]-2*temp[i][j][k][0]+temp[i][j][k-1][0]); 173 } 174 z = z + dz; 175 } 176 y = y + dy; 177 } 178 x = x + dz; 179 } 180 181 for (i = 0; i <= imax; i++) { 182 for (j = 0; j <= jmax - 1; j++) { 183 for (k = 0; k <= kmax; k++) { 184 temp[i][j][k][0] = temp[i][j][k][1]; 185 } 186 } 187 } 188 189 /************************************************ 190 caluculate thrust 191 ************************************************/ 192 //_chdir("C:\Users\Public\ThrustFolder"); 193 //fopen_s(&fl, "thrust.csv", "w"); 194 195 if (n == 0) fprintf(fl, "Time (s),Thrust (N),LaserPower (W)\n"); 196 x = 0; 197 for (i = 0; i <= imax; i++) { 198 y = 0; 199 for (j = 0; j <= jmax; j++) { 200 z = 0; 201 pl[i][j] = 0; 202 ph[i][j] = 0; 203 for (k = 0; k <= kmax; k++) { 204 if (z >= z1 - sqrt(pow(rs, 2) - pow((y - y1), 2) - pow((x - x1), 2)) && z <= z1 + sqrt(pow(rs, 2) - pow((y - y1), 2) - pow((x - x1), 2))) { 205 if (k == lmin[i][j] || k == lmax[i][j]) { 206 if (z > z1) ph[i][j] = fun_p(temp[i][j][k][1]); 207 else if(z<z1) pl[i][j] = fun_p(temp[i][j][k][1]); 208 } 209 } 210 z = z + dz; 211 } 212 y = y + dy; 213 } 214 x = x + dx; 215 } 216 if (n % 100==0){ 217 x=0; 218 ff = 0; 219 for (i = 0; i <= imax; i++) { 220 y=0; 221 for (j = 0; j <= jmax; j++) { 222 pres = ph[i][j] - pl[i][j]; 223 ff = ff + ver(x,y) * pres * dx * dy; 224 y = y + dy ; 225 } 226 x = x + dx; 227 } 228 fprintf(fl, "%8.3lf,", time); 229 fprintf(fl, "%5.3e,", ff); 230 fprintf(fl, "%8.3lf\n", lp); 231 } 232 if(n % 1000 == 0) printf("Time: %8.3lf, Temp: %8.3lf\n", time, temp[imax/2][jmax/2][kmax/2][1]); 233 time = time + dt; 234 fclose(fp); 235 } 236 fclose(fl); 237 238 return 0; 239}
補足情報（FW/ツールのバージョンなど）
gccのバージョン: 13.0.0
M1 MacBook Air
macOS Monterey
何が nan になるのでしょうか。
私の環境(Windows/cygwin)では (S_IRWXU はともかく) "tm" が無いと言われますが、これは?
main の中の処理をもっと関数に分けて、入力・出力をはっきりさせてください。
また、(コンパイラに因りますが)変数の宣言は各々の使用開始時に行うとか、定数と変数の見分けが付くように定数名は大文字・変数名は小文字にする等の命名規則をご利用くださると、他人にも見易くなると思います。
警告レベルを上げるとtmがないと言われました．
最初にtmを宣言すると3回連続で数値が出力されました．（おそらくこれで解決しました）
ご教示いただいた皆様ありがとうございます．
未定義変数がエラーでなく警告だというのは謎ですね。
お使いの環境ではどこかの include ファイルの中に定義があって cygwin には無いモノなのかと思ったのですが、どうなんでしょうね。
あなたの回答
tips
プレビュー
