#include <GJPStabilisation.h>

Public Member Functions
	GJPStabilisation (ExpListSharedPtr field)

	~GJPStabilisation ()

void	Apply (const Array< OneD, NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const Array< OneD, NekDouble > &pUnorm=NullNekDouble1DArray, const NekDouble scale=1.0) const

Array< OneD, Array< OneD, NekDouble > > &	GetTraceNormals (void)

int	GetNumTracePts (void) const

bool	IsSemiImplicit () const

bool	IsExplicit () const

bool	IsImplicit () const

Array< OneD, NekDouble >	GetTraceWeightVarFactors (void)

Private Member Functions
void	ConstructLocalTraceJump (const int dir, const Array< OneD, const NekDouble > &in, Array< OneD, NekDouble > &store) const

void	ConstructLocalTraceJumpSI (const int dir, const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &store) const

void	IProductwrtStdDerivBaseOnTraceMat (int i, Array< OneD, NekDouble > &in, Array< OneD, NekDouble > &out) const

void	StdDerivOnTraceFromModes (int i, Array< OneD, NekDouble > &in, Array< OneD, NekDouble > &out) const

void	TraceJumpFromLocTraceNormDeriv (Array< OneD, NekDouble > &normderiv, Array< OneD, NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd) const

Private Attributes
unsigned	m_coordDim

unsigned	m_traceDim

unsigned	m_nLocTracePts

GJPFormulation	m_formulation = eGJPNoFormulation

Array< OneD, Array< OneD, NekDouble > >	m_traceNormals

MultiRegions::ExpListSharedPtr	m_dgfield
	DG expansion for projection evalaution along trace.

Array< OneD, Array< OneD, NekDouble > >	m_scalTrace
	Scale factor for phys values along trace involving the local normals and geometric factors.

Array< OneD, NekDouble >	m_locTraceWeights

std::vector< unsigned >	m_traceOffset
	phys offset in trace expannsion of each trace as we loop over elmts

std::vector< unsigned >	m_locTracePts0
	npoints in local trace expannsion in dir 0

std::vector< unsigned >	m_locTracePts1
	npoints in local trace expannsion in dir 1

std::vector< unsigned >	m_ntrace
	list of the number of traces over an element;

std::map< int, std::pair< LocalRegions::ExpansionSharedPtr, LocalRegions::ExpansionSharedPtr > >	m_interpTrace
	local trace and multiregion dg trace if different (i.e variable p and BC trace)

std::vector< bool >	m_traceFwd

std::vector< std::pair< int, Array< OneD, DNekMatSharedPtr > > >	m_StdDBaseOnTraceMat

std::vector< NekDouble >	m_locEdgeScale

Static Private Attributes
static std::string	GJPStabilisationLookupIds []

Detailed Description

Definition at line 54 of file GJPStabilisation.h.

Constructor & Destructor Documentation

◆ GJPStabilisation()

Nektar::MultiRegions::GJPStabilisation::GJPStabilisation ( ExpListSharedPtr field )

Definition at line 47 of file GJPStabilisation.cpp.

