csimple2d.c

CSIMPLE2d CFD C源程序

/*
* SIMPLE--Semi-Implicit Method for Pressure-Linked Equations--program
*								by Suhas V. Patankar
*
* This is a two demensional SIMPLE program in ANSI C language rewritten from its fortran versoin
*
* Copyright May 2000, rewritten by xu xianghua([e-mail])
*								mailto: [e-mail]
*								http://166.111.37.31
*
* This code is freedom, may be used in any form in any way you desire. This  
* file may be redistributed unmodified by any means PROVIDING it is   
* not sold for profit without the authors written consent, and   
* providing that this notice and the authors name and all copyright   
* notices remains intact.  
*  
* I will be very pleased if you find any bug or improve this code and inform me.
* Thank you! Good Luck!
*/
/*
* CSimple2D.c : main program of the SIMPLE 
*/
/*
* 注意C与FORTRAN的差别：
* 1. C中的数组下标是从0开始的,而FORTRAN则是从1开始的
*   因此本程序中的一些量(如L1,L2,L3,M1,M2,M3,NFMAX,nGam,nP,nRho,IPref,JPref)比FORTRAN程序中的小1
*   所有数组的下标从0开始,程序中数组下标中出现的数字比FORTRAN中的小1
* 2. C的数组是按行优先读的,而FORTRAN的则是按列读的
*   因此本程序中的F的定义不同,是F[NFMAX][I][J]
* 3. C用指针来实现FORTRAN的EQUIVALENCE功能。
*   如定义double (*U)[NSIZE],U=F[0],则U[i][j]等同于F[0][i][j];
*/
#include "staincl.h"
#include "CSimple2D.h"

/*********************************************************************************
* the main func
* input param:	output file name
*
*********************************************************************************/
int main(int argc,char* argv[])
{
	int i;

/* do some initialize */
	if(argc==2)
		szFileName=argv[1];
	else if(argc>2)
	{
		printf("too many parameters!\n");
		return 1;
	}
	else
		szFileName="Out.dat";

	nFMax=9;
	nP=10;
	nRho=11;
	nGam=12;

	bStop=FALSE;
	for(i=0;i<10;i++)
		bSolve[i]=FALSE;
	for(i=0;i<13;i++)
		bPrint[i]=FALSE;
	for(i=0;i<10;i++)
		bBlk[i]=TRUE;

	nMode=1;
	iLast=5;
	Time=0;
	iIter=0;
	for(i=0;i<13;i++)
		Relax[i]=1.0;
	for(i=0;i<10;i++)
		nTimes[i]=1;

	Dt=1e10;
	IPRef=0;
	JPRef=0;
	RhoCon=1.0;

	U=F[0];			//U[i][j]==F[0][i][j]
	V=F[1];			//V[i][j]==F[1][i][j]
	PC=F[2];		//PC[i][j]==F[2][i][j]
	P=F[nP];		//P[i][j]==F[nP][i][j]
	Rho=F[nRho];	//Rho[i][j]==F[nRho][i][j]
	Gam=F[nGam];	//Gam[i][j]==F[nGam][i][j]

	fnGrid();
	fnBegin();
	fnSetUp1();
	fnStart();
	/* do the main computation loop */
	while(1)
	{
		fnDense();
		fnBound();
		fnOutPut();
		if(bStop)
			break;
		fnSetUp2();
	}
	fnEnd();

	return 0;
}

/*********************************************************************************
* function name:	fnBegin
* description:		open the output file and get/print the current time
* input param:		-none-
* out param:		-none-
*********************************************************************************/
void fnBegin()
{
	if((pFileOut=fopen(szFileName,"wb"))==NULL)
	{
		printf("Cannot open the file %s.\n",szFileName);
		exit(0);
	}
	
	time(&start_time);

	printf("Computation begins at: \t\t%s", ctime( &start_time ));
	printf("Output file is %s.\n",szFileName);
	fprintf(pFileOut,"Computation begins at:\t\t%s", ctime( &start_time ));
}

