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

Generated on Wed May 12 21:25:42 2010 for Rythmos - Transient Integration for Differential Equations by  doxygen 1.4.7