Outflow characteristic boundary conditions for compressible flow problems. More...

#include <EnforceEntropyTotalEnthalpy.h>

Inheritance diagram for Nektar::EnforceEntropyTotalEnthalpy:

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 Array< OneD, Array< OneD, NekDouble > > &pGridVelocity, 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

void	GenerateRotationMatrices (const Array< OneD, const Array< OneD, NekDouble > > &normals)

void	FromToRotation (Array< OneD, const NekDouble > &from, Array< OneD, const NekDouble > &to, NekDouble *mat)

SOLVER_UTILS_EXPORT void	rotateToNormal (const Array< OneD, const Array< OneD, NekDouble > > &inarray, const Array< OneD, const Array< OneD, NekDouble > > &normals, const Array< OneD, const Array< OneD, NekDouble > > &vecLocs, Array< OneD, Array< OneD, NekDouble > > &outarray)

SOLVER_UTILS_EXPORT void	rotateFromNormal (const Array< OneD, const Array< OneD, NekDouble > > &inarray, const Array< OneD, const Array< OneD, NekDouble > > &normals, const Array< OneD, const Array< OneD, NekDouble > > &vecLocs, Array< OneD, Array< OneD, NekDouble > > &outarray)

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 Array< OneD, Array< OneD, NekDouble > > &pGridVelocity, 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 ()

Protected Attributes
int	m_npts

Array< OneD, NekDouble >	m_rhoBC

Array< OneD, Array< OneD, NekDouble > >	m_velBC

Array< OneD, NekDouble >	m_pBC

Array< OneD, NekDouble >	m_VnInf
	Reference normal velocity. More...

Array< OneD, int >	m_bndToTraceMap

Array< OneD, Array< OneD, NekDouble > >	m_bndPhys

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

Array< OneD, Array< OneD, NekDouble > >	m_gridVelocityTrace
	Grid Velocity. 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

NekDouble	m_angVel

int	m_bcRegion
	Id of the boundary region. More...

int	m_offset
	Offset. More...

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

	~EnforceEntropyTotalEnthalpy (void) override

Friends
class	MemoryManager< EnforceEntropyTotalEnthalpy >

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

Detailed Description

Outflow characteristic boundary conditions for compressible flow problems.

Definition at line 51 of file EnforceEntropyTotalEnthalpy.h.

Constructor & Destructor Documentation

◆ EnforceEntropyTotalEnthalpy()

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

private

Definition at line 51 of file EnforceEntropyTotalEnthalpy.cpp.

    : CFSBndCond(pSession, pFields, pTraceNormals, pGridVelocity, pSpaceDim,
                 bcRegion, cnt)
{
 
    const MultiRegions::ExpListSharedPtr bndexp =
        m_fields[0]->GetBndCondExpansions()[m_bcRegion];
 
    //-> Gather a list of index from trace to this boundary
    m_npts = bndexp->GetTotPoints();
 
    m_bndToTraceMap = Array<OneD, int>(m_npts, -1);
 
    const Array<OneD, const int> &traceBndMap = m_fields[0]->GetTraceBndMap();
 
    // Construct a map for the boundary to trace map for easy acess to
    // phys space points
    int cnt1 = 0;
    for (int e = 0; e < bndexp->GetNumElmts(); ++e)
    {
        int nTracePts = bndexp->GetExp(e)->GetTotPoints();
 
        int id =
            m_fields[0]->GetTrace()->GetPhys_Offset(traceBndMap[m_offset + e]);
 
        // Loop on the points of the m_bcRegion
        for (int i = 0; i < nTracePts; i++)
        {
            // the ith point in region e
            m_bndToTraceMap[cnt1++] = id + i;
        }
    }
 
    Array<OneD, Array<OneD, NekDouble>> BCvals(m_fields.size());
    m_bndPhys = Array<OneD, Array<OneD, NekDouble>>(m_fields.size());
 
    for (int i = 0; i < m_fields.size(); ++i)
    {
        m_bndPhys[i] =
            m_fields[i]->GetBndCondExpansions()[m_bcRegion]->UpdatePhys();
 
        BCvals[i] = Array<OneD, NekDouble>(m_npts);
        Vmath::Vcopy(m_npts, m_bndPhys[i], 1, BCvals[i], 1);
    }
 
    // Set up boudnary required BCs
    m_rhoBC = Array<OneD, NekDouble>(m_npts);
    Vmath::Vcopy(m_npts, BCvals[0], 1, m_rhoBC, 1);
 
    m_velBC = Array<OneD, Array<OneD, NekDouble>>(m_spacedim);
    // Evaluate velocity on boundary
    for (int i = 0; i < m_spacedim; ++i)
    {
        m_velBC[i] = Array<OneD, NekDouble>(m_npts);
        Vmath::Vcopy(m_npts, BCvals[i + 1], 1, m_velBC[i], 1);
        Vmath::Vdiv(m_npts, m_velBC[i], 1, m_rhoBC, 1, m_velBC[i], 1);
    }
    m_pBC = Array<OneD, NekDouble>(m_npts);
    m_varConv->GetPressure(BCvals, m_pBC);
 
    // Computing the normal velocity for characteristics coming
    // from outside the computational domain
    m_VnInf = Array<OneD, NekDouble>(m_npts, 0.0);
    for (int i = 0; i < m_spacedim; i++)
    {
        for (int j = 0; j < m_npts; ++j)
        {
            m_VnInf[j] += m_traceNormals[i][m_bndToTraceMap[j]] * m_velBC[i][j];
        }
    }
}

