20 #include "libmesh/libmesh_common.h" 22 #if defined(LIBMESH_TRILINOS_HAVE_NOX) && defined(LIBMESH_TRILINOS_HAVE_EPETRA) 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" 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" 46 public NOX::Epetra::Interface::Jacobian,
47 public NOX::Epetra::Interface::Preconditioner
56 bool computeF(
const Epetra_Vector & x,
58 NOX::Epetra::Interface::Required::FillType fillType);
62 Epetra_Operator & Jac);
71 Epetra_RowMatrix & M);
76 Epetra_Operator & Prec,
77 Teuchos::ParameterList * p);
96 NOX::Epetra::Interface::Required::FillType )
98 LOG_SCOPE(
"residual()",
"TrilinosNoxNonlinearSolver");
107 X_global.
swap(X_sys);
115 X_global.swap(X_sys);
127 "ERROR: cannot specify both a function and object to compute the combined Residual & Jacobian!");
153 Epetra_Operator & jac)
155 LOG_SCOPE(
"jacobian()",
"TrilinosNoxNonlinearSolver");
167 X_global.
swap(X_sys);
173 X_global.swap(X_sys);
179 "ERROR: cannot specify both a function and object to compute the Jacobian!");
182 "ERROR: cannot specify both a function and object to compute the combined Residual & Jacobian!");
197 libmesh_error_msg(
"Error! Unable to compute residual and/or Jacobian!");
216 Epetra_Operator & prec,
217 Teuchos::ParameterList * )
219 LOG_SCOPE(
"preconditioner()",
"TrilinosNoxNonlinearSolver");
232 X_global.
swap(X_sys);
234 if (this->_exact_constraint_enforcement)
238 X_global.swap(X_sys);
244 "ERROR: cannot specify both a function and object to compute the Jacobian!");
247 "ERROR: cannot specify both a function and object to compute the combined Residual & Jacobian!");
262 libmesh_error_msg(
"Error! Unable to compute residual and/or Jacobian!");
276 template <
typename T>
283 template <
typename T>
292 #include <libmesh/ignore_warnings.h> 294 template <
typename T>
295 std::pair<unsigned int, Real>
304 if (this->user_presolve)
305 this->user_presolve(this->system());
310 NOX::Epetra::Vector x(*x_epetra->
vec());
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");
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);
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");
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);
342 Teuchos::RCP<NOX::Epetra::Interface::Required> iReq(_interface);
343 Teuchos::RCP<NOX::Epetra::LinearSystemAztecOO> linSys;
344 Teuchos::RCP<Epetra_Operator> pc;
346 if (this->jacobian || this->jacobian_object || this->residual_and_jacobian_object)
348 if (this->_preconditioner)
351 lsParams.set(
"Preconditioner",
"User Defined");
354 cast_ptr<TrilinosPreconditioner<Number> *>(this->_preconditioner);
355 pc = Teuchos::rcp(trilinos_pc);
357 Teuchos::RCP<NOX::Epetra::Interface::Preconditioner> iPrec(_interface);
358 linSys = Teuchos::rcp(
new NOX::Epetra::LinearSystemAztecOO(printParams, lsParams, iReq, iPrec, pc, x));
362 lsParams.set(
"Preconditioner",
"None");
369 Teuchos::RCP<Epetra_RowMatrix> jacMat = Teuchos::rcp(jacSys.
mat());
371 Teuchos::RCP<NOX::Epetra::Interface::Jacobian> iJac(_interface);
372 linSys = Teuchos::rcp(
new NOX::Epetra::LinearSystemAztecOO(printParams, lsParams, iReq, iJac, jacMat, x));
378 Teuchos::RCP<NOX::Epetra::MatrixFree> MF = Teuchos::rcp(
new NOX::Epetra::MatrixFree(printParams, iReq, x));
380 Teuchos::RCP<NOX::Epetra::Interface::Jacobian> iJac(MF);
381 linSys = Teuchos::rcp(
new NOX::Epetra::LinearSystemAztecOO(printParams, lsParams, iReq, iJac, MF, x));
385 Teuchos::RCP<NOX::Epetra::Group> grpPtr = Teuchos::rcp(
new NOX::Epetra::Group(printParams, iReq, x, linSys));
386 NOX::Epetra::Group & grp = *(grpPtr.get());
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);
406 Teuchos::RCP<Teuchos::ParameterList> finalPars = nlParamsPtr;
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);
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());
415 *x_epetra->
vec() = noxFinalSln.getEpetraVector();
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);
426 return std::make_pair(total_iters, residual_norm);
429 #include <libmesh/restore_warnings.h> 433 template <
typename T>
437 return _n_linear_iterations;
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...
Epetra_FECrsMatrix * mat()
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...
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.
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.
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'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.
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.
bool initialized()
Checks that library initialization has been done.
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
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. ...
bool computeF(const Epetra_Vector &x, Epetra_Vector &FVec, NOX::Epetra::Interface::Required::FillType fillType)
Compute and return F.