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