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 "Teuchos_RCP.hpp"
00034 #include "Teuchos_ParameterList.hpp"
00035 #include "Thyra_VectorBase.hpp"
00036 #include "Thyra_ModelEvaluator.hpp"
00037 #include "Thyra_ModelEvaluatorHelpers.hpp"
00038 #include "RTOpPack_RTOpTHelpers.hpp"
00039 
00040 namespace Rythmos {
00041 
00139   template<class Scalar>
00140   class ExplicitTaylorPolynomialStepper : virtual public StepperBase<Scalar>
00141   {
00142     public:
00143 
00145     typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType ScalarMag;
00146     
00148     ExplicitTaylorPolynomialStepper();
00149     
00151     ~ExplicitTaylorPolynomialStepper();
00152 
00154     RCP<const Thyra::VectorSpaceBase<Scalar> > get_x_space() const;
00155 
00157     void setModel(const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > &model);
00158 
00160     Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >
00161     getModel() const;
00162     
00164     Scalar takeStep(Scalar dt, StepSizeType flag);
00165 
00167     const StepStatus<Scalar> getStepStatus() const;
00168 
00170 
00171     void setParameterList(Teuchos::RCP<Teuchos::ParameterList> const& paramList);
00172 
00174     Teuchos::RCP<Teuchos::ParameterList> getParameterList();
00175 
00177     Teuchos::RCP<Teuchos::ParameterList> unsetParameterList();
00178 
00180     std::string description() const;
00181 
00183     std::ostream& describe(
00184       std::ostream                &out
00185       ,const Teuchos::EVerbosityLevel      verbLevel
00186       ,const std::string          leadingIndent
00187       ,const std::string          indentSpacer
00188       ) const;
00189 
00192     void addPoints(
00193       const Array<Scalar>& time_vec
00194       ,const Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >& x_vec
00195       ,const Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >& xdot_vec
00196       );
00197     
00199     void getPoints(
00200       const Array<Scalar>& time_vec
00201       ,Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >* x_vec
00202       ,Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec
00203       ,Array<ScalarMag>* accuracy_vec) const;
00204 
00206     void setRange(
00207       const TimeRange<Scalar>& range,
00208       const InterpolationBufferBase<Scalar> & IB
00209       );
00210 
00212     TimeRange<Scalar> getTimeRange() const;
00213 
00215     void getNodes(Array<Scalar>* time_vec) const;
00216 
00218     void removeNodes(Array<Scalar>& time_vec);
00219 
00221     int getOrder() const;
00222 
00223     private:
00224 
00226     void computeTaylorSeriesSolution_();
00227 
00232     ScalarMag estimateLogRadius_();
00233 
00235     Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > model_;
00236 
00238     Teuchos::RCP<Teuchos::ParameterList> parameterList_;
00239 
00241     Teuchos::RCP<Thyra::VectorBase<Scalar> > x_vector_;
00242 
00244     Teuchos::RCP<Thyra::VectorBase<Scalar> > x_dot_vector_;
00245 
00247     Teuchos::RCP<Thyra::VectorBase<Scalar> > f_vector_;
00248 
00250     Teuchos::RCP<Teuchos::Polynomial<Thyra::VectorBase<Scalar> > > x_poly_;
00251 
00253     Teuchos::RCP<Teuchos::Polynomial<Thyra::VectorBase<Scalar> > > f_poly_;
00254 
00256     Scalar t_;
00257 
00259     Scalar dt_;
00260 
00262     Scalar t_initial_;
00263 
00265     Scalar t_final_;
00266 
00268     ScalarMag local_error_tolerance_;
00269 
00271     Scalar min_step_size_;
00272 
00274     Scalar max_step_size_;
00275 
00277     unsigned int degree_;
00278 
00280     Scalar linc_;
00281   };
00282 
00284 
00291   template <typename Scalar>
00292     class ROpLogNormInf : public RTOpPack::ROpScalarReductionBase<typename Teuchos::ScalarTraits<Scalar>::magnitudeType> {
00293       public:
00294 
00296         ROpLogNormInf() : RTOpPack::RTOpT<Scalar>("ROpLogInfNorm"){}
00297 
00299         typename Teuchos::ScalarTraits<Scalar>::magnitudeType 
00300           operator() (const RTOpPack::ReductTarget& reduct_obj) const { 
00301             return this->getRawVal(reduct_obj); 
00302           }
00303 
00306 
00308         void reduce_reduct_objs(const RTOpPack::ReductTarget& _in_reduct_obj, 
00309             RTOpPack::ReductTarget* _inout_reduct_obj) const
00310         {
00311           const RTOpPack::ReductTargetScalar<Scalar>& in_reduct_obj = 
00312             Teuchos::dyn_cast<const RTOpPack::ReductTargetScalar<Scalar> >(_in_reduct_obj);
00313           RTOpPack::ReductTargetScalar<Scalar>& inout_reduct_obj = 
00314             Teuchos::dyn_cast<RTOpPack::ReductTargetScalar<Scalar> >(*_inout_reduct_obj);
00315           if(in_reduct_obj.get() > inout_reduct_obj.get()) 
00316             inout_reduct_obj.set(in_reduct_obj.get());
00317         }
00318 
00320         void apply_op(const int num_vecs, 
00321             const RTOpPack::ConstSubVectorView<Scalar> sub_vecs[],
00322             const int num_targ_vecs, 
00323             const RTOpPack::SubVectorView<Scalar> targ_sub_vecs[], 
00324             RTOpPack::ReductTarget *_reduct_obj) const
00325         {
00326           RTOpPack::ReductTargetScalar<Scalar>& reduct_obj = 
00327             Teuchos::dyn_cast<RTOpPack::ReductTargetScalar<Scalar> >(*_reduct_obj);
00328 
00329           RTOP_APPLY_OP_1_0(num_vecs,sub_vecs,num_targ_vecs,targ_sub_vecs);
00330 
00331           typename Teuchos::ScalarTraits<Scalar>::magnitudeType norm_inf = 
00332             reduct_obj.get();
00333 
00334           // unit stride
00335           if(v0_s == 1) {
00336             for(Teuchos_Index i=0; i<subDim; i++) {
00337               typename Teuchos::ScalarTraits<Scalar>::magnitudeType
00338                 mag = Teuchos::ScalarTraits<Scalar>::magnitude(*v0_val++);
00339               mag = std::log(Teuchos::ScalarTraits<Scalar>::one()*1e-100 + mag);
00340               norm_inf = mag > norm_inf ? mag : norm_inf;
00341             }
00342           } else {
00343             for(Teuchos_Index i=0; i<subDim; i++, v0_val+=v0_s) {
00344               typename Teuchos::ScalarTraits<Scalar>::magnitudeType
00345                 mag = Teuchos::ScalarTraits<Scalar>::magnitude(*v0_val);
00346               mag = std::log(Teuchos::ScalarTraits<Scalar>::one()*1e-100 + mag);
00347               norm_inf = mag > norm_inf ? mag : norm_inf;
00348             }
00349           }
00350           reduct_obj.set(norm_inf);
00351         }
00353 
00354     }; // class ROpLogNormInf
00355 
00357   template <typename Scalar>
00358   typename Teuchos::ScalarTraits<Scalar>::magnitudeType
00359   log_norm_inf(const Thyra::VectorBase<Scalar>& x)
00360   {
00361     ROpLogNormInf<Scalar> log_norm_inf_op;
00362     Teuchos::RCP<RTOpPack::ReductTarget> log_norm_inf_targ = 
00363       log_norm_inf_op.reduct_obj_create();
00364     const Thyra::VectorBase<Scalar>* vecs[] = { &x };
00365     Thyra::applyOp<Scalar>(log_norm_inf_op,1,vecs,0,
00366          (Thyra::VectorBase<Scalar>**)NULL,
00367          log_norm_inf_targ.get());
00368     
00369     return log_norm_inf_op(*log_norm_inf_targ);
00370   }
00371 
00372   template<class Scalar>
00373   ExplicitTaylorPolynomialStepper<Scalar>::ExplicitTaylorPolynomialStepper()
00374   {
00375   }
00376 
00377   template<class Scalar>
00378   ExplicitTaylorPolynomialStepper<Scalar>::~ExplicitTaylorPolynomialStepper()
00379   {
00380   }
00381 
00382   template<class Scalar>
00383   void ExplicitTaylorPolynomialStepper<Scalar>::setModel(
00384     const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > &model
00385     )
00386   {
00387     TEST_FOR_EXCEPT(model == Teuchos::null)
00388     
00389     model_ = model;
00390     x_vector_ = model_->getNominalValues().get_x()->clone_v();
00391     x_dot_vector_ = Thyra::createMember(model_->get_x_space());
00392     f_vector_ = Thyra::createMember(model_->get_f_space());
00393   }
00394 
00395 
00396   template<class Scalar>
00397   Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >
00398   ExplicitTaylorPolynomialStepper<Scalar>::getModel() const
00399   {
00400     return model_;
00401   }
00402 
00403   template<class Scalar>
00404   Scalar 
00405   ExplicitTaylorPolynomialStepper<Scalar>::takeStep(Scalar dt, StepSizeType flag)
00406   {
00407     if (x_poly_ == Teuchos::null)
00408       x_poly_ = 
00409   Teuchos::rcp(new Teuchos::Polynomial<
00410          Thyra::VectorBase<Scalar> >(0,*x_vector_,degree_));
00411 
00412     if (f_poly_ == Teuchos::null)
00413       f_poly_ = 
00414   Teuchos::rcp(new Teuchos::Polynomial<
00415          Thyra::VectorBase<Scalar> >(0, *f_vector_, degree_));
00416     if (flag == STEP_TYPE_VARIABLE) {
00417       // If t_ > t_final_, we're done
00418       if (t_ > t_final_) {
00419         dt_ = Teuchos::ScalarTraits<Scalar>::zero();
00420         return dt_;
00421       }
00422 
00423       // Compute a local truncated Taylor series solution to system
00424       computeTaylorSeriesSolution_();
00425 
00426       // Estimate log of radius of convergence of Taylor series
00427       Scalar rho = estimateLogRadius_();
00428 
00429       // Set step size
00430       Scalar dt = std::exp(linc_ - rho);
00431 
00432       // If step size is too big, reduce
00433       if (dt > max_step_size_) {
00434         dt = max_step_size_;
00435       }
00436 
00437       // If step goes past t_final_, reduce
00438       if (t_+dt > t_final_) {
00439         dt = t_final_-t_;
00440       }
00441 
00442       ScalarMag local_error;
00443 
00444       do {
00445 
00446         // compute x(t_+dt), xdot(t_+dt)
00447         x_poly_->evaluate(dt, x_vector_.get(), x_dot_vector_.get());
00448 
00449         // compute f( x(t_+dt), t_+dt )
00450         Thyra::eval_f(*model_, *x_vector_, t_+dt, f_vector_.get());
00451 
00452         // compute || xdot(t_+dt) - f( x(t_+dt), t_+dt ) ||
00453         Thyra::Vp_StV(x_dot_vector_.get(), -Teuchos::ScalarTraits<Scalar>::one(),
00454           *f_vector_);
00455         local_error = norm_inf(*x_dot_vector_);
00456 
00457         if (local_error > local_error_tolerance_) {
00458           dt *= 0.7;
00459         }
00460 
00461       } while (local_error > local_error_tolerance_ && dt > min_step_size_);
00462 
00463       // Check if minimum step size was reached
00464       TEST_FOR_EXCEPTION(dt < min_step_size_, 
00465             std::runtime_error,
00466             "ExplicitTaylorPolynomialStepper<Scalar>::takeStep(): " 
00467             << "Step size reached minimum step size " 
00468             << min_step_size_ << ".  Failing step." );
00469 
00470       // Increment t_
00471       t_ += dt;
00472 
00473       dt_ = dt;
00474 
00475       return dt;
00476 
00477     } else {
00478 
00479       // If t_ > t_final_, we're done
00480       if (t_ > t_final_) {
00481         dt_ = Teuchos::ScalarTraits<Scalar>::zero();
00482         return dt_;
00483       }
00484 
00485       // Compute a local truncated Taylor series solution to system
00486       computeTaylorSeriesSolution_();
00487 
00488       // If step size is too big, reduce
00489       if (dt > max_step_size_) {
00490         dt = max_step_size_;
00491       }
00492 
00493       // If step goes past t_final_, reduce
00494       if (t_+dt > t_final_) {
00495         dt = t_final_-t_;
00496       }
00497 
00498       // compute x(t_+dt)
00499       x_poly_->evaluate(dt, x_vector_.get());
00500 
00501       // Increment t_
00502       t_ += dt;
00503 
00504       dt_ = dt;
00505 
00506       return dt;
00507     }
00508   }
00509 
00510 
00511   template<class Scalar>
00512   const StepStatus<Scalar>
00513   ExplicitTaylorPolynomialStepper<Scalar>::getStepStatus() const
00514   {
00515     typedef Teuchos::ScalarTraits<Scalar> ST;
00516     StepStatus<Scalar> stepStatus;
00517 
00518     stepStatus.stepSize = dt_;
00519     stepStatus.order = this->getOrder();
00520     stepStatus.time = t_;
00521     stepStatus.solution = x_vector_;
00522     stepStatus.solutionDot = x_dot_vector_;
00523     stepStatus.residual = f_vector_;
00524 
00525     return(stepStatus);
00526   }
00527 
00528   template<class Scalar>
00529   void ExplicitTaylorPolynomialStepper<Scalar>::setParameterList(Teuchos::RCP<Teuchos::ParameterList> const& paramList)
00530   {
00531     typedef Teuchos::ScalarTraits<Scalar> ST;
00532 
00533     parameterList_ = paramList;
00534 
00535     // Get initial time
00536     t_initial_ = parameterList_->get("Initial Time", ST::zero());
00537 
00538     // Get final time
00539     t_final_ = parameterList_->get("Final Time", ST::one());
00540 
00541     // Get local error tolerance
00542     local_error_tolerance_ = 
00543       parameterList_->get("Local Error Tolerance", ScalarMag(1.0e-10));
00544 
00545     // Get minimum step size
00546     min_step_size_ = parameterList_->get("Minimum Step Size", Scalar(1.0e-10));
00547 
00548     // Get maximum step size
00549     max_step_size_ = parameterList_->get("Maximum Step Size", Scalar(1.0));
00550 
00551     // Get degree_ of Taylor polynomial expansion
00552     degree_ = parameterList_->get("Taylor Polynomial Degree", 40);
00553 
00554     linc_ = Scalar(-16.0*std::log(10.0)/degree_);
00555     t_ = t_initial_;
00556   }
00557 
00558   template<class Scalar>
00559   Teuchos::RCP<Teuchos::ParameterList>
00560   ExplicitTaylorPolynomialStepper<Scalar>::getParameterList()
00561   {
00562     return parameterList_;
00563   }
00564 
00565   template<class Scalar>
00566   Teuchos::RCP<Teuchos::ParameterList>
00567   ExplicitTaylorPolynomialStepper<Scalar>:: unsetParameterList()
00568   {
00569     Teuchos::RCP<Teuchos::ParameterList> temp_param_list = parameterList_;
00570     parameterList_ = Teuchos::null;
00571     return temp_param_list;
00572   }
00573 
00574   template<class Scalar>
00575   std::string ExplicitTaylorPolynomialStepper<Scalar>::description() const
00576   {
00577     std::string name = "Rythmos::ExplicitTaylorPolynomialStepper";
00578     return name;
00579   }
00580 
00581   template<class Scalar>
00582   std::ostream& ExplicitTaylorPolynomialStepper<Scalar>::describe(
00583         std::ostream                &out
00584         ,const Teuchos::EVerbosityLevel      verbLevel
00585         ,const std::string          leadingIndent
00586         ,const std::string          indentSpacer
00587         ) const
00588   {
00589     if (verbLevel == Teuchos::VERB_EXTREME) {
00590       out << description() << "::describe" << std::endl;
00591       out << "model_ = " << std::endl;
00592       out << model_->describe(out,verbLevel,leadingIndent,indentSpacer) << std::endl;
00593       out << "x_vector_ = " << std::endl;
00594       out << x_vector_->describe(out,verbLevel,leadingIndent,indentSpacer) << std::endl;
00595       out << "x_dot_vector_ = " << std::endl;
00596       out << x_dot_vector_->describe(out,verbLevel,leadingIndent,indentSpacer) << std::endl;
00597       out << "f_vector_ = " << std::endl;
00598       out << f_vector_->describe(out,verbLevel,leadingIndent,indentSpacer) << std::endl;
00599       out << "x_poly_ = " << std::endl;
00600       out << x_poly_->describe(out,verbLevel,leadingIndent,indentSpacer) << std::endl;
00601       out << "f_poly_ = " << std::endl;
00602       out << f_poly_->describe(out,verbLevel,leadingIndent,indentSpacer) << std::endl;
00603       out << "t_ = " << t_ << std::endl;
00604       out << "t_initial_ = " << t_initial_ << std::endl;
00605       out << "t_final_ = " << t_final_ << std::endl;
00606       out << "local_error_tolerance_ = " << local_error_tolerance_ << std::endl;
00607       out << "min_step_size_ = " << min_step_size_ << std::endl;
00608       out << "max_step_size_ = " << max_step_size_ << std::endl;
00609       out << "degree_ = " << degree_ << std::endl;
00610       out << "linc_ = " << linc_ << std::endl;
00611     }
00612     return(out);
00613   }
00614 
00615 
00616   template<class Scalar>
00617   void ExplicitTaylorPolynomialStepper<Scalar>::addPoints(
00618     const Array<Scalar>& time_vec
00619     ,const Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >& x_vec
00620     ,const Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >& xdot_vec
00621     )
00622   {
00623     TEST_FOR_EXCEPTION(true,std::logic_error,"Error, addPoints is not implemented for the ExplicitTaylorPolynomialStepper.\n");
00624   }
00625 
00626   template<class Scalar>
00627   void ExplicitTaylorPolynomialStepper<Scalar>::getPoints(
00628     const Array<Scalar>& time_vec
00629     ,Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >* x_vec
00630     ,Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec
00631     ,Array<ScalarMag>* accuracy_vec) const
00632   {
00633     TEST_FOR_EXCEPTION(true,std::logic_error,"Error, getPoints is not implemented for the ExplicitTaylorPolynomialStepper.\n");
00634   }
00635 
00636   template<class Scalar>
00637   TimeRange<Scalar> ExplicitTaylorPolynomialStepper<Scalar>::getTimeRange() const
00638   {
00639     return invalidTimeRange<Scalar>();
00640   }
00641 
00642   template<class Scalar>
00643   void ExplicitTaylorPolynomialStepper<Scalar>::getNodes(Array<Scalar>* time_list) const
00644   {
00645     TEST_FOR_EXCEPTION(true,std::logic_error,"Error, getNodes is not implemented for the ExplicitTaylorPolynomialStepper.\n");
00646   }
00647 
00648   template<class Scalar>
00649   void ExplicitTaylorPolynomialStepper<Scalar>::removeNodes(Array<Scalar>& time_vec)
00650   {
00651     TEST_FOR_EXCEPTION(true,std::logic_error,"Error, removeNodes is not implemented for the ExplicitTaylorPolynomialStepper.\n");
00652   }
00653 
00654 
00655   template<class Scalar>
00656   int ExplicitTaylorPolynomialStepper<Scalar>::getOrder() const
00657   {
00658     return degree_;
00659   }
00660 
00661   //
00662   // Definitions of protected methods
00663   //
00664 
00665   template<class Scalar>
00666   void
00667   ExplicitTaylorPolynomialStepper<Scalar>::computeTaylorSeriesSolution_()
00668   {
00669     Teuchos::RCP<Thyra::VectorBase<Scalar> > tmp;
00670 
00671     // Set degree_ of polynomials to 0
00672     x_poly_->setDegree(0);
00673     f_poly_->setDegree(0);
00674 
00675     // Set degree_ 0 coefficient
00676     x_poly_->setCoefficient(0, *x_vector_);
00677 
00678     for (unsigned int k=1; k<=degree_; k++) {
00679 
00680       // compute [f] = f([x])
00681       Thyra::eval_f_poly(*model_, *x_poly_, t_, f_poly_.get());
00682 
00683       x_poly_->setDegree(k);
00684       f_poly_->setDegree(k);
00685       
00686       // x[k] = f[k-1] / k
00687       tmp = x_poly_->getCoefficient(k);
00688       copy(*(f_poly_->getCoefficient(k-1)), tmp.get());
00689       scale(Scalar(1.0)/Scalar(k), tmp.get());
00690     }
00691 
00692   }
00693 
00694   template<class Scalar>
00695   typename ExplicitTaylorPolynomialStepper<Scalar>::ScalarMag
00696   ExplicitTaylorPolynomialStepper<Scalar>::estimateLogRadius_()
00697   {
00698     ScalarMag rho = 0;
00699     ScalarMag tmp;
00700     for (unsigned int k=degree_/2; k<=degree_; k++) {
00701       tmp = log_norm_inf(*(x_poly_->getCoefficient(k))) / k;
00702       if (tmp > rho) {
00703         rho = tmp;
00704       }
00705     }
00706     return rho;
00707   }
00708 
00709   template<class Scalar>
00710   RCP<const Thyra::VectorSpaceBase<Scalar> > ExplicitTaylorPolynomialStepper<Scalar>::get_x_space() const
00711   {
00712     return(x_vector_->space());
00713   }
00714 
00715 } // namespace Rythmos
00716 
00717 #endif // RYTHMOS_EXPLICIT_TAYLOR_POLYNOMIAL_STEPPER_H

Generated on Tue Oct 20 12:46:07 2009 for Rythmos - Transient Integration for Differential Equations by doxygen 1.4.7