A global linear system. More...

#include <GlobalLinSysIterativeStaticCond.h>

Inheritance diagram for Nektar::MultiRegions::GlobalLinSysIterativeStaticCond:

Public Types
typedef NekSparseDiagBlkMatrix< StorageSmvBsr< NekDouble > >	DNekSmvBsrDiagBlkMat

typedef std::shared_ptr< DNekSmvBsrDiagBlkMat >	DNekSmvBsrDiagBlkMatSharedPtr

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

	GlobalLinSysIterativeStaticCond (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, const PreconditionerSharedPtr pPrecon)
	Constructor for full direct matrix solve. More...

	~GlobalLinSysIterativeStaticCond () override

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

	~GlobalLinSysIterative () override

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

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

Public Member Functions inherited from Nektar::MultiRegions::GlobalLinSysStaticCond
	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

Static Public Member Functions
static GlobalLinSysSharedPtr	create (const GlobalLinSysKey &pLinSysKey, const std::weak_ptr< ExpList > &pExpList, const std::shared_ptr< AssemblyMap > &pLocToGloMap)
	Creates an instance of this class. More...

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

static std::string	className2

Protected Member Functions
void	v_InitObject () override

void	v_AssembleSchurComplement (const std::shared_ptr< AssemblyMap > locToGloMap) override
	Assemble the Schur complement matrix. More...

void	v_DoMatrixMultiply (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput) override
	Perform a Shur-complement matrix multiply operation. More...

void	v_UniqueMap () override

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

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) override

void	v_PreSolve (int scLevel, Array< OneD, NekDouble > &F_bnd) override

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

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

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

Protected Member Functions inherited from Nektar::MultiRegions::GlobalLinSysIterative
void	DoProjection (const int pNumRows, const Array< OneD, const NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const int pNumDir, const NekDouble tol, const bool isAconjugate)
	projection technique More...

void	Set_Rhs_Magnitude (const NekVector< NekDouble > &pIn)

virtual void	v_UniqueMap ()=0

virtual void	v_DoMatrixMultiply (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)=0

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 Member Functions inherited from Nektar::MultiRegions::GlobalLinSysStaticCond
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)

Private Member Functions
void	PrepareLocalSchurComplement ()
	Prepares local representation of Schur complement stored as a sparse block-diagonal matrix. More...

void	v_SolveLinearSystem (const int pNumRows, const Array< OneD, const NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const AssemblyMapSharedPtr &locToGloMap, const int pNumDir) override
	Solve the matrix system. More...

Private Attributes
std::vector< double >	m_storage
	Dense storage for block Schur complement matrix. More...

std::vector< const double * >	m_denseBlocks
	Vector of pointers to local matrix data. More...

Array< OneD, unsigned int >	m_rows
	Ranks of local matrices. More...

Array< OneD, NekDouble >	m_scale
	Scaling factors for local matrices. More...

DNekSmvBsrDiagBlkMatSharedPtr	m_sparseSchurCompl
	Sparse representation of Schur complement matrix at this level. More...

Static Private Attributes
static std::string	storagedef
	Utility strings. More...

static std::string	storagelookupIds []

Additional Inherited Members
Protected Attributes inherited from Nektar::MultiRegions::GlobalLinSysIterative
Array< OneD, int >	m_map
	Global to universal unique map. More...

NekDouble	m_tolerance
	Tolerance of iterative solver. More...

NekDouble	m_rhs_magnitude
	dot product of rhs to normalise stopping criterion More...

NekDouble	m_rhs_mag_sm
	cnt to how many times rhs_magnitude is called More...

PreconditionerSharedPtr	m_precon

std::string	m_precontype

int	m_totalIterations

bool	m_useProjection
	Whether to apply projection technique. More...

bool	m_root
	Root if parallel. More...

std::string	m_linSysIterSolver
	Iterative solver: Conjugate Gradient, GMRES. More...

std::vector< Array< OneD, NekDouble > >	m_prevLinSol
	Storage for solutions to previous linear problems. More...

std::vector< Array< OneD, NekDouble > >	m_prevBasis

DNekMatSharedPtr	m_coeffMatrix

Array< OneD, NekDouble >	m_coeffMatrixFactor

Array< OneD, int >	m_ipivot

int	m_numSuccessiveRHS

bool	m_isAconjugate

std::string	m_matrixType

bool	m_isNonSymmetricLinSys

int	m_numPrevSols

bool	m_isAbsoluteTolerance

LibUtilities::NekSysOperators	m_NekSysOp

LibUtilities::NekLinSysIterSharedPtr	m_linsol

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

Protected Attributes inherited from Nektar::MultiRegions::GlobalLinSysStaticCond
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

Static Protected Attributes inherited from Nektar::MultiRegions::GlobalLinSysIterative
static std::string	IteratSolverlookupIds []

