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