Nektar::MultiRegions::PreconditionerLowEnergy Class Reference

#include <PreconditionerLowEnergy.h>

Inheritance diagram for Nektar::MultiRegions::PreconditionerLowEnergy:

Public Member Functions
	PreconditionerLowEnergy (const std::shared_ptr< GlobalLinSys > &plinsys, const AssemblyMapSharedPtr &pLocToGloMap)
	~PreconditionerLowEnergy () override
Public Member Functions inherited from Nektar::MultiRegions::Preconditioner
	Preconditioner (const std::shared_ptr< GlobalLinSys > &plinsys, const AssemblyMapSharedPtr &pLocToGloMap)
virtual	~Preconditioner ()
void	DoPreconditioner (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const bool &IsLocal=false)
void	DoAssembleLoc (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const bool &ZeroDir)
	Apply an assembly and scatter back to lcoal array.
void	DoPreconditionerWithNonVertOutput (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const Array< OneD, NekDouble > &pNonVertOutput, Array< OneD, NekDouble > &pVertForce=NullNekDouble1DArray)
void	DoTransformBasisToLowEnergy (Array< OneD, NekDouble > &pInOut)
void	DoTransformCoeffsFromLowEnergy (Array< OneD, NekDouble > &pInOut)
void	DoTransformCoeffsToLowEnergy (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)
void	DoTransformBasisFromLowEnergy (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)
void	BuildPreconditioner ()
void	InitObject ()
Array< OneD, NekDouble >	AssembleStaticCondGlobalDiagonals ()
	Performs global assembly of diagonal entries to global Schur complement matrix.
const DNekScalBlkMatSharedPtr &	GetBlockTransformedSchurCompl () const
const DNekScalBlkMatSharedPtr &	GetBlockCMatrix () const
const DNekScalBlkMatSharedPtr &	GetBlockInvDMatrix () const
const DNekScalBlkMatSharedPtr &	GetBlockSchurCompl () const
const DNekScalBlkMatSharedPtr &	GetBlockTransformationMatrix () const
const DNekScalBlkMatSharedPtr &	GetBlockTransposedTransformationMatrix () const
DNekScalMatSharedPtr	TransformedSchurCompl (int offset, int bndoffset, const std::shared_ptr< DNekScalMat > &loc_mat)

Static Public Member Functions
static PreconditionerSharedPtr	create (const std::shared_ptr< GlobalLinSys > &plinsys, const std::shared_ptr< AssemblyMap > &pLocToGloMap)
	Creates an instance of this class.

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

Protected Member Functions
void	v_InitObject () override
void	v_DoPreconditioner (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const bool &isLocal=false) override
void	v_BuildPreconditioner () override
	Construct the low energy preconditioner from \(\mathbf{S}_{2}\).
DNekScalMatSharedPtr	v_TransformedSchurCompl (int n, int offset, const std::shared_ptr< DNekScalMat > &loc_mat) override
	Set up the transformed block matrix system.
void	v_DoTransformBasisToLowEnergy (Array< OneD, NekDouble > &pInOut) override
	Transform the basis vector to low energy space.
void	v_DoTransformCoeffsFromLowEnergy (Array< OneD, NekDouble > &pInOut) override
	transform the solution coeffiicents from low energy back to the original coefficient space.
void	v_DoTransformBasisFromLowEnergy (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput) override
	Multiply by the block inverse transformation matrix This transforms the bassi from Low Energy to original basis.
void	v_DoTransformCoeffsToLowEnergy (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput) override
	Multiply by the block tranposed inverse transformation matrix (R^T)^{-1} which is equivlaent to transforming coefficients to LowEnergy space.
Protected Member Functions inherited from Nektar::MultiRegions::Preconditioner
virtual void	v_DoPreconditionerWithNonVertOutput (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const Array< OneD, NekDouble > &pNonVertOutput, Array< OneD, NekDouble > &pVertForce)
	Apply a preconditioner to the conjugate gradient method with an output for non-vertex degrees of freedom.

Protected Attributes
DNekBlkMatSharedPtr	m_BlkMat
DNekBlkMatSharedPtr	m_RBlk
DNekBlkMatSharedPtr	m_InvRBlk
SpatialDomains::EntityHolder	m_holder
Array< OneD, NekDouble >	m_variablePmask
bool	m_signChange
std::vector< std::pair< int, int > >	m_sameBlock
Protected Attributes inherited from Nektar::MultiRegions::Preconditioner
const std::weak_ptr< GlobalLinSys >	m_linsys
std::string	m_preconType
DNekMatSharedPtr	m_preconditioner
std::weak_ptr< AssemblyMap >	m_locToGloMap

Private Types
typedef std::map< LibUtilities::ShapeType, DNekScalMatSharedPtr >	ShapeToDNekMap
typedef std::map< LibUtilities::ShapeType, LocalRegions::Expansion3DSharedPtr >	ShapeToExpMap
typedef std::map< LibUtilities::ShapeType, Array< OneD, unsigned int > >	ShapeToIntArrayMap
typedef std::map< LibUtilities::ShapeType, Array< OneD, Array< OneD, unsigned int > > >	ShapeToIntArrayArrayMap

Private Member Functions
void	SetupBlockTransformationMatrix (void)
void	SetUpReferenceElements (ShapeToDNekMap &maxRmat, ShapeToExpMap &maxElmt, ShapeToIntArrayMap &vertMapMaxR, ShapeToIntArrayArrayMap &edgeMapMaxR)
	Loop expansion and determine different variants of the transformation matrix.
void	SetUpPyrMaxRMat (int nummodesmax, LocalRegions::PyrExpSharedPtr &PyrExp, ShapeToDNekMap &maxRmat, ShapeToIntArrayMap &vertMapMaxR, ShapeToIntArrayArrayMap &edgeMapMaxR, ShapeToIntArrayArrayMap &faceMapMaxR)
void	ReSetTetMaxRMat (int nummodesmax, LocalRegions::TetExpSharedPtr &TetExp, ShapeToDNekMap &maxRmat, ShapeToIntArrayMap &vertMapMaxR, ShapeToIntArrayArrayMap &edgeMapMaxR, ShapeToIntArrayArrayMap &faceMapMaxR)
void	ReSetPrismMaxRMat (int nummodesmax, LocalRegions::PrismExpSharedPtr &PirsmExp, ShapeToDNekMap &maxRmat, ShapeToIntArrayMap &vertMapMaxR, ShapeToIntArrayArrayMap &edgeMapMaxR, ShapeToIntArrayArrayMap &faceMapMaxR, bool UseTetOnly)
DNekMatSharedPtr	ExtractLocMat (LocalRegions::Expansion3DSharedPtr &locExp, DNekScalMatSharedPtr &maxRmat, LocalRegions::Expansion3DSharedPtr &expMax, Array< OneD, unsigned int > &vertMapMaxR, Array< OneD, Array< OneD, unsigned int > > &edgeMapMaxR)
void	CreateVariablePMask (void)
SpatialDomains::TetGeomUniquePtr	CreateRefTetGeom (void)
	Sets up the reference tretrahedral element needed to construct a low energy basis.
SpatialDomains::PyrGeomUniquePtr	CreateRefPyrGeom (void)
	Sets up the reference prismatic element needed to construct a low energy basis mapping arrays.
SpatialDomains::PrismGeomUniquePtr	CreateRefPrismGeom (void)
	Sets up the reference prismatic element needed to construct a low energy basis.
SpatialDomains::HexGeomUniquePtr	CreateRefHexGeom (void)
	Sets up the reference hexahedral element needed to construct a low energy basis.

Detailed Description

This class implements low energy preconditioning for the conjugate gradient matrix solver.

Definition at line 51 of file PreconditionerLowEnergy.h.

Member Typedef Documentation

ShapeToDNekMap

typedef std::map<LibUtilities::ShapeType, DNekScalMatSharedPtr> Nektar::MultiRegions::PreconditionerLowEnergy::ShapeToDNekMap

private

Definition at line 121 of file PreconditionerLowEnergy.h.

ShapeToExpMap

typedef std::map<LibUtilities::ShapeType, LocalRegions::Expansion3DSharedPtr> Nektar::MultiRegions::PreconditionerLowEnergy::ShapeToExpMap

private

Definition at line 124 of file PreconditionerLowEnergy.h.

ShapeToIntArrayArrayMap

typedef std::map<LibUtilities::ShapeType, Array<OneD, Array<OneD, unsigned int> > > Nektar::MultiRegions::PreconditionerLowEnergy::ShapeToIntArrayArrayMap

private

Definition at line 129 of file PreconditionerLowEnergy.h.

ShapeToIntArrayMap

typedef std::map<LibUtilities::ShapeType, Array<OneD, unsigned int> > Nektar::MultiRegions::PreconditionerLowEnergy::ShapeToIntArrayMap

private

Definition at line 126 of file PreconditionerLowEnergy.h.

Constructor & Destructor Documentation

PreconditionerLowEnergy()

Nektar::MultiRegions::PreconditionerLowEnergy::PreconditionerLowEnergy	(	const std::shared_ptr< GlobalLinSys > &	plinsys,
		const AssemblyMapSharedPtr &	pLocToGloMap
	)

Definition at line 65 of file PreconditionerLowEnergy.cpp.

    : Preconditioner(plinsys, pLocToGloMap)
{
}

~PreconditionerLowEnergy()

Nektar::MultiRegions::PreconditionerLowEnergy::~PreconditionerLowEnergy ( )

inlineoverride

Definition at line 74 of file PreconditionerLowEnergy.h.

    {
    }

Member Function Documentation

create()

static PreconditionerSharedPtr Nektar::MultiRegions::PreconditionerLowEnergy::create ( const std::shared_ptr< GlobalLinSys > &  plinsys, const std::shared_ptr< AssemblyMap > &  pLocToGloMap  )

inlinestatic

Creates an instance of this class.

Definition at line 55 of file PreconditionerLowEnergy.h.

    {
        PreconditionerSharedPtr p =
            MemoryManager<PreconditionerLowEnergy>::AllocateSharedPtr(
                plinsys, pLocToGloMap);
        p->InitObject();
        return p;
    }

References Nektar::MemoryManager< DataType >::AllocateSharedPtr().

CreateRefHexGeom()

SpatialDomains::HexGeomUniquePtr Nektar::MultiRegions::PreconditionerLowEnergy::CreateRefHexGeom ( void  )

private

Sets up the reference hexahedral element needed to construct a low energy basis.

Definition at line 1575 of file PreconditionerLowEnergy.cpp.

{
    ////////////////////////////////
    // Set up Hexahedron vertices //
    ////////////////////////////////
 
    const int three = 3;
 
    const int nVerts        = 8;
    const double point[][3] = {{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0},
                               {0, 0, 1}, {1, 0, 1}, {1, 1, 1}, {0, 1, 1}};
 
    // Populate the list of verts
    SpatialDomains::PointGeomUniquePtr verts[8];
    for (int i = 0; i < nVerts; ++i)
    {
        verts[i] = ObjPoolManager<SpatialDomains::PointGeom>::AllocateUniquePtr(
            three, i, point[i][0], point[i][1], point[i][2]);
    }
 
    /////////////////////////////
    // Set up Hexahedron Edges //
    /////////////////////////////
 
    // SegGeom (int id, const int coordim), EdgeComponent(id, coordim)
    const int nEdges                  = 12;
    const int vertexConnectivity[][2] = {{0, 1}, {1, 2}, {2, 3}, {0, 3},
                                         {0, 4}, {1, 5}, {2, 6}, {3, 7},
                                         {4, 5}, {5, 6}, {6, 7}, {4, 7}};
 
    // Populate the list of edges
    SpatialDomains::SegGeomUniquePtr edges[nEdges];
    for (int i = 0; i < nEdges; ++i)
    {
        std::array<SpatialDomains::PointGeom *, 2> vertsArray;
        for (int j = 0; j < 2; ++j)
        {
            vertsArray[j] = verts[vertexConnectivity[i][j]].get();
        }
        edges[i] = ObjPoolManager<SpatialDomains::SegGeom>::AllocateUniquePtr(
            i, three, vertsArray);
    }
    for (int i = 0; i < nVerts; ++i)
    {
        m_holder.m_pointVec.push_back(std::move(verts[i]));
    }
 
    /////////////////////////////
    // Set up Hexahedron faces //
    /////////////////////////////
 
    constexpr int nFaces            = 6;
    const int edgeConnectivity[][4] = {{0, 1, 2, 3},  {0, 5, 8, 4},
                                       {1, 6, 9, 5},  {2, 7, 10, 6},
                                       {3, 7, 11, 4}, {8, 9, 10, 11}};
 
    // Populate the list of faces
    std::array<SpatialDomains::QuadGeom *, nFaces> faces;
    for (int i = 0; i < nFaces; ++i)
    {
        std::array<SpatialDomains::SegGeom *, 4> edgeArray;
        for (int j = 0; j < 4; ++j)
        {
            edgeArray[j] = edges[edgeConnectivity[i][j]].get();
        }
        auto tmp = ObjPoolManager<SpatialDomains::QuadGeom>::AllocateUniquePtr(
            i, edgeArray);
        faces[i] = tmp.get();
        m_holder.m_quadVec.push_back(std::move(tmp));
    }
    for (int i = 0; i < nEdges; ++i)
    {
        m_holder.m_segVec.push_back(std::move(edges[i]));
    }
 
    SpatialDomains::HexGeomUniquePtr geom =
        ObjPoolManager<SpatialDomains::HexGeom>::AllocateUniquePtr(0, faces);
 
    return geom;
}

References Nektar::ObjPoolManager< DataType >::AllocateUniquePtr(), m_holder, Nektar::SpatialDomains::EntityHolder::m_pointVec, Nektar::SpatialDomains::EntityHolder::m_quadVec, and Nektar::SpatialDomains::EntityHolder::m_segVec.

Referenced by SetUpReferenceElements().

CreateRefPrismGeom()

SpatialDomains::PrismGeomUniquePtr Nektar::MultiRegions::PreconditionerLowEnergy::CreateRefPrismGeom ( void  )

private

Sets up the reference prismatic element needed to construct a low energy basis.

Definition at line 1282 of file PreconditionerLowEnergy.cpp.

{
    //////////////////////////
    // Set up Prism element //
    //////////////////////////
 
    const int three         = 3;
    const int nVerts        = 6;
    const double point[][3] = {
        {-1, -1, 0},
        {1, -1, 0},
        {1, 1, 0},
        {-1, 1, 0},
        {0, -1, sqrt(double(3))},
        {0, 1, sqrt(double(3))},
    };
 
    // std::shared_ptr<SpatialDomains::PointGeom> verts[6];
    SpatialDomains::PointGeomUniquePtr verts[6];
    for (int i = 0; i < nVerts; ++i)
    {
        verts[i] = ObjPoolManager<SpatialDomains::PointGeom>::AllocateUniquePtr(
            three, i, point[i][0], point[i][1], point[i][2]);
    }
    const int nEdges                  = 9;
    const int vertexConnectivity[][2] = {{0, 1}, {1, 2}, {3, 2}, {0, 3}, {0, 4},
                                         {1, 4}, {2, 5}, {3, 5}, {4, 5}};
 
    // Populate the list of edges
    SpatialDomains::SegGeomUniquePtr edges[nEdges];
    for (int i = 0; i < nEdges; ++i)
    {
        std::array<SpatialDomains::PointGeom *, 2> vertsArray;
        for (int j = 0; j < 2; ++j)
        {
            vertsArray[j] = verts[vertexConnectivity[i][j]].get();
        }
        edges[i] = ObjPoolManager<SpatialDomains::SegGeom>::AllocateUniquePtr(
            i, three, vertsArray);
    }
    for (int i = 0; i < nVerts; ++i)
    {
        m_holder.m_pointVec.push_back(std::move(verts[i]));
    }
 
    ////////////////////////
    // Set up Prism faces //
    ////////////////////////
 
    constexpr int nFaces = 5;
    // quad-edge connectivity base-face0, vertical-quadface2, vertical-quadface4
    const int quadEdgeConnectivity[][4] = {
        {0, 1, 2, 3}, {1, 6, 8, 5}, {3, 7, 8, 4}};
    // QuadId ordered as 0, 1, 2, otherwise return false
    const int quadId[] = {0, -1, 1, -1, 2};
 
    // triangle-edge connectivity side-triface-1, side triface-3
    const int triEdgeConnectivity[][3] = {{0, 5, 4}, {2, 6, 7}};
    // TriId ordered as 0, 1, otherwise return false
    const int triId[] = {-1, 0, -1, 1, -1};
 
    // Populate the list of faces
    std::array<SpatialDomains::Geometry2D *, nFaces> faces;
    for (int f = 0; f < nFaces; ++f)
    {
        if (f == 1 || f == 3)
        {
            int i = triId[f];
            std::array<SpatialDomains::SegGeom *, 3> edgeArray;
            for (int j = 0; j < 3; ++j)
            {
                edgeArray[j] = edges[triEdgeConnectivity[i][j]].get();
            }
            auto tmp =
                ObjPoolManager<SpatialDomains::TriGeom>::AllocateUniquePtr(
                    i, edgeArray);
            faces[f] = tmp.get();
            m_holder.m_triVec.push_back(std::move(tmp));
        }
        else
        {
            int i = quadId[f];
            std::array<SpatialDomains::SegGeom *, 4> edgeArray;
            for (int j = 0; j < 4; ++j)
            {
                edgeArray[j] = edges[quadEdgeConnectivity[i][j]].get();
            }
            auto tmp =
                ObjPoolManager<SpatialDomains::QuadGeom>::AllocateUniquePtr(
                    i, edgeArray);
            faces[f] = tmp.get();
            m_holder.m_quadVec.push_back(std::move(tmp));
        }
    }
    for (int i = 0; i < nEdges; ++i)
    {
        m_holder.m_segVec.push_back(std::move(edges[i]));
    }
 
    SpatialDomains::PrismGeomUniquePtr geom =
        ObjPoolManager<SpatialDomains::PrismGeom>::AllocateUniquePtr(0, faces);
 
    return geom;
}

References Nektar::ObjPoolManager< DataType >::AllocateUniquePtr(), m_holder, Nektar::SpatialDomains::EntityHolder::m_pointVec, Nektar::SpatialDomains::EntityHolder::m_quadVec, Nektar::SpatialDomains::EntityHolder::m_segVec, Nektar::SpatialDomains::EntityHolder::m_triVec, and tinysimd::sqrt().

Referenced by SetUpReferenceElements().

CreateRefPyrGeom()

SpatialDomains::PyrGeomUniquePtr Nektar::MultiRegions::PreconditionerLowEnergy::CreateRefPyrGeom ( void  )

private

Sets up the reference prismatic element needed to construct a low energy basis mapping arrays.

Definition at line 1391 of file PreconditionerLowEnergy.cpp.

{
    //////////////////////////
    // Set up Pyramid element //
    //////////////////////////
 
    const int nVerts        = 5;
    const double point[][3] = {{-1, -1, 0},
                               {1, -1, 0},
                               {1, 1, 0},
                               {-1, 1, 0},
                               {0, 0, sqrt(double(2))}};
 
    // boost::shared_ptr<SpatialDomains::PointGeom> verts[6];
    const int three = 3;
    SpatialDomains::PointGeomUniquePtr verts[5];
    for (int i = 0; i < nVerts; ++i)
    {
        verts[i] = ObjPoolManager<SpatialDomains::PointGeom>::AllocateUniquePtr(
            three, i, point[i][0], point[i][1], point[i][2]);
    }
    const int nEdges                  = 8;
    const int vertexConnectivity[][2] = {{0, 1}, {1, 2}, {2, 3}, {3, 0},
                                         {0, 4}, {1, 4}, {2, 4}, {3, 4}};
 
    // Populate the list of edges
    SpatialDomains::SegGeomUniquePtr edges[nEdges];
    for (int i = 0; i < nEdges; ++i)
    {
        std::array<SpatialDomains::PointGeom *, 2> vertsArray;
        for (int j = 0; j < 2; ++j)
        {
            vertsArray[j] = verts[vertexConnectivity[i][j]].get();
        }
        edges[i] = ObjPoolManager<SpatialDomains::SegGeom>::AllocateUniquePtr(
            i, three, vertsArray);
    }
 
    ////////////////////////
    // Set up Pyramid faces //
    ////////////////////////
 
    constexpr int nFaces = 5;
    // quad-edge connectivity base-face0,
    const int quadEdgeConnectivity[][4] = {{0, 1, 2, 3}};
 
    // triangle-edge connectivity side-triface-1, side triface-2
    const int triEdgeConnectivity[][3] = {
        {0, 5, 4}, {1, 6, 5}, {2, 7, 6}, {3, 4, 7}};
 
    // Populate the list of faces
    std::array<SpatialDomains::Geometry2D *, nFaces> faces;
    for (int f = 0; f < nFaces; ++f)
    {
        if (f == 0)
        {
            std::array<SpatialDomains::SegGeom *, 4> edgeArray;
            for (int j = 0; j < 4; ++j)
            {
                edgeArray[j] = edges[quadEdgeConnectivity[f][j]].get();
            }
 
            auto tmp =
                ObjPoolManager<SpatialDomains::QuadGeom>::AllocateUniquePtr(
                    f, edgeArray);
            faces[f] = tmp.get();
            m_holder.m_quadVec.push_back(std::move(tmp));
        }
        else
        {
            int i = f - 1;
            std::array<SpatialDomains::SegGeom *, 3> edgeArray;
            for (int j = 0; j < 3; ++j)
            {
                edgeArray[j] = edges[triEdgeConnectivity[i][j]].get();
            }
            auto tmp =
                ObjPoolManager<SpatialDomains::TriGeom>::AllocateUniquePtr(
                    f, edgeArray);
            faces[f] = tmp.get();
            m_holder.m_triVec.push_back(std::move(tmp));
        }
    }
    for (int i = 0; i < nEdges; ++i)
    {
        m_holder.m_segVec.push_back(std::move(edges[i]));
    }
 
    SpatialDomains::PyrGeomUniquePtr geom =
        ObjPoolManager<SpatialDomains::PyrGeom>::AllocateUniquePtr(0, faces);
 
    return geom;
}

References Nektar::ObjPoolManager< DataType >::AllocateUniquePtr(), m_holder, Nektar::SpatialDomains::EntityHolder::m_quadVec, Nektar::SpatialDomains::EntityHolder::m_segVec, Nektar::SpatialDomains::EntityHolder::m_triVec, and tinysimd::sqrt().

Referenced by SetUpReferenceElements().

CreateRefTetGeom()

SpatialDomains::TetGeomUniquePtr Nektar::MultiRegions::PreconditionerLowEnergy::CreateRefTetGeom ( void  )

private

Sets up the reference tretrahedral element needed to construct a low energy basis.

Definition at line 1489 of file PreconditionerLowEnergy.cpp.

{
    /////////////////////////////////
    // Set up Tetrahedron vertices //
    /////////////////////////////////
 
    int i, j;
    const int three         = 3;
    const int nVerts        = 4;
    const double point[][3] = {{-1, -1 / sqrt(double(3)), -1 / sqrt(double(6))},
                               {1, -1 / sqrt(double(3)), -1 / sqrt(double(6))},
                               {0, 2 / sqrt(double(3)), -1 / sqrt(double(6))},
                               {0, 0, 3 / sqrt(double(6))}};
 
    SpatialDomains::PointGeomUniquePtr verts[4];
    for (i = 0; i < nVerts; ++i)
    {
        verts[i] = ObjPoolManager<SpatialDomains::PointGeom>::AllocateUniquePtr(
            three, i, point[i][0], point[i][1], point[i][2]);
    }
 
    //////////////////////////////
    // Set up Tetrahedron Edges //
    //////////////////////////////
 
    // SegGeom (int id, const int coordim), EdgeComponent(id, coordim)
    const int nEdges                  = 6;
    const int vertexConnectivity[][2] = {{0, 1}, {1, 2}, {0, 2},
                                         {0, 3}, {1, 3}, {2, 3}};
 
    // Populate the list of edges
    SpatialDomains::SegGeomUniquePtr edges[nEdges];
    for (i = 0; i < nEdges; ++i)
    {
        std::array<SpatialDomains::PointGeom *, 2> vertsArray;
        for (j = 0; j < 2; ++j)
        {
            vertsArray[j] = verts[vertexConnectivity[i][j]].get();
        }
 
        edges[i] = ObjPoolManager<SpatialDomains::SegGeom>::AllocateUniquePtr(
            i, three, vertsArray);
    }
    for (i = 0; i < nVerts; ++i)
    {
        m_holder.m_pointVec.push_back(std::move(verts[i]));
    }
 
    //////////////////////////////
    // Set up Tetrahedron faces //
    //////////////////////////////
 
    constexpr int nFaces            = 4;
    const int edgeConnectivity[][3] = {
        {0, 1, 2}, {0, 4, 3}, {1, 5, 4}, {2, 5, 3}};
 
    // Populate the list of faces
    std::array<SpatialDomains::TriGeom *, nFaces> faces;
    for (i = 0; i < nFaces; ++i)
    {
        std::array<SpatialDomains::SegGeom *, 3> edgeArray;
        for (j = 0; j < 3; ++j)
        {
            edgeArray[j] = edges[edgeConnectivity[i][j]].get();
        }
 
        auto tmp = ObjPoolManager<SpatialDomains::TriGeom>::AllocateUniquePtr(
            i, edgeArray);
        faces[i] = tmp.get();
        m_holder.m_triVec.push_back(std::move(tmp));
    }
    for (i = 0; i < nEdges; ++i)
    {
        m_holder.m_segVec.push_back(std::move(edges[i]));
    }
 
    SpatialDomains::TetGeomUniquePtr geom =
        ObjPoolManager<SpatialDomains::TetGeom>::AllocateUniquePtr(0, faces);
 
    return geom;
}

References Nektar::ObjPoolManager< DataType >::AllocateUniquePtr(), m_holder, Nektar::SpatialDomains::EntityHolder::m_pointVec, Nektar::SpatialDomains::EntityHolder::m_segVec, Nektar::SpatialDomains::EntityHolder::m_triVec, and tinysimd::sqrt().

Referenced by SetUpReferenceElements().

CreateVariablePMask()

void Nektar::MultiRegions::PreconditionerLowEnergy::CreateVariablePMask ( void  )

private

Create the mask array

Definition at line 1255 of file PreconditionerLowEnergy.cpp.

{
    auto asmMap          = m_locToGloMap.lock();
    unsigned int nLocBnd = asmMap->GetNumLocalBndCoeffs();
 
    m_signChange = asmMap->GetSignChange();
 
    // Construct a mask array
    m_variablePmask = Array<OneD, NekDouble>(nLocBnd, 1.0);
 
    if (m_signChange)
    {
        unsigned int i;
        const Array<OneD, const NekDouble> &sign =
            asmMap->GetLocalToGlobalBndSign();
 
        for (i = 0; i < nLocBnd; ++i)
        {
            m_variablePmask[i] = fabs(sign[i]);
        }
    }
}

References Nektar::MultiRegions::Preconditioner::m_locToGloMap, m_signChange, m_variablePmask, and sign.

Referenced by v_InitObject().

ExtractLocMat()

DNekMatSharedPtr Nektar::MultiRegions::PreconditionerLowEnergy::ExtractLocMat ( LocalRegions::Expansion3DSharedPtr &  locExp, DNekScalMatSharedPtr &  maxRmat, LocalRegions::Expansion3DSharedPtr &  expMax, Array< OneD, unsigned int > &  vertMapMaxR, Array< OneD, Array< OneD, unsigned int > > &  edgeMapMaxR  )

private

Definition at line 2317 of file PreconditionerLowEnergy.cpp.