static std::string	IteratSolverdef

Detailed Description

A global linear system.

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

Definition at line 63 of file GlobalLinSysIterativeStaticCond.h.

Member Typedef Documentation

◆ DNekSmvBsrDiagBlkMat

typedef NekSparseDiagBlkMatrix<StorageSmvBsr<NekDouble> > Nektar::MultiRegions::GlobalLinSysIterativeStaticCond::DNekSmvBsrDiagBlkMat

Definition at line 68 of file GlobalLinSysIterativeStaticCond.h.

◆ DNekSmvBsrDiagBlkMatSharedPtr

typedef std::shared_ptr<DNekSmvBsrDiagBlkMat> Nektar::MultiRegions::GlobalLinSysIterativeStaticCond::DNekSmvBsrDiagBlkMatSharedPtr

Definition at line 69 of file GlobalLinSysIterativeStaticCond.h.

Constructor & Destructor Documentation

◆ GlobalLinSysIterativeStaticCond() [1/2]

Nektar::MultiRegions::GlobalLinSysIterativeStaticCond::GlobalLinSysIterativeStaticCond	(	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 99 of file GlobalLinSysIterativeStaticCond.cpp.

    : GlobalLinSys(pKey, pExpList, pLocToGloMap),
      GlobalLinSysIterative(pKey, pExpList, pLocToGloMap),
      GlobalLinSysStaticCond(pKey, pExpList, pLocToGloMap)
{
    ASSERTL1(
        (pKey.GetGlobalSysSolnType() == eIterativeStaticCond) ||
            (pKey.GetGlobalSysSolnType() == eIterativeMultiLevelStaticCond),
        "This constructor is only valid when using static "
        "condensation");
    ASSERTL1(pKey.GetGlobalSysSolnType() ==
                 pLocToGloMap->GetGlobalSysSolnType(),
             "The local to global map is not set up for the requested "
             "solution type");
}

References ASSERTL1, Nektar::MultiRegions::eIterativeMultiLevelStaticCond, Nektar::MultiRegions::eIterativeStaticCond, and Nektar::MultiRegions::GlobalLinSysKey::GetGlobalSysSolnType().

◆ GlobalLinSysIterativeStaticCond() [2/2]

Nektar::MultiRegions::GlobalLinSysIterativeStaticCond::GlobalLinSysIterativeStaticCond	(	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,
		const PreconditionerSharedPtr	pPrecon
	)

Constructor for full direct matrix solve.

Definition at line 120 of file GlobalLinSysIterativeStaticCond.cpp.

    : GlobalLinSys(pKey, pExpList, pLocToGloMap),
      GlobalLinSysIterative(pKey, pExpList, pLocToGloMap),
      GlobalLinSysStaticCond(pKey, pExpList, pLocToGloMap)
{
    m_schurCompl = pSchurCompl;
    m_BinvD      = pBinvD;
    m_C          = pC;
    m_invD       = pInvD;
    m_precon     = pPrecon;
}

References Nektar::MultiRegions::GlobalLinSysStaticCond::m_BinvD, Nektar::MultiRegions::GlobalLinSysStaticCond::m_C, Nektar::MultiRegions::GlobalLinSysStaticCond::m_invD, Nektar::MultiRegions::GlobalLinSysIterative::m_precon, and Nektar::MultiRegions::GlobalLinSysStaticCond::m_schurCompl.

◆ ~GlobalLinSysIterativeStaticCond()

Nektar::MultiRegions::GlobalLinSysIterativeStaticCond::~GlobalLinSysIterativeStaticCond ( )

override

Definition at line 174 of file GlobalLinSysIterativeStaticCond.cpp.

175{

176}

Member Function Documentation

◆ create()

static GlobalLinSysSharedPtr Nektar::MultiRegions::GlobalLinSysIterativeStaticCond::create	(	const GlobalLinSysKey &	pLinSysKey,
		const std::weak_ptr< ExpList > &	pExpList,
		const std::shared_ptr< AssemblyMap > &	pLocToGloMap
	)

inlinestatic

Creates an instance of this class.

Definition at line 72 of file GlobalLinSysIterativeStaticCond.h.

    {
        GlobalLinSysSharedPtr p =
            MemoryManager<GlobalLinSysIterativeStaticCond>::AllocateSharedPtr(
                pLinSysKey, pExpList, pLocToGloMap);
        p->InitObject();
        return p;
    }

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

◆ PrepareLocalSchurComplement()

void Nektar::MultiRegions::GlobalLinSysIterativeStaticCond::PrepareLocalSchurComplement ( )

private

Prepares local representation of Schur complement stored as a sparse block-diagonal matrix.

Populates sparse block-diagonal schur complement matrix from the block matrices stored in #m_blkMatrices.

Definition at line 220 of file GlobalLinSysIterativeStaticCond.cpp.

{
    LocalMatrixStorageStrategy storageStrategy =
        m_expList.lock()
            ->GetSession()
            ->GetSolverInfoAsEnum<LocalMatrixStorageStrategy>(
                "LocalMatrixStorageStrategy");
 
    switch (storageStrategy)
    {
        case MultiRegions::eContiguous:
        case MultiRegions::eNonContiguous:
        {
            size_t storageSize = 0;
            int nBlk           = m_schurCompl->GetNumberOfBlockRows();
 
            m_scale = Array<OneD, NekDouble>(nBlk, 1.0);
            m_rows  = Array<OneD, unsigned int>(nBlk, 0U);
 
            // Determine storage requirements for dense blocks.
            for (int i = 0; i < nBlk; ++i)
            {
                m_rows[i]  = m_schurCompl->GetBlock(i, i)->GetRows();
                m_scale[i] = m_schurCompl->GetBlock(i, i)->Scale();
                storageSize += m_rows[i] * m_rows[i];
            }
 
            // Assemble dense storage blocks.
            DNekScalMatSharedPtr loc_mat;
            m_denseBlocks.resize(nBlk);
            double *ptr = nullptr;
 
            if (MultiRegions::eContiguous == storageStrategy)
            {
                m_storage.resize(storageSize);
                ptr = &m_storage[0];
            }
 
            for (unsigned int n = 0; n < nBlk; ++n)
            {
                loc_mat = m_schurCompl->GetBlock(n, n);
 
                if (MultiRegions::eContiguous == storageStrategy)
                {
                    int loc_lda      = loc_mat->GetRows();
                    int blockSize    = loc_lda * loc_lda;
                    m_denseBlocks[n] = ptr;
                    for (int i = 0; i < loc_lda; ++i)
                    {
                        for (int j = 0; j < loc_lda; ++j)
                        {
                            ptr[j * loc_lda + i] = (*loc_mat)(i, j);
                        }
                    }
                    ptr += blockSize;
                    GlobalLinSys::v_DropStaticCondBlock(n);
                }
                else
                {
                    m_denseBlocks[n] = loc_mat->GetRawPtr();
                }
            }
            break;
        }
        case MultiRegions::eSparse:
        {
            DNekScalMatSharedPtr loc_mat;
            int loc_lda;
            int blockSize = 0;
 
            // First run through to split the set of local matrices into
            // partitions of fixed block size, and count number of local
            // matrices that belong to each partition.
            std::vector<std::pair<int, int>> partitions;
            for (int n = 0; n < m_schurCompl->GetNumberOfBlockRows(); ++n)
            {
                loc_mat = m_schurCompl->GetBlock(n, n);
                loc_lda = loc_mat->GetRows();
 
                ASSERTL1(loc_lda >= 0,
                         std::to_string(n) +
                             "-th "
                             "matrix block in Schur complement has "
                             "rank 0!");
 
                if (blockSize == loc_lda)
                {
                    partitions[partitions.size() - 1].first++;
                }
                else
                {
                    blockSize = loc_lda;
                    partitions.push_back(make_pair(1, loc_lda));
                }
            }
 
            MatrixStorage matStorage = eFULL;
 
            // Create a vector of sparse storage holders
            DNekSmvBsrDiagBlkMat::SparseStorageSharedPtrVector sparseStorage(
                partitions.size());
 
            for (int part = 0, n = 0; part < partitions.size(); ++part)
            {
                BCOMatType partMat;
 
                for (int k = 0; k < partitions[part].first; ++k, ++n)
                {
                    loc_mat = m_schurCompl->GetBlock(n, n);
                    loc_lda = loc_mat->GetRows();
 
                    ASSERTL1(loc_lda == partitions[part].second,
                             std::to_string(n) +
                                 "-th"
                                 " matrix block in Schur complement has "
                                 "unexpected rank");
 
                    NekDouble scale = loc_mat->Scale();
                    if (fabs(scale - 1.0) > NekConstants::kNekZeroTol)
                    {
                        Array<OneD, NekDouble> matarray(loc_lda * loc_lda);
                        Vmath::Smul(loc_lda * loc_lda, scale,
                                    loc_mat->GetRawPtr(), 1, &matarray[0], 1);
                        partMat[make_pair(k, k)] = BCOEntryType(matarray);
                    }
                    else // scale factor is 1.0
                    {
                        partMat[make_pair(k, k)] = BCOEntryType(
                            loc_lda * loc_lda, loc_mat->GetRawPtr());
                    }
 
                    GlobalLinSys::v_DropStaticCondBlock(n);
                }
 
                sparseStorage[part] =
                    MemoryManager<DNekSmvBsrDiagBlkMat::StorageType>::
                        AllocateSharedPtr(
                            partitions[part].first, partitions[part].first,
                            partitions[part].second, partMat, matStorage);
            }
 
            // Create block diagonal matrix
            m_sparseSchurCompl =
                MemoryManager<DNekSmvBsrDiagBlkMat>::AllocateSharedPtr(
                    sparseStorage);
 
            break;
        }
        default:
            ErrorUtil::NekError("Solver info property \
                        LocalMatrixStorageStrategy takes values \
                        Contiguous, Non-contiguous and Sparse");
    }
}

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL1, Nektar::MultiRegions::eContiguous, Nektar::eFULL, Nektar::MultiRegions::eNonContiguous, Nektar::MultiRegions::eSparse, Nektar::NekConstants::kNekZeroTol, m_denseBlocks, Nektar::MultiRegions::GlobalLinSys::m_expList, m_rows, m_scale, Nektar::MultiRegions::GlobalLinSysStaticCond::m_schurCompl, m_sparseSchurCompl, m_storage, Vmath::Smul(), and Nektar::MultiRegions::GlobalLinSys::v_DropStaticCondBlock().

Referenced by v_AssembleSchurComplement().

◆ v_AssembleSchurComplement()

void Nektar::MultiRegions::GlobalLinSysIterativeStaticCond::v_AssembleSchurComplement ( const std::shared_ptr< AssemblyMap > locToGloMap )

overrideprotectedvirtual

Assemble the Schur complement matrix.

Assemble the schur complement matrix from the block matrices stored in #m_blkMatrices and the given local to global mapping information.

Parameters

locToGloMap Local to global mapping information.

Reimplemented from Nektar::MultiRegions::GlobalLinSysStaticCond.

Definition at line 199 of file GlobalLinSysIterativeStaticCond.cpp.

{
    // Set up unique map
    v_UniqueMap();
 
    // Build precon again if we in multi-level static condensation (a
    // bit of a hack)
    if (m_linSysKey.GetGlobalSysSolnType() == eIterativeMultiLevelStaticCond)
    {
        m_precon = CreatePrecon(m_locToGloMap.lock());
        m_precon->BuildPreconditioner();
    }
 
    PrepareLocalSchurComplement();
}

References Nektar::MultiRegions::GlobalLinSys::CreatePrecon(), Nektar::MultiRegions::eIterativeMultiLevelStaticCond, Nektar::MultiRegions::GlobalLinSysKey::GetGlobalSysSolnType(), Nektar::MultiRegions::GlobalLinSys::m_linSysKey, Nektar::MultiRegions::GlobalLinSysStaticCond::m_locToGloMap, Nektar::MultiRegions::GlobalLinSysIterative::m_precon, PrepareLocalSchurComplement(), and v_UniqueMap().

◆ v_BasisFwdTransform()

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

overrideprotectedvirtual

Reimplemented from Nektar::MultiRegions::GlobalLinSysStaticCond.

Definition at line 493 of file GlobalLinSysIterativeStaticCond.cpp.

{
    m_precon->DoTransformBasisToLowEnergy(pInOut);
}

References Nektar::MultiRegions::GlobalLinSysIterative::m_precon.

◆ v_CoeffsBwdTransform()

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

overrideprotectedvirtual

Reimplemented from Nektar::MultiRegions::GlobalLinSysStaticCond.

Definition at line 499 of file GlobalLinSysIterativeStaticCond.cpp.

{
    m_precon->DoTransformCoeffsFromLowEnergy(pInOut);
}

References Nektar::MultiRegions::GlobalLinSysIterative::m_precon.

◆ v_CoeffsFwdTransform()

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

overrideprotectedvirtual

Reimplemented from Nektar::MultiRegions::GlobalLinSysStaticCond.

Definition at line 505 of file GlobalLinSysIterativeStaticCond.cpp.

{
    m_precon->DoTransformCoeffsToLowEnergy(pInput, pOutput);
}

References Nektar::MultiRegions::GlobalLinSysIterative::m_precon.

◆ v_DoMatrixMultiply()

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

overrideprotectedvirtual

Perform a Shur-complement matrix multiply operation.

Implements Nektar::MultiRegions::GlobalLinSysIterative.

Definition at line 378 of file GlobalLinSysIterativeStaticCond.cpp.

{
    if (m_linsol->IsLocal())
    {
        if (m_sparseSchurCompl)
        {
            m_sparseSchurCompl->Multiply(pInput, pOutput);
        }
        else
        {
            // Do matrix multiply locally, using direct BLAS calls
            int i, cnt;
            for (i = cnt = 0; i < m_denseBlocks.size(); cnt += m_rows[i], ++i)
            {
                const int rows = m_rows[i];
                Blas::Dgemv('N', rows, rows, m_scale[i], m_denseBlocks[i], rows,
                            pInput.get() + cnt, 1, 0.0, pOutput.get() + cnt, 1);
            }
        }
    }
    else
    {
        int nLocal = m_locToGloMap.lock()->GetNumLocalBndCoeffs();
        AssemblyMapSharedPtr asmMap = m_locToGloMap.lock();
 
        if (m_sparseSchurCompl)
        {
            // Do matrix multiply locally using block-diagonal sparse matrix
            Array<OneD, NekDouble> tmp = m_wsp + nLocal;
 
            asmMap->GlobalToLocalBnd(pInput, m_wsp);
            m_sparseSchurCompl->Multiply(m_wsp, tmp);
            asmMap->AssembleBnd(tmp, pOutput);
        }
        else
        {
            // Do matrix multiply locally, using direct BLAS calls
            asmMap->GlobalToLocalBnd(pInput, m_wsp);
            int i, cnt;
            Array<OneD, NekDouble> tmpout = m_wsp + nLocal;
            for (i = cnt = 0; i < m_denseBlocks.size(); cnt += m_rows[i], ++i)
            {
                const int rows = m_rows[i];
                Blas::Dgemv('N', rows, rows, m_scale[i], m_denseBlocks[i], rows,
                            m_wsp.get() + cnt, 1, 0.0, tmpout.get() + cnt, 1);
            }
            asmMap->AssembleBnd(tmpout, pOutput);
        }
    }
}

References Blas::Dgemv(), m_denseBlocks, Nektar::MultiRegions::GlobalLinSysIterative::m_linsol, Nektar::MultiRegions::GlobalLinSysStaticCond::m_locToGloMap, m_rows, m_scale, m_sparseSchurCompl, and Nektar::MultiRegions::GlobalLinSysStaticCond::m_wsp.

◆ v_GetStaticCondBlock()

DNekScalBlkMatSharedPtr Nektar::MultiRegions::GlobalLinSysIterativeStaticCond::v_GetStaticCondBlock ( unsigned int n )

overrideprotectedvirtual

Retrieves a the static condensation block matrices from n-th expansion using the matrix key provided by the m_linSysKey.

Parameters

n	Number of the expansion

Returns: 2x2 Block matrix holding the static condensation matrices for the n-th expansion.

Reimplemented from Nektar::MultiRegions::GlobalLinSys.

Definition at line 178 of file GlobalLinSysIterativeStaticCond.cpp.

{
    DNekScalBlkMatSharedPtr schurComplBlock;
    DNekScalMatSharedPtr localMat = m_schurCompl->GetBlock(n, n);
    unsigned int nbdry            = localMat->GetRows();
    unsigned int nblks            = 1;
    unsigned int esize[1]         = {nbdry};
 
    schurComplBlock = MemoryManager<DNekScalBlkMat>::AllocateSharedPtr(
        nblks, nblks, esize, esize);
    schurComplBlock->SetBlock(0, 0, localMat);
 
    return schurComplBlock;
}

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), and Nektar::MultiRegions::GlobalLinSysStaticCond::m_schurCompl.