/*********************************************************************************
* function name:	fnEnd
* description:		close the output file and get/print the current time,total 
*																	consumed time
* input param:		-none-
* out param:		-none-
*********************************************************************************/
void fnEnd()
{
	time_t stop_time;
	time_t span_time;

	time(&stop_time);

	span_time=stop_time-start_time;

	printf("\n\nComputation stops at: \t\t%s", ctime( &stop_time ));
	printf("\nComputation time is: \t\t%ld second.\n", span_time );
	fprintf(pFileOut,"\n\nComputation stops at: \t\t%s", ctime( &stop_time ));
	fprintf(pFileOut,"\nComputation time is: \t\t%ld second.\n", span_time );

	fclose(pFileOut);
}

/*********************************************************************************
* function name:	fnDiFlow
* description:		
* input param:		-none-
* out param:		-none-
*********************************************************************************/
void fnDiFlow()
{
	double temp;
	ACof=Diff;
	if(Flow!=0)
	{
		temp=Diff-fabs(Flow)*0.1;
		ACof=0;
		if(temp>0)
		{
			temp/=Diff;
			ACof=Diff*power(temp,5);
		}
	}
}

/*********************************************************************************
* function name:	fnSolve
* description:		solve the eqs
* input param:		-none-
* out param:		-none-
*********************************************************************************/
void fnSolve()
{
	int i,j,nT;
	double temp,denom;
	short ISTF,JSTF;
	double PT[NSIZE],QT[NSIZE];
	double BL,BLM,BLP,BLC;

	ISTF=IST-1;
	JSTF=JST-1;

	for(nT=0;nT<nTimes[nF];nT++)
	{
		if(bBlk[nF])
		{
//垂直方向块修正

			PT[ISTF]=0;
			QT[ISTF]=0;

			for(i=IST;i<=L2;i++)
			{
				BL=0;
				BLP=0;
				BLM=0;
				BLC=0;
				
				for(j=JST;j<=M2;j++)
				{
					BL+=AP[i][j];
					if(j!=M2)
						BL-=AJP[i][j];
					if(j!=JST)
						BL-=AJM[i][j];
					BLP+=AIP[i][j];
					BLM+=AIM[i][j];
					BLC+=AIP[i][j]*F[nF][i+1][j]+AIM[i][j]*F[nF][i-1][j]
						+AJP[i][j]*F[nF][i][j+1]+AJM[i][j]*F[nF][i][j-1]
						+Con[i][j]-AP[i][j]*F[nF][i][j];
				}
				denom=BL-PT[i-1]*BLM;
				if(fabs(denom/BL)<1e-10)
					denom=1e38;
				PT[i]=BLP/denom;
				QT[i]=(BLC+BLM*QT[i-1])/denom;
			}
			BL=0;
			for(i=L2;i>=IST;i--)
			{
				BL=BL*PT[i]+QT[i];
				for(j=JST;j<=M2;j++)
					F[nF][i][j]+=BL;
			}
//水平方向块修正
			PT[JSTF]=0;
			QT[JSTF]=0;
			for(j=JST;j<=M2;j++)
			{
				BL=0;
				BLP=0;
				BLM=0;
				BLC=0;
				for(i=IST;i<=L2;i++)
				{
					BL+=AP[i][j];
					if(i!=L2)
						BL-=AIP[i][j];
					if(i!=IST)
						BL-=AIM[i][j];
					BLP+=AJP[i][j];
					BLM+=AJM[i][j];
					BLC+=AIP[i][j]*F[nF][i+1][j]+AIM[i][j]*F[nF][i-1][j]
						+AJP[i][j]*F[nF][i][j+1]+AJM[i][j]*F[nF][i][j-1]
						+Con[i][j]-AP[i][j]*F[nF][i][j];
				}
				denom=BL-PT[j-1]*BLM;
			    if (fabs(denom/BL)<1e-10)
					denom=1e38;
				PT[j]=BLP/denom;
				QT[j]=(BLC+BLM*QT[j-1])/denom;
			}
			BL=0;
			for(j=M2;j>=JST;j--)
			{
				BL=BL*PT[j]+QT[j];
				for(i=IST;i<=L2;i++)
					F[nF][i][j]+=BL;
			}
		}
 //j=JST to M2 TDMA
		for(j=JST;j<=M2;j++)
		{
			PT[ISTF]=0;
			QT[ISTF]=F[nF][ISTF][j];
			for(i=IST;i<=L2;i++)
			{
				denom=AP[i][j]-PT[i-1]*AIM[i][j];
				PT[i]=AIP[i][j]/denom;
				temp=Con[i][j]+AJP[i][j]*F[nF][i][j+1]
					+AJM[i][j]*F[nF][i][j-1];
				QT[i]=(temp+AIM[i][j]*QT[i-1])/denom;
			}
			for(i=L2;i>=IST;i--)
				F[nF][i][j]=F[nF][i+1][j]*PT[i]+QT[i];
		}
 //j=M3 to JST TDMA
		for(j=M3;j>=JST;j--)
		{
			PT[ISTF]=0;
			QT[ISTF]=F[nF][ISTF][j];
			for(i=IST;i<=L2;i++)
			{
				denom=AP[i][j]-PT[i-1]*AIM[i][j];
				PT[i]=AIP[i][j]/denom;
				temp=Con[i][j]+AJP[i][j]*F[nF][i][j+1]
					+AJM[i][j]*F[nF][i][j-1];
				QT[i]=(temp+AIM[i][j]*QT[i-1])/denom;
			}
			for(i=L2;i>=IST;i--)
				F[nF][i][j]=F[nF][i+1][j]*PT[i]+QT[i];
		}
 //i=IST to L2 TDMA
		for(i=IST;i<=L2;i++)
		{
			PT[JSTF]=0;
			QT[JSTF]=F[nF][i][JSTF];
			for(j=JST;j<=M2;j++)
			{
				denom=AP[i][j]-PT[j-1]*AJM[i][j];
				PT[j]=AJP[i][j]/denom;
				temp=Con[i][j]+AIP[i][j]*F[nF][i+1][j]
					+AIM[i][j]*F[nF][i-1][j];
				QT[j]=(temp+AJM[i][j]*QT[j-1])/denom;
			}
			for(j=M2;j>=JST;j--)
				F[nF][i][j]=F[nF][i][j+1]*PT[j]+QT[j];
		}
 //i=L3 to IST TDMA
		for(i=L3;i>=IST;i--)
		{
			PT[JSTF]=0;
			QT[JSTF]=F[nF][i][JSTF];
			for(j=JST;j<=M2;j++)
			{
				denom=AP[i][j]-PT[j-1]*AJM[i][j];
				PT[j]=AJP[i][j]/denom;
				temp=Con[i][j]+AIP[i][j]*F[nF][i+1][j]
					+AIM[i][j]*F[nF][i-1][j];
				QT[j]=(temp+AJM[i][j]*QT[j-1])/denom;
			}
			for(j=M2;j>=JST;j--)
				F[nF][i][j]=F[nF][i][j+1]*PT[j]+QT[j];
		}
	}

	for(i=1;i<=L2;i++)
	{
		for(j=1;j<=M2;j++)
		{
			Con[i][j]=0;
			AP[i][j]=0;
		}
	}
}

