trilinos_nox_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_NOX
23 
24 
25 // C++ includes
26 
27 // Local Includes
29 #include "libmesh/dof_map.h"
32 #include "libmesh/system.h"
36 
37 // ---------- Standard Includes ----------
38 #include <iostream>
39 #include "Epetra_Vector.h"
40 #include "Epetra_Operator.h"
41 #include "Epetra_RowMatrix.h"
42 #include "NOX_Epetra_Interface_Required.H" // base class
43 #include "NOX_Epetra_Interface_Jacobian.H" // base class
44 #include "NOX_Epetra_Interface_Preconditioner.H" // base class
45 #include "NOX_Epetra_MatrixFree.H"
46 #include "NOX_Epetra_LinearSystem_AztecOO.H"
47 #include "NOX_Epetra_Group.H" // class definition
48 #include "NOX_Epetra_Vector.H"
49 
50 namespace libMesh
51 {
52 
53 class Problem_Interface : public NOX::Epetra::Interface::Required,
54  public NOX::Epetra::Interface::Jacobian,
55  public NOX::Epetra::Interface::Preconditioner
56 {
57 public:
58  explicit
60 
62 
64  bool computeF(const Epetra_Vector& x, Epetra_Vector& FVec,
65  NOX::Epetra::Interface::Required::FillType fillType);
66 
68  bool computeJacobian(const Epetra_Vector& x, Epetra_Operator& Jac);
69 
71  bool computePrecMatrix(const Epetra_Vector& x, Epetra_RowMatrix& M);
72 
74  bool computePreconditioner(const Epetra_Vector& x, Epetra_Operator& Prec,
75  Teuchos::ParameterList* p);
76 
78 };
79 
80 //-----------------------------------------------------------------------------
82  _solver(solver)
83 { }
84 
86 { }
87 
88 bool Problem_Interface::computeF(const Epetra_Vector& x, Epetra_Vector& r,
89  NOX::Epetra::Interface::Required::FillType /*fillType*/)
90 {
91  START_LOG("residual()", "TrilinosNoxNonlinearSolver");
92 
94 
95  EpetraVector<Number> X_global(*const_cast<Epetra_Vector *>(&x), sys.comm()), R(r, sys.comm());
96  EpetraVector<Number>& X_sys = *libmesh_cast_ptr<EpetraVector<Number>*>(sys.solution.get());
97  EpetraVector<Number>& R_sys = *libmesh_cast_ptr<EpetraVector<Number>*>(sys.rhs);
98 
99  // Use the systems update() to get a good local version of the parallel solution
100  X_global.swap(X_sys);
101  R.swap(R_sys);
102 
104  sys.update();
105 
106  // Swap back
107  X_global.swap(X_sys);
108  R.swap(R_sys);
109 
110  R.zero();
111 
112  //-----------------------------------------------------------------------------
113  // if the user has provided both function pointers and objects only the pointer
114  // will be used, so catch that as an error
115 
117  {
118  libMesh::err << "ERROR: cannot specifiy both a function and object to compute the Residual!" << std::endl;
119  libmesh_error();
120  }
121 
123  {
124  libMesh::err << "ERROR: cannot specifiy both a function and object to compute the combined Residual & Jacobian!" << std::endl;
125  libmesh_error();
126  }
127  //-----------------------------------------------------------------------------
128 
129  if (_solver->residual != NULL) _solver->residual (*sys.current_local_solution.get(), R, sys);
130  else if (_solver->residual_object != NULL) _solver->residual_object->residual (*sys.current_local_solution.get(), R, sys);
131  else if (_solver->matvec != NULL) _solver->matvec (*sys.current_local_solution.get(), &R, NULL, sys);
133  else return false;
134 
135  R.close();
136  X_global.close();
137 
138  STOP_LOG("residual()", "TrilinosNoxNonlinearSolver");
139 
140  return true;
141 }
142 
143 bool Problem_Interface::computeJacobian(const Epetra_Vector & x,
144  Epetra_Operator & jac)
145 {
146  START_LOG("jacobian()", "TrilinosNoxNonlinearSolver");
147 
149 
150  EpetraMatrix<Number> Jac(&dynamic_cast<Epetra_FECrsMatrix &>(jac), sys.comm());
151  EpetraVector<Number>& X_sys = *libmesh_cast_ptr<EpetraVector<Number>*>(sys.solution.get());
152  EpetraVector<Number> X_global(*const_cast<Epetra_Vector *>(&x), sys.comm());
153 
154  // Set the dof maps
155  Jac.attach_dof_map(sys.get_dof_map());
156 
157  // Use the systems update() to get a good local version of the parallel solution
158  X_global.swap(X_sys);
159 
161  sys.update();
162 
163  X_global.swap(X_sys);
164 
165  //-----------------------------------------------------------------------------
166  // if the user has provided both function pointers and objects only the pointer
167  // will be used, so catch that as an error
169  {
170  libMesh::err << "ERROR: cannot specify both a function and object to compute the Jacobian!" << std::endl;
171  libmesh_error();
172  }
173 
175  {
176  libMesh::err << "ERROR: cannot specify both a function and object to compute the combined Residual & Jacobian!" << std::endl;
177  libmesh_error();
178  }
179  //-----------------------------------------------------------------------------
180 
181  if (_solver->jacobian != NULL) _solver->jacobian (*sys.current_local_solution.get(), Jac, sys);
182  else if (_solver->jacobian_object != NULL) _solver->jacobian_object->jacobian (*sys.current_local_solution.get(), Jac, sys);
183  else if (_solver->matvec != NULL) _solver->matvec (*sys.current_local_solution.get(), NULL, &Jac, sys);
185  else libmesh_error();
186 
187  Jac.close();
188  X_global.close();
189 
190  STOP_LOG("jacobian()", "TrilinosNoxNonlinearSolver");
191 
192  return true;
193 }
194 
195 bool Problem_Interface::computePrecMatrix(const Epetra_Vector & /*x*/, Epetra_RowMatrix & /*M*/)
196 {
197 // libMesh::out << "ERROR: Problem_Interface::preconditionVector() - Use Explicit Jacobian only for this test problem!" << endl;
198  throw 1;
199 }
200 
201 bool Problem_Interface::computePreconditioner(const Epetra_Vector & x,
202  Epetra_Operator & prec,
203  Teuchos::ParameterList * /*p*/)
204 {
205  START_LOG("preconditioner()", "TrilinosNoxNonlinearSolver");
206 
209 
210  EpetraMatrix<Number> Jac(dynamic_cast<Epetra_FECrsMatrix *>(tpc.mat()),sys.comm());
211  EpetraVector<Number>& X_sys = *libmesh_cast_ptr<EpetraVector<Number>*>(sys.solution.get());
212  EpetraVector<Number> X_global(*const_cast<Epetra_Vector *>(&x), sys.comm());
213 
214  // Set the dof maps
215  Jac.attach_dof_map(sys.get_dof_map());
216 
217  // Use the systems update() to get a good local version of the parallel solution
218  X_global.swap(X_sys);
219 
221  sys.update();
222 
223  X_global.swap(X_sys);
224 
225  //-----------------------------------------------------------------------------
226  // if the user has provided both function pointers and objects only the pointer
227  // will be used, so catch that as an error
229  {
230  libMesh::err << "ERROR: cannot specify both a function and object to compute the Jacobian!" << std::endl;
231  libmesh_error();
232  }
233 
235  {
236  libMesh::err << "ERROR: cannot specify both a function and object to compute the combined Residual & Jacobian!" << std::endl;
237  libmesh_error();
238  }
239  //-----------------------------------------------------------------------------
240 
241  if (_solver->jacobian != NULL) _solver->jacobian (*sys.current_local_solution.get(), Jac, sys);
242  else if (_solver->jacobian_object != NULL) _solver->jacobian_object->jacobian (*sys.current_local_solution.get(), Jac, sys);
243  else if (_solver->matvec != NULL) _solver->matvec (*sys.current_local_solution.get(), NULL, &Jac, sys);
245  else libmesh_error();
246 
247  Jac.close();
248  X_global.close();
249 
250  tpc.compute();
251 
252  STOP_LOG("preconditioner()", "TrilinosNoxNonlinearSolver");
253 
254  return true;
255 }
256 
257 
258 //---------------------------------------------------------------------
259 // NoxNonlinearSolver<> methods
260 template <typename T>
262 {
263 }
264 
265 template <typename T>
267 {
268  if (!this->initialized())
269  _interface = new Problem_Interface(this);
270 }
271 
272 template <typename T>
273 std::pair<unsigned int, Real>
274 NoxNonlinearSolver<T>::solve (SparseMatrix<T>& /* jac_in */, // System Jacobian Matrix
275  NumericVector<T>& x_in, // Solution vector
276  NumericVector<T>& /* r_in */, // Residual vector
277  const double, // Stopping tolerance
278  const unsigned int)
279 {
280  this->init ();
281 
282  if (this->user_presolve)
283  this->user_presolve(this->system());
284 
285  EpetraVector<T> * x_epetra = libmesh_cast_ptr<EpetraVector<T>*>(&x_in);
286  // Creating a Teuchos::RCP as they do in NOX examples does not work here - we get some invalid memory references
287  // thus we make a local copy
288  NOX::Epetra::Vector x(*x_epetra->vec());
289 
290  Teuchos::RCP<Teuchos::ParameterList> nlParamsPtr = Teuchos::rcp(new Teuchos::ParameterList);
291  Teuchos::ParameterList& nlParams = *(nlParamsPtr.get());
292  nlParams.set("Nonlinear Solver", "Line Search Based");
293 
294  //print params
295  Teuchos::ParameterList& printParams = nlParams.sublist("Printing");
296  printParams.set("Output Precision", 3);
297  printParams.set("Output Processor", 0);
298  printParams.set("Output Information",
299  NOX::Utils::OuterIteration +
300  NOX::Utils::OuterIterationStatusTest +
301  NOX::Utils::InnerIteration +
302  NOX::Utils::LinearSolverDetails +
303  NOX::Utils::Parameters +
304  NOX::Utils::Details +
305  NOX::Utils::Warning);
306 
307  Teuchos::ParameterList& dirParams = nlParams.sublist("Direction");
308  dirParams.set("Method", "Newton");
309  Teuchos::ParameterList& newtonParams = dirParams.sublist("Newton");
310  newtonParams.set("Forcing Term Method", "Constant");
311 
312  Teuchos::ParameterList& lsParams = newtonParams.sublist("Linear Solver");
313  lsParams.set("Aztec Solver", "GMRES");
314  lsParams.set("Max Iterations", static_cast<int>(this->max_linear_iterations));
315  lsParams.set("Tolerance", this->initial_linear_tolerance);
316  lsParams.set("Output Frequency", 1);
317  lsParams.set("Size of Krylov Subspace", 1000);
318 
319  //create linear system
320  Teuchos::RCP<NOX::Epetra::Interface::Required> iReq(_interface);
321  Teuchos::RCP<NOX::Epetra::LinearSystemAztecOO> linSys;
322  Teuchos::RCP<Epetra_Operator> pc;
323 
324  if (this->jacobian || this->jacobian_object || this->residual_and_jacobian_object)
325  {
326  if(this->_preconditioner)
327  {
328  // PJNFK
329  lsParams.set("Preconditioner", "User Defined");
330 
331  TrilinosPreconditioner<Number> * trilinos_pc = libmesh_cast_ptr<TrilinosPreconditioner<Number> *>(this->_preconditioner);
332  pc = Teuchos::rcp(trilinos_pc);
333 
334  Teuchos::RCP<NOX::Epetra::Interface::Preconditioner> iPrec(_interface);
335  linSys = Teuchos::rcp(new NOX::Epetra::LinearSystemAztecOO(printParams, lsParams, iReq, iPrec, pc, x));
336  }
337  else
338  {
339  lsParams.set("Preconditioner", "None");
340 // lsParams.set("Preconditioner", "Ifpack");
341 // lsParams.set("Preconditioner", "AztecOO");
342 
343  // full jacobian
344  NonlinearImplicitSystem & sys = _interface->_solver->system();
345  EpetraMatrix<Number> & jacSys = *libmesh_cast_ptr<EpetraMatrix<Number>*>(sys.matrix);
346  Teuchos::RCP<Epetra_RowMatrix> jacMat = Teuchos::rcp(jacSys.mat());
347 
348  Teuchos::RCP<NOX::Epetra::Interface::Jacobian> iJac(_interface);
349  linSys = Teuchos::rcp(new NOX::Epetra::LinearSystemAztecOO(printParams, lsParams, iReq, iJac, jacMat, x));
350  }
351  }
352  else
353  {
354  // matrix free
355  Teuchos::RCP<NOX::Epetra::MatrixFree> MF = Teuchos::rcp(new NOX::Epetra::MatrixFree(printParams, iReq, x));
356 
357  Teuchos::RCP<NOX::Epetra::Interface::Jacobian> iJac(MF);
358  linSys = Teuchos::rcp(new NOX::Epetra::LinearSystemAztecOO(printParams, lsParams, iReq, iJac, MF, x));
359  }
360 
361  //create group
362  Teuchos::RCP<NOX::Epetra::Group> grpPtr = Teuchos::rcp(new NOX::Epetra::Group(printParams, iReq, x, linSys));
363  NOX::Epetra::Group& grp = *(grpPtr.get());
364 
365  Teuchos::RCP<NOX::StatusTest::NormF> absresid =
366  Teuchos::rcp(new NOX::StatusTest::NormF(this->absolute_residual_tolerance, NOX::StatusTest::NormF::Unscaled));
367  Teuchos::RCP<NOX::StatusTest::NormF> relresid =
368  Teuchos::rcp(new NOX::StatusTest::NormF(grp, this->relative_residual_tolerance));
369  Teuchos::RCP<NOX::StatusTest::MaxIters> maxiters =
370  Teuchos::rcp(new NOX::StatusTest::MaxIters(this->max_nonlinear_iterations));
371  Teuchos::RCP<NOX::StatusTest::FiniteValue> finiteval =
372  Teuchos::rcp(new NOX::StatusTest::FiniteValue());
373  Teuchos::RCP<NOX::StatusTest::NormUpdate> normupdate =
374  Teuchos::rcp(new NOX::StatusTest::NormUpdate(this->absolute_step_tolerance));
375  Teuchos::RCP<NOX::StatusTest::Combo> combo =
376  Teuchos::rcp(new NOX::StatusTest::Combo(NOX::StatusTest::Combo::OR));
377  combo->addStatusTest(absresid);
378  combo->addStatusTest(relresid);
379  combo->addStatusTest(maxiters);
380  combo->addStatusTest(finiteval);
381  combo->addStatusTest(normupdate);
382 
383  Teuchos::RCP<Teuchos::ParameterList> finalPars = nlParamsPtr;
384 
385  Teuchos::RCP<NOX::Solver::Generic> solver = NOX::Solver::buildSolver(grpPtr, combo, nlParamsPtr);
386  NOX::StatusTest::StatusType status = solver->solve();
387  this->converged = (status == NOX::StatusTest::Converged);
388 
389  const NOX::Epetra::Group& finalGroup = dynamic_cast<const NOX::Epetra::Group&>(solver->getSolutionGroup());
390  const NOX::Epetra::Vector& noxFinalSln = dynamic_cast<const NOX::Epetra::Vector&>(finalGroup.getX());
391 
392  *x_epetra->vec() = noxFinalSln.getEpetraVector();
393  x_in.close();
394 
395  Real residual_norm = finalGroup.getNormF();
396  unsigned int total_iters = solver->getNumIterations();
397  _n_linear_iterations = finalPars->sublist("Direction").sublist("Newton").sublist("Linear Solver").sublist("Output").get("Total Number of Linear Iterations", -1);
398 
399  // do not let Trilinos to deallocate what we allocated
400  pc.release();
401  iReq.release();
402 
403  return std::make_pair(total_iters, residual_norm);
404 }
405 
406 template <typename T>
407 int
409 {
410  return _n_linear_iterations;
411 }
412 
413 
414 
415 //------------------------------------------------------------------
416 // Explicit instantiations
417 template class NoxNonlinearSolver<Number>;
418 
419 } // namespace libMesh
420 
421 
422 
423 #endif // #ifdef LIBMESH_HAVE_NOX

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

Hosted By:
SourceForge.net Logo