Nektar++
Functions
FldAddFalknerSkanBL.cpp File Reference
#include <cstdio>
#include <cstdlib>
#include <iomanip>
#include <iostream>
#include <LibUtilities/Memory/NekMemoryManager.hpp>
#include <MultiRegions/ContField.h>
#include <MultiRegions/ExpList.h>
#include <SpatialDomains/MeshGraph.h>

Go to the source code of this file.

Functions

int main (int argc, char *argv[])
 Main. More...
 

Function Documentation

◆ main()

int main ( int  argc,
char *  argv[] 
)

Main.

Setting up the decimal precision to machine precision

Auxiliary counters for the x and y directions

Usage check

Reading the session file

Loading the parameters to define the BL

Check on the physical parameters

Computation of the kinematic viscosity

Read in mesh from input file and create an object of class MeshGraph2D

Feed our spatial discretisation object

Get the total number of elements

Get the total number of quadrature points (depends on n. modes)

Coordinates of the quadrature points

Reading the .txt file with eta, f(eta) and f'(eta)

Open the .txt file with the BL data

Saving eta, f(eta) and f'(eta) into separate arrays

Inizialisation of the arrays for computations on the Quadrature points

PARALLEL CASE

NON-PARALLEL CASE

INFLOW SECTION: X = 0; Y > 0.

SINGULARITY POINT: X = 0; Y = 0.

Inspection of the interpolation

Definition of the 2D expansion using the mesh data specified on the session file –

Filling the 2D expansion using a recursive algorithm based on the mesh ordering ---------—

Copying the ukGlobal vector (with the same pattern of m_phys) in m_phys

Initialisation of the ExpList Exp

Expansion coefficient extraction (necessary to write the .fld file)

Generation .FLD file with one field only (at the moment) --------------------------------— Definition of the name of the .fld file

Definition of the number of the fields

Definition of the Field

Definition at line 32 of file FldAddFalknerSkanBL.cpp.

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
49  LibUtilities::SessionReaderSharedPtr vSession =
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
136  SpatialDomains::MeshGraphSharedPtr graphShPt;
137  graphShPt = SpatialDomains::MeshGraph::Read(vSession);
138 
139  //! Feed our spatial discretisation object
140  MultiRegions::ContFieldSharedPtr Domain;
141  Domain = MemoryManager<MultiRegions::ContField>::AllocateSharedPtr(
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  //! ----------------------------------------
196  Array<OneD, NekDouble> eta;
197  Array<OneD, NekDouble> f;
198  Array<OneD, NekDouble> df;
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 --
379  MultiRegions::ExpListSharedPtr Exp2D_uk;
380  Exp2D_uk = MemoryManager<MultiRegions::ExpList>::AllocateSharedPtr(
381  vSession, graphShPt);
382 
383  MultiRegions::ExpListSharedPtr Exp2D_vk;
384  Exp2D_vk = MemoryManager<MultiRegions::ExpList>::AllocateSharedPtr(
385  vSession, graphShPt);
386 
387  MultiRegions::ExpListSharedPtr Exp2D_pk;
388  Exp2D_pk = MemoryManager<MultiRegions::ExpList>::AllocateSharedPtr(
389  vSession, graphShPt);
390  //! -----------------------------------------------------------------------------------
391 
392  //! Filling the 2D expansion using a recursive algorithm based on the mesh
393  //! ordering ------------
394  LibUtilities::BasisSharedPtr Basis;
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
407  Array<OneD, MultiRegions::ExpListSharedPtr> Exp(3);
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 }
L
NekDouble L
Definition: CompressibleBL.cpp:74
Nektar::LibUtilities::Basis
Represents a basis of a given type.
Definition: Basis.h:213
Nektar::LibUtilities::Basis::GetNumModes
int GetNumModes() const
Return order of basis from the basis specification.
Definition: Basis.h:222
Nektar::MemoryManager
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
Definition: NekMemoryManager.hpp:83
Nektar::LibUtilities::Write
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
Nektar::LibUtilities::BasisSharedPtr
std::shared_ptr< Basis > BasisSharedPtr
Definition: FoundationsFwd.hpp:71
Nektar::LibUtilities::SessionReaderSharedPtr
std::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: SessionReader.h:115
Nektar::MultiRegions::ExpListSharedPtr
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
Definition: LocTraceToTraceMap.h:55
Nektar::MultiRegions::ContFieldSharedPtr
std::shared_ptr< ContField > ContFieldSharedPtr
Definition: ContField.h:277
Nektar::SpatialDomains::MeshGraphSharedPtr
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:172
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1255
tinysimd::sqrt
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:291

References Nektar::LibUtilities::Basis::GetNumModes(), L, tinysimd::sqrt(), Vmath::Vcopy(), and Nektar::LibUtilities::Write().

Generated on Tue Aug 23 2022 07:27:16 for Nektar++ by doxygen 1.9.1