libMesh
trilinos_nox_nonlinear_solver.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2025 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 bool Problem_Interface::computeF(const Epetra_Vector & x,
95  Epetra_Vector & r,
96  NOX::Epetra::Interface::Required::FillType /*fillType*/)
97 {
98  LOG_SCOPE("residual()", "TrilinosNoxNonlinearSolver");
99 
101 
102  EpetraVector<Number> X_global(*const_cast<Epetra_Vector *>(&x), sys.comm()), R(r, sys.comm());
103  EpetraVector<Number> & X_sys = *cast_ptr<EpetraVector<Number> *>(sys.solution.get());
104  EpetraVector<Number> & R_sys = *cast_ptr<EpetraVector<Number> *>(sys.rhs);
105 
106  // Use the systems update() to get a good local version of the parallel solution
107  X_global.swap(X_sys);
108  R.swap(R_sys);
109 
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 
124  libmesh_error_msg_if(_solver->residual && _solver->residual_object, "ERROR: cannot specify both a function and object to compute the Residual!");
125 
126  libmesh_error_msg_if(_solver->matvec && _solver->residual_and_jacobian_object,
127  "ERROR: cannot specify both a function and object to compute the combined Residual & Jacobian!");
128 
129  if (_solver->residual != nullptr)
130  _solver->residual(*sys.current_local_solution.get(), R, sys);
131 
132  else if (_solver->residual_object != nullptr)
134 
135  else if (_solver->matvec != nullptr)
136  _solver->matvec(*sys.current_local_solution.get(), &R, nullptr, sys);
137 
138  else if (_solver->residual_and_jacobian_object != nullptr)
140 
141  else
142  return false;
143 
144  R.close();
145  X_global.close();
146 
147  return true;
148 }
149 
150 
151 
152 bool Problem_Interface::computeJacobian(const Epetra_Vector & x,
153  Epetra_Operator & jac)
154 {
155  LOG_SCOPE("jacobian()", "TrilinosNoxNonlinearSolver");
156 
158 
159  EpetraMatrix<Number> J(&dynamic_cast<Epetra_FECrsMatrix &>(jac), sys.comm());
160  EpetraVector<Number> & X_sys = *cast_ptr<EpetraVector<Number> *>(sys.solution.get());
161  EpetraVector<Number> X_global(*const_cast<Epetra_Vector *>(&x), sys.comm());
162 
163  // Set the dof maps
164  J.attach_dof_map(sys.get_dof_map());
165 
166  // Use the systems update() to get a good local version of the parallel solution
167  X_global.swap(X_sys);
168 
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
178  libmesh_error_msg_if(_solver->jacobian && _solver->jacobian_object,
179  "ERROR: cannot specify both a function and object to compute the Jacobian!");
180 
181  libmesh_error_msg_if(_solver->matvec && _solver->residual_and_jacobian_object,
182  "ERROR: cannot specify both a function and object to compute the combined Residual & Jacobian!");
183 
184  if (_solver->jacobian != nullptr)
185  _solver->jacobian(*sys.current_local_solution.get(), J, sys);
186 
187  else if (_solver->jacobian_object != nullptr)
189 
190  else if (_solver->matvec != nullptr)
191  _solver->matvec(*sys.current_local_solution.get(), nullptr, &J, sys);
192 
193  else if (_solver->residual_and_jacobian_object != nullptr)
195 
196  else
197  libmesh_error_msg("Error! Unable to compute residual and/or Jacobian!");
198 
199  J.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> J(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  J.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 
234  if (this->_exact_constraint_enforcement)
236  sys.update();
237 
238  X_global.swap(X_sys);
239 
240  //-----------------------------------------------------------------------------
241  // if the user has provided both function pointers and objects only the pointer
242  // will be used, so catch that as an error
243  libmesh_error_msg_if(_solver->jacobian && _solver->jacobian_object,
244  "ERROR: cannot specify both a function and object to compute the Jacobian!");
245 
246  libmesh_error_msg_if(_solver->matvec && _solver->residual_and_jacobian_object,
247  "ERROR: cannot specify both a function and object to compute the combined Residual & Jacobian!");
248 
249  if (_solver->jacobian != nullptr)
250  _solver->jacobian(*sys.current_local_solution.get(), J, sys);
251 
252  else if (_solver->jacobian_object != nullptr)
254 
255  else if (_solver->matvec != nullptr)
256  _solver->matvec(*sys.current_local_solution.get(), nullptr, &J, sys);
257 
258  else if (_solver->residual_and_jacobian_object != nullptr)
260 
261  else
262  libmesh_error_msg("Error! Unable to compute residual and/or Jacobian!");
263 
264  J.close();
265  X_global.close();
266 
267  tpc.compute();
268 
269  return true;
270 }
271 
272 
273 
274 //---------------------------------------------------------------------
275 // NoxNonlinearSolver<> methods
276 template <typename T>
278 {
279 }
280 
281 
282 
283 template <typename T>
284 void NoxNonlinearSolver<T>::init (const char * /*name*/)
285 {
286  if (!this->initialized())
287  _interface = new Problem_Interface(this);
288 }
289 
290 
291 
292 #include <libmesh/ignore_warnings.h> // deprecated-copy in Epetra_Vector
293 
294 template <typename T>
295 std::pair<unsigned int, Real>
296 NoxNonlinearSolver<T>::solve (SparseMatrix<T> & /* jac_in */, // System Jacobian Matrix
297  NumericVector<T> & x_in, // Solution vector
298  NumericVector<T> & /* r_in */, // Residual vector
299  const double, // Stopping tolerance
300  const unsigned int)
301 {
302  this->init ();
303 
304  if (this->user_presolve)
305  this->user_presolve(this->system());
306 
307  EpetraVector<T> * x_epetra = cast_ptr<EpetraVector<T> *>(&x_in);
308  // Creating a Teuchos::RCP as they do in NOX examples does not work here - we get some invalid memory references
309  // thus we make a local copy
310  NOX::Epetra::Vector x(*x_epetra->vec());
311 
312  Teuchos::RCP<Teuchos::ParameterList> nlParamsPtr = Teuchos::rcp(new Teuchos::ParameterList);
313  Teuchos::ParameterList & nlParams = *(nlParamsPtr.get());
314  nlParams.set("Nonlinear Solver", "Line Search Based");
315 
316  //print params
317  Teuchos::ParameterList & printParams = nlParams.sublist("Printing");
318  printParams.set("Output Precision", 3);
319  printParams.set("Output Processor", 0);
320  printParams.set("Output Information",
321  NOX::Utils::OuterIteration +
322  NOX::Utils::OuterIterationStatusTest +
323  NOX::Utils::InnerIteration +
324  NOX::Utils::LinearSolverDetails +
325  NOX::Utils::Parameters +
326  NOX::Utils::Details +
327  NOX::Utils::Warning);
328 
329  Teuchos::ParameterList & dirParams = nlParams.sublist("Direction");
330  dirParams.set("Method", "Newton");
331  Teuchos::ParameterList & newtonParams = dirParams.sublist("Newton");
332  newtonParams.set("Forcing Term Method", "Constant");
333 
334  Teuchos::ParameterList & lsParams = newtonParams.sublist("Linear Solver");
335  lsParams.set("Aztec Solver", "GMRES");
336  lsParams.set("Max Iterations", static_cast<int>(this->max_linear_iterations));
337  lsParams.set("Tolerance", this->initial_linear_tolerance);
338  lsParams.set("Output Frequency", 1);
339  lsParams.set("Size of Krylov Subspace", 1000);
340 
341  //create linear system
342  Teuchos::RCP<NOX::Epetra::Interface::Required> iReq(_interface);
343  Teuchos::RCP<NOX::Epetra::LinearSystemAztecOO> linSys;
344  Teuchos::RCP<Epetra_Operator> pc;
345 
346  if (this->jacobian || this->jacobian_object || this->residual_and_jacobian_object)
347  {
348  if (this->_preconditioner)
349  {
350  // PJNFK
351  lsParams.set("Preconditioner", "User Defined");
352 
353  TrilinosPreconditioner<Number> * trilinos_pc =
354  cast_ptr<TrilinosPreconditioner<Number> *>(this->_preconditioner);
355  pc = Teuchos::rcp(trilinos_pc);
356 
357  Teuchos::RCP<NOX::Epetra::Interface::Preconditioner> iPrec(_interface);
358  linSys = Teuchos::rcp(new NOX::Epetra::LinearSystemAztecOO(printParams, lsParams, iReq, iPrec, pc, x));
359  }
360  else
361  {
362  lsParams.set("Preconditioner", "None");
363  // lsParams.set("Preconditioner", "Ifpack");
364  // lsParams.set("Preconditioner", "AztecOO");
365 
366  // full jacobian
367  NonlinearImplicitSystem & sys = _interface->_solver->system();
368  EpetraMatrix<Number> & jacSys = *cast_ptr<EpetraMatrix<Number> *>(sys.matrix);
369  Teuchos::RCP<Epetra_RowMatrix> jacMat = Teuchos::rcp(jacSys.mat());
370 
371  Teuchos::RCP<NOX::Epetra::Interface::Jacobian> iJac(_interface);
372  linSys = Teuchos::rcp(new NOX::Epetra::LinearSystemAztecOO(printParams, lsParams, iReq, iJac, jacMat, x));
373  }
374  }
375  else
376  {
377  // matrix free
378  Teuchos::RCP<NOX::Epetra::MatrixFree> MF = Teuchos::rcp(new NOX::Epetra::MatrixFree(printParams, iReq, x));
379 
380  Teuchos::RCP<NOX::Epetra::Interface::Jacobian> iJac(MF);
381  linSys = Teuchos::rcp(new NOX::Epetra::LinearSystemAztecOO(printParams, lsParams, iReq, iJac, MF, x));
382  }
383 
384  //create group
385  Teuchos::RCP<NOX::Epetra::Group> grpPtr = Teuchos::rcp(new NOX::Epetra::Group(printParams, iReq, x, linSys));
386  NOX::Epetra::Group & grp = *(grpPtr.get());
387 
388  Teuchos::RCP<NOX::StatusTest::NormF> absresid =
389  Teuchos::rcp(new NOX::StatusTest::NormF(this->absolute_residual_tolerance, NOX::StatusTest::NormF::Unscaled));
390  Teuchos::RCP<NOX::StatusTest::NormF> relresid =
391  Teuchos::rcp(new NOX::StatusTest::NormF(grp, this->relative_residual_tolerance));
392  Teuchos::RCP<NOX::StatusTest::MaxIters> maxiters =
393  Teuchos::rcp(new NOX::StatusTest::MaxIters(this->max_nonlinear_iterations));
394  Teuchos::RCP<NOX::StatusTest::FiniteValue> finiteval =
395  Teuchos::rcp(new NOX::StatusTest::FiniteValue());
396  Teuchos::RCP<NOX::StatusTest::NormUpdate> normupdate =
397  Teuchos::rcp(new NOX::StatusTest::NormUpdate(this->absolute_step_tolerance));
398  Teuchos::RCP<NOX::StatusTest::Combo> combo =
399  Teuchos::rcp(new NOX::StatusTest::Combo(NOX::StatusTest::Combo::OR));
400  combo->addStatusTest(absresid);
401  combo->addStatusTest(relresid);
402  combo->addStatusTest(maxiters);
403  combo->addStatusTest(finiteval);
404  combo->addStatusTest(normupdate);
405 
406  Teuchos::RCP<Teuchos::ParameterList> finalPars = nlParamsPtr;
407 
408  Teuchos::RCP<NOX::Solver::Generic> solver = NOX::Solver::buildSolver(grpPtr, combo, nlParamsPtr);
409  NOX::StatusTest::StatusType status = solver->solve();
410  this->converged = (status == NOX::StatusTest::Converged);
411 
412  const NOX::Epetra::Group & finalGroup = dynamic_cast<const NOX::Epetra::Group &>(solver->getSolutionGroup());
413  const NOX::Epetra::Vector & noxFinalSln = dynamic_cast<const NOX::Epetra::Vector &>(finalGroup.getX());
414 
415  *x_epetra->vec() = noxFinalSln.getEpetraVector();
416  x_in.close();
417 
418  Real residual_norm = finalGroup.getNormF();
419  unsigned int total_iters = solver->getNumIterations();
420  _n_linear_iterations = finalPars->sublist("Direction").sublist("Newton").sublist("Linear Solver").sublist("Output").get("Total Number of Linear Iterations", -1);
421 
422  // do not let Trilinos to deallocate what we allocated
423  pc.release();
424  iReq.release();
425 
426  return std::make_pair(total_iters, residual_norm);
427 }
428 
429 #include <libmesh/restore_warnings.h>
430 
431 
432 
433 template <typename T>
434 int
436 {
437  return _n_linear_iterations;
438 }
439 
440 
441 
442 //------------------------------------------------------------------
443 // Explicit instantiations
444 template class LIBMESH_EXPORT NoxNonlinearSolver<Number>;
445 
446 } // namespace libMesh
447 
448 
449 
450 #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...
virtual int get_total_linear_iterations() override
Get the total number of linear iterations done in the last solve.
This class provides a nice interface to the Trilinos Epetra_Vector object.
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...
virtual void init(const char *name=nullptr) override
Initialize data structures if not done so already.
NumericVector< Number > * rhs
The system matrix.
Provides a uniform interface to vector storage schemes for different linear algebra libraries...
Definition: vector_fe_ex5.C:44
const Parallel::Communicator & comm() const
NoxNonlinearSolver< Number > * _solver
virtual std::pair< unsigned int, Real > solve(SparseMatrix< T > &, NumericVector< T > &, NumericVector< T > &, const double, const unsigned int) override
Call the Nox solver.
The libMesh namespace provides an interface to certain functionality in the library.
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.
virtual void clear() override
Release all memory and clear data structures.
Generic sparse matrix.
Definition: vector_fe_ex5.C:46
MPI_Status status
const sys_type & system() const
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)
bool _exact_constraint_enforcement
Whether we should enforce exact constraints globally during a solve.
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.
std::unique_ptr< NumericVector< Number > > solution
Data structure to hold solution values.
Definition: system.h:1593
void init(triangulateio &t)
Initializes the fields of t to nullptr/0 as necessary.
Manages consistently variables, degrees of freedom, coefficient vectors, matrices and non-linear solv...
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:510
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
SparseMatrix< Number > * matrix
The system matrix.
virtual void swap(NumericVector< T > &v) override
Swaps the contents of this with v.
std::unique_ptr< NumericVector< Number > > current_local_solution
All the values I need to compute my contribution to the simulation at hand.
Definition: system.h:1605
bool initialized()
Checks that library initialization has been done.
Definition: libmesh.C:276
bool computeJacobian(const Epetra_Vector &x, Epetra_Operator &Jac)
Compute an explicit Jacobian.
Epetra_FECrsMatrix * mat()
This class provides a nice interface to the Epetra data structures for parallel, sparse matrices...
const DofMap & get_dof_map() const
Definition: system.h:2374
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...
void enforce_constraints_exactly(const System &system, NumericVector< Number > *v=nullptr, bool homogeneous=false) const
Constrains the numeric vector v, which represents a solution defined on the mesh. ...
Definition: dof_map.h:2350
bool computeF(const Epetra_Vector &x, Epetra_Vector &FVec, NOX::Epetra::Interface::Required::FillType fillType)
Compute and return F.