Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator 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.

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

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::LibUtilities::SessionReader::CreateInstance(), Vmath::Vcopy(), and Nektar::LibUtilities::Write().

{
//! 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))
{
fprintf(stderr,"Usage: ./FldAddFalknerSkanBL sessionFile [SysSolnType]\n");exit(1);
}
//! Reading the session file
LibUtilities::SessionReaderSharedPtr vSession = LibUtilities::SessionReader::CreateInstance(argc, argv);
//! Loading the parameters to define the BL
NekDouble Re;
NekDouble L;
NekDouble U_inf;
NekDouble x;
NekDouble x_0;
NekDouble nu;
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);
}
vSession->LoadParameter("Re", Re, 1.0);
vSession->LoadParameter("L", L, 1.0);
vSession->LoadParameter("U_inf", U_inf, 1.0);
vSession->LoadParameter("x", x, 1.0);
vSession->LoadParameter("x_0", x_0, 1.0);
vSession->LoadParameter("NumberLines_txt", numLines, 1.0);
//! 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
SpatialDomains::MeshGraphSharedPtr graphShPt;
graphShPt = MemoryManager<SpatialDomains::MeshGraph2D>::AllocateSharedPtr(vSession);
//! Feed our spatial discretisation object
MultiRegions::ContField2DSharedPtr Domain;
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)
int nQuadraturePts;
nQuadraturePts = Domain->GetTotPoints();
std::cout << "Number of quadrature points = " << nQuadraturePts << std::endl;
//! Coordinates of the quadrature points
Array<OneD,NekDouble> x_QuadraturePts;
Array<OneD,NekDouble> y_QuadraturePts;
Array<OneD,NekDouble> z_QuadraturePts;
x_QuadraturePts = Array<OneD,NekDouble>(nQuadraturePts);
y_QuadraturePts = Array<OneD,NekDouble>(nQuadraturePts);
z_QuadraturePts = Array<OneD,NekDouble>(nQuadraturePts);
Domain->GetCoords(x_QuadraturePts,y_QuadraturePts,z_QuadraturePts);
//! 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;
NekDouble d;
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 ---------------------
Array<OneD,NekDouble> eta_QuadraturePts;
eta_QuadraturePts = Array<OneD,NekDouble>(nQuadraturePts);
Array<OneD,NekDouble> f_QuadraturePts;
f_QuadraturePts = Array<OneD,NekDouble>(nQuadraturePts);
Array<OneD,NekDouble> df_QuadraturePts;
df_QuadraturePts = Array<OneD,NekDouble>(nQuadraturePts);
Array<OneD,NekDouble> u_QuadraturePts;
u_QuadraturePts = Array<OneD,NekDouble>(nQuadraturePts);
Array<OneD,NekDouble> v_QuadraturePts;
v_QuadraturePts = Array<OneD,NekDouble>(nQuadraturePts);
Array<OneD,NekDouble> p_QuadraturePts;
p_QuadraturePts = Array<OneD,NekDouble>(nQuadraturePts);
//! --------------------------------------------------------------------------------------------
//! PARALLEL CASE ------------------------------------------------------------------------------
if(BL_type == "Parallel")
{
for(i=0; i<nQuadraturePts; i++)
{
eta_QuadraturePts[i] = y_QuadraturePts[i] * sqrt(U_inf / (2 * x * nu));
for(j=0; j<numLines-1; j++)
{
if(eta_QuadraturePts[i] >= eta[j] && eta_QuadraturePts[i] <= eta[j+1])
{
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];
}
else if(eta_QuadraturePts[i] == 1000000)
{
f_QuadraturePts[i] = f[numLines-1];
df_QuadraturePts[i] = df[numLines-1];
}
else if(eta_QuadraturePts[i] > eta[numLines-1])
{
f_QuadraturePts[i] = f[numLines-1] + df[numLines-1] * (eta_QuadraturePts[i] - eta[numLines-1]);
df_QuadraturePts[i] = df[numLines-1];
}
}
u_QuadraturePts[i] = U_inf * df_QuadraturePts[i];
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]);
p_QuadraturePts[i] = 0.0;
}
}
//! --------------------------------------------------------------------------------------------
//! NON-PARALLEL CASE --------------------------------------------------------------------------
if(BL_type == "nonParallel")
{
for(i=0; i<nQuadraturePts; i++)
{
eta_QuadraturePts[i] = y_QuadraturePts[i] * sqrt(U_inf / (2 * (x_QuadraturePts[i] + x_0) * nu));
if((x_QuadraturePts[i] + x_0) == 0)
{
eta_QuadraturePts[i] = 1000000;
}
for(j=0; j<numLines-1; j++)
{
if(eta_QuadraturePts[i] >= eta[j] && eta_QuadraturePts[i] <= eta[j+1])
{
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];
}
else if(eta_QuadraturePts[i] == 1000000)
{
f_QuadraturePts[i] = f[numLines-1];
df_QuadraturePts[i] = df[numLines-1];
}
else if(eta_QuadraturePts[i] > eta[numLines-1])
{
f_QuadraturePts[i] = f[numLines-1] + df[numLines-1] * (eta_QuadraturePts[i] - eta[numLines-1]);
df_QuadraturePts[i] = df[numLines-1];
}
}
u_QuadraturePts[i] = U_inf * df_QuadraturePts[i];
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))) *
df_QuadraturePts[i] - f_QuadraturePts[i]);
//! INFLOW SECTION: X = 0; Y > 0.
if((x_QuadraturePts[i] + x_0) == 0)
{
u_QuadraturePts[i] = U_inf;
v_QuadraturePts[i] = 0.0;
}
//! SINGULARITY POINT: X = 0; Y = 0.
if((x_QuadraturePts[i] + x_0) == 0 && y_QuadraturePts[i] == 0)
{
u_QuadraturePts[i] = 0.0;
v_QuadraturePts[i] = 0.0;
}
p_QuadraturePts[i] = 0.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+");
for(i=0; i<nQuadraturePts; i++)
{
fprintf(yQ,"%f %f %f %f %f %f %f\n", x_QuadraturePts[i], y_QuadraturePts[i], u_QuadraturePts[i],
v_QuadraturePts[i], eta_QuadraturePts[i], f_QuadraturePts[i], df_QuadraturePts[i]);
}
fclose(yQ);
//! --------------------------------------------------------------------------------------------
//! Definition of the 2D expansion using the mesh data specified on the session file -----------
MultiRegions::ExpList2DSharedPtr Exp2D_uk;
Exp2D_uk = MemoryManager<MultiRegions::ExpList2D>::AllocateSharedPtr(vSession,graphShPt);
MultiRegions::ExpList2DSharedPtr Exp2D_vk;
Exp2D_vk = MemoryManager<MultiRegions::ExpList2D>::AllocateSharedPtr(vSession,graphShPt);
MultiRegions::ExpList2DSharedPtr Exp2D_pk;
Exp2D_pk = MemoryManager<MultiRegions::ExpList2D>::AllocateSharedPtr(vSession,graphShPt);
//! --------------------------------------------------------------------------------------------
//! Filling the 2D expansion using a recursive algorithm based on the mesh ordering ------------
LibUtilities::BasisSharedPtr Basis;
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
Vmath::Vcopy(nQuadraturePts, u_QuadraturePts , 1, Exp2D_uk->UpdatePhys(), 1);
Vmath::Vcopy(nQuadraturePts, v_QuadraturePts, 1, Exp2D_vk->UpdatePhys(), 1);
Vmath::Vcopy(nQuadraturePts, p_QuadraturePts , 1, Exp2D_pk->UpdatePhys(), 1);
//! 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;
}
Generated on Sun Mar 15 2015 00:11:39 for Nektar++ by   doxygen 1.8.1.2