{
    NekDouble val;
    NekDouble zero = 0.0;
 
    int nRows = locExp->NumBndryCoeffs();
    DNekMatSharedPtr newmat =
        MemoryManager<DNekMat>::AllocateSharedPtr(nRows, nRows, zero, eFULL);
 
    Array<OneD, unsigned int> vlocmap;
    Array<OneD, Array<OneD, unsigned int>> elocmap;
    Array<OneD, Array<OneD, unsigned int>> flocmap;
 
    locExp->GetInverseBoundaryMaps(vlocmap, elocmap, flocmap);
 
    // fill diagonal
    for (int i = 0; i < nRows; ++i)
    {
        val = 1.0;
        newmat->SetValue(i, i, val);
    }
 
    int nverts = locExp->GetNverts();
    int nedges = locExp->GetNedges();
    int nfaces = locExp->GetNtraces();
 
    // fill vertex to edge coupling
    for (int e = 0; e < nedges; ++e)
    {
        int nEdgeInteriorCoeffs = locExp->GetEdgeNcoeffs(e) - 2;
 
        for (int v = 0; v < nverts; ++v)
        {
            for (int i = 0; i < nEdgeInteriorCoeffs; ++i)
            {
                val = (*maxRmat)(vmap[v], emap[e][i]);
                newmat->SetValue(vlocmap[v], elocmap[e][i], val);
            }
        }
    }
 
    for (int f = 0; f < nfaces; ++f)
    {
        // Get details to extrac this face from max reference matrix
        StdRegions::Orientation FwdOrient =
            StdRegions::eDir1FwdDir1_Dir2FwdDir2;
        int m0, m1; // Local face expansion orders.
 
        int nFaceInteriorCoeffs = locExp->GetTraceIntNcoeffs(f);
 
        locExp->GetTraceNumModes(f, m0, m1, FwdOrient);
 
        Array<OneD, unsigned int> fmapRmat =
            maxExp->GetTraceInverseBoundaryMap(f, FwdOrient, m0, m1);
 
        // fill in vertex to face coupling
        for (int v = 0; v < nverts; ++v)
        {
            for (int i = 0; i < nFaceInteriorCoeffs; ++i)
            {
                val = (*maxRmat)(vmap[v], fmapRmat[i]);
                newmat->SetValue(vlocmap[v], flocmap[f][i], val);
            }
        }
 
        // fill in edges to face coupling
        for (int e = 0; e < nedges; ++e)
        {
            int nEdgeInteriorCoeffs = locExp->GetEdgeNcoeffs(e) - 2;
 
            for (int j = 0; j < nEdgeInteriorCoeffs; ++j)
            {
                for (int i = 0; i < nFaceInteriorCoeffs; ++i)
                {
                    val = (*maxRmat)(emap[e][j], fmapRmat[i]);
                    newmat->SetValue(elocmap[e][j], flocmap[f][i], val);
                }
            }
        }
    }
 
    return newmat;
}

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::StdRegions::eDir1FwdDir1_Dir2FwdDir2, and Nektar::eFULL.

Referenced by SetupBlockTransformationMatrix().

ReSetPrismMaxRMat()

void Nektar::MultiRegions::PreconditionerLowEnergy::ReSetPrismMaxRMat ( int  nummodesmax, LocalRegions::PrismExpSharedPtr &  PirsmExp, ShapeToDNekMap &  maxRmat, ShapeToIntArrayMap &  vertMapMaxR, ShapeToIntArrayArrayMap &  edgeMapMaxR, ShapeToIntArrayArrayMap &  faceMapMaxR, bool  UseTetOnly  )

private

Definition at line 2101 of file PreconditionerLowEnergy.cpp.

{
    int nRows = PrismExp->NumBndryCoeffs();
    NekDouble val;
    NekDouble zero = 0.0;
    DNekMatSharedPtr newmat =
        MemoryManager<DNekMat>::AllocateSharedPtr(nRows, nRows, zero, eFULL);
 
    int nfacemodes;
 
    if (UseTetOnly)
    {
        // copy existing system
        for (int i = 0; i < nRows; ++i)
        {
            for (int j = 0; j < nRows; ++j)
            {
                val = (*maxRmat[ePrism])(i, j);
                newmat->SetValue(i, j, val);
            }
        }
 
        // Reset vertex to edge mapping from tet.
        const int VETetVert[][2]   = {{0, 0}, {1, 1}, {1, 1},
                                      {0, 0}, {3, 3}, {3, 3}};
        const int VETetEdge[][2]   = {{0, 3}, {0, 4}, {0, 4},
                                      {0, 3}, {3, 4}, {4, 3}};
        const int VEPrismEdge[][2] = {{0, 4}, {0, 5}, {2, 6},
                                      {2, 7}, {4, 5}, {6, 7}};
 
        // fill vertex to edge coupling
        for (int v = 0; v < 6; ++v)
        {
            for (int e = 0; e < 2; ++e)
            {
                for (int i = 0; i < nummodesmax - 2; ++i)
                {
                    // Note this is using wrong shape but gives
                    // answer that seems to have correct error!
                    val = (*maxRmat[eTetrahedron])(
                        vertMapMaxR[eTetrahedron][VETetVert[v][e]],
                        edgeMapMaxR[eTetrahedron][VETetEdge[v][e]][i]);
                    newmat->SetValue(vertMapMaxR[ePrism][v],
                                     edgeMapMaxR[ePrism][VEPrismEdge[v][e]][i],
                                     val);
                }
            }
        }
    }
    else
    {
        // set diagonal to 1
        for (int i = 0; i < nRows; ++i)
        {
            newmat->SetValue(i, i, 1.0);
        }
 
        // Set vertex to edge mapping from Hex.
 
        // The following lists specify the number of adjacent
        // edges to each vertex (nadj) then the Hex vert to
        // use for each prism ver in the vert-edge map (VEVert)
        // followed by the hex edge to use for each prism edge
        // in the vert-edge map (VEEdge)
        const int VEHexVert[][3]   = {{0, 0, 0}, {1, 1, 1}, {2, 2, 2},
                                      {3, 3, 3}, {4, 5, 5}, {6, 7, 7}};
        const int VEHexEdge[][3]   = {{0, 3, 4}, {0, 1, 5}, {1, 2, 6},
                                      {2, 3, 7}, {4, 5, 9}, {6, 7, 11}};
        const int VEPrismEdge[][3] = {{0, 3, 4}, {0, 1, 5}, {1, 2, 6},
                                      {2, 3, 7}, {4, 5, 8}, {6, 7, 8}};
 
        // fill vertex to edge coupling
        for (int v = 0; v < 6; ++v)
        {
            for (int e = 0; e < 3; ++e)
            {
                for (int i = 0; i < nummodesmax - 2; ++i)
                {
                    // Note this is using wrong shape but gives
                    // answer that seems to have correct error!
                    val = (*maxRmat[eHexahedron])(
                        vertMapMaxR[eHexahedron][VEHexVert[v][e]],
                        edgeMapMaxR[eHexahedron][VEHexEdge[v][e]][i]);
                    newmat->SetValue(vertMapMaxR[ePrism][v],
                                     edgeMapMaxR[ePrism][VEPrismEdge[v][e]][i],
                                     val);
                }
            }
        }
 
        // Setup vertex to face mapping from Hex
        const int VFHexVert[][2] = {{0, 0}, {1, 1}, {4, 5},
                                    {2, 2}, {3, 3}, {6, 7}};
        const int VFHexFace[][2] = {{0, 4}, {0, 2}, {4, 2},
                                    {0, 2}, {0, 4}, {2, 4}};
 
        const int VQFPrismVert[][2] = {{0, 0}, {1, 1}, {4, 4},
                                       {2, 2}, {3, 3}, {5, 5}};
        const int VQFPrismFace[][2] = {{0, 4}, {0, 2}, {4, 2},
                                       {0, 2}, {0, 4}, {2, 4}};
 
        nfacemodes = (nummodesmax - 2) * (nummodesmax - 2);
        // Replace two Quad faces  on every vertex
        for (int v = 0; v < 6; ++v)
        {
            for (int f = 0; f < 2; ++f)
            {
                for (int i = 0; i < nfacemodes; ++i)
                {
                    val = (*maxRmat[eHexahedron])(
                        vertMapMaxR[eHexahedron][VFHexVert[v][f]],
                        faceMapMaxR[eHexahedron][VFHexFace[v][f]][i]);
                    newmat->SetValue(vertMapMaxR[ePrism][VQFPrismVert[v][f]],
                                     faceMapMaxR[ePrism][VQFPrismFace[v][f]][i],
                                     val);
                }
            }
        }
 
        // Mapping of Hex Edge-Face mappings to Prism Edge-Face Mappings
        const int nadjface[]     = {1, 2, 1, 2, 1, 1, 1, 1, 2};
        const int EFHexEdge[][2] = {{0, -1}, {1, 1},  {2, -1}, {3, 3}, {4, -1},
                                    {5, -1}, {6, -1}, {7, -1}, {9, 11}};
        const int EFHexFace[][2] = {{0, -1}, {0, 2},  {0, -1}, {0, 4}, {4, -1},
                                    {2, -1}, {2, -1}, {4, -1}, {2, 4}};
        const int EQFPrismEdge[][2] = {{0, -1}, {1, 1},  {2, -1},
                                       {3, 3},  {4, -1}, {5, -1},
                                       {6, -1}, {7, -1}, {8, 8}};
        const int EQFPrismFace[][2] = {{0, -1}, {0, 2},  {0, -1},
                                       {0, 4},  {4, -1}, {2, -1},
                                       {2, -1}, {4, -1}, {2, 4}};
 
        // all base edges are coupled to face 0
        nfacemodes = (nummodesmax - 2) * (nummodesmax - 2);
        for (int e = 0; e < 9; ++e)
        {
            for (int f = 0; f < nadjface[e]; ++f)
            {
                for (int i = 0; i < nummodesmax - 2; ++i)
                {
                    for (int j = 0; j < nfacemodes; ++j)
                    {
                        int edgemapid =
                            edgeMapMaxR[eHexahedron][EFHexEdge[e][f]][i];
                        int facemapid =
                            faceMapMaxR[eHexahedron][EFHexFace[e][f]][j];
 
                        val = (*maxRmat[eHexahedron])(edgemapid, facemapid);
 
                        int edgemapid1 =
                            edgeMapMaxR[ePrism][EQFPrismEdge[e][f]][i];
                        int facemapid1 =
                            faceMapMaxR[ePrism][EQFPrismFace[e][f]][j];
                        newmat->SetValue(edgemapid1, facemapid1, val);
                    }
                }
            }
        }
    }
 
    const int VFTetVert[]    = {0, 1, 3, 1, 0, 3};
    const int VFTetFace[]    = {1, 1, 1, 1, 1, 1};
    const int VTFPrismVert[] = {0, 1, 4, 2, 3, 5};
    const int VTFPrismFace[] = {1, 1, 1, 3, 3, 3};
 
    //  Handle all triangular faces from tetrahedron
    nfacemodes = (nummodesmax - 3) * (nummodesmax - 2) / 2;
    for (int v = 0; v < 6; ++v)
    {
        for (int i = 0; i < nfacemodes; ++i)
        {
            val = (*maxRmat[eTetrahedron])(
                vertMapMaxR[eTetrahedron][VFTetVert[v]],
                faceMapMaxR[eTetrahedron][VFTetFace[v]][i]);
 
            newmat->SetValue(vertMapMaxR[ePrism][VTFPrismVert[v]],
                             faceMapMaxR[ePrism][VTFPrismFace[v]][i], val);
        }
    }
 
    // Mapping of Tet Edge-Face mappings to Prism Edge-Face Mappings
    const int EFTetEdge[]    = {0, 3, 4, 0, 4, 3};
    const int EFTetFace[]    = {1, 1, 1, 1, 1, 1};
    const int ETFPrismEdge[] = {0, 4, 5, 2, 6, 7};
    const int ETFPrismFace[] = {1, 1, 1, 3, 3, 3};
 
    // handle all edge to triangular faces from tetrahedron
    // (only 6 this time)
    nfacemodes = (nummodesmax - 3) * (nummodesmax - 2) / 2;
    for (int e = 0; e < 6; ++e)
    {
        for (int i = 0; i < nummodesmax - 2; ++i)
        {
            for (int j = 0; j < nfacemodes; ++j)
            {
                int edgemapid = edgeMapMaxR[eTetrahedron][EFTetEdge[e]][i];
                int facemapid = faceMapMaxR[eTetrahedron][EFTetFace[e]][j];
                val           = (*maxRmat[eTetrahedron])(edgemapid, facemapid);
 
                newmat->SetValue(edgeMapMaxR[ePrism][ETFPrismEdge[e]][i],
                                 faceMapMaxR[ePrism][ETFPrismFace[e]][j], val);
            }
        }
    }
 
    DNekScalMatSharedPtr PrismR =
        MemoryManager<DNekScalMat>::AllocateSharedPtr(1.0, newmat);
    maxRmat[ePrism] = PrismR;
}

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::eFULL, Nektar::LibUtilities::eHexahedron, Nektar::LibUtilities::ePrism, and Nektar::LibUtilities::eTetrahedron.

Referenced by SetUpReferenceElements().

ReSetTetMaxRMat()

void Nektar::MultiRegions::PreconditionerLowEnergy::ReSetTetMaxRMat ( int  nummodesmax, LocalRegions::TetExpSharedPtr &  TetExp, ShapeToDNekMap &  maxRmat, ShapeToIntArrayMap &  vertMapMaxR, ShapeToIntArrayArrayMap &  edgeMapMaxR, ShapeToIntArrayArrayMap &  faceMapMaxR  )

private

Definition at line 2043 of file PreconditionerLowEnergy.cpp.