References Nektar::CFSBndCond::m_bcRegion, m_bndPhys, m_bndToTraceMap, Nektar::CFSBndCond::m_fields, m_npts, Nektar::CFSBndCond::m_offset, m_pBC, m_rhoBC, Nektar::CFSBndCond::m_spacedim, Nektar::CFSBndCond::m_traceNormals, Nektar::CFSBndCond::m_varConv, m_velBC, m_VnInf, Vmath::Vcopy(), and Vmath::Vdiv().

◆ ~EnforceEntropyTotalEnthalpy()

Nektar::EnforceEntropyTotalEnthalpy::~EnforceEntropyTotalEnthalpy ( void )

inlineoverrideprivate

Definition at line 121 of file EnforceEntropyTotalEnthalpy.h.

121{};

Member Function Documentation

◆ create()

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

inlinestatic

Creates an instance of this class.

Definition at line 57 of file EnforceEntropyTotalEnthalpy.h.

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

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

◆ FromToRotation()

void Nektar::EnforceEntropyTotalEnthalpy::FromToRotation	(	Array< OneD, const NekDouble > &	from,
		Array< OneD, const NekDouble > &	to,
		NekDouble *	mat
	)

protected

◆ GenerateRotationMatrices()

void Nektar::EnforceEntropyTotalEnthalpy::GenerateRotationMatrices ( const Array< OneD, const Array< OneD, NekDouble > > & normals )

protected

◆ rotateFromNormal()

SOLVER_UTILS_EXPORT void Nektar::EnforceEntropyTotalEnthalpy::rotateFromNormal	(	const Array< OneD, const Array< OneD, NekDouble > > &	inarray,
		const Array< OneD, const Array< OneD, NekDouble > > &	normals,
		const Array< OneD, const Array< OneD, NekDouble > > &	vecLocs,
		Array< OneD, Array< OneD, NekDouble > > &	outarray
	)

protected

◆ rotateToNormal()

SOLVER_UTILS_EXPORT void Nektar::EnforceEntropyTotalEnthalpy::rotateToNormal	(	const Array< OneD, const Array< OneD, NekDouble > > &	inarray,
		const Array< OneD, const Array< OneD, NekDouble > > &	normals,
		const Array< OneD, const Array< OneD, NekDouble > > &	vecLocs,
		Array< OneD, Array< OneD, NekDouble > > &	outarray
	)

protected

◆ v_Apply()

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

overrideprotectedvirtual

Implements Nektar::CFSBndCond.

Definition at line 128 of file EnforceEntropyTotalEnthalpy.cpp.

