Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 :
19 :
20 : #include "libmesh/libmesh_common.h"
21 :
22 : #ifdef LIBMESH_HAVE_PETSC
23 :
24 : // Local Includes
25 : #include "libmesh/libmesh_logging.h"
26 : #include "libmesh/nonlinear_implicit_system.h"
27 : #include "libmesh/petsc_nonlinear_solver.h"
28 : #include "libmesh/petsc_linear_solver.h"
29 : #include "libmesh/petsc_vector.h"
30 : #include "libmesh/petsc_mffd_matrix.h"
31 : #include "libmesh/dof_map.h"
32 : #include "libmesh/preconditioner.h"
33 : #include "libmesh/solver_configuration.h"
34 : #include "libmesh/petscdmlibmesh.h"
35 : #include "libmesh/petsc_mffd_matrix.h"
36 :
37 : #if defined(LIBMESH_HAVE_PETSC_HYPRE) && PETSC_VERSION_LESS_THAN(3, 23, 0) && \
38 : !PETSC_VERSION_LESS_THAN(3, 12, 0) && defined(PETSC_HAVE_HYPRE_DEVICE)
39 : #include <HYPRE_utilities.h>
40 : #endif
41 :
42 : namespace libMesh
43 : {
44 : class ResidualContext
45 : {
46 : public:
47 5660 : ResidualContext(PetscNonlinearSolver<Number> * solver_in, NonlinearImplicitSystem & sys_in) :
48 : solver(solver_in),
49 5660 : sys(sys_in)
50 5660 : {}
51 :
52 : PetscNonlinearSolver<Number> * solver;
53 : NonlinearImplicitSystem & sys;
54 : };
55 :
56 : ResidualContext
57 281117 : libmesh_petsc_snes_residual_helper (SNES snes, Vec x, void * ctx)
58 : {
59 11320 : LOG_SCOPE("residual()", "PetscNonlinearSolver");
60 :
61 5660 : libmesh_assert(x);
62 5660 : libmesh_assert(ctx);
63 :
64 : // No way to safety-check this cast, since we got a void *...
65 5660 : PetscNonlinearSolver<Number> * solver =
66 : static_cast<PetscNonlinearSolver<Number> *> (ctx);
67 :
68 5660 : libmesh_parallel_only(solver->comm());
69 :
70 : // Get the current iteration number from the snes object,
71 : // store it in the PetscNonlinearSolver object for possible use
72 : // by the user's residual function.
73 : {
74 281117 : PetscInt n_iterations = 0;
75 281117 : LibmeshPetscCall2(solver->comm(), SNESGetIterationNumber(snes, &n_iterations));
76 281117 : solver->_current_nonlinear_iteration_number = cast_int<unsigned>(n_iterations);
77 : }
78 :
79 10912 : NonlinearImplicitSystem & sys = solver->system();
80 :
81 5660 : PetscVector<Number> & X_sys = *cast_ptr<PetscVector<Number> *>(sys.solution.get());
82 :
83 281117 : PetscVector<Number> X_global(x, sys.comm());
84 :
85 : // Use the system's update() to get a good local version of the
86 : // parallel solution. This operation does not modify the incoming
87 : // "x" vector, it only localizes information from "x" into
88 : // sys.current_local_solution.
89 281117 : X_global.swap(X_sys);
90 281117 : sys.update();
91 281117 : X_global.swap(X_sys);
92 :
93 : // Enforce constraints (if any) exactly on the
94 : // current_local_solution. This is the solution vector that is
95 : // actually used in the computation of the residual below, and is
96 : // not locked by debug-enabled PETSc the way that "x" is.
97 281117 : if (solver->_exact_constraint_enforcement)
98 281117 : sys.get_dof_map().enforce_constraints_exactly(sys, sys.current_local_solution.get());
99 :
100 292029 : return ResidualContext(solver, sys);
101 270205 : }
102 :
103 : //--------------------------------------------------------------------
104 : // Functions with C linkage to pass to PETSc. PETSc will call these
105 : // methods as needed.
106 : //
107 : // Since they must have C linkage they have no knowledge of a namespace.
108 : // Give them an obscure name to avoid namespace pollution.
109 : extern "C"
110 : {
111 : // -----------------------------------------------------------------
112 : // this function monitors the nonlinear solve and checks to see
113 : // if we want to recalculate the preconditioner. It only gets
114 : // added to the SNES instance if we're reusing the preconditioner
115 : PetscErrorCode
116 0 : libmesh_petsc_recalculate_monitor(SNES snes, PetscInt, PetscReal, void* ctx)
117 : {
118 : PetscFunctionBegin;
119 :
120 : // No way to safety-check this cast, since we got a void *...
121 0 : PetscNonlinearSolver<Number> * solver =
122 : static_cast<PetscNonlinearSolver<Number> *> (ctx);
123 :
124 : KSP ksp;
125 0 : LibmeshPetscCall2(solver->comm(), SNESGetKSP(snes, &ksp));
126 :
127 : PetscInt niter;
128 0 : LibmeshPetscCall2(solver->comm(), KSPGetIterationNumber(ksp, &niter));
129 :
130 0 : if (niter > cast_int<PetscInt>(solver->reuse_preconditioner_max_linear_its()))
131 : {
132 : // -2 is a magic number for "recalculate next time you need it
133 : // and then not again"
134 0 : LibmeshPetscCall2(solver->comm(), SNESSetLagPreconditioner(snes, -2));
135 : }
136 0 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
137 : }
138 :
139 : //-------------------------------------------------------------------
140 : // this function is called by PETSc at the end of each nonlinear step
141 : PetscErrorCode
142 105241 : libmesh_petsc_snes_monitor (SNES, PetscInt its, PetscReal fnorm, void *)
143 : {
144 : PetscFunctionBegin;
145 2896 : libMesh::out << " NL step "
146 2896 : << std::setw(2) << its
147 2896 : << std::scientific
148 2896 : << ", |residual|_2 = " << fnorm
149 2896 : << std::endl;
150 105241 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
151 : }
152 :
153 : //---------------------------------------------------------------
154 : // this function is called by PETSc to evaluate the residual at X
155 : PetscErrorCode
156 159137 : libmesh_petsc_snes_residual (SNES snes, Vec x, Vec r, void * ctx)
157 : {
158 : PetscFunctionBegin;
159 :
160 159137 : 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 165661 : PetscVector<Number> R(r, rc.sys.comm());
166 :
167 159137 : if (rc.solver->_zero_out_residual)
168 159137 : 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 159137 : 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 159137 : 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 159137 : if (rc.solver->residual != nullptr)
180 280 : rc.solver->residual(*rc.sys.current_local_solution.get(), R, rc.sys);
181 :
182 158857 : else if (rc.solver->residual_object != nullptr)
183 3277 : 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 162399 : PetscVector<Number> X_global(x, rc.sys.comm());
213 :
214 159137 : X_global.swap(X_sys);
215 159137 : rc.sys.update();
216 159137 : X_global.swap(X_sys);
217 :
218 159137 : R.close();
219 :
220 159137 : if (rc.solver->_exact_constraint_enforcement)
221 : {
222 159137 : rc.sys.get_dof_map().enforce_constraints_on_residual(rc.sys, &R, rc.sys.current_local_solution.get());
223 159137 : R.close();
224 : }
225 :
226 162399 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
227 152613 : }
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 76064 : 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 76064 : PetscInt n_iterations = 0;
420 76064 : LibmeshPetscCall2(solver->comm(), SNESGetIterationNumber(snes, &n_iterations));
421 76064 : 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 76064 : libmesh_error_msg_if(solver->jacobian && solver->jacobian_object,
428 : "ERROR: cannot specify both a function and object to compute the Jacobian!");
429 :
430 76064 : 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 76064 : PetscMatrixBase<Number> * const PC = pc ? PetscMatrixBase<Number>::get_context(pc, sys.comm()) : nullptr;
436 76064 : 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 80184 : PetscVector<Number> X_global(x, sys.comm());
439 :
440 6180 : PetscMFFDMatrix<Number> mffd_jac(sys.comm());
441 76064 : PetscBool p_is_shell = PETSC_FALSE;
442 76064 : PetscBool j_is_mffd = PETSC_FALSE;
443 76064 : PetscBool j_is_shell = PETSC_FALSE;
444 76064 : if (pc)
445 76064 : LibmeshPetscCall2(sys.comm(), PetscObjectTypeCompare((PetscObject)pc, MATSHELL, &p_is_shell));
446 2060 : libmesh_assert(jac);
447 76064 : LibmeshPetscCall2(sys.comm(), PetscObjectTypeCompare((PetscObject)jac, MATMFFD, &j_is_mffd));
448 76064 : LibmeshPetscCall2(sys.comm(), PetscObjectTypeCompare((PetscObject)jac, MATSHELL, &j_is_shell));
449 76064 : if (j_is_mffd == PETSC_TRUE)
450 : {
451 408 : libmesh_assert(!Jac);
452 408 : Jac = &mffd_jac;
453 408 : mffd_jac = jac;
454 : }
455 :
456 : // We already computed the Jacobian during the residual evaluation
457 76064 : if (solver->residual_and_jacobian_object)
458 : {
459 : // We could be doing matrix-free in which case we cannot rely on closing of explicit matrices
460 : // that occurs during the PETSc residual callback
461 73312 : if ((j_is_shell == PETSC_TRUE) || (j_is_mffd == PETSC_TRUE))
462 14154 : Jac->close();
463 :
464 73312 : if (pc && (p_is_shell == PETSC_TRUE))
465 0 : PC->close();
466 :
467 73312 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
468 : }
469 :
470 : // Set the dof maps
471 2836 : PC->attach_dof_map(sys.get_dof_map());
472 2836 : Jac->attach_dof_map(sys.get_dof_map());
473 :
474 : // Use the systems update() to get a good local version of the parallel solution
475 2752 : X_global.swap(X_sys);
476 2752 : sys.update();
477 2752 : X_global.swap(X_sys);
478 :
479 : // Enforce constraints (if any) exactly on the
480 : // current_local_solution. This is the solution vector that is
481 : // actually used in the computation of the residual below, and is
482 : // not locked by debug-enabled PETSc the way that "x" is.
483 2752 : if (solver->_exact_constraint_enforcement)
484 2752 : sys.get_dof_map().enforce_constraints_exactly(sys, sys.current_local_solution.get());
485 :
486 2752 : if (solver->_zero_out_jacobian)
487 2752 : PC->zero();
488 :
489 :
490 2752 : if (solver->jacobian != nullptr)
491 210 : solver->jacobian(*sys.current_local_solution.get(), *PC, sys);
492 :
493 2542 : else if (solver->jacobian_object != nullptr)
494 2620 : solver->jacobian_object->jacobian(*sys.current_local_solution.get(), *PC, sys);
495 :
496 0 : else if (solver->matvec != nullptr)
497 0 : solver->matvec(*sys.current_local_solution.get(), nullptr, PC, sys);
498 :
499 : else
500 0 : libmesh_error_msg("Error! Unable to compute residual and/or Jacobian!");
501 :
502 2752 : PC->close();
503 2752 : if (solver->_exact_constraint_enforcement)
504 : {
505 2752 : sys.get_dof_map().enforce_constraints_on_jacobian(sys, PC);
506 2752 : PC->close();
507 : }
508 :
509 2752 : if (Jac != PC)
510 : {
511 : // Assume that shells know what they're doing
512 0 : libmesh_assert(!solver->_exact_constraint_enforcement || (j_is_mffd == PETSC_TRUE) ||
513 : (j_is_shell == PETSC_TRUE));
514 0 : Jac->close();
515 : }
516 :
517 84 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
518 71944 : }
519 :
520 : // This function gets called by PETSc in place of the standard Petsc line searches
521 : // if a linesearch object is supplied to the PetscNonlinearSolver class. It wraps
522 : // the linesearch algorithm implemented on the linesearch object.
523 : // * "linesearch" is an object that can be used to access the non-linear and linear solution
524 : // vectors as well as the residual and SNES object
525 : // * "ctx" is the PetscNonlinearSolver context
526 0 : PetscErrorCode libmesh_petsc_linesearch_shellfunc (SNESLineSearch linesearch, void * ctx)
527 : {
528 : PetscFunctionBegin;
529 :
530 : // No way to safety-check this cast, since we got a void *...
531 0 : PetscNonlinearSolver<Number> * solver =
532 : static_cast<PetscNonlinearSolver<Number> *> (ctx);
533 :
534 0 : libmesh_parallel_only(solver->comm());
535 :
536 0 : solver->linesearch_object->linesearch(linesearch);
537 0 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
538 : }
539 :
540 : // This function gets called by PETSc after the SNES linesearch is
541 : // complete. We use it to exactly enforce any constraints on the
542 : // solution which may have drifted during the linear solve. In the
543 : // PETSc nomenclature:
544 : // * "x" is the old solution vector,
545 : // * "y" is the search direction (Newton step) vector,
546 : // * "w" is the candidate solution vector, and
547 : // the user is responsible for setting changed_y and changed_w
548 : // appropriately, depending on whether or not the search
549 : // direction or solution vector was changed, respectively.
550 210 : PetscErrorCode libmesh_petsc_snes_postcheck(SNESLineSearch, Vec x, Vec y, Vec w, PetscBool * changed_y, PetscBool * changed_w, void * context)
551 : {
552 : PetscFunctionBegin;
553 :
554 12 : LOG_SCOPE("postcheck()", "PetscNonlinearSolver");
555 :
556 : // PETSc almost certainly initializes these to false already, but
557 : // it doesn't hurt to be explicit.
558 210 : *changed_w = PETSC_FALSE;
559 210 : *changed_y = PETSC_FALSE;
560 :
561 6 : libmesh_assert(context);
562 :
563 : // Cast the context to a NonlinearSolver object.
564 6 : PetscNonlinearSolver<Number> * solver =
565 : static_cast<PetscNonlinearSolver<Number> *> (context);
566 :
567 6 : libmesh_parallel_only(solver->comm());
568 :
569 : // If the user has provided both postcheck function pointer and
570 : // object, this is ambiguous, so throw an error.
571 210 : libmesh_error_msg_if(solver->postcheck && solver->postcheck_object,
572 : "ERROR: cannot specify both a function and object for performing the solve postcheck!");
573 :
574 : // It's also possible that we don't need to do anything at all, in
575 : // that case return early...
576 12 : NonlinearImplicitSystem & sys = solver->system();
577 :
578 210 : if (!solver->postcheck && !solver->postcheck_object)
579 0 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
580 :
581 : // We definitely need to wrap at least "w"
582 216 : PetscVector<Number> petsc_w(w, sys.comm());
583 :
584 : // The user sets these flags in his/her postcheck function to
585 : // indicate whether they changed something.
586 : bool
587 210 : changed_search_direction = false,
588 210 : changed_new_soln = false;
589 :
590 210 : if (solver->postcheck || solver->postcheck_object)
591 : {
592 222 : PetscVector<Number> petsc_x(x, sys.comm());
593 222 : PetscVector<Number> petsc_y(y, sys.comm());
594 :
595 210 : if (solver->postcheck)
596 0 : solver->postcheck(petsc_x,
597 : petsc_y,
598 : petsc_w,
599 : changed_search_direction,
600 : changed_new_soln,
601 : sys);
602 :
603 210 : else if (solver->postcheck_object)
604 216 : solver->postcheck_object->postcheck(petsc_x,
605 : petsc_y,
606 : petsc_w,
607 : changed_search_direction,
608 : changed_new_soln,
609 12 : sys);
610 198 : }
611 :
612 : // Record whether the user changed the solution or the search direction.
613 210 : if (changed_search_direction)
614 0 : *changed_y = PETSC_TRUE;
615 :
616 210 : if (changed_new_soln)
617 0 : *changed_w = PETSC_TRUE;
618 :
619 6 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
620 198 : }
621 :
622 0 : PetscErrorCode libmesh_petsc_snes_precheck(SNESLineSearch, Vec X, Vec Y, PetscBool * changed, void * context)
623 : {
624 : PetscFunctionBegin;
625 :
626 0 : LOG_SCOPE("precheck()", "PetscNonlinearSolver");
627 :
628 : // PETSc almost certainly initializes these to false already, but
629 : // it doesn't hurt to be explicit.
630 0 : *changed = PETSC_FALSE;
631 :
632 0 : libmesh_assert(context);
633 :
634 : // Cast the context to a NonlinearSolver object.
635 0 : PetscNonlinearSolver<Number> * solver =
636 : static_cast<PetscNonlinearSolver<Number> *> (context);
637 :
638 0 : libmesh_parallel_only(solver->comm());
639 :
640 : // It's possible that we don't need to do anything at all, in
641 : // that case return early...
642 0 : if (!solver->precheck_object)
643 0 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
644 :
645 : // The user sets these flags in his/her postcheck function to
646 : // indicate whether they changed something.
647 : bool
648 0 : petsc_changed = false;
649 :
650 0 : auto & sys = solver->system();
651 0 : auto & x_sys = *cast_ptr<PetscVector<Number> *>(sys.solution.get());
652 0 : PetscVector<Number> petsc_x(X, sys.comm());
653 0 : PetscVector<Number> petsc_y(Y, sys.comm());
654 :
655 : // Use the systems update() to get a good local version of the parallel solution
656 0 : petsc_x.swap(x_sys);
657 0 : sys.update();
658 0 : petsc_x.swap(x_sys);
659 :
660 : // Enforce constraints (if any) exactly on the
661 : // current_local_solution. This is the solution vector that is
662 : // actually used in the computation of residuals and Jacobians, and is
663 : // not locked by debug-enabled PETSc the way that "x" is.
664 0 : libmesh_assert(sys.current_local_solution.get());
665 0 : auto & local_soln = *sys.current_local_solution.get();
666 0 : if (solver->_exact_constraint_enforcement)
667 0 : sys.get_dof_map().enforce_constraints_exactly(sys, &local_soln);
668 :
669 0 : solver->precheck_object->precheck(local_soln,
670 : petsc_y,
671 : petsc_changed,
672 0 : sys);
673 :
674 : // Record whether the user changed the solution or the search direction.
675 0 : if (petsc_changed)
676 0 : *changed = PETSC_TRUE;
677 :
678 0 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
679 0 : }
680 :
681 : } // end extern "C"
682 :
683 :
684 :
685 : //---------------------------------------------------------------------
686 : // PetscNonlinearSolver<> methods
687 : template <typename T>
688 1470 : PetscNonlinearSolver<T>::PetscNonlinearSolver (sys_type & system_in) :
689 : NonlinearSolver<T>(system_in),
690 1386 : linesearch_object(nullptr),
691 1386 : _reason(SNES_CONVERGED_ITERATING/*==0*/), // Arbitrary initial value...
692 1386 : _n_linear_iterations(0),
693 1386 : _current_nonlinear_iteration_number(0),
694 1386 : _zero_out_residual(true),
695 1386 : _zero_out_jacobian(true),
696 1386 : _default_monitor(true),
697 1386 : _snesmf_reuse_base(true),
698 1386 : _computing_base_vector(true),
699 1386 : _setup_reuse(false),
700 1470 : _mffd_jac(this->_communicator)
701 : {
702 1470 : }
703 :
704 :
705 :
706 : template <typename T>
707 2772 : PetscNonlinearSolver<T>::~PetscNonlinearSolver () = default;
708 :
709 :
710 :
711 : template <typename T>
712 29680 : void PetscNonlinearSolver<T>::clear ()
713 : {
714 29680 : if (this->initialized())
715 : {
716 29400 : this->_is_initialized = false;
717 :
718 : // If we don't need the preconditioner next time
719 : // retain the original behavior of clearing the data
720 : // between solves.
721 29400 : if (!(reuse_preconditioner()))
722 : {
723 : // SNESReset really ought to work but replacing destroy() with
724 : // SNESReset causes a very slight change in behavior that
725 : // manifests as two failed MOOSE tests...
726 29400 : _snes.destroy();
727 : }
728 :
729 : // Reset the nonlinear iteration counter. This information is only relevant
730 : // *during* the solve(). After the solve is completed it should return to
731 : // the default value of 0.
732 29400 : _current_nonlinear_iteration_number = 0;
733 : }
734 29680 : }
735 :
736 : template <typename T>
737 180780 : void PetscNonlinearSolver<T>::init (const char * name)
738 : {
739 4078 : parallel_object_only();
740 :
741 : // Initialize the data structures if not done so already.
742 180780 : if (!this->initialized())
743 : {
744 29400 : this->_is_initialized = true;
745 :
746 : // Make only if we don't already have a retained snes
747 : // hanging around from the last solve
748 29400 : if (!_snes)
749 29400 : LibmeshPetscCall(SNESCreate(this->comm().get(), _snes.get()));
750 :
751 : // I believe all of the following can be safely repeated
752 : // even on an old snes instance from the last solve
753 :
754 29400 : if (name)
755 : {
756 0 : libmesh_assert(std::string(name).front() != '-');
757 0 : libmesh_assert(std::string(name).back() == '_');
758 0 : LibmeshPetscCall(SNESSetOptionsPrefix(_snes, name));
759 : }
760 :
761 : // Attaching a DM to SNES.
762 : #if defined(LIBMESH_ENABLE_AMR) && defined(LIBMESH_HAVE_METAPHYSICL)
763 57960 : bool use_petsc_dm = libMesh::on_command_line(
764 27720 : "--" + (name ? std::string(name) : std::string("")) + "use_petsc_dm");
765 :
766 : // This needs to be called before SNESSetFromOptions
767 29400 : if (use_petsc_dm)
768 0 : this->_dm_wrapper.init_and_attach_petscdm(this->system(), _snes);
769 : else
770 : #endif
771 : {
772 1680 : WrappedPetsc<DM> dm;
773 29400 : LibmeshPetscCall(DMCreate(this->comm().get(), dm.get()));
774 29400 : LibmeshPetscCall(DMSetType(dm, DMLIBMESH));
775 29400 : LibmeshPetscCall(DMlibMeshSetSystem(dm, this->system()));
776 :
777 29400 : if (name)
778 0 : LibmeshPetscCall(DMSetOptionsPrefix(dm, name));
779 :
780 29400 : LibmeshPetscCall(DMSetFromOptions(dm));
781 29400 : LibmeshPetscCall(DMSetUp(dm));
782 29400 : LibmeshPetscCall(SNESSetDM(_snes, dm));
783 : // SNES now owns the reference to dm.
784 : }
785 :
786 29400 : setup_default_monitor();
787 :
788 : // If the SolverConfiguration object is provided, use it to set
789 : // options during solver initialization.
790 29400 : if (this->_solver_configuration)
791 : {
792 0 : this->_solver_configuration->set_options_during_init();
793 : }
794 :
795 29400 : if (this->_preconditioner)
796 : {
797 : KSP ksp;
798 280 : LibmeshPetscCall(SNESGetKSP (_snes, &ksp));
799 : PC pc;
800 280 : LibmeshPetscCall(KSPGetPC(ksp,&pc));
801 :
802 280 : this->_preconditioner->init();
803 :
804 280 : LibmeshPetscCall(PCSetType(pc, PCSHELL));
805 280 : LibmeshPetscCall(PCShellSetContext(pc,(void *)this->_preconditioner));
806 :
807 : //Re-Use the shell functions from petsc_linear_solver
808 280 : LibmeshPetscCall(PCShellSetSetUp(pc,libmesh_petsc_preconditioner_setup));
809 280 : LibmeshPetscCall(PCShellSetApply(pc,libmesh_petsc_preconditioner_apply));
810 : }
811 : }
812 :
813 :
814 : // Tell PETSc about our linesearch "post-check" function, but only
815 : // if the user has provided one. There seem to be extra,
816 : // unnecessary residual calculations if a postcheck function is
817 : // attached for no reason.
818 180780 : if (this->postcheck || this->postcheck_object)
819 : {
820 : SNESLineSearch linesearch;
821 140 : LibmeshPetscCall(SNESGetLineSearch(_snes, &linesearch));
822 :
823 140 : LibmeshPetscCall(SNESLineSearchSetPostCheck(linesearch, libmesh_petsc_snes_postcheck, this));
824 : }
825 :
826 180780 : if (this->precheck_object)
827 : {
828 : SNESLineSearch linesearch;
829 0 : LibmeshPetscCall(SNESGetLineSearch(_snes, &linesearch));
830 :
831 0 : LibmeshPetscCall(SNESLineSearchSetPreCheck(linesearch, libmesh_petsc_snes_precheck, this));
832 : }
833 180780 : }
834 :
835 :
836 : template <typename T>
837 4388 : SNES PetscNonlinearSolver<T>::snes(const char * name)
838 : {
839 121980 : this->init(name);
840 4388 : return _snes;
841 : }
842 :
843 :
844 :
845 : template <typename T>
846 : void
847 0 : PetscNonlinearSolver<T>::build_mat_null_space(NonlinearImplicitSystem::ComputeVectorSubspace * computeSubspaceObject,
848 : void (*computeSubspace)(std::vector<NumericVector<Number> *> &, sys_type &),
849 : MatNullSpace * msp)
850 : {
851 0 : parallel_object_only();
852 :
853 0 : std::vector<NumericVector<Number> *> sp;
854 0 : if (computeSubspaceObject)
855 0 : (*computeSubspaceObject)(sp, this->system());
856 : else
857 0 : (*computeSubspace)(sp, this->system());
858 :
859 0 : *msp = LIBMESH_PETSC_NULLPTR;
860 0 : if (sp.size())
861 : {
862 0 : PetscInt nmodes = cast_int<PetscInt>(sp.size());
863 :
864 0 : std::vector<Vec> modes(nmodes);
865 0 : std::vector<PetscScalar> dots(nmodes);
866 :
867 0 : for (PetscInt i=0; i<nmodes; ++i)
868 : {
869 0 : auto pv = cast_ptr<PetscVector<T> *>(sp[i]);
870 :
871 0 : LibmeshPetscCall(VecDuplicate(pv->vec(), &modes[i]));
872 :
873 0 : LibmeshPetscCall(VecCopy(pv->vec(), modes[i]));
874 : }
875 :
876 : // Normalize.
877 0 : LibmeshPetscCall(VecNormalize(modes[0], LIBMESH_PETSC_NULLPTR));
878 :
879 0 : for (PetscInt i=1; i<nmodes; i++)
880 : {
881 : // Orthonormalize vec[i] against vec[0:i-1]
882 0 : LibmeshPetscCall(VecMDot(modes[i], i, modes.data(), dots.data()));
883 :
884 0 : for (PetscInt j=0; j<i; j++)
885 0 : dots[j] *= -1.;
886 :
887 0 : LibmeshPetscCall(VecMAXPY(modes[i], i, dots.data(), modes.data()));
888 :
889 0 : LibmeshPetscCall(VecNormalize(modes[i], LIBMESH_PETSC_NULLPTR));
890 : }
891 :
892 0 : LibmeshPetscCall(MatNullSpaceCreate(this->comm().get(), PETSC_FALSE, nmodes, modes.data(), msp));
893 :
894 0 : for (PetscInt i=0; i<nmodes; ++i)
895 0 : LibmeshPetscCall(VecDestroy(&modes[i]));
896 : }
897 0 : }
898 :
899 : template <typename T>
900 : std::pair<unsigned int, Real>
901 29400 : PetscNonlinearSolver<T>::solve (SparseMatrix<T> & pre_in, // System Preconditioning Matrix
902 : NumericVector<T> & x_in, // Solution vector
903 : NumericVector<T> & r_in, // Residual vector
904 : const double, // Stopping tolerance
905 : const unsigned int)
906 : {
907 840 : parallel_object_only();
908 :
909 840 : LOG_SCOPE("solve()", "PetscNonlinearSolver");
910 29400 : this->init ();
911 :
912 : // Make sure the data passed in are really of Petsc types
913 840 : PetscMatrixBase<T> * pre = cast_ptr<PetscMatrixBase<T> *>(&pre_in);
914 840 : PetscVector<T> * x = cast_ptr<PetscVector<T> *>(&x_in);
915 840 : PetscVector<T> * r = cast_ptr<PetscVector<T> *>(&r_in);
916 :
917 29400 : PetscInt n_iterations =0;
918 : // Should actually be a PetscReal, but I don't know which version of PETSc first introduced PetscReal
919 29400 : Real final_residual_norm=0.;
920 :
921 : // We don't want to do this twice because it resets
922 : // SNESSetLagPreconditioner
923 29400 : if ((reuse_preconditioner()) && (!_setup_reuse))
924 : {
925 0 : _setup_reuse = true;
926 0 : LibmeshPetscCall(SNESSetLagPreconditionerPersists(_snes, PETSC_TRUE));
927 : // According to the PETSC 3.16.5 docs -2 is a magic number which
928 : // means "recalculate the next time you need it and then not again"
929 0 : LibmeshPetscCall(SNESSetLagPreconditioner(_snes, -2));
930 : // Add in our callback which will trigger recalculating
931 : // the preconditioner when we hit reuse_preconditioner_max_linear_its
932 0 : LibmeshPetscCall(SNESMonitorSet(_snes, &libmesh_petsc_recalculate_monitor,
933 : this,
934 : NULL));
935 : }
936 29400 : else if (!(reuse_preconditioner()))
937 : // This covers the case where it was enabled but was then disabled
938 : {
939 29400 : LibmeshPetscCall(SNESSetLagPreconditionerPersists(_snes, PETSC_FALSE));
940 29400 : if (_setup_reuse)
941 : {
942 0 : _setup_reuse = false;
943 0 : LibmeshPetscCall(SNESMonitorCancel(_snes));
944 : // Readd default monitor
945 0 : setup_default_monitor();
946 : }
947 : }
948 :
949 29400 : LibmeshPetscCall(SNESSetFunction (_snes, r->vec(), libmesh_petsc_snes_residual, this));
950 :
951 : // Only set the jacobian function if we've been provided with something to call.
952 : // This allows a user to set their own jacobian function if they want to
953 29400 : if (this->jacobian || this->jacobian_object || this->residual_and_jacobian_object)
954 29400 : LibmeshPetscCall(SNESSetJacobian (_snes, pre->mat(), pre->mat(), libmesh_petsc_snes_jacobian, this));
955 :
956 : // Have the Krylov subspace method use our good initial guess rather than 0
957 : KSP ksp;
958 29400 : LibmeshPetscCall(SNESGetKSP (_snes, &ksp));
959 :
960 : // Set the tolerances for the iterative solver. Use the user-supplied
961 : // tolerance for the relative residual & leave the others at default values
962 29400 : LibmeshPetscCall(KSPSetTolerances (ksp, this->initial_linear_tolerance, PETSC_DEFAULT,
963 : PETSC_DEFAULT, this->max_linear_iterations));
964 :
965 : // Set the tolerances for the non-linear solver.
966 29400 : LibmeshPetscCall(SNESSetTolerances(_snes,
967 : this->absolute_residual_tolerance,
968 : this->relative_residual_tolerance,
969 : this->relative_step_tolerance,
970 : this->max_nonlinear_iterations,
971 : this->max_function_evaluations));
972 :
973 : // Set the divergence tolerance for the non-linear solver
974 : #if !PETSC_VERSION_LESS_THAN(3,8,0)
975 29400 : LibmeshPetscCall(SNESSetDivergenceTolerance(_snes, this->divergence_tolerance));
976 : #endif
977 :
978 : //Pull in command-line options
979 : #if PETSC_VERSION_LESS_THAN(3,7,0)
980 : LibmeshPetscCall(KSPSetFromOptions(ksp));
981 : #endif
982 29400 : LibmeshPetscCall(SNESSetFromOptions(_snes));
983 :
984 : #if defined(LIBMESH_HAVE_PETSC_HYPRE) && PETSC_VERSION_LESS_THAN(3, 23, 0) && \
985 : !PETSC_VERSION_LESS_THAN(3, 12, 0) && defined(PETSC_HAVE_HYPRE_DEVICE)
986 : {
987 : // Make sure hypre has been initialized
988 : LibmeshPetscCallExternal(HYPRE_Initialize);
989 : PetscScalar * dummyarray;
990 : PetscMemType mtype;
991 : LibmeshPetscCall(VecGetArrayAndMemType(x->vec(), &dummyarray, &mtype));
992 : LibmeshPetscCall(VecRestoreArrayAndMemType(x->vec(), &dummyarray));
993 : if (PetscMemTypeHost(mtype))
994 : LibmeshPetscCallExternal(HYPRE_SetMemoryLocation, HYPRE_MEMORY_HOST);
995 : }
996 : #endif
997 :
998 29400 : if (this->user_presolve)
999 0 : this->user_presolve(this->system());
1000 :
1001 : //Set the preconditioning matrix
1002 29400 : if (this->_preconditioner)
1003 : {
1004 8 : this->_preconditioner->set_matrix(pre_in);
1005 280 : this->_preconditioner->init();
1006 : }
1007 :
1008 : // If the SolverConfiguration object is provided, use it to override
1009 : // solver options.
1010 29400 : if (this->_solver_configuration)
1011 0 : this->_solver_configuration->configure_solver();
1012 :
1013 : // In PETSc versions before 3.5.0, it is not possible to call
1014 : // SNESSetUp() before the solution and rhs vectors are initialized, as
1015 : // this triggers the
1016 : //
1017 : // "Solution vector cannot be right hand side vector!"
1018 : //
1019 : // error message. It is also not possible to call SNESSetSolution()
1020 : // in those versions of PETSc to work around the problem, since that
1021 : // API was removed in 3.0.0 and only restored in 3.6.0. The
1022 : // overzealous check was moved out of SNESSetUp in PETSc 3.5.0
1023 : // (petsc/petsc@154060b), so this code block should be safe to use
1024 : // in 3.5.0 and later.
1025 : #if !PETSC_VERSION_LESS_THAN(3,6,0)
1026 29400 : LibmeshPetscCall(SNESSetSolution(_snes, x->vec()));
1027 : #endif
1028 29400 : LibmeshPetscCall(SNESSetUp(_snes));
1029 :
1030 : Mat J, P;
1031 29400 : LibmeshPetscCall(SNESGetJacobian(_snes, &J, &P,
1032 : LIBMESH_PETSC_NULLPTR,
1033 : LIBMESH_PETSC_NULLPTR));
1034 29400 : LibmeshPetscCall(MatMFFDSetFunction(J, libmesh_petsc_snes_mffd_interface, this));
1035 : #if !PETSC_VERSION_LESS_THAN(3,8,4)
1036 : #ifndef NDEBUG
1037 : // If we're in debug mode, do not reuse the nonlinear function evaluation as the base for doing
1038 : // matrix-free approximations of the Jacobian action. Instead if the user requested that we reuse
1039 : // the base, we'll check the base function evaluation and compare it to the nonlinear residual
1040 : // evaluation. If they are different, then we'll error and inform the user that it's unsafe to
1041 : // reuse the base
1042 840 : LibmeshPetscCall(MatSNESMFSetReuseBase(J, PETSC_FALSE));
1043 : #else
1044 : // Resue the residual vector from SNES
1045 28560 : LibmeshPetscCall(MatSNESMFSetReuseBase(J, static_cast<PetscBool>(_snesmf_reuse_base)));
1046 : #endif
1047 : #endif
1048 :
1049 : // Only set the nullspace if we have a way of computing it and the result is non-empty.
1050 29400 : if (this->nullspace || this->nullspace_object)
1051 : {
1052 0 : WrappedPetsc<MatNullSpace> msp;
1053 0 : this->build_mat_null_space(this->nullspace_object, this->nullspace, msp.get());
1054 0 : if (msp)
1055 : {
1056 0 : LibmeshPetscCall(MatSetNullSpace(J, msp));
1057 0 : if (P != J)
1058 0 : LibmeshPetscCall(MatSetNullSpace(P, msp));
1059 : }
1060 : }
1061 :
1062 : // Only set the transpose nullspace if we have a way of computing it and the result is non-empty.
1063 29400 : if (this->transpose_nullspace || this->transpose_nullspace_object)
1064 : {
1065 : #if PETSC_VERSION_LESS_THAN(3,6,0)
1066 : libmesh_warning("MatSetTransposeNullSpace is only supported for PETSc >= 3.6, transpose nullspace will be ignored.");
1067 : #else
1068 0 : WrappedPetsc<MatNullSpace> msp;
1069 0 : this->build_mat_null_space(this->transpose_nullspace_object, this->transpose_nullspace, msp.get());
1070 0 : if (msp)
1071 : {
1072 0 : LibmeshPetscCall(MatSetTransposeNullSpace(J, msp));
1073 0 : if (P != J)
1074 0 : LibmeshPetscCall(MatSetTransposeNullSpace(P, msp));
1075 : }
1076 : #endif
1077 : }
1078 :
1079 : // Only set the nearnullspace if we have a way of computing it and the result is non-empty.
1080 29400 : if (this->nearnullspace || this->nearnullspace_object)
1081 : {
1082 0 : WrappedPetsc<MatNullSpace> msp;
1083 0 : this->build_mat_null_space(this->nearnullspace_object, this->nearnullspace, msp.get());
1084 :
1085 0 : if (msp)
1086 : {
1087 0 : LibmeshPetscCall(MatSetNearNullSpace(J, msp));
1088 0 : if (P != J)
1089 0 : LibmeshPetscCall(MatSetNearNullSpace(P, msp));
1090 : }
1091 : }
1092 :
1093 : SNESLineSearch linesearch;
1094 29400 : if (linesearch_object)
1095 : {
1096 0 : LibmeshPetscCall(SNESGetLineSearch(_snes, &linesearch));
1097 0 : LibmeshPetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHSHELL));
1098 : #if PETSC_RELEASE_GREATER_EQUALS(3, 21, 0)
1099 : LibmeshPetscCall(SNESLineSearchShellSetApply(linesearch, libmesh_petsc_linesearch_shellfunc, this));
1100 : #else
1101 0 : LibmeshPetscCall(SNESLineSearchShellSetUserFunc(linesearch, libmesh_petsc_linesearch_shellfunc, this));
1102 : #endif
1103 : }
1104 :
1105 29400 : LibmeshPetscCall(SNESSolve (_snes, LIBMESH_PETSC_NULLPTR, x->vec()));
1106 :
1107 29400 : LibmeshPetscCall(SNESGetIterationNumber(_snes, &n_iterations));
1108 :
1109 29400 : LibmeshPetscCall(SNESGetLinearSolveIterations(_snes, &_n_linear_iterations));
1110 :
1111 : // SNESGetFunction has been around forever and should work on all
1112 : // versions of PETSc. This is also now the recommended approach
1113 : // according to the documentation for the PETSc 3.5.1 release:
1114 : // http://www.mcs.anl.gov/petsc/documentation/changes/35.html
1115 : Vec f;
1116 29400 : LibmeshPetscCall(SNESGetFunction(_snes, &f, 0, 0));
1117 29400 : LibmeshPetscCall(VecNorm(f, NORM_2, pPR(&final_residual_norm)));
1118 :
1119 : // Get and store the reason for convergence
1120 29400 : LibmeshPetscCall(SNESGetConvergedReason(_snes, &_reason));
1121 :
1122 : //Based on Petsc 2.3.3 documentation all diverged reasons are negative
1123 29400 : this->converged = (_reason >= 0);
1124 :
1125 : // Reset data structure
1126 29400 : this->clear();
1127 :
1128 : // return the # of its. and the final residual norm.
1129 31080 : return std::make_pair(n_iterations, final_residual_norm);
1130 : }
1131 :
1132 :
1133 :
1134 : template <typename T>
1135 280 : void PetscNonlinearSolver<T>::print_converged_reason()
1136 : {
1137 :
1138 8 : libMesh::out << "Nonlinear solver convergence/divergence reason: "
1139 280 : << SNESConvergedReasons[this->get_converged_reason()] << std::endl;
1140 280 : }
1141 :
1142 :
1143 :
1144 : template <typename T>
1145 280 : SNESConvergedReason PetscNonlinearSolver<T>::get_converged_reason()
1146 : {
1147 280 : if (this->initialized())
1148 0 : LibmeshPetscCall(SNESGetConvergedReason(_snes, &_reason));
1149 :
1150 280 : return _reason;
1151 : }
1152 :
1153 : template <typename T>
1154 0 : int PetscNonlinearSolver<T>::get_total_linear_iterations()
1155 : {
1156 0 : return _n_linear_iterations;
1157 : }
1158 :
1159 : template <typename T>
1160 29400 : void PetscNonlinearSolver<T>::setup_default_monitor()
1161 : {
1162 29400 : if (_default_monitor)
1163 29400 : LibmeshPetscCall(
1164 : SNESMonitorSet(_snes, libmesh_petsc_snes_monitor, this, LIBMESH_PETSC_NULLPTR));
1165 29400 : }
1166 :
1167 : template <typename T>
1168 88200 : bool PetscNonlinearSolver<T>::reuse_preconditioner() const
1169 : {
1170 88200 : return this->_reuse_preconditioner;
1171 : }
1172 :
1173 : template <typename T>
1174 0 : unsigned int PetscNonlinearSolver<T>::reuse_preconditioner_max_linear_its() const
1175 : {
1176 0 : return this->_reuse_preconditioner_max_linear_its;
1177 : }
1178 :
1179 : template <typename T>
1180 0 : void PetscNonlinearSolver<T>::force_new_preconditioner()
1181 : {
1182 : // Easiest way is just to clear everything out
1183 0 : this->_is_initialized = false;
1184 0 : _snes.destroy();
1185 0 : _setup_reuse = false;
1186 0 : }
1187 :
1188 : //------------------------------------------------------------------
1189 : // Explicit instantiations
1190 : template class LIBMESH_EXPORT PetscNonlinearSolver<Number>;
1191 :
1192 : } // namespace libMesh
1193 :
1194 :
1195 :
1196 : #endif // #ifdef LIBMESH_HAVE_PETSC
|
