A global linear system. More...

#include <GlobalLinSysStaticCond.h>

Inheritance diagram for Nektar::MultiRegions::GlobalLinSysStaticCond:

Public Member Functions
	GlobalLinSysStaticCond (const GlobalLinSysKey &mkey, const std::weak_ptr< ExpList > &pExpList, const std::shared_ptr< AssemblyMap > &locToGloMap)
	Constructor for full direct matrix solve. More...

	~GlobalLinSysStaticCond () override

Public Member Functions inherited from Nektar::MultiRegions::GlobalLinSys
	GlobalLinSys (const GlobalLinSysKey &pKey, const std::weak_ptr< ExpList > &pExpList, const std::shared_ptr< AssemblyMap > &pLocToGloMap)
	Constructor for full direct matrix solve. More...

virtual	~GlobalLinSys ()

const GlobalLinSysKey &	GetKey (void) const
	Returns the key associated with the system. More...

const std::weak_ptr< ExpList > &	GetLocMat (void) const

void	InitObject ()

void	Initialise (const std::shared_ptr< AssemblyMap > &pLocToGloMap)

void	Solve (const Array< OneD, const NekDouble > &in, Array< OneD, NekDouble > &out, const AssemblyMapSharedPtr &locToGloMap, const Array< OneD, const NekDouble > &dirForcing=NullNekDouble1DArray)
	Solve the linear system for given input and output vectors using a specified local to global map. More...

std::shared_ptr< GlobalLinSys >	GetSharedThisPtr ()
	Returns a shared pointer to the current object. More...

int	GetNumBlocks ()

DNekScalMatSharedPtr	GetBlock (unsigned int n)

void	DropBlock (unsigned int n)

DNekScalBlkMatSharedPtr	GetStaticCondBlock (unsigned int n)

void	DropStaticCondBlock (unsigned int n)

void	SolveLinearSystem (const int pNumRows, const Array< OneD, const NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const AssemblyMapSharedPtr &locToGloMap, const int pNumDir=0)
	Solve the linear system for given input and output vectors. More...

Protected Member Functions
virtual void	v_PreSolve (int scLevel, Array< OneD, NekDouble > &F_bnd)

virtual void	v_BasisFwdTransform (Array< OneD, NekDouble > &pInOut)

virtual void	v_CoeffsBwdTransform (Array< OneD, NekDouble > &pInOut)

virtual void	v_CoeffsFwdTransform (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)

virtual void	v_AssembleSchurComplement (std::shared_ptr< AssemblyMap > pLoctoGloMap)

int	v_GetNumBlocks () override
	Get the number of blocks in this system. More...

virtual GlobalLinSysStaticCondSharedPtr	v_Recurse (const GlobalLinSysKey &mkey, const std::weak_ptr< ExpList > &pExpList, const DNekScalBlkMatSharedPtr pSchurCompl, const DNekScalBlkMatSharedPtr pBinvD, const DNekScalBlkMatSharedPtr pC, const DNekScalBlkMatSharedPtr pInvD, const std::shared_ptr< AssemblyMap > &locToGloMap)=0

void	v_Solve (const Array< OneD, const NekDouble > &in, Array< OneD, NekDouble > &out, const AssemblyMapSharedPtr &locToGloMap, const Array< OneD, const NekDouble > &dirForcing=NullNekDouble1DArray) override
	Solve the linear system for given input and output vectors using a specified local to global map. More...

void	v_InitObject () override

void	v_Initialise (const std::shared_ptr< AssemblyMap > &locToGloMap) override
	Initialise this object. More...

void	SetupTopLevel (const std::shared_ptr< AssemblyMap > &locToGloMap)
	Set up the storage for the Schur complement or the top level of the multi-level Schur complement. More...

void	ConstructNextLevelCondensedSystem (const std::shared_ptr< AssemblyMap > &locToGloMap)

Protected Member Functions inherited from Nektar::MultiRegions::GlobalLinSys
virtual void	v_Solve (const Array< OneD, const NekDouble > &in, Array< OneD, NekDouble > &out, const AssemblyMapSharedPtr &locToGloMap, const Array< OneD, const NekDouble > &dirForcing=NullNekDouble1DArray)=0
	Solve a linear system based on mapping. More...

virtual void	v_SolveLinearSystem (const int pNumRows, const Array< OneD, const NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const AssemblyMapSharedPtr &locToGloMap, const int pNumDir)=0
	Solve a basic matrix system. More...

virtual void	v_InitObject ()

virtual void	v_Initialise (const std::shared_ptr< AssemblyMap > &pLocToGloMap)

virtual int	v_GetNumBlocks ()
	Get the number of blocks in this system. More...

virtual DNekScalMatSharedPtr	v_GetBlock (unsigned int n)
	Retrieves the block matrix from n-th expansion using the matrix key provided by the m_linSysKey. More...

virtual void	v_DropBlock (unsigned int n)
	Releases the local block matrix from NekManager of n-th expansion using the matrix key provided by the m_linSysKey. More...

virtual DNekScalBlkMatSharedPtr	v_GetStaticCondBlock (unsigned int n)
	Retrieves a the static condensation block matrices from n-th expansion using the matrix key provided by the m_linSysKey. More...

virtual void	v_DropStaticCondBlock (unsigned int n)
	Releases the static condensation block matrices from NekManager of n-th expansion using the matrix key provided by the m_linSysKey. More...

