Rythmos_ExplicitRKStepper_def.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_DEF_H
00030 #define Rythmos_ExplicitRK_STEPPER_DEF_H
00031 
00032 #include "Rythmos_ExplicitRKStepper_decl.hpp"
00033 
00034 #include "Rythmos_RKButcherTableau.hpp"
00035 #include "Rythmos_RKButcherTableauHelpers.hpp"
00036 #include "Rythmos_RKButcherTableauBuilder.hpp"
00037 #include "Rythmos_StepperHelpers.hpp"
00038 #include "Rythmos_LinearInterpolator.hpp"
00039 #include "Rythmos_InterpolatorBaseHelpers.hpp"
00040 
00041 #include "Teuchos_VerboseObjectParameterListHelpers.hpp"
00042 
00043 #include "Thyra_ModelEvaluatorHelpers.hpp"
00044 
00045 namespace Rythmos {
00046 
00047 // Non-member constructors
00048 template<class Scalar>
00049 RCP<ExplicitRKStepper<Scalar> > explicitRKStepper()
00050 {
00051   RCP<ExplicitRKStepper<Scalar> > stepper = rcp(new ExplicitRKStepper<Scalar>());
00052   return stepper;
00053 }
00054 
00055 template<class Scalar>
00056 RCP<ExplicitRKStepper<Scalar> > explicitRKStepper(
00057     const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > &model 
00058     )
00059 {
00060   RCP<RKButcherTableauBase<Scalar> > rkbt = createRKBT<Scalar>("Explicit 4 Stage");
00061   //RCP<RKButcherTableauBase<Scalar> > rkbt = rcp(new Explicit4Stage4thOrder_RKBT<Scalar>());
00062   RCP<ExplicitRKStepper<Scalar> > stepper = rcp(new ExplicitRKStepper<Scalar>());
00063   stepper->setModel(model);
00064   stepper->setRKButcherTableau(rkbt);
00065   return stepper;
00066 }
00067 
00068 template<class Scalar>
00069 RCP<ExplicitRKStepper<Scalar> > explicitRKStepper(
00070     const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > &model,
00071     const RCP<const RKButcherTableauBase<Scalar> > &rkbt 
00072     )
00073 {
00074   RCP<ExplicitRKStepper<Scalar> > stepper = rcp(new ExplicitRKStepper<Scalar>());
00075   stepper->setModel(model);
00076   stepper->setRKButcherTableau(rkbt);
00077   return stepper;
00078 }
00079 
00080 template<class Scalar>
00081 ExplicitRKStepper<Scalar>::ExplicitRKStepper()
00082 {
00083   this->defaultInitializeAll_();
00084   Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
00085   out->precision(15);
00086   erkButcherTableau_ = rKButcherTableau<Scalar>();
00087   numSteps_ = 0;
00088 }
00089 
00090 template<class Scalar>
00091 void ExplicitRKStepper<Scalar>::defaultInitializeAll_()
00092 {
00093   model_ = Teuchos::null;
00094   solution_vector_ = Teuchos::null;
00095   solution_vector_old_ = Teuchos::null;
00096   //k_vector_;
00097   ktemp_vector_ = Teuchos::null;
00098   //basePoint_;
00099   erkButcherTableau_ = Teuchos::null;
00100   t_ = ST::nan();
00101   t_old_ = ST::nan();
00102   dt_ = ST::nan();
00103   numSteps_ = -1;
00104   parameterList_ = Teuchos::null;
00105   isInitialized_ = false;
00106   haveInitialCondition_ = false;
00107 }
00108 
00109 template<class Scalar>
00110 void ExplicitRKStepper<Scalar>::setRKButcherTableau(const RCP<const RKButcherTableauBase<Scalar> > &rkbt)
00111 {
00112   TEUCHOS_ASSERT( !is_null(rkbt) );
00113   validateERKButcherTableau(*rkbt);
00114   int numStages_old = erkButcherTableau_->numStages();
00115   int numStages_new = rkbt->numStages();
00116   TEST_FOR_EXCEPTION( numStages_new == 0, std::logic_error,
00117       "Error!  The Runge-Kutta Butcher tableau has no stages!"
00118       );
00119   if (!is_null(model_)) {
00120     int numNewStages = numStages_new - numStages_old;
00121     if ( numNewStages > 0 ) {
00122       k_vector_.reserve(numStages_new);
00123       for (int i=0 ; i<numNewStages ; ++i) {
00124         k_vector_.push_back(Thyra::createMember(model_->get_f_space()));
00125       }
00126     }
00127   }
00128   erkButcherTableau_ = rkbt;
00129 }
00130 
00131 template<class Scalar>
00132 RCP<const RKButcherTableauBase<Scalar> > ExplicitRKStepper<Scalar>::getRKButcherTableau() const
00133 {
00134   return erkButcherTableau_;
00135 }
00136 
00137 template<class Scalar>
00138 void ExplicitRKStepper<Scalar>::initialize_()
00139 {
00140   if (!isInitialized_) {
00141     TEUCHOS_ASSERT( !is_null(model_) );
00142     TEUCHOS_ASSERT( !is_null(erkButcherTableau_) );
00143     TEUCHOS_ASSERT( haveInitialCondition_ );
00144     TEST_FOR_EXCEPTION( erkButcherTableau_->numStages() == 0, std::logic_error,
00145         "Error!  The Runge-Kutta Butcher tableau has no stages!"
00146         );
00147     ktemp_vector_ = Thyra::createMember(model_->get_f_space());
00148     // Initialize the stage vectors
00149     int numStages = erkButcherTableau_->numStages();
00150     k_vector_.reserve(numStages);
00151     for (int i=0 ; i<numStages ; ++i) {
00152       k_vector_.push_back(Thyra::createMember(model_->get_f_space()));
00153     }
00154   }
00155 #ifdef RYTHMOS_DEBUG
00156   THYRA_ASSERT_VEC_SPACES(
00157     "Rythmos::ExplicitRKStepper::initialize_(...)",
00158     *solution_vector_->space(), *model_->get_x_space() );
00159 #endif // RYTHMOS_DEBUG
00160   isInitialized_ = true;
00161 }
00162 
00163 
00164 template<class Scalar>
00165 ExplicitRKStepper<Scalar>::~ExplicitRKStepper()
00166 {
00167 }
00168 
00169 template<class Scalar>
00170 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > ExplicitRKStepper<Scalar>::get_x_space() const
00171 {
00172   TEUCHOS_ASSERT( !is_null(model_) );
00173   return(model_->get_x_space());
00174 }
00175 
00176 template<class Scalar>
00177 Scalar ExplicitRKStepper<Scalar>::takeStep(Scalar dt, StepSizeType flag)
00178 {
00179   typedef typename Thyra::ModelEvaluatorBase::InArgs<Scalar>::ScalarMag ScalarMag;
00180   this->initialize_();
00181 #ifdef RYTHMOS_DEBUG
00182     TEST_FOR_EXCEPTION( flag == STEP_TYPE_VARIABLE, std::logic_error,
00183         "Error!  ExplicitRKStepper does not support variable time steps at this time."
00184         );
00185 #endif // RYTHMOS_DEBUG
00186   if ((flag == STEP_TYPE_VARIABLE) || (dt == ST::zero())) {
00187     return(Scalar(-ST::one()));
00188   }
00189   // Store old solution & old time
00190   V_V(&*solution_vector_old_, *solution_vector_);
00191   t_old_ = t_;
00192 
00193   dt_ = dt;
00194 
00195   int stages = erkButcherTableau_->numStages();
00196   Teuchos::SerialDenseMatrix<int,Scalar> A = erkButcherTableau_->A();
00197   Teuchos::SerialDenseVector<int,Scalar> b = erkButcherTableau_->b();
00198   Teuchos::SerialDenseVector<int,Scalar> c = erkButcherTableau_->c();
00199   // Compute stage solutions
00200   for (int s=0 ; s < stages ; ++s) {
00201     Thyra::assign(&*ktemp_vector_, *solution_vector_); // ktemp = solution_vector
00202     for (int j=0 ; j < s ; ++j) { // assuming Butcher matix is strictly lower triangular
00203       if (A(s,j) != ST::zero()) {
00204         Thyra::Vp_StV(&*ktemp_vector_, A(s,j), *k_vector_[j]); // ktemp = ktemp + a_{s+1,j+1}*k_{j+1}
00205       }
00206     }
00207     ScalarMag ts = t_ + c(s)*dt;
00208     eval_model_explicit<Scalar>(*model_,basePoint_,*ktemp_vector_,ts,Teuchos::outArg(*k_vector_[s]));
00209     Thyra::Vt_S(&*k_vector_[s],dt); // k_s = k_s*dt
00210   } 
00211   // Sum for solution:
00212   for (int s=0 ; s < stages ; ++s) {
00213     if (b(s) != ST::zero()) {
00214       Thyra::Vp_StV(&*solution_vector_, b(s), *k_vector_[s]); // solution_vector += b_{s+1}*k_{s+1}
00215     }
00216   }
00217 
00218   // update current time:
00219   t_ = t_ + dt;
00220 
00221   numSteps_++;
00222 
00223   return(dt);
00224 }
00225 
00226 template<class Scalar>
00227 const StepStatus<Scalar> ExplicitRKStepper<Scalar>::getStepStatus() const
00228 {
00229   StepStatus<Scalar> stepStatus;
00230 
00231   if (!haveInitialCondition_) {
00232     stepStatus.stepStatus = STEP_STATUS_UNINITIALIZED;
00233   } else if (numSteps_ == 0) {
00234     stepStatus.stepStatus = STEP_STATUS_UNKNOWN;
00235     stepStatus.order = erkButcherTableau_->order();
00236     stepStatus.time = t_;
00237     stepStatus.solution = solution_vector_;
00238   } else {
00239     stepStatus.stepStatus = STEP_STATUS_CONVERGED;
00240     stepStatus.stepSize = dt_;
00241     stepStatus.order = erkButcherTableau_->order();
00242     stepStatus.time = t_;
00243     stepStatus.solution = solution_vector_;
00244   }
00245 
00246   return(stepStatus);
00247 }
00248 
00249 template<class Scalar>
00250 void ExplicitRKStepper<Scalar>::describe(
00251     Teuchos::FancyOStream &out,
00252     const Teuchos::EVerbosityLevel verbLevel
00253     ) const
00254 {
00255   if ( (static_cast<int>(verbLevel) == static_cast<int>(Teuchos::VERB_DEFAULT) ) ||
00256        (static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW)     )
00257      ) {
00258     out << this->description() << "::describe" << std::endl;
00259     out << "model = " << model_->description() << std::endl;
00260     out << erkButcherTableau_->numStages() << " stage Explicit RK method" << std::endl;
00261   } else if (static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW)) {
00262     out << "solution_vector = " << std::endl;
00263     out << Teuchos::describe(*solution_vector_,verbLevel);
00264   } else if (static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_MEDIUM)) {
00265   } else if (static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_HIGH)) {
00266     out << "model = " << std::endl;
00267     out << Teuchos::describe(*model_,verbLevel);
00268     int stages = erkButcherTableau_->numStages();
00269     for (int i=0 ; i<stages ; ++i) {
00270       out << "k_vector[" << i << "] = " << std::endl;
00271       out << Teuchos::describe(*k_vector_[i],verbLevel);
00272     }
00273     out << "ktemp_vector = " << std::endl;
00274     out << Teuchos::describe(*ktemp_vector_,verbLevel);
00275     out << "ERK Butcher Tableau A matrix: " << erkButcherTableau_->A() << std::endl;
00276     out << "ERK Butcher Tableau b vector: " << erkButcherTableau_->b() << std::endl;
00277     out << "ERK Butcher Tableau c vector: " << erkButcherTableau_->c() << std::endl;
00278     out << "t = " << t_ << std::endl;
00279   }
00280 }
00281 
00282 template<class Scalar>
00283 void ExplicitRKStepper<Scalar>::addPoints(
00284     const Array<Scalar>& time_vec
00285     ,const Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >& x_vec
00286     ,const Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >& xdot_vec
00287     )
00288 {
00289   TEST_FOR_EXCEPTION(true,std::logic_error,"Error, addPoints is not implemented for ExplicitRKStepper at this time.\n");
00290 }
00291 
00292 template<class Scalar>
00293 TimeRange<Scalar> ExplicitRKStepper<Scalar>::getTimeRange() const
00294 {
00295   if (!haveInitialCondition_) {
00296     return(invalidTimeRange<Scalar>());
00297   } else {
00298     return(TimeRange<Scalar>(t_old_,t_));
00299   }
00300 }
00301 
00302 template<class Scalar>
00303 void ExplicitRKStepper<Scalar>::getPoints(
00304   const Array<Scalar>& time_vec,
00305   Array<RCP<const Thyra::VectorBase<Scalar> > >* x_vec,
00306   Array<RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec,
00307   Array<ScalarMag>* accuracy_vec
00308   ) const
00309 {
00310   TEUCHOS_ASSERT( haveInitialCondition_ );
00311   using Teuchos::constOptInArg;
00312   using Teuchos::null;
00313   defaultGetPoints<Scalar>(
00314       t_old_,constOptInArg(*solution_vector_old_),null,
00315       t_,constOptInArg(*solution_vector_),null,
00316       time_vec,ptr(x_vec),ptr(xdot_vec),ptr(accuracy_vec),
00317       null
00318       );
00319 }
00320 
00321 template<class Scalar>
00322 void ExplicitRKStepper<Scalar>::getNodes(Array<Scalar>* time_vec) const
00323 {
00324   TEUCHOS_ASSERT( time_vec != NULL );
00325   time_vec->clear();
00326   if (!haveInitialCondition_) {
00327     return;
00328   }
00329   time_vec->push_back(t_old_);
00330   if (t_ != t_old_) {
00331     time_vec->push_back(t_);
00332   }
00333 }
00334 
00335 template<class Scalar>
00336 void ExplicitRKStepper<Scalar>::removeNodes(Array<Scalar>& time_vec) 
00337 {
00338   TEST_FOR_EXCEPTION(true,std::logic_error,"Error, removeNodes is not implemented for ExplicitRKStepper at this time.\n");
00339 }
00340 
00341 template<class Scalar>
00342 int ExplicitRKStepper<Scalar>::getOrder() const
00343 {
00344   return(erkButcherTableau_->order());
00345 }
00346 
00347 template <class Scalar>
00348 void ExplicitRKStepper<Scalar>::setParameterList(Teuchos::RCP<Teuchos::ParameterList> const& paramList)
00349 {
00350   TEST_FOR_EXCEPT(is_null(paramList));
00351   paramList->validateParametersAndSetDefaults(*this->getValidParameters());
00352   parameterList_ = paramList;
00353   Teuchos::readVerboseObjectSublist(&*parameterList_,this);
00354 }
00355 
00356 template <class Scalar>
00357 Teuchos::RCP<Teuchos::ParameterList> ExplicitRKStepper<Scalar>::getNonconstParameterList()
00358 {
00359   return(parameterList_);
00360 }
00361 
00362 template <class Scalar>
00363 Teuchos::RCP<Teuchos::ParameterList> ExplicitRKStepper<Scalar>::unsetParameterList()
00364 {
00365   Teuchos::RCP<Teuchos::ParameterList> temp_param_list = parameterList_;
00366   parameterList_ = Teuchos::null;
00367   return(temp_param_list);
00368 }
00369 
00370 template<class Scalar>
00371 RCP<const Teuchos::ParameterList>
00372 ExplicitRKStepper<Scalar>::getValidParameters() const
00373 {
00374   using Teuchos::ParameterList;
00375   static RCP<const ParameterList> validPL;
00376   if (is_null(validPL)) {
00377     RCP<ParameterList> pl = Teuchos::parameterList();
00378     Teuchos::setupVerboseObjectSublist(&*pl);
00379     validPL = pl;
00380   }
00381   return validPL;
00382 }
00383 
00384 template<class Scalar>
00385 void ExplicitRKStepper<Scalar>::setModel(const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > &model)
00386 {
00387   TEST_FOR_EXCEPT( is_null(model) );
00388   TEST_FOR_EXCEPT( !is_null(model_) ); // For now you can only call this once.
00389   assertValidModel( *this, *model );
00390   model_ = model;
00391 }
00392 
00393 template<class Scalar>
00394 Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >
00395 ExplicitRKStepper<Scalar>::getModel() const
00396 {
00397   return model_;
00398 }
00399 
00400 template<class Scalar>
00401 void ExplicitRKStepper<Scalar>::setInitialCondition(
00402   const Thyra::ModelEvaluatorBase::InArgs<Scalar> &initialCondition
00403   )
00404 {
00405 
00406   typedef Teuchos::ScalarTraits<Scalar> ST;
00407   typedef Thyra::ModelEvaluatorBase MEB;
00408 
00409   basePoint_ = initialCondition;
00410 
00411   // x
00412 
00413   RCP<const Thyra::VectorBase<Scalar> >
00414     x_init = initialCondition.get_x();
00415 
00416 #ifdef RYTHMOS_DEBUG
00417   TEST_FOR_EXCEPTION(
00418     is_null(x_init), std::logic_error,
00419     "Error, if the client passes in an intial condition to setInitialCondition(...),\n"
00420     "then x can not be null!" );
00421 #endif
00422 
00423   solution_vector_ = x_init->clone_v();
00424   solution_vector_old_ = x_init->clone_v();
00425 
00426   // t
00427   
00428   t_ = initialCondition.get_t();
00429   t_old_ = t_;
00430 
00431   haveInitialCondition_ = true;
00432 
00433 }
00434 
00435 template<class Scalar>
00436 bool ExplicitRKStepper<Scalar>::supportsCloning() const
00437 {
00438   return true;
00439 }
00440 
00441 template<class Scalar>
00442 RCP<StepperBase<Scalar> > ExplicitRKStepper<Scalar>::cloneStepperAlgorithm() const
00443 {
00444   // Just use the interface to clone the algorithm in a basically
00445   // uninitialized state
00446   RCP<ExplicitRKStepper<Scalar> >
00447     stepper = Teuchos::rcp(new ExplicitRKStepper<Scalar>());
00448 
00449   if (!is_null(model_)) {
00450     stepper->setModel(model_); // Shallow copy is okay!
00451   }
00452 
00453   if (!is_null(erkButcherTableau_)) {
00454     // 06/16/09 tscoffe:  should we clone the RKBT here?
00455     stepper->setRKButcherTableau(erkButcherTableau_);
00456   }
00457 
00458   if (!is_null(parameterList_)) {
00459     stepper->setParameterList(Teuchos::parameterList(*parameterList_));
00460   }
00461 
00462   return stepper;
00463 
00464 }
00465     
00466 
00467 
00468 // 
00469 // Explicit Instantiation macro
00470 //
00471 // Must be expanded from within the Rythmos namespace!
00472 //
00473 
00474 #define RYTHMOS_EXPLICIT_RK_STEPPER_INSTANT(SCALAR) \
00475   \
00476   template class ExplicitRKStepper< SCALAR >; \
00477   \
00478   template RCP< ExplicitRKStepper< SCALAR > > \
00479   explicitRKStepper();  \
00480   \
00481   template RCP< ExplicitRKStepper< SCALAR > > \
00482   explicitRKStepper( \
00483     const RCP<const Thyra::ModelEvaluator< SCALAR > > &model \
00484       ); \
00485   \
00486   template RCP< ExplicitRKStepper< SCALAR > > \
00487   explicitRKStepper( \
00488     const RCP<const Thyra::ModelEvaluator< SCALAR > > &model, \
00489     const RCP<const RKButcherTableauBase< SCALAR > > &rkbt \
00490       ); \
00491    
00492 } // namespace Rythmos
00493 
00494 #endif //Rythmos_ExplicitRK_STEPPER_DEF_H
00495 

Generated on Wed May 12 21:25:42 2010 for Rythmos - Transient Integration for Differential Equations by  doxygen 1.4.7