{
    LibUtilities::SessionReaderSharedPtr session = pField->GetSession();
 
    bool test;
 
    session->MatchSolverInfo("GJPStabilisation", "Explicit", test, false);
    if (test)
    {
        m_formulation = eGJPExplicit;
    }
 
    session->MatchSolverInfo("GJPStabilisation", "Implicit", test, false);
    if (test)
    {
        m_formulation = eGJPImplicit;
 
        ASSERTL0(session->MatchSolverInfo("GlobalSysSoln", "IterativeFull"),
                 "To use GJP Fully implicit stabilisation you must use a "
                 "Iterative Full solver");
    }
 
    session->MatchSolverInfo("GJPStabilisation", "SemiImplicit", test, false);
    if (test)
    {
        m_formulation = eGJPSemiImplicit;
    }
    ASSERTL0(m_formulation != eGJPNoFormulation,
             "Need a valid formualtion type for GradientJumpStabilisation: "
             "Explicit, Implicit, SemiImplicit");
 
    // Call GetTrace on the initialising field will set up
    // DG. Store a copy so that if we make a soft copy of
    // this class we can re-used this field for operators.
    pField->GetTrace();
    m_dgfield = pField;
 
    m_coordDim = m_dgfield->GetCoordim(0);
    m_traceDim = m_dgfield->GetShapeDimension() - 1;
 
    // set up trace normals but would be better if could use
    // equation system definition
    m_traceNormals = Array<OneD, Array<OneD, NekDouble>>(m_coordDim);
    for (int i = 0; i < m_coordDim; ++i)
    {
        m_traceNormals[i] =
            Array<OneD, NekDouble>(m_dgfield->GetTrace()->GetNpoints());
    }
    m_dgfield->GetTrace()->GetNormals(m_traceNormals);
 
    m_scalTrace = Array<OneD, Array<OneD, NekDouble>>(m_traceDim + 1);
 
    MultiRegions::ExpListSharedPtr dgtrace = m_dgfield->GetTrace();
 
    const std::shared_ptr<LocalRegions::ExpansionVector> exp =
        m_dgfield->GetExp();
 
    Array<OneD, Array<OneD, NekDouble>> dfactors[3];
    Array<OneD, NekDouble> e_tmp;
 
    m_nLocTracePts = 0;
    for (unsigned e = 0; e < (*exp).size(); ++e)
    {
        for (unsigned t = 0; t < (*exp)[e]->GetNtraces(); ++t)
        {
            m_nLocTracePts += (*exp)[e]->GetLocTraceExp(t)->GetTotPoints();
        }
    }
 
    m_scalTrace = Array<OneD, Array<OneD, NekDouble>>(m_traceDim + 1);
    for (int i = 0; i < m_traceDim + 1; ++i)
    {
        m_scalTrace[i] = Array<OneD, NekDouble>(m_nLocTracePts);
    }
 
    int cnt         = 0;
    int offset_phys = 0;
    Array<OneD, Array<OneD, Array<OneD, NekDouble>>> dbasis;
    Array<OneD, Array<OneD, Array<OneD, unsigned int>>> traceToCoeffMap;
 
    Array<OneD, int> sign, sign1;
    NekDouble h, p;
    std::map<unsigned, std::pair<NekDouble, unsigned>> hpscale;
 
    MultiRegions::DisContFieldSharedPtr dgfield =
        std::dynamic_pointer_cast<MultiRegions::DisContField>(m_dgfield);
 
    for (int e = 0; e < m_dgfield->GetExpSize(); ++e)
    {
        LocalRegions::ExpansionSharedPtr elmt = (*exp)[e];
 
        elmt->NormalTraceDerivFactors(dfactors[0], dfactors[1], dfactors[2]);
 
        for (int n = 0; n < elmt->GetNtraces(); ++n, ++cnt)
        {
            // collect offset for traces in dgtrace
            unsigned eid = dgfield->GetTraceElmtId(e, n);
            m_traceOffset.push_back(dgtrace->GetPhys_Offset(eid));
            LocalRegions::ExpansionSharedPtr LocTraceExp =
                elmt->GetLocTraceExp(n);
            unsigned LocTracepts = LocTraceExp->GetTotPoints();
            if (LocTracepts != dgtrace->GetExp(eid)->GetTotPoints())
            {
                m_interpTrace[cnt] =
                    std::make_pair(LocTraceExp, dgtrace->GetExp(eid));
            }
 
            m_locTracePts0.push_back(LocTraceExp->GetNumPoints(0));
            if (m_traceDim == 1)
            {
                m_locTracePts1.push_back(1);
            }
            else
            {
                m_locTracePts1.push_back(LocTraceExp->GetNumPoints(1));
            }
            m_traceFwd.push_back(dgfield->IsLeftAdjacentTrace(e, n));
 
            // gather h and p average informatoin
            elmt->TraceNormLen(n, h, p);
            unsigned edgeid = elmt->GetTraceExp(n)->GetGeom()->GetGlobalID();
            if (hpscale.count(edgeid))
            {
                auto hp         = hpscale[edgeid];
                hpscale[edgeid] = std::pair<NekDouble, unsigned>(
                    0.5 * (hp.first + h), (int)(0.5 * (hp.second + p)));
            }
            else
            {
                hpscale[edgeid] = std::pair<NekDouble, unsigned>(h, p);
            }
 
            int nptrace = dfactors[0][n].size();
            for (int i = 0; i < m_traceDim + 1; ++i)
            {
                Vmath::Smul(nptrace, 1.0, dfactors[i][n], 1,
                            e_tmp = m_scalTrace[i] + offset_phys, 1);
            }
            offset_phys += nptrace;
        }
        m_ntrace.push_back(elmt->GetNtraces());
    }
    m_nLocTracePts = offset_phys;
 
    m_locTraceWeights = Array<OneD, NekDouble>(m_nLocTracePts);
    // set up scale factors
    for (int e = 0; e < m_dgfield->GetExpSize(); ++e)
    {
        LocalRegions::ExpansionSharedPtr elmt = (*exp)[e];
        for (int n = 0; n < elmt->GetNtraces(); ++n, ++cnt)
        {
            unsigned edgeid = elmt->GetTraceExp(n)->GetGeom()->GetGlobalID();
 
            ASSERTL1(hpscale.count(edgeid), "Scale has not been defined");
            auto hp     = hpscale[edgeid];
            NekDouble h = hp.first;
            unsigned p  = hp.second;
            NekDouble jumpScal =
                (p == 1) ? 0.02 * h * h : 0.8 * pow(p + 1, -4.0) * h * h;
 
            m_locEdgeScale.push_back(jumpScal);
        }
    }
 
    m_locTraceWeights = Array<OneD, NekDouble>(m_nLocTracePts, 1.0);
 
    //  Generate array of quadrature and Jacobian for loc Trace
    unsigned offset = 0;
    cnt             = 0;
    for (int e = 0; e < m_dgfield->GetExpSize(); ++e)
    {
        for (unsigned j = 0; j < m_ntrace[e]; ++j)
        {
            unsigned eid =
                std::dynamic_pointer_cast<MultiRegions::DisContField>(m_dgfield)
                    ->GetTraceElmtId(e, j);
 
            if (m_interpTrace.count(cnt)) // manage interpolated case
            {
                m_interpTrace[cnt].first->MultiplyByQuadratureMetric(
                    m_locTraceWeights + offset,
                    e_tmp = m_locTraceWeights + offset);
            }
            else
            {
                dgtrace->GetExp(eid)->MultiplyByQuadratureMetric(
                    m_locTraceWeights + offset,
                    e_tmp = m_locTraceWeights + offset);
 
                // reverse data if necessary
                m_dgfield->GetExp(e)->ReOrientTracePhysVals(
                    m_dgfield->GetExp(e)->GetTraceOrient(j), e_tmp, e_tmp,
                    m_locTracePts0[cnt], m_locTracePts1[cnt], false);
            }
            offset += m_locTracePts0[cnt] * m_locTracePts1[cnt];
            cnt++;
        }
    }
    // Assemble list of Matrix Product
    Array<OneD, DNekMatSharedPtr> TraceMat;
 
    m_dgfield->GetExp(0)->StdDerivBaseOnTraceMat(TraceMat);
 
    int nelmt = 1;
    Array<OneD, const LibUtilities::BasisSharedPtr> base_sav =
        m_dgfield->GetExp(0)->GetBase();
 
    for (int n = 1; n < m_dgfield->GetExpSize(); ++n)
    {
        const Array<OneD, const LibUtilities::BasisSharedPtr> &base =
            m_dgfield->GetExp(n)->GetBase();
 
        // check to see if same element expansion as previous matrix and if
        // so can reused
        int i;
        for (i = 0; i < base.size(); ++i)
        {
            if (base[i] != base_sav[i])
            {
                break;
            }
        }
 
        if (i == base.size())
        {
            nelmt++;
        }
        else
        {
            // save previous block of data.
            m_StdDBaseOnTraceMat.push_back(
                std::pair<int, Array<OneD, DNekMatSharedPtr>>(nelmt, TraceMat));
 
            // start new block
            m_dgfield->GetExp(n)->StdDerivBaseOnTraceMat(TraceMat);
            nelmt    = 1;
            base_sav = m_dgfield->GetExp(n)->GetBase();
        }
    }
 
    // save latest block of data.
    m_StdDBaseOnTraceMat.push_back(
        std::pair<int, Array<OneD, DNekMatSharedPtr>>(nelmt, TraceMat));
}