PreconditionerSharedPtr	CreatePrecon (AssemblyMapSharedPtr asmMap)
	Create a preconditioner object from the parameters defined in the supplied assembly map. More...

Protected Attributes
GlobalLinSysStaticCondSharedPtr	m_recursiveSchurCompl
	Schur complement for Direct Static Condensation. More...

DNekScalBlkMatSharedPtr	m_schurCompl
	Block Schur complement matrix. More...

DNekScalBlkMatSharedPtr	m_BinvD
	Block \( BD^{-1} \) matrix. More...

DNekScalBlkMatSharedPtr	m_C
	Block \( C \) matrix. More...

DNekScalBlkMatSharedPtr	m_invD
	Block \( D^{-1} \) matrix. More...

std::weak_ptr< AssemblyMap >	m_locToGloMap
	Local to global map. More...

Array< OneD, NekDouble >	m_wsp
	Workspace array for matrix multiplication. More...

Array< OneD, const NekDouble >	m_sign

Protected Attributes inherited from Nektar::MultiRegions::GlobalLinSys
const GlobalLinSysKey	m_linSysKey
	Key associated with this linear system. More...

const std::weak_ptr< ExpList >	m_expList
	Local Matrix System. More...

const std::map< int, RobinBCInfoSharedPtr >	m_robinBCInfo
	Robin boundary info. More...

bool	m_verbose

Detailed Description

A global linear system.

Solves a linear system using single- or multi-level static condensation.

Definition at line 51 of file GlobalLinSysStaticCond.h.

Constructor & Destructor Documentation

◆ GlobalLinSysStaticCond()

Nektar::MultiRegions::GlobalLinSysStaticCond::GlobalLinSysStaticCond	(	const GlobalLinSysKey &	pKey,
		const std::weak_ptr< ExpList > &	pExpList,
		const std::shared_ptr< AssemblyMap > &	pLocToGloMap
	)

Constructor for full direct matrix solve.

For a matrix system of the form

\[ \left[ \begin{array}{cc} \boldsymbol{A} & \boldsymbol{B}\\ \boldsymbol{C} & \boldsymbol{D} \end{array} \right] \left[ \begin{array}{c} \boldsymbol{x_1}\\ \boldsymbol{x_2} \end{array}\right] = \left[ \begin{array}{c} \boldsymbol{y_1}\\ \boldsymbol{y_2} \end{array}\right], \]

where \(\boldsymbol{D}\) and \((\boldsymbol{A-BD^{-1}C})\) are invertible, store and assemble a static condensation system, according to a given local to global mapping. #m_linSys is constructed by AssembleSchurComplement().

Parameters

mKey	Associated matrix key.
pLocMatSys	LocalMatrixSystem
locToGloMap	Local to global mapping.

Definition at line 72 of file GlobalLinSysStaticCond.cpp.

    : GlobalLinSys(pKey, pExpList, pLocToGloMap), m_locToGloMap(pLocToGloMap)
{
}

◆ ~GlobalLinSysStaticCond()

Nektar::MultiRegions::GlobalLinSysStaticCond::~GlobalLinSysStaticCond ( )

override

Definition at line 91 of file GlobalLinSysStaticCond.cpp.

92{

93}

Member Function Documentation

◆ ConstructNextLevelCondensedSystem()

void Nektar::MultiRegions::GlobalLinSysStaticCond::ConstructNextLevelCondensedSystem ( const std::shared_ptr< AssemblyMap > & locToGloMap )

protected

Definition at line 311 of file GlobalLinSysStaticCond.cpp.

