Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 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 :
19 :
20 : #ifndef LIBMESH_PETSC_NONLINEAR_SOLVER_H
21 : #define LIBMESH_PETSC_NONLINEAR_SOLVER_H
22 :
23 : #include "libmesh/libmesh_config.h"
24 :
25 : // Petsc include files.
26 : #ifdef LIBMESH_HAVE_PETSC
27 :
28 : // Local includes
29 : #include "libmesh/nonlinear_solver.h"
30 : #include "libmesh/petsc_macro.h"
31 : #include "libmesh/wrapped_petsc.h"
32 : #include "libmesh/petsc_dm_wrapper.h"
33 : #include "libmesh/petsc_mffd_matrix.h"
34 :
35 : // PETSc includes
36 : #ifdef I
37 : # define LIBMESH_SAW_I
38 : #endif
39 : #include <petscsnes.h>
40 : #ifndef LIBMESH_SAW_I
41 : # undef I // Avoid complex.h contamination
42 : #endif
43 :
44 : namespace libMesh
45 : {
46 : class ResidualContext;
47 :
48 : // Allow users access to these functions in case they want to reuse them. Users shouldn't
49 : // need access to these most of the time as they are used internally by this object.
50 : extern "C"
51 : {
52 : PetscErrorCode libmesh_petsc_recalculate_monitor(SNES snes, PetscInt it, PetscReal norm, void* mctx);
53 : PetscErrorCode libmesh_petsc_snes_monitor (SNES, PetscInt its, PetscReal fnorm, void *);
54 : PetscErrorCode libmesh_petsc_snes_residual (SNES, Vec x, Vec r, void * ctx);
55 : PetscErrorCode libmesh_petsc_snes_fd_residual (SNES, Vec x, Vec r, void * ctx);
56 : PetscErrorCode libmesh_petsc_snes_mffd_residual (SNES snes, Vec x, Vec r, void * ctx);
57 : PetscErrorCode libmesh_petsc_snes_mffd_interface (void * ctx, Vec x, Vec r);
58 : PetscErrorCode libmesh_petsc_snes_jacobian (SNES, Vec x, Mat jac, Mat pc, void * ctx);
59 : PetscErrorCode libmesh_petsc_snes_precheck(SNESLineSearch, Vec X, Vec Y, PetscBool * changed, void * context);
60 : PetscErrorCode libmesh_petsc_snes_postcheck(SNESLineSearch, Vec x, Vec y, Vec w, PetscBool * changed_y, PetscBool * changed_w, void * context);
61 : PetscErrorCode libmesh_petsc_linesearch_shellfunc(SNESLineSearch linesearch, void * ctx);
62 : }
63 :
64 : /**
65 : * This class provides an interface to PETSc
66 : * iterative solvers that is compatible with the \p libMesh
67 : * \p NonlinearSolver<>
68 : *
69 : * \author Benjamin Kirk
70 : * \date 2002-2007
71 : */
72 : template <typename T>
73 210 : class PetscNonlinearSolver : public NonlinearSolver<T>
74 : {
75 : public:
76 : /**
77 : * The type of system
78 : */
79 : typedef NonlinearImplicitSystem sys_type;
80 :
81 : /**
82 : * Constructor. Initializes Petsc data structures
83 : */
84 : explicit
85 : PetscNonlinearSolver (sys_type & system);
86 :
87 : /**
88 : * Destructor.
89 : */
90 : ~PetscNonlinearSolver ();
91 :
92 : /**
93 : * Release all memory and clear data structures.
94 : */
95 : virtual void clear () override;
96 :
97 : /**
98 : * Initialize data structures if not done so already.
99 : * May assign a name to the solver in some implementations
100 : */
101 : virtual void init (const char * name = nullptr) override;
102 :
103 : /**
104 : * \returns The raw PETSc snes context pointer. This method calls init() so in that vein, we have
105 : * an optional name prefix argument that if provided will be given to the SNES context
106 : */
107 : SNES snes(const char * name = nullptr);
108 :
109 : /**
110 : * Call the Petsc solver, using the same matrix for the system and preconditioner matrices. Calls
111 : * the two-matrix overload below with \p pre_in for both.
112 : */
113 : virtual std::pair<unsigned int, Real>
114 : solve (SparseMatrix<T> & pre_in, // System Preconditioning Matrix
115 : NumericVector<T> &, // Solution vector
116 : NumericVector<T> &, // Residual vector
117 : const double, // Stopping tolerance
118 : const unsigned int) override; // N. Iterations
119 :
120 : /**
121 : * Call the Petsc solver, using \p jac_in as the actual SNES Jacobian operator (Amat) and
122 : * \p pre_in as the preconditioning matrix (Pmat).
123 : */
124 : virtual std::pair<unsigned int, Real>
125 : solve (SparseMatrix<T> & jac_in, // Jacobian operator matrix (Amat)
126 : SparseMatrix<T> & pre_in, // Preconditioning matrix (Pmat)
127 : NumericVector<T> &, // Solution vector
128 : NumericVector<T> &, // Residual vector
129 : const double, // Stopping tolerance
130 : const unsigned int) override; // N. Iterations
131 :
132 : /**
133 : * Prints a useful message about why the latest nonlinear solve
134 : * con(di)verged.
135 : */
136 : virtual void print_converged_reason() override;
137 :
138 : /**
139 : * \returns The currently-available (or most recently obtained, if
140 : * the SNES object has been destroyed) convergence reason.
141 : *
142 : * Refer to PETSc docs for the meaning of different
143 : * SNESConvergedReasons.
144 : */
145 : SNESConvergedReason get_converged_reason();
146 :
147 : /**
148 : * Get the total number of linear iterations done in the last solve
149 : */
150 : virtual int get_total_linear_iterations() override;
151 :
152 : /**
153 : * \returns The current nonlinear iteration number if called
154 : * *during* the solve(), for example by the user-specified residual
155 : * or Jacobian function.
156 : */
157 0 : virtual unsigned get_current_nonlinear_iteration_number() const override
158 0 : { return _current_nonlinear_iteration_number; }
159 :
160 : /**
161 : * Set if the residual should be zeroed out in the callback
162 : */
163 0 : void set_residual_zero_out(bool state) { _zero_out_residual = state; }
164 :
165 : /**
166 : * Set if the jacobian should be zeroed out in the callback
167 : */
168 0 : void set_jacobian_zero_out(bool state) { _zero_out_jacobian = state; }
169 :
170 : /**
171 : * Set to true to use the libMesh's default monitor, set to false to use your own
172 : */
173 0 : void use_default_monitor(bool state) { _default_monitor = state; }
174 :
175 : /**
176 : * Set to true to let PETSc reuse the base vector
177 : */
178 0 : void set_snesmf_reuse_base(bool state) { _snesmf_reuse_base = state; }
179 :
180 : /**
181 : * @return Whether we are reusing the nonlinear function evaluation as the base for doing
182 : * matrix-free approximation of the Jacobian action
183 : */
184 2398 : bool snes_mf_reuse_base() const { return _snesmf_reuse_base; }
185 :
186 : /**
187 : * Set whether we are computing the base vector for matrix-free finite-differencing
188 : */
189 0 : void set_computing_base_vector(bool computing_base_vector) { _computing_base_vector = computing_base_vector; }
190 :
191 : /**
192 : * @return whether we are computing the base vector for matrix-free finite-differencing
193 : */
194 0 : bool computing_base_vector() const { return _computing_base_vector; }
195 :
196 : /**
197 : * Abstract base class to be used to implement a custom line-search algorithm
198 : */
199 : class ComputeLineSearchObject
200 : {
201 : public:
202 : virtual ~ComputeLineSearchObject () = default;
203 :
204 : virtual void linesearch (SNESLineSearch linesearch) = 0;
205 : };
206 :
207 : /**
208 : * A callable object that can be used to specify a custom line-search
209 : */
210 : std::unique_ptr<ComputeLineSearchObject> linesearch_object;
211 :
212 : /**
213 : * Setup the default monitor if required
214 : */
215 : void setup_default_monitor();
216 :
217 : /**
218 : * Getter for preconditioner reuse
219 : */
220 : virtual bool reuse_preconditioner() const override;
221 :
222 : /**
223 : * Getter for the maximum iterations flag for preconditioner reuse
224 : */
225 : virtual unsigned int reuse_preconditioner_max_linear_its() const override;
226 :
227 : /**
228 : * Immediately force a new preconditioner, even if reuse is set
229 : */
230 : virtual void force_new_preconditioner() override;
231 :
232 : protected:
233 :
234 : /**
235 : * Nonlinear solver context
236 : */
237 : WrappedPetsc<SNES> _snes;
238 :
239 : /**
240 : * Store the reason for SNES convergence/divergence for use even after the _snes
241 : * has been cleared.
242 : *
243 : * \note \p print_converged_reason() will always \e try to get the
244 : * current reason with SNESGetConvergedReason(), but if the SNES
245 : * object has already been cleared, it will fall back on this stored
246 : * value. This value is therefore necessarily \e not cleared by the
247 : * \p clear() function.
248 : */
249 : SNESConvergedReason _reason;
250 :
251 : /**
252 : * Stores the total number of linear iterations from the last solve.
253 : */
254 : PetscInt _n_linear_iterations;
255 :
256 : /**
257 : * Stores the current nonlinear iteration number
258 : */
259 : unsigned _current_nonlinear_iteration_number;
260 :
261 : /**
262 : * true to zero out residual before going into application level call-back, otherwise false
263 : */
264 : bool _zero_out_residual;
265 :
266 : /**
267 : * true to zero out jacobian before going into application level call-back, otherwise false
268 : */
269 : bool _zero_out_jacobian;
270 :
271 : /**
272 : * true if we want the default monitor to be set, false for no monitor (i.e. user code can use their own)
273 : */
274 : bool _default_monitor;
275 :
276 : /**
277 : * True, If we want the base vector to be used for differencing even if the function provided to SNESSetFunction()
278 : * is not the same as that provided to MatMFFDSetFunction().
279 : * https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/SNES/MatSNESMFSetReuseBase.html
280 : */
281 : bool _snesmf_reuse_base;
282 :
283 : void build_mat_null_space(NonlinearImplicitSystem::ComputeVectorSubspace * computeSubspaceObject,
284 : void (*)(std::vector<NumericVector<Number> *> &, sys_type &),
285 : MatNullSpace *);
286 :
287 : /**
288 : * Whether we are computing the base vector for matrix-free finite differencing
289 : */
290 : bool _computing_base_vector;
291 :
292 : /**
293 : * Whether we've triggered the preconditioner reuse
294 : */
295 : bool _setup_reuse;
296 :
297 : #if defined(LIBMESH_ENABLE_AMR) && defined(LIBMESH_HAVE_METAPHYSICL)
298 : /**
299 : * Wrapper object for interacting with the "new" libMesh PETSc DM. The new libMesh PETSc DM
300 : * implementation is capable of geometric multigrid while the old implementation is not. The new
301 : * implementation can be activated from the command line with --use_petsc_dm
302 : */
303 : PetscDMWrapper _dm_wrapper;
304 : #endif
305 :
306 : private:
307 : friend ResidualContext libmesh_petsc_snes_residual_helper (SNES snes, Vec x, void * ctx);
308 : friend PetscErrorCode libmesh_petsc_snes_residual (SNES snes, Vec x, Vec r, void * ctx);
309 : friend PetscErrorCode libmesh_petsc_snes_fd_residual (SNES snes, Vec x, Vec r, void * ctx);
310 : friend PetscErrorCode libmesh_petsc_snes_mffd_residual (SNES snes, Vec x, Vec r, void * ctx);
311 : friend PetscErrorCode libmesh_petsc_snes_jacobian (SNES snes, Vec x, Mat jac, Mat pc, void * ctx);
312 : friend PetscErrorCode libmesh_petsc_snes_precheck (SNESLineSearch,
313 : Vec X,
314 : Vec Y,
315 : PetscBool * changed,
316 : void * context);
317 : };
318 :
319 : } // namespace libMesh
320 :
321 :
322 : #endif // #ifdef LIBMESH_HAVE_PETSC
323 : #endif // LIBMESH_PETSC_NONLINEAR_SOLVER_H
|
