nsesolver.cpp

PDE simulator on GPU.

// NSESolver.cpp: implementation of the NSESolver class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "shadow.h"
#include "NSESolver.h"

#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif

//////////////////////////////////////////////////////////////////////
// Try to solver impressible NS Equations ----2D
// Date: 2004/1/10
//////////////////////////////////////////////////////////////////////
#define XWidth 20
#define YHeight 20
#define TimeSteps  10
void CalculateNSEWithSimple()
{
	int i,j,k;
	float eps = 1.0e-5;
	float maxerror; 
	float u[TimeSteps][XWidth][YHeight];
	float v[TimeSteps][XWidth][YHeight];
	float p[TimeSteps][XWidth][YHeight];
	float dp[XWidth][YHeight];
	
	float dx = 1;
	float dy = 1;
	float dt = 0.1;
	float a1,a2,a3,a4;
	float b1,b2,b3,b4;
	float u1,u2;
	float v1,v2;
	float AA,BB;
	float Re = 500;
	
	float a = 2.0*dt*(1.0/(dx*dx) + 1.0/(dy*dy));
	float b = -dt/(dx*dx);
	float c = -dt/(dy*dy);
	float d;
	//////////////////////////////////////////////////////////////////////////
	/// initialize the field 
	memset(u, 0, sizeof(float)*TimeSteps*XWidth*YHeight);
	memset(v, 0, sizeof(float)*TimeSteps*XWidth*YHeight);
	memset(p, 0, sizeof(float)*TimeSteps*XWidth*YHeight);
	
	FILE *fp = fopen("Data/accurate&time/test.txt", "w");
	int m,n;
	for(k=0; k<TimeSteps; k++)
	{
		n=0;
loop:
	    n++;
		//////////////////////////////////////////////////////////////////////////
		/// Set the boundary condition here
		for(j=1; j<YHeight-1; j++)
		{
			p[k][0][j] = p[k][1][j];
			p[k][XWidth-1][j] = p[k][XWidth-2][j];

			//Left
			u[k][0][j] = -2.0*u[k][0][j-1]/3.0;
			//Right
			u[k][XWidth-1][j] = -2.0*u[k][XWidth-1][j-1]/3.0;

			v[k][0][j] = 0.0;
			v[k][XWidth-1][j] = 0.0;
		}
		for(i=1; i<XWidth-1; i++)
		{
			p[k][i][0] = p[k][i][1];
			p[k][i][YHeight-1] = p[k][i][YHeight-2];

			//Top
			u[k][i][0] = 1.0;
			//Bottom
			u[k][i][YHeight-1] = 0.0;

			v[k][i][0] = 0.0;
			v[k][i][YHeight-1] = 0.0;
		}
		p[k][0][0]					= p[k][1][1];
		p[k][0][YHeight-1]			= p[k][1][YHeight-2];
		p[k][XWidth-1][0]			= p[k][XWidth-2][1];
		p[k][XWidth-1][YHeight-1]	= p[k][XWidth-2][YHeight-2];

		u[k][0][0]					= 0.0;//u[k][1][1];
		u[k][0][YHeight-1]			= 0.0;//u[k][1][YHeight-2];
		u[k][XWidth-1][0]			= 0.0;//u[k][XWidth-2][1];
		u[k][XWidth-1][YHeight-1]	= 0.0;//u[k][XWidth-2][YHeight-2];

		v[k][0][0]					= 0.0;//v[k][1][1];
		v[k][0][YHeight-1]			= 0.0;//v[k][1][YHeight-2];
		v[k][XWidth-1][0]			= 0.0;//v[k][XWidth-2][1];
		v[k][XWidth-1][YHeight-1]	= 0.0;//v[k][XWidth-2][YHeight-2];
		
		for(j=1; j<YHeight-1; j++)
		{
			for(i=1; i<XWidth-1; i++)
			{
				u1 = 0.5*(u[k][i-1][j] + u[k][i-1][j+1]);
				u2 = 0.5*(u[k][i][j] + u[k][i][j+1]);

				v1 = 0.5*(v[k][i][j] + v[k][i+1][j]);
				v2 = 0.5*(v[k][i][j-1] + v[k][i+1][j-1]);

				a1 = 0.5*(u[k][i+1][j]*u[k][i+1][j] - u[k][i-1][j]*u[k][i-1][j])/dx - 0.5*(u[k][i][j+1]*v1 - u[k][i][j-1]*v2)/dy;
				b1 = 0.5*(v[k][i][j+1]*v[k][i][j+1] - v[k][i][j-1]*v[k][i][j-1])/dy - 0.5*(v[k][i+1][j]*u1 - v[k][i-1][j]*u2)/dx;

				a3 = (u[k][i+1][j] - 2.0*u[k][i][j] + u[k][i-1][j])/(dx*dx);
				a4 = (u[k][i][j+1] - 2.0*u[k][i][j] + u[k][i][j-1])/(dy*dy);

				b3 = (v[k][i][j+1] - 2.0*v[k][i][j] + v[k][i][j-1])/(dy*dy);
				b4 = (v[k][i+1][j] - 2.0*v[k][i][j] + v[k][i-1][j])/(dx*dx);

				AA = -a1 + (a3+a4)/Re;
				BB = -b1 + (b3+b4)/Re;
				//////////////////////////////////////////////////////////////////////////
				u[k+1][i][j] = u[k][i][j] + dt*(AA - (p[k][i+1][j] - p[k][i][j])/dx);
				v[k+1][i][j] = v[k][i][j] + dt*(BB - (p[k][i][j+1] - p[k][i][j])/dy);
			}
		}	
		
		for(j=0; j<YHeight; j++)
		{
			for(i=0; i<XWidth; i++)
			{
				//////////////////////////////////////////////////////////////////////////
				fprintf(fp, "k=%d,i=%d,j=%d(%4.2f,%4.2f,%4.2f)\n", k,i,j, u[k][i][j], v[k][i][j], p[k][i][j]);
			}
		}	
		
		//////////////////////////////////////////////////////////////////////////
		/// Check Pressure Equation Here
		//////////////////////////////////////////////////////////////////////////
		memset(dp, 0, sizeof(float)*XWidth*YHeight);
		float tmp = 0.0;
		float sum = 0.0;
		do{
			for(j=1; j<YHeight-1; j++)
			{
				for(i=1; i<XWidth-1; i++)
				{
					tmp = dp[i][j];
					d = (u[k][i][j] - u[k][i-1][j])/dx + (v[k][i][j] - v[k][i][j-1])/dy;
					//////////////////////////////////////////////////////////////////////////
					dp[i][j] = -(b*(dp[i+1][j] + dp[i-1][j]) + c*(dp[i][j+1] + dp[i][j-1]) + d)/a;
					sum += (dp[i][j] - tmp)*(dp[i][j] - tmp);
				}
			}	
			fprintf(fp,"sum = %f\n",sum);
		}while(sum>0.001); 

		for(j=0; j<YHeight; j++)
		{
			for(i=0; i<XWidth; i++)
			{
				//////////////////////////////////////////////////////////////////////////
				fprintf(fp, "%4.2f  ", dp[i][j]);
			}
			fprintf(fp,"\n");
		}	
		
		//////////////////////////////////////////////////////////////////////////
		/// adjust the result of pressure
		maxerror = 0;
		for(j=1; j<YHeight-1; j++)
		{
			for(i=1; i<XWidth-1; i++)
			{
				//////////////////////////////////////////////////////////////////////////
				p[k+1][i][j] = p[k][i][j] + dp[i][j];
				if(maxerror < dp[i][j])
					maxerror = dp[i][j];
			}
		}	

		fprintf(fp, "\nIteration%d***********************************Max Error=%f\n", n, maxerror);

		
		//////////////////////////////////////////////////////////////////////////
		/// adjust velocity
		for(j=1; j<YHeight-1; j++)
		{
			for(i=1; i<XWidth-1; i++)
			{
				//////////////////////////////////////////////////////////////////////////
				u[k+1][i][j] = u[k][i][j] + (p[k][i+1][j] - p[k][i][j])/dx;
				v[k+1][i][j] = v[k][i][j] + (p[k][i][j+1] - p[k][i][j])/dy;
			}
		}	
		for(j=0; j<YHeight; j++)
		{
			for(i=0; i<XWidth; i++)
			{
				//////////////////////////////////////////////////////////////////////////
				fprintf(fp, "k=%d,i=%d,j=%d(%4.2f,%4.2f,%4.2f)\n", k,i,j, u[k][i][j], v[k][i][j], p[k][i][j]);
			}
		}	
		//////////////////////////////////////////////////////////////////////////
		/// check convergence
		if(maxerror > eps)
			goto loop;
	}
	
	//////////////////////////////////////////////////////////////////////////
	fprintf(fp, "-------------------Final Result----------------\n");
	for(k=0; k<TimeSteps; k++)
	{
		fprintf(fp, "Time=%d\n", k);
		for(j=0; j<YHeight; j++)
		{
			for(i=0; i<XWidth; i++)
			{
				//////////////////////////////////////////////////////////////////////////
				fprintf(fp, "%4.2f  ", u[k][i][j], v[k][i][j], p[k][i][j]);
			}
			fprintf(fp,"\n");
		}		
	}
	fclose(fp);

	return;
}