solve.c

这是一个用SIMPLE法模拟的矩形流场的程序

#include "global.h"
void solve()
{	
	double temp[maxnum],at[maxnum],bt[maxnum],ct[maxnum],con[maxnum],pt[maxnum],qt[maxnum];//三对角方程系数
	int k=0;
	//从左到右扫描
		
	do{
	//从左到右扫描
	for(i=ist;i<=l-1;i++)
	{	
		pt[jst-1]=0;
		qt[jst-1]=0;
		//消元	
		for(j=jst;j<=m-1;j++)
		{	
			at[j]=ap[i][j][n];
			bt[j]=an[i][j][n];
			ct[j]=as[i][j][n];
			con[j]=b[i][j][n]+ae[i][j][n]*f[i+1][j][n]+aw[i][j][n]*f[i-1][j][n];
			if(j==jst){con[j]=con[j]+ct[j]*f[i][j-1][n];ct[j]=0;}
			if(j==m-1){con[j]=con[j]+bt[j]*f[i][j+1][n];bt[j]=0;}
			temp[j]=at[j]-ct[j]*pt[j-1];
			pt[j]=bt[j]/temp[j];
			qt[j]=(ct[j]*qt[j-1]+con[j])/temp[j];
			//test
	
//				printf("\nj=%d,%5.2f,%5.2f, %5.3f, %5.3f\n",j,ct[j],at[j],bt[j],con[j]);
			//test over
		}
		//回代
		f[i][m-1][n]=qt[m-1];
		for(j=m-2;j>=jst;j--)
			f[i][j][n]=pt[j]*f[i][j+1][n]+qt[j];
	
	}

	//从右到左扫描		
	for(i=l-1;i>=ist;i--)
	{	
		pt[jst-1]=0;
		qt[jst-1]=0;
		//消元	
		for(j=jst;j<=m-1;j++)
		{	
			at[j]=ap[i][j][n];
			bt[j]=an[i][j][n];
			ct[j]=as[i][j][n];
			con[j]=b[i][j][n]+ae[i][j][n]*f[i+1][j][n]+aw[i][j][n]*f[i-1][j][n];
			if(j==jst){con[j]=con[j]+ct[j]*f[i][j-1][n];ct[j]=0;}
			if(j==m-1){con[j]=con[j]+bt[j]*f[i][j+1][n];bt[j]=0;}
			temp[j]=at[j]-ct[j]*pt[j-1];
			pt[j]=bt[j]/temp[j];
			qt[j]=(ct[j]*qt[j-1]+con[j])/temp[j];
			//test
	
//				printf("\nj=%d,%5.2f,%5.2f, %5.3f, %5.3f\n",j,ct[j],at[j],bt[j],con[j]);
			//test over
		}
		//回代
		f[i][m-1][n]=qt[m-1];
		for(j=m-2;j>=jst;j--)
			f[i][j][n]=pt[j]*f[i][j+1][n]+qt[j];
	
	}

	//从下到上扫描

	for(j=jst;j<=m-1;j++)
	{	
		pt[ist-1]=0;
		qt[ist-1]=0;
		
		//消元
		for(i=ist;i<=l-1;i++)
		{	
			at[i]=ap[i][j][n];
			bt[i]=ae[i][j][n];
			ct[i]=aw[i][j][n];
			con[i]=b[i][j][n]+an[i][j][n]*f[i][j+1][n]+as[i][j][n]*f[i][j-1][n];
			if(i==ist){con[i]=con[i]+ct[i]*f[i-1][j][n];ct[i]=0;}
			if(i==l-1){con[i]=con[i]+bt[i]*f[i+1][j][n];bt[i]=0;}
			temp[i]=at[i]-ct[i]*pt[i-1];
			pt[i]=bt[i]/temp[i];
			qt[i]=(ct[i]*qt[i-1]+con[i])/temp[i];
			//test
	
//			printf("\ni=%d,%5.2f,%5.2f,%5.2f, %5.3f\n",i,ct[i],at[i],bt[i],f[i][j+1][n]);
			//test over
		}
		//回代
		f[l-1][j][n]=qt[l-1];
		for(i=l-2;i>=ist;i--)
			f[i][j][n]=pt[i]*f[i+1][j][n]+qt[i];
	
	}

	//从上到下扫描
	for(j=m-1;j>=jst;j--)
	{	
		pt[ist-1]=0;
		qt[ist-1]=0;
		
		//消元
		for(i=ist;i<=l-1;i++)
		{	
			at[i]=ap[i][j][n];
			bt[i]=ae[i][j][n];
			ct[i]=aw[i][j][n];
			con[i]=b[i][j][n]+an[i][j][n]*f[i][j+1][n]+as[i][j][n]*f[i][j-1][n];
			if(i==ist){con[i]=con[i]+ct[i]*f[i-1][j][n];ct[i]=0;}
			if(i==l-1){con[i]=con[i]+bt[i]*f[i+1][j][n];bt[i]=0;}
			temp[i]=at[i]-ct[i]*pt[i-1];
			pt[i]=bt[i]/temp[i];
			qt[i]=(ct[i]*qt[i-1]+con[i])/temp[i];
			//test
	
//				printf("\ni=%d,%5.2f,%5.2f, %5.3f, %5.3f\n",i,ct[i],at[i],bt[i],con[i]);
			//test over
		}
		//回代
		f[l-1][j][n]=qt[l-1];
		for(i=l-2;i>=ist;i--)
			f[i][j][n]=pt[i]*f[i+1][j][n]+qt[i];
	
	}
	
	k=k+1;
	}while(k<10);



}