libMesh
petsc_preconditioner.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2024 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 // Local Includes
23 #include "libmesh/petsc_preconditioner.h"
24 #include "libmesh/petsc_macro.h"
25 #include "libmesh/petsc_matrix.h"
26 #include "libmesh/petsc_vector.h"
27 #include "libmesh/libmesh_common.h"
28 #include "libmesh/enum_preconditioner_type.h"
29 
30 namespace libMesh
31 {
32 
33 template <typename T>
35  Preconditioner<T>(comm_in)
36 {}
37 
38 
39 
40 template <typename T>
42 {
43  PetscVector<T> & x_pvec = cast_ref<PetscVector<T> &>(const_cast<NumericVector<T> &>(x));
44  PetscVector<T> & y_pvec = cast_ref<PetscVector<T> &>(const_cast<NumericVector<T> &>(y));
45 
46  Vec x_vec = x_pvec.vec();
47  Vec y_vec = y_pvec.vec();
48 
49  PetscErrorCode ierr = PCApply(_pc, x_vec, y_vec);
50  LIBMESH_CHKERR(ierr);
51 }
52 
53 
54 
55 
56 template <typename T>
58 {
59  libmesh_error_msg_if(!this->_matrix, "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  _pc.destroy();
67 
68  PetscErrorCode ierr = PCCreate(this->comm().get(), _pc.get());
69  LIBMESH_CHKERR(ierr);
70 
71  auto pmatrix = cast_ptr<PetscMatrix<T> *>(this->_matrix);
72  _mat = pmatrix->mat();
73  }
74 
75  PetscErrorCode ierr = PCSetOperators(_pc, _mat, _mat);
76  LIBMESH_CHKERR(ierr);
77 
78  // Set the PCType. Note: this used to be done *before* the call to
79  // PCSetOperators(), and only when !_is_initialized, but
80  // 1.) Some preconditioners (those employing sub-preconditioners,
81  // for example) have to call PCSetUp(), and can only do this after
82  // the operators have been set.
83  // 2.) It should be safe to call set_petsc_preconditioner_type()
84  // multiple times.
85  set_petsc_preconditioner_type(this->_preconditioner_type, *_pc);
86 
87  this->_is_initialized = true;
88 }
89 
90 
91 
92 template <typename T>
94 {
95  // Calls custom deleter
96  _pc.destroy();
97 }
98 
99 
100 
101 template <typename T>
103 {
104  return _pc;
105 }
106 
107 
108 
109 template <typename T>
111 {
112  PetscErrorCode ierr = 0;
113 
114  // get the communicator from the PETSc object
115  Parallel::communicator comm;
116  PetscObjectGetComm((PetscObject)pc, & comm);
118 
119  switch (preconditioner_type)
120  {
121  case IDENTITY_PRECOND:
122  ierr = PCSetType (pc, const_cast<KSPType>(PCNONE));
123  CHKERRABORT(comm,ierr);
124  break;
125 
126  case CHOLESKY_PRECOND:
127  ierr = PCSetType (pc, const_cast<KSPType>(PCCHOLESKY));
128  CHKERRABORT(comm,ierr);
129  break;
130 
131  case ICC_PRECOND:
132  ierr = PCSetType (pc, const_cast<KSPType>(PCICC));
133  CHKERRABORT(comm,ierr);
134  break;
135 
136  case ILU_PRECOND:
137  {
138  // In serial, just set the ILU preconditioner type
139  if (communicator.size())
140  {
141  ierr = PCSetType (pc, const_cast<KSPType>(PCILU));
142  CHKERRABORT(comm,ierr);
143  }
144  else
145  {
146  // But PETSc has no truly parallel ILU, instead you have to set
147  // an actual parallel preconditioner (e.g. block Jacobi) and then
148  // assign ILU sub-preconditioners.
149  ierr = PCSetType (pc, const_cast<KSPType>(PCBJACOBI));
150  CHKERRABORT(comm,ierr);
151 
152  // Set ILU as the sub preconditioner type
153  set_petsc_subpreconditioner_type(PCILU, pc);
154  }
155  break;
156  }
157 
158  case LU_PRECOND:
159  {
160  // In serial, just set the LU preconditioner type
161  if (communicator.size())
162  {
163  ierr = PCSetType (pc, const_cast<KSPType>(PCLU));
164  CHKERRABORT(comm,ierr);
165  }
166  else
167  {
168  // But PETSc has no truly parallel LU, instead you have to set
169  // an actual parallel preconditioner (e.g. block Jacobi) and then
170  // assign LU sub-preconditioners.
171  ierr = PCSetType (pc, const_cast<KSPType>(PCBJACOBI));
172  CHKERRABORT(comm,ierr);
173 
174  // Set ILU as the sub preconditioner type
175  set_petsc_subpreconditioner_type(PCLU, pc);
176  }
177  break;
178  }
179 
180  case ASM_PRECOND:
181  {
182  // In parallel, I think ASM uses ILU by default as the sub-preconditioner...
183  // I tried setting a different sub-preconditioner here, but apparently the matrix
184  // is not in the correct state (at this point) to call PCSetUp().
185  ierr = PCSetType (pc, const_cast<KSPType>(PCASM));
186  CHKERRABORT(comm,ierr);
187  break;
188  }
189 
190  case JACOBI_PRECOND:
191  ierr = PCSetType (pc, const_cast<KSPType>(PCJACOBI));
192  CHKERRABORT(comm,ierr);
193  break;
194 
196  ierr = PCSetType (pc, const_cast<KSPType>(PCBJACOBI));
197  CHKERRABORT(comm,ierr);
198  break;
199 
200  case SOR_PRECOND:
201  ierr = PCSetType (pc, const_cast<KSPType>(PCSOR));
202  CHKERRABORT(comm,ierr);
203  break;
204 
205  case EISENSTAT_PRECOND:
206  ierr = PCSetType (pc, const_cast<KSPType>(PCEISENSTAT));
207  CHKERRABORT(comm,ierr);
208  break;
209 
210  case AMG_PRECOND:
211  ierr = PCSetType (pc, const_cast<KSPType>(PCHYPRE));
212  CHKERRABORT(comm,ierr);
213  break;
214 
215  case SVD_PRECOND:
216  ierr = PCSetType (pc, const_cast<KSPType>(PCSVD));
217  CHKERRABORT(comm,ierr);
218  break;
219 
220  case USER_PRECOND:
221  ierr = PCSetType (pc, const_cast<KSPType>(PCMAT));
222  CHKERRABORT(comm,ierr);
223  break;
224 
225  case SHELL_PRECOND:
226  ierr = PCSetType (pc, const_cast<KSPType>(PCSHELL));
227  CHKERRABORT(comm,ierr);
228  break;
229 
230  default:
231  libMesh::err << "ERROR: Unsupported PETSC Preconditioner: "
232  << preconditioner_type << std::endl
233  << "Continuing with PETSC defaults" << std::endl;
234  }
235 
236  // Set additional options if we are doing AMG and
237  // HYPRE is available
238 #ifdef LIBMESH_HAVE_PETSC_HYPRE
239  if (preconditioner_type == AMG_PRECOND)
240  {
241  ierr = PCHYPRESetType(pc, "boomeramg");
242  CHKERRABORT(comm,ierr);
243  }
244 #endif
245 
246  // Let the commandline override stuff
247  ierr = PCSetFromOptions(pc);
248  CHKERRABORT(comm,ierr);
249 }
250 
251 
252 template <typename T>
254 {
255  // For catching PETSc error return codes
256  PetscErrorCode ierr = 0;
257 
258  // get the communicator from the PETSc object
259  Parallel::communicator comm;
260  PetscObjectGetComm((PetscObject)pc, & comm);
262 
263  // All docs say must call KSPSetUp or PCSetUp before calling PCBJacobiGetSubKSP.
264  // You must call PCSetUp after the preconditioner operators have been set, otherwise you get the:
265  //
266  // "Object is in wrong state!"
267  // "Matrix must be set first."
268  //
269  // error messages...
270  ierr = PCSetUp(pc);
271  CHKERRABORT(comm,ierr);
272 
273  // To store array of local KSP contexts on this processor
274  KSP * subksps;
275 
276  // the number of blocks on this processor
277  PetscInt n_local;
278 
279  // The global number of the first block on this processor.
280  // This is not used, so we just pass null instead.
281  // int first_local;
282 
283  // Fill array of local KSP contexts
284  ierr = PCBJacobiGetSubKSP(pc, &n_local, LIBMESH_PETSC_NULLPTR,
285  &subksps);
286  CHKERRABORT(comm,ierr);
287 
288  // Loop over sub-ksp objects, set ILU preconditioner
289  for (PetscInt i=0; i<n_local; ++i)
290  {
291  // Get pointer to sub KSP object's PC
292  PC subpc;
293  ierr = KSPGetPC(subksps[i], &subpc);
294  CHKERRABORT(comm,ierr);
295 
296  // Set requested type on the sub PC
297  ierr = PCSetType(subpc, type);
298  CHKERRABORT(comm,ierr);
299  }
300 }
301 
302 
303 
304 
305 //------------------------------------------------------------------
306 // Explicit instantiations
307 template class LIBMESH_EXPORT PetscPreconditioner<Number>;
308 
309 } // namespace libMesh
310 
311 #endif // #ifdef LIBMESH_HAVE_PETSC
OStreamProxy err
This class provides a nice interface to PETSc&#39;s Vec object.
Definition: petsc_vector.h:73
This class provides an interface to the suite of preconditioners available from PETSc.
virtual void init() override
Initialize data structures if not done so already.
Provides a uniform interface to vector storage schemes for different linear algebra libraries...
Definition: vector_fe_ex5.C:43
static void set_petsc_subpreconditioner_type(const PCType type, PC &pc)
Some PETSc preconditioners (ILU, LU) don&#39;t work in parallel.
The libMesh namespace provides an interface to certain functionality in the library.
virtual void apply(const NumericVector< T > &x, NumericVector< T > &y) override
Computes the preconditioned vector y based on input vector x.
This class provides a uniform interface for preconditioners.
bool _is_initialized
Flag that tells if init() has been called.
Definition: libmesh.C:247
PetscPreconditioner(const libMesh::Parallel::Communicator &comm_in)
Constructor.
virtual void clear() override
Release all memory and clear data structures.
static void set_petsc_preconditioner_type(const PreconditionerType &preconditioner_type, PC &pc)
Tells PETSc to use the user-specified preconditioner.
PreconditionerType
Defines an enum for preconditioner types.
DIE A HORRIBLE DEATH HERE typedef MPI_Comm communicator