petsc_dm_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 #include "libmesh/petsc_macro.h"
20 
21 // This only works with petsc-3.3 and above.
22 #if !PETSC_VERSION_LESS_THAN(3,3,0)
23 
24 // Local Includes
25 #include "libmesh/libmesh_common.h"
29 #include "libmesh/petsc_vector.h"
30 #include "libmesh/petsc_matrix.h"
31 #include "libmesh/dof_map.h"
32 #include "libmesh/preconditioner.h"
33 
34 extern "C"
35 {
36 #if PETSC_VERSION_LESS_THAN(2,2,1)
37  typedef int PetscErrorCode;
38  typedef int PetscInt;
39 #endif
40 }
41 
42 // PETSc extern definition
43 EXTERN_C_BEGIN
45 EXTERN_C_END
46 
47 #include <libmesh/petscdmlibmesh.h>
48 
49 namespace libMesh {
52  {
54  return;
55 
57 #if PETSC_RELEASE_LESS_THAN(3,4,0)
58  ierr = DMRegister(DMLIBMESH, PETSC_NULL, "DMCreate_libMesh", DMCreate_libMesh); CHKERRABORT(libMesh::COMM_WORLD,ierr);
59 #else
60  ierr = DMRegister(DMLIBMESH, DMCreate_libMesh); CHKERRABORT(libMesh::COMM_WORLD,ierr);
61 #endif
62  PetscDMRegistered = PETSC_TRUE;
63  }
64 
65 
66  template <typename T>
68  PetscNonlinearSolver<T>(system_in)
69  {
71  }
72 
73  template <typename T>
74  inline
76  {
77  this->clear ();
78  }
79 
80 
81 
82  template <typename T>
84  {
86  DM dm;
88 
89  // Attaching a DM with the function and Jacobian callbacks to SNES.
90  ierr = DMCreateLibMesh(this->comm().get(), this->system(), &dm); LIBMESH_CHKERRABORT(ierr);
91  ierr = DMSetFromOptions(dm); LIBMESH_CHKERRABORT(ierr);
92  ierr = DMSetUp(dm); LIBMESH_CHKERRABORT(ierr);
93  ierr = SNESSetDM(this->_snes, dm); LIBMESH_CHKERRABORT(ierr);
94  // SNES now owns the reference to dm.
95  ierr = DMDestroy(&dm); LIBMESH_CHKERRABORT(ierr);
96  KSP ksp;
97  ierr = SNESGetKSP (this->_snes, &ksp); LIBMESH_CHKERRABORT(ierr);
98 
99  // Set the tolerances for the iterative solver. Use the user-supplied
100  // tolerance for the relative residual & leave the others at default values
101  ierr = KSPSetTolerances (ksp, this->initial_linear_tolerance, PETSC_DEFAULT,PETSC_DEFAULT, this->max_linear_iterations); LIBMESH_CHKERRABORT(ierr);
102 
103  // Set the tolerances for the non-linear solver.
104  ierr = SNESSetTolerances(this->_snes,
105  this->absolute_residual_tolerance,
106  this->relative_residual_tolerance,
107  this->absolute_step_tolerance,
108  this->max_nonlinear_iterations,
109  this->max_function_evaluations);
110  LIBMESH_CHKERRABORT(ierr);
111 
112  //Pull in command-line options
113  KSPSetFromOptions(ksp);
114  SNESSetFromOptions(this->_snes);
115  }
116 
117 
118  template <typename T>
119  std::pair<unsigned int, Real>
120  PetscDMNonlinearSolver<T>::solve (SparseMatrix<T>& jac_in, // System Jacobian Matrix
121  NumericVector<T>& x_in, // Solution vector
122  NumericVector<T>& r_in, // Residual vector
123  const double, // Stopping tolerance
124  const unsigned int)
125  {
126  START_LOG("solve()", "PetscNonlinearSolver");
127  this->init ();
128 
129  // Make sure the data passed in are really of Petsc types
130  libmesh_cast_ptr<PetscMatrix<T>*>(&jac_in);
131  libmesh_cast_ptr<PetscVector<T>*>(&r_in);
132 
133  // Extract solution vector
134  PetscVector<T>* x = libmesh_cast_ptr<PetscVector<T>*>(&x_in);
135 
137  PetscInt n_iterations =0;
138 
139 
140  PetscReal final_residual_norm=0.;
141 
142  if (this->user_presolve)
143  this->user_presolve(this->system());
144 
145  //Set the preconditioning matrix
146  if (this->_preconditioner)
147  this->_preconditioner->set_matrix(jac_in);
148 
149  ierr = SNESSolve (this->_snes, PETSC_NULL, x->vec());
150  LIBMESH_CHKERRABORT(ierr);
151 
152  ierr = SNESGetIterationNumber(this->_snes,&n_iterations);
153  LIBMESH_CHKERRABORT(ierr);
154 
155  ierr = SNESGetLinearSolveIterations(this->_snes, &this->_n_linear_iterations);
156  LIBMESH_CHKERRABORT(ierr);
157 
158  ierr = SNESGetFunctionNorm(this->_snes,&final_residual_norm);
159  LIBMESH_CHKERRABORT(ierr);
160 
161  // Get and store the reason for convergence
162  SNESGetConvergedReason(this->_snes, &this->_reason);
163 
164  //Based on Petsc 2.3.3 documentation all diverged reasons are negative
165  this->converged = (this->_reason >= 0);
166 
167  this->clear();
168 
169  STOP_LOG("solve()", "PetscNonlinearSolver");
170 
171  // return the # of its. and the final residual norm.
172  return std::make_pair(n_iterations, final_residual_norm);
173  }
174 
175 
176  //------------------------------------------------------------------
177  // Explicit instantiations
178  template class PetscDMNonlinearSolver<Number>;
179 
180 } // namespace libMesh
181 
182 
183 #endif // #if !PETSC_VERSION_LESS_THAN(3,3,0)

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

Hosted By:
SourceForge.net Logo