◆ v_InitObject()

void Nektar::MultiRegions::GlobalLinSysIterativeStaticCond::v_InitObject ( )

overrideprotectedvirtual

Reimplemented from Nektar::MultiRegions::GlobalLinSys.

Definition at line 138 of file GlobalLinSysIterativeStaticCond.cpp.

{
    auto asmMap = m_locToGloMap.lock();
 
    m_precon = CreatePrecon(asmMap);
 
    // Allocate memory for top-level structure
    SetupTopLevel(asmMap);
 
    // Setup Block Matrix systems
    int n, n_exp = m_expList.lock()->GetNumElmts();
 
    // Build preconditioner
    m_precon->BuildPreconditioner();
 
    // Do transform of Schur complement matrix
    int cnt = 0;
    for (n = 0; n < n_exp; ++n)
    {
        if (m_linSysKey.GetMatrixType() != StdRegions::eHybridDGHelmBndLam)
        {
            DNekScalMatSharedPtr mat = m_schurCompl->GetBlock(n, n);
            DNekScalMatSharedPtr t =
                m_precon->TransformedSchurCompl(n, cnt, mat);
            m_schurCompl->SetBlock(n, n, t);
            cnt += mat->GetRows();
        }
    }
 
    // Construct this level
    Initialise(asmMap);
}

