Nektar::RinglebFlowBC Class Reference

Wall boundary conditions for compressible flow problems. More...

#include <RinglebFlowBC.h>

Inheritance diagram for Nektar::RinglebFlowBC:

Static Public Member Functions
static CFSBndCondSharedPtr	create (const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const Array< OneD, Array< OneD, NekDouble > > &pTraceNormals, const int pSpaceDim, const int bcRegion, const int cnt)
	Creates an instance of this class. More...

Static Public Attributes
static std::string	className
	Name of the class. More...

Protected Member Functions
void	v_Apply (Array< OneD, Array< OneD, NekDouble > > &Fwd, Array< OneD, Array< OneD, NekDouble > > &physarray, const NekDouble &time) override
Protected Member Functions inherited from Nektar::CFSBndCond
	CFSBndCond (const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const Array< OneD, Array< OneD, NekDouble > > &pTraceNormals, const int pSpaceDim, const int bcRegion, const int cnt)
	Constructor. More...
virtual void	v_Apply (Array< OneD, Array< OneD, NekDouble > > &Fwd, Array< OneD, Array< OneD, NekDouble > > &physarray, const NekDouble &time)=0
virtual void	v_ApplyBwdWeight ()

Private Member Functions
	RinglebFlowBC (const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const Array< OneD, Array< OneD, NekDouble > > &pTraceNormals, const int pSpaceDim, const int bcRegion, const int cnt)
	~RinglebFlowBC (void) override

Private Attributes
int	m_expdim
bool	m_homo1D

Friends
class	MemoryManager< RinglebFlowBC >

Additional Inherited Members
Public Member Functions inherited from Nektar::CFSBndCond
virtual	~CFSBndCond ()
void	Apply (Array< OneD, Array< OneD, NekDouble > > &Fwd, Array< OneD, Array< OneD, NekDouble > > &physarray, const NekDouble &time=0)
	Apply the boundary condition. More...
void	ApplyBwdWeight ()
	Apply the Weight of boundary condition. More...
Protected Attributes inherited from Nektar::CFSBndCond
LibUtilities::SessionReaderSharedPtr	m_session
	Session reader. More...
Array< OneD, MultiRegions::ExpListSharedPtr >	m_fields
	Array of fields. More...
Array< OneD, Array< OneD, NekDouble > >	m_traceNormals
	Trace normals. More...
int	m_spacedim
	Space dimension. More...
VariableConverterSharedPtr	m_varConv
	Auxiliary object to convert variables. More...
NekDouble	m_diffusionAveWeight
	Weight for average calculation of diffusion term. More...
NekDouble	m_gamma
	Parameters of the flow. More...
NekDouble	m_rhoInf
NekDouble	m_pInf
NekDouble	m_pOut
Array< OneD, NekDouble >	m_velInf
int	m_bcRegion
	Id of the boundary region. More...
int	m_offset
	Offset. More...

Detailed Description

Wall boundary conditions for compressible flow problems.

Definition at line 46 of file RinglebFlowBC.h.

Constructor & Destructor Documentation

RinglebFlowBC()

Nektar::RinglebFlowBC::RinglebFlowBC ( const LibUtilities::SessionReaderSharedPtr &  pSession, const Array< OneD, MultiRegions::ExpListSharedPtr > &  pFields, const Array< OneD, Array< OneD, NekDouble > > &  pTraceNormals, const int  pSpaceDim, const int  bcRegion, const int  cnt  )

private

Definition at line 48 of file RinglebFlowBC.cpp.

    : CFSBndCond(pSession, pFields, pTraceNormals, pSpaceDim, bcRegion, cnt)
{
    m_expdim = pFields[0]->GetGraph()->GetMeshDimension();
 
    m_homo1D = false;
    if (m_session->DefinesSolverInfo("HOMOGENEOUS"))
    {
        std::string HomoStr = m_session->GetSolverInfo("HOMOGENEOUS");
        if ((HomoStr == "HOMOGENEOUS1D") || (HomoStr == "Homogeneous1D") ||
            (HomoStr == "1D") || (HomoStr == "Homo1D"))
        {
            m_homo1D = true;
        }
    }
}