{
    int i, j, n, cnt;
    DNekScalBlkMatSharedPtr blkMatrices[4];
 
    // Create temporary matrices within an inner-local scope to ensure
    // any references to the intermediate storage is lost before
    // the recursive step, rather than at the end of the routine.
    // This allows the schur complement matrix from this level to be
    // disposed of in the next level after use without being retained
    // due to lingering shared pointers.
    {
 
        const Array<OneD, const unsigned int> &nBndDofsPerPatch =
            pLocToGloMap->GetNumLocalBndCoeffsPerPatch();
        const Array<OneD, const unsigned int> &nIntDofsPerPatch =
            pLocToGloMap->GetNumLocalIntCoeffsPerPatch();
 
        // STEP 1:
        // Based upon the schur complement of the the current level we
        // will substructure this matrix in the form
        //      --     --
        //      | A   B |
        //      | C   D |
        //      --     --
        // All matrices A,B,C and D are (diagonal) blockmatrices.
        // However, as a start we will use an array of DNekMatrices as
        // it is too hard to change the individual entries of a
        // DNekScalBlkMatSharedPtr.
 
        // In addition, we will also try to ensure that the memory of
        // the blockmatrices will be contiguous. This will probably
        // enhance the efficiency
        // - Calculate the total number of entries in the blockmatrices
        int nPatches  = pLocToGloMap->GetNumPatches();
        int nEntriesA = 0;
        int nEntriesB = 0;
        int nEntriesC = 0;
        int nEntriesD = 0;
 
        for (i = 0; i < nPatches; i++)
        {
            nEntriesA += nBndDofsPerPatch[i] * nBndDofsPerPatch[i];
            nEntriesB += nBndDofsPerPatch[i] * nIntDofsPerPatch[i];
            nEntriesC += nIntDofsPerPatch[i] * nBndDofsPerPatch[i];
            nEntriesD += nIntDofsPerPatch[i] * nIntDofsPerPatch[i];
        }
 
        // Now create the DNekMatrices and link them to the memory
        // allocated above
        Array<OneD, DNekMatSharedPtr> substructuredMat[4] = {
            Array<OneD, DNekMatSharedPtr>(nPatches),  // Matrix A
            Array<OneD, DNekMatSharedPtr>(nPatches),  // Matrix B
            Array<OneD, DNekMatSharedPtr>(nPatches),  // Matrix C
            Array<OneD, DNekMatSharedPtr>(nPatches)}; // Matrix D
 
        // Initialise storage for the matrices. We do this separately
        // for each matrix so the matrices may be independently
        // deallocated when no longer required.
        Array<OneD, NekDouble> storageA(nEntriesA, 0.0);
        Array<OneD, NekDouble> storageB(nEntriesB, 0.0);
        Array<OneD, NekDouble> storageC(nEntriesC, 0.0);
        Array<OneD, NekDouble> storageD(nEntriesD, 0.0);
 
        Array<OneD, NekDouble> tmparray;
        PointerWrapper wType = eWrapper;
        int cntA             = 0;
        int cntB             = 0;
        int cntC             = 0;
        int cntD             = 0;
 
        // Use symmetric storage for invD if possible
        MatrixStorage storageTypeD = eFULL;
        if ((m_linSysKey.GetMatrixType() == StdRegions::eMass) ||
            (m_linSysKey.GetMatrixType() == StdRegions::eLaplacian) ||
            (m_linSysKey.GetMatrixType() == StdRegions::eHelmholtz))
        {
            storageTypeD = eSYMMETRIC;
        }
 
        for (i = 0; i < nPatches; i++)
        {
            // Matrix A
            tmparray               = storageA + cntA;
            substructuredMat[0][i] = MemoryManager<DNekMat>::AllocateSharedPtr(
                nBndDofsPerPatch[i], nBndDofsPerPatch[i], tmparray, wType);
            // Matrix B
            tmparray               = storageB + cntB;
            substructuredMat[1][i] = MemoryManager<DNekMat>::AllocateSharedPtr(
                nBndDofsPerPatch[i], nIntDofsPerPatch[i], tmparray, wType);
            // Matrix C
            tmparray               = storageC + cntC;
            substructuredMat[2][i] = MemoryManager<DNekMat>::AllocateSharedPtr(
                nIntDofsPerPatch[i], nBndDofsPerPatch[i], tmparray, wType);
            // Matrix D
            tmparray               = storageD + cntD;
            substructuredMat[3][i] = MemoryManager<DNekMat>::AllocateSharedPtr(
                nIntDofsPerPatch[i], nIntDofsPerPatch[i], tmparray, wType,
                storageTypeD);
 
            cntA += nBndDofsPerPatch[i] * nBndDofsPerPatch[i];
            cntB += nBndDofsPerPatch[i] * nIntDofsPerPatch[i];
            cntC += nIntDofsPerPatch[i] * nBndDofsPerPatch[i];
            cntD += nIntDofsPerPatch[i] * nIntDofsPerPatch[i];
        }
 
        // Then, project SchurComplement onto
        // the substructured matrices of the next level
        DNekScalBlkMatSharedPtr SchurCompl = m_schurCompl;
        DNekScalMatSharedPtr schurComplSubMat;
        int schurComplSubMatnRows;
        Array<OneD, const int> patchId, dofId;
        Array<OneD, const unsigned int> isBndDof;
        Array<OneD, const NekDouble> sign;
 
        for (n = cnt = 0; n < SchurCompl->GetNumberOfBlockRows(); ++n)
        {
            schurComplSubMat      = SchurCompl->GetBlock(n, n);
            schurComplSubMatnRows = schurComplSubMat->GetRows();
 
            patchId =
                pLocToGloMap->GetPatchMapFromPrevLevel()->GetPatchId() + cnt;
            dofId = pLocToGloMap->GetPatchMapFromPrevLevel()->GetDofId() + cnt;
            isBndDof =
                pLocToGloMap->GetPatchMapFromPrevLevel()->IsBndDof() + cnt;
            sign = pLocToGloMap->GetPatchMapFromPrevLevel()->GetSign() + cnt;
 
            // Set up  Matrix;
            for (i = 0; i < schurComplSubMatnRows; ++i)
            {
                int pId = patchId[i];
                Array<OneD, NekDouble> subMat0 =
                    substructuredMat[0][pId]->GetPtr();
                Array<OneD, NekDouble> subMat1 =
                    substructuredMat[1][patchId[i]]->GetPtr();
                Array<OneD, NekDouble> subMat2 =
                    substructuredMat[2][patchId[i]]->GetPtr();
                DNekMatSharedPtr subMat3 = substructuredMat[3][patchId[i]];
                int subMat0rows          = substructuredMat[0][pId]->GetRows();
                int subMat1rows          = substructuredMat[1][pId]->GetRows();
                int subMat2rows          = substructuredMat[2][pId]->GetRows();
 
                if (isBndDof[i])
                {
                    for (j = 0; j < schurComplSubMatnRows; ++j)
                    {
                        ASSERTL0(patchId[i] == patchId[j],
                                 "These values should be equal");
 
                        if (isBndDof[j])
                        {
                            subMat0[dofId[i] + dofId[j] * subMat0rows] +=
                                sign[i] * sign[j] * (*schurComplSubMat)(i, j);
                        }
                        else
                        {
                            subMat1[dofId[i] + dofId[j] * subMat1rows] +=
                                sign[i] * sign[j] * (*schurComplSubMat)(i, j);
                        }
                    }
                }
                else
                {
                    for (j = 0; j < schurComplSubMatnRows; ++j)
                    {
                        ASSERTL0(patchId[i] == patchId[j],
                                 "These values should be equal");
 
                        if (isBndDof[j])
                        {
                            subMat2[dofId[i] + dofId[j] * subMat2rows] +=
                                sign[i] * sign[j] * (*schurComplSubMat)(i, j);
                        }
                        else
                        {
                            if (storageTypeD == eSYMMETRIC)
                            {
                                if (dofId[i] <= dofId[j])
                                {
                                    (*subMat3)(dofId[i], dofId[j]) +=
                                        sign[i] * sign[j] *
                                        (*schurComplSubMat)(i, j);
                                }
                            }
                            else
                            {
                                (*subMat3)(dofId[i], dofId[j]) +=
                                    sign[i] * sign[j] *
                                    (*schurComplSubMat)(i, j);
                            }
                        }
                    }
                }
            }
            cnt += schurComplSubMatnRows;
        }
 
        // STEP 2: condense the system
        // This can be done elementally (i.e. patch per patch)
        for (i = 0; i < nPatches; i++)
        {
            if (nIntDofsPerPatch[i])
            {
                Array<OneD, NekDouble> subMat0 =
                    substructuredMat[0][i]->GetPtr();
                Array<OneD, NekDouble> subMat1 =
                    substructuredMat[1][i]->GetPtr();
                Array<OneD, NekDouble> subMat2 =
                    substructuredMat[2][i]->GetPtr();
                int subMat0rows = substructuredMat[0][i]->GetRows();
                int subMat1rows = substructuredMat[1][i]->GetRows();
                int subMat2rows = substructuredMat[2][i]->GetRows();
                int subMat2cols = substructuredMat[2][i]->GetColumns();
 
                // 1. D -> InvD
                substructuredMat[3][i]->Invert();
                // 2. B -> BInvD
                (*substructuredMat[1][i]) =
                    (*substructuredMat[1][i]) * (*substructuredMat[3][i]);
                // 3. A -> A - BInvD*C (= schurcomplement)
                // (*substructuredMat[0][i]) =(*substructuredMat[0][i])-
                // (*substructuredMat[1][i])*(*substructuredMat[2][i]);
                // Note: faster to use blas directly
                if (subMat1rows && subMat2cols)
                {
                    Blas::Dgemm('N', 'N', subMat1rows, subMat2cols, subMat2rows,
                                -1.0, &subMat1[0], subMat1rows, &subMat2[0],
                                subMat2rows, 1.0, &subMat0[0], subMat0rows);
                }
            }
        }
 
        // STEP 3: fill the block matrices. However, do note that we
        // first have to convert them to a DNekScalMat in order to be
        // compatible with the first level of static condensation
        const Array<OneD, const unsigned int> &nbdry_size =
            pLocToGloMap->GetNumLocalBndCoeffsPerPatch();
        const Array<OneD, const unsigned int> &nint_size =
            pLocToGloMap->GetNumLocalIntCoeffsPerPatch();
        MatrixStorage blkmatStorage = eDIAGONAL;
 
        blkMatrices[0] = MemoryManager<DNekScalBlkMat>::AllocateSharedPtr(
            nbdry_size, nbdry_size, blkmatStorage);
        blkMatrices[1] = MemoryManager<DNekScalBlkMat>::AllocateSharedPtr(
            nbdry_size, nint_size, blkmatStorage);
        blkMatrices[2] = MemoryManager<DNekScalBlkMat>::AllocateSharedPtr(
            nint_size, nbdry_size, blkmatStorage);
        blkMatrices[3] = MemoryManager<DNekScalBlkMat>::AllocateSharedPtr(
            nint_size, nint_size, blkmatStorage);
 
        DNekScalMatSharedPtr tmpscalmat;
        for (i = 0; i < nPatches; i++)
        {
            for (j = 0; j < 4; j++)
            {
                tmpscalmat = MemoryManager<DNekScalMat>::AllocateSharedPtr(
                    1.0, substructuredMat[j][i]);
                blkMatrices[j]->SetBlock(i, i, tmpscalmat);
            }
        }
    }
 
    // We've finished with the Schur complement matrix passed to this
    // level, so return the memory to the system. The Schur complement
    // matrix need only be retained at the last level. Save the other
    // matrices at this level though.
    m_schurCompl.reset();
 
    m_recursiveSchurCompl =
        v_Recurse(m_linSysKey, m_expList, blkMatrices[0], blkMatrices[1],
                  blkMatrices[2], blkMatrices[3], pLocToGloMap);
}

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL0, Blas::Dgemm(), Nektar::eDIAGONAL, Nektar::eFULL, Nektar::StdRegions::eHelmholtz, Nektar::StdRegions::eLaplacian, Nektar::StdRegions::eMass, Nektar::eSYMMETRIC, Nektar::eWrapper, Nektar::MultiRegions::GlobalMatrixKey::GetMatrixType(), Nektar::MultiRegions::GlobalLinSys::m_expList, Nektar::MultiRegions::GlobalLinSys::m_linSysKey, m_recursiveSchurCompl, m_schurCompl, sign, and v_Recurse().

