libMesh
petsc_preconditioner.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 #include "libmesh/libmesh_common.h"
19 
20 #ifdef LIBMESH_HAVE_PETSC
21 
22 // C++ includes
23 
24 // Local Includes
25 #include "libmesh/petsc_preconditioner.h"
26 #include "libmesh/petsc_macro.h"
27 #include "libmesh/petsc_matrix.h"
28 #include "libmesh/petsc_vector.h"
29 #include "libmesh/libmesh_common.h"
30 
31 // PCBJacobiGetSubKSP was defined in petscksp.h in PETSc 2.3.3, 3.1.0
32 #if PETSC_VERSION_LESS_THAN(3,1,0)
33 # include "petscksp.h"
34 #endif
35 
36 namespace libMesh
37 {
38 
39 template <typename T>
41 {
42  PetscVector<T> & x_pvec = cast_ref<PetscVector<T> &>(const_cast<NumericVector<T> &>(x));
43  PetscVector<T> & y_pvec = cast_ref<PetscVector<T> &>(const_cast<NumericVector<T> &>(y));
44 
45  Vec x_vec = x_pvec.vec();
46  Vec y_vec = y_pvec.vec();
47 
48  int ierr = PCApply(_pc,x_vec,y_vec);
49  LIBMESH_CHKERR(ierr);
50 }
51 
52 
53 
54 
55 template <typename T>
57 {
58  if (!this->_matrix)
59  libmesh_error_msg("ERROR: No matrix set for PetscPreconditioner, but init() called");
60 
61  // Clear the preconditioner in case it has been created in the past
62  if (!this->_is_initialized)
63  {
64  // Should probably use PCReset(), but it's not working at the moment so we'll destroy instead
65  if (_pc)
66  {
67  int ierr = LibMeshPCDestroy(&_pc);
68  LIBMESH_CHKERR(ierr);
69  }
70 
71  int ierr = PCCreate(this->comm().get(),&_pc);
72  LIBMESH_CHKERR(ierr);
73 
74  PetscMatrix<T> * pmatrix = cast_ptr<PetscMatrix<T> *, SparseMatrix<T>>(this->_matrix);
75 
76  _mat = pmatrix->mat();
77  }
78 
79 #if PETSC_RELEASE_LESS_THAN(3,5,0)
80  int ierr = PCSetOperators(_pc,_mat,_mat,SAME_NONZERO_PATTERN);
81 #else
82  int ierr = PCSetOperators(_pc,_mat,_mat);
83 #endif
84  LIBMESH_CHKERR(ierr);
85 
86  // Set the PCType. Note: this used to be done *before* the call to
87  // PCSetOperators(), and only when !_is_initialized, but
88  // 1.) Some preconditioners (those employing sub-preconditioners,
89  // for example) have to call PCSetUp(), and can only do this after
90  // the operators have been set.
91  // 2.) It should be safe to call set_petsc_preconditioner_type()
92  // multiple times.
93  set_petsc_preconditioner_type(this->_preconditioner_type, _pc);
94 
95  this->_is_initialized = true;
96 }
97 
98 
99 
100 template <typename T>
102 {
103  if (_pc)
104  {
105  int ierr = LibMeshPCDestroy(&_pc);
106  LIBMESH_CHKERR(ierr);
107  }
108 }
109 
110 
111 
112 
113 template <typename T>
115 {
116  int ierr = 0;
117 
118  // get the communicator from the PETSc object
120  PetscObjectGetComm((PetscObject)pc, & comm);
122 
123  switch (preconditioner_type)
124  {
125  case IDENTITY_PRECOND:
126  ierr = PCSetType (pc, const_cast<KSPType>(PCNONE));
127  CHKERRABORT(comm,ierr);
128  break;
129 
130  case CHOLESKY_PRECOND:
131  ierr = PCSetType (pc, const_cast<KSPType>(PCCHOLESKY));
132  CHKERRABORT(comm,ierr);
133  break;
134 
135  case ICC_PRECOND:
136  ierr = PCSetType (pc, const_cast<KSPType>(PCICC));
137  CHKERRABORT(comm,ierr);
138  break;
139 
140  case ILU_PRECOND:
141  {
142  // In serial, just set the ILU preconditioner type
143  if (communicator.size())
144  {
145  ierr = PCSetType (pc, const_cast<KSPType>(PCILU));
146  CHKERRABORT(comm,ierr);
147  }
148  else
149  {
150  // But PETSc has no truly parallel ILU, instead you have to set
151  // an actual parallel preconditioner (e.g. block Jacobi) and then
152  // assign ILU sub-preconditioners.
153  ierr = PCSetType (pc, const_cast<KSPType>(PCBJACOBI));
154  CHKERRABORT(comm,ierr);
155 
156  // Set ILU as the sub preconditioner type
157  set_petsc_subpreconditioner_type(PCILU, pc);
158  }
159  break;
160  }
161 
162  case LU_PRECOND:
163  {
164  // In serial, just set the LU preconditioner type
165  if (communicator.size())
166  {
167  ierr = PCSetType (pc, const_cast<KSPType>(PCLU));
168  CHKERRABORT(comm,ierr);
169  }
170  else
171  {
172  // But PETSc has no truly parallel LU, instead you have to set
173  // an actual parallel preconditioner (e.g. block Jacobi) and then
174  // assign LU sub-preconditioners.
175  ierr = PCSetType (pc, const_cast<KSPType>(PCBJACOBI));
176  CHKERRABORT(comm,ierr);
177 
178  // Set ILU as the sub preconditioner type
179  set_petsc_subpreconditioner_type(PCLU, pc);
180  }
181  break;
182  }
183 
184  case ASM_PRECOND:
185  {
186  // In parallel, I think ASM uses ILU by default as the sub-preconditioner...
187  // I tried setting a different sub-preconditioner here, but apparently the matrix
188  // is not in the correct state (at this point) to call PCSetUp().
189  ierr = PCSetType (pc, const_cast<KSPType>(PCASM));
190  CHKERRABORT(comm,ierr);
191  break;
192  }
193 
194  case JACOBI_PRECOND:
195  ierr = PCSetType (pc, const_cast<KSPType>(PCJACOBI));
196  CHKERRABORT(comm,ierr);
197  break;
198 
200  ierr = PCSetType (pc, const_cast<KSPType>(PCBJACOBI));
201  CHKERRABORT(comm,ierr);
202  break;
203 
204  case SOR_PRECOND:
205  ierr = PCSetType (pc, const_cast<KSPType>(PCSOR));
206  CHKERRABORT(comm,ierr);
207  break;
208 
209  case EISENSTAT_PRECOND:
210  ierr = PCSetType (pc, const_cast<KSPType>(PCEISENSTAT));
211  CHKERRABORT(comm,ierr);
212  break;
213 
214  case AMG_PRECOND:
215  ierr = PCSetType (pc, const_cast<KSPType>(PCHYPRE));
216  CHKERRABORT(comm,ierr);
217  break;
218 
219  case USER_PRECOND:
220  ierr = PCSetType (pc, const_cast<KSPType>(PCMAT));
221  CHKERRABORT(comm,ierr);
222  break;
223 
224  case SHELL_PRECOND:
225  ierr = PCSetType (pc, const_cast<KSPType>(PCSHELL));
226  CHKERRABORT(comm,ierr);
227  break;
228 
229  default:
230  libMesh::err << "ERROR: Unsupported PETSC Preconditioner: "
231  << preconditioner_type << std::endl
232  << "Continuing with PETSC defaults" << std::endl;
233  }
234 
235  // Set additional options if we are doing AMG and
236  // HYPRE is available
237 #ifdef LIBMESH_HAVE_PETSC_HYPRE
238  if (preconditioner_type == AMG_PRECOND)
239  {
240  ierr = PCHYPRESetType(pc, "boomeramg");
241  CHKERRABORT(comm,ierr);
242  }
243 #endif
244 
245  // Let the commandline override stuff
246  ierr = PCSetFromOptions(pc);
247  CHKERRABORT(comm,ierr);
248 }
249 
250 
251 
252 template <typename T>
253 #if PETSC_VERSION_LESS_THAN(3,0,0)
255 #else
256  void PetscPreconditioner<T>::set_petsc_subpreconditioner_type(const PCType type, PC & pc)
257 #endif
258 {
259  // For catching PETSc error return codes
260  int ierr = 0;
261 
262  // get the communicator from the PETSc object
264  PetscObjectGetComm((PetscObject)pc, & comm);
266 
267  // All docs say must call KSPSetUp or PCSetUp before calling PCBJacobiGetSubKSP.
268  // You must call PCSetUp after the preconditioner operators have been set, otherwise you get the:
269  //
270  // "Object is in wrong state!"
271  // "Matrix must be set first."
272  //
273  // error messages...
274  ierr = PCSetUp(pc);
275  CHKERRABORT(comm,ierr);
276 
277  // To store array of local KSP contexts on this processor
278  KSP * subksps;
279 
280  // the number of blocks on this processor
281  PetscInt n_local;
282 
283  // The global number of the first block on this processor.
284  // This is not used, so we just pass PETSC_NULL instead.
285  // int first_local;
286 
287  // Fill array of local KSP contexts
288  ierr = PCBJacobiGetSubKSP(pc, &n_local, PETSC_NULL, &subksps);
289  CHKERRABORT(comm,ierr);
290 
291  // Loop over sub-ksp objects, set ILU preconditioner
292  for (PetscInt i=0; i<n_local; ++i)
293  {
294  // Get pointer to sub KSP object's PC
295  PC subpc;
296  ierr = KSPGetPC(subksps[i], &subpc);
297  CHKERRABORT(comm,ierr);
298 
299  // Set requested type on the sub PC
300  ierr = PCSetType(subpc, type);
301  CHKERRABORT(comm,ierr);
302  }
303 
304 }
305 
306 
307 
308 
309 //------------------------------------------------------------------
310 // Explicit instantiations
311 template class PetscPreconditioner<Number>;
312 
313 } // namespace libMesh
314 
315 #endif // #ifdef LIBMESH_HAVE_PETSC
OStreamProxy err
Encapsulates the MPI_Comm object.
Definition: parallel.h:657
This class provides an interface to the suite of preconditioners available from PETSc.
unsigned int size() const
Definition: parallel.h:726
MPI_Comm communicator
Communicator object for talking with subsets of processors.
Definition: parallel.h:181
Numeric vector.
Definition: dof_map.h:66
The libMesh namespace provides an interface to certain functionality in the library.
PetscErrorCode Vec Mat Mat pc
Generic sparse matrix.
Definition: dof_map.h:65
PreconditionerType
Defines an enum for preconditioner types.
PetscErrorCode Vec x
bool _is_initialized
Flag that tells if init() has been called.
Definition: libmesh.C:255
static void set_petsc_subpreconditioner_type(PCType type, PC &pc)
Some PETSc preconditioners (ILU, LU) don&#39;t work in parallel.
static void set_petsc_preconditioner_type(const PreconditionerType &preconditioner_type, PC &pc)
Tells PETSc to use the user-specified preconditioner.
virtual void init() libmesh_override
Initialize data structures if not done so already.
PetscErrorCode ierr
virtual void apply(const NumericVector< T > &x, NumericVector< T > &y) libmesh_override
Computes the preconditioned vector y based on input vector x.
virtual void clear() libmesh_override
Release all memory and clear data structures.