Untitled
raw download clone
TEXT
views 24
,
size 5235 b
#include <iostream>
#include <fstream>
#include <sstream>
#include <string>
#include <cmath>
#include <iomanip>
using namespace std;
/////////////////////////////////////////////


//
int const Ni = 100; // So diem nut
int const NPHI = 2;			 // So ptr can giai
//

double PHI[NPHI][Ni], PHIN[NPHI][Ni], G[NPHI][Ni], EX[NPHI][Ni];
double ROU[Ni], U[Ni], P[Ni], T[Ni], ON[Ni], X[Ni], RN1[Ni], UN1[Ni], PN1[Ni], TN1[Ni], RT1[Ni], UT1[Ni], PT1[Ni], TT1[Ni];

void InputData();
void metric();
void grid();
void Maccormack();
void outPut(int t);
void contcal();
void phical();
void phincal();
void invisx();
void varical();
void boundio();
void phasedt();
void Output(int t);
int NA, NAF, ITOUT, ITCLD, ITPRI;
double DT, P0, T0, U0;
double DH = 1, H1 = 0.5;
double CR = 287.22, CP = 1007;

int NOUT = 0;
int I0   = Ni - 1;
int I1   = I0 - 1;
int I2   = I1 - 1;

int main()
{
	float time_run = 0;
    float time_limit = 0.01;
	InputData();
	grid();
	double NTIMST = NAF - NA;
	double H    = 1.0;
	for(int i = 1; i <= Ni;i++)
        T[i] = T0;

    cout << "Start at: " << time_run << endl;
    int i=1;
    while(time_run < time_limit)
    {
        Maccormack();
        time_run += DT;
        cout << i << "\t" << DT << "\t" << fixed << setprecision(5) << time_run << endl;
        ++i;
    }
    Output(time_run);
    cout << "End at: " << time_run << endl;
	return 0;
}

void InputData()
{
	ifstream PARAMETER(".//DATA//INPUT.dat");
	string temp;
	PARAMETER >> temp >> NA;
	PARAMETER >> temp >> NAF;
	PARAMETER >> temp >> ITOUT;
	PARAMETER >> temp >> ITCLD;
	PARAMETER >> temp >> ITPRI;
	PARAMETER >> temp >> DT;
	PARAMETER >> temp >> P0;
	PARAMETER >> temp >> T0;
	PARAMETER >> temp >> U0;
	PARAMETER.close();
	
	for (int i = 1; i <= Ni; i++)
	{
		if (X[i] <= 10)
		{
			ROU[i] = 1;
			P[i] = 100000;
			U[i] = 0;
		}
		else
		{
			ROU[i] = 0.1;
			P[i] = 10000;
			U[i] = 0;
		}
	}
}

void metric() // Tinh cac nut luoi
{
	for (int i = 1; i <= Ni; i++)
	{
		DH = 1.0/(20.0/(Ni - 1));
		H1 = 0.5*DH;
		X[i] = 20.0/(Ni - 1)*(i - 1);
	}
}

void grid()
{
	double d_xi = (double)1/(Ni - 1);
	for (int i = 1; i <= Ni; i++)
	{
		double xi = (double)d_xi*(i - 1);
		X[i] = 20*xi;
	}
}

/*void residual()
{
	Umax = 0;
	for (int i = 1; i <= Ni; i++)
	{
		HMIN = DH;
		Umax = max(Umax, ON[i] + abs(U[i]));
	}
	
	CFL = Umax*DT/HMIN; //Courant number
	DRMAX = 0;
	DUMAX = 0;
    DPMAX = 0;
    for (int i = 1; i <= Ni; i++)
    {
		DRMAX = max(DRMAX, abs(RN1[i] - ROU[i]));
    	DUMAX = max(DUMAX, abs(UN1[i] -  U[i]));
    	DPMAX = max(DPMAX, abs(PN1[i] -  P[i]));	
	}
	
	
}*/

void contcal()
{
	for (int i = 1; i <= Ni; i++)
	{
		ON[i] = sqrt(1.4*CR*T0);
	}
}