Referenced by v_Initialise().

◆ SetupTopLevel()

void Nektar::MultiRegions::GlobalLinSysStaticCond::SetupTopLevel ( const std::shared_ptr< AssemblyMap > & pLocToGloMap )

protected

Set up the storage for the Schur complement or the top level of the multi-level Schur complement.

For the first level in multi-level static condensation, or the only level in the case of single-level static condensation, allocate the condensed matrices and populate them with the local matrices retrieved from the expansion list.

Parameters

Definition at line 266 of file GlobalLinSysStaticCond.cpp.

{
    int n;
    int n_exp = m_expList.lock()->GetNumElmts();
 
    const Array<OneD, const unsigned int> &nbdry_size =
        pLocToGloMap->GetNumLocalBndCoeffsPerPatch();
    const Array<OneD, const unsigned int> &nint_size =
        pLocToGloMap->GetNumLocalIntCoeffsPerPatch();
 
    // Setup Block Matrix systems
    MatrixStorage blkmatStorage = eDIAGONAL;
    m_schurCompl = MemoryManager<DNekScalBlkMat>::AllocateSharedPtr(
        nbdry_size, nbdry_size, blkmatStorage);
    m_BinvD = MemoryManager<DNekScalBlkMat>::AllocateSharedPtr(
        nbdry_size, nint_size, blkmatStorage);
    m_C = MemoryManager<DNekScalBlkMat>::AllocateSharedPtr(
        nint_size, nbdry_size, blkmatStorage);
    m_invD = MemoryManager<DNekScalBlkMat>::AllocateSharedPtr(
        nint_size, nint_size, blkmatStorage);
 
    for (n = 0; n < n_exp; ++n)
    {
        if (m_linSysKey.GetMatrixType() == StdRegions::eHybridDGHelmBndLam)
        {
            DNekScalMatSharedPtr loc_mat = GlobalLinSys::v_GetBlock(n);
            m_schurCompl->SetBlock(n, n, loc_mat);
        }
        else
        {
            DNekScalBlkMatSharedPtr loc_schur =
                GlobalLinSys::v_GetStaticCondBlock(n);
            DNekScalMatSharedPtr t;
            m_schurCompl->SetBlock(n, n, t = loc_schur->GetBlock(0, 0));
            m_BinvD->SetBlock(n, n, t = loc_schur->GetBlock(0, 1));
            m_C->SetBlock(n, n, t = loc_schur->GetBlock(1, 0));
            m_invD->SetBlock(n, n, t = loc_schur->GetBlock(1, 1));
        }
    }
}

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::eDIAGONAL, Nektar::StdRegions::eHybridDGHelmBndLam, Nektar::MultiRegions::GlobalMatrixKey::GetMatrixType(), m_BinvD, m_C, Nektar::MultiRegions::GlobalLinSys::m_expList, m_invD, Nektar::MultiRegions::GlobalLinSys::m_linSysKey, m_schurCompl, Nektar::MultiRegions::GlobalLinSys::v_GetBlock(), and Nektar::MultiRegions::GlobalLinSys::v_GetStaticCondBlock().

