A global linear system. More...

#include <GlobalLinSysPETSc.h>

Inheritance diagram for Nektar::MultiRegions::GlobalLinSysPETSc:

Collaboration diagram for Nektar::MultiRegions::GlobalLinSysPETSc:

Public Member Functions
	GlobalLinSysPETSc (const GlobalLinSysKey &pKey, const boost::weak_ptr< ExpList > &pExp, const boost::shared_ptr< AssemblyMap > &pLocToGloMap)
	Constructor for full direct matrix solve.
virtual	~GlobalLinSysPETSc ()
virtual void	v_SolveLinearSystem (const int pNumRows, const Array< OneD, const NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const AssemblyMapSharedPtr &locToGloMap, const int pNumDir)
	Solve a basic matrix system.
Public Member Functions inherited from Nektar::MultiRegions::GlobalLinSys
	GlobalLinSys (const GlobalLinSysKey &pKey, const boost::weak_ptr< ExpList > &pExpList, const boost::shared_ptr< AssemblyMap > &pLocToGloMap)
	Constructor for full direct matrix solve.
virtual	~GlobalLinSys ()
const GlobalLinSysKey &	GetKey (void) const
	Returns the key associated with the system.
const boost::weak_ptr< ExpList > &	GetLocMat (void) const
void	InitObject ()
void	Initialise (const boost::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.
boost::shared_ptr< GlobalLinSys >	GetSharedThisPtr ()
	Returns a shared pointer to the current object.
int	GetNumBlocks ()
DNekScalMatSharedPtr	GetBlock (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.

Protected Member Functions
void	SetUpScatter ()
void	SetUpMatVec ()
void	SetUpSolver (NekDouble tolerance)
void	CalculateReordering (const Array< OneD, const int > &glo2uniMap, const Array< OneD, const int > &glo2unique, const AssemblyMapSharedPtr &pLocToGloMap)
Protected Member Functions inherited from Nektar::MultiRegions::GlobalLinSys
virtual int	v_GetNumBlocks ()
	Get the number of blocks in this system.
virtual DNekScalMatSharedPtr	v_GetBlock (unsigned int n)
	Retrieves the block matrix from n-th expansion using the matrix key provided by the m_linSysKey.
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.
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.

Protected Attributes
Mat	m_matrix
Vec	m_x
Vec	m_b
Vec	m_locVec
KSP	m_ksp
vector< int >	m_reorderedMap
VecScatter	m_ctx
int	m_nLocal
Protected Attributes inherited from Nektar::MultiRegions::GlobalLinSys
const GlobalLinSysKey	m_linSysKey
	Key associated with this linear system.
const boost::weak_ptr< ExpList >	m_expList
	Local Matrix System.
const map< int, RobinBCInfoSharedPtr >	m_robinBCInfo
	Robin boundary info.

Detailed Description

A global linear system.

Solves a linear system using direct methods.

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

Definition at line 52 of file GlobalLinSysPETSc.h.

Constructor & Destructor Documentation

Nektar::MultiRegions::GlobalLinSysPETSc::GlobalLinSysPETSc	(	const GlobalLinSysKey &	pKey,
		const boost::weak_ptr< ExpList > &	pExp,
		const boost::shared_ptr< AssemblyMap > &	pLocToGloMap
	)

Constructor for full direct matrix solve.

Definition at line 53 of file GlobalLinSysPETSc.cpp.

References m_matrix.

            : GlobalLinSys(pKey, pExp, pLocToGloMap)
        {
            // Initialise PETSc
            PetscInitialize(0, NULL, NULL, NULL);
            // Create matrix
            MatCreate(PETSC_COMM_WORLD, &m_matrix);
        }

Nektar::MultiRegions::GlobalLinSysPETSc::~GlobalLinSysPETSc ( )

virtual

Definition at line 66 of file GlobalLinSysPETSc.cpp.

{

}

Member Function Documentation

void Nektar::MultiRegions::GlobalLinSysPETSc::CalculateReordering	(	const Array< OneD, const int > &	glo2uniMap,
		const Array< OneD, const int > &	glo2unique,
		const AssemblyMapSharedPtr &	pLocToGloMap
	)

protected

Definition at line 128 of file GlobalLinSysPETSc.cpp.

References ASSERTL0, Nektar::MultiRegions::GlobalLinSys::m_expList, m_nLocal, m_reorderedMap, Nektar::LibUtilities::ReduceSum, and Vmath::Vsum().

Referenced by Nektar::MultiRegions::GlobalLinSysPETScFull::GlobalLinSysPETScFull(), and Nektar::MultiRegions::GlobalLinSysPETScStaticCond::v_AssembleSchurComplement().

        {
            LibUtilities::CommSharedPtr vComm
                = m_expList.lock()->GetSession()->GetComm();
            const int nDirDofs = pLocToGloMap->GetNumGlobalDirBndCoeffs();
            const int nHomDofs = glo2uniMap.num_elements() - nDirDofs;
            const int nProc    = vComm->GetSize();
            const int rank     = vComm->GetRank();
            int n, cnt;
            m_nLocal = Vmath::Vsum(nHomDofs, glo2unique + nDirDofs, 1);
            m_reorderedMap.resize(nHomDofs);
            Array<OneD, int> localCounts(nProc, 0), localOffset(nProc, 0);
            localCounts[rank] = nHomDofs;
            vComm->AllReduce(localCounts, LibUtilities::ReduceSum);
            for (n = 1; n < nProc; ++n)
            {
                localOffset[n] = localOffset[n-1] + localCounts[n-1];
            }
            int totHomDofs = Vmath::Vsum(nProc, localCounts, 1);
            vector<unsigned int> allUniIds(totHomDofs, 0);
            for (n = 0; n < nHomDofs; ++n)
            {
                int gid = n + nDirDofs;
                allUniIds[n + localOffset[rank]] = glo2uniMap[gid];
            }
            vComm->AllReduce(allUniIds, LibUtilities::ReduceSum);
            std::sort(allUniIds.begin(), allUniIds.end());
            map<int,int> uniIdReorder;
            for (cnt = n = 0; n < allUniIds.size(); ++n)
            {
                if (uniIdReorder.count(allUniIds[n]) > 0)
                {
                    continue;
                }
                uniIdReorder[allUniIds[n]] = cnt++;
            }
            for (n = 0; n < nHomDofs; ++n)
            {
                int gid = n + nDirDofs;
                int uniId = glo2uniMap[gid];
                ASSERTL0(uniIdReorder.count(uniId) > 0, "Error in ordering");
                m_reorderedMap[n] = uniIdReorder[uniId];
            }
        }

void Nektar::MultiRegions::GlobalLinSysPETSc::SetUpMatVec ( )

protected

Definition at line 188 of file GlobalLinSysPETSc.cpp.

References m_b, m_matrix, m_nLocal, and m_x.

Referenced by Nektar::MultiRegions::GlobalLinSysPETScFull::GlobalLinSysPETScFull(), and Nektar::MultiRegions::GlobalLinSysPETScStaticCond::v_AssembleSchurComplement().

        {
            // CREATE VECTORS
            VecCreate        (PETSC_COMM_WORLD, &m_x);
            VecSetSizes      (m_x, m_nLocal, PETSC_DECIDE);
            VecSetFromOptions(m_x);
            VecDuplicate     (m_x, &m_b);
            // CREATE MATRICES
            MatCreate        (PETSC_COMM_WORLD, &m_matrix);
            MatSetType       (m_matrix, MATMPIAIJ);
            MatSetSizes      (m_matrix, m_nLocal, m_nLocal,
                              PETSC_DETERMINE, PETSC_DETERMINE);
            MatSetFromOptions(m_matrix);
            MatSetUp         (m_matrix);
        }

void Nektar::MultiRegions::GlobalLinSysPETSc::SetUpScatter ( )

protected

Definition at line 105 of file GlobalLinSysPETSc.cpp.

References m_ctx, m_locVec, m_reorderedMap, and m_x.

Referenced by Nektar::MultiRegions::GlobalLinSysPETScFull::GlobalLinSysPETScFull(), and Nektar::MultiRegions::GlobalLinSysPETScStaticCond::v_AssembleSchurComplement().

        {
            const int nHomDofs = m_reorderedMap.size();
            // Create local and global numbering systems for vector
            IS isGlobal, isLocal;
            ISCreateGeneral(PETSC_COMM_SELF, nHomDofs, &m_reorderedMap[0],
                            PETSC_COPY_VALUES, &isGlobal);
            ISCreateStride (PETSC_COMM_SELF, nHomDofs, 0, 1, &isLocal);
            // Create local vector for output
            VecCreate        (PETSC_COMM_SELF, &m_locVec);
            VecSetSizes      (m_locVec, m_reorderedMap.size(), PETSC_DECIDE);
            VecSetFromOptions(m_locVec);
            // Create scatter context
            VecScatterCreate (m_x, isGlobal, m_locVec, isLocal, &m_ctx);
            // Clean up
            ISDestroy(&isGlobal);
            ISDestroy(&isLocal);
        }

void Nektar::MultiRegions::GlobalLinSysPETSc::SetUpSolver ( NekDouble tolerance )

protected

Definition at line 205 of file GlobalLinSysPETSc.cpp.

References m_ksp, and m_matrix.

Referenced by Nektar::MultiRegions::GlobalLinSysPETScFull::GlobalLinSysPETScFull(), and Nektar::MultiRegions::GlobalLinSysPETScStaticCond::v_AssembleSchurComplement().

        {
            KSPCreate(PETSC_COMM_WORLD, &m_ksp);
            KSPSetTolerances(
                m_ksp, tolerance, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT);
            KSPSetFromOptions(m_ksp);
#if PETSC_VERSION_GE(3,5,0)
            KSPSetOperators(m_ksp, m_matrix, m_matrix);
#else
            KSPSetOperators(m_ksp, m_matrix, m_matrix, SAME_NONZERO_PATTERN);
#endif
        }

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

virtual

Solve a basic matrix system.

Implements Nektar::MultiRegions::GlobalLinSys.

Definition at line 70 of file GlobalLinSysPETSc.cpp.

References m_b, m_ctx, m_ksp, m_locVec, m_reorderedMap, m_x, and Vmath::Vcopy().

        {
            const int nHomDofs = pNumRows - pNumDir;
            // Populate RHS vector from input
            VecSetValues(m_b, nHomDofs, &m_reorderedMap[0], &pInput[pNumDir], INSERT_VALUES);
            // Assemble RHS vector
            VecAssemblyBegin(m_b);
            VecAssemblyEnd  (m_b);
            // Do system solve
            PetscErrorCode ierr = KSPSolve(m_ksp, m_b, m_x);
            // Grab number of iterations taken
            PetscInt its;
            KSPGetIterationNumber(m_ksp, &its);
            cout << "iterations = " << its << endl;
            // Scatter results to local vector
            VecScatterBegin(m_ctx, m_x, m_locVec, INSERT_VALUES, SCATTER_FORWARD);
            VecScatterEnd  (m_ctx, m_x, m_locVec, INSERT_VALUES, SCATTER_FORWARD);
            // Copy results into output vector
            PetscScalar *avec;
            VecGetArray    (m_locVec, &avec);
            Vmath::Vcopy   (nHomDofs, avec, 1, &pOutput[pNumDir], 1);
            VecRestoreArray(m_locVec, &avec);
        }