◆ v_PreSolve()

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

overrideprotectedvirtual

Reimplemented from Nektar::MultiRegions::GlobalLinSysStaticCond.

Definition at line 435 of file GlobalLinSysIterativeStaticCond.cpp.

{
    if (m_isAbsoluteTolerance)
    {
        m_rhs_magnitude = 1.0;
        return;
    }
 
    if (scLevel == 0)
    {
        // When matrices are supplied to the constructor at the top
        // level, the preconditioner is never set up.
        if (!m_precon)
        {
            m_precon = CreatePrecon(m_locToGloMap.lock());
            m_precon->BuildPreconditioner();
        }
 
        int nGloBndDofs = m_locToGloMap.lock()->GetNumGlobalBndCoeffs();
        Array<OneD, NekDouble> F_gloBnd(nGloBndDofs);
        NekVector<NekDouble> F_GloBnd(nGloBndDofs, F_gloBnd, eWrapper);
        m_locToGloMap.lock()->AssembleBnd(F_bnd, F_gloBnd);
 
        NekDouble vExchange(0.0);
        if (m_map.size() > 0)
        {
            vExchange = Vmath::Dot2(F_GloBnd.GetDimension(), &F_GloBnd[0],
                                    &F_GloBnd[0], &m_map[0]);
        }
 
        m_expList.lock()->GetComm()->GetRowComm()->AllReduce(
            vExchange, Nektar::LibUtilities::ReduceSum);
 
        // To ensure that very different rhs values are not being
        // used in subsequent solvers such as the velocity solve in
        // INC NS. If this works we then need to work out a better
        // way to control this.
        NekDouble new_rhs_mag = (vExchange > 1e-6) ? vExchange : 1.0;
 
        if (m_rhs_magnitude == NekConstants::kNekUnsetDouble)
        {
            m_rhs_magnitude = new_rhs_mag;
        }
        else
        {
            m_rhs_magnitude = (m_rhs_mag_sm * (m_rhs_magnitude) +
                               (1.0 - m_rhs_mag_sm) * new_rhs_mag);
        }
    }
    else
    {
        // for multilevel iterative solver always use rhs
        // vector value with no weighting
        m_rhs_magnitude = NekConstants::kNekUnsetDouble;
    }
}

References Nektar::MultiRegions::GlobalLinSys::CreatePrecon(), Vmath::Dot2(), Nektar::eWrapper, Nektar::NekVector< DataType >::GetDimension(), Nektar::NekConstants::kNekUnsetDouble, Nektar::MultiRegions::GlobalLinSys::m_expList, Nektar::MultiRegions::GlobalLinSysIterative::m_isAbsoluteTolerance, Nektar::MultiRegions::GlobalLinSysStaticCond::m_locToGloMap, Nektar::MultiRegions::GlobalLinSysIterative::m_map, Nektar::MultiRegions::GlobalLinSysIterative::m_precon, Nektar::MultiRegions::GlobalLinSysIterative::m_rhs_mag_sm, Nektar::MultiRegions::GlobalLinSysIterative::m_rhs_magnitude, and Nektar::LibUtilities::ReduceSum.

◆ v_Recurse()

GlobalLinSysStaticCondSharedPtr Nektar::MultiRegions::GlobalLinSysIterativeStaticCond::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
	)

overrideprotectedvirtual

Implements Nektar::MultiRegions::GlobalLinSysStaticCond.

Definition at line 511 of file GlobalLinSysIterativeStaticCond.cpp.

{
    GlobalLinSysIterativeStaticCondSharedPtr sys =
        MemoryManager<GlobalLinSysIterativeStaticCond>::AllocateSharedPtr(
            mkey, pExpList, pSchurCompl, pBinvD, pC, pInvD, l2gMap, m_precon);
    sys->Initialise(l2gMap);
    return sys;
}

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), and Nektar::MultiRegions::GlobalLinSysIterative::m_precon.

