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