petsc_nonlinear_solver.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 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 
25 // C++ includes
26 
27 // Local Includes
31 #include "libmesh/petsc_vector.h"
32 #include "libmesh/petsc_matrix.h"
33 #include "libmesh/dof_map.h"
34 #include "libmesh/preconditioner.h"
36 
37 /* DMlibMesh include. */
38 #include "libmesh/petscdmlibmesh.h"
39 
40 namespace libMesh
41 {
42 
43 //--------------------------------------------------------------------
44 // Functions with C linkage to pass to PETSc. PETSc will call these
45 // methods as needed.
46 //
47 // Since they must have C linkage they have no knowledge of a namespace.
48 // Give them an obscure name to avoid namespace pollution.
49 extern "C"
50 {
51  //-------------------------------------------------------------------
52  // this function is called by PETSc at the end of each nonlinear step
53  PetscErrorCode
54  __libmesh_petsc_snes_monitor (SNES, PetscInt its, PetscReal fnorm, void *)
55  {
56  //PetscErrorCode ierr=0;
57 
58  //if (its > 0)
59  libMesh::out << " NL step "
60  << std::setw(2) << its
61  << std::scientific
62  << ", |residual|_2 = " << fnorm
63  << std::endl;
64 
65  //return ierr;
66  return 0;
67  }
68 
69 
70 
71  //---------------------------------------------------------------
72  // this function is called by PETSc to evaluate the residual at X
73  PetscErrorCode
74  __libmesh_petsc_snes_residual (SNES snes, Vec x, Vec r, void * ctx)
75  {
76  LOG_SCOPE("residual()", "PetscNonlinearSolver");
77 
78  PetscErrorCode ierr=0;
79 
80  libmesh_assert(x);
81  libmesh_assert(r);
82  libmesh_assert(ctx);
83 
84  // No way to safety-check this cast, since we got a void *...
86  static_cast<PetscNonlinearSolver<Number> *> (ctx);
87 
88  // Get the current iteration number from the snes object,
89  // store it in the PetscNonlinearSolver object for possible use
90  // by the user's residual function.
91  {
92  PetscInt n_iterations = 0;
93  ierr = SNESGetIterationNumber(snes, &n_iterations);
94  CHKERRABORT(solver->comm().get(),ierr);
95  solver->_current_nonlinear_iteration_number = cast_int<unsigned>(n_iterations);
96  }
97 
98  NonlinearImplicitSystem & sys = solver->system();
99 
100  PetscVector<Number> & X_sys = *cast_ptr<PetscVector<Number> *>(sys.solution.get());
101  PetscVector<Number> X_global(x, sys.comm()), R(r, sys.comm());
102 
103  // Use the system's update() to get a good local version of the
104  // parallel solution. This operation does not modify the incoming
105  // "x" vector, it only localizes information from "x" into
106  // sys.current_local_solution.
107  X_global.swap(X_sys);
108  sys.update();
109  X_global.swap(X_sys);
110 
111  // Enforce constraints (if any) exactly on the
112  // current_local_solution. This is the solution vector that is
113  // actually used in the computation of the residual below, and is
114  // not locked by debug-enabled PETSc the way that "x" is.
115  sys.get_dof_map().enforce_constraints_exactly(sys, sys.current_local_solution.get());
116 
117  if (solver->_zero_out_residual)
118  R.zero();
119 
120  //-----------------------------------------------------------------------------
121  // if the user has provided both function pointers and objects only the pointer
122  // will be used, so catch that as an error
123  if (solver->residual && solver->residual_object)
124  libmesh_error_msg("ERROR: cannot specify both a function and object to compute the Residual!");
125 
126  if (solver->matvec && solver->residual_and_jacobian_object)
127  libmesh_error_msg("ERROR: cannot specify both a function and object to compute the combined Residual & Jacobian!");
128 
129  if (solver->residual != libmesh_nullptr)
130  solver->residual(*sys.current_local_solution.get(), R, sys);
131 
132  else if (solver->residual_object != libmesh_nullptr)
133  solver->residual_object->residual(*sys.current_local_solution.get(), R, sys);
134 
135  else if (solver->matvec != libmesh_nullptr)
136  solver->matvec (*sys.current_local_solution.get(), &R, libmesh_nullptr, sys);
137 
138  else if (solver->residual_and_jacobian_object != libmesh_nullptr)
139  solver->residual_and_jacobian_object->residual_and_jacobian (*sys.current_local_solution.get(), &R, libmesh_nullptr, sys);
140 
141  else
142  libmesh_error_msg("Error! Unable to compute residual and/or Jacobian!");
143 
144  R.close();
145 
146  return ierr;
147  }
148 
149 
150 
151  //---------------------------------------------------------------
152  // this function is called by PETSc to evaluate the Jacobian at X
153  PetscErrorCode
155 #if PETSC_RELEASE_LESS_THAN(3,5,0)
156  SNES snes, Vec x, Mat * jac, Mat * pc, MatStructure * msflag, void * ctx
157 #else
158  SNES snes, Vec x, Mat jac, Mat pc, void * ctx
159 #endif
160  )
161  {
162  LOG_SCOPE("jacobian()", "PetscNonlinearSolver");
163 
164  PetscErrorCode ierr=0;
165 
167 
168  // No way to safety-check this cast, since we got a void *...
170  static_cast<PetscNonlinearSolver<Number> *> (ctx);
171 
172  // Get the current iteration number from the snes object,
173  // store it in the PetscNonlinearSolver object for possible use
174  // by the user's Jacobian function.
175  {
176  PetscInt n_iterations = 0;
177  ierr = SNESGetIterationNumber(snes, &n_iterations);
178  CHKERRABORT(solver->comm().get(),ierr);
179  solver->_current_nonlinear_iteration_number = cast_int<unsigned>(n_iterations);
180  }
181 
182  NonlinearImplicitSystem & sys = solver->system();
183 #if PETSC_RELEASE_LESS_THAN(3,5,0)
184  PetscMatrix<Number> PC(*pc, sys.comm());
185  PetscMatrix<Number> Jac(*jac, sys.comm());
186 #else
187  PetscMatrix<Number> PC(pc, sys.comm());
188  PetscMatrix<Number> Jac(jac, sys.comm());
189 #endif
190  PetscVector<Number> & X_sys = *cast_ptr<PetscVector<Number> *>(sys.solution.get());
191  PetscVector<Number> X_global(x, sys.comm());
192 
193  // Set the dof maps
194  PC.attach_dof_map(sys.get_dof_map());
195  Jac.attach_dof_map(sys.get_dof_map());
196 
197  // Use the systems update() to get a good local version of the parallel solution
198  X_global.swap(X_sys);
199  sys.update();
200  X_global.swap(X_sys);
201 
202  // Enforce constraints (if any) exactly on the
203  // current_local_solution. This is the solution vector that is
204  // actually used in the computation of the residual below, and is
205  // not locked by debug-enabled PETSc the way that "x" is.
206  sys.get_dof_map().enforce_constraints_exactly(sys, sys.current_local_solution.get());
207 
208  if (solver->_zero_out_jacobian)
209  PC.zero();
210 
211  //-----------------------------------------------------------------------------
212  // if the user has provided both function pointers and objects only the pointer
213  // will be used, so catch that as an error
214  if (solver->jacobian && solver->jacobian_object)
215  libmesh_error_msg("ERROR: cannot specify both a function and object to compute the Jacobian!");
216 
217  if (solver->matvec && solver->residual_and_jacobian_object)
218  libmesh_error_msg("ERROR: cannot specify both a function and object to compute the combined Residual & Jacobian!");
219 
220  if (solver->jacobian != libmesh_nullptr)
221  solver->jacobian(*sys.current_local_solution.get(), PC, sys);
222 
223  else if (solver->jacobian_object != libmesh_nullptr)
224  solver->jacobian_object->jacobian(*sys.current_local_solution.get(), PC, sys);
225 
226  else if (solver->matvec != libmesh_nullptr)
227  solver->matvec(*sys.current_local_solution.get(), libmesh_nullptr, &PC, sys);
228 
229  else if (solver->residual_and_jacobian_object != libmesh_nullptr)
230  solver->residual_and_jacobian_object->residual_and_jacobian (*sys.current_local_solution.get(), libmesh_nullptr, &PC, sys);
231 
232  else
233  libmesh_error_msg("Error! Unable to compute residual and/or Jacobian!");
234 
235  PC.close();
236  Jac.close();
237 #if PETSC_RELEASE_LESS_THAN(3,5,0)
238  *msflag = SAME_NONZERO_PATTERN;
239 #endif
240 
241  return ierr;
242  }
243 
244  // This function gets called by PETSc after the SNES linesearch is
245  // complete. We use it to exactly enforce any constraints on the
246  // solution which may have drifted during the linear solve. In the
247  // PETSc nomenclature:
248  // * "x" is the old solution vector,
249  // * "y" is the search direction (Newton step) vector,
250  // * "w" is the candidate solution vector, and
251  // the user is responsible for setting changed_y and changed_w
252  // appropriately, depending on whether or not the search
253  // direction or solution vector was changed, respectively.
255 #if PETSC_VERSION_LESS_THAN(3,3,0)
256  SNES, Vec x, Vec y, Vec w, void * context, PetscBool * changed_y, PetscBool * changed_w
257 #else
258  SNESLineSearch, Vec x, Vec y, Vec w, PetscBool * changed_y, PetscBool * changed_w, void * context
259 #endif
260  )
261  {
262  LOG_SCOPE("postcheck()", "PetscNonlinearSolver");
263 
264  PetscErrorCode ierr = 0;
265 
266  // PETSc almost certainly initializes these to false already, but
267  // it doesn't hurt to be explicit.
268  *changed_w = PETSC_FALSE;
269  *changed_y = PETSC_FALSE;
270 
271  libmesh_assert(context);
272 
273  // Cast the context to a NonlinearSolver object.
275  static_cast<PetscNonlinearSolver<Number> *> (context);
276 
277  // If the user has provided both postcheck function pointer and
278  // object, this is ambiguous, so throw an error.
279  if (solver->postcheck && solver->postcheck_object)
280  libmesh_error_msg("ERROR: cannot specify both a function and object for performing the solve postcheck!");
281 
282  // It's also possible that we don't need to do anything at all, in
283  // that case return early...
284  NonlinearImplicitSystem & sys = solver->system();
285  DofMap & dof_map = sys.get_dof_map();
286 
287  if (!dof_map.n_constrained_dofs() && !solver->postcheck && !solver->postcheck_object)
288  return ierr;
289 
290  // We definitely need to wrap at least "w"
291  PetscVector<Number> petsc_w(w, sys.comm());
292 
293  // The user sets these flags in his/her postcheck function to
294  // indicate whether they changed something.
295  bool
296  changed_search_direction = false,
297  changed_new_soln = false;
298 
299  if (solver->postcheck || solver->postcheck_object)
300  {
301  PetscVector<Number> petsc_x(x, sys.comm());
302  PetscVector<Number> petsc_y(y, sys.comm());
303 
304  if (solver->postcheck)
305  solver->postcheck(petsc_x,
306  petsc_y,
307  petsc_w,
308  changed_search_direction,
309  changed_new_soln,
310  sys);
311 
312  else if (solver->postcheck_object)
313  solver->postcheck_object->postcheck(petsc_x,
314  petsc_y,
315  petsc_w,
316  changed_search_direction,
317  changed_new_soln,
318  sys);
319  }
320 
321  // Record whether the user changed the solution or the search direction.
322  if (changed_search_direction)
323  *changed_y = PETSC_TRUE;
324 
325  if (changed_new_soln)
326  *changed_w = PETSC_TRUE;
327 
328  if (dof_map.n_constrained_dofs())
329  {
330  PetscVector<Number> & system_soln = *cast_ptr<PetscVector<Number> *>(sys.solution.get());
331 
332  // ... and swap it in before enforcing the constraints.
333  petsc_w.swap(system_soln);
334 
335  dof_map.enforce_constraints_exactly(sys);
336 
337  // If we have constraints, we'll assume that we did change the
338  // solution w (hopefully slightly). Enforcing constraints
339  // does not change the search direction, y, but the user may
340  // have, so we leave it alone.
341  *changed_w = PETSC_TRUE;
342 
343  // Swap back
344  petsc_w.swap(system_soln);
345  }
346 
347  return ierr;
348  }
349 
350 } // end extern "C"
351 
352 
353 
354 //---------------------------------------------------------------------
355 // PetscNonlinearSolver<> methods
356 template <typename T>
358  NonlinearSolver<T>(system_in),
359  _reason(SNES_CONVERGED_ITERATING/*==0*/), // Arbitrary initial value...
360  _n_linear_iterations(0),
361  _current_nonlinear_iteration_number(0),
362  _zero_out_residual(true),
363  _zero_out_jacobian(true),
364  _default_monitor(true)
365 {
366 }
367 
368 
369 
370 template <typename T>
372 {
373  this->clear ();
374 }
375 
376 
377 
378 template <typename T>
380 {
381  if (this->initialized())
382  {
383  this->_is_initialized = false;
384 
385  PetscErrorCode ierr=0;
386 
387  ierr = LibMeshSNESDestroy(&_snes);
388  LIBMESH_CHKERR(ierr);
389 
390  // Reset the nonlinear iteration counter. This information is only relevant
391  // *during* the solve(). After the solve is completed it should return to
392  // the default value of 0.
394  }
395 }
396 
397 
398 
399 template <typename T>
401 {
402  // Initialize the data structures if not done so already.
403  if (!this->initialized())
404  {
405  this->_is_initialized = true;
406 
407  PetscErrorCode ierr=0;
408 
409  ierr = SNESCreate(this->comm().get(),&_snes);
410  LIBMESH_CHKERR(ierr);
411 
412  if (name)
413  {
414  ierr = SNESSetOptionsPrefix(_snes, name);
415  LIBMESH_CHKERR(ierr);
416  }
417 
418 #if !PETSC_RELEASE_LESS_THAN(3,3,0)
419  // Attaching a DM to SNES.
420  DM dm;
421  ierr = DMCreate(this->comm().get(), &dm);LIBMESH_CHKERR(ierr);
422  ierr = DMSetType(dm,DMLIBMESH);LIBMESH_CHKERR(ierr);
423  ierr = DMlibMeshSetSystem(dm,this->system());LIBMESH_CHKERR(ierr);
424  if (name)
425  {
426  ierr = DMSetOptionsPrefix(dm,name); LIBMESH_CHKERR(ierr);
427  }
428  ierr = DMSetFromOptions(dm); LIBMESH_CHKERR(ierr);
429  ierr = DMSetUp(dm); LIBMESH_CHKERR(ierr);
430  ierr = SNESSetDM(this->_snes, dm); LIBMESH_CHKERR(ierr);
431  // SNES now owns the reference to dm.
432  ierr = DMDestroy(&dm); LIBMESH_CHKERR(ierr);
433 
434 #endif
435 
436  if (_default_monitor)
437  {
438  ierr = SNESMonitorSet (_snes, __libmesh_petsc_snes_monitor,
439  this, PETSC_NULL);
440  LIBMESH_CHKERR(ierr);
441  }
442 
443  // If the SolverConfiguration object is provided, use it to set
444  // options during solver initialization.
445  if (this->_solver_configuration)
446  {
448  }
449 
450 #if PETSC_VERSION_LESS_THAN(3,1,0)
451  // Cannot call SNESSetOptions before SNESSetFunction when using
452  // any matrix free options with PETSc 3.1.0+
453  ierr = SNESSetFromOptions(_snes);
454  LIBMESH_CHKERR(ierr);
455 #endif
456 
457  if (this->_preconditioner)
458  {
459  KSP ksp;
460  ierr = SNESGetKSP (_snes, &ksp);
461  LIBMESH_CHKERR(ierr);
462  PC pc;
463  ierr = KSPGetPC(ksp,&pc);
464  LIBMESH_CHKERR(ierr);
465 
466  this->_preconditioner->init();
467 
468  PCSetType(pc, PCSHELL);
469  PCShellSetContext(pc,(void *)this->_preconditioner);
470 
471  //Re-Use the shell functions from petsc_linear_solver
472  PCShellSetSetUp(pc,__libmesh_petsc_preconditioner_setup);
473  PCShellSetApply(pc,__libmesh_petsc_preconditioner_apply);
474  }
475  }
476 
477 
478  // Tell PETSc about our linesearch "post-check" function, but only
479  // if the user has provided one. There seem to be extra,
480  // unnecessary residual calculations if a postcheck function is
481  // attached for no reason.
482  if (this->postcheck || this->postcheck_object)
483  {
484 #if PETSC_VERSION_LESS_THAN(3,3,0)
485  PetscErrorCode ierr = SNESLineSearchSetPostCheck(_snes, __libmesh_petsc_snes_postcheck, this);
486  LIBMESH_CHKERR(ierr);
487 
488 #else
489  // Pick the right version of the function name
490 #if PETSC_VERSION_LESS_THAN(3,4,0)
491 # define SNESGETLINESEARCH SNESGetSNESLineSearch
492 #else
493 # define SNESGETLINESEARCH SNESGetLineSearch
494 #endif
495 
496  SNESLineSearch linesearch;
497  PetscErrorCode ierr = SNESGETLINESEARCH(_snes, &linesearch);
498  LIBMESH_CHKERR(ierr);
499 
500  ierr = SNESLineSearchSetPostCheck(linesearch, __libmesh_petsc_snes_postcheck, this);
501  LIBMESH_CHKERR(ierr);
502 #endif
503  }
504 }
505 
506 #if !PETSC_VERSION_LESS_THAN(3,3,0)
507 template <typename T>
508 void
510  void (*computeSubspace)(std::vector<NumericVector<Number> *> &, sys_type &),
511  MatNullSpace * msp)
512 {
513  PetscErrorCode ierr;
514  std::vector<NumericVector<Number> * > sp;
515  if (computeSubspaceObject)
516  (*computeSubspaceObject)(sp, this->system());
517  else
518  (*computeSubspace)(sp, this->system());
519 
520  *msp = PETSC_NULL;
521  if (sp.size())
522  {
523  Vec * modes;
524  PetscScalar * dots;
525  PetscInt nmodes = cast_int<PetscInt>(sp.size());
526 
527 #if PETSC_RELEASE_LESS_THAN(3,5,0)
528  ierr = PetscMalloc2(nmodes,Vec,&modes,nmodes,PetscScalar,&dots);
529 #else
530  ierr = PetscMalloc2(nmodes,&modes,nmodes,&dots);
531 #endif
532  LIBMESH_CHKERR(ierr);
533 
534  for (PetscInt i=0; i<nmodes; ++i)
535  {
536  PetscVector<T> * pv = cast_ptr<PetscVector<T> *>(sp[i]);
537  Vec v = pv->vec();
538 
539  ierr = VecDuplicate(v, modes+i);
540  LIBMESH_CHKERR(ierr);
541 
542  ierr = VecCopy(v,modes[i]);
543  LIBMESH_CHKERR(ierr);
544  }
545 
546  // Normalize.
547  ierr = VecNormalize(modes[0],PETSC_NULL);
548  LIBMESH_CHKERR(ierr);
549 
550  for (PetscInt i=1; i<nmodes; i++)
551  {
552  // Orthonormalize vec[i] against vec[0:i-1]
553  ierr = VecMDot(modes[i],i,modes,dots);
554  LIBMESH_CHKERR(ierr);
555 
556  for (PetscInt j=0; j<i; j++)
557  dots[j] *= -1.;
558 
559  ierr = VecMAXPY(modes[i],i,dots,modes);
560  LIBMESH_CHKERR(ierr);
561 
562  ierr = VecNormalize(modes[i],PETSC_NULL);
563  LIBMESH_CHKERR(ierr);
564  }
565 
566  ierr = MatNullSpaceCreate(this->comm().get(), PETSC_FALSE, nmodes, modes, msp);
567  LIBMESH_CHKERR(ierr);
568 
569  for (PetscInt i=0; i<nmodes; ++i)
570  {
571  ierr = VecDestroy(modes+i);
572  LIBMESH_CHKERR(ierr);
573  }
574 
575  ierr = PetscFree2(modes,dots);
576  LIBMESH_CHKERR(ierr);
577  }
578 }
579 #endif
580 
581 template <typename T>
582 std::pair<unsigned int, Real>
583 PetscNonlinearSolver<T>::solve (SparseMatrix<T> & jac_in, // System Jacobian Matrix
584  NumericVector<T> & x_in, // Solution vector
585  NumericVector<T> & r_in, // Residual vector
586  const double, // Stopping tolerance
587  const unsigned int)
588 {
589  LOG_SCOPE("solve()", "PetscNonlinearSolver");
590  this->init ();
591 
592  // Make sure the data passed in are really of Petsc types
593  PetscMatrix<T> * jac = cast_ptr<PetscMatrix<T> *>(&jac_in);
594  PetscVector<T> * x = cast_ptr<PetscVector<T> *>(&x_in);
595  PetscVector<T> * r = cast_ptr<PetscVector<T> *>(&r_in);
596 
597  PetscErrorCode ierr=0;
598  PetscInt n_iterations =0;
599  // Should actually be a PetscReal, but I don't know which version of PETSc first introduced PetscReal
600  Real final_residual_norm=0.;
601 
602  ierr = SNESSetFunction (_snes, r->vec(), __libmesh_petsc_snes_residual, this);
603  LIBMESH_CHKERR(ierr);
604 
605  // Only set the jacobian function if we've been provided with something to call.
606  // This allows a user to set their own jacobian function if they want to
607  if (this->jacobian || this->jacobian_object || this->residual_and_jacobian_object)
608  {
609  ierr = SNESSetJacobian (_snes, jac->mat(), jac->mat(), __libmesh_petsc_snes_jacobian, this);
610  LIBMESH_CHKERR(ierr);
611  }
612 #if !PETSC_VERSION_LESS_THAN(3,3,0)
613  // Only set the nullspace if we have a way of computing it and the result is non-empty.
614  if (this->nullspace || this->nullspace_object)
615  {
616  MatNullSpace msp;
617  this->build_mat_null_space(this->nullspace_object, this->nullspace, &msp);
618  if (msp)
619  {
620  ierr = MatSetNullSpace(jac->mat(), msp);
621  LIBMESH_CHKERR(ierr);
622 
623  ierr = MatNullSpaceDestroy(&msp);
624  LIBMESH_CHKERR(ierr);
625  }
626  }
627 
628  // Only set the transpose nullspace if we have a way of computing it and the result is non-empty.
630  {
631 #if PETSC_VERSION_LESS_THAN(3,6,0)
632  libmesh_warning("MatSetTransposeNullSpace is only supported for PETSc >= 3.6, transpose nullspace will be ignored.");
633 #else
634  MatNullSpace msp = PETSC_NULL;
636  if (msp)
637  {
638  ierr = MatSetTransposeNullSpace(jac->mat(), msp);
639  LIBMESH_CHKERR(ierr);
640 
641  ierr = MatNullSpaceDestroy(&msp);
642  LIBMESH_CHKERR(ierr);
643  }
644 #endif
645  }
646 
647  // Only set the nearnullspace if we have a way of computing it and the result is non-empty.
648  if (this->nearnullspace || this->nearnullspace_object)
649  {
650  MatNullSpace msp = PETSC_NULL;
651  this->build_mat_null_space(this->nearnullspace_object, this->nearnullspace, &msp);
652 
653  if (msp)
654  {
655  ierr = MatSetNearNullSpace(jac->mat(), msp);
656  LIBMESH_CHKERR(ierr);
657 
658  ierr = MatNullSpaceDestroy(&msp);
659  LIBMESH_CHKERR(ierr);
660  }
661  }
662 #endif
663  // Have the Krylov subspace method use our good initial guess rather than 0
664  KSP ksp;
665  ierr = SNESGetKSP (_snes, &ksp);
666  LIBMESH_CHKERR(ierr);
667 
668  // Set the tolerances for the iterative solver. Use the user-supplied
669  // tolerance for the relative residual & leave the others at default values
670  ierr = KSPSetTolerances (ksp, this->initial_linear_tolerance, PETSC_DEFAULT,
671  PETSC_DEFAULT, this->max_linear_iterations);
672  LIBMESH_CHKERR(ierr);
673 
674  // Set the tolerances for the non-linear solver.
675  ierr = SNESSetTolerances(_snes, this->absolute_residual_tolerance, this->relative_residual_tolerance,
677  LIBMESH_CHKERR(ierr);
678 
679  //Pull in command-line options
680 #if PETSC_VERSION_LESS_THAN(3,7,0)
681  KSPSetFromOptions(ksp);
682 #endif
683  SNESSetFromOptions(_snes);
684 
685  if (this->user_presolve)
686  this->user_presolve(this->system());
687 
688  //Set the preconditioning matrix
689  if (this->_preconditioner)
690  {
691  this->_preconditioner->set_matrix(jac_in);
692  this->_preconditioner->init();
693  }
694 
695  // If the SolverConfiguration object is provided, use it to override
696  // solver options.
697  if (this->_solver_configuration)
698  {
700  }
701 
702  ierr = SNESSolve (_snes, PETSC_NULL, x->vec());
703  LIBMESH_CHKERR(ierr);
704 
705  ierr = SNESGetIterationNumber(_snes,&n_iterations);
706  LIBMESH_CHKERR(ierr);
707 
708  ierr = SNESGetLinearSolveIterations(_snes, &_n_linear_iterations);
709  LIBMESH_CHKERR(ierr);
710 
711  // Enforce constraints exactly now that the solve is done. We have
712  // been enforcing them on the current_local_solution during the
713  // solve, but now need to be sure they are enforced on the parallel
714  // solution vector as well.
716 
717  // SNESGetFunction has been around forever and should work on all
718  // versions of PETSc. This is also now the recommended approach
719  // according to the documentation for the PETSc 3.5.1 release:
720  // http://www.mcs.anl.gov/petsc/documentation/changes/35.html
721  Vec f;
722  ierr = SNESGetFunction(_snes, &f, 0, 0);
723  LIBMESH_CHKERR(ierr);
724  ierr = VecNorm(f, NORM_2, &final_residual_norm);
725  LIBMESH_CHKERR(ierr);
726 
727  // Get and store the reason for convergence
728  SNESGetConvergedReason(_snes, &_reason);
729 
730  //Based on Petsc 2.3.3 documentation all diverged reasons are negative
731  this->converged = (_reason >= 0);
732 
733  this->clear();
734 
735  // return the # of its. and the final residual norm.
736  return std::make_pair(n_iterations, final_residual_norm);
737 }
738 
739 
740 
741 template <typename T>
743 {
744 
745  libMesh::out << "Nonlinear solver convergence/divergence reason: "
746  << SNESConvergedReasons[this->get_converged_reason()] << std::endl;
747 }
748 
749 
750 
751 template <typename T>
753 {
754  PetscErrorCode ierr=0;
755 
756  if (this->initialized())
757  {
758  ierr = SNESGetConvergedReason(_snes, &_reason);
759  LIBMESH_CHKERR(ierr);
760  }
761 
762  return _reason;
763 }
764 
765 template <typename T>
767 {
768  return _n_linear_iterations;
769 }
770 
771 
772 //------------------------------------------------------------------
773 // Explicit instantiations
774 template class PetscNonlinearSolver<Number>;
775 
776 } // namespace libMesh
777 
778 
779 
780 #endif // #ifdef LIBMESH_HAVE_PETSC
std::string name(const ElemQuality q)
Definition: elem_quality.C:39
virtual void clear() libmesh_override
NumericVector interface to PETSc Vec.
Definition: petsc_vector.h:64
void(* residual)(const NumericVector< Number > &X, NumericVector< Number > &R, sys_type &S)
NonlinearImplicitSystem::ComputeResidualandJacobian * residual_and_jacobian_object
virtual void residual(const NumericVector< Number > &X, NumericVector< Number > &R, sys_type &S)=0
void(* transpose_nullspace)(std::vector< NumericVector< Number > * > &sp, sys_type &S)
PetscErrorCode __libmesh_petsc_snes_residual(SNES snes, Vec x, Vec r, void *ctx)
dof_id_type n_constrained_dofs() const
ImplicitSystem & sys
unsigned int max_function_evaluations
Preconditioner< T > * _preconditioner
PetscErrorCode __libmesh_petsc_snes_monitor(SNES, PetscInt its, PetscReal fnorm, void *)
void(* user_presolve)(sys_type &S)
void(* nullspace)(std::vector< NumericVector< Number > * > &sp, sys_type &S)
const class libmesh_nullptr_t libmesh_nullptr
PetscErrorCode __libmesh_petsc_snes_postcheck(#if PETSC_VERSION_LESS_THAN(3, 3, 0) SNES, Vec x, Vec y, Vec w, void *context, PetscBool *changed_y, PetscBool *changed_w#else SNESLineSearch, Vec x, Vec y, Vec w, PetscBool *changed_y, PetscBool *changed_w, void *context#endif)
virtual void configure_solver()=0
unsigned int max_linear_iterations
PetscErrorCode Vec Mat Mat pc
SolverConfiguration * _solver_configuration
NonlinearImplicitSystem::ComputeResidual * residual_object
void(* jacobian)(const NumericVector< Number > &X, SparseMatrix< Number > &J, sys_type &S)
virtual void postcheck(const NumericVector< Number > &old_soln, NumericVector< Number > &search_direction, NumericVector< Number > &new_soln, bool &changed_search_direction, bool &changed_new_soln, sys_type &S)=0
libmesh_assert(j)
PetscDiffSolver & solver
Manages the degrees of freedom (DOFs) in a simulation.
Definition: dof_map.h:167
NonlinearImplicitSystem::ComputeJacobian * jacobian_object
PetscErrorCode Vec x
const DofMap & get_dof_map() const
Definition: system.h:2028
virtual void jacobian(const NumericVector< Number > &X, SparseMatrix< Number > &J, sys_type &S)=0
PetscErrorCode __libmesh_petsc_snes_jacobian(#if PETSC_RELEASE_LESS_THAN(3, 5, 0) SNES snes, Vec x, Mat *jac, Mat *pc, MatStructure *msflag, void *ctx#else SNES snes, Vec x, Mat jac, Mat pc, void *ctx#endif)
Used for solving nonlinear implicit systems of equations.
PetscErrorCode Vec Mat Mat void * ctx
UniquePtr< NumericVector< Number > > solution
Definition: system.h:1521
virtual void init(const char *name=libmesh_nullptr) libmesh_override
void enforce_constraints_exactly(const System &system, NumericVector< Number > *v=libmesh_nullptr, bool homogeneous=false) const
PetscErrorCode __libmesh_petsc_preconditioner_setup(void *ctx)
friend PetscErrorCode __libmesh_petsc_snes_residual(SNES snes, Vec x, Vec r, void *ctx)
void build_mat_null_space(NonlinearImplicitSystem::ComputeVectorSubspace *computeSubspaceObject, void(*)(std::vector< NumericVector< Number > * > &, sys_type &), MatNullSpace *)
PetscErrorCode DMlibMeshSetSystem(DM dm, libMesh::NonlinearImplicitSystem &sys)
SNESConvergedReason get_converged_reason()
friend PetscErrorCode __libmesh_petsc_snes_jacobian(SNES snes, Vec x, Mat *jac, Mat *pc, MatStructure *msflag, void *ctx)
PetscTruth PetscBool
Definition: petsc_macro.h:64
const sys_type & system() const
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
void(* nearnullspace)(std::vector< NumericVector< Number > * > &sp, sys_type &S)
NonlinearImplicitSystem::ComputePostCheck * postcheck_object
virtual void print_converged_reason() libmesh_override
virtual std::pair< unsigned int, Real > solve(SparseMatrix< T > &, NumericVector< T > &, NumericVector< T > &, const double, const unsigned int) libmesh_override
PetscErrorCode ierr
const Parallel::Communicator & comm() const
SparseMatrix interface to PETSc Mat.
void(* postcheck)(const NumericVector< Number > &old_soln, NumericVector< Number > &search_direction, NumericVector< Number > &new_soln, bool &changed_search_direction, bool &changed_new_soln, sys_type &S)
NonlinearImplicitSystem::ComputeVectorSubspace * nearnullspace_object
virtual void swap(NumericVector< T > &v) libmesh_override
OStreamProxy out(std::cout)
PetscErrorCode __libmesh_petsc_preconditioner_apply(void *ctx, Vec x, Vec y)
virtual int get_total_linear_iterations() libmesh_override
unsigned int max_nonlinear_iterations
NonlinearImplicitSystem::ComputeVectorSubspace * transpose_nullspace_object
virtual void residual_and_jacobian(const NumericVector< Number > &X, NumericVector< Number > *R, SparseMatrix< Number > *J, sys_type &S)=0
void(* matvec)(const NumericVector< Number > &X, NumericVector< Number > *R, SparseMatrix< Number > *J, sys_type &S)
NonlinearImplicitSystem::ComputeVectorSubspace * nullspace_object