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., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
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   Thyra::applyOp<Scalar>(log_norm_inf_op,
00348     Teuchos::tuple(Teuchos::ptrFromRef(x))(), Teuchos::null,
00349     log_norm_inf_targ.ptr());
00350   return log_norm_inf_op(*log_norm_inf_targ);
00351 }
00352 
00353 
00354 // Non-member constructor
00355 template<class Scalar>
00356 RCP<ExplicitTaylorPolynomialStepper<Scalar> > explicitTaylorPolynomialStepper()
00357 {
00358   RCP<ExplicitTaylorPolynomialStepper<Scalar> > stepper = rcp(new ExplicitTaylorPolynomialStepper<Scalar>());
00359   return stepper;
00360 }
00361 
00362 
00363 template<class Scalar>
00364 ExplicitTaylorPolynomialStepper<Scalar>::ExplicitTaylorPolynomialStepper()
00365 {
00366   this->defaultInitializAll_();
00367   numSteps_ = 0;
00368 }
00369 
00370 
00371 template<class Scalar>
00372 ExplicitTaylorPolynomialStepper<Scalar>::~ExplicitTaylorPolynomialStepper()
00373 {
00374 }
00375 
00376 
00377 template<class Scalar>
00378 void ExplicitTaylorPolynomialStepper<Scalar>::defaultInitializAll_()
00379 {
00380   typedef Teuchos::ScalarTraits<Scalar> ST;
00381   Scalar nan = ST::nan();
00382   model_ = Teuchos::null;
00383   parameterList_ = Teuchos::null;
00384   x_vector_ = Teuchos::null;
00385   x_vector_old_ = Teuchos::null;
00386   x_dot_vector_ = Teuchos::null;
00387   x_dot_vector_old_ = Teuchos::null;
00388   f_vector_ = Teuchos::null;
00389   x_poly_ = Teuchos::null;
00390   f_poly_ = Teuchos::null;
00391   haveInitialCondition_ = false;
00392   numSteps_ = -1;
00393   t_ = nan;
00394   dt_ = nan;
00395   t_initial_ = nan;
00396   t_final_ = nan;
00397   local_error_tolerance_ = nan;
00398   min_step_size_ = nan;
00399   max_step_size_ = nan;
00400   degree_ = 0;
00401   linc_ = nan;
00402 }
00403 
00404 
00405 template<class Scalar>
00406 void ExplicitTaylorPolynomialStepper<Scalar>::setModel(
00407   const RCP<const Thyra::ModelEvaluator<Scalar> >& model
00408   )
00409 {
00410   TEUCHOS_TEST_FOR_EXCEPT( is_null(model) );
00411   assertValidModel( *this, *model );
00412     
00413   model_ = model;
00414   f_vector_ = Thyra::createMember(model_->get_f_space());
00415 }
00416 
00417 
00418 template<class Scalar>
00419 void ExplicitTaylorPolynomialStepper<Scalar>::setNonconstModel(
00420   const RCP<Thyra::ModelEvaluator<Scalar> >& model
00421   )
00422 {
00423   this->setModel(model); // TODO 09/09/09 tscoffe:  use ConstNonconstObjectContainer!
00424 }
00425 
00426 
00427 template<class Scalar>
00428 RCP<const Thyra::ModelEvaluator<Scalar> >
00429 ExplicitTaylorPolynomialStepper<Scalar>::getModel() const
00430 {
00431   return model_;
00432 }
00433 
00434 
00435 template<class Scalar>
00436 RCP<Thyra::ModelEvaluator<Scalar> >
00437 ExplicitTaylorPolynomialStepper<Scalar>::getNonconstModel() 
00438 {
00439   return Teuchos::null;
00440 }
00441 
00442 
00443 template<class Scalar>
00444 void ExplicitTaylorPolynomialStepper<Scalar>::setInitialCondition(
00445   const Thyra::ModelEvaluatorBase::InArgs<Scalar> &initialCondition
00446   )
00447 {
00448   typedef Teuchos::ScalarTraits<Scalar> ST;
00449   typedef Thyra::ModelEvaluatorBase MEB;
00450   basePoint_ = initialCondition;
00451   if (initialCondition.supports(MEB::IN_ARG_t)) {
00452     t_ = initialCondition.get_t();
00453   } else {
00454     t_ = ST::zero();
00455   }
00456   dt_ = ST::zero();
00457   x_vector_ = initialCondition.get_x()->clone_v();
00458   x_dot_vector_ = x_vector_->clone_v();
00459   x_vector_old_ = x_vector_->clone_v();
00460   x_dot_vector_old_ = x_dot_vector_->clone_v();
00461   haveInitialCondition_ = true;
00462 }
00463 
00464 
00465 template<class Scalar>
00466 Thyra::ModelEvaluatorBase::InArgs<Scalar> 
00467 ExplicitTaylorPolynomialStepper<Scalar>::getInitialCondition() const
00468 {
00469   return basePoint_;
00470 }
00471 
00472 
00473 template<class Scalar>
00474 Scalar 
00475 ExplicitTaylorPolynomialStepper<Scalar>::takeStep(Scalar dt, StepSizeType flag)
00476 {
00477   typedef Teuchos::ScalarTraits<Scalar> ST;
00478   TEUCHOS_ASSERT( haveInitialCondition_ );
00479   TEUCHOS_ASSERT( !is_null(model_) );
00480   TEUCHOS_ASSERT( !is_null(parameterList_) ); // parameters are nan otherwise
00481 
00482   V_V(outArg(*x_vector_old_),*x_vector_); // x_vector_old = x_vector
00483   V_V(outArg(*x_dot_vector_old_),*x_dot_vector_); // x_dot_vector_old = x_dot_vector
00484 
00485   if (x_poly_ == Teuchos::null) {
00486     x_poly_ = Teuchos::rcp(new Teuchos::Polynomial<Thyra::VectorBase<Scalar> >(0,*x_vector_,degree_));
00487   }
00488 
00489   if (f_poly_ == Teuchos::null) {
00490     f_poly_ = Teuchos::rcp(new Teuchos::Polynomial<Thyra::VectorBase<Scalar> >(0, *f_vector_, degree_));
00491   }
00492   if (flag == STEP_TYPE_VARIABLE) {
00493     // If t_ > t_final_, we're done
00494     if (t_ > t_final_) {
00495       dt_ = ST::zero();
00496       return dt_;
00497     }
00498 
00499     // Compute a local truncated Taylor series solution to system
00500     computeTaylorSeriesSolution_();
00501 
00502     // Estimate log of radius of convergence of Taylor series
00503     Scalar rho = estimateLogRadius_();
00504 
00505     // Set step size
00506     Scalar shadowed_dt = std::exp(linc_ - rho);
00507 
00508     // If step size is too big, reduce
00509     if (shadowed_dt > max_step_size_) {
00510       shadowed_dt = max_step_size_;
00511     }
00512 
00513     // If step goes past t_final_, reduce
00514     if (t_+shadowed_dt > t_final_) {
00515       shadowed_dt = t_final_-t_;
00516     }
00517 
00518     ScalarMag local_error;
00519 
00520     do {
00521 
00522       // compute x(t_+shadowed_dt), xdot(t_+shadowed_dt)
00523       x_poly_->evaluate(shadowed_dt, x_vector_.get(), x_dot_vector_.get());
00524 
00525       // compute f( x(t_+shadowed_dt), t_+shadowed_dt )
00526       eval_model_explicit<Scalar>(*model_,basePoint_,*x_vector_,t_+shadowed_dt,Teuchos::outArg(*f_vector_));
00527 
00528       // compute || xdot(t_+shadowed_dt) - f( x(t_+shadowed_dt), t_+shadowed_dt ) ||
00529       Thyra::Vp_StV(x_dot_vector_.ptr(), -ST::one(),
00530         *f_vector_);
00531       local_error = norm_inf(*x_dot_vector_);
00532 
00533       if (local_error > local_error_tolerance_) {
00534         shadowed_dt *= 0.7;
00535       }
00536 
00537     } while (local_error > local_error_tolerance_ && shadowed_dt > min_step_size_);
00538 
00539     // Check if minimum step size was reached
00540     TEUCHOS_TEST_FOR_EXCEPTION(shadowed_dt < min_step_size_, 
00541       std::runtime_error,
00542       "ExplicitTaylorPolynomialStepper<Scalar>::takeStep(): " 
00543       << "Step size reached minimum step size " 
00544       << min_step_size_ << ".  Failing step." );
00545 
00546     // Increment t_
00547     t_ += shadowed_dt;
00548 
00549     numSteps_++;
00550 
00551     dt_ = shadowed_dt;
00552 
00553     return shadowed_dt;
00554 
00555   } else {
00556 
00557     // If t_ > t_final_, we're done
00558     if (t_ > t_final_) {
00559       dt_ = Teuchos::ScalarTraits<Scalar>::zero();
00560       return dt_;
00561     }
00562 
00563     // Compute a local truncated Taylor series solution to system
00564     computeTaylorSeriesSolution_();
00565 
00566     // If step size is too big, reduce
00567     if (dt > max_step_size_) {
00568       dt = max_step_size_;
00569     }
00570 
00571     // If step goes past t_final_, reduce
00572     if (t_+dt > t_final_) {
00573       dt = t_final_-t_;
00574     }
00575 
00576     // compute x(t_+dt)
00577     x_poly_->evaluate(dt, x_vector_.get());
00578 
00579     // Increment t_
00580     t_ += dt;
00581 
00582     numSteps_++;
00583 
00584     dt_ = dt;
00585 
00586     return dt;
00587   }
00588 }
00589 
00590 
00591 template<class Scalar>
00592 const StepStatus<Scalar>
00593 ExplicitTaylorPolynomialStepper<Scalar>::getStepStatus() const
00594 {
00595   typedef Teuchos::ScalarTraits<Scalar> ST;
00596   StepStatus<Scalar> stepStatus;
00597 
00598   if (!haveInitialCondition_) {
00599     stepStatus.stepStatus = STEP_STATUS_UNINITIALIZED;
00600   } 
00601   else if (numSteps_ == 0) {
00602     stepStatus.stepStatus = STEP_STATUS_UNKNOWN;
00603     stepStatus.stepSize = dt_;
00604     stepStatus.order = this->getOrder();
00605     stepStatus.time = t_;
00606     stepStatus.solution = x_vector_;
00607     stepStatus.solutionDot = x_dot_vector_;
00608     if (!is_null(model_)) {
00609       stepStatus.residual = f_vector_;
00610     }
00611   } 
00612   else  {
00613     stepStatus.stepStatus = STEP_STATUS_CONVERGED;
00614     stepStatus.stepSize = dt_;
00615     stepStatus.order = this->getOrder();
00616     stepStatus.time = t_;
00617     stepStatus.solution = x_vector_;
00618     stepStatus.solutionDot = x_dot_vector_;
00619     stepStatus.residual = f_vector_;
00620   }
00621   return(stepStatus);
00622 }
00623 
00624 
00625 template<class Scalar>
00626 void ExplicitTaylorPolynomialStepper<Scalar>::setParameterList(RCP<Teuchos::ParameterList> const& paramList)
00627 {
00628   typedef Teuchos::ScalarTraits<Scalar> ST;
00629 
00630   TEUCHOS_TEST_FOR_EXCEPT(is_null(paramList));
00631   paramList->validateParameters(*this->getValidParameters());
00632   parameterList_ = paramList;
00633   Teuchos::readVerboseObjectSublist(&*parameterList_,this);
00634 
00635   // Get initial time
00636   t_initial_ = parameterList_->get("Initial Time", ST::zero());
00637 
00638   // Get final time
00639   t_final_ = parameterList_->get("Final Time", ST::one());
00640 
00641   // Get local error tolerance
00642   local_error_tolerance_ = 
00643     parameterList_->get("Local Error Tolerance", ScalarMag(1.0e-10));
00644 
00645   // Get minimum step size
00646   min_step_size_ = parameterList_->get("Minimum Step Size", Scalar(1.0e-10));
00647 
00648   // Get maximum step size
00649   max_step_size_ = parameterList_->get("Maximum Step Size", Scalar(1.0));
00650 
00651   // Get degree_ of Taylor polynomial expansion
00652   degree_ = parameterList_->get("Taylor Polynomial Degree", Teuchos::as<unsigned int>(40));
00653 
00654   linc_ = Scalar(-16.0*std::log(10.0)/degree_);
00655   t_ = t_initial_;
00656 }
00657 
00658 
00659 template<class Scalar>
00660 RCP<Teuchos::ParameterList>
00661 ExplicitTaylorPolynomialStepper<Scalar>::getNonconstParameterList()
00662 {
00663   return parameterList_;
00664 }
00665 
00666 
00667 template<class Scalar>
00668 RCP<Teuchos::ParameterList>
00669 ExplicitTaylorPolynomialStepper<Scalar>:: unsetParameterList()
00670 {
00671   RCP<Teuchos::ParameterList> temp_param_list = parameterList_;
00672   parameterList_ = Teuchos::null;
00673   return temp_param_list;
00674 }
00675 
00676 
00677 template<class Scalar>
00678 RCP<const Teuchos::ParameterList>
00679 ExplicitTaylorPolynomialStepper<Scalar>::getValidParameters() const
00680 {
00681   typedef ScalarTraits<Scalar> ST;
00682   static RCP<const ParameterList> validPL;
00683   if (is_null(validPL)) {
00684     RCP<ParameterList> pl = Teuchos::parameterList();
00685 
00686     pl->set<Scalar>("Initial Time", ST::zero());
00687     pl->set<Scalar>("Final Time", ST::one());
00688     pl->set<ScalarMag>("Local Error Tolerance", ScalarMag(1.0e-10));
00689     pl->set<Scalar>("Minimum Step Size", Scalar(1.0e-10));
00690     pl->set<Scalar>("Maximum Step Size", Scalar(1.0));
00691     pl->set<unsigned int>("Taylor Polynomial Degree", 40);
00692 
00693     Teuchos::setupVerboseObjectSublist(&*pl);
00694     validPL = pl;
00695   }
00696   return validPL;
00697 }
00698 
00699 
00700 template<class Scalar>
00701 std::string ExplicitTaylorPolynomialStepper<Scalar>::description() const
00702 {
00703   std::string name = "Rythmos::ExplicitTaylorPolynomialStepper";
00704   return name;
00705 }
00706 
00707 
00708 template<class Scalar>
00709 void ExplicitTaylorPolynomialStepper<Scalar>::describe(
00710   Teuchos::FancyOStream &out,
00711   const Teuchos::EVerbosityLevel verbLevel
00712   ) const
00713 {
00714   if (verbLevel == Teuchos::VERB_EXTREME) {
00715     out << description() << "::describe" << std::endl;
00716     out << "model_ = " << std::endl;
00717     out << Teuchos::describe(*model_, verbLevel) << std::endl;
00718     out << "x_vector_ = " << std::endl;
00719     out << Teuchos::describe(*x_vector_, verbLevel) << std::endl;
00720     out << "x_dot_vector_ = " << std::endl;
00721     out << Teuchos::describe(*x_dot_vector_, verbLevel) << std::endl;
00722     out << "f_vector_ = " << std::endl;
00723     out << Teuchos::describe(*f_vector_, verbLevel) << std::endl;
00724     out << "x_poly_ = " << std::endl;
00725     out << Teuchos::describe(*x_poly_, verbLevel) << std::endl;
00726     out << "f_poly_ = " << std::endl;
00727     out << Teuchos::describe(*f_poly_, verbLevel) << std::endl;
00728     out << "t_ = " << t_ << std::endl;
00729     out << "t_initial_ = " << t_initial_ << std::endl;
00730     out << "t_final_ = " << t_final_ << std::endl;
00731     out << "local_error_tolerance_ = " << local_error_tolerance_ << std::endl;
00732     out << "min_step_size_ = " << min_step_size_ << std::endl;
00733     out << "max_step_size_ = " << max_step_size_ << std::endl;
00734     out << "degree_ = " << degree_ << std::endl;
00735     out << "linc_ = " << linc_ << std::endl;
00736   }
00737 }
00738 
00739 
00740 template<class Scalar>
00741 void ExplicitTaylorPolynomialStepper<Scalar>::addPoints(
00742   const Array<Scalar>& time_vec
00743   ,const Array<RCP<const Thyra::VectorBase<Scalar> > >& x_vec
00744   ,const Array<RCP<const Thyra::VectorBase<Scalar> > >& xdot_vec
00745   )
00746 {
00747   TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Error, addPoints is not implemented for the ExplicitTaylorPolynomialStepper.\n");
00748 }
00749 
00750 
00751 template<class Scalar>
00752 void ExplicitTaylorPolynomialStepper<Scalar>::getPoints(
00753   const Array<Scalar>& time_vec
00754   ,Array<RCP<const Thyra::VectorBase<Scalar> > >* x_vec
00755   ,Array<RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec
00756   ,Array<ScalarMag>* accuracy_vec) const
00757 {
00758   TEUCHOS_ASSERT( haveInitialCondition_ );
00759   using Teuchos::constOptInArg;
00760   using Teuchos::null;
00761   defaultGetPoints<Scalar>(
00762     t_-dt_,constOptInArg(*x_vector_old_),constOptInArg(*x_dot_vector_old_),
00763     t_,constOptInArg(*x_vector_),constOptInArg(*x_dot_vector_),
00764     time_vec,ptr(x_vec),ptr(xdot_vec),ptr(accuracy_vec),
00765     Ptr<InterpolatorBase<Scalar> >(null)
00766     );
00767 }
00768 
00769 
00770 template<class Scalar>
00771 TimeRange<Scalar> ExplicitTaylorPolynomialStepper<Scalar>::getTimeRange() const
00772 {
00773   if (!haveInitialCondition_) {
00774     return invalidTimeRange<Scalar>();
00775   } else {
00776     return(TimeRange<Scalar>(t_-dt_,t_));
00777   }
00778 }
00779 
00780 
00781 template<class Scalar>
00782 void ExplicitTaylorPolynomialStepper<Scalar>::getNodes(Array<Scalar>* time_vec) const
00783 {
00784   TEUCHOS_ASSERT( time_vec != NULL );
00785   time_vec->clear();
00786   if (!haveInitialCondition_) {
00787     return; 
00788   } else {
00789     time_vec->push_back(t_);
00790   }
00791   if (numSteps_ > 0) {
00792     time_vec->push_back(t_-dt_);
00793   }
00794 }
00795 
00796 
00797 template<class Scalar>
00798 void ExplicitTaylorPolynomialStepper<Scalar>::removeNodes(Array<Scalar>& time_vec)
00799 {
00800   TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Error, removeNodes is not implemented for the ExplicitTaylorPolynomialStepper.\n");
00801 }
00802 
00803 
00804 template<class Scalar>
00805 int ExplicitTaylorPolynomialStepper<Scalar>::getOrder() const
00806 {
00807   return degree_;
00808 }
00809 
00810 
00811 //
00812 // Definitions of protected methods
00813 //
00814 
00815 
00816 template<class Scalar>
00817 void
00818 ExplicitTaylorPolynomialStepper<Scalar>::computeTaylorSeriesSolution_()
00819 {
00820   RCP<Thyra::VectorBase<Scalar> > tmp;
00821 
00822   // Set degree_ of polynomials to 0
00823   x_poly_->setDegree(0);
00824   f_poly_->setDegree(0);
00825 
00826   // Set degree_ 0 coefficient
00827   x_poly_->setCoefficient(0, *x_vector_);
00828 
00829   for (unsigned int k=1; k<=degree_; k++) {
00830 
00831     // compute [f] = f([x])
00832     eval_model_explicit_poly(*model_, basePoint_, *x_poly_, t_, Teuchos::outArg(*f_poly_));
00833 
00834     x_poly_->setDegree(k);
00835     f_poly_->setDegree(k);
00836       
00837     // x[k] = f[k-1] / k
00838     tmp = x_poly_->getCoefficient(k);
00839     copy(*(f_poly_->getCoefficient(k-1)), tmp.ptr());
00840     scale(Scalar(1.0)/Scalar(k), tmp.ptr());
00841   }
00842 
00843 }
00844 
00845 
00846 template<class Scalar>
00847 typename ExplicitTaylorPolynomialStepper<Scalar>::ScalarMag
00848 ExplicitTaylorPolynomialStepper<Scalar>::estimateLogRadius_()
00849 {
00850   ScalarMag rho = 0;
00851   ScalarMag tmp;
00852   for (unsigned int k=degree_/2; k<=degree_; k++) {
00853     tmp = log_norm_inf(*(x_poly_->getCoefficient(k))) / k;
00854     if (tmp > rho) {
00855       rho = tmp;
00856     }
00857   }
00858   return rho;
00859 }
00860 
00861 
00862 template<class Scalar>
00863 RCP<const Thyra::VectorSpaceBase<Scalar> > ExplicitTaylorPolynomialStepper<Scalar>::get_x_space() const
00864 {
00865   if (haveInitialCondition_) {
00866     return(x_vector_->space());
00867   } else {
00868     return Teuchos::null;
00869   }
00870 }
00871 
00872 
00873 } // namespace Rythmos
00874 
00875 #endif // RYTHMOS_EXPLICIT_TAYLOR_POLYNOMIAL_STEPPER_H
 All Classes Functions Variables Typedefs Friends