Rythmos_ExplicitRKStepper.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_ExplicitRK_STEPPER_H
00030 #define Rythmos_ExplicitRK_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 ExplicitRKStepper : virtual public StepperBase<Scalar>
00043 {
00044   public:
00045     typedef Teuchos::ScalarTraits<Scalar> ST;
00046     typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType ScalarMag;
00047     
00049     ExplicitRKStepper();
00050     ExplicitRKStepper(const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > &model_);
00051 
00053     Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > get_x_space() const;
00054 
00056     void setModel(const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > &model);
00057 
00059     Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >
00060     getModel() const;
00061     
00063     ~ExplicitRKStepper();
00064 
00066     Scalar takeStep(Scalar dt, StepSizeType flag);
00067 
00069     const StepStatus<Scalar> getStepStatus() const;
00070 
00072     std::ostream& describe(
00073       std::ostream                &out
00074       ,const Teuchos::EVerbosityLevel      verbLevel
00075       ,const std::string          leadingIndent
00076       ,const std::string          indentSpacer
00077       ) const;
00078 
00081     void addPoints(
00082       const Array<Scalar>& time_vec
00083       ,const Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >& x_vec
00084       ,const Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >& xdot_vec
00085       );
00086 
00088     void getPoints(
00089       const Array<Scalar>& time_vec
00090       ,Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >* x_vec
00091       ,Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec
00092       ,Array<ScalarMag>* accuracy_vec) const;
00093 
00095     TimeRange<Scalar> getTimeRange() const;
00096 
00098     void getNodes(Array<Scalar>* time_vec) const;
00099 
00101     void removeNodes(Array<Scalar>& time_vec);
00102 
00104     int getOrder() const;
00105 
00107 
00108     void setParameterList(Teuchos::RCP<Teuchos::ParameterList> const& paramList);
00109 
00111     Teuchos::RCP<Teuchos::ParameterList> getNonconstParameterList();
00112 
00114     Teuchos::RCP<Teuchos::ParameterList> unsetParameterList();
00115     
00116 
00117   private:
00118 
00119     Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > model_;
00120     Teuchos::RCP<Thyra::VectorBase<Scalar> > solution_vector_;
00121     Array<Teuchos::RCP<Thyra::VectorBase<Scalar> > > k_vector_;
00122     Teuchos::RCP<Thyra::VectorBase<Scalar> > ktemp_vector_;
00123 
00124     int stages_; // Number of stages of RK
00125     Array<Array<Scalar> > b_A_; // Butcher tableau A matrix
00126     Array<Scalar> b_b_; // Butcher b vector
00127     Array<Scalar> b_c_; // Butcher c vector
00128 
00129     Scalar t_;
00130     Scalar dt_;
00131 
00132     Teuchos::RCP<Teuchos::ParameterList> parameterList_;
00133 
00134     bool isInitialized_;
00135 
00136     // Private member functions:
00137     void initialize_();
00138 
00139 };
00140 
00141 template<class Scalar>
00142 ExplicitRKStepper<Scalar>::ExplicitRKStepper(const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > &model)
00143 {
00144   Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
00145   out->precision(15);
00146 
00147   this->setModel(model);
00148   initialize_();
00149 }
00150 
00151 template<class Scalar>
00152 void ExplicitRKStepper<Scalar>::initialize_()
00153 {
00154   t_ = ST::zero();
00155   solution_vector_ = model_->getNominalValues().get_x()->clone_v();
00156   stages_ = 4; // 4 stage ERK
00157   k_vector_.reserve(stages_);
00158   Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
00159     f_space = model_->get_f_space();
00160   for (int i=0 ; i<stages_ ; ++i) {
00161     k_vector_.push_back(Thyra::createMember(f_space));
00162   }
00163   ktemp_vector_ = Thyra::createMember(f_space);
00164 
00165   // initialize the Butcher tableau and its vectors
00166   Scalar zero = ST::zero();
00167   b_A_.reserve(stages_);
00168   b_b_.reserve(stages_);
00169   b_c_.reserve(stages_);
00170 
00171   // fill b with zeros
00172   b_b_.push_back(zero); 
00173   b_b_.push_back(zero); 
00174   b_b_.push_back(zero);
00175   b_b_.push_back(zero);
00176   // fill c with zeros
00177   b_c_.push_back(zero); 
00178   b_c_.push_back(zero); 
00179   b_c_.push_back(zero);
00180   b_c_.push_back(zero);
00181   // fill A with zeros
00182   for (int i=0 ; i<stages_ ; ++i) {
00183     b_A_.push_back(b_b_);
00184   }
00185 
00186   // Runge-Kutta methods in Butcher tableau form:
00187   //  c | A    A = sxs matrix
00188   //   -|---   b = s vector
00189   //      b'   c = s vector
00190   
00191 /*
00192   // 3/8 Rule Runge-Kutta Method: (not implemented yet)
00193   // c = [  0  1/3 2/3  1  ]'
00194   // A = [  0              ]
00195   //     [ 1/3  0          ]
00196   //     [-1/3  1   0      ]
00197   //     [  1  -1   1   0  ]
00198   // b = [ 1/8 3/8 3/8 1/8 ]'
00199   Scalar one = ST::one();
00200   Scalar one_third    = Scalar(ST::one()/(3*ST::one()));
00201   Scalar two_third    = Scalar(2*ST::one()/(3*ST::one()));
00202   Scalar one_eighth   = Scalar(ST::one()/(8*ST::one()));
00203   Scalar three_eighth = Scalar(3*ST::one()/(8*ST::one()));
00204 
00205   // fill b_A_
00206   b_A_[0][0] = zero;
00207   b_A_[0][1] = zero;
00208   b_A_[0][2] = zero;
00209   b_A_[0][3] = zero;
00210 
00211   b_A_[1][0] = one_third;
00212   b_A_[1][1] = zero;
00213   b_A_[1][2] = zero;
00214   b_A_[1][3] = zero;
00215 
00216   b_A_[2][0] = Scalar(-one_third);
00217   b_A_[2][1] = one;
00218   b_A_[2][2] = zero;
00219   b_A_[2][3] = zero;
00220 
00221   b_A_[3][0] = one;
00222   b_A_[3][1] = Scalar(-one);
00223   b_A_[3][2] = one;
00224   b_A_[3][3] = zero;
00225 
00226   // fill b_b_
00227   b_b_[0] = one_eighth;
00228   b_b_[1] = three_eighth;
00229   b_b_[2] = three_eighth;
00230   b_b_[3] = one_eighth;
00231   
00232   // fill b_c_
00233   b_c_[0] = zero;
00234   b_c_[1] = one_third;
00235   b_c_[2] = two_third;
00236   b_c_[3] = one;
00237 */
00238 
00239 
00240   // "The" Runge-Kutta Method: (implemented below)
00241   // c = [  0  1/2 1/2  1  ]'
00242   // A = [  0              ] 
00243   //     [ 1/2  0          ]
00244   //     [  0  1/2  0      ]
00245   //     [  0   0   1   0  ]
00246   // b = [ 1/6 1/3 1/3 1/6 ]'
00247   Scalar one = ST::one();
00248   Scalar onehalf = ST::one()/(2*ST::one());
00249   Scalar onesixth = ST::one()/(6*ST::one());
00250   Scalar onethird = ST::one()/(3*ST::one());
00251 
00252   // fill b_A_
00253   b_A_[0][0] = zero;
00254   b_A_[0][1] = zero;
00255   b_A_[0][2] = zero;
00256   b_A_[0][3] = zero;
00257 
00258   b_A_[1][0] = onehalf;
00259   b_A_[1][1] = zero;
00260   b_A_[1][2] = zero;
00261   b_A_[1][3] = zero;
00262 
00263   b_A_[2][0] = zero;
00264   b_A_[2][1] = onehalf;
00265   b_A_[2][2] = zero;
00266   b_A_[2][3] = zero;
00267 
00268   b_A_[3][0] = zero;
00269   b_A_[3][1] = zero;
00270   b_A_[3][2] = one;
00271   b_A_[3][3] = zero;
00272 
00273   // fill b_b_
00274   b_b_[0] = onesixth;
00275   b_b_[1] = onethird;
00276   b_b_[2] = onethird;
00277   b_b_[3] = onesixth;
00278   
00279   // fill b_c_
00280   b_c_[0] = zero;
00281   b_c_[1] = onehalf;
00282   b_c_[2] = onehalf;
00283   b_c_[3] = one;
00284 
00285   isInitialized_ = true;
00286 }
00287 
00288 template<class Scalar>
00289 ExplicitRKStepper<Scalar>::ExplicitRKStepper()
00290   : isInitialized_(false)
00291 {
00292 }
00293 
00294 template<class Scalar>
00295 ExplicitRKStepper<Scalar>::~ExplicitRKStepper()
00296 {
00297 }
00298 
00299 template<class Scalar>
00300 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > ExplicitRKStepper<Scalar>::get_x_space() const
00301 {
00302   TEST_FOR_EXCEPTION(!isInitialized_,std::logic_error,"Error, attempting to call get_x_space before initialization!\n");
00303   return(solution_vector_->space());
00304 }
00305 
00306 template<class Scalar>
00307 Scalar ExplicitRKStepper<Scalar>::takeStep(Scalar dt, StepSizeType flag)
00308 {
00309   typedef typename Thyra::ModelEvaluatorBase::InArgs<Scalar>::ScalarMag ScalarMag;
00310   if ((flag == STEP_TYPE_VARIABLE) || (dt == ST::zero())) {
00311     return(Scalar(-ST::one()));
00312   }
00313   dt_ = dt;
00314 
00315   // Compute stage solutions
00316   for (int s=0 ; s < stages_ ; ++s) {
00317     Thyra::assign(&*ktemp_vector_, *solution_vector_); // ktemp = solution_vector
00318     for (int j=0 ; j < s ; ++j) { // assuming Butcher matix is strictly lower triangular
00319       if (b_A_[s][j] != ST::zero())
00320         Thyra::Vp_StV(&*ktemp_vector_, b_A_[s][j], *k_vector_[j]); // ktemp = ktemp + a_{s+1,j+1}*k_{j+1}
00321     }
00322     ScalarMag ts = t_ + b_c_[s]*dt;
00323     Thyra::eval_f<Scalar>(*model_,*ktemp_vector_,ts,&*k_vector_[s]);
00324     Thyra::Vt_S(&*k_vector_[s],dt); // k_s = k_s*dt
00325   } 
00326   // Sum for solution:
00327   for (int s=0 ; s < stages_ ; ++s) {
00328     if (b_b_[s] != ST::zero()) {
00329       Thyra::Vp_StV(&*solution_vector_, b_b_[s], *k_vector_[s]); // solution_vector += b_{s+1}*k_{s+1}
00330     }
00331   }
00332 
00333   // update current time:
00334   t_ = t_ + dt;
00335 
00336   return(dt);
00337 }
00338 
00339 template<class Scalar>
00340 const StepStatus<Scalar> ExplicitRKStepper<Scalar>::getStepStatus() const
00341 {
00342   StepStatus<Scalar> stepStatus;
00343 
00344   stepStatus.stepSize = dt_;
00345   stepStatus.order = -1;
00346   stepStatus.time = t_;
00347   stepStatus.solution = solution_vector_;
00348 
00349   return(stepStatus);
00350 }
00351 
00352 template<class Scalar>
00353 std::ostream& ExplicitRKStepper<Scalar>::describe(
00354       std::ostream                &out
00355       ,const Teuchos::EVerbosityLevel      verbLevel
00356       ,const std::string          leadingIndent
00357       ,const std::string          indentSpacer
00358       ) const
00359 {
00360   if ( (static_cast<int>(verbLevel) == static_cast<int>(Teuchos::VERB_DEFAULT) ) ||
00361        (static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW)     )
00362      ) {
00363     out << this->description() << "::describe" << std::endl;
00364     out << "model = " << model_->description() << std::endl;
00365     out << stages_ << " stage Explicit RK method" << std::endl;
00366   } else if (static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW)) {
00367     out << "solution_vector = " << std::endl;
00368     solution_vector_->describe(out,verbLevel,leadingIndent,indentSpacer); 
00369   } else if (static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_MEDIUM)) {
00370   } else if (static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_HIGH)) {
00371     out << "model = " << std::endl;
00372     model_->describe(out,verbLevel,leadingIndent,indentSpacer); 
00373     for (int i=0 ; i<stages_ ; ++i) {
00374       out << "k_vector[" << i << "] = " << std::endl;
00375       out << k_vector_[i]->describe(out,verbLevel,leadingIndent,indentSpacer) << std::endl;
00376     }
00377     out << "ktemp_vector = " << std::endl;
00378     ktemp_vector_->describe(out,verbLevel,leadingIndent,indentSpacer); 
00379     for (int i=0 ; i<stages_ ; ++i) {
00380       for (int j=0 ; j<stages_ ; ++j) {
00381         out << "b_A_[" << i << "][" << j << "] = " << b_A_[i][j] << std::endl;
00382       }
00383     }
00384     for (int i=0 ; i<stages_ ; ++i) {
00385       out << "b_b_[" << i << "] = " << b_b_[i] << std::endl;
00386     }
00387     for (int i=0 ; i<stages_ ; ++i) {
00388       out << "b_c_[" << i << "] = " << b_c_[i] << std::endl;
00389     }
00390     out << "t = " << t_ << std::endl;
00391   }
00392   return(out);
00393 }
00394 
00395 template<class Scalar>
00396 void ExplicitRKStepper<Scalar>::addPoints(
00397     const Array<Scalar>& time_vec
00398     ,const Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >& x_vec
00399     ,const Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >& xdot_vec
00400     )
00401 {
00402   TEST_FOR_EXCEPTION(true,std::logic_error,"Error, addPoints is not implemented for ExplicitRKStepper at this time.\n");
00403 }
00404 
00405 template<class Scalar>
00406 void ExplicitRKStepper<Scalar>::getPoints(
00407     const Array<Scalar>& time_vec
00408     ,Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >* x_vec
00409     ,Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec
00410     ,Array<ScalarMag>* accuracy_vec) const
00411 {
00412   TEST_FOR_EXCEPTION(true,std::logic_error,"Error, getPoints is not implemented for ExplicitRKStepper at this time.\n");
00413 }
00414 
00415 template<class Scalar>
00416 TimeRange<Scalar> ExplicitRKStepper<Scalar>::getTimeRange() const
00417 {
00418   return invalidTimeRange<Scalar>();
00419 }
00420 
00421 template<class Scalar>
00422 void ExplicitRKStepper<Scalar>::getNodes(Array<Scalar>* time_vec) const
00423 {
00424   TEST_FOR_EXCEPTION(true,std::logic_error,"Error, getNodes is not implemented for ExplicitRKStepper at this time.\n");
00425 }
00426 
00427 template<class Scalar>
00428 void ExplicitRKStepper<Scalar>::removeNodes(Array<Scalar>& time_vec) 
00429 {
00430   TEST_FOR_EXCEPTION(true,std::logic_error,"Error, removeNodes is not implemented for ExplicitRKStepper at this time.\n");
00431 }
00432 
00433 template<class Scalar>
00434 int ExplicitRKStepper<Scalar>::getOrder() const
00435 {
00436   return(4);
00437 }
00438 
00439 template <class Scalar>
00440 void ExplicitRKStepper<Scalar>::setParameterList(Teuchos::RCP<Teuchos::ParameterList> const& paramList)
00441 {
00442   parameterList_ = paramList;
00443   int outputLevel = parameterList_->get( "outputLevel", int(-1) );
00444   outputLevel = std::min(std::max(outputLevel,-1),4);
00445   this->setVerbLevel(static_cast<Teuchos::EVerbosityLevel>(outputLevel));
00446 }
00447 
00448 template <class Scalar>
00449 Teuchos::RCP<Teuchos::ParameterList> ExplicitRKStepper<Scalar>::getNonconstParameterList()
00450 {
00451   return(parameterList_);
00452 }
00453 
00454 template <class Scalar>
00455 Teuchos::RCP<Teuchos::ParameterList> ExplicitRKStepper<Scalar>::unsetParameterList()
00456 {
00457   Teuchos::RCP<Teuchos::ParameterList> temp_param_list = parameterList_;
00458   parameterList_ = Teuchos::null;
00459   return(temp_param_list);
00460 }
00461 
00462 template<class Scalar>
00463 void ExplicitRKStepper<Scalar>::setModel(const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > &model)
00464 {
00465   TEST_FOR_EXCEPT(model == Teuchos::null)
00466   model_ = model;
00467 }
00468 
00469 template<class Scalar>
00470 Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >
00471 ExplicitRKStepper<Scalar>::getModel() const
00472 {
00473   return model_;
00474 }
00475 
00476 } // namespace Rythmos
00477 
00478 #endif //Rythmos_ExplicitRK_STEPPER_H
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends

Generated on Tue Oct 20 10:24:09 2009 for Rythmos - Transient Integration for Differential Equations by  doxygen 1.6.1