References ASSERTL0, ASSERTL1, Nektar::MultiRegions::eGJPExplicit, Nektar::MultiRegions::eGJPImplicit, Nektar::MultiRegions::eGJPNoFormulation, Nektar::MultiRegions::eGJPSemiImplicit, m_coordDim, m_dgfield, m_formulation, m_interpTrace, m_locEdgeScale, m_locTracePts0, m_locTracePts1, m_locTraceWeights, m_nLocTracePts, m_ntrace, m_scalTrace, m_StdDBaseOnTraceMat, m_traceDim, m_traceFwd, m_traceNormals, m_traceOffset, sign, and Vmath::Smul().

◆ ~GJPStabilisation()

Nektar::MultiRegions::GJPStabilisation::~GJPStabilisation ( )

inline

Definition at line 59 of file GJPStabilisation.h.

59{};

Member Function Documentation

◆ Apply()

void Nektar::MultiRegions::GJPStabilisation::Apply	(	const Array< OneD, NekDouble > &	inarray,
		Array< OneD, NekDouble > &	outarray,
		const Array< OneD, NekDouble > &	pUnorm = `NullNekDouble1DArray`,
		const NekDouble	scale = `1.0`
	)		const

Definition at line 292 of file GJPStabilisation.cpp.

{
    int ncoeffs   = m_dgfield->GetNcoeffs();
    int nTracePts = m_dgfield->GetTrace()->GetTotPoints();
    LibUtilities::Timer timer, timer1;
    timer.Start();
    timer1.Start();
 
    Array<OneD, NekDouble> FilterCoeffs(ncoeffs), CoeffsTmp(ncoeffs);
    Array<OneD, NekDouble> Fwd(nTracePts, 0.0), Bwd(nTracePts, 0.0);
    Array<OneD, NekDouble> Store(m_nLocTracePts);
    Array<OneD, NekDouble> unorm;
 
    if (pUnorm == NullNekDouble1DArray)
    {
        unorm = Array<OneD, NekDouble>(nTracePts, 1.0);
    }
    else
    {
        unorm = pUnorm;
    }
 
    Array<OneD, NekDouble> dudn(m_nLocTracePts, 0.0);
 
    timer1.Stop();
    timer1.AccumulateRegion("GJP:Init", 10);
    timer1.Start();
 
    // Fwd Trans
    m_dgfield->FwdTransLocalElmt(inarray, FilterCoeffs);
 
    // Construct local derivaitves on trace
    for (int n = 0; n < m_coordDim; ++n)
    {
        StdDerivOnTraceFromModes(n, FilterCoeffs, Store);
        Vmath::Vvtvp(m_nLocTracePts, Store, 1, m_scalTrace[n], 1, dudn, 1, dudn,
                     1);
    }
 
    timer1.Stop();
    timer1.AccumulateRegion("GJP:FwdTrans + StdDerivFromModes", 10);
    timer1.Start();
 
    if (m_formulation == eGJPSemiImplicit)
    {
        // want to put Fwd vals on bwd trace and vice versa
        TraceJumpFromLocTraceNormDeriv(dudn, Bwd, Fwd);
    }
    else
    {
        TraceJumpFromLocTraceNormDeriv(dudn, Fwd, Bwd);
        Vmath::Vadd(nTracePts, Fwd, 1, Bwd, 1, Fwd, 1);
    }
    timer1.Stop();
    timer1.AccumulateRegion("GJP:TransJumpFromLocDeriv", 10);
    timer1.Start();
 
    if (m_formulation == eGJPSemiImplicit)
    {
        Vmath::Vmul(nTracePts, unorm, 1, Fwd, 1, Fwd, 1);
        Vmath::Vmul(nTracePts, unorm, 1, Bwd, 1, Bwd, 1);
        // Evaluate trace inner product with respect to derivative of the
        // basis
        ConstructLocalTraceJumpSI(0, Fwd, Bwd, Store);
        timer1.Stop();
        timer1.AccumulateRegion("GJP:Construct Trace ", 10);
        timer1.Start();
        IProductwrtStdDerivBaseOnTraceMat(0, Store, FilterCoeffs);
        timer1.Stop();
        timer1.AccumulateRegion("GJP:Deriv on Trace", 10);
 
        timer1.Start();
 
        for (int i = 0; i < m_traceDim; ++i)
        {
            ConstructLocalTraceJumpSI(i + 1, Fwd, Bwd, Store);
            timer1.Stop();
            timer1.AccumulateRegion("GJP:Construct Trace ", 10);
 
            timer1.Start();
            IProductwrtStdDerivBaseOnTraceMat(i + 1, Store, CoeffsTmp);
            timer1.Stop();
            timer1.AccumulateRegion("GJP:Deriv on Trace", 10);
 
            timer1.Start();
            Vmath::Vadd(ncoeffs, CoeffsTmp, 1, FilterCoeffs, 1, FilterCoeffs,
                        1);
        }
    }
    else
    {
        Vmath::Vmul(nTracePts, unorm, 1, Fwd, 1, Fwd, 1);
 
        // Evaluate trace inner product with respect to derivitive of the
        // basis
        ConstructLocalTraceJump(0, Fwd, Store);
        IProductwrtStdDerivBaseOnTraceMat(0, Store, FilterCoeffs);
 
        for (int i = 0; i < m_traceDim; ++i)
        {
            ConstructLocalTraceJump(i + 1, Fwd, Store);
            IProductwrtStdDerivBaseOnTraceMat(i + 1, Store, CoeffsTmp);
            Vmath::Vadd(ncoeffs, CoeffsTmp, 1, FilterCoeffs, 1, FilterCoeffs,
                        1);
        }
    }
 
    Vmath::Svtvp(ncoeffs, scale, FilterCoeffs, 1, outarray, 1, outarray, 1);
    timer.Stop();
    // Elapsed time
    timer.AccumulateRegion("GJP:Total", 10);
}

