cfd eno格式.cpp

CFD各种差分格式,ENO/Lax-Friedriches/Beam-Worming/Lax实例

#include<iostream.h>
#include<stdio.h>
#include<math.h>
#define T 0.5
#define h 1
void main()
{  
	int j,n;
    int	P=1000,M=100,G=80;
	int i;
	int k=2;
	double ub[1000][2];
	int is[1000];
	double u[1000];
	double uz[1000][100];
                    
	for(j=0;j<501;j++)         uz[j][0]=u[j]=1.0;     //初始条件离散;
	for(j=501;j<1000;j++)      uz[j][0]=u[j]=0.0;

    for(n=1;n<M-1;n++)
	{
		for(j=2*n-2;j<P-n+1;j++)
		 {
			 ub[j][0]=u[j];
		 }
		 for(i=1;i<=k-1;i++)
			 for(j=2*n-2;j<P-n;j++)
			 {
				 ub[j][i]=ub[j+1][i-1]-ub[j][i-1];
			 }
		 for(j=2*n;j<P-n;j++)
		 {
			 is[j]=j;
			 for(i=1;i<=k-1;i++)
			 {
				 if(fabs(ub[is[j]-1][i])<=fabs(ub[is[j]][i]))
					 is[j]=is[j]-1;
			 }
			 if(j-1==is[j])
			 {
    			 uz[j][n]=uz[j][n-1]-0.5*(0.5*uz[j-2][n-1]-2*uz[j-1][n-1]+1.5*uz[j][n-1]);
			 }
			 if(j==is[j])
			 {
				 uz[j][n]=uz[j][n-1]-0.5*(-0.5*uz[j-1][n-1]+0.5*uz[j+1][n-1]);
			 }
			 u[j]=uz[j][n];
		 }
	}
 
	FILE *fp;
	fp=fopen("d:\ENOTEST.plt","w");
	
	for(j=2*G;j<1000-G;j++)
	{   
	   fprintf(fp,"%d%c",j,' ');
	   fprintf(fp,"%f%c",uz[j][G],'\n');
	}
	fclose(fp);
}