/*********************************************************************************
* function name:	fnSetUp1
* description:		solve the eqs
* input param:		-none-
* out param:		-none-
*********************************************************************************/
void fnSetUp1()
{
	int i,j;

	L2=L1-1;
    L3=L2-1;
    M2=M1-1;
    M3=M2-1;

    X[0]=XU[1];
    for(i=1;i<=L2;i++)
		X[i]=0.5*(XU[i+1]+XU[i]);
    X[L1]=XU[L1];

    Y[0]=YV[1];
    for(j=1;j<=M2;j++)
		Y[j]=0.5*(YV[j+1]+YV[j]);
    Y[M1]=YV[M1];

    for(i=1;i<=L1;i++)
		XDif[i]=X[i]-X[i-1];
    for(i=1;i<=L2;i++)
		XCV[i]=XU[i+1]-XU[i];
    for(i=2;i<=L2;i++)
		XCVS[i]=XDif[i];
    XCVS[2]=XCVS[2]+XDif[1];
    XCVS[L2]=XCVS[L2]+XDif[L1];
    for(i=2;i<=L3;i++)
	{
		XCVI[i]=0.5*XCV[i];
		XCVIP[i]=XCVI[i];
	}
    XCVIP[1]=XCV[1];
    XCVI[L2]=XCV[L2];

    for(j=1;j<=M1;j++)
		YDif[j]=Y[j]-Y[j-1];
    for(j=1;j<=M2;j++)
		YCV[j]=YV[j+1]-YV[j];
    for(j=2;j<=M2;j++)
		YCVS[j]=YDif[j];
    YCVS[2]+=YDif[1];
    YCVS[M2]+=YDif[M1];

    if(nMode==1)
	{
		for(j=0;j<=M1;j++)
		{
			RMN[j]=1;
			R[j]=1;
		}
	}
    else
	{
		for(j=1;j<=M1;j++)
			R[j]=R[j-1]+YDif[j];
		RMN[1]=R[0];
		for(j=2;j<=M2;j++)
			RMN[j]=RMN[j-1]+YCV[j-1];
		RMN[M1]=R[M1];
	}

    for(j=0;j<=M1;j++)
	{
		SX[j]=1;
		SXMN[j]=1;
		if(nMode==3)
		{
			SX[j]=R[j];
			if(j!=0)
				SXMN[j]=RMN[j];
		}
	}

    for(j=1;j<=M2;j++)
	{
		YCVR[j]=R[j]*YCV[j];
		ARX[j]=YCVR[j];
		if(nMode==3)
			ARX[j]=YCV[j];
	}

    for(j=3;j<=M3;j++)
		YCVRS[j]=0.5*(R[j]+R[j-1])*YDif[j];

    YCVRS[2]=0.5*(R[2]+R[0])*YCVS[2];
    YCVRS[M2]=0.5*(R[M1]+R[M3])*YCVS[M2];

    if(nMode==2)
	{
		for(j=2;j<=M3;j++)
		{
			ARXJ[j]=0.25*(1+RMN[j]/R[j])*ARX[j];
			ARXJP[j]=ARX[j]-ARXJ[j];
		}
	}
    else
	{
		for(j=2;j<=M3;j++)
		{
			ARXJ[j]=0.5*ARX[j];
			ARXJP[j]=ARXJ[j];
		}
	}

	ARXJP[1]=ARX[1];
	ARXJ[M2]=ARX[M2];

	for(j=2;j<=M3;j++)
	{
		FV[j]=ARXJP[j]/ARX[j];
		FVP[j]=1-FV[j];
	}

	for(i=2;i<=L2;i++)
	{
		FX[i]=0.5*XCV[i-1]/XDif[i];
		FXM[i]=1-FX[i];
	}

	FX[1]=0;
    FXM[1]=1;
    FX[L1]=1;
    FXM[L1]=0;

    for(j=2;j<=M2;j++)
	{
		FY[j]=0.5*YCV[j-1]/YDif[j];
		FYM[j]=1-FY[j];
	}

    FY[1]=0;
    FYM[1]=1;
    FY[M1]=1;
    FYM[M1]=0;

//CON,AP,U,V,RHO,PC AND P ARRAYS ARE INITIALIZED HERE
	for(i=0;i<=L1;i++)
	{
		for(j=0;j<=M1;j++)
		{
			PC[i][j]=0;
			U[i][j]=0;
			V[i][j]=0;
			Con[i][j]=0;
			AP[i][j]=0;
			Rho[i][j]=RhoCon;
			P[i][j]=0.0;
		}
	}
	
	switch(nMode)
	{
	case 1:
		fprintf(pFileOut,"               computation in cartesian coordinates\n");
		break;
	case 2: