libMesh
trilinos_nox_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 #if defined(LIBMESH_TRILINOS_HAVE_NOX) && defined(LIBMESH_TRILINOS_HAVE_EPETRA)
23 
24 // Local Includes
25 #include "libmesh/libmesh_logging.h"
26 #include "libmesh/dof_map.h"
27 #include "libmesh/nonlinear_implicit_system.h"
28 #include "libmesh/trilinos_nox_nonlinear_solver.h"
29 #include "libmesh/system.h"
30 #include "libmesh/trilinos_epetra_vector.h"
31 #include "libmesh/trilinos_epetra_matrix.h"
32 #include "libmesh/trilinos_preconditioner.h"
33 
34 // Trilinos includes
35 #include "libmesh/ignore_warnings.h"
36 #include "NOX_Epetra_MatrixFree.H"
37 #include "NOX_Epetra_LinearSystem_AztecOO.H"
38 #include "NOX_Epetra_Group.H"
39 #include "NOX_Epetra_Vector.H"
40 #include "libmesh/restore_warnings.h"
41 
42 namespace libMesh
43 {
44 
45 class Problem_Interface : public NOX::Epetra::Interface::Required,
46  public NOX::Epetra::Interface::Jacobian,
47  public NOX::Epetra::Interface::Preconditioner
48 {
49 public:
50  explicit
52 
54 
56  bool computeF(const Epetra_Vector & x,
57  Epetra_Vector & FVec,
58  NOX::Epetra::Interface::Required::FillType fillType);
59 
61  bool computeJacobian(const Epetra_Vector & x,
62  Epetra_Operator & Jac);
63 
70  bool computePrecMatrix(const Epetra_Vector & x,
71  Epetra_RowMatrix & M);
72 
75  bool computePreconditioner(const Epetra_Vector & x,
76  Epetra_Operator & Prec,
77  Teuchos::ParameterList * p);
78 
80 };
81 
82 
83 
85  _solver(solver)
86 {}
87 
88 
89 
91 {}
92 
93 
94 
95 bool Problem_Interface::computeF(const Epetra_Vector & x,
96  Epetra_Vector & r,
97  NOX::Epetra::Interface::Required::FillType /*fillType*/)
98 {
99  LOG_SCOPE("residual()", "TrilinosNoxNonlinearSolver");
100 
102 
103  EpetraVector<Number> X_global(*const_cast<Epetra_Vector *>(&x), sys.comm()), R(r, sys.comm());
104  EpetraVector<Number> & X_sys = *cast_ptr<EpetraVector<Number> *>(sys.solution.get());
105  EpetraVector<Number> & R_sys = *cast_ptr<EpetraVector<Number> *>(sys.rhs);
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  sys.update();
113 
114  // Swap back
115  X_global.swap(X_sys);
116  R.swap(R_sys);
117 
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 
125  libmesh_error_msg("ERROR: cannot specify both a function and object to compute the Residual!");
126 
128  libmesh_error_msg("ERROR: cannot specify both a function and object to compute the combined Residual & Jacobian!");
129 
131  _solver->residual(*sys.current_local_solution.get(), R, sys);
132 
135 
136  else if (_solver->matvec != libmesh_nullptr)
138 
139  else if (_solver->residual_and_jacobian_object != libmesh_nullptr)
141 
142  else
143  return false;
144 
145  R.close();
146  X_global.close();
147 
148  return true;
149 }
150 
151 
152 
153 bool Problem_Interface::computeJacobian(const Epetra_Vector & x,
154  Epetra_Operator & jac)
155 {
156  LOG_SCOPE("jacobian()", "TrilinosNoxNonlinearSolver");
157 
159 
160  EpetraMatrix<Number> Jac(&dynamic_cast<Epetra_FECrsMatrix &>(jac), sys.comm());
161  EpetraVector<Number> & X_sys = *cast_ptr<EpetraVector<Number> *>(sys.solution.get());
162  EpetraVector<Number> X_global(*const_cast<Epetra_Vector *>(&x), sys.comm());
163 
164  // Set the dof maps
165  Jac.attach_dof_map(sys.get_dof_map());
166 
167  // Use the systems update() to get a good local version of the parallel solution
168  X_global.swap(X_sys);
169 
171  sys.update();
172 
173  X_global.swap(X_sys);
174 
175  //-----------------------------------------------------------------------------
176  // if the user has provided both function pointers and objects only the pointer
177  // will be used, so catch that as an error
179  libmesh_error_msg("ERROR: cannot specify both a function and object to compute the Jacobian!");
180 
182  libmesh_error_msg("ERROR: cannot specify both a function and object to compute the combined Residual & Jacobian!");
183 
185  _solver->jacobian(*sys.current_local_solution.get(), Jac, sys);
186 
189 
190  else if (_solver->matvec != libmesh_nullptr)
192 
195 
196  else
197  libmesh_error_msg("Error! Unable to compute residual and/or Jacobian!");
198 
199  Jac.close();
200  X_global.close();
201 
202  return true;
203 }
204 
205 
206 
207 bool Problem_Interface::computePrecMatrix(const Epetra_Vector & /*x*/, Epetra_RowMatrix & /*M*/)
208 {
209  // libMesh::out << "ERROR: Problem_Interface::preconditionVector() - Use Explicit Jacobian only for this test problem!" << endl;
210  throw 1;
211 }
212 
213 
214 
215 bool Problem_Interface::computePreconditioner(const Epetra_Vector & x,
216  Epetra_Operator & prec,
217  Teuchos::ParameterList * /*p*/)
218 {
219  LOG_SCOPE("preconditioner()", "TrilinosNoxNonlinearSolver");
220 
223 
224  EpetraMatrix<Number> Jac(dynamic_cast<Epetra_FECrsMatrix *>(tpc.mat()),sys.comm());
225  EpetraVector<Number> & X_sys = *cast_ptr<EpetraVector<Number> *>(sys.solution.get());
226  EpetraVector<Number> X_global(*const_cast<Epetra_Vector *>(&x), sys.comm());
227 
228  // Set the dof maps
229  Jac.attach_dof_map(sys.get_dof_map());
230 
231  // Use the systems update() to get a good local version of the parallel solution
232  X_global.swap(X_sys);
233 
235  sys.update();
236 
237  X_global.swap(X_sys);
238 
239  //-----------------------------------------------------------------------------
240  // if the user has provided both function pointers and objects only the pointer
241  // will be used, so catch that as an error
243  libmesh_error_msg("ERROR: cannot specify both a function and object to compute the Jacobian!");
244 
246  libmesh_error_msg("ERROR: cannot specify both a function and object to compute the combined Residual & Jacobian!");
247 
249  _solver->jacobian(*sys.current_local_solution.get(), Jac, sys);
250 
253 
254  else if (_solver->matvec != libmesh_nullptr)
256 
259 
260  else
261  libmesh_error_msg("Error! Unable to compute residual and/or Jacobian!");
262 
263  Jac.close();
264  X_global.close();
265 
266  tpc.compute();
267 
268  return true;
269 }
270 
271 
272 
273 //---------------------------------------------------------------------
274 // NoxNonlinearSolver<> methods
275 template <typename T>
277 {
278 }
279 
280 
281 
282 template <typename T>
283 void NoxNonlinearSolver<T>::init (const char * /*name*/)
284 {
285  if (!this->initialized())
286  _interface = new Problem_Interface(this);
287 }
288 
289 
290 
291 template <typename T>
292 std::pair<unsigned int, Real>
293 NoxNonlinearSolver<T>::solve (SparseMatrix<T> & /* jac_in */, // System Jacobian Matrix
294  NumericVector<T> & x_in, // Solution vector
295  NumericVector<T> & /* r_in */, // Residual vector
296  const double, // Stopping tolerance
297  const unsigned int)
298 {
299  this->init ();
300 
301  if (this->user_presolve)
302  this->user_presolve(this->system());
303 
304  EpetraVector<T> * x_epetra = cast_ptr<EpetraVector<T> *>(&x_in);
305  // Creating a Teuchos::RCP as they do in NOX examples does not work here - we get some invalid memory references
306  // thus we make a local copy
307  NOX::Epetra::Vector x(*x_epetra->vec());
308 
309  Teuchos::RCP<Teuchos::ParameterList> nlParamsPtr = Teuchos::rcp(new Teuchos::ParameterList);
310  Teuchos::ParameterList & nlParams = *(nlParamsPtr.get());
311  nlParams.set("Nonlinear Solver", "Line Search Based");
312 
313  //print params
314  Teuchos::ParameterList & printParams = nlParams.sublist("Printing");
315  printParams.set("Output Precision", 3);
316  printParams.set("Output Processor", 0);
317  printParams.set("Output Information",
318  NOX::Utils::OuterIteration +
319  NOX::Utils::OuterIterationStatusTest +
320  NOX::Utils::InnerIteration +
321  NOX::Utils::LinearSolverDetails +
322  NOX::Utils::Parameters +
323  NOX::Utils::Details +
324  NOX::Utils::Warning);
325 
326  Teuchos::ParameterList & dirParams = nlParams.sublist("Direction");
327  dirParams.set("Method", "Newton");
328  Teuchos::ParameterList & newtonParams = dirParams.sublist("Newton");
329  newtonParams.set("Forcing Term Method", "Constant");
330 
331  Teuchos::ParameterList & lsParams = newtonParams.sublist("Linear Solver");
332  lsParams.set("Aztec Solver", "GMRES");
333  lsParams.set("Max Iterations", static_cast<int>(this->max_linear_iterations));
334  lsParams.set("Tolerance", this->initial_linear_tolerance);
335  lsParams.set("Output Frequency", 1);
336  lsParams.set("Size of Krylov Subspace", 1000);
337 
338  //create linear system
339  Teuchos::RCP<NOX::Epetra::Interface::Required> iReq(_interface);
340  Teuchos::RCP<NOX::Epetra::LinearSystemAztecOO> linSys;
341  Teuchos::RCP<Epetra_Operator> pc;
342 
343  if (this->jacobian || this->jacobian_object || this->residual_and_jacobian_object)
344  {
345  if (this->_preconditioner)
346  {
347  // PJNFK
348  lsParams.set("Preconditioner", "User Defined");
349 
350  TrilinosPreconditioner<Number> * trilinos_pc =
351  cast_ptr<TrilinosPreconditioner<Number> *>(this->_preconditioner);
352  pc = Teuchos::rcp(trilinos_pc);
353 
354  Teuchos::RCP<NOX::Epetra::Interface::Preconditioner> iPrec(_interface);
355  linSys = Teuchos::rcp(new NOX::Epetra::LinearSystemAztecOO(printParams, lsParams, iReq, iPrec, pc, x));
356  }
357  else
358  {
359  lsParams.set("Preconditioner", "None");
360  // lsParams.set("Preconditioner", "Ifpack");
361  // lsParams.set("Preconditioner", "AztecOO");
362 
363  // full jacobian
364  NonlinearImplicitSystem & sys = _interface->_solver->system();
365  EpetraMatrix<Number> & jacSys = *cast_ptr<EpetraMatrix<Number> *>(sys.matrix);
366  Teuchos::RCP<Epetra_RowMatrix> jacMat = Teuchos::rcp(jacSys.mat());
367 
368  Teuchos::RCP<NOX::Epetra::Interface::Jacobian> iJac(_interface);
369  linSys = Teuchos::rcp(new NOX::Epetra::LinearSystemAztecOO(printParams, lsParams, iReq, iJac, jacMat, x));
370  }
371  }
372  else
373  {
374  // matrix free
375  Teuchos::RCP<NOX::Epetra::MatrixFree> MF = Teuchos::rcp(new NOX::Epetra::MatrixFree(printParams, iReq, x));
376 
377  Teuchos::RCP<NOX::Epetra::Interface::Jacobian> iJac(MF);
378  linSys = Teuchos::rcp(new NOX::Epetra::LinearSystemAztecOO(printParams, lsParams, iReq, iJac, MF, x));
379  }
380 
381  //create group
382  Teuchos::RCP<NOX::Epetra::Group> grpPtr = Teuchos::rcp(new NOX::Epetra::Group(printParams, iReq, x, linSys));
383  NOX::Epetra::Group & grp = *(grpPtr.get());
384 
385  Teuchos::RCP<NOX::StatusTest::NormF> absresid =
386  Teuchos::rcp(new NOX::StatusTest::NormF(this->absolute_residual_tolerance, NOX::StatusTest::NormF::Unscaled));
387  Teuchos::RCP<NOX::StatusTest::NormF> relresid =
388  Teuchos::rcp(new NOX::StatusTest::NormF(grp, this->relative_residual_tolerance));
389  Teuchos::RCP<NOX::StatusTest::MaxIters> maxiters =
390  Teuchos::rcp(new NOX::StatusTest::MaxIters(this->max_nonlinear_iterations));
391  Teuchos::RCP<NOX::StatusTest::FiniteValue> finiteval =
392  Teuchos::rcp(new NOX::StatusTest::FiniteValue());
393  Teuchos::RCP<NOX::StatusTest::NormUpdate> normupdate =
394  Teuchos::rcp(new NOX::StatusTest::NormUpdate(this->absolute_step_tolerance));
395  Teuchos::RCP<NOX::StatusTest::Combo> combo =
396  Teuchos::rcp(new NOX::StatusTest::Combo(NOX::StatusTest::Combo::OR));
397  combo->addStatusTest(absresid);
398  combo->addStatusTest(relresid);
399  combo->addStatusTest(maxiters);
400  combo->addStatusTest(finiteval);
401  combo->addStatusTest(normupdate);
402 
403  Teuchos::RCP<Teuchos::ParameterList> finalPars = nlParamsPtr;
404 
405  Teuchos::RCP<NOX::Solver::Generic> solver = NOX::Solver::buildSolver(grpPtr, combo, nlParamsPtr);
406  NOX::StatusTest::StatusType status = solver->solve();
407  this->converged = (status == NOX::StatusTest::Converged);
408 
409  const NOX::Epetra::Group & finalGroup = dynamic_cast<const NOX::Epetra::Group &>(solver->getSolutionGroup());
410  const NOX::Epetra::Vector & noxFinalSln = dynamic_cast<const NOX::Epetra::Vector &>(finalGroup.getX());
411 
412  *x_epetra->vec() = noxFinalSln.getEpetraVector();
413  x_in.close();
414 
415  Real residual_norm = finalGroup.getNormF();
416  unsigned int total_iters = solver->getNumIterations();
417  _n_linear_iterations = finalPars->sublist("Direction").sublist("Newton").sublist("Linear Solver").sublist("Output").get("Total Number of Linear Iterations", -1);
418 
419  // do not let Trilinos to deallocate what we allocated
420  pc.release();
421  iReq.release();
422 
423  return std::make_pair(total_iters, residual_norm);
424 }
425 
426 
427 
428 template <typename T>
429 int
431 {
432  return _n_linear_iterations;
433 }
434 
435 
436 
437 //------------------------------------------------------------------
438 // Explicit instantiations
439 template class NoxNonlinearSolver<Number>;
440 
441 } // namespace libMesh
442 
443 
444 
445 #endif // LIBMESH_TRILINOS_HAVE_NOX && LIBMESH_TRILINOS_HAVE_EPETRA
void(* residual)(const NumericVector< Number > &X, NumericVector< Number > &R, sys_type &S)
Function that computes the residual R(X) of the nonlinear system at the input iterate X...
NonlinearImplicitSystem::ComputeResidualandJacobian * residual_and_jacobian_object
Object that computes either the residual or the Jacobian of the nonlinear system at the input itera...
virtual void residual(const NumericVector< Number > &X, NumericVector< Number > &R, sys_type &S)=0
Residual function.
bool computePrecMatrix(const Epetra_Vector &x, Epetra_RowMatrix &M)
Compute the Epetra_RowMatrix M, that will be used by the Aztec preconditioner instead of the Jacobian...
ImplicitSystem & sys
UniquePtr< NumericVector< Number > > current_local_solution
All the values I need to compute my contribution to the simulation at hand.
Definition: system.h:1535
const class libmesh_nullptr_t libmesh_nullptr
NumericVector< Number > * rhs
The system matrix.
Numeric vector.
Definition: dof_map.h:66
NoxNonlinearSolver< Number > * _solver
The libMesh namespace provides an interface to certain functionality in the library.
PetscErrorCode Vec Mat Mat pc
NonlinearImplicitSystem::ComputeResidual * residual_object
Object that computes the residual R(X) of the nonlinear system at the input iterate X...
void(* jacobian)(const NumericVector< Number > &X, SparseMatrix< Number > &J, sys_type &S)
Function that computes the Jacobian J(X) of the nonlinear system at the input iterate X...
This class provides an interface to the suite of preconditioners available from Trilinos.
PetscDiffSolver & solver
Generic sparse matrix.
Definition: dof_map.h:65
NonlinearImplicitSystem::ComputeJacobian * jacobian_object
Object that computes the Jacobian J(X) of the nonlinear system at the input iterate X...
Problem_Interface(NoxNonlinearSolver< Number > *solver)
PetscErrorCode Vec x
const DofMap & get_dof_map() const
Definition: system.h:2030
void compute()
Compute the preconditioner.
virtual void jacobian(const NumericVector< Number > &X, SparseMatrix< Number > &J, sys_type &S)=0
Jacobian function.
bool computePreconditioner(const Epetra_Vector &x, Epetra_Operator &Prec, Teuchos::ParameterList *p)
Computes a user supplied preconditioner based on input vector x.
virtual void clear() libmesh_override
Release all memory and clear data structures.
void init(triangulateio &t)
Initializes the fields of t to NULL/0 as necessary.
This class provides a specific system class.
MPI_Status status
Status object for querying messages.
Definition: parallel.h:176
virtual void init(const char *name=libmesh_nullptr) libmesh_override
Initialize data structures if not done so already.
UniquePtr< NumericVector< Number > > solution
Data structure to hold solution values.
Definition: system.h:1523
void enforce_constraints_exactly(const System &system, NumericVector< Number > *v=libmesh_nullptr, bool homogeneous=false) const
Constrains the numeric vector v, which represents a solution defined on the mesh. ...
Definition: dof_map.h:1816
virtual void close()=0
Calls the NumericVector&#39;s internal assembly routines, ensuring that the values are consistent across ...
virtual void update()
Update the local values to reflect the solution on neighboring processors.
Definition: system.C:425
const sys_type & system() const
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual int get_total_linear_iterations() libmesh_override
Get the total number of linear iterations done in the last solve.
const Parallel::Communicator & comm() const
SparseMatrix< Number > * matrix
The system matrix.
bool initialized()
Checks that library initialization has been done.
Definition: libmesh.C:274
bool computeJacobian(const Epetra_Vector &x, Epetra_Operator &Jac)
Compute an explicit Jacobian.
virtual std::pair< unsigned int, Real > solve(SparseMatrix< T > &, NumericVector< T > &, NumericVector< T > &, const double, const unsigned int) libmesh_override
Call the Nox solver.
virtual void residual_and_jacobian(const NumericVector< Number > &X, NumericVector< Number > *R, SparseMatrix< Number > *J, sys_type &S)=0
Residual & Jacobian function, calculated simultaneously.
void(* matvec)(const NumericVector< Number > &X, NumericVector< Number > *R, SparseMatrix< Number > *J, sys_type &S)
Function that computes either the residual or the Jacobian of the nonlinear system at the input ite...
bool computeF(const Epetra_Vector &x, Epetra_Vector &FVec, NOX::Epetra::Interface::Required::FillType fillType)
Compute and return F.