Untitled
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();
}``````