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 :
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 : #ifdef LIBMESH_ENABLE_DEPRECATED
64 : PetscErrorCode __libmesh_petsc_snes_monitor (SNES, PetscInt its, PetscReal fnorm, void *);
65 : PetscErrorCode __libmesh_petsc_snes_residual (SNES, Vec x, Vec r, void * ctx);
66 : PetscErrorCode __libmesh_petsc_snes_fd_residual (SNES, Vec x, Vec r, void * ctx);
67 : PetscErrorCode __libmesh_petsc_snes_mffd_interface (void * ctx, Vec x, Vec r);
68 : PetscErrorCode __libmesh_petsc_snes_jacobian (SNES, Vec x, Mat jac, Mat pc, void * ctx);
69 : PetscErrorCode __libmesh_petsc_snes_postcheck(SNESLineSearch, Vec x, Vec y, Vec w, PetscBool * changed_y, PetscBool * changed_w, void * context);
70 : #endif
71 : }
72 :
73 : /**
74 : * This class provides an interface to PETSc
75 : * iterative solvers that is compatible with the \p libMesh
76 : * \p NonlinearSolver<>
77 : *
78 : * \author Benjamin Kirk
79 : * \date 2002-2007
80 : */
81 : template <typename T>
82 210 : class PetscNonlinearSolver : public NonlinearSolver<T>
83 : {
84 : public:
85 : /**
86 : * The type of system
87 : */
88 : typedef NonlinearImplicitSystem sys_type;
89 :
90 : /**
91 : * Constructor. Initializes Petsc data structures
92 : */
93 : explicit
94 : PetscNonlinearSolver (sys_type & system);
95 :
96 : /**
97 : * Destructor.
98 : */
99 : ~PetscNonlinearSolver ();
100 :
101 : /**
102 : * Release all memory and clear data structures.
103 : */
104 : virtual void clear () override;
105 :
106 : /**
107 : * Initialize data structures if not done so already.
108 : * May assign a name to the solver in some implementations
109 : */
110 : virtual void init (const char * name = nullptr) override;
111 :
112 : /**
113 : * \returns The raw PETSc snes context pointer. This method calls init() so in that vein, we have
114 : * an optional name prefix argument that if provided will be given to the SNES context
115 : */
116 : SNES snes(const char * name = nullptr);
117 :
118 : /**
119 : * Call the Petsc solver. It calls the method below, using the
120 : * same matrix for the system and preconditioner matrices.
121 : */
122 : virtual std::pair<unsigned int, Real>
123 : solve (SparseMatrix<T> &, // System Jacobian Matrix
124 : NumericVector<T> &, // Solution vector
125 : NumericVector<T> &, // Residual vector
126 : const double, // Stopping tolerance
127 : const unsigned int) override; // N. Iterations
128 :
129 : /**
130 : * Prints a useful message about why the latest nonlinear solve
131 : * con(di)verged.
132 : */
133 : virtual void print_converged_reason() override;
134 :
135 : /**
136 : * \returns The currently-available (or most recently obtained, if
137 : * the SNES object has been destroyed) convergence reason.
138 : *
139 : * Refer to PETSc docs for the meaning of different
140 : * SNESConvergedReasons.
141 : */
142 : SNESConvergedReason get_converged_reason();
143 :
144 : /**
145 : * Get the total number of linear iterations done in the last solve
146 : */
147 : virtual int get_total_linear_iterations() override;
148 :
149 : /**
150 : * \returns The current nonlinear iteration number if called
151 : * *during* the solve(), for example by the user-specified residual
152 : * or Jacobian function.
153 : */
154 0 : virtual unsigned get_current_nonlinear_iteration_number() const override
155 0 : { return _current_nonlinear_iteration_number; }
156 :
157 : /**
158 : * Set if the residual should be zeroed out in the callback
159 : */
160 0 : void set_residual_zero_out(bool state) { _zero_out_residual = state; }
161 :
162 : /**
163 : * Set if the jacobian should be zeroed out in the callback
164 : */
165 0 : void set_jacobian_zero_out(bool state) { _zero_out_jacobian = state; }
166 :
167 : /**
168 : * Set to true to use the libMesh's default monitor, set to false to use your own
169 : */
170 0 : void use_default_monitor(bool state) { _default_monitor = state; }
171 :
172 : /**
173 : * Set to true to let PETSc reuse the base vector
174 : */
175 0 : void set_snesmf_reuse_base(bool state) { _snesmf_reuse_base = state; }
176 :
177 : /**
178 : * @return Whether we are reusing the nonlinear function evaluation as the base for doing
179 : * matrix-free approximation of the Jacobian action
180 : */
181 2398 : bool snes_mf_reuse_base() const { return _snesmf_reuse_base; }
182 :
183 : /**
184 : * Set whether we are computing the base vector for matrix-free finite-differencing
185 : */
186 0 : void set_computing_base_vector(bool computing_base_vector) { _computing_base_vector = computing_base_vector; }
187 :
188 : /**
189 : * @return whether we are computing the base vector for matrix-free finite-differencing
190 : */
191 0 : bool computing_base_vector() const { return _computing_base_vector; }
192 :
193 : /**
194 : * Abstract base class to be used to implement a custom line-search algorithm
195 : */
196 : class ComputeLineSearchObject
197 : {
198 : public:
199 : virtual ~ComputeLineSearchObject () = default;
200 :
201 : virtual void linesearch (SNESLineSearch linesearch) = 0;
202 : };
203 :
204 : /**
205 : * A callable object that can be used to specify a custom line-search
206 : */
207 : std::unique_ptr<ComputeLineSearchObject> linesearch_object;
208 :
209 : /**
210 : * Setup the default monitor if required
211 : */
212 : void setup_default_monitor();
213 :
214 : /**
215 : * Getter for preconditioner reuse
216 : */
217 : virtual bool reuse_preconditioner() const override;
218 :
219 : /**
220 : * Getter for the maximum iterations flag for preconditioner reuse
221 : */
222 : virtual unsigned int reuse_preconditioner_max_linear_its() const override;
223 :
224 : /**
225 : * Immediately force a new preconditioner, even if reuse is set
226 : */
227 : virtual void force_new_preconditioner() override;
228 :
229 : protected:
230 :
231 : /**
232 : * Nonlinear solver context
233 : */
234 : WrappedPetsc<SNES> _snes;
235 :
236 : /**
237 : * Store the reason for SNES convergence/divergence for use even after the _snes
238 : * has been cleared.
239 : *
240 : * \note \p print_converged_reason() will always \e try to get the
241 : * current reason with SNESGetConvergedReason(), but if the SNES
242 : * object has already been cleared, it will fall back on this stored
243 : * value. This value is therefore necessarily \e not cleared by the
244 : * \p clear() function.
245 : */
246 : SNESConvergedReason _reason;
247 :
248 : /**
249 : * Stores the total number of linear iterations from the last solve.
250 : */
251 : PetscInt _n_linear_iterations;
252 :
253 : /**
254 : * Stores the current nonlinear iteration number
255 : */
256 : unsigned _current_nonlinear_iteration_number;
257 :
258 : /**
259 : * true to zero out residual before going into application level call-back, otherwise false
260 : */
261 : bool _zero_out_residual;
262 :
263 : /**
264 : * true to zero out jacobian before going into application level call-back, otherwise false
265 : */
266 : bool _zero_out_jacobian;
267 :
268 : /**
269 : * true if we want the default monitor to be set, false for no monitor (i.e. user code can use their own)
270 : */
271 : bool _default_monitor;
272 :
273 : /**
274 : * True, If we want the base vector to be used for differencing even if the function provided to SNESSetFunction()
275 : * is not the same as that provided to MatMFFDSetFunction().
276 : * https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/SNES/MatSNESMFSetReuseBase.html
277 : */
278 : bool _snesmf_reuse_base;
279 :
280 : void build_mat_null_space(NonlinearImplicitSystem::ComputeVectorSubspace * computeSubspaceObject,
281 : void (*)(std::vector<NumericVector<Number> *> &, sys_type &),
282 : MatNullSpace *);
283 :
284 : /**
285 : * Whether we are computing the base vector for matrix-free finite differencing
286 : */
287 : bool _computing_base_vector;
288 :
289 : /**
290 : * Whether we've triggered the preconditioner reuse
291 : */
292 : bool _setup_reuse;
293 :
294 : #if defined(LIBMESH_ENABLE_AMR) && defined(LIBMESH_HAVE_METAPHYSICL)
295 : /**
296 : * Wrapper object for interacting with the "new" libMesh PETSc DM. The new libMesh PETSc DM
297 : * implementation is capable of geometric multigrid while the old implementation is not. The new
298 : * implementation can be activated from the command line with --use_petsc_dm
299 : */
300 : PetscDMWrapper _dm_wrapper;
301 : #endif
302 :
303 : /// Wrapper for matrix-free finite-difference Jacobians
304 : PetscMFFDMatrix<Number> _mffd_jac;
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
|
