Rythmos - Transient Integration for Differential Equations  Version of the Day
Rythmos_TimeStepNonlinearSolver_def.hpp
1 //@HEADER
2 // ***********************************************************************
3 //
4 // Rythmos Package
5 // Copyright (2006) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // This library is free software; you can redistribute it and/or modify
11 // it under the terms of the GNU Lesser General Public License as
12 // published by the Free Software Foundation; either version 2.1 of the
13 // License, or (at your option) any later version.
14 //
15 // This library is distributed in the hope that it will be useful, but
16 // WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 // Lesser General Public License for more details.
19 //
20 // You should have received a copy of the GNU Lesser General Public
21 // License along with this library; if not, write to the Free Software
22 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
23 // USA
24 // Questions? Contact Todd S. Coffey (tscoffe@sandia.gov)
25 //
26 // ***********************************************************************
27 //@HEADER
28 
29 
30 #ifndef RYTHMOS_TIME_STEP_NONLINEAR_SOLVER_DEF_HPP
31 #define RYTHMOS_TIME_STEP_NONLINEAR_SOLVER_DEF_HPP
32 
33 #include "Rythmos_TimeStepNonlinearSolver_decl.hpp"
34 
35 #include "Thyra_TestingTools.hpp"
36 #include "Thyra_ModelEvaluatorHelpers.hpp"
37 #include "Teuchos_VerboseObjectParameterListHelpers.hpp"
38 #include "Teuchos_StandardParameterEntryValidators.hpp"
39 #include "Teuchos_as.hpp"
40 
41 namespace Rythmos {
42 
43 
44 // ////////////////////////
45 // Defintions
46 
47 
48 // Static members
49 
50 
51 template<class Scalar>
52 const std::string
53 TimeStepNonlinearSolver<Scalar>::DefaultTol_name_ = "Default Tol";
54 
55 template<class Scalar>
56 const double
57 TimeStepNonlinearSolver<Scalar>::DefaultTol_default_ = 1e-2;
58 
59 
60 template<class Scalar>
61 const std::string
62 TimeStepNonlinearSolver<Scalar>::DefaultMaxIters_name_ = "Default Max Iters";
63 
64 template<class Scalar>
65 const int
66 TimeStepNonlinearSolver<Scalar>::DefaultMaxIters_default_ = 3;
67 
68 
69 template<class Scalar>
70 const std::string
71 TimeStepNonlinearSolver<Scalar>::NonlinearSafetyFactor_name_
72 = "Nonlinear Safety Factor";
73 
74 template<class Scalar>
75 const double
76 TimeStepNonlinearSolver<Scalar>::NonlinearSafetyFactor_default_ = 0.1;
77 
78 
79 template<class Scalar>
80 const std::string
81 TimeStepNonlinearSolver<Scalar>::LinearSafetyFactor_name_ = "Linear Safety Factor";
82 
83 template<class Scalar>
84 const double
85 TimeStepNonlinearSolver<Scalar>::LinearSafetyFactor_default_ = 0.05;
86 
87 
88 template<class Scalar>
89 const std::string
90 TimeStepNonlinearSolver<Scalar>::RMinFraction_name_ = "R Min Fraction";
91 
92 template<class Scalar>
93 const double
94 TimeStepNonlinearSolver<Scalar>::RMinFraction_default_ = 0.3;
95 
96 
97 template<class Scalar>
98 const std::string
99 TimeStepNonlinearSolver<Scalar>::ThrownOnLinearSolveFailure_name_
100 = "Thrown on Linear Solve Failure";
101 
102 template<class Scalar>
103 const bool
104 TimeStepNonlinearSolver<Scalar>::ThrownOnLinearSolveFailure_default_ = false;
105 
106 
107 // Constructors/Intializers/Misc
108 
109 
110 template <class Scalar>
112  :J_is_current_(false),
113  defaultTol_(DefaultTol_default_),
114  defaultMaxIters_(DefaultMaxIters_default_),
115  nonlinearSafetyFactor_(NonlinearSafetyFactor_default_),
116  linearSafetyFactor_(LinearSafetyFactor_default_),
117  RMinFraction_(RMinFraction_default_),
118  throwOnLinearSolveFailure_(ThrownOnLinearSolveFailure_default_)
119 {}
120 
121 
122 // Overridden from ParameterListAcceptor
123 
124 
125 template<class Scalar>
127  RCP<ParameterList> const& paramList
128  )
129 {
130  using Teuchos::get;
131  TEUCHOS_TEST_FOR_EXCEPT(is_null(paramList));
132  paramList->validateParametersAndSetDefaults(*getValidParameters(),0);
133  paramList_ = paramList;
134  defaultTol_ = get<double>(*paramList_,DefaultTol_name_);
135  defaultMaxIters_ = get<int>(*paramList_,DefaultMaxIters_name_);
136  nonlinearSafetyFactor_ = get<double>(*paramList_,NonlinearSafetyFactor_name_);
137  linearSafetyFactor_ = get<double>(*paramList_,LinearSafetyFactor_name_);
138  RMinFraction_ = get<double>(*paramList_,RMinFraction_name_);
139  throwOnLinearSolveFailure_ = get<bool>(
140  *paramList_,ThrownOnLinearSolveFailure_name_);
141  Teuchos::readVerboseObjectSublist(&*paramList_,this);
142 #ifdef HAVE_RYTHMOS_DEBUG
143  paramList_->validateParameters(*getValidParameters(),0);
144 #endif // HAVE_RYTHMOS_DEBUG
145 }
146 
147 
148 template<class Scalar>
149 RCP<ParameterList>
151 {
152  return paramList_;
153 }
154 
155 
156 template<class Scalar>
157 RCP<ParameterList>
159 {
160  RCP<ParameterList> _paramList = paramList_;
161  paramList_ = Teuchos::null;
162  return _paramList;
163 }
164 
165 
166 template<class Scalar>
167 RCP<const ParameterList>
169 {
170  return paramList_;
171 }
172 
173 
174 template<class Scalar>
175 RCP<const ParameterList>
177 {
178  using Teuchos::setDoubleParameter; using Teuchos::setIntParameter;
179  static RCP<const ParameterList> validPL;
180  if (is_null(validPL)) {
181  RCP<ParameterList> pl = Teuchos::parameterList();
182  setDoubleParameter(
183  DefaultTol_name_, DefaultTol_default_,
184  "The default base tolerance for the nonlinear timestep solve.\n"
185  "This tolerance can be overridden ???",
186  &*pl );
187  setIntParameter(
188  DefaultMaxIters_name_, DefaultMaxIters_default_,
189  "The default maximum number of Newton iterations to perform.\n"
190  "This default can be overridden ???",
191  &*pl );
192  setDoubleParameter(
193  NonlinearSafetyFactor_name_, NonlinearSafetyFactor_default_,
194  "The factor (< 1.0) to multiply tol to bound R*||dx|||.\n"
195  "The exact nonlinear convergence test is:\n"
196  " R*||dx|| <= \"" + NonlinearSafetyFactor_name_ + "\" * tol.",
197  &*pl );
198  setDoubleParameter(
199  LinearSafetyFactor_name_, LinearSafetyFactor_default_,
200  "This factor multiplies the nonlinear safety factor which multiplies\n"
201  "tol when determining the linear solve tolerence.\n"
202  "The exact linear convergence tolerance is:\n"
203  " ||J*dx+f||/||f|| <= \"" + LinearSafetyFactor_name_ + "\" * "
204  "\"" + NonlinearSafetyFactor_name_ + "\" * tol.",
205  &*pl );
206  setDoubleParameter(
207  RMinFraction_name_, RMinFraction_default_,
208  "The faction below which the R factor is not allowed to drop\n"
209  "below each Newton iteration. The R factor is related to the\n"
210  "ratio of ||dx||/||dx_last|| between nonlinear iterations.",
211  &*pl );
212  pl->set(
213  ThrownOnLinearSolveFailure_name_, ThrownOnLinearSolveFailure_default_,
214  "If set to true (\"1\"), then an Thyra::CatastrophicSolveFailure\n"
215  "exception will be thrown when a linear solve fails to meet it's tolerance."
216  );
217  Teuchos::setupVerboseObjectSublist(&*pl);
218  validPL = pl;
219  }
220  return validPL;
221 }
222 
223 
224 // Overridden from NonlinearSolverBase
225 
226 
227 template <class Scalar>
229  const RCP<const Thyra::ModelEvaluator<Scalar> > &model
230  )
231 {
232  TEUCHOS_TEST_FOR_EXCEPT(model.get()==NULL);
233  model_ = model;
234  J_ = Teuchos::null;
235  current_x_ = Teuchos::null;
236  J_is_current_ = false;
237 }
238 
239 
240 template <class Scalar>
241 RCP<const Thyra::ModelEvaluator<Scalar> >
243 {
244  return model_;
245 }
246 
247 template <class Scalar>
248 Thyra::SolveStatus<Scalar>
250  Thyra::VectorBase<Scalar> *x,
251  const Thyra::SolveCriteria<Scalar> *solveCriteria,
252  Thyra::VectorBase<Scalar> *delta
253  )
254 {
255 
256  RYTHMOS_FUNC_TIME_MONITOR("Rythmos:TimeStepNonlinearSolver::solve");
257 
258  using std::endl;
259  using Teuchos::incrVerbLevel;
260  using Teuchos::describe;
261  using Teuchos::as;
262  using Teuchos::rcp;
263  using Teuchos::OSTab;
264  using Teuchos::getFancyOStream;
265  typedef Thyra::ModelEvaluatorBase MEB;
266  typedef Teuchos::VerboseObjectTempState<MEB> VOTSME;
267  typedef Thyra::LinearOpWithSolveBase<Scalar> LOWSB;
268  typedef Teuchos::VerboseObjectTempState<LOWSB> VOTSLOWSB;
269 
270 #ifdef HAVE_RYTHMOS_DEBUG
271  TEUCHOS_TEST_FOR_EXCEPT(0==x);
272  THYRA_ASSERT_VEC_SPACES(
273  "TimeStepNonlinearSolver<Scalar>::solve(...)",
274  *x->space(),*model_->get_x_space() );
275  TEUCHOS_TEST_FOR_EXCEPT(
276  0!=solveCriteria && "ToDo: Support passed in solve criteria!" );
277 #endif
278 
279  const RCP<Teuchos::FancyOStream> out = this->getOStream();
280  const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
281  const bool showNewtonDetails =
282  (!is_null(out) && (as<int>(verbLevel) >= as<int>(Teuchos::VERB_MEDIUM)));
283  const bool dumpAll =
284  (!is_null(out) && (as<int>(verbLevel) == as<int>(Teuchos::VERB_EXTREME)));
285  TEUCHOS_OSTAB;
286  VOTSME stateModel_outputTempState(model_,out,incrVerbLevel(verbLevel,-1));
287 
288  if (showNewtonDetails)
289  *out
290  << "\nEntering TimeStepNonlinearSolver::solve(...) ...\n"
291  << "\nmodel = " << Teuchos::describe(*model_,verbLevel);
292 
293  if(dumpAll) {
294  *out << "\nInitial guess:\n";
295  *out << "\nx = " << *x;
296  }
297 
298  // Initialize storage for algorithm
299  if(!J_.get()) J_ = model_->create_W();
300  TEUCHOS_TEST_FOR_EXCEPTION( Teuchos::is_null(J_), std::logic_error,
301  "Error! model->create_W() returned a null pointer!\n"
302  );
303  RCP<Thyra::VectorBase<Scalar> > f = createMember(model_->get_f_space());
304  RCP<Thyra::VectorBase<Scalar> > dx = createMember(model_->get_x_space());
305  RCP<Thyra::VectorBase<Scalar> > dx_last = createMember(model_->get_x_space());
306  RCP<Thyra::VectorBase<Scalar> > x_curr = createMember(model_->get_x_space());
307  if (delta != NULL)
308  Thyra::V_S(ptr(delta),ST::zero()); // delta stores the cumulative update to x over the whole Newton solve.
309  Thyra::assign(x_curr.ptr(),*x);
310  J_is_current_ = false;
311  current_x_ = Teuchos::null;
312 
313  // Initialize convergence criteria
314  ScalarMag R = SMT::one();
315  ScalarMag linearTolSafety = linearSafetyFactor_ * nonlinearSafetyFactor_;
316  int maxIters = defaultMaxIters_;
317  ScalarMag tol = defaultTol_;
318  // ToDo: Get above from solveCriteria!
319 
320  // Do the undampened Newton iterations
321  bool converged = false;
322  bool sawFailedLinearSolve = false;
323  Thyra::SolveStatus<Scalar> failedLinearSolveStatus;
324  ScalarMag nrm_dx = SMT::nan();
325  ScalarMag nrm_dx_last = SMT::nan();
326  int iter = 1;
327  for( ; iter <= maxIters; ++iter ) {
328  if (showNewtonDetails)
329  *out << "\n*** newtonIter = " << iter << endl;
330  if (showNewtonDetails)
331  *out << "\nEvaluating the model f and W ...\n";
332  Thyra::eval_f_W( *model_, *x_curr, &*f, &*J_ );
333  if (showNewtonDetails)
334  *out << "\nSolving the system J*dx = -f ...\n";
335  Thyra::V_S(dx.ptr(),ST::zero()); // Initial guess is needed!
336  Thyra::SolveCriteria<Scalar>
337  linearSolveCriteria(
338  Thyra::SolveMeasureType(
339  Thyra::SOLVE_MEASURE_NORM_RESIDUAL, Thyra::SOLVE_MEASURE_NORM_RHS
340  ),
341  linearTolSafety*tol
342  );
343  VOTSLOWSB J_outputTempState(J_,out,incrVerbLevel(verbLevel,-1));
344  Thyra::SolveStatus<Scalar> linearSolveStatus
345  = J_->solve(Thyra::NOTRANS, *f, dx.ptr(), Teuchos::ptr(&linearSolveCriteria) );
346  if (showNewtonDetails)
347  *out << "\nLinear solve status:\n" << linearSolveStatus;
348  Thyra::Vt_S(dx.ptr(),Scalar(-ST::one()));
349  if (dumpAll)
350  *out << "\ndx = " << Teuchos::describe(*dx,verbLevel);
351  if (delta != NULL) {
352  Thyra::Vp_V(ptr(delta),*dx);
353  if (dumpAll)
354  *out << "\ndelta = " << Teuchos::describe(*delta,verbLevel);
355  }
356  // Check the linear solve
357  if(linearSolveStatus.solveStatus != Thyra::SOLVE_STATUS_CONVERGED) {
358  sawFailedLinearSolve = true;
359  failedLinearSolveStatus = linearSolveStatus;
360  if (throwOnLinearSolveFailure_) {
361  TEUCHOS_TEST_FOR_EXCEPTION(
362  throwOnLinearSolveFailure_, Thyra::CatastrophicSolveFailure,
363  "Error, the linear solver did not converge!"
364  );
365  }
366  if (showNewtonDetails)
367  *out << "\nWarning, linear solve did not converge! Continuing anyway :-)\n";
368  }
369  // Update the solution: x_curr = x_curr + dx
370  Vp_V( x_curr.ptr(), *dx );
371  if (dumpAll)
372  *out << "\nUpdated solution x = " << Teuchos::describe(*x_curr,verbLevel);
373  // Convergence test
374  nrm_dx = Thyra::norm(*dx);
375  if ( R*nrm_dx <= nonlinearSafetyFactor_*tol )
376  converged = true;
377  if (showNewtonDetails)
378  *out
379  << "\nConvergence test:\n"
380  << " R*||dx|| = " << R << "*" << nrm_dx
381  << " = " << (R*nrm_dx) << "\n"
382  << " <= nonlinearSafetyFactor*tol = " << nonlinearSafetyFactor_ << "*" << tol
383  << " = " << (nonlinearSafetyFactor_*tol)
384  << " : " << ( converged ? "converged!" : " unconverged" )
385  << endl;
386  if(converged)
387  break; // We have converged!!!
388  // Update convergence criteria for the next iteration ...
389  if(iter > 1) {
390  const Scalar
391  MinR = RMinFraction_*R,
392  nrm_dx_ratio = nrm_dx/nrm_dx_last;
393  R = std::max(MinR,nrm_dx_ratio);
394  if (showNewtonDetails)
395  *out
396  << "\nUpdated R\n"
397  << " = max(RMinFraction*R,||dx||/||dx_last||)\n"
398  << " = max("<<RMinFraction_<<"*"<<R<<","<<nrm_dx<<"/"<<nrm_dx_last<<")\n"
399  << " = max("<<MinR<<","<<nrm_dx_ratio<<")\n"
400  << " = " << R << endl;
401  }
402  // Save to old
403  std::swap(dx_last,dx);
404  nrm_dx_last = nrm_dx;
405  }
406 
407  // Set the solution
408  Thyra::assign(ptr(x),*x_curr);
409 
410  if (dumpAll)
411  *out << "\nFinal solution x = " << Teuchos::describe(*x,verbLevel);
412 
413  // Check the status
414 
415  Thyra::SolveStatus<Scalar> solveStatus;
416 
417  std::ostringstream oss;
418  Teuchos::FancyOStream omsg(rcp(&oss,false));
419 
420  omsg << "Solver: " << this->description() << endl;
421 
422  if(converged) {
423  solveStatus.solveStatus = Thyra::SOLVE_STATUS_CONVERGED;
424  omsg << "CVODE status test converged!\n";
425  }
426  else {
427  solveStatus.solveStatus = Thyra::SOLVE_STATUS_UNCONVERGED;
428  omsg << "CVODE status test failed!\n";
429  }
430 
431  if (sawFailedLinearSolve) {
432  omsg << "Warning! A failed linear solve was encountered with status:\n";
433  OSTab tab(omsg);
434  omsg << failedLinearSolveStatus;
435  }
436 
437  omsg
438  << "R*||dx|| = " << R << "*" << nrm_dx
439  << " <= nonlinearSafetyFactor*tol = " << nonlinearSafetyFactor_ << "*" << tol << " : "
440  << ( converged ? "converged!" : " unconverged" ) << endl;
441 
442  omsg
443  << "Iterations = " << iter;
444  // Above, we leave off the last newline since this is the convention for the
445  // SolveStatus::message string!
446 
447  solveStatus.message = oss.str();
448 
449  // Update the solution state for external clients
450  current_x_ = x->clone_v();
451  J_is_current_ = false;
452  // 2007/09/04: rabartl: Note, above the Jacobian J is always going to be out
453  // of date since this algorithm computes x_curr = x_curr + dx for at least
454  // one solve for dx = -inv(J)*f. Therefore, J is never at the updated
455  // x_curr, only the old x_curr!
456 
457  if (showNewtonDetails)
458  *out << "\nLeaving TimeStepNonlinearSolver::solve(...) ...\n";
459 
460  return solveStatus;
461 
462 }
463 
464 
465 template <class Scalar>
467 {
468  return true;
469 }
470 
471 
472 template <class Scalar>
473 RCP<Thyra::NonlinearSolverBase<Scalar> >
475 {
476  RCP<TimeStepNonlinearSolver<Scalar> >
477  nonlinearSolver = Teuchos::rcp(new TimeStepNonlinearSolver<Scalar>);
478  nonlinearSolver->model_ = model_; // Shallow copy is okay, model is stateless
479  nonlinearSolver->defaultTol_ = defaultTol_;
480  nonlinearSolver->defaultMaxIters_ = defaultMaxIters_;
481  nonlinearSolver->nonlinearSafetyFactor_ = nonlinearSafetyFactor_;
482  nonlinearSolver->linearSafetyFactor_ = linearSafetyFactor_;
483  nonlinearSolver->RMinFraction_ = RMinFraction_;
484  nonlinearSolver->throwOnLinearSolveFailure_ = throwOnLinearSolveFailure_;
485  // Note: The specification of this virtual function in the interface class
486  // allows us to just copy the algorithm, not the entire state so we are
487  // done!
488  return nonlinearSolver;
489 }
490 
491 
492 template <class Scalar>
493 RCP<const Thyra::VectorBase<Scalar> >
495 {
496  return current_x_;
497 }
498 
499 
500 template <class Scalar>
502 {
503  return J_is_current_;
504 }
505 
506 
507 template <class Scalar>
508 RCP<Thyra::LinearOpWithSolveBase<Scalar> >
510 {
511  if (is_null(J_))
512  return Teuchos::null;
513  if (forceUpToDate) {
514 #ifdef HAVE_RYTHMOS_DEBUG
515  TEUCHOS_TEST_FOR_EXCEPT(is_null(current_x_));
516 #endif
517  Thyra::eval_f_W<Scalar>( *model_, *current_x_, 0, &*J_ );
518  J_is_current_ = true;
519  }
520  return J_;
521 }
522 
523 
524 template <class Scalar>
525 RCP<const Thyra::LinearOpWithSolveBase<Scalar> >
527 {
528  return J_;
529 }
530 
531 
532 template <class Scalar>
534 {
535  J_is_current_ = W_is_current;
536 }
537 
538 
539 } // namespace Rythmos
540 
541 
542 // Nonmember constructors
543 
544 
545 template <class Scalar>
546 Teuchos::RCP<Rythmos::TimeStepNonlinearSolver<Scalar> >
547 Rythmos::timeStepNonlinearSolver()
548 {
549  return Teuchos::rcp(new TimeStepNonlinearSolver<Scalar>);
550 }
551 
552 
553 template <class Scalar>
554 Teuchos::RCP<Rythmos::TimeStepNonlinearSolver<Scalar> >
555 Rythmos::timeStepNonlinearSolver(const RCP<ParameterList> &pl)
556 {
557  const RCP<Rythmos::TimeStepNonlinearSolver<Scalar> >
558  solver = timeStepNonlinearSolver<Scalar>();
559  solver->setParameterList(pl);
560  return solver;
561 }
562 
563 
564 //
565 // Explicit Instantiation macro
566 //
567 // Must be expanded from within the Rythmos namespace!
568 //
569 
570 #define RYTHMOS_TIME_STEP_NONLINEAR_SOLVER_INSTANT(SCALAR) \
571  \
572  template class TimeStepNonlinearSolver< SCALAR >; \
573  \
574  template RCP<TimeStepNonlinearSolver< SCALAR > > timeStepNonlinearSolver(); \
575  \
576  template RCP<TimeStepNonlinearSolver<SCALAR > > \
577  timeStepNonlinearSolver(const RCP<ParameterList> &pl);
578 
579 
580 
581 #endif // RYTHMOS_TIME_STEP_NONLINEAR_SOLVER_DEF_HPP
Simple undampended Newton solver designed to solve time step equations in accurate times-tepping meth...
RCP< const Thyra::LinearOpWithSolveBase< Scalar > > get_W() const
void setModel(const RCP< const Thyra::ModelEvaluator< Scalar > > &model)
RCP< const Thyra::VectorBase< Scalar > > get_current_x() const
RCP< const Thyra::ModelEvaluator< Scalar > > getModel() const
RCP< Thyra::NonlinearSolverBase< Scalar > > cloneNonlinearSolver() const
Thyra::SolveStatus< Scalar > solve(Thyra::VectorBase< Scalar > *x, const Thyra::SolveCriteria< Scalar > *solveCriteria, Thyra::VectorBase< Scalar > *delta=NULL)
RCP< const ParameterList > getParameterList() const
RCP< Thyra::LinearOpWithSolveBase< Scalar > > get_nonconst_W(const bool forceUpToDate)
RCP< const ParameterList > getValidParameters() const
void setParameterList(RCP< ParameterList > const &paramList)