Referenced by Nektar::MultiRegions::GlobalLinSysIterativeStaticCond::v_InitObject(), Nektar::MultiRegions::GlobalLinSysPETScStaticCond::v_InitObject(), and v_InitObject().

◆ v_AssembleSchurComplement()

virtual void Nektar::MultiRegions::GlobalLinSysStaticCond::v_AssembleSchurComplement ( std::shared_ptr< AssemblyMap > pLoctoGloMap )

inlineprotectedvirtual

Reimplemented in Nektar::MultiRegions::GlobalLinSysIterativeStaticCond, Nektar::MultiRegions::GlobalLinSysPETScStaticCond, Nektar::MultiRegions::GlobalLinSysXxtStaticCond, and Nektar::MultiRegions::GlobalLinSysDirectStaticCond.

Definition at line 83 of file GlobalLinSysStaticCond.h.

85 {

86 }

Referenced by v_Initialise().

◆ v_BasisFwdTransform()

virtual void Nektar::MultiRegions::GlobalLinSysStaticCond::v_BasisFwdTransform ( Array< OneD, NekDouble > & pInOut )

inlineprotectedvirtual

Reimplemented in Nektar::MultiRegions::GlobalLinSysIterativeStaticCond, and Nektar::MultiRegions::GlobalLinSysPETScStaticCond.

Definition at line 67 of file GlobalLinSysStaticCond.h.

69 {

70 }

Referenced by v_Solve().

◆ v_CoeffsBwdTransform()

virtual void Nektar::MultiRegions::GlobalLinSysStaticCond::v_CoeffsBwdTransform ( Array< OneD, NekDouble > & pInOut )

inlineprotectedvirtual

Reimplemented in Nektar::MultiRegions::GlobalLinSysIterativeStaticCond, and Nektar::MultiRegions::GlobalLinSysPETScStaticCond.

Definition at line 72 of file GlobalLinSysStaticCond.h.

74 {

75 }

Referenced by v_Solve().

◆ v_CoeffsFwdTransform()