{
    int i, j;
    int nDimensions = m_spacedim;
 
    Array<OneD, Array<OneD, NekDouble>> FwdBnd(Fwd.size());
    Array<OneD, Array<OneD, NekDouble>> bndPhys(Fwd.size());
 
    // make a local copy of Fwd along boundary of interest
    for (i = 0; i < Fwd.size(); ++i)
    {
        FwdBnd[i] = Array<OneD, NekDouble>(m_npts);
        for (j = 0; j < m_npts; ++j)
        {
            FwdBnd[i][j] = Fwd[i][m_bndToTraceMap[j]];
        }
    }
 
    // Computing the normal velocity for characteristics coming
    // from inside the computational domain
    Array<OneD, NekDouble> Vn(m_npts, 0.0);
 
    for (i = 0; i < nDimensions; ++i)
    {
        for (j = 0; j < m_npts; ++j)
        {
            Vn[j] += m_traceNormals[i][m_bndToTraceMap[j]] * FwdBnd[i + 1][j];
        }
    }
    // divide by density.
    Vmath::Vdiv(m_npts, Vn, 1, FwdBnd[0], 1, Vn, 1);
 
    // Get speed of sound
    Array<OneD, NekDouble> pressure(m_npts);
    Array<OneD, NekDouble> soundSpeed(m_npts);
 
    m_varConv->GetPressure(FwdBnd, pressure);
    m_varConv->GetSoundSpeed(FwdBnd, soundSpeed);
 
    // Get Mach. Note: it is computed by Vn/c
    Array<OneD, NekDouble> Mach(m_npts, 0.0);
    Vmath::Vdiv(m_npts, Vn, 1, soundSpeed, 1, Mach, 1);
    Vmath::Vabs(m_npts, Mach, 1, Mach, 1);
 
    // Auxiliary variables
    Array<OneD, NekDouble> velBC(nDimensions, 0.0);
 
    // L represents properties outside boundary
    // R represents properties inside boundary (numerical state)
    NekDouble rhoL, uL, pL;
    NekDouble EBC, rR, cstar, pstar, rhostar, ustar, Sstar, vn;
 
    NekDouble gamMinOne        = m_gamma - 1.0;
    NekDouble twoOverGamMinOne = 2.0 / gamMinOne;
    NekDouble gamInv           = 1.0 / m_gamma;
 
    NekDouble tmp1, tmp2;
    // Loop on m_bcRegions
    for (int pnt = 0; pnt < m_npts; ++pnt)
    {
        // Impose inflow Riemann invariant
        if (Vn[pnt] <= 0.0)
        {
            // Subsonic flows
            if (Mach[pnt] < 1.00)
            {
                // right characteristic
                rR = -Vn[pnt] - sqrt(m_gamma * pressure[pnt] / FwdBnd[0][pnt]) *
                                    twoOverGamMinOne;
                vn = -m_VnInf[pnt]; // vn BC
 
                // fix total entropy and entropy to be input values
                // compute cstar and ustar using left-pointint characteristic
                // IR^-
                tmp1 = twoOverGamMinOne * rR;
                tmp1 = tmp1 * tmp1;
                tmp2 = rR * rR - vn * vn -
                       twoOverGamMinOne * m_gamma * m_pBC[pnt] / m_rhoBC[pnt];
                tmp2 =
                    (twoOverGamMinOne * twoOverGamMinOne + twoOverGamMinOne) *
                    tmp2;
                cstar   = -twoOverGamMinOne * rR + sqrt(tmp1 - tmp2);
                cstar   = cstar / (twoOverGamMinOne * twoOverGamMinOne +
                                 twoOverGamMinOne);
                ustar   = rR + twoOverGamMinOne * cstar;
                Sstar   = m_pBC[pnt] / pow(m_rhoBC[pnt], m_gamma);
                rhostar = pow(cstar * cstar / (m_gamma * Sstar),
                              0.5 * twoOverGamMinOne);
                pstar   = rhostar * cstar * cstar * gamInv;
 
                // add supplement equation that rhoL=rhostar
                // then pL=pstar, according to IL^0
                // and  uL=ustar, according to IL^+
                rhoL = rhostar;
                pL   = pstar;
                uL   = ustar;
            }
            else // Supersonic inflow
            {
                // all characteristics are from left so just impose
                // star state to left values
                // Note: m_vnInf is the negative of the normal velocity
                // across boundary
                rhoL = m_rhoBC[pnt];
                uL   = -m_VnInf[pnt];
                pL   = m_pBC[pnt];
            }
 
            // Boundary energy
            EBC = pL * twoOverGamMinOne * 0.5;
 
            // evaluate the different between the left state normal
            // velocity and that from the desired condition (note
            // m_VnInf is using an outwards normal definition.
            NekDouble VnDiff = uL + m_VnInf[pnt];
 
            // Boundary velocities & Kinite energy
            // Note: normals are negated since they point outwards in
            // the domain
 
            // Note: Can just use the BC values directly!!
            for (j = 0; j < nDimensions; ++j)
            {
                // Set velocity to the desired conditions modified to
                // take account of the normal state for Riemann
                // problem. (Negative accounts for outwards normal definition)
                // velBC[j] = m_velBC[j][pnt];
                velBC[j] = m_velBC[j][pnt] -
                           VnDiff * m_traceNormals[j][m_bndToTraceMap[pnt]];
 
                EBC += 0.5 * rhoL * velBC[j] * velBC[j];
            }
 
            //-------------------------------------------------------------------------
            // Impose Left hand Riemann Invariant boundary conditions
            m_bndPhys[0][pnt] = rhoL;
            for (j = 0; j < nDimensions; ++j)
            {
                m_bndPhys[j + 1][pnt] = rhoL * velBC[j];
            }
            m_bndPhys[nDimensions + 1][pnt] = EBC;
        }
        else // Outflow
        {
 
            // Note: Allowing the switch can cause worse convergence in this
            // type BC.
            //       So improve it later.
            if (Mach[pnt] < 1.00)
            {
                // subsonic outflow: fix pstar
                rR = -Vn[pnt] - sqrt(m_gamma * pressure[pnt] / FwdBnd[0][pnt]) *
                                    twoOverGamMinOne;
                // vn = -m_VnInf[pnt];
 
                pstar   = m_pBC[pnt];
                rhostar = FwdBnd[0][pnt] * pow((pstar / pressure[pnt]), gamInv);
                cstar   = sqrt(m_gamma * pstar / rhostar);
                ustar   = rR + cstar * twoOverGamMinOne;
 
                rhoL = rhostar;
                uL   = ustar;
                pL   = pstar;
            }
            else
            {
                // supersonic outflow
                // Just set to imposed state and let Riemann BC dictate values
                rhoL = m_rhoBC[pnt];
                uL   = -m_VnInf[pnt];
                pL   = m_pBC[pnt];
            }
 
            // Boundary energy
            EBC = pL * twoOverGamMinOne * 0.5;
 
            // Boundary velocities & Kinite energy
            // Note: normals are negated since they point outwards in
            // the domain
            for (j = 0; j < nDimensions; ++j)
            {
                velBC[j] = -1.0 * uL * m_traceNormals[j][m_bndToTraceMap[pnt]];
                EBC += 0.5 * rhoL * velBC[j] * velBC[j];
            }
 
            // Impose Left hand Riemann Invariant boundary conditions
            m_bndPhys[0][pnt] = rhoL;
            for (j = 0; j < nDimensions; ++j)
            {
                m_bndPhys[j + 1][pnt] = rhoL * velBC[j];
            }
            m_bndPhys[nDimensions + 1][pnt] = EBC;
        }
    }
}

References m_bndPhys, m_bndToTraceMap, Nektar::CFSBndCond::m_gamma, m_npts, m_pBC, m_rhoBC, Nektar::CFSBndCond::m_spacedim, Nektar::CFSBndCond::m_traceNormals, Nektar::CFSBndCond::m_varConv, m_velBC, m_VnInf, CG_Iterations::pressure, tinysimd::sqrt(), Vmath::Vabs(), and Vmath::Vdiv().

Friends And Related Function Documentation

◆ MemoryManager< EnforceEntropyTotalEnthalpy >

friend class MemoryManager< EnforceEntropyTotalEnthalpy >

friend

Definition at line 1 of file EnforceEntropyTotalEnthalpy.h.

Member Data Documentation

◆ className

std::string Nektar::EnforceEntropyTotalEnthalpy::className

static

Initial value:

=
    GetCFSBndCondFactory().RegisterCreatorFunction(
        "EnforceEntropyTotalEnthalpy", EnforceEntropyTotalEnthalpy::create,
        "Riemann invariant boundary condition, \
                        fixing H and S.")