petsc_nonlinear_solver.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 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"
35 
36 namespace libMesh
37 {
38 
39 //--------------------------------------------------------------------
40 // Functions with C linkage to pass to PETSc. PETSc will call these
41 // methods as needed.
42 //
43 // Since they must have C linkage they have no knowledge of a namespace.
44 // Give them an obscure name to avoid namespace pollution.
45 extern "C"
46 {
47  // Older versions of PETSc do not have the different int typedefs.
48  // On 64-bit machines, PetscInt may actually be a long long int.
49  // This change occurred in Petsc-2.2.1.
50 #if PETSC_VERSION_LESS_THAN(2,2,1)
51  typedef int PetscErrorCode;
52  typedef int PetscInt;
53 #endif
54 
55  //-------------------------------------------------------------------
56  // this function is called by PETSc at the end of each nonlinear step
57  PetscErrorCode
58  __libmesh_petsc_snes_monitor (SNES, PetscInt its, PetscReal fnorm, void *)
59  {
60  //PetscErrorCode ierr=0;
61 
62  //if (its > 0)
63  libMesh::out << " NL step "
64  << std::setw(2) << its
65  << std::scientific
66  << ", |residual|_2 = " << fnorm
67  << std::endl;
68 
69  //return ierr;
70  return 0;
71  }
72 
73 
74 
75  //---------------------------------------------------------------
76  // this function is called by PETSc to evaluate the residual at X
77  PetscErrorCode
78  __libmesh_petsc_snes_residual (SNES snes, Vec x, Vec r, void *ctx)
79  {
80  START_LOG("residual()", "PetscNonlinearSolver");
81 
82  PetscErrorCode ierr=0;
83 
84  libmesh_assert(x);
85  libmesh_assert(r);
86  libmesh_assert(ctx);
87 
89  static_cast<PetscNonlinearSolver<Number>*> (ctx);
90 
91  // Get the current iteration number from the snes object,
92  // store it in the PetscNonlinearSolver object for possible use
93  // by the user's residual function.
94  {
95  PetscInt n_iterations = 0;
96  ierr = SNESGetIterationNumber(snes, &n_iterations);
97  CHKERRABORT(solver->comm().get(),ierr);
98  solver->_current_nonlinear_iteration_number = static_cast<unsigned>(n_iterations);
99  }
100 
101  NonlinearImplicitSystem &sys = solver->system();
102 
103  PetscVector<Number>& X_sys = *libmesh_cast_ptr<PetscVector<Number>*>(sys.solution.get());
104  PetscVector<Number>& R_sys = *libmesh_cast_ptr<PetscVector<Number>*>(sys.rhs);
105  PetscVector<Number> X_global(x, sys.comm()), R(r, sys.comm());
106 
107  // Use the systems update() to get a good local version of the parallel solution
108  X_global.swap(X_sys);
109  R.swap(R_sys);
110 
112 
113  sys.update();
114 
115  //Swap back
116  X_global.swap(X_sys);
117  R.swap(R_sys);
118 
119  if (solver->_zero_out_residual)
120  R.zero();
121 
122  //-----------------------------------------------------------------------------
123  // if the user has provided both function pointers and objects only the pointer
124  // will be used, so catch that as an error
125  if (solver->residual && solver->residual_object)
126  {
127  libMesh::err << "ERROR: cannot specifiy both a function and object to compute the Residual!" << std::endl;
128  libmesh_error();
129  }
130 
131  if (solver->matvec && solver->residual_and_jacobian_object)
132  {
133  libMesh::err << "ERROR: cannot specifiy both a function and object to compute the combined Residual & Jacobian!" << std::endl;
134  libmesh_error();
135  }
136  //-----------------------------------------------------------------------------
137 
138  if (solver->residual != NULL) solver->residual (*sys.current_local_solution.get(), R, sys);
139  else if (solver->residual_object != NULL) solver->residual_object->residual (*sys.current_local_solution.get(), R, sys);
140  else if (solver->matvec != NULL) solver->matvec (*sys.current_local_solution.get(), &R, NULL, sys);
141  else if (solver->residual_and_jacobian_object != NULL) solver->residual_and_jacobian_object->residual_and_jacobian (*sys.current_local_solution.get(), &R, NULL, sys);
142  else libmesh_error();
143 
144  R.close();
145  X_global.close();
146 
147  STOP_LOG("residual()", "PetscNonlinearSolver");
148 
149  return ierr;
150  }
151 
152 
153 
154  //---------------------------------------------------------------
155  // this function is called by PETSc to evaluate the Jacobian at X
156  PetscErrorCode
157  __libmesh_petsc_snes_jacobian (SNES snes, Vec x, Mat *jac, Mat *pc, MatStructure *msflag, void *ctx)
158  {
159  START_LOG("jacobian()", "PetscNonlinearSolver");
160 
161  PetscErrorCode ierr=0;
162 
163  libmesh_assert(ctx);
164 
166  static_cast<PetscNonlinearSolver<Number>*> (ctx);
167 
168  // Get the current iteration number from the snes object,
169  // store it in the PetscNonlinearSolver object for possible use
170  // by the user's Jacobian function.
171  {
172  PetscInt n_iterations = 0;
173  ierr = SNESGetIterationNumber(snes, &n_iterations);
174  CHKERRABORT(solver->comm().get(),ierr);
175  solver->_current_nonlinear_iteration_number = static_cast<unsigned>(n_iterations);
176  }
177 
178  NonlinearImplicitSystem &sys = solver->system();
179 
180  PetscMatrix<Number> PC(*pc, sys.comm());
181  PetscMatrix<Number> Jac(*jac, sys.comm());
182  PetscVector<Number>& X_sys = *libmesh_cast_ptr<PetscVector<Number>*>(sys.solution.get());
183  PetscMatrix<Number>& Jac_sys = *libmesh_cast_ptr<PetscMatrix<Number>*>(sys.matrix);
184  PetscVector<Number> X_global(x, sys.comm());
185 
186  // Set the dof maps
187  PC.attach_dof_map(sys.get_dof_map());
188  Jac.attach_dof_map(sys.get_dof_map());
189 
190  // Use the systems update() to get a good local version of the parallel solution
191  X_global.swap(X_sys);
192  Jac.swap(Jac_sys);
193 
195  sys.update();
196 
197  X_global.swap(X_sys);
198  Jac.swap(Jac_sys);
199 
200  if (solver->_zero_out_jacobian)
201  PC.zero();
202 
203  //-----------------------------------------------------------------------------
204  // if the user has provided both function pointers and objects only the pointer
205  // will be used, so catch that as an error
206  if (solver->jacobian && solver->jacobian_object)
207  {
208  libMesh::err << "ERROR: cannot specify both a function and object to compute the Jacobian!" << std::endl;
209  libmesh_error();
210  }
211 
212  if (solver->matvec && solver->residual_and_jacobian_object)
213  {
214  libMesh::err << "ERROR: cannot specify both a function and object to compute the combined Residual & Jacobian!" << std::endl;
215  libmesh_error();
216  }
217  //-----------------------------------------------------------------------------
218 
219  if (solver->jacobian != NULL) solver->jacobian (*sys.current_local_solution.get(), PC, sys);
220  else if (solver->jacobian_object != NULL) solver->jacobian_object->jacobian (*sys.current_local_solution.get(), PC, sys);
221  else if (solver->matvec != NULL) solver->matvec (*sys.current_local_solution.get(), NULL, &PC, sys);
222  else if (solver->residual_and_jacobian_object != NULL) solver->residual_and_jacobian_object->residual_and_jacobian (*sys.current_local_solution.get(), NULL, &PC, sys);
223  else libmesh_error();
224 
225  PC.close();
226  Jac.close();
227  X_global.close();
228 
229  *msflag = SAME_NONZERO_PATTERN;
230 
231  STOP_LOG("jacobian()", "PetscNonlinearSolver");
232 
233  return ierr;
234  }
235 
236 } // end extern "C"
237 //---------------------------------------------------------------------
238 
239 
240 
241 //---------------------------------------------------------------------
242 // PetscNonlinearSolver<> methods
243 template <typename T>
245  NonlinearSolver<T>(system_in),
246  _reason(SNES_CONVERGED_ITERATING/*==0*/), // Arbitrary initial value...
247  _n_linear_iterations(0),
248  _current_nonlinear_iteration_number(0),
249  _zero_out_residual(true),
250  _zero_out_jacobian(true),
251  _default_monitor(true)
252 {
253 }
254 
255 
256 
257 template <typename T>
259 {
260  this->clear ();
261 }
262 
263 
264 
265 template <typename T>
267 {
268  if (this->initialized())
269  {
270  this->_is_initialized = false;
271 
273 
274  ierr = LibMeshSNESDestroy(&_snes);
275  LIBMESH_CHKERRABORT(ierr);
276 
277  // Reset the nonlinear iteration counter. This information is only relevant
278  // *during* the solve(). After the solve is completed it should return to
279  // the default value of 0.
280  _current_nonlinear_iteration_number = 0;
281  }
282 }
283 
284 
285 
286 template <typename T>
288 {
289  // Initialize the data structures if not done so already.
290  if (!this->initialized())
291  {
292  this->_is_initialized = true;
293 
295 
296 #if PETSC_VERSION_LESS_THAN(2,1,2)
297  // At least until Petsc 2.1.1, the SNESCreate had a different calling syntax.
298  // The second argument was of type SNESProblemType, and could have a value of
299  // either SNES_NONLINEAR_EQUATIONS or SNES_UNCONSTRAINED_MINIMIZATION.
300  ierr = SNESCreate(this->comm().get(), SNES_NONLINEAR_EQUATIONS, &_snes);
301  LIBMESH_CHKERRABORT(ierr);
302 
303 #else
304 
305  ierr = SNESCreate(this->comm().get(),&_snes);
306  LIBMESH_CHKERRABORT(ierr);
307 
308 #endif
309 
310  if (_default_monitor)
311  {
312 #if PETSC_VERSION_LESS_THAN(2,3,3)
313  ierr = SNESSetMonitor (_snes, __libmesh_petsc_snes_monitor,
314  this, PETSC_NULL);
315 #else
316  // API name change in PETSc 2.3.3
317  ierr = SNESMonitorSet (_snes, __libmesh_petsc_snes_monitor,
318  this, PETSC_NULL);
319 #endif
320  LIBMESH_CHKERRABORT(ierr);
321  }
322 
323 #if PETSC_VERSION_LESS_THAN(3,1,0)
324  // Cannot call SNESSetOptions before SNESSetFunction when using
325  // any matrix free options with PETSc 3.1.0+
326  ierr = SNESSetFromOptions(_snes);
327  LIBMESH_CHKERRABORT(ierr);
328 #endif
329 
330  if(this->_preconditioner)
331  {
332  KSP ksp;
333  ierr = SNESGetKSP (_snes, &ksp);
334  LIBMESH_CHKERRABORT(ierr);
335  PC pc;
336  ierr = KSPGetPC(ksp,&pc);
337  LIBMESH_CHKERRABORT(ierr);
338 
339  this->_preconditioner->init();
340 
341  PCSetType(pc, PCSHELL);
342  PCShellSetContext(pc,(void*)this->_preconditioner);
343 
344  //Re-Use the shell functions from petsc_linear_solver
345  PCShellSetSetUp(pc,__libmesh_petsc_preconditioner_setup);
346  PCShellSetApply(pc,__libmesh_petsc_preconditioner_apply);
347  }
348  }
349 }
350 
351 #if !PETSC_VERSION_LESS_THAN(3,3,0)
352 template <typename T>
353 void
355  void (*computeSubspace)(std::vector<NumericVector<Number>*>&, sys_type&),
356  MatNullSpace *msp)
357 {
359  std::vector<NumericVector<Number>* > sp;
360  if (computeSubspaceObject)
361  (*computeSubspaceObject)(sp, this->system());
362  else
363  (*computeSubspace)(sp, this->system());
364 
365  *msp = PETSC_NULL;
366  if (sp.size())
367  {
368  Vec *modes;
369  PetscScalar *dots;
370  PetscInt nmodes = sp.size();
371 
372 #if PETSC_RELEASE_LESS_THAN(3,5,0)
373  ierr = PetscMalloc2(nmodes,Vec,&modes,nmodes,PetscScalar,&dots);
374 #else
375  ierr = PetscMalloc2(nmodes,&modes,nmodes,&dots);
376 #endif
377  LIBMESH_CHKERRABORT(ierr);
378 
379  for (PetscInt i=0; i<nmodes; ++i)
380  {
381  PetscVector<T>* pv = libmesh_cast_ptr<PetscVector<T>*>(sp[i]);
382  Vec v = pv->vec();
383 
384  ierr = VecDuplicate(v, modes+i);
385  LIBMESH_CHKERRABORT(ierr);
386 
387  ierr = VecCopy(v,modes[i]);
388  LIBMESH_CHKERRABORT(ierr);
389  }
390 
391  // Normalize.
392  ierr = VecNormalize(modes[0],PETSC_NULL);
393  LIBMESH_CHKERRABORT(ierr);
394 
395  for (PetscInt i=1; i<nmodes; i++)
396  {
397  // Orthonormalize vec[i] against vec[0:i-1]
398  ierr = VecMDot(modes[i],i,modes,dots);
399  LIBMESH_CHKERRABORT(ierr);
400 
401  for (PetscInt j=0; j<i; j++)
402  dots[j] *= -1.;
403 
404  ierr = VecMAXPY(modes[i],i,dots,modes);
405  LIBMESH_CHKERRABORT(ierr);
406 
407  ierr = VecNormalize(modes[i],PETSC_NULL);
408  LIBMESH_CHKERRABORT(ierr);
409  }
410 
411  ierr = MatNullSpaceCreate(this->comm().get(), PETSC_FALSE, nmodes, modes, msp);
412  LIBMESH_CHKERRABORT(ierr);
413 
414  for (PetscInt i=0; i<nmodes; ++i)
415  {
416  ierr = VecDestroy(modes+i);
417  LIBMESH_CHKERRABORT(ierr);
418  }
419 
420  ierr = PetscFree2(modes,dots);
421  LIBMESH_CHKERRABORT(ierr);
422  }
423 }
424 #endif
425 
426 template <typename T>
427 std::pair<unsigned int, Real>
428 PetscNonlinearSolver<T>::solve (SparseMatrix<T>& jac_in, // System Jacobian Matrix
429  NumericVector<T>& x_in, // Solution vector
430  NumericVector<T>& r_in, // Residual vector
431  const double, // Stopping tolerance
432  const unsigned int)
433 {
434  START_LOG("solve()", "PetscNonlinearSolver");
435  this->init ();
436 
437  // Make sure the data passed in are really of Petsc types
438  PetscMatrix<T>* jac = libmesh_cast_ptr<PetscMatrix<T>*>(&jac_in);
439  PetscVector<T>* x = libmesh_cast_ptr<PetscVector<T>*>(&x_in);
440  PetscVector<T>* r = libmesh_cast_ptr<PetscVector<T>*>(&r_in);
441 
443  PetscInt n_iterations =0;
444  // Should actually be a PetscReal, but I don't know which version of PETSc first introduced PetscReal
445  Real final_residual_norm=0.;
446 
447  ierr = SNESSetFunction (_snes, r->vec(), __libmesh_petsc_snes_residual, this);
448  LIBMESH_CHKERRABORT(ierr);
449 
450  // Only set the jacobian function if we've been provided with something to call.
451  // This allows a user to set their own jacobian function if they want to
452  if (this->jacobian || this->jacobian_object || this->residual_and_jacobian_object)
453  {
454  ierr = SNESSetJacobian (_snes, jac->mat(), jac->mat(), __libmesh_petsc_snes_jacobian, this);
455  LIBMESH_CHKERRABORT(ierr);
456  }
457 #if !PETSC_VERSION_LESS_THAN(3,3,0)
458  // Only set the nullspace if we have a way of computing it and the result is non-empty.
459  if (this->nullspace || this->nullspace_object)
460  {
461  MatNullSpace msp;
462  this->build_mat_null_space(this->nullspace_object, this->nullspace, &msp);
463  if (msp)
464  {
465  ierr = MatSetNullSpace(jac->mat(), msp);
466  LIBMESH_CHKERRABORT(ierr);
467 
468  ierr = MatNullSpaceDestroy(&msp);
469  LIBMESH_CHKERRABORT(ierr);
470  }
471  }
472 
473  // Only set the nearnullspace if we have a way of computing it and the result is non-empty.
474  if (this->nearnullspace || this->nearnullspace_object)
475  {
476  MatNullSpace msp = PETSC_NULL;
477  this->build_mat_null_space(this->nearnullspace_object, this->nearnullspace, &msp);
478 
479  if(msp) {
480  ierr = MatSetNearNullSpace(jac->mat(), msp);
481  LIBMESH_CHKERRABORT(ierr);
482 
483  ierr = MatNullSpaceDestroy(&msp);
484  LIBMESH_CHKERRABORT(ierr);
485  }
486  }
487 #endif
488  // Have the Krylov subspace method use our good initial guess rather than 0
489  KSP ksp;
490  ierr = SNESGetKSP (_snes, &ksp);
491  LIBMESH_CHKERRABORT(ierr);
492 
493  // Set the tolerances for the iterative solver. Use the user-supplied
494  // tolerance for the relative residual & leave the others at default values
495  ierr = KSPSetTolerances (ksp, this->initial_linear_tolerance, PETSC_DEFAULT,
496  PETSC_DEFAULT, this->max_linear_iterations);
497  LIBMESH_CHKERRABORT(ierr);
498 
499  // Set the tolerances for the non-linear solver.
500  ierr = SNESSetTolerances(_snes, this->absolute_residual_tolerance, this->relative_residual_tolerance,
501  this->relative_step_tolerance, this->max_nonlinear_iterations, this->max_function_evaluations);
502  LIBMESH_CHKERRABORT(ierr);
503 
504  //Pull in command-line options
505  KSPSetFromOptions(ksp);
506  SNESSetFromOptions(_snes);
507 
508  if (this->user_presolve)
509  this->user_presolve(this->system());
510 
511  //Set the preconditioning matrix
512  if(this->_preconditioner)
513  {
514  this->_preconditioner->set_matrix(jac_in);
515  this->_preconditioner->init();
516  }
517 
518 // ierr = KSPSetInitialGuessNonzero (ksp, PETSC_TRUE);
519 // LIBMESH_CHKERRABORT(ierr);
520 
521 // Older versions (at least up to 2.1.5) of SNESSolve took 3 arguments,
522 // the last one being a pointer to an int to hold the number of iterations required.
523 # if PETSC_VERSION_LESS_THAN(2,2,0)
524 
525  ierr = SNESSolve (_snes, x->vec(), &n_iterations);
526  LIBMESH_CHKERRABORT(ierr);
527 
528 // 2.2.x style
529 #elif PETSC_VERSION_LESS_THAN(2,3,0)
530 
531  ierr = SNESSolve (_snes, x->vec());
532  LIBMESH_CHKERRABORT(ierr);
533 
534 // 2.3.x & newer style
535 #else
536 
537  ierr = SNESSolve (_snes, PETSC_NULL, x->vec());
538  LIBMESH_CHKERRABORT(ierr);
539 
540  ierr = SNESGetIterationNumber(_snes,&n_iterations);
541  LIBMESH_CHKERRABORT(ierr);
542 
543  ierr = SNESGetLinearSolveIterations(_snes, &_n_linear_iterations);
544  LIBMESH_CHKERRABORT(ierr);
545 
546  ierr = SNESGetFunctionNorm(_snes,&final_residual_norm);
547  LIBMESH_CHKERRABORT(ierr);
548 
549 #endif
550 
551  // Get and store the reason for convergence
552  SNESGetConvergedReason(_snes, &_reason);
553 
554  //Based on Petsc 2.3.3 documentation all diverged reasons are negative
555  this->converged = (_reason >= 0);
556 
557  this->clear();
558 
559  STOP_LOG("solve()", "PetscNonlinearSolver");
560 
561  // return the # of its. and the final residual norm.
562  return std::make_pair(n_iterations, final_residual_norm);
563 }
564 
565 
566 
567 template <typename T>
569 {
570 
571  libMesh::out << "Nonlinear solver convergence/divergence reason: "
572  << SNESConvergedReasons[this->get_converged_reason()] << std::endl;
573 }
574 
575 
576 
577 template <typename T>
579 {
581 
582  if (this->initialized())
583  {
584  ierr = SNESGetConvergedReason(_snes, &_reason);
585  LIBMESH_CHKERRABORT(ierr);
586  }
587 
588  return _reason;
589 }
590 
591 template <typename T>
593 {
594  return _n_linear_iterations;
595 }
596 
597 
598 //------------------------------------------------------------------
599 // Explicit instantiations
600 template class PetscNonlinearSolver<Number>;
601 
602 } // namespace libMesh
603 
604 
605 
606 #endif // #ifdef LIBMESH_HAVE_PETSC

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:06 UTC

Hosted By:
SourceForge.net Logo