Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
Functions
FldAddFalknerSkanBL.cpp File Reference
#include <cstdio>
#include <cstdlib>
#include <iostream>
#include <iomanip>
#include <LibUtilities/Memory/NekMemoryManager.hpp>
#include <MultiRegions/ExpList.h>
#include <MultiRegions/ExpList2D.h>
#include <MultiRegions/ContField2D.h>
#include <SpatialDomains/MeshGraph2D.h>
Include dependency graph for FldAddFalknerSkanBL.cpp:

Go to the source code of this file.

Functions

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

Function Documentation

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

References L, Vmath::Vcopy(), and Nektar::LibUtilities::Write().

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