◆ v_SolveLinearSystem()

void Nektar::MultiRegions::GlobalLinSysIterativeStaticCond::v_SolveLinearSystem	(	const int	pNumRows,
		const Array< OneD, const NekDouble > &	pInput,
		Array< OneD, NekDouble > &	pOutput,
		const AssemblyMapSharedPtr &	locToGloMap,
		const int	pNumDir
	)

overrideprivatevirtual

Solve the matrix system.

Implements Nektar::MultiRegions::GlobalLinSys.

Definition at line 528 of file GlobalLinSysIterativeStaticCond.cpp.

{
 
    if (!m_linsol)
    {
        LibUtilities::CommSharedPtr vRowComm =
            m_expList.lock()->GetComm()->GetRowComm();
        LibUtilities::SessionReaderSharedPtr pSession =
            m_expList.lock()->GetSession();
 
        // Check such a module exists for this equation.
        ASSERTL0(LibUtilities::GetNekLinSysIterFactory().ModuleExists(
                     m_linSysIterSolver),
                 "NekLinSysIter '" + m_linSysIterSolver +
                     "' is not defined.\n");
 
        // Create the key to hold solver settings
        auto sysKey     = LibUtilities::NekSysKey();
        string variable = plocToGloMap->GetVariable();
 
        // Either get the solnInfo from <GlobalSysSolInfo> or from
        // <Parameters>
        if (pSession->DefinesGlobalSysSolnInfo(variable,
                                               "NekLinSysMaxIterations"))
        {
            sysKey.m_NekLinSysMaxIterations = boost::lexical_cast<int>(
                pSession
                    ->GetGlobalSysSolnInfo(variable, "NekLinSysMaxIterations")
                    .c_str());
        }
        else
        {
            pSession->LoadParameter("NekLinSysMaxIterations",
                                    sysKey.m_NekLinSysMaxIterations, 5000);
        }
 
        if (pSession->DefinesGlobalSysSolnInfo(variable, "LinSysMaxStorage"))
        {
            sysKey.m_LinSysMaxStorage = boost::lexical_cast<int>(
                pSession->GetGlobalSysSolnInfo(variable, "LinSysMaxStorage")
                    .c_str());
        }
        else
        {
            pSession->LoadParameter("LinSysMaxStorage",
                                    sysKey.m_LinSysMaxStorage, 100);
        }
 
        if (pSession->DefinesGlobalSysSolnInfo(variable, "GMRESMaxHessMatBand"))
        {
            sysKey.m_KrylovMaxHessMatBand = boost::lexical_cast<int>(
                pSession->GetGlobalSysSolnInfo(variable, "GMRESMaxHessMatBand")
                    .c_str());
        }
        else
        {
            pSession->LoadParameter("GMRESMaxHessMatBand",
                                    sysKey.m_KrylovMaxHessMatBand,
                                    sysKey.m_LinSysMaxStorage + 1);
        }
 
        // The following settings have no correponding tests and are rarely
        // used.
        pSession->MatchSolverInfo("GMRESLeftPrecon", "True",
                                  sysKey.m_NekLinSysLeftPrecon, false);
        pSession->MatchSolverInfo("GMRESRightPrecon", "True",
                                  sysKey.m_NekLinSysRightPrecon, true);
 
        m_linsol = LibUtilities::GetNekLinSysIterFactory().CreateInstance(
            m_linSysIterSolver, pSession, vRowComm, nGlobal - nDir, sysKey);
 
        m_linsol->SetSysOperators(m_NekSysOp);
        v_UniqueMap();
        m_linsol->SetUniversalUniqueMap(m_map);
    }
 
    if (!m_precon)
    {
        m_precon = CreatePrecon(plocToGloMap);
        m_precon->BuildPreconditioner();
    }
 
    m_linsol->setRhsMagnitude(m_rhs_magnitude);
 
    if (m_useProjection)
    {
        Array<OneD, NekDouble> gloIn(nGlobal);
        Array<OneD, NekDouble> gloOut(nGlobal, 0.0);
        plocToGloMap->AssembleBnd(pInput, gloIn);
        DoProjection(nGlobal, gloIn, gloOut, nDir, m_tolerance, m_isAconjugate);
        plocToGloMap->GlobalToLocalBnd(gloOut, pOutput);
    }
    else
    {
        if (m_linsol->IsLocal())
        {
            int nLocDofs = plocToGloMap->GetNumLocalBndCoeffs();
            Vmath::Zero(nLocDofs, pOutput, 1);
            m_linsol->SolveSystem(nLocDofs, pInput, pOutput, nDir, m_tolerance);
        }
        else
        {
            Array<OneD, NekDouble> gloIn(nGlobal);
            Array<OneD, NekDouble> gloOut(nGlobal, 0.0);
            plocToGloMap->AssembleBnd(pInput, gloIn);
            m_linsol->SolveSystem(nGlobal, gloIn, gloOut, nDir, m_tolerance);
            plocToGloMap->GlobalToLocalBnd(gloOut, pOutput);
        }
    }
}

