Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
FldAddFalknerSkanBL.cpp
Go to the documentation of this file.
1 /* ==========================================================================
2  * Generation of a .fld file for the Falkner-Skan boundary layer starting
3  * from a session file with some physical data for the definition of the
4  * BL and a .txt file with the similarity solution.
5  * ======================================================================== */
6 
7 /* =====================================
8  * Author: Gianmarco Mengaldo
9  * Generation: dd/mm/aa = 08/03/12
10  * Mantainer: Gianmarco Mengaldo
11 ===================================== */
12 
13 //! Loading cc libraries
14 #include <cstdio>
15 #include <cstdlib>
16 #include <iostream>
17 #include <iomanip>
18 
19 //! Loading Nektar++ libraries
21 #include <MultiRegions/ExpList.h>
22 #include <MultiRegions/ExpList2D.h>
25 
26 //! STL namespace
27 using namespace std;
28 
29 //! Nektar++ namespace
30 using namespace Nektar;
31 
32 //! Main
33 int main(int argc, char *argv[])
34 {
35  //! Setting up the decimal precision to machine precision
36  setprecision (16);
37 
38  //! Auxiliary counters for the x and y directions
39  int i, j, numModes, nFields;
40 
41  //! Usage check
42  if((argc != 2) && (argc != 3))
43  {
44  fprintf(stderr,"Usage: ./FldAddFalknerSkanBL sessionFile [SysSolnType]\n");exit(1);
45  }
46 
47  //! Reading the session file
48  LibUtilities::SessionReaderSharedPtr vSession = LibUtilities::SessionReader::CreateInstance(argc, argv);
49 
50  //! Loading the parameters to define the BL
51  NekDouble Re;
52  NekDouble L;
53  NekDouble U_inf;
54  NekDouble x;
55  NekDouble x_0;
56  NekDouble nu;
57  string BL_type;
58  string txt_file;
59  string stability_solver;
60  int numLines;
61 
62 
63  BL_type = vSession->GetSolverInfo("BL_type");
64  txt_file = vSession->GetSolverInfo("txt_file");
65  stability_solver = vSession->GetSolverInfo("stability_solver");
66 
67  if(stability_solver != "velocity_correction_scheme" &&
68  stability_solver != "coupled_scheme")
69  {
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");
72  exit(1);
73  }
74 
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);
81 
82  //! Check on the physical parameters
83  if(x <= 0)
84  {
85  fprintf(stderr,"Error: x must be positive ===> CHECK the session file\n");
86  exit(1);
87  }
88 
89  if(x_0 < 0)
90  {
91  fprintf(stderr,"Error: x_0 must be positive or at least equal to 0 ===> CHECK the session file\n");
92  exit(1);
93  }
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";
111 
112  //! Computation of the kinematic viscosity
113  nu = U_inf * L / Re;
114 
115  //! Read in mesh from input file and create an object of class MeshGraph2D
118 
119  //! Feed our spatial discretisation object
121  Domain = MemoryManager<MultiRegions::ContField2D>::AllocateSharedPtr(vSession,graphShPt,vSession->GetVariable(0));
122 
123  //! Get the total number of elements
124  int nElements;
125  nElements = Domain->GetExpSize();
126  std::cout << "Number of elements = " << nElements << std::endl;
127 
128  //! Get the total number of quadrature points (depends on n. modes)
129  int nQuadraturePts;
130  nQuadraturePts = Domain->GetTotPoints();
131  std::cout << "Number of quadrature points = " << nQuadraturePts << std::endl;
132 
133  //! Coordinates of the quadrature points
134  Array<OneD,NekDouble> x_QuadraturePts;
135  Array<OneD,NekDouble> y_QuadraturePts;
136  Array<OneD,NekDouble> z_QuadraturePts;
137  x_QuadraturePts = Array<OneD,NekDouble>(nQuadraturePts);
138  y_QuadraturePts = Array<OneD,NekDouble>(nQuadraturePts);
139  z_QuadraturePts = Array<OneD,NekDouble>(nQuadraturePts);
140  Domain->GetCoords(x_QuadraturePts,y_QuadraturePts,z_QuadraturePts);
141 
142  //! Reading the .txt file with eta, f(eta) and f'(eta) -----------------------------------------
143  const char *txt_file_char;
144  //string txt_file(argv[argc-1]);
145  txt_file_char = txt_file.c_str();
146 
147  //! Open the .txt file with the BL data
148  ifstream pFile(txt_file_char);
149  if (!pFile)
150  {
151  cout << "Errro: Unable to open the .txt file with the BL data\n";
152  exit(1);
153  }
154 
155  numLines = numLines/3;
156  NekDouble d;
157  std::vector<std::vector<NekDouble> > GlobalArray (3);
158 
159  for (j=0; j<=2; j++)
160  {
161  GlobalArray[j].resize(numLines);
162  for (i=0; i<=numLines-1; i++)
163  {
164  pFile >> d;
165  GlobalArray[j][i] = d;
166  }
167  }
168  //! --------------------------------------------------------------------------------------------
169 
170 
171  //! Saving eta, f(eta) and f'(eta) into separate arrays ----------------------------------------
175 
176  eta = Array<OneD,NekDouble>(numLines);
177  f = Array<OneD,NekDouble>(numLines);
178  df = Array<OneD,NekDouble>(numLines);
179 
180  for (i=0; i<numLines; i++)
181  {
182  eta[i] = GlobalArray[0][i];
183  f[i] = GlobalArray[1][i];
184  df[i] = GlobalArray[2][i];
185  }
186  //! --------------------------------------------------------------------------------------------
187 
188 
189  //! Inizialisation of the arrays for computations on the Quadrature points ---------------------
190  Array<OneD,NekDouble> eta_QuadraturePts;
191  eta_QuadraturePts = Array<OneD,NekDouble>(nQuadraturePts);
192 
193  Array<OneD,NekDouble> f_QuadraturePts;
194  f_QuadraturePts = Array<OneD,NekDouble>(nQuadraturePts);
195 
196  Array<OneD,NekDouble> df_QuadraturePts;
197  df_QuadraturePts = Array<OneD,NekDouble>(nQuadraturePts);
198 
199  Array<OneD,NekDouble> u_QuadraturePts;
200  u_QuadraturePts = Array<OneD,NekDouble>(nQuadraturePts);
201 
202  Array<OneD,NekDouble> v_QuadraturePts;
203  v_QuadraturePts = Array<OneD,NekDouble>(nQuadraturePts);
204 
205  Array<OneD,NekDouble> p_QuadraturePts;
206  p_QuadraturePts = Array<OneD,NekDouble>(nQuadraturePts);
207  //! --------------------------------------------------------------------------------------------
208 
209 
210 
211  //! PARALLEL CASE ------------------------------------------------------------------------------
212  if(BL_type == "Parallel")
213  {
214  for(i=0; i<nQuadraturePts; i++)
215  {
216  eta_QuadraturePts[i] = y_QuadraturePts[i] * sqrt(U_inf / (2 * x * nu));
217  for(j=0; j<numLines-1; j++)
218  {
219  if(eta_QuadraturePts[i] >= eta[j] && eta_QuadraturePts[i] <= eta[j+1])
220  {
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];
223  }
224 
225  else if(eta_QuadraturePts[i] == 1000000)
226  {
227  f_QuadraturePts[i] = f[numLines-1];
228  df_QuadraturePts[i] = df[numLines-1];
229  }
230 
231  else if(eta_QuadraturePts[i] > eta[numLines-1])
232  {
233  f_QuadraturePts[i] = f[numLines-1] + df[numLines-1] * (eta_QuadraturePts[i] - eta[numLines-1]);
234  df_QuadraturePts[i] = df[numLines-1];
235  }
236 
237  }
238 
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;
242  }
243  }
244  //! --------------------------------------------------------------------------------------------
245 
246 
247 
248  //! NON-PARALLEL CASE --------------------------------------------------------------------------
249  if(BL_type == "nonParallel")
250  {
251  for(i=0; i<nQuadraturePts; i++)
252  {
253  eta_QuadraturePts[i] = y_QuadraturePts[i] * sqrt(U_inf / (2 * (x_QuadraturePts[i] + x_0) * nu));
254 
255  if((x_QuadraturePts[i] + x_0) == 0)
256  {
257  eta_QuadraturePts[i] = 1000000;
258  }
259 
260  for(j=0; j<numLines-1; j++)
261  {
262  if(eta_QuadraturePts[i] >= eta[j] && eta_QuadraturePts[i] <= eta[j+1])
263  {
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];
266  }
267 
268  else if(eta_QuadraturePts[i] == 1000000)
269  {
270  f_QuadraturePts[i] = f[numLines-1];
271  df_QuadraturePts[i] = df[numLines-1];
272  }
273 
274  else if(eta_QuadraturePts[i] > eta[numLines-1])
275  {
276  f_QuadraturePts[i] = f[numLines-1] + df[numLines-1] * (eta_QuadraturePts[i] - eta[numLines-1]);
277  df_QuadraturePts[i] = df[numLines-1];
278  }
279  }
280 
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]);
285 
286  //! INFLOW SECTION: X = 0; Y > 0.
287  if((x_QuadraturePts[i] + x_0) == 0)
288  {
289  u_QuadraturePts[i] = U_inf;
290  v_QuadraturePts[i] = 0.0;
291  }
292 
293  //! SINGULARITY POINT: X = 0; Y = 0.
294  if((x_QuadraturePts[i] + x_0) == 0 && y_QuadraturePts[i] == 0)
295  {
296  u_QuadraturePts[i] = 0.0;
297  v_QuadraturePts[i] = 0.0;
298  }
299 
300  p_QuadraturePts[i] = 0.0;
301  }
302  }
303  //! --------------------------------------------------------------------------------------------
304 
305 
306 
307  //! Inspection of the interpolation ------------------------------------------------------------
308  FILE *etaOriginal;
309  etaOriginal = fopen("eta.txt","w+");
310  for(i=0; i<numLines; i++)
311  {
312  fprintf(etaOriginal,"%f %f %f \n", eta[i], f[i], df[i]);
313  }
314  fclose(etaOriginal);
315 
316 
317  FILE *yQ;
318  yQ = fopen("yQ.txt","w+");
319  for(i=0; i<nQuadraturePts; i++)
320  {
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]);
323  }
324  fclose(yQ);
325  //! --------------------------------------------------------------------------------------------
326 
327 
328 
329  //! Definition of the 2D expansion using the mesh data specified on the session file -----------
331  Exp2D_uk = MemoryManager<MultiRegions::ExpList2D>::AllocateSharedPtr(vSession,graphShPt);
332 
334  Exp2D_vk = MemoryManager<MultiRegions::ExpList2D>::AllocateSharedPtr(vSession,graphShPt);
335 
337  Exp2D_pk = MemoryManager<MultiRegions::ExpList2D>::AllocateSharedPtr(vSession,graphShPt);
338  //! --------------------------------------------------------------------------------------------
339 
340 
341 
342  //! Filling the 2D expansion using a recursive algorithm based on the mesh ordering ------------
344  Basis = Domain->GetExp(0)->GetBasis(0);
345  numModes = Basis->GetNumModes();
346 
347  std::cout<< "Number of modes = " << numModes << std::endl;
348 
349  //! Copying the ukGlobal vector (with the same pattern of m_phys) in m_phys
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);
353 
354  //! Initialisation of the ExpList Exp
356  Exp[0] = Exp2D_uk;
357  Exp[1] = Exp2D_vk;
358 
359  if(stability_solver == "velocity_correction_scheme")
360  Exp[2] = Exp2D_pk;
361 
362  //! Expansion coefficient extraction (necessary to write the .fld file)
363  Exp[0]->FwdTrans(Exp2D_uk->GetPhys(),Exp[0]->UpdateCoeffs());
364  Exp[1]->FwdTrans(Exp2D_vk->GetPhys(),Exp[1]->UpdateCoeffs());
365 
366  if(stability_solver == "velocity_correction_scheme")
367  Exp[2]->FwdTrans(Exp2D_pk->GetPhys(),Exp[2]->UpdateCoeffs());
368  //! --------------------------------------------------------------------------------------------
369 
370 
371 
372  //! Generation .FLD file with one field only (at the moment) -----------------------------------
373  //! Definition of the name of the .fld file
374  string FalknerSkan = "FalknerSkan.fld";
375 
376  //! Definition of the number of the fields
377  if(stability_solver == "coupled_scheme")
378  nFields = 2;
379 
380  if(stability_solver == "velocity_correction_scheme")
381  nFields = 3;
382 
383 
384  //! Definition of the Field
385  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef = Exp[0]->GetFieldDefinitions();
386  std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
387 
388  for(j = 0; j < nFields; ++j)
389  {
390  for(i = 0; i < FieldDef.size(); ++i)
391  {
392  if(j == 0)
393  {
394  FieldDef[i]->m_fields.push_back("u");
395  }
396  else if(j == 1)
397  {
398  FieldDef[i]->m_fields.push_back("v");
399  }
400  else if(j == 2 && stability_solver == "velocity_correction_scheme")
401  {
402  FieldDef[i]->m_fields.push_back("p");
403  }
404  Exp[j]->AppendFieldData(FieldDef[i], FieldData[i]);
405  }
406  }
407 
408  LibUtilities::Write(FalknerSkan, FieldDef, FieldData);
409 
410  std::cout<<"-------------------------------------------------------------------\n";
411  //! ----------------------------------------------------------------------------
412 
413  return 0;
414 }
415 
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
STL namespace.
boost::shared_ptr< ContField2D > ContField2DSharedPtr
Definition: ContField2D.h:287
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
int main(int argc, char *argv[])
Main.
double NekDouble
boost::shared_ptr< ExpList2D > ExpList2DSharedPtr
Shared pointer to an ExpList2D object.
Definition: ExpList2D.h:49
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 ...
Definition: FieldIO.cpp:235
NekDouble L
boost::shared_ptr< Basis > BasisSharedPtr
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:442
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1061