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