Rythmos_IntegratorDefault.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_INTEGRATOR_DEFAULT_H
00030 #define Rythmos_INTEGRATOR_DEFAULT_H
00031 
00032 #include "Rythmos_IntegratorBase.hpp"
00033 #include "Rythmos_InterpolationBufferBase.hpp"
00034 #include "Rythmos_InterpolationBufferAppender.hpp"
00035 #include "Rythmos_BreakPointInformerBase.hpp"
00036 #include "Rythmos_IntegrationObserverBase.hpp"
00037 #include "Rythmos_StepperBase.hpp"
00038 #include "Rythmos_InterpolationBufferHelpers.hpp"
00039 
00040 #include "Teuchos_VerboseObjectParameterListHelpers.hpp"
00041 #include "Teuchos_as.hpp"
00042 
00043 
00044 namespace Rythmos {
00045 
00046 
00050 template<class Scalar> 
00051 class IntegratorDefault : virtual public IntegratorBase<Scalar>
00052 {
00053 public:
00054   
00056   typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType ScalarMag;
00057 
00060   
00062   ~IntegratorDefault() {};
00063  
00065   IntegratorDefault();
00066   
00068   IntegratorDefault(
00069     const Teuchos::RCP<StepperBase<Scalar> > &stepper,
00070     const Teuchos::RCP<InterpolationBufferBase<Scalar> > &trailingInterpBuffer,
00071     const Teuchos::RCP<Teuchos::ParameterList> &paramList = Teuchos::null
00072     );
00073   
00075   void setTrailingInterpolationBuffer(
00076     const Teuchos::RCP<InterpolationBufferBase<Scalar> > &trailingInterpBuffer
00077     );
00078 
00080   bool acceptsTrailingInterpolationBuffer() const;
00081 
00083   Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > get_x_space() const;
00084 
00086   Teuchos::RCP<const InterpolationBufferBase<Scalar> > getInterpolationBuffer();
00087 
00089   Teuchos::RCP<InterpolationBufferBase<Scalar> > unSetInterpolationBuffer();
00090 
00092   void setInterpolationBufferAppender(
00093     const Teuchos::RCP<InterpolationBufferAppenderBase<Scalar> > &interpolationBufferAppender
00094     );
00095 
00097   Teuchos::RCP<const InterpolationBufferAppenderBase<Scalar> > getInterpolationBufferAppender();
00098 
00100   Teuchos::RCP<InterpolationBufferAppenderBase<Scalar> > unSetInterpolationBufferAppender();
00101 
00103   void setStepper(
00104     const Teuchos::RCP<StepperBase<Scalar> > &stepper,
00105     const Scalar &finalTime,
00106     const bool landOnFinalTime
00107     );
00108 
00110   Teuchos::RCP<const Rythmos::StepperBase<Scalar> > getStepper() const;
00111 
00113   Teuchos::RCP<StepperBase<Scalar> > unSetStepper();
00114   
00116   void setBreakPointInformer(Teuchos::RCP<BreakPointInformerBase<Scalar> >& breakPointInformer);
00117 
00119   Teuchos::RCP<const BreakPointInformerBase<Scalar> > getBreakPointInformer();
00120 
00122   Teuchos::RCP<BreakPointInformerBase<Scalar> > unSetBreakPointInformer();
00123   
00125   void setObserver(Teuchos::RCP<IntegrationObserverBase<Scalar> >& observer);
00126 
00128   Teuchos::RCP<const IntegrationObserverBase<Scalar> > getObserver();
00129 
00131   Teuchos::RCP<IntegrationObserverBase<Scalar> > unSetObserver();
00132   
00134   
00135   void getFwdPoints(
00136     const Array<Scalar>& time_vec,
00137     Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >* x_vec,
00138     Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec,
00139     Array<ScalarMag>* accuracy_vec
00140     );
00141 
00143   TimeRange<Scalar> getFwdTimeRange() const;
00144   
00145 
00148     
00150   void addPoints(
00151     const Array<Scalar>& time_vec,
00152     const Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >& x_vec,
00153     const Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >& xdot_vec
00154     );
00155 
00157   void getPoints(
00158     const Array<Scalar>& time_vec,
00159     Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >* x_vec,
00160     Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec,
00161     Array<ScalarMag>* accuracy_vec
00162     ) const;
00163 
00165   TimeRange<Scalar> getTimeRange() const;
00166 
00168   void getNodes(Array<Scalar>* time_vec) const;
00169 
00171   int getOrder() const;
00172 
00174 
00176   virtual void removeNodes(Array<Scalar>& time_vec);
00177 
00180 
00182   void describe(
00183     Teuchos::FancyOStream &out,
00184     const Teuchos::EVerbosityLevel verbLevel
00185     ) const;
00186 
00188 
00191 
00193   void setParameterList(Teuchos::RCP<Teuchos::ParameterList> const& paramList);
00194 
00196   Teuchos::RCP<Teuchos::ParameterList> getParameterList();
00197 
00199   Teuchos::RCP<Teuchos::ParameterList> unsetParameterList();
00200 
00202   Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const;
00203 
00205 
00206 private:
00207 
00208 
00209   // Interpolation Buffer used to store past results
00210   Teuchos::RCP<InterpolationBufferBase<Scalar> > trailingInterpBuffer_;
00211 
00212   // Interpolation Buffer Appender strategy object for transfering data from
00213   // the stepper to the traling interpolation buffer
00214   Teuchos::RCP<InterpolationBufferAppenderBase<Scalar> > interpolationBufferAppender_;
00215 
00216   // Stepper used to fill interpolation buffer.
00217   Teuchos::RCP<StepperBase<Scalar> > stepper_;
00218 
00219   // Stepper Observer 
00220   Teuchos::RCP<IntegrationObserverBase<Scalar> > observer_;
00221 
00222   // ParameterList to control behavior
00223   Teuchos::RCP<Teuchos::ParameterList> paramList_;
00224 
00225   // BreakPoint informer strategy object
00226   // If this is null, then it has no effect
00227   Teuchos::RCP<BreakPointInformerBase<Scalar> > breakPointInformer_;
00228 
00229   // Take variable steps or not
00230   bool takeVariableSteps_;
00231 
00232   // The fixed step size to take
00233   Scalar fixed_dt_;
00234 
00235   // The final time to integrate to
00236   Scalar finalTime_;
00237 
00238   // Flag for whether the integrator is initialized or not.
00239   bool isInitialized_;
00240 
00241   // Private member functions:
00242 
00244   void initialize_();
00245 
00246 };
00247 
00248 
00249 // ////////////////////////////
00250 // Defintions
00251 
00252 
00253 // Constructors, Initializers, Misc
00254 
00255 
00256 template<class Scalar>
00257 IntegratorDefault<Scalar>::IntegratorDefault()
00258   : takeVariableSteps_(true), fixed_dt_(-1.0), isInitialized_(false)
00259 {}
00260 
00261 
00262 template<class Scalar>
00263 IntegratorDefault<Scalar>::IntegratorDefault(
00264   const Teuchos::RCP<StepperBase<Scalar> > &stepper,
00265   const Teuchos::RCP<InterpolationBufferBase<Scalar> > &trailingInterpBuffer,
00266   const Teuchos::RCP<Teuchos::ParameterList> &paramList 
00267   )
00268   : takeVariableSteps_(true), fixed_dt_(-1.0), isInitialized_(false)
00269 {
00270   using Teuchos::as;
00271   if (!is_null(paramList))
00272     setParameterList(paramList);
00273   const Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
00274   const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00275   Teuchos::OSTab ostab(out,1,"IntegratorDefault::constructor");
00276   *out << "Initializing IntegratorDefault" << std::endl;
00277   if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH)) {
00278     *out << "Calling setStepper..." << std::endl;
00279   }
00280   setStepper(stepper,0.0,BREAK_POINT_TYPE_HARD);
00281   if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH)) {
00282     *out << "Calling setInterpolationBuffer..." << std::endl;
00283   }
00284   setTrailingInterpolationBuffer(trailingInterpBuffer);
00285   initialize_();
00286 }
00287 
00288 template<class Scalar>
00289 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > IntegratorDefault<Scalar>::get_x_space() const
00290 {
00291   if (!isInitialized_) { 
00292     Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > space;
00293     return(space);
00294   }
00295   return(stepper_->get_x_space());
00296 }
00297 
00298 template<class Scalar>
00299 void IntegratorDefault<Scalar>::initialize_()
00300 {
00301   TEST_FOR_EXCEPTION(
00302       stepper_ == Teuchos::null, std::logic_error,
00303       "Error, Attempting to intialize integrator without setting a stepper!\n"
00304       );
00305   TEST_FOR_EXCEPTION(
00306       trailingInterpBuffer_ == Teuchos::null, std::logic_error,
00307       "Error, Attempting to intialize integrator without setting an interpolation buffer!\n"
00308       );
00309   if (paramList_ == Teuchos::null) {
00310     Teuchos::RCP<Teuchos::ParameterList> params;
00311     setParameterList(params);
00312   }
00313   if (interpolationBufferAppender_ == Teuchos::null) {
00314     Teuchos::RCP<InterpolationBufferAppenderDefault<Scalar> > 
00315       defaultInterpolationBufferAppender = Teuchos::rcp(new InterpolationBufferAppenderDefault<Scalar>);
00316     interpolationBufferAppender_ = defaultInterpolationBufferAppender;
00317   }
00318   isInitialized_ = true;
00319 }
00320 
00321 template<class Scalar>
00322 void IntegratorDefault<Scalar>::setTrailingInterpolationBuffer(
00323   const Teuchos::RCP<InterpolationBufferBase<Scalar> > &trailingInterpBuffer
00324   )
00325 {
00326 #ifdef TEUCHOS_DEBUG
00327   // Check preconditions:
00328   TEST_FOR_EXCEPTION(trailingInterpBuffer == Teuchos::null,std::logic_error,"Error, trailingInterpBuffer == Teuchos::null!\n");
00329   if (stepper_ != Teuchos::null) {
00330     TEST_FOR_EXCEPTION(
00331         trailingInterpBuffer->getTimeRange().upper() != stepper_->getTimeRange().lower(),
00332         std::logic_error,
00333         "Error, specified interpolation buffer's upper time range = " <<
00334           trailingInterpBuffer->getTimeRange().upper() << 
00335           " does not match internal stepper's lower time range = " <<
00336           stepper_->getTimeRange().lower() << "!\n"
00337         );
00338   }
00339 #endif // TEUCHOS_DEBUG
00340   using Teuchos::as;
00341   Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
00342   const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00343   trailingInterpBuffer_ = trailingInterpBuffer;
00344   Teuchos::OSTab ostab(out,1,"IntegratorDefault::setInterpolationBuffer");
00345   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00346     *out << "trailingInterpBuffer_ = " << trailingInterpBuffer_->description() << std::endl;
00347   }
00348   if (stepper_ != Teuchos::null) {
00349     initialize_();
00350   }
00351 }
00352 
00353 template<class Scalar>
00354 bool IntegratorDefault<Scalar>::acceptsTrailingInterpolationBuffer() const
00355 {
00356   return(true);
00357 }
00358 
00359 template<class Scalar>
00360 void IntegratorDefault<Scalar>::setStepper(
00361   const Teuchos::RCP<StepperBase<Scalar> > &stepper,
00362   const Scalar &finalTime,
00363   const bool landOnFinalTime // Not using this yet!
00364   )
00365 {
00366 #ifdef TEUCHOS_DEBUG
00367   // Check preconditions:
00368   TEST_FOR_EXCEPTION(stepper == Teuchos::null, std::logic_error, "Error, steppeer == Teuchos::null!\n");
00369   if (trailingInterpBuffer_ != Teuchos::null) {
00370     TEST_FOR_EXCEPTION(
00371         trailingInterpBuffer_->getTimeRange().upper() != stepper->getTimeRange().lower(),
00372         std::logic_error,
00373         "Error, specified stepper's lower time range = " <<
00374           stepper->getTimeRange().lower() <<
00375           " does not match internal trailing interpolation buffer's upper time range = " <<
00376           trailingInterpBuffer_->getTimeRange().upper() << "!\n"
00377         );
00378   }
00379 #endif // TEUCHOS_DEBUG
00380   using Teuchos::as;
00381   Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
00382   const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00383   stepper_ = stepper;
00384   Teuchos::OSTab ostab(out,1,"IntegratorDefault::setStepper");
00385   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00386     *out << "stepper_ = " << stepper_->description() << std::endl;
00387   }
00388   if (trailingInterpBuffer_ != Teuchos::null) {
00389     initialize_();
00390   }
00391 }
00392 
00393 template<class Scalar>
00394 Teuchos::RCP<StepperBase<Scalar> > IntegratorDefault<Scalar>::unSetStepper()
00395 {
00396   Teuchos::RCP<StepperBase<Scalar> > stepper_old = stepper_;
00397   stepper_ = Teuchos::null;
00398   isInitialized_ = false;
00399   return(stepper_old);
00400 }
00401 
00402 template<class Scalar>
00403 void IntegratorDefault<Scalar>::setBreakPointInformer(
00404     Teuchos::RCP<BreakPointInformerBase<Scalar> >& breakPointInformer
00405     )
00406 {
00407   breakPointInformer_ = breakPointInformer;
00408 }
00409 
00410 
00411 template<class Scalar>
00412 Teuchos::RCP<const BreakPointInformerBase<Scalar> > IntegratorDefault<Scalar>::getBreakPointInformer()
00413 {
00414   return(breakPointInformer_);
00415 }
00416 
00417 
00418 template<class Scalar>
00419 Teuchos::RCP<BreakPointInformerBase<Scalar> > IntegratorDefault<Scalar>::unSetBreakPointInformer()
00420 {
00421   Teuchos::RCP<BreakPointInformerBase<Scalar> > breakPointInformer_old = breakPointInformer_;
00422   breakPointInformer_ = Teuchos::null;
00423   return(breakPointInformer_old);
00424 }
00425 
00426 template<class Scalar>
00427 Teuchos::RCP<const Rythmos::StepperBase<Scalar> > IntegratorDefault<Scalar>::getStepper() const
00428 {
00429   return(stepper_);
00430 }
00431 
00432 template<class Scalar>
00433 Teuchos::RCP<const InterpolationBufferBase<Scalar> > IntegratorDefault<Scalar>::getInterpolationBuffer()
00434 {
00435   return(trailingInterpBuffer_);
00436 }
00437 
00438 template<class Scalar>
00439 Teuchos::RCP<InterpolationBufferBase<Scalar> > IntegratorDefault<Scalar>::unSetInterpolationBuffer()
00440 {
00441   Teuchos::RCP<InterpolationBufferBase<Scalar> > interpolationBuffer_old = trailingInterpBuffer_;
00442   trailingInterpBuffer_ = Teuchos::null;
00443   isInitialized_ = false;
00444   return(interpolationBuffer_old);
00445 }
00446 
00447 template<class Scalar>
00448 void IntegratorDefault<Scalar>::setInterpolationBufferAppender(
00449     const Teuchos::RCP<InterpolationBufferAppenderBase<Scalar> > &interpolationBufferAppender
00450     )
00451 {
00452   using Teuchos::as;
00453   Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
00454   const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00455   interpolationBufferAppender_ = interpolationBufferAppender;
00456   Teuchos::OSTab ostab(out,1,"IntegratorDefault::setiInterpolationBufferAppender");
00457   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00458     *out << "interpolationBufferAppender_ = " << interpolationBufferAppender_->description() << std::endl;
00459   }
00460 }
00461 
00462 template<class Scalar>
00463 Teuchos::RCP<const InterpolationBufferAppenderBase<Scalar> > IntegratorDefault<Scalar>::getInterpolationBufferAppender()
00464 {
00465   return(interpolationBufferAppender_);
00466 }
00467 
00468 template<class Scalar>
00469 Teuchos::RCP<InterpolationBufferAppenderBase<Scalar> > IntegratorDefault<Scalar>::unSetInterpolationBufferAppender()
00470 {
00471   Teuchos::RCP<InterpolationBufferAppenderBase<Scalar> > interpolationBufferAppender_old = interpolationBufferAppender_;
00472   interpolationBufferAppender_ = Teuchos::null;
00473   isInitialized_ = false;
00474   return(interpolationBufferAppender_old);
00475 }
00476 
00477 template<class Scalar>
00478 void IntegratorDefault<Scalar>::setObserver(Teuchos::RCP<IntegrationObserverBase<Scalar> >& observer)
00479 {
00480   observer_ = observer;
00481 }
00482 
00483 template<class Scalar>
00484 Teuchos::RCP<const IntegrationObserverBase<Scalar> > IntegratorDefault<Scalar>::getObserver()
00485 {
00486   return(observer_);
00487 }
00488 
00489 template<class Scalar>
00490 Teuchos::RCP<IntegrationObserverBase<Scalar> > IntegratorDefault<Scalar>::unSetObserver()
00491 {
00492   Teuchos::RCP<IntegrationObserverBase<Scalar> > observer_old = observer_;
00493   observer_ = Teuchos::null;
00494   return(observer_old);
00495 }
00496   
00497 
00498 // Overridden from InterpolationBufferBase
00499 
00500 
00501 template<class Scalar>
00502 void IntegratorDefault<Scalar>::addPoints(
00503   const Array<Scalar>& time_vec
00504   ,const Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >& x_vec
00505   ,const Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >& xdot_vec
00506   ) 
00507 {
00508   TEST_FOR_EXCEPTION(
00509       true, std::logic_error,
00510       "Error, addPoints is not defined for IntegratorDefault.\n"
00511       );
00512 }
00513 
00514 
00515 template<class Scalar>
00516 void IntegratorDefault<Scalar>::getPoints(
00517   const Array<Scalar>& time_vec,
00518   Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >* x_vec_in,
00519   Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec_in,
00520   Array<ScalarMag>* accuracy_vec_in
00521   ) const
00522 {
00523   TEST_FOR_EXCEPTION(!isInitialized_,std::logic_error,"Error, attemping to call getPoints (const) before initialization.\n");
00524 #ifdef TEUCHOS_DEBUG
00525   // check preconditions:
00526   for (int i=0; i<Teuchos::as<int>(time_vec.size()) ; ++i) {
00527     TEST_FOR_EXCEPTION(
00528         ~(trailingInterpBuffer_->getTimeRange().isInRange(time_vec[i]) || stepper_->getTimeRange().isInRange(time_vec[i])),
00529         std::logic_error,
00530         "Error, time_vec[" << i << "] is not in TimeRange of trailing interpolation buffer = [" <<
00531           trailingInterpBuffer_->getTimeRange().lower() << "," <<
00532           trailingInterpBuffer_->getTimeRange().upper() << "] nor in TimeRange of stepper = [" << 
00533           stepper_->getTimeRange().lower() << "," << stepper_->getTimeRange().upper() << "]!\n"
00534         );
00535   }
00536 #endif // TEUCHOS_DEBUG
00537 
00538   using Teuchos::as;
00539   typedef Teuchos::ScalarTraits<Scalar> ST;
00540 
00541   Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
00542   const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00543   Teuchos::OSTab ostab(out,1,"IntegratorDefault::getPoints");
00544 
00545   // Dump the input
00546 
00547   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00548     *out << "time_vec = " << std::endl;
00549     for (int i=0 ; i<Teuchos::as<int>(time_vec.size()) ; ++i) {
00550       *out << "time_vec[" << i << "] = " << time_vec[i] << std::endl;
00551     }
00552     if (x_vec_in == NULL) {
00553       *out << "x_vec_in = NULL" << std::endl;
00554     }
00555     else if (x_vec_in->size() == 0) {
00556       *out << "x_vec_in = ptr to empty vector" << std::endl;
00557     }
00558     else {
00559       *out << "x_vec_in = " << std::endl;
00560       for (int i=0 ; i<Teuchos::as<int>(x_vec_in->size()) ; ++i) {
00561         *out << "x_vec[" << i << "] = " << std::endl;
00562         (*x_vec_in)[i]->describe(*out,Teuchos::VERB_EXTREME);
00563       }
00564     }
00565     if (xdot_vec_in == NULL) {
00566       *out << "xdot_vec_in = NULL" << std::endl;
00567     }
00568     else if (xdot_vec_in->size() == 0) {
00569       *out << "xdot_vec_in = ptr to empty vector" << std::endl;
00570     }
00571     else {
00572       *out << "xdot_vec = " << std::endl;
00573       for (int i=0 ; i<Teuchos::as<int>(xdot_vec_in->size()) ; ++i) {
00574         if ((*xdot_vec_in)[i] == Teuchos::null) {
00575           *out << "xdot_vec[" << i << "] = Teuchos::null" << std::endl;
00576         } else {
00577           *out << "xdot_vec[" << i << "] = " << std::endl;
00578           (*xdot_vec_in)[i]->describe(*out,Teuchos::VERB_EXTREME);
00579         }
00580       }
00581     }
00582     if (accuracy_vec_in == NULL) {
00583       *out << "accuracy_vec_in = NULL" << std::endl;
00584     }
00585     else if (accuracy_vec_in->size() == 0) {
00586       *out << "accuracy_vec_in = ptr to empty vector" << std::endl;
00587     }
00588     else {
00589       *out << "accuracy_vec = " << std::endl;
00590       for (int i=0 ; i<Teuchos::as<int>(accuracy_vec_in->size()) ; ++i) {
00591         *out << "accuracy_vec[" << i << "] = " << (*accuracy_vec_in)[i] << std::endl;
00592       }
00593     }
00594   }
00595 
00596   Array<Scalar> time_vec_local = time_vec;
00597 
00598   // First we get points from the stepper_.
00599   Array<Scalar> stepper_time_vec;
00600   Array<RCP<const Thyra::VectorBase<Scalar> > > stepper_x_vec, stepper_xdot_vec;
00601   Array<ScalarMag> stepper_accuracy_vec;
00602 
00603   selectPointsInTimeRange<Scalar>(&stepper_time_vec,time_vec_local,stepper_->getTimeRange());
00604   removePointsInTimeRange<Scalar>(&time_vec_local,stepper_->getTimeRange());
00605   stepper_->getPoints(
00606       stepper_time_vec, 
00607       &stepper_x_vec, 
00608       &stepper_xdot_vec, 
00609       &stepper_accuracy_vec
00610       );
00611   typename DataStore<Scalar>::DataStoreVector_t stepper_dsv;
00612   vectorToDataStoreVector(
00613       stepper_time_vec, 
00614       stepper_x_vec, 
00615       stepper_xdot_vec, 
00616       stepper_accuracy_vec, 
00617       &stepper_dsv
00618       );
00619   
00620   // Next we get points from the trailingInterpBuffer_
00621   Array<Scalar> IB_time_vec;
00622   Array<RCP<const Thyra::VectorBase<Scalar> > > IB_x_vec, IB_xdot_vec;
00623   Array<ScalarMag> IB_accuracy_vec;
00624 
00625   selectPointsInTimeRange<Scalar>(&IB_time_vec,time_vec_local,trailingInterpBuffer_->getTimeRange());
00626   removePointsInTimeRange<Scalar>(&time_vec_local,trailingInterpBuffer_->getTimeRange());
00627   trailingInterpBuffer_->getPoints(
00628       IB_time_vec, 
00629       &IB_x_vec, 
00630       &IB_xdot_vec, 
00631       &IB_accuracy_vec
00632       );
00633   typename DataStore<Scalar>::DataStoreVector_t IB_dsv;
00634   vectorToDataStoreVector(
00635       IB_time_vec, 
00636       IB_x_vec, 
00637       IB_xdot_vec, 
00638       IB_accuracy_vec, 
00639       &IB_dsv
00640       );
00641 
00642   TEST_FOR_EXCEPTION(
00643       time_vec_local.size() != 0, std::logic_error,
00644       "Error, there are " << time_vec_local.size() << 
00645       " points in time_vec that were not found in either the stepper_ or the trailingInterpBuffer_!\n"
00646       );
00647 
00648   int IB_N = IB_dsv.size();
00649   for (int i=0 ; i < IB_N ; ++i) {
00650     stepper_dsv.push_back(IB_dsv[i]);
00651   }
00652 
00653   std::sort(stepper_dsv.begin(),stepper_dsv.end());
00654 
00655   Array<Scalar> time_vec_out;
00656   dataStoreVectorToVector(
00657       stepper_dsv,
00658       &time_vec_out,
00659       x_vec_in,
00660       xdot_vec_in,
00661       accuracy_vec_in
00662       );
00663 
00664   TEST_FOR_EXCEPTION(
00665      time_vec.size() != time_vec_out.size(),
00666      std::logic_error,
00667      "Error, number of output points = " << 
00668        time_vec_out.size()  << " != " << 
00669        time_vec.size() << " = number of time point requested!\n"
00670      ); 
00671 }
00672 
00673 template<class Scalar>
00674 TimeRange<Scalar> IntegratorDefault<Scalar>::getTimeRange() const
00675 {
00676   if (!isInitialized_) {
00677     TimeRange<Scalar> range; // return invalid time range
00678     return(range);
00679   }
00680   TimeRange<Scalar> timerange(trailingInterpBuffer_->getTimeRange().lower(),stepper_->getTimeRange().upper());
00681   return(timerange);
00682 }
00683 
00684 
00685 template<class Scalar>
00686 void IntegratorDefault<Scalar>::getNodes(
00687   Array<Scalar>* time_vec
00688   ) const
00689 {
00690   TEST_FOR_EXCEPTION(!isInitialized_,std::logic_error,"Error, attempting to call getNodes before initialized!\n");
00691   using Teuchos::as;
00692   Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
00693   const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00694   Teuchos::OSTab ostab(out,1,"IntegratorDefault::getNodes");
00695   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00696     *out << this->description() << std::endl;
00697   }
00698   return(trailingInterpBuffer_->getNodes(time_vec));
00699 }
00700 
00701 
00702 template<class Scalar>
00703 void IntegratorDefault<Scalar>::removeNodes(
00704   Array<Scalar>& time_vec
00705   ) 
00706 {
00707   TEST_FOR_EXCEPTION(!isInitialized_,std::logic_error,"Error, attempting to call removeNodes before initialized!\n");
00708   return(trailingInterpBuffer_->removeNodes(time_vec));
00709 }
00710 
00711 
00712 template<class Scalar>
00713 int IntegratorDefault<Scalar>::getOrder() const
00714 {
00715   TEST_FOR_EXCEPTION(!isInitialized_,std::logic_error,"Error, attempting to call getOrder before initialized!\n");
00716   return(trailingInterpBuffer_->getOrder());
00717 }
00718 
00719 
00720 // Overridden from Teuchos::Describable
00721 
00722 
00723 template<class Scalar>
00724 void IntegratorDefault<Scalar>::describe(
00725   Teuchos::FancyOStream &out,
00726   const Teuchos::EVerbosityLevel verbLevel
00727   ) const
00728 {
00729 
00730   using Teuchos::as;
00731   using Teuchos::describe;
00732 
00733   const bool printSomething =
00734     (
00735       as<int>(verbLevel) == as<int>(Teuchos::VERB_DEFAULT)
00736       ||
00737       as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW)
00738       );
00739 
00740   Teuchos::OSTab tab(out);
00741 
00742   if (printSomething) {
00743     out << this->description() << endl;
00744   }
00745 
00746   Teuchos::OSTab tab2(out);
00747 
00748   if (printSomething) {
00749     out << "trailing interpolation buffer = ";
00750     if (!is_null(trailingInterpBuffer_)) {
00751       out << describe(*trailingInterpBuffer_,verbLevel);
00752     } else {
00753       out << "NULL\n";
00754     }
00755     out << "stepper = ";
00756     if (!is_null(stepper_)) {
00757       out << describe(*stepper_,verbLevel);
00758     } else {
00759       out << "NULL\n";
00760     }
00761     out << "interpolationBufferAppender = ";
00762     if (!is_null(interpolationBufferAppender_)) {
00763       out << describe(*interpolationBufferAppender_,verbLevel);
00764     } else {
00765       out << "NULL\n";
00766     }
00767   }
00768   out << "paramList = ";
00769   if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_MEDIUM)) {
00770     if (!is_null(paramList_)) {
00771       out << paramList_->print(out);
00772     } else {
00773       out << "NULL\n";
00774     }
00775   }
00776   
00777 }
00778 
00779 
00780 // Overridden from Teuchos::ParameterListAcceptor
00781 
00782 
00783 template <class Scalar>
00784 void IntegratorDefault<Scalar>::setParameterList(
00785   Teuchos::RCP<Teuchos::ParameterList> const& paramList
00786   )
00787 {
00788   using Teuchos::as;
00789   using Teuchos::get;
00790   typedef Teuchos::ScalarTraits<Scalar> ST;
00791 
00792   TEST_FOR_EXCEPT(is_null(paramList));
00793 
00794   Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
00795   const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00796   Teuchos::OSTab ostab(out,1,"IntegratorDefault::setParameterList");
00797 
00798   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00799     *out << "paramList = " << paramList->print(*out) << std::endl;
00800   }
00801 
00802   paramList->validateParametersAndSetDefaults(*getValidParameters());
00803   paramList_ = paramList;
00804 
00805   takeVariableSteps_ = get<bool>(*paramList_,"Take Variable Steps");
00806   fixed_dt_ = get<Scalar>(*paramList_,"fixed_dt");
00807   finalTime_ = get<Scalar>(*paramList_,"Final Time");
00808 
00809   TEST_FOR_EXCEPTION(
00810     !takeVariableSteps_ && fixed_dt_ <= ST::zero(),
00811     Teuchos::Exceptions::InvalidParameterValue,
00812     "Error, if we are taking fixed steps then \"fixed_dt\" must be set to"
00813     " a positive number!"
00814     );
00815 
00816   Teuchos::readVerboseObjectSublist(&*paramList_,this);
00817 
00818 #ifdef TEUCHOS_DEBUG
00819   paramList->validateParameters(*getValidParameters());
00820 #endif
00821 }
00822 
00823 
00824 template <class Scalar>
00825 Teuchos::RCP<Teuchos::ParameterList>
00826 IntegratorDefault<Scalar>::getParameterList()
00827 {
00828   return(paramList_);
00829 }
00830 
00831 
00832 template <class Scalar>
00833 Teuchos::RCP<Teuchos::ParameterList>
00834 IntegratorDefault<Scalar>::unsetParameterList()
00835 {
00836   Teuchos::RCP<Teuchos::ParameterList> temp_param_list = paramList_;
00837   paramList_ = Teuchos::null;
00838   return(temp_param_list);
00839 }
00840 
00841 
00842 template <class Scalar>
00843 Teuchos::RCP<const Teuchos::ParameterList>
00844 IntegratorDefault<Scalar>::getValidParameters() const
00845 {
00846   using Teuchos::RCP;
00847   using Teuchos::ParameterList;
00848   static RCP<const ParameterList> validPL;
00849   if (is_null(validPL)) {
00850     RCP<ParameterList> pl = Teuchos::parameterList();
00851     pl->set(
00852       "Take Variable Steps", bool(true),
00853       "If set to true, then the stepper will take variable steps.\n"
00854       "If set to false, then fixed stepps will be taken of size\n"
00855       "\"fixed_dt\".  Warning, you must set \"fixed_dt\" in this case!"
00856       );
00857     pl->set(
00858       "fixed_dt", Scalar(-1.0),
00859       "If fixed steps are being taken, then this is the step size."
00860       );
00861     pl->set(
00862       "Final Time", Scalar(0.0),
00863       "This specifies the final time to which the integrator will integrate.\n"
00864       );
00865     Teuchos::setupVerboseObjectSublist(&*pl);
00866     validPL = pl;
00867   }
00868   return (validPL);
00869 }
00870 
00871 /* 08/13/07 tscoffe:  To simplify getPoints and allow the time points to be
00872  * unsorted & to come from current values and forward values, I can write a
00873  * helper function to do all of this.
00874  */
00875 template <class Scalar>
00876 void IntegratorDefault<Scalar>::getFwdPoints(
00877     const Array<Scalar>& time_vec,
00878     Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >* x_vec,
00879     Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec,
00880     Array<ScalarMag>* accuracy_vec
00881     )
00882 {
00883   typedef Teuchos::ScalarTraits<Scalar> ST;
00884   if (!isInitialized_) {
00885     initialize_();
00886   }
00887 #ifdef TEUCHOS_DEBUG
00888   // Check preconditions:
00889   for (int i=0 ; i<Teuchos::as<int>(time_vec.size()) ; ++i) { 
00890     TEST_FOR_EXCEPTION(
00891         isInRange_oc(this->getFwdTimeRange(),time_vec[i]),
00892         std::logic_error,
00893         "Error, time_vec[" << i << "] = " << time_vec[i] << 
00894           " is not in this->getFwdTimeRange:  (" << 
00895           this->getFwdTimeRange().lower() << "," <<
00896           this->getFwdTimeRange().upper() << "]!\n"
00897         );
00898   }
00899   // Check that time_vec is sorted:
00900   assertTimePointsAreSorted(time_vec);
00901 #endif // TEUCHOS_DEBUG
00902 
00903   Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
00904   const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00905   Teuchos::OSTab ostab(out,1,"IntegratorDefault::getFwdPoints");
00906 
00907   Array<Scalar> time_vec_local = time_vec;
00908   Array<Scalar> temp_time_vec;
00909   
00910   // First we use this->getPoints to recover points in the trailingInterpBuffer_ and the stepper_.
00911   selectPointsInTimeRange<Scalar>(&temp_time_vec,time_vec_local,this->getTimeRange());
00912   removePointsInTimeRange<Scalar>(&time_vec_local,this->getTimeRange());
00913   this->getPoints(
00914       temp_time_vec, 
00915       x_vec, 
00916       xdot_vec, 
00917       accuracy_vec
00918       );
00919    
00920   // Then we integrate forward for new points
00921   
00922   // Assumptions:
00923   // stepper_ has initial condition and is ready to take steps.
00924   // trailingInterpBuffer_ is initialized and is ready for points to be imported to it from stepper_.
00925   //
00926   // 08/13/07 tscoffe:  Question:  We're getting RCPs from the Stepper here and
00927   // stuffing them into the trailing interpolation buffer, what if the stepper
00928   // decides to re-use the vector sitting bethind the RCP?  It will corrupt the
00929   // data in the trailing interpolation buffer.  Oof.
00930   // This is okay for ImplicitBDFStepper, but we need to define this behavior
00931   // at the StepperBase level.  TODO
00932 
00933   RCP<Array<RCP<const Thyra::VectorBase<Scalar> > > > local_x_vec, local_xdot_vec;
00934   RCP<Array<ScalarMag> > local_accuracy_vec;
00935 
00936   if (x_vec == 0) { 
00937     local_x_vec = Teuchos::null; 
00938   } else {
00939     local_x_vec = rcp(new Array<RCP<const Thyra::VectorBase<Scalar> > >());
00940   }
00941   if (xdot_vec == 0) { 
00942     local_xdot_vec = Teuchos::null; 
00943   } else {
00944     local_xdot_vec = rcp(new Array<RCP<const Thyra::VectorBase<Scalar> > >());
00945   }
00946   if (accuracy_vec == 0) { 
00947     local_accuracy_vec = Teuchos::null; 
00948   } else {
00949     local_accuracy_vec = rcp(new Array<ScalarMag>());
00950   }
00951 
00952   for (int i=0 ; i<Teuchos::as<int>(time_vec_local.size()) ; ++i) {
00953     // 08/13/07 tscoffe:  we're going to grab points one at a time from the stepper.
00954 
00955     while (!(stepper_->getTimeRange().isInRange(time_vec_local[i]))) {
00956       if (takeVariableSteps_) {
00957         Scalar step_taken = stepper_->takeStep(ST::zero(),STEP_TYPE_VARIABLE);
00958         if (Teuchos::as<int>(verbLevel) >= Teuchos::as<int>(Teuchos::VERB_HIGH)) {
00959           *out << "step_taken = " << step_taken << std::endl;
00960         }
00961       } 
00962       else {
00963         Scalar step_taken = stepper_->takeStep(fixed_dt_,STEP_TYPE_FIXED);
00964         if (Teuchos::as<int>(verbLevel) >= Teuchos::as<int>(Teuchos::VERB_HIGH)) {
00965           *out << "step_taken = " << step_taken << std::endl;
00966         }
00967       }
00968       interpolationBufferAppender_->import(&*trailingInterpBuffer_,*stepper_,stepper_->getTimeRange());
00969       // After we've stepped forward & transfered the data into the trailing interpolation buffer, call the observer.
00970       if (observer_ != Teuchos::null)
00971       {
00972         TEST_FOR_EXCEPT(true); // ToDo: Call the new interface!
00973         //observer_->notify(*stepper_);
00974       } 
00975     }
00976     Array<Scalar> temp_time_vec;
00977     temp_time_vec.push_back(time_vec_local[i]);
00978 
00979     stepper_->getPoints(temp_time_vec,&*local_x_vec,&*local_xdot_vec,&*local_accuracy_vec);
00980 
00981     if (x_vec) {
00982       x_vec->push_back((*local_x_vec)[0]);
00983     }
00984     if (xdot_vec) {
00985       xdot_vec->push_back((*local_xdot_vec)[0]);
00986     }
00987     if (accuracy_vec) {
00988       accuracy_vec->push_back((*local_accuracy_vec)[0]);
00989     }
00990   }
00991 
00992 }
00993 
00994 template <class Scalar>
00995 TimeRange<Scalar> IntegratorDefault<Scalar>::getFwdTimeRange() const
00996 {
00997   if (!isInitialized_) {
00998     TimeRange<Scalar> range; // invalid time range.
00999     return(range);
01000   }
01001   TimeRange<Scalar> timerange(trailingInterpBuffer_->getTimeRange().lower(),finalTime_);
01002   return(timerange);
01003 }
01004 
01005 } // namespace Rythmos
01006 
01007 
01008 #endif //Rythmos_INTEGRATOR_DEFAULT_H

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