References m_expdim, m_homo1D, and Nektar::CFSBndCond::m_session.

~RinglebFlowBC()

Nektar::RinglebFlowBC::~RinglebFlowBC ( void  )

inlineoverrideprivate

Definition at line 77 of file RinglebFlowBC.h.

{};

Member Function Documentation

create()

static CFSBndCondSharedPtr Nektar::RinglebFlowBC::create ( const LibUtilities::SessionReaderSharedPtr &  pSession, const Array< OneD, MultiRegions::ExpListSharedPtr > &  pFields, const Array< OneD, Array< OneD, NekDouble > > &  pTraceNormals, const int  pSpaceDim, const int  bcRegion, const int  cnt  )

inlinestatic

Creates an instance of this class.

Definition at line 52 of file RinglebFlowBC.h.

    {
        CFSBndCondSharedPtr p = MemoryManager<RinglebFlowBC>::AllocateSharedPtr(
            pSession, pFields, pTraceNormals, pSpaceDim, bcRegion, cnt);
        return p;
    }

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), and CellMLToNektar.cellml_metadata::p.

v_Apply()

void Nektar::RinglebFlowBC::v_Apply ( Array< OneD, Array< OneD, NekDouble > > &  Fwd, Array< OneD, Array< OneD, NekDouble > > &  physarray, const NekDouble &  time  )

overrideprotectedvirtual

Implements Nektar::CFSBndCond.

Definition at line 69 of file RinglebFlowBC.cpp.

