A PETSc global linear system. More...

#include <GlobalLinSysPETSc.h>

Inheritance diagram for Nektar::MultiRegions::GlobalLinSysPETSc:

Collaboration diagram for Nektar::MultiRegions::GlobalLinSysPETSc:

Classes
struct	ShellCtx
	Internal struct for MatShell and PCShell calls to store current context for callback. More...

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 ()
	Clean up PETSc objects. More...

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 (int nGlobal, int nDir)
	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...

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

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

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

Protected Attributes
Mat	m_matrix
	PETSc matrix object. More...

Vec	m_x
	PETSc vector objects used for local storage. More...

Vec	m_b

Vec	m_locVec

KSP	m_ksp
	KSP object that represents solver system. More...

PC	m_pc
	PCShell for preconditioner. More...

PETScMatMult	m_matMult
	Enumerator to select matrix multiplication type. More...

std::vector< int >	m_reorderedMap
	Reordering that takes universal IDs to a unique row in the PETSc matrix. More...

VecScatter	m_ctx
	PETSc scatter context that takes us between Nektar++ global ordering and PETSc vector ordering. More...

int	m_nLocal
	Number of unique degrees of freedom on this process. More...

PreconditionerSharedPtr	m_precon

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 std::map< int, RobinBCInfoSharedPtr >	m_robinBCInfo
	Robin boundary info. More...

bool	m_verbose

Static Private Member Functions
static PetscErrorCode	DoMatrixMultiply (Mat M, Vec in, Vec out)
	Perform matrix multiplication using Nektar++ routines. More...

static PetscErrorCode	DoPreconditioner (PC pc, Vec in, Vec out)
	Apply preconditioning using Nektar++ routines. More...

static void	DoNekppOperation (Vec &in, Vec &out, ShellCtx *ctx, bool precon)
	Perform either matrix multiplication or preconditioning using Nektar++ routines. More...

static PetscErrorCode	DoDestroyMatCtx (Mat M)
	Destroy matrix shell context object. More...

static PetscErrorCode	DoDestroyPCCtx (PC pc)
	Destroy preconditioner context object. More...

Static Private Attributes
static std::string	matMult

static std::string	matMultIds []

Detailed Description

A PETSc global linear system.

Solves a linear system using PETSc.

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

Definition at line 59 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 65 of file GlobalLinSysPETSc.cpp.

References Nektar::LibUtilities::CommMpi::GetComm(), Nektar::MultiRegions::GlobalLinSys::m_expList, m_matMult, and m_matrix.

             : GlobalLinSys(pKey, pExp, pLocToGloMap)
         {
             // Determine whether to use standard sparse matrix approach or
             // shell.
             m_matMult = pExp.lock()->GetSession()
                 ->GetSolverInfoAsEnum<PETScMatMult>(
                     "PETScMatMult");
 
             // Check PETSc is initialized. For some reason, this is needed on
             // OS X as logging is not initialized properly in the call within
             // CommMpi.
             PetscBool isInitialized;
             PetscInitialized(&isInitialized);
             if (!isInitialized)
             {
 #ifdef NEKTAR_USE_MPI
                 std::string commType =
                     m_expList.lock()->GetSession()->GetComm()->GetType();
                 if (commType.find("MPI") != std::string::npos)
                 {
                     LibUtilities::CommMpiSharedPtr comm =
                         boost::static_pointer_cast<LibUtilities::CommMpi>(
                             m_expList.lock()->GetSession()->GetComm());
                     PETSC_COMM_WORLD = comm->GetComm();
                 }
 #endif
                 PetscInitializeNoArguments();
             }
 
             // Create matrix
             MatCreate(PETSC_COMM_WORLD, &m_matrix);
         }

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

virtual

Clean up PETSc objects.

Note that if SessionReader::Finalize is called before the end of the program, PETSc may have been finalized already, at which point we cannot deallocate our objects. If that's the case we do nothing and let the kernel clear up after us.

Definition at line 110 of file GlobalLinSysPETSc.cpp.