virtual void Nektar::MultiRegions::GlobalLinSysStaticCond::v_CoeffsFwdTransform	(	const Array< OneD, NekDouble > &	pInput,
		Array< OneD, NekDouble > &	pOutput
	)

inlineprotectedvirtual

Reimplemented in Nektar::MultiRegions::GlobalLinSysIterativeStaticCond, and Nektar::MultiRegions::GlobalLinSysPETScStaticCond.

Definition at line 77 of file GlobalLinSysStaticCond.h.

80 {

81 }

Referenced by v_Solve().

◆ v_GetNumBlocks()

int Nektar::MultiRegions::GlobalLinSysStaticCond::v_GetNumBlocks ( )

overrideprotectedvirtual

Get the number of blocks in this system.

At the top level this corresponds to the number of elements in the expansion list.

Reimplemented from Nektar::MultiRegions::GlobalLinSys.

Definition at line 254 of file GlobalLinSysStaticCond.cpp.

{
    return m_schurCompl->GetNumberOfBlockRows();
}

References m_schurCompl.

◆ v_Initialise()

void Nektar::MultiRegions::GlobalLinSysStaticCond::v_Initialise ( const std::shared_ptr< AssemblyMap > & pLocToGloMap )

overrideprotectedvirtual

Initialise this object.

If at the last level of recursion (or the only level in the case of single-level static condensation), assemble the Schur complement. For other levels, in the case of multi-level static condensation, the next level of the condensed system is computed.

Parameters

pLocToGloMap Local to global mapping.

Reimplemented from Nektar::MultiRegions::GlobalLinSys.

Definition at line 236 of file GlobalLinSysStaticCond.cpp.

{
    int nLocalBnd = m_locToGloMap.lock()->GetNumLocalBndCoeffs();
    int nIntDofs  = m_locToGloMap.lock()->GetNumLocalCoeffs() - nLocalBnd;
    m_wsp         = Array<OneD, NekDouble>(4 * nLocalBnd + nIntDofs, 0.0);
 
    if (pLocToGloMap->AtLastLevel())
    {
        v_AssembleSchurComplement(pLocToGloMap);
    }
    else
    {
        ConstructNextLevelCondensedSystem(
            pLocToGloMap->GetNextLevelLocalToGlobalMap());
    }
}

References ConstructNextLevelCondensedSystem(), m_locToGloMap, m_wsp, and v_AssembleSchurComplement().

◆ v_InitObject()

void Nektar::MultiRegions::GlobalLinSysStaticCond::v_InitObject ( )

overrideprotectedvirtual

Reimplemented from Nektar::MultiRegions::GlobalLinSys.

Definition at line 79 of file GlobalLinSysStaticCond.cpp.

{
    // Allocate memory for top-level structure
    SetupTopLevel(m_locToGloMap.lock());
 
    // Construct this level
    Initialise(m_locToGloMap.lock());
}

References Nektar::MultiRegions::GlobalLinSys::Initialise(), m_locToGloMap, and SetupTopLevel().

◆ v_PreSolve()

virtual void Nektar::MultiRegions::GlobalLinSysStaticCond::v_PreSolve	(	int	scLevel,
		Array< OneD, NekDouble > &	F_bnd
	)

inlineprotectedvirtual

Reimplemented in Nektar::MultiRegions::GlobalLinSysPETScStaticCond, and Nektar::MultiRegions::GlobalLinSysIterativeStaticCond.

Definition at line 62 of file GlobalLinSysStaticCond.h.

64 {

65 }

Referenced by v_Solve().

◆ v_Recurse()

virtual GlobalLinSysStaticCondSharedPtr Nektar::MultiRegions::GlobalLinSysStaticCond::v_Recurse	(	const GlobalLinSysKey &	mkey,
		const std::weak_ptr< ExpList > &	pExpList,
		const DNekScalBlkMatSharedPtr	pSchurCompl,
		const DNekScalBlkMatSharedPtr	pBinvD,
		const DNekScalBlkMatSharedPtr	pC,
		const DNekScalBlkMatSharedPtr	pInvD,
		const std::shared_ptr< AssemblyMap > &	locToGloMap
	)

protectedpure virtual

Implemented in Nektar::MultiRegions::GlobalLinSysDirectStaticCond, Nektar::MultiRegions::GlobalLinSysIterativeStaticCond, Nektar::MultiRegions::GlobalLinSysPETScStaticCond, and Nektar::MultiRegions::GlobalLinSysXxtStaticCond.

Referenced by ConstructNextLevelCondensedSystem().

◆ v_Solve()

void Nektar::MultiRegions::GlobalLinSysStaticCond::v_Solve	(	const Array< OneD, const NekDouble > &	in,
		Array< OneD, NekDouble > &	out,
		const AssemblyMapSharedPtr &	locToGloMap,
		const Array< OneD, const NekDouble > &	dirForcing = `NullNekDouble1DArray`
	)

overrideprotectedvirtual

Solve the linear system for given input and output vectors using a specified local to global map.

Implements Nektar::MultiRegions::GlobalLinSys.

Definition at line 98 of file GlobalLinSysStaticCond.cpp.

