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

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 linear system using PETSc. More...

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

virtual	~GlobalLinSys ()

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

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

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

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

Protected Member Functions
void	SetUpScatter ()
	Set up PETSc local (equivalent to Nektar++ global) and global (equivalent to universal) scatter maps. More...

void	SetUpMatVec ()
	Construct PETSc matrix and vector handles. More...

void	SetUpSolver (NekDouble tolerance)
	Set up KSP solver object. More...

void	CalculateReordering (const Array< OneD, const int > &glo2uniMap, const Array< OneD, const int > &glo2unique, const AssemblyMapSharedPtr &pLocToGloMap)
	Calculate a reordering of universal IDs for PETSc. More...

Protected Member Functions inherited from Nektar::MultiRegions::GlobalLinSys
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 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...

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

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

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

Detailed Description

A global linear system.

Solves a linear system using PETSc.

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 52 of file GlobalLinSysPETSc.cpp.

References m_matrix.

             : GlobalLinSys(pKey, pExp, pLocToGloMap)
         {
             // Check PETSc is initialized
             // For some reason, this is needed on OS X as logging is not
             // initialized properly in the call within CommMpi.
             PetscInitializeNoArguments();
 
             // Create matrix
             MatCreate(PETSC_COMM_WORLD, &m_matrix);
         }

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

virtual

Definition at line 67 of file GlobalLinSysPETSc.cpp.

68 {

69 }

Member Function Documentation

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

protected

Calculate a reordering of universal IDs for PETSc.

PETSc requires a unique, contiguous index of all global and universal degrees of freedom which represents its position inside the matrix. Presently Gs does not guarantee this, so this routine constructs a new universal mapping.

Definition at line 152 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;
 
             // Count number of unique degrees of freedom on each process.
             m_nLocal = Vmath::Vsum(nHomDofs, glo2unique + nDirDofs, 1);
             m_reorderedMap.resize(nHomDofs);
 
             // Reduce coefficient counts across all processors.
             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);
 
             // Assemble list of universal IDs
             for (n = 0; n < nHomDofs; ++n)
             {
                 int gid = n + nDirDofs;
                 allUniIds[n + localOffset[rank]] = glo2uniMap[gid];
             }
 
             // Reduce this across processors so that each process has a list of
             // all universal IDs.
             vComm->AllReduce(allUniIds, LibUtilities::ReduceSum);
             std::sort(allUniIds.begin(), allUniIds.end());
             map<int,int> uniIdReorder;
 
             // Renumber starting from 0.
             for (cnt = n = 0; n < allUniIds.size(); ++n)
             {
                 if (uniIdReorder.count(allUniIds[n]) > 0)
                 {
                     continue;
                 }
 
                 uniIdReorder[allUniIds[n]] = cnt++;
             }
 
             // Populate reordering map.
             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

Construct PETSc matrix and vector handles.

Todo:: Preallocation should be done at this point, since presently matrix allocation takes a significant amount of time.

Definition at line 224 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, MATAIJ);
             MatSetSizes      (m_matrix, m_nLocal, m_nLocal,
                               PETSC_DETERMINE, PETSC_DETERMINE);
             MatSetFromOptions(m_matrix);
             MatSetUp         (m_matrix);
         }

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

protected

Set up PETSc local (equivalent to Nektar++ global) and global (equivalent to universal) scatter maps.

These maps are used in GlobalLinSysPETSc::v_SolveLinearSystem to scatter the solution vector back to each process.

Definition at line 121 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, nHomDofs, 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

Set up KSP solver object.

This is reasonably generic setup – most solver types can be changed using the .petscrc file.

Definition at line 247 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 linear system using PETSc.

The general strategy being a PETSc solve is to:

Copy values into the PETSc vector m_b
Solve the system m_ksp and place result into m_x.
Scatter results back into m_locVec using m_ctx scatter object.
Copy from m_locVec to output array #pOutput.

Implements Nektar::MultiRegions::GlobalLinSys.

Definition at line 81 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
             KSPSolve(m_ksp, m_b, m_x);
 
             // 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 *tmp;
             VecGetArray    (m_locVec, &tmp);
             Vmath::Vcopy   (nHomDofs, tmp, 1, &pOutput[pNumDir], 1);
             VecRestoreArray(m_locVec, &tmp);
         }