Rythmos_BackwardEulerStepper.hpp

Go to the documentation of this file.
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 #ifndef Rythmos_BACKWARD_EULER_STEPPER_H
00030 #define Rythmos_BACKWARD_EULER_STEPPER_H
00031 
00032 #include "Rythmos_Stepper.hpp"
00033 #include "Teuchos_RefCountPtr.hpp"
00034 #include "Thyra_VectorBase.hpp"
00035 #include "Thyra_ModelEvaluator.hpp"
00036 #include "Thyra_ModelEvaluatorHelpers.hpp"
00037 #include "Thyra_NonlinearSolverBase.hpp"
00038 #include "Thyra_SingleResidSSDAEModelEvaluator.hpp"
00039 
00040 namespace Rythmos {
00041 
00043 template<class Scalar>
00044 class BackwardEulerStepper : virtual public Stepper<Scalar>
00045 {
00046   public:
00047 
00048     typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType ScalarMag;
00049 
00051     BackwardEulerStepper(
00052       const Teuchos::RefCountPtr<const Thyra::ModelEvaluator<Scalar> >  &model
00053       ,const Teuchos::RefCountPtr<Thyra::NonlinearSolverBase<Scalar> >  &solver
00054       );
00055 
00057     void setModel(const Teuchos::RefCountPtr<const Thyra::ModelEvaluator<Scalar> > &model);
00058 
00060     void setSolver(const Teuchos::RefCountPtr<Thyra::NonlinearSolverBase<Scalar> > &solver);
00061 
00063     Scalar TakeStep(Scalar dt);
00064    
00066     Scalar TakeStep();
00067 
00069     Teuchos::RefCountPtr<const Thyra::VectorBase<Scalar> > get_solution() const;
00070 
00072     Teuchos::RefCountPtr<const Thyra::VectorBase<Scalar> > get_residual() const;
00073 
00075 
00076     std::string description() const;
00077 
00079     void describe(
00080       Teuchos::FancyOStream                &out
00081       ,const Teuchos::EVerbosityLevel      verbLevel
00082       ) const;
00083 
00086     bool SetPoints(
00087       const std::vector<Scalar>& time_list
00088       ,const std::vector<Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > >& x_list
00089       ,const std::vector<Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > >& xdot_list);
00090     
00092     bool GetPoints(
00093       const std::vector<Scalar>& time_list
00094       ,std::vector<Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > >* x_list
00095       ,std::vector<Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > >* xdot_list
00096       ,std::vector<ScalarMag>* accuracy_list) const;
00097 
00099     bool SetRange(
00100       const Scalar& time_lower
00101       ,const Scalar& time_upper
00102       ,const InterpolationBuffer<Scalar> & IB);
00103 
00105     bool GetNodes(std::vector<Scalar>* time_list) const;
00106 
00108     bool RemoveNodes(std::vector<Scalar>& time_list) const;
00109 
00111     int GetOrder() const;
00112 
00113   private:
00114 
00115     Teuchos::RefCountPtr<const Thyra::ModelEvaluator<Scalar> > model_;
00116     Teuchos::RefCountPtr<Thyra::NonlinearSolverBase<Scalar> > solver_;
00117     Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > x_;
00118     Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > scaled_x_old_;
00119     Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > f_;
00120     Scalar t_;
00121     Scalar t_old_;
00122 
00123     Teuchos::RefCountPtr<Thyra::SingleResidSSDAEModelEvaluator<Scalar> >  neModel_;
00124 
00125 };
00126 
00127 // ////////////////////////////
00128 // Defintions
00129 
00130 template<class Scalar>
00131 BackwardEulerStepper<Scalar>::BackwardEulerStepper(
00132   const Teuchos::RefCountPtr<const Thyra::ModelEvaluator<Scalar> > &model
00133   ,const Teuchos::RefCountPtr<Thyra::NonlinearSolverBase<Scalar> > &solver
00134   )
00135 {
00136   setModel(model);
00137   setSolver(solver);
00138 }
00139 
00140 template<class Scalar>
00141 void BackwardEulerStepper<Scalar>::setModel(const Teuchos::RefCountPtr<const Thyra::ModelEvaluator<Scalar> > &model)
00142 {
00143   typedef Teuchos::ScalarTraits<Scalar> ST;
00144   model_ = model;
00145   t_ = ST::zero();
00146   x_ = model_->getNominalValues().get_x()->clone_v();
00147   f_ = Thyra::createMember(model_->get_f_space());
00148 
00149   scaled_x_old_ = x_->clone_v();
00150   
00151 }
00152 
00153 template<class Scalar>
00154 void BackwardEulerStepper<Scalar>::setSolver(const Teuchos::RefCountPtr<Thyra::NonlinearSolverBase<Scalar> > &solver)
00155 {
00156   solver_ = solver;
00157 }
00158 
00159 template<class Scalar>
00160 Scalar BackwardEulerStepper<Scalar>::TakeStep()
00161 {
00162   // print something out about this method not supporting automatic variable step-size
00163   typedef Teuchos::ScalarTraits<Scalar> ST;
00164   return(-ST::one());
00165 }
00166 
00167 template<class Scalar>
00168 Scalar BackwardEulerStepper<Scalar>::TakeStep(Scalar dt)
00169 {
00170   typedef Teuchos::ScalarTraits<Scalar> ST;
00171   //
00172   // Setup the nonlinear equations:
00173   //
00174   //   f( (1/dt)* x + (-1/dt)*x_old), x, t ) = 0
00175   //
00176   V_StV( &*scaled_x_old_, Scalar(-ST::one()/dt), *x_ );
00177   t_old_ = t_;
00178   if(!neModel_.get())
00179     neModel_ = Teuchos::rcp(new Thyra::SingleResidSSDAEModelEvaluator<Scalar>());
00180   neModel_->initialize(model_,Scalar(ST::one()/dt),scaled_x_old_,ST::one(),Teuchos::null,t_old_+dt,Teuchos::null);
00181   if(solver_->getModel().get()!=neModel_.get())
00182     solver_->setModel(neModel_);
00183   //
00184   // Solve the implicit nonlinear system to a tolerance of ???
00185   //
00186   Thyra::assign(&*x_,ST::zero());
00187   solver_->solve(&*x_); // Note that x in input is x_old and on output is the solved x!
00188   //
00189   // Update the step
00190   //
00191   t_ += dt;
00192 
00193   return(dt);
00194 }
00195 
00196 template<class Scalar>
00197 Teuchos::RefCountPtr<const Thyra::VectorBase<Scalar> > BackwardEulerStepper<Scalar>::get_solution() const
00198 {
00199   return(x_);
00200 }
00201 
00202 template<class Scalar>
00203 Teuchos::RefCountPtr<const Thyra::VectorBase<Scalar> > BackwardEulerStepper<Scalar>::get_residual() const
00204 {
00205   return(f_);
00206 }
00207 
00208 template<class Scalar>
00209 std::string BackwardEulerStepper<Scalar>::description() const
00210 {
00211   std::string name = "Rythmos::BackwardEulerStepper";
00212   return(name);
00213 }
00214 
00215 template<class Scalar>
00216 void BackwardEulerStepper<Scalar>::describe(
00217       Teuchos::FancyOStream                &out
00218       ,const Teuchos::EVerbosityLevel      verbLevel
00219       ) const
00220 {
00221   if (verbLevel == Teuchos::VERB_EXTREME)
00222   {
00223     out << description() << "::describe:";
00224     out << "\nmodel_ = " << std::endl;
00225     model_->describe(out,verbLevel);
00226     out << "\nsolver_ = " << std::endl;
00227     solver_->describe(out,verbLevel);
00228     out << "\nx_ = " << std::endl;
00229     x_->describe(out,verbLevel);
00230     out << "\nscaled_x_old_ = " << std::endl;
00231     scaled_x_old_->describe(out,verbLevel);
00232     out << "\nf_ = " << std::endl;
00233     f_->describe(out,verbLevel);
00234     out << "\nt_ = " << t_;
00235     out << "\nt_old_ = " << t_old_;
00236     out << std::endl;
00237 //    out << "neModel_ = " << std::endl;
00238 //    out << neModel_->describe(out,verbLevel) << std::endl;
00239   }
00240 }
00241 
00242 template<class Scalar>
00243 bool BackwardEulerStepper<Scalar>::SetPoints(
00244     const std::vector<Scalar>& time_list
00245     ,const std::vector<Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > >& x_list
00246     ,const std::vector<Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > >& xdot_list)
00247 {
00248   return(false);
00249 }
00250 
00251 template<class Scalar>
00252 bool BackwardEulerStepper<Scalar>::GetPoints(
00253     const std::vector<Scalar>& time_list
00254     ,std::vector<Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > >* x_list
00255     ,std::vector<Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > >* xdot_list
00256     ,std::vector<ScalarMag>* accuracy_list) const
00257 {
00258   return(false);
00259 }
00260 
00261 template<class Scalar>
00262 bool BackwardEulerStepper<Scalar>::SetRange(
00263     const Scalar& time_lower
00264     ,const Scalar& time_upper
00265     ,const InterpolationBuffer<Scalar>& IB)
00266 {
00267   return(false);
00268 }
00269 
00270 template<class Scalar>
00271 bool BackwardEulerStepper<Scalar>::GetNodes(std::vector<Scalar>* time_list) const
00272 {
00273   return(false);
00274 }
00275 
00276 template<class Scalar>
00277 bool BackwardEulerStepper<Scalar>::RemoveNodes(std::vector<Scalar>& time_list) const
00278 {
00279   return(false);
00280 }
00281 
00282 template<class Scalar>
00283 int BackwardEulerStepper<Scalar>::GetOrder() const
00284 {
00285   return(1);
00286 }
00287 
00288 
00289 } // namespace Rythmos
00290 
00291 #endif //Rythmos_BACKWARD_EULER_STEPPER_H

Generated on Thu Sep 18 12:30:05 2008 for Rythmos - Transient Integration for Differential Equations by doxygen 1.3.9.1