void phical()
{
	for (int i = 1; i <= Ni; i++)
	{
		PHI[1][i] = ROU[i];
		PHI[2][i] = PHI[1][i]*U[i];
	}

	for (int i =2; i <= I0; i++)
	{
		for (int k = 1; k <= NPHI; k++)
		{
			G[k][i] = 1.0*fabs(P[i + 1] - 2.0*P[i] + P[i - 1]) / (P[i + 1] + 2.0*P[i] + P[i - 1]) * (PHI[k][i + 1] - 2.0*PHI[k][i] + PHI[k][i - 1]);
		}
	}
	for ( int k = 1; k <= NPHI; k++)
	{
		G[k][1] = G[k][2];
		G[k][Ni] = G[k][I0];
	}
}

void phincal()
{
	for (int i = 1; i <= Ni; i++)
	{
		PHIN[1][i] = ROU[i];
		PHIN[2][i] = PHIN[1][i]*U[i];
	}
	
	for (int i = 2; i <= I0; i++)
	{
		for (int k = 1; k <= NPHI; k++)
		{
			G[k][i] = 1.0*fabs(P[i + 1] - 2.0*P[i] + P[i - 1]) / (P[i + 1] + 2.0*P[i] + P[i - 1]) * (PHIN[k][i + 1] - 2.0*PHIN[k][i] + PHIN[k][i - 1]);
		}
	}
	
	for (int k = 1; k <= NPHI; k++)
	{
		G[k][1] = G[k][2];
		G[k][Ni] = G[k][I0];
	}
}

void Maccormack()
{
	//////////////// Predictor Step ///////////////
	contcal();
	phical();
	invisx();
	
	for (int i = 2; i <= I0; i++)
	{
		for (int k = 1; k <= NPHI; k++)
		{
			PHIN[k][i] = PHI[k][i] - DT*(EX[k][i] - EX[k][i - 1])*DH + G[k][i];
		}
	}
	for (int k = 1; k <= NPHI; k++)
	{
		PHIN[k][Ni] = PHIN[k][I0];
		PHIN[k][1] = PHIN[k][2];
	}
	
	varical();
	boundio();
	
	/////////////// Corrector Step ///////////////
	contcal();
	phincal();
	invisx();
	
	for (int i = 2; i <= I0; i++)
	{
		for (int k = 1; k <= NPHI; k++)
		{
			PHIN[k][i] = 0.5*(PHI[k][i] + PHIN[k][i]) - 0.5*DT*(EX[k][i + 1] - EX[k][i])*DH + G[k][i];
		}
	}
	for (int k = 1; k <= NPHI; k++)
	{
		PHIN[k][Ni] = PHIN[k][I0];
		PHIN[k][1] = PHIN[k][2];
	}
	
	varical();
	boundio();
}

void invisx()
{
	for (int i = 1; i <= Ni; i++)
	{
		EX[1][i] = ROU[i]*U[i];
		EX[2][i] = ROU[i]*U[i]*U[i] + P[i];
	}
}

void varical()
{
	for (int i = 1; i <= Ni; i++)
	{
		ROU[i] = PHIN[1][i];
		U[i] = PHIN[2][i]/PHIN[1][i];
	}
	
	for (int i = 1; i <= Ni; i++)
	{
		P[i] = CR*T[i]*ROU[i];
	}
}

void boundio()
{
	ROU[1] = ROU[2];
	U[1] = U[2];
	P[1] = P[2];
 
	ROU[Ni] = ROU[I0];
	U[Ni] = U[I0];
	P[Ni] = P[I0];
}

void Output(int t)
{
	stringstream nameFile;
	nameFile << NOUT;
	string name;
	ofstream outData;
	name = ".//DATA//" + nameFile.str() + ".dat";
	outData.open(name.c_str());
	for (int i = 1; i <= Ni; i++)
	{
		outData << ROU[i] << "\t" << U[i] << "\t" << P[i] << endl;
	}
	outData.close();
}
close fullscreen
Login or Register to edit or fork this paste. It's free.