Nektar++
FldAddFalknerSkanBL.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: FldAddFalknerSkanBL.cpp
4//
5// For more information, please see: http://www.nektar.info
6//
7// The MIT License
8//
9// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10// Department of Aeronautics, Imperial College London (UK), and Scientific
11// Computing and Imaging Institute, University of Utah (USA).
12//
13// Permission is hereby granted, free of charge, to any person obtaining a
14// copy of this software and associated documentation files (the "Software"),
15// to deal in the Software without restriction, including without limitation
16// the rights to use, copy, modify, merge, publish, distribute, sublicense,
17// and/or sell copies of the Software, and to permit persons to whom the
18// Software is furnished to do so, subject to the following conditions:
19//
20// The above copyright notice and this permission notice shall be included
21// in all copies or substantial portions of the Software.
22//
23// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29// DEALINGS IN THE SOFTWARE.
30//
31// Description:
32//
33///////////////////////////////////////////////////////////////////////////////
34
35/* ==========================================================================
36 * Generation of a .fld file for the Falkner-Skan boundary layer starting
37 * from a session file with some physical data for the definition of the
38 * BL and a .txt file with the similarity solution.
39 * ======================================================================== */
40
41/* =====================================
42 * Author: Gianmarco Mengaldo
43 * Generation: dd/mm/aa = 08/03/12
44 * Mantainer: Gianmarco Mengaldo
45===================================== */
46
47//! Loading cc libraries
48#include <cstdio>
49#include <cstdlib>
50#include <iomanip>
51#include <iostream>
52
53//! Loading Nektar++ libraries
58
59//! STL namespace
60using namespace std;
61
62//! Nektar++ namespace
63using namespace Nektar;
64
65//! Main
66int main(int argc, char *argv[])
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
84 LibUtilities::SessionReader::CreateInstance(argc, argv);
85
86 //! Loading the parameters to define the BL
87 NekDouble Re;
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
171 graphShPt = SpatialDomains::MeshGraphIO::Read(vSession);
172
173 //! Feed our spatial discretisation object
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 //! ----------------------------------------
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 --
415 vSession, graphShPt);
416
419 vSession, graphShPt);
420
423 vSession, graphShPt);
424 //! -----------------------------------------------------------------------------------
425
426 //! Filling the 2D expansion using a recursive algorithm based on the mesh
427 //! ordering ------------
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
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}
NekDouble L
int main(int argc, char *argv[])
Main.
Represents a basis of a given type.
Definition: Basis.h:198
int GetNumModes() const
Return order of basis from the basis specification.
Definition: Basis.h:207
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
std::shared_ptr< Basis > BasisSharedPtr
std::shared_ptr< SessionReader > SessionReaderSharedPtr
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
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
std::shared_ptr< ContField > ContFieldSharedPtr
Definition: ContField.h:268
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:174
std::vector< double > d(NPUPPER *NPUPPER)
double NekDouble
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:825
STL namespace.
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:294