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 : #include "libmesh/libmesh_common.h"
21 :
22 : #ifdef LIBMESH_HAVE_PETSC
23 :
24 : // Local Includes
25 : #include "libmesh/libmesh_logging.h"
26 : #include "libmesh/nonlinear_implicit_system.h"
27 : #include "libmesh/petsc_nonlinear_solver.h"
28 : #include "libmesh/petsc_linear_solver.h"
29 : #include "libmesh/petsc_vector.h"
30 : #include "libmesh/petsc_mffd_matrix.h"
31 : #include "libmesh/dof_map.h"
32 : #include "libmesh/preconditioner.h"
33 : #include "libmesh/solver_configuration.h"
34 : #include "libmesh/petscdmlibmesh.h"
35 : #include "libmesh/petsc_mffd_matrix.h"
36 :
37 : #if defined(LIBMESH_HAVE_PETSC_HYPRE) && PETSC_VERSION_LESS_THAN(3, 23, 0) && \
38 : !PETSC_VERSION_LESS_THAN(3, 12, 0) && defined(PETSC_HAVE_HYPRE_DEVICE)
39 : #include <HYPRE_utilities.h>
40 : #endif
41 :
42 : namespace libMesh
43 : {
44 : class ResidualContext
45 : {
46 : public:
47 5660 : ResidualContext(PetscNonlinearSolver<Number> * solver_in, NonlinearImplicitSystem & sys_in) :
48 : solver(solver_in),
49 5660 : sys(sys_in)
50 5660 : {}
51 :
52 : PetscNonlinearSolver<Number> * solver;
53 : NonlinearImplicitSystem & sys;
54 : };
55 :
56 : ResidualContext
57 281117 : libmesh_petsc_snes_residual_helper (SNES snes, Vec x, void * ctx)
58 : {
59 11320 : LOG_SCOPE("residual()", "PetscNonlinearSolver");
60 :
61 5660 : libmesh_assert(x);
62 5660 : libmesh_assert(ctx);
63 :
64 : // No way to safety-check this cast, since we got a void *...
65 5660 : PetscNonlinearSolver<Number> * solver =
66 : static_cast<PetscNonlinearSolver<Number> *> (ctx);
67 :
68 5660 : libmesh_parallel_only(solver->comm());
69 :
70 : // Get the current iteration number from the snes object,
71 : // store it in the PetscNonlinearSolver object for possible use
72 : // by the user's residual function.
73 : {
74 281117 : PetscInt n_iterations = 0;
75 281117 : LibmeshPetscCall2(solver->comm(), SNESGetIterationNumber(snes, &n_iterations));
76 281117 : solver->_current_nonlinear_iteration_number = cast_int<unsigned>(n_iterations);
77 : }
78 :
79 10912 : NonlinearImplicitSystem & sys = solver->system();
80 :
81 5660 : PetscVector<Number> & X_sys = *cast_ptr<PetscVector<Number> *>(sys.solution.get());
82 :
83 281117 : PetscVector<Number> X_global(x, sys.comm());
84 :
85 : // Use the system's update() to get a good local version of the
86 : // parallel solution. This operation does not modify the incoming
87 : // "x" vector, it only localizes information from "x" into
88 : // sys.current_local_solution.
89 281117 : X_global.swap(X_sys);
90 281117 : sys.update();
91 281117 : X_global.swap(X_sys);
92 :
93 : // Enforce constraints (if any) exactly on the
94 : // current_local_solution. This is the solution vector that is
95 : // actually used in the computation of the residual below, and is
96 : // not locked by debug-enabled PETSc the way that "x" is.
97 281117 : if (solver->_exact_constraint_enforcement)
98 281117 : sys.get_dof_map().enforce_constraints_exactly(sys, sys.current_local_solution.get());
99 :
100 292029 : return ResidualContext(solver, sys);
101 270205 : }
102 :
103 : //--------------------------------------------------------------------
104 : // Functions with C linkage to pass to PETSc. PETSc will call these
105 : // methods as needed.
106 : //
107 : // Since they must have C linkage they have no knowledge of a namespace.
108 : // Give them an obscure name to avoid namespace pollution.
109 : extern "C"
110 : {
111 : // -----------------------------------------------------------------
112 : // this function monitors the nonlinear solve and checks to see
113 : // if we want to recalculate the preconditioner. It only gets
114 : // added to the SNES instance if we're reusing the preconditioner
115 : PetscErrorCode
116 0 : libmesh_petsc_recalculate_monitor(SNES snes, PetscInt, PetscReal, void* ctx)
117 : {
118 : PetscFunctionBegin;
119 :
120 : // No way to safety-check this cast, since we got a void *...
121 0 : PetscNonlinearSolver<Number> * solver =
122 : static_cast<PetscNonlinearSolver<Number> *> (ctx);
123 :
124 : KSP ksp;
125 0 : LibmeshPetscCall2(solver->comm(), SNESGetKSP(snes, &ksp));
126 :
127 : PetscInt niter;
128 0 : LibmeshPetscCall2(solver->comm(), KSPGetIterationNumber(ksp, &niter));
129 :
130 0 : if (niter > cast_int<PetscInt>(solver->reuse_preconditioner_max_linear_its()))
131 : {
132 : // -2 is a magic number for "recalculate next time you need it
133 : // and then not again"
134 0 : LibmeshPetscCall2(solver->comm(), SNESSetLagPreconditioner(snes, -2));
135 : }
136 0 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
137 : }
138 :
139 : //-------------------------------------------------------------------
140 : // this function is called by PETSc at the end of each nonlinear step
141 : PetscErrorCode
142 105241 : libmesh_petsc_snes_monitor (SNES, PetscInt its, PetscReal fnorm, void *)
143 : {
144 : PetscFunctionBegin;
145 2896 : libMesh::out << " NL step "
146 2896 : << std::setw(2) << its
147 2896 : << std::scientific
148 2896 : << ", |residual|_2 = " << fnorm
149 2896 : << std::endl;
150 105241 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
151 : }
152 :
153 : #ifdef LIBMESH_ENABLE_DEPRECATED
154 : PetscErrorCode
155 0 : __libmesh_petsc_snes_monitor (SNES, PetscInt its, PetscReal fnorm, void *)
156 : {
157 : PetscFunctionBegin;
158 : libmesh_deprecated();
159 0 : PetscFunctionReturn(libmesh_petsc_snes_monitor(nullptr, its, fnorm, nullptr));
160 : }
161 : #endif
162 :
163 :
164 : //---------------------------------------------------------------
165 : // this function is called by PETSc to evaluate the residual at X
166 : PetscErrorCode
167 159137 : libmesh_petsc_snes_residual (SNES snes, Vec x, Vec r, void * ctx)
168 : {
169 : PetscFunctionBegin;
170 :
171 159137 : ResidualContext rc = libmesh_petsc_snes_residual_helper(snes, x, ctx);
172 :
173 3262 : libmesh_parallel_only(rc.sys.comm());
174 :
175 3262 : libmesh_assert(r);
176 165661 : PetscVector<Number> R(r, rc.sys.comm());
177 :
178 159137 : if (rc.solver->_zero_out_residual)
179 159137 : R.zero();
180 :
181 : //-----------------------------------------------------------------------------
182 : // if the user has provided both function pointers and objects only the pointer
183 : // will be used, so catch that as an error
184 159137 : libmesh_error_msg_if(rc.solver->residual && rc.solver->residual_object,
185 : "ERROR: cannot specify both a function and object to compute the Residual!");
186 :
187 159137 : libmesh_error_msg_if(rc.solver->matvec && rc.solver->residual_and_jacobian_object,
188 : "ERROR: cannot specify both a function and object to compute the combined Residual & Jacobian!");
189 :
190 159137 : if (rc.solver->residual != nullptr)
191 16 : rc.solver->residual(*rc.sys.current_local_solution.get(), R, rc.sys);
192 :
193 158857 : else if (rc.solver->residual_object != nullptr)
194 3277 : rc.solver->residual_object->residual(*rc.sys.current_local_solution.get(), R, rc.sys);
195 :
196 155678 : else if (rc.solver->matvec != nullptr)
197 0 : rc.solver->matvec (*rc.sys.current_local_solution.get(), &R, nullptr, rc.sys);
198 :
199 155678 : else if (rc.solver->residual_and_jacobian_object != nullptr)
200 : {
201 155678 : auto & jac = rc.sys.get_system_matrix();
202 :
203 155678 : if (rc.solver->_zero_out_jacobian)
204 155678 : jac.zero();
205 :
206 155678 : rc.solver->residual_and_jacobian_object->residual_and_jacobian(
207 6312 : *rc.sys.current_local_solution.get(), &R, &jac, rc.sys);
208 :
209 155678 : jac.close();
210 155678 : if (rc.solver->_exact_constraint_enforcement)
211 : {
212 155678 : rc.sys.get_dof_map().enforce_constraints_on_jacobian(rc.sys, &jac);
213 155678 : jac.close();
214 : }
215 : }
216 :
217 : else
218 0 : libmesh_error_msg("Error! Unable to compute residual and/or Jacobian!");
219 :
220 :
221 : // Synchronize PETSc x to local solution since the local solution may be changed due to the constraints
222 3262 : PetscVector<Number> & X_sys = *cast_ptr<PetscVector<Number> *>(rc.sys.solution.get());
223 162399 : PetscVector<Number> X_global(x, rc.sys.comm());
224 :
225 159137 : X_global.swap(X_sys);
226 159137 : rc.sys.update();
227 159137 : X_global.swap(X_sys);
228 :
229 159137 : R.close();
230 :
231 159137 : if (rc.solver->_exact_constraint_enforcement)
232 : {
233 159137 : rc.sys.get_dof_map().enforce_constraints_on_residual(rc.sys, &R, rc.sys.current_local_solution.get());
234 159137 : R.close();
235 : }
236 :
237 162399 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
238 152613 : }
239 :
240 : #ifdef LIBMESH_ENABLE_DEPRECATED
241 : PetscErrorCode
242 0 : __libmesh_petsc_snes_residual (SNES snes, Vec x, Vec r, void * ctx)
243 : {
244 : PetscFunctionBegin;
245 : libmesh_deprecated();
246 0 : PetscFunctionReturn(libmesh_petsc_snes_residual(snes, x, r, ctx));
247 : }
248 : #endif
249 :
250 : //-----------------------------------------------------------------------------------------
251 : // this function is called by PETSc to approximate the Jacobian at X via finite differences
252 : PetscErrorCode
253 0 : libmesh_petsc_snes_fd_residual (SNES snes, Vec x, Vec r, void * ctx)
254 : {
255 : PetscFunctionBegin;
256 :
257 0 : ResidualContext rc = libmesh_petsc_snes_residual_helper(snes, x, ctx);
258 :
259 0 : libmesh_parallel_only(rc.sys.comm());
260 :
261 0 : libmesh_assert(r);
262 0 : PetscVector<Number> R(r, rc.sys.comm());
263 :
264 0 : if (rc.solver->_zero_out_residual)
265 0 : R.zero();
266 :
267 0 : if (rc.solver->fd_residual_object != nullptr)
268 0 : rc.solver->fd_residual_object->residual(*rc.sys.current_local_solution.get(), R, rc.sys);
269 :
270 0 : else if (rc.solver->residual_object != nullptr)
271 0 : rc.solver->residual_object->residual(*rc.sys.current_local_solution.get(), R, rc.sys);
272 :
273 : else
274 0 : libmesh_error_msg("Error! Unable to compute residual for forming finite difference Jacobian!");
275 :
276 : // Synchronize PETSc x to local solution since the local solution may be changed due to the constraints
277 0 : PetscVector<Number> & X_sys = *cast_ptr<PetscVector<Number> *>(rc.sys.solution.get());
278 0 : PetscVector<Number> X_global(x, rc.sys.comm());
279 :
280 0 : X_global.swap(X_sys);
281 0 : rc.sys.update();
282 0 : X_global.swap(X_sys);
283 :
284 0 : R.close();
285 :
286 0 : if (rc.solver->_exact_constraint_enforcement)
287 : {
288 0 : rc.sys.get_dof_map().enforce_constraints_on_residual(rc.sys, &R, rc.sys.current_local_solution.get());
289 0 : R.close();
290 : }
291 :
292 0 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
293 0 : }
294 :
295 : #ifdef LIBMESH_ENABLE_DEPRECATED
296 : PetscErrorCode
297 0 : __libmesh_petsc_snes_fd_residual (SNES snes, Vec x, Vec r, void * ctx)
298 : {
299 : PetscFunctionBegin;
300 : libmesh_deprecated();
301 0 : PetscFunctionReturn(libmesh_petsc_snes_fd_residual(snes, x, r, ctx));
302 : }
303 : #endif
304 :
305 : //----------------------------------------------------------------
306 : // this function is called by PETSc to approximate Jacobian-vector
307 : // products at X via finite differences
308 : PetscErrorCode
309 121980 : libmesh_petsc_snes_mffd_residual (SNES snes, Vec x, Vec r, void * ctx)
310 : {
311 : PetscFunctionBegin;
312 :
313 121980 : ResidualContext rc = libmesh_petsc_snes_residual_helper(snes, x, ctx);
314 :
315 2398 : libmesh_parallel_only(rc.sys.comm());
316 :
317 2398 : libmesh_assert(r);
318 126368 : PetscVector<Number> R(r, rc.sys.comm());
319 :
320 121980 : if (rc.solver->_zero_out_residual)
321 121980 : R.zero();
322 :
323 121980 : if (rc.solver->mffd_residual_object != nullptr)
324 123970 : rc.solver->mffd_residual_object->residual(*rc.sys.current_local_solution.get(), R, rc.sys);
325 :
326 0 : else if (rc.solver->residual_object != nullptr)
327 0 : rc.solver->residual_object->residual(*rc.sys.current_local_solution.get(), R, rc.sys);
328 :
329 : else
330 0 : libmesh_error_msg("Error! Unable to compute residual for forming finite differenced"
331 : "Jacobian-vector products!");
332 :
333 : // Synchronize PETSc x to local solution since the local solution may be changed due to the constraints
334 2398 : PetscVector<Number> & X_sys = *cast_ptr<PetscVector<Number> *>(rc.sys.solution.get());
335 123970 : PetscVector<Number> X_global(x, rc.sys.comm());
336 :
337 121980 : X_global.swap(X_sys);
338 121980 : rc.sys.update();
339 121980 : X_global.swap(X_sys);
340 :
341 121980 : R.close();
342 :
343 121980 : if (rc.solver->_exact_constraint_enforcement)
344 : {
345 121980 : rc.sys.get_dof_map().enforce_constraints_on_residual(rc.sys, &R, rc.sys.current_local_solution.get());
346 121980 : R.close();
347 : }
348 :
349 124378 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
350 117592 : }
351 :
352 : //----------------------------------------------------------
353 : // this function serves an interface between the petsc layer
354 : // and the actual mffd residual computing routine
355 : PetscErrorCode
356 121980 : libmesh_petsc_snes_mffd_interface (void * ctx, Vec x, Vec r)
357 : {
358 : PetscFunctionBegin;
359 :
360 : // No way to safety-check this cast, since we got a void *...
361 2398 : PetscNonlinearSolver<Number> * solver =
362 : static_cast<PetscNonlinearSolver<Number> *> (ctx);
363 :
364 121980 : LibmeshPetscCall2(solver->comm(), libmesh_petsc_snes_mffd_residual(solver->snes(), x, r, ctx));
365 :
366 : #if !PETSC_VERSION_LESS_THAN(3,8,4)
367 : #ifndef NDEBUG
368 :
369 : // When the user requested to reuse the nonlinear residual as the base for doing matrix-free
370 : // approximation of the Jacobian, we'll do a sanity check to make sure that that was safe to do
371 2398 : if (solver->snes_mf_reuse_base() && (solver->comm().size() == 1) && (libMesh::n_threads() == 1))
372 : {
373 0 : SNES snes = solver->snes();
374 :
375 : KSP ksp;
376 0 : LibmeshPetscCall2(solver->comm(), SNESGetKSP(snes, &ksp));
377 :
378 : PetscInt ksp_it;
379 0 : LibmeshPetscCall2(solver->comm(), KSPGetIterationNumber(ksp, &ksp_it));
380 :
381 : SNESType snes_type;
382 0 : LibmeshPetscCall2(solver->comm(), SNESGetType(snes, &snes_type));
383 :
384 0 : libmesh_assert_msg(snes_type, "We're being called from SNES; snes_type should be non-null");
385 :
386 : Mat J;
387 0 : LibmeshPetscCall2(solver->comm(), SNESGetJacobian(snes, &J, NULL, NULL, NULL));
388 0 : libmesh_assert_msg(J, "We're being called from SNES; J should be non-null");
389 :
390 : MatType mat_type;
391 0 : LibmeshPetscCall2(solver->comm(), MatGetType(J, &mat_type));
392 0 : libmesh_assert_msg(mat_type, "We're being called from SNES; mat_type should be non-null");
393 :
394 0 : bool is_operator_mffd = strcmp(mat_type, MATMFFD) == 0;
395 :
396 0 : if ((ksp_it == PetscInt(0)) && is_operator_mffd)
397 : {
398 0 : bool computing_base_vector = solver->computing_base_vector();
399 :
400 0 : if (computing_base_vector)
401 : {
402 : Vec nonlinear_residual;
403 :
404 0 : LibmeshPetscCall2(solver->comm(), SNESGetFunction(snes, &nonlinear_residual, NULL, NULL));
405 :
406 : PetscBool vecs_equal;
407 0 : LibmeshPetscCall2(solver->comm(), VecEqual(r, nonlinear_residual, &vecs_equal));
408 :
409 0 : libmesh_error_msg_if(!(vecs_equal == PETSC_TRUE),
410 : "You requested to reuse the nonlinear residual vector as the base vector for "
411 : "computing the action of the matrix-free Jacobian, but the vectors are not "
412 : "the same. Your physics must have states; either remove the states "
413 : "from your code or make sure that you set_mf_reuse_base(false)");
414 : }
415 :
416 : // There are always exactly two function evaluations for the zeroth ksp iteration when doing
417 : // matrix-free approximation of the Jacobian action: one corresponding to the evaluation of
418 : // the base vector, and the other corresponding to evaluation of the perturbed vector. So we
419 : // toggle back and forth between states
420 0 : solver->set_computing_base_vector(!computing_base_vector);
421 : }
422 : }
423 : #endif
424 : #endif
425 :
426 121980 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
427 : }
428 :
429 : #ifdef LIBMESH_ENABLE_DEPRECATED
430 : PetscErrorCode
431 0 : __libmesh_petsc_snes_mffd_interface (void * ctx, Vec x, Vec r)
432 : {
433 : PetscFunctionBegin;
434 : libmesh_deprecated();
435 0 : PetscFunctionReturn(libmesh_petsc_snes_mffd_interface(ctx, x, r));
436 : }
437 : #endif
438 :
439 : //---------------------------------------------------------------
440 : // this function is called by PETSc to evaluate the Jacobian at X
441 : PetscErrorCode
442 76064 : libmesh_petsc_snes_jacobian(SNES snes, Vec x, Mat jac, Mat pc, void * ctx)
443 : {
444 : PetscFunctionBegin;
445 :
446 4120 : LOG_SCOPE("jacobian()", "PetscNonlinearSolver");
447 :
448 2060 : libmesh_assert(ctx);
449 :
450 : // No way to safety-check this cast, since we got a void *...
451 2060 : PetscNonlinearSolver<Number> * solver =
452 : static_cast<PetscNonlinearSolver<Number> *> (ctx);
453 :
454 2060 : libmesh_parallel_only(solver->comm());
455 :
456 : // Get the current iteration number from the snes object,
457 : // store it in the PetscNonlinearSolver object for possible use
458 : // by the user's Jacobian function.
459 : {
460 76064 : PetscInt n_iterations = 0;
461 76064 : LibmeshPetscCall2(solver->comm(), SNESGetIterationNumber(snes, &n_iterations));
462 76064 : solver->_current_nonlinear_iteration_number = cast_int<unsigned>(n_iterations);
463 : }
464 :
465 : //-----------------------------------------------------------------------------
466 : // if the user has provided both function pointers and objects only the pointer
467 : // will be used, so catch that as an error
468 76064 : libmesh_error_msg_if(solver->jacobian && solver->jacobian_object,
469 : "ERROR: cannot specify both a function and object to compute the Jacobian!");
470 :
471 76064 : libmesh_error_msg_if(solver->matvec && solver->residual_and_jacobian_object,
472 : "ERROR: cannot specify both a function and object to compute the combined Residual & Jacobian!");
473 :
474 4120 : NonlinearImplicitSystem & sys = solver->system();
475 :
476 76064 : PetscMatrixBase<Number> * const PC = pc ? PetscMatrixBase<Number>::get_context(pc, sys.comm()) : nullptr;
477 76064 : PetscMatrixBase<Number> * Jac = jac ? PetscMatrixBase<Number>::get_context(jac, sys.comm()) : nullptr;
478 2060 : PetscVector<Number> & X_sys = *cast_ptr<PetscVector<Number> *>(sys.solution.get());
479 80184 : PetscVector<Number> X_global(x, sys.comm());
480 :
481 6180 : PetscMFFDMatrix<Number> mffd_jac(sys.comm());
482 76064 : PetscBool p_is_shell = PETSC_FALSE;
483 76064 : PetscBool j_is_mffd = PETSC_FALSE;
484 76064 : PetscBool j_is_shell = PETSC_FALSE;
485 76064 : if (pc)
486 76064 : LibmeshPetscCall2(sys.comm(), PetscObjectTypeCompare((PetscObject)pc, MATSHELL, &p_is_shell));
487 2060 : libmesh_assert(jac);
488 76064 : LibmeshPetscCall2(sys.comm(), PetscObjectTypeCompare((PetscObject)jac, MATMFFD, &j_is_mffd));
489 76064 : LibmeshPetscCall2(sys.comm(), PetscObjectTypeCompare((PetscObject)jac, MATSHELL, &j_is_shell));
490 76064 : if (j_is_mffd == PETSC_TRUE)
491 : {
492 408 : libmesh_assert(!Jac);
493 408 : Jac = &mffd_jac;
494 408 : mffd_jac = jac;
495 : }
496 :
497 : // We already computed the Jacobian during the residual evaluation
498 76064 : if (solver->residual_and_jacobian_object)
499 : {
500 : // We could be doing matrix-free in which case we cannot rely on closing of explicit matrices
501 : // that occurs during the PETSc residual callback
502 73312 : if ((j_is_shell == PETSC_TRUE) || (j_is_mffd == PETSC_TRUE))
503 14154 : Jac->close();
504 :
505 73312 : if (pc && (p_is_shell == PETSC_TRUE))
506 0 : PC->close();
507 :
508 73312 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
509 : }
510 :
511 : // Set the dof maps
512 2836 : PC->attach_dof_map(sys.get_dof_map());
513 2836 : Jac->attach_dof_map(sys.get_dof_map());
514 :
515 : // Use the systems update() to get a good local version of the parallel solution
516 2752 : X_global.swap(X_sys);
517 2752 : sys.update();
518 2752 : X_global.swap(X_sys);
519 :
520 : // Enforce constraints (if any) exactly on the
521 : // current_local_solution. This is the solution vector that is
522 : // actually used in the computation of the residual below, and is
523 : // not locked by debug-enabled PETSc the way that "x" is.
524 2752 : if (solver->_exact_constraint_enforcement)
525 2752 : sys.get_dof_map().enforce_constraints_exactly(sys, sys.current_local_solution.get());
526 :
527 2752 : if (solver->_zero_out_jacobian)
528 2752 : PC->zero();
529 :
530 :
531 2752 : if (solver->jacobian != nullptr)
532 210 : solver->jacobian(*sys.current_local_solution.get(), *PC, sys);
533 :
534 2542 : else if (solver->jacobian_object != nullptr)
535 2620 : solver->jacobian_object->jacobian(*sys.current_local_solution.get(), *PC, sys);
536 :
537 0 : else if (solver->matvec != nullptr)
538 0 : solver->matvec(*sys.current_local_solution.get(), nullptr, PC, sys);
539 :
540 : else
541 0 : libmesh_error_msg("Error! Unable to compute residual and/or Jacobian!");
542 :
543 2752 : PC->close();
544 2752 : if (solver->_exact_constraint_enforcement)
545 : {
546 2752 : sys.get_dof_map().enforce_constraints_on_jacobian(sys, PC);
547 2752 : PC->close();
548 : }
549 :
550 2752 : if (Jac != PC)
551 : {
552 : // Assume that shells know what they're doing
553 0 : libmesh_assert(!solver->_exact_constraint_enforcement || (j_is_mffd == PETSC_TRUE) ||
554 : (j_is_shell == PETSC_TRUE));
555 0 : Jac->close();
556 : }
557 :
558 84 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
559 71944 : }
560 :
561 : // This function gets called by PETSc in place of the standard Petsc line searches
562 : // if a linesearch object is supplied to the PetscNonlinearSolver class. It wraps
563 : // the linesearch algorithm implemented on the linesearch object.
564 : // * "linesearch" is an object that can be used to access the non-linear and linear solution
565 : // vectors as well as the residual and SNES object
566 : // * "ctx" is the PetscNonlinearSolver context
567 0 : PetscErrorCode libmesh_petsc_linesearch_shellfunc (SNESLineSearch linesearch, void * ctx)
568 : {
569 : PetscFunctionBegin;
570 :
571 : // No way to safety-check this cast, since we got a void *...
572 0 : PetscNonlinearSolver<Number> * solver =
573 : static_cast<PetscNonlinearSolver<Number> *> (ctx);
574 :
575 0 : libmesh_parallel_only(solver->comm());
576 :
577 0 : solver->linesearch_object->linesearch(linesearch);
578 0 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
579 : }
580 :
581 : #ifdef LIBMESH_ENABLE_DEPRECATED
582 : PetscErrorCode
583 0 : __libmesh_petsc_snes_jacobian(SNES snes, Vec x, Mat jac, Mat pc, void * ctx)
584 : {
585 : PetscFunctionBegin;
586 : libmesh_deprecated();
587 0 : PetscFunctionReturn(libmesh_petsc_snes_jacobian(snes, x, jac, pc, ctx));
588 : }
589 : #endif
590 :
591 : // This function gets called by PETSc after the SNES linesearch is
592 : // complete. We use it to exactly enforce any constraints on the
593 : // solution which may have drifted during the linear solve. In the
594 : // PETSc nomenclature:
595 : // * "x" is the old solution vector,
596 : // * "y" is the search direction (Newton step) vector,
597 : // * "w" is the candidate solution vector, and
598 : // the user is responsible for setting changed_y and changed_w
599 : // appropriately, depending on whether or not the search
600 : // direction or solution vector was changed, respectively.
601 210 : PetscErrorCode libmesh_petsc_snes_postcheck(SNESLineSearch, Vec x, Vec y, Vec w, PetscBool * changed_y, PetscBool * changed_w, void * context)
602 : {
603 : PetscFunctionBegin;
604 :
605 12 : LOG_SCOPE("postcheck()", "PetscNonlinearSolver");
606 :
607 : // PETSc almost certainly initializes these to false already, but
608 : // it doesn't hurt to be explicit.
609 210 : *changed_w = PETSC_FALSE;
610 210 : *changed_y = PETSC_FALSE;
611 :
612 6 : libmesh_assert(context);
613 :
614 : // Cast the context to a NonlinearSolver object.
615 6 : PetscNonlinearSolver<Number> * solver =
616 : static_cast<PetscNonlinearSolver<Number> *> (context);
617 :
618 6 : libmesh_parallel_only(solver->comm());
619 :
620 : // If the user has provided both postcheck function pointer and
621 : // object, this is ambiguous, so throw an error.
622 210 : libmesh_error_msg_if(solver->postcheck && solver->postcheck_object,
623 : "ERROR: cannot specify both a function and object for performing the solve postcheck!");
624 :
625 : // It's also possible that we don't need to do anything at all, in
626 : // that case return early...
627 12 : NonlinearImplicitSystem & sys = solver->system();
628 :
629 210 : if (!solver->postcheck && !solver->postcheck_object)
630 0 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
631 :
632 : // We definitely need to wrap at least "w"
633 216 : PetscVector<Number> petsc_w(w, sys.comm());
634 :
635 : // The user sets these flags in his/her postcheck function to
636 : // indicate whether they changed something.
637 : bool
638 210 : changed_search_direction = false,
639 210 : changed_new_soln = false;
640 :
641 210 : if (solver->postcheck || solver->postcheck_object)
642 : {
643 222 : PetscVector<Number> petsc_x(x, sys.comm());
644 222 : PetscVector<Number> petsc_y(y, sys.comm());
645 :
646 210 : if (solver->postcheck)
647 0 : solver->postcheck(petsc_x,
648 : petsc_y,
649 : petsc_w,
650 : changed_search_direction,
651 : changed_new_soln,
652 : sys);
653 :
654 210 : else if (solver->postcheck_object)
655 216 : solver->postcheck_object->postcheck(petsc_x,
656 : petsc_y,
657 : petsc_w,
658 : changed_search_direction,
659 : changed_new_soln,
660 12 : sys);
661 198 : }
662 :
663 : // Record whether the user changed the solution or the search direction.
664 210 : if (changed_search_direction)
665 0 : *changed_y = PETSC_TRUE;
666 :
667 210 : if (changed_new_soln)
668 0 : *changed_w = PETSC_TRUE;
669 :
670 6 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
671 198 : }
672 :
673 : #ifdef LIBMESH_ENABLE_DEPRECATED
674 0 : PetscErrorCode __libmesh_petsc_snes_postcheck(SNESLineSearch, Vec x, Vec y, Vec w, PetscBool * changed_y, PetscBool * changed_w, void * context)
675 : {
676 : PetscFunctionBegin;
677 : libmesh_deprecated();
678 0 : PetscFunctionReturn(libmesh_petsc_snes_postcheck(nullptr, x, y, w, changed_y, changed_w, context));
679 : }
680 : #endif
681 :
682 0 : PetscErrorCode libmesh_petsc_snes_precheck(SNESLineSearch, Vec X, Vec Y, PetscBool * changed, void * context)
683 : {
684 : PetscFunctionBegin;
685 :
686 0 : LOG_SCOPE("precheck()", "PetscNonlinearSolver");
687 :
688 : // PETSc almost certainly initializes these to false already, but
689 : // it doesn't hurt to be explicit.
690 0 : *changed = PETSC_FALSE;
691 :
692 0 : libmesh_assert(context);
693 :
694 : // Cast the context to a NonlinearSolver object.
695 0 : PetscNonlinearSolver<Number> * solver =
696 : static_cast<PetscNonlinearSolver<Number> *> (context);
697 :
698 0 : libmesh_parallel_only(solver->comm());
699 :
700 : // It's possible that we don't need to do anything at all, in
701 : // that case return early...
702 0 : if (!solver->precheck_object)
703 0 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
704 :
705 : // The user sets these flags in his/her postcheck function to
706 : // indicate whether they changed something.
707 : bool
708 0 : petsc_changed = false;
709 :
710 0 : auto & sys = solver->system();
711 0 : auto & x_sys = *cast_ptr<PetscVector<Number> *>(sys.solution.get());
712 0 : PetscVector<Number> petsc_x(X, sys.comm());
713 0 : PetscVector<Number> petsc_y(Y, sys.comm());
714 :
715 : // Use the systems update() to get a good local version of the parallel solution
716 0 : petsc_x.swap(x_sys);
717 0 : sys.update();
718 0 : petsc_x.swap(x_sys);
719 :
720 : // Enforce constraints (if any) exactly on the
721 : // current_local_solution. This is the solution vector that is
722 : // actually used in the computation of residuals and Jacobians, and is
723 : // not locked by debug-enabled PETSc the way that "x" is.
724 0 : libmesh_assert(sys.current_local_solution.get());
725 0 : auto & local_soln = *sys.current_local_solution.get();
726 0 : if (solver->_exact_constraint_enforcement)
727 0 : sys.get_dof_map().enforce_constraints_exactly(sys, &local_soln);
728 :
729 0 : solver->precheck_object->precheck(local_soln,
730 : petsc_y,
731 : petsc_changed,
732 0 : sys);
733 :
734 : // Record whether the user changed the solution or the search direction.
735 0 : if (petsc_changed)
736 0 : *changed = PETSC_TRUE;
737 :
738 0 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
739 0 : }
740 :
741 : } // end extern "C"
742 :
743 :
744 :
745 : //---------------------------------------------------------------------
746 : // PetscNonlinearSolver<> methods
747 : template <typename T>
748 1470 : PetscNonlinearSolver<T>::PetscNonlinearSolver (sys_type & system_in) :
749 : NonlinearSolver<T>(system_in),
750 1386 : linesearch_object(nullptr),
751 1386 : _reason(SNES_CONVERGED_ITERATING/*==0*/), // Arbitrary initial value...
752 1386 : _n_linear_iterations(0),
753 1386 : _current_nonlinear_iteration_number(0),
754 1386 : _zero_out_residual(true),
755 1386 : _zero_out_jacobian(true),
756 1386 : _default_monitor(true),
757 1386 : _snesmf_reuse_base(true),
758 1386 : _computing_base_vector(true),
759 1386 : _setup_reuse(false),
760 1470 : _mffd_jac(this->_communicator)
761 : {
762 1470 : }
763 :
764 :
765 :
766 : template <typename T>
767 2772 : PetscNonlinearSolver<T>::~PetscNonlinearSolver () = default;
768 :
769 :
770 :
771 : template <typename T>
772 29680 : void PetscNonlinearSolver<T>::clear ()
773 : {
774 29680 : if (this->initialized())
775 : {
776 29400 : this->_is_initialized = false;
777 :
778 : // If we don't need the preconditioner next time
779 : // retain the original behavior of clearing the data
780 : // between solves.
781 29400 : if (!(reuse_preconditioner()))
782 : {
783 : // SNESReset really ought to work but replacing destroy() with
784 : // SNESReset causes a very slight change in behavior that
785 : // manifests as two failed MOOSE tests...
786 29400 : _snes.destroy();
787 : }
788 :
789 : // Reset the nonlinear iteration counter. This information is only relevant
790 : // *during* the solve(). After the solve is completed it should return to
791 : // the default value of 0.
792 29400 : _current_nonlinear_iteration_number = 0;
793 : }
794 29680 : }
795 :
796 : template <typename T>
797 180780 : void PetscNonlinearSolver<T>::init (const char * name)
798 : {
799 4078 : parallel_object_only();
800 :
801 : // Initialize the data structures if not done so already.
802 180780 : if (!this->initialized())
803 : {
804 29400 : this->_is_initialized = true;
805 :
806 : // Make only if we don't already have a retained snes
807 : // hanging around from the last solve
808 29400 : if (!_snes)
809 29400 : LibmeshPetscCall(SNESCreate(this->comm().get(), _snes.get()));
810 :
811 : // I believe all of the following can be safely repeated
812 : // even on an old snes instance from the last solve
813 :
814 29400 : if (name)
815 : {
816 0 : libmesh_assert(std::string(name).front() != '-');
817 0 : libmesh_assert(std::string(name).back() == '_');
818 0 : LibmeshPetscCall(SNESSetOptionsPrefix(_snes, name));
819 : }
820 :
821 : // Attaching a DM to SNES.
822 : #if defined(LIBMESH_ENABLE_AMR) && defined(LIBMESH_HAVE_METAPHYSICL)
823 57960 : bool use_petsc_dm = libMesh::on_command_line(
824 27720 : "--" + (name ? std::string(name) : std::string("")) + "use_petsc_dm");
825 :
826 : // This needs to be called before SNESSetFromOptions
827 29400 : if (use_petsc_dm)
828 0 : this->_dm_wrapper.init_and_attach_petscdm(this->system(), _snes);
829 : else
830 : #endif
831 : {
832 1680 : WrappedPetsc<DM> dm;
833 29400 : LibmeshPetscCall(DMCreate(this->comm().get(), dm.get()));
834 29400 : LibmeshPetscCall(DMSetType(dm, DMLIBMESH));
835 29400 : LibmeshPetscCall(DMlibMeshSetSystem(dm, this->system()));
836 :
837 29400 : if (name)
838 0 : LibmeshPetscCall(DMSetOptionsPrefix(dm, name));
839 :
840 29400 : LibmeshPetscCall(DMSetFromOptions(dm));
841 29400 : LibmeshPetscCall(DMSetUp(dm));
842 29400 : LibmeshPetscCall(SNESSetDM(_snes, dm));
843 : // SNES now owns the reference to dm.
844 : }
845 :
846 29400 : setup_default_monitor();
847 :
848 : // If the SolverConfiguration object is provided, use it to set
849 : // options during solver initialization.
850 29400 : if (this->_solver_configuration)
851 : {
852 0 : this->_solver_configuration->set_options_during_init();
853 : }
854 :
855 29400 : if (this->_preconditioner)
856 : {
857 : KSP ksp;
858 280 : LibmeshPetscCall(SNESGetKSP (_snes, &ksp));
859 : PC pc;
860 280 : LibmeshPetscCall(KSPGetPC(ksp,&pc));
861 :
862 280 : this->_preconditioner->init();
863 :
864 280 : LibmeshPetscCall(PCSetType(pc, PCSHELL));
865 280 : LibmeshPetscCall(PCShellSetContext(pc,(void *)this->_preconditioner));
866 :
867 : //Re-Use the shell functions from petsc_linear_solver
868 280 : LibmeshPetscCall(PCShellSetSetUp(pc,libmesh_petsc_preconditioner_setup));
869 280 : LibmeshPetscCall(PCShellSetApply(pc,libmesh_petsc_preconditioner_apply));
870 : }
871 : }
872 :
873 :
874 : // Tell PETSc about our linesearch "post-check" function, but only
875 : // if the user has provided one. There seem to be extra,
876 : // unnecessary residual calculations if a postcheck function is
877 : // attached for no reason.
878 180780 : if (this->postcheck || this->postcheck_object)
879 : {
880 : SNESLineSearch linesearch;
881 140 : LibmeshPetscCall(SNESGetLineSearch(_snes, &linesearch));
882 :
883 140 : LibmeshPetscCall(SNESLineSearchSetPostCheck(linesearch, libmesh_petsc_snes_postcheck, this));
884 : }
885 :
886 180780 : if (this->precheck_object)
887 : {
888 : SNESLineSearch linesearch;
889 0 : LibmeshPetscCall(SNESGetLineSearch(_snes, &linesearch));
890 :
891 0 : LibmeshPetscCall(SNESLineSearchSetPreCheck(linesearch, libmesh_petsc_snes_precheck, this));
892 : }
893 180780 : }
894 :
895 :
896 : template <typename T>
897 4388 : SNES PetscNonlinearSolver<T>::snes(const char * name)
898 : {
899 121980 : this->init(name);
900 4388 : return _snes;
901 : }
902 :
903 :
904 :
905 : template <typename T>
906 : void
907 0 : PetscNonlinearSolver<T>::build_mat_null_space(NonlinearImplicitSystem::ComputeVectorSubspace * computeSubspaceObject,
908 : void (*computeSubspace)(std::vector<NumericVector<Number> *> &, sys_type &),
909 : MatNullSpace * msp)
910 : {
911 0 : parallel_object_only();
912 :
913 0 : std::vector<NumericVector<Number> *> sp;
914 0 : if (computeSubspaceObject)
915 0 : (*computeSubspaceObject)(sp, this->system());
916 : else
917 0 : (*computeSubspace)(sp, this->system());
918 :
919 0 : *msp = LIBMESH_PETSC_NULLPTR;
920 0 : if (sp.size())
921 : {
922 0 : PetscInt nmodes = cast_int<PetscInt>(sp.size());
923 :
924 0 : std::vector<Vec> modes(nmodes);
925 0 : std::vector<PetscScalar> dots(nmodes);
926 :
927 0 : for (PetscInt i=0; i<nmodes; ++i)
928 : {
929 0 : auto pv = cast_ptr<PetscVector<T> *>(sp[i]);
930 :
931 0 : LibmeshPetscCall(VecDuplicate(pv->vec(), &modes[i]));
932 :
933 0 : LibmeshPetscCall(VecCopy(pv->vec(), modes[i]));
934 : }
935 :
936 : // Normalize.
937 0 : LibmeshPetscCall(VecNormalize(modes[0], LIBMESH_PETSC_NULLPTR));
938 :
939 0 : for (PetscInt i=1; i<nmodes; i++)
940 : {
941 : // Orthonormalize vec[i] against vec[0:i-1]
942 0 : LibmeshPetscCall(VecMDot(modes[i], i, modes.data(), dots.data()));
943 :
944 0 : for (PetscInt j=0; j<i; j++)
945 0 : dots[j] *= -1.;
946 :
947 0 : LibmeshPetscCall(VecMAXPY(modes[i], i, dots.data(), modes.data()));
948 :
949 0 : LibmeshPetscCall(VecNormalize(modes[i], LIBMESH_PETSC_NULLPTR));
950 : }
951 :
952 0 : LibmeshPetscCall(MatNullSpaceCreate(this->comm().get(), PETSC_FALSE, nmodes, modes.data(), msp));
953 :
954 0 : for (PetscInt i=0; i<nmodes; ++i)
955 0 : LibmeshPetscCall(VecDestroy(&modes[i]));
956 : }
957 0 : }
958 :
959 : template <typename T>
960 : std::pair<unsigned int, Real>
961 29400 : PetscNonlinearSolver<T>::solve (SparseMatrix<T> & pre_in, // System Preconditioning Matrix
962 : NumericVector<T> & x_in, // Solution vector
963 : NumericVector<T> & r_in, // Residual vector
964 : const double, // Stopping tolerance
965 : const unsigned int)
966 : {
967 840 : parallel_object_only();
968 :
969 840 : LOG_SCOPE("solve()", "PetscNonlinearSolver");
970 29400 : this->init ();
971 :
972 : // Make sure the data passed in are really of Petsc types
973 840 : PetscMatrixBase<T> * pre = cast_ptr<PetscMatrixBase<T> *>(&pre_in);
974 840 : PetscVector<T> * x = cast_ptr<PetscVector<T> *>(&x_in);
975 840 : PetscVector<T> * r = cast_ptr<PetscVector<T> *>(&r_in);
976 :
977 29400 : PetscInt n_iterations =0;
978 : // Should actually be a PetscReal, but I don't know which version of PETSc first introduced PetscReal
979 29400 : Real final_residual_norm=0.;
980 :
981 : // We don't want to do this twice because it resets
982 : // SNESSetLagPreconditioner
983 29400 : if ((reuse_preconditioner()) && (!_setup_reuse))
984 : {
985 0 : _setup_reuse = true;
986 0 : LibmeshPetscCall(SNESSetLagPreconditionerPersists(_snes, PETSC_TRUE));
987 : // According to the PETSC 3.16.5 docs -2 is a magic number which
988 : // means "recalculate the next time you need it and then not again"
989 0 : LibmeshPetscCall(SNESSetLagPreconditioner(_snes, -2));
990 : // Add in our callback which will trigger recalculating
991 : // the preconditioner when we hit reuse_preconditioner_max_linear_its
992 0 : LibmeshPetscCall(SNESMonitorSet(_snes, &libmesh_petsc_recalculate_monitor,
993 : this,
994 : NULL));
995 : }
996 29400 : else if (!(reuse_preconditioner()))
997 : // This covers the case where it was enabled but was then disabled
998 : {
999 29400 : LibmeshPetscCall(SNESSetLagPreconditionerPersists(_snes, PETSC_FALSE));
1000 29400 : if (_setup_reuse)
1001 : {
1002 0 : _setup_reuse = false;
1003 0 : LibmeshPetscCall(SNESMonitorCancel(_snes));
1004 : // Readd default monitor
1005 0 : setup_default_monitor();
1006 : }
1007 : }
1008 :
1009 29400 : LibmeshPetscCall(SNESSetFunction (_snes, r->vec(), libmesh_petsc_snes_residual, this));
1010 :
1011 : // Only set the jacobian function if we've been provided with something to call.
1012 : // This allows a user to set their own jacobian function if they want to
1013 29400 : if (this->jacobian || this->jacobian_object || this->residual_and_jacobian_object)
1014 29400 : LibmeshPetscCall(SNESSetJacobian (_snes, pre->mat(), pre->mat(), libmesh_petsc_snes_jacobian, this));
1015 :
1016 : // Have the Krylov subspace method use our good initial guess rather than 0
1017 : KSP ksp;
1018 29400 : LibmeshPetscCall(SNESGetKSP (_snes, &ksp));
1019 :
1020 : // Set the tolerances for the iterative solver. Use the user-supplied
1021 : // tolerance for the relative residual & leave the others at default values
1022 29400 : LibmeshPetscCall(KSPSetTolerances (ksp, this->initial_linear_tolerance, PETSC_DEFAULT,
1023 : PETSC_DEFAULT, this->max_linear_iterations));
1024 :
1025 : // Set the tolerances for the non-linear solver.
1026 29400 : LibmeshPetscCall(SNESSetTolerances(_snes,
1027 : this->absolute_residual_tolerance,
1028 : this->relative_residual_tolerance,
1029 : this->relative_step_tolerance,
1030 : this->max_nonlinear_iterations,
1031 : this->max_function_evaluations));
1032 :
1033 : // Set the divergence tolerance for the non-linear solver
1034 : #if !PETSC_VERSION_LESS_THAN(3,8,0)
1035 29400 : LibmeshPetscCall(SNESSetDivergenceTolerance(_snes, this->divergence_tolerance));
1036 : #endif
1037 :
1038 : //Pull in command-line options
1039 : #if PETSC_VERSION_LESS_THAN(3,7,0)
1040 : LibmeshPetscCall(KSPSetFromOptions(ksp));
1041 : #endif
1042 29400 : LibmeshPetscCall(SNESSetFromOptions(_snes));
1043 :
1044 : #if defined(LIBMESH_HAVE_PETSC_HYPRE) && PETSC_VERSION_LESS_THAN(3, 23, 0) && \
1045 : !PETSC_VERSION_LESS_THAN(3, 12, 0) && defined(PETSC_HAVE_HYPRE_DEVICE)
1046 : {
1047 : // Make sure hypre has been initialized
1048 : LibmeshPetscCallExternal(HYPRE_Initialize);
1049 : PetscScalar * dummyarray;
1050 : PetscMemType mtype;
1051 : LibmeshPetscCall(VecGetArrayAndMemType(x->vec(), &dummyarray, &mtype));
1052 : LibmeshPetscCall(VecRestoreArrayAndMemType(x->vec(), &dummyarray));
1053 : if (PetscMemTypeHost(mtype))
1054 : LibmeshPetscCallExternal(HYPRE_SetMemoryLocation, HYPRE_MEMORY_HOST);
1055 : }
1056 : #endif
1057 :
1058 29400 : if (this->user_presolve)
1059 0 : this->user_presolve(this->system());
1060 :
1061 : //Set the preconditioning matrix
1062 29400 : if (this->_preconditioner)
1063 : {
1064 8 : this->_preconditioner->set_matrix(pre_in);
1065 280 : this->_preconditioner->init();
1066 : }
1067 :
1068 : // If the SolverConfiguration object is provided, use it to override
1069 : // solver options.
1070 29400 : if (this->_solver_configuration)
1071 0 : this->_solver_configuration->configure_solver();
1072 :
1073 : // In PETSc versions before 3.5.0, it is not possible to call
1074 : // SNESSetUp() before the solution and rhs vectors are initialized, as
1075 : // this triggers the
1076 : //
1077 : // "Solution vector cannot be right hand side vector!"
1078 : //
1079 : // error message. It is also not possible to call SNESSetSolution()
1080 : // in those versions of PETSc to work around the problem, since that
1081 : // API was removed in 3.0.0 and only restored in 3.6.0. The
1082 : // overzealous check was moved out of SNESSetUp in PETSc 3.5.0
1083 : // (petsc/petsc@154060b), so this code block should be safe to use
1084 : // in 3.5.0 and later.
1085 : #if !PETSC_VERSION_LESS_THAN(3,6,0)
1086 29400 : LibmeshPetscCall(SNESSetSolution(_snes, x->vec()));
1087 : #endif
1088 29400 : LibmeshPetscCall(SNESSetUp(_snes));
1089 :
1090 : Mat J, P;
1091 29400 : LibmeshPetscCall(SNESGetJacobian(_snes, &J, &P,
1092 : LIBMESH_PETSC_NULLPTR,
1093 : LIBMESH_PETSC_NULLPTR));
1094 29400 : LibmeshPetscCall(MatMFFDSetFunction(J, libmesh_petsc_snes_mffd_interface, this));
1095 : #if !PETSC_VERSION_LESS_THAN(3,8,4)
1096 : #ifndef NDEBUG
1097 : // If we're in debug mode, do not reuse the nonlinear function evaluation as the base for doing
1098 : // matrix-free approximations of the Jacobian action. Instead if the user requested that we reuse
1099 : // the base, we'll check the base function evaluation and compare it to the nonlinear residual
1100 : // evaluation. If they are different, then we'll error and inform the user that it's unsafe to
1101 : // reuse the base
1102 840 : LibmeshPetscCall(MatSNESMFSetReuseBase(J, PETSC_FALSE));
1103 : #else
1104 : // Resue the residual vector from SNES
1105 28560 : LibmeshPetscCall(MatSNESMFSetReuseBase(J, static_cast<PetscBool>(_snesmf_reuse_base)));
1106 : #endif
1107 : #endif
1108 :
1109 : // Only set the nullspace if we have a way of computing it and the result is non-empty.
1110 29400 : if (this->nullspace || this->nullspace_object)
1111 : {
1112 0 : WrappedPetsc<MatNullSpace> msp;
1113 0 : this->build_mat_null_space(this->nullspace_object, this->nullspace, msp.get());
1114 0 : if (msp)
1115 : {
1116 0 : LibmeshPetscCall(MatSetNullSpace(J, msp));
1117 0 : if (P != J)
1118 0 : LibmeshPetscCall(MatSetNullSpace(P, msp));
1119 : }
1120 : }
1121 :
1122 : // Only set the transpose nullspace if we have a way of computing it and the result is non-empty.
1123 29400 : if (this->transpose_nullspace || this->transpose_nullspace_object)
1124 : {
1125 : #if PETSC_VERSION_LESS_THAN(3,6,0)
1126 : libmesh_warning("MatSetTransposeNullSpace is only supported for PETSc >= 3.6, transpose nullspace will be ignored.");
1127 : #else
1128 0 : WrappedPetsc<MatNullSpace> msp;
1129 0 : this->build_mat_null_space(this->transpose_nullspace_object, this->transpose_nullspace, msp.get());
1130 0 : if (msp)
1131 : {
1132 0 : LibmeshPetscCall(MatSetTransposeNullSpace(J, msp));
1133 0 : if (P != J)
1134 0 : LibmeshPetscCall(MatSetTransposeNullSpace(P, msp));
1135 : }
1136 : #endif
1137 : }
1138 :
1139 : // Only set the nearnullspace if we have a way of computing it and the result is non-empty.
1140 29400 : if (this->nearnullspace || this->nearnullspace_object)
1141 : {
1142 0 : WrappedPetsc<MatNullSpace> msp;
1143 0 : this->build_mat_null_space(this->nearnullspace_object, this->nearnullspace, msp.get());
1144 :
1145 0 : if (msp)
1146 : {
1147 0 : LibmeshPetscCall(MatSetNearNullSpace(J, msp));
1148 0 : if (P != J)
1149 0 : LibmeshPetscCall(MatSetNearNullSpace(P, msp));
1150 : }
1151 : }
1152 :
1153 : SNESLineSearch linesearch;
1154 29400 : if (linesearch_object)
1155 : {
1156 0 : LibmeshPetscCall(SNESGetLineSearch(_snes, &linesearch));
1157 0 : LibmeshPetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHSHELL));
1158 : #if PETSC_RELEASE_GREATER_EQUALS(3, 21, 0)
1159 : LibmeshPetscCall(SNESLineSearchShellSetApply(linesearch, libmesh_petsc_linesearch_shellfunc, this));
1160 : #else
1161 0 : LibmeshPetscCall(SNESLineSearchShellSetUserFunc(linesearch, libmesh_petsc_linesearch_shellfunc, this));
1162 : #endif
1163 : }
1164 :
1165 29400 : LibmeshPetscCall(SNESSolve (_snes, LIBMESH_PETSC_NULLPTR, x->vec()));
1166 :
1167 29400 : LibmeshPetscCall(SNESGetIterationNumber(_snes, &n_iterations));
1168 :
1169 29400 : LibmeshPetscCall(SNESGetLinearSolveIterations(_snes, &_n_linear_iterations));
1170 :
1171 : // SNESGetFunction has been around forever and should work on all
1172 : // versions of PETSc. This is also now the recommended approach
1173 : // according to the documentation for the PETSc 3.5.1 release:
1174 : // http://www.mcs.anl.gov/petsc/documentation/changes/35.html
1175 : Vec f;
1176 29400 : LibmeshPetscCall(SNESGetFunction(_snes, &f, 0, 0));
1177 29400 : LibmeshPetscCall(VecNorm(f, NORM_2, pPR(&final_residual_norm)));
1178 :
1179 : // Get and store the reason for convergence
1180 29400 : LibmeshPetscCall(SNESGetConvergedReason(_snes, &_reason));
1181 :
1182 : //Based on Petsc 2.3.3 documentation all diverged reasons are negative
1183 29400 : this->converged = (_reason >= 0);
1184 :
1185 : // Reset data structure
1186 29400 : this->clear();
1187 :
1188 : // return the # of its. and the final residual norm.
1189 31080 : return std::make_pair(n_iterations, final_residual_norm);
1190 : }
1191 :
1192 :
1193 :
1194 : template <typename T>
1195 280 : void PetscNonlinearSolver<T>::print_converged_reason()
1196 : {
1197 :
1198 8 : libMesh::out << "Nonlinear solver convergence/divergence reason: "
1199 280 : << SNESConvergedReasons[this->get_converged_reason()] << std::endl;
1200 280 : }
1201 :
1202 :
1203 :
1204 : template <typename T>
1205 280 : SNESConvergedReason PetscNonlinearSolver<T>::get_converged_reason()
1206 : {
1207 280 : if (this->initialized())
1208 0 : LibmeshPetscCall(SNESGetConvergedReason(_snes, &_reason));
1209 :
1210 280 : return _reason;
1211 : }
1212 :
1213 : template <typename T>
1214 0 : int PetscNonlinearSolver<T>::get_total_linear_iterations()
1215 : {
1216 0 : return _n_linear_iterations;
1217 : }
1218 :
1219 : template <typename T>
1220 29400 : void PetscNonlinearSolver<T>::setup_default_monitor()
1221 : {
1222 29400 : if (_default_monitor)
1223 29400 : LibmeshPetscCall(
1224 : SNESMonitorSet(_snes, libmesh_petsc_snes_monitor, this, LIBMESH_PETSC_NULLPTR));
1225 29400 : }
1226 :
1227 : template <typename T>
1228 88200 : bool PetscNonlinearSolver<T>::reuse_preconditioner() const
1229 : {
1230 88200 : return this->_reuse_preconditioner;
1231 : }
1232 :
1233 : template <typename T>
1234 0 : unsigned int PetscNonlinearSolver<T>::reuse_preconditioner_max_linear_its() const
1235 : {
1236 0 : return this->_reuse_preconditioner_max_linear_its;
1237 : }
1238 :
1239 : template <typename T>
1240 0 : void PetscNonlinearSolver<T>::force_new_preconditioner()
1241 : {
1242 : // Easiest way is just to clear everything out
1243 0 : this->_is_initialized = false;
1244 0 : _snes.destroy();
1245 0 : _setup_reuse = false;
1246 0 : }
1247 :
1248 : //------------------------------------------------------------------
1249 : // Explicit instantiations
1250 : template class LIBMESH_EXPORT PetscNonlinearSolver<Number>;
1251 :
1252 : } // namespace libMesh
1253 :
1254 :
1255 :
1256 : #endif // #ifdef LIBMESH_HAVE_PETSC
|