References m_b, m_ksp, m_locVec, m_matrix, m_pc, and m_x.

         {
             PetscBool isFinalized;
             PetscFinalized(&isFinalized);
 
             // Sometimes, PetscFinalized returns false when (in fact) CommMpi's
             // Finalise routine has been called. We therefore also need to check
             // whether MPI has been finalised. This might arise from the
             // additional call to PetscInitializeNoArguments in the constructor
             // above.
 #ifdef NEKTAR_USE_MPI
             int mpiFinal = 0;
             MPI_Finalized(&mpiFinal);
             isFinalized = isFinalized || mpiFinal ? PETSC_TRUE : PETSC_FALSE;
 #endif
 
             if (!isFinalized)
             {
                 KSPDestroy(&m_ksp);
                 PCDestroy (&m_pc);
                 MatDestroy(&m_matrix);
                 VecDestroy(&m_x);
                 VecDestroy(&m_b);
                 VecDestroy(&m_locVec);
             }
         }

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.

Parameters

glo2uniMap	Global to universal map
glo2unique	Global to unique map
pLocToGloMap	Assembly map for this system

Definition at line 235 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];
             }
         }

PetscErrorCode Nektar::MultiRegions::GlobalLinSysPETSc::DoDestroyMatCtx ( Mat M )

staticprivate

Destroy matrix shell context object.

Note the matrix shell and preconditioner share a common context so this might have already been deallocated below, in which case we do nothing.

Parameters

M	Matrix shell object

Definition at line 486 of file GlobalLinSysPETSc.cpp.

Referenced by SetUpMatVec().

         {
             void *ptr;
             MatShellGetContext(M, &ptr);
             ShellCtx *ctx = (ShellCtx *)ptr;
             delete ctx;
             return 0;
         }

PetscErrorCode Nektar::MultiRegions::GlobalLinSysPETSc::DoDestroyPCCtx ( PC pc )

staticprivate

Destroy preconditioner context object.

Note the matrix shell and preconditioner share a common context so this might have already been deallocated above, in which case we do nothing.

Parameters

pc	Preconditioner object

Definition at line 504 of file GlobalLinSysPETSc.cpp.

Referenced by SetUpMatVec().

         {
             void *ptr;
             PCShellGetContext(pc, &ptr);
             ShellCtx *ctx = (ShellCtx *)ptr;
             delete ctx;
             return 0;
         }

PetscErrorCode Nektar::MultiRegions::GlobalLinSysPETSc::DoMatrixMultiply	(	Mat	M,
		Vec	in,
		Vec	out
	)

staticprivate

Perform matrix multiplication using Nektar++ routines.

This static function uses Nektar++ routines to calculate the matrix-vector product of M with in, storing the output in out.

Parameters

M	Original MatShell matrix, which stores the ShellCtx object.
in	Input vector.
out	Output vector.

Definition at line 438 of file GlobalLinSysPETSc.cpp.

References DoNekppOperation().

Referenced by SetUpMatVec().

         {
             // Grab our shell context from M.
             void *ptr;
             MatShellGetContext(M, &ptr);
             ShellCtx *ctx = (ShellCtx *)ptr;
 
             DoNekppOperation(in, out, ctx, false);
 
             // Must return 0, otherwise PETSc complains.
             return 0;
         }

void Nektar::MultiRegions::GlobalLinSysPETSc::DoNekppOperation	(	Vec &	in,
		Vec &	out,
		ShellCtx *	ctx,
		bool	precon
	)

staticprivate

Perform either matrix multiplication or preconditioning using Nektar++ routines.

This static function uses Nektar++ routines to calculate the matrix-vector product of M with in, storing the output in out.

Todo:: There's a lot of scatters and copies that might possibly be eliminated to make this more efficient.

Parameters

in	Input vector.
out	Output vector.
ctx	ShellCtx object that points to our instance of GlobalLinSysPETSc.
precon	If true, we apply a preconditioner, if false, we perform a matrix multiplication.

Definition at line 384 of file GlobalLinSysPETSc.cpp.

References Nektar::MultiRegions::GlobalLinSysPETSc::ShellCtx::linSys, m_ctx, m_locVec, m_precon, m_reorderedMap, Nektar::MultiRegions::GlobalLinSysPETSc::ShellCtx::nDir, Nektar::MultiRegions::GlobalLinSysPETSc::ShellCtx::nGlobal, v_DoMatrixMultiply(), and Vmath::Vcopy().

