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