{
    int nRows = TetExp->NumBndryCoeffs();
    NekDouble val;
    NekDouble zero = 0.0;
    DNekMatSharedPtr newmat =
        MemoryManager<DNekMat>::AllocateSharedPtr(nRows, nRows, zero, eFULL);
 
    // copy existing system
    for (int i = 0; i < nRows; ++i)
    {
        for (int j = 0; j < nRows; ++j)
        {
            val = (*maxRmat[eTetrahedron])(i, j);
            newmat->SetValue(i, j, val);
        }
    }
 
    // The following lists specify the number of adjacent
    // edges to each vertex (nadj) then the Hex vert to
    // use for each pyramid ver in the vert-edge map (VEVert)
    // followed by the hex edge to use for each Tet edge
    // in the vert-edge map (VEEdge)
    const int VEHexVert[][4] = {{0, 0, 0}, {1, 1, 1}, {2, 2, 2}, {4, 5, 6}};
    const int VEHexEdge[][4] = {{0, 3, 4}, {0, 1, 5}, {1, 2, 6}, {4, 5, 6}};
    const int VETetEdge[][4] = {{0, 2, 3}, {0, 1, 4}, {1, 2, 5}, {3, 4, 5}};
 
    // fill vertex to edge coupling
    for (int v = 0; v < 4; ++v)
    {
        for (int e = 0; e < 3; ++e)
        {
            for (int i = 0; i < nummodesmax - 2; ++i)
            {
                // Note this is using wrong shape but gives
                // answer that seems to have correct error!
                val = (*maxRmat[eHexahedron])(
                    vertMapMaxR[eHexahedron][VEHexVert[v][e]],
                    edgeMapMaxR[eHexahedron][VEHexEdge[v][e]][i]);
                newmat->SetValue(vertMapMaxR[eTetrahedron][v],
                                 edgeMapMaxR[eTetrahedron][VETetEdge[v][e]][i],
                                 val);
            }
        }
    }
 
    DNekScalMatSharedPtr TetR =
        MemoryManager<DNekScalMat>::AllocateSharedPtr(1.0, newmat);
 
    maxRmat[eTetrahedron] = TetR;
}

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::eFULL, Nektar::LibUtilities::eHexahedron, and Nektar::LibUtilities::eTetrahedron.

Referenced by SetUpReferenceElements().

SetupBlockTransformationMatrix()

void Nektar::MultiRegions::PreconditionerLowEnergy::SetupBlockTransformationMatrix ( void  )

private

Set a block transformation matrices for each element type. These are needed in routines that transform the schur complement matrix to and from the low energy basis.

Definition at line 939 of file PreconditionerLowEnergy.cpp.

{
    std::shared_ptr<MultiRegions::ExpList> expList =
        ((m_linsys.lock())->GetLocMat()).lock();
    LocalRegions::Expansion3DSharedPtr locExp;
    LocalRegions::Expansion3DSharedPtr locExpSav;
    map<int, int> EdgeSize;
 
    int n;
 
    std::map<ShapeType, DNekScalMatSharedPtr> maxRmat;
    map<ShapeType, LocalRegions::Expansion3DSharedPtr> maxElmt;
    map<ShapeType, Array<OneD, unsigned int>> vertMapMaxR;
    map<ShapeType, Array<OneD, Array<OneD, unsigned int>>> edgeMapMaxR;
 
    // Sets up reference element and builds transformation matrix for
    // maximum polynomial order meshes
    SetUpReferenceElements(maxRmat, maxElmt, vertMapMaxR, edgeMapMaxR);
 
    const Array<OneD, const unsigned int> &nbdry_size =
        m_locToGloMap.lock()->GetNumLocalBndCoeffsPerPatch();
 
    int n_exp = expList->GetNumElmts();
 
    MatrixStorage blkmatStorage = eDIAGONAL;
 
    // Variants of R matrices required for low energy preconditioning
    m_RBlk = MemoryManager<DNekBlkMat>::AllocateSharedPtr(
        nbdry_size, nbdry_size, blkmatStorage);
    m_InvRBlk = MemoryManager<DNekBlkMat>::AllocateSharedPtr(
        nbdry_size, nbdry_size, blkmatStorage);
 
    DNekMatSharedPtr rmat, invrmat;
 
    int offset = 0;
 
    // Set up transformation matrices whilst checking to see if
    // consecutive matrices are the same and if so reuse the
    // matrices and store how many consecutive offsets there
    // are
    for (n = 0; n < n_exp; ++n)
    {
        locExp = expList->GetExp(n)->as<LocalRegions::Expansion3D>();
 
        ShapeType eltype = locExp->DetShapeType();
 
        int nbndcoeffs = locExp->NumBndryCoeffs();
 
        if (m_sameBlock.size() == 0)
        {
            rmat = ExtractLocMat(locExp, maxRmat[eltype], maxElmt[eltype],
                                 vertMapMaxR[eltype], edgeMapMaxR[eltype]);
            // Block R matrix
            m_RBlk->SetBlock(n, n, rmat);
 
            invrmat = MemoryManager<DNekMat>::AllocateSharedPtr(*rmat);
            invrmat->Invert();
 
            // Block inverse R matrix
            m_InvRBlk->SetBlock(n, n, invrmat);
 
            m_sameBlock.push_back(pair<int, int>(1, nbndcoeffs));
            locExpSav = locExp;
        }
        else
        {
            bool reuse = true;
 
            // check to see if same as previous matrix and
            // reuse if we can
            for (int i = 0; i < 3; ++i)
            {
                if (locExpSav->GetBasis(i) != locExp->GetBasis(i))
                {
                    reuse = false;
                    break;
                }
            }
 
            if (reuse)
            {
                m_RBlk->SetBlock(n, n, rmat);
                m_InvRBlk->SetBlock(n, n, invrmat);
 
                m_sameBlock[offset] =
                    (pair<int, int>(m_sameBlock[offset].first + 1, nbndcoeffs));
            }
            else
            {
                rmat = ExtractLocMat(locExp, maxRmat[eltype], maxElmt[eltype],
                                     vertMapMaxR[eltype], edgeMapMaxR[eltype]);
 
                // Block R matrix
                m_RBlk->SetBlock(n, n, rmat);
 
                invrmat = MemoryManager<DNekMat>::AllocateSharedPtr(*rmat);
                invrmat->Invert();
                // Block inverse R matrix
                m_InvRBlk->SetBlock(n, n, invrmat);
 
                m_sameBlock.push_back(pair<int, int>(1, nbndcoeffs));
                offset++;
                locExpSav = locExp;
            }
        }
    }
}

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::StdRegions::StdExpansion::DetShapeType(), Nektar::eDIAGONAL, ExtractLocMat(), m_InvRBlk, Nektar::MultiRegions::Preconditioner::m_linsys, Nektar::MultiRegions::Preconditioner::m_locToGloMap, m_RBlk, m_sameBlock, and SetUpReferenceElements().

Referenced by v_InitObject().

SetUpPyrMaxRMat()

void Nektar::MultiRegions::PreconditionerLowEnergy::SetUpPyrMaxRMat ( int  nummodesmax, LocalRegions::PyrExpSharedPtr &  PyrExp, ShapeToDNekMap &  maxRmat, ShapeToIntArrayMap &  vertMapMaxR, ShapeToIntArrayArrayMap &  edgeMapMaxR, ShapeToIntArrayArrayMap &  faceMapMaxR  )

private

Definition at line 1878 of file PreconditionerLowEnergy.cpp.

{
    int nRows = PyrExp->NumBndryCoeffs();
    NekDouble val;
    NekDouble zero = 0.0;
    DNekMatSharedPtr newmat =
        MemoryManager<DNekMat>::AllocateSharedPtr(nRows, nRows, zero, eFULL);
 
    // set diagonal to 1
    for (int i = 0; i < nRows; ++i)
    {
        newmat->SetValue(i, i, 1.0);
    }
 
    // The following lists specify the number of adjacent
    // edges to each vertex (nadj) then the Hex vert to
    // use for each pyramid ver in the vert-edge map (VEVert)
    // followed by the hex edge to use for each pyramid edge
    // in the vert-edge map (VEEdge)
    const int nadjedge[]     = {3, 3, 3, 3, 4};
    const int VEHexVert[][4] = {{0, 0, 0, -1},
                                {1, 1, 1, -1},
                                {2, 2, 2, -1},
                                {3, 3, 3, -1},
                                {4, 5, 6, 7}};
    const int VEHexEdge[][4] = {{0, 3, 4, -1},
                                {0, 1, 5, -1},
                                {1, 2, 6, -1},
                                {2, 3, 7, -1},
                                {4, 5, 6, 7}};
    const int VEPyrEdge[][4] = {{0, 3, 4, -1},
                                {0, 1, 5, -1},
                                {1, 2, 6, -1},
                                {2, 3, 7, -1},
                                {4, 5, 6, 7}};
 
    // fill vertex to edge coupling
    for (int v = 0; v < 5; ++v)
    {
        for (int e = 0; e < nadjedge[v]; ++e)
        {
            for (int i = 0; i < nummodesmax - 2; ++i)
            {
                // Note this is using wrong shape but gives
                // answer that seems to have correct error!
                val = (*maxRmat[eHexahedron])(
                    vertMapMaxR[eHexahedron][VEHexVert[v][e]],
                    edgeMapMaxR[eHexahedron][VEHexEdge[v][e]][i]);
                newmat->SetValue(vertMapMaxR[ePyramid][v],
                                 edgeMapMaxR[ePyramid][VEPyrEdge[v][e]][i],
                                 val);
            }
        }
    }
 
    int nfacemodes;
    nfacemodes = (nummodesmax - 2) * (nummodesmax - 2);
    // First four verties are all adjacent to base face
    for (int v = 0; v < 4; ++v)
    {
        for (int i = 0; i < nfacemodes; ++i)
        {
            val = (*maxRmat[eHexahedron])(vertMapMaxR[eHexahedron][v],
                                          faceMapMaxR[eHexahedron][0][i]);
            newmat->SetValue(vertMapMaxR[ePyramid][v],
                             faceMapMaxR[ePyramid][0][i], val);
        }
    }
 
    const int nadjface[]     = {2, 2, 2, 2, 4};
    const int VFTetVert[][4] = {{0, 0, -1, -1},
                                {1, 1, -1, -1},
                                {2, 2, -1, -1},
                                {0, 2, -1, -1},
                                {3, 3, 3, 3}};
    const int VFTetFace[][4] = {{1, 3, -1, -1},
                                {1, 2, -1, -1},
                                {2, 3, -1, -1},
                                {1, 3, -1, -1},
                                {1, 2, 1, 2}};
    const int VFPyrFace[][4] = {{1, 4, -1, -1},
                                {1, 2, -1, -1},
                                {2, 3, -1, -1},
                                {3, 4, -1, -1},
                                {1, 2, 3, 4}};
 
    // next handle all triangular faces from tetrahedron
    nfacemodes = (nummodesmax - 3) * (nummodesmax - 2) / 2;
    for (int v = 0; v < 5; ++v)
    {
        for (int f = 0; f < nadjface[v]; ++f)
        {
            for (int i = 0; i < nfacemodes; ++i)
            {
                val = (*maxRmat[eTetrahedron])(
                    vertMapMaxR[eTetrahedron][VFTetVert[v][f]],
                    faceMapMaxR[eTetrahedron][VFTetFace[v][f]][i]);
                newmat->SetValue(vertMapMaxR[ePyramid][v],
                                 faceMapMaxR[ePyramid][VFPyrFace[v][f]][i],
                                 val);
            }
        }
    }
 
    // Edge to face coupling
    // all base edges are coupled to face 0
    nfacemodes = (nummodesmax - 2) * (nummodesmax - 2);
    for (int e = 0; e < 4; ++e)
    {
        for (int i = 0; i < nummodesmax - 2; ++i)
        {
            for (int j = 0; j < nfacemodes; ++j)
            {
                int edgemapid = edgeMapMaxR[eHexahedron][e][i];
                int facemapid = faceMapMaxR[eHexahedron][0][j];
 
                val = (*maxRmat[eHexahedron])(edgemapid, facemapid);
                newmat->SetValue(edgeMapMaxR[ePyramid][e][i],
                                 faceMapMaxR[ePyramid][0][j], val);
            }
        }
    }
 
    const int nadjface1[]    = {1, 1, 1, 1, 2, 2, 2, 2};
    const int EFTetEdge[][2] = {{0, -1}, {1, -1}, {0, -1}, {2, -1},
                                {3, 3},  {4, 4},  {5, 5},  {3, 5}};
    const int EFTetFace[][2] = {{1, -1}, {2, -1}, {1, -1}, {3, -1},
                                {1, 3},  {1, 2},  {2, 3},  {1, 3}};
    const int EFPyrFace[][2] = {{1, -1}, {2, -1}, {3, -1}, {4, -1},
                                {1, 4},  {1, 2},  {2, 3},  {3, 4}};
 
    // next handle all triangular faces from tetrahedron
    nfacemodes = (nummodesmax - 3) * (nummodesmax - 2) / 2;
    for (int e = 0; e < 8; ++e)
    {
        for (int f = 0; f < nadjface1[e]; ++f)
        {
            for (int i = 0; i < nummodesmax - 2; ++i)
            {
                for (int j = 0; j < nfacemodes; ++j)
                {
                    int edgemapid =
                        edgeMapMaxR[eTetrahedron][EFTetEdge[e][f]][i];
                    int facemapid =
                        faceMapMaxR[eTetrahedron][EFTetFace[e][f]][j];
 
                    val = (*maxRmat[eTetrahedron])(edgemapid, facemapid);
                    newmat->SetValue(edgeMapMaxR[ePyramid][e][i],
                                     faceMapMaxR[ePyramid][EFPyrFace[e][f]][j],
                                     val);
                }
            }
        }
    }
 
    DNekScalMatSharedPtr PyrR;
    PyrR = MemoryManager<DNekScalMat>::AllocateSharedPtr(1.0, newmat);
    maxRmat[ePyramid] = PyrR;
}

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::eFULL, Nektar::LibUtilities::eHexahedron, Nektar::LibUtilities::ePyramid, and Nektar::LibUtilities::eTetrahedron.

Referenced by SetUpReferenceElements().

SetUpReferenceElements()

void Nektar::MultiRegions::PreconditionerLowEnergy::SetUpReferenceElements ( ShapeToDNekMap &  maxRmat, ShapeToExpMap &  maxElmt, ShapeToIntArrayMap &  vertMapMaxR, ShapeToIntArrayArrayMap &  edgeMapMaxR  )

private

Loop expansion and determine different variants of the transformation matrix.

Sets up multiple reference elements based on the element expansion.

Definition at line 1662 of file PreconditionerLowEnergy.cpp.