Referenced by DoMatrixMultiply(), and DoPreconditioner().

         {
             const int          nGlobal  = ctx->nGlobal;
             const int          nDir     = ctx->nDir;
             const int          nHomDofs = nGlobal - nDir;
             GlobalLinSysPETSc *linSys   = ctx->linSys;
 
             // Scatter from PETSc ordering to our local ordering. It's actually
             // unclear whether this step might also do some communication in
             // parallel, which is probably not ideal.
             VecScatterBegin(linSys->m_ctx, in, linSys->m_locVec,
                             INSERT_VALUES, SCATTER_FORWARD);
             VecScatterEnd  (linSys->m_ctx, in, linSys->m_locVec,
                             INSERT_VALUES, SCATTER_FORWARD);
 
             // Temporary storage to pass to Nektar++
             Array<OneD, NekDouble> tmpIn(nHomDofs), tmpOut(nHomDofs);
 
             // Get values from input vector and copy to tmpIn.
             PetscScalar *tmpLocIn;
             VecGetArray    (linSys->m_locVec, &tmpLocIn);
             Vmath::Vcopy   (nHomDofs, tmpLocIn, 1, &tmpIn[0], 1);
             VecRestoreArray(linSys->m_locVec, &tmpLocIn);
 
             // Do matrix multiply in Nektar++, store in tmpOut.
             if (precon)
             {
                 linSys->m_precon->DoPreconditioner(tmpIn, tmpOut);
             }
             else
             {
                 linSys->v_DoMatrixMultiply(tmpIn, tmpOut);
             }
 
             // Scatter back to PETSc ordering and put in out.
             VecSetValues(out, nHomDofs, &linSys->m_reorderedMap[0],
                          &tmpOut[0], INSERT_VALUES);
             VecAssemblyBegin(out);
             VecAssemblyEnd  (out);
         }

PetscErrorCode Nektar::MultiRegions::GlobalLinSysPETSc::DoPreconditioner	(	PC	pc,
		Vec	in,
		Vec	out
	)

staticprivate

Apply preconditioning using Nektar++ routines.

This static function uses Nektar++ routines to apply the preconditioner stored in GlobalLinSysPETSc::m_precon from the context of pc to the vector in, storing the output in out.

Parameters

pc	Preconditioner object that stores the ShellCtx.
in	Input vector.
out	Output vector.

Definition at line 463 of file GlobalLinSysPETSc.cpp.

References DoNekppOperation().

Referenced by SetUpMatVec().

         {
             // Grab our PCShell context from pc.
             void *ptr;
             PCShellGetContext(pc, &ptr);
             ShellCtx *ctx = (ShellCtx *)ptr;
 
             DoNekppOperation(in, out, ctx, true);
 
             // Must return 0, otherwise PETSc complains.
             return 0;
         }

void Nektar::MultiRegions::GlobalLinSysPETSc::SetUpMatVec	(	int	nGlobal,
		int	nDir
	)

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.

Parameters

nGlobal	Number of global degrees of freedom in the system (on this processor)
nDir	Number of Dirichlet degrees of freedom (on this processor).

Definition at line 312 of file GlobalLinSysPETSc.cpp.