◆ v_UniqueMap()

void Nektar::MultiRegions::GlobalLinSysIterativeStaticCond::v_UniqueMap ( )

overrideprotectedvirtual

Implements Nektar::MultiRegions::GlobalLinSysIterative.

Definition at line 430 of file GlobalLinSysIterativeStaticCond.cpp.

{
    m_map = m_locToGloMap.lock()->GetGlobalToUniversalBndMapUnique();
}

References Nektar::MultiRegions::GlobalLinSysStaticCond::m_locToGloMap, and Nektar::MultiRegions::GlobalLinSysIterative::m_map.

Referenced by v_AssembleSchurComplement(), and v_SolveLinearSystem().

Member Data Documentation

◆ className

string Nektar::MultiRegions::GlobalLinSysIterativeStaticCond::className

static

Initial value:

=
    GetGlobalLinSysFactory().RegisterCreatorFunction(
        "IterativeStaticCond", GlobalLinSysIterativeStaticCond::create,
        "Iterative static condensation.")

Name of class.

Registers the class with the Factory.

Definition at line 85 of file GlobalLinSysIterativeStaticCond.h.

◆ className2

string Nektar::MultiRegions::GlobalLinSysIterativeStaticCond::className2

static

Initial value:

=
    GetGlobalLinSysFactory().RegisterCreatorFunction(
        "IterativeMultiLevelStaticCond",
        GlobalLinSysIterativeStaticCond::create,
        "Iterative multi-level static condensation.")

