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 template<class Scalar>
00125 const std::string
00126 DefaultIntegrator<Scalar>::observeInitCond_name_ = "Observe Initial Condition";
00127 
00128 // replaced in code with std::numeric_limits<int>::max()
00129 // template<class Scalar>
00130 // const int
00131 // DefaultIntegrator<Scalar>::maxNumTimeSteps_default_ = 10000;
00132 
00133 
00134 // Constructors, Initializers, Misc
00135 
00136 
00137 template<class Scalar>
00138 DefaultIntegrator<Scalar>::DefaultIntegrator()
00139   :landOnFinalTime_(true),
00140    maxNumTimeSteps_(std::numeric_limits<int>::max()),
00141    observeInitCond_(false),
00142    currTimeStepIndex_(-1)
00143 {}
00144 
00145 
00146 template<class Scalar>
00147 void DefaultIntegrator<Scalar>::setIntegrationControlStrategy(
00148   const RCP<IntegrationControlStrategyBase<Scalar> > &integrationControlStrategy
00149   )
00150 {
00151 #ifdef HAVE_RYTHMOS_DEBUG
00152   TEUCHOS_TEST_FOR_EXCEPT(is_null(integrationControlStrategy));
00153 #endif
00154   integrationControlStrategy_ = integrationControlStrategy;
00155 }
00156 
00157 template<class Scalar>
00158 RCP<IntegrationControlStrategyBase<Scalar> > 
00159   DefaultIntegrator<Scalar>::getNonconstIntegrationControlStrategy()
00160 {
00161   return integrationControlStrategy_;
00162 }
00163 
00164 template<class Scalar>
00165 RCP<const IntegrationControlStrategyBase<Scalar> > 
00166   DefaultIntegrator<Scalar>::getIntegrationControlStrategy() const
00167 {
00168   return integrationControlStrategy_;
00169 }
00170 
00171 
00172 template<class Scalar>
00173 void DefaultIntegrator<Scalar>::setIntegrationObserver(
00174   const RCP<IntegrationObserverBase<Scalar> > &integrationObserver
00175   )
00176 {
00177 #ifdef HAVE_RYTHMOS_DEBUG
00178   TEUCHOS_TEST_FOR_EXCEPT(is_null(integrationObserver));
00179 #endif
00180   integrationObserver_ = integrationObserver;
00181 }
00182 
00183 
00184 template<class Scalar> 
00185 void DefaultIntegrator<Scalar>::setInterpolationBufferAppender(
00186   const RCP<InterpolationBufferAppenderBase<Scalar> > &interpBufferAppender
00187   )
00188 {
00189   interpBufferAppender_ = interpBufferAppender.assert_not_null();
00190 }
00191 
00192 
00193 template<class Scalar> 
00194 RCP<const InterpolationBufferAppenderBase<Scalar> >
00195 DefaultIntegrator<Scalar>::getInterpolationBufferAppender()
00196 {
00197   return interpBufferAppender_;
00198 }
00199 
00200 template<class Scalar> 
00201 RCP<InterpolationBufferAppenderBase<Scalar> >
00202 DefaultIntegrator<Scalar>::getNonconstInterpolationBufferAppender()
00203 {
00204   return interpBufferAppender_;
00205 }
00206 
00207 template<class Scalar> 
00208 RCP<InterpolationBufferAppenderBase<Scalar> >
00209 DefaultIntegrator<Scalar>::unSetInterpolationBufferAppender()
00210 {
00211   RCP<InterpolationBufferAppenderBase<Scalar> > InterpBufferAppender;
00212   std::swap( InterpBufferAppender, interpBufferAppender_ );
00213   return InterpBufferAppender;
00214 }
00215 
00216 
00217 // Overridden from ParameterListAcceptor
00218 
00219 
00220 template<class Scalar> 
00221 void DefaultIntegrator<Scalar>::setParameterList(
00222   RCP<ParameterList> const& paramList
00223   )
00224 {
00225   TEUCHOS_TEST_FOR_EXCEPT(is_null(paramList));
00226   paramList->validateParameters(*getValidParameters());
00227   this->setMyParamList(paramList);
00228   maxNumTimeSteps_ = paramList->get(
00229     maxNumTimeSteps_name_, std::numeric_limits<int>::max());
00230   observeInitCond_ = paramList->get(observeInitCond_name_, false);
00231   Teuchos::readVerboseObjectSublist(&*paramList,this);
00232 }
00233 
00234 
00235 template<class Scalar> 
00236 RCP<const ParameterList>
00237 DefaultIntegrator<Scalar>::getValidParameters() const
00238 {
00239   static RCP<const ParameterList> validPL;
00240   if (is_null(validPL) ) {
00241     RCP<ParameterList> pl = Teuchos::parameterList();
00242     pl->set(maxNumTimeSteps_name_, std::numeric_limits<int>::max(),
00243       "Set the maximum number of integration time-steps allowed."
00244       );
00245     pl->set(observeInitCond_name_, false,
00246       "Flag to call Observer with initial conditions beofre integration."
00247       );
00248     Teuchos::setupVerboseObjectSublist(&*pl);
00249     validPL = pl;
00250   }
00251   return validPL;
00252 }
00253 
00254 
00255 // Overridden from IntegratorBase
00256 
00257 
00258 template<class Scalar>
00259 RCP<IntegratorBase<Scalar> >
00260 DefaultIntegrator<Scalar>::cloneIntegrator() const
00261 {
00262 
00263   using Teuchos::null;
00264   RCP<DefaultIntegrator<Scalar> >
00265     newIntegrator = Teuchos::rcp(new DefaultIntegrator<Scalar>());
00266   // Only copy control information, not the stepper or the model it contains!
00267   newIntegrator->stepper_ = Teuchos::null;
00268   const RCP<const ParameterList> paramList = this->getParameterList();
00269   if (!is_null(paramList)) {
00270     newIntegrator->setParameterList(Teuchos::parameterList(*paramList));
00271   }
00272   if (!is_null(integrationControlStrategy_)) {
00273     newIntegrator->integrationControlStrategy_ =
00274       integrationControlStrategy_->cloneIntegrationControlStrategy().assert_not_null();
00275   }
00276   if (!is_null(integrationObserver_)) {
00277     newIntegrator->integrationObserver_ =
00278       integrationObserver_->cloneIntegrationObserver().assert_not_null();
00279   }
00280   if (!is_null(trailingInterpBuffer_)) {
00281     // ToDo: implement the clone!
00282     newIntegrator->trailingInterpBuffer_ = null;
00283     //newIntegrator->trailingInterpBuffer_ =
00284     //  trailingInterpBuffer_->cloneInterploationBuffer().assert_not_null();
00285   }
00286   if (!is_null(interpBufferAppender_)) {
00287     // ToDo: implement the clone!
00288     newIntegrator->interpBufferAppender_ = null;
00289     //newIntegrator->interpBufferAppender_ =
00290     //  interpBufferAppender_->cloneInterpolationBufferAppender().assert_not_null();
00291   }
00292   return newIntegrator;
00293 }
00294 
00295 
00296 template<class Scalar> 
00297 void DefaultIntegrator<Scalar>::setStepper(
00298   const RCP<StepperBase<Scalar> > &stepper,
00299   const Scalar &finalTime,
00300   const bool landOnFinalTime
00301   )
00302 {
00303   typedef Teuchos::ScalarTraits<Scalar> ST;
00304   TEUCHOS_TEST_FOR_EXCEPT(is_null(stepper));
00305   TEUCHOS_TEST_FOR_EXCEPT( finalTime <= stepper->getTimeRange().lower() );
00306   TEUCHOS_ASSERT( stepper->getTimeRange().length() == ST::zero() );
00307   // 2007/07/25: rabartl: ToDo: Validate state of the stepper!
00308   stepper_ = stepper;
00309   integrationTimeDomain_ = timeRange(stepper_->getTimeRange().lower(), finalTime);
00310   landOnFinalTime_ = landOnFinalTime;
00311   currTimeStepIndex_ = 0;
00312   stepCtrlInfoLast_ = StepControlInfo<Scalar>();
00313   if (!is_null(integrationControlStrategy_))
00314     integrationControlStrategy_->resetIntegrationControlStrategy(
00315       integrationTimeDomain_
00316       );
00317   if (!is_null(integrationObserver_))
00318     integrationObserver_->resetIntegrationObserver(
00319       integrationTimeDomain_
00320       );
00321 }
00322 
00323 
00324 template<class Scalar>
00325 RCP<StepperBase<Scalar> > DefaultIntegrator<Scalar>::unSetStepper()
00326 {
00327   RCP<StepperBase<Scalar> > stepper_temp = stepper_;
00328   stepper_ = Teuchos::null;
00329   return(stepper_temp);
00330 }
00331 
00332 
00333 template<class Scalar>
00334 RCP<const StepperBase<Scalar> > DefaultIntegrator<Scalar>::getStepper() const
00335 {
00336   return(stepper_);
00337 }
00338 
00339 
00340 template<class Scalar>
00341 RCP<StepperBase<Scalar> > DefaultIntegrator<Scalar>::getNonconstStepper() const
00342 {
00343   return(stepper_);
00344 }
00345 
00346 
00347 template<class Scalar>
00348 void DefaultIntegrator<Scalar>::setTrailingInterpolationBuffer(
00349   const RCP<InterpolationBufferBase<Scalar> > &trailingInterpBuffer
00350   )
00351 {
00352   trailingInterpBuffer_ = trailingInterpBuffer.assert_not_null();
00353 }
00354 
00355 
00356 template<class Scalar>
00357 RCP<InterpolationBufferBase<Scalar> >
00358 DefaultIntegrator<Scalar>::getNonconstTrailingInterpolationBuffer()
00359 {
00360   return trailingInterpBuffer_;
00361 }
00362 
00363 
00364 template<class Scalar>
00365 RCP<const InterpolationBufferBase<Scalar> >
00366 DefaultIntegrator<Scalar>::getTrailingInterpolationBuffer() const
00367 {
00368   return trailingInterpBuffer_;
00369 }
00370 
00371 template<class Scalar>
00372 RCP<InterpolationBufferBase<Scalar> >
00373 DefaultIntegrator<Scalar>::unSetTrailingInterpolationBuffer()
00374 {
00375   RCP<InterpolationBufferBase<Scalar> > trailingInterpBuffer;
00376   std::swap( trailingInterpBuffer, trailingInterpBuffer_ );
00377   return trailingInterpBuffer;
00378 }
00379 
00380 
00381 template<class Scalar>
00382 void DefaultIntegrator<Scalar>::getFwdPoints(
00383   const Array<Scalar>& time_vec,
00384   Array<RCP<const Thyra::VectorBase<Scalar> > >* x_vec,
00385   Array<RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec,
00386   Array<ScalarMag>* accuracy_vec
00387   )
00388 {
00389 
00390   RYTHMOS_FUNC_TIME_MONITOR_DIFF("Rythmos:DefaultIntegrator::getFwdPoints",
00391     TopLevel);
00392 
00393   using Teuchos::incrVerbLevel;
00394 #ifndef _MSC_VER
00395   using Teuchos::Describable;
00396 #endif
00397   typedef Teuchos::ScalarTraits<Scalar> ST;
00398   typedef InterpolationBufferBase<Scalar> IBB;
00399   typedef Teuchos::VerboseObjectTempState<IBB> VOTSIBB;
00400 
00401   finalizeSetup();
00402 
00403   RCP<Teuchos::FancyOStream> out = this->getOStream();
00404   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00405   Teuchos::OSTab tab(out);
00406   VOTSIBB stepper_outputTempState(stepper_,out,incrVerbLevel(verbLevel,-1));
00407 
00408   if ( includesVerbLevel(verbLevel,Teuchos::VERB_LOW) )
00409     *out << "\nEntering " << this->Describable::description() << "::getFwdPoints(...) ...\n"
00410          << "\nStepper: " << Teuchos::describe(*stepper_,verbLevel);
00411 
00412   if ( includesVerbLevel(verbLevel,Teuchos::VERB_MEDIUM) )
00413     *out << "\nRequested time points: " << Teuchos::toString(time_vec) << "\n";
00414 
00415   // Observe start of a time integration
00416   if (!is_null(integrationObserver_)) {
00417     integrationObserver_->setOStream(out);
00418     integrationObserver_->setVerbLevel(incrVerbLevel(verbLevel,-1));
00419     integrationObserver_->observeStartTimeIntegration(*stepper_);
00420 
00421   // AGS: attempt to write observe init conditions
00422   if ( observeInitCond_ )
00423      integrationObserver_->observeCompletedTimeStep(
00424       *stepper_, stepCtrlInfoLast_, 0);
00425   }
00426 
00427   //
00428   // 0) Initial setup
00429   //
00430 
00431   const int numTimePoints = time_vec.size();
00432 
00433   // Assert preconditions
00434   assertTimePointsAreSorted(time_vec);
00435   TEUCHOS_TEST_FOR_EXCEPT(accuracy_vec!=0); // ToDo: Remove accuracy_vec!
00436 
00437   // Resize the storage for the output arrays
00438   if (x_vec)
00439     x_vec->resize(numTimePoints);
00440   if (xdot_vec)
00441     xdot_vec->resize(numTimePoints);
00442 
00443   // This int records the next time point offset in time_vec[timePointIndex]
00444   // that needs to be handled.  This gets updated as the time points are
00445   // filled below.
00446   int nextTimePointIndex = 0;
00447   
00448   assertNoTimePointsBeforeCurrentTimeRange(*this,time_vec,nextTimePointIndex);
00449 
00450   //
00451   // 1) First, get all time points that fall within the current time range
00452   //
00453 
00454   {
00455     RYTHMOS_FUNC_TIME_MONITOR("Rythmos:DefaultIntegrator::getFwdPoints: getPoints");
00456     // 2007/10/05: rabartl: ToDo: Get points from trailingInterpBuffer_ first!
00457     getCurrentPoints(*stepper_,time_vec,x_vec,xdot_vec,&nextTimePointIndex);
00458   }
00459 
00460   //
00461   // 2) Advance the stepper to satisfy time points in time_vec that fall
00462   // before the current time.
00463   //
00464 
00465   while ( nextTimePointIndex < numTimePoints ) {
00466     
00467     // Use the time stepping algorithm to step up to or past the next
00468     // requested time but not so far as to step past the point entirely.
00469     const Scalar t = time_vec[nextTimePointIndex];
00470     bool advanceStepperToTimeSucceeded = false;
00471     {
00472       RYTHMOS_FUNC_TIME_MONITOR("Rythmos:DefaultIntegrator::getFwdPoints: advanceStepperToTime");
00473       advanceStepperToTimeSucceeded= advanceStepperToTime(t);
00474     }
00475     if (!advanceStepperToTimeSucceeded) {
00476       bool reachedMaxNumTimeSteps = (currTimeStepIndex_ >= maxNumTimeSteps_);
00477       if (reachedMaxNumTimeSteps) {
00478         // Break out of the while loop and attempt to exit gracefully.
00479         break;
00480       }
00481       TEUCHOS_TEST_FOR_EXCEPTION(
00482           !advanceStepperToTimeSucceeded, Exceptions::GetFwdPointsFailed,
00483           this->description() << "\n\n"
00484           "Error:  The integration failed to get to time " << t << " and only achieved\n"
00485           "getting to " << stepper_->getTimeRange().upper() << "!"
00486           );
00487     }
00488     
00489     // Extract the next set of points (perhaps just one) from the stepper
00490     {
00491       RYTHMOS_FUNC_TIME_MONITOR("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() << "::getFwdPoints(...) ...\n";
00504   
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("Rythmos:DefaultIntegrator::advanceStepperToTime",
00615     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 we the requested time is reached (or passed)
00647 
00648   TimeRange<Scalar> currStepperTimeRange = stepper_->getTimeRange();
00649 
00650   // Start by assume we can reach the time advance_to_t
00651   bool return_val = true;
00652 
00653   
00654   while ( !currStepperTimeRange.isInRange(advance_to_t) ) {
00655 
00656     // Halt immediatly if exceeded max iterations
00657     if (currTimeStepIndex_ >= maxNumTimeSteps_) {
00658       if ( includesVerbLevel(verbLevel,Teuchos::VERB_LOW) )
00659         *out
00660           << "\n***"
00661           << "\n*** NOTICE: currTimeStepIndex = "<<currTimeStepIndex_
00662           << " >= maxNumTimeSteps = "<<maxNumTimeSteps_<< ", 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 = " << currStepperTimeRange.upper()
00670            << ", currTimeStepIndex = " << currTimeStepIndex_ << endl;
00671     }
00672 
00673     OSTab tab(out);
00674 
00675     //
00676     // A) Reinitialize if a hard breakpoint was reached on the last time step
00677     //
00678 
00679     if (stepCtrlInfoLast_.limitedByBreakPoint) {
00680       if ( stepCtrlInfoLast_.breakPointType == BREAK_POINT_TYPE_HARD ) {
00681         RYTHMOS_FUNC_TIME_MONITOR("Rythmos:DefaultIntegrator::restart");
00682         if ( includesVerbLevel(verbLevel,Teuchos::VERB_LOW) )
00683           *out << "\nAt a hard-breakpoint, restarting time integrator ...\n";
00684         restart(&*stepper_);
00685       }
00686       else  {
00687         if ( includesVerbLevel(verbLevel,Teuchos::VERB_LOW) )
00688           *out << "\nAt a soft-breakpoint, NOT restarting time integrator ...\n";
00689       }
00690     }
00691 
00692     //
00693     // B) Find an acceptable time step in a loop
00694     //
00695     // NOTE: Look for continue statements to iterate the loop!
00696     //
00697 
00698     bool foundAcceptableTimeStep = false;
00699     StepControlInfo<Scalar> stepCtrlInfo;
00700 
00701     // \todo Limit the maximum number of trial time steps to avoid an infinite
00702     // loop!
00703 
00704     while (!foundAcceptableTimeStep) {
00705 
00706       //
00707       // B.1) Get the trial step control info
00708       //
00709 
00710       StepControlInfo<Scalar> trialStepCtrlInfo;
00711       {
00712         RYTHMOS_FUNC_TIME_MONITOR("Rythmos:DefaultIntegrator::advanceStepperToTime: getStepCtrl");
00713         if (!is_null(integrationControlStrategy_)) {
00714           // Let an external strategy object determine the step size and type.
00715           // Note that any breakpoint info is also related through this call.
00716           trialStepCtrlInfo = integrationControlStrategy_->getNextStepControlInfo(
00717             *stepper_, stepCtrlInfoLast_, currTimeStepIndex_
00718             );
00719         }
00720         else {
00721           // Take a variable step if we have no control strategy
00722           trialStepCtrlInfo.stepType = STEP_TYPE_VARIABLE;
00723           trialStepCtrlInfo.stepSize = NL::max();
00724         }
00725       }
00726 
00727       // Print the initial trial step
00728       if ( includesVerbLevel(verbLevel,Teuchos::VERB_MEDIUM) ) {
00729         *out << "\nTrial step:\n";
00730         OSTab tab2(out);
00731         *out << trialStepCtrlInfo;
00732       }
00733 
00734       // Halt immediately if we where told to do so
00735       if (trialStepCtrlInfo.stepSize < ST::zero()) {
00736         if ( includesVerbLevel(verbLevel,Teuchos::VERB_MEDIUM) )
00737           *out
00738             << "\n***"
00739             << "\n*** NOTICE: The IntegrationControlStrategy object return stepSize < 0.0, halting time integration!"
00740             << "\n***\n";
00741         return_val = false;
00742         break; // Exit the loop immediately!
00743       }
00744 
00745       // Make sure we don't step past the final time if asked not to
00746       bool updatedTrialStepCtrlInfo = false;
00747       {
00748         const Scalar finalTime = integrationTimeDomain_.upper();
00749         if (landOnFinalTime_ && trialStepCtrlInfo.stepSize + currStepperTimeRange.upper() > finalTime) {
00750           if ( includesVerbLevel(verbLevel,Teuchos::VERB_LOW) )
00751             *out << "\nCutting trial step to avoid stepping past final time ...\n";
00752           trialStepCtrlInfo.stepSize = finalTime - currStepperTimeRange.upper();
00753           updatedTrialStepCtrlInfo = true;
00754         }
00755       }
00756     
00757       // Print the modified trial step
00758       if ( updatedTrialStepCtrlInfo
00759         && includesVerbLevel(verbLevel,Teuchos::VERB_MEDIUM) )
00760       {
00761         *out << "\nUpdated trial step:\n";
00762         OSTab tab2(out);
00763         *out << trialStepCtrlInfo;
00764       }
00765 
00766       //
00767       // B.2) Take the step
00768       //
00769 
00770       // Output to the observer we are starting a step
00771       if (!is_null(integrationObserver_))
00772   integrationObserver_->observeStartTimeStep(
00773             *stepper_, trialStepCtrlInfo, currTimeStepIndex_
00774             );
00775 
00776       // Print step type and size
00777       if ( includesVerbLevel(verbLevel,Teuchos::VERB_MEDIUM) ) {
00778         if (trialStepCtrlInfo.stepType == STEP_TYPE_VARIABLE)
00779           *out << "\nTaking a variable time step with max step size = "
00780                << trialStepCtrlInfo.stepSize << " ....\n";
00781         else
00782           *out << "\nTaking a fixed time step of size = "
00783                << trialStepCtrlInfo.stepSize << " ....\n";
00784       }
00785 
00786       // Take step
00787       Scalar stepSizeTaken;
00788       {
00789         RYTHMOS_FUNC_TIME_MONITOR("Rythmos:DefaultIntegrator::advanceStepperToTime: takeStep");
00790         stepSizeTaken = stepper_->takeStep(
00791           trialStepCtrlInfo.stepSize, trialStepCtrlInfo.stepType
00792           );
00793       }
00794 
00795       // Update info about this step
00796       currStepperTimeRange = stepper_->getTimeRange();
00797       stepCtrlInfo = stepCtrlInfoTaken(trialStepCtrlInfo,stepSizeTaken);
00798 
00799       // Print the step actually taken 
00800       if ( includesVerbLevel(verbLevel,Teuchos::VERB_MEDIUM) ) {
00801         *out << "\nStep actually taken:\n";
00802         OSTab tab2(out);
00803         *out << stepCtrlInfo;
00804       }
00805 
00806       // Determine if the timestep failed
00807       const bool timeStepFailed = (stepCtrlInfo.stepSize <= ST::zero());
00808       if (timeStepFailed && includesVerbLevel(verbLevel,Teuchos::VERB_MEDIUM)) {
00809         *out << "\nWARNING: timeStep = "<<trialStepCtrlInfo.stepSize<<" failed!\n";
00810       }
00811 
00812       // Notify observer of a failed time step
00813       if (timeStepFailed) {
00814         if (!is_null(integrationObserver_))
00815           integrationObserver_->observeFailedTimeStep(
00816             *stepper_, stepCtrlInfo, currTimeStepIndex_
00817             );
00818       }
00819 
00820       // Allow the IntegrationControlStrategy object to suggest another
00821       // timestep when a timestep fails.
00822       if (timeStepFailed && integrationControlStrategy_->handlesFailedTimeSteps())
00823       {
00824         // See if a new timestep can be suggested
00825         if (integrationControlStrategy_->resetForFailedTimeStep(
00826               *stepper_, stepCtrlInfoLast_, currTimeStepIndex_, trialStepCtrlInfo)
00827           )
00828         {
00829           if ( includesVerbLevel(verbLevel,Teuchos::VERB_MEDIUM) ) {
00830             *out << "\nThe IntegrationControlStrategy object indicated that"
00831                  << " it would like to suggest another timestep!\n";
00832           }
00833           // Skip the rest of the code in the loop and back to the top to try
00834           // another timestep!  Note: By doing this we skip the statement that
00835           // sets
00836           continue;
00837         }
00838         else
00839         {
00840           if ( includesVerbLevel(verbLevel,Teuchos::VERB_MEDIUM) ) {
00841             *out << "\nThe IntegrationControlStrategy object could not suggest"
00842                  << " a better time step!  Allowing to fail the time step!\n";
00843           }
00844           // Fall through to the failure checking!
00845         }
00846       }
00847 
00848       // Validate step taken
00849       if (trialStepCtrlInfo.stepType == STEP_TYPE_VARIABLE) {
00850         TEUCHOS_TEST_FOR_EXCEPTION(
00851           stepSizeTaken < ST::zero(), std::logic_error,
00852           "Error, stepper took negative step of dt = " << stepSizeTaken << "!\n"
00853           );
00854         TEUCHOS_TEST_FOR_EXCEPTION(
00855           stepSizeTaken > trialStepCtrlInfo.stepSize, std::logic_error,
00856           "Error, stepper took step of dt = " << stepSizeTaken
00857           << " > max step size of = " << trialStepCtrlInfo.stepSize << "!\n"
00858           );
00859       }
00860       else { // STEP_TYPE_FIXED
00861         TEUCHOS_TEST_FOR_EXCEPTION(
00862           stepSizeTaken != trialStepCtrlInfo.stepSize, std::logic_error,
00863           "Error, stepper took step of dt = " << stepSizeTaken 
00864           << " when asked to take step of dt = " << trialStepCtrlInfo.stepSize << "\n"
00865           );
00866       }
00867 
00868       // If we get here, the timestep is fine and is accepted!
00869       foundAcceptableTimeStep = true;
00870 
00871       // Append the trailing interpolation buffer (if defined)
00872       if (!is_null(trailingInterpBuffer_)) {
00873         interpBufferAppender_->append(*stepper_,currStepperTimeRange,
00874           trailingInterpBuffer_.ptr() );
00875       }
00876 
00877       //
00878       // B.3) Output info about step
00879       //
00880 
00881       {
00882 
00883         RYTHMOS_FUNC_TIME_MONITOR("Rythmos:DefaultIntegrator::advanceStepperToTime: output");
00884       
00885         // Print our own brief output
00886         if ( includesVerbLevel(verbLevel,Teuchos::VERB_MEDIUM) ) {
00887           StepStatus<Scalar> stepStatus = stepper_->getStepStatus();
00888           *out << "\nTime point reached = " << stepStatus.time << endl;
00889           *out << "\nstepStatus:\n" << stepStatus;
00890           if ( includesVerbLevel(verbLevel,Teuchos::VERB_EXTREME) ) {
00891             RCP<const Thyra::VectorBase<Scalar> >
00892               solution = stepStatus.solution,
00893               solutionDot = stepStatus.solutionDot;
00894             if (!is_null(solution))
00895               *out << "\nsolution = \n" << Teuchos::describe(*solution,verbLevel);
00896             if (!is_null(solutionDot))
00897               *out << "\nsolutionDot = \n" << Teuchos::describe(*solutionDot,verbLevel);
00898           }
00899         }
00900       
00901         // Output to the observer
00902         if (!is_null(integrationObserver_))
00903           integrationObserver_->observeCompletedTimeStep(
00904             *stepper_, stepCtrlInfo, currTimeStepIndex_
00905             );
00906 
00907       }
00908 
00909     } // end loop to find a valid time step
00910 
00911     //
00912     // C) Update info for next time step
00913     //
00914 
00915     stepCtrlInfoLast_ = stepCtrlInfo;
00916     ++currTimeStepIndex_;
00917     
00918   }
00919 
00920   if ( includesVerbLevel(verbLevel,Teuchos::VERB_LOW) )
00921     *out << "\nNumber of steps taken in this call to advanceStepperToTime(...) = "
00922          << (currTimeStepIndex_ - initCurrTimeStepIndex) << endl
00923          << "\nLeaving" << this->Describable::description()
00924          << "::advanceStepperToTime("<<advance_to_t<<") ...\n";
00925 
00926   return return_val;
00927   
00928 }
00929 
00930 // 
00931 // Explicit Instantiation macro
00932 //
00933 // Must be expanded from within the Rythmos namespace!
00934 //
00935 
00936 #define RYTHMOS_DEFAULT_INTEGRATOR_INSTANT(SCALAR) \
00937   \
00938   template class DefaultIntegrator< SCALAR >; \
00939   \
00940   template RCP<DefaultIntegrator< SCALAR > > \
00941   defaultIntegrator(); \
00942   \
00943   template RCP<DefaultIntegrator< SCALAR > > \
00944   defaultIntegrator( \
00945     const RCP<IntegrationControlStrategyBase< SCALAR > > &integrationControlStrategy, \
00946     const RCP<IntegrationObserverBase< SCALAR > > &integrationObserver \
00947     ); \
00948   \
00949   template RCP<DefaultIntegrator< SCALAR > > \
00950   controlledDefaultIntegrator( \
00951     const RCP<IntegrationControlStrategyBase< SCALAR > > &integrationControlStrategy \
00952     ); 
00953 
00954 
00955 } // namespace Rythmos
00956 
00957 
00958 #endif //RYTHMOS_DEFAULT_INTEGRATOR_DEF_HPP
 All Classes Functions Variables Typedefs Friends