?? test.cpp
字號:
// test.cpp : 定義控制臺應用程序的入口點。
//
#include <errno.h> //讀文件用的
#include "stdafx.h"
#include<stdio.h>
#include<math.h>
#define N1 150 //斷面數
#define N2 20 //斷面的分層數
#define N3 50 //計算次數
#define G 9.8f
int IZ;
float A1[N1],B1[N1],C1[N1],D1[N1],E1[N1],A2[N1];
float B2[N1],DKDZJ[N1];
float C2[N1],D2[N1],E2[N1],H1[N1],H2[N1],H3[N1];
float DKDZS[N1][N2],DBDZS[N1][N2],DBDZJ[N1];
float TT;//后加上去的
int JS,JE,IK[2];//斷數的起始個數
//*********************************共享去了*****************************************
float Z0[N1],Q0[N1]; //初始水深和流量
float DZ[N1],DQ[N1]; //水深和流量的增量
float BCUP[N2][2],BCDN[N2][2]; //上下游邊界條件
float FJ[N1],GJ[N1];
//**********************************************************************************
//計算水面積
//*****************************************************************
// 各段面分層水深 AS對應過水面積 Z0初始水深 過水面積(要求的量) 斷面起始 結束 每個斷面的層數
void SufaceArea( float A[][N2], float B[][N2], float C[N1], float D[N1], int JS, int JE, int IZJ[N1])
{
//int N1=150,N2=20,N3=50;
// float A[N1][N2],
// float B[N1][N2],C[N1],D[N1],IZJ[N1];
int LevelNum,j;
float InitWaterHigh;
for(j=JS;j<=JE;j++) //求每個斷面的過水面積，斷面循環
{
LevelNum=IZJ[j]; //每個斷面的層數
InitWaterHigh=C[j]; //斷面初始水位高度
for(int i=1;i<LevelNum;i++) //求每層的過水面積
{
if(InitWaterHigh<=A[j][i]) //當前水位低于該層水位時，計算該層過水面積
break;
}
D[j]=B[j][i-1]+(B[j][i]-B[j][i-1])*(InitWaterHigh-A[j][i-1])/(A[j][i]-A[j][i-1]);//求過水面積
//計算理論參見PPT
}
}
//計算水位流量
//*****************************************************************
// 上游邊界條件 要求的量 實際分層數
void WaterHigh( float A[][2], float A1, float* A2, int IA)
{
//float A[N2][2];
for(int i=1;i<IA;i++) //每一層的流量
{
if(A1<A[i][0]) //如果A1<上游邊界流量條件
break;
}
*A2=A[i-1][1]+(A[i][1]-A[i-1][1])*(A1-A[i-1][0])/(A[i][0]-A[i-1][0]);//計算流量
//流量=上游入水流量+（本段流量-上段流量）*（當前-上）/（本斷流量-）
}
//輸入邊界
//*****************************************************************
void BOUND(float DATA[][2], float A1, int IA,int J,int IGG,int KD)
{
float ZT=0.0f;
float QT,QZD,HIQ,DFQ,QIQ;
// J=J-1;
if(IGG==1)
{
WaterHigh(DATA,A1,&ZT,IA);
DZ[J]=ZT-Z0[J];
if(KD==1)
{
FJ[J]=1000000.0;
GJ[J]=-FJ[J]*DZ[J];
}
else
DQ[J]=FJ[J]*DZ[J]+GJ[J];
}
if(IGG==2)
{
WaterHigh(DATA,A1,&QT,IA);
DQ[J]=QT-Q0[J];
if(KD==1)
{
FJ[J]=0.0;
GJ[J]=DQ[J];
}
else
DZ[J]=(DQ[J]-GJ[J])/FJ[J];
}
if(IGG==3)
{
WaterHigh(DATA,Z0[J],&QZD,IA);
HIQ=Z0[J]+0.1f;
WaterHigh(DATA,HIQ,&QIQ,IA);
DFQ=(QIQ-QZD)/0.1f;
if(KD==1)
{
FJ[J]=DFQ;
GJ[J]=QZD-Q0[J];
}
else
{
DZ[J]=(QZD-GJ[J]-Q0[J])/(FJ[J]-DFQ);
DQ[J]=FJ[J]*DZ[J]+GJ[J];
}
}
}
int _tmain(int argc, _TCHAR* argv[])
{
float ZN[N3][N1],QN[N3][N1]; //計算的水深和流量
float T,temp,DZS,DBS; //
// AS面積, BS面寬， PS濕周 各斷面分層水深 流量模數
float AS[N1][N2],BS[N1][N2], PS[N1][N2], ZS[N1][N2], RKS[N1][N2];
//對應ZJ[]的過水面面積，水面寬度，濕周和流量模數，用于計算流量模數，斷面坐標值
float AJ[N1],BJ[N1],RKJ[N1],RNJ[N1],XL[N1];
int IZJ[N1]; //各斷面對應ZS[][]的分層數
float ZETA,DT;
int NX;
FILE * fp1;
if(NULL==(fp1=fopen("E:\\WangStu\\Water Simulation\\N當前成果代碼\\test\\test\\input.txt","rb+")))
{
printf("Open file 1 error:\n");
return 1;
}
printf("Open file(s) successfull!\n");
printf("Get data from file 1 : \n");
fscanf(fp1,"%f %d %f %f", &TT, &NX, &DT, &ZETA);//時間間隔，段數，時間間隔，代表權值,應該>0且<1 一般取0.5
printf("\n時間長度(HOURS)\t河道段數\t時間間隔(HOURS)\t代表權值\n%f\t%f\t%f\t%f\t\n",TT,NX,DT,ZETA);
float NT,RS; //后加上去的,未聲明的標識符
int IBU,IBD; //后加上去的,未聲明的標識符
NT=int(TT/DT)+1.001f;
JS=0;
JE=NX;
fscanf(fp1,"%f %f %f %f %f %f %f %f %f %f %f", &XL[0],&XL[1],&XL[2],&XL[3],&XL[4],&XL[5],&XL[6],&XL[7],&XL[8],&XL[9],&XL[10]); //READ(1,*) (XL(J),J=JS,JE);//各斷面坐標值
fscanf(fp1,"%f %f %f %f %f %f %f %f %f %f %f", &RNJ[0],&RNJ[1],&RNJ[2],&RNJ[3],&RNJ[4],&RNJ[5],&RNJ[6],&RNJ[7],&RNJ[8],&RNJ[9],&RNJ[10]); //READ(1,*) (RNJ(J),J=JS,JE);
fscanf(fp1,"%d %d %d %d %d %d %d %d %d %d %d", &IZJ[0],&IZJ[1],&IZJ[2],&IZJ[3],&IZJ[4],&IZJ[5],&IZJ[6],&IZJ[7],&IZJ[8],&IZJ[9],&IZJ[10]); //READ(1,*) (IZJ(J),J=JS,JE);各斷面對應zs分層數
int i;
printf("\n各斷面分層水深值，水面寬度 一層 二層\n");
for(i=JS;i<JE;i++)//段循環
{
IZ=IZJ[i];//各段的分層數
fscanf(fp1,"%f %f %f %f",&ZS[i][0],&BS[i][0],&ZS[i][1],&BS[i][1]);//各斷面分層水深值，水面寬度
printf("%d:\t%f\t%f\t%f\t%f\t\n",i+1,ZS[i][0],BS[i][0],ZS[i][1],BS[i][1]);
}
fscanf(fp1,"%d %d", &IK[0], &IK[1]);//READ(1,*) (IK(I),I=1,2);//邊界條件指示數，IK[1]上游邊界條件IK[2]下游邊界條件，IK[*]=1給定水位時間過程Z(t) =2給定流量時間過程Q(t) =3給定流量水位關系Q(Z)
fscanf(fp1,"%d %d", &IBU, &IBD); //READ(1,*) IBU,IBD;上下游分層數
printf("\n邊界條件類型值\n[%d]\t\t[%d]\t\t\n上下游條件的個數\n[%d]\t\t[%d]\t\t\n",IK[0],IK[1],IBU,IBD);
printf("\n讀取 上 游各層的邊界條件（水位時間過程） 一層\t二層\n");
for(i=0;i<IBU;i++)//讀取 上 游各層的邊界條件//以列表形式給出時，數據的個數
{
fscanf(fp1,"%f %f", &BCUP[i][0],&BCUP[i][1] );//READ(1,*) ((BCUP(J,I),I=1,2),J=1,IBU); //上游邊界條件
printf("[%f]\t[%f]\t\n",BCUP[i][0],BCUP[i][1]);
}
printf("\n讀取 下 游各層的邊界條件（水位時間過程） 一層\t二層\n");
for(i=0;i<IBD;i++)//讀取 下 游各層的邊界條件
{
fscanf(fp1,"%f %f", &BCDN[i][0],&BCDN[i][1] );//下游邊界條件
printf("[%f]\t[%f]\t\n",BCDN[i][0],BCDN[i][1]);
}
fscanf(fp1,"%f %f %f %f %f %f %f %f %f %f %f",&Q0[0], &Q0[1], &Q0[2], &Q0[3], &Q0[4], &Q0[5], &Q0[6], &Q0[7], &Q0[8], &Q0[9], &Q0[10]);//READ(1,*) (Q0(J),J=JS,JE); //各斷面初識流量
fscanf(fp1,"%f %f %f %f %f %f %f %f %f %f %f",&Z0[0], &Z0[1], &Z0[2], &Z0[3], &Z0[4], &Z0[5], &Z0[6], &Z0[7], &Z0[8], &Z0[9], &Z0[10]);//READ(1,*) (Z0(J),J=JS,JE); //各斷面初識水深
printf("\n各斷面初始流量\n%f\t%f\t%f\t%f\t%f\t%f\t\n",Q0[0], Q0[1], Q0[2], Q0[3], Q0[4], Q0[5]);
printf("\n各斷面初始水深\n%f\t%f\t%f\t%f\t%f\t%f\t\n",Z0[0], Z0[1], Z0[2], Z0[3], Z0[4], Z0[5]);
printf("--------------------文件讀取結束-----------------\n");
fclose(fp1);
int j;
//未知變量賦初值
for(j=JS;j<JE;j++)//DO 30 J=JS,JE
{
AS[j][0]=0.0; //過水面積
PS[j][0]=0.0; //濕周
RKS[j][0]=0.0; //流量模數
DQ[j]=0.0; //水深增量
DZ[j]=0.0; //流量增量
}
//計算過水面面積，水面寬，流量模數 dbdz dkdz 與z的關系
for(j=JS;j<JE;j++)
{
IZ=IZJ[j];
for(int i=1;i<IZ;i++)//DO 20 i=2,IZ
{
DZS=ZS[j][i]-ZS[j][i-1]; //各斷分層水深
DBS=(BS[j][i]-BS[j][i-1])/2.0f; //水面寬
AS[j][i]=AS[j][i-1]+(BS[j][i-1]+BS[j][i])*DZS/2.0f; //計算面積=下面的面積+（上底+下底）*高/2
PS[j][i]=PS[j][i-1]+2.0f*sqrtf(DBS*DBS+DZS*DZS); //計算濕周=可以理解為水與土壤的接觸長度
RS=AS[j][i]/PS[j][i]; //面積/濕周
RKS[j][i]=AS[j][i]*powf(float(RS),float(2.0/3.0))/RNJ[j];//流量模數
DKDZS[j][i]=(RKS[j][i]-RKS[j][i-1])/DZS;
DBDZS[j][i]=(BS[j][i]-BS[j][i-1])/DZS;
}
DKDZS[j][0]=DKDZS[j][1];
DBDZS[j][0]=DBDZS[j][1];
}
//計算系數 A,B,C,D,E
float DX,CR1,BJM,AJM,DQJ,DZJ; //后加上去的,未聲明的標識符
int JEM1,JSP1;
JEM1=JE-1;//斷面數-1
JSP1=JS+1;//斷面數起始+1
NT=(TT/DT)+1;//TT為讀進來的值．DT為讀進來的值 時間間隔，TT可以假設是時間總長度
for(int N=0;N<NT;N++)//模擬次數循環
{
T=float(N+1)*DT;
//Z0初始水深，AJ=對應ZJ過水斷面積，JS=0，JE=最大斷面數，IZJ=各斷面對應ZS的分層數
SufaceArea(ZS,AS,Z0,AJ,JS,JE,IZJ); //ZS=各段面分層水深值，AS=對應ZS過水斷面積,求對應ZJ過水斷面積AJ
SufaceArea(ZS,BS,Z0,BJ,JS,JE,IZJ); //BS=對應ZS過水面寬度，BJ=對應ZJ水面寬度,求對應ZJ過水面寬度BJ
SufaceArea(ZS,RKS,Z0,RKJ,JS,JE,IZJ); //RKS=對應ZS流量模數，RKJ=對應ZJ流量模數,求對應ZJ流量模數RKJ
SufaceArea(ZS,DKDZS,Z0,DKDZJ,JS,JE,IZJ);//DKDZS=各斷對應不同水深的dk/dz值，DKDZJ=（dk/dz）f值,求對應ZJ（dk/dz）f值DKDZJ
SufaceArea(ZS,DBDZS,Z0,DBDZJ,JS,JE,IZJ);//DBDZS=各斷對應不同水深的db/dz值，DKDZJ=（db/dz）f值,求對應ZJ（db/dz）f值DBDZJ
for(int o=0;o<5;o++)
// printf("##### AJ=%.2f,BJ=%.2f,RKJ=%.2f,DKDZJ=%.2f,DBDZJ=%.2f\n",AJ[o],BJ[o],RKJ[o],DKDZJ[o],DBDZJ[o]);
for(j=JS;j<JEM1;j++)//斷數-1個循環
{
DX = XL[j+1]-XL[j];//各斷面坐標值，這里理解為段河道的長度
CR1 = ZETA*DT/DX;//ZETA讀進來的
BJM = BJ[j]+BJ[j+1];//水面寬度
AJM = AJ[j]+AJ[j+1];//水斷面面積
DQJ = Q0[j+1]-Q0[j];
DZJ = Z0[j+1]-Z0[j];
//計算各個參數或系數
A1[j] = -4.0f*CR1/BJM;
B1[j] = 1.0f-4.0f*CR1*ZETA*DQJ*DBDZJ[j]/BJM/BJM;
C1[j] = 4.0f/BJM*CR1;
D1[j] = 1.0f-4.0f*CR1*DQJ*DBDZJ[j+1]/BJM/BJM;
E1[j] = -4*CR1*DQJ/BJM;
A2[j] = 1.0f-4.0f*CR1*Q0[j]/AJ[j]+2.0f*CR1*DX*G*AJ[j]*fabsf(Q0[j])/RKJ[j]/RKJ[j];
B2[j] = CR1*(2.0f*Q0[j]*Q0[j]*BJ[j]/AJ[j]/AJ[j]-G*AJM+G*DZJ*BJ[j]);
temp = G*CR1*DX*Q0[j]*fabsf(Q0[j])/RKJ[j]/RKJ[j]*(BJ[j]-2*AJ[j]*DKDZJ[j]/RKJ[j]);
B2[j] = B2[j]+temp;
C2[j] = 1.0f+4.0f*CR1*Q0[j+1]/AJ[j+1]+2.0f*G*CR1*DX*AJ[j+1]*fabsf(Q0[j+1])/RKJ[j+1]/RKJ[j+1];
temp = -2.0f*Q0[j+1]*Q0[j+1]*BJ[j+1]/AJ[j+1]/AJ[j+1]+G*AJM+G*DZJ*BJ[j+1];
D2[j] = CR1*temp;
temp = BJ[j+1]-2.0f*AJ[j+1]*DKDZJ[j+1]/RKJ[j+1];
D2[j] = D2[j]+G*CR1*DX*Q0[j+1]*fabsf(Q0[j+1])/RKJ[j+1]/RKJ[j+1]*temp;
E2[j] = DT/DX*(-2.0f*Q0[j+1]*Q0[j+1]/AJ[j+1]+2.0f*Q0[j]*Q0[j]/AJ[j]-G*AJM*DZJ);
E2[j]=E2[j]-G*DT*(AJ[j+1]*Q0[j+1]*fabsf(Q0[j+1])/RKJ[j+1]/RKJ[j+1]+Q0[j]*fabsf(Q0[j])*AJ[j]/RKJ[j]/RKJ[j]);
}
float AFB1,AFB2,AHC2;
int JR,IGG;
//輸入邊界條件 上游
IGG=IK[0];
BOUND(BCUP,T,IBU,JS,IGG,1);
//THE FIRST SWEEP
for(j=JS;j<JEM1;j++)
{
AFB1=A1[j]*FJ[j]+B1[j];
H1[j]=-C1[j]/AFB1;
H2[j]=-D1[j]/AFB1;
H3[j]=(E1[j]-A1[j]*GJ[j])/AFB1;
AFB2=A2[j]*FJ[j]+B2[j];
AHC2=AFB2*H1[j]+C2[j];
FJ[j+1]=-(AFB2*H2[j]+D2[j])/AHC2;
GJ[j+1]=(E2[j]-AFB2*H3[j]-A2[j]*GJ[j])/AHC2;
}
//輸入邊界條件 下游
IGG=IK[1];
BOUND(BCDN,T,IBD,JE-1,IGG,2);
//THE SECOND SWEEP
for(j=JS;j<JEM1;j++)
{
JR=JE-j-2;
DZ[JR]=H1[JR]*DQ[JR+1]+H2[JR]*DZ[JR+1]+H3[JR]; //水深的增量
DQ[JR]=FJ[JR]*DZ[JR]+GJ[JR]; //流量的增量
}
//JE個斷數循環，重置參數
for(j=JS;j<JE;j++) //RESET INITIAL CONDITION
{
Q0[j]=Q0[j]+DQ[j]; //初識流量+流量增量
Z0[j]=Z0[j]+DZ[j]; //初始水深+水深增量
QN[N][j]=Q0[j]; //計算的流量
ZN[N][j]=Z0[j]; //計算的水深
}
}
//輸出 計算次數，時間點，水深，流量
for(N=0;N<NT;N++)
{
float TIME=float(N+1)*DT;
printf("\n%d\t%f\t\n\t\t",N+1,TIME);
for(i=0;i<JE;i++)
{
printf("斷面%d ",i);
}
printf("\n計算的流量\t");
for(i=0;i<JE;i++)
{
printf("%.3f ",QN[N][i]);
}
printf("\n計算的水深\t");
for(i=0;i<JE;i++)
{
printf("%.3f ",ZN[N][i]);
}
printf("\n");
}
return 0;
}?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -