Rythmos_SimpleIntegrator.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_SIMPLE_INTEGRATOR_H
00030 #define Rythmos_SIMPLE_INTEGRATOR_H
00031 
00032 
00033 #include "Rythmos_IntegratorBase.hpp"
00034 #include "Rythmos_InterpolationBufferHelpers.hpp"
00035 #include "Rythmos_IntegrationControlStrategyBase.hpp"
00036 #include "Rythmos_IntegrationObserverBase.hpp"
00037 #include "Rythmos_StepperHelpers.hpp"
00038 #include "Teuchos_ParameterListAcceptorDefaultBase.hpp"
00039 #include "Teuchos_VerboseObjectParameterListHelpers.hpp"
00040 #include "Teuchos_Assert.hpp"
00041 #include "Teuchos_as.hpp"
00042 
00043 
00044 namespace Rythmos {
00045 
00046 
00050 template<class Scalar> 
00051 class SimpleIntegrator
00052   : virtual public IntegratorBase<Scalar>,
00053     virtual public Teuchos::ParameterListAcceptorDefaultBase
00054 {
00055 public:
00056   
00058   typedef typename ScalarTraits<Scalar>::magnitudeType ScalarMag;
00059 
00062   
00064   SimpleIntegrator();
00065 
00067   void setIntegrationControlStrategy(
00068     const RCP<IntegrationControlStrategyBase<Scalar> > &integrationControlStrategy
00069     );
00070 
00072   void setIntegrationObserver(
00073     const RCP<IntegrationObserverBase<Scalar> > &integrationObserver
00074     );
00075 
00077 
00080 
00082   void setParameterList(RCP<ParameterList> const& paramList);
00083 
00085   RCP<const ParameterList> getValidParameters() const;
00086 
00088 
00091 
00093   RCP<IntegratorBase<Scalar> > cloneIntegrator() const;
00094   
00096   void setStepper(
00097     const RCP<StepperBase<Scalar> > &stepper,
00098     const Scalar &finalTime,
00099     const bool landOnFinalTime
00100     );
00101 
00103   RCP<StepperBase<Scalar> > unSetStepper();
00104 
00106   Teuchos::RCP<const StepperBase<Scalar> > getStepper() const;
00107 
00109   void getFwdPoints(
00110     const Array<Scalar>& time_vec,
00111     Array<RCP<const Thyra::VectorBase<Scalar> > >* x_vec,
00112     Array<RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec,
00113     Array<ScalarMag>* accuracy_vec
00114     );
00115 
00117   TimeRange<Scalar> getFwdTimeRange() const;
00118 
00120 
00123 
00125   RCP<const Thyra::VectorSpaceBase<Scalar> > get_x_space() const;
00126     
00128   void addPoints(
00129     const Array<Scalar>& time_vec,
00130     const Array<RCP<const Thyra::VectorBase<Scalar> > >& x_vec,
00131     const Array<RCP<const Thyra::VectorBase<Scalar> > >& xdot_vec
00132     );
00133 
00135   void getPoints(
00136     const Array<Scalar>& time_vec,
00137     Array<RCP<const Thyra::VectorBase<Scalar> > >* x_vec,
00138     Array<RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec,
00139     Array<ScalarMag>* accuracy_vec
00140     ) const;
00141 
00143   TimeRange<Scalar> getTimeRange() const;
00144 
00146   void getNodes(Array<Scalar>* time_vec) const;
00147 
00149   void removeNodes(Array<Scalar>& time_vec);
00150 
00152   int getOrder() const;
00153 
00155 
00156 private:
00157 
00158   // ////////////////////////
00159   // Private data members
00160 
00161   RCP<IntegrationControlStrategyBase<Scalar> > integrationControlStrategy_;
00162   RCP<IntegrationObserverBase<Scalar> > integrationObserver_;
00163   RCP<StepperBase<Scalar> > stepper_;
00164   TimeRange<Scalar> integrationTimeDomain_;
00165   bool landOnFinalTime_;
00166 
00167   int maxNumTimeSteps_;
00168 
00169   int currTimeStepIndex_;
00170   StepControlInfo<Scalar> stepCtrlInfoLast_;
00171 
00172   static const std::string maxNumTimeSteps_name_;
00173   static const int maxNumTimeSteps_default_;
00174 
00175   // /////////////////////////
00176   // Private member functions
00177 
00178   bool advanceStepperToTime( const Scalar& t );
00179 
00180 };
00181 
00182 
00187 template<class Scalar> 
00188 RCP<SimpleIntegrator<Scalar> >
00189 simpleIntegrator()
00190 {
00191   RCP<SimpleIntegrator<Scalar> >
00192     integrator = Teuchos::rcp(new SimpleIntegrator<Scalar>());
00193   return integrator;
00194 }
00195 
00196 
00201 template<class Scalar> 
00202 RCP<SimpleIntegrator<Scalar> >
00203 controlledSimpleIntegrator(
00204   const RCP<IntegrationControlStrategyBase<Scalar> > &integrationControlStrategy
00205   )
00206 {
00207   RCP<SimpleIntegrator<Scalar> >
00208     integrator = Teuchos::rcp(new SimpleIntegrator<Scalar>());
00209   integrator->setIntegrationControlStrategy(integrationControlStrategy);
00210   return integrator;
00211 }
00212 
00213 
00218 template<class Scalar> 
00219 RCP<SimpleIntegrator<Scalar> >
00220 observedSimpleIntegrator(
00221   const RCP<IntegrationObserverBase<Scalar> > &integrationObserver
00222   )
00223 {
00224   RCP<SimpleIntegrator<Scalar> >
00225     integrator = Teuchos::rcp(new SimpleIntegrator<Scalar>());
00226   integrator->setIntegrationObserver(integrationObserver);
00227   return integrator;
00228 }
00229 
00230 
00235 template<class Scalar> 
00236 RCP<SimpleIntegrator<Scalar> >
00237 simpleIntegrator(
00238   const RCP<IntegrationControlStrategyBase<Scalar> > &integrationControlStrategy,
00239   const RCP<IntegrationObserverBase<Scalar> > &integrationObserver
00240   )
00241 {
00242   RCP<SimpleIntegrator<Scalar> >
00243     integrator = Teuchos::rcp(new SimpleIntegrator<Scalar>());
00244   integrator->setIntegrationControlStrategy(integrationControlStrategy);
00245   integrator->setIntegrationObserver(integrationObserver);
00246   return integrator;
00247 }
00248 
00249 
00250 // 2007/08/30: rabartl: Above, note that I had to name the nonmember
00251 // constructors taking an single RCP argument different names from each other
00252 // in order to get around the classic ambiguity problem with implicit
00253 // conversions of smart pointers.
00254 
00255 
00256 // ////////////////////////////
00257 // Defintions
00258 
00259 
00260 // Static data members
00261 
00262 
00263 template<class Scalar>
00264 const std::string
00265 SimpleIntegrator<Scalar>::maxNumTimeSteps_name_ = "Max Number Time Steps";
00266 
00267 template<class Scalar>
00268 const int
00269 SimpleIntegrator<Scalar>::maxNumTimeSteps_default_ = 10000;
00270 
00271 
00272 // Constructors, Initializers, Misc
00273 
00274 
00275 template<class Scalar>
00276 SimpleIntegrator<Scalar>::SimpleIntegrator()
00277   :landOnFinalTime_(true),
00278    maxNumTimeSteps_(maxNumTimeSteps_default_),
00279    currTimeStepIndex_(-1)
00280 {}
00281 
00282 
00283 template<class Scalar>
00284 void SimpleIntegrator<Scalar>::setIntegrationControlStrategy(
00285   const RCP<IntegrationControlStrategyBase<Scalar> > &integrationControlStrategy
00286   )
00287 {
00288 #ifdef TEUCHOS_DEBUG
00289   TEST_FOR_EXCEPT(is_null(integrationControlStrategy));
00290 #endif
00291   integrationControlStrategy_ = integrationControlStrategy;
00292 }
00293 
00294 
00295 template<class Scalar>
00296 void SimpleIntegrator<Scalar>::setIntegrationObserver(
00297   const RCP<IntegrationObserverBase<Scalar> > &integrationObserver
00298   )
00299 {
00300 #ifdef TEUCHOS_DEBUG
00301   TEST_FOR_EXCEPT(is_null(integrationObserver));
00302 #endif
00303   integrationObserver_ = integrationObserver;
00304 }
00305 
00306 
00307 // Overridden from ParameterListAcceptor
00308 
00309 
00310 template<class Scalar> 
00311 void SimpleIntegrator<Scalar>::setParameterList(
00312   RCP<ParameterList> const& paramList
00313   )
00314 {
00315   TEST_FOR_EXCEPT(is_null(paramList));
00316   paramList->validateParameters(*getValidParameters());
00317   this->setMyParamList(paramList);
00318   maxNumTimeSteps_ = paramList->get(
00319     maxNumTimeSteps_name_, maxNumTimeSteps_default_);
00320   Teuchos::readVerboseObjectSublist(&*paramList,this);
00321 }
00322 
00323 
00324 template<class Scalar> 
00325 RCP<const ParameterList>
00326 SimpleIntegrator<Scalar>::getValidParameters() const
00327 {
00328   static RCP<const ParameterList> validPL;
00329   if (is_null(validPL) ) {
00330     RCP<ParameterList> pl = Teuchos::parameterList();
00331     pl->set(maxNumTimeSteps_name_, maxNumTimeSteps_default_,
00332       "Set the maximum number of integration time-steps allowed."
00333       );
00334     Teuchos::setupVerboseObjectSublist(&*pl);
00335     validPL = pl;
00336   }
00337   return validPL;
00338 }
00339 
00340 
00341 // Overridden from IntegratorBase
00342 
00343 
00344 template<class Scalar>
00345 RCP<IntegratorBase<Scalar> >
00346 SimpleIntegrator<Scalar>::cloneIntegrator() const
00347 {
00348   RCP<SimpleIntegrator<Scalar> >
00349     newIntegrator = Teuchos::rcp(new SimpleIntegrator<Scalar>());
00350   // Only copy control information, not the state of an existing integration!
00351   newIntegrator->stepper_ = Teuchos::null;
00352   const RCP<const ParameterList> paramList = this->getParameterList();
00353   if (!is_null(paramList))
00354     newIntegrator->setParameterList(Teuchos::parameterList(*paramList));
00355   if(!is_null(integrationControlStrategy_))
00356     newIntegrator->integrationControlStrategy_ =
00357       integrationControlStrategy_->cloneIntegrationControlStrategy().assert_not_null();
00358   if(!is_null(integrationObserver_))
00359     newIntegrator->integrationObserver_ =
00360       integrationObserver_->cloneIntegrationObserver().assert_not_null();
00361   return newIntegrator;
00362 }
00363 
00364 
00365 template<class Scalar> 
00366 void SimpleIntegrator<Scalar>::setStepper(
00367   const RCP<StepperBase<Scalar> > &stepper,
00368   const Scalar &finalTime,
00369   const bool landOnFinalTime
00370   )
00371 {
00372   typedef Teuchos::ScalarTraits<Scalar> ST;
00373   TEST_FOR_EXCEPT(is_null(stepper));
00374   TEST_FOR_EXCEPT( finalTime <= stepper->getTimeRange().lower() );
00375   TEUCHOS_ASSERT( stepper->getTimeRange().length() == ST::zero() );
00376   // 2007/07/25: rabartl: ToDo: Validate state of the stepper!
00377   stepper_ = stepper;
00378   integrationTimeDomain_ = timeRange(stepper_->getTimeRange().lower(), finalTime);
00379   landOnFinalTime_ = landOnFinalTime;
00380   currTimeStepIndex_ = 0;
00381   stepCtrlInfoLast_ = StepControlInfo<Scalar>();
00382   if (!is_null(integrationControlStrategy_))
00383     integrationControlStrategy_->resetIntegrationControlStrategy(
00384       integrationTimeDomain_
00385       );
00386   if (!is_null(integrationObserver_))
00387     integrationObserver_->resetIntegrationObserver(
00388       integrationTimeDomain_
00389       );
00390 }
00391 
00392 template<class Scalar>
00393 RCP<StepperBase<Scalar> > SimpleIntegrator<Scalar>::unSetStepper()
00394 {
00395   RCP<StepperBase<Scalar> > stepper_temp = stepper_;
00396   stepper_ = Teuchos::null;
00397   return(stepper_temp);
00398 }
00399 
00400 template<class Scalar>
00401 Teuchos::RCP<const StepperBase<Scalar> > SimpleIntegrator<Scalar>::getStepper() const
00402 {
00403   return(stepper_);
00404 }
00405 
00406 template<class Scalar>
00407 void SimpleIntegrator<Scalar>::getFwdPoints(
00408   const Array<Scalar>& time_vec,
00409   Array<RCP<const Thyra::VectorBase<Scalar> > >* x_vec,
00410   Array<RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec,
00411   Array<ScalarMag>* accuracy_vec
00412   )
00413 {
00414 
00415 #ifdef ENABLE_RYTHMOS_TIMERS
00416   TEUCHOS_FUNC_TIME_MONITOR("Rythmos:SimpleIntegrator::getFwdPoints");
00417 #endif
00418 
00419   using Teuchos::incrVerbLevel;
00420   using Teuchos::Describable;
00421   typedef Teuchos::ScalarTraits<Scalar> ST;
00422   typedef InterpolationBufferBase<Scalar> IBB;
00423   typedef Teuchos::VerboseObjectTempState<IBB> VOTSIBB;
00424 
00425   Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
00426   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00427   Teuchos::OSTab tab(out);
00428   VOTSIBB stepper_outputTempState(stepper_,out,incrVerbLevel(verbLevel,-1));
00429 
00430   if ( includesVerbLevel(verbLevel,Teuchos::VERB_LOW) )
00431     *out << "\nEntering " << this->Describable::description() << "::getFwdPoints(...) ...\n"
00432          << "\nStepper: " << Teuchos::describe(*stepper_,verbLevel);
00433 
00434   if ( includesVerbLevel(verbLevel,Teuchos::VERB_MEDIUM) )
00435     *out << "\nRequested time points: " << Teuchos::toString(time_vec) << "\n";
00436 
00437   //
00438   // 1) Initial setup
00439   //
00440 
00441   const int numTimePoints = time_vec.size();
00442 
00443   // Assert preconditions
00444   assertTimePointsAreSorted(time_vec);
00445   TEST_FOR_EXCEPT(accuracy_vec!=0); // ToDo: Remove accuracy_vec!
00446 
00447   // Resize the storage for the output arrays
00448   if (x_vec)
00449     x_vec->resize(numTimePoints);
00450   if (xdot_vec)
00451     xdot_vec->resize(numTimePoints);
00452 
00453   // This int records the next time point offset in time_vec[timePointIndex]
00454   // that needs to be handled.  This gets updated as the time points are
00455   // filled below.
00456   int nextTimePointIndex = 0;
00457   
00458   assertNoTimePointsBeforeCurrentTimeRange(*this,time_vec,nextTimePointIndex);
00459 
00460   //
00461   // 1) First, get all time points that fall within the current time range
00462   //
00463 
00464   {
00465 #ifdef ENABLE_RYTHMOS_TIMERS
00466     TEUCHOS_FUNC_TIME_MONITOR("Rythmos:SimpleIntegrator::getFwdPoints: getPoints");
00467 #endif
00468     getCurrentPoints(*stepper_,time_vec,x_vec,xdot_vec,&nextTimePointIndex);
00469   }
00470 
00471   //
00472   // 2) Advance the stepper to satisfy time points in time_vec that fall
00473   // before the current time.
00474   //
00475 
00476   while ( nextTimePointIndex < numTimePoints ) {
00477     
00478     // Use the time stepping algorithm to step up to or past the next
00479     // requested time but not so far as to step past the point entirely.
00480     const Scalar t = time_vec[nextTimePointIndex];
00481     bool advanceStepperToTimeSucceeded = false;
00482     {
00483 #ifdef ENABLE_RYTHMOS_TIMERS
00484       TEUCHOS_FUNC_TIME_MONITOR("Rythmos:SimpleIntegrator::getFwdPoints: advanceStepperToTime");
00485 #endif
00486       advanceStepperToTimeSucceeded= advanceStepperToTime(t);
00487     }
00488     TEST_FOR_EXCEPTION(
00489       !advanceStepperToTimeSucceeded, Exceptions::GetFwdPointsFailed,
00490       this->description() << "\n\n"
00491       "Error:  The integration failed to get to time " << t << " and only achieved\n"
00492       "getting to " << stepper_->getTimeRange().upper() << "!"
00493       );
00494     
00495     // Extract the next set of points (perhaps just one) from the stepper
00496     {
00497 #ifdef ENABLE_RYTHMOS_TIMERS
00498       TEUCHOS_FUNC_TIME_MONITOR("Rythmos:SimpleIntegrator::getFwdPoints: getPoints (fwd)");
00499 #endif
00500       getCurrentPoints(*stepper_,time_vec,x_vec,xdot_vec,&nextTimePointIndex);
00501     }
00502     
00503   }
00504 
00505   if ( includesVerbLevel(verbLevel,Teuchos::VERB_LOW) )
00506     *out << "\nLeaving " << this->Describable::description() << "::getFwdPoints(...) ...\n";
00507   
00508 }
00509 
00510 
00511 template<class Scalar> 
00512 TimeRange<Scalar>
00513 SimpleIntegrator<Scalar>::getFwdTimeRange() const
00514 {
00515   return timeRange(
00516     stepper_->getTimeRange().lower(),
00517     integrationTimeDomain_.upper()
00518     );
00519 }
00520 
00521 
00522 // Overridden from InterpolationBufferBase
00523 
00524 template<class Scalar>
00525 RCP<const Thyra::VectorSpaceBase<Scalar> > SimpleIntegrator<Scalar>::get_x_space() const
00526 {
00527   if (stepper_ == Teuchos::null) {
00528     RCP<const Thyra::VectorSpaceBase<Scalar> > space;
00529     return(space);
00530   } else {
00531     return(stepper_->get_x_space());
00532   }
00533 }
00534 
00535 template<class Scalar> 
00536 void SimpleIntegrator<Scalar>::addPoints(
00537   const Array<Scalar>& time_vec,
00538   const Array<RCP<const Thyra::VectorBase<Scalar> > >& x_vec,
00539   const Array<RCP<const Thyra::VectorBase<Scalar> > >& xdot_vec
00540   )
00541 {
00542   stepper_->addPoints(time_vec,x_vec,xdot_vec);
00543 }
00544 
00545 
00546 template<class Scalar> 
00547 void SimpleIntegrator<Scalar>::getPoints(
00548   const Array<Scalar>& time_vec,
00549   Array<RCP<const Thyra::VectorBase<Scalar> > >* x_vec,
00550   Array<RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec,
00551   Array<ScalarMag>* accuracy_vec
00552   ) const
00553 {
00554   stepper_->getPoints(time_vec,x_vec,xdot_vec,accuracy_vec);
00555 }
00556 
00557 
00558 template<class Scalar> 
00559 TimeRange<Scalar> SimpleIntegrator<Scalar>::getTimeRange() const
00560 {
00561   return stepper_->getTimeRange();
00562 }
00563 
00564 
00565 template<class Scalar> 
00566 void SimpleIntegrator<Scalar>::getNodes(Array<Scalar>* time_vec) const
00567 {
00568   stepper_->getNodes(time_vec);
00569 }
00570 
00571 
00572 template<class Scalar> 
00573 void SimpleIntegrator<Scalar>::removeNodes(Array<Scalar>& time_vec)
00574 {
00575   stepper_->removeNodes(time_vec);
00576 }
00577 
00578 
00579 template<class Scalar> 
00580 int SimpleIntegrator<Scalar>::getOrder() const
00581 {
00582   return stepper_->getOrder();
00583 }
00584 
00585 
00586 // private
00587 
00588 
00589 template<class Scalar> 
00590 bool SimpleIntegrator<Scalar>::advanceStepperToTime( const Scalar& advance_to_t )
00591 {
00592 
00593 #ifdef ENABLE_RYTHMOS_TIMERS
00594   TEUCHOS_FUNC_TIME_MONITOR("Rythmos:SimpleIntegrator::advanceStepperToTime");
00595 #endif
00596 
00597   using std::endl;
00598   typedef std::numeric_limits<Scalar> NL;
00599   using Teuchos::incrVerbLevel;
00600   using Teuchos::Describable;
00601   using Teuchos::OSTab;
00602   typedef Teuchos::ScalarTraits<Scalar> ST;
00603 
00604   Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
00605   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00606   Teuchos::OSTab tab(out);
00607 
00608   if (!is_null(integrationControlStrategy_)) {
00609     integrationControlStrategy_->setOStream(out);
00610     integrationControlStrategy_->setVerbLevel(incrVerbLevel(verbLevel,-1));
00611   }
00612 
00613   if (!is_null(integrationObserver_)) {
00614     integrationObserver_->setOStream(out);
00615     integrationObserver_->setVerbLevel(incrVerbLevel(verbLevel,-1));
00616   }
00617 
00618   if ( includesVerbLevel(verbLevel,Teuchos::VERB_LOW) )
00619     *out << "\nEntering " << this->Describable::description()
00620          << "::advanceStepperToTime("<<advance_to_t<<") ...\n";
00621 
00622   // Remember what timestep index we are on so we can report it later
00623   const int initCurrTimeStepIndex = currTimeStepIndex_;
00624 
00625   // Take steps until we the requested time is reached (or passed)
00626 
00627   TimeRange<Scalar> currStepperTimeRange = stepper_->getTimeRange();
00628 
00629   // Start by assume we can reach the time advance_to_t
00630   bool return_val = true;
00631   
00632   while ( !currStepperTimeRange.isInRange(advance_to_t) ) {
00633 
00634     // Halt immediatly if exceeded max iterations
00635     if (currTimeStepIndex_ >= maxNumTimeSteps_) {
00636       if ( includesVerbLevel(verbLevel,Teuchos::VERB_LOW) )
00637         *out
00638           << "\n***"
00639           << "\n*** NOTICE: currTimeStepIndex = "<<currTimeStepIndex_
00640           << " >= maxNumTimeSteps = "<<maxNumTimeSteps_<< ", halting time integration!"
00641           << "\n***\n";
00642       return_val = false;
00643       break; // Exit the loop immediately!
00644     }
00645 
00646     if ( includesVerbLevel(verbLevel,Teuchos::VERB_LOW) )
00647       *out << "\nTake step:  current_stepper_t = " << currStepperTimeRange.upper()
00648            << ", currTimeStepIndex = " << currTimeStepIndex_ << endl;
00649     Teuchos::OSTab tab(out);
00650 
00651     //
00652     // A) Reinitialize if a hard breakpoint was reached on the last time step
00653     //
00654 
00655     if (stepCtrlInfoLast_.limitedByBreakPoint) {
00656       if ( stepCtrlInfoLast_.breakPointType == BREAK_POINT_TYPE_HARD ) {
00657 #ifdef ENABLE_RYTHMOS_TIMERS
00658         TEUCHOS_FUNC_TIME_MONITOR("Rythmos:SimpleIntegrator::restart");
00659 #endif
00660         if ( includesVerbLevel(verbLevel,Teuchos::VERB_LOW) )
00661           *out << "\nAt a hard-breakpoint, restarting time integrator ...\n";
00662         restart(&*stepper_);
00663       }
00664       else  {
00665         if ( includesVerbLevel(verbLevel,Teuchos::VERB_LOW) )
00666           *out << "\nAt a soft-breakpoint, NOT restarting time integrator ...\n";
00667       }
00668     }
00669 
00670     //
00671     // B) Get the trial step control info
00672     //
00673 
00674     StepControlInfo<Scalar> trialStepCtrlInfo;
00675     {
00676 #ifdef ENABLE_RYTHMOS_TIMERS
00677       TEUCHOS_FUNC_TIME_MONITOR("Rythmos:SimpleIntegrator::advanceStepperToTime: getStepCtrl");
00678 #endif
00679       if (!is_null(integrationControlStrategy_)) {
00680         // Let an external strategy object determine the step size and type.
00681         // Note that any breakpoint info is also related through this call.
00682         trialStepCtrlInfo = integrationControlStrategy_->getNextStepControlInfo(
00683           *stepper_, stepCtrlInfoLast_, currTimeStepIndex_
00684           );
00685       }
00686       else {
00687         // Take a variable step if we have no control strategy
00688         trialStepCtrlInfo.stepType = STEP_TYPE_VARIABLE;
00689         trialStepCtrlInfo.stepSize = NL::max();
00690       }
00691     }
00692 
00693     // Print the initial trial step
00694     if ( includesVerbLevel(verbLevel,Teuchos::VERB_MEDIUM) ) {
00695       *out << "\nTrial step:\n";
00696       OSTab tab(out);
00697       *out << trialStepCtrlInfo;
00698     }
00699 
00700     // Halt immediately if we where told to do so
00701     if (trialStepCtrlInfo.stepSize < ST::zero()) {
00702       if ( includesVerbLevel(verbLevel,Teuchos::VERB_MEDIUM) )
00703         *out
00704           << "\n***"
00705           << "\n*** NOTICE: The IntegrationControlStrategy object return stepSize < 0.0, halting time integration!"
00706           << "\n***\n";
00707       return_val = false;
00708       break; // Exit the loop immediately!
00709     }
00710 
00711     // Make sure we don't step past the final time if asked not to
00712     bool updatedTrialStepCtrlInfo = false;
00713     {
00714       const Scalar finalTime = integrationTimeDomain_.upper();
00715       if (landOnFinalTime_ && trialStepCtrlInfo.stepSize + currStepperTimeRange.upper() > finalTime) {
00716         if ( includesVerbLevel(verbLevel,Teuchos::VERB_LOW) )
00717           *out << "\nCutting trial step to avoid stepping past final time ...\n";
00718         trialStepCtrlInfo.stepSize = finalTime - currStepperTimeRange.upper();
00719         updatedTrialStepCtrlInfo = true;
00720       }
00721     }
00722     
00723     // Print the modified trial step
00724     if ( updatedTrialStepCtrlInfo
00725       && includesVerbLevel(verbLevel,Teuchos::VERB_MEDIUM) )
00726     {
00727       *out << "\nUpdated trial step:\n";
00728       OSTab tab(out);
00729       *out << trialStepCtrlInfo;
00730     }
00731 
00732     //
00733     // C) Take the step
00734     //
00735 
00736     // Print step type and size
00737     if ( includesVerbLevel(verbLevel,Teuchos::VERB_MEDIUM) ) {
00738       if (trialStepCtrlInfo.stepType == STEP_TYPE_VARIABLE)
00739         *out << "\nTaking a variable time step with max step size = "
00740              << trialStepCtrlInfo.stepSize << " ....\n";
00741       else
00742         *out << "\nTaking a fixed time step of size = "
00743              << trialStepCtrlInfo.stepSize << " ....\n";
00744     }
00745 
00746     // Take step
00747     Scalar stepSizeTaken;
00748     {
00749 #ifdef ENABLE_RYTHMOS_TIMERS
00750       TEUCHOS_FUNC_TIME_MONITOR("Rythmos:SimpleIntegrator::advanceStepperToTime: takeStep");
00751 #endif
00752       stepSizeTaken = stepper_->takeStep(
00753         trialStepCtrlInfo.stepSize, trialStepCtrlInfo.stepType
00754         );
00755     }
00756 
00757     // Validate step taken
00758     if (trialStepCtrlInfo.stepType == STEP_TYPE_VARIABLE) {
00759       TEST_FOR_EXCEPTION(
00760         stepSizeTaken < ST::zero(), std::logic_error,
00761         "Error, stepper took negative step of dt = " << stepSizeTaken << "!\n"
00762         );
00763       TEST_FOR_EXCEPTION(
00764         stepSizeTaken > trialStepCtrlInfo.stepSize, std::logic_error,
00765         "Error, stepper took step of dt = " << stepSizeTaken
00766         << " > max step size of = " << trialStepCtrlInfo.stepSize << "!\n"
00767         );
00768     }
00769     else { // STEP_TYPE_FIXED
00770       TEST_FOR_EXCEPTION(
00771         stepSizeTaken != trialStepCtrlInfo.stepSize, std::logic_error,
00772         "Error, stepper took step of dt = " << stepSizeTaken 
00773         << " when asked to take step of dt = " << trialStepCtrlInfo.stepSize << "\n"
00774         );
00775     }
00776 
00777     // Update info about this step
00778     currStepperTimeRange = stepper_->getTimeRange();
00779     const StepControlInfo<Scalar> stepCtrlInfo =
00780       stepCtrlInfoTaken(trialStepCtrlInfo,stepSizeTaken);
00781 
00782     // Print the step actually taken 
00783     if ( includesVerbLevel(verbLevel,Teuchos::VERB_MEDIUM) ) {
00784       *out << "\nStep actually taken:\n";
00785       OSTab tab(out);
00786       *out << stepCtrlInfo;
00787     }
00788 
00789     //
00790     // D) Output info about step
00791     //
00792 
00793     {
00794 
00795 #ifdef ENABLE_RYTHMOS_TIMERS
00796       TEUCHOS_FUNC_TIME_MONITOR("Rythmos:SimpleIntegrator::advanceStepperToTime: output");
00797 #endif
00798       
00799       // Print our own brief output
00800       if ( includesVerbLevel(verbLevel,Teuchos::VERB_MEDIUM) ) {
00801         StepStatus<Scalar> stepStatus = stepper_->getStepStatus();
00802         *out << "\nTime point reached = " << stepStatus.time << endl;
00803         *out << "\nstepStatus:\n" << stepStatus;
00804         if ( includesVerbLevel(verbLevel,Teuchos::VERB_EXTREME) ) {
00805           RCP<const Thyra::VectorBase<Scalar> >
00806             solution = stepStatus.solution,
00807             solutionDot = stepStatus.solutionDot;
00808           if (!is_null(solution))
00809             *out << "\nsolution = \n" << Teuchos::describe(*solution,verbLevel);
00810           if (!is_null(solutionDot))
00811             *out << "\nsolutionDot = \n" << Teuchos::describe(*solutionDot,verbLevel);
00812         }
00813       }
00814       
00815       // Output to the observer
00816       if (!is_null(integrationObserver_))
00817         integrationObserver_->observeCompletedTimeStep(
00818           *stepper_, stepCtrlInfo, currTimeStepIndex_
00819           );
00820 
00821     }
00822 
00823     //
00824     // E) Update info for next time step
00825     //
00826 
00827     stepCtrlInfoLast_ = stepCtrlInfo;
00828     ++currTimeStepIndex_;
00829     
00830   }
00831 
00832   if ( includesVerbLevel(verbLevel,Teuchos::VERB_LOW) )
00833     *out << "\nNumber of steps taken in this call to advanceStepperToTime(...) = "
00834          << (currTimeStepIndex_ - initCurrTimeStepIndex) << endl
00835          << "\nLeaving" << this->Describable::description()
00836          << "::advanceStepperToTime("<<advance_to_t<<") ...\n";
00837 
00838   return return_val;
00839   
00840 }
00841 
00842 
00843 } // namespace Rythmos
00844 
00845 
00846 #endif //Rythmos_SIMPLE_INTEGRATOR_H

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