Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 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_base.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>
34 0 : PetscPreconditioner<T>::PetscPreconditioner (const libMesh::Parallel::Communicator & comm_in) :
35 0 : Preconditioner<T>(comm_in)
36 0 : {}
37 :
38 :
39 :
40 : template <typename T>
41 0 : void PetscPreconditioner<T>::apply(const NumericVector<T> & x, NumericVector<T> & y)
42 : {
43 0 : PetscVector<T> & x_pvec = cast_ref<PetscVector<T> &>(const_cast<NumericVector<T> &>(x));
44 0 : PetscVector<T> & y_pvec = cast_ref<PetscVector<T> &>(const_cast<NumericVector<T> &>(y));
45 :
46 0 : Vec x_vec = x_pvec.vec();
47 0 : Vec y_vec = y_pvec.vec();
48 :
49 0 : LibmeshPetscCall(PCApply(_pc, x_vec, y_vec));
50 0 : }
51 :
52 :
53 :
54 :
55 : template <typename T>
56 0 : void PetscPreconditioner<T>::init ()
57 : {
58 0 : libmesh_error_msg_if(!this->_matrix, "ERROR: No matrix set for PetscPreconditioner, but init() called");
59 :
60 : // Clear the preconditioner in case it has been created in the past
61 0 : if (!this->_is_initialized)
62 : {
63 : // Should probably use PCReset(), but it's not working at the moment so we'll destroy instead
64 0 : if (_pc)
65 0 : _pc.destroy();
66 :
67 0 : LibmeshPetscCall(PCCreate(this->comm().get(), _pc.get()));
68 :
69 0 : auto pmatrix = cast_ptr<PetscMatrixBase<T> *>(this->_matrix);
70 0 : _mat = pmatrix->mat();
71 : }
72 :
73 0 : LibmeshPetscCall(PCSetOperators(_pc, _mat, _mat));
74 :
75 : // Set the PCType. Note: this used to be done *before* the call to
76 : // PCSetOperators(), and only when !_is_initialized, but
77 : // 1.) Some preconditioners (those employing sub-preconditioners,
78 : // for example) have to call PCSetUp(), and can only do this after
79 : // the operators have been set.
80 : // 2.) It should be safe to call set_petsc_preconditioner_type()
81 : // multiple times.
82 0 : set_petsc_preconditioner_type(this->_preconditioner_type, *_pc);
83 :
84 0 : this->_is_initialized = true;
85 0 : }
86 :
87 :
88 :
89 : template <typename T>
90 0 : void PetscPreconditioner<T>::clear()
91 : {
92 : // Calls custom deleter
93 0 : _pc.destroy();
94 0 : }
95 :
96 :
97 :
98 : template <typename T>
99 0 : PC PetscPreconditioner<T>::pc()
100 : {
101 0 : return _pc;
102 : }
103 :
104 :
105 :
106 : template <typename T>
107 44861 : void PetscPreconditioner<T>::set_petsc_preconditioner_type (const PreconditionerType & preconditioner_type, PC & pc)
108 : {
109 : // get the communicator from the PETSc object
110 : Parallel::communicator comm;
111 44861 : PetscErrorCode ierr = PetscObjectGetComm((PetscObject)pc, & comm);
112 44861 : if (ierr != LIBMESH_PETSC_SUCCESS)
113 0 : libmesh_error_msg("Error retrieving communicator");
114 :
115 : #define CasePCSetType(PreconditionerType, PCType) \
116 : case PreconditionerType: \
117 : LibmeshPetscCallA(comm, PCSetType (pc, const_cast<KSPType>(PCType))); \
118 : break;
119 :
120 44861 : switch (preconditioner_type)
121 : {
122 140 : CasePCSetType(IDENTITY_PRECOND, PCNONE)
123 0 : CasePCSetType(CHOLESKY_PRECOND, PCCHOLESKY)
124 0 : CasePCSetType(ICC_PRECOND, PCICC)
125 633 : CasePCSetType(ILU_PRECOND, PCILU)
126 0 : CasePCSetType(LU_PRECOND, PCLU)
127 0 : CasePCSetType(ASM_PRECOND, PCASM)
128 0 : CasePCSetType(JACOBI_PRECOND, PCJACOBI)
129 43878 : CasePCSetType(BLOCK_JACOBI_PRECOND, PCBJACOBI)
130 0 : CasePCSetType(SOR_PRECOND, PCSOR)
131 0 : CasePCSetType(EISENSTAT_PRECOND, PCEISENSTAT)
132 0 : CasePCSetType(AMG_PRECOND, PCHYPRE)
133 0 : CasePCSetType(SVD_PRECOND, PCSVD)
134 0 : CasePCSetType(USER_PRECOND, PCMAT)
135 210 : CasePCSetType(SHELL_PRECOND, PCSHELL)
136 :
137 0 : default:
138 0 : libMesh::err << "ERROR: Unsupported PETSC Preconditioner: "
139 0 : << Utility::enum_to_string(preconditioner_type) << std::endl
140 0 : << "Continuing with PETSC defaults" << std::endl;
141 : }
142 :
143 : // Set additional options if we are doing AMG and
144 : // HYPRE is available
145 : #ifdef LIBMESH_HAVE_PETSC_HYPRE
146 44861 : if (preconditioner_type == AMG_PRECOND)
147 0 : LibmeshPetscCallA(comm, PCHYPRESetType(pc, "boomeramg"));
148 : #endif
149 :
150 : // Let the commandline override stuff
151 44861 : LibmeshPetscCallA(comm, PCSetFromOptions(pc));
152 44861 : }
153 :
154 :
155 :
156 : //------------------------------------------------------------------
157 : // Explicit instantiations
158 : template class LIBMESH_EXPORT PetscPreconditioner<Number>;
159 :
160 : } // namespace libMesh
161 :
162 : #endif // #ifdef LIBMESH_HAVE_PETSC
|
