Nektar++
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 
27 //! Nektar++ namespace
28 using namespace Nektar;
29 
30 //! Main
31 int main(int argc, char *argv[])
32 {
33  //! Setting up the decimal precision to machine precision
34  setprecision (16);
35 
36  //! Auxiliary counters for the x and y directions
37  int i, j, numModes, nFields;
38 
39  //! Usage check
40  if((argc != 2) && (argc != 3))
41  {
42  fprintf(stderr,"Usage: ./FldAddFalknerSkanBL sessionFile [SysSolnType]\n");exit(1);
43  }
44 
45  //! Reading the session file
47 
48  //! Loading the parameters to define the BL
49  NekDouble Re;
50  NekDouble L;
51  NekDouble U_inf;
52  NekDouble x;
53  NekDouble x_0;
54  NekDouble nu;
55  string BL_type;
56  string txt_file;
57  string stability_solver;
58  int numLines;
59 
60 
61  BL_type = vSession->GetSolverInfo("BL_type");
62  txt_file = vSession->GetSolverInfo("txt_file");
63  stability_solver = vSession->GetSolverInfo("stability_solver");
64 
65  if(stability_solver != "velocity_correction_scheme" &&
66  stability_solver != "coupled_scheme")
67  {
68  fprintf(stderr,"Error: You must specify the stability solver in the session file properly.\n");
69  fprintf(stderr,"Options: 'velocity_correction_scheme' [output ===> (u,v,p)]; 'coupled_scheme' [output ===>(u,v)]\n");
70  exit(1);
71  }
72 
73  vSession->LoadParameter("Re", Re, 1.0);
74  vSession->LoadParameter("L", L, 1.0);
75  vSession->LoadParameter("U_inf", U_inf, 1.0);
76  vSession->LoadParameter("x", x, 1.0);
77  vSession->LoadParameter("x_0", x_0, 1.0);
78  vSession->LoadParameter("NumberLines_txt", numLines, 1.0);
79 
80  //! Check on the physical parameters
81  if(x <= 0)
82  {
83  fprintf(stderr,"Error: x must be positive ===> CHECK the session file\n");
84  exit(1);
85  }
86 
87  if(x_0 < 0)
88  {
89  fprintf(stderr,"Error: x_0 must be positive or at least equal to 0 ===> CHECK the session file\n");
90  exit(1);
91  }
92  std::cout<<"\n=========================================================================\n";
93  std::cout<<"Falkner-Skan Boundary Layer Generation (version of July 12th 2012)\n";
94  std::cout<<"=========================================================================\n";
95  std::cout<<"*************************************************************************\n";
96  std::cout<<"DATA FROM THE SESSION FILE:\n";
97  std::cout << "Reynolds number = " << Re << "\t[-]" << std::endl;
98  std::cout << "Characteristic length = " << L << "\t\t[m]" << std::endl;
99  std::cout << "U_infinity = " << U_inf << "\t[m/s]" << std::endl;
100  std::cout << "Position x (parallel case only) = " << x << "\t\t[m]" << std::endl;
101  std::cout << "Position x_0 to start the BL [m] = " << x_0 << "\t\t[m]" << std::endl;
102  std::cout << "Number of lines of the .txt file = " << numLines << "\t[-]" << std::endl;
103  std::cout << "BL type = " << BL_type << std::endl;
104  std::cout << ".txt file read = " << txt_file << std::endl;
105  std::cout << "Stability solver = " << stability_solver << std::endl;
106  std::cout<<"*************************************************************************\n";
107  std::cout<<"-------------------------------------------------------------------------\n";
108  std::cout<<"MESH and EXPANSION DATA:\n";
109 
110  //! Computation of the kinematic viscosity
111  nu = U_inf * L / Re;
112 
113  //! Read in mesh from input file and create an object of class MeshGraph2D
116 
117  //! Feed our spatial discretisation object
119  Domain = MemoryManager<MultiRegions::ContField2D>::AllocateSharedPtr(vSession,graphShPt,vSession->GetVariable(0));
120 
121  //! Get the total number of elements
122  int nElements;
123  nElements = Domain->GetExpSize();
124  std::cout << "Number of elements = " << nElements << std::endl;
125 
126  //! Get the total number of quadrature points (depends on n. modes)
127  int nQuadraturePts;
128  nQuadraturePts = Domain->GetTotPoints();
129  std::cout << "Number of quadrature points = " << nQuadraturePts << std::endl;
130 
131  //! Coordinates of the quadrature points
132  Array<OneD,NekDouble> x_QuadraturePts;
133  Array<OneD,NekDouble> y_QuadraturePts;
134  Array<OneD,NekDouble> z_QuadraturePts;
135  x_QuadraturePts = Array<OneD,NekDouble>(nQuadraturePts);
136  y_QuadraturePts = Array<OneD,NekDouble>(nQuadraturePts);
137  z_QuadraturePts = Array<OneD,NekDouble>(nQuadraturePts);
138  Domain->GetCoords(x_QuadraturePts,y_QuadraturePts,z_QuadraturePts);
139 
140  //! Reading the .txt file with eta, f(eta) and f'(eta) -----------------------------------------
141  const char *txt_file_char;
142  //string txt_file(argv[argc-1]);
143  txt_file_char = txt_file.c_str();
144 
145  //! Open the .txt file with the BL data
146  ifstream pFile(txt_file_char);
147  if (!pFile)
148  {
149  cout << "Errro: Unable to open the .txt file with the BL data\n";
150  exit(1);
151  }
152 
153  numLines = numLines/3;
154  NekDouble d;
155  std::vector<std::vector<NekDouble> > GlobalArray (3);
156 
157  for (j=0; j<=2; j++)
158  {
159  GlobalArray[j].resize(numLines);
160  for (i=0; i<=numLines-1; i++)
161  {
162  pFile >> d;
163  GlobalArray[j][i] = d;
164  }
165  }
166  //! --------------------------------------------------------------------------------------------
167 
168 
169  //! Saving eta, f(eta) and f'(eta) into separate arrays ----------------------------------------
173 
174  eta = Array<OneD,NekDouble>(numLines);
175  f = Array<OneD,NekDouble>(numLines);
176  df = Array<OneD,NekDouble>(numLines);
177 
178  for (i=0; i<numLines; i++)
179  {
180  eta[i] = GlobalArray[0][i];
181  f[i] = GlobalArray[1][i];
182  df[i] = GlobalArray[2][i];
183  }
184  //! --------------------------------------------------------------------------------------------
185 
186 
187  //! Inizialisation of the arrays for computations on the Quadrature points ---------------------
188  Array<OneD,NekDouble> eta_QuadraturePts;
189  eta_QuadraturePts = Array<OneD,NekDouble>(nQuadraturePts);
190 
191  Array<OneD,NekDouble> f_QuadraturePts;
192  f_QuadraturePts = Array<OneD,NekDouble>(nQuadraturePts);
193 
194  Array<OneD,NekDouble> df_QuadraturePts;
195  df_QuadraturePts = Array<OneD,NekDouble>(nQuadraturePts);
196 
197  Array<OneD,NekDouble> u_QuadraturePts;
198  u_QuadraturePts = Array<OneD,NekDouble>(nQuadraturePts);
199 
200  Array<OneD,NekDouble> v_QuadraturePts;
201  v_QuadraturePts = Array<OneD,NekDouble>(nQuadraturePts);
202 
203  Array<OneD,NekDouble> p_QuadraturePts;
204  p_QuadraturePts = Array<OneD,NekDouble>(nQuadraturePts);
205  //! --------------------------------------------------------------------------------------------
206 
207 
208 
209  //! PARALLEL CASE ------------------------------------------------------------------------------
210  if(BL_type == "Parallel")
211  {
212  for(i=0; i<nQuadraturePts; i++)
213  {
214  eta_QuadraturePts[i] = y_QuadraturePts[i] * sqrt(U_inf / (2 * x * nu));
215  for(j=0; j<numLines-1; j++)
216  {
217  if(eta_QuadraturePts[i] >= eta[j] && eta_QuadraturePts[i] <= eta[j+1])
218  {
219  f_QuadraturePts[i] = (eta_QuadraturePts[i] - eta[j]) * (f[j+1] - f[j]) / (eta[j+1] - eta[j]) + f[j];
220  df_QuadraturePts[i] = (eta_QuadraturePts[i] - eta[j]) * (df[j+1] - df[j]) / (eta[j+1] - eta[j]) + df[j];
221  }
222 
223  else if(eta_QuadraturePts[i] == 1000000)
224  {
225  f_QuadraturePts[i] = f[numLines-1];
226  df_QuadraturePts[i] = df[numLines-1];
227  }
228 
229  else if(eta_QuadraturePts[i] > eta[numLines-1])
230  {
231  f_QuadraturePts[i] = f[numLines-1] + df[numLines-1] * (eta_QuadraturePts[i] - eta[numLines-1]);
232  df_QuadraturePts[i] = df[numLines-1];
233  }
234 
235  }
236 
237  u_QuadraturePts[i] = U_inf * df_QuadraturePts[i];
238  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]);
239  p_QuadraturePts[i] = 0.0;
240  }
241  }
242  //! --------------------------------------------------------------------------------------------
243 
244 
245 
246  //! NON-PARALLEL CASE --------------------------------------------------------------------------
247  if(BL_type == "nonParallel")
248  {
249  for(i=0; i<nQuadraturePts; i++)
250  {
251  eta_QuadraturePts[i] = y_QuadraturePts[i] * sqrt(U_inf / (2 * (x_QuadraturePts[i] + x_0) * nu));
252 
253  if((x_QuadraturePts[i] + x_0) == 0)
254  {
255  eta_QuadraturePts[i] = 1000000;
256  }
257 
258  for(j=0; j<numLines-1; j++)
259  {
260  if(eta_QuadraturePts[i] >= eta[j] && eta_QuadraturePts[i] <= eta[j+1])
261  {
262  f_QuadraturePts[i] = (eta_QuadraturePts[i] - eta[j]) * (f[j+1] - f[j]) / (eta[j+1] - eta[j]) + f[j];
263  df_QuadraturePts[i] = (eta_QuadraturePts[i] - eta[j]) * (df[j+1] - df[j]) / (eta[j+1] - eta[j]) + df[j];
264  }
265 
266  else if(eta_QuadraturePts[i] == 1000000)
267  {
268  f_QuadraturePts[i] = f[numLines-1];
269  df_QuadraturePts[i] = df[numLines-1];
270  }
271 
272  else if(eta_QuadraturePts[i] > eta[numLines-1])
273  {
274  f_QuadraturePts[i] = f[numLines-1] + df[numLines-1] * (eta_QuadraturePts[i] - eta[numLines-1]);
275  df_QuadraturePts[i] = df[numLines-1];
276  }
277  }
278 
279  u_QuadraturePts[i] = U_inf * df_QuadraturePts[i];
280  v_QuadraturePts[i] = nu * sqrt(U_inf / (2.0 * nu * (x_QuadraturePts[i] + x_0))) *
281  (y_QuadraturePts[i] * sqrt(U_inf / (2.0 * nu * (x_QuadraturePts[i] + x_0))) *
282  df_QuadraturePts[i] - f_QuadraturePts[i]);
283 
284  //! INFLOW SECTION: X = 0; Y > 0.
285  if((x_QuadraturePts[i] + x_0) == 0)
286  {
287  u_QuadraturePts[i] = U_inf;
288  v_QuadraturePts[i] = 0.0;
289  }
290 
291  //! SINGULARITY POINT: X = 0; Y = 0.
292  if((x_QuadraturePts[i] + x_0) == 0 && y_QuadraturePts[i] == 0)
293  {
294  u_QuadraturePts[i] = 0.0;
295  v_QuadraturePts[i] = 0.0;
296  }
297 
298  p_QuadraturePts[i] = 0.0;
299  }
300  }
301  //! --------------------------------------------------------------------------------------------
302 
303 
304 
305  //! Inspection of the interpolation ------------------------------------------------------------
306  FILE *etaOriginal;
307  etaOriginal = fopen("eta.txt","w+");
308  for(i=0; i<numLines; i++)
309  {
310  fprintf(etaOriginal,"%f %f %f \n", eta[i], f[i], df[i]);
311  }
312  fclose(etaOriginal);
313 
314 
315  FILE *yQ;
316  yQ = fopen("yQ.txt","w+");
317  for(i=0; i<nQuadraturePts; i++)
318  {
319  fprintf(yQ,"%f %f %f %f %f %f %f\n", x_QuadraturePts[i], y_QuadraturePts[i], u_QuadraturePts[i],
320  v_QuadraturePts[i], eta_QuadraturePts[i], f_QuadraturePts[i], df_QuadraturePts[i]);
321  }
322  fclose(yQ);
323  //! --------------------------------------------------------------------------------------------
324 
325 
326 
327  //! Definition of the 2D expansion using the mesh data specified on the session file -----------
329  Exp2D_uk = MemoryManager<MultiRegions::ExpList2D>::AllocateSharedPtr(vSession,graphShPt);
330 
332  Exp2D_vk = MemoryManager<MultiRegions::ExpList2D>::AllocateSharedPtr(vSession,graphShPt);
333 
335  Exp2D_pk = MemoryManager<MultiRegions::ExpList2D>::AllocateSharedPtr(vSession,graphShPt);
336  //! --------------------------------------------------------------------------------------------
337 
338 
339 
340  //! Filling the 2D expansion using a recursive algorithm based on the mesh ordering ------------
342  Basis = Domain->GetExp(0)->GetBasis(0);
343  numModes = Basis->GetNumModes();
344 
345  std::cout<< "Number of modes = " << numModes << std::endl;
346 
347  //! Copying the ukGlobal vector (with the same pattern of m_phys) in m_phys
348  Vmath::Vcopy(nQuadraturePts, u_QuadraturePts , 1, Exp2D_uk->UpdatePhys(), 1);
349  Vmath::Vcopy(nQuadraturePts, v_QuadraturePts, 1, Exp2D_vk->UpdatePhys(), 1);
350  Vmath::Vcopy(nQuadraturePts, p_QuadraturePts , 1, Exp2D_pk->UpdatePhys(), 1);
351 
352  //! Initialisation of the ExpList Exp
354  Exp[0] = Exp2D_uk;
355  Exp[1] = Exp2D_vk;
356 
357  if(stability_solver == "velocity_correction_scheme")
358  Exp[2] = Exp2D_pk;
359 
360  //! Expansion coefficient extraction (necessary to write the .fld file)
361  Exp[0]->FwdTrans(Exp2D_uk->GetPhys(),Exp[0]->UpdateCoeffs());
362  Exp[1]->FwdTrans(Exp2D_vk->GetPhys(),Exp[1]->UpdateCoeffs());
363 
364  if(stability_solver == "velocity_correction_scheme")
365  Exp[2]->FwdTrans(Exp2D_pk->GetPhys(),Exp[2]->UpdateCoeffs());
366  //! --------------------------------------------------------------------------------------------
367 
368 
369 
370  //! Generation .FLD file with one field only (at the moment) -----------------------------------
371  //! Definition of the name of the .fld file
372  string FalknerSkan = "FalknerSkan.fld";
373 
374  //! Definition of the number of the fields
375  if(stability_solver == "coupled_scheme")
376  nFields = 2;
377 
378  if(stability_solver == "velocity_correction_scheme")
379  nFields = 3;
380 
381 
382  //! Definition of the Field
383  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef = Exp[0]->GetFieldDefinitions();
384  std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
385 
386  for(j = 0; j < nFields; ++j)
387  {
388  for(i = 0; i < FieldDef.size(); ++i)
389  {
390  if(j == 0)
391  {
392  FieldDef[i]->m_fields.push_back("u");
393  }
394  else if(j == 1)
395  {
396  FieldDef[i]->m_fields.push_back("v");
397  }
398  else if(j == 2 && stability_solver == "velocity_correction_scheme")
399  {
400  FieldDef[i]->m_fields.push_back("p");
401  }
402  Exp[j]->AppendFieldData(FieldDef[i], FieldData[i]);
403  }
404  }
405 
406  LibUtilities::Write(FalknerSkan, FieldDef, FieldData);
407 
408  std::cout<<"-------------------------------------------------------------------\n";
409  //! ----------------------------------------------------------------------------
410 
411  return 0;
412 }
413 
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
boost::shared_ptr< ContField2D > ContField2DSharedPtr
Definition: ContField2D.h:291
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:50
static SessionReaderSharedPtr CreateInstance(int argc, char *argv[])
Creates an instance of the SessionReader class.
int main(int argc, char *argv[])
Main.
double NekDouble
boost::shared_ptr< ExpList2D > ExpList2DSharedPtr
Shared pointer to an ExpList2D object.
Definition: ExpList2D.h:49
NekDouble L
void Write(const std::string &outFile, std::vector< FieldDefinitionsSharedPtr > &fielddefs, std::vector< std::vector< NekDouble > > &fielddata, const FieldMetaDataMap &fieldinfomap)
Write a field file in serial only.
Definition: FieldIO.cpp:72
boost::shared_ptr< Basis > BasisSharedPtr
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:432
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1038