References Nektar::LibUtilities::Timer::AccumulateRegion(), ConstructLocalTraceJump(), ConstructLocalTraceJumpSI(), Nektar::MultiRegions::eGJPSemiImplicit, IProductwrtStdDerivBaseOnTraceMat(), m_coordDim, m_dgfield, m_formulation, m_nLocTracePts, m_scalTrace, m_traceDim, Nektar::NullNekDouble1DArray, Nektar::LibUtilities::Timer::Start(), StdDerivOnTraceFromModes(), Nektar::LibUtilities::Timer::Stop(), Vmath::Svtvp(), TraceJumpFromLocTraceNormDeriv(), Vmath::Vadd(), Vmath::Vmul(), and Vmath::Vvtvp().

◆ ConstructLocalTraceJump()

void Nektar::MultiRegions::GJPStabilisation::ConstructLocalTraceJump	(	const int	dir,
		const Array< OneD, const NekDouble > &	in,
		Array< OneD, NekDouble > &	store
	)		const

private

construct the gradient jump on local trace using input 'in' which is in trace orientation

Definition at line 412 of file GJPStabilisation.cpp.

{
    LibUtilities::Timer timer;
    unsigned offset = 0;
    unsigned cnt    = 0;
    Array<OneD, NekDouble> tmp;
 
    for (unsigned e = 0; e < m_dgfield->GetExpSize(); ++e)
    {
        for (unsigned j = 0; j < m_ntrace[e]; ++j)
        {
            unsigned locTracePts = m_locTracePts0[cnt] * m_locTracePts1[cnt];
 
            NekDouble jumpScal = m_locEdgeScale[cnt];
 
            if (m_interpTrace.count(cnt)) // variable p and BCs
            {
                auto it = m_interpTrace.find(cnt);
 
                StdRegions::Orientation orient =
                    m_dgfield->GetExp(e)->GetTraceOrient(j);
 
                // interpolate to new space
                it->second.first->PhysInterp(
                    it->second.second, in + m_traceOffset[cnt],
                    tmp = store + offset,
                    orient >= StdRegions::eDir1FwdDir2_Dir2FwdDir1);
 
                // reorientate
                m_dgfield->GetExp(e)->ReOrientTracePhysVals(
                    orient, tmp, tmp, m_locTracePts0[cnt], m_locTracePts1[cnt],
                    false);
            }
            else
            {
                tmp = store + offset;
                timer.Start();
                // Reverse data if necessary - using mapping from trace to local
                // trace
                m_dgfield->GetExp(e)->ReOrientTracePhysVals(
                    m_dgfield->GetExp(e)->GetTraceOrient(j),
                    in + m_traceOffset[cnt], tmp, m_locTracePts0[cnt],
                    m_locTracePts1[cnt], false);
                timer.Stop();
                timer.AccumulateRegion("GJPStab::Construct Trace:: Reorient",
                                       10);
            }
 
            // multiply by local factors in local trace space
            for (unsigned i = 0; i < locTracePts; ++i)
            {
                tmp[i] *= jumpScal * m_scalTrace[dir][offset + i] *
                          m_locTraceWeights[offset + i];
            }
 
            offset += locTracePts;
            cnt++;
        }
    }
}

