Nektar++
`#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>`

Go to the source code of this file.

Functions

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

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

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)

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

{
//! Setting up the decimal precision to machine precision
setprecision (16);
//! Auxiliary counters for the x and y directions
int i, j, numModes, nFields;
//! Usage check
if((argc != 2) && (argc != 3))
{
}
NekDouble U_inf;
NekDouble x_0;
string BL_type;
string txt_file;
string stability_solver;
int numLines;
BL_type = vSession->GetSolverInfo("BL_type");
txt_file = vSession->GetSolverInfo("txt_file");
stability_solver = vSession->GetSolverInfo("stability_solver");
if(stability_solver != "velocity_correction_scheme" &&
stability_solver != "coupled_scheme")
{
fprintf(stderr,"Error: You must specify the stability solver in the session file properly.\n");
fprintf(stderr,"Options: 'velocity_correction_scheme' [output ===> (u,v,p)]; 'coupled_scheme' [output ===>(u,v)]\n");
exit(1);
}
//! Check on the physical parameters
if(x <= 0)
{
fprintf(stderr,"Error: x must be positive ===> CHECK the session file\n");
exit(1);
}
if(x_0 < 0)
{
fprintf(stderr,"Error: x_0 must be positive or at least equal to 0 ===> CHECK the session file\n");
exit(1);
}
std::cout<<"\n=========================================================================\n";
std::cout<<"Falkner-Skan Boundary Layer Generation (version of July 12th 2012)\n";
std::cout<<"=========================================================================\n";
std::cout<<"*************************************************************************\n";
std::cout<<"DATA FROM THE SESSION FILE:\n";
std::cout << "Reynolds number = " << Re << "\t[-]" << std::endl;
std::cout << "Characteristic length = " << L << "\t\t[m]" << std::endl;
std::cout << "U_infinity = " << U_inf << "\t[m/s]" << std::endl;
std::cout << "Position x (parallel case only) = " << x << "\t\t[m]" << std::endl;
std::cout << "Position x_0 to start the BL [m] = " << x_0 << "\t\t[m]" << std::endl;
std::cout << "Number of lines of the .txt file = " << numLines << "\t[-]" << std::endl;
std::cout << "BL type = " << BL_type << std::endl;
std::cout << ".txt file read = " << txt_file << std::endl;
std::cout << "Stability solver = " << stability_solver << std::endl;
std::cout<<"*************************************************************************\n";
std::cout<<"-------------------------------------------------------------------------\n";
std::cout<<"MESH and EXPANSION DATA:\n";
//! Computation of the kinematic viscosity
nu = U_inf * L / Re;
//! Read in mesh from input file and create an object of class MeshGraph2D
//! Feed our spatial discretisation object
Domain = MemoryManager<MultiRegions::ContField2D>::AllocateSharedPtr(vSession,graphShPt,vSession->GetVariable(0));
//! Get the total number of elements
int nElements;
nElements = Domain->GetExpSize();
std::cout << "Number of elements = " << nElements << std::endl;
//! Get the total number of quadrature points (depends on n. modes)
std::cout << "Number of quadrature points = " << nQuadraturePts << std::endl;
//! Coordinates of the quadrature points
//! Reading the .txt file with eta, f(eta) and f'(eta) -----------------------------------------
const char *txt_file_char;
//string txt_file(argv[argc-1]);
txt_file_char = txt_file.c_str();
//! Open the .txt file with the BL data
ifstream pFile(txt_file_char);
if (!pFile)
{
cout << "Errro: Unable to open the .txt file with the BL data\n";
exit(1);
}
numLines = numLines/3;
std::vector<std::vector<NekDouble> > GlobalArray (3);
for (j=0; j<=2; j++)
{
GlobalArray[j].resize(numLines);
for (i=0; i<=numLines-1; i++)
{
pFile >> d;
GlobalArray[j][i] = d;
}
}
//! --------------------------------------------------------------------------------------------
//! Saving eta, f(eta) and f'(eta) into separate arrays ----------------------------------------
Array<OneD,NekDouble> eta;
Array<OneD,NekDouble> f;
Array<OneD,NekDouble> df;
eta = Array<OneD,NekDouble>(numLines);
f = Array<OneD,NekDouble>(numLines);
df = Array<OneD,NekDouble>(numLines);
for (i=0; i<numLines; i++)
{
eta[i] = GlobalArray[0][i];
f[i] = GlobalArray[1][i];
df[i] = GlobalArray[2][i];
}
//! --------------------------------------------------------------------------------------------
//! Inizialisation of the arrays for computations on the Quadrature points ---------------------
//! --------------------------------------------------------------------------------------------
//! PARALLEL CASE ------------------------------------------------------------------------------
if(BL_type == "Parallel")
{
{
for(j=0; j<numLines-1; j++)
{
{
f_QuadraturePts[i] = (eta_QuadraturePts[i] - eta[j]) * (f[j+1] - f[j]) / (eta[j+1] - eta[j]) + f[j];
df_QuadraturePts[i] = (eta_QuadraturePts[i] - eta[j]) * (df[j+1] - df[j]) / (eta[j+1] - eta[j]) + df[j];
}
{
}
{
}
}
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]);
}
}
//! --------------------------------------------------------------------------------------------
//! NON-PARALLEL CASE --------------------------------------------------------------------------
if(BL_type == "nonParallel")
{
{
{
}
for(j=0; j<numLines-1; j++)
{
{
f_QuadraturePts[i] = (eta_QuadraturePts[i] - eta[j]) * (f[j+1] - f[j]) / (eta[j+1] - eta[j]) + f[j];
df_QuadraturePts[i] = (eta_QuadraturePts[i] - eta[j]) * (df[j+1] - df[j]) / (eta[j+1] - eta[j]) + df[j];
}
{
}
{
}
}
v_QuadraturePts[i] = nu * sqrt(U_inf / (2.0 * nu * (x_QuadraturePts[i] + x_0))) *
(y_QuadraturePts[i] * sqrt(U_inf / (2.0 * nu * (x_QuadraturePts[i] + x_0))) *
//! INFLOW SECTION: X = 0; Y > 0.
{
}
//! SINGULARITY POINT: X = 0; Y = 0.
{
}
}
}
//! --------------------------------------------------------------------------------------------
//! Inspection of the interpolation ------------------------------------------------------------
FILE *etaOriginal;
etaOriginal = fopen("eta.txt","w+");
for(i=0; i<numLines; i++)
{
fprintf(etaOriginal,"%f %f %f \n", eta[i], f[i], df[i]);
}
fclose(etaOriginal);
FILE *yQ;
yQ = fopen("yQ.txt","w+");
{
}
fclose(yQ);
//! --------------------------------------------------------------------------------------------
//! 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 ------------
Basis = Domain->GetExp(0)->GetBasis(0);
numModes = Basis->GetNumModes();
std::cout<< "Number of modes = " << numModes << std::endl;
//! Copying the ukGlobal vector (with the same pattern of m_phys) in m_phys
//! Initialisation of the ExpList Exp
Array<OneD, MultiRegions::ExpListSharedPtr> Exp(3);
Exp[0] = Exp2D_uk;
Exp[1] = Exp2D_vk;
if(stability_solver == "velocity_correction_scheme")
Exp[2] = Exp2D_pk;
//! Expansion coefficient extraction (necessary to write the .fld file)
Exp[0]->FwdTrans(Exp2D_uk->GetPhys(),Exp[0]->UpdateCoeffs());
Exp[1]->FwdTrans(Exp2D_vk->GetPhys(),Exp[1]->UpdateCoeffs());
if(stability_solver == "velocity_correction_scheme")
Exp[2]->FwdTrans(Exp2D_pk->GetPhys(),Exp[2]->UpdateCoeffs());
//! --------------------------------------------------------------------------------------------
//! Generation .FLD file with one field only (at the moment) -----------------------------------
//! Definition of the name of the .fld file
string FalknerSkan = "FalknerSkan.fld";
//! Definition of the number of the fields
if(stability_solver == "coupled_scheme")
nFields = 2;
if(stability_solver == "velocity_correction_scheme")
nFields = 3;
//! Definition of the Field
std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef = Exp[0]->GetFieldDefinitions();
std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
for(j = 0; j < nFields; ++j)
{
for(i = 0; i < FieldDef.size(); ++i)
{
if(j == 0)
{
FieldDef[i]->m_fields.push_back("u");
}
else if(j == 1)
{
FieldDef[i]->m_fields.push_back("v");
}
else if(j == 2 && stability_solver == "velocity_correction_scheme")
{
FieldDef[i]->m_fields.push_back("p");
}
Exp[j]->AppendFieldData(FieldDef[i], FieldData[i]);
}
}
LibUtilities::Write(FalknerSkan, FieldDef, FieldData);
std::cout<<"-------------------------------------------------------------------\n";
//! ----------------------------------------------------------------------------
return 0;
}