Rythmos_ForwardEulerStepper.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 #ifndef Rythmos_FORWARDEULER_STEPPER_H
00030 #define Rythmos_FORWARDEULER_STEPPER_H
00031 
00032 #include "Rythmos_StepperBase.hpp"
00033 #include "Teuchos_RCP.hpp"
00034 #include "Thyra_VectorBase.hpp"
00035 #include "Thyra_ModelEvaluator.hpp"
00036 #include "Thyra_ModelEvaluatorHelpers.hpp"
00037 
00038 namespace Rythmos {
00039 
00041 template<class Scalar>
00042 class ForwardEulerStepper : virtual public StepperBase<Scalar>
00043 {
00044   public:
00045 
00046     typedef Teuchos::ScalarTraits<Scalar> ST;
00047     typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType ScalarMag;
00048     
00050     ForwardEulerStepper();
00051     ForwardEulerStepper(const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > &model);
00052 
00054     void setModel(const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > &model);
00055 
00057     Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >
00058     getModel() const;
00059 
00061     RCP<const Thyra::VectorSpaceBase<Scalar> > get_x_space() const;
00062     
00064     ~ForwardEulerStepper();
00065 
00067     Scalar takeStep(Scalar dt, StepSizeType flag);
00068 
00070     const StepStatus<Scalar> getStepStatus() const;
00071 
00073     std::string description() const;
00074 
00076     std::ostream& describe(
00077       std::ostream                &out
00078       ,const Teuchos::EVerbosityLevel      verbLevel
00079       ,const std::string          leadingIndent
00080       ,const std::string          indentSpacer
00081       ) const;
00082     
00085     void addPoints(
00086       const Array<Scalar>& time_vec
00087       ,const Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >& x_vec
00088       ,const Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >& xdot_vec
00089       );
00090     
00092     void getPoints(
00093       const Array<Scalar>& time_vec
00094       ,Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >* x_vec
00095       ,Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec
00096       ,Array<ScalarMag>* accuracy_vec
00097       ) const;
00098 
00100     void setRange(
00101       const TimeRange<Scalar>& range,
00102       const InterpolationBufferBase<Scalar> & IB
00103       );
00104 
00106     TimeRange<Scalar> getTimeRange() const;
00107 
00109     void getNodes(Array<Scalar>* time_vec) const;
00110 
00112     void removeNodes(Array<Scalar>& time_vec);
00113 
00115     int getOrder() const;
00116 
00118 
00119     void setParameterList(Teuchos::RCP<Teuchos::ParameterList> const& paramList);
00120 
00122     Teuchos::RCP<Teuchos::ParameterList> getParameterList();
00123 
00125     Teuchos::RCP<Teuchos::ParameterList> unsetParameterList();
00126 
00127   private:
00128 
00129     Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > model_;
00130     Teuchos::RCP<Thyra::VectorBase<Scalar> > solution_vector_;
00131     Teuchos::RCP<Thyra::VectorBase<Scalar> > residual_vector_;
00132     Scalar t_;
00133     Scalar dt_;
00134 
00135     Teuchos::RCP<Teuchos::ParameterList> parameterList_;
00136     bool isInitialized_;
00137 
00138     // Private member functions:
00139     void initialize_();
00140 
00141 };
00142 
00143 template<class Scalar>
00144 ForwardEulerStepper<Scalar>::ForwardEulerStepper(const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > &model)
00145   : isInitialized_(false)
00146 {
00147   this->setModel(model);
00148   initialize_();
00149 }
00150 
00151 template<class Scalar>
00152 void ForwardEulerStepper<Scalar>::initialize_()
00153 {
00154   t_ = ST::zero();
00155   solution_vector_ = model_->getNominalValues().get_x()->clone_v();
00156   residual_vector_ = Thyra::createMember(model_->get_f_space());
00157   isInitialized_ = true;
00158 }
00159 
00160 template<class Scalar>
00161 ForwardEulerStepper<Scalar>::ForwardEulerStepper()
00162   : isInitialized_(false)
00163 {
00164 }
00165 
00166 template<class Scalar>
00167 ForwardEulerStepper<Scalar>::~ForwardEulerStepper()
00168 {
00169 }
00170 
00171 template<class Scalar>
00172 RCP<const Thyra::VectorSpaceBase<Scalar> > ForwardEulerStepper<Scalar>::get_x_space() const
00173 {
00174   TEST_FOR_EXCEPTION(!isInitialized_,std::logic_error,"Error, attempting to call get_x_space before initialization!\n");
00175   return(solution_vector_->space());
00176 }
00177 
00178 template<class Scalar>
00179 Scalar ForwardEulerStepper<Scalar>::takeStep(Scalar dt, StepSizeType flag)
00180 {
00181   if (flag == STEP_TYPE_VARIABLE) { 
00182     // print something out about this method not supporting automatic variable step-size
00183     typedef Teuchos::ScalarTraits<Scalar> ST;
00184     return(-ST::one());
00185   }
00186 /*
00187   Thyra::ModelEvaluatorBase::InArgs<Scalar>   inArgs  = model_->createInArgs();
00188   Thyra::ModelEvaluatorBase::OutArgs<Scalar>  outArgs = model_->createOutArgs();
00189 
00190   inArgs.set_x(solution_vector_);
00191   inArgs.set_t(t_+dt);
00192 
00193   outArgs.set_f(residual_vector_);
00194 
00195   model_->evalModel(inArgs,outArgs);
00196 */
00197   Thyra::eval_f<Scalar>(*model_,*solution_vector_,t_+dt,&*residual_vector_);
00198 
00199   // solution_vector = solution_vector + dt*residual_vector
00200   Thyra::Vp_StV(&*solution_vector_,dt,*residual_vector_); 
00201   t_ += dt;
00202 
00203   dt_ = dt;
00204 
00205   return(dt);
00206 }
00207 
00208 template<class Scalar>
00209 const StepStatus<Scalar> ForwardEulerStepper<Scalar>::getStepStatus() const
00210 {
00211   typedef Teuchos::ScalarTraits<Scalar> ST;
00212   StepStatus<Scalar> stepStatus;
00213 
00214   stepStatus.stepSize = dt_; 
00215   stepStatus.order = 1;
00216   stepStatus.time = t_;
00217   stepStatus.stepLETValue = Scalar(-ST::one()); 
00218   stepStatus.solution = solution_vector_;
00219   stepStatus.residual = residual_vector_;
00220 
00221   return(stepStatus);
00222 }
00223 
00224 template<class Scalar>
00225 std::string ForwardEulerStepper<Scalar>::description() const
00226 {
00227   std::string name = "Rythmos::ForwardEulerStepper";
00228   return(name);
00229 }
00230 
00231 template<class Scalar>
00232 std::ostream& ForwardEulerStepper<Scalar>::describe(
00233       std::ostream                &out
00234       ,const Teuchos::EVerbosityLevel      verbLevel
00235       ,const std::string          leadingIndent
00236       ,const std::string          indentSpacer
00237       ) const
00238 {
00239   if ( (static_cast<int>(verbLevel) == static_cast<int>(Teuchos::VERB_DEFAULT) ) ||
00240        (static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW)     )
00241      ) {
00242     out << description() << "::describe" << std::endl;
00243     out << "model = " << model_->description() << std::endl;
00244   } else if (static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW)) {
00245   } else if (static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_MEDIUM)) {
00246   } else if (static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_HIGH)) {
00247     out << "model = " << std::endl;
00248     model_->describe(out,verbLevel,leadingIndent,indentSpacer);
00249     out << "solution_vector = " << std::endl;
00250     solution_vector_->describe(out,verbLevel,leadingIndent,indentSpacer);
00251     out << "residual_vector = " << std::endl;
00252     residual_vector_->describe(out,verbLevel,leadingIndent,indentSpacer);
00253   }
00254   return(out);
00255 }
00256 
00257 template<class Scalar>
00258 void ForwardEulerStepper<Scalar>::addPoints(
00259     const Array<Scalar>& time_vec
00260     ,const Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >& x_vec
00261     ,const Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >& xdot_vec
00262     )
00263 {
00264   TEST_FOR_EXCEPTION(true,std::logic_error,"Error, addPoints is not implemented for ForwardEulerStepper.\n");
00265 }
00266 
00267 template<class Scalar>
00268 void ForwardEulerStepper<Scalar>::getPoints(
00269     const Array<Scalar>& time_vec
00270     ,Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >* x_vec
00271     ,Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec
00272     ,Array<ScalarMag>* accuracy_vec) const
00273 {
00274   TEST_FOR_EXCEPTION(true,std::logic_error,"Error, getPoints is not implemented for ForwardEulerStepper.\n");
00275 }
00276 
00277 template<class Scalar>
00278 TimeRange<Scalar> ForwardEulerStepper<Scalar>::getTimeRange() const
00279 {
00280   return invalidTimeRange<Scalar>();
00281 }
00282 
00283 template<class Scalar>
00284 void ForwardEulerStepper<Scalar>::getNodes(Array<Scalar>* time_vec) const
00285 {
00286   TEST_FOR_EXCEPTION(true,std::logic_error,"Error, getNodes is not implemented for ForwardEulerStepper.\n");
00287 }
00288 
00289 template<class Scalar>
00290 void ForwardEulerStepper<Scalar>::removeNodes(Array<Scalar>& time_vec) 
00291 {
00292   TEST_FOR_EXCEPTION(true,std::logic_error,"Error, removeNodes is not implemented for ForwardEulerStepper.\n");
00293 }
00294 
00295 template<class Scalar>
00296 int ForwardEulerStepper<Scalar>::getOrder() const
00297 {
00298   return(1);
00299 }
00300 
00301 template <class Scalar>
00302 void ForwardEulerStepper<Scalar>::setParameterList(Teuchos::RCP<Teuchos::ParameterList> const& paramList)
00303 {
00304   parameterList_ = paramList;
00305   int outputLevel = parameterList_->get( "outputLevel", int(-1) );
00306   outputLevel = std::min(std::max(outputLevel,-1),4);
00307   this->setVerbLevel(static_cast<Teuchos::EVerbosityLevel>(outputLevel));
00308 }
00309 
00310 template <class Scalar>
00311 Teuchos::RCP<Teuchos::ParameterList> ForwardEulerStepper<Scalar>::getParameterList()
00312 {
00313   return(parameterList_);
00314 }
00315 
00316 template <class Scalar>
00317 Teuchos::RCP<Teuchos::ParameterList> ForwardEulerStepper<Scalar>::unsetParameterList()
00318 {
00319   Teuchos::RCP<Teuchos::ParameterList> temp_param_list = parameterList_;
00320   parameterList_ = Teuchos::null;
00321   return(temp_param_list);
00322 }
00323 
00324 template<class Scalar>
00325 void ForwardEulerStepper<Scalar>::setModel(const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > &model)
00326 {
00327   TEST_FOR_EXCEPT(model == Teuchos::null)
00328   model_ = model;
00329 }
00330 
00331 template<class Scalar>
00332 Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >
00333 ForwardEulerStepper<Scalar>::getModel() const
00334 {
00335   return model_;
00336 }
00337 
00338 } // namespace Rythmos
00339 
00340 #endif //Rythmos_FORWARDEULER_STEPPER_H

Generated on Tue Oct 20 12:46:07 2009 for Rythmos - Transient Integration for Differential Equations by doxygen 1.4.7