2013-04-11之前的版本
求解压力poisson方程
View Code
#include <iomanip> #include <fstream> #include "matvec.h" #include"GS.h" #include"GMRES.h" #include"ggje.h" #include "callumfpack.h" //2013-2-22 核查了A和B是否于Poisson离散出来形式一致 //加入了左边界条件为强制边界(Dirichlet条件) //#define Temp(i,j) Temp[(i)*(N+1)*(M+1)+j] //使用宏定义将二维转换成一维 //2013-4-2 添加了umfpack计算压力 void Sol_P(int N,int M,double delta_x,double delta_y,double delta_t ,double **ustar ,double **vstar ,double **upie ,double **vpie,double **Pres) { //void (*p)(double * ,int , double [],double []); //定义一个函数指针变量 p //cout<<"N="<<N<<" "<<"M="<<M<<endl; double **A=new double *[(N+1)*(M+1)]; for(int k=0;k<(N+1)*(M+1);k++) A[k]= new double [(N+1)*(M+1)]; for(int i=0;i<(N+1)*(M+1);i++) for(int j=0;j<(N+1)*(M+1);j++) A[i][j]=0; //cout<<endl; // 输出初始化的 A矩阵 // cout<<" 初始化A矩阵 "<<endl; /* for(int i=0;i<(N+1)*(M+1);i++) { for(int j=0;j<(N+1)*(M+1);j++) { cout<<A[i][j]<<" "; } cout<<endl; }*/ double *B =new double [(N+1)*(M+1)]; for(int bk=0;bk<(N+1)*(M+1);bk++) { B[bk]=0; // cout<<B[bk]<<" "; } //cout<<endl; double *sol_P =new double [(N+1)*(M+1)]; for(int bk=0;bk<(N+1)*(M+1);bk++) { sol_P[bk]=0; // cout<<sol_P[bk]<<" "; } //cout<<endl; /*---------------------------------------------------------------------------------------------*/ /*点(0,0)处的离散方程*/ A[0*(M+1)+0][0*(M+1)+0]=-2*(1/pow(delta_x,2)+1/pow(delta_y,2)); A[0*(M+1)+0][0*(M+1)+1]=2*(1/pow(delta_y,2)); A[0*(M+1)+0][1*(M+1)+0]=2*(1/pow(delta_x,2)); B[0*(M+1)+0]=1/delta_t*(ustar[0+1][0]-ustar[0][0])/delta_x+1/delta_t*(vstar[0][0+1]-vstar[0][0])/delta_y+2/delta_x*(-(upie[0][0]-ustar[0][0])/delta_t)+2/delta_y*(-(vpie[0][0]-vstar[0][0])/delta_t); // 向后差分 /*i=0上的点（除去(0,0),(0,M)）*/ for(int j=1;j<=M-1;j++) { int i=0;//第一行 A[i*(M+1)+j][i*(M+1)+j-1]=1/pow(delta_y,2); A[i*(M+1)+j][i*(M+1)+j]=-2*(1/pow(delta_x,2)+1/pow(delta_y,2)); A[i*(M+1)+j][i*(M+1)+j+1]=1/pow(delta_y,2); A[i*(M+1)+j][(i+1)*(M+1)+j]=2.0/pow(delta_x,2); B[i*(M+1)+j]=1/delta_t*(ustar[i+1][j]-ustar[i][j])/delta_x+1/delta_t*(vstar[i][j+1]-vstar[i][j-1])/(2*delta_y)+2/delta_x*(-(upie[i][j]-ustar[i][j])/delta_t); // 向后差分和中心差分（能用中心差分就用中心差分） } /*点(0,M)处的离散方程*/ A[0*(M+1)+M][0*(M+1)+M-1]=2*(1/pow(delta_y,2)); A[0*(M+1)+M][0*(M+1)+M]=-2*(1/pow(delta_x,2)+1/pow(delta_y,2)); A[0*(M+1)+M][1*(M+1)+M]=2*(1/pow(delta_x,2)); B[0*(M+1)+M]=1/delta_t*(ustar[0+1][M]-ustar[0][M])/delta_x+1/delta_t*(vstar[0][M]-vstar[0][M-1])/delta_y+2/delta_x*(-(upie[0][M]-ustar[0][M])/delta_t)-2/delta_y*(-(vpie[0][M]-vstar[0][M])/delta_t); /*i=1~N-1,j=1~M-1 之间的离散方程*/ for(int i=1;i<N;i++) { for(int j=1;j<M;j++) { A[i*(M+1)+j][i*(M+1)+j-1]=1*(1/pow(delta_y,2)); A[i*(M+1)+j][i*(M+1)+j]=-2*(1/pow(delta_x,2)+1/pow(delta_y,2)); A[i*(M+1)+j][i*(M+1)+j+1]=1*(1/pow(delta_y,2)); A[i*(M+1)+j][(i-1)*(M+1)+j]=1*(1/pow(delta_x,2)); A[i*(M+1)+j][(i+1)*(M+1)+j]=1*(1/pow(delta_x,2)); B[i*(M+1)+j]=1/delta_t*((ustar[i+1][j]-ustar[i-1][j])/(2*delta_x)+(vstar[i][j+1]-vstar[i][j-1])/(2*delta_y));//对一阶导数采用中心差分 } } /*------------------------------------------------------------------------------------------------------*/ /*点(N,0)处的离散方程*/ A[N*(M+1)+0][N*(M+1)+0]=-2*(1/pow(delta_x,2)+1/pow(delta_y,2)); A[N*(M+1)+0][N*(M+1)+1]=2*(1/pow(delta_y,2)); A[N*(M+1)+0][(N-1)*(M+1)+0]=2*(1/pow(delta_x,2)); B[N*(M+1)+0]=1/delta_t*(ustar[N][0]-ustar[N-1][0])/delta_x+1/delta_t*(vstar[N][0+1]-vstar[N][0])/delta_y-2/delta_x*(-(upie[N][0]-ustar[N][0])/delta_t)+2/delta_y*(-(vpie[N][0]-vstar[N][0])/delta_t); /*i=N上的点（除去(0,0),(0,M)）*/ for(int j=1;j<=M-1;j++) { int i=N;//第N行 A[i*(M+1)+j][i*(M+1)+j-1]=1/pow(delta_y,2); A[i*(M+1)+j][i*(M+1)+j]=-2*(1/pow(delta_x,2)+1/pow(delta_y,2)); A[i*(M+1)+j][i*(M+1)+j+1]=1/pow(delta_y,2); A[i*(M+1)+j][(i-1)*(M+1)+j]=2.0/pow(delta_x,2); B[i*(M+1)+j]=1/delta_t*(ustar[i][j]-ustar[i-1][j])/delta_x+1/delta_t*(vstar[i][j+1]-vstar[i][j-1])/(2*delta_y)-2/delta_x*(-(upie[i][j]-ustar[i][j])/delta_t); } /*点(N,M)处的离散方程*/ A[N*(M+1)+M][N*(M+1)+M-1]=2*(1/pow(delta_y,2)); A[N*(M+1)+M][N*(M+1)+M]=-2*(1/pow(delta_x,2)+1/pow(delta_y,2)); A[N*(M+1)+M][(N-1)*(M+1)+M]=2*(1/pow(delta_x,2)); B[N*(M+1)+M]=1/delta_t*(ustar[N][M]-ustar[N-1][M])/delta_x+1/delta_t*(vstar[N][M]-vstar[N][M-1])/delta_y-2/delta_x*(-(upie[N][M]-ustar[N][M])/delta_t)-2/delta_y*(-(vpie[N][M]-vstar[N][M])/delta_t); //cout<<"B[N*(M+1)+M]="<<B[N*(M+1)+M]<<endl; /*---------------------------------------------------------------------------------------------------------------------*/ /*j=M上的点（除去(0,M),(N,M)）*/ for(int i=1;i<=N-1;i++) { int j=M;//第N行 A[i*(M+1)+j][(i-1)*(M+1)+j]=1/pow(delta_y,2); A[i*(M+1)+j][i*(M+1)+j]=-2*(1/pow(delta_x,2)+1/pow(delta_y,2)); A[i*(M+1)+j][i*(M+1)+j-1]=2*(1/pow(delta_y,2)); A[i*(M+1)+j][(i+1)*(M+1)+j]=1.0/pow(delta_x,2); B[i*(M+1)+j]=1/delta_t*(ustar[i+1][j]-ustar[i-1][j])/(2*delta_x)+1/delta_t*(vstar[i][j]-vstar[i][j-1])/delta_y-2/delta_y*(-(vpie[i][j]-vstar[i][j])/delta_t); } /*---------------------------------------------------------------------------------------------------------------------*/ /*j=0上的点（除去(0,0),(N,0)）*/ for(int i=1;i<=N-1;i++) { int j=0;//第0行 A[i*(M+1)+j][(i-1)*(M+1)+j]=1/pow(delta_x,2); A[i*(M+1)+j][i*(M+1)+j]=-2*(1/pow(delta_x,2)+1/pow(delta_y,2)); A[i*(M+1)+j][i*(M+1)+j+1]=2*(1/pow(delta_y,2)); A[i*(M+1)+j][(i+1)*(M+1)+j]=1.0/pow(delta_x,2); B[i*(M+1)+j]=1/delta_t*(ustar[i+1][j]-ustar[i-1][j])/(2*delta_x)+1/delta_t*(vstar[i][j+1]-vstar[i][j])/delta_y+2/delta_y*(-(vpie[i][j]-vstar[i][j])/delta_t); } /*---------------------------------------------------------------------------------------------------------------------*/ //double *Temp=new double[((N+1)*(M+1))*((N+1)*(M+1))]; //动态分配内存空间用于存储临时矩阵 ; //将矩阵Temp全部元素初始化为 //注：row为行号， col为列号 /*-----------------------处理Ax=b中A的形式让其利用Gauss或者稀疏可解-----------------*/ //???????????????????????????????? for(int i=0;i<M+1;i++) { for(int j=0;j<(N+1)*(M+1);j++) { A[i][j]=0; } A[i][i]=1; B[i]=10;//左边界赋值为常数10 } //ofstream outAMatrix; //outAMatrix.open("F:\\fluid\\C-V-IBM\\AMatrix.txt"); ////速度太慢 //for(int row=0;row<(N+1)*(M+1);row++) //{ // for(int col=0;col<(N+1)*(M+1);col++) // { // //Temp(row,col)=A[row][col]; // outAMatrix<<setprecision(5)<<A[row][col]<<" "; // } // outAMatrix<<endl; //} /*ofstream outBVector; outBVector.open("F:\\fluid\\C-V-IBM\\BVector.txt"); for(int i=0;i<(N+1)*(M+1);i++) { outBVector<<setprecision(5)<<B[i]<<" "; outBVector<<endl; }*/ double *testx =new double [(N+1)*(M+1)]; for(int bk=0;bk<(N+1)*(M+1);bk++) { testx[bk]=0; // cout<<testx[bk]<<" "; } //cout<<endl; /*---------------------调用UMFPACK求解 A*testx=B ---------------------------*/ umf(A,B,testx,(N+1)*(M+1)); for(int row=0;row<=N;row++) { for(int col=0;col<=M;col++) { //Pres[row][col]=sol_P[(row-1)*M+col-1]; //Pres represents Pressure Pres[row][col]=testx[row*M+col];// 调用umf 之前写成row*(M)+col 出现问题 } } std::ofstream outPres; outPres.open("F:\\fluid\\C-V-IBM\\Pres.txt"); //cout<< "输出求解的矩阵Pres"<<endl; for(int i=0;i<N+1;i++) { for(int j=0;j<M+1;j++) { //cout<<Pres[i][j]<< " "; outPres<<setprecision(5)<<Pres[i][j]<<" "; } //cout<<endl; outPres<<endl; } for(int t=0;t<(N+1)*(M+1);t++) delete []A[t]; delete []A; /*撤销u*空间*/ delete [] B; delete [] sol_P; outPres.close(); //outAMatrix.close(); }
这个有问题 主要有
1 压力不对称
2 压力输出有问题
Pres[row][col]=testx[row*M+col];的错误 困扰了很久
修改版本
View Code
#include <iomanip> #include <fstream> #include "matvec.h" #include"GS.h" #include"GMRES.h" #include"ggje.h" #include "callumfpack.h" //2013-2-22 核查了A和B是否于Poisson离散出来形式一致 //加入了左边界条件为强制边界(Dirichlet条件) //#define Temp(i,j) Temp[(i)*(N+1)*(M+1)+j] //使用宏定义将二维转换成一维 //2013-4-2 添加了umfpack计算压力 //2013-4-8 修改了在边界处离散的错误，以前为在每个边界按照一个方向差分 pian_u/pian_x=(u_n+1-u_n-1)/delta_x //在边界处的差分都是由边界内部节点到外部节点 //2013-4-9 修改了pian_ustar/pian_x的迎风格式 按照具体物理模型进行离散 //2013-4-11 修改了压力输出的错误 以前导致有两条斜着的带子出现，错位导致 void Sol_P(int N,int M,double delta_x,double delta_y,double delta_t ,double **ustar ,double **vstar ,double **upie ,double **vpie,double **Pres) { //void (*p)(double * ,int , double [],double []); //定义一个函数指针变量 p //cout<<"N="<<N<<" "<<"M="<<M<<endl; double **A=new double *[(N+1)*(M+1)]; for(int k=0;k<(N+1)*(M+1);k++) A[k]= new double [(N+1)*(M+1)]; for(int i=0;i<(N+1)*(M+1);i++) for(int j=0;j<(N+1)*(M+1);j++) A[i][j]=0; //cout<<endl; // 输出初始化的 A矩阵 // cout<<" 初始化A矩阵 "<<endl; /* for(int i=0;i<(N+1)*(M+1);i++) { for(int j=0;j<(N+1)*(M+1);j++) { cout<<A[i][j]<<" "; } cout<<endl; }*/ double *B =new double [(N+1)*(M+1)]; for(int bk=0;bk<(N+1)*(M+1);bk++) { B[bk]=0; // cout<<B[bk]<<" "; } //cout<<endl; double *sol_P =new double [(N+1)*(M+1)]; for(int bk=0;bk<(N+1)*(M+1);bk++) { sol_P[bk]=0; // cout<<sol_P[bk]<<" "; } //cout<<endl; /*---------------------------------------------------------------------------------------------*/ /*点(0,0)处的离散方程*/ A[0*(M+1)+0][0*(M+1)+0]=-2*(1/pow(delta_x,2)+1/pow(delta_y,2)); A[0*(M+1)+0][0*(M+1)+1]=2*(1/pow(delta_y,2)); A[0*(M+1)+0][1*(M+1)+0]=2*(1/pow(delta_x,2)); B[0*(M+1)+0]=1/delta_t*(ustar[0+1][0]-ustar[0][0])/delta_x+1/delta_t*(vstar[0][0+1]-vstar[0][0])/delta_y+2/delta_x*(-(upie[0][0]-ustar[0][0])/delta_t)+2/delta_y*(-(vpie[0][0]-vstar[0][0])/delta_t); // 向后差分 /*i=0上的点（除去(0,0),(0,M)）*/ for(int j=1;j<=M-1;j++) { int i=0;//第一行 A[i*(M+1)+j][i*(M+1)+j-1]=1/pow(delta_y,2); A[i*(M+1)+j][i*(M+1)+j]=-2*(1/pow(delta_x,2)+1/pow(delta_y,2)); A[i*(M+1)+j][i*(M+1)+j+1]=1/pow(delta_y,2); A[i*(M+1)+j][(i+1)*(M+1)+j]=2.0/pow(delta_x,2); B[i*(M+1)+j]=1/delta_t*(ustar[i+1][j]-ustar[i][j])/delta_x+1/delta_t*(vstar[i][j+1]-vstar[i][j-1])/(2*delta_y)+2/delta_x*(-(upie[i][j]-ustar[i][j])/delta_t); // 向后差分和中心差分（能用中心差分就用中心差分） } /*点(0,M)处的离散方程*/ A[0*(M+1)+M][0*(M+1)+M-1]=2*(1/pow(delta_y,2)); A[0*(M+1)+M][0*(M+1)+M]=-2*(1/pow(delta_x,2)+1/pow(delta_y,2)); A[0*(M+1)+M][1*(M+1)+M]=2*(1/pow(delta_x,2)); B[0*(M+1)+M]=1/delta_t*(ustar[0+1][M]-ustar[0][M])/delta_x+1/delta_t*(vstar[0][M]-vstar[0][M-1])/delta_y+2/delta_x*(-(upie[0][M]-ustar[0][M])/delta_t)-2/delta_y*(-(vpie[0][M]-vstar[0][M])/delta_t); /*i=1~N-1,j=1~M-1 之间的离散方程*/ for(int i=1;i<N;i++) { for(int j=1;j<M;j++) { A[i*(M+1)+j][i*(M+1)+j-1]=1*(1/pow(delta_y,2)); A[i*(M+1)+j][i*(M+1)+j]=-2*(1/pow(delta_x,2)+1/pow(delta_y,2)); A[i*(M+1)+j][i*(M+1)+j+1]=1*(1/pow(delta_y,2)); A[i*(M+1)+j][(i-1)*(M+1)+j]=1*(1/pow(delta_x,2)); A[i*(M+1)+j][(i+1)*(M+1)+j]=1*(1/pow(delta_x,2)); B[i*(M+1)+j]=1/delta_t*((ustar[i+1][j]-ustar[i-1][j])/(2*delta_x)+(vstar[i][j+1]-vstar[i][j-1])/(2*delta_y));//对一阶导数采用中心差分 } } /*------------------------------------------------------------------------------------------------------*/ /*点(N,0)处的离散方程*/ A[N*(M+1)+0][N*(M+1)+0]=-2*(1/pow(delta_x,2)+1/pow(delta_y,2)); A[N*(M+1)+0][N*(M+1)+1]=2*(1/pow(delta_y,2)); A[N*(M+1)+0][(N-1)*(M+1)+0]=2*(1/pow(delta_x,2)); B[N*(M+1)+0]=1/delta_t*(ustar[N][0]-ustar[N-1][0])/delta_x+1/delta_t*(vstar[N][0+1]-vstar[N][0])/delta_y-2/delta_x*(-(upie[N][0]-ustar[N][0])/delta_t)+2/delta_y*(-(vpie[N][0]-vstar[N][0])/delta_t);//改变下边界符号 /*i=N上的点（除去(0,0),(0,M)）*/ for(int j=1;j<=M-1;j++) { int i=N;//第N行 A[i*(M+1)+j][i*(M+1)+j-1]=1/pow(delta_y,2); A[i*(M+1)+j][i*(M+1)+j]=-2*(1/pow(delta_x,2)+1/pow(delta_y,2)); A[i*(M+1)+j][i*(M+1)+j+1]=1/pow(delta_y,2); A[i*(M+1)+j][(i-1)*(M+1)+j]=2.0/pow(delta_x,2); B[i*(M+1)+j]=1/delta_t*(ustar[i][j]-ustar[i-1][j])/delta_x+1/delta_t*(vstar[i][j+1]-vstar[i][j-1])/(2*delta_y)-2/delta_x*(-(upie[i][j]-ustar[i][j])/delta_t); } /*点(N,M)处的离散方程*/ A[N*(M+1)+M][N*(M+1)+M-1]=2*(1/pow(delta_y,2)); A[N*(M+1)+M][N*(M+1)+M]=-2*(1/pow(delta_x,2)+1/pow(delta_y,2)); A[N*(M+1)+M][(N-1)*(M+1)+M]=2*(1/pow(delta_x,2)); B[N*(M+1)+M]=1/delta_t*(ustar[N][M]-ustar[N-1][M])/delta_x+1/delta_t*(vstar[N][M]-vstar[N][M-1])/delta_y-2/delta_x*(-(upie[N][M]-ustar[N][M])/delta_t)-2/delta_y*(-(vpie[N][M]-vstar[N][M])/delta_t); //cout<<"B[N*(M+1)+M]="<<B[N*(M+1)+M]<<endl; /*---------------------------------------------------------------------------------------------------------------------*/ /*j=M上的点（除去(0,M),(N,M)）*/ for(int i=1;i<=N-1;i++) { int j=M;//第N行 A[i*(M+1)+j][(i-1)*(M+1)+j]=1/pow(delta_y,2); A[i*(M+1)+j][i*(M+1)+j]=-2*(1/pow(delta_x,2)+1/pow(delta_y,2)); A[i*(M+1)+j][i*(M+1)+j-1]=2*(1/pow(delta_y,2)); A[i*(M+1)+j][(i+1)*(M+1)+j]=1.0/pow(delta_x,2); B[i*(M+1)+j]=1/delta_t*(ustar[i+1][j]-ustar[i-1][j])/(2*delta_x)+1/delta_t*(vstar[i][j]-vstar[i][j-1])/delta_y-2/delta_y*(-(vpie[i][j]-vstar[i][j])/delta_t); } /*---------------------------------------------------------------------------------------------------------------------*/ /*j=0上的点（除去(0,0),(N,0)）*/ for(int i=1;i<=N-1;i++) { int j=0;//第0行 A[i*(M+1)+j][(i-1)*(M+1)+j]=1/pow(delta_x,2); A[i*(M+1)+j][i*(M+1)+j]=-2*(1/pow(delta_x,2)+1/pow(delta_y,2)); A[i*(M+1)+j][i*(M+1)+j+1]=2*(1/pow(delta_y,2)); A[i*(M+1)+j][(i+1)*(M+1)+j]=1.0/pow(delta_x,2); B[i*(M+1)+j]=1/delta_t*(ustar[i+1][j]-ustar[i-1][j])/(2*delta_x)+1/delta_t*(vstar[i][j+1]-vstar[i][j])/delta_y+2/delta_y*(-(vpie[i][j]-vstar[i][j])/delta_t);//改变下边界 } /*---------------------------------------------------------------------------------------------------------------------*/ //double *Temp=new double[((N+1)*(M+1))*((N+1)*(M+1))]; //动态分配内存空间用于存储临时矩阵 ; //将矩阵Temp全部元素初始化为 //注：row为行号， col为列号 //下面八行需要注意---需满足delta_x=delta_y for(int i=0;i<=N;i++) { for(int j=0;j<=M;j++) { A[i][j]=pow(delta_x,2.0)*A[i][j]; B[i*(M+1)+j]=pow(delta_x,2.0)*B[i*(M+1)+j]; } } /*-----------------------处理Ax=b中A的形式让其利用Gauss或者稀疏可解-----------------*/ //???????????????????????????????? for(int i=0;i<M+1;i++) { for(int j=0;j<(N+1)*(M+1);j++) { A[i][j]=0; } A[i][i]=1; B[i]=10;//左边界赋值为常数10 } //ofstream outAMatrix; //outAMatrix.open("F:\\fluid\\cC-V-IBM\\AMatrix.txt"); ////速度太慢 //for(int row=0;row<(N+1)*(M+1);row++) //{ // for(int col=0;col<(N+1)*(M+1);col++) // { // //Temp(row,col)=A[row][col]; // outAMatrix<<setprecision(5)<<A[row][col]<<" "; // } // outAMatrix<<endl; //} ofstream outBVector; outBVector.open("F:\\fluid\\cC-V-IBM\\BVector.txt"); /*for(int i=0;i<(N+1)*(M+1);i++) { outBVector<<setprecision(5)<<B[i]<<" "; outBVector<<endl; }*/ for(int j=M;j>=0;j--) { for(int i=0;i<N+1;i++) { //cout<<" "<<u[i][j]; outBVector<<" "<<B[i*(M+1)+j]; //outBVector<<endl; } //cout<<endl; outBVector<<endl; } outBVector.close(); double *testx =new double [(N+1)*(M+1)]; for(int bk=0;bk<(N+1)*(M+1);bk++) { testx[bk]=0; // cout<<testx[bk]<<" "; } //cout<<endl; /*---------------------调用UMFPACK求解 A*testx=B ---------------------------*/ umf(A,B,testx,(N+1)*(M+1)); for(int row=0;row<=N;row++) { for(int col=0;col<=M;col++) { //Pres[row][col]=sol_P[(row-1)*M+col-1]; //Pres represents Pressure Pres[row][col]=testx[row*(M+1)+col];// 调用umf } } std::ofstream outPres; outPres.open("F:\\fluid\\cC-V-IBM\\Pres.txt"); //cout<< "输出求解的矩阵Pres"<<endl; for(int i=0;i<N+1;i++) { for(int j=0;j<M+1;j++) { //cout<<Pres[i][j]<< " "; outPres<<setprecision(5)<<Pres[i][j]<<" "; } //cout<<endl; outPres<<endl; } for(int t=0;t<(N+1)*(M+1);t++) delete []A[t]; delete []A; /*撤销u*空间*/ delete [] B; delete [] sol_P; outPres.close(); //outAMatrix.close(); }