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