#include<stdio.h>
#include<math.h>
int main(){
int i,j,k,l,m,n;
double gamma=1.4,dltt=0.1,dltx=0.1;
double q[3][100],qnew[3][100];
double E1[100],E2[100],E3[100];
double count[100],rho[100],u[100],p[100],ene[100],H[100];
double rhoave[99],uave[99],cave[99],Have[99];
double detA[99],R[3][3],Ri[3][3];
double lambda[3],b1,b2,A[3][3]={0.0};
FILE *fp;
fp=fopen("FDS.csv","w");
//初期条件//
for(j=0;j<51;j++){
count[j]=j;
rho[j] = 1.0;
u[j] = 0.0;
p[j] = 0.8;
}
for(j=51;j<100;j++){
count[j]=j;
rho[j] = 0.125;
u[j] = 0.0;
p[j] = 0.1;
}
for(j=0;j<100;j++){
ene[j] = p[j]/(gamma - 1.0) + 0.5*rho[j]*u[j]*u[j];
H[j] = (ene[j]+p[j])/rho[j];
}
for(j=0;j<100;j++){
q[0][j] = rho[j];
q[1][j] = rho[j]*u[j];
q[2][j] = ene[j];
}
//メインループ//
for(m=0;m<10;m++){
//平均値//
for(j=0;j<99;j++){
rhoave[j] = sqrt(rho[j]*rho[j+1]);
}
for(j=0;j<99;j++){
uave[j] = (sqrt(rho[j])*u[j] + sqrt(rho[j+1])*u[j+1])/(sqrt(rho[j]) + sqrt(rho[j+1]));
}
for(j=0;j<99;j++){
Have[j] = (sqrt(rho[j])*H[j] + sqrt(rho[j+1])*H[j+1])/(sqrt(rho[j])+sqrt(rho[j+1]));
}
for(j=0;j<99;j++){
cave[j] = sqrt((gamma - 1.0)*(Have[j] - 0.5*uave[j]*uave[j]));
}
//行列//
//定義//
for(j=0;j<100;j++){
lambda[0]=fabs(uave[j] - cave[j]);
lambda[1]=fabs(uave[j]);
lambda[2]=fabs(uave[j] + cave[j]);
R[0][0]=1.0;
R[0][1]=1.0;
R[0][2]=1.0;
R[1][0]=uave[j] - cave[j];
R[1][1]=uave[j];
R[1][2]=uave[j] + cave[j];
R[2][0]=Have[j] - uave[j]*cave[j];
R[2][1]=0.5*uave[j]*uave[j];
R[2][2]=Have[j] + uave[j]*cave[j];
b1 = 0.5*(gamma - 1.0)*uave[j]*uave[j]/(cave[j]*cave[j]);
b2 = (gamma - 1.0)/(cave[j]*cave[j]);
Ri[0][0]= 0.5*(b1 + uave[j]/cave[j]);
Ri[0][1]= -0.5*(1.0/cave[j] + b2*uave[j]);
Ri[0][2]= 0.5*b2;
Ri[1][0]= 1.0 - b1;
Ri[1][1]= b2*uave[j];
Ri[1][2]= -b2;
Ri[2][0]= 0.5*(b1 - uave[j]/cave[j]);
Ri[2][1]= 0.5*(1.0/cave[j] - b2*uave[j]);
Ri[2][2]= 0.5*b2;
//行列計算//
for(l=0;l<3;l++){
for(m=0;m<3;m++){
A[l][m]=0.0;
for(n=0;n<3;n++){
A[l][m] += R[l][n]*lambda[n]*Ri[n][m];
}
}
}
E1[j]=0.5*(q[1][j+1]+q[1][j]
-A[0][0]*(q[0][j+1]-q[0][j])
-A[0][1]*(q[1][j+1]-q[1][j])
-A[0][2]*(q[2][j+1]-q[2][j]));
E2[j]=0.5*((gamma - 1.0)*(q[2][j+1]+q[2][j] - 0.5*(q[1][j+1]*q[1][j+1]/(q[0][j+1]*q[0][j+1]) + q[1][j]*q[1][j]/(q[0][j]*q[0][j]))) + (q[1][j+1]*q[1][j+1]/q[0][j+1]+q[1][j]*q[1][j]/q[0][j])
-A[1][0]*(q[0][j+1]-q[0][j])
-A[1][1]*(q[1][j+1]-q[1][j])
-A[1][2]*(q[2][j+1]-q[2][j]));
E3[j]=0.5*((q[2][j+1] + (gamma - 1.0)*(q[2][j+1] - 0.5*q[1][j+1]*q[1][j+1]/q[0][j+1]))*(q[1][j+1]/q[0][j+1]) + (q[2][j] + (gamma - 1.0)*(q[2][j] - 0.5*q[1][j]*q[1][j]/q[0][j]))*(q[1][j]/q[0][j])
-A[2][0]*(q[0][j+1]-q[0][j])
-A[2][1]*(q[1][j+1]-q[1][j])
-A[2][2]*(q[2][j+1]-q[2][j]));
}
//次ステップのQの計算//
for(j=1;j<99;j++){
qnew[0][j]=q[0][j]-(E1[j]-E1[j-1]);
qnew[1][j]=q[1][j]-(E2[j]-E2[j-1]);
qnew[2][j]=q[2][j]-(E3[j]-E3[j-1]);
}
//端点//
for(i=0;i<3;i++){
qnew[i][0]=q[i][0];
qnew[i][99]=q[i][99];
}
//qnewの中身をrho,u,pなどに変換//
for(j=0;j<100;j++){
rho[j]=qnew[0][j];
u[j] =qnew[1][j]/qnew[0][j];
p[j] =(gamma - 1.0)*(qnew[2][j] - 0.5*qnew[1][j]*qnew[1][j]/qnew[0][j]);
H[j] =(qnew[2][j] + (gamma - 1.0)*(qnew[2][j] - 0.5*qnew[1][j]*qnew[1][j]/qnew[0][j]))/qnew[0][j];
}
}
//ファイル出力//
for(j=0;j<100;j++){
fprintf(fp,"%.4f ", count[j]);
}
fprintf(fp,"\n");
for(j=0;j<100;j++){
fprintf(fp,"%.4f ", rho[j]);
}
fprintf(fp,"\n");
for(j=0;j<100;j++){
fprintf(fp,"%.4f ", u[j]);
}
fprintf(fp,"\n");
for(j=0;j<100;j++){
fprintf(fp,"%.4f ", p[j]);
}
fclose(fp);
return 0;
}
I2luY2x1ZGU8c3RkaW8uaD4KI2luY2x1ZGU8bWF0aC5oPgoKCmludCBtYWluKCl7CgoJaW50IGksaixrLGwsbSxuOwoJZG91YmxlIGdhbW1hPTEuNCxkbHR0PTAuMSxkbHR4PTAuMTsKCWRvdWJsZSBxWzNdWzEwMF0scW5ld1szXVsxMDBdOwoJZG91YmxlIEUxWzEwMF0sRTJbMTAwXSxFM1sxMDBdOwoJZG91YmxlIGNvdW50WzEwMF0scmhvWzEwMF0sdVsxMDBdLHBbMTAwXSxlbmVbMTAwXSxIWzEwMF07Cglkb3VibGUgcmhvYXZlWzk5XSx1YXZlWzk5XSxjYXZlWzk5XSxIYXZlWzk5XTsKCWRvdWJsZSBkZXRBWzk5XSxSWzNdWzNdLFJpWzNdWzNdOwoJZG91YmxlIGxhbWJkYVszXSxiMSxiMixBWzNdWzNdPXswLjB9OwoKCUZJTEUgKmZwOwoJZnA9Zm9wZW4oIkZEUy5jc3YiLCJ3Iik7CgoJLy/liJ3mnJ/mnaHku7YvLwoJZm9yKGo9MDtqPDUxO2orKyl7CgkJY291bnRbal09ajsKCQlyaG9bal0gPSAxLjA7CgkJdVtqXSAgID0gMC4wOwoJCXBbal0gICA9IDAuODsKCX0KCWZvcihqPTUxO2o8MTAwO2orKyl7CgkJY291bnRbal09ajsKCQlyaG9bal0gPSAwLjEyNTsKCQl1W2pdICAgPSAwLjA7CgkJcFtqXSAgID0gMC4xOwoJfQoJZm9yKGo9MDtqPDEwMDtqKyspewoJCWVuZVtqXSA9IHBbal0vKGdhbW1hIC0gMS4wKSArIDAuNSpyaG9bal0qdVtqXSp1W2pdOwoJCUhbal0gPSAoZW5lW2pdK3Bbal0pL3Job1tqXTsKCX0KCWZvcihqPTA7ajwxMDA7aisrKXsKCQlxWzBdW2pdID0gcmhvW2pdOwoJCXFbMV1bal0gPSByaG9bal0qdVtqXTsKCQlxWzJdW2pdID0gZW5lW2pdOwoJfQoKCgoKCS8v44Oh44Kk44Oz44Or44O844OXLy8KCWZvcihtPTA7bTwxMDttKyspewoJCS8v5bmz5Z2H5YCkLy8KCQlmb3Ioaj0wO2o8OTk7aisrKXsKCQkJcmhvYXZlW2pdID0gc3FydChyaG9bal0qcmhvW2orMV0pOwoJCX0KCQlmb3Ioaj0wO2o8OTk7aisrKXsKCQkJdWF2ZVtqXSA9IChzcXJ0KHJob1tqXSkqdVtqXSArIHNxcnQocmhvW2orMV0pKnVbaisxXSkvKHNxcnQocmhvW2pdKSArIHNxcnQocmhvW2orMV0pKTsKCQl9CgkJZm9yKGo9MDtqPDk5O2orKyl7CgkJCUhhdmVbal0gPSAoc3FydChyaG9bal0pKkhbal0gKyBzcXJ0KHJob1tqKzFdKSpIW2orMV0pLyhzcXJ0KHJob1tqXSkrc3FydChyaG9baisxXSkpOwoJCX0KCQlmb3Ioaj0wO2o8OTk7aisrKXsKCQkJY2F2ZVtqXSA9IHNxcnQoKGdhbW1hIC0gMS4wKSooSGF2ZVtqXSAtIDAuNSp1YXZlW2pdKnVhdmVbal0pKTsKCQl9CgoJCS8v6KGM5YiXLy8KCQkvL+Wumue+qS8vCgkJZm9yKGo9MDtqPDEwMDtqKyspewoKCQkJbGFtYmRhWzBdPWZhYnModWF2ZVtqXSAtIGNhdmVbal0pOwoJCQlsYW1iZGFbMV09ZmFicyh1YXZlW2pdKTsKCQkJbGFtYmRhWzJdPWZhYnModWF2ZVtqXSArIGNhdmVbal0pOwoKCQkJUlswXVswXT0xLjA7CgkJCVJbMF1bMV09MS4wOwoJCQlSWzBdWzJdPTEuMDsKCQkJUlsxXVswXT11YXZlW2pdIC0gY2F2ZVtqXTsKCQkJUlsxXVsxXT11YXZlW2pdOwoJCQlSWzFdWzJdPXVhdmVbal0gKyBjYXZlW2pdOwoJCQlSWzJdWzBdPUhhdmVbal0gLSB1YXZlW2pdKmNhdmVbal07CgkJCVJbMl1bMV09MC41KnVhdmVbal0qdWF2ZVtqXTsKCQkJUlsyXVsyXT1IYXZlW2pdICsgdWF2ZVtqXSpjYXZlW2pdOwoKCQkJYjEgPSAwLjUqKGdhbW1hIC0gMS4wKSp1YXZlW2pdKnVhdmVbal0vKGNhdmVbal0qY2F2ZVtqXSk7CgkJCWIyID0gKGdhbW1hIC0gMS4wKS8oY2F2ZVtqXSpjYXZlW2pdKTsKCgoJCQlSaVswXVswXT0gIDAuNSooYjEgKyB1YXZlW2pdL2NhdmVbal0pOwoJCQlSaVswXVsxXT0gLTAuNSooMS4wL2NhdmVbal0gKyBiMip1YXZlW2pdKTsKCQkJUmlbMF1bMl09ICAwLjUqYjI7CgkJCVJpWzFdWzBdPSAgMS4wIC0gYjE7CgkJCVJpWzFdWzFdPSAgYjIqdWF2ZVtqXTsKCQkJUmlbMV1bMl09IC1iMjsKCQkJUmlbMl1bMF09ICAwLjUqKGIxIC0gdWF2ZVtqXS9jYXZlW2pdKTsKCQkJUmlbMl1bMV09ICAwLjUqKDEuMC9jYXZlW2pdIC0gYjIqdWF2ZVtqXSk7CgkJCVJpWzJdWzJdPSAgMC41KmIyOwoJCQkvL+ihjOWIl+ioiOeuly8vCgkJCWZvcihsPTA7bDwzO2wrKyl7CgkJCQlmb3IobT0wO208MzttKyspewoJCQkJCUFbbF1bbV09MC4wOwoJCQkJCWZvcihuPTA7bjwzO24rKyl7CgkJCQkJCUFbbF1bbV0gKz0gUltsXVtuXSpsYW1iZGFbbl0qUmlbbl1bbV07CgkJCQkJfQoJCQkJfQoJCQl9CgoJCQlFMVtqXT0wLjUqKHFbMV1baisxXStxWzFdW2pdCgkJCS1BWzBdWzBdKihxWzBdW2orMV0tcVswXVtqXSkKCQkJCS1BWzBdWzFdKihxWzFdW2orMV0tcVsxXVtqXSkKCQkJCS1BWzBdWzJdKihxWzJdW2orMV0tcVsyXVtqXSkpOwoKCQkJRTJbal09MC41KigoZ2FtbWEgLSAxLjApKihxWzJdW2orMV0rcVsyXVtqXSAtIDAuNSoocVsxXVtqKzFdKnFbMV1baisxXS8ocVswXVtqKzFdKnFbMF1baisxXSkgKyBxWzFdW2pdKnFbMV1bal0vKHFbMF1bal0qcVswXVtqXSkpKSArIChxWzFdW2orMV0qcVsxXVtqKzFdL3FbMF1baisxXStxWzFdW2pdKnFbMV1bal0vcVswXVtqXSkKCQkJCS1BWzFdWzBdKihxWzBdW2orMV0tcVswXVtqXSkKCQkJCS1BWzFdWzFdKihxWzFdW2orMV0tcVsxXVtqXSkKCQkJCS1BWzFdWzJdKihxWzJdW2orMV0tcVsyXVtqXSkpOwoKCQkJRTNbal09MC41KigocVsyXVtqKzFdICsgKGdhbW1hIC0gMS4wKSoocVsyXVtqKzFdIC0gMC41KnFbMV1baisxXSpxWzFdW2orMV0vcVswXVtqKzFdKSkqKHFbMV1baisxXS9xWzBdW2orMV0pICsgKHFbMl1bal0gKyAoZ2FtbWEgLSAxLjApKihxWzJdW2pdIC0gMC41KnFbMV1bal0qcVsxXVtqXS9xWzBdW2pdKSkqKHFbMV1bal0vcVswXVtqXSkKCQkJCS1BWzJdWzBdKihxWzBdW2orMV0tcVswXVtqXSkKCQkJCS1BWzJdWzFdKihxWzFdW2orMV0tcVsxXVtqXSkKCQkJCS1BWzJdWzJdKihxWzJdW2orMV0tcVsyXVtqXSkpOwoKCQl9CgkJLy/mrKHjgrnjg4bjg4Pjg5fjga5R44Gu6KiI566XLy8KCQlmb3Ioaj0xO2o8OTk7aisrKXsKCQkJcW5ld1swXVtqXT1xWzBdW2pdLShFMVtqXS1FMVtqLTFdKTsKCQkJcW5ld1sxXVtqXT1xWzFdW2pdLShFMltqXS1FMltqLTFdKTsKCQkJcW5ld1syXVtqXT1xWzJdW2pdLShFM1tqXS1FM1tqLTFdKTsKCQl9CgoJCS8v56uv54K5Ly8KCQlmb3IoaT0wO2k8MztpKyspewoJCQlxbmV3W2ldWzBdPXFbaV1bMF07CgkJCXFuZXdbaV1bOTldPXFbaV1bOTldOwoJCX0KCgkJLy9xbmV344Gu5Lit6Lqr44KScmhvLHUscOOBquOBqeOBq+WkieaPmy8vCgkJZm9yKGo9MDtqPDEwMDtqKyspewoJCQlyaG9bal09cW5ld1swXVtqXTsKCQkJdVtqXSAgPXFuZXdbMV1bal0vcW5ld1swXVtqXTsKCQkJcFtqXSAgPShnYW1tYSAtIDEuMCkqKHFuZXdbMl1bal0gLSAwLjUqcW5ld1sxXVtqXSpxbmV3WzFdW2pdL3FuZXdbMF1bal0pOwoJCQlIW2pdICA9KHFuZXdbMl1bal0gKyAoZ2FtbWEgLSAxLjApKihxbmV3WzJdW2pdIC0gMC41KnFuZXdbMV1bal0qcW5ld1sxXVtqXS9xbmV3WzBdW2pdKSkvcW5ld1swXVtqXTsKCQl9CgoJfQoKCS8v44OV44Kh44Kk44Or5Ye65YqbLy8KCWZvcihqPTA7ajwxMDA7aisrKXsKCQlmcHJpbnRmKGZwLCIlLjRmICIsIGNvdW50W2pdKTsKCX0KCWZwcmludGYoZnAsIlxuIik7Cglmb3Ioaj0wO2o8MTAwO2orKyl7CgkJZnByaW50ZihmcCwiJS40ZiAiLCByaG9bal0pOwoJfQoJZnByaW50ZihmcCwiXG4iKTsKCWZvcihqPTA7ajwxMDA7aisrKXsKCQlmcHJpbnRmKGZwLCIlLjRmICIsIHVbal0pOwoJfQoJZnByaW50ZihmcCwiXG4iKTsKCWZvcihqPTA7ajwxMDA7aisrKXsKCQlmcHJpbnRmKGZwLCIlLjRmICIsIHBbal0pOwoJfQoJZmNsb3NlKGZwKTsKCglyZXR1cm4gMDsKfQ==