33 int main(
int argc,
char *argv[])
39 int i, j, numModes, nFields;
42 if((argc != 2) && (argc != 3))
44 fprintf(stderr,
"Usage: ./FldAddFalknerSkanBL sessionFile [SysSolnType]\n");exit(1);
59 string stability_solver;
63 BL_type = vSession->GetSolverInfo(
"BL_type");
64 txt_file = vSession->GetSolverInfo(
"txt_file");
65 stability_solver = vSession->GetSolverInfo(
"stability_solver");
67 if(stability_solver !=
"velocity_correction_scheme" &&
68 stability_solver !=
"coupled_scheme")
70 fprintf(stderr,
"Error: You must specify the stability solver in the session file properly.\n");
71 fprintf(stderr,
"Options: 'velocity_correction_scheme' [output ===> (u,v,p)]; 'coupled_scheme' [output ===>(u,v)]\n");
75 vSession->LoadParameter(
"Re", Re, 1.0);
76 vSession->LoadParameter(
"L", L, 1.0);
77 vSession->LoadParameter(
"U_inf", U_inf, 1.0);
78 vSession->LoadParameter(
"x", x, 1.0);
79 vSession->LoadParameter(
"x_0", x_0, 1.0);
80 vSession->LoadParameter(
"NumberLines_txt", numLines, 1.0);
85 fprintf(stderr,
"Error: x must be positive ===> CHECK the session file\n");
91 fprintf(stderr,
"Error: x_0 must be positive or at least equal to 0 ===> CHECK the session file\n");
94 std::cout<<
"\n=========================================================================\n";
95 std::cout<<
"Falkner-Skan Boundary Layer Generation (version of July 12th 2012)\n";
96 std::cout<<
"=========================================================================\n";
97 std::cout<<
"*************************************************************************\n";
98 std::cout<<
"DATA FROM THE SESSION FILE:\n";
99 std::cout <<
"Reynolds number = " << Re <<
"\t[-]" << std::endl;
100 std::cout <<
"Characteristic length = " << L <<
"\t\t[m]" << std::endl;
101 std::cout <<
"U_infinity = " << U_inf <<
"\t[m/s]" << std::endl;
102 std::cout <<
"Position x (parallel case only) = " << x <<
"\t\t[m]" << std::endl;
103 std::cout <<
"Position x_0 to start the BL [m] = " << x_0 <<
"\t\t[m]" << std::endl;
104 std::cout <<
"Number of lines of the .txt file = " << numLines <<
"\t[-]" << std::endl;
105 std::cout <<
"BL type = " << BL_type << std::endl;
106 std::cout <<
".txt file read = " << txt_file << std::endl;
107 std::cout <<
"Stability solver = " << stability_solver << std::endl;
108 std::cout<<
"*************************************************************************\n";
109 std::cout<<
"-------------------------------------------------------------------------\n";
110 std::cout<<
"MESH and EXPANSION DATA:\n";
125 nElements = Domain->GetExpSize();
126 std::cout <<
"Number of elements = " << nElements << std::endl;
130 nQuadraturePts = Domain->GetTotPoints();
131 std::cout <<
"Number of quadrature points = " << nQuadraturePts << std::endl;
140 Domain->GetCoords(x_QuadraturePts,y_QuadraturePts,z_QuadraturePts);
143 const char *txt_file_char;
145 txt_file_char = txt_file.c_str();
148 ifstream pFile(txt_file_char);
151 cout <<
"Errro: Unable to open the .txt file with the BL data\n";
155 numLines = numLines/3;
157 std::vector<std::vector<NekDouble> > GlobalArray (3);
161 GlobalArray[j].resize(numLines);
162 for (i=0; i<=numLines-1; i++)
165 GlobalArray[j][i] = d;
180 for (i=0; i<numLines; i++)
182 eta[i] = GlobalArray[0][i];
183 f[i] = GlobalArray[1][i];
184 df[i] = GlobalArray[2][i];
212 if(BL_type ==
"Parallel")
214 for(i=0; i<nQuadraturePts; i++)
216 eta_QuadraturePts[i] = y_QuadraturePts[i] * sqrt(U_inf / (2 * x * nu));
217 for(j=0; j<numLines-1; j++)
219 if(eta_QuadraturePts[i] >= eta[j] && eta_QuadraturePts[i] <= eta[j+1])
221 f_QuadraturePts[i] = (eta_QuadraturePts[i] - eta[j]) * (f[j+1] - f[j]) / (eta[j+1] - eta[j]) + f[j];
222 df_QuadraturePts[i] = (eta_QuadraturePts[i] - eta[j]) * (df[j+1] - df[j]) / (eta[j+1] - eta[j]) + df[j];
225 else if(eta_QuadraturePts[i] == 1000000)
227 f_QuadraturePts[i] = f[numLines-1];
228 df_QuadraturePts[i] = df[numLines-1];
231 else if(eta_QuadraturePts[i] > eta[numLines-1])
233 f_QuadraturePts[i] = f[numLines-1] + df[numLines-1] * (eta_QuadraturePts[i] - eta[numLines-1]);
234 df_QuadraturePts[i] = df[numLines-1];
239 u_QuadraturePts[i] = U_inf * df_QuadraturePts[i];
240 v_QuadraturePts[i] = nu * sqrt(U_inf / (2.0 * nu * x)) * (y_QuadraturePts[i] * sqrt(U_inf / (2.0 * nu * x)) * df_QuadraturePts[i] - f_QuadraturePts[i]);
241 p_QuadraturePts[i] = 0.0;
249 if(BL_type ==
"nonParallel")
251 for(i=0; i<nQuadraturePts; i++)
253 eta_QuadraturePts[i] = y_QuadraturePts[i] * sqrt(U_inf / (2 * (x_QuadraturePts[i] + x_0) * nu));
255 if((x_QuadraturePts[i] + x_0) == 0)
257 eta_QuadraturePts[i] = 1000000;
260 for(j=0; j<numLines-1; j++)
262 if(eta_QuadraturePts[i] >= eta[j] && eta_QuadraturePts[i] <= eta[j+1])
264 f_QuadraturePts[i] = (eta_QuadraturePts[i] - eta[j]) * (f[j+1] - f[j]) / (eta[j+1] - eta[j]) + f[j];
265 df_QuadraturePts[i] = (eta_QuadraturePts[i] - eta[j]) * (df[j+1] - df[j]) / (eta[j+1] - eta[j]) + df[j];
268 else if(eta_QuadraturePts[i] == 1000000)
270 f_QuadraturePts[i] = f[numLines-1];
271 df_QuadraturePts[i] = df[numLines-1];
274 else if(eta_QuadraturePts[i] > eta[numLines-1])
276 f_QuadraturePts[i] = f[numLines-1] + df[numLines-1] * (eta_QuadraturePts[i] - eta[numLines-1]);
277 df_QuadraturePts[i] = df[numLines-1];
281 u_QuadraturePts[i] = U_inf * df_QuadraturePts[i];
282 v_QuadraturePts[i] = nu * sqrt(U_inf / (2.0 * nu * (x_QuadraturePts[i] + x_0))) *
283 (y_QuadraturePts[i] * sqrt(U_inf / (2.0 * nu * (x_QuadraturePts[i] + x_0))) *
284 df_QuadraturePts[i] - f_QuadraturePts[i]);
287 if((x_QuadraturePts[i] + x_0) == 0)
289 u_QuadraturePts[i] = U_inf;
290 v_QuadraturePts[i] = 0.0;
294 if((x_QuadraturePts[i] + x_0) == 0 && y_QuadraturePts[i] == 0)
296 u_QuadraturePts[i] = 0.0;
297 v_QuadraturePts[i] = 0.0;
300 p_QuadraturePts[i] = 0.0;
309 etaOriginal = fopen(
"eta.txt",
"w+");
310 for(i=0; i<numLines; i++)
312 fprintf(etaOriginal,
"%f %f %f \n", eta[i], f[i], df[i]);
318 yQ = fopen(
"yQ.txt",
"w+");
319 for(i=0; i<nQuadraturePts; i++)
321 fprintf(yQ,
"%f %f %f %f %f %f %f\n", x_QuadraturePts[i], y_QuadraturePts[i], u_QuadraturePts[i],
322 v_QuadraturePts[i], eta_QuadraturePts[i], f_QuadraturePts[i], df_QuadraturePts[i]);
344 Basis = Domain->GetExp(0)->GetBasis(0);
345 numModes = Basis->GetNumModes();
347 std::cout<<
"Number of modes = " << numModes << std::endl;
350 Vmath::Vcopy(nQuadraturePts, u_QuadraturePts , 1, Exp2D_uk->UpdatePhys(), 1);
351 Vmath::Vcopy(nQuadraturePts, v_QuadraturePts, 1, Exp2D_vk->UpdatePhys(), 1);
352 Vmath::Vcopy(nQuadraturePts, p_QuadraturePts , 1, Exp2D_pk->UpdatePhys(), 1);
359 if(stability_solver ==
"velocity_correction_scheme")
363 Exp[0]->FwdTrans(Exp2D_uk->GetPhys(),Exp[0]->UpdateCoeffs());
364 Exp[1]->FwdTrans(Exp2D_vk->GetPhys(),Exp[1]->UpdateCoeffs());
366 if(stability_solver ==
"velocity_correction_scheme")
367 Exp[2]->FwdTrans(Exp2D_pk->GetPhys(),Exp[2]->UpdateCoeffs());
374 string FalknerSkan =
"FalknerSkan.fld";
377 if(stability_solver ==
"coupled_scheme")
380 if(stability_solver ==
"velocity_correction_scheme")
385 std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef = Exp[0]->GetFieldDefinitions();
386 std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
388 for(j = 0; j < nFields; ++j)
390 for(i = 0; i < FieldDef.size(); ++i)
394 FieldDef[i]->m_fields.push_back(
"u");
398 FieldDef[i]->m_fields.push_back(
"v");
400 else if(j == 2 && stability_solver ==
"velocity_correction_scheme")
402 FieldDef[i]->m_fields.push_back(
"p");
404 Exp[j]->AppendFieldData(FieldDef[i], FieldData[i]);
410 std::cout<<
"-------------------------------------------------------------------\n";
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
int main(int argc, char *argv[])
Main.
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< Basis > BasisSharedPtr
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)