Rythmos_ExplicitRKStepper.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_ExplicitRK_STEPPER_H
00030 #define Rythmos_ExplicitRK_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 
00038 namespace Rythmos {
00039 
00041 template<class Scalar>
00042 class ExplicitRKStepper : virtual public Stepper<Scalar>
00043 {
00044   public:
00045 
00046     typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType ScalarMag;
00047     
00049     ExplicitRKStepper();
00050     ExplicitRKStepper(const Teuchos::RefCountPtr<const Thyra::ModelEvaluator<Scalar> > &model_);
00051     
00053     ~ExplicitRKStepper();
00054 
00056     Scalar TakeStep(Scalar dt);
00057    
00059     Scalar TakeStep();
00060 
00062     Teuchos::RefCountPtr<const Thyra::VectorBase<Scalar> > get_solution() const;
00063     
00065     std::string description() const;
00066 
00068     std::ostream& describe(
00069       std::ostream                &out
00070       ,const Teuchos::EVerbosityLevel      verbLevel
00071       ,const std::string          leadingIndent
00072       ,const std::string          indentSpacer
00073       ) const;
00074 
00077     bool SetPoints(
00078       const std::vector<Scalar>& time_list
00079       ,const std::vector<Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > >& x_list
00080       ,const std::vector<Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > >& xdot_list);
00081     
00083     bool GetPoints(
00084       const std::vector<Scalar>& time_list
00085       ,std::vector<Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > >* x_list
00086       ,std::vector<Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > >* xdot_list
00087       ,std::vector<ScalarMag>* accuracy_list) const;
00088 
00090     bool SetRange(
00091       const Scalar& time_lower
00092       ,const Scalar& time_upper
00093       ,const InterpolationBuffer<Scalar> & IB);
00094 
00096     bool GetNodes(std::vector<Scalar>* time_list) const;
00097 
00099     bool RemoveNodes(std::vector<Scalar>& time_list) const;
00100 
00102     int GetOrder() const;
00103 
00104   private:
00105 
00106     Teuchos::RefCountPtr<const Thyra::ModelEvaluator<Scalar> > model_;
00107     Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > solution_vector_;
00108     std::vector<Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > > k_vector_;
00109     Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > ktemp_vector_;
00110 
00111     int stages_; // Number of stages of RK
00112     std::vector<std::vector<Scalar> > b_A; // Butcher tableau A matrix
00113     std::vector<Scalar> b_b; // Butcher b vector
00114     std::vector<Scalar> b_c; // Butcher c vector
00115 
00116     Scalar t_;
00117 
00118 };
00119 
00120 template<class Scalar>
00121 ExplicitRKStepper<Scalar>::ExplicitRKStepper(const Teuchos::RefCountPtr<const Thyra::ModelEvaluator<Scalar> > &model)
00122 {
00123   typedef Teuchos::ScalarTraits<Scalar> ST;
00124   model_ = model;
00125   t_ = ST::zero();
00126   solution_vector_ = model_->getNominalValues().get_x()->clone_v();
00127   stages_ = 4; // 4 stage ERK
00128   k_vector_.reserve(stages_);
00129   Teuchos::RefCountPtr<const Thyra::VectorSpaceBase<Scalar> >
00130     f_space = model_->get_f_space();
00131   for (int i=0 ; i<stages_ ; ++i)
00132   {
00133     k_vector_.push_back(Thyra::createMember(f_space));
00134   }
00135   ktemp_vector_ = Thyra::createMember(f_space);
00136 
00137   // initialize the Butcher tableau and its vectors
00138   Scalar zero = ST::zero();
00139   b_A.reserve(stages_);
00140   b_b.reserve(stages_);
00141   b_c.reserve(stages_);
00142 
00143   // fill b with zeros
00144   b_b.push_back(zero); 
00145   b_b.push_back(zero); 
00146   b_b.push_back(zero);
00147   b_b.push_back(zero);
00148   // fill c with zeros
00149   b_c.push_back(zero); 
00150   b_c.push_back(zero); 
00151   b_c.push_back(zero);
00152   b_c.push_back(zero);
00153   // fill A with zeros
00154   for (int i=0 ; i<stages_ ; ++i)
00155   {
00156     b_A.push_back(b_b);
00157   }
00158 
00159   // Runge-Kutta methods in Butcher tableau form:
00160   //  c | A    A = sxs matrix
00161   //   -|---   b = s vector
00162   //      b'   c = s vector
00163   
00164 /*
00165   // 3/8 Rule Runge-Kutta Method: (not implemented yet)
00166   // c = [  0  1/3 2/3  1  ]'
00167   // A = [  0              ]
00168   //     [ 1/3  0          ]
00169   //     [-1/3  1   0      ]
00170   //     [  1  -1   1   0  ]
00171   // b = [ 1/8 3/8 3/8 1/8 ]'
00172   Scalar one = ST::one();
00173   Scalar one_third    = Scalar(ST::one()/(3*ST::one()));
00174   Scalar two_third    = Scalar(2*ST::one()/(3*ST::one()));
00175   Scalar one_eighth   = Scalar(ST::one()/(8*ST::one()));
00176   Scalar three_eighth = Scalar(3*ST::one()/(8*ST::one()));
00177 
00178   // fill b_A
00179   b_A[0][0] = zero;
00180   b_A[0][1] = zero;
00181   b_A[0][2] = zero;
00182   b_A[0][3] = zero;
00183 
00184   b_A[1][0] = one_third;
00185   b_A[1][1] = zero;
00186   b_A[1][2] = zero;
00187   b_A[1][3] = zero;
00188 
00189   b_A[2][0] = Scalar(-one_third);
00190   b_A[2][1] = one;
00191   b_A[2][2] = zero;
00192   b_A[2][3] = zero;
00193 
00194   b_A[3][0] = one;
00195   b_A[3][1] = Scalar(-one);
00196   b_A[3][2] = one;
00197   b_A[3][3] = zero;
00198 
00199   // fill b_b
00200   b_b[0] = one_eighth;
00201   b_b[1] = three_eighth;
00202   b_b[2] = three_eighth;
00203   b_b[3] = one_eighth;
00204   
00205   // fill b_c
00206   b_c[0] = zero;
00207   b_c[1] = one_third;
00208   b_c[2] = two_third;
00209   b_c[3] = one;
00210 */
00211 
00212 
00213   // "The" Runge-Kutta Method: (implemented below)
00214   // c = [  0  1/2 1/2  1  ]'
00215   // A = [  0              ] 
00216   //     [ 1/2  0          ]
00217   //     [  0  1/2  0      ]
00218   //     [  0   0   1   0  ]
00219   // b = [ 1/6 1/3 1/3 1/6 ]'
00220   Scalar one = ST::one();
00221   Scalar onehalf = ST::one()/(2*ST::one());
00222   Scalar onesixth = ST::one()/(6*ST::one());
00223   Scalar onethird = ST::one()/(3*ST::one());
00224 
00225   // fill b_A
00226   b_A[0][0] = zero;
00227   b_A[0][1] = zero;
00228   b_A[0][2] = zero;
00229   b_A[0][3] = zero;
00230 
00231   b_A[1][0] = onehalf;
00232   b_A[1][1] = zero;
00233   b_A[1][2] = zero;
00234   b_A[1][3] = zero;
00235 
00236   b_A[2][0] = zero;
00237   b_A[2][1] = onehalf;
00238   b_A[2][2] = zero;
00239   b_A[2][3] = zero;
00240 
00241   b_A[3][0] = zero;
00242   b_A[3][1] = zero;
00243   b_A[3][2] = one;
00244   b_A[3][3] = zero;
00245 
00246   // fill b_b
00247   b_b[0] = onesixth;
00248   b_b[1] = onethird;
00249   b_b[2] = onethird;
00250   b_b[3] = onesixth;
00251   
00252   // fill b_c
00253   b_c[0] = zero;
00254   b_c[1] = onehalf;
00255   b_c[2] = onehalf;
00256   b_c[3] = one;
00257 
00258 }
00259 
00260 template<class Scalar>
00261 ExplicitRKStepper<Scalar>::ExplicitRKStepper()
00262 {
00263 }
00264 
00265 template<class Scalar>
00266 ExplicitRKStepper<Scalar>::~ExplicitRKStepper()
00267 {
00268 }
00269 
00270 template<class Scalar>
00271 Scalar ExplicitRKStepper<Scalar>::TakeStep()
00272 {
00273   // print something out about this method not supporting automatic variable step-size
00274   typedef Teuchos::ScalarTraits<Scalar> ST;
00275   return(-ST::one());
00276 }
00277 
00278 template<class Scalar>
00279 Scalar ExplicitRKStepper<Scalar>::TakeStep(Scalar dt)
00280 {
00281   typedef Teuchos::ScalarTraits<Scalar> ST;
00282   typedef typename Thyra::ModelEvaluatorBase::InArgs<Scalar>::ScalarMag ScalarMag;
00283 
00284   // Compute stage solutions
00285   for (int s=0 ; s < stages_ ; ++s)
00286   {
00287     Thyra::assign(&*ktemp_vector_, *solution_vector_); // ktemp = solution_vector
00288     for (int j=0 ; j < s ; ++j) // assuming Butcher matix is strictly lower triangular
00289     {
00290       if (b_A[s][j] != ST::zero())
00291         Thyra::Vp_StV(&*ktemp_vector_, b_A[s][j], *k_vector_[j]); // ktemp = ktemp + a_{s+1,j+1}*k_{j+1}
00292     }
00293     ScalarMag ts = t_ + b_c[s]*dt;
00294     Thyra::eval_f<Scalar>(*model_,*ktemp_vector_,ts,&*k_vector_[s]);
00295     Thyra::Vt_S(&*k_vector_[s],dt); // k_s = k_s*dt
00296   
00297   } 
00298   // Sum for solution:
00299   for (int s=0 ; s < stages_ ; ++s)
00300   {
00301     if (b_b[s] != ST::zero())
00302       Thyra::Vp_StV(&*solution_vector_, b_b[s], *k_vector_[s]); // solution_vector += b_{s+1}*k_{s+1}
00303   }
00304 
00305   // update current time:
00306   t_ = t_ + dt;
00307 
00308   return(dt);
00309 }
00310 
00311 template<class Scalar>
00312 Teuchos::RefCountPtr<const Thyra::VectorBase<Scalar> > ExplicitRKStepper<Scalar>::get_solution() const
00313 {
00314   return(solution_vector_);
00315 }
00316 
00317 template<class Scalar>
00318 std::string ExplicitRKStepper<Scalar>::description() const
00319 {
00320   std::string name = "Rythmos::ExplicitRKStepper";
00321   return(name);
00322 }
00323 
00324 template<class Scalar>
00325 std::ostream& ExplicitRKStepper<Scalar>::describe(
00326       std::ostream                &out
00327       ,const Teuchos::EVerbosityLevel      verbLevel
00328       ,const std::string          leadingIndent
00329       ,const std::string          indentSpacer
00330       ) const
00331 {
00332   if (verbLevel == Teuchos::VERB_EXTREME)
00333   {
00334     out << description() << "::describe" << std::endl;
00335     out << "model_ = " << std::endl;
00336     out << model_->describe(out,verbLevel,leadingIndent,indentSpacer) << std::endl;
00337     out << "solution_vector_ = " << std::endl;
00338     out << solution_vector_->describe(out,verbLevel,leadingIndent,indentSpacer) << std::endl;
00339     for (int i=0 ; i<stages_ ; ++i)
00340     {
00341       out << "k_vector_[" << i << "] = " << std::endl;
00342       out << k_vector_[i]->describe(out,verbLevel,leadingIndent,indentSpacer) << std::endl;
00343     }
00344     out << "ktemp_vector_ = " << std::endl;
00345     out << ktemp_vector_->describe(out,verbLevel,leadingIndent,indentSpacer) << std::endl;
00346     for (int i=0 ; i<stages_ ; ++i)
00347       for (int j=0 ; j<stages_ ; ++j)
00348         out << "b_A[" << i << "][" << j << "] = " << b_A[i][j] << std::endl;
00349     for (int i=0 ; i<stages_ ; ++i)
00350       out << "b_b[" << i << "] = " << b_b[i] << std::endl;
00351     for (int i=0 ; i<stages_ ; ++i)
00352       out << "b_c[" << i << "] = " << b_c[i] << std::endl;
00353     out << "t_ = " << t_ << std::endl;
00354   }
00355   return(out);
00356 }
00357 
00358 template<class Scalar>
00359 bool ExplicitRKStepper<Scalar>::SetPoints(
00360     const std::vector<Scalar>& time_list
00361     ,const std::vector<Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > >& x_list
00362     ,const std::vector<Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > >& xdot_list)
00363 {
00364   return(false);
00365 }
00366 
00367 template<class Scalar>
00368 bool ExplicitRKStepper<Scalar>::GetPoints(
00369     const std::vector<Scalar>& time_list
00370     ,std::vector<Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > >* x_list
00371     ,std::vector<Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > >* xdot_list
00372     ,std::vector<ScalarMag>* accuracy_list) const
00373 {
00374   return(false);
00375 }
00376 
00377 template<class Scalar>
00378 bool ExplicitRKStepper<Scalar>::SetRange(
00379     const Scalar& time_lower
00380     ,const Scalar& time_upper
00381     ,const InterpolationBuffer<Scalar>& IB)
00382 {
00383   return(false);
00384 }
00385 
00386 template<class Scalar>
00387 bool ExplicitRKStepper<Scalar>::GetNodes(std::vector<Scalar>* time_list) const
00388 {
00389   return(false);
00390 }
00391 
00392 template<class Scalar>
00393 bool ExplicitRKStepper<Scalar>::RemoveNodes(std::vector<Scalar>& time_list) const
00394 {
00395   return(false);
00396 }
00397 
00398 template<class Scalar>
00399 int ExplicitRKStepper<Scalar>::GetOrder() const
00400 {
00401   return(4);
00402 }
00403 
00404 } // namespace Rythmos
00405 
00406 #endif //Rythmos_ExplicitRK_STEPPER_H

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