A PETSc global linear system. More...

#include <GlobalLinSysPETSc.h>

Inheritance 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 std::weak_ptr< ExpList > &pExp, const std::shared_ptr< AssemblyMap > &pLocToGloMap)
	Constructor for full direct matrix solve. More...

	~GlobalLinSysPETSc () override
	Clean up PETSc objects. More...

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

virtual	~GlobalLinSys ()

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

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

void	InitObject ()

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

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

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

int	GetNumBlocks ()

DNekScalMatSharedPtr	GetBlock (unsigned int n)

void	DropBlock (unsigned int n)

DNekScalBlkMatSharedPtr	GetStaticCondBlock (unsigned int n)

void	DropStaticCondBlock (unsigned int n)

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

Protected Member Functions
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...

void	v_SolveLinearSystem (const int pNumRows, const Array< OneD, const NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const AssemblyMapSharedPtr &locToGloMap, const int pNumDir) override
	Solve linear system using 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 void	v_Solve (const Array< OneD, const NekDouble > &in, Array< OneD, NekDouble > &out, const AssemblyMapSharedPtr &locToGloMap, const Array< OneD, const NekDouble > &dirForcing=NullNekDouble1DArray)=0
	Solve a linear system based on mapping. More...

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

virtual void	v_InitObject ()

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

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

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

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

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

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

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

Protected Attributes
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 std::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 56 of file GlobalLinSysPETSc.h.

Constructor & Destructor Documentation

◆ GlobalLinSysPETSc()

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

Constructor for full direct matrix solve.

Definition at line 59 of file GlobalLinSysPETSc.cpp.

    : 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 =
                std::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);
}

References Nektar::MultiRegions::GlobalLinSys::m_expList, m_matMult, and m_matrix.

◆ ~GlobalLinSysPETSc()

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

override

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

{
    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);
    }
}

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

Member Function Documentation

◆ CalculateReordering()

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

{
    LibUtilities::CommSharedPtr vComm =
        m_expList.lock()->GetSession()->GetComm();
 
    const int nDirDofs = pLocToGloMap->GetNumGlobalDirBndCoeffs();
    const int nHomDofs = glo2uniMap.size() - 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];
    }
}

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

◆ DoDestroyMatCtx()

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

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

Referenced by SetUpMatVec().

◆ DoDestroyPCCtx()

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

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

Referenced by SetUpMatVec().

◆ DoMatrixMultiply()

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

{
    // 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;
}

References DoNekppOperation().

Referenced by SetUpMatVec().

◆ DoNekppOperation()

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

{
    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);
}

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

◆ DoPreconditioner()

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

{
    // 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;
}

References DoNekppOperation().

Referenced by SetUpMatVec().

◆ SetUpMatVec()

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

{
    // 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);
    }
}

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

◆ SetUpScatter()

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

{
    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);
}

References m_ctx, m_locVec, m_reorderedMap, and m_x.

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

◆ SetUpSolver()

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

{
    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);
    }
}

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

◆ v_DoMatrixMultiply()

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

protectedpure virtual

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

Referenced by DoNekppOperation().

◆ v_SolveLinearSystem()

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
	)

overrideprotectedvirtual

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.

Reimplemented in Nektar::MultiRegions::GlobalLinSysPETScStaticCond.

Definition at line 139 of file GlobalLinSysPETSc.cpp.

{
    const int nHomDofs = pNumRows - pNumDir;
 
    if (!m_precon && m_matMult == ePETScMatMultShell)
    {
        m_precon = CreatePrecon(locToGloMap);
        m_precon->BuildPreconditioner();
    }
 
    Array<OneD, NekDouble> Glo(pNumRows);
    locToGloMap->Assemble(pInput, Glo);
 
    // Populate RHS vector from input
    VecSetValues(m_b, nHomDofs, &m_reorderedMap[0], &Glo[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, &Glo[pNumDir], 1);
    Vmath::Zero(pNumDir, Glo, 1);
    locToGloMap->GlobalToLocal(Glo, pOutput);
    VecRestoreArray(m_locVec, &tmp);
}

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, Vmath::Vcopy(), and Vmath::Zero().

Member Data Documentation

◆ m_b

Vec Nektar::MultiRegions::GlobalLinSysPETSc::m_b

protected

Definition at line 70 of file GlobalLinSysPETSc.h.

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

◆ m_ctx

VecScatter Nektar::MultiRegions::GlobalLinSysPETSc::m_ctx

protected

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

Definition at line 82 of file GlobalLinSysPETSc.h.

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

◆ m_ksp

KSP Nektar::MultiRegions::GlobalLinSysPETSc::m_ksp

protected

KSP object that represents solver system.

Definition at line 72 of file GlobalLinSysPETSc.h.

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

◆ m_locVec

Vec Nektar::MultiRegions::GlobalLinSysPETSc::m_locVec

protected

Definition at line 70 of file GlobalLinSysPETSc.h.

Referenced by DoNekppOperation(), SetUpScatter(), v_SolveLinearSystem(), Nektar::MultiRegions::GlobalLinSysPETScStaticCond::v_SolveLinearSystem(), and ~GlobalLinSysPETSc().

◆ m_matMult

PETScMatMult Nektar::MultiRegions::GlobalLinSysPETSc::m_matMult

protected

Enumerator to select matrix multiplication type.

Definition at line 76 of file GlobalLinSysPETSc.h.

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

◆ m_matrix

Mat Nektar::MultiRegions::GlobalLinSysPETSc::m_matrix

protected

PETSc matrix object.

Definition at line 68 of file GlobalLinSysPETSc.h.

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

◆ m_nLocal

int Nektar::MultiRegions::GlobalLinSysPETSc::m_nLocal

protected

Number of unique degrees of freedom on this process.

Definition at line 84 of file GlobalLinSysPETSc.h.

Referenced by CalculateReordering(), and SetUpMatVec().

◆ m_pc

PC Nektar::MultiRegions::GlobalLinSysPETSc::m_pc

protected

PCShell for preconditioner.

Definition at line 74 of file GlobalLinSysPETSc.h.

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

◆ m_precon

PreconditionerSharedPtr Nektar::MultiRegions::GlobalLinSysPETSc::m_precon

protected

Definition at line 86 of file GlobalLinSysPETSc.h.

◆ m_reorderedMap

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 79 of file GlobalLinSysPETSc.h.

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

◆ m_x

Vec Nektar::MultiRegions::GlobalLinSysPETSc::m_x

protected

PETSc vector objects used for local storage.

Definition at line 70 of file GlobalLinSysPETSc.h.

Referenced by SetUpMatVec(), SetUpScatter(), v_SolveLinearSystem(), Nektar::MultiRegions::GlobalLinSysPETScStaticCond::v_SolveLinearSystem(), and ~GlobalLinSysPETSc().

◆ matMult

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

staticprivate

Initial value:

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

Definition at line 128 of file GlobalLinSysPETSc.h.

◆ matMultIds

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

◆ GlobalLinSysPETSc()

◆ ~GlobalLinSysPETSc()

Member Function Documentation

◆ CalculateReordering()

◆ DoDestroyMatCtx()

◆ DoDestroyPCCtx()

◆ DoMatrixMultiply()

◆ DoNekppOperation()

◆ DoPreconditioner()

◆ SetUpMatVec()

◆ SetUpScatter()

◆ SetUpSolver()

◆ v_DoMatrixMultiply()

◆ v_SolveLinearSystem()

Member Data Documentation

◆ m_b

◆ m_ctx

◆ m_ksp

◆ m_locVec

◆ m_matMult

◆ m_matrix

◆ m_nLocal

◆ m_pc

◆ m_precon

◆ m_reorderedMap

◆ m_x

◆ matMult

◆ matMultIds