{
    std::shared_ptr<MultiRegions::ExpList> expList =
        ((m_linsys.lock())->GetLocMat()).lock();
    GlobalLinSysKey linSysKey = (m_linsys.lock())->GetKey();
 
    LibUtilities::CommSharedPtr vComm = expList->GetComm()->GetSpaceComm();
 
    LocalRegions::ExpansionSharedPtr locExp;
 
    // face maps for pyramid and hybrid meshes - not need to return.
    map<ShapeType, Array<OneD, Array<OneD, unsigned int>>> faceMapMaxR;
 
    /* Determine the maximum expansion order for all elements */
    int nummodesmax = 0;
    Array<OneD, int> Shapes(LibUtilities::SIZE_ShapeType, 0);
 
    for (int n = 0; n < expList->GetNumElmts(); ++n)
    {
        locExp = expList->GetExp(n);
 
        nummodesmax = max(nummodesmax, locExp->GetBasisNumModes(0));
        nummodesmax = max(nummodesmax, locExp->GetBasisNumModes(1));
        nummodesmax = max(nummodesmax, locExp->GetBasisNumModes(2));
 
        Shapes[locExp->DetShapeType()] = 1;
    }
 
    vComm->AllReduce(nummodesmax, ReduceMax);
    vComm->AllReduce(Shapes, ReduceMax);
 
    if (Shapes[ePyramid] || Shapes[ePrism]) // if Pyramids or Prisms used then
                                            // need Tet and Hex expansion
    {
        Shapes[eTetrahedron] = 1;
        Shapes[eHexahedron]  = 1;
    }
 
    StdRegions::MatrixType PreconR;
    if (linSysKey.GetMatrixType() == StdRegions::eMass)
    {
        PreconR = StdRegions::ePreconRMass;
    }
    else
    {
        PreconR = StdRegions::ePreconR;
    }
 
    Array<OneD, unsigned int> vmap;
    Array<OneD, Array<OneD, unsigned int>> emap;
    Array<OneD, Array<OneD, unsigned int>> fmap;
 
    /*
     * Set up a transformation matrices for equal max order
     * polynomial meshes
     */
 
    if (Shapes[eHexahedron])
    {
        SpatialDomains::HexGeomUniquePtr hexgeom = CreateRefHexGeom();
        // Bases for Hex element
        const BasisKey HexBa(eModified_A, nummodesmax,
                             PointsKey(nummodesmax + 1, eGaussLobattoLegendre));
        const BasisKey HexBb(eModified_A, nummodesmax,
                             PointsKey(nummodesmax + 1, eGaussLobattoLegendre));
        const BasisKey HexBc(eModified_A, nummodesmax,
                             PointsKey(nummodesmax + 1, eGaussLobattoLegendre));
 
        // Create reference Hexahdedral expansion
        LocalRegions::HexExpSharedPtr HexExp;
 
        HexExp = MemoryManager<LocalRegions::HexExp>::AllocateSharedPtr(
            HexBa, HexBb, HexBc, hexgeom.get());
        m_holder.m_hexVec.push_back(std::move(hexgeom));
 
        maxElmt[eHexahedron] = HexExp;
 
        // Hex:
        HexExp->GetInverseBoundaryMaps(vmap, emap, fmap);
        vertMapMaxR[eHexahedron] = vmap;
        edgeMapMaxR[eHexahedron] = emap;
        faceMapMaxR[eHexahedron] = fmap;
 
        // Get hexahedral transformation matrix
        LocalRegions::MatrixKey HexR(PreconR, eHexahedron, *HexExp,
                                     linSysKey.GetConstFactors());
        maxRmat[eHexahedron] = HexExp->GetLocMatrix(HexR);
        HexExp->DropLocMatrix(HexR); // Need to delete reference from manager
    }
 
    if (Shapes[eTetrahedron])
    {
        SpatialDomains::TetGeomUniquePtr tetgeom = CreateRefTetGeom();
        // Bases for Tetrahedral element
        const BasisKey TetBa(eModified_A, nummodesmax,
                             PointsKey(nummodesmax + 1, eGaussLobattoLegendre));
        const BasisKey TetBb(eModified_B, nummodesmax,
                             PointsKey(nummodesmax, eGaussRadauMAlpha1Beta0));
        const BasisKey TetBc(eModified_C, nummodesmax,
                             PointsKey(nummodesmax, eGaussRadauMAlpha2Beta0));
 
        // Create reference tetrahedral expansion
        LocalRegions::TetExpSharedPtr TetExp;
 
        TetExp = MemoryManager<LocalRegions::TetExp>::AllocateSharedPtr(
            TetBa, TetBb, TetBc, tetgeom.get());
        m_holder.m_tetVec.push_back(std::move(tetgeom));
 
        maxElmt[eTetrahedron] = TetExp;
 
        TetExp->GetInverseBoundaryMaps(vmap, emap, fmap);
        vertMapMaxR[eTetrahedron] = vmap;
        edgeMapMaxR[eTetrahedron] = emap;
        faceMapMaxR[eTetrahedron] = fmap;
 
        // Get tetrahedral transformation matrix
        LocalRegions::MatrixKey TetR(PreconR, eTetrahedron, *TetExp,
                                     linSysKey.GetConstFactors());
        maxRmat[eTetrahedron] = TetExp->GetLocMatrix(TetR);
        TetExp->DropLocMatrix(TetR); // Need to delete reference from manager
 
        if ((Shapes[ePyramid]) || (Shapes[eHexahedron]))
        {
            ReSetTetMaxRMat(nummodesmax, TetExp, maxRmat, vertMapMaxR,
                            edgeMapMaxR, faceMapMaxR);
        }
    }
 
    if (Shapes[ePyramid])
    {
        SpatialDomains::PyrGeomUniquePtr pyrgeom = CreateRefPyrGeom();
 
        // Bases for Pyramid element
        const BasisKey PyrBa(eModified_A, nummodesmax,
                             PointsKey(nummodesmax + 1, eGaussLobattoLegendre));
        const BasisKey PyrBb(eModified_A, nummodesmax,
                             PointsKey(nummodesmax + 1, eGaussLobattoLegendre));
        const BasisKey PyrBc(eModifiedPyr_C, nummodesmax,
                             PointsKey(nummodesmax, eGaussRadauMAlpha2Beta0));
 
        // Create reference pyramid expansion
        LocalRegions::PyrExpSharedPtr PyrExp;
 
        PyrExp = MemoryManager<LocalRegions::PyrExp>::AllocateSharedPtr(
            PyrBa, PyrBb, PyrBc, pyrgeom.get());
        m_holder.m_pyrVec.push_back(std::move(pyrgeom));
 
        maxElmt[ePyramid] = PyrExp;
 
        // Pyramid:
        PyrExp->GetInverseBoundaryMaps(vmap, emap, fmap);
        vertMapMaxR[ePyramid] = vmap;
        edgeMapMaxR[ePyramid] = emap;
        faceMapMaxR[ePyramid] = fmap;
 
        // Set up Pyramid Transformation Matrix based on Tet
        // and Hex expansion
        SetUpPyrMaxRMat(nummodesmax, PyrExp, maxRmat, vertMapMaxR, edgeMapMaxR,
                        faceMapMaxR);
    }
 
    if (Shapes[ePrism])
    {
        SpatialDomains::PrismGeomUniquePtr prismgeom = CreateRefPrismGeom();
        // Bases for Prism element
        const BasisKey PrismBa(
            eModified_A, nummodesmax,
            PointsKey(nummodesmax + 1, eGaussLobattoLegendre));
        const BasisKey PrismBb(
            eModified_A, nummodesmax,
            PointsKey(nummodesmax + 1, eGaussLobattoLegendre));
        const BasisKey PrismBc(eModified_B, nummodesmax,
                               PointsKey(nummodesmax, eGaussRadauMAlpha1Beta0));
 
        // Create reference prismatic expansion
        LocalRegions::PrismExpSharedPtr PrismExp;
 
        PrismExp = MemoryManager<LocalRegions::PrismExp>::AllocateSharedPtr(
            PrismBa, PrismBb, PrismBc, prismgeom.get());
        m_holder.m_prismVec.push_back(std::move(prismgeom));
        maxElmt[ePrism] = PrismExp;
 
        // Prism:
        PrismExp->GetInverseBoundaryMaps(vmap, emap, fmap);
        vertMapMaxR[ePrism] = vmap;
        edgeMapMaxR[ePrism] = emap;
 
        faceMapMaxR[ePrism] = fmap;
 
        if ((Shapes[ePyramid]) || (Shapes[eHexahedron]))
        {
            ReSetPrismMaxRMat(nummodesmax, PrismExp, maxRmat, vertMapMaxR,
                              edgeMapMaxR, faceMapMaxR, false);
        }
        else
        {
            // Get prismatic transformation matrix
            LocalRegions::MatrixKey PrismR(PreconR, ePrism, *PrismExp,
                                           linSysKey.GetConstFactors());
            maxRmat[ePrism] = PrismExp->GetLocMatrix(PrismR);
            PrismExp->DropLocMatrix(
                PrismR); // Need to delete reference from manager
 
            if (Shapes[eTetrahedron]) // reset triangular faces from Tet
            {
                ReSetPrismMaxRMat(nummodesmax, PrismExp, maxRmat, vertMapMaxR,
                                  edgeMapMaxR, faceMapMaxR, true);
            }
        }
    }
}

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), CreateRefHexGeom(), CreateRefPrismGeom(), CreateRefPyrGeom(), CreateRefTetGeom(), Nektar::LibUtilities::eGaussLobattoLegendre, Nektar::LibUtilities::eHexahedron, Nektar::StdRegions::eMass, Nektar::LibUtilities::eModified_A, Nektar::LibUtilities::eModified_B, Nektar::LibUtilities::eModified_C, Nektar::LibUtilities::eModifiedPyr_C, Nektar::StdRegions::ePreconR, Nektar::StdRegions::ePreconRMass, Nektar::LibUtilities::ePrism, Nektar::LibUtilities::ePyramid, Nektar::LibUtilities::eTetrahedron, Nektar::MultiRegions::GlobalMatrixKey::GetConstFactors(), Nektar::MultiRegions::GlobalMatrixKey::GetMatrixType(), Nektar::SpatialDomains::EntityHolder::m_hexVec, m_holder, Nektar::MultiRegions::Preconditioner::m_linsys, Nektar::SpatialDomains::EntityHolder::m_prismVec, Nektar::SpatialDomains::EntityHolder::m_pyrVec, Nektar::SpatialDomains::EntityHolder::m_tetVec, tinysimd::max(), Nektar::LibUtilities::ReduceMax, ReSetPrismMaxRMat(), ReSetTetMaxRMat(), SetUpPyrMaxRMat(), and Nektar::LibUtilities::SIZE_ShapeType.

Referenced by SetupBlockTransformationMatrix().

v_BuildPreconditioner()

void Nektar::MultiRegions::PreconditionerLowEnergy::v_BuildPreconditioner ( )

overrideprotectedvirtual

Construct the low energy preconditioner from \(\mathbf{S}_{2}\).

\[\mathbf{M}^{-1}=\left[\begin{array}{ccc} Diag[(\mathbf{S_{2}})_{vv}] & & \\ & (\mathbf{S}_{2})_{eb} & \\ & & (\mathbf{S}_{2})_{fb} \end{array}\right] \]

where \(\mathbf{R}\) is the transformation matrix and \(\mathbf{S}_{2}\) the Schur complement of the modified basis, given by

\[\mathbf{S}_{2}=\mathbf{R}\mathbf{S}_{1}\mathbf{R}^{T}\]

where \(\mathbf{S}_{1}\) is the local schur complement matrix for each element.

• Count edges, face and add up min edges and min face sizes

Reimplemented from Nektar::MultiRegions::Preconditioner.

Definition at line 115 of file PreconditionerLowEnergy.cpp.