{
    int nvariables = physarray.size();
 
    // For 3DHomogenoeus1D
    int n_planes = 1;
    if (m_expdim == 2 && m_homo1D)
    {
        int nPointsTot       = m_fields[0]->GetTotPoints();
        int nPointsTot_plane = m_fields[0]->GetPlane(0)->GetTotPoints();
        n_planes             = nPointsTot / nPointsTot_plane;
    }
 
    int id2, id2_plane, e_max;
 
    e_max = m_fields[0]->GetBndCondExpansions()[m_bcRegion]->GetExpSize();
 
    for (int e = 0; e < e_max; ++e)
    {
        int npoints = m_fields[0]
                          ->GetBndCondExpansions()[m_bcRegion]
                          ->GetExp(e)
                          ->GetTotPoints();
        int id1 =
            m_fields[0]->GetBndCondExpansions()[m_bcRegion]->GetPhys_Offset(e);
 
        // For 3DHomogenoeus1D
        if (m_expdim == 2 && m_homo1D)
        {
            int m_offset_plane = m_offset / n_planes;
            int e_plane;
            int e_max_plane     = e_max / n_planes;
            int nTracePts_plane = m_fields[0]->GetTrace()->GetNpoints();
 
            int planeID = floor((e + 0.5) / e_max_plane);
            e_plane     = e - e_max_plane * planeID;
 
            id2_plane = m_fields[0]->GetTrace()->GetPhys_Offset(
                m_fields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(
                    m_offset_plane + e_plane));
            id2 = id2_plane + planeID * nTracePts_plane;
        }
        else // For general case
        {
            id2 = m_fields[0]->GetTrace()->GetPhys_Offset(
                m_fields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(
                    m_offset + e));
        }
 
        Array<OneD, NekDouble> x0(npoints, 0.0);
        Array<OneD, NekDouble> x1(npoints, 0.0);
        Array<OneD, NekDouble> x2(npoints, 0.0);
 
        m_fields[0]->GetBndCondExpansions()[m_bcRegion]->GetExp(e)->GetCoords(
            x0, x1, x2);
 
        // Flow parameters
        NekDouble c, k, phi, r, J, VV, pp, sint, P, ss;
        NekDouble J11, J12, J21, J22, det;
        NekDouble Fx, Fy;
        NekDouble xi, yi;
        NekDouble dV;
        NekDouble dtheta;
        NekDouble par1;
        NekDouble theta     = M_PI / 4.0;
        NekDouble kExt      = 0.7;
        NekDouble V         = kExt * sin(theta);
        NekDouble toll      = 1.0e-8;
        NekDouble errV      = 1.0;
        NekDouble errTheta  = 1.0;
        NekDouble gamma     = m_gamma;
        NekDouble gamma_1_2 = (gamma - 1.0) / 2.0;
 
        // Loop on all the points of that edge
        for (int j = 0; j < npoints; j++)
        {
 
            while ((abs(errV) > toll) || (abs(errTheta) > toll))
            {
                VV   = V * V;
                sint = sin(theta);
                c    = sqrt(1.0 - gamma_1_2 * VV);
                k    = V / sint;
                phi  = 1.0 / k;
                pp   = phi * phi;
                J    = 1.0 / c + 1.0 / (3.0 * c * c * c) +
                    1.0 / (5.0 * c * c * c * c * c) -
                    0.5 * log((1.0 + c) / (1.0 - c));
 
                r    = pow(c, 1.0 / gamma_1_2);
                xi   = 1.0 / (2.0 * r) * (1.0 / VV - 2.0 * pp) + J / 2.0;
                yi   = phi / (r * V) * sqrt(1.0 - VV * pp);
                par1 = 25.0 - 5.0 * VV;
                ss   = sint * sint;
 
                Fx = xi - x0[j];
                Fy = yi - x1[j];
 
                J11 =
                    39062.5 / pow(par1, 3.5) * (1.0 / VV - 2.0 / VV * ss) * V +
                    1562.5 / pow(par1, 2.5) *
                        (-2.0 / (VV * V) + 4.0 / (VV * V) * ss) +
                    12.5 / pow(par1, 1.5) * V + 312.5 / pow(par1, 2.5) * V +
                    7812.5 / pow(par1, 3.5) * V -
                    0.25 *
                        (-1.0 / pow(par1, 0.5) * V /
                             (1.0 - 0.2 * pow(par1, 0.5)) -
                         (1.0 + 0.2 * pow(par1, 0.5)) /
                             pow((1.0 - 0.2 * pow(par1, 0.5)), 2.0) /
                             pow(par1, 0.5) * V) /
                        (1.0 + 0.2 * pow(par1, 0.5)) *
                        (1.0 - 0.2 * pow(par1, 0.5));
 
                J12 = -6250.0 / pow(par1, 2.5) / VV * sint * cos(theta);
                J21 =
                    -6250.0 / (VV * V) * sint / pow(par1, 2.5) *
                        pow((1.0 - ss), 0.5) +
                    78125.0 / V * sint / pow(par1, 3.5) * pow((1.0 - ss), 0.5);
 
                // the matrix is singular when theta = pi/2
                if (abs(x1[j]) < toll && abs(cos(theta)) < toll)
                {
                    J22 = -39062.5 / pow(par1, 3.5) / V +
                          3125 / pow(par1, 2.5) / (VV * V) +
                          12.5 / pow(par1, 1.5) * V +
                          312.5 / pow(par1, 2.5) * V +
                          7812.5 / pow(par1, 3.5) * V -
                          0.25 *
                              (-1.0 / pow(par1, 0.5) * V /
                                   (1.0 - 0.2 * pow(par1, 0.5)) -
                               (1.0 + 0.2 * pow(par1, 0.5)) /
                                   pow((1.0 - 0.2 * pow(par1, 0.5)), 2.0) /
                                   pow(par1, 0.5) * V) /
                              (1.0 + 0.2 * pow(par1, 0.5)) *
                              (1.0 - 0.2 * pow(par1, 0.5));
 
                    // dV = -dV/dx * Fx
                    dV     = -1.0 / J22 * Fx;
                    dtheta = 0.0;
                    theta  = M_PI / 2.0;
                }
                else
                {
                    J22 = 3125.0 / VV * cos(theta) / pow(par1, 2.5) *
                              pow((1.0 - ss), 0.5) -
                          3125.0 / VV * ss / pow(par1, 2.5) /
                              pow((1.0 - ss), 0.5) * cos(theta);
 
                    det = -1.0 / (J11 * J22 - J12 * J21);
 
                    // [dV dtheta]' = -[invJ]*[Fx Fy]'
                    dV     = det * (J22 * Fx - J12 * Fy);
                    dtheta = det * (-J21 * Fx + J11 * Fy);
                }
 
                V     = V + dV;
                theta = theta + dtheta;
 
                errV     = abs(dV);
                errTheta = abs(dtheta);
            }
 
            c                  = sqrt(1.0 - gamma_1_2 * VV);
            int kk             = id2 + j;
            NekDouble timeramp = 200.0;
            ;
            if (time < timeramp &&
                !(m_session->DefinesFunction("InitialConditions") &&
                  m_session->GetFunctionType("InitialConditions", 0) ==
                      LibUtilities::eFunctionTypeFile))
            {
                Fwd[0][kk] =
                    pow(c, 1.0 / gamma_1_2) * exp(-1.0 + time / timeramp);
 
                Fwd[1][kk] =
                    Fwd[0][kk] * V * cos(theta) * exp(-1 + time / timeramp);
 
                Fwd[2][kk] =
                    Fwd[0][kk] * V * sin(theta) * exp(-1 + time / timeramp);
            }
            else
            {
                Fwd[0][kk] = pow(c, 1.0 / gamma_1_2);
                Fwd[1][kk] = Fwd[0][kk] * V * cos(theta);
                Fwd[2][kk] = Fwd[0][kk] * V * sin(theta);
            }
 
            P          = (c * c) * Fwd[0][kk] / gamma;
            Fwd[3][kk] = P / (gamma - 1.0) +
                         0.5 * (Fwd[1][kk] * Fwd[1][kk] / Fwd[0][kk] +
                                Fwd[2][kk] * Fwd[2][kk] / Fwd[0][kk]);
 
            errV     = 1.0;
            errTheta = 1.0;
            theta    = M_PI / 4.0;
            V        = kExt * sin(theta);
        }
 
        for (int i = 0; i < nvariables; ++i)
        {
            Vmath::Vcopy(npoints, &Fwd[i][id2], 1,
                         &(m_fields[i]
                               ->GetBndCondExpansions()[m_bcRegion]
                               ->UpdatePhys())[id1],
                         1);
        }
    }
}

