petsc_preconditioner.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 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
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 
34 EXTERN_C_FOR_PETSC_BEGIN
35 #include "petscksp.h"
36 EXTERN_C_FOR_PETSC_END
37 
38 #endif
39 
40 namespace libMesh
41 {
42 
43 template <typename T>
45 {
46  PetscVector<T> & x_pvec = libmesh_cast_ref<PetscVector<T>&>(const_cast<NumericVector<T>&>(x));
47  PetscVector<T> & y_pvec = libmesh_cast_ref<PetscVector<T>&>(const_cast<NumericVector<T>&>(y));
48 
49  Vec x_vec = x_pvec.vec();
50  Vec y_vec = y_pvec.vec();
51 
52  int ierr = PCApply(_pc,x_vec,y_vec);
53  LIBMESH_CHKERRABORT(ierr);
54 }
55 
56 
57 
58 
59 template <typename T>
61 {
62  if(!this->_matrix)
63  {
64  libMesh::err << "ERROR: No matrix set for PetscPreconditioner, but init() called" << std::endl;
65  libmesh_error();
66  }
67 
68  // Clear the preconditioner in case it has been created in the past
69  if (!this->_is_initialized)
70  {
71  // Should probably use PCReset(), but it's not working at the moment so we'll destroy instead
72  if (_pc)
73  {
74  int ierr = LibMeshPCDestroy(&_pc);
75  LIBMESH_CHKERRABORT(ierr);
76  }
77 
78  int ierr = PCCreate(this->comm().get(),&_pc);
79  LIBMESH_CHKERRABORT(ierr);
80 
81  PetscMatrix<T> * pmatrix = libmesh_cast_ptr<PetscMatrix<T>*, SparseMatrix<T> >(this->_matrix);
82 
83  _mat = pmatrix->mat();
84  }
85 
86  int ierr = PCSetOperators(_pc,_mat,_mat,SAME_NONZERO_PATTERN);
87  LIBMESH_CHKERRABORT(ierr);
88 
89  // Set the PCType. Note: this used to be done *before* the call to
90  // PCSetOperators(), and only when !_is_initialized, but
91  // 1.) Some preconditioners (those employing sub-preconditioners,
92  // for example) have to call PCSetUp(), and can only do this after
93  // the operators have been set.
94  // 2.) It should be safe to call set_petsc_preconditioner_type()
95  // multiple times.
96  set_petsc_preconditioner_type(this->_preconditioner_type, _pc);
97 
98  this->_is_initialized = true;
99 }
100 
101 
102 
103 template <typename T>
105 {
106  if (_pc)
107  {
108  int ierr = LibMeshPCDestroy(&_pc);
109  LIBMESH_CHKERRABORT(ierr);
110  }
111 }
112 
113 
114 
115 
116 template <typename T>
118 {
119  int ierr = 0;
120 
121  // get the communicator from the PETSc object
123  PetscObjectGetComm((PetscObject)pc, &comm);
125 
126  switch (preconditioner_type)
127  {
128  case IDENTITY_PRECOND:
129  ierr = PCSetType (pc, (char*) PCNONE); CHKERRABORT(comm,ierr); break;
130 
131  case CHOLESKY_PRECOND:
132  ierr = PCSetType (pc, (char*) PCCHOLESKY); CHKERRABORT(comm,ierr); break;
133 
134  case ICC_PRECOND:
135  ierr = PCSetType (pc, (char*) PCICC); CHKERRABORT(comm,ierr); break;
136 
137  case ILU_PRECOND:
138  {
139  // In serial, just set the ILU preconditioner type
140  if (communicator.size())
141  {
142  ierr = PCSetType (pc, (char*) PCILU);
143  CHKERRABORT(comm,ierr);
144  }
145  else
146  {
147  // But PETSc has no truly parallel ILU, instead you have to set
148  // an actual parallel preconditioner (e.g. block Jacobi) and then
149  // assign ILU sub-preconditioners.
150  ierr = PCSetType (pc, (char*) PCBJACOBI);
151  CHKERRABORT(comm,ierr);
152 
153  // Set ILU as the sub preconditioner type
154  set_petsc_subpreconditioner_type(PCILU, pc);
155  }
156  break;
157  }
158 
159  case LU_PRECOND:
160  {
161  // In serial, just set the LU preconditioner type
162  if (communicator.size())
163  {
164  ierr = PCSetType (pc, (char*) PCLU);
165  CHKERRABORT(comm,ierr);
166  }
167  else
168  {
169  // But PETSc has no truly parallel LU, instead you have to set
170  // an actual parallel preconditioner (e.g. block Jacobi) and then
171  // assign LU sub-preconditioners.
172  ierr = PCSetType (pc, (char*) PCBJACOBI);
173  CHKERRABORT(comm,ierr);
174 
175  // Set ILU as the sub preconditioner type
176  set_petsc_subpreconditioner_type(PCLU, pc);
177  }
178  break;
179  }
180 
181  case ASM_PRECOND:
182  {
183  // In parallel, I think ASM uses ILU by default as the sub-preconditioner...
184  // I tried setting a different sub-preconditioner here, but apparently the matrix
185  // is not in the correct state (at this point) to call PCSetUp().
186  ierr = PCSetType (pc, (char*) PCASM);
187  CHKERRABORT(comm,ierr);
188  break;
189  }
190 
191  case JACOBI_PRECOND:
192  ierr = PCSetType (pc, (char*) PCJACOBI); CHKERRABORT(comm,ierr); break;
193 
195  ierr = PCSetType (pc, (char*) PCBJACOBI); CHKERRABORT(comm,ierr); break;
196 
197  case SOR_PRECOND:
198  ierr = PCSetType (pc, (char*) PCSOR); CHKERRABORT(comm,ierr); break;
199 
200  case EISENSTAT_PRECOND:
201  ierr = PCSetType (pc, (char*) PCEISENSTAT); CHKERRABORT(comm,ierr); break;
202 
203  case AMG_PRECOND:
204  ierr = PCSetType (pc, (char*) PCHYPRE); CHKERRABORT(comm,ierr); break;
205 
206 #if !(PETSC_VERSION_LESS_THAN(2,1,2))
207  // Only available for PETSC >= 2.1.2
208  case USER_PRECOND:
209  ierr = PCSetType (pc, (char*) PCMAT); CHKERRABORT(comm,ierr); break;
210 #endif
211 
212  case SHELL_PRECOND:
213  ierr = PCSetType (pc, (char*) PCSHELL); CHKERRABORT(comm,ierr); break;
214 
215  default:
216  libMesh::err << "ERROR: Unsupported PETSC Preconditioner: "
217  << preconditioner_type << std::endl
218  << "Continuing with PETSC defaults" << std::endl;
219  }
220 
221  // Set additional options if we are doing AMG and
222  // HYPRE is available
223 #ifdef LIBMESH_HAVE_PETSC_HYPRE
224  if (preconditioner_type == AMG_PRECOND)
225  {
226  ierr = PCHYPRESetType(pc, "boomeramg");
227  CHKERRABORT(comm,ierr);
228  }
229 #endif
230 
231  // Let the commandline override stuff
232  // FIXME: Unless we are doing AMG???
233  if (preconditioner_type != AMG_PRECOND)
234  {
235  ierr = PCSetFromOptions(pc);
236  CHKERRABORT(comm,ierr);
237  }
238 }
239 
240 
241 
242 template <typename T>
243 #if PETSC_VERSION_LESS_THAN(3,0,0)
245 #else
246  void PetscPreconditioner<T>::set_petsc_subpreconditioner_type(const PCType type, PC& pc)
247 #endif
248 {
249  // For catching PETSc error return codes
250  int ierr = 0;
251 
252  // get the communicator from the PETSc object
254  PetscObjectGetComm((PetscObject)pc, &comm);
256 
257  // All docs say must call KSPSetUp or PCSetUp before calling PCBJacobiGetSubKSP.
258  // You must call PCSetUp after the preconditioner operators have been set, otherwise you get the:
259  //
260  // "Object is in wrong state!"
261  // "Matrix must be set first."
262  //
263  // error messages...
264  ierr = PCSetUp(pc);
265  CHKERRABORT(comm,ierr);
266 
267  // To store array of local KSP contexts on this processor
268  KSP* subksps;
269 
270  // the number of blocks on this processor
272 
273  // The global number of the first block on this processor.
274  // This is not used, so we just pass PETSC_NULL instead.
275  // int first_local;
276 
277  // Fill array of local KSP contexts
278  ierr = PCBJacobiGetSubKSP(pc, &n_local, PETSC_NULL, &subksps);
279  CHKERRABORT(comm,ierr);
280 
281  // Loop over sub-ksp objects, set ILU preconditioner
282  for (PetscInt i=0; i<n_local; ++i)
283  {
284  // Get pointer to sub KSP object's PC
285  PC subpc;
286  ierr = KSPGetPC(subksps[i], &subpc);
287  CHKERRABORT(comm,ierr);
288 
289  // Set requested type on the sub PC
290  ierr = PCSetType(subpc, type);
291  CHKERRABORT(comm,ierr);
292  }
293 
294 }
295 
296 
297 
298 
299 //------------------------------------------------------------------
300 // Explicit instantiations
301 template class PetscPreconditioner<Number>;
302 
303 } // namespace libMesh
304 
305 #endif // #ifdef LIBMESH_HAVE_PETSC

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:06 UTC

Hosted By:
SourceForge.net Logo