Rythmos_TimeStepNonlinearSolver_def.hpp

00001 //@HEADER
00002 // ***********************************************************************
00003 //
00004 //                           Rythmos Package
00005 //                 Copyright (2006) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Todd S. Coffey (tscoffe@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 //@HEADER
00028 
00029 
00030 #ifndef RYTHMOS_TIME_STEP_NONLINEAR_SOLVER_DEF_HPP
00031 #define RYTHMOS_TIME_STEP_NONLINEAR_SOLVER_DEF_HPP
00032 
00033 #include "Rythmos_TimeStepNonlinearSolver_decl.hpp"
00034 
00035 #include "Thyra_TestingTools.hpp"
00036 #include "Thyra_ModelEvaluatorHelpers.hpp"
00037 #include "Teuchos_VerboseObjectParameterListHelpers.hpp"
00038 #include "Teuchos_StandardParameterEntryValidators.hpp"
00039 #include "Teuchos_as.hpp"
00040 
00041 namespace Rythmos {
00042 
00043 
00044 template <class Scalar>
00045 RCP<TimeStepNonlinearSolver<Scalar> > timeStepNonlinearSolver()
00046 {
00047   return Teuchos::rcp(new TimeStepNonlinearSolver<Scalar>);
00048 }
00049 
00050 
00051 // ////////////////////////
00052 // Defintions
00053 
00054 
00055 // Static members
00056 
00057 
00058 template<class Scalar>
00059 const std::string
00060 TimeStepNonlinearSolver<Scalar>::DefaultTol_name_ = "Default Tol";
00061 
00062 template<class Scalar>
00063 const double
00064 TimeStepNonlinearSolver<Scalar>::DefaultTol_default_ = 1e-2;
00065 
00066 
00067 template<class Scalar>
00068 const std::string
00069 TimeStepNonlinearSolver<Scalar>::DefaultMaxIters_name_ = "Default Max Iters";
00070 
00071 template<class Scalar>
00072 const int
00073 TimeStepNonlinearSolver<Scalar>::DefaultMaxIters_default_ = 3;
00074 
00075 
00076 template<class Scalar>
00077 const std::string
00078 TimeStepNonlinearSolver<Scalar>::NonlinearSafetyFactor_name_
00079 = "Nonlinear Safety Factor";
00080 
00081 template<class Scalar>
00082 const double
00083 TimeStepNonlinearSolver<Scalar>::NonlinearSafetyFactor_default_ = 0.1;
00084 
00085 
00086 template<class Scalar>
00087 const std::string
00088 TimeStepNonlinearSolver<Scalar>::LinearSafetyFactor_name_ = "Linear Safety Factor";
00089 
00090 template<class Scalar>
00091 const double
00092 TimeStepNonlinearSolver<Scalar>::LinearSafetyFactor_default_ = 0.05;
00093 
00094 
00095 template<class Scalar>
00096 const std::string
00097 TimeStepNonlinearSolver<Scalar>::RMinFraction_name_ = "R Min Fraction";
00098 
00099 template<class Scalar>
00100 const double
00101 TimeStepNonlinearSolver<Scalar>::RMinFraction_default_ = 0.3;
00102 
00103 
00104 template<class Scalar>
00105 const std::string
00106 TimeStepNonlinearSolver<Scalar>::ThrownOnLinearSolveFailure_name_
00107 = "Thrown on Linear Solve Failure";
00108 
00109 template<class Scalar>
00110 const bool
00111 TimeStepNonlinearSolver<Scalar>::ThrownOnLinearSolveFailure_default_ = false;
00112 
00113 
00114 // Constructors/Intializers/Misc
00115 
00116 
00117 template <class Scalar>
00118 TimeStepNonlinearSolver<Scalar>::TimeStepNonlinearSolver()
00119   :J_is_current_(false),
00120    defaultTol_(DefaultTol_default_),
00121    defaultMaxIters_(DefaultMaxIters_default_),
00122    nonlinearSafetyFactor_(NonlinearSafetyFactor_default_),
00123    linearSafetyFactor_(LinearSafetyFactor_default_),
00124    RMinFraction_(RMinFraction_default_),
00125    throwOnLinearSolveFailure_(ThrownOnLinearSolveFailure_default_)
00126 {}
00127 
00128 
00129 // Overridden from ParameterListAcceptor
00130 
00131 
00132 template<class Scalar>
00133 void TimeStepNonlinearSolver<Scalar>::setParameterList(
00134   RCP<ParameterList> const& paramList
00135   )
00136 {
00137   using Teuchos::get;
00138   TEST_FOR_EXCEPT(is_null(paramList));
00139   paramList->validateParametersAndSetDefaults(*getValidParameters(),0);
00140   paramList_ = paramList;
00141   defaultTol_ = get<double>(*paramList_,DefaultTol_name_);
00142   defaultMaxIters_ = get<int>(*paramList_,DefaultMaxIters_name_);
00143   nonlinearSafetyFactor_ = get<double>(*paramList_,NonlinearSafetyFactor_name_);
00144   linearSafetyFactor_ = get<double>(*paramList_,LinearSafetyFactor_name_);
00145   RMinFraction_ = get<double>(*paramList_,RMinFraction_name_);
00146   throwOnLinearSolveFailure_ = get<bool>(
00147     *paramList_,ThrownOnLinearSolveFailure_name_);
00148   Teuchos::readVerboseObjectSublist(&*paramList_,this);
00149 #ifdef RYTHMOS_DEBUG
00150   paramList_->validateParameters(*getValidParameters(),0);
00151 #endif // RYTHMOS_DEBUG
00152 }
00153 
00154 
00155 template<class Scalar>
00156 RCP<ParameterList>
00157 TimeStepNonlinearSolver<Scalar>::getNonconstParameterList()
00158 {
00159   return paramList_;
00160 }
00161 
00162 
00163 template<class Scalar>
00164 RCP<ParameterList>
00165 TimeStepNonlinearSolver<Scalar>::unsetParameterList()
00166 {
00167   RCP<ParameterList> _paramList = paramList_;
00168   paramList_ = Teuchos::null;
00169   return _paramList;
00170 }
00171 
00172 
00173 template<class Scalar>
00174 RCP<const ParameterList>
00175 TimeStepNonlinearSolver<Scalar>::getParameterList() const
00176 {
00177   return paramList_;
00178 }
00179 
00180 
00181 template<class Scalar>
00182 RCP<const ParameterList>
00183 TimeStepNonlinearSolver<Scalar>::getValidParameters() const
00184 {
00185   using Teuchos::setDoubleParameter; using Teuchos::setIntParameter;
00186   static RCP<const ParameterList> validPL;
00187   if (is_null(validPL)) {
00188     RCP<ParameterList> pl = Teuchos::parameterList();
00189     setDoubleParameter(
00190       DefaultTol_name_, DefaultTol_default_,
00191       "The default base tolerance for the nonlinear timestep solve.\n"
00192       "This tolerance can be overridden ???",
00193       &*pl );
00194     setIntParameter(
00195       DefaultMaxIters_name_, DefaultMaxIters_default_,
00196       "The default maximum number of Newton iterations to perform.\n"
00197       "This default can be overridden ???",
00198       &*pl );
00199     setDoubleParameter(
00200       NonlinearSafetyFactor_name_, NonlinearSafetyFactor_default_,
00201       "The factor (< 1.0) to multiply tol to bound R*||dx|||.\n"
00202       "The exact nonlinear convergence test is:\n"
00203       "  R*||dx|| <= \"" + NonlinearSafetyFactor_name_ + "\" * tol.",
00204       &*pl );
00205     setDoubleParameter(
00206       LinearSafetyFactor_name_, LinearSafetyFactor_default_,
00207       "This factor multiplies the nonlinear safety factor which multiplies\n"
00208       "tol when determining the linear solve tolerence.\n"
00209       "The exact linear convergence tolerance is:\n"
00210       "  ||J*dx+f||/||f|| <= \"" + LinearSafetyFactor_name_ + "\" * "
00211       "\"" + NonlinearSafetyFactor_name_ + "\" * tol.",
00212       &*pl );
00213     setDoubleParameter(
00214       RMinFraction_name_, RMinFraction_default_,
00215       "The faction below which the R factor is not allowed to drop\n"
00216       "below each Newton iteration.  The R factor is related to the\n"
00217       "ratio of ||dx||/||dx_last|| between nonlinear iterations.",
00218       &*pl );
00219     pl->set(
00220       ThrownOnLinearSolveFailure_name_, ThrownOnLinearSolveFailure_default_,
00221       "If set to true (\"1\"), then an Thyra::CatastrophicSolveFailure\n"
00222       "exception will be thrown when a linear solve fails to meet it's tolerance."
00223       );
00224     Teuchos::setupVerboseObjectSublist(&*pl);
00225     validPL = pl;
00226   }
00227   return validPL;
00228 }
00229 
00230 
00231 // Overridden from NonlinearSolverBase
00232 
00233 
00234 template <class Scalar>
00235 void TimeStepNonlinearSolver<Scalar>::setModel(
00236   const RCP<const Thyra::ModelEvaluator<Scalar> > &model
00237   )
00238 {
00239   TEST_FOR_EXCEPT(model.get()==NULL);
00240   model_ = model;
00241   J_ = Teuchos::null;
00242   current_x_ = Teuchos::null;
00243   J_is_current_ = false;
00244 }
00245 
00246 
00247 template <class Scalar>
00248 RCP<const Thyra::ModelEvaluator<Scalar> >
00249 TimeStepNonlinearSolver<Scalar>::getModel() const
00250 {
00251   return model_;
00252 }
00253 
00254 template <class Scalar>
00255 Thyra::SolveStatus<Scalar>
00256 TimeStepNonlinearSolver<Scalar>::solve(
00257   Thyra::VectorBase<Scalar> *x,
00258   const Thyra::SolveCriteria<Scalar> *solveCriteria,
00259   Thyra::VectorBase<Scalar> *delta
00260   )
00261 {
00262   
00263 #ifdef ENABLE_RYTHMOS_TIMERS
00264   TEUCHOS_FUNC_TIME_MONITOR("Rythmos:TimeStepNonlinearSolver::solve");
00265 #endif
00266 
00267   using std::endl;
00268   using Teuchos::incrVerbLevel;
00269   using Teuchos::describe;
00270   using Teuchos::as;
00271   using Teuchos::rcp;
00272   using Teuchos::OSTab;
00273   using Teuchos::getFancyOStream;
00274   typedef Thyra::ModelEvaluatorBase MEB;
00275   typedef Teuchos::VerboseObjectTempState<MEB> VOTSME;
00276   typedef Thyra::LinearOpWithSolveBase<Scalar> LOWSB;
00277   typedef Teuchos::VerboseObjectTempState<LOWSB> VOTSLOWSB;
00278 
00279 #ifdef RYTHMOS_DEBUG
00280   TEST_FOR_EXCEPT(0==x);
00281   THYRA_ASSERT_VEC_SPACES(
00282     "TimeStepNonlinearSolver<Scalar>::solve(...)",
00283     *x->space(),*model_->get_x_space() );
00284   TEST_FOR_EXCEPT(
00285     0!=solveCriteria && "ToDo: Support passed in solve criteria!" );
00286 #endif
00287 
00288   const RCP<Teuchos::FancyOStream> out = this->getOStream();
00289   const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00290   const bool showNewtonDetails =
00291     (!is_null(out) && (as<int>(verbLevel) >= as<int>(Teuchos::VERB_MEDIUM)));
00292   const bool dumpAll =
00293     (!is_null(out) && (as<int>(verbLevel) == as<int>(Teuchos::VERB_EXTREME))); 
00294   TEUCHOS_OSTAB;
00295   VOTSME stateModel_outputTempState(model_,out,incrVerbLevel(verbLevel,-1));
00296 
00297   if (showNewtonDetails)
00298     *out
00299       << "\nEntering TimeStepNonlinearSolver::solve(...) ...\n"
00300       << "\nmodel = " << Teuchos::describe(*model_,verbLevel);
00301 
00302   if(dumpAll) {
00303     *out << "\nInitial guess:\n";
00304     *out << "\nx = " << *x;
00305   }
00306 
00307   // Initialize storage for algorithm
00308   if(!J_.get()) J_ = model_->create_W();
00309   TEST_FOR_EXCEPTION( Teuchos::is_null(J_), std::logic_error,
00310       "Error!  model->create_W() returned a null pointer!\n"
00311       );
00312   RCP<Thyra::VectorBase<Scalar> > f = createMember(model_->get_f_space());
00313   RCP<Thyra::VectorBase<Scalar> > dx = createMember(model_->get_x_space());
00314   RCP<Thyra::VectorBase<Scalar> > dx_last = createMember(model_->get_x_space());
00315   RCP<Thyra::VectorBase<Scalar> > x_curr = createMember(model_->get_x_space());
00316   if (delta != NULL)
00317     Thyra::V_S(delta,ST::zero()); // delta stores the cumulative update to x over the whole Newton solve.
00318   Thyra::assign(&*x_curr,*x);
00319   J_is_current_ = false;
00320   current_x_ = Teuchos::null;
00321 
00322   // Initialize convergence criteria
00323   ScalarMag R = SMT::one();
00324   ScalarMag linearTolSafety = linearSafetyFactor_ * nonlinearSafetyFactor_;
00325   int maxIters = defaultMaxIters_;
00326   ScalarMag tol = defaultTol_;
00327   // ToDo: Get above from solveCriteria!
00328 
00329   // Do the undampened Newton iterations
00330   bool converged = false;
00331   bool sawFailedLinearSolve = false;
00332   Thyra::SolveStatus<Scalar> failedLinearSolveStatus;
00333   ScalarMag nrm_dx = SMT::nan();
00334   ScalarMag nrm_dx_last = SMT::nan();
00335   int iter = 1;
00336   for( ; iter <= maxIters; ++iter ) {
00337     if (showNewtonDetails)
00338       *out << "\n*** newtonIter = " << iter << endl;
00339     if (showNewtonDetails)
00340       *out << "\nEvaluating the model f and W ...\n";
00341     Thyra::eval_f_W( *model_, *x_curr, &*f, &*J_ );
00342     if (showNewtonDetails)
00343       *out << "\nSolving the system J*dx = -f ...\n";
00344     Thyra::V_S(&*dx,ST::zero()); // Initial guess is needed!
00345     Thyra::SolveCriteria<Scalar>
00346       linearSolveCriteria(
00347         Thyra::SolveMeasureType(
00348           Thyra::SOLVE_MEASURE_NORM_RESIDUAL, Thyra::SOLVE_MEASURE_NORM_RHS
00349           ),
00350         linearTolSafety*tol
00351         );
00352     VOTSLOWSB J_outputTempState(J_,out,incrVerbLevel(verbLevel,-1));
00353     Thyra::SolveStatus<Scalar> linearSolveStatus
00354       = Thyra::solve( *J_, Thyra::NOTRANS, *f, &*dx, &linearSolveCriteria );
00355     if (showNewtonDetails)
00356       *out << "\nLinear solve status:\n" << linearSolveStatus;
00357     Thyra::Vt_S(&*dx,Scalar(-ST::one()));
00358     if (dumpAll)
00359       *out << "\ndx = " << Teuchos::describe(*dx,verbLevel);
00360     if (delta != NULL) {
00361       Thyra::Vp_V(delta,*dx);
00362       if (dumpAll)
00363         *out << "\ndelta = " << Teuchos::describe(*delta,verbLevel);
00364     }
00365     // Check the linear solve
00366     if(linearSolveStatus.solveStatus != Thyra::SOLVE_STATUS_CONVERGED) {
00367       sawFailedLinearSolve = true;
00368       failedLinearSolveStatus = linearSolveStatus;
00369       if (throwOnLinearSolveFailure_) {
00370         TEST_FOR_EXCEPTION(
00371           throwOnLinearSolveFailure_, Thyra::CatastrophicSolveFailure,
00372           "Error, the linear solver did not converge!"
00373           );
00374       }
00375       if (showNewtonDetails)
00376         *out << "\nWarning, linear solve did not converge!  Continuing anyway :-)\n";
00377     }
00378     // Update the solution: x_curr = x_curr + dx
00379     Vp_V( &*x_curr, *dx );
00380     if (dumpAll)
00381       *out << "\nUpdated solution x = " << Teuchos::describe(*x_curr,verbLevel);
00382     // Convergence test
00383     nrm_dx = Thyra::norm(*dx);
00384     if ( R*nrm_dx <= nonlinearSafetyFactor_*tol )
00385       converged = true;
00386     if (showNewtonDetails)
00387       *out
00388         << "\nConvergence test:\n"
00389         << "  R*||dx|| = " << R << "*" << nrm_dx
00390         << " = " << (R*nrm_dx) << "\n"
00391         << "    <= nonlinearSafetyFactor*tol = " << nonlinearSafetyFactor_ << "*" << tol
00392         << " = " << (nonlinearSafetyFactor_*tol)
00393         << " : " << ( converged ? "converged!" : " unconverged" )
00394         << endl;
00395     if(converged)
00396       break; // We have converged!!!
00397     // Update convergence criteria for the next iteration ...
00398     if(iter > 1) {
00399       const Scalar
00400         MinR = RMinFraction_*R,
00401         nrm_dx_ratio = nrm_dx/nrm_dx_last;
00402       R = std::max(MinR,nrm_dx_ratio);
00403       if (showNewtonDetails)
00404       *out
00405         << "\nUpdated R\n"
00406         << "  = max(RMinFraction*R,||dx||/||dx_last||)\n"
00407         << "  = max("<<RMinFraction_<<"*"<<R<<","<<nrm_dx<<"/"<<nrm_dx_last<<")\n"
00408         << "  = max("<<MinR<<","<<nrm_dx_ratio<<")\n"
00409         << "  = " << R << endl;
00410     }
00411     // Save to old
00412     std::swap(dx_last,dx);
00413     nrm_dx_last = nrm_dx;
00414   }
00415 
00416   // Set the solution
00417   Thyra::assign(x,*x_curr);
00418   
00419   if (dumpAll)
00420     *out << "\nFinal solution x = " << Teuchos::describe(*x,verbLevel);
00421 
00422   // Check the status
00423 
00424   Thyra::SolveStatus<Scalar> solveStatus;
00425 
00426   std::ostringstream oss;
00427   Teuchos::FancyOStream omsg(rcp(&oss,false));
00428 
00429   omsg << "Solver: " << this->description() << endl;
00430 
00431   if(converged) {
00432     solveStatus.solveStatus = Thyra::SOLVE_STATUS_CONVERGED;
00433     omsg << "CVODE status test converged!\n";
00434   }
00435   else {
00436     solveStatus.solveStatus = Thyra::SOLVE_STATUS_UNCONVERGED;
00437     omsg << "CVODE status test failed!\n";
00438   }
00439 
00440   if (sawFailedLinearSolve) {
00441     omsg << "Warning!  A failed linear solve was encountered with status:\n";
00442     OSTab tab(omsg);
00443     omsg << failedLinearSolveStatus;
00444   }
00445 
00446   omsg
00447     << "R*||dx|| = " << R << "*" << nrm_dx
00448     << " <= nonlinearSafetyFactor*tol = " << nonlinearSafetyFactor_ << "*" << tol << " : "
00449     << ( converged ? "converged!" : " unconverged" ) << endl;
00450 
00451   omsg
00452     << "Iterations = " << iter;
00453   // Above, we leave off the last newline since this is the convention for the
00454   // SolveStatus::message string!
00455 
00456   solveStatus.message = oss.str();
00457 
00458   // Update the solution state for external clients
00459   current_x_ = x->clone_v();
00460   J_is_current_ = false;
00461   // 2007/09/04: rabartl: Note, above the Jacobian J is always going to be out
00462   // of date since this algorithm computes x_curr = x_curr + dx for at least
00463   // one solve for dx = -inv(J)*f.  Therefore, J is never at the updated
00464   // x_curr, only the old x_curr!
00465 
00466   if (showNewtonDetails)
00467     *out << "\nLeaving TimeStepNonlinearSolver::solve(...) ...\n";
00468 
00469   return solveStatus;
00470 
00471 }
00472 
00473 
00474 template <class Scalar>
00475 bool TimeStepNonlinearSolver<Scalar>::supportsCloning() const
00476 {
00477   return true;
00478 }
00479 
00480 
00481 template <class Scalar>
00482 RCP<Thyra::NonlinearSolverBase<Scalar> >
00483 TimeStepNonlinearSolver<Scalar>::cloneNonlinearSolver() const
00484 {
00485   RCP<TimeStepNonlinearSolver<Scalar> >
00486     nonlinearSolver = Teuchos::rcp(new TimeStepNonlinearSolver<Scalar>);
00487   nonlinearSolver->model_ = model_; // Shallow copy is okay, model is stateless
00488   nonlinearSolver->defaultTol_ = defaultTol_;
00489   nonlinearSolver->defaultMaxIters_ = defaultMaxIters_;
00490   nonlinearSolver->nonlinearSafetyFactor_ = nonlinearSafetyFactor_;
00491   nonlinearSolver->linearSafetyFactor_ = linearSafetyFactor_;
00492   nonlinearSolver->RMinFraction_ = RMinFraction_;
00493   nonlinearSolver->throwOnLinearSolveFailure_ = throwOnLinearSolveFailure_;
00494   // Note: The specification of this virtual function in the interface class
00495   // allows us to just copy the algorithm, not the entire state so we are
00496   // done!
00497   return nonlinearSolver;
00498 }
00499 
00500 
00501 template <class Scalar>
00502 RCP<const Thyra::VectorBase<Scalar> >
00503 TimeStepNonlinearSolver<Scalar>::get_current_x() const
00504 {
00505   return current_x_;
00506 }
00507 
00508 
00509 template <class Scalar>
00510 bool TimeStepNonlinearSolver<Scalar>::is_W_current() const
00511 {
00512   return J_is_current_;
00513 }
00514 
00515 
00516 template <class Scalar>
00517 RCP<Thyra::LinearOpWithSolveBase<Scalar> >
00518 TimeStepNonlinearSolver<Scalar>::get_nonconst_W(const bool forceUpToDate)
00519 {
00520   if (is_null(J_))
00521     return Teuchos::null;
00522   if (forceUpToDate) {
00523 #ifdef RYTHMOS_DEBUG
00524     TEST_FOR_EXCEPT(is_null(current_x_));
00525 #endif
00526     Thyra::eval_f_W<Scalar>( *model_, *current_x_, 0, &*J_ );
00527     J_is_current_ = true;
00528   }
00529   return J_;
00530 }
00531 
00532 
00533 template <class Scalar>
00534 RCP<const Thyra::LinearOpWithSolveBase<Scalar> >
00535 TimeStepNonlinearSolver<Scalar>::get_W() const
00536 {
00537   return J_;
00538 }
00539 
00540 
00541 template <class Scalar>
00542 void TimeStepNonlinearSolver<Scalar>::set_W_is_current(bool W_is_current)
00543 {
00544   J_is_current_ = W_is_current;
00545 }
00546 
00547 // 
00548 // Explicit Instantiation macro
00549 //
00550 // Must be expanded from within the Rythmos namespace!
00551 //
00552 
00553 #define RYTHMOS_TIME_STEP_NONLINEAR_SOLVER_INSTANT(SCALAR) \
00554   \
00555   template class TimeStepNonlinearSolver< SCALAR >; \
00556   \
00557   template RCP<TimeStepNonlinearSolver< SCALAR > > timeStepNonlinearSolver(); 
00558    
00559 
00560 
00561 
00562 } // namespace Rythmos
00563 
00564 
00565 #endif // RYTHMOS_TIME_STEP_NONLINEAR_SOLVER_DEF_HPP

Generated on Wed May 12 21:25:44 2010 for Rythmos - Transient Integration for Differential Equations by  doxygen 1.4.7