Rythmos - Transient Integration for Differential Equations Version of the Day
Rythmos_StackedStepper.hpp
00001 //@HEADER
00002 // ***********************************************************************
00003 //
00004 //                           Rythmos Package
00005 //                 Copyright (2009) 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., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
00023 // USA
00024 // Questions? Contact Todd S. Coffey (tscoffe@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 //@HEADER
00028 
00029 #ifndef RYTHMOS_STACKED_STEPPER_HPP
00030 #define RYTHMOS_STACKED_STEPPER_HPP
00031 
00032 
00033 #include "Rythmos_StepperBase.hpp"
00034 #include "Thyra_ModelEvaluatorHelpers.hpp"
00035 #include "Teuchos_ParameterListAcceptorDefaultBase.hpp"
00036 #include "Teuchos_VerboseObjectParameterListHelpers.hpp"
00037 #include "Teuchos_Assert.hpp"
00038 #include "Teuchos_as.hpp"
00039 
00040 // Includes for ForwardSensitivityStackedStepper
00041 #include "Rythmos_ForwardSensitivityModelEvaluatorBase.hpp"
00042 
00043 namespace Rythmos {
00044 
00045 template<class Scalar>
00046 class StackedStepperStepStrategyBase
00047 {
00048   public:
00049     virtual ~StackedStepperStepStrategyBase() {}
00050     virtual void setupNextStepper(
00051         int i, 
00052         Array<RCP<StepperBase<Scalar> > > &stepperArray
00053         ) = 0;
00054     virtual Scalar evaluateStep(
00055         Scalar local_dt_taken, 
00056         int i, 
00057         Array<RCP<StepperBase<Scalar> > > &stepperArray
00058         ) = 0;
00059 };
00060 
00061 template<class Scalar>
00062 class DefaultStackedStepperStepStrategy
00063  : virtual public StackedStepperStepStrategyBase<Scalar>
00064 {
00065   public:
00066     DefaultStackedStepperStepStrategy();
00067     virtual ~DefaultStackedStepperStepStrategy();
00068     void setupNextStepper(
00069         int i, 
00070         Array<RCP<StepperBase<Scalar> > > &stepperArray
00071         );
00072     Scalar evaluateStep(
00073         Scalar local_dt_taken, 
00074         int i, 
00075         Array<RCP<StepperBase<Scalar> > > &stepperArray
00076         );
00077   private:
00078     Scalar dt_taken_;
00079 };
00080 
00081 // Nonmember constructor declaration
00082 template<class Scalar>
00083 RCP<DefaultStackedStepperStepStrategy<Scalar> > 
00084 defaultStackedStepperStepStrategy();
00085 
00086 // Nonmember constructor definition
00087 template<class Scalar>
00088 RCP<DefaultStackedStepperStepStrategy<Scalar> > 
00089 defaultStackedStepperStepStrategy()
00090 {
00091   return Teuchos::rcp(new DefaultStackedStepperStepStrategy<Scalar>());
00092 }
00093 
00094 template<class Scalar>
00095 DefaultStackedStepperStepStrategy<Scalar>::DefaultStackedStepperStepStrategy()
00096 {
00097   typedef Teuchos::ScalarTraits<Scalar> ST;
00098   dt_taken_ = ST::nan();
00099 }
00100 
00101 
00102 template<class Scalar>
00103 DefaultStackedStepperStepStrategy<Scalar>::~DefaultStackedStepperStepStrategy()
00104 {
00105 }
00106 
00107 
00108 template<class Scalar>
00109 void DefaultStackedStepperStepStrategy<Scalar>::setupNextStepper(
00110         int i, 
00111         Array<RCP<StepperBase<Scalar> > > &stepperArray
00112         )
00113 {
00114   if (i > 0) {
00115     RCP<StepperBase<Scalar> > baseStepper = stepperArray[0];
00116     stepperArray[i]->setStepControlData(*baseStepper);
00117   }
00118 }
00119 
00120 
00121 template<class Scalar>
00122 Scalar DefaultStackedStepperStepStrategy<Scalar>::evaluateStep(
00123         Scalar local_dt_taken, 
00124         int i, 
00125         Array<RCP<StepperBase<Scalar> > > &stepperArray
00126         )
00127 {
00128   if (i == 0) {
00129     dt_taken_ = local_dt_taken;
00130   }
00131   else {
00132     TEUCHOS_TEST_FOR_EXCEPTION( local_dt_taken != dt_taken_, std::logic_error,
00133         "Error!  sub-stepper["<<i<<"] did not take the same "
00134         "size step as the base stepper!"
00135         );
00136   }
00137   return dt_taken_;
00138 }
00139 
00140 
00141 
00142 template<class Scalar>
00143 class ForwardSensitivityStackedStepperStepStrategy
00144  : virtual public StackedStepperStepStrategyBase<Scalar>
00145 {
00146   public:
00147     ForwardSensitivityStackedStepperStepStrategy();
00148     virtual ~ForwardSensitivityStackedStepperStepStrategy();
00149     void setupNextStepper(
00150         int i, 
00151         Array<RCP<StepperBase<Scalar> > > &stepperArray
00152         );
00153     Scalar evaluateStep(
00154         Scalar local_dt_taken, 
00155         int i, 
00156         Array<RCP<StepperBase<Scalar> > > &stepperArray
00157         );
00158     void setStateBasePoint( 
00159         const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint
00160         );
00161   private:
00162     Scalar dt_taken_;
00163     Thyra::ModelEvaluatorBase::InArgs<Scalar> stateBasePoint_;
00164 };
00165 
00166 // Nonmember constructor declaration
00167 template<class Scalar>
00168 RCP<ForwardSensitivityStackedStepperStepStrategy<Scalar> >
00169 forwardSensitivityStackedStepperStepStrategy();
00170 
00171 // Nonmember constructor definition
00172 template<class Scalar>
00173 RCP<ForwardSensitivityStackedStepperStepStrategy<Scalar> >
00174 forwardSensitivityStackedStepperStepStrategy()
00175 {
00176   return Teuchos::rcp(new ForwardSensitivityStackedStepperStepStrategy<Scalar>());
00177 }
00178 
00179 template<class Scalar>
00180 ForwardSensitivityStackedStepperStepStrategy<Scalar>::ForwardSensitivityStackedStepperStepStrategy()
00181 {
00182   typedef Teuchos::ScalarTraits<Scalar> ST;
00183   dt_taken_ = ST::nan();
00184 }
00185 
00186 
00187 template<class Scalar>
00188 ForwardSensitivityStackedStepperStepStrategy<Scalar>::~ForwardSensitivityStackedStepperStepStrategy()
00189 {
00190 }
00191 
00192 
00193 template<class Scalar>
00194 void ForwardSensitivityStackedStepperStepStrategy<Scalar>::setupNextStepper(
00195         int i, 
00196         Array<RCP<StepperBase<Scalar> > > &stepperArray
00197         )
00198 {
00199   typedef Thyra::ModelEvaluatorBase MEB;
00200   if (i > 0) {
00201     RCP<StepperBase<Scalar> > baseStepper = stepperArray[0];
00202     RCP<Thyra::ModelEvaluator<Scalar> > modelI = 
00203       // TODO:  09/09/09 tscoffe:  This should be replaced by
00204       // getNonconstModel() when it is implemented correctly. 
00205       Teuchos::rcp_const_cast<Thyra::ModelEvaluator<Scalar> >(
00206           stepperArray[i]->getModel()
00207           );
00208     RCP<ForwardSensitivityModelEvaluatorBase<Scalar> > fwdSensME = 
00209       Teuchos::rcp_dynamic_cast<ForwardSensitivityModelEvaluatorBase<Scalar> > (
00210           modelI,false
00211           );
00212     if (Teuchos::nonnull(fwdSensME)) {
00213       bool forceUpToDateW = true;
00214       fwdSensME->initializePointState( Teuchos::inOutArg(*baseStepper), forceUpToDateW);
00215     }
00216     stepperArray[i]->setStepControlData(*baseStepper);
00217   }
00218 }
00219 
00220 
00221 template<class Scalar>
00222 Scalar ForwardSensitivityStackedStepperStepStrategy<Scalar>::evaluateStep(
00223         Scalar local_dt_taken, 
00224         int i, 
00225         Array<RCP<StepperBase<Scalar> > > &stepperArray
00226         )
00227 {
00228   if (i == 0) {
00229     dt_taken_ = local_dt_taken;
00230   }
00231   else {
00232     TEUCHOS_TEST_FOR_EXCEPTION( local_dt_taken != dt_taken_, std::logic_error,
00233         "Error!  sub-stepper["<<i<<"] did not take the same "
00234         "size step as the base stepper!"
00235         );
00236   }
00237   return dt_taken_;
00238 }
00239 
00240 
00241 template<class Scalar>
00242 void ForwardSensitivityStackedStepperStepStrategy<Scalar>::setStateBasePoint( 
00243     const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint
00244     )
00245 {
00246   stateBasePoint_ = stateBasePoint;
00247 }
00248 
00249 
00250 template<class Scalar> 
00251 class StackedStepper
00252   : virtual public StepperBase<Scalar>,
00253     virtual public Teuchos::ParameterListAcceptorDefaultBase
00254 {
00255 public:
00256 
00258   typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType ScalarMag;
00259 
00262 
00264   StackedStepper();
00265 
00268   void addStepper(const RCP<StepperBase<Scalar> >& stepper);
00269 
00272   ArrayView<const RCP<StepperBase<Scalar> > > getNonconstSteppers() const;
00273 
00276   void setStackedStepperStepControlStrategy( 
00277       const RCP<StackedStepperStepStrategyBase<Scalar> >& stepControl
00278       );
00279 
00282   RCP<const StackedStepperStepStrategyBase<Scalar> > 
00283     getStackedStepperStepCongrolStrategy() const;
00284 
00286 
00289 
00291   void setParameterList(RCP<Teuchos::ParameterList> const& paramList);
00293   RCP<const Teuchos::ParameterList> getValidParameters() const;
00294 
00296 
00299 
00301   bool acceptsModel() const;
00302 
00304   void setModel(
00305     const RCP<const Thyra::ModelEvaluator<Scalar> >& model
00306     );
00307 
00309   void setNonconstModel(
00310     const RCP<Thyra::ModelEvaluator<Scalar> >& model
00311     );
00312 
00319   RCP<const Thyra::ModelEvaluator<Scalar> > getModel() const;
00320 
00322   RCP<Thyra::ModelEvaluator<Scalar> > getNonconstModel();
00323 
00336   void setInitialCondition(
00337     const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stacked_ic
00338     );
00339 
00340   /* \brief Get the initial condigion
00341    */
00342   Thyra::ModelEvaluatorBase::InArgs<Scalar> getInitialCondition() const;
00343 
00345   Scalar takeStep( Scalar dt, StepSizeType stepType );
00346 
00348   const StepStatus<Scalar> getStepStatus() const;
00349 
00351 
00354 
00361   RCP<const Thyra::VectorSpaceBase<Scalar> >
00362   get_x_space() const;
00363 
00365   void addPoints(
00366     const Array<Scalar>& time_vec,
00367     const Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >& x_vec,
00368     const Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >& xdot_vec
00369     );
00370 
00372   TimeRange<Scalar> getTimeRange() const;
00373 
00375   void getPoints(
00376     const Array<Scalar>& time_vec,
00377     Array<RCP<const Thyra::VectorBase<Scalar> > >* x_vec,
00378     Array<RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec,
00379     Array<ScalarMag>* accuracy_vec
00380     ) const;
00381 
00383   void getNodes(Array<Scalar>* time_vec) const;
00384 
00386   void removeNodes(Array<Scalar>& time_vec);
00387 
00389   int getOrder() const;
00390 
00392 
00393 private:
00394 
00395   mutable bool isInitialized_;
00396 
00397   Array<RCP<StepperBase<Scalar> > > stepperArray_;
00398 
00399   mutable Array<RCP<const Thyra::VectorSpaceBase<Scalar> > > xSpaceArray_;
00400   //Array<RCP<const Thyra::VectorSpaceBase<Scalar> > > fSpaceArray_;
00401 
00402   mutable Teuchos::RCP<const Thyra::ProductVectorSpaceBase<Scalar> > 
00403     x_bar_space_;
00404   //Teuchos::RCP<const Thyra::ProductVectorSpaceBase<Scalar> > 
00405   //  f_bar_space_;
00406 
00407   RCP<StackedStepperStepStrategyBase<Scalar> > stackedStepperStepStrategy_;
00408   
00409   // Private functions:
00410   void setupSpaces_() const;
00411 
00412 };
00413 
00414 
00419 template<class Scalar> 
00420 inline
00421 RCP<StackedStepper<Scalar> >
00422 stackedStepper()
00423 {
00424   return Teuchos::rcp(new StackedStepper<Scalar>());
00425 }
00426 
00427 
00428 // Constructors, Intializers, Misc.
00429 
00430 
00431 template<class Scalar> 
00432 StackedStepper<Scalar>::StackedStepper()
00433   : isInitialized_(false)
00434 {}
00435 
00436 template<class Scalar> 
00437 void StackedStepper<Scalar>::setupSpaces_() const
00438 {
00439   using Teuchos::as;
00440   if (!isInitialized_) {
00441     for (int i=0 ; i < as<int>(stepperArray_.size()) ; ++i) {
00442       xSpaceArray_.push_back(stepperArray_[i]->get_x_space());
00443       //fSpaceArray_.push_back(stepperArray_[i]->get_f_space());
00444     }
00445 
00446     x_bar_space_ = Thyra::productVectorSpace<Scalar>(xSpaceArray_);
00447     //f_bar_space_ = Thyra::productVectorSpace<Scalar>(fSpaceArray_);
00448 
00449     isInitialized_ = true;
00450   }
00451 }
00452 
00453 template<class Scalar> 
00454 void StackedStepper<Scalar>::addStepper(
00455     const RCP<StepperBase<Scalar> >& stepper
00456     )
00457 {
00458   TEUCHOS_ASSERT(!is_null(stepper));
00459   stepperArray_.push_back(stepper);
00460   isInitialized_ = false;
00461 }
00462 
00463 
00464 
00465 template<class Scalar> 
00466 ArrayView<const RCP<StepperBase<Scalar> > >
00467 StackedStepper<Scalar>::getNonconstSteppers() const
00468 {
00469   return stepperArray_();
00470 }
00471 
00472 
00473 template<class Scalar> 
00474 void StackedStepper<Scalar>::setStackedStepperStepControlStrategy( 
00475     const RCP<StackedStepperStepStrategyBase<Scalar> >& stepControl
00476     )
00477 {
00478   TEUCHOS_ASSERT( Teuchos::nonnull(stepControl) );
00479   stackedStepperStepStrategy_ = stepControl;
00480 }
00481 
00482 
00483 template<class Scalar> 
00484 RCP<const StackedStepperStepStrategyBase<Scalar> > 
00485 StackedStepper<Scalar>::getStackedStepperStepCongrolStrategy() const
00486 {
00487   return stackedStepperStepStrategy_;
00488 }
00489 
00490 
00491 template<class Scalar> 
00492 void StackedStepper<Scalar>::setParameterList(
00493   RCP<Teuchos::ParameterList> const& paramList
00494   )
00495 {
00496   TEUCHOS_TEST_FOR_EXCEPT(is_null(paramList));
00497   paramList->validateParameters(*getValidParameters());
00498   this->setMyParamList(paramList);
00499   Teuchos::readVerboseObjectSublist(&*paramList,this);
00500 }
00501 
00502 
00503 template<class Scalar> 
00504 RCP<const Teuchos::ParameterList>
00505 StackedStepper<Scalar>::getValidParameters() const
00506 {
00507   static RCP<const ParameterList> validPL;
00508   if (is_null(validPL) ) {
00509     RCP<ParameterList> pl = Teuchos::parameterList();
00510     Teuchos::setupVerboseObjectSublist(&*pl);
00511     validPL = pl;
00512   }
00513   return validPL;
00514 }
00515 
00516 
00517 // Overridden from StepperBase
00518 
00519 template<class Scalar> 
00520 bool StackedStepper<Scalar>::acceptsModel() const
00521 {
00522   return false;
00523 }
00524 
00525 template<class Scalar> 
00526 void StackedStepper<Scalar>::setModel(
00527   const RCP<const Thyra::ModelEvaluator<Scalar> >& model
00528   )
00529 {
00530   TEUCHOS_TEST_FOR_EXCEPT_MSG( true,
00531     "Error, this stepper subclass does not accept a model"
00532     " as defined by the StepperBase interface!");
00533 }
00534 
00535 
00536 template<class Scalar> 
00537 void StackedStepper<Scalar>::setNonconstModel(
00538   const RCP<Thyra::ModelEvaluator<Scalar> >& model
00539   )
00540 {
00541   TEUCHOS_TEST_FOR_EXCEPT_MSG( true,
00542     "Error, this stepper subclass does not accept a model"
00543     " as defined by the StepperBase interface!");
00544 }
00545 
00546 
00547 template<class Scalar> 
00548 RCP<const Thyra::ModelEvaluator<Scalar> >
00549 StackedStepper<Scalar>::getModel() const
00550 {
00551   TEUCHOS_TEST_FOR_EXCEPT(true);
00552   return Teuchos::null;
00553 }
00554 
00555 
00556 template<class Scalar> 
00557 RCP<Thyra::ModelEvaluator<Scalar> >
00558 StackedStepper<Scalar>::getNonconstModel() 
00559 {
00560   TEUCHOS_TEST_FOR_EXCEPT(true);
00561   return Teuchos::null;
00562 }
00563 
00564 
00565 template<class Scalar> 
00566 void StackedStepper<Scalar>::setInitialCondition(
00567   const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stacked_ic
00568   )
00569 {
00570   TEUCHOS_TEST_FOR_EXCEPT(true);
00571 }
00572  
00573 
00574 template<class Scalar> 
00575 Thyra::ModelEvaluatorBase::InArgs<Scalar> 
00576 StackedStepper<Scalar>::getInitialCondition() const
00577 {
00578   TEUCHOS_TEST_FOR_EXCEPT(true);
00579   Thyra::ModelEvaluatorBase::InArgs<Scalar> inArgs;
00580   return inArgs;
00581 }
00582  
00583 
00584 template<class Scalar> 
00585 Scalar
00586 StackedStepper<Scalar>::takeStep(
00587   Scalar dt, StepSizeType stepType
00588   )
00589 {
00590   typedef Teuchos::ScalarTraits<Scalar> ST;
00591 
00592   TEUCHOS_TEST_FOR_EXCEPTION( stepperArray_.size() == 0, std::logic_error,
00593      "Error!  Rythmos::StackedStepper::takeStep: "
00594      "addStepper must be called at least once before takeStep!" 
00595      );
00596   this->setupSpaces_();
00597 
00598   RYTHMOS_FUNC_TIME_MONITOR("Rythmos:StackedStepper::takeStep");
00599   
00600   if (Teuchos::is_null(stackedStepperStepStrategy_)) {
00601     stackedStepperStepStrategy_ = defaultStackedStepperStepStrategy<Scalar>();
00602   }
00603 
00604   Scalar dt_taken = Scalar(-ST::one());
00605   for (int i=0 ; i<Teuchos::as<int>(stepperArray_.size()) ; ++i) {
00606     stackedStepperStepStrategy_->setupNextStepper(i,stepperArray_);
00607     Scalar local_dt_taken = stepperArray_[i]->takeStep(dt,stepType);
00608     dt_taken = stackedStepperStepStrategy_->evaluateStep(
00609           local_dt_taken,
00610           i,
00611           stepperArray_
00612           );
00613   }
00615   //RCP<StepperBase<Scalar> > baseStepper = stepperArray_[0];
00616   //Scalar dt_taken = baseStepper->takeStep(dt,stepType);
00617   //for (int i=1 ; i<Teuchos::as<int>(stepperArray_.size()) ; ++i) {
00618   //  // 07/27/09 tscoffe:  This line should be handled by a strategy object
00619   //  stepperArray_[i]->setStepControlData(*baseStepper);
00620   //  Scalar local_dt_taken = stepperArray_[i]->takeStep(dt,stepType);
00621   //  TEUCHOS_TEST_FOR_EXCEPTION( local_dt_taken != dt_taken, std::logic_error,
00622   //      "Error!  sub-stepper["<<i<<"] did not take the same "
00623   //      "size step as the base stepper!"
00624   //      );
00625   //}
00626   return dt_taken;
00627 }
00628 
00629 
00630 template<class Scalar> 
00631 const StepStatus<Scalar>
00632 StackedStepper<Scalar>::getStepStatus() const
00633 {
00634   TEUCHOS_TEST_FOR_EXCEPT(true);
00635   const StepStatus<Scalar> stepStatus;
00636   return stepStatus;
00637 
00638 }
00639 
00640 
00641 // Overridden from InterpolationBufferBase
00642 
00643 
00644 template<class Scalar> 
00645 RCP<const Thyra::VectorSpaceBase<Scalar> >
00646 StackedStepper<Scalar>::get_x_space() const
00647 {
00648   this->setupSpaces_();
00649   TEUCHOS_TEST_FOR_EXCEPTION( is_null(x_bar_space_), std::logic_error,
00650       "Rythmos::StackedStepper::get_x_space():  "
00651       "addStepper must be called at least once before get_x_space()!"
00652       );
00653   return(x_bar_space_);
00654 }
00655 
00656 
00657 template<class Scalar> 
00658 void StackedStepper<Scalar>::addPoints(
00659   const Array<Scalar>& time_vec,
00660   const Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >& x_vec,
00661   const Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >& xdot_vec
00662   )
00663 {
00664   TEUCHOS_TEST_FOR_EXCEPT(
00665       "Not implemented addPoints(...) yet but we could if we wanted!"
00666       );
00667 }
00668 
00669 
00670 template<class Scalar> 
00671 TimeRange<Scalar>
00672 StackedStepper<Scalar>::getTimeRange() const
00673 {
00674   TEUCHOS_TEST_FOR_EXCEPT(true);
00675   TimeRange<Scalar> tr;
00676   return tr;
00677 }
00678 
00679 
00680 template<class Scalar> 
00681 void StackedStepper<Scalar>::getPoints(
00682   const Array<Scalar>& time_vec,
00683   Array<RCP<const Thyra::VectorBase<Scalar> > >* x_bar_vec,
00684   Array<RCP<const Thyra::VectorBase<Scalar> > >* x_bar_dot_vec,
00685   Array<ScalarMag>* accuracy_vec
00686   ) const
00687 {
00688   using Teuchos::as;
00689   this->setupSpaces_();
00690   TEUCHOS_TEST_FOR_EXCEPTION( stepperArray_.size() == 0, std::logic_error,
00691       "Rythmos::StackedStepper:getPoints:  Error!  "
00692       "addStepper must be called at least once before getPoints!"
00693       );
00694   int numSteppers = as<int>(stepperArray_.size());
00695   int numTimePoints = as<int>(time_vec.size());
00696 
00697   // Set up local x_bar_vec for filling:
00698   Array<RCP<Thyra::VectorBase<Scalar> > > x_bar_vec_local;
00699   if (x_bar_vec != NULL) {
00700     x_bar_vec_local.clear();
00701     x_bar_vec->clear();
00702     for (int n = 0 ; n < numTimePoints ; ++n) {
00703       x_bar_vec_local.push_back(Thyra::createMember(this->get_x_space()));
00704       x_bar_vec->push_back(x_bar_vec_local[n]);
00705     }
00706   }
00707 
00708   // Set up local x_bar_dot_vec for filling:
00709   Array<RCP<Thyra::VectorBase<Scalar> > > x_bar_dot_vec_local;
00710   if (x_bar_dot_vec != NULL) {
00711     x_bar_dot_vec_local.clear();
00712     x_bar_dot_vec->clear();
00713     for (int n = 0 ; n < numTimePoints ; ++n) {
00714       x_bar_dot_vec_local.push_back(Thyra::createMember(this->get_x_space()));
00715       x_bar_dot_vec->push_back(x_bar_dot_vec_local[n]);
00716     }
00717   }
00718 
00719   // Set up accuracy_vec
00720   if (accuracy_vec) {
00721     accuracy_vec->clear();
00722     accuracy_vec->resize(numTimePoints);
00723   }
00724 
00725   for (int i=0 ; i < numSteppers; ++i) {
00726     Array<RCP<const Thyra::VectorBase<Scalar> > > sub_x_bar_vec;
00727     Array<RCP<const Thyra::VectorBase<Scalar> > > sub_x_bar_dot_vec;
00728     Array<ScalarMag> sub_accuracy_vec;
00729     stepperArray_[i]->getPoints(
00730         time_vec,
00731         x_bar_vec ? &sub_x_bar_vec : 0,
00732         x_bar_dot_vec ? &sub_x_bar_dot_vec : 0,
00733         accuracy_vec ? &sub_accuracy_vec: 0
00734         );
00735     // Copy vectors into the sub-blocks of the 
00736     // x_bar_vec, x_bar_dot_vec, and accuracy_vec vectors.
00737     for (int n=0; n < numTimePoints ; ++n ) {
00738       if (x_bar_vec) {
00739         RCP<Thyra::ProductVectorBase<Scalar> > x_bar_pv = 
00740           Teuchos::rcp_dynamic_cast<Thyra::ProductVectorBase<Scalar> >(
00741               x_bar_vec_local[n],
00742               true
00743               );
00744         RCP<Thyra::VectorBase<Scalar> > xi = 
00745           x_bar_pv->getNonconstVectorBlock(i);
00746         V_V(Teuchos::outArg(*xi),*sub_x_bar_vec[n]);
00747       }
00748       if (x_bar_dot_vec) {
00749         RCP<Thyra::ProductVectorBase<Scalar> > x_bar_dot_pv = 
00750           Teuchos::rcp_dynamic_cast<Thyra::ProductVectorBase<Scalar> >(
00751               x_bar_dot_vec_local[n],
00752               true
00753               );
00754         RCP<Thyra::VectorBase<Scalar> > xdoti = 
00755           x_bar_dot_pv->getNonconstVectorBlock(i);
00756         V_V(Teuchos::outArg(*xdoti),*sub_x_bar_dot_vec[n]);
00757       }
00758       // Just copy the accuracy from the first stepper:
00759       if ( (i == 0) && accuracy_vec) {
00760           (*accuracy_vec)[n] = sub_accuracy_vec[n];
00761       }
00762     }
00763   }
00764 }
00765 
00766 
00767 template<class Scalar>
00768 void StackedStepper<Scalar>::getNodes(
00769   Array<Scalar>* time_vec
00770   ) const
00771 {
00772   TEUCHOS_TEST_FOR_EXCEPT(true);
00773 }
00774 
00775 
00776 template<class Scalar> 
00777 void StackedStepper<Scalar>::removeNodes(
00778   Array<Scalar>& time_vec
00779   )
00780 {
00781   TEUCHOS_TEST_FOR_EXCEPT("Not implemented yet but we can!");
00782 }
00783 
00784 
00785 template<class Scalar> 
00786 int StackedStepper<Scalar>::getOrder() const
00787 {
00788   TEUCHOS_TEST_FOR_EXCEPT(true);
00789   return -1;
00790 }
00791 
00792 
00793 // private
00794 
00795 
00796 } // namespace Rythmos
00797 
00798 
00799 #endif //RYTHMOS_STACKED_STEPPER_HPP
 All Classes Functions Variables Typedefs Friends