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 240 : 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 evaluating the residual twice at the same solution gives exactly the same result.
188 : * This is used to decide whether to run the debug-mode check of matrix-free base vector reuse,
189 : * which requires an exact match. The check is already skipped for MPI-parallel and threaded runs;
190 : * applications whose assembly is non-deterministic for other reasons (e.g. due to device
191 : * parallelism) should set this to false. This only disables the check; the base vector is still
192 : * reused.
193 : */
194 0 : void set_reproducible_residual(bool state) { _reproducible_residual = state; }
195 :
196 : /**
197 : * @return Whether evaluating the residual twice at the same solution gives exactly the same
198 : * result
199 : */
200 2398 : bool reproducible_residual() const { return _reproducible_residual; }
201 :
202 : /**
203 : * Set whether we are computing the base vector for matrix-free finite-differencing
204 : */
205 0 : void set_computing_base_vector(bool computing_base_vector) { _computing_base_vector = computing_base_vector; }
206 :
207 : /**
208 : * @return whether we are computing the base vector for matrix-free finite-differencing
209 : */
210 0 : bool computing_base_vector() const { return _computing_base_vector; }
211 :
212 : /**
213 : * Abstract base class to be used to implement a custom line-search algorithm
214 : */
215 : class ComputeLineSearchObject
216 : {
217 : public:
218 : virtual ~ComputeLineSearchObject () = default;
219 :
220 : virtual void linesearch (SNESLineSearch linesearch) = 0;
221 : };
222 :
223 : /**
224 : * A callable object that can be used to specify a custom line-search
225 : */
226 : std::unique_ptr<ComputeLineSearchObject> linesearch_object;
227 :
228 : /**
229 : * Setup the default monitor if required
230 : */
231 : void setup_default_monitor();
232 :
233 : /**
234 : * Getter for preconditioner reuse
235 : */
236 : virtual bool reuse_preconditioner() const override;
237 :
238 : /**
239 : * Getter for the maximum iterations flag for preconditioner reuse
240 : */
241 : virtual unsigned int reuse_preconditioner_max_linear_its() const override;
242 :
243 : /**
244 : * Immediately force a new preconditioner, even if reuse is set
245 : */
246 : virtual void force_new_preconditioner() override;
247 :
248 : protected:
249 :
250 : /**
251 : * Nonlinear solver context
252 : */
253 : WrappedPetsc<SNES> _snes;
254 :
255 : /**
256 : * Store the reason for SNES convergence/divergence for use even after the _snes
257 : * has been cleared.
258 : *
259 : * \note \p print_converged_reason() will always \e try to get the
260 : * current reason with SNESGetConvergedReason(), but if the SNES
261 : * object has already been cleared, it will fall back on this stored
262 : * value. This value is therefore necessarily \e not cleared by the
263 : * \p clear() function.
264 : */
265 : SNESConvergedReason _reason;
266 :
267 : /**
268 : * Stores the total number of linear iterations from the last solve.
269 : */
270 : PetscInt _n_linear_iterations;
271 :
272 : /**
273 : * Stores the current nonlinear iteration number
274 : */
275 : unsigned _current_nonlinear_iteration_number;
276 :
277 : /**
278 : * true to zero out residual before going into application level call-back, otherwise false
279 : */
280 : bool _zero_out_residual;
281 :
282 : /**
283 : * true to zero out jacobian before going into application level call-back, otherwise false
284 : */
285 : bool _zero_out_jacobian;
286 :
287 : /**
288 : * true if we want the default monitor to be set, false for no monitor (i.e. user code can use their own)
289 : */
290 : bool _default_monitor;
291 :
292 : /**
293 : * True, If we want the base vector to be used for differencing even if the function provided to SNESSetFunction()
294 : * is not the same as that provided to MatMFFDSetFunction().
295 : * https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/SNES/MatSNESMFSetReuseBase.html
296 : */
297 : bool _snesmf_reuse_base;
298 :
299 : /**
300 : * Whether residual evaluation is reproducible. See set_reproducible_residual()
301 : */
302 : bool _reproducible_residual;
303 :
304 : void build_mat_null_space(NonlinearImplicitSystem::ComputeVectorSubspace * computeSubspaceObject,
305 : void (*)(std::vector<NumericVector<Number> *> &, sys_type &),
306 : MatNullSpace *);
307 :
308 : /**
309 : * Whether we are computing the base vector for matrix-free finite differencing
310 : */
311 : bool _computing_base_vector;
312 :
313 : /**
314 : * Whether we've triggered the preconditioner reuse
315 : */
316 : bool _setup_reuse;
317 :
318 : #if defined(LIBMESH_ENABLE_AMR) && defined(LIBMESH_HAVE_METAPHYSICL)
319 : /**
320 : * Wrapper object for interacting with the "new" libMesh PETSc DM. The new libMesh PETSc DM
321 : * implementation is capable of geometric multigrid while the old implementation is not. The new
322 : * implementation can be activated from the command line with --use_petsc_dm
323 : */
324 : PetscDMWrapper _dm_wrapper;
325 : #endif
326 :
327 : private:
328 : friend ResidualContext libmesh_petsc_snes_residual_helper (SNES snes, Vec x, void * ctx);
329 : friend PetscErrorCode libmesh_petsc_snes_residual (SNES snes, Vec x, Vec r, void * ctx);
330 : friend PetscErrorCode libmesh_petsc_snes_fd_residual (SNES snes, Vec x, Vec r, void * ctx);
331 : friend PetscErrorCode libmesh_petsc_snes_mffd_residual (SNES snes, Vec x, Vec r, void * ctx);
332 : friend PetscErrorCode libmesh_petsc_snes_jacobian (SNES snes, Vec x, Mat jac, Mat pc, void * ctx);
333 : friend PetscErrorCode libmesh_petsc_snes_precheck (SNESLineSearch,
334 : Vec X,
335 : Vec Y,
336 : PetscBool * changed,
337 : void * context);
338 : };
339 :
340 : } // namespace libMesh
341 :
342 :
343 : #endif // #ifdef LIBMESH_HAVE_PETSC
344 : #endif // LIBMESH_PETSC_NONLINEAR_SOLVER_H
|
