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 {
447 Exp[2] = Exp2D_pk;
448 }
449
450 //! Expansion coefficient extraction (necessary to write the .fld file)
451 Exp[0]->FwdTrans(Exp2D_uk->GetPhys(), Exp[0]->UpdateCoeffs());
452 Exp[1]->FwdTrans(Exp2D_vk->GetPhys(), Exp[1]->UpdateCoeffs());
453
454 if (stability_solver == "velocity_correction_scheme")
455 {
456 Exp[2]->FwdTrans(Exp2D_pk->GetPhys(), Exp[2]->UpdateCoeffs());
457 }
458 //! --------------------------------------------------------------------------------------------
459
460 //! Generation .FLD file with one field only (at the moment)
461 //! ----------------------------------- Definition of the name of the .fld
462 //! file
463 string FalknerSkan = "FalknerSkan.fld";
464
465 //! Definition of the number of the fields
466 if (stability_solver == "coupled_scheme")
467 {
468 nFields = 2;
469 }
470
471 if (stability_solver == "velocity_correction_scheme")
472 {
473 nFields = 3;
474 }
475
476 //! Definition of the Field
477 std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef =
478 Exp[0]->GetFieldDefinitions();
479 std::vector<std::vector<NekDouble>> FieldData(FieldDef.size());
480
481 for (j = 0; j < nFields; ++j)
482 {
483 for (i = 0; i < FieldDef.size(); ++i)
484 {
485 if (j == 0)
486 {
487 FieldDef[i]->m_fields.push_back("u");
488 }
489 else if (j == 1)
490 {
491 FieldDef[i]->m_fields.push_back("v");
492 }
493 else if (j == 2 && stability_solver == "velocity_correction_scheme")
494 {
495 FieldDef[i]->m_fields.push_back("p");
496 }
497 Exp[j]->AppendFieldData(FieldDef[i], FieldData[i]);
498 }
499 }
500
501 LibUtilities::Write(FalknerSkan, FieldDef, FieldData);
502
503 std::cout << "-------------------------------------------------------------"
504 "------\n";
505 //! ----------------------------------------------------------------------------
506
507 return 0;
508}
L
NekDouble L
Definition: CompressibleBL.cpp:69
Nektar::LibUtilities::Basis
Represents a basis of a given type.
Definition: Basis.h:198
Nektar::LibUtilities::Basis::GetNumModes
int GetNumModes() const
Return order of basis from the basis specification.
Definition: Basis.h:207
Nektar::MemoryManager
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
Definition: NekMemoryManager.hpp:81
Nektar::LibUtilities::BasisSharedPtr
std::shared_ptr< Basis > BasisSharedPtr
Definition: FoundationsFwd.hpp:69
Nektar::LibUtilities::SessionReaderSharedPtr
std::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: SessionReader.h:113
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:245
Nektar::MultiRegions::ExpListSharedPtr
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
Definition: LocTraceToTraceMap.h:56
Nektar::MultiRegions::ContFieldSharedPtr
std::shared_ptr< ContField > ContFieldSharedPtr
Definition: ContField.h:268
Nektar::SpatialDomains::MeshGraphSharedPtr
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:174
Nektar::UnitTests::d
std::vector< double > d(NPUPPER *NPUPPER)
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.hpp:825
tinysimd::sqrt
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:294

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

Generated on Wed Jan 31 2024 09:35:58 for Nektar++ by doxygen 1.9.4