113 dcdg = 1.0 / (2. * sqrt(v[3])) - sqrt(v[3]) / (v[3]+
m_Suth);
119 c = pow(v[3], (
Omega-1.0));
120 dcdg = (
Omega - 1.0) * pow(v[3], (
Omega - 2.0));
131 dv[2] = - v[2] * (cp + v[0]) / c;
133 dv[4] = - v[4] * (cp +
m_Pr * v[0]) / c -
156 for (
int i = 0; i < n ; i++)
158 yt[i] = y[i] + hh * dydx[i];
163 for (
int i = 0; i < n; i++)
165 yt[i] = y[i] + hh * dyt[i];
170 for (
int i = 0; i < n; i++)
172 yt[i] = y[i] + h * dym[i];
173 dym[i] = dyt[i] + dym[i];
178 for (
int i = 0; i < n; i++)
180 yout[i] = y[i] + h6 * (dydx[i] + dyt[i] + 2 * dym[i]);
201 for (
int i = 0; i < nvar; i++)
209 h = (x2-x1) / m_xpoints;
214 RK4 (v, dv, nvar, x, h, v);
218 cout <<
"bug" << endl;
224 for (
int i = 0; i < nvar; i++)
266 z[i] = z[i-1] + 0.5 * (xx[i] - xx[i-1]) * (ff[3][i] + ff[3][i-1]);
267 dm = ff[3][i-1] - ff[1][i-1];
268 dd = ff[3][i] - ff[1][i];
269 sumd = sumd + 0.5 * (xx[i] - xx[i-1]) * (dd + dm);
271 if ((ff[1][i] > 0.999) && (flg < 1.0))
280 file3.open(
"physical_data.dat");
286 for (
int k = 0; k < 5; k++)
293 rho[i] = (1.0 / ff[3][i]);
294 vv[i] = -ff[0][i]/sqrt(
m_uInf);
296 velocity[i] = ff[0][i] ;
301 for (
int i = 0; i < nQuadraturePts; i++)
305 cout <<
"i" <<
" " << i <<
"/" << nQuadraturePts << endl;
308 xcher = x_QuadraturePts[i];
309 ycher = y_QuadraturePts[i];
313 rex = 0.5 * pow(((
m_Re) / scale), 2) + (
m_Re) * xin;
314 delsx = sqrt(2.0 / rex) * scale * (xin)*
m_Pr;
315 scale = scale / delsx;
317 scale2 = ycher * (scale * delta) / sqrt(
etamax) ;
318 coeff = 0.5 * sqrt( 2 / (xcher*
m_Re)) ;
320 if (scale2 > z[m_xpoints-3])
322 u_QuadraturePts[i] = 1;
323 rho_QuadraturePts[i] = 1;
324 T_QuadraturePts[i] = 1.0 / rho_QuadraturePts[i];
325 v_QuadraturePts[i] = coeff * (z[m_xpoints-3] -
326 velocity[m_xpoints-3]);
328 file3 << xcher <<
" "
330 << velocity[m_xpoints-3] <<
" "
331 << z[m_xpoints-3] <<
" "
337 for (
int j = 0 ; j< m_xpoints-1; j++)
339 if ((z[j] <= scale2) && (z[j+1] > scale2))
347 ASSERTL0(
false,
"Could not determine index in CompressibleBL");
350 u_QuadraturePts[i] = u[index];
351 rho_QuadraturePts[i] = rho[index];
352 T_QuadraturePts[i] = 1.0/rho_QuadraturePts[i];
353 v_QuadraturePts[i] = coeff * (u[index]*scale2 - velocity[index]);
363 int main(
int argc,
char *argv[])
371 int i, j, k, numModes;
377 for (i = 0; i < nmax; i++)
393 = LibUtilities::SessionReader::CreateInstance(argc, argv);
397 = SpatialDomains::MeshGraph::Read(vSession);
399 int expdim = graphShPt->GetMeshDimension();
401 int nElements, nQuadraturePts = 0;
414 vSession->GetVariable(0));
417 nElements = Domain->GetExpSize();
418 std::cout <<
"Number of elements = "
419 << nElements << std::endl;
422 nQuadraturePts = Domain->GetTotPoints();
423 std::cout <<
"Number of quadrature points = "
424 << nQuadraturePts << std::endl;
430 Domain->GetCoords(x_QuadraturePts, y_QuadraturePts, z_QuadraturePts);
432 else if (expdim == 3)
442 nElements = Domain->GetExpSize();
443 std::cout <<
"Number of elements = "
444 << nElements << std::endl;
447 nQuadraturePts = Domain->GetTotPoints();
448 std::cout <<
"Number of quadrature points = "
449 << nQuadraturePts << std::endl;
455 Domain->GetCoords(x_QuadraturePts, y_QuadraturePts, z_QuadraturePts);
459 ASSERTL0(
false,
"Routine available for 2D and 3D problem only.")
463 vSession->LoadParameter(
"Re",
m_Re, 1.0);
464 vSession->LoadParameter(
"Mach",
m_Mach, 1.0);
465 vSession->LoadParameter(
"TInf",
m_Tinf, 1.0);
466 vSession->LoadParameter(
"Twall",
m_Twall, 1.0);
467 vSession->LoadParameter(
"Gamma",
m_Gamma, 1.0);
468 vSession->LoadParameter(
"Pr",
m_Pr, 1.0);
469 vSession->LoadParameter(
"L",
m_long, 1.0);
470 vSession->LoadParameter(
"rhoInf",
m_rhoInf, 1.0);
471 vSession->LoadParameter(
"uInf",
m_uInf, 1.0);
472 vSession->LoadParameter(
"GasConstant",
m_R, 1.0);
473 vSession->LoadParameter(
"vInf",
m_vInf, 1.0);
474 vSession->LoadParameter(
"mu",
m_mu, 1.0);
481 cout <<
"Number of points" <<
" " <<
m_xpoints << endl;
497 v[0] = 0.47 * pow(v[1], 0.21);
501 v[1] = 0.062 * pow(
m_Mach, 2) - 0.1 * (
m_Tw - 1.0) *
503 v[0] = 0.45 - 0.01 *
m_Mach + (
m_Tw - 1.0) * 0.06;
532 for (k = 0; k < maxit; k++)
551 cout <<
"err" << scientific << setprecision(9) <<
" " << err << endl;
557 cout <<
"Calculating" << endl;
559 y_QuadraturePts, u_QuadraturePts, v_QuadraturePts,
560 rho_QuadraturePts, T_QuadraturePts);
567 vstart[2] = v[0] + dv[0];
588 vstart[3] = v[1] + dv[1];
593 vstart[4] = v[1] + dv[1];
601 al11 = (f1[0] - f[0]) / dv[0];
602 al21 = (f1[1] - f[1]) / dv[0];
603 al12 = (f2[0] - f[0]) / dv[1];
604 al22 = (f2[1] - f[1]) / dv[1];
605 det = al11 * al22 - al21 * al12;
607 dv[0] = ( - al22 * f[0] + al12 * f[1]) / det;
608 dv[1] = (al21 * f[0] - al11 * f[1]) / det;
613 else if (expdim == 3)
618 cout <<
"Calculating" << endl;
620 z_QuadraturePts, u_QuadraturePts, v_QuadraturePts,
621 rho_QuadraturePts, T_QuadraturePts);
628 vstart[2] = v[0] + dv[0];
649 vstart[3] = v[1] + dv[1];
654 vstart[4] = v[1] + dv[1];
662 al11 = (f1[0] - f[0]) / dv[0];
663 al21 = (f1[1] - f[1]) / dv[0];
664 al12 = (f2[0] - f[0]) / dv[1];
665 al22 = (f2[1] - f[1]) / dv[1];
666 det = al11 * al22 - al21 * al12;
668 dv[0] = ( - al22 * f[0] + al12 * f[1]) / det;
669 dv[1] = (al21 * f[0] - al11 * f[1]) / det;
679 verif.open(
"similarity_solution.dat");
680 for (i=0; i< nQuadraturePts; i++)
682 verif << scientific << setprecision(9) << x_QuadraturePts[i]
683 <<
" \t " << y_QuadraturePts[i] <<
" \t " ;
684 verif << scientific << setprecision(9) << u_QuadraturePts[i]
685 <<
" \t " << v_QuadraturePts[i] <<
" \t " ;
686 verif << scientific << setprecision(9) << rho_QuadraturePts[i]
687 <<
" \t " << T_QuadraturePts[i] << endl;
692 for (i = 0; i < nQuadraturePts; i++)
694 rho_QuadraturePts[i] = rho_QuadraturePts[i] *
m_rhoInf;
695 u_QuadraturePts[i] = u_QuadraturePts[i] *
m_uInf;
696 v_QuadraturePts[i] = v_QuadraturePts[i] *
m_uInf;
697 T_QuadraturePts[i] = T_QuadraturePts[i] *
m_Tinf;
699 T_QuadraturePts[i] = T_QuadraturePts[i] * rho_QuadraturePts[i] *
m_R;
700 T_QuadraturePts[i] = T_QuadraturePts[i] / (
m_Gamma-1);
701 T_QuadraturePts[i] = T_QuadraturePts[i] + 0.5 * rho_QuadraturePts[i] * (
702 pow(u_QuadraturePts[i], 2.0) + pow(v_QuadraturePts[i], 2.0));
704 u_QuadraturePts[i] = u_QuadraturePts[i] * rho_QuadraturePts[i];
705 v_QuadraturePts[i] = v_QuadraturePts[i] * rho_QuadraturePts[i];
737 Basis = Domain->GetExp(0)->GetBasis(0);
738 numModes = Basis->GetNumModes();
740 std::cout <<
"Number of modes = " << numModes << std::endl;
745 Exp2D_uk->UpdatePhys(), 1);
747 Exp2D_vk->UpdatePhys(), 1);
749 Exp2D_rhok->UpdatePhys(), 1);
751 Exp2D_Tk->UpdatePhys(), 1);
760 Exp[0]->FwdTrans(Exp2D_rhok->GetPhys(), Exp[0]->UpdateCoeffs());
761 Exp[1]->FwdTrans(Exp2D_uk->GetPhys(), Exp[1]->UpdateCoeffs());
762 Exp[2]->FwdTrans(Exp2D_vk->GetPhys(), Exp[2]->UpdateCoeffs());
763 Exp[3]->FwdTrans(Exp2D_Tk->GetPhys(), Exp[3]->UpdateCoeffs());
766 cout << argv[1] << endl;
767 string tmp = argv[1];
768 int len = tmp.size();
769 for (i = 0; i < len-4; ++i)
771 file_name += argv[1][i];
773 file_name = file_name+
".rst";
776 std::vector<LibUtilities::FieldDefinitionsSharedPtr>
777 FieldDef = Exp[0]->GetFieldDefinitions();
778 std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
780 for (j = 0; j < 4; j++)
782 for (i = 0; i < FieldDef.size(); i++)
786 FieldDef[i]->m_fields.push_back(
"rho");
790 FieldDef[i]->m_fields.push_back(
"rhou");
794 FieldDef[i]->m_fields.push_back(
"rhov");
798 FieldDef[i]->m_fields.push_back(
"E");
800 Exp[j]->AppendFieldData(FieldDef[i], FieldData[i]);
806 else if (expdim == 3)
842 Basis = Domain->GetExp(0)->GetBasis(0);
843 numModes = Basis->GetNumModes();
845 std::cout<<
"Number of modes = " << numModes << std::endl;
850 Exp3D_rhok->UpdatePhys(), 1);
852 Exp3D_uk->UpdatePhys(), 1);
854 Exp3D_vk->UpdatePhys(), 1);
856 Exp3D_wk->UpdatePhys(), 1);
858 Exp3D_Tk->UpdatePhys(), 1);
868 Exp[0]->FwdTrans(Exp3D_rhok->GetPhys(), Exp[0]->UpdateCoeffs());
869 Exp[1]->FwdTrans(Exp3D_uk->GetPhys(), Exp[1]->UpdateCoeffs());
870 Exp[2]->FwdTrans(Exp3D_vk->GetPhys(), Exp[2]->UpdateCoeffs());
871 Exp[3]->FwdTrans(Exp3D_wk->GetPhys(), Exp[3]->UpdateCoeffs());
872 Exp[4]->FwdTrans(Exp3D_Tk->GetPhys(), Exp[4]->UpdateCoeffs());
875 cout << argv[1] << endl;
876 string tmp = argv[1];
877 int len = tmp.size();
878 for (i = 0; i < len-4; ++i)
880 file_name += argv[1][i];
882 file_name = file_name+
".rst";
885 std::vector<LibUtilities::FieldDefinitionsSharedPtr>
886 FieldDef = Exp[0]->GetFieldDefinitions();
887 std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
889 for (j = 0; j < 5; j++)
891 for (i = 0; i < FieldDef.size(); i++)
895 FieldDef[i]->m_fields.push_back(
"rho");
899 FieldDef[i]->m_fields.push_back(
"rhou");
903 FieldDef[i]->m_fields.push_back(
"rhov");
907 FieldDef[i]->m_fields.push_back(
"rhow");
911 FieldDef[i]->m_fields.push_back(
"E");
913 Exp[j]->AppendFieldData(FieldDef[i], FieldData[i]);
920 std::cout <<
"----------------------------------------------------\n";
921 std::cout <<
"\n=================================================\n";
922 std::cout <<
"Similarity solution \n";
923 std::cout <<
"===================================================\n";
924 std::cout <<
"***************************************************\n";
925 std::cout <<
"DATA FROM THE SESSION FILE:\n";
926 std::cout <<
"Reynolds number = " <<
m_Re
927 <<
"\t[-]" << std::endl;
928 std::cout <<
"Mach number = " <<
m_Mach
929 <<
"\t[-]" << std::endl;
930 std::cout <<
"Characteristic length = " <<
m_long
931 <<
"\t[m]" << std::endl;
932 std::cout <<
"U_infinity = " <<
m_uInf
933 <<
"\t[m/s]" << std::endl;
934 std::cout <<
"***************************************************\n";
935 std::cout <<
"---------------------------------------------------\n";
936 std::cout <<
"MESH and EXPANSION DATA:\n";
937 std::cout <<
"Done." << std::endl;
#define ASSERTL0(condition, msg)
void OUTPUT(int m_xpoints, Array< OneD, NekDouble > xx, Array< OneD, Array< OneD, NekDouble > > ff, int nQuadraturePts, Array< OneD, NekDouble > x_QuadraturePts, Array< OneD, NekDouble > y_QuadraturePts, Array< OneD, NekDouble > u_QuadraturePts, Array< OneD, NekDouble > v_QuadraturePts, Array< OneD, NekDouble > rho_QuadraturePts, Array< OneD, NekDouble > T_QuadraturePts)
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
boost::shared_ptr< ContField2D > ContField2DSharedPtr
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
void RKDUMB(Array< OneD, NekDouble > vstart, int nvar, NekDouble x1, NekDouble x2, int m_xpoints, Array< OneD, NekDouble > xx, Array< OneD, Array< OneD, NekDouble > > y)
void RK4(Array< OneD, NekDouble > y, Array< OneD, NekDouble > dydx, int n, NekDouble x, NekDouble h, Array< OneD, NekDouble > yout)
boost::shared_ptr< ExpList2D > ExpList2DSharedPtr
Shared pointer to an ExpList2D object.
void Write(const std::string &outFile, std::vector< FieldDefinitionsSharedPtr > &fielddefs, std::vector< std::vector< NekDouble > > &fielddata, const FieldMetaDataMap &fieldinfomap, const bool backup)
This function allows for data to be written to an FLD file when a session and/or communicator is not ...
boost::shared_ptr< ExpList3D > ExpList3DSharedPtr
Shared pointer to an ExpList3D object.
int main(int argc, char *argv[])
boost::shared_ptr< Basis > BasisSharedPtr
boost::shared_ptr< ContField3D > ContField3DSharedPtr
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
void COMPBL(Array< OneD, NekDouble > v, Array< OneD, NekDouble > dv)