References Nektar::LibUtilities::Timer::AccumulateRegion(), Nektar::StdRegions::eDir1FwdDir2_Dir2FwdDir1, m_dgfield, m_interpTrace, m_locEdgeScale, m_locTracePts0, m_locTracePts1, m_locTraceWeights, m_ntrace, m_scalTrace, m_traceOffset, Nektar::LibUtilities::Timer::Start(), and Nektar::LibUtilities::Timer::Stop().

Referenced by Apply().

◆ ConstructLocalTraceJumpSI()

void Nektar::MultiRegions::GJPStabilisation::ConstructLocalTraceJumpSI	(	const int	dir,
		const Array< OneD, const NekDouble > &	Fwd,
		const Array< OneD, const NekDouble > &	Bwd,
		Array< OneD, NekDouble > &	store
	)		const

private

Definition at line 475 of file GJPStabilisation.cpp.

{
    unsigned offset = 0;
    unsigned cnt    = 0;
    Array<OneD, NekDouble> tmp;
    LibUtilities::Timer timer;
 
    for (unsigned e = 0; e < m_dgfield->GetExpSize(); ++e)
    {
        for (unsigned j = 0; j < m_ntrace[e]; ++j)
        {
            unsigned locTracePts = m_locTracePts0[cnt] * m_locTracePts1[cnt];
 
            NekDouble jumpScal = m_locEdgeScale[cnt];
 
            if (m_traceFwd[cnt])
            {
                if (m_interpTrace.count(cnt)) // variable p and BCs
                {
                    auto it = m_interpTrace.find(cnt);
 
                    StdRegions::Orientation orient =
                        m_dgfield->GetExp(e)->GetTraceOrient(j);
 
                    // interpolate to new space
                    it->second.first->PhysInterp(
                        it->second.second, Fwd + m_traceOffset[cnt],
                        tmp = store + offset,
                        orient >= StdRegions::eDir1FwdDir2_Dir2FwdDir1);
 
                    // reorientate
                    m_dgfield->GetExp(e)->ReOrientTracePhysVals(
                        orient, tmp, tmp, m_locTracePts0[cnt],
                        m_locTracePts1[cnt], false);
                }
                else
                {
                    // Reverse data if necessary - using mapping from trace to
                    // local trace
                    timer.Start();
                    m_dgfield->GetExp(e)->ReOrientTracePhysVals(
                        m_dgfield->GetExp(e)->GetTraceOrient(j),
                        Fwd + m_traceOffset[cnt], tmp = store + offset,
                        m_locTracePts0[cnt], m_locTracePts1[cnt], false);
                    timer.Stop();
                    timer.AccumulateRegion(
                        "GJPStab::Construct Trace:: Reorient", 10);
                }
            }
            else
            {
                if (m_interpTrace.count(cnt)) // variable p and BCs
                {
                    auto it = m_interpTrace.find(cnt);
 
                    StdRegions::Orientation orient =
                        m_dgfield->GetExp(e)->GetTraceOrient(j);
 
                    // interpolate to new space
                    it->second.first->PhysInterp(
                        it->second.second, Fwd + m_traceOffset[cnt],
                        tmp = store + offset,
                        orient >= StdRegions::eDir1FwdDir2_Dir2FwdDir1);
 
                    // reorientate
                    m_dgfield->GetExp(e)->ReOrientTracePhysVals(
                        orient, tmp, tmp, m_locTracePts0[cnt],
                        m_locTracePts1[cnt], false);
                }
                else
                {
                    // Reverse data if necessary - using mapping from trace to
                    // local trace
                    timer.Start();
                    m_dgfield->GetExp(e)->ReOrientTracePhysVals(
                        m_dgfield->GetExp(e)->GetTraceOrient(j),
                        Bwd + m_traceOffset[cnt], tmp = store + offset,
                        m_locTracePts0[cnt], m_locTracePts1[cnt], false);
                    timer.Stop();
                    timer.AccumulateRegion(
                        "GJPStab::Construct Trace:: Reorient", 10);
                }
            }
            // multiply by local factors in local trace space
            for (unsigned i = 0; i < locTracePts; ++i)
            {
                tmp[i] *= jumpScal * m_scalTrace[dir][offset + i] *
                          m_locTraceWeights[offset + i];
            }
 
            offset += locTracePts;
            cnt++;
        }
    }
    ASSERTL1(offset <= 2 * m_dgfield->GetTrace()->GetTotPoints(),
             "Stroage is not large enough");
}

