Rythmos - Transient Integration for Differential Equations Version of the Day
Rythmos_DefaultIntegrator_def.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_DEFAULT_INTEGRATOR_DEF_HPP
00030 #define RYTHMOS_DEFAULT_INTEGRATOR_DEF_HPP
00031 
00032 #include "Rythmos_DefaultIntegrator_decl.hpp"
00033 #include "Rythmos_InterpolationBufferHelpers.hpp"
00034 #include "Rythmos_IntegrationControlStrategyBase.hpp"
00035 #include "Rythmos_IntegrationObserverBase.hpp"
00036 #include "Rythmos_InterpolationBufferAppenderBase.hpp"
00037 #include "Rythmos_PointwiseInterpolationBufferAppender.hpp"
00038 #include "Rythmos_StepperHelpers.hpp"
00039 #include "Teuchos_VerboseObjectParameterListHelpers.hpp"
00040 #include "Teuchos_Assert.hpp"
00041 #include "Teuchos_as.hpp"
00042 #include <limits>
00043 
00044 namespace Rythmos {
00045 
00050 template<class Scalar> 
00051 RCP<DefaultIntegrator<Scalar> >
00052 defaultIntegrator()
00053 {
00054   RCP<DefaultIntegrator<Scalar> >
00055     integrator = Teuchos::rcp(new DefaultIntegrator<Scalar>());
00056   return integrator;
00057 }
00058 
00059 
00064 template<class Scalar> 
00065 RCP<DefaultIntegrator<Scalar> >
00066 defaultIntegrator(
00067   const RCP<IntegrationControlStrategyBase<Scalar> > &integrationControlStrategy,
00068   const RCP<IntegrationObserverBase<Scalar> > &integrationObserver
00069   )
00070 {
00071   RCP<DefaultIntegrator<Scalar> >
00072     integrator = Teuchos::rcp(new DefaultIntegrator<Scalar>());
00073   integrator->setIntegrationControlStrategy(integrationControlStrategy);
00074   integrator->setIntegrationObserver(integrationObserver);
00075   return integrator;
00076 }
00077 
00078 
00083 template<class Scalar> 
00084 RCP<DefaultIntegrator<Scalar> >
00085 controlledDefaultIntegrator(
00086   const RCP<IntegrationControlStrategyBase<Scalar> > &integrationControlStrategy
00087   )
00088 {
00089   RCP<DefaultIntegrator<Scalar> >
00090     integrator = Teuchos::rcp(new DefaultIntegrator<Scalar>());
00091   integrator->setIntegrationControlStrategy(integrationControlStrategy);
00092   return integrator;
00093 }
00094 
00095 
00100 template<class Scalar> 
00101 RCP<DefaultIntegrator<Scalar> >
00102 observedDefaultIntegrator(
00103   const RCP<IntegrationObserverBase<Scalar> > &integrationObserver
00104   )
00105 {
00106   RCP<DefaultIntegrator<Scalar> >
00107     integrator = Teuchos::rcp(new DefaultIntegrator<Scalar>());
00108   integrator->setIntegrationObserver(integrationObserver);
00109   return integrator;
00110 }
00111 
00112 
00113 
00114 // ////////////////////////////
00115 // Defintions
00116 
00117 
00118 // Static data members
00119 
00120 
00121 template<class Scalar>
00122 const std::string
00123 DefaultIntegrator<Scalar>::maxNumTimeSteps_name_ = "Max Number Time Steps";
00124 
00125 // replaced in code with std::numeric_limits<int>::max()
00126 // template<class Scalar>
00127 // const int
00128 // DefaultIntegrator<Scalar>::maxNumTimeSteps_default_ = 10000;
00129 
00130 
00131 // Constructors, Initializers, Misc
00132 
00133 
00134 template<class Scalar>
00135 DefaultIntegrator<Scalar>::DefaultIntegrator()
00136   :landOnFinalTime_(true),
00137    maxNumTimeSteps_(std::numeric_limits<int>::max()),
00138    currTimeStepIndex_(-1)
00139 {}
00140 
00141 
00142 template<class Scalar>
00143 void DefaultIntegrator<Scalar>::setIntegrationControlStrategy(
00144   const RCP<IntegrationControlStrategyBase<Scalar> > &integrationControlStrategy
00145   )
00146 {
00147 #ifdef HAVE_RYTHMOS_DEBUG
00148   TEUCHOS_TEST_FOR_EXCEPT(is_null(integrationControlStrategy));
00149 #endif
00150   integrationControlStrategy_ = integrationControlStrategy;
00151 }
00152 
00153 template<class Scalar>
00154 RCP<IntegrationControlStrategyBase<Scalar> > 
00155   DefaultIntegrator<Scalar>::getNonconstIntegrationControlStrategy()
00156 {
00157   return integrationControlStrategy_;
00158 }
00159 
00160 template<class Scalar>
00161 RCP<const IntegrationControlStrategyBase<Scalar> > 
00162   DefaultIntegrator<Scalar>::getIntegrationControlStrategy() const
00163 {
00164   return integrationControlStrategy_;
00165 }
00166 
00167 
00168 template<class Scalar>
00169 void DefaultIntegrator<Scalar>::setIntegrationObserver(
00170   const RCP<IntegrationObserverBase<Scalar> > &integrationObserver
00171   )
00172 {
00173 #ifdef HAVE_RYTHMOS_DEBUG
00174   TEUCHOS_TEST_FOR_EXCEPT(is_null(integrationObserver));
00175 #endif
00176   integrationObserver_ = integrationObserver;
00177 }
00178 
00179 
00180 template<class Scalar> 
00181 void DefaultIntegrator<Scalar>::setInterpolationBufferAppender(
00182   const RCP<InterpolationBufferAppenderBase<Scalar> > &interpBufferAppender
00183   )
00184 {
00185   interpBufferAppender_ = interpBufferAppender.assert_not_null();
00186 }
00187 
00188 
00189 template<class Scalar> 
00190 RCP<const InterpolationBufferAppenderBase<Scalar> >
00191 DefaultIntegrator<Scalar>::getInterpolationBufferAppender()
00192 {
00193   return interpBufferAppender_;
00194 }
00195 
00196 template<class Scalar> 
00197 RCP<InterpolationBufferAppenderBase<Scalar> >
00198 DefaultIntegrator<Scalar>::getNonconstInterpolationBufferAppender()
00199 {
00200   return interpBufferAppender_;
00201 }
00202 
00203 template<class Scalar> 
00204 RCP<InterpolationBufferAppenderBase<Scalar> >
00205 DefaultIntegrator<Scalar>::unSetInterpolationBufferAppender()
00206 {
00207   RCP<InterpolationBufferAppenderBase<Scalar> > InterpBufferAppender;
00208   std::swap( InterpBufferAppender, interpBufferAppender_ );
00209   return InterpBufferAppender;
00210 }
00211 
00212 
00213 // Overridden from ParameterListAcceptor
00214 
00215 
00216 template<class Scalar> 
00217 void DefaultIntegrator<Scalar>::setParameterList(
00218   RCP<ParameterList> const& paramList
00219   )
00220 {
00221   TEUCHOS_TEST_FOR_EXCEPT(is_null(paramList));
00222   paramList->validateParameters(*getValidParameters());
00223   this->setMyParamList(paramList);
00224   maxNumTimeSteps_ = paramList->get(
00225     maxNumTimeSteps_name_, std::numeric_limits<int>::max());
00226   Teuchos::readVerboseObjectSublist(&*paramList,this);
00227 }
00228 
00229 
00230 template<class Scalar> 
00231 RCP<const ParameterList>
00232 DefaultIntegrator<Scalar>::getValidParameters() const
00233 {
00234   static RCP<const ParameterList> validPL;
00235   if (is_null(validPL) ) {
00236     RCP<ParameterList> pl = Teuchos::parameterList();
00237     pl->set(maxNumTimeSteps_name_, std::numeric_limits<int>::max(),
00238       "Set the maximum number of integration time-steps allowed."
00239       );
00240     Teuchos::setupVerboseObjectSublist(&*pl);
00241     validPL = pl;
00242   }
00243   return validPL;
00244 }
00245 
00246 
00247 // Overridden from IntegratorBase
00248 
00249 
00250 template<class Scalar>
00251 RCP<IntegratorBase<Scalar> >
00252 DefaultIntegrator<Scalar>::cloneIntegrator() const
00253 {
00254 
00255   using Teuchos::null;
00256   RCP<DefaultIntegrator<Scalar> >
00257     newIntegrator = Teuchos::rcp(new DefaultIntegrator<Scalar>());
00258   // Only copy control information, not the stepper or the model it contains!
00259   newIntegrator->stepper_ = Teuchos::null;
00260   const RCP<const ParameterList> paramList = this->getParameterList();
00261   if (!is_null(paramList)) {
00262     newIntegrator->setParameterList(Teuchos::parameterList(*paramList));
00263   }
00264   if (!is_null(integrationControlStrategy_)) {
00265     newIntegrator->integrationControlStrategy_ =
00266       integrationControlStrategy_->cloneIntegrationControlStrategy().assert_not_null();
00267   }
00268   if (!is_null(integrationObserver_)) {
00269     newIntegrator->integrationObserver_ =
00270       integrationObserver_->cloneIntegrationObserver().assert_not_null();
00271   }
00272   if (!is_null(trailingInterpBuffer_)) {
00273     // ToDo: implement the clone!
00274     newIntegrator->trailingInterpBuffer_ = null;
00275     //newIntegrator->trailingInterpBuffer_ =
00276     //  trailingInterpBuffer_->cloneInterploationBuffer().assert_not_null();
00277   }
00278   if (!is_null(interpBufferAppender_)) {
00279     // ToDo: implement the clone!
00280     newIntegrator->interpBufferAppender_ = null;
00281     //newIntegrator->interpBufferAppender_ =
00282     //  interpBufferAppender_->cloneInterpolationBufferAppender().assert_not_null();
00283   }
00284   return newIntegrator;
00285 }
00286 
00287 
00288 template<class Scalar> 
00289 void DefaultIntegrator<Scalar>::setStepper(
00290   const RCP<StepperBase<Scalar> > &stepper,
00291   const Scalar &finalTime,
00292   const bool landOnFinalTime
00293   )
00294 {
00295   typedef Teuchos::ScalarTraits<Scalar> ST;
00296   TEUCHOS_TEST_FOR_EXCEPT(is_null(stepper));
00297   TEUCHOS_TEST_FOR_EXCEPT( finalTime <= stepper->getTimeRange().lower() );
00298   TEUCHOS_ASSERT( stepper->getTimeRange().length() == ST::zero() );
00299   // 2007/07/25: rabartl: ToDo: Validate state of the stepper!
00300   stepper_ = stepper;
00301   integrationTimeDomain_ = timeRange(stepper_->getTimeRange().lower(), finalTime);
00302   landOnFinalTime_ = landOnFinalTime;
00303   currTimeStepIndex_ = 0;
00304   stepCtrlInfoLast_ = StepControlInfo<Scalar>();
00305   if (!is_null(integrationControlStrategy_))
00306     integrationControlStrategy_->resetIntegrationControlStrategy(
00307       integrationTimeDomain_
00308       );
00309   if (!is_null(integrationObserver_))
00310     integrationObserver_->resetIntegrationObserver(
00311       integrationTimeDomain_
00312       );
00313 }
00314 
00315 
00316 template<class Scalar>
00317 RCP<StepperBase<Scalar> > DefaultIntegrator<Scalar>::unSetStepper()
00318 {
00319   RCP<StepperBase<Scalar> > stepper_temp = stepper_;
00320   stepper_ = Teuchos::null;
00321   return(stepper_temp);
00322 }
00323 
00324 
00325 template<class Scalar>
00326 RCP<const StepperBase<Scalar> > DefaultIntegrator<Scalar>::getStepper() const
00327 {
00328   return(stepper_);
00329 }
00330 
00331 
00332 template<class Scalar>
00333 RCP<StepperBase<Scalar> > DefaultIntegrator<Scalar>::getNonconstStepper() const
00334 {
00335   return(stepper_);
00336 }
00337 
00338 
00339 template<class Scalar>
00340 void DefaultIntegrator<Scalar>::setTrailingInterpolationBuffer(
00341   const RCP<InterpolationBufferBase<Scalar> > &trailingInterpBuffer
00342   )
00343 {
00344   trailingInterpBuffer_ = trailingInterpBuffer.assert_not_null();
00345 }
00346 
00347 
00348 template<class Scalar>
00349 RCP<InterpolationBufferBase<Scalar> >
00350 DefaultIntegrator<Scalar>::getNonconstTrailingInterpolationBuffer()
00351 {
00352   return trailingInterpBuffer_;
00353 }
00354 
00355 
00356 template<class Scalar>
00357 RCP<const InterpolationBufferBase<Scalar> >
00358 DefaultIntegrator<Scalar>::getTrailingInterpolationBuffer() const
00359 {
00360   return trailingInterpBuffer_;
00361 }
00362 
00363 template<class Scalar>
00364 RCP<InterpolationBufferBase<Scalar> >
00365 DefaultIntegrator<Scalar>::unSetTrailingInterpolationBuffer()
00366 {
00367   RCP<InterpolationBufferBase<Scalar> > trailingInterpBuffer;
00368   std::swap( trailingInterpBuffer, trailingInterpBuffer_ );
00369   return trailingInterpBuffer;
00370 }
00371 
00372 
00373 template<class Scalar>
00374 void DefaultIntegrator<Scalar>::getFwdPoints(
00375   const Array<Scalar>& time_vec,
00376   Array<RCP<const Thyra::VectorBase<Scalar> > >* x_vec,
00377   Array<RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec,
00378   Array<ScalarMag>* accuracy_vec
00379   )
00380 {
00381 
00382   RYTHMOS_FUNC_TIME_MONITOR_DIFF("Rythmos:DefaultIntegrator::getFwdPoints",
00383     TopLevel);
00384 
00385   using Teuchos::incrVerbLevel;
00386 #ifndef _MSC_VER
00387   using Teuchos::Describable;
00388 #endif
00389   typedef Teuchos::ScalarTraits<Scalar> ST;
00390   typedef InterpolationBufferBase<Scalar> IBB;
00391   typedef Teuchos::VerboseObjectTempState<IBB> VOTSIBB;
00392 
00393   finalizeSetup();
00394 
00395   RCP<Teuchos::FancyOStream> out = this->getOStream();
00396   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00397   Teuchos::OSTab tab(out);
00398   VOTSIBB stepper_outputTempState(stepper_,out,incrVerbLevel(verbLevel,-1));
00399 
00400   if ( includesVerbLevel(verbLevel,Teuchos::VERB_LOW) )
00401     *out << "\nEntering " << this->Describable::description() << "::getFwdPoints(...) ...\n"
00402          << "\nStepper: " << Teuchos::describe(*stepper_,verbLevel);
00403 
00404   if ( includesVerbLevel(verbLevel,Teuchos::VERB_MEDIUM) )
00405     *out << "\nRequested time points: " << Teuchos::toString(time_vec) << "\n";
00406 
00407   // Observe start of a time integration
00408   if (!is_null(integrationObserver_)) {
00409     integrationObserver_->setOStream(out);
00410     integrationObserver_->setVerbLevel(incrVerbLevel(verbLevel,-1));
00411     integrationObserver_->observeStartTimeIntegration(*stepper_);
00412   }
00413 
00414   //
00415   // 0) Initial setup
00416   //
00417 
00418   const int numTimePoints = time_vec.size();
00419 
00420   // Assert preconditions
00421   assertTimePointsAreSorted(time_vec);
00422   TEUCHOS_TEST_FOR_EXCEPT(accuracy_vec!=0); // ToDo: Remove accuracy_vec!
00423 
00424   // Resize the storage for the output arrays
00425   if (x_vec)
00426     x_vec->resize(numTimePoints);
00427   if (xdot_vec)
00428     xdot_vec->resize(numTimePoints);
00429 
00430   // This int records the next time point offset in time_vec[timePointIndex]
00431   // that needs to be handled.  This gets updated as the time points are
00432   // filled below.
00433   int nextTimePointIndex = 0;
00434   
00435   assertNoTimePointsBeforeCurrentTimeRange(*this,time_vec,nextTimePointIndex);
00436 
00437   //
00438   // 1) First, get all time points that fall within the current time range
00439   //
00440 
00441   {
00442     RYTHMOS_FUNC_TIME_MONITOR("Rythmos:DefaultIntegrator::getFwdPoints: getPoints");
00443     // 2007/10/05: rabartl: ToDo: Get points from trailingInterpBuffer_ first!
00444     getCurrentPoints(*stepper_,time_vec,x_vec,xdot_vec,&nextTimePointIndex);
00445   }
00446 
00447   //
00448   // 2) Advance the stepper to satisfy time points in time_vec that fall
00449   // before the current time.
00450   //
00451 
00452   while ( nextTimePointIndex < numTimePoints ) {
00453     
00454     // Use the time stepping algorithm to step up to or past the next
00455     // requested time but not so far as to step past the point entirely.
00456     const Scalar t = time_vec[nextTimePointIndex];
00457     bool advanceStepperToTimeSucceeded = false;
00458     {
00459       RYTHMOS_FUNC_TIME_MONITOR("Rythmos:DefaultIntegrator::getFwdPoints: advanceStepperToTime");
00460       advanceStepperToTimeSucceeded= advanceStepperToTime(t);
00461     }
00462     if (!advanceStepperToTimeSucceeded) {
00463       bool reachedMaxNumTimeSteps = (currTimeStepIndex_ >= maxNumTimeSteps_);
00464       if (reachedMaxNumTimeSteps) {
00465         // Break out of the while loop and attempt to exit gracefully.
00466         break;
00467       }
00468       TEUCHOS_TEST_FOR_EXCEPTION(
00469           !advanceStepperToTimeSucceeded, Exceptions::GetFwdPointsFailed,
00470           this->description() << "\n\n"
00471           "Error:  The integration failed to get to time " << t << " and only achieved\n"
00472           "getting to " << stepper_->getTimeRange().upper() << "!"
00473           );
00474     }
00475     
00476     // Extract the next set of points (perhaps just one) from the stepper
00477     {
00478       RYTHMOS_FUNC_TIME_MONITOR("Rythmos:DefaultIntegrator::getFwdPoints: getPoints (fwd)");
00479       getCurrentPoints(*stepper_,time_vec,x_vec,xdot_vec,&nextTimePointIndex);
00480     }
00481     
00482   }
00483 
00484   // Observe end of a time integration
00485   if (!is_null(integrationObserver_)) {
00486     integrationObserver_->observeEndTimeIntegration(*stepper_);
00487   }
00488 
00489   if ( includesVerbLevel(verbLevel,Teuchos::VERB_LOW) )
00490     *out << "\nLeaving " << this->Describable::description() << "::getFwdPoints(...) ...\n";
00491   
00492 }
00493 
00494 
00495 template<class Scalar> 
00496 TimeRange<Scalar>
00497 DefaultIntegrator<Scalar>::getFwdTimeRange() const
00498 {
00499   return timeRange(
00500     stepper_->getTimeRange().lower(),
00501     integrationTimeDomain_.upper()
00502     );
00503 }
00504 
00505 
00506 // Overridden from InterpolationBufferBase
00507 
00508 
00509 template<class Scalar> 
00510 RCP<const Thyra::VectorSpaceBase<Scalar> >
00511 DefaultIntegrator<Scalar>::get_x_space() const
00512 {
00513   return stepper_->get_x_space();
00514 }
00515 
00516 
00517 template<class Scalar> 
00518 void DefaultIntegrator<Scalar>::addPoints(
00519   const Array<Scalar>& time_vec,
00520   const Array<RCP<const Thyra::VectorBase<Scalar> > >& x_vec,
00521   const Array<RCP<const Thyra::VectorBase<Scalar> > >& xdot_vec
00522   )
00523 {
00524   stepper_->addPoints(time_vec,x_vec,xdot_vec);
00525 }
00526 
00527 
00528 template<class Scalar> 
00529 void DefaultIntegrator<Scalar>::getPoints(
00530   const Array<Scalar>& time_vec,
00531   Array<RCP<const Thyra::VectorBase<Scalar> > >* x_vec,
00532   Array<RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec,
00533   Array<ScalarMag>* accuracy_vec
00534   ) const
00535 {
00536 //  if (nonnull(trailingInterpBuffer_)) {
00537 //    int nextTimePointIndex = 0;
00538 //    getCurrentPoints(*trailingInterpBuffer_, time_vec, x_vec, xdot_vec, &nextTimePointIndex);
00539 //    getCurrentPoints(*stepper_, time_vec, x_vec, xdot_vec, &nextTimePointIndex);
00540 //    TEUCHOS_TEST_FOR_EXCEPTION( nextTimePointIndex < Teuchos::as<int>(time_vec.size()),
00541 //      std::out_of_range,
00542 //      "Error, the time point time_vec["<<nextTimePointIndex<<"] = "
00543 //      << time_vec[nextTimePointIndex] << " falls outside of the time range "
00544 //      << getTimeRange() << "!"
00545 //      );
00546 //  }
00547 //  else {
00548   stepper_->getPoints(time_vec, x_vec, xdot_vec, accuracy_vec);
00549 //  }
00550 }
00551 
00552 
00553 template<class Scalar> 
00554 TimeRange<Scalar> DefaultIntegrator<Scalar>::getTimeRange() const
00555 {
00556   if (nonnull(trailingInterpBuffer_)) {
00557     return timeRange(trailingInterpBuffer_->getTimeRange().lower(),
00558       stepper_->getTimeRange().upper());
00559   }
00560   return stepper_->getTimeRange();
00561 }
00562 
00563 
00564 template<class Scalar> 
00565 void DefaultIntegrator<Scalar>::getNodes(Array<Scalar>* time_vec) const
00566 {
00567   stepper_->getNodes(time_vec);
00568 }
00569 
00570 
00571 template<class Scalar> 
00572 void DefaultIntegrator<Scalar>::removeNodes(Array<Scalar>& time_vec)
00573 {
00574   stepper_->removeNodes(time_vec);
00575 }
00576 
00577 
00578 template<class Scalar> 
00579 int DefaultIntegrator<Scalar>::getOrder() const
00580 {
00581   return stepper_->getOrder();
00582 }
00583 
00584 
00585 // private
00586 
00587 
00588 template<class Scalar> 
00589 void DefaultIntegrator<Scalar>::finalizeSetup()
00590 {
00591   if (!is_null(trailingInterpBuffer_) && is_null(interpBufferAppender_))
00592     interpBufferAppender_ = pointwiseInterpolationBufferAppender<Scalar>();
00593   // ToDo: Do other setup?
00594 }
00595 
00596 
00597 template<class Scalar> 
00598 bool DefaultIntegrator<Scalar>::advanceStepperToTime( const Scalar& advance_to_t )
00599 {
00600 
00601   RYTHMOS_FUNC_TIME_MONITOR_DIFF("Rythmos:DefaultIntegrator::advanceStepperToTime",
00602     TopLevel);
00603 
00604   using std::endl;
00605   typedef std::numeric_limits<Scalar> NL;
00606   using Teuchos::incrVerbLevel; 
00607 #ifndef _MSC_VER
00608   using Teuchos::Describable;
00609 #endif
00610   using Teuchos::OSTab;
00611   typedef Teuchos::ScalarTraits<Scalar> ST;
00612 
00613   RCP<Teuchos::FancyOStream> out = this->getOStream();
00614   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00615 
00616   if (!is_null(integrationControlStrategy_)) {
00617     integrationControlStrategy_->setOStream(out);
00618     integrationControlStrategy_->setVerbLevel(incrVerbLevel(verbLevel,-1));
00619   }
00620 
00621   if (!is_null(integrationObserver_)) {
00622     integrationObserver_->setOStream(out);
00623     integrationObserver_->setVerbLevel(incrVerbLevel(verbLevel,-1));
00624   }
00625 
00626   if ( includesVerbLevel(verbLevel,Teuchos::VERB_LOW) )
00627     *out << "\nEntering " << this->Describable::description()
00628          << "::advanceStepperToTime("<<advance_to_t<<") ...\n";
00629 
00630   // Remember what timestep index we are on so we can report it later
00631   const int initCurrTimeStepIndex = currTimeStepIndex_;
00632 
00633   // Take steps until we the requested time is reached (or passed)
00634 
00635   TimeRange<Scalar> currStepperTimeRange = stepper_->getTimeRange();
00636 
00637   // Start by assume we can reach the time advance_to_t
00638   bool return_val = true;
00639   
00640   while ( !currStepperTimeRange.isInRange(advance_to_t) ) {
00641 
00642     // Halt immediatly if exceeded max iterations
00643     if (currTimeStepIndex_ >= maxNumTimeSteps_) {
00644       if ( includesVerbLevel(verbLevel,Teuchos::VERB_LOW) )
00645         *out
00646           << "\n***"
00647           << "\n*** NOTICE: currTimeStepIndex = "<<currTimeStepIndex_
00648           << " >= maxNumTimeSteps = "<<maxNumTimeSteps_<< ", halting time integration!"
00649           << "\n***\n";
00650       return_val = false;
00651       break; // Exit the loop immediately!
00652     }
00653 
00654     if ( includesVerbLevel(verbLevel,Teuchos::VERB_LOW) ) {
00655       *out << "\nTake step:  current_stepper_t = " << currStepperTimeRange.upper()
00656            << ", currTimeStepIndex = " << currTimeStepIndex_ << endl;
00657     }
00658 
00659     OSTab tab(out);
00660 
00661     //
00662     // A) Reinitialize if a hard breakpoint was reached on the last time step
00663     //
00664 
00665     if (stepCtrlInfoLast_.limitedByBreakPoint) {
00666       if ( stepCtrlInfoLast_.breakPointType == BREAK_POINT_TYPE_HARD ) {
00667         RYTHMOS_FUNC_TIME_MONITOR("Rythmos:DefaultIntegrator::restart");
00668         if ( includesVerbLevel(verbLevel,Teuchos::VERB_LOW) )
00669           *out << "\nAt a hard-breakpoint, restarting time integrator ...\n";
00670         restart(&*stepper_);
00671       }
00672       else  {
00673         if ( includesVerbLevel(verbLevel,Teuchos::VERB_LOW) )
00674           *out << "\nAt a soft-breakpoint, NOT restarting time integrator ...\n";
00675       }
00676     }
00677 
00678     //
00679     // B) Find an acceptable time step in a loop
00680     //
00681     // NOTE: Look for continue statements to iterate the loop!
00682     //
00683 
00684     bool foundAcceptableTimeStep = false;
00685     StepControlInfo<Scalar> stepCtrlInfo;
00686 
00687     // \todo Limit the maximum number of trial time steps to avoid an infinite
00688     // loop!
00689 
00690     while (!foundAcceptableTimeStep) {
00691 
00692       //
00693       // B.1) Get the trial step control info
00694       //
00695 
00696       StepControlInfo<Scalar> trialStepCtrlInfo;
00697       {
00698         RYTHMOS_FUNC_TIME_MONITOR("Rythmos:DefaultIntegrator::advanceStepperToTime: getStepCtrl");
00699         if (!is_null(integrationControlStrategy_)) {
00700           // Let an external strategy object determine the step size and type.
00701           // Note that any breakpoint info is also related through this call.
00702           trialStepCtrlInfo = integrationControlStrategy_->getNextStepControlInfo(
00703             *stepper_, stepCtrlInfoLast_, currTimeStepIndex_
00704             );
00705         }
00706         else {
00707           // Take a variable step if we have no control strategy
00708           trialStepCtrlInfo.stepType = STEP_TYPE_VARIABLE;
00709           trialStepCtrlInfo.stepSize = NL::max();
00710         }
00711       }
00712 
00713       // Print the initial trial step
00714       if ( includesVerbLevel(verbLevel,Teuchos::VERB_MEDIUM) ) {
00715         *out << "\nTrial step:\n";
00716         OSTab tab2(out);
00717         *out << trialStepCtrlInfo;
00718       }
00719 
00720       // Halt immediately if we where told to do so
00721       if (trialStepCtrlInfo.stepSize < ST::zero()) {
00722         if ( includesVerbLevel(verbLevel,Teuchos::VERB_MEDIUM) )
00723           *out
00724             << "\n***"
00725             << "\n*** NOTICE: The IntegrationControlStrategy object return stepSize < 0.0, halting time integration!"
00726             << "\n***\n";
00727         return_val = false;
00728         break; // Exit the loop immediately!
00729       }
00730 
00731       // Make sure we don't step past the final time if asked not to
00732       bool updatedTrialStepCtrlInfo = false;
00733       {
00734         const Scalar finalTime = integrationTimeDomain_.upper();
00735         if (landOnFinalTime_ && trialStepCtrlInfo.stepSize + currStepperTimeRange.upper() > finalTime) {
00736           if ( includesVerbLevel(verbLevel,Teuchos::VERB_LOW) )
00737             *out << "\nCutting trial step to avoid stepping past final time ...\n";
00738           trialStepCtrlInfo.stepSize = finalTime - currStepperTimeRange.upper();
00739           updatedTrialStepCtrlInfo = true;
00740         }
00741       }
00742     
00743       // Print the modified trial step
00744       if ( updatedTrialStepCtrlInfo
00745         && includesVerbLevel(verbLevel,Teuchos::VERB_MEDIUM) )
00746       {
00747         *out << "\nUpdated trial step:\n";
00748         OSTab tab2(out);
00749         *out << trialStepCtrlInfo;
00750       }
00751 
00752       //
00753       // B.2) Take the step
00754       //
00755 
00756       // Output to the observer we are starting a step
00757       if (!is_null(integrationObserver_))
00758   integrationObserver_->observeStartTimeStep(
00759             *stepper_, trialStepCtrlInfo, currTimeStepIndex_
00760             );
00761 
00762       // Print step type and size
00763       if ( includesVerbLevel(verbLevel,Teuchos::VERB_MEDIUM) ) {
00764         if (trialStepCtrlInfo.stepType == STEP_TYPE_VARIABLE)
00765           *out << "\nTaking a variable time step with max step size = "
00766                << trialStepCtrlInfo.stepSize << " ....\n";
00767         else
00768           *out << "\nTaking a fixed time step of size = "
00769                << trialStepCtrlInfo.stepSize << " ....\n";
00770       }
00771 
00772       // Take step
00773       Scalar stepSizeTaken;
00774       {
00775         RYTHMOS_FUNC_TIME_MONITOR("Rythmos:DefaultIntegrator::advanceStepperToTime: takeStep");
00776         stepSizeTaken = stepper_->takeStep(
00777           trialStepCtrlInfo.stepSize, trialStepCtrlInfo.stepType
00778           );
00779       }
00780 
00781       // Update info about this step
00782       currStepperTimeRange = stepper_->getTimeRange();
00783       stepCtrlInfo = stepCtrlInfoTaken(trialStepCtrlInfo,stepSizeTaken);
00784 
00785       // Print the step actually taken 
00786       if ( includesVerbLevel(verbLevel,Teuchos::VERB_MEDIUM) ) {
00787         *out << "\nStep actually taken:\n";
00788         OSTab tab2(out);
00789         *out << stepCtrlInfo;
00790       }
00791 
00792       // Determine if the timestep failed
00793       const bool timeStepFailed = (stepCtrlInfo.stepSize <= ST::zero());
00794       if (timeStepFailed && includesVerbLevel(verbLevel,Teuchos::VERB_MEDIUM)) {
00795         *out << "\nWARNING: timeStep = "<<trialStepCtrlInfo.stepSize<<" failed!\n";
00796       }
00797 
00798       // Notify observer of a failed time step
00799       if (timeStepFailed) {
00800         if (!is_null(integrationObserver_))
00801           integrationObserver_->observeFailedTimeStep(
00802             *stepper_, stepCtrlInfo, currTimeStepIndex_
00803             );
00804       }
00805 
00806       // Allow the IntegrationControlStrategy object to suggest another
00807       // timestep when a timestep fails.
00808       if (timeStepFailed && integrationControlStrategy_->handlesFailedTimeSteps())
00809       {
00810         // See if a new timestep can be suggested
00811         if (integrationControlStrategy_->resetForFailedTimeStep(
00812               *stepper_, stepCtrlInfoLast_, currTimeStepIndex_, trialStepCtrlInfo)
00813           )
00814         {
00815           if ( includesVerbLevel(verbLevel,Teuchos::VERB_MEDIUM) ) {
00816             *out << "\nThe IntegrationControlStrategy object indicated that"
00817                  << " it would like to suggest another timestep!\n";
00818           }
00819           // Skip the rest of the code in the loop and back to the top to try
00820           // another timestep!  Note: By doing this we skip the statement that
00821           // sets
00822           continue;
00823         }
00824         else
00825         {
00826           if ( includesVerbLevel(verbLevel,Teuchos::VERB_MEDIUM) ) {
00827             *out << "\nThe IntegrationControlStrategy object could not suggest"
00828                  << " a better time step!  Allowing to fail the time step!\n";
00829           }
00830           // Fall through to the failure checking!
00831         }
00832       }
00833 
00834       // Validate step taken
00835       if (trialStepCtrlInfo.stepType == STEP_TYPE_VARIABLE) {
00836         TEUCHOS_TEST_FOR_EXCEPTION(
00837           stepSizeTaken < ST::zero(), std::logic_error,
00838           "Error, stepper took negative step of dt = " << stepSizeTaken << "!\n"
00839           );
00840         TEUCHOS_TEST_FOR_EXCEPTION(
00841           stepSizeTaken > trialStepCtrlInfo.stepSize, std::logic_error,
00842           "Error, stepper took step of dt = " << stepSizeTaken
00843           << " > max step size of = " << trialStepCtrlInfo.stepSize << "!\n"
00844           );
00845       }
00846       else { // STEP_TYPE_FIXED
00847         TEUCHOS_TEST_FOR_EXCEPTION(
00848           stepSizeTaken != trialStepCtrlInfo.stepSize, std::logic_error,
00849           "Error, stepper took step of dt = " << stepSizeTaken 
00850           << " when asked to take step of dt = " << trialStepCtrlInfo.stepSize << "\n"
00851           );
00852       }
00853 
00854       // If we get here, the timestep is fine and is accepted!
00855       foundAcceptableTimeStep = true;
00856 
00857       // Append the trailing interpolation buffer (if defined)
00858       if (!is_null(trailingInterpBuffer_)) {
00859         interpBufferAppender_->append(*stepper_,currStepperTimeRange,
00860           trailingInterpBuffer_.ptr() );
00861       }
00862 
00863       //
00864       // B.3) Output info about step
00865       //
00866 
00867       {
00868 
00869         RYTHMOS_FUNC_TIME_MONITOR("Rythmos:DefaultIntegrator::advanceStepperToTime: output");
00870       
00871         // Print our own brief output
00872         if ( includesVerbLevel(verbLevel,Teuchos::VERB_MEDIUM) ) {
00873           StepStatus<Scalar> stepStatus = stepper_->getStepStatus();
00874           *out << "\nTime point reached = " << stepStatus.time << endl;
00875           *out << "\nstepStatus:\n" << stepStatus;
00876           if ( includesVerbLevel(verbLevel,Teuchos::VERB_EXTREME) ) {
00877             RCP<const Thyra::VectorBase<Scalar> >
00878               solution = stepStatus.solution,
00879               solutionDot = stepStatus.solutionDot;
00880             if (!is_null(solution))
00881               *out << "\nsolution = \n" << Teuchos::describe(*solution,verbLevel);
00882             if (!is_null(solutionDot))
00883               *out << "\nsolutionDot = \n" << Teuchos::describe(*solutionDot,verbLevel);
00884           }
00885         }
00886       
00887         // Output to the observer
00888         if (!is_null(integrationObserver_))
00889           integrationObserver_->observeCompletedTimeStep(
00890             *stepper_, stepCtrlInfo, currTimeStepIndex_
00891             );
00892 
00893       }
00894 
00895     } // end loop to find a valid time step
00896 
00897     //
00898     // C) Update info for next time step
00899     //
00900 
00901     stepCtrlInfoLast_ = stepCtrlInfo;
00902     ++currTimeStepIndex_;
00903     
00904   }
00905 
00906   if ( includesVerbLevel(verbLevel,Teuchos::VERB_LOW) )
00907     *out << "\nNumber of steps taken in this call to advanceStepperToTime(...) = "
00908          << (currTimeStepIndex_ - initCurrTimeStepIndex) << endl
00909          << "\nLeaving" << this->Describable::description()
00910          << "::advanceStepperToTime("<<advance_to_t<<") ...\n";
00911 
00912   return return_val;
00913   
00914 }
00915 
00916 // 
00917 // Explicit Instantiation macro
00918 //
00919 // Must be expanded from within the Rythmos namespace!
00920 //
00921 
00922 #define RYTHMOS_DEFAULT_INTEGRATOR_INSTANT(SCALAR) \
00923   \
00924   template class DefaultIntegrator< SCALAR >; \
00925   \
00926   template RCP<DefaultIntegrator< SCALAR > > \
00927   defaultIntegrator(); \
00928   \
00929   template RCP<DefaultIntegrator< SCALAR > > \
00930   defaultIntegrator( \
00931     const RCP<IntegrationControlStrategyBase< SCALAR > > &integrationControlStrategy, \
00932     const RCP<IntegrationObserverBase< SCALAR > > &integrationObserver \
00933     ); \
00934   \
00935   template RCP<DefaultIntegrator< SCALAR > > \
00936   controlledDefaultIntegrator( \
00937     const RCP<IntegrationControlStrategyBase< SCALAR > > &integrationControlStrategy \
00938     ); 
00939 
00940 
00941 } // namespace Rythmos
00942 
00943 
00944 #endif //RYTHMOS_DEFAULT_INTEGRATOR_DEF_HPP
 All Classes Functions Variables Typedefs Friends