Rythmos_InterpolationBuffer.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_INTERPOLATION_BUFFER_H
00030 #define Rythmos_INTERPOLATION_BUFFER_H
00031 
00032 #include "Rythmos_InterpolationBufferBase.hpp"
00033 #include "Rythmos_InterpolatorBase.hpp"
00034 #include "Rythmos_LinearInterpolator.hpp"
00035 #include "Rythmos_DataStore.hpp"
00036 
00037 #include "Thyra_VectorBase.hpp"
00038 
00039 
00040 namespace Rythmos {
00041 
00043 template<class Scalar> 
00044 class InterpolationBuffer : virtual public InterpolationBufferBase<Scalar>
00045 {
00046   public:
00047 
00048     typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType ScalarMag;
00049 
00051 
00052     RCP<const Thyra::VectorSpaceBase<Scalar> > get_x_space() const;
00053     
00055     InterpolationBuffer();
00056     InterpolationBuffer( const RCP<InterpolatorBase<Scalar> >& interpolator_, int storage_ );
00057 
00059     void initialize( const RCP<InterpolatorBase<Scalar> >& interpolator_, int storage_ );
00060 
00062     void setInterpolator(const RCP<InterpolatorBase<Scalar> >& interpolator_);
00063     
00065     RCP<InterpolatorBase<Scalar> >& unSetInterpolator();
00066 
00068     void setStorage( int storage );
00069     
00071     int getStorage() const;
00072         
00074     ~InterpolationBuffer() {};
00075 
00077     void addPoints(
00078       const Array<Scalar>& time_vec
00079       ,const Array<RCP<const Thyra::VectorBase<Scalar> > >& x_vec
00080       ,const Array<RCP<const Thyra::VectorBase<Scalar> > >& xdot_vec);
00081 
00083     void getPoints(
00084       const Array<Scalar>& time_vec
00085       ,Array<RCP<const Thyra::VectorBase<Scalar> > >* x_vec
00086       ,Array<RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec
00087       ,Array<ScalarMag>* accuracy_vec
00088       ) const;
00089 
00091     TimeRange<Scalar> getTimeRange() const;
00092 
00094     void getNodes(Array<Scalar>* time_vec) const;
00095 
00097     int getOrder() const;
00098     
00100     void removeNodes(Array<Scalar>& time_vec);
00101 
00103 
00104     std::string description() const;
00105 
00107     void describe(
00108       Teuchos::FancyOStream       &out
00109       ,const Teuchos::EVerbosityLevel      verbLevel
00110       ) const;
00111 
00113 
00114     void setParameterList(RCP<Teuchos::ParameterList> const& paramList);
00115 
00117     RCP<Teuchos::ParameterList> getParameterList();
00118 
00120     RCP<Teuchos::ParameterList> unsetParameterList();
00121 
00122     enum IBPolicy {
00123       BUFFER_STATIC = 0
00124       ,BUFFER_KEEP_NEWEST = 1
00125     };
00126     
00127   private:
00128 
00129     RCP<InterpolatorBase<Scalar> > interpolator_;
00130     int storage_limit_;
00131     typename DataStore<Scalar>::DataStoreVector_t data_vec_;
00132 
00133     RCP<Teuchos::ParameterList> parameterList_;
00134 
00135     IBPolicy policy_;
00136 
00137 };
00138 
00139 // ////////////////////////////
00140 // Defintions
00141 template<class Scalar>
00142 InterpolationBuffer<Scalar>::InterpolationBuffer()
00143 {
00144   initialize(Teuchos::null,0);
00145 }
00146 
00147 template<class Scalar>
00148 InterpolationBuffer<Scalar>::InterpolationBuffer( 
00149     const RCP<InterpolatorBase<Scalar> >& interpolator_
00150     ,int storage_ 
00151     )
00152 {
00153   initialize(interpolator_,storage_);
00154 }
00155 
00156 template<class Scalar>
00157 RCP<const Thyra::VectorSpaceBase<Scalar> > InterpolationBuffer<Scalar>::get_x_space() const
00158 {
00159   if (data_vec_.size() == 0) {
00160     RCP<const Thyra::VectorSpaceBase<Scalar> > space;
00161     return(space);
00162   } else {
00163     return(data_vec_[0].x->space());
00164   }
00165 }
00166 
00167 template<class Scalar>
00168 void InterpolationBuffer<Scalar>::initialize( 
00169     const RCP<InterpolatorBase<Scalar> >& interpolator_
00170     ,int storage_
00171     )
00172 {
00173   RCP<Teuchos::FancyOStream> out = this->getOStream();
00174   Teuchos::OSTab ostab(out,1,"IB::initialize");
00175   *out << "Initializing InterpolationBuffer" << std::endl;
00176   if ( Teuchos::as<int>(this->getVerbLevel()) >= Teuchos::as<int>(Teuchos::VERB_HIGH) ) {
00177     *out << "Calling setInterpolator..." << std::endl;
00178   }
00179   setInterpolator(interpolator_);
00180   if ( Teuchos::as<int>(this->getVerbLevel()) >= Teuchos::as<int>(Teuchos::VERB_HIGH) ) {
00181     *out << "Calling setStorage..." << std::endl;
00182   }
00183   setStorage(storage_);
00184   policy_ = BUFFER_KEEP_NEWEST;
00185 }
00186 
00187 template<class Scalar>
00188 void InterpolationBuffer<Scalar>::setStorage( int storage )
00189 {
00190   int storage_limit = std::max(2,storage); // Minimum of two points so interpolation is possible
00191   TEST_FOR_EXCEPTION(
00192       Teuchos::as<int>(data_vec_.size()) > storage_limit, std::logic_error,
00193       "Error, specified storage = " << storage_limit << " is below current number of vectors stored = " << data_vec_.size() << "!\n"
00194       );
00195   storage_limit_ = storage_limit;
00196   RCP<Teuchos::FancyOStream> out = this->getOStream();
00197   Teuchos::OSTab ostab(out,1,"IB::setStorage");
00198   if ( Teuchos::as<int>(this->getVerbLevel()) >= Teuchos::as<int>(Teuchos::VERB_HIGH) ) {
00199     *out << "storage_limit = " << storage_limit_ << std::endl;
00200   }
00201 }
00202 
00203 template<class Scalar>
00204 int InterpolationBuffer<Scalar>::getStorage() const
00205 {
00206   return(storage_limit_);
00207 }
00208 
00209 template<class Scalar>
00210 void InterpolationBuffer<Scalar>::setInterpolator(
00211     const RCP<InterpolatorBase<Scalar> >& interpolator
00212     )
00213 {
00214   if (interpolator_ == Teuchos::null) {
00215     interpolator_ = Teuchos::rcp(new LinearInterpolator<Scalar>);
00216   } else {
00217     interpolator_ = interpolator;
00218   }
00219   RCP<Teuchos::FancyOStream> out = this->getOStream();
00220   Teuchos::OSTab ostab(out,1,"IB::setInterpolator");
00221   if ( Teuchos::as<int>(this->getVerbLevel()) >= Teuchos::as<int>(Teuchos::VERB_HIGH) ) {
00222     *out << "interpolator = " << interpolator_->description() << std::endl;
00223   }
00224 }
00225 
00226 template<class Scalar>
00227 void InterpolationBuffer<Scalar>::addPoints( 
00228     const Array<Scalar>& time_vec
00229     ,const Array<RCP<const Thyra::VectorBase<Scalar> > >& x_vec
00230     ,const Array<RCP<const Thyra::VectorBase<Scalar> > >& xdot_vec 
00231     )
00232 {
00233 #ifdef TEUCHOS_DEBUG
00234   // Check preconditions
00235   assertTimePointsAreSorted(time_vec);
00236   int tsize = time_vec.size();
00237   TEST_FOR_EXCEPTION(
00238       tsize == 0, std::logic_error,
00239       "Error, time_vec is empty!"
00240       );
00241   TEST_FOR_EXCEPTION(
00242       Teuchos::as<int>(x_vec.size()) == tsize, std::logic_error,
00243       "Error, size of x_vec = " << x_vec.size() << " != " << tsize << " = size of time_vec!\n"
00244       );
00245   TEST_FOR_EXCEPTION(
00246       Teuchos::as<int>(xdot_vec.size()) == tsize, std::logic_error,
00247       "Error, size of xdot_vec = " << x_vec.size() << " != " << tsize << " = size of time_vec!\n"
00248       );
00249   for (int i=0; i<tsize ; ++i) {
00250     TEST_FOR_EXCEPTION(
00251         x_vec[i] == Teuchos::null, std::logic_error,
00252         "Error, x_vec[" << i << "] == null!\n"
00253         );
00254     TEST_FOR_EXCEPTION(
00255         xdot_vec[i] == Teuchos::null, std::logic_error,
00256         "Error, xdot_vec[" << i << "] == null!\n"
00257         );
00258   }
00259   assertNoTimePointsInsideCurrentTimeRange(*this,time_vec);
00260 #endif // TEUCHOS_DEBUG
00261   RCP<Teuchos::FancyOStream> out = this->getOStream();
00262   Teuchos::OSTab ostab(out,1,"IB::setPoints");
00263   if ( Teuchos::as<int>(this->getVerbLevel()) >= Teuchos::as<int>(Teuchos::VERB_HIGH) ) {
00264     *out << "time_vec = " << std::endl;
00265     for (unsigned int i=0 ; i<time_vec.size() ; ++i) {
00266       *out << "time_vec[" << i << "] = " << time_vec[i] << std::endl;
00267     }
00268     *out << "x_vec = " << std::endl;
00269     for (unsigned int i=0 ; i<x_vec.size() ; ++i) {
00270       *out << "x_vec[" << i << "] = " << std::endl;
00271       x_vec[i]->describe(*out,Teuchos::VERB_EXTREME);
00272     }
00273     *out << "xdot_vec = " << std::endl;
00274     for (unsigned int i=0 ; i<xdot_vec.size() ; ++i) {
00275       *out << "xdot_vec[" << i << "] = " << std::endl;
00276       xdot_vec[i]->describe(*out,Teuchos::VERB_EXTREME);
00277     }
00278   }
00279   typename DataStore<Scalar>::DataStoreList_t input_data_list;
00280   vectorToDataStoreList<Scalar>(time_vec,x_vec,xdot_vec,&input_data_list);
00281   // Check that we're not going to exceed our storage limit:
00282   if ((data_vec_.size()+input_data_list.size()) > Teuchos::as<unsigned int>(storage_limit_)) { 
00283     if (policy_ == BUFFER_STATIC) {
00284       TEST_FOR_EXCEPTION(
00285           true, std::logic_error,
00286           "Error, buffer is full and buffer policy is BUFFER_STATIC, no points can be added\n"
00287           );
00288     } else if (policy_ == BUFFER_KEEP_NEWEST) {
00289       if (input_data_list.front() > data_vec_.back()) {
00290         // Case:  all of new points are past end of existing points
00291         // Remove points from the beginning of data_vec, then add new points
00292         int num_extra_points = input_data_list.size();
00293         typename DataStore<Scalar>::DataStoreVector_t::iterator 
00294           data_it = data_vec_.begin();
00295         for (int i=0 ; i < num_extra_points ; ++i) {
00296           data_it++;
00297         }
00298         if ( Teuchos::as<int>(this->getVerbLevel()) >= Teuchos::as<int>(Teuchos::VERB_HIGH) ) {
00299           *out << "Removing " << num_extra_points 
00300             << " from beginning of data_vec to make room for new points." << std::endl;
00301         }
00302         data_vec_.erase(data_vec_.begin(),data_it);
00303       } else if (input_data_list.back() < data_vec_.front()) {
00304         // Case:  all of new points are before beginning of existing points
00305         // Remove points from end of data_vec, then add new points
00306         int num_extra_points = input_data_list.size();
00307         typename DataStore<Scalar>::DataStoreVector_t::iterator 
00308           data_it = data_vec_.end();
00309         for (int i=0 ; i < num_extra_points ; ++i) {
00310           data_it--;
00311         }
00312         if ( Teuchos::as<int>(this->getVerbLevel()) >= Teuchos::as<int>(Teuchos::VERB_HIGH) ) {
00313           *out << "Removing " << num_extra_points 
00314             << " from end of data_vec to make room for new points." << std::endl;
00315         }
00316         data_vec_.erase(data_it,data_vec_.end());
00317       } else {
00318         // Case:  Some points are before beginning of data_vec and some points are after end of data_vec
00319         TEST_FOR_EXCEPTION(
00320             true, std::logic_error,
00321             "Error, incoming points are on both sides of TimeRange, this feature not implemented yet.\n"
00322             );
00323       }
00324     } else {
00325       // Unknown Buffer policy:
00326       TEST_FOR_EXCEPTION(
00327           true, std::logic_error,
00328           "Error, unknown buffer policy.\n"
00329           );
00330     }
00331   }
00332   // Now add all the remaining points to data_vec
00333   data_vec_.insert(data_vec_.end(),input_data_list.begin(),input_data_list.end());
00334   // And sort data_vec:
00335   std::sort(data_vec_.begin(),data_vec_.end());
00336   if ( Teuchos::as<int>(this->getVerbLevel()) >= Teuchos::as<int>(Teuchos::VERB_HIGH) ) {
00337     *out << "data_vec at end of setPoints:" << std::endl;
00338     for (unsigned int i=0 ; i<data_vec_.size() ; ++i) {
00339       *out << "data_vec[" << i << "] = " << std::endl;
00340       data_vec_[i].describe(*out,Teuchos::VERB_EXTREME);
00341     }
00342   }
00343 }
00344 
00345 template<class Scalar>
00346 void InterpolationBuffer<Scalar>::getPoints(
00347     const Array<Scalar>& time_vec
00348     ,Array<RCP<const Thyra::VectorBase<Scalar> > >* x_vec
00349     ,Array<RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec
00350     ,Array<ScalarMag>* accuracy_vec
00351     ) const
00352 {
00353   RCP<Teuchos::FancyOStream> out = this->getOStream();
00354   Teuchos::OSTab ostab(out,1,"IB::getPoints");
00355   typename DataStore<Scalar>::DataStoreVector_t data_out;
00356   interpolator_->interpolate(data_vec_, time_vec, &data_out);
00357   Array<Scalar> time_out_vec;
00358   dataStoreVectorToVector<Scalar>(data_out, &time_out_vec, x_vec, xdot_vec, accuracy_vec);
00359   TEST_FOR_EXCEPTION(
00360     (time_vec.size() != time_out_vec.size()), std::logic_error,
00361     "Error, number of time points returned from interpolator = " <<
00362       time_out_vec.size() << " != " << time_vec.size() << 
00363       " = number of time points requested\n"
00364     );
00365   if ( Teuchos::as<int>(this->getVerbLevel()) >= Teuchos::as<int>(Teuchos::VERB_HIGH) ) {
00366     *out << "Conversion of DataStoreVector to Vector successful" << std::endl;
00367   }
00368 }
00369 
00370 
00371 template<class Scalar>
00372 TimeRange<Scalar> InterpolationBuffer<Scalar>::getTimeRange() const
00373 {
00374   TimeRange<Scalar> timerange(data_vec_.front().time,data_vec_.back().time);
00375   return(timerange);
00376 }
00377 
00378 template<class Scalar>
00379 void InterpolationBuffer<Scalar>::getNodes( Array<Scalar>* time_vec ) const
00380 {
00381   int N = data_vec_.size();
00382   time_vec->clear();
00383   time_vec->reserve(N);
00384   for (int i=0 ; i<N ; ++i) {
00385     time_vec->push_back(data_vec_[i].time);
00386   }
00387   RCP<Teuchos::FancyOStream> out = this->getOStream();
00388   Teuchos::OSTab ostab(out,1,"IB::getNodes");
00389   if ( Teuchos::as<int>(this->getVerbLevel()) >= Teuchos::as<int>(Teuchos::VERB_HIGH) ) {
00390     *out << this->description() << std::endl;
00391     for (unsigned int i=0 ; i<time_vec->size() ; ++i) {
00392       *out << "time_vec[" << i << "] = " << (*time_vec)[i] << std::endl;
00393     }
00394   }
00395 }
00396 
00397 template<class Scalar>
00398 void InterpolationBuffer<Scalar>::removeNodes( Array<Scalar>& time_vec ) 
00399 {
00400   typedef Teuchos::ScalarTraits<Scalar> ST;
00401   int N = time_vec.size();
00402 #ifdef TEUCHOS_DEBUG
00403   // Check preconditions:
00404   TimeRange<Scalar> range = this->getTimeRange();
00405   for (int i=0; i<N ; ++i) {
00406     TEST_FOR_EXCEPTION(
00407         ~(range.lower() <= time_vec[i]) && (time_vec[i] <= range.upper()),
00408         std::logic_error,
00409         "Error, time_vec[" << i << "] = " << time_vec[i] << 
00410           "is not in range of this interpolation buffer = [" << 
00411           range.lower() << "," << range.upper() << "]!\n"
00412         );
00413   }
00414 #endif // TEUCHOS_DEBUG
00415   RCP<Thyra::VectorBase<Scalar> > vec_temp;
00416   ScalarMag z = ST::zero();
00417   for (int i=0; i<N ; ++i) {
00418     DataStore<Scalar> ds_temp(time_vec[i],vec_temp,vec_temp,z);
00419     typename DataStore<Scalar>::DataStoreVector_t::iterator 
00420       data_it = std::find(data_vec_.begin(),data_vec_.end(),ds_temp);
00421     TEST_FOR_EXCEPTION(
00422         data_it == data_vec_.end(), std::logic_error,
00423         "Error, time_vec[" << i << "] = " << time_vec[i] << "is not a node in the interpolation buffer!\n"
00424         );
00425     data_vec_.erase(data_it);
00426   }
00427 }
00428 
00429 template<class Scalar>
00430 int InterpolationBuffer<Scalar>::getOrder() const
00431 {
00432   return(interpolator_->order());
00433 }
00434 
00435 template<class Scalar>
00436 std::string InterpolationBuffer<Scalar>::description() const
00437 {
00438   std::string name = "Rythmos::InterpolationBuffer";
00439   return(name);
00440 }
00441 
00442 template<class Scalar>
00443 void InterpolationBuffer<Scalar>::describe(
00444       Teuchos::FancyOStream                &out
00445       ,const Teuchos::EVerbosityLevel      verbLevel
00446       ) const
00447 {
00448   if ( (Teuchos::as<int>(verbLevel) == Teuchos::as<int>(Teuchos::VERB_DEFAULT) ) ||
00449        (Teuchos::as<int>(verbLevel) >= Teuchos::as<int>(Teuchos::VERB_LOW)     )
00450      ) {
00451     out << description() << "::describe" << std::endl;
00452     out << "interpolator = " << interpolator_->description() << std::endl;
00453     out << "storage_limit = " << storage_limit_ << std::endl;
00454   } else if (Teuchos::as<int>(verbLevel) >= Teuchos::as<int>(Teuchos::VERB_LOW)) {
00455   } else if (Teuchos::as<int>(verbLevel) >= Teuchos::as<int>(Teuchos::VERB_MEDIUM)) {
00456   } else if (Teuchos::as<int>(verbLevel) >= Teuchos::as<int>(Teuchos::VERB_HIGH)) {
00457     out << "data_vec = " << std::endl;
00458     for (unsigned int i=0; i<data_vec_.size() ; ++i) {
00459       out << "data_vec[" << i << "] = " << std::endl;
00460       data_vec_[i].describe(out,this->getVerbLevel());
00461     }
00462   }
00463 }
00464 
00465 template <class Scalar>
00466 void InterpolationBuffer<Scalar>::setParameterList(RCP<Teuchos::ParameterList> const& paramList)
00467 {
00468   parameterList_ = paramList;
00469   int outputLevel = parameterList_->get( "outputLevel", int(-1) );
00470   outputLevel = std::min(std::max(outputLevel,-1),4);
00471   this->setVerbLevel(static_cast<Teuchos::EVerbosityLevel>(outputLevel));
00472   int policyLevel = parameterList_->get( "InterpolationBufferPolicy", int(1) );
00473   policyLevel = std::min(std::max(policyLevel,0),1);
00474   policy_ = static_cast<IBPolicy>(policyLevel);
00475   int storage_limit = parameterList_->get( "StorageLimit", storage_limit_ );
00476   if (storage_limit != storage_limit_) {
00477     this->setStorage(storage_limit);
00478   }
00479 }
00480 
00481 template <class Scalar>
00482 RCP<Teuchos::ParameterList> InterpolationBuffer<Scalar>::getParameterList()
00483 {
00484   return(parameterList_);
00485 }
00486 
00487 template <class Scalar>
00488 RCP<Teuchos::ParameterList> InterpolationBuffer<Scalar>::unsetParameterList()
00489 {
00490   RCP<Teuchos::ParameterList> temp_param_list = parameterList_;
00491   parameterList_ = Teuchos::null;
00492   return(temp_param_list);
00493 }
00494 
00495 } // namespace Rythmos
00496 
00497 #endif // Rythmos_INTERPOLATION_BUFFER_H

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