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 66 of file FldAddFalknerSkanBL.cpp.

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

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

Generated on Sat Mar 25 2023 12:35:53 for Nektar++ by doxygen 1.9.1