{
    std::shared_ptr<MultiRegions::ExpList> expList =
        ((m_linsys.lock())->GetLocMat()).lock();
    LocalRegions::Expansion3DSharedPtr locExpansion;
    GlobalLinSysKey linSysKey = (m_linsys.lock())->GetKey();
 
    LibUtilities::CommSharedPtr vComm = expList->GetComm()->GetSpaceComm();
 
    int i, j, k;
    int nVerts, nEdges, nFaces;
    int eid, fid, n, cnt, nmodes, nedgemodes, nfacemodes;
    int nedgemodesloc;
    NekDouble zero = 0.0;
 
    int vMap1, vMap2, sign1, sign2;
    int m, v, eMap1, eMap2, fMap1, fMap2;
    int offset, globalrow, globalcol, nCoeffs;
 
    // Periodic information
    PeriodicMap periodicVerts;
    PeriodicMap periodicEdges;
    PeriodicMap periodicFaces;
    expList->GetPeriodicEntities(periodicVerts, periodicEdges, periodicFaces);
 
    // matrix storage
    MatrixStorage storage       = eFULL;
    MatrixStorage vertstorage   = eDIAGONAL;
    MatrixStorage blkmatStorage = eDIAGONAL;
 
    // local element static condensed matrices
    DNekScalBlkMatSharedPtr loc_mat;
    DNekScalMatSharedPtr bnd_mat;
 
    DNekMatSharedPtr pRSRT;
 
    DNekMat RS;
    DNekMat RSRT;
 
    auto asmMap = m_locToGloMap.lock();
 
    int nDirBnd      = asmMap->GetNumGlobalDirBndCoeffs();
    int nNonDirVerts = asmMap->GetNumNonDirVertexModes();
 
    // Vertex, edge and face preconditioner matrices
    DNekMatSharedPtr VertBlk = MemoryManager<DNekMat>::AllocateSharedPtr(
        nNonDirVerts, nNonDirVerts, zero, vertstorage);
 
    Array<OneD, NekDouble> vertArray(nNonDirVerts, 0.0);
    Array<OneD, long> VertBlockToUniversalMap(nNonDirVerts, -1);
 
    // maps for different element types
    int n_exp          = expList->GetNumElmts();
    int nNonDirEdgeIDs = asmMap->GetNumNonDirEdges();
    int nNonDirFaceIDs = asmMap->GetNumNonDirFaces();
 
    set<int> edgeDirMap;
    set<int> faceDirMap;
    map<int, int> uniqueEdgeMap;
    map<int, int> uniqueFaceMap;
 
    // this should be of size total number of local edges + faces
    Array<OneD, int> modeoffset(1 + nNonDirEdgeIDs + nNonDirFaceIDs, 0);
    Array<OneD, int> globaloffset(1 + nNonDirEdgeIDs + nNonDirFaceIDs, 0);
 
    const Array<OneD, const ExpListSharedPtr> &bndCondExp =
        expList->GetBndCondExpansions();
    LocalRegions::Expansion2DSharedPtr bndCondFaceExp;
    const Array<OneD, const SpatialDomains::BoundaryConditionShPtr>
        &bndConditions = expList->GetBndConditions();
 
    int meshVertId;
    int meshEdgeId;
    int meshFaceId;
 
    const Array<OneD, const int> &extradiredges = asmMap->GetExtraDirEdges();
    for (i = 0; i < extradiredges.size(); ++i)
    {
        meshEdgeId = extradiredges[i];
        edgeDirMap.insert(meshEdgeId);
    }
 
    // Determine which boundary edges and faces have dirichlet values
    for (i = 0; i < bndCondExp.size(); i++)
    {
        cnt = 0;
        for (j = 0; j < bndCondExp[i]->GetNumElmts(); j++)
        {
            bndCondFaceExp =
                std::dynamic_pointer_cast<LocalRegions::Expansion2D>(
                    bndCondExp[i]->GetExp(j));
            if (bndConditions[i]->GetBoundaryConditionType() ==
                SpatialDomains::eDirichlet)
            {
                for (k = 0; k < bndCondFaceExp->GetNtraces(); k++)
                {
                    meshEdgeId = bndCondFaceExp->GetGeom()->GetEid(k);
                    if (edgeDirMap.count(meshEdgeId) == 0)
                    {
                        edgeDirMap.insert(meshEdgeId);
                    }
                }
                meshFaceId = bndCondFaceExp->GetGeom()->GetGlobalID();
                faceDirMap.insert(meshFaceId);
            }
        }
    }
 
    int dof               = 0;
    int maxFaceDof        = 0;
    int maxEdgeDof        = 0;
    int nlocalNonDirEdges = 0;
    int nlocalNonDirFaces = 0;
    int ntotalentries     = 0;
 
    map<int, int> EdgeSize;
    map<int, int> FaceSize;
    map<int, pair<int, int>> FaceModes;
 
    /// -  Count  edges, face and add up min edges and min face sizes
    for (n = 0; n < n_exp; ++n)
    {
        locExpansion = expList->GetExp(n)->as<LocalRegions::Expansion3D>();
 
        nEdges = locExpansion->GetNedges();
        for (j = 0; j < nEdges; ++j)
        {
            int nEdgeInteriorCoeffs = locExpansion->GetEdgeNcoeffs(j) - 2;
            meshEdgeId = locExpansion->as<LocalRegions::Expansion3D>()
                             ->GetGeom3D()
                             ->GetEid(j);
            if (EdgeSize.count(meshEdgeId) == 0)
            {
                EdgeSize[meshEdgeId] = nEdgeInteriorCoeffs;
            }
            else
            {
                EdgeSize[meshEdgeId] =
                    min(EdgeSize[meshEdgeId], nEdgeInteriorCoeffs);
            }
        }
 
        nFaces = locExpansion->GetNtraces();
        for (j = 0; j < nFaces; ++j)
        {
            int nFaceInteriorCoeffs = locExpansion->GetTraceIntNcoeffs(j);
            meshFaceId              = locExpansion->GetGeom3D()->GetFid(j);
 
            if (FaceSize.count(meshFaceId) == 0)
            {
                FaceSize[meshFaceId] = nFaceInteriorCoeffs;
 
                int m0, m1;
                locExpansion->GetTraceNumModes(j, m0, m1,
                                               locExpansion->GetTraceOrient(j));
                FaceModes[meshFaceId] = pair<int, int>(m0, m1);
            }
            else
            {
                if (nFaceInteriorCoeffs < FaceSize[meshFaceId])
                {
                    FaceSize[meshFaceId] = nFaceInteriorCoeffs;
                    int m0, m1;
                    locExpansion->GetTraceNumModes(
                        j, m0, m1, locExpansion->GetTraceOrient(j));
                    FaceModes[meshFaceId] = pair<int, int>(m0, m1);
                }
            }
        }
    }
 
    bool verbose = expList->GetSession()->DefinesCmdLineArgument("verbose");
 
    // For parallel runs need to check have minimum of edges and faces over
    // partition boundaries
    if (vComm->GetSize() > 1)
    {
        int EdgeSizeLen = EdgeSize.size();
        int FaceSizeLen = FaceSize.size();
        Array<OneD, long> FacetMap(EdgeSizeLen + FaceSizeLen, -1);
        Array<OneD, NekDouble> FacetLen(EdgeSizeLen + FaceSizeLen, -1);
 
        map<int, int>::iterator it;
 
        cnt       = 0;
        int maxid = 0;
        for (it = EdgeSize.begin(); it != EdgeSize.end(); ++it, ++cnt)
        {
            FacetMap[cnt] = it->first;
            maxid         = max(it->first, maxid);
            FacetLen[cnt] = it->second;
        }
        maxid++;
 
        vComm->AllReduce(maxid, ReduceMax);
 
        for (it = FaceSize.begin(); it != FaceSize.end(); ++it, ++cnt)
        {
            FacetMap[cnt] = it->first + maxid;
            FacetLen[cnt] = it->second;
        }
 
        // Exchange vertex data over different processes
        Gs::gs_data *tmp = Gs::Init(FacetMap, vComm, verbose);
        Gs::Gather(FacetLen, Gs::gs_min, tmp);
        Gs::Free(tmp);
 
        cnt = 0;
        for (it = EdgeSize.begin(); it != EdgeSize.end(); ++it, ++cnt)
        {
            it->second = (int)FacetLen[cnt];
        }
 
        for (it = FaceSize.begin(); it != FaceSize.end(); ++it, ++cnt)
        {
            it->second = (int)FacetLen[cnt];
        }
    }
 
    // Loop over all the elements in the domain and compute total edge
    // DOF and set up unique ordering.
    map<int, int> nblks;
    int matrixlocation = 0;
 
    // First do periodic edges
    for (auto &pIt : periodicEdges)
    {
        meshEdgeId = pIt.first;
 
        if (edgeDirMap.count(meshEdgeId) == 0)
        {
            dof = EdgeSize[meshEdgeId];
            if (uniqueEdgeMap.count(meshEdgeId) == 0 && dof > 0)
            {
                bool SetUpNewEdge = true;
 
                for (i = 0; i < pIt.second.size(); ++i)
                {
                    if (!pIt.second[i].isLocal)
                    {
                        continue;
                    }
 
                    int meshEdgeId2 = pIt.second[i].id;
                    if (edgeDirMap.count(meshEdgeId2) == 0)
                    {
                        if (uniqueEdgeMap.count(meshEdgeId2) != 0)
                        {
                            // set unique map to same location
                            uniqueEdgeMap[meshEdgeId] =
                                uniqueEdgeMap[meshEdgeId2];
                            SetUpNewEdge = false;
                        }
                    }
                    else
                    {
                        edgeDirMap.insert(meshEdgeId);
                        SetUpNewEdge = false;
                    }
                }
 
                if (SetUpNewEdge)
                {
                    uniqueEdgeMap[meshEdgeId] = matrixlocation;
                    globaloffset[matrixlocation] += ntotalentries;
                    modeoffset[matrixlocation] = dof * dof;
                    ntotalentries += dof * dof;
                    nblks[matrixlocation++] = dof;
                }
            }
        }
    }
 
    for (cnt = n = 0; n < n_exp; ++n)
    {
        locExpansion = expList->GetExp(n)->as<LocalRegions::Expansion3D>();
 
        for (j = 0; j < locExpansion->GetNedges(); ++j)
        {
            meshEdgeId = locExpansion->GetGeom()->GetEid(j);
            dof        = EdgeSize[meshEdgeId];
            maxEdgeDof = (dof > maxEdgeDof ? dof : maxEdgeDof);
 
            if (edgeDirMap.count(meshEdgeId) == 0)
            {
                if (uniqueEdgeMap.count(meshEdgeId) == 0 && dof > 0)
 
                {
                    uniqueEdgeMap[meshEdgeId] = matrixlocation;
 
                    globaloffset[matrixlocation] += ntotalentries;
                    modeoffset[matrixlocation] = dof * dof;
                    ntotalentries += dof * dof;
                    nblks[matrixlocation++] = dof;
                }
                nlocalNonDirEdges += dof * dof;
            }
        }
    }
 
    // Loop over all the elements in the domain and compute max face
    // DOF. Reduce across all processes to get universal maximum.
    // - Periodic faces
    for (auto &pIt : periodicFaces)
    {
        meshFaceId = pIt.first;
 
        if (faceDirMap.count(meshFaceId) == 0)
        {
            dof = FaceSize[meshFaceId];
 
            if (uniqueFaceMap.count(meshFaceId) == 0 && dof > 0)
            {
                bool SetUpNewFace = true;
 
                if (pIt.second[0].isLocal)
                {
                    int meshFaceId2 = pIt.second[0].id;
 
                    if (faceDirMap.count(meshFaceId2) == 0)
                    {
                        if (uniqueFaceMap.count(meshFaceId2) != 0)
                        {
                            // set unique map to same location
                            uniqueFaceMap[meshFaceId] =
                                uniqueFaceMap[meshFaceId2];
                            SetUpNewFace = false;
                        }
                    }
                    else // set face to be a Dirichlet face
                    {
                        faceDirMap.insert(meshFaceId);
                        SetUpNewFace = false;
                    }
                }
 
                if (SetUpNewFace)
                {
                    uniqueFaceMap[meshFaceId] = matrixlocation;
 
                    modeoffset[matrixlocation] = dof * dof;
                    globaloffset[matrixlocation] += ntotalentries;
                    ntotalentries += dof * dof;
                    nblks[matrixlocation++] = dof;
                }
            }
        }
    }
 
    for (cnt = n = 0; n < n_exp; ++n)
    {
        locExpansion = expList->GetExp(n)->as<LocalRegions::Expansion3D>();
 
        for (j = 0; j < locExpansion->GetNtraces(); ++j)
        {
            meshFaceId = locExpansion->GetGeom()->GetFid(j);
 
            dof        = FaceSize[meshFaceId];
            maxFaceDof = (dof > maxFaceDof ? dof : maxFaceDof);
 
            if (faceDirMap.count(meshFaceId) == 0)
            {
                if (uniqueFaceMap.count(meshFaceId) == 0 && dof > 0)
                {
                    uniqueFaceMap[meshFaceId]  = matrixlocation;
                    modeoffset[matrixlocation] = dof * dof;
                    globaloffset[matrixlocation] += ntotalentries;
                    ntotalentries += dof * dof;
                    nblks[matrixlocation++] = dof;
                }
                nlocalNonDirFaces += dof * dof;
            }
        }
    }
 
    vComm->AllReduce(maxEdgeDof, ReduceMax);
    vComm->AllReduce(maxFaceDof, ReduceMax);
 
    // Allocate arrays for block to universal map (number of expansions * p^2)
    Array<OneD, long> BlockToUniversalMap(ntotalentries, -1);
    Array<OneD, int> localToGlobalMatrixMap(
        nlocalNonDirEdges + nlocalNonDirFaces, -1);
 
    // Allocate arrays to store matrices (number of expansions * p^2)
    Array<OneD, NekDouble> BlockArray(nlocalNonDirEdges + nlocalNonDirFaces,
                                      0.0);
 
    int matrixoffset = 0;
    int vGlobal;
    int uniEdgeOffset = 0;
 
    // Need to obtain a fixed offset for the universal number
    // of the faces which come after the edge numbering
    for (n = 0; n < n_exp; ++n)
    {
        for (j = 0; j < locExpansion->GetNedges(); ++j)
        {
            meshEdgeId = locExpansion->as<LocalRegions::Expansion3D>()
                             ->GetGeom3D()
                             ->GetEid(j);
 
            uniEdgeOffset = max(meshEdgeId, uniEdgeOffset);
        }
    }
    uniEdgeOffset++;
 
    vComm->AllReduce(uniEdgeOffset, ReduceMax);
    uniEdgeOffset = uniEdgeOffset * maxEdgeDof * maxEdgeDof;
 
    for (n = 0; n < n_exp; ++n)
    {
        locExpansion = expList->GetExp(n)->as<LocalRegions::Expansion3D>();
 
        // loop over the edges of the expansion
        for (j = 0; j < locExpansion->GetNedges(); ++j)
        {
            // get mesh edge id
            meshEdgeId = locExpansion->GetGeom()->GetEid(j);
 
            nedgemodes = EdgeSize[meshEdgeId];
 
            if (edgeDirMap.count(meshEdgeId) == 0 && nedgemodes > 0)
            {
                // Determine the Global edge offset
                int edgeOffset = globaloffset[uniqueEdgeMap[meshEdgeId]];
 
                // Determine a universal map offset
                int uniOffset = meshEdgeId;
                auto pIt      = periodicEdges.find(meshEdgeId);
                if (pIt != periodicEdges.end())
                {
                    for (int l = 0; l < pIt->second.size(); ++l)
                    {
                        uniOffset = min(uniOffset, pIt->second[l].id);
                    }
                }
                uniOffset = uniOffset * maxEdgeDof * maxEdgeDof;
 
                for (k = 0; k < nedgemodes * nedgemodes; ++k)
                {
                    vGlobal                                  = edgeOffset + k;
                    localToGlobalMatrixMap[matrixoffset + k] = vGlobal;
                    BlockToUniversalMap[vGlobal] = uniOffset + k + 1;
                }
                matrixoffset += nedgemodes * nedgemodes;
            }
        }
 
        Array<OneD, unsigned int> faceInteriorMap;
        Array<OneD, int> faceInteriorSign;
        // loop over the faces of the expansion
        for (j = 0; j < locExpansion->GetNtraces(); ++j)
        {
            // get mesh face id
            meshFaceId = locExpansion->GetGeom()->GetFid(j);
 
            nfacemodes = FaceSize[meshFaceId];
 
            // Check if face has dirichlet values
            if (faceDirMap.count(meshFaceId) == 0 && nfacemodes > 0)
            {
                // Determine the Global edge offset
                int faceOffset = globaloffset[uniqueFaceMap[meshFaceId]];
                // Determine a universal map offset
                int uniOffset = meshFaceId;
                // use minimum face edge when periodic
                auto pIt = periodicFaces.find(meshFaceId);
                if (pIt != periodicFaces.end())
                {
                    uniOffset = min(uniOffset, pIt->second[0].id);
                }
                uniOffset = uniOffset * maxFaceDof * maxFaceDof;
 
                for (k = 0; k < nfacemodes * nfacemodes; ++k)
                {
                    vGlobal = faceOffset + k;
 
                    localToGlobalMatrixMap[matrixoffset + k] = vGlobal;
 
                    BlockToUniversalMap[vGlobal] =
                        uniOffset + uniEdgeOffset + k + 1;
                }
                matrixoffset += nfacemodes * nfacemodes;
            }
        }
    }
 
    matrixoffset = 0;
 
    map<int, int>::iterator it;
    Array<OneD, unsigned int> n_blks(nblks.size() + 1);
    n_blks[0] = nNonDirVerts;
    for (i = 1, it = nblks.begin(); it != nblks.end(); ++it)
    {
        n_blks[i++] = it->second;
    }
 
    m_BlkMat = MemoryManager<DNekBlkMat>::AllocateSharedPtr(n_blks, n_blks,
                                                            blkmatStorage);
 
    // Here we loop over the expansion and build the block low energy
    // preconditioner as well as the block versions of the transformation
    // matrices.
    for (cnt = n = 0; n < n_exp; ++n)
    {
        locExpansion = expList->GetExp(n)->as<LocalRegions::Expansion3D>();
        nCoeffs      = locExpansion->NumBndryCoeffs();
 
        // Get correct transformation matrix for element type
        DNekMat &R = (*m_RBlk->GetBlock(n, n));
 
        pRSRT = MemoryManager<DNekMat>::AllocateSharedPtr(nCoeffs, nCoeffs,
                                                          zero, storage);
        RSRT  = (*pRSRT);
 
        nVerts = locExpansion->GetGeom()->GetNumVerts();
        nEdges = locExpansion->GetGeom()->GetNumEdges();
        nFaces = locExpansion->GetGeom()->GetNumFaces();
 
        // Get statically condensed matrix
        loc_mat = (m_linsys.lock())->GetStaticCondBlock(n);
 
        // Extract boundary block (elemental S1)
        bnd_mat = loc_mat->GetBlock(0, 0);
 
        // offset by number of rows
        offset = bnd_mat->GetRows();
 
        DNekScalMat &S = (*bnd_mat);
 
        DNekMat Sloc(nCoeffs, nCoeffs);
 
        // For variable p we need to just use the active modes
        NekDouble val;
 
        for (int i = 0; i < nCoeffs; ++i)
        {
            for (int j = 0; j < nCoeffs; ++j)
            {
                val = S(i, j) * m_variablePmask[cnt + i] *
                      m_variablePmask[cnt + j];
                Sloc.SetValue(i, j, val);
            }
        }
 
        // Calculate R*S*trans(R)
        RSRT = R * Sloc * Transpose(R);
 
        // loop over vertices of the element and return the vertex map
        // for each vertex
        for (v = 0; v < nVerts; ++v)
        {
            vMap1 = locExpansion->GetVertexMap(v);
 
            // Get vertex map
            globalrow = asmMap->GetLocalToGlobalBndMap(cnt + vMap1) - nDirBnd;
 
            if (globalrow >= 0)
            {
                for (m = 0; m < nVerts; ++m)
                {
                    vMap2 = locExpansion->GetVertexMap(m);
 
                    // global matrix location (without offset due to
                    // dirichlet values)
                    globalcol =
                        asmMap->GetLocalToGlobalBndMap(cnt + vMap2) - nDirBnd;
 
                    // offset for dirichlet conditions
                    if (globalcol == globalrow)
                    {
                        // modal connectivity between elements
                        sign1 = asmMap->GetLocalToGlobalBndSign(cnt + vMap1);
                        sign2 = asmMap->GetLocalToGlobalBndSign(cnt + vMap2);
 
                        vertArray[globalrow] +=
                            sign1 * sign2 * RSRT(vMap1, vMap2);
 
                        meshVertId = locExpansion->GetGeom()->GetVid(v);
 
                        auto pIt = periodicVerts.find(meshVertId);
                        if (pIt != periodicVerts.end())
                        {
                            for (k = 0; k < pIt->second.size(); ++k)
                            {
                                meshVertId = min(meshVertId, pIt->second[k].id);
                            }
                        }
 
                        VertBlockToUniversalMap[globalrow] = meshVertId + 1;
                    }
                }
            }
        }
 
        // loop over edges of the element and return the edge map
        for (eid = 0; eid < nEdges; ++eid)
        {
            meshEdgeId = locExpansion->as<LocalRegions::Expansion3D>()
                             ->GetGeom3D()
                             ->GetEid(eid);
 
            nedgemodes = EdgeSize[meshEdgeId];
            if (nedgemodes)
            {
                nedgemodesloc = locExpansion->GetEdgeNcoeffs(eid) - 2;
                DNekMatSharedPtr m_locMat =
                    MemoryManager<DNekMat>::AllocateSharedPtr(
                        nedgemodes, nedgemodes, zero, storage);
 
                Array<OneD, unsigned int> edgemodearray =
                    locExpansion->GetEdgeInverseBoundaryMap(eid);
 
                if (edgeDirMap.count(meshEdgeId) == 0)
                {
                    for (v = 0; v < nedgemodesloc; ++v)
                    {
                        eMap1 = edgemodearray[v];
                        sign1 = asmMap->GetLocalToGlobalBndSign(cnt + eMap1);
 
                        if (sign1 == 0)
                        {
                            continue;
                        }
 
                        for (m = 0; m < nedgemodesloc; ++m)
                        {
                            eMap2 = edgemodearray[m];
 
                            // modal connectivity between elements
                            sign2 =
                                asmMap->GetLocalToGlobalBndSign(cnt + eMap2);
 
                            NekDouble globalEdgeValue =
                                sign1 * sign2 * RSRT(eMap1, eMap2);
 
                            if (sign2 != 0)
                            {
                                // if(eMap1 == eMap2)
                                BlockArray[matrixoffset + v * nedgemodes + m] =
                                    globalEdgeValue;
                            }
                        }
                    }
                    matrixoffset += nedgemodes * nedgemodes;
                }
            }
        }
 
        // loop over faces of the element and return the face map
        for (fid = 0; fid < nFaces; ++fid)
        {
            meshFaceId = locExpansion->as<LocalRegions::Expansion3D>()
                             ->GetGeom3D()
                             ->GetFid(fid);
 
            nfacemodes = FaceSize[meshFaceId];
            if (nfacemodes > 0)
            {
                DNekMatSharedPtr m_locMat =
                    MemoryManager<DNekMat>::AllocateSharedPtr(
                        nfacemodes, nfacemodes, zero, storage);
 
                if (faceDirMap.count(meshFaceId) == 0)
                {
                    Array<OneD, unsigned int> facemodearray;
                    StdRegions::Orientation faceOrient =
                        locExpansion->GetTraceOrient(fid);
 
                    auto pIt = periodicFaces.find(meshFaceId);
                    if (pIt != periodicFaces.end())
                    {
                        if (meshFaceId == min(meshFaceId, pIt->second[0].id))
                        {
                            faceOrient = DeterminePeriodicFaceOrient(
                                faceOrient, pIt->second[0].orient);
                        }
                    }
 
                    facemodearray = locExpansion->GetTraceInverseBoundaryMap(
                        fid, faceOrient, FaceModes[meshFaceId].first,
                        FaceModes[meshFaceId].second);
 
                    for (v = 0; v < nfacemodes; ++v)
                    {
                        fMap1 = facemodearray[v];
 
                        sign1 = asmMap->GetLocalToGlobalBndSign(cnt + fMap1);
 
                        ASSERTL1(sign1 != 0,
                                 "Something is wrong since we "
                                 "shoudl only be extracting modes for "
                                 "lowest order expansion");
 
                        for (m = 0; m < nfacemodes; ++m)
                        {
                            fMap2 = facemodearray[m];
 
                            // modal connectivity between elements
                            sign2 =
                                asmMap->GetLocalToGlobalBndSign(cnt + fMap2);
 
                            ASSERTL1(sign2 != 0,
                                     "Something is wrong since "
                                     "we shoudl only be extracting modes for "
                                     "lowest order expansion");
 
                            // Get the face-face value from the
                            // low energy matrix (S2)
                            NekDouble globalFaceValue =
                                sign1 * sign2 * RSRT(fMap1, fMap2);
 
                            // local face value to global face value
                            // if(fMap1 == fMap2)
                            BlockArray[matrixoffset + v * nfacemodes + m] =
                                globalFaceValue;
                        }
                    }
                    matrixoffset += nfacemodes * nfacemodes;
                }
            }
        }
 
        // offset for the expansion
        cnt += nCoeffs;
    }
 
    // Exchange vertex data over different processes
    Gs::gs_data *tmp = Gs::Init(VertBlockToUniversalMap, vComm, verbose);
    Gs::Gather(vertArray, Gs::gs_add, tmp);
    Gs::Free(tmp);
 
    Array<OneD, NekDouble> GlobalBlock(ntotalentries, 0.0);
    if (ntotalentries)
    {
        // Assemble edge matrices of each process
        Vmath::Assmb(BlockArray.size(), BlockArray, localToGlobalMatrixMap,
                     GlobalBlock);
    }
 
    // Exchange edge & face data over different processes
    Gs::gs_data *tmp1 = Gs::Init(BlockToUniversalMap, vComm, verbose);
    Gs::Gather(GlobalBlock, Gs::gs_add, tmp1);
    Gs::Free(tmp1);
 
    // Populate vertex block
    for (int i = 0; i < nNonDirVerts; ++i)
    {
        VertBlk->SetValue(i, i, 1.0 / vertArray[i]);
    }
 
    // Set the first block to be the diagonal of the vertex space
    m_BlkMat->SetBlock(0, 0, VertBlk);
 
    // Build the edge and face matrices from the vector
    DNekMatSharedPtr gmat;
 
    offset = 0;
    // -1 since we ignore vert block
    for (int loc = 0; loc < n_blks.size() - 1; ++loc)
    {
        nmodes = n_blks[1 + loc];
        gmat   = MemoryManager<DNekMat>::AllocateSharedPtr(nmodes, nmodes, zero,
                                                           storage);
 
        for (v = 0; v < nmodes; ++v)
        {
            for (m = 0; m < nmodes; ++m)
            {
                NekDouble Value = GlobalBlock[offset + v * nmodes + m];
                gmat->SetValue(v, m, Value);
            }
        }
        m_BlkMat->SetBlock(1 + loc, 1 + loc, gmat);
        offset += modeoffset[loc];
    }
 
    // invert blocks.
    int totblks = m_BlkMat->GetNumberOfBlockRows();
    for (i = 1; i < totblks; ++i)
    {
        unsigned int nmodes = m_BlkMat->GetNumberOfRowsInBlockRow(i);
        if (nmodes)
        {
            DNekMatSharedPtr tmp_mat =
                MemoryManager<DNekMat>::AllocateSharedPtr(nmodes, nmodes, zero,
                                                          storage);
 
            tmp_mat = m_BlkMat->GetBlock(i, i);
            tmp_mat->Invert();
 
            m_BlkMat->SetBlock(i, i, tmp_mat);
        }
    }
}

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL1, Vmath::Assmb(), Nektar::MultiRegions::DeterminePeriodicFaceOrient(), Nektar::eDIAGONAL, Nektar::SpatialDomains::eDirichlet, Nektar::eFULL, Gs::Free(), Gs::Gather(), Nektar::StdRegions::StdExpansion3D::GetEdgeNcoeffs(), Nektar::SpatialDomains::Geometry::GetEid(), Nektar::SpatialDomains::Geometry::GetFid(), Nektar::LocalRegions::Expansion::GetGeom(), Nektar::StdRegions::StdExpansion3D::GetNedges(), Gs::gs_add, Gs::gs_min, Gs::Init(), m_BlkMat, Nektar::MultiRegions::Preconditioner::m_linsys, Nektar::MultiRegions::Preconditioner::m_locToGloMap, m_RBlk, m_variablePmask, tinysimd::max(), tinysimd::min(), Nektar::StdRegions::StdExpansion::NumBndryCoeffs(), Nektar::LibUtilities::ReduceMax, and Nektar::Transpose().