References DoDestroyMatCtx(), DoDestroyPCCtx(), DoMatrixMultiply(), DoPreconditioner(), Nektar::MultiRegions::ePETScMatMultShell, Nektar::MultiRegions::GlobalLinSysPETSc::ShellCtx::linSys, m_b, m_matMult, m_matrix, m_nLocal, m_pc, m_x, Nektar::MultiRegions::GlobalLinSysPETSc::ShellCtx::nDir, and Nektar::MultiRegions::GlobalLinSysPETSc::ShellCtx::nGlobal.

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
             if (m_matMult == ePETScMatMultShell)
             {
                 // Create ShellCtx context object which will store the matrix
                 // size and a pointer to the linear system. We do this so that
                 // we can call a member function to the matrix-vector and
                 // preconditioning multiplication in a subclass.
                 ShellCtx *ctx1 = new ShellCtx(), *ctx2 = new ShellCtx();
                 ctx1->nGlobal = ctx2->nGlobal = nGlobal;
                 ctx1->nDir    = ctx2->nDir    = nDir;
                 ctx1->linSys  = ctx2->linSys  = this;
 
                 // Set up MatShell object.
                 MatCreateShell      (PETSC_COMM_WORLD, m_nLocal, m_nLocal,
                                      PETSC_DETERMINE, PETSC_DETERMINE,
                                      (void *)ctx1, &m_matrix);
                 MatShellSetOperation(m_matrix, MATOP_MULT,
                                      (void(*)(void))DoMatrixMultiply);
                 MatShellSetOperation(m_matrix, MATOP_DESTROY,
                                      (void(*)(void))DoDestroyMatCtx);
 
                 // Create a PCShell to go alongside the MatShell.
                 PCCreate         (PETSC_COMM_WORLD, &m_pc);
 #if PETSC_VERSION_GE(3,5,0)
                 PCSetOperators   (m_pc, m_matrix, m_matrix);
 #else
                 PCSetOperators   (m_pc, m_matrix, m_matrix, SAME_NONZERO_PATTERN);
 #endif
                 PCSetType        (m_pc, PCSHELL);
                 PCShellSetApply  (m_pc, DoPreconditioner);
                 PCShellSetDestroy(m_pc, DoDestroyPCCtx);
                 PCShellSetContext(m_pc, ctx2);
             }
             else
             {
                 // Otherwise we create a PETSc matrix and use MatSetFromOptions
                 // so that we can set various options on the command line.
                 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 200 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.

Parameters

tolerance Residual tolerance to converge to.

Definition at line 521 of file GlobalLinSysPETSc.cpp.

References Nektar::MultiRegions::ePETScMatMultShell, m_ksp, m_matMult, m_matrix, and m_pc.

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
 
             if (m_matMult == ePETScMatMultShell)
             {
                 KSPSetPC(m_ksp, m_pc);
             }
         }

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

protectedpure virtual

Implemented in Nektar::MultiRegions::GlobalLinSysPETScStaticCond, and Nektar::MultiRegions::GlobalLinSysPETScFull.

Referenced by DoNekppOperation().

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

References ASSERTL0, Nektar::MultiRegions::GlobalLinSys::CreatePrecon(), Nektar::MultiRegions::ePETScMatMultShell, m_b, m_ctx, m_ksp, m_locVec, m_matMult, m_precon, m_reorderedMap, m_x, and Vmath::Vcopy().

         {
             const int nHomDofs = pNumRows - pNumDir;
 
             if (!m_precon && m_matMult == ePETScMatMultShell)
             {
                 m_precon = CreatePrecon(locToGloMap);
                 m_precon->BuildPreconditioner();
             }
 
             // 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);
 
             KSPConvergedReason reason;
             KSPGetConvergedReason(m_ksp, &reason);
             ASSERTL0(reason > 0,
                     "PETSc solver diverged, reason is: " +
                         std::string(KSPConvergedReasons[reason]));
 
 
             // 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);
         }

Member Data Documentation

Vec Nektar::MultiRegions::GlobalLinSysPETSc::m_b

protected

Definition at line 81 of file GlobalLinSysPETSc.h.

Referenced by SetUpMatVec(), v_SolveLinearSystem(), and ~GlobalLinSysPETSc().

VecScatter Nektar::MultiRegions::GlobalLinSysPETSc::m_ctx

protected

PETSc scatter context that takes us between Nektar++ global ordering and PETSc vector ordering.

Definition at line 93 of file GlobalLinSysPETSc.h.

Referenced by DoNekppOperation(), SetUpScatter(), and v_SolveLinearSystem().

KSP Nektar::MultiRegions::GlobalLinSysPETSc::m_ksp

protected

KSP object that represents solver system.

Definition at line 83 of file GlobalLinSysPETSc.h.

Referenced by SetUpSolver(), v_SolveLinearSystem(), and ~GlobalLinSysPETSc().

Vec Nektar::MultiRegions::GlobalLinSysPETSc::m_locVec

protected

Definition at line 81 of file GlobalLinSysPETSc.h.

Referenced by DoNekppOperation(), SetUpScatter(), v_SolveLinearSystem(), and ~GlobalLinSysPETSc().

PETScMatMult Nektar::MultiRegions::GlobalLinSysPETSc::m_matMult

