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 : template <typename T>
48 : class NumericVector;
49 : void libmesh_petsc_snes_finish_residual (ResidualContext & rc,
50 : Vec x,
51 : Vec r,
52 : NumericVector<Number> & R);
53 :
54 : // Allow users access to these functions in case they want to reuse them. Users shouldn't
55 : // need access to these most of the time as they are used internally by this object.
56 : extern "C"
57 : {
58 : PetscErrorCode libmesh_petsc_recalculate_monitor(SNES snes, PetscInt it, PetscReal norm, void* mctx);
59 : PetscErrorCode libmesh_petsc_snes_monitor (SNES, PetscInt its, PetscReal fnorm, void *);
60 : PetscErrorCode libmesh_petsc_snes_residual (SNES, Vec x, Vec r, void * ctx);
61 : PetscErrorCode libmesh_petsc_snes_fd_residual (SNES, Vec x, Vec r, void * ctx);
62 : PetscErrorCode libmesh_petsc_snes_mffd_residual (SNES snes, Vec x, Vec r, void * ctx);
63 : PetscErrorCode libmesh_petsc_snes_mffd_interface (void * ctx, Vec x, Vec r);
64 : PetscErrorCode libmesh_petsc_snes_jacobian (SNES, Vec x, Mat jac, Mat pc, void * ctx);
65 : PetscErrorCode libmesh_petsc_snes_precheck(SNESLineSearch, Vec X, Vec Y, PetscBool * changed, void * context);
66 : PetscErrorCode libmesh_petsc_snes_postcheck(SNESLineSearch, Vec x, Vec y, Vec w, PetscBool * changed_y, PetscBool * changed_w, void * context);
67 : PetscErrorCode libmesh_petsc_linesearch_shellfunc(SNESLineSearch linesearch, void * ctx);
68 : }
69 :
70 : /**
71 : * This class provides an interface to PETSc
72 : * iterative solvers that is compatible with the \p libMesh
73 : * \p NonlinearSolver<>
74 : *
75 : * \author Benjamin Kirk
76 : * \date 2002-2007
77 : */
78 : template <typename T>
79 330 : class PetscNonlinearSolver : public NonlinearSolver<T>
80 : {
81 : public:
82 : /**
83 : * The type of system
84 : */
85 : typedef NonlinearImplicitSystem sys_type;
86 :
87 : /**
88 : * Constructor. Initializes Petsc data structures
89 : */
90 : explicit
91 : PetscNonlinearSolver (sys_type & system);
92 :
93 : /**
94 : * Destructor.
95 : */
96 : ~PetscNonlinearSolver ();
97 :
98 : /**
99 : * Release all memory and clear data structures.
100 : */
101 : virtual void clear () override;
102 :
103 : /**
104 : * Initialize data structures if not done so already.
105 : * May assign a name to the solver in some implementations
106 : */
107 : virtual void init (const char * name = nullptr) override;
108 :
109 : /**
110 : * \returns The raw PETSc snes context pointer. This method calls init() so in that vein, we have
111 : * an optional name prefix argument that if provided will be given to the SNES context
112 : */
113 : SNES snes(const char * name = nullptr);
114 :
115 : /**
116 : * Call the Petsc solver, using the same matrix for the system and preconditioner matrices. Calls
117 : * the two-matrix overload below with \p pre_in for both.
118 : */
119 : virtual std::pair<unsigned int, Real>
120 : solve (SparseMatrix<T> & pre_in, // System Preconditioning Matrix
121 : NumericVector<T> &, // Solution vector
122 : NumericVector<T> &, // Residual vector
123 : const double, // Stopping tolerance
124 : const unsigned int) override; // N. Iterations
125 :
126 : /**
127 : * Call the Petsc solver, using \p jac_in as the actual SNES Jacobian operator (Amat) and
128 : * \p pre_in as the preconditioning matrix (Pmat).
129 : */
130 : virtual std::pair<unsigned int, Real>
131 : solve (SparseMatrix<T> & jac_in, // Jacobian operator matrix (Amat)
132 : SparseMatrix<T> & pre_in, // Preconditioning matrix (Pmat)
133 : NumericVector<T> &, // Solution vector
134 : NumericVector<T> &, // Residual vector
135 : const double, // Stopping tolerance
136 : const unsigned int) override; // N. Iterations
137 :
138 : /**
139 : * Prints a useful message about why the latest nonlinear solve
140 : * con(di)verged.
141 : */
142 : virtual void print_converged_reason() override;
143 :
144 : /**
145 : * \returns The currently-available (or most recently obtained, if
146 : * the SNES object has been destroyed) convergence reason.
147 : *
148 : * Refer to PETSc docs for the meaning of different
149 : * SNESConvergedReasons.
150 : */
151 : SNESConvergedReason get_converged_reason();
152 :
153 : /**
154 : * Get the total number of linear iterations done in the last solve
155 : */
156 : virtual int get_total_linear_iterations() override;
157 :
158 : /**
159 : * \returns The current nonlinear iteration number if called
160 : * *during* the solve(), for example by the user-specified residual
161 : * or Jacobian function.
162 : */
163 0 : virtual unsigned get_current_nonlinear_iteration_number() const override
164 0 : { return _current_nonlinear_iteration_number; }
165 :
166 : /**
167 : * Set if the residual should be zeroed out in the callback
168 : */
169 0 : void set_residual_zero_out(bool state) { _zero_out_residual = state; }
170 :
171 : /**
172 : * Set if the jacobian should be zeroed out in the callback
173 : */
174 0 : void set_jacobian_zero_out(bool state) { _zero_out_jacobian = state; }
175 :
176 : /**
177 : * Set to true to use the libMesh's default monitor, set to false to use your own
178 : */
179 0 : void use_default_monitor(bool state) { _default_monitor = state; }
180 :
181 : /**
182 : * Set to true to let PETSc reuse the base vector
183 : */
184 0 : void set_snesmf_reuse_base(bool state) { _snesmf_reuse_base = state; }
185 :
186 : /**
187 : * @return Whether we are reusing the nonlinear function evaluation as the base for doing
188 : * matrix-free approximation of the Jacobian action
189 : */
190 2398 : bool snes_mf_reuse_base() const { return _snesmf_reuse_base; }
191 :
192 : /**
193 : * Set whether evaluating the residual twice at the same solution gives exactly the same result.
194 : * This is used to decide whether to run the debug-mode check of matrix-free base vector reuse,
195 : * which requires an exact match. The check is already skipped for MPI-parallel and threaded runs;
196 : * applications whose assembly is non-deterministic for other reasons (e.g. due to device
197 : * parallelism) should set this to false. This only disables the check; the base vector is still
198 : * reused.
199 : */
200 0 : void set_reproducible_residual(bool state) { _reproducible_residual = state; }
201 :
202 : /**
203 : * @return Whether evaluating the residual twice at the same solution gives exactly the same
204 : * result
205 : */
206 2398 : bool reproducible_residual() const { return _reproducible_residual; }
207 :
208 : /**
209 : * Set whether we are computing the base vector for matrix-free finite-differencing
210 : */
211 0 : void set_computing_base_vector(bool computing_base_vector) { _computing_base_vector = computing_base_vector; }
212 :
213 : /**
214 : * @return whether we are computing the base vector for matrix-free finite-differencing
215 : */
216 0 : bool computing_base_vector() const { return _computing_base_vector; }
217 :
218 : /**
219 : * Abstract base class to be used to implement a custom line-search algorithm. When the system
220 : * eliminates constrained degrees of freedom, the line search's vectors are reduced; \p prolong()
221 : * forms full vectors from them
222 : */
223 : class ComputeLineSearchObject
224 : {
225 : public:
226 : virtual ~ComputeLineSearchObject () = default;
227 :
228 : virtual void linesearch (SNESLineSearch linesearch) = 0;
229 : };
230 :
231 : /**
232 : * A callable object that can be used to specify a custom line-search
233 : */
234 : std::unique_ptr<ComputeLineSearchObject> linesearch_object;
235 :
236 : /**
237 : * Setup the default monitor if required
238 : */
239 : void setup_default_monitor();
240 :
241 : /**
242 : * Getter for preconditioner reuse
243 : */
244 : virtual bool reuse_preconditioner() const override;
245 :
246 : /**
247 : * Getter for the maximum iterations flag for preconditioner reuse
248 : */
249 : virtual unsigned int reuse_preconditioner_max_linear_its() const override;
250 :
251 : /**
252 : * Immediately force a new preconditioner, even if reuse is set
253 : */
254 : virtual void force_new_preconditioner() override;
255 :
256 : /**
257 : * Builds \p _retained_is from the system's \p ConstraintEliminationDofMap, or clears it when
258 : * the system does not eliminate constraints
259 : */
260 : void setup_constraint_elimination();
261 :
262 : /**
263 : * Whether the solve eliminates constrained degrees of freedom
264 : */
265 45732 : bool eliminating_constraints() const { return static_cast<IS>(_retained_is) != nullptr; }
266 :
267 : /**
268 : * Makes the system's \p current_local_solution the full iterate that \p x represents: \p x
269 : * itself, with constraints enforced exactly if requested, or \f$ C x + h \f$ when eliminating
270 : * constraints. When eliminating constraints the system's \p solution is set as well
271 : */
272 : void localize_solution(Vec x);
273 :
274 : /**
275 : * Prolongs the reduced vector \p y into the full vector \p x, which must be parallel. With \p
276 : * homogeneous false this forms \f$ x = C y + h \f$; with \p homogeneous true, as for a step or
277 : * a search direction, \f$ x = C y \f$
278 : */
279 : void prolong(Vec y, NumericVector<T> & x, bool homogeneous) const;
280 :
281 : /**
282 : * Restricts the full vector \p x into the reduced vector \p y by taking its retained entries
283 : */
284 : void restrict_vector(const NumericVector<T> & x, Vec y) const;
285 :
286 :
287 : protected:
288 :
289 : /**
290 : * Nonlinear solver context
291 : */
292 : WrappedPetsc<SNES> _snes;
293 :
294 : /**
295 : * Store the reason for SNES convergence/divergence for use even after the _snes
296 : * has been cleared.
297 : *
298 : * \note \p print_converged_reason() will always \e try to get the
299 : * current reason with SNESGetConvergedReason(), but if the SNES
300 : * object has already been cleared, it will fall back on this stored
301 : * value. This value is therefore necessarily \e not cleared by the
302 : * \p clear() function.
303 : */
304 : SNESConvergedReason _reason;
305 :
306 : /**
307 : * Stores the total number of linear iterations from the last solve.
308 : */
309 : PetscInt _n_linear_iterations;
310 :
311 : /**
312 : * Stores the current nonlinear iteration number
313 : */
314 : unsigned _current_nonlinear_iteration_number;
315 :
316 : /**
317 : * true to zero out residual before going into application level call-back, otherwise false
318 : */
319 : bool _zero_out_residual;
320 :
321 : /**
322 : * true to zero out jacobian before going into application level call-back, otherwise false
323 : */
324 : bool _zero_out_jacobian;
325 :
326 : /**
327 : * true if we want the default monitor to be set, false for no monitor (i.e. user code can use their own)
328 : */
329 : bool _default_monitor;
330 :
331 : /**
332 : * True, If we want the base vector to be used for differencing even if the function provided to SNESSetFunction()
333 : * is not the same as that provided to MatMFFDSetFunction().
334 : * https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/SNES/MatSNESMFSetReuseBase.html
335 : */
336 : bool _snesmf_reuse_base;
337 :
338 : /**
339 : * Whether residual evaluation is reproducible. See set_reproducible_residual()
340 : */
341 : bool _reproducible_residual;
342 :
343 : void build_mat_null_space(NonlinearImplicitSystem::ComputeVectorSubspace * computeSubspaceObject,
344 : void (*)(std::vector<NumericVector<Number> *> &, sys_type &),
345 : MatNullSpace *);
346 :
347 : /**
348 : * Whether we are computing the base vector for matrix-free finite differencing
349 : */
350 : bool _computing_base_vector;
351 :
352 : /**
353 : * Whether we've triggered the preconditioner reuse
354 : */
355 : bool _setup_reuse;
356 :
357 : /**
358 : * The local retained degrees of freedom of the system's \p ConstraintEliminationDofMap, in the
359 : * full numbering, when the solve eliminates constrained degrees of freedom. Entry \p i is the
360 : * full index of the \p i-th local entry of a reduced vector. Null otherwise
361 : */
362 : WrappedPetsc<IS> _retained_is;
363 :
364 : #if defined(LIBMESH_ENABLE_AMR) && defined(LIBMESH_HAVE_METAPHYSICL)
365 : /**
366 : * Wrapper object for interacting with the "new" libMesh PETSc DM. The new libMesh PETSc DM
367 : * implementation is capable of geometric multigrid while the old implementation is not. The new
368 : * implementation can be activated from the command line with --use_petsc_dm
369 : */
370 : PetscDMWrapper _dm_wrapper;
371 : #endif
372 :
373 : private:
374 11160 : friend ResidualContext libmesh_petsc_snes_residual_helper (SNES snes, Vec x, void * ctx);
375 : friend void libmesh_petsc_snes_finish_residual (ResidualContext & rc,
376 : Vec x,
377 : Vec r,
378 : NumericVector<Number> & R);
379 : friend PetscErrorCode libmesh_petsc_snes_residual (SNES snes, Vec x, Vec r, void * ctx);
380 : friend PetscErrorCode libmesh_petsc_snes_fd_residual (SNES snes, Vec x, Vec r, void * ctx);
381 : friend PetscErrorCode libmesh_petsc_snes_mffd_residual (SNES snes, Vec x, Vec r, void * ctx);
382 : friend PetscErrorCode libmesh_petsc_snes_jacobian (SNES snes, Vec x, Mat jac, Mat pc, void * ctx);
383 : friend PetscErrorCode libmesh_petsc_snes_postcheck (SNESLineSearch,
384 : Vec x,
385 : Vec y,
386 : Vec w,
387 : PetscBool * changed_y,
388 : PetscBool * changed_w,
389 : void * context);
390 : friend PetscErrorCode libmesh_petsc_snes_precheck (SNESLineSearch,
391 : Vec X,
392 : Vec Y,
393 : PetscBool * changed,
394 : void * context);
395 : };
396 :
397 : } // namespace libMesh
398 :
399 :
400 : #endif // #ifdef LIBMESH_HAVE_PETSC
401 : #endif // LIBMESH_PETSC_NONLINEAR_SOLVER_H
|
