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 
00253   using Teuchos::null;
00254   RCP<DefaultIntegrator<Scalar> >
00255     newIntegrator = Teuchos::rcp(new DefaultIntegrator<Scalar>());
00256   // Only copy control information, not the stepper or the model it contains!
00257   newIntegrator->stepper_ = Teuchos::null;
00258   const RCP<const ParameterList> paramList = this->getParameterList();
00259   if (!is_null(paramList)) {
00260     newIntegrator->setParameterList(Teuchos::parameterList(*paramList));
00261   }
00262   if (!is_null(integrationControlStrategy_)) {
00263     newIntegrator->integrationControlStrategy_ =
00264       integrationControlStrategy_->cloneIntegrationControlStrategy().assert_not_null();
00265   }
00266   if (!is_null(integrationObserver_)) {
00267     newIntegrator->integrationObserver_ =
00268       integrationObserver_->cloneIntegrationObserver().assert_not_null();
00269   }
00270   if (!is_null(trailingInterpBuffer_)) {
00271     // ToDo: implement the clone!
00272     newIntegrator->trailingInterpBuffer_ = null;
00273     //newIntegrator->trailingInterpBuffer_ =
00274     //  trailingInterpBuffer_->cloneInterploationBuffer().assert_not_null();
00275   }
00276   if (!is_null(interpBufferAppender_)) {
00277     // ToDo: implement the clone!
00278     newIntegrator->interpBufferAppender_ = null;
00279     //newIntegrator->interpBufferAppender_ =
00280     //  interpBufferAppender_->cloneInterpolationBufferAppender().assert_not_null();
00281   }
00282   return newIntegrator;
00283 }
00284 
00285 
00286 template<class Scalar> 
00287 void DefaultIntegrator<Scalar>::setStepper(
00288   const RCP<StepperBase<Scalar> > &stepper,
00289   const Scalar &finalTime,
00290   const bool landOnFinalTime
00291   )
00292 {
00293   typedef Teuchos::ScalarTraits<Scalar> ST;
00294   TEST_FOR_EXCEPT(is_null(stepper));
00295   TEST_FOR_EXCEPT( finalTime <= stepper->getTimeRange().lower() );
00296   TEUCHOS_ASSERT( stepper->getTimeRange().length() == ST::zero() );
00297   // 2007/07/25: rabartl: ToDo: Validate state of the stepper!
00298   stepper_ = stepper;
00299   integrationTimeDomain_ = timeRange(stepper_->getTimeRange().lower(), finalTime);
00300   landOnFinalTime_ = landOnFinalTime;
00301   currTimeStepIndex_ = 0;
00302   stepCtrlInfoLast_ = StepControlInfo<Scalar>();
00303   if (!is_null(integrationControlStrategy_))
00304     integrationControlStrategy_->resetIntegrationControlStrategy(
00305       integrationTimeDomain_
00306       );
00307   if (!is_null(integrationObserver_))
00308     integrationObserver_->resetIntegrationObserver(
00309       integrationTimeDomain_
00310       );
00311 }
00312 
00313 
00314 template<class Scalar>
00315 RCP<StepperBase<Scalar> > DefaultIntegrator<Scalar>::unSetStepper()
00316 {
00317   RCP<StepperBase<Scalar> > stepper_temp = stepper_;
00318   stepper_ = Teuchos::null;
00319   return(stepper_temp);
00320 }
00321 
00322 
00323 template<class Scalar>
00324 RCP<const StepperBase<Scalar> > DefaultIntegrator<Scalar>::getStepper() const
00325 {
00326   return(stepper_);
00327 }
00328 
00329 
00330 template<class Scalar>
00331 RCP<StepperBase<Scalar> > DefaultIntegrator<Scalar>::getNonconstStepper() const
00332 {
00333   return(stepper_);
00334 }
00335 
00336 
00337 template<class Scalar>
00338 void DefaultIntegrator<Scalar>::setTrailingInterpolationBuffer(
00339   const RCP<InterpolationBufferBase<Scalar> > &trailingInterpBuffer
00340   )
00341 {
00342   trailingInterpBuffer_ = trailingInterpBuffer.assert_not_null();
00343 }
00344 
00345 
00346 template<class Scalar>
00347 RCP<InterpolationBufferBase<Scalar> >
00348 DefaultIntegrator<Scalar>::getNonconstTrailingInterpolationBuffer()
00349 {
00350   return trailingInterpBuffer_;
00351 }
00352 
00353 
00354 template<class Scalar>
00355 RCP<const InterpolationBufferBase<Scalar> >
00356 DefaultIntegrator<Scalar>::getTrailingInterpolationBuffer() const
00357 {
00358   return trailingInterpBuffer_;
00359 }
00360 
00361 template<class Scalar>
00362 RCP<InterpolationBufferBase<Scalar> >
00363 DefaultIntegrator<Scalar>::unSetTrailingInterpolationBuffer()
00364 {
00365   RCP<InterpolationBufferBase<Scalar> > trailingInterpBuffer;
00366   std::swap( trailingInterpBuffer, trailingInterpBuffer_ );
00367   return trailingInterpBuffer;
00368 }
00369 
00370 
00371 template<class Scalar>
00372 void DefaultIntegrator<Scalar>::getFwdPoints(
00373   const Array<Scalar>& time_vec,
00374   Array<RCP<const Thyra::VectorBase<Scalar> > >* x_vec,
00375   Array<RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec,
00376   Array<ScalarMag>* accuracy_vec
00377   )
00378 {
00379 
00380 #ifdef ENABLE_RYTHMOS_TIMERS
00381   TEUCHOS_FUNC_TIME_MONITOR_DIFF("Rythmos:DefaultIntegrator::getFwdPoints",
00382     TopLevel);
00383 #endif
00384 
00385   using Teuchos::incrVerbLevel;
00386 #ifndef _MSC_VER
00387   using Teuchos::Describable;
00388 #endif
00389   typedef Teuchos::ScalarTraits<Scalar> ST;
00390   typedef InterpolationBufferBase<Scalar> IBB;
00391   typedef Teuchos::VerboseObjectTempState<IBB> VOTSIBB;
00392 
00393   finalizeSetup();
00394 
00395   RCP<Teuchos::FancyOStream> out = this->getOStream();
00396   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00397   Teuchos::OSTab tab(out);
00398   VOTSIBB stepper_outputTempState(stepper_,out,incrVerbLevel(verbLevel,-1));
00399 
00400   if ( includesVerbLevel(verbLevel,Teuchos::VERB_LOW) )
00401     *out << "\nEntering " << this->Describable::description() << "::getFwdPoints(...) ...\n"
00402          << "\nStepper: " << Teuchos::describe(*stepper_,verbLevel);
00403 
00404   if ( includesVerbLevel(verbLevel,Teuchos::VERB_MEDIUM) )
00405     *out << "\nRequested time points: " << Teuchos::toString(time_vec) << "\n";
00406 
00407   //
00408   // 0) Initial setup
00409   //
00410 
00411   const int numTimePoints = time_vec.size();
00412 
00413   // Assert preconditions
00414   assertTimePointsAreSorted(time_vec);
00415   TEST_FOR_EXCEPT(accuracy_vec!=0); // ToDo: Remove accuracy_vec!
00416 
00417   // Resize the storage for the output arrays
00418   if (x_vec)
00419     x_vec->resize(numTimePoints);
00420   if (xdot_vec)
00421     xdot_vec->resize(numTimePoints);
00422 
00423   // This int records the next time point offset in time_vec[timePointIndex]
00424   // that needs to be handled.  This gets updated as the time points are
00425   // filled below.
00426   int nextTimePointIndex = 0;
00427   
00428   assertNoTimePointsBeforeCurrentTimeRange(*this,time_vec,nextTimePointIndex);
00429 
00430   //
00431   // 1) First, get all time points that fall within the current time range
00432   //
00433 
00434   {
00435 #ifdef ENABLE_RYTHMOS_TIMERS
00436     TEUCHOS_FUNC_TIME_MONITOR("Rythmos:DefaultIntegrator::getFwdPoints: getPoints");
00437 #endif
00438     // 2007/10/05: rabartl: ToDo: Get points from trailingInterpBuffer_ first!
00439     getCurrentPoints(*stepper_,time_vec,x_vec,xdot_vec,&nextTimePointIndex);
00440   }
00441 
00442   //
00443   // 2) Advance the stepper to satisfy time points in time_vec that fall
00444   // before the current time.
00445   //
00446 
00447   while ( nextTimePointIndex < numTimePoints ) {
00448     
00449     // Use the time stepping algorithm to step up to or past the next
00450     // requested time but not so far as to step past the point entirely.
00451     const Scalar t = time_vec[nextTimePointIndex];
00452     bool advanceStepperToTimeSucceeded = false;
00453     {
00454 #ifdef ENABLE_RYTHMOS_TIMERS
00455       TEUCHOS_FUNC_TIME_MONITOR("Rythmos:DefaultIntegrator::getFwdPoints: advanceStepperToTime");
00456 #endif
00457       advanceStepperToTimeSucceeded= advanceStepperToTime(t);
00458     }
00459     if (!advanceStepperToTimeSucceeded) {
00460       bool reachedMaxNumTimeSteps = (currTimeStepIndex_ >= maxNumTimeSteps_);
00461       if (reachedMaxNumTimeSteps) {
00462         // Break out of the while loop and attempt to exit gracefully.
00463         break;
00464       }
00465       TEST_FOR_EXCEPTION(
00466           !advanceStepperToTimeSucceeded, Exceptions::GetFwdPointsFailed,
00467           this->description() << "\n\n"
00468           "Error:  The integration failed to get to time " << t << " and only achieved\n"
00469           "getting to " << stepper_->getTimeRange().upper() << "!"
00470           );
00471     }
00472     
00473     // Extract the next set of points (perhaps just one) from the stepper
00474     {
00475 #ifdef ENABLE_RYTHMOS_TIMERS
00476       TEUCHOS_FUNC_TIME_MONITOR("Rythmos:DefaultIntegrator::getFwdPoints: getPoints (fwd)");
00477 #endif
00478       getCurrentPoints(*stepper_,time_vec,x_vec,xdot_vec,&nextTimePointIndex);
00479     }
00480     
00481   }
00482 
00483   if ( includesVerbLevel(verbLevel,Teuchos::VERB_LOW) )
00484     *out << "\nLeaving " << this->Describable::description() << "::getFwdPoints(...) ...\n";
00485   
00486 }
00487 
00488 
00489 template<class Scalar> 
00490 TimeRange<Scalar>
00491 DefaultIntegrator<Scalar>::getFwdTimeRange() const
00492 {
00493   return timeRange(
00494     stepper_->getTimeRange().lower(),
00495     integrationTimeDomain_.upper()
00496     );
00497 }
00498 
00499 
00500 // Overridden from InterpolationBufferBase
00501 
00502 
00503 template<class Scalar> 
00504 RCP<const Thyra::VectorSpaceBase<Scalar> >
00505 DefaultIntegrator<Scalar>::get_x_space() const
00506 {
00507   return stepper_->get_x_space();
00508 }
00509 
00510 
00511 template<class Scalar> 
00512 void DefaultIntegrator<Scalar>::addPoints(
00513   const Array<Scalar>& time_vec,
00514   const Array<RCP<const Thyra::VectorBase<Scalar> > >& x_vec,
00515   const Array<RCP<const Thyra::VectorBase<Scalar> > >& xdot_vec
00516   )
00517 {
00518   stepper_->addPoints(time_vec,x_vec,xdot_vec);
00519 }
00520 
00521 
00522 template<class Scalar> 
00523 void DefaultIntegrator<Scalar>::getPoints(
00524   const Array<Scalar>& time_vec,
00525   Array<RCP<const Thyra::VectorBase<Scalar> > >* x_vec,
00526   Array<RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec,
00527   Array<ScalarMag>* accuracy_vec
00528   ) const
00529 {
00530 //  if (nonnull(trailingInterpBuffer_)) {
00531 //    int nextTimePointIndex = 0;
00532 //    getCurrentPoints(*trailingInterpBuffer_, time_vec, x_vec, xdot_vec, &nextTimePointIndex);
00533 //    getCurrentPoints(*stepper_, time_vec, x_vec, xdot_vec, &nextTimePointIndex);
00534 //    TEST_FOR_EXCEPTION( nextTimePointIndex < Teuchos::as<int>(time_vec.size()),
00535 //      std::out_of_range,
00536 //      "Error, the time point time_vec["<<nextTimePointIndex<<"] = "
00537 //      << time_vec[nextTimePointIndex] << " falls outside of the time range "
00538 //      << getTimeRange() << "!"
00539 //      );
00540 //  }
00541 //  else {
00542   stepper_->getPoints(time_vec, x_vec, xdot_vec, accuracy_vec);
00543 //  }
00544 }
00545 
00546 
00547 template<class Scalar> 
00548 TimeRange<Scalar> DefaultIntegrator<Scalar>::getTimeRange() const
00549 {
00550   if (nonnull(trailingInterpBuffer_)) {
00551     return timeRange(trailingInterpBuffer_->getTimeRange().lower(),
00552       stepper_->getTimeRange().upper());
00553   }
00554   return stepper_->getTimeRange();
00555 }
00556 
00557 
00558 template<class Scalar> 
00559 void DefaultIntegrator<Scalar>::getNodes(Array<Scalar>* time_vec) const
00560 {
00561   stepper_->getNodes(time_vec);
00562 }
00563 
00564 
00565 template<class Scalar> 
00566 void DefaultIntegrator<Scalar>::removeNodes(Array<Scalar>& time_vec)
00567 {
00568   stepper_->removeNodes(time_vec);
00569 }
00570 
00571 
00572 template<class Scalar> 
00573 int DefaultIntegrator<Scalar>::getOrder() const
00574 {
00575   return stepper_->getOrder();
00576 }
00577 
00578 
00579 // private
00580 
00581 
00582 template<class Scalar> 
00583 void DefaultIntegrator<Scalar>::finalizeSetup()
00584 {
00585   if (!is_null(trailingInterpBuffer_) && is_null(interpBufferAppender_))
00586     interpBufferAppender_ = pointwiseInterpolationBufferAppender<Scalar>();
00587   // ToDo: Do other setup?
00588 }
00589 
00590 
00591 template<class Scalar> 
00592 bool DefaultIntegrator<Scalar>::advanceStepperToTime( const Scalar& advance_to_t )
00593 {
00594 
00595 #ifdef ENABLE_RYTHMOS_TIMERS
00596   TEUCHOS_FUNC_TIME_MONITOR_DIFF("Rythmos:DefaultIntegrator::advanceStepperToTime",
00597     TopLevel);
00598 #endif
00599 
00600   using std::endl;
00601   typedef std::numeric_limits<Scalar> NL;
00602   using Teuchos::incrVerbLevel; 
00603 #ifndef _MSC_VER
00604   using Teuchos::Describable;
00605 #endif
00606   using Teuchos::OSTab;
00607   typedef Teuchos::ScalarTraits<Scalar> ST;
00608 
00609   RCP<Teuchos::FancyOStream> out = this->getOStream();
00610   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00611 
00612   if (!is_null(integrationControlStrategy_)) {
00613     integrationControlStrategy_->setOStream(out);
00614     integrationControlStrategy_->setVerbLevel(incrVerbLevel(verbLevel,-1));
00615   }
00616 
00617   if (!is_null(integrationObserver_)) {
00618     integrationObserver_->setOStream(out);
00619     integrationObserver_->setVerbLevel(incrVerbLevel(verbLevel,-1));
00620   }
00621 
00622   if ( includesVerbLevel(verbLevel,Teuchos::VERB_LOW) )
00623     *out << "\nEntering " << this->Describable::description()
00624          << "::advanceStepperToTime("<<advance_to_t<<") ...\n";
00625 
00626   // Remember what timestep index we are on so we can report it later
00627   const int initCurrTimeStepIndex = currTimeStepIndex_;
00628 
00629   // Take steps until we the requested time is reached (or passed)
00630 
00631   TimeRange<Scalar> currStepperTimeRange = stepper_->getTimeRange();
00632 
00633   // Start by assume we can reach the time advance_to_t
00634   bool return_val = true;
00635   
00636   while ( !currStepperTimeRange.isInRange(advance_to_t) ) {
00637 
00638     // Halt immediatly if exceeded max iterations
00639     if (currTimeStepIndex_ >= maxNumTimeSteps_) {
00640       if ( includesVerbLevel(verbLevel,Teuchos::VERB_LOW) )
00641         *out
00642           << "\n***"
00643           << "\n*** NOTICE: currTimeStepIndex = "<<currTimeStepIndex_
00644           << " >= maxNumTimeSteps = "<<maxNumTimeSteps_<< ", halting time integration!"
00645           << "\n***\n";
00646       return_val = false;
00647       break; // Exit the loop immediately!
00648     }
00649 
00650     if ( includesVerbLevel(verbLevel,Teuchos::VERB_LOW) )
00651       *out << "\nTake step:  current_stepper_t = " << currStepperTimeRange.upper()
00652            << ", currTimeStepIndex = " << currTimeStepIndex_ << endl;
00653 
00654     //
00655     // A) Reinitialize if a hard breakpoint was reached on the last time step
00656     //
00657 
00658     if (stepCtrlInfoLast_.limitedByBreakPoint) {
00659       if ( stepCtrlInfoLast_.breakPointType == BREAK_POINT_TYPE_HARD ) {
00660 #ifdef ENABLE_RYTHMOS_TIMERS
00661         TEUCHOS_FUNC_TIME_MONITOR("Rythmos:DefaultIntegrator::restart");
00662 #endif
00663         if ( includesVerbLevel(verbLevel,Teuchos::VERB_LOW) )
00664           *out << "\nAt a hard-breakpoint, restarting time integrator ...\n";
00665         restart(&*stepper_);
00666       }
00667       else  {
00668         if ( includesVerbLevel(verbLevel,Teuchos::VERB_LOW) )
00669           *out << "\nAt a soft-breakpoint, NOT restarting time integrator ...\n";
00670       }
00671     }
00672 
00673     //
00674     // B) Get the trial step control info
00675     //
00676 
00677     StepControlInfo<Scalar> trialStepCtrlInfo;
00678     {
00679 #ifdef ENABLE_RYTHMOS_TIMERS
00680       TEUCHOS_FUNC_TIME_MONITOR("Rythmos:DefaultIntegrator::advanceStepperToTime: getStepCtrl");
00681 #endif
00682       if (!is_null(integrationControlStrategy_)) {
00683         // Let an external strategy object determine the step size and type.
00684         // Note that any breakpoint info is also related through this call.
00685         trialStepCtrlInfo = integrationControlStrategy_->getNextStepControlInfo(
00686           *stepper_, stepCtrlInfoLast_, currTimeStepIndex_
00687           );
00688       }
00689       else {
00690         // Take a variable step if we have no control strategy
00691         trialStepCtrlInfo.stepType = STEP_TYPE_VARIABLE;
00692         trialStepCtrlInfo.stepSize = NL::max();
00693       }
00694     }
00695 
00696     // Print the initial trial step
00697     if ( includesVerbLevel(verbLevel,Teuchos::VERB_MEDIUM) ) {
00698       *out << "\nTrial step:\n";
00699       OSTab tab(out);
00700       *out << trialStepCtrlInfo;
00701     }
00702 
00703     // Halt immediately if we where told to do so
00704     if (trialStepCtrlInfo.stepSize < ST::zero()) {
00705       if ( includesVerbLevel(verbLevel,Teuchos::VERB_MEDIUM) )
00706         *out
00707           << "\n***"
00708           << "\n*** NOTICE: The IntegrationControlStrategy object return stepSize < 0.0, halting time integration!"
00709           << "\n***\n";
00710       return_val = false;
00711       break; // Exit the loop immediately!
00712     }
00713 
00714     // Make sure we don't step past the final time if asked not to
00715     bool updatedTrialStepCtrlInfo = false;
00716     {
00717       const Scalar finalTime = integrationTimeDomain_.upper();
00718       if (landOnFinalTime_ && trialStepCtrlInfo.stepSize + currStepperTimeRange.upper() > finalTime) {
00719         if ( includesVerbLevel(verbLevel,Teuchos::VERB_LOW) )
00720           *out << "\nCutting trial step to avoid stepping past final time ...\n";
00721         trialStepCtrlInfo.stepSize = finalTime - currStepperTimeRange.upper();
00722         updatedTrialStepCtrlInfo = true;
00723       }
00724     }
00725     
00726     // Print the modified trial step
00727     if ( updatedTrialStepCtrlInfo
00728       && includesVerbLevel(verbLevel,Teuchos::VERB_MEDIUM) )
00729     {
00730       *out << "\nUpdated trial step:\n";
00731       OSTab tab(out);
00732       *out << trialStepCtrlInfo;
00733     }
00734 
00735     //
00736     // C) Take the step
00737     //
00738 
00739     // Print step type and size
00740     if ( includesVerbLevel(verbLevel,Teuchos::VERB_MEDIUM) ) {
00741       if (trialStepCtrlInfo.stepType == STEP_TYPE_VARIABLE)
00742         *out << "\nTaking a variable time step with max step size = "
00743              << trialStepCtrlInfo.stepSize << " ....\n";
00744       else
00745         *out << "\nTaking a fixed time step of size = "
00746              << trialStepCtrlInfo.stepSize << " ....\n";
00747     }
00748 
00749     // Take step
00750     Scalar stepSizeTaken;
00751     {
00752 #ifdef ENABLE_RYTHMOS_TIMERS
00753       TEUCHOS_FUNC_TIME_MONITOR("Rythmos:DefaultIntegrator::advanceStepperToTime: takeStep");
00754 #endif
00755       stepSizeTaken = stepper_->takeStep(
00756         trialStepCtrlInfo.stepSize, trialStepCtrlInfo.stepType
00757         );
00758     }
00759 
00760     // Validate step taken
00761     if (trialStepCtrlInfo.stepType == STEP_TYPE_VARIABLE) {
00762       TEST_FOR_EXCEPTION(
00763         stepSizeTaken < ST::zero(), std::logic_error,
00764         "Error, stepper took negative step of dt = " << stepSizeTaken << "!\n"
00765         );
00766       TEST_FOR_EXCEPTION(
00767         stepSizeTaken > trialStepCtrlInfo.stepSize, std::logic_error,
00768         "Error, stepper took step of dt = " << stepSizeTaken
00769         << " > max step size of = " << trialStepCtrlInfo.stepSize << "!\n"
00770         );
00771     }
00772     else { // STEP_TYPE_FIXED
00773       TEST_FOR_EXCEPTION(
00774         stepSizeTaken != trialStepCtrlInfo.stepSize, std::logic_error,
00775         "Error, stepper took step of dt = " << stepSizeTaken 
00776         << " when asked to take step of dt = " << trialStepCtrlInfo.stepSize << "\n"
00777         );
00778     }
00779 
00780     // Update info about this step
00781     currStepperTimeRange = stepper_->getTimeRange();
00782     const StepControlInfo<Scalar> stepCtrlInfo =
00783       stepCtrlInfoTaken(trialStepCtrlInfo,stepSizeTaken);
00784 
00785     // Print the step actually taken 
00786     if ( includesVerbLevel(verbLevel,Teuchos::VERB_MEDIUM) ) {
00787       *out << "\nStep actually taken:\n";
00788       OSTab tab(out);
00789       *out << stepCtrlInfo;
00790     }
00791 
00792     // Append the trailing interpolation buffer (if defined)
00793     if (!is_null(trailingInterpBuffer_)) {
00794       interpBufferAppender_->append(*stepper_,currStepperTimeRange,
00795         trailingInterpBuffer_.ptr() );
00796     }
00797 
00798     //
00799     // D) Output info about step
00800     //
00801 
00802     {
00803 
00804 #ifdef ENABLE_RYTHMOS_TIMERS
00805       TEUCHOS_FUNC_TIME_MONITOR("Rythmos:DefaultIntegrator::advanceStepperToTime: output");
00806 #endif
00807       
00808       // Print our own brief output
00809       if ( includesVerbLevel(verbLevel,Teuchos::VERB_MEDIUM) ) {
00810         StepStatus<Scalar> stepStatus = stepper_->getStepStatus();
00811         *out << "\nTime point reached = " << stepStatus.time << endl;
00812         *out << "\nstepStatus:\n" << stepStatus;
00813         if ( includesVerbLevel(verbLevel,Teuchos::VERB_EXTREME) ) {
00814           RCP<const Thyra::VectorBase<Scalar> >
00815             solution = stepStatus.solution,
00816             solutionDot = stepStatus.solutionDot;
00817           if (!is_null(solution))
00818             *out << "\nsolution = \n" << Teuchos::describe(*solution,verbLevel);
00819           if (!is_null(solutionDot))
00820             *out << "\nsolutionDot = \n" << Teuchos::describe(*solutionDot,verbLevel);
00821         }
00822       }
00823       
00824       // Output to the observer
00825       if (!is_null(integrationObserver_))
00826         integrationObserver_->observeCompletedTimeStep(
00827           *stepper_, stepCtrlInfo, currTimeStepIndex_
00828           );
00829 
00830     }
00831 
00832     //
00833     // E) Update info for next time step
00834     //
00835 
00836     stepCtrlInfoLast_ = stepCtrlInfo;
00837     ++currTimeStepIndex_;
00838     
00839   }
00840 
00841   if ( includesVerbLevel(verbLevel,Teuchos::VERB_LOW) )
00842     *out << "\nNumber of steps taken in this call to advanceStepperToTime(...) = "
00843          << (currTimeStepIndex_ - initCurrTimeStepIndex) << endl
00844          << "\nLeaving" << this->Describable::description()
00845          << "::advanceStepperToTime("<<advance_to_t<<") ...\n";
00846 
00847   return return_val;
00848   
00849 }
00850 
00851 // 
00852 // Explicit Instantiation macro
00853 //
00854 // Must be expanded from within the Rythmos namespace!
00855 //
00856 
00857 #define RYTHMOS_DEFAULT_INTEGRATOR_INSTANT(SCALAR) \
00858   \
00859   template class DefaultIntegrator< SCALAR >; \
00860   \
00861   template RCP<DefaultIntegrator< SCALAR > > \
00862   defaultIntegrator(); \
00863   \
00864   template RCP<DefaultIntegrator< SCALAR > > \
00865   defaultIntegrator( \
00866     const RCP<IntegrationControlStrategyBase< SCALAR > > &integrationControlStrategy, \
00867     const RCP<IntegrationObserverBase< SCALAR > > &integrationObserver \
00868     ); \
00869   \
00870   template RCP<DefaultIntegrator< SCALAR > > \
00871   controlledDefaultIntegrator( \
00872     const RCP<IntegrationControlStrategyBase< SCALAR > > &integrationControlStrategy \
00873     ); 
00874 
00875 
00876 } // namespace Rythmos
00877 
00878 
00879 #endif //RYTHMOS_DEFAULT_INTEGRATOR_DEF_HPP
 All Classes Functions Variables Typedefs Friends
Generated on Wed Apr 13 09:58:54 2011 for Rythmos - Transient Integration for Differential Equations by  doxygen 1.6.3