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

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