Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
GlobalLinSysPETSc.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File GlobalLinSys.cpp
4 //
5 // For more information, please see: http://www.nektar.info
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: GlobalLinSys definition
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
38 
39 #include "petscis.h"
40 #include "petscversion.h"
41 
42 using namespace std;
43 
44 namespace Nektar
45 {
46  namespace MultiRegions
47  {
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",
58  };
59 
60  /**
61  * @class GlobalLinSysPETSc
62  *
63  * Solves a linear system using PETSc.
64  */
65  GlobalLinSysPETSc::GlobalLinSysPETSc(
66  const GlobalLinSysKey &pKey,
67  const boost::weak_ptr<ExpList> &pExp,
68  const boost::shared_ptr<AssemblyMap> &pLocToGloMap)
69  : GlobalLinSys(pKey, pExp, pLocToGloMap)
70  {
71  // Determine whether to use standard sparse matrix approach or
72  // shell.
73  m_matMult = pExp.lock()->GetSession()
74  ->GetSolverInfoAsEnum<PETScMatMult>(
75  "PETScMatMult");
76 
77  // Check PETSc is initialized. For some reason, this is needed on
78  // OS X as logging is not initialized properly in the call within
79  // CommMpi.
80  PetscBool isInitialized;
81  PetscInitialized(&isInitialized);
82  if (!isInitialized)
83  {
84 #ifdef NEKTAR_USE_MPI
85  std::string commType =
86  m_expList.lock()->GetSession()->GetComm()->GetType();
87  if (commType.find("MPI") != std::string::npos)
88  {
90  boost::static_pointer_cast<LibUtilities::CommMpi>(
91  m_expList.lock()->GetSession()->GetComm());
92  PETSC_COMM_WORLD = comm->GetComm();
93  }
94 #endif
95  PetscInitializeNoArguments();
96  }
97 
98  // Create matrix
99  MatCreate(PETSC_COMM_WORLD, &m_matrix);
100  }
101 
102  /**
103  * @brief Clean up PETSc objects.
104  *
105  * Note that if SessionReader::Finalize is called before the end of the
106  * program, PETSc may have been finalized already, at which point we
107  * cannot deallocate our objects. If that's the case we do nothing and
108  * let the kernel clear up after us.
109  */
111  {
112  PetscBool isFinalized;
113  PetscFinalized(&isFinalized);
114 
115  // Sometimes, PetscFinalized returns false when (in fact) CommMpi's
116  // Finalise routine has been called. We therefore also need to check
117  // whether MPI has been finalised. This might arise from the
118  // additional call to PetscInitializeNoArguments in the constructor
119  // above.
120 #ifdef NEKTAR_USE_MPI
121  int mpiFinal = 0;
122  MPI_Finalized(&mpiFinal);
123  isFinalized = isFinalized || mpiFinal ? PETSC_TRUE : PETSC_FALSE;
124 #endif
125 
126  if (!isFinalized)
127  {
128  KSPDestroy(&m_ksp);
129  PCDestroy (&m_pc);
130  MatDestroy(&m_matrix);
131  VecDestroy(&m_x);
132  VecDestroy(&m_b);
133  VecDestroy(&m_locVec);
134  }
135  }
136 
137  /**
138  * @brief Solve linear system using PETSc.
139  *
140  * The general strategy being a PETSc solve is to:
141  *
142  * - Copy values into the PETSc vector #m_b
143  * - Solve the system #m_ksp and place result into #m_x.
144  * - Scatter results back into #m_locVec using #m_ctx scatter object.
145  * - Copy from #m_locVec to output array #pOutput.
146  */
148  const int pNumRows,
149  const Array<OneD,const NekDouble> &pInput,
150  Array<OneD, NekDouble> &pOutput,
151  const AssemblyMapSharedPtr &locToGloMap,
152  const int pNumDir)
153  {
154  const int nHomDofs = pNumRows - pNumDir;
155 
157  {
158  m_precon = CreatePrecon(locToGloMap);
159  m_precon->BuildPreconditioner();
160  }
161 
162  // Populate RHS vector from input
163  VecSetValues(m_b, nHomDofs, &m_reorderedMap[0],
164  &pInput[pNumDir], INSERT_VALUES);
165 
166  // Assemble RHS vector
167  VecAssemblyBegin(m_b);
168  VecAssemblyEnd (m_b);
169 
170  // Do system solve
171  KSPSolve(m_ksp, m_b, m_x);
172 
173  KSPConvergedReason reason;
174  KSPGetConvergedReason(m_ksp, &reason);
175  ASSERTL0(reason > 0,
176  "PETSc solver diverged, reason is: " +
177  std::string(KSPConvergedReasons[reason]));
178 
179 
180  // Scatter results to local vector
181  VecScatterBegin(m_ctx, m_x, m_locVec,
182  INSERT_VALUES, SCATTER_FORWARD);
183  VecScatterEnd (m_ctx, m_x, m_locVec,
184  INSERT_VALUES, SCATTER_FORWARD);
185 
186  // Copy results into output vector
187  PetscScalar *tmp;
188  VecGetArray (m_locVec, &tmp);
189  Vmath::Vcopy (nHomDofs, tmp, 1, &pOutput[pNumDir], 1);
190  VecRestoreArray(m_locVec, &tmp);
191  }
192 
193  /**
194  * @brief Set up PETSc local (equivalent to Nektar++ global) and global
195  * (equivalent to universal) scatter maps.
196  *
197  * These maps are used in GlobalLinSysPETSc::v_SolveLinearSystem to
198  * scatter the solution vector back to each process.
199  */
201  {
202  const int nHomDofs = m_reorderedMap.size();
203 
204  // Create local and global numbering systems for vector
205  IS isGlobal, isLocal;
206  ISCreateGeneral (PETSC_COMM_SELF, nHomDofs, &m_reorderedMap[0],
207  PETSC_COPY_VALUES, &isGlobal);
208  ISCreateStride (PETSC_COMM_SELF, nHomDofs, 0, 1, &isLocal);
209 
210  // Create local vector for output
211  VecCreate (PETSC_COMM_SELF, &m_locVec);
212  VecSetSizes (m_locVec, nHomDofs, PETSC_DECIDE);
213  VecSetFromOptions(m_locVec);
214 
215  // Create scatter context
216  VecScatterCreate (m_x, isGlobal, m_locVec, isLocal, &m_ctx);
217 
218  // Clean up
219  ISDestroy(&isGlobal);
220  ISDestroy(&isLocal);
221  }
222 
223  /**
224  * @brief Calculate a reordering of universal IDs for PETSc.
225  *
226  * PETSc requires a unique, contiguous index of all global and universal
227  * degrees of freedom which represents its position inside the
228  * matrix. Presently Gs does not guarantee this, so this routine
229  * constructs a new universal mapping.
230  *
231  * @param glo2uniMap Global to universal map
232  * @param glo2unique Global to unique map
233  * @param pLocToGloMap Assembly map for this system
234  */
236  const Array<OneD, const int> &glo2uniMap,
237  const Array<OneD, const int> &glo2unique,
238  const AssemblyMapSharedPtr &pLocToGloMap)
239  {
241  = m_expList.lock()->GetSession()->GetComm();
242 
243  const int nDirDofs = pLocToGloMap->GetNumGlobalDirBndCoeffs();
244  const int nHomDofs = glo2uniMap.num_elements() - nDirDofs;
245  const int nProc = vComm->GetSize();
246  const int rank = vComm->GetRank();
247 
248  int n, cnt;
249 
250  // Count number of unique degrees of freedom on each process.
251  m_nLocal = Vmath::Vsum(nHomDofs, glo2unique + nDirDofs, 1);
252  m_reorderedMap.resize(nHomDofs);
253 
254  // Reduce coefficient counts across all processors.
255  Array<OneD, int> localCounts(nProc, 0), localOffset(nProc, 0);
256  localCounts[rank] = nHomDofs;
257  vComm->AllReduce(localCounts, LibUtilities::ReduceSum);
258 
259  for (n = 1; n < nProc; ++n)
260  {
261  localOffset[n] = localOffset[n-1] + localCounts[n-1];
262  }
263 
264  int totHomDofs = Vmath::Vsum(nProc, localCounts, 1);
265  vector<unsigned int> allUniIds(totHomDofs, 0);
266 
267  // Assemble list of universal IDs
268  for (n = 0; n < nHomDofs; ++n)
269  {
270  int gid = n + nDirDofs;
271  allUniIds[n + localOffset[rank]] = glo2uniMap[gid];
272  }
273 
274  // Reduce this across processors so that each process has a list of
275  // all universal IDs.
276  vComm->AllReduce(allUniIds, LibUtilities::ReduceSum);
277  std::sort(allUniIds.begin(), allUniIds.end());
278  map<int,int> uniIdReorder;
279 
280  // Renumber starting from 0.
281  for (cnt = n = 0; n < allUniIds.size(); ++n)
282  {
283  if (uniIdReorder.count(allUniIds[n]) > 0)
284  {
285  continue;
286  }
287 
288  uniIdReorder[allUniIds[n]] = cnt++;
289  }
290 
291  // Populate reordering map.
292  for (n = 0; n < nHomDofs; ++n)
293  {
294  int gid = n + nDirDofs;
295  int uniId = glo2uniMap[gid];
296  ASSERTL0(uniIdReorder.count(uniId) > 0, "Error in ordering");
297  m_reorderedMap[n] = uniIdReorder[uniId];
298  }
299  }
300 
301  /**
302  * @brief Construct PETSc matrix and vector handles.
303  *
304  * @todo Preallocation should be done at this point, since presently
305  * matrix allocation takes a significant amount of time.
306  *
307  * @param nGlobal Number of global degrees of freedom in the system (on
308  * this processor)
309  * @param nDir Number of Dirichlet degrees of freedom (on this
310  * processor).
311  */
312  void GlobalLinSysPETSc::SetUpMatVec(int nGlobal, int nDir)
313  {
314  // CREATE VECTORS
315  VecCreate (PETSC_COMM_WORLD, &m_x);
316  VecSetSizes (m_x, m_nLocal, PETSC_DECIDE);
317  VecSetFromOptions(m_x);
318  VecDuplicate (m_x, &m_b);
319 
320  // CREATE MATRICES
322  {
323  // Create ShellCtx context object which will store the matrix
324  // size and a pointer to the linear system. We do this so that
325  // we can call a member function to the matrix-vector and
326  // preconditioning multiplication in a subclass.
327  ShellCtx *ctx1 = new ShellCtx(), *ctx2 = new ShellCtx();
328  ctx1->nGlobal = ctx2->nGlobal = nGlobal;
329  ctx1->nDir = ctx2->nDir = nDir;
330  ctx1->linSys = ctx2->linSys = this;
331 
332  // Set up MatShell object.
333  MatCreateShell (PETSC_COMM_WORLD, m_nLocal, m_nLocal,
334  PETSC_DETERMINE, PETSC_DETERMINE,
335  (void *)ctx1, &m_matrix);
336  MatShellSetOperation(m_matrix, MATOP_MULT,
337  (void(*)(void))DoMatrixMultiply);
338  MatShellSetOperation(m_matrix, MATOP_DESTROY,
339  (void(*)(void))DoDestroyMatCtx);
340 
341  // Create a PCShell to go alongside the MatShell.
342  PCCreate (PETSC_COMM_WORLD, &m_pc);
343 #if PETSC_VERSION_GE(3,5,0)
344  PCSetOperators (m_pc, m_matrix, m_matrix);
345 #else
346  PCSetOperators (m_pc, m_matrix, m_matrix, SAME_NONZERO_PATTERN);
347 #endif
348  PCSetType (m_pc, PCSHELL);
349  PCShellSetApply (m_pc, DoPreconditioner);
350  PCShellSetDestroy(m_pc, DoDestroyPCCtx);
351  PCShellSetContext(m_pc, ctx2);
352  }
353  else
354  {
355  // Otherwise we create a PETSc matrix and use MatSetFromOptions
356  // so that we can set various options on the command line.
357  MatCreate (PETSC_COMM_WORLD, &m_matrix);
358  MatSetType (m_matrix, MATAIJ);
359  MatSetSizes (m_matrix, m_nLocal, m_nLocal,
360  PETSC_DETERMINE, PETSC_DETERMINE);
361  MatSetFromOptions(m_matrix);
362  MatSetUp (m_matrix);
363  }
364  }
365 
366  /**
367  * @brief Perform either matrix multiplication or preconditioning using
368  * Nektar++ routines.
369  *
370  * This static function uses Nektar++ routines to calculate the
371  * matrix-vector product of @p M with @p in, storing the output in @p
372  * out.
373  *
374  * @todo There's a lot of scatters and copies that might possibly be
375  * eliminated to make this more efficient.
376  *
377  * @param in Input vector.
378  * @param out Output vector.
379  * @param ctx ShellCtx object that points to our instance of
380  * GlobalLinSysPETSc.
381  * @param precon If true, we apply a preconditioner, if false, we
382  * perform a matrix multiplication.
383  */
385  Vec &in, Vec &out, ShellCtx *ctx, bool precon)
386  {
387  const int nGlobal = ctx->nGlobal;
388  const int nDir = ctx->nDir;
389  const int nHomDofs = nGlobal - nDir;
390  GlobalLinSysPETSc *linSys = ctx->linSys;
391 
392  // Scatter from PETSc ordering to our local ordering. It's actually
393  // unclear whether this step might also do some communication in
394  // parallel, which is probably not ideal.
395  VecScatterBegin(linSys->m_ctx, in, linSys->m_locVec,
396  INSERT_VALUES, SCATTER_FORWARD);
397  VecScatterEnd (linSys->m_ctx, in, linSys->m_locVec,
398  INSERT_VALUES, SCATTER_FORWARD);
399 
400  // Temporary storage to pass to Nektar++
401  Array<OneD, NekDouble> tmpIn(nHomDofs), tmpOut(nHomDofs);
402 
403  // Get values from input vector and copy to tmpIn.
404  PetscScalar *tmpLocIn;
405  VecGetArray (linSys->m_locVec, &tmpLocIn);
406  Vmath::Vcopy (nHomDofs, tmpLocIn, 1, &tmpIn[0], 1);
407  VecRestoreArray(linSys->m_locVec, &tmpLocIn);
408 
409  // Do matrix multiply in Nektar++, store in tmpOut.
410  if (precon)
411  {
412  linSys->m_precon->DoPreconditioner(tmpIn, tmpOut);
413  }
414  else
415  {
416  linSys->v_DoMatrixMultiply(tmpIn, tmpOut);
417  }
418 
419  // Scatter back to PETSc ordering and put in out.
420  VecSetValues(out, nHomDofs, &linSys->m_reorderedMap[0],
421  &tmpOut[0], INSERT_VALUES);
422  VecAssemblyBegin(out);
423  VecAssemblyEnd (out);
424  }
425 
426  /**
427  * @brief Perform matrix multiplication using Nektar++ routines.
428  *
429  * This static function uses Nektar++ routines to calculate the
430  * matrix-vector product of @p M with @p in, storing the output in @p
431  * out.
432  *
433  * @param M Original MatShell matrix, which stores the ShellCtx
434  * object.
435  * @param in Input vector.
436  * @param out Output vector.
437  */
439  Mat M, Vec in, Vec out)
440  {
441  // Grab our shell context from M.
442  void *ptr;
443  MatShellGetContext(M, &ptr);
444  ShellCtx *ctx = (ShellCtx *)ptr;
445 
446  DoNekppOperation(in, out, ctx, false);
447 
448  // Must return 0, otherwise PETSc complains.
449  return 0;
450  }
451 
452  /**
453  * @brief Apply preconditioning using Nektar++ routines.
454  *
455  * This static function uses Nektar++ routines to apply the
456  * preconditioner stored in GlobalLinSysPETSc::m_precon from the context
457  * of @p pc to the vector @p in, storing the output in @p out.
458  *
459  * @param pc Preconditioner object that stores the ShellCtx.
460  * @param in Input vector.
461  * @param out Output vector.
462  */
464  PC pc, Vec in, Vec out)
465  {
466  // Grab our PCShell context from pc.
467  void *ptr;
468  PCShellGetContext(pc, &ptr);
469  ShellCtx *ctx = (ShellCtx *)ptr;
470 
471  DoNekppOperation(in, out, ctx, true);
472 
473  // Must return 0, otherwise PETSc complains.
474  return 0;
475  }
476 
477  /**
478  * @brief Destroy matrix shell context object.
479  *
480  * Note the matrix shell and preconditioner share a common context so
481  * this might have already been deallocated below, in which case we do
482  * nothing.
483  *
484  * @param M Matrix shell object
485  */
487  {
488  void *ptr;
489  MatShellGetContext(M, &ptr);
490  ShellCtx *ctx = (ShellCtx *)ptr;
491  delete ctx;
492  return 0;
493  }
494 
495  /**
496  * @brief Destroy preconditioner context object.
497  *
498  * Note the matrix shell and preconditioner share a common context so
499  * this might have already been deallocated above, in which case we do
500  * nothing.
501  *
502  * @param pc Preconditioner object
503  */
504  PetscErrorCode GlobalLinSysPETSc::DoDestroyPCCtx(PC pc)
505  {
506  void *ptr;
507  PCShellGetContext(pc, &ptr);
508  ShellCtx *ctx = (ShellCtx *)ptr;
509  delete ctx;
510  return 0;
511  }
512 
513  /**
514  * @brief Set up KSP solver object.
515  *
516  * This is reasonably generic setup -- most solver types can be changed
517  * using the .petscrc file.
518  *
519  * @param tolerance Residual tolerance to converge to.
520  */
522  {
523  KSPCreate(PETSC_COMM_WORLD, &m_ksp);
524  KSPSetTolerances(
525  m_ksp, tolerance, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT);
526  KSPSetFromOptions(m_ksp);
527 #if PETSC_VERSION_GE(3,5,0)
528  KSPSetOperators(m_ksp, m_matrix, m_matrix);
529 #else
530  KSPSetOperators(m_ksp, m_matrix, m_matrix, SAME_NONZERO_PATTERN);
531 #endif
532 
534  {
535  KSPSetPC(m_ksp, m_pc);
536  }
537  }
538  }
539 }
VecScatter m_ctx
PETSc scatter context that takes us between Nektar++ global ordering and PETSc vector ordering...
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
GlobalLinSysPETSc * linSys
Pointer to the original calling object.
boost::shared_ptr< AssemblyMap > AssemblyMapSharedPtr
Definition: AssemblyMap.h:53
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.
STL namespace.
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.
Definition: Comm.h:55
boost::shared_ptr< CommMpi > CommMpiSharedPtr
Pointer to a Communicator object.
Definition: CommMpi.h:55
int m_nLocal
Number of unique degrees of freedom on this process.
KSP m_ksp
KSP object that represents solver system.
double NekDouble
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 global linear system.
Definition: CommMpi.h:61
A PETSc global linear system.
void SetUpScatter()
Set up PETSc local (equivalent to Nektar++ global) and global (equivalent to universal) scatter maps...
A global linear system.
Definition: GlobalLinSys.h:74
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)
Definition: Vmath.cpp:737
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)
Definition: Vmath.cpp:1061
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.
Definition: GlobalLinSys.h:129