protected

Enumerator to select matrix multiplication type.

Definition at line 87 of file GlobalLinSysPETSc.h.

Referenced by GlobalLinSysPETSc(), SetUpMatVec(), SetUpSolver(), Nektar::MultiRegions::GlobalLinSysPETScStaticCond::v_AssembleSchurComplement(), and v_SolveLinearSystem().

Mat Nektar::MultiRegions::GlobalLinSysPETSc::m_matrix

protected

PETSc matrix object.

Definition at line 79 of file GlobalLinSysPETSc.h.

Referenced by GlobalLinSysPETSc(), Nektar::MultiRegions::GlobalLinSysPETScFull::GlobalLinSysPETScFull(), SetUpMatVec(), SetUpSolver(), Nektar::MultiRegions::GlobalLinSysPETScStaticCond::v_AssembleSchurComplement(), and ~GlobalLinSysPETSc().

int Nektar::MultiRegions::GlobalLinSysPETSc::m_nLocal

protected

Number of unique degrees of freedom on this process.

Definition at line 95 of file GlobalLinSysPETSc.h.

Referenced by CalculateReordering(), and SetUpMatVec().

PC Nektar::MultiRegions::GlobalLinSysPETSc::m_pc

protected

PCShell for preconditioner.

Definition at line 85 of file GlobalLinSysPETSc.h.

Referenced by SetUpMatVec(), SetUpSolver(), and ~GlobalLinSysPETSc().

PreconditionerSharedPtr Nektar::MultiRegions::GlobalLinSysPETSc::m_precon

protected

Definition at line 97 of file GlobalLinSysPETSc.h.

Referenced by DoNekppOperation(), Nektar::MultiRegions::GlobalLinSysPETScStaticCond::GlobalLinSysPETScStaticCond(), Nektar::MultiRegions::GlobalLinSysPETScStaticCond::v_AssembleSchurComplement(), Nektar::MultiRegions::GlobalLinSysPETScStaticCond::v_BasisInvTransform(), Nektar::MultiRegions::GlobalLinSysPETScStaticCond::v_BasisTransform(), Nektar::MultiRegions::GlobalLinSysPETScStaticCond::v_InitObject(), Nektar::MultiRegions::GlobalLinSysPETScStaticCond::v_PreSolve(), Nektar::MultiRegions::GlobalLinSysPETScStaticCond::v_Recurse(), and v_SolveLinearSystem().

std::vector<int> Nektar::MultiRegions::GlobalLinSysPETSc::m_reorderedMap

protected

Reordering that takes universal IDs to a unique row in the PETSc matrix.

See also: GlobalLinSysPETSc::CalculateReordering

Definition at line 90 of file GlobalLinSysPETSc.h.

Referenced by CalculateReordering(), DoNekppOperation(), Nektar::MultiRegions::GlobalLinSysPETScFull::GlobalLinSysPETScFull(), SetUpScatter(), Nektar::MultiRegions::GlobalLinSysPETScStaticCond::v_AssembleSchurComplement(), and v_SolveLinearSystem().

Vec Nektar::MultiRegions::GlobalLinSysPETSc::m_x

protected

PETSc vector objects used for local storage.

Definition at line 81 of file GlobalLinSysPETSc.h.

Referenced by SetUpMatVec(), SetUpScatter(), v_SolveLinearSystem(), and ~GlobalLinSysPETSc().

std::string Nektar::MultiRegions::GlobalLinSysPETSc::matMult

staticprivate

Initial value:

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

Definition at line 134 of file GlobalLinSysPETSc.h.

std::string Nektar::MultiRegions::GlobalLinSysPETSc::matMultIds

staticprivate

Initial value:

= {
            LibUtilities::SessionReader::RegisterEnumValue(
                "PETScMatMult", "Sparse",
                MultiRegions::ePETScMatMultSparse),
            LibUtilities::SessionReader::RegisterEnumValue(
                "PETScMatMult", "Shell",
                MultiRegions::ePETScMatMultShell)
        }

Definition at line 135 of file GlobalLinSysPETSc.h.

Classes

Public Member Functions

Protected Member Functions

Protected Attributes

Static Private Member Functions

Static Private Attributes

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation