40 #include "petscversion.h"
46 namespace MultiRegions
48 std::string GlobalLinSysPETSc::matMult =
49 LibUtilities::SessionReader::RegisterDefaultSolverInfo(
50 "PETScMatMult",
"Sparse");
51 std::string GlobalLinSysPETSc::matMultIds[] = {
52 LibUtilities::SessionReader::RegisterEnumValue(
53 "PETScMatMult",
"Sparse",
55 LibUtilities::SessionReader::RegisterEnumValue(
56 "PETScMatMult",
"Shell",
65 GlobalLinSysPETSc::GlobalLinSysPETSc(
67 const boost::weak_ptr<ExpList> &pExp,
68 const boost::shared_ptr<AssemblyMap> &pLocToGloMap)
80 PetscBool isInitialized;
81 PetscInitialized(&isInitialized);
84 PetscInitializeNoArguments();
88 MatCreate(PETSC_COMM_WORLD, &
m_matrix);
101 PetscBool isFinalized;
102 PetscFinalized(&isFinalized);
109 #ifdef NEKTAR_USE_MPI
111 MPI_Finalized(&mpiFinal);
112 isFinalized = isFinalized || mpiFinal ? PETSC_TRUE : PETSC_FALSE;
143 const int nHomDofs = pNumRows - pNumDir;
153 &pInput[pNumDir], INSERT_VALUES);
156 VecAssemblyBegin(
m_b);
157 VecAssemblyEnd (
m_b);
164 INSERT_VALUES, SCATTER_FORWARD);
166 INSERT_VALUES, SCATTER_FORWARD);
187 IS isGlobal, isLocal;
189 PETSC_COPY_VALUES, &isGlobal);
190 ISCreateStride (PETSC_COMM_SELF, nHomDofs, 0, 1, &isLocal);
193 VecCreate (PETSC_COMM_SELF, &
m_locVec);
194 VecSetSizes (
m_locVec, nHomDofs, PETSC_DECIDE);
201 ISDestroy(&isGlobal);
223 =
m_expList.lock()->GetSession()->GetComm();
225 const int nDirDofs = pLocToGloMap->GetNumGlobalDirBndCoeffs();
226 const int nHomDofs = glo2uniMap.num_elements() - nDirDofs;
227 const int nProc = vComm->GetSize();
228 const int rank = vComm->GetRank();
238 localCounts[rank] = nHomDofs;
241 for (n = 1; n < nProc; ++n)
243 localOffset[n] = localOffset[n-1] + localCounts[n-1];
246 int totHomDofs =
Vmath::Vsum(nProc, localCounts, 1);
247 vector<unsigned int> allUniIds(totHomDofs, 0);
250 for (n = 0; n < nHomDofs; ++n)
252 int gid = n + nDirDofs;
253 allUniIds[n + localOffset[rank]] = glo2uniMap[gid];
259 std::sort(allUniIds.begin(), allUniIds.end());
260 map<int,int> uniIdReorder;
263 for (cnt = n = 0; n < allUniIds.size(); ++n)
265 if (uniIdReorder.count(allUniIds[n]) > 0)
270 uniIdReorder[allUniIds[n]] = cnt++;
274 for (n = 0; n < nHomDofs; ++n)
276 int gid = n + nDirDofs;
277 int uniId = glo2uniMap[gid];
278 ASSERTL0(uniIdReorder.count(uniId) > 0,
"Error in ordering");
297 VecCreate (PETSC_COMM_WORLD, &
m_x);
299 VecSetFromOptions(
m_x);
310 ctx1->
nGlobal = ctx2->nGlobal = nGlobal;
311 ctx1->
nDir = ctx2->nDir = nDir;
312 ctx1->
linSys = ctx2->linSys =
this;
316 PETSC_DETERMINE, PETSC_DETERMINE,
318 MatShellSetOperation(
m_matrix, MATOP_MULT,
320 MatShellSetOperation(
m_matrix, MATOP_DESTROY,
324 PCCreate (PETSC_COMM_WORLD, &
m_pc);
325 #if PETSC_VERSION_GE(3,5,0)
330 PCSetType (
m_pc, PCSHELL);
333 PCShellSetContext(
m_pc, ctx2);
339 MatCreate (PETSC_COMM_WORLD, &
m_matrix);
342 PETSC_DETERMINE, PETSC_DETERMINE);
367 Vec &in, Vec &out,
ShellCtx *ctx,
bool precon)
369 const int nGlobal = ctx->
nGlobal;
370 const int nDir = ctx->
nDir;
371 const int nHomDofs = nGlobal - nDir;
378 INSERT_VALUES, SCATTER_FORWARD);
380 INSERT_VALUES, SCATTER_FORWARD);
386 PetscScalar *tmpLocIn;
387 VecGetArray (linSys->
m_locVec, &tmpLocIn);
389 VecRestoreArray(linSys->
m_locVec, &tmpLocIn);
394 linSys->
m_precon->DoPreconditioner(tmpIn, tmpOut);
403 &tmpOut[0], INSERT_VALUES);
404 VecAssemblyBegin(out);
405 VecAssemblyEnd (out);
421 Mat M, Vec in, Vec out)
425 MatShellGetContext(M, &ptr);
446 PC pc, Vec in, Vec out)
450 PCShellGetContext(pc, &ptr);
471 MatShellGetContext(M, &ptr);
489 PCShellGetContext(pc, &ptr);
505 KSPCreate(PETSC_COMM_WORLD, &
m_ksp);
507 m_ksp, tolerance, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT);
508 KSPSetFromOptions(
m_ksp);
509 #if PETSC_VERSION_GE(3,5,0)
VecScatter m_ctx
PETSc scatter context that takes us between Nektar++ global ordering and PETSc vector ordering...
#define ASSERTL0(condition, msg)
GlobalLinSysPETSc * linSys
Pointer to the original calling object.
PreconditionerSharedPtr m_precon
boost::shared_ptr< AssemblyMap > AssemblyMapSharedPtr
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.
PC m_pc
PCShell for preconditioner.
void SetUpSolver(NekDouble tolerance)
Set up KSP solver object.
Internal struct for MatShell and PCShell calls to store current context for callback.
Vec m_x
PETSc vector objects used for local storage.
static PetscErrorCode DoPreconditioner(PC pc, Vec in, Vec out)
Apply preconditioning using Nektar++ routines.
int nGlobal
Number of global degrees of freedom.
boost::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Mat m_matrix
PETSc matrix object.
int m_nLocal
Number of unique degrees of freedom on this process.
KSP m_ksp
KSP object that represents solver system.
virtual void v_DoMatrixMultiply(const Array< OneD, const NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)=0
virtual ~GlobalLinSysPETSc()
Clean up PETSc objects.
static PetscErrorCode DoMatrixMultiply(Mat M, Vec in, Vec out)
Perform matrix multiplication using Nektar++ routines.
Describe a linear system.
A PETSc global linear system.
void SetUpScatter()
Set up PETSc local (equivalent to Nektar++ global) and global (equivalent to universal) scatter maps...
int nDir
Number of Dirichlet degrees of freedom.
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.
static PetscErrorCode DoDestroyMatCtx(Mat M)
Destroy matrix shell context object.
void SetUpMatVec(int nGlobal, int nDir)
Construct PETSc matrix and vector handles.
static void DoNekppOperation(Vec &in, Vec &out, ShellCtx *ctx, bool precon)
Perform either matrix multiplication or preconditioning using Nektar++ routines.
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
std::vector< int > m_reorderedMap
Reordering that takes universal IDs to a unique row in the PETSc matrix.
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
static PetscErrorCode DoDestroyPCCtx(PC pc)
Destroy preconditioner context object.
PETScMatMult m_matMult
Enumerator to select matrix multiplication type.
PreconditionerSharedPtr CreatePrecon(AssemblyMapSharedPtr asmMap)
Create a preconditioner object from the parameters defined in the supplied assembly map...
const boost::weak_ptr< ExpList > m_expList
Local Matrix System.