References Nektar::LibUtilities::Timer::AccumulateRegion(), ASSERTL1, Nektar::StdRegions::eDir1FwdDir2_Dir2FwdDir1, m_dgfield, m_interpTrace, m_locEdgeScale, m_locTracePts0, m_locTracePts1, m_locTraceWeights, m_ntrace, m_scalTrace, m_traceFwd, m_traceOffset, Nektar::LibUtilities::Timer::Start(), and Nektar::LibUtilities::Timer::Stop().

Referenced by Apply().

◆ GetNumTracePts()

int Nektar::MultiRegions::GJPStabilisation::GetNumTracePts ( void ) const

inline

Definition at line 71 of file GJPStabilisation.h.

    {
        return m_dgfield->GetTrace()->GetTotPoints();
    }

References m_dgfield.

◆ GetTraceNormals()

Array< OneD, Array< OneD, NekDouble > > & Nektar::MultiRegions::GJPStabilisation::GetTraceNormals ( void )

inline

Definition at line 66 of file GJPStabilisation.h.

    {
        return m_traceNormals;
    }

References m_traceNormals.

◆ GetTraceWeightVarFactors()

Array< OneD, NekDouble > Nektar::MultiRegions::GJPStabilisation::GetTraceWeightVarFactors ( void )

Definition at line 673 of file GJPStabilisation.cpp.

{
 
    Array<OneD, NekDouble> returnval(m_dgfield->GetNumElmts() * 6, 0.0);
 
    unsigned cnt = 0;
    for (unsigned e = 0; e < m_dgfield->GetExpSize(); ++e)
    {
        for (unsigned n = 0; n < m_dgfield->GetExp(e)->GetNtraces(); ++n)
        {
            returnval[6 * e + n] = m_locEdgeScale[cnt++];
        }
    }
    return returnval;
}

References m_dgfield, and m_locEdgeScale.

◆ IProductwrtStdDerivBaseOnTraceMat()

void Nektar::MultiRegions::GJPStabilisation::IProductwrtStdDerivBaseOnTraceMat	(	int	i,
		Array< OneD, NekDouble > &	in,
		Array< OneD, NekDouble > &	out
	)		const

private

Definition at line 579 of file GJPStabilisation.cpp.

{
    unsigned cnt  = 0;
    unsigned cnt1 = 0;
 
    for (auto &it : m_StdDBaseOnTraceMat)
    {
        unsigned modes    = it.second[dir]->GetRows();
        unsigned tracepts = it.second[dir]->GetColumns();
 
        Blas::Dgemm('N', 'N', modes, it.first, tracepts, 1.0,
                    &(it.second[dir]->GetPtr())[0], modes, &in[0] + cnt,
                    tracepts, 0.0, &out[0] + cnt1, modes);
 
        cnt += tracepts * it.first;
        cnt1 += modes * it.first;
    }
}

References Blas::Dgemm(), and m_StdDBaseOnTraceMat.

Referenced by Apply().

◆ IsExplicit()

bool Nektar::MultiRegions::GJPStabilisation::IsExplicit ( ) const

inline

Definition at line 81 of file GJPStabilisation.h.

    {
        return (m_formulation == eGJPExplicit);
    }

References Nektar::MultiRegions::eGJPExplicit, and m_formulation.

◆ IsImplicit()

bool Nektar::MultiRegions::GJPStabilisation::IsImplicit ( ) const

inline

Definition at line 86 of file GJPStabilisation.h.

    {
        return (m_formulation == eGJPImplicit);
    }

References Nektar::MultiRegions::eGJPImplicit, and m_formulation.

◆ IsSemiImplicit()

bool Nektar::MultiRegions::GJPStabilisation::IsSemiImplicit ( ) const

inline

Definition at line 76 of file GJPStabilisation.h.

    {
        return (m_formulation == eGJPSemiImplicit);
    }

References Nektar::MultiRegions::eGJPSemiImplicit, and m_formulation.

◆ StdDerivOnTraceFromModes()

void Nektar::MultiRegions::GJPStabilisation::StdDerivOnTraceFromModes	(	int	i,
		Array< OneD, NekDouble > &	in,
		Array< OneD, NekDouble > &	out
	)		const

private

Definition at line 653 of file GJPStabilisation.cpp.

{
    unsigned cnt  = 0;
    unsigned cnt1 = 0;
 
    for (auto &it : m_StdDBaseOnTraceMat)
    {
        unsigned modes    = it.second[dir]->GetRows();
        unsigned tracepts = it.second[dir]->GetColumns();
 
        Blas::Dgemm('T', 'N', tracepts, it.first, modes, 1.0,
                    &(it.second[dir]->GetPtr())[0], modes, &in[0] + cnt, modes,
                    0.0, &out[0] + cnt1, tracepts);
 
        cnt += modes * it.first;
        cnt1 += tracepts * it.first;
    }
}

References Blas::Dgemm(), and m_StdDBaseOnTraceMat.

Referenced by Apply().

◆ TraceJumpFromLocTraceNormDeriv()

void Nektar::MultiRegions::GJPStabilisation::TraceJumpFromLocTraceNormDeriv	(	Array< OneD, NekDouble > &	normderiv,
		Array< OneD, NekDouble > &	Fwd,
		Array< OneD, NekDouble > &	Bwd
	)		const

private

Definition at line 599 of file GJPStabilisation.cpp.

{
    unsigned offset = 0;
    unsigned cnt    = 0;
    Array<OneD, NekDouble> tmp, tmp1;
 
    for (unsigned e = 0; e < m_dgfield->GetExpSize(); ++e)
    {
        for (unsigned j = 0; j < m_ntrace[e]; ++j)
        {
            unsigned locTracePts = m_locTracePts0[cnt] * m_locTracePts1[cnt];
 
            tmp = normderiv + offset;
            ASSERTL0(offset < normderiv.size(), "Issue");
            // Reverse data if necessary
            m_dgfield->GetExp(e)->ReOrientTracePhysVals(
                m_dgfield->GetExp(e)->GetTraceOrient(j), tmp, tmp,
                m_locTracePts0[cnt], m_locTracePts1[cnt]);
 
            if (m_traceFwd[cnt])
            {
                tmp1 = Fwd + m_traceOffset[cnt];
            }
            else
            {
                tmp1 = Bwd + m_traceOffset[cnt];
            }
 
            if (m_interpTrace.count(cnt)) // variable p and BCs
            {
                auto it = m_interpTrace.find(cnt);
                it->second.second->PhysInterp(it->second.first, tmp, tmp1);
            }
            else
            {
                Vmath::Vcopy(locTracePts, tmp, 1, tmp1, 1);
            }
 
            offset += locTracePts;
            cnt++;
        }
    }
 
    // globally assemble Fwd and Bwd traces;
    m_dgfield->PeriodicBwdCopy(Fwd, Bwd);
    m_dgfield->FillBwdWithBoundCond(Fwd, Bwd, true);
    m_dgfield->GetTraceMap()->GetAssemblyCommDG()->PerformExchange(Fwd, Bwd);
}

References ASSERTL0, m_dgfield, m_interpTrace, m_locTracePts0, m_locTracePts1, m_ntrace, m_traceFwd, m_traceOffset, and Vmath::Vcopy().

Referenced by Apply().

Member Data Documentation

◆ GJPStabilisationLookupIds

std::string Nektar::MultiRegions::GJPStabilisation::GJPStabilisationLookupIds

staticprivate

Initial value:

= {
    LibUtilities::SessionReader::RegisterEnumValue(
        "GJPStabilisation", "Explicit", eExplicitGJPStabilisation),
    LibUtilities::SessionReader::RegisterEnumValue(
        "GJPStabilisation", "SemiImplicit", eSemiImplicitGJPStabilisation),
    LibUtilities::SessionReader::RegisterEnumValue(
        "GJPStabilisation", "Implicit", eFullImplicitGJPStabilisation),
}

Definition at line 98 of file GJPStabilisation.h.

◆ m_coordDim

unsigned Nektar::MultiRegions::GJPStabilisation::m_coordDim

private

Definition at line 94 of file GJPStabilisation.h.

Referenced by Apply(), and GJPStabilisation().

◆ m_dgfield

MultiRegions::ExpListSharedPtr Nektar::MultiRegions::GJPStabilisation::m_dgfield

private

DG expansion for projection evalaution along trace.

Definition at line 104 of file GJPStabilisation.h.

Referenced by Apply(), ConstructLocalTraceJump(), ConstructLocalTraceJumpSI(), GetNumTracePts(), GetTraceWeightVarFactors(), GJPStabilisation(), and TraceJumpFromLocTraceNormDeriv().

◆ m_formulation

GJPFormulation Nektar::MultiRegions::GJPStabilisation::m_formulation = eGJPNoFormulation

private

Definition at line 97 of file GJPStabilisation.h.

Referenced by Apply(), GJPStabilisation(), IsExplicit(), IsImplicit(), and IsSemiImplicit().

◆ m_interpTrace

std::map<int, std::pair<LocalRegions::ExpansionSharedPtr, LocalRegions::ExpansionSharedPtr> > Nektar::MultiRegions::GJPStabilisation::m_interpTrace

private

local trace and multiregion dg trace if different (i.e variable p and BC trace)

Definition at line 124 of file GJPStabilisation.h.

Referenced by ConstructLocalTraceJump(), ConstructLocalTraceJumpSI(), GJPStabilisation(), and TraceJumpFromLocTraceNormDeriv().

◆ m_locEdgeScale

std::vector<NekDouble> Nektar::MultiRegions::GJPStabilisation::m_locEdgeScale

private

Definition at line 131 of file GJPStabilisation.h.

Referenced by ConstructLocalTraceJump(), ConstructLocalTraceJumpSI(), GetTraceWeightVarFactors(), and GJPStabilisation().

◆ m_locTracePts0

std::vector<unsigned> Nektar::MultiRegions::GJPStabilisation::m_locTracePts0

private

npoints in local trace expannsion in dir 0

Definition at line 115 of file GJPStabilisation.h.

Referenced by ConstructLocalTraceJump(), ConstructLocalTraceJumpSI(), GJPStabilisation(), and TraceJumpFromLocTraceNormDeriv().

◆ m_locTracePts1

std::vector<unsigned> Nektar::MultiRegions::GJPStabilisation::m_locTracePts1

private

npoints in local trace expannsion in dir 1

Definition at line 117 of file GJPStabilisation.h.

Referenced by ConstructLocalTraceJump(), ConstructLocalTraceJumpSI(), GJPStabilisation(), and TraceJumpFromLocTraceNormDeriv().

◆ m_locTraceWeights

Array<OneD, NekDouble> Nektar::MultiRegions::GJPStabilisation::m_locTraceWeights

private

Definition at line 110 of file GJPStabilisation.h.

Referenced by ConstructLocalTraceJump(), ConstructLocalTraceJumpSI(), and GJPStabilisation().

◆ m_nLocTracePts

unsigned Nektar::MultiRegions::GJPStabilisation::m_nLocTracePts

private

Definition at line 96 of file GJPStabilisation.h.

Referenced by Apply(), and GJPStabilisation().

◆ m_ntrace

std::vector<unsigned> Nektar::MultiRegions::GJPStabilisation::m_ntrace

private

list of the number of traces over an element;

Definition at line 119 of file GJPStabilisation.h.

Referenced by ConstructLocalTraceJump(), ConstructLocalTraceJumpSI(), GJPStabilisation(), and TraceJumpFromLocTraceNormDeriv().

◆ m_scalTrace

Array<OneD, Array<OneD, NekDouble> > Nektar::MultiRegions::GJPStabilisation::m_scalTrace

private

Scale factor for phys values along trace involving the local normals and geometric factors.

Definition at line 108 of file GJPStabilisation.h.

Referenced by Apply(), ConstructLocalTraceJump(), ConstructLocalTraceJumpSI(), and GJPStabilisation().

◆ m_StdDBaseOnTraceMat

std::vector<std::pair<int, Array<OneD, DNekMatSharedPtr> > > Nektar::MultiRegions::GJPStabilisation::m_StdDBaseOnTraceMat

private

Definition at line 129 of file GJPStabilisation.h.

Referenced by GJPStabilisation(), IProductwrtStdDerivBaseOnTraceMat(), and StdDerivOnTraceFromModes().

◆ m_traceDim

unsigned Nektar::MultiRegions::GJPStabilisation::m_traceDim

private

Definition at line 95 of file GJPStabilisation.h.

Referenced by Apply(), and GJPStabilisation().

◆ m_traceFwd

std::vector<bool> Nektar::MultiRegions::GJPStabilisation::m_traceFwd

private

Definition at line 126 of file GJPStabilisation.h.

Referenced by ConstructLocalTraceJumpSI(), GJPStabilisation(), and TraceJumpFromLocTraceNormDeriv().

◆ m_traceNormals

Array<OneD, Array<OneD, NekDouble> > Nektar::MultiRegions::GJPStabilisation::m_traceNormals

private

Definition at line 101 of file GJPStabilisation.h.

Referenced by GetTraceNormals(), and GJPStabilisation().

◆ m_traceOffset

std::vector<unsigned> Nektar::MultiRegions::GJPStabilisation::m_traceOffset

private

phys offset in trace expannsion of each trace as we loop over elmts

Definition at line 113 of file GJPStabilisation.h.

Referenced by ConstructLocalTraceJump(), ConstructLocalTraceJumpSI(), GJPStabilisation(), and TraceJumpFromLocTraceNormDeriv().

Public Member Functions

Private Member Functions

Private Attributes

Static Private Attributes

Detailed Description

Constructor & Destructor Documentation

◆ GJPStabilisation()

◆ ~GJPStabilisation()

Member Function Documentation

◆ Apply()

◆ ConstructLocalTraceJump()

◆ ConstructLocalTraceJumpSI()

◆ GetNumTracePts()

◆ GetTraceNormals()

◆ GetTraceWeightVarFactors()

◆ IProductwrtStdDerivBaseOnTraceMat()

◆ IsExplicit()

◆ IsImplicit()

◆ IsSemiImplicit()

◆ StdDerivOnTraceFromModes()

◆ TraceJumpFromLocTraceNormDeriv()

Member Data Documentation

◆ GJPStabilisationLookupIds

◆ m_coordDim

◆ m_dgfield

◆ m_formulation

◆ m_interpTrace

◆ m_locEdgeScale

◆ m_locTracePts0

◆ m_locTracePts1

◆ m_locTraceWeights

◆ m_nLocTracePts

◆ m_ntrace

◆ m_scalTrace

◆ m_StdDBaseOnTraceMat

◆ m_traceDim

◆ m_traceFwd

◆ m_traceNormals

◆ m_traceOffset