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