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

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

Generated on Wed Nov 24 2021 10:33:29 for Nektar++ by doxygen 1.9.1