Definition at line 86 of file GlobalLinSysIterativeStaticCond.h.

◆ m_denseBlocks

std::vector<const double *> Nektar::MultiRegions::GlobalLinSysIterativeStaticCond::m_denseBlocks

private

Vector of pointers to local matrix data.

Definition at line 135 of file GlobalLinSysIterativeStaticCond.h.

Referenced by PrepareLocalSchurComplement(), and v_DoMatrixMultiply().

◆ m_rows

Array<OneD, unsigned int> Nektar::MultiRegions::GlobalLinSysIterativeStaticCond::m_rows

private

Ranks of local matrices.

Definition at line 137 of file GlobalLinSysIterativeStaticCond.h.

Referenced by PrepareLocalSchurComplement(), and v_DoMatrixMultiply().

◆ m_scale

Array<OneD, NekDouble> Nektar::MultiRegions::GlobalLinSysIterativeStaticCond::m_scale

private

Scaling factors for local matrices.

Definition at line 139 of file GlobalLinSysIterativeStaticCond.h.

Referenced by PrepareLocalSchurComplement(), and v_DoMatrixMultiply().

◆ m_sparseSchurCompl

DNekSmvBsrDiagBlkMatSharedPtr Nektar::MultiRegions::GlobalLinSysIterativeStaticCond::m_sparseSchurCompl

private

Sparse representation of Schur complement matrix at this level.

Definition at line 141 of file GlobalLinSysIterativeStaticCond.h.

Referenced by PrepareLocalSchurComplement(), and v_DoMatrixMultiply().

◆ m_storage

std::vector<double> Nektar::MultiRegions::GlobalLinSysIterativeStaticCond::m_storage

private

Dense storage for block Schur complement matrix.

Definition at line 133 of file GlobalLinSysIterativeStaticCond.h.

Referenced by PrepareLocalSchurComplement().

◆ storagedef

std::string Nektar::MultiRegions::GlobalLinSysIterativeStaticCond::storagedef

staticprivate

Initial value:

=
    LibUtilities::SessionReader::RegisterDefaultSolverInfo(
        "LocalMatrixStorageStrategy", "Sparse")

Utility strings.

Definition at line 143 of file GlobalLinSysIterativeStaticCond.h.

◆ storagelookupIds

std::string Nektar::MultiRegions::GlobalLinSysIterativeStaticCond::storagelookupIds

staticprivate

Initial value:

= {
    LibUtilities::SessionReader::RegisterEnumValue(
        "LocalMatrixStorageStrategy", "Contiguous", MultiRegions::eContiguous),
    LibUtilities::SessionReader::RegisterEnumValue(
        "LocalMatrixStorageStrategy", "Non-contiguous",
        MultiRegions::eNonContiguous),
    LibUtilities::SessionReader::RegisterEnumValue(
        "LocalMatrixStorageStrategy", "Sparse", MultiRegions::eSparse),
}

Definition at line 144 of file GlobalLinSysIterativeStaticCond.h.

Public Types

Public Member Functions

Static Public Member Functions

Static Public Attributes

Protected Member Functions

Private Member Functions

Private Attributes

Static Private Attributes

Additional Inherited Members

Detailed Description

Member Typedef Documentation

◆ DNekSmvBsrDiagBlkMat

◆ DNekSmvBsrDiagBlkMatSharedPtr

Constructor & Destructor Documentation

◆ GlobalLinSysIterativeStaticCond() [1/2]

◆ GlobalLinSysIterativeStaticCond() [2/2]

◆ ~GlobalLinSysIterativeStaticCond()

Member Function Documentation

◆ create()

◆ PrepareLocalSchurComplement()

◆ v_AssembleSchurComplement()

◆ v_BasisFwdTransform()

◆ v_CoeffsBwdTransform()

◆ v_CoeffsFwdTransform()

◆ v_DoMatrixMultiply()

◆ v_GetStaticCondBlock()

◆ v_InitObject()

◆ v_PreSolve()

◆ v_Recurse()

◆ v_SolveLinearSystem()

◆ v_UniqueMap()

Member Data Documentation

◆ className

◆ className2

◆ m_denseBlocks

◆ m_rows

◆ m_scale

◆ m_sparseSchurCompl

◆ m_storage

◆ storagedef

◆ storagelookupIds