{
    ASSERTL1(dirForcing.size() == 0,
             "GlobalLinSysStaticCond: Not setup for dirForcing");
 
    bool atLastLevel = pLocToGloMap->AtLastLevel();
    int scLevel      = pLocToGloMap->GetStaticCondLevel();
 
    int nGlobDofs    = pLocToGloMap->GetNumGlobalCoeffs();
    int nLocBndDofs  = pLocToGloMap->GetNumLocalBndCoeffs();
    int nGlobBndDofs = pLocToGloMap->GetNumGlobalBndCoeffs();
    int nDirBndDofs  = pLocToGloMap->GetNumGlobalDirBndCoeffs();
    int nIntDofs     = nGlobDofs - nGlobBndDofs;
 
    if ((nGlobDofs - nDirBndDofs) == 0)
    {
        return; // nothing to solve;
    }
 
    Array<OneD, NekDouble> F_bnd, F_bnd1, F_int, V_bnd;
    Array<OneD, NekDouble> tmp, unused;
 
    F_bnd = m_wsp;
    // unsued = Temporary storage required for DoMatrixMultiply, not used here
    unused = m_wsp + 1 * nLocBndDofs;
    F_bnd1 = m_wsp + 2 * nLocBndDofs;
    V_bnd  = m_wsp + 3 * nLocBndDofs;
    F_int  = m_wsp + 4 * nLocBndDofs;
 
    pLocToGloMap->LocalToLocalBnd(pLocOutput, V_bnd);
 
    NekVector<NekDouble> F_Int(nIntDofs, F_int, eWrapper);
    NekVector<NekDouble> V_Bnd(nLocBndDofs, V_bnd, eWrapper);
 
    if (nIntDofs)
    {
        m_locToGloMap.lock()->LocalToLocalInt(pLocInput, F_int);
    }
 
    // Boundary system solution
    if (nGlobBndDofs - nDirBndDofs)
    {
        pLocToGloMap->LocalToLocalBnd(pLocInput, F_bnd);
 
        // set up normalisation factor for right hand side on first SC level
        v_PreSolve(scLevel, F_bnd);
 
        // Gather boundary expansison into locbnd
        NekVector<NekDouble> F_Bnd(nLocBndDofs, F_bnd1, eWrapper);
 
        // construct boundary forcing
        if (nIntDofs)
        {
            DNekScalBlkMat &BinvD = *m_BinvD;
 
            F_Bnd = BinvD * F_Int;
 
            Vmath::Vsub(nLocBndDofs, F_bnd, 1, F_bnd1, 1, F_bnd, 1);
        }
 
        if (atLastLevel)
        {
            // Transform to new basis if required
            v_BasisFwdTransform(F_bnd);
 
            DNekScalBlkMat &SchurCompl = *m_schurCompl;
 
            v_CoeffsFwdTransform(V_bnd, V_bnd);
 
            // subtract dirichlet boundary forcing
            F_Bnd = SchurCompl * V_Bnd;
 
            Vmath::Vsub(nLocBndDofs, F_bnd, 1, F_bnd1, 1, F_bnd, 1);
 
            Array<OneD, NekDouble> F_hom;
 
            // Solve for difference from initial solution given inout;
            SolveLinearSystem(nGlobBndDofs, F_bnd, F_bnd1, pLocToGloMap,
                              nDirBndDofs);
 
            // Add back initial conditions onto difference
            Vmath::Vadd(nLocBndDofs, V_bnd, 1, F_bnd1, 1, V_bnd, 1);
 
            // Transform back to original basis
            v_CoeffsBwdTransform(V_bnd);
 
            // put final bnd solution back in output array
            m_locToGloMap.lock()->LocalBndToLocal(V_bnd, pLocOutput);
        }
        else // Process next level of recursion for multi level SC
        {
            AssemblyMapSharedPtr nextLevLocToGloMap =
                pLocToGloMap->GetNextLevelLocalToGlobalMap();
 
            // partially assemble F for next level and
            // reshuffle V_bnd
            nextLevLocToGloMap->PatchAssemble(F_bnd, F_bnd);
            nextLevLocToGloMap->PatchLocalToGlobal(V_bnd, V_bnd);
 
            m_recursiveSchurCompl->Solve(F_bnd, V_bnd, nextLevLocToGloMap);
 
            // unpack V_bnd
            nextLevLocToGloMap->PatchGlobalToLocal(V_bnd, V_bnd);
 
            // place V_bnd in output array
            m_locToGloMap.lock()->LocalBndToLocal(V_bnd, pLocOutput);
        }
    }
 
    // solve interior system
    if (nIntDofs)
    {
        Array<OneD, NekDouble> V_int(nIntDofs);
        NekVector<NekDouble> V_Int(nIntDofs, V_int, eWrapper);
 
        // get array of local solutions
        DNekScalBlkMat &invD = *m_invD;
        DNekScalBlkMat &C    = *m_C;
 
        F_Int = F_Int - C * V_Bnd;
 
        Multiply(V_Int, invD, F_Int);
 
        m_locToGloMap.lock()->LocalIntToLocal(V_int, pLocOutput);
    }
}

References ASSERTL1, Nektar::eWrapper, m_BinvD, m_C, m_invD, m_locToGloMap, m_recursiveSchurCompl, m_schurCompl, m_wsp, Nektar::Multiply(), Nektar::MultiRegions::GlobalLinSys::SolveLinearSystem(), v_BasisFwdTransform(), v_CoeffsBwdTransform(), v_CoeffsFwdTransform(), v_PreSolve(), Vmath::Vadd(), and Vmath::Vsub().

Member Data Documentation

◆ m_BinvD

