66int main(
int argc,
char *argv[])
72 int i, j, numModes, nFields;
75 if ((argc != 2) && (argc != 3))
78 "Usage: ./FldAddFalknerSkanBL sessionFile [SysSolnType]\n");
84 LibUtilities::SessionReader::CreateInstance(argc, argv);
95 string stability_solver;
98 BL_type = vSession->GetSolverInfo(
"BL_type");
99 txt_file = vSession->GetSolverInfo(
"txt_file");
100 stability_solver = vSession->GetSolverInfo(
"stability_solver");
102 if (stability_solver !=
"velocity_correction_scheme" &&
103 stability_solver !=
"coupled_scheme")
105 fprintf(stderr,
"Error: You must specify the stability solver in the "
106 "session file properly.\n");
107 fprintf(stderr,
"Options: 'velocity_correction_scheme' [output ===> "
108 "(u,v,p)]; 'coupled_scheme' [output ===>(u,v)]\n");
112 vSession->LoadParameter(
"Re", Re, 1.0);
113 vSession->LoadParameter(
"L",
L, 1.0);
114 vSession->LoadParameter(
"U_inf", U_inf, 1.0);
115 vSession->LoadParameter(
"x", x, 1.0);
116 vSession->LoadParameter(
"x_0", x_0, 1.0);
117 vSession->LoadParameter(
"NumberLines_txt", numLines, 1.0);
123 "Error: x must be positive ===> CHECK the session file\n");
129 fprintf(stderr,
"Error: x_0 must be positive or at least equal to 0 "
130 "===> CHECK the session file\n");
133 std::cout <<
"\n==========================================================="
135 std::cout <<
"Falkner-Skan Boundary Layer Generation (version of July 12th "
137 std::cout <<
"============================================================="
139 std::cout <<
"*************************************************************"
141 std::cout <<
"DATA FROM THE SESSION FILE:\n";
142 std::cout <<
"Reynolds number = " << Re <<
"\t[-]"
144 std::cout <<
"Characteristic length = " <<
L <<
"\t\t[m]"
146 std::cout <<
"U_infinity = " << U_inf
147 <<
"\t[m/s]" << std::endl;
148 std::cout <<
"Position x (parallel case only) = " << x <<
"\t\t[m]"
150 std::cout <<
"Position x_0 to start the BL [m] = " << x_0
151 <<
"\t\t[m]" << std::endl;
152 std::cout <<
"Number of lines of the .txt file = " << numLines
153 <<
"\t[-]" << std::endl;
154 std::cout <<
"BL type = " << BL_type
156 std::cout <<
".txt file read = " << txt_file
158 std::cout <<
"Stability solver = "
159 << stability_solver << std::endl;
160 std::cout <<
"*************************************************************"
162 std::cout <<
"-------------------------------------------------------------"
164 std::cout <<
"MESH and EXPANSION DATA:\n";
171 graphShPt = SpatialDomains::MeshGraphIO::Read(vSession);
176 vSession, graphShPt, vSession->GetVariable(0));
180 nElements = Domain->GetExpSize();
181 std::cout <<
"Number of elements = " << nElements
186 nQuadraturePts = Domain->GetTotPoints();
187 std::cout <<
"Number of quadrature points = " << nQuadraturePts
197 Domain->GetCoords(x_QuadraturePts, y_QuadraturePts, z_QuadraturePts);
201 const char *txt_file_char;
203 txt_file_char = txt_file.c_str();
206 ifstream pFile(txt_file_char);
209 cout <<
"Errro: Unable to open the .txt file with the BL data\n";
213 numLines = numLines / 3;
215 std::vector<std::vector<NekDouble>> GlobalArray(3);
217 for (j = 0; j <= 2; j++)
219 GlobalArray[j].resize(numLines);
220 for (i = 0; i <= numLines - 1; i++)
223 GlobalArray[j][i] =
d;
238 for (i = 0; i < numLines; i++)
240 eta[i] = GlobalArray[0][i];
241 f[i] = GlobalArray[1][i];
242 df[i] = GlobalArray[2][i];
269 if (BL_type ==
"Parallel")
271 for (i = 0; i < nQuadraturePts; i++)
273 eta_QuadraturePts[i] =
274 y_QuadraturePts[i] *
sqrt(U_inf / (2 * x * nu));
275 for (j = 0; j < numLines - 1; j++)
277 if (eta_QuadraturePts[i] >= eta[j] &&
278 eta_QuadraturePts[i] <= eta[j + 1])
280 f_QuadraturePts[i] = (eta_QuadraturePts[i] - eta[j]) *
282 (eta[j + 1] - eta[j]) +
284 df_QuadraturePts[i] = (eta_QuadraturePts[i] - eta[j]) *
285 (df[j + 1] - df[j]) /
286 (eta[j + 1] - eta[j]) +
290 else if (eta_QuadraturePts[i] == 1000000)
292 f_QuadraturePts[i] = f[numLines - 1];
293 df_QuadraturePts[i] = df[numLines - 1];
296 else if (eta_QuadraturePts[i] > eta[numLines - 1])
301 (eta_QuadraturePts[i] - eta[numLines - 1]);
302 df_QuadraturePts[i] = df[numLines - 1];
306 u_QuadraturePts[i] = U_inf * df_QuadraturePts[i];
308 nu *
sqrt(U_inf / (2.0 * nu * x)) *
309 (y_QuadraturePts[i] *
sqrt(U_inf / (2.0 * nu * x)) *
310 df_QuadraturePts[i] -
312 p_QuadraturePts[i] = 0.0;
319 if (BL_type ==
"nonParallel")
321 for (i = 0; i < nQuadraturePts; i++)
323 eta_QuadraturePts[i] =
325 sqrt(U_inf / (2 * (x_QuadraturePts[i] + x_0) * nu));
327 if ((x_QuadraturePts[i] + x_0) == 0)
329 eta_QuadraturePts[i] = 1000000;
332 for (j = 0; j < numLines - 1; j++)
334 if (eta_QuadraturePts[i] >= eta[j] &&
335 eta_QuadraturePts[i] <= eta[j + 1])
337 f_QuadraturePts[i] = (eta_QuadraturePts[i] - eta[j]) *
339 (eta[j + 1] - eta[j]) +
341 df_QuadraturePts[i] = (eta_QuadraturePts[i] - eta[j]) *
342 (df[j + 1] - df[j]) /
343 (eta[j + 1] - eta[j]) +
347 else if (eta_QuadraturePts[i] == 1000000)
349 f_QuadraturePts[i] = f[numLines - 1];
350 df_QuadraturePts[i] = df[numLines - 1];
353 else if (eta_QuadraturePts[i] > eta[numLines - 1])
358 (eta_QuadraturePts[i] - eta[numLines - 1]);
359 df_QuadraturePts[i] = df[numLines - 1];
363 u_QuadraturePts[i] = U_inf * df_QuadraturePts[i];
365 nu *
sqrt(U_inf / (2.0 * nu * (x_QuadraturePts[i] + x_0))) *
366 (y_QuadraturePts[i] *
367 sqrt(U_inf / (2.0 * nu * (x_QuadraturePts[i] + x_0))) *
368 df_QuadraturePts[i] -
372 if ((x_QuadraturePts[i] + x_0) == 0)
374 u_QuadraturePts[i] = U_inf;
375 v_QuadraturePts[i] = 0.0;
379 if ((x_QuadraturePts[i] + x_0) == 0 && y_QuadraturePts[i] == 0)
381 u_QuadraturePts[i] = 0.0;
382 v_QuadraturePts[i] = 0.0;
385 p_QuadraturePts[i] = 0.0;
393 etaOriginal = fopen(
"eta.txt",
"w+");
394 for (i = 0; i < numLines; i++)
396 fprintf(etaOriginal,
"%f %f %f \n", eta[i], f[i], df[i]);
401 yQ = fopen(
"yQ.txt",
"w+");
402 for (i = 0; i < nQuadraturePts; i++)
404 fprintf(yQ,
"%f %f %f %f %f %f %f\n", x_QuadraturePts[i],
405 y_QuadraturePts[i], u_QuadraturePts[i], v_QuadraturePts[i],
406 eta_QuadraturePts[i], f_QuadraturePts[i], df_QuadraturePts[i]);
415 vSession, graphShPt);
419 vSession, graphShPt);
423 vSession, graphShPt);
429 Basis = Domain->GetExp(0)->GetBasis(0);
432 std::cout <<
"Number of modes = " << numModes
436 Vmath::Vcopy(nQuadraturePts, u_QuadraturePts, 1, Exp2D_uk->UpdatePhys(), 1);
437 Vmath::Vcopy(nQuadraturePts, v_QuadraturePts, 1, Exp2D_vk->UpdatePhys(), 1);
438 Vmath::Vcopy(nQuadraturePts, p_QuadraturePts, 1, Exp2D_pk->UpdatePhys(), 1);
445 if (stability_solver ==
"velocity_correction_scheme")
451 Exp[0]->FwdTrans(Exp2D_uk->GetPhys(), Exp[0]->UpdateCoeffs());
452 Exp[1]->FwdTrans(Exp2D_vk->GetPhys(), Exp[1]->UpdateCoeffs());
454 if (stability_solver ==
"velocity_correction_scheme")
456 Exp[2]->FwdTrans(Exp2D_pk->GetPhys(), Exp[2]->UpdateCoeffs());
463 string FalknerSkan =
"FalknerSkan.fld";
466 if (stability_solver ==
"coupled_scheme")
471 if (stability_solver ==
"velocity_correction_scheme")
477 std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef =
478 Exp[0]->GetFieldDefinitions();
479 std::vector<std::vector<NekDouble>> FieldData(FieldDef.size());
481 for (j = 0; j < nFields; ++j)
483 for (i = 0; i < FieldDef.size(); ++i)
487 FieldDef[i]->m_fields.push_back(
"u");
491 FieldDef[i]->m_fields.push_back(
"v");
493 else if (j == 2 && stability_solver ==
"velocity_correction_scheme")
495 FieldDef[i]->m_fields.push_back(
"p");
497 Exp[j]->AppendFieldData(FieldDef[i], FieldData[i]);
503 std::cout <<
"-------------------------------------------------------------"
int main(int argc, char *argv[])
Main.
Represents a basis of a given type.
int GetNumModes() const
Return order of basis from the basis specification.
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
std::shared_ptr< Basis > BasisSharedPtr
std::shared_ptr< SessionReader > SessionReaderSharedPtr
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 ...
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
std::shared_ptr< ContField > ContFieldSharedPtr
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
std::vector< double > d(NPUPPER *NPUPPER)
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
scalarT< T > sqrt(scalarT< T > in)