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 :
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. It calls the method below, using the
111 : * same matrix for the system and preconditioner matrices.
112 : */
113 : virtual std::pair<unsigned int, Real>
114 : solve (SparseMatrix<T> &, // System Jacobian 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 : * Prints a useful message about why the latest nonlinear solve
122 : * con(di)verged.
123 : */
124 : virtual void print_converged_reason() override;
125 :
126 : /**
127 : * \returns The currently-available (or most recently obtained, if
128 : * the SNES object has been destroyed) convergence reason.
129 : *
130 : * Refer to PETSc docs for the meaning of different
131 : * SNESConvergedReasons.
132 : */
133 : SNESConvergedReason get_converged_reason();
134 :
135 : /**
136 : * Get the total number of linear iterations done in the last solve
137 : */
138 : virtual int get_total_linear_iterations() override;
139 :
140 : /**
141 : * \returns The current nonlinear iteration number if called
142 : * *during* the solve(), for example by the user-specified residual
143 : * or Jacobian function.
144 : */
145 0 : virtual unsigned get_current_nonlinear_iteration_number() const override
146 0 : { return _current_nonlinear_iteration_number; }
147 :
148 : /**
149 : * Set if the residual should be zeroed out in the callback
150 : */
151 0 : void set_residual_zero_out(bool state) { _zero_out_residual = state; }
152 :
153 : /**
154 : * Set if the jacobian should be zeroed out in the callback
155 : */
156 0 : void set_jacobian_zero_out(bool state) { _zero_out_jacobian = state; }
157 :
158 : /**
159 : * Set to true to use the libMesh's default monitor, set to false to use your own
160 : */
161 0 : void use_default_monitor(bool state) { _default_monitor = state; }
162 :
163 : /**
164 : * Set to true to let PETSc reuse the base vector
165 : */
166 0 : void set_snesmf_reuse_base(bool state) { _snesmf_reuse_base = state; }
167 :
168 : /**
169 : * @return Whether we are reusing the nonlinear function evaluation as the base for doing
170 : * matrix-free approximation of the Jacobian action
171 : */
172 2398 : bool snes_mf_reuse_base() const { return _snesmf_reuse_base; }
173 :
174 : /**
175 : * Set whether we are computing the base vector for matrix-free finite-differencing
176 : */
177 0 : void set_computing_base_vector(bool computing_base_vector) { _computing_base_vector = computing_base_vector; }
178 :
179 : /**
180 : * @return whether we are computing the base vector for matrix-free finite-differencing
181 : */
182 0 : bool computing_base_vector() const { return _computing_base_vector; }
183 :
184 : /**
185 : * Abstract base class to be used to implement a custom line-search algorithm
186 : */
187 : class ComputeLineSearchObject
188 : {
189 : public:
190 : virtual ~ComputeLineSearchObject () = default;
191 :
192 : virtual void linesearch (SNESLineSearch linesearch) = 0;
193 : };
194 :
195 : /**
196 : * A callable object that can be used to specify a custom line-search
197 : */
198 : std::unique_ptr<ComputeLineSearchObject> linesearch_object;
199 :
200 : /**
201 : * Setup the default monitor if required
202 : */
203 : void setup_default_monitor();
204 :
205 : /**
206 : * Getter for preconditioner reuse
207 : */
208 : virtual bool reuse_preconditioner() const override;
209 :
210 : /**
211 : * Getter for the maximum iterations flag for preconditioner reuse
212 : */
213 : virtual unsigned int reuse_preconditioner_max_linear_its() const override;
214 :
215 : /**
216 : * Immediately force a new preconditioner, even if reuse is set
217 : */
218 : virtual void force_new_preconditioner() override;
219 :
220 : protected:
221 :
222 : /**
223 : * Nonlinear solver context
224 : */
225 : WrappedPetsc<SNES> _snes;
226 :
227 : /**
228 : * Store the reason for SNES convergence/divergence for use even after the _snes
229 : * has been cleared.
230 : *
231 : * \note \p print_converged_reason() will always \e try to get the
232 : * current reason with SNESGetConvergedReason(), but if the SNES
233 : * object has already been cleared, it will fall back on this stored
234 : * value. This value is therefore necessarily \e not cleared by the
235 : * \p clear() function.
236 : */
237 : SNESConvergedReason _reason;
238 :
239 : /**
240 : * Stores the total number of linear iterations from the last solve.
241 : */
242 : PetscInt _n_linear_iterations;
243 :
244 : /**
245 : * Stores the current nonlinear iteration number
246 : */
247 : unsigned _current_nonlinear_iteration_number;
248 :
249 : /**
250 : * true to zero out residual before going into application level call-back, otherwise false
251 : */
252 : bool _zero_out_residual;
253 :
254 : /**
255 : * true to zero out jacobian before going into application level call-back, otherwise false
256 : */
257 : bool _zero_out_jacobian;
258 :
259 : /**
260 : * true if we want the default monitor to be set, false for no monitor (i.e. user code can use their own)
261 : */
262 : bool _default_monitor;
263 :
264 : /**
265 : * True, If we want the base vector to be used for differencing even if the function provided to SNESSetFunction()
266 : * is not the same as that provided to MatMFFDSetFunction().
267 : * https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/SNES/MatSNESMFSetReuseBase.html
268 : */
269 : bool _snesmf_reuse_base;
270 :
271 : void build_mat_null_space(NonlinearImplicitSystem::ComputeVectorSubspace * computeSubspaceObject,
272 : void (*)(std::vector<NumericVector<Number> *> &, sys_type &),
273 : MatNullSpace *);
274 :
275 : /**
276 : * Whether we are computing the base vector for matrix-free finite differencing
277 : */
278 : bool _computing_base_vector;
279 :
280 : /**
281 : * Whether we've triggered the preconditioner reuse
282 : */
283 : bool _setup_reuse;
284 :
285 : #if defined(LIBMESH_ENABLE_AMR) && defined(LIBMESH_HAVE_METAPHYSICL)
286 : /**
287 : * Wrapper object for interacting with the "new" libMesh PETSc DM. The new libMesh PETSc DM
288 : * implementation is capable of geometric multigrid while the old implementation is not. The new
289 : * implementation can be activated from the command line with --use_petsc_dm
290 : */
291 : PetscDMWrapper _dm_wrapper;
292 : #endif
293 :
294 : /// Wrapper for matrix-free finite-difference Jacobians
295 : PetscMFFDMatrix<Number> _mffd_jac;
296 :
297 : private:
298 : friend ResidualContext libmesh_petsc_snes_residual_helper (SNES snes, Vec x, void * ctx);
299 : friend PetscErrorCode libmesh_petsc_snes_residual (SNES snes, Vec x, Vec r, void * ctx);
300 : friend PetscErrorCode libmesh_petsc_snes_fd_residual (SNES snes, Vec x, Vec r, void * ctx);
301 : friend PetscErrorCode libmesh_petsc_snes_mffd_residual (SNES snes, Vec x, Vec r, void * ctx);
302 : friend PetscErrorCode libmesh_petsc_snes_jacobian (SNES snes, Vec x, Mat jac, Mat pc, void * ctx);
303 : friend PetscErrorCode libmesh_petsc_snes_precheck (SNESLineSearch,
304 : Vec X,
305 : Vec Y,
306 : PetscBool * changed,
307 : void * context);
308 : };
309 :
310 : } // namespace libMesh
311 :
312 :
313 : #endif // #ifdef LIBMESH_HAVE_PETSC
314 : #endif // LIBMESH_PETSC_NONLINEAR_SOLVER_H
|