References tinysimd::abs(), Nektar::LibUtilities::eFunctionTypeFile, tinysimd::log(), Nektar::CFSBndCond::m_bcRegion, m_expdim, Nektar::CFSBndCond::m_fields, Nektar::CFSBndCond::m_gamma, m_homo1D, Nektar::CFSBndCond::m_offset, Nektar::CFSBndCond::m_session, Nektar::LibUtilities::P, tinysimd::sqrt(), and Vmath::Vcopy().

Friends And Related Function Documentation

MemoryManager< RinglebFlowBC >

friend class MemoryManager< RinglebFlowBC >

friend

Definition at line 1 of file RinglebFlowBC.h.

Member Data Documentation

className

std::string Nektar::RinglebFlowBC::className

static

Initial value:

=
    GetCFSBndCondFactory().RegisterCreatorFunction(
        "RinglebFlow", RinglebFlowBC::create,
        "Ringleb flow boundary condition.")

Name of the class.

Definition at line 64 of file RinglebFlowBC.h.

m_expdim

int Nektar::RinglebFlowBC::m_expdim

private

Definition at line 79 of file RinglebFlowBC.h.

Referenced by RinglebFlowBC(), and v_Apply().

m_homo1D

bool Nektar::RinglebFlowBC::m_homo1D

private

Definition at line 80 of file RinglebFlowBC.h.

Referenced by RinglebFlowBC(), and v_Apply().