v_DoPreconditioner()

void Nektar::MultiRegions::PreconditionerLowEnergy::v_DoPreconditioner ( const Array< OneD, NekDouble > &  pInput, Array< OneD, NekDouble > &  pOutput, const bool &  isLocal = false  )

overrideprotectedvirtual

Apply the low energy preconditioner during the conjugate gradient routine

Reimplemented from Nektar::MultiRegions::Preconditioner.

Definition at line 915 of file PreconditionerLowEnergy.cpp.

{
    ASSERTL0(isLocal == false, "PreconditionerLowEnergy is only currently "
                               "set up for Global iteratives sovles");
    int nDir      = m_locToGloMap.lock()->GetNumGlobalDirBndCoeffs();
    int nGlobal   = m_locToGloMap.lock()->GetNumGlobalBndCoeffs();
    int nNonDir   = nGlobal - nDir;
    DNekBlkMat &M = (*m_BlkMat);
 
    NekVector<NekDouble> r(nNonDir, pInput, eWrapper);
 
    NekVector<NekDouble> z(nNonDir, pOutput, eWrapper);
 
    z = M * r;
}

References ASSERTL0, Nektar::eWrapper, and Nektar::MultiRegions::Preconditioner::m_locToGloMap.

v_DoTransformBasisFromLowEnergy()

void Nektar::MultiRegions::PreconditionerLowEnergy::v_DoTransformBasisFromLowEnergy ( const Array< OneD, NekDouble > &  pInput, Array< OneD, NekDouble > &  pOutput  )

overrideprotectedvirtual

Multiply by the block inverse transformation matrix This transforms the bassi from Low Energy to original basis.

Note; not typically required

Reimplemented from Nektar::MultiRegions::Preconditioner.

Definition at line 1130 of file PreconditionerLowEnergy.cpp.

