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().

{
    //! 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;
}