Rythmos - Transient Integration for Differential Equations Version of the Day
Rythmos_ExplicitTaylorPolynomialStepper.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_EXPLICIT_TAYLOR_POLYNOMIAL_STEPPER_H
00030 #define RYTHMOS_EXPLICIT_TAYLOR_POLYNOMIAL_STEPPER_H
00031 
00032 #include "Rythmos_StepperBase.hpp"
00033 #include "Rythmos_StepperHelpers.hpp"
00034 #include "Teuchos_RCP.hpp"
00035 #include "Teuchos_ParameterList.hpp"
00036 #include "Teuchos_VerboseObjectParameterListHelpers.hpp"
00037 #include "Thyra_VectorBase.hpp"
00038 #include "Thyra_ModelEvaluator.hpp"
00039 #include "Thyra_ModelEvaluatorHelpers.hpp"
00040 #include "Thyra_PolynomialVectorTraits.hpp"
00041 #include "RTOpPack_RTOpTHelpers.hpp"
00042 
00043 namespace Rythmos {
00044 
00046 
00053 RTOP_ROP_1_REDUCT_SCALAR( ROpLogNormInf,
00054   typename ScalarTraits<Scalar>::magnitudeType, // Reduction object type
00055   RTOpPack::REDUCT_TYPE_MAX // Reduction object reduction
00056   )
00057 {
00058   using Teuchos::as;
00059   typedef ScalarTraits<Scalar> ST;
00060   typedef typename ST::magnitudeType ScalarMag;
00061   const ScalarMag mag = std::log(as<ScalarMag>(1e-100) + ST::magnitude(v0));
00062   reduct = TEUCHOS_MAX( mag, reduct );
00063 }
00064 
00162 template<class Scalar>
00163 class ExplicitTaylorPolynomialStepper : virtual public StepperBase<Scalar>
00164 {
00165 public:
00166 
00168   typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType ScalarMag;
00169     
00171   ExplicitTaylorPolynomialStepper();
00172     
00174   ~ExplicitTaylorPolynomialStepper();
00175 
00177   RCP<const Thyra::VectorSpaceBase<Scalar> > get_x_space() const;
00178 
00180   void setModel(const RCP<const Thyra::ModelEvaluator<Scalar> >& model);
00181 
00183   void setNonconstModel(const RCP<Thyra::ModelEvaluator<Scalar> >& model);
00184 
00186   RCP<const Thyra::ModelEvaluator<Scalar> > getModel() const;
00187 
00189   RCP<Thyra::ModelEvaluator<Scalar> > getNonconstModel();
00190 
00192   void setInitialCondition(
00193     const Thyra::ModelEvaluatorBase::InArgs<Scalar> &initialCondition
00194     );
00195     
00197   Thyra::ModelEvaluatorBase::InArgs<Scalar> getInitialCondition() const;
00198     
00200   Scalar takeStep(Scalar dt, StepSizeType flag);
00201 
00203   const StepStatus<Scalar> getStepStatus() const;
00204 
00206 
00207   void setParameterList(RCP<Teuchos::ParameterList> const& paramList);
00208 
00210   RCP<Teuchos::ParameterList> getNonconstParameterList();
00211 
00213   RCP<Teuchos::ParameterList> unsetParameterList();
00214 
00216   RCP<const Teuchos::ParameterList> getValidParameters() const;
00217 
00219   std::string description() const;
00220 
00222   void describe(
00223     Teuchos::FancyOStream &out,
00224     const Teuchos::EVerbosityLevel verbLevel = Teuchos::Describable::verbLevel_default
00225     ) const;
00226 
00229   void addPoints(
00230     const Array<Scalar>& time_vec
00231     ,const Array<RCP<const Thyra::VectorBase<Scalar> > >& x_vec
00232     ,const Array<RCP<const Thyra::VectorBase<Scalar> > >& xdot_vec
00233     );
00234     
00236   void getPoints(
00237     const Array<Scalar>& time_vec
00238     ,Array<RCP<const Thyra::VectorBase<Scalar> > >* x_vec
00239     ,Array<RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec
00240     ,Array<ScalarMag>* accuracy_vec) const;
00241 
00243   void setRange(
00244     const TimeRange<Scalar>& range,
00245     const InterpolationBufferBase<Scalar> & IB
00246     );
00247 
00249   TimeRange<Scalar> getTimeRange() const;
00250 
00252   void getNodes(Array<Scalar>* time_vec) const;
00253 
00255   void removeNodes(Array<Scalar>& time_vec);
00256 
00258   int getOrder() const;
00259 
00260 private:
00261 
00263   void defaultInitializAll_();
00264 
00266   void computeTaylorSeriesSolution_();
00267 
00272   ScalarMag estimateLogRadius_();
00273 
00275   RCP<const Thyra::ModelEvaluator<Scalar> > model_;
00276 
00278   RCP<Teuchos::ParameterList> parameterList_;
00279 
00281   RCP<Thyra::VectorBase<Scalar> > x_vector_;
00282 
00284   RCP<Thyra::VectorBase<Scalar> > x_vector_old_;
00285 
00287   RCP<Thyra::VectorBase<Scalar> > x_dot_vector_;
00288 
00290   RCP<Thyra::VectorBase<Scalar> > x_dot_vector_old_;
00291 
00293   RCP<Thyra::VectorBase<Scalar> > f_vector_;
00294 
00296   RCP<Teuchos::Polynomial<Thyra::VectorBase<Scalar> > > x_poly_;
00297 
00299   RCP<Teuchos::Polynomial<Thyra::VectorBase<Scalar> > > f_poly_;
00300 
00302   Thyra::ModelEvaluatorBase::InArgs<Scalar> basePoint_;
00303 
00305   bool haveInitialCondition_;
00306 
00308   int numSteps_;
00309 
00311   Scalar t_;
00312 
00314   Scalar dt_;
00315 
00317   Scalar t_initial_;
00318 
00320   Scalar t_final_;
00321 
00323   ScalarMag local_error_tolerance_;
00324 
00326   Scalar min_step_size_;
00327 
00329   Scalar max_step_size_;
00330 
00332   unsigned int degree_;
00333 
00335   Scalar linc_;
00336 };
00337 
00338 
00340 template <typename Scalar>
00341 typename Teuchos::ScalarTraits<Scalar>::magnitudeType
00342 log_norm_inf(const Thyra::VectorBase<Scalar>& x)
00343 {
00344   ROpLogNormInf<Scalar> log_norm_inf_op;
00345   RCP<RTOpPack::ReductTarget> log_norm_inf_targ = 
00346     log_norm_inf_op.reduct_obj_create();
00347   const Thyra::VectorBase<Scalar>* vecs[] = { &x };
00348   Thyra::applyOp<Scalar>(log_norm_inf_op,1,vecs,0,
00349     (Thyra::VectorBase<Scalar>**)NULL,
00350     log_norm_inf_targ.get());
00351     
00352   return log_norm_inf_op(*log_norm_inf_targ);
00353 }
00354 
00355 
00356 // Non-member constructor
00357 template<class Scalar>
00358 RCP<ExplicitTaylorPolynomialStepper<Scalar> > explicitTaylorPolynomialStepper()
00359 {
00360   RCP<ExplicitTaylorPolynomialStepper<Scalar> > stepper = rcp(new ExplicitTaylorPolynomialStepper<Scalar>());
00361   return stepper;
00362 }
00363 
00364 
00365 template<class Scalar>
00366 ExplicitTaylorPolynomialStepper<Scalar>::ExplicitTaylorPolynomialStepper()
00367 {
00368   this->defaultInitializAll_();
00369   numSteps_ = 0;
00370 }
00371 
00372 
00373 template<class Scalar>
00374 ExplicitTaylorPolynomialStepper<Scalar>::~ExplicitTaylorPolynomialStepper()
00375 {
00376 }
00377 
00378 
00379 template<class Scalar>
00380 void ExplicitTaylorPolynomialStepper<Scalar>::defaultInitializAll_()
00381 {
00382   typedef Teuchos::ScalarTraits<Scalar> ST;
00383   Scalar nan = ST::nan();
00384   model_ = Teuchos::null;
00385   parameterList_ = Teuchos::null;
00386   x_vector_ = Teuchos::null;
00387   x_vector_old_ = Teuchos::null;
00388   x_dot_vector_ = Teuchos::null;
00389   x_dot_vector_old_ = Teuchos::null;
00390   f_vector_ = Teuchos::null;
00391   x_poly_ = Teuchos::null;
00392   f_poly_ = Teuchos::null;
00393   haveInitialCondition_ = false;
00394   numSteps_ = -1;
00395   t_ = nan;
00396   dt_ = nan;
00397   t_initial_ = nan;
00398   t_final_ = nan;
00399   local_error_tolerance_ = nan;
00400   min_step_size_ = nan;
00401   max_step_size_ = nan;
00402   degree_ = 0;
00403   linc_ = nan;
00404 }
00405 
00406 
00407 template<class Scalar>
00408 void ExplicitTaylorPolynomialStepper<Scalar>::setModel(
00409   const RCP<const Thyra::ModelEvaluator<Scalar> >& model
00410   )
00411 {
00412   TEUCHOS_TEST_FOR_EXCEPT( is_null(model) );
00413   assertValidModel( *this, *model );
00414     
00415   model_ = model;
00416   f_vector_ = Thyra::createMember(model_->get_f_space());
00417 }
00418 
00419 
00420 template<class Scalar>
00421 void ExplicitTaylorPolynomialStepper<Scalar>::setNonconstModel(
00422   const RCP<Thyra::ModelEvaluator<Scalar> >& model
00423   )
00424 {
00425   this->setModel(model); // TODO 09/09/09 tscoffe:  use ConstNonconstObjectContainer!
00426 }
00427 
00428 
00429 template<class Scalar>
00430 RCP<const Thyra::ModelEvaluator<Scalar> >
00431 ExplicitTaylorPolynomialStepper<Scalar>::getModel() const
00432 {
00433   return model_;
00434 }
00435 
00436 
00437 template<class Scalar>
00438 RCP<Thyra::ModelEvaluator<Scalar> >
00439 ExplicitTaylorPolynomialStepper<Scalar>::getNonconstModel() 
00440 {
00441   return Teuchos::null;
00442 }
00443 
00444 
00445 template<class Scalar>
00446 void ExplicitTaylorPolynomialStepper<Scalar>::setInitialCondition(
00447   const Thyra::ModelEvaluatorBase::InArgs<Scalar> &initialCondition
00448   )
00449 {
00450   typedef Teuchos::ScalarTraits<Scalar> ST;
00451   typedef Thyra::ModelEvaluatorBase MEB;
00452   basePoint_ = initialCondition;
00453   if (initialCondition.supports(MEB::IN_ARG_t)) {
00454     t_ = initialCondition.get_t();
00455   } else {
00456     t_ = ST::zero();
00457   }
00458   dt_ = ST::zero();
00459   x_vector_ = initialCondition.get_x()->clone_v();
00460   x_dot_vector_ = x_vector_->clone_v();
00461   x_vector_old_ = x_vector_->clone_v();
00462   x_dot_vector_old_ = x_dot_vector_->clone_v();
00463   haveInitialCondition_ = true;
00464 }
00465 
00466 
00467 template<class Scalar>
00468 Thyra::ModelEvaluatorBase::InArgs<Scalar> 
00469 ExplicitTaylorPolynomialStepper<Scalar>::getInitialCondition() const
00470 {
00471   return basePoint_;
00472 }
00473 
00474 
00475 template<class Scalar>
00476 Scalar 
00477 ExplicitTaylorPolynomialStepper<Scalar>::takeStep(Scalar dt, StepSizeType flag)
00478 {
00479   typedef Teuchos::ScalarTraits<Scalar> ST;
00480   TEUCHOS_ASSERT( haveInitialCondition_ );
00481   TEUCHOS_ASSERT( !is_null(model_) );
00482   TEUCHOS_ASSERT( !is_null(parameterList_) ); // parameters are nan otherwise
00483 
00484   V_V(outArg(*x_vector_old_),*x_vector_); // x_vector_old = x_vector
00485   V_V(outArg(*x_dot_vector_old_),*x_dot_vector_); // x_dot_vector_old = x_dot_vector
00486 
00487   if (x_poly_ == Teuchos::null) {
00488     x_poly_ = Teuchos::rcp(new Teuchos::Polynomial<Thyra::VectorBase<Scalar> >(0,*x_vector_,degree_));
00489   }
00490 
00491   if (f_poly_ == Teuchos::null) {
00492     f_poly_ = Teuchos::rcp(new Teuchos::Polynomial<Thyra::VectorBase<Scalar> >(0, *f_vector_, degree_));
00493   }
00494   if (flag == STEP_TYPE_VARIABLE) {
00495     // If t_ > t_final_, we're done
00496     if (t_ > t_final_) {
00497       dt_ = ST::zero();
00498       return dt_;
00499     }
00500 
00501     // Compute a local truncated Taylor series solution to system
00502     computeTaylorSeriesSolution_();
00503 
00504     // Estimate log of radius of convergence of Taylor series
00505     Scalar rho = estimateLogRadius_();
00506 
00507     // Set step size
00508     Scalar shadowed_dt = std::exp(linc_ - rho);
00509 
00510     // If step size is too big, reduce
00511     if (shadowed_dt > max_step_size_) {
00512       shadowed_dt = max_step_size_;
00513     }
00514 
00515     // If step goes past t_final_, reduce
00516     if (t_+shadowed_dt > t_final_) {
00517       shadowed_dt = t_final_-t_;
00518     }
00519 
00520     ScalarMag local_error;
00521 
00522     do {
00523 
00524       // compute x(t_+shadowed_dt), xdot(t_+shadowed_dt)
00525       x_poly_->evaluate(shadowed_dt, x_vector_.get(), x_dot_vector_.get());
00526 
00527       // compute f( x(t_+shadowed_dt), t_+shadowed_dt )
00528       eval_model_explicit<Scalar>(*model_,basePoint_,*x_vector_,t_+shadowed_dt,Teuchos::outArg(*f_vector_));
00529 
00530       // compute || xdot(t_+shadowed_dt) - f( x(t_+shadowed_dt), t_+shadowed_dt ) ||
00531       Thyra::Vp_StV(x_dot_vector_.ptr(), -ST::one(),
00532         *f_vector_);
00533       local_error = norm_inf(*x_dot_vector_);
00534 
00535       if (local_error > local_error_tolerance_) {
00536         shadowed_dt *= 0.7;
00537       }
00538 
00539     } while (local_error > local_error_tolerance_ && shadowed_dt > min_step_size_);
00540 
00541     // Check if minimum step size was reached
00542     TEUCHOS_TEST_FOR_EXCEPTION(shadowed_dt < min_step_size_, 
00543       std::runtime_error,
00544       "ExplicitTaylorPolynomialStepper<Scalar>::takeStep(): " 
00545       << "Step size reached minimum step size " 
00546       << min_step_size_ << ".  Failing step." );
00547 
00548     // Increment t_
00549     t_ += shadowed_dt;
00550 
00551     numSteps_++;
00552 
00553     dt_ = shadowed_dt;
00554 
00555     return shadowed_dt;
00556 
00557   } else {
00558 
00559     // If t_ > t_final_, we're done
00560     if (t_ > t_final_) {
00561       dt_ = Teuchos::ScalarTraits<Scalar>::zero();
00562       return dt_;
00563     }
00564 
00565     // Compute a local truncated Taylor series solution to system
00566     computeTaylorSeriesSolution_();
00567 
00568     // If step size is too big, reduce
00569     if (dt > max_step_size_) {
00570       dt = max_step_size_;
00571     }
00572 
00573     // If step goes past t_final_, reduce
00574     if (t_+dt > t_final_) {
00575       dt = t_final_-t_;
00576     }
00577 
00578     // compute x(t_+dt)
00579     x_poly_->evaluate(dt, x_vector_.get());
00580 
00581     // Increment t_
00582     t_ += dt;
00583 
00584     numSteps_++;
00585 
00586     dt_ = dt;
00587 
00588     return dt;
00589   }
00590 }
00591 
00592 
00593 template<class Scalar>
00594 const StepStatus<Scalar>
00595 ExplicitTaylorPolynomialStepper<Scalar>::getStepStatus() const
00596 {
00597   typedef Teuchos::ScalarTraits<Scalar> ST;
00598   StepStatus<Scalar> stepStatus;
00599 
00600   if (!haveInitialCondition_) {
00601     stepStatus.stepStatus = STEP_STATUS_UNINITIALIZED;
00602   } 
00603   else if (numSteps_ == 0) {
00604     stepStatus.stepStatus = STEP_STATUS_UNKNOWN;
00605     stepStatus.stepSize = dt_;
00606     stepStatus.order = this->getOrder();
00607     stepStatus.time = t_;
00608     stepStatus.solution = x_vector_;
00609     stepStatus.solutionDot = x_dot_vector_;
00610     if (!is_null(model_)) {
00611       stepStatus.residual = f_vector_;
00612     }
00613   } 
00614   else  {
00615     stepStatus.stepStatus = STEP_STATUS_CONVERGED;
00616     stepStatus.stepSize = dt_;
00617     stepStatus.order = this->getOrder();
00618     stepStatus.time = t_;
00619     stepStatus.solution = x_vector_;
00620     stepStatus.solutionDot = x_dot_vector_;
00621     stepStatus.residual = f_vector_;
00622   }
00623   return(stepStatus);
00624 }
00625 
00626 
00627 template<class Scalar>
00628 void ExplicitTaylorPolynomialStepper<Scalar>::setParameterList(RCP<Teuchos::ParameterList> const& paramList)
00629 {
00630   typedef Teuchos::ScalarTraits<Scalar> ST;
00631 
00632   TEUCHOS_TEST_FOR_EXCEPT(is_null(paramList));
00633   paramList->validateParameters(*this->getValidParameters());
00634   parameterList_ = paramList;
00635   Teuchos::readVerboseObjectSublist(&*parameterList_,this);
00636 
00637   // Get initial time
00638   t_initial_ = parameterList_->get("Initial Time", ST::zero());
00639 
00640   // Get final time
00641   t_final_ = parameterList_->get("Final Time", ST::one());
00642 
00643   // Get local error tolerance
00644   local_error_tolerance_ = 
00645     parameterList_->get("Local Error Tolerance", ScalarMag(1.0e-10));
00646 
00647   // Get minimum step size
00648   min_step_size_ = parameterList_->get("Minimum Step Size", Scalar(1.0e-10));
00649 
00650   // Get maximum step size
00651   max_step_size_ = parameterList_->get("Maximum Step Size", Scalar(1.0));
00652 
00653   // Get degree_ of Taylor polynomial expansion
00654   degree_ = parameterList_->get("Taylor Polynomial Degree", Teuchos::as<unsigned int>(40));
00655 
00656   linc_ = Scalar(-16.0*std::log(10.0)/degree_);
00657   t_ = t_initial_;
00658 }
00659 
00660 
00661 template<class Scalar>
00662 RCP<Teuchos::ParameterList>
00663 ExplicitTaylorPolynomialStepper<Scalar>::getNonconstParameterList()
00664 {
00665   return parameterList_;
00666 }
00667 
00668 
00669 template<class Scalar>
00670 RCP<Teuchos::ParameterList>
00671 ExplicitTaylorPolynomialStepper<Scalar>:: unsetParameterList()
00672 {
00673   RCP<Teuchos::ParameterList> temp_param_list = parameterList_;
00674   parameterList_ = Teuchos::null;
00675   return temp_param_list;
00676 }
00677 
00678 
00679 template<class Scalar>
00680 RCP<const Teuchos::ParameterList>
00681 ExplicitTaylorPolynomialStepper<Scalar>::getValidParameters() const
00682 {
00683   typedef ScalarTraits<Scalar> ST;
00684   static RCP<const ParameterList> validPL;
00685   if (is_null(validPL)) {
00686     RCP<ParameterList> pl = Teuchos::parameterList();
00687 
00688     pl->set<Scalar>("Initial Time", ST::zero());
00689     pl->set<Scalar>("Final Time", ST::one());
00690     pl->set<ScalarMag>("Local Error Tolerance", ScalarMag(1.0e-10));
00691     pl->set<Scalar>("Minimum Step Size", Scalar(1.0e-10));
00692     pl->set<Scalar>("Maximum Step Size", Scalar(1.0));
00693     pl->set<unsigned int>("Taylor Polynomial Degree", 40);
00694 
00695     Teuchos::setupVerboseObjectSublist(&*pl);
00696     validPL = pl;
00697   }
00698   return validPL;
00699 }
00700 
00701 
00702 template<class Scalar>
00703 std::string ExplicitTaylorPolynomialStepper<Scalar>::description() const
00704 {
00705   std::string name = "Rythmos::ExplicitTaylorPolynomialStepper";
00706   return name;
00707 }
00708 
00709 
00710 template<class Scalar>
00711 void ExplicitTaylorPolynomialStepper<Scalar>::describe(
00712   Teuchos::FancyOStream &out,
00713   const Teuchos::EVerbosityLevel verbLevel
00714   ) const
00715 {
00716   if (verbLevel == Teuchos::VERB_EXTREME) {
00717     out << description() << "::describe" << std::endl;
00718     out << "model_ = " << std::endl;
00719     out << Teuchos::describe(*model_, verbLevel) << std::endl;
00720     out << "x_vector_ = " << std::endl;
00721     out << Teuchos::describe(*x_vector_, verbLevel) << std::endl;
00722     out << "x_dot_vector_ = " << std::endl;
00723     out << Teuchos::describe(*x_dot_vector_, verbLevel) << std::endl;
00724     out << "f_vector_ = " << std::endl;
00725     out << Teuchos::describe(*f_vector_, verbLevel) << std::endl;
00726     out << "x_poly_ = " << std::endl;
00727     out << Teuchos::describe(*x_poly_, verbLevel) << std::endl;
00728     out << "f_poly_ = " << std::endl;
00729     out << Teuchos::describe(*f_poly_, verbLevel) << std::endl;
00730     out << "t_ = " << t_ << std::endl;
00731     out << "t_initial_ = " << t_initial_ << std::endl;
00732     out << "t_final_ = " << t_final_ << std::endl;
00733     out << "local_error_tolerance_ = " << local_error_tolerance_ << std::endl;
00734     out << "min_step_size_ = " << min_step_size_ << std::endl;
00735     out << "max_step_size_ = " << max_step_size_ << std::endl;
00736     out << "degree_ = " << degree_ << std::endl;
00737     out << "linc_ = " << linc_ << std::endl;
00738   }
00739 }
00740 
00741 
00742 template<class Scalar>
00743 void ExplicitTaylorPolynomialStepper<Scalar>::addPoints(
00744   const Array<Scalar>& time_vec
00745   ,const Array<RCP<const Thyra::VectorBase<Scalar> > >& x_vec
00746   ,const Array<RCP<const Thyra::VectorBase<Scalar> > >& xdot_vec
00747   )
00748 {
00749   TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Error, addPoints is not implemented for the ExplicitTaylorPolynomialStepper.\n");
00750 }
00751 
00752 
00753 template<class Scalar>
00754 void ExplicitTaylorPolynomialStepper<Scalar>::getPoints(
00755   const Array<Scalar>& time_vec
00756   ,Array<RCP<const Thyra::VectorBase<Scalar> > >* x_vec
00757   ,Array<RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec
00758   ,Array<ScalarMag>* accuracy_vec) const
00759 {
00760   TEUCHOS_ASSERT( haveInitialCondition_ );
00761   using Teuchos::constOptInArg;
00762   using Teuchos::null;
00763   defaultGetPoints<Scalar>(
00764     t_-dt_,constOptInArg(*x_vector_old_),constOptInArg(*x_dot_vector_old_),
00765     t_,constOptInArg(*x_vector_),constOptInArg(*x_dot_vector_),
00766     time_vec,ptr(x_vec),ptr(xdot_vec),ptr(accuracy_vec),
00767     Ptr<InterpolatorBase<Scalar> >(null)
00768     );
00769 }
00770 
00771 
00772 template<class Scalar>
00773 TimeRange<Scalar> ExplicitTaylorPolynomialStepper<Scalar>::getTimeRange() const
00774 {
00775   if (!haveInitialCondition_) {
00776     return invalidTimeRange<Scalar>();
00777   } else {
00778     return(TimeRange<Scalar>(t_-dt_,t_));
00779   }
00780 }
00781 
00782 
00783 template<class Scalar>
00784 void ExplicitTaylorPolynomialStepper<Scalar>::getNodes(Array<Scalar>* time_vec) const
00785 {
00786   TEUCHOS_ASSERT( time_vec != NULL );
00787   time_vec->clear();
00788   if (!haveInitialCondition_) {
00789     return; 
00790   } else {
00791     time_vec->push_back(t_);
00792   }
00793   if (numSteps_ > 0) {
00794     time_vec->push_back(t_-dt_);
00795   }
00796 }
00797 
00798 
00799 template<class Scalar>
00800 void ExplicitTaylorPolynomialStepper<Scalar>::removeNodes(Array<Scalar>& time_vec)
00801 {
00802   TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Error, removeNodes is not implemented for the ExplicitTaylorPolynomialStepper.\n");
00803 }
00804 
00805 
00806 template<class Scalar>
00807 int ExplicitTaylorPolynomialStepper<Scalar>::getOrder() const
00808 {
00809   return degree_;
00810 }
00811 
00812 
00813 //
00814 // Definitions of protected methods
00815 //
00816 
00817 
00818 template<class Scalar>
00819 void
00820 ExplicitTaylorPolynomialStepper<Scalar>::computeTaylorSeriesSolution_()
00821 {
00822   RCP<Thyra::VectorBase<Scalar> > tmp;
00823 
00824   // Set degree_ of polynomials to 0
00825   x_poly_->setDegree(0);
00826   f_poly_->setDegree(0);
00827 
00828   // Set degree_ 0 coefficient
00829   x_poly_->setCoefficient(0, *x_vector_);
00830 
00831   for (unsigned int k=1; k<=degree_; k++) {
00832 
00833     // compute [f] = f([x])
00834     eval_model_explicit_poly(*model_, basePoint_, *x_poly_, t_, Teuchos::outArg(*f_poly_));
00835 
00836     x_poly_->setDegree(k);
00837     f_poly_->setDegree(k);
00838       
00839     // x[k] = f[k-1] / k
00840     tmp = x_poly_->getCoefficient(k);
00841     copy(*(f_poly_->getCoefficient(k-1)), tmp.ptr());
00842     scale(Scalar(1.0)/Scalar(k), tmp.ptr());
00843   }
00844 
00845 }
00846 
00847 
00848 template<class Scalar>
00849 typename ExplicitTaylorPolynomialStepper<Scalar>::ScalarMag
00850 ExplicitTaylorPolynomialStepper<Scalar>::estimateLogRadius_()
00851 {
00852   ScalarMag rho = 0;
00853   ScalarMag tmp;
00854   for (unsigned int k=degree_/2; k<=degree_; k++) {
00855     tmp = log_norm_inf(*(x_poly_->getCoefficient(k))) / k;
00856     if (tmp > rho) {
00857       rho = tmp;
00858     }
00859   }
00860   return rho;
00861 }
00862 
00863 
00864 template<class Scalar>
00865 RCP<const Thyra::VectorSpaceBase<Scalar> > ExplicitTaylorPolynomialStepper<Scalar>::get_x_space() const
00866 {
00867   if (haveInitialCondition_) {
00868     return(x_vector_->space());
00869   } else {
00870     return Teuchos::null;
00871   }
00872 }
00873 
00874 
00875 } // namespace Rythmos
00876 
00877 #endif // RYTHMOS_EXPLICIT_TAYLOR_POLYNOMIAL_STEPPER_H
 All Classes Functions Variables Typedefs Friends