{
    int nLocBndDofs = m_locToGloMap.lock()->GetNumLocalBndCoeffs();
 
    ASSERTL1(pInput.size() >= nLocBndDofs,
             "Input array is smaller than nLocBndDofs");
    ASSERTL1(pOutput.size() >= nLocBndDofs,
             "Output array is smaller than nLocBndDofs");
 
    // Block inverse transformation matrix
    DNekBlkMat &invR = *m_InvRBlk;
 
    Array<OneD, NekDouble> pLocalIn(nLocBndDofs, pInput.data());
 
    // Multiply by the inverse transformation matrix
    int cnt  = 0;
    int cnt1 = 0;
    for (int i = 0; i < m_sameBlock.size(); ++i)
    {
        int nexp       = m_sameBlock[i].first;
        int nbndcoeffs = m_sameBlock[i].second;
        Blas::Dgemm('N', 'N', nbndcoeffs, nexp, nbndcoeffs, 1.0,
                    &(invR.GetBlock(cnt1, cnt1)->GetPtr()[0]), nbndcoeffs,
                    pLocalIn.data() + cnt, nbndcoeffs, 0.0,
                    pOutput.data() + cnt, nbndcoeffs);
        cnt += nbndcoeffs * nexp;
        cnt1 += nexp;
    }
}

References ASSERTL1, Blas::Dgemm(), m_InvRBlk, Nektar::MultiRegions::Preconditioner::m_locToGloMap, and m_sameBlock.

v_DoTransformBasisToLowEnergy()

void Nektar::MultiRegions::PreconditionerLowEnergy::v_DoTransformBasisToLowEnergy ( Array< OneD, NekDouble > &  pInOut )

overrideprotectedvirtual

Transform the basis vector to low energy space.

As the conjugate gradient system is solved for the low energy basis, the solution vector \(\mathbf{x}\) must be transformed to the low energy basis i.e. \(\overline{\mathbf{x}}=\mathbf{R}\mathbf{x}\).

Reimplemented from Nektar::MultiRegions::Preconditioner.

Definition at line 1054 of file PreconditionerLowEnergy.cpp.

{
    int nLocBndDofs = m_locToGloMap.lock()->GetNumLocalBndCoeffs();
 
    // Block transformation matrix
    DNekBlkMat &R = *m_RBlk;
 
    Array<OneD, NekDouble> pLocalIn(nLocBndDofs, pInOut.data());
 
    // Apply mask in case of variable P
    Vmath::Vmul(nLocBndDofs, pLocalIn, 1, m_variablePmask, 1, pLocalIn, 1);
 
    // Multiply by the block transformation matrix
    int cnt  = 0;
    int cnt1 = 0;
    for (int i = 0; i < m_sameBlock.size(); ++i)
    {
        int nexp       = m_sameBlock[i].first;
        int nbndcoeffs = m_sameBlock[i].second;
        Blas::Dgemm('N', 'N', nbndcoeffs, nexp, nbndcoeffs, 1.0,
                    &(R.GetBlock(cnt1, cnt1)->GetPtr()[0]), nbndcoeffs,
                    pLocalIn.data() + cnt, nbndcoeffs, 0.0, pInOut.data() + cnt,
                    nbndcoeffs);
        cnt += nbndcoeffs * nexp;
        cnt1 += nexp;
    }
}

References Blas::Dgemm(), Nektar::MultiRegions::Preconditioner::m_locToGloMap, m_RBlk, m_sameBlock, m_variablePmask, and Vmath::Vmul().

v_DoTransformCoeffsFromLowEnergy()

void Nektar::MultiRegions::PreconditionerLowEnergy::v_DoTransformCoeffsFromLowEnergy ( Array< OneD, NekDouble > &  pInOut )

overrideprotectedvirtual

transform the solution coeffiicents from low energy back to the original coefficient space.

After the conjugate gradient routine the output vector is in the low energy basis and must be trasnformed back to the original basis in order to get the correct solution out. the solution vector i.e. \(\mathbf{x}=\mathbf{R^{T}}\mathbf{\overline{x}}\).

Reimplemented from Nektar::MultiRegions::Preconditioner.

Definition at line 1092 of file PreconditionerLowEnergy.cpp.

{
    int nLocBndDofs = m_locToGloMap.lock()->GetNumLocalBndCoeffs();
 
    ASSERTL1(pInOut.size() >= nLocBndDofs,
             "Output array is not greater than the nLocBndDofs");
 
    // Block transposed transformation matrix
    DNekBlkMat &R = *m_RBlk;
 
    Array<OneD, NekDouble> pLocalIn(nLocBndDofs, pInOut.data());
 
    // Multiply by the transpose of block transformation matrix
    int cnt  = 0;
    int cnt1 = 0;
    for (int i = 0; i < m_sameBlock.size(); ++i)
    {
        int nexp       = m_sameBlock[i].first;
        int nbndcoeffs = m_sameBlock[i].second;
        Blas::Dgemm('T', 'N', nbndcoeffs, nexp, nbndcoeffs, 1.0,
                    &(R.GetBlock(cnt1, cnt1)->GetPtr()[0]), nbndcoeffs,
                    pLocalIn.data() + cnt, nbndcoeffs, 0.0, pInOut.data() + cnt,
                    nbndcoeffs);
        cnt += nbndcoeffs * nexp;
        cnt1 += nexp;
    }
 
    Vmath::Vmul(nLocBndDofs, pInOut, 1, m_variablePmask, 1, pInOut, 1);
}

References ASSERTL1, Blas::Dgemm(), Nektar::MultiRegions::Preconditioner::m_locToGloMap, m_RBlk, m_sameBlock, m_variablePmask, and Vmath::Vmul().

v_DoTransformCoeffsToLowEnergy()

void Nektar::MultiRegions::PreconditionerLowEnergy::v_DoTransformCoeffsToLowEnergy ( const Array< OneD, NekDouble > &  pInput, Array< OneD, NekDouble > &  pOutput  )

overrideprotectedvirtual

Multiply by the block tranposed inverse transformation matrix (R^T)^{-1} which is equivlaent to transforming coefficients to LowEnergy space.

In JCP 2001 paper on low energy this is seen as (C^T)^{-1}

Reimplemented from Nektar::MultiRegions::Preconditioner.

Definition at line 1168 of file PreconditionerLowEnergy.cpp.

{
    int nLocBndDofs = m_locToGloMap.lock()->GetNumLocalBndCoeffs();
 
    ASSERTL1(pInput.size() >= nLocBndDofs,
             "Input array is less than nLocBndDofs");
    ASSERTL1(pOutput.size() >= nLocBndDofs,
             "Output array is less than nLocBndDofs");
 
    // Block inverse transformation matrix
    DNekBlkMat &invR = *m_InvRBlk;
 
    Array<OneD, NekDouble> pLocalIn(nLocBndDofs, pInput.data());
 
    // Multiply by the transpose of block transformation matrix
    int cnt  = 0;
    int cnt1 = 0;
    for (int i = 0; i < m_sameBlock.size(); ++i)
    {
        int nexp       = m_sameBlock[i].first;
        int nbndcoeffs = m_sameBlock[i].second;
        Blas::Dgemm('T', 'N', nbndcoeffs, nexp, nbndcoeffs, 1.0,
                    &(invR.GetBlock(cnt1, cnt1)->GetPtr()[0]), nbndcoeffs,
                    pLocalIn.data() + cnt, nbndcoeffs, 0.0,
                    pOutput.data() + cnt, nbndcoeffs);
        cnt += nbndcoeffs * nexp;
        cnt1 += nexp;
    }
}

References ASSERTL1, Blas::Dgemm(), m_InvRBlk, Nektar::MultiRegions::Preconditioner::m_locToGloMap, and m_sameBlock.

v_InitObject()

void Nektar::MultiRegions::PreconditionerLowEnergy::v_InitObject ( )

overrideprotectedvirtual

Reimplemented from Nektar::MultiRegions::Preconditioner.

Definition at line 72 of file PreconditionerLowEnergy.cpp.

{
    GlobalSysSolnType solvertype = m_locToGloMap.lock()->GetGlobalSysSolnType();
 
    ASSERTL0(solvertype == eIterativeStaticCond ||
                 solvertype == ePETScStaticCond,
             "Solver type not valid");
 
    std::shared_ptr<MultiRegions::ExpList> expList =
        ((m_linsys.lock())->GetLocMat()).lock();
 
    LocalRegions::ExpansionSharedPtr locExpansion;
 
    locExpansion = expList->GetExp(0);
 
    int nDim = locExpansion->GetShapeDimension();
 
    ASSERTL0(nDim == 3, "Preconditioner type only valid in 3D");
 
    // Set up block transformation matrix
    SetupBlockTransformationMatrix();
 
    // Sets up variable p mask
    CreateVariablePMask();
}

References ASSERTL0, CreateVariablePMask(), Nektar::MultiRegions::eIterativeStaticCond, Nektar::MultiRegions::ePETScStaticCond, Nektar::MultiRegions::Preconditioner::m_linsys, Nektar::MultiRegions::Preconditioner::m_locToGloMap, and SetupBlockTransformationMatrix().

v_TransformedSchurCompl()

DNekScalMatSharedPtr Nektar::MultiRegions::PreconditionerLowEnergy::v_TransformedSchurCompl ( int  n, int  bndoffset, const std::shared_ptr< DNekScalMat > &  loc_mat  )

overrideprotectedvirtual

Set up the transformed block matrix system.

Sets up a block elemental matrix in which each of the block matrix is the low energy equivalent i.e. \(\mathbf{S}_{2}=\mathbf{R}\mathbf{S}_{1}\mathbf{R}^{T}\)

Reimplemented from Nektar::MultiRegions::Preconditioner.

Definition at line 1206 of file PreconditionerLowEnergy.cpp.

{
    std::shared_ptr<MultiRegions::ExpList> expList =
        ((m_linsys.lock())->GetLocMat()).lock();
 
    LocalRegions::ExpansionSharedPtr locExpansion;
    locExpansion      = expList->GetExp(n);
    unsigned int nbnd = locExpansion->NumBndryCoeffs();
 
    MatrixStorage storage = eFULL;
    DNekMatSharedPtr pS2 =
        MemoryManager<DNekMat>::AllocateSharedPtr(nbnd, nbnd, 0.0, storage);
 
    // transformation matrices
    DNekMat &R = (*m_RBlk->GetBlock(n, n));
 
    // Original Schur Complement
    DNekScalMat &S1 = (*loc_mat);
 
    DNekMat Sloc(nbnd, nbnd);
 
    // For variable p we need to just use the active modes
    NekDouble val;
 
    for (int i = 0; i < nbnd; ++i)
    {
        for (int j = 0; j < nbnd; ++j)
        {
            val = S1(i, j) * m_variablePmask[bndoffset + i] *
                  m_variablePmask[bndoffset + j];
            Sloc.SetValue(i, j, val);
        }
    }
 
    // create low energy matrix
    DNekMat &S2 = (*pS2);
 
    S2 = R * Sloc * Transpose(R);
 
    DNekScalMatSharedPtr return_val;
    return_val = MemoryManager<DNekScalMat>::AllocateSharedPtr(1.0, pS2);
 
    return return_val;
}

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::eFULL, Nektar::MultiRegions::Preconditioner::m_linsys, m_RBlk, m_variablePmask, and Nektar::Transpose().

Member Data Documentation

className

string Nektar::MultiRegions::PreconditionerLowEnergy::className

static

Initial value:

=
    GetPreconFactory().RegisterCreatorFunction("LowEnergyBlock",
                                               PreconditionerLowEnergy::create,
                                               "LowEnergy Preconditioning")

Name of class.

Registers the class with the Factory.

Definition at line 67 of file PreconditionerLowEnergy.h.

m_BlkMat

DNekBlkMatSharedPtr Nektar::MultiRegions::PreconditionerLowEnergy::m_BlkMat

protected

Definition at line 79 of file PreconditionerLowEnergy.h.

Referenced by v_BuildPreconditioner().

m_holder

SpatialDomains::EntityHolder Nektar::MultiRegions::PreconditionerLowEnergy::m_holder

protected

Definition at line 82 of file PreconditionerLowEnergy.h.

Referenced by CreateRefHexGeom(), CreateRefPrismGeom(), CreateRefPyrGeom(), CreateRefTetGeom(), and SetUpReferenceElements().

m_InvRBlk

DNekBlkMatSharedPtr Nektar::MultiRegions::PreconditionerLowEnergy::m_InvRBlk

protected

Definition at line 81 of file PreconditionerLowEnergy.h.

Referenced by SetupBlockTransformationMatrix(), v_DoTransformBasisFromLowEnergy(), and v_DoTransformCoeffsToLowEnergy().

m_RBlk

DNekBlkMatSharedPtr Nektar::MultiRegions::PreconditionerLowEnergy::m_RBlk

protected

Definition at line 80 of file PreconditionerLowEnergy.h.

Referenced by SetupBlockTransformationMatrix(), v_BuildPreconditioner(), v_DoTransformBasisToLowEnergy(), v_DoTransformCoeffsFromLowEnergy(), and v_TransformedSchurCompl().

m_sameBlock

std::vector<std::pair<int, int> > Nektar::MultiRegions::PreconditionerLowEnergy::m_sameBlock

protected

Definition at line 90 of file PreconditionerLowEnergy.h.

Referenced by SetupBlockTransformationMatrix(), v_DoTransformBasisFromLowEnergy(), v_DoTransformBasisToLowEnergy(), v_DoTransformCoeffsFromLowEnergy(), and v_DoTransformCoeffsToLowEnergy().

m_signChange

bool Nektar::MultiRegions::PreconditionerLowEnergy::m_signChange

protected

Definition at line 86 of file PreconditionerLowEnergy.h.

Referenced by CreateVariablePMask().

m_variablePmask

Array<OneD, NekDouble> Nektar::MultiRegions::PreconditionerLowEnergy::m_variablePmask

protected

Definition at line 84 of file PreconditionerLowEnergy.h.

Referenced by CreateVariablePMask(), v_BuildPreconditioner(), v_DoTransformBasisToLowEnergy(), v_DoTransformCoeffsFromLowEnergy(), and v_TransformedSchurCompl().