DNekScalBlkMatSharedPtr Nektar::MultiRegions::GlobalLinSysStaticCond::m_BinvD

protected

Block \( BD^{-1} \) matrix.

Definition at line 102 of file GlobalLinSysStaticCond.h.

Referenced by Nektar::MultiRegions::GlobalLinSysDirectStaticCond::GlobalLinSysDirectStaticCond(), Nektar::MultiRegions::GlobalLinSysIterativeStaticCond::GlobalLinSysIterativeStaticCond(), Nektar::MultiRegions::GlobalLinSysPETScStaticCond::GlobalLinSysPETScStaticCond(), Nektar::MultiRegions::GlobalLinSysXxtStaticCond::GlobalLinSysXxtStaticCond(), SetupTopLevel(), Nektar::MultiRegions::GlobalLinSysPETScStaticCond::v_AssembleSchurComplement(), Nektar::MultiRegions::GlobalLinSysDirectStaticCond::v_AssembleSchurComplement(), and v_Solve().

◆ m_C

DNekScalBlkMatSharedPtr Nektar::MultiRegions::GlobalLinSysStaticCond::m_C

protected

Block \( C \) matrix.

Definition at line 104 of file GlobalLinSysStaticCond.h.

Referenced by Nektar::MultiRegions::GlobalLinSysDirectStaticCond::GlobalLinSysDirectStaticCond(), Nektar::MultiRegions::GlobalLinSysIterativeStaticCond::GlobalLinSysIterativeStaticCond(), Nektar::MultiRegions::GlobalLinSysPETScStaticCond::GlobalLinSysPETScStaticCond(), Nektar::MultiRegions::GlobalLinSysXxtStaticCond::GlobalLinSysXxtStaticCond(), SetupTopLevel(), Nektar::MultiRegions::GlobalLinSysPETScStaticCond::v_AssembleSchurComplement(), Nektar::MultiRegions::GlobalLinSysDirectStaticCond::v_AssembleSchurComplement(), and v_Solve().

◆ m_invD

DNekScalBlkMatSharedPtr Nektar::MultiRegions::GlobalLinSysStaticCond::m_invD

protected

Block \( D^{-1} \) matrix.

Definition at line 106 of file GlobalLinSysStaticCond.h.

Referenced by Nektar::MultiRegions::GlobalLinSysDirectStaticCond::GlobalLinSysDirectStaticCond(), Nektar::MultiRegions::GlobalLinSysIterativeStaticCond::GlobalLinSysIterativeStaticCond(), Nektar::MultiRegions::GlobalLinSysPETScStaticCond::GlobalLinSysPETScStaticCond(), Nektar::MultiRegions::GlobalLinSysXxtStaticCond::GlobalLinSysXxtStaticCond(), SetupTopLevel(), Nektar::MultiRegions::GlobalLinSysPETScStaticCond::v_AssembleSchurComplement(), Nektar::MultiRegions::GlobalLinSysDirectStaticCond::v_AssembleSchurComplement(), and v_Solve().

◆ m_locToGloMap

std::weak_ptr<AssemblyMap> Nektar::MultiRegions::GlobalLinSysStaticCond::m_locToGloMap

protected

Local to global map.

Definition at line 108 of file GlobalLinSysStaticCond.h.

◆ m_recursiveSchurCompl

GlobalLinSysStaticCondSharedPtr Nektar::MultiRegions::GlobalLinSysStaticCond::m_recursiveSchurCompl

protected

Schur complement for Direct Static Condensation.

Definition at line 98 of file GlobalLinSysStaticCond.h.

Referenced by ConstructNextLevelCondensedSystem(), and v_Solve().

◆ m_schurCompl

DNekScalBlkMatSharedPtr Nektar::MultiRegions::GlobalLinSysStaticCond::m_schurCompl

protected

Block Schur complement matrix.

Definition at line 100 of file GlobalLinSysStaticCond.h.

◆ m_sign

Array<OneD, const NekDouble> Nektar::MultiRegions::GlobalLinSysStaticCond::m_sign

protected

Definition at line 112 of file GlobalLinSysStaticCond.h.

◆ m_wsp

Array<OneD, NekDouble> Nektar::MultiRegions::GlobalLinSysStaticCond::m_wsp

protected

Workspace array for matrix multiplication.

Definition at line 110 of file GlobalLinSysStaticCond.h.

Referenced by Nektar::MultiRegions::GlobalLinSysIterativeStaticCond::v_DoMatrixMultiply(), v_Initialise(), and v_Solve().

Public Member Functions

Protected Member Functions

Protected Attributes

Detailed Description

Constructor & Destructor Documentation

◆ GlobalLinSysStaticCond()

◆ ~GlobalLinSysStaticCond()

Member Function Documentation

◆ ConstructNextLevelCondensedSystem()

◆ SetupTopLevel()

◆ v_AssembleSchurComplement()

◆ v_BasisFwdTransform()

◆ v_CoeffsBwdTransform()

◆ v_CoeffsFwdTransform()

◆ v_GetNumBlocks()

◆ v_Initialise()

◆ v_InitObject()

◆ v_PreSolve()

◆ v_Recurse()

◆ v_Solve()

Member Data Documentation

◆ m_BinvD

◆ m_C

◆ m_invD

◆ m_locToGloMap

◆ m_recursiveSchurCompl

◆ m_schurCompl

◆ m_sign

◆ m_wsp