Thyra_TimeStepNewtonNonlinearSolver.hpp

00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 //    Thyra: Interfaces and Support for Abstract Numerical Algorithms
00005 //                 Copyright (2004) 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 Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #ifndef THYRA_NEWTON_NONLINEAR_SOLVER_HPP
00030 #define THYRA_NEWTON_NONLINEAR_SOLVER_HPP
00031 
00032 #include "Thyra_NonlinearSolverBase.hpp"
00033 #include "Teuchos_StandardMemberCompositionMacros.hpp"
00034 #include "Teuchos_StandardCompositionMacros.hpp"
00035 
00036 namespace Thyra {
00037 
00042 template <class Scalar>
00043 class TimeStepNewtonNonlinearSolver : public NonlinearSolverBase<Scalar> {
00044 public:
00045 
00047   typedef Teuchos::ScalarTraits<Scalar> ST;
00049   typedef typename ST::magnitudeType ScalarMag;
00051   typedef Teuchos::ScalarTraits<ScalarMag> SMT;
00052 
00054   STANDARD_MEMBER_COMPOSITION_MEMBERS( int, defaultMaxIterations )
00055 
00056   
00057    STANDARD_MEMBER_COMPOSITION_MEMBERS( ScalarMag, defaultTol )
00058 
00060    STANDARD_NONCONST_COMPOSITION_MEMBERS( std::ostream, warningOut )
00061 
00062   TimeStepNewtonNonlinearSolver(
00063     const int                 defaultMaxIterations = 3
00064     ,const ScalarMag          defaultTol           = 1e-2
00065     ,const warningOut_ptr_t   &warningOut          = Teuchos::rcp(&std::cerr,false)
00066     );
00067 
00070 
00071 
00073   void setModel(
00074     const Teuchos::RefCountPtr<const ModelEvaluator<Scalar> > &model
00075     );
00077   Teuchos::RefCountPtr<const ModelEvaluator<Scalar> > getModel() const;
00079   SolveStatus<Scalar> solve(
00080     VectorBase<Scalar>              *x
00081     ,const SolveCriteria<Scalar>    *solveCriteria
00082     ,VectorBase<Scalar>             *delta = NULL
00083     );
00085   Teuchos::RefCountPtr<const VectorBase<Scalar> > get_current_x() const;
00087   bool is_W_current() const;
00089   Teuchos::RefCountPtr<LinearOpWithSolveBase<Scalar> > get_nonconst_W();
00091   Teuchos::RefCountPtr<const LinearOpWithSolveBase<Scalar> > get_W() const;
00093   void set_W_is_current(bool W_is_current);
00094 
00096 
00097 private:
00098 
00099   Teuchos::RefCountPtr<const ModelEvaluator<Scalar> >    model_;
00100   Teuchos::RefCountPtr<LinearOpWithSolveBase<Scalar> >   J_;
00101   Teuchos::RefCountPtr<VectorBase<Scalar> >              current_x_;
00102   bool                                                   J_is_current_;
00103 
00104 };
00105 
00106 // ////////////////////////
00107 // Defintions
00108 
00109 template <class Scalar>
00110 TimeStepNewtonNonlinearSolver<Scalar>::TimeStepNewtonNonlinearSolver(
00111   const int                 defaultMaxIterations
00112   ,const ScalarMag          defaultTol
00113   ,const warningOut_ptr_t   &warningOut
00114   )
00115   :defaultMaxIterations_(defaultMaxIterations)
00116   ,defaultTol_(defaultTol)
00117   ,warningOut_(warningOut)
00118   ,J_is_current_(false)
00119 {}
00120 
00121 // Overridden from NonlinearSolverBase
00122 
00123 template <class Scalar>
00124 void TimeStepNewtonNonlinearSolver<Scalar>::setModel(
00125   const Teuchos::RefCountPtr<const ModelEvaluator<Scalar> > &model
00126   )
00127 {
00128   TEST_FOR_EXCEPT(model.get()==NULL);
00129   model_ = model;
00130   J_ = Teuchos::null;
00131   current_x_ = Teuchos::null;
00132   J_is_current_ = false;
00133 }
00134 
00135 template <class Scalar>
00136 Teuchos::RefCountPtr<const ModelEvaluator<Scalar> >
00137 TimeStepNewtonNonlinearSolver<Scalar>::getModel() const
00138 {
00139   return model_;
00140 }
00141 
00142 template <class Scalar>
00143 SolveStatus<Scalar> TimeStepNewtonNonlinearSolver<Scalar>::solve(
00144   VectorBase<Scalar>             *x
00145   ,const SolveCriteria<Scalar>   *solveCriteria
00146   ,VectorBase<Scalar>            *delta
00147   )
00148 {
00149   using std::endl;
00150   using Teuchos::rcp;
00151   using Teuchos::OSTab;
00152   using Teuchos::getFancyOStream;
00153   TEST_FOR_EXCEPT(solveCriteria!=NULL); // ToDo: Pass to linear solver?
00154   // Initialize storage for algorithm
00155   if(!J_.get())                                 J_ = model_->create_W();
00156   Teuchos::RefCountPtr<VectorBase<Scalar> >     f = createMember(model_->get_f_space());
00157   Teuchos::RefCountPtr<VectorBase<Scalar> >     dx = createMember(model_->get_x_space());
00158   Teuchos::RefCountPtr<VectorBase<Scalar> >     dx_last = createMember(model_->get_x_space());
00159   Teuchos::RefCountPtr<VectorBase<Scalar> >     x_curr = createMember(model_->get_x_space());
00160   if (delta != NULL)
00161     Thyra::V_S(delta,ST::zero()); // delta stores the cumulative update to x over the whole Newton solve.
00162   Thyra::assign(&*x_curr,*x);
00163   // Initialize convergence criteria
00164   ScalarMag R = SMT::one();
00165   ScalarMag RMin = 0.3; // ToDo: Make this adjustable
00166   ScalarMag tolSafety = 0.1; // ToDo: Make this adjustable
00167   ScalarMag linearTolSafety = 0.05 * tolSafety; // ToDo: Make this adjustable
00168   int maxIters = this->defaultMaxIterations();
00169   ScalarMag tol = this->defaultTol();
00170   // ToDo: Get these from solveCriteria
00171   bool converged = false;
00172   ScalarMag nrm_dx_last;
00173   int iter = 1;
00174   for( ; iter <= maxIters; ++iter ) {
00175     // Evaluate f and W
00176     eval_f_W( *model_, *x_curr, ST::one(), &*f, &*J_ );
00177     // Zero out dx to deal with iterative linear solvers
00178     Thyra::V_S(&*dx,ST::zero());
00179     // Solve the system: J*dx = -f
00180     SolveCriteria<Scalar>
00181       linearSolveCriteria(
00182         SolveMeasureType(
00183           SOLVE_MEASURE_NORM_RESIDUAL,SOLVE_MEASURE_NORM_RHS
00184           )
00185         ,linearTolSafety*tol
00186         );
00187     SolveStatus<Scalar>
00188       linearSolveStatus = Thyra::solve( *J_, NOTRANS, *f, &*dx, &linearSolveCriteria );
00189     Thyra::Vt_S(&*dx,Scalar(-ST::one()));
00190     if (delta != NULL)
00191       Thyra::Vp_V(delta,*dx); 
00192     // Check the linear solve
00193     if(linearSolveStatus.solveStatus != SOLVE_STATUS_CONVERGED) {
00194       warningOut()
00195         << "Thyra::TimeStepNewtonNonlinearSolver<Scalar>::solve(...), Warning, linear solve did not converge with solve status:\n\n";
00196       *OSTab(getFancyOStream(rcp(&warningOut(),false))).getOStream() << linearSolveStatus;
00197       warningOut()
00198         <<endl<< "Continuing anyway :-)\n";
00199       // ToDo: Add option to throw exception failure
00200     }
00201     // Update the solution: x_curr = x_curr + dx
00202     Vp_V( &*x_curr, *dx );
00203     // Convergence test
00204     const ScalarMag nrm_dx = Thyra::norm(*dx);
00205     if(R*nrm_dx < tolSafety*tol) {
00206       converged = true;
00207       break;
00208     }
00209     // Update convergence criteria
00210     if(iter > 1) {
00211       R = std::max(RMin*R,nrm_dx/nrm_dx_last);
00212     }
00213     // Save to old
00214     std::swap(dx_last,dx);
00215     nrm_dx_last = nrm_dx;
00216   }
00217   // Set the solution
00218   Thyra::assign(x,*x_curr);
00219   // Check the status
00220   SolveStatus<Scalar> solveStatus;
00221   std::ostringstream oss;
00222   if(converged) {
00223     solveStatus.solveStatus = SOLVE_STATUS_CONVERGED;
00224     oss << "CVODE status test converged!  Iterations = " << iter << ".";
00225   }
00226   else {
00227     solveStatus.solveStatus = SOLVE_STATUS_UNCONVERGED;
00228     oss << "CVODE status test failed!  Iterations = " << iter << ".";
00229   }
00230   solveStatus.message = oss.str();
00231   //
00232   current_x_ = x->clone_v();
00233   J_is_current_ = true;
00234   //
00235   return solveStatus;
00236 }
00237 
00238 template <class Scalar>
00239 Teuchos::RefCountPtr<const VectorBase<Scalar> >
00240 TimeStepNewtonNonlinearSolver<Scalar>::get_current_x() const
00241 {
00242   return current_x_;
00243 }
00244 
00245 template <class Scalar>
00246 bool TimeStepNewtonNonlinearSolver<Scalar>::is_W_current() const
00247 {
00248   return J_is_current_;
00249 }
00250 
00251 template <class Scalar>
00252 Teuchos::RefCountPtr<LinearOpWithSolveBase<Scalar> >
00253 TimeStepNewtonNonlinearSolver<Scalar>::get_nonconst_W()
00254 {
00255   return J_;
00256 }
00257 
00258 template <class Scalar>
00259 Teuchos::RefCountPtr<const LinearOpWithSolveBase<Scalar> >
00260 TimeStepNewtonNonlinearSolver<Scalar>::get_W() const
00261 {
00262   return J_;
00263 }
00264 
00265 template <class Scalar>
00266 void TimeStepNewtonNonlinearSolver<Scalar>::set_W_is_current(bool W_is_current)
00267 {
00268   J_is_current_ = W_is_current;
00269 }
00270 
00271 } // namespace Thyra
00272 
00273 #endif // THYRA_NEWTON_NONLINEAR_SOLVER_HPP

Generated on Thu Sep 18 12:32:51 2008 for Thyra Nonlinear Solver Support by doxygen 1.3.9.1