Rythmos_ExplicitTaylorPolynomialStepper.hpp

Go to the documentation of this file.
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_Stepper.hpp"
00033 #include "Teuchos_RefCountPtr.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   typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType ScalarMag;
00140 
00141   template<class Scalar>
00142   class ExplicitTaylorPolynomialStepper : public Stepper<Scalar>
00143   {
00144   public:
00145 
00147     typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitude_type;
00148     
00153     ExplicitTaylorPolynomialStepper(const Teuchos::RefCountPtr<const Thyra::ModelEvaluator<Scalar> > &model, Teuchos::ParameterList& params);
00154     
00156     ~ExplicitTaylorPolynomialStepper();
00157     
00159     Scalar TakeStep(Scalar dt);
00160    
00162     Scalar TakeStep();
00163 
00165     Teuchos::RefCountPtr<const Thyra::VectorBase<Scalar> > get_solution() const;
00166 
00168     std::string description() const;
00169 
00171     std::ostream& describe(
00172       std::ostream                &out
00173       ,const Teuchos::EVerbosityLevel      verbLevel
00174       ,const std::string          leadingIndent
00175       ,const std::string          indentSpacer
00176       ) const;
00177 
00180     bool SetPoints(
00181       const std::vector<Scalar>& time_list
00182       ,const std::vector<Thyra::VectorBase<Scalar> >& x_list
00183       ,const std::vector<Thyra::VectorBase<Scalar> >& xdot_list);
00184 
00186     bool GetPoints(
00187       const std::vector<Scalar>& time_list
00188       ,std::vector<Thyra::VectorBase<Scalar> >* x_list
00189       ,std::vector<Thyra::VectorBase<Scalar> >* xdot_list
00190       ,std::vector<ScalarMag>* accuracy_list) const;
00191 
00193     bool SetRange(
00194       const Scalar& time_lower
00195       ,const Scalar& time_upper
00196       ,const InterpolationBuffer<Scalar> & IB);
00197 
00199     bool GetNodes(std::vector<Scalar>* time_list) const;
00200 
00202     bool RemoveNodes(std::vector<Scalar>* time_list) const;
00203 
00205     int GetOrder() const;
00206 
00207   protected:
00208 
00210     void computeTaylorSeriesSolution();
00211 
00216     magnitude_type estimateLogRadius();
00217 
00218   protected:
00219 
00221     Teuchos::RefCountPtr<const Thyra::ModelEvaluator<Scalar> > model;
00222 
00224     Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > x_vector;
00225 
00227     Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > x_dot_vector;
00228 
00230     Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > f_vector;
00231 
00233     Teuchos::RefCountPtr<Teuchos::Polynomial<Thyra::VectorBase<Scalar> > > x_poly;
00234 
00236     Teuchos::RefCountPtr<Teuchos::Polynomial<Thyra::VectorBase<Scalar> > > f_poly;
00237 
00239     Scalar t;
00240 
00242     Scalar t_initial;
00243 
00245     Scalar t_final;
00246 
00248     magnitude_type local_error_tolerance;
00249 
00251     Scalar min_step_size;
00252 
00254     Scalar max_step_size;
00255 
00257     unsigned int degree;
00258 
00260     Scalar linc;
00261   };
00262 
00264 
00271   template <typename Scalar>
00272   class ROpLogNormInf : public RTOpPack::ROpScalarReductionBase<typename Teuchos::ScalarTraits<Scalar>::magnitudeType> {
00273   public:
00274 
00276     ROpLogNormInf() : RTOpPack::RTOpT<Scalar>("ROpLogInfNorm"){}
00277 
00279     typename Teuchos::ScalarTraits<Scalar>::magnitudeType 
00280     operator() (const RTOpPack::ReductTarget& reduct_obj) const { 
00281       return this->getRawVal(reduct_obj); 
00282     }
00283 
00286 
00288     void reduce_reduct_objs(const RTOpPack::ReductTarget& _in_reduct_obj, 
00289           RTOpPack::ReductTarget* _inout_reduct_obj) const
00290     {
00291       const RTOpPack::ReductTargetScalar<Scalar>& in_reduct_obj = 
00292   Teuchos::dyn_cast<const RTOpPack::ReductTargetScalar<Scalar> >(_in_reduct_obj);
00293       RTOpPack::ReductTargetScalar<Scalar>& inout_reduct_obj = 
00294   Teuchos::dyn_cast<RTOpPack::ReductTargetScalar<Scalar> >(*_inout_reduct_obj);
00295       if(in_reduct_obj.get() > inout_reduct_obj.get()) 
00296   inout_reduct_obj.set(in_reduct_obj.get());
00297     }
00298 
00300     void apply_op(const int num_vecs, 
00301       const RTOpPack::SubVectorT<Scalar> sub_vecs[],
00302       const int num_targ_vecs, 
00303       const RTOpPack::MutableSubVectorT<Scalar> targ_sub_vecs[], 
00304       RTOpPack::ReductTarget *_reduct_obj) const
00305     {
00306       RTOpPack::ReductTargetScalar<Scalar>& reduct_obj = 
00307   Teuchos::dyn_cast<RTOpPack::ReductTargetScalar<Scalar> >(*_reduct_obj);
00308  
00309       RTOP_APPLY_OP_1_0(num_vecs,sub_vecs,num_targ_vecs,targ_sub_vecs);
00310 
00311       typename Teuchos::ScalarTraits<Scalar>::magnitudeType norm_inf = 
00312   reduct_obj.get();
00313 
00314       // unit stride
00315       if(v0_s == 1)
00316         for(RTOp_index_type i=0; i<subDim; i++) {
00317           typename Teuchos::ScalarTraits<Scalar>::magnitudeType
00318             mag = Teuchos::ScalarTraits<Scalar>::magnitude(*v0_val++);
00319     mag = std::log(Teuchos::ScalarTraits<Scalar>::one() + mag);
00320           norm_inf = mag > norm_inf ? mag : norm_inf;
00321         }
00322       else 
00323         for(RTOp_index_type i=0; i<subDim; i++, v0_val+=v0_s) {
00324           typename Teuchos::ScalarTraits<Scalar>::magnitudeType
00325             mag = Teuchos::ScalarTraits<Scalar>::magnitude(*v0_val);
00326     mag = std::log(Teuchos::ScalarTraits<Scalar>::one() + mag);
00327           norm_inf = mag > norm_inf ? mag : norm_inf;
00328         }
00329 
00330       reduct_obj.set(norm_inf);
00331     }
00333 
00334   }; // class ROpLogNormInf
00335 
00337   template <typename Scalar>
00338   typename Teuchos::ScalarTraits<Scalar>::magnitudeType
00339   log_norm_inf(const Thyra::VectorBase<Scalar>& x)
00340   {
00341     ROpLogNormInf<Scalar> log_norm_inf_op;
00342     Teuchos::RefCountPtr<RTOpPack::ReductTarget> log_norm_inf_targ = 
00343       log_norm_inf_op.reduct_obj_create();
00344     const Thyra::VectorBase<Scalar>* vecs[] = { &x };
00345     Thyra::applyOp<Scalar>(log_norm_inf_op,1,vecs,0,
00346          (Thyra::VectorBase<Scalar>**)NULL,
00347          log_norm_inf_targ.get());
00348     
00349     return log_norm_inf_op(*log_norm_inf_targ);
00350   }
00351 
00352   template<class Scalar>
00353   ExplicitTaylorPolynomialStepper<Scalar>::ExplicitTaylorPolynomialStepper(
00354   const Teuchos::RefCountPtr<const Thyra::ModelEvaluator<Scalar> > &m,
00355   Teuchos::ParameterList& params)
00356   {
00357     typedef Teuchos::ScalarTraits<Scalar> ST;
00358 
00359     // Get initial time
00360     t_initial = params.get("Initial Time", ST::zero());
00361 
00362     // Get final time
00363     t_final = params.get("Final Time", ST::one());
00364 
00365     // Get local error tolerance
00366     local_error_tolerance = 
00367       params.get("Local Error Tolerance", magnitude_type(1.0e-10));
00368 
00369     // Get minimum step size
00370     min_step_size = params.get("Minimum Step Size", Scalar(1.0e-10));
00371 
00372     // Get maximum step size
00373     max_step_size = params.get("Maximum Step Size", Scalar(1.0));
00374 
00375     // Get degree of Taylor polynomial expansion
00376     degree = params.get("Taylor Polynomial Degree", 40);
00377 
00378     linc = Scalar(-16.0*std::log(10.0)/degree);
00379   
00380     model = m;
00381     t = t_initial;
00382     x_vector = model->getNominalValues().get_x()->clone_v();
00383     x_dot_vector = Thyra::createMember(model->get_x_space());
00384     f_vector = Thyra::createMember(model->get_f_space());
00385     x_poly = 
00386       Teuchos::rcp(new Teuchos::Polynomial<Thyra::VectorBase<Scalar> >(0, 
00387                      *x_vector,
00388                      degree));
00389     f_poly = 
00390       Teuchos::rcp(new Teuchos::Polynomial<Thyra::VectorBase<Scalar> >(0, 
00391                      *f_vector,
00392                      degree));
00393   }
00394 
00395   template<class Scalar>
00396   ExplicitTaylorPolynomialStepper<Scalar>::~ExplicitTaylorPolynomialStepper()
00397   {
00398   }
00399 
00400   template<class Scalar>
00401   Scalar 
00402   ExplicitTaylorPolynomialStepper<Scalar>::TakeStep()
00403   {
00404     // If t > t_final, we're done
00405     if (t > t_final)
00406       return Teuchos::ScalarTraits<Scalar>::zero();
00407 
00408     // Compute a local truncated Taylor series solution to system
00409     computeTaylorSeriesSolution();
00410 
00411     // Estimate log of radius of convergence of Taylor series
00412     Scalar rho = estimateLogRadius();
00413 
00414     // Set step size
00415     Scalar dt = std::exp(linc - rho);
00416 
00417     // If step size is too big, reduce
00418     if (dt > max_step_size)
00419       dt = max_step_size;
00420 
00421     // If step goes past t_final, reduce
00422     if (t+dt > t_final)
00423       dt = t_final-t;
00424 
00425     magnitude_type local_error;
00426 
00427     do {
00428 
00429       // compute x(t+dt), xdot(t+dt)
00430       x_poly->evaluate(dt, x_vector.get(), x_dot_vector.get());
00431 
00432       // compute f( x(t+dt), t+dt )
00433       Thyra::eval_f(*model, *x_vector, t+dt, f_vector.get());
00434 
00435       // compute || xdot(t+dt) - f( x(t+dt), t+dt ) ||
00436       Thyra::Vp_StV(x_dot_vector.get(), -Teuchos::ScalarTraits<Scalar>::one(),
00437         *f_vector);
00438       local_error = norm_inf(*x_dot_vector);
00439 
00440       if (local_error > local_error_tolerance)
00441   dt *= 0.7;
00442 
00443     }
00444     while (local_error > local_error_tolerance && dt > min_step_size);
00445 
00446     // Check if minimum step size was reached
00447     TEST_FOR_EXCEPTION(dt < min_step_size, 
00448            std::runtime_error,
00449            "ExplicitTaylorPolynomialStepper<Scalar>::takeStep(): " 
00450            << "Step size reached minimum step size " 
00451            << min_step_size << ".  Failing step." );
00452 
00453     // Increment t
00454     t += dt;
00455 
00456     return dt;
00457   }
00458 
00459   template<class Scalar>
00460   Scalar 
00461   ExplicitTaylorPolynomialStepper<Scalar>::TakeStep(Scalar dt)
00462   {
00463     // If t > t_final, we're done
00464     if (t > t_final)
00465       return Teuchos::ScalarTraits<Scalar>::zero();
00466 
00467     // Compute a local truncated Taylor series solution to system
00468     computeTaylorSeriesSolution();
00469 
00470     // If step size is too big, reduce
00471     if (dt > max_step_size)
00472       dt = max_step_size;
00473 
00474     // If step goes past t_final, reduce
00475     if (t+dt > t_final)
00476       dt = t_final-t;
00477 
00478     // compute x(t+dt)
00479     x_poly->evaluate(dt, x_vector.get());
00480 
00481     // Increment t
00482     t += dt;
00483 
00484     return dt;
00485   }
00486 
00487   template<class Scalar>
00488   Teuchos::RefCountPtr<const Thyra::VectorBase<Scalar> >
00489   ExplicitTaylorPolynomialStepper<Scalar>::get_solution() const
00490   {
00491     return x_vector;
00492   }
00493 
00494   template<class Scalar>
00495   void
00496   ExplicitTaylorPolynomialStepper<Scalar>::computeTaylorSeriesSolution()
00497   {
00498     Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > tmp;
00499 
00500     // Set degree of polynomials to 0
00501     x_poly->setDegree(0);
00502     f_poly->setDegree(0);
00503 
00504     // Set degree 0 coefficient
00505     x_poly->setCoefficient(0, *x_vector);
00506 
00507     for (unsigned int k=1; k<=degree; k++) {
00508 
00509       // compute [f] = f([x])
00510       Thyra::eval_f_poly(*model, *x_poly, t, f_poly.get());
00511 
00512       x_poly->setDegree(k);
00513       f_poly->setDegree(k);
00514       
00515       // x[k] = f[k-1] / k
00516       tmp = x_poly->getCoefficient(k);
00517       copy(*(f_poly->getCoefficient(k-1)), tmp.get());
00518       scale(Scalar(1.0)/Scalar(k), tmp.get());
00519     }
00520 
00521   }
00522 
00523   template<class Scalar>
00524   typename ExplicitTaylorPolynomialStepper<Scalar>::magnitude_type
00525   ExplicitTaylorPolynomialStepper<Scalar>::estimateLogRadius()
00526   {
00527     magnitude_type rho = 0;
00528     magnitude_type tmp;
00529     for (unsigned int k=degree/2; k<=degree; k++) {
00530       tmp = log_norm_inf(*(x_poly->getCoefficient(k))) / k;
00531       if (tmp > rho)
00532   rho = tmp;
00533     }
00534     return rho;
00535   }
00536 
00537   template<class Scalar>
00538   std::string ExplicitTaylorPolynomialStepper<Scalar>::description() const
00539   {
00540     std::string name = "Rythmos::ExplicitTaylorPolynomialStepper";
00541     return(name);
00542   }
00543 
00544   template<class Scalar>
00545   std::ostream& ExplicitTaylorPolynomialStepper<Scalar>::describe(
00546         std::ostream                &out
00547         ,const Teuchos::EVerbosityLevel      verbLevel
00548         ,const std::string          leadingIndent
00549         ,const std::string          indentSpacer
00550         ) const
00551   {
00552     if (verbLevel == Teuchos::VERB_EXTREME)
00553     {
00554       out << description() << "::describe" << std::endl;
00555       out << "model = " << std::endl;
00556       out << model->describe(out,verbLevel,leadingIndent,indentSpacer) << std::endl;
00557       out << "x_vector = " << std::endl;
00558       out << x_vector->describe(out,verbLevel,leadingIndent,indentSpacer) << std::endl;
00559       out << "x_dot_vector = " << std::endl;
00560       out << x_dot_vector->describe(out,verbLevel,leadingIndent,indentSpacer) << std::endl;
00561       out << "f_vector = " << std::endl;
00562       out << f_vector->describe(out,verbLevel,leadingIndent,indentSpacer) << std::endl;
00563       out << "x_poly = " << std::endl;
00564       out << x_poly->describe(out,verbLevel,leadingIndent,indentSpacer) << std::endl;
00565       out << "f_poly = " << std::endl;
00566       out << f_poly->describe(out,verbLevel,leadingIndent,indentSpacer) << std::endl;
00567       out << "t = " << t << std::endl;
00568       out << "t_initial = " << t_initial << std::endl;
00569       out << "t_final = " << t_final << std::endl;
00570       out << "local_error_tolerance = " << local_error_tolerance << std::endl;
00571       out << "min_step_size = " << min_step_size << std::endl;
00572       out << "max_step_size = " << max_step_size << std::endl;
00573       out << "degree = " << degree << std::endl;
00574       out << "linc = " << linc << std::endl;
00575     }
00576     return(out);
00577   }
00578 
00579 
00580 template<class Scalar>
00581 bool ExplicitTaylorPolynomialStepper<Scalar>::SetPoints(
00582     const std::vector<Scalar>& time_list
00583     ,const std::vector<Thyra::VectorBase<Scalar> >& x_list
00584     ,const std::vector<Thyra::VectorBase<Scalar> >& xdot_list)
00585 {
00586   return(false);
00587 }
00588 
00589 template<class Scalar>
00590 bool ExplicitTaylorPolynomialStepper<Scalar>::GetPoints(
00591     const std::vector<Scalar>& time_list
00592     ,std::vector<Thyra::VectorBase<Scalar> >* x_list
00593     ,std::vector<Thyra::VectorBase<Scalar> >* xdot_list
00594     ,std::vector<ScalarMag>* accuracy_list) const
00595 {
00596   return(false);
00597 }
00598 
00599 template<class Scalar>
00600 bool ExplicitTaylorPolynomialStepper<Scalar>::SetRange(
00601     const Scalar& time_lower
00602     ,const Scalar& time_upper
00603     ,const InterpolationBuffer<Scalar>& IB)
00604 {
00605   return(false);
00606 }
00607 
00608 template<class Scalar>
00609 bool ExplicitTaylorPolynomialStepper<Scalar>::GetNodes(std::vector<Scalar>* time_list) const
00610 {
00611   return(false);
00612 }
00613 
00614 template<class Scalar>
00615 bool ExplicitTaylorPolynomialStepper<Scalar>::RemoveNodes(std::vector<Scalar>* time_list) const
00616 {
00617   return(false);
00618 }
00619 
00620 
00621 template<class Scalar>
00622 int ExplicitTaylorPolynomialStepper<Scalar>::GetOrder() const
00623 {
00624   return(4);
00625 }
00626 
00627 
00628 } // namespace Rythmos
00629 
00630 #endif // RYTHMOS_EXPLICIT_TAYLOR_POLYNOMIAL_STEPPER_H

Generated on Thu Sep 18 12:30:05 2008 for Rythmos - Transient Integration for Differential Equations by doxygen 1.3.9.1