diff --git a/.gitignore b/.gitignore index ddbc873..4843333 100644 --- a/.gitignore +++ b/.gitignore @@ -16,6 +16,3 @@ obj/ /.dropbox src/.vscode/ .vscode -grid/*.vtk -grid/*.dat -a.out diff --git a/Makefile b/Makefile index d91a2b1..bcdb846 100644 --- a/Makefile +++ b/Makefile @@ -2,7 +2,7 @@ CC := g++ #The Target Binary Program -TARGET := 2D_NS_Solver +TARGET := 2D_Euler_FFS #The Directories, Source, Includes, Objects, Binary and Resources SRCDIR := src diff --git a/grid/FlatPlate.cpp b/grid/FlatPlate.cpp deleted file mode 100644 index 7253d59..0000000 --- a/grid/FlatPlate.cpp +++ /dev/null @@ -1,287 +0,0 @@ -#include -#include -#include -#include -#include -#include - - -using namespace std; - -enum stretch_flag{str_of_line, end_of_line}; - -int Nx, Ny; -double **x, **y; -double dx, dy, xll, xul, yll, yul, len, D, A, num, den; - -template -T **Allocate_2D(T ** &m, int t1, int t2) { - m=new T* [t1]; - for (int i=0; i -T *linePartition(T *&t_line_to_part, const T t_start_of_line, const T t_end_of_line, const int t_part_size){ - - T part_len = (t_end_of_line - t_start_of_line)/t_part_size; - - for(int i=0;i<=t_part_size; i++){ - t_line_to_part[i] = t_start_of_line + i*part_len; - } - - return t_line_to_part; -} - -void WriteToFile(const int Nx, const int Ny, double **&x, double **&y){ - - //ofstream GridFFS("SWBLIGrid"+std::to_string(Nx)+std::to_string(Ny)+".dat"); - ofstream GridFFS("Blasius"+std::to_string(Nx)+std::to_string(Ny)+".dat"); - - ofstream GridFFSTecPlot("SWBLIGridTecPlot"+std::to_string(Nx)+std::to_string(Ny)+".dat"); - //GridFFS12040.flags(ios::dec | ios::scientific); - //GridFFS12040.precision(5); - - if(!GridFFS){ - cerr<<"File couldn't be opened to write the solution"< -T *strechedGrid(T *&vec, T vec_max, T vec_min, T vec_size){ - - T *eta, *A, B; - const float beta = 1.01; - - B = (beta + 1.0)/(beta - 1.0); - eta = new T [vec_size]; - A = new T [vec_size]; - - for(int j=0;j<=vec_size;j++){ - - A[j] = (beta + (1 - vec[j]/vec_max))/(beta - (1 - vec[j]/vec_max)); - std::cout<<"val of A="<0) and (fabs(delS)<=1*emax)) {mnb = 1.0*(-ric);} else mnb =0; @@ -343,4 +341,4 @@ T rounding(T dec, T n){ sp = round(pow(10,n)*dec); return sp/pow(10,n); -} +} \ No newline at end of file diff --git a/src/diss_nkfds_movers_2o.cpp b/src/diss_nkfds_movers_2o.cpp index 853fd29..609029c 100644 --- a/src/diss_nkfds_movers_2o.cpp +++ b/src/diss_nkfds_movers_2o.cpp @@ -323,9 +323,7 @@ T Mahalanobis(T rhoL, T rhoR, T unL, T unR, T pL, T pR, T maxa, T emax) { SR = cv*log(pR/(rhoR*(Gamma-1))) - RR*log(rhoR); delS = SL - SR; // if (fabs(delS/emax)<1) {delS = 0;} - //ric = (0.5*(fabs(unL)+fabs(unR))); // maxa; // - ric = maxa; // - + ric = (0.5*(fabs(unL)+fabs(unR))); // maxa; // /* if ric = maxa it is LLF. If ric = (0.5*(fabs(unL)+fabs(unR))) then it is RICCA type */ if ((DD >0) and (fabs(delS)<=1*emax)) {mnb = 1.0*(-ric);} else mnb =0; @@ -338,4 +336,4 @@ T rounding(T dec, T n){ sp = round(pow(10,n)*dec); return sp/pow(10,n); -} +} \ No newline at end of file diff --git a/src/diss_roe2_prim.cpp b/src/diss_roe2_prim.cpp index 9b8f66f..828c84f 100755 --- a/src/diss_roe2_prim.cpp +++ b/src/diss_roe2_prim.cpp @@ -52,11 +52,11 @@ void Diss_Roe2_Prim(int ib, int id1, int id2, int jb, int jd1, int jd2, double * for (int k = 0; k < 4; k++) { - // deltl[k] = 0.5 * dui[k][im1][j] * Minmod_Flux(dui[k][i][j], dui[k][im1][j]); - // deltr[k] = 0.5 * dui[k][i][j] * Minmod_Flux(dui[k][ip1][j], dui[k][i][j]); + deltl[k] = 0.5 * dui[k][im1][j] * Minmod_Flux(dui[k][i][j], dui[k][im1][j]); + deltr[k] = 0.5 * dui[k][i][j] * Minmod_Flux(dui[k][ip1][j], dui[k][i][j]); - deltl[k] = 0.5*dui[k][im1][j]*VanAlbada_Limiter(dui[k][i][j], dui[k][im1][j]); - deltr[k] = 0.5*dui[k][i][j]*VanAlbada_Limiter(dui[k][ip1][j], dui[k][i][j]); + /*deltl[k] = 0.5*dui[k][im1][j]*VanAlbada_Limiter(dui[k][i][j], dui[k][im1][j]); + deltr[k] = 0.5*dui[k][i][j]*VanAlbada_Limiter(dui[k][ip1][j], dui[k][i][j]);*/ /*deltl[k] = 0.5 * Venki_Limiter(dui[k][i][j], dui[k][im1][j], dx); deltr[k] = 0.5 * Venki_Limiter(dui[k][ip1][j], dui[k][i][j], dx);*/ @@ -167,11 +167,11 @@ void Diss_Roe2_Prim(int ib, int id1, int id2, int jb, int jd1, int jd2, double * for (int k = 0; k < 4; k++) { - // deltl[k] = 0.5 * duj[k][i][jm1] * Minmod_Flux(duj[k][i][j], duj[k][i][jm1]); - // deltr[k] = 0.5 * duj[k][i][j] * Minmod_Flux(duj[k][i][jp1], duj[k][i][j]); + deltl[k] = 0.5 * duj[k][i][jm1] * Minmod_Flux(duj[k][i][j], duj[k][i][jm1]); + deltr[k] = 0.5 * duj[k][i][j] * Minmod_Flux(duj[k][i][jp1], duj[k][i][j]); - deltl[k] = 0.5*duj[k][i][jm1]*VanAlbada_Limiter(duj[k][i][j], duj[k][i][jm1]); - deltr[k] = 0.5*duj[k][i][j]*VanAlbada_Limiter(duj[k][i][jp1], duj[k][i][j]); + /*deltl[k] = 0.5*duj[k][i][jm1]*VanAlbada_Limiter(duj[k][i][j], duj[k][i][jm1]); + deltr[k] = 0.5*duj[k][i][j]*VanAlbada_Limiter(duj[k][i][jp1], duj[k][i][j]);*/ /*deltl[k] = 0.5 * Venki_Limiter(duj[k][i][j], duj[k][i][jm1], dy); deltr[k] = 0.5 * Venki_Limiter(duj[k][i][jp1], duj[k][i][j], dy);*/ diff --git a/src/ini_flow_dimensional.cpp b/src/ini_flow_dimensional.cpp index 3224227..3859ff8 100755 --- a/src/ini_flow_dimensional.cpp +++ b/src/ini_flow_dimensional.cpp @@ -6,9 +6,6 @@ void InitFlowDimensional(int id2, int jd2, double Re_inf, double Machinf, double Lref, double alpha, double pinf, double ***&cv, double ***&dv) { - const double s1 = 110.0, s2 = 288.16; - double s12; - s12 = 1.0 + s1/s2; Cp = Rgas*GamaByGamaMin1; rhoinf = pinf/(Rgas*tinf); qinf = Machinf*sqrt((Gamma-1.0)*Cp*tinf); @@ -47,15 +44,12 @@ void InitFlowDimensional(int id2, int jd2, double Re_inf, double Machinf, double dv[3][i][j] = pinf; //p dv[4][i][j] = tinf; //T dv[5][i][j] = sqrt(Gamma*pinf/rhoinf); //a - // dv[6][i][j] = ((C0+C1)/(C1+dv[4][i][j]))*pow((dv[4][i][j]/C0),1.5)*refvisc; //Mu - dv[6][i][j] = (sqrt(dv[4][i][j]/s2)*s12/(1.0+s1/dv[4][i][j]))*refvisc; //Mu //both are same - + dv[6][i][j] = ((C0+C1)/(C1+dv[4][i][j]))*pow((dv[4][i][j]/C0),1.5)*refvisc; //Mu dv[7][i][j] = dv[6][i][j]*Cp/Pr; //K - /*std::cout <<"I.Cs: Dep Var"<<"\t"<>tot_time; infile.ignore(std::numeric_limits::max(), '\n'); - std::cout<<"Maximum time at which unsteady computations will be stopped="<49 && j<18) continue; + adtv = tstep[i][j] / area[i][j]; rhs[0][i][j] = adtv * rhs[0][i][j]; rhs[1][i][j] = adtv * rhs[1][i][j]; @@ -95,6 +103,9 @@ void Solver(int ib, int id1, int id2, int jb, int jd1, int jd2, int nconv, int n for (int j = 2; j <= jb; j++) { for (int i = 2; i <= ib; i++) { + + if(i>49 && j<18) continue; + cv[0][i][j] = ConstSSPRK2a[k] * cvold[0][i][j] + ConstSSPRK2b[k] * (cv[0][i][j] - rhs[0][i][j]); cv[1][i][j] = ConstSSPRK2a[k] * cvold[1][i][j] + ConstSSPRK2b[k] * (cv[1][i][j] - rhs[1][i][j]); cv[2][i][j] = ConstSSPRK2a[k] * cvold[2][i][j] + ConstSSPRK2b[k] * (cv[2][i][j] - rhs[2][i][j]); @@ -126,6 +137,9 @@ void Solver(int ib, int id1, int id2, int jb, int jd1, int jd2, int nconv, int n for (int j = 2; j <= jb; j++) { for (int i = 2; i <= ib; i++) { + + if(i>49 && j<18) continue; + adtv = tstep[i][j] / area[i][j]; rhs[0][i][j] = adtv * rhs[0][i][j]; rhs[1][i][j] = adtv * rhs[1][i][j]; @@ -138,6 +152,9 @@ void Solver(int ib, int id1, int id2, int jb, int jd1, int jd2, int nconv, int n for (int j = 2; j <= jb; j++) { for (int i = 2; i <= ib; i++) { + + if(i>49 && j<18) continue; + cv[0][i][j] = ConstRK1[k] * cvold[0][i][j] + ConstRK2[k] * (cv[0][i][j] - rhs[0][i][j]); cv[1][i][j] = ConstRK1[k] * cvold[1][i][j] + ConstRK2[k] * (cv[1][i][j] - rhs[1][i][j]); cv[2][i][j] = ConstRK1[k] * cvold[2][i][j] + ConstRK2[k] * (cv[2][i][j] - rhs[2][i][j]);