29#ifndef RYTHMOS_STACKED_STEPPER_HPP
30#define RYTHMOS_STACKED_STEPPER_HPP
33#include "Rythmos_StepperBase.hpp"
34#include "Thyra_ModelEvaluatorHelpers.hpp"
35#include "Teuchos_ParameterListAcceptorDefaultBase.hpp"
36#include "Teuchos_VerboseObjectParameterListHelpers.hpp"
37#include "Teuchos_Assert.hpp"
38#include "Teuchos_as.hpp"
41#include "Rythmos_ForwardSensitivityModelEvaluatorBase.hpp"
46class StackedStepperStepStrategyBase
49 virtual ~StackedStepperStepStrategyBase() {}
50 virtual void setupNextStepper(
52 Array<RCP<StepperBase<Scalar> > > &stepperArray
54 virtual Scalar evaluateStep(
55 Scalar local_dt_taken,
57 Array<RCP<StepperBase<Scalar> > > &stepperArray
62class DefaultStackedStepperStepStrategy
63 :
virtual public StackedStepperStepStrategyBase<Scalar>
66 DefaultStackedStepperStepStrategy();
67 virtual ~DefaultStackedStepperStepStrategy();
68 void setupNextStepper(
70 Array<RCP<StepperBase<Scalar> > > &stepperArray
73 Scalar local_dt_taken,
75 Array<RCP<StepperBase<Scalar> > > &stepperArray
83RCP<DefaultStackedStepperStepStrategy<Scalar> >
84defaultStackedStepperStepStrategy();
88RCP<DefaultStackedStepperStepStrategy<Scalar> >
89defaultStackedStepperStepStrategy()
91 return Teuchos::rcp(
new DefaultStackedStepperStepStrategy<Scalar>());
95DefaultStackedStepperStepStrategy<Scalar>::DefaultStackedStepperStepStrategy()
97 typedef Teuchos::ScalarTraits<Scalar> ST;
98 dt_taken_ = ST::nan();
102template<
class Scalar>
103DefaultStackedStepperStepStrategy<Scalar>::~DefaultStackedStepperStepStrategy()
108template<
class Scalar>
109void DefaultStackedStepperStepStrategy<Scalar>::setupNextStepper(
111 Array<RCP<StepperBase<Scalar> > > &stepperArray
115 RCP<StepperBase<Scalar> > baseStepper = stepperArray[0];
116 stepperArray[i]->setStepControlData(*baseStepper);
121template<
class Scalar>
122Scalar DefaultStackedStepperStepStrategy<Scalar>::evaluateStep(
123 Scalar local_dt_taken,
125 Array<RCP<StepperBase<Scalar> > > &stepperArray
129 dt_taken_ = local_dt_taken;
132 TEUCHOS_TEST_FOR_EXCEPTION( local_dt_taken != dt_taken_, std::logic_error,
133 "Error! sub-stepper["<<i<<
"] did not take the same "
134 "size step as the base stepper!"
142template<
class Scalar>
143class ForwardSensitivityStackedStepperStepStrategy
144 :
virtual public StackedStepperStepStrategyBase<Scalar>
147 ForwardSensitivityStackedStepperStepStrategy();
148 virtual ~ForwardSensitivityStackedStepperStepStrategy();
149 void setupNextStepper(
151 Array<RCP<StepperBase<Scalar> > > &stepperArray
154 Scalar local_dt_taken,
156 Array<RCP<StepperBase<Scalar> > > &stepperArray
158 void setStateBasePoint(
159 const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint
163 Thyra::ModelEvaluatorBase::InArgs<Scalar> stateBasePoint_;
167template<
class Scalar>
168RCP<ForwardSensitivityStackedStepperStepStrategy<Scalar> >
169forwardSensitivityStackedStepperStepStrategy();
172template<
class Scalar>
173RCP<ForwardSensitivityStackedStepperStepStrategy<Scalar> >
174forwardSensitivityStackedStepperStepStrategy()
176 return Teuchos::rcp(
new ForwardSensitivityStackedStepperStepStrategy<Scalar>());
179template<
class Scalar>
180ForwardSensitivityStackedStepperStepStrategy<Scalar>::ForwardSensitivityStackedStepperStepStrategy()
182 typedef Teuchos::ScalarTraits<Scalar> ST;
183 dt_taken_ = ST::nan();
187template<
class Scalar>
188ForwardSensitivityStackedStepperStepStrategy<Scalar>::~ForwardSensitivityStackedStepperStepStrategy()
193template<
class Scalar>
194void ForwardSensitivityStackedStepperStepStrategy<Scalar>::setupNextStepper(
196 Array<RCP<StepperBase<Scalar> > > &stepperArray
199 typedef Thyra::ModelEvaluatorBase MEB;
201 RCP<StepperBase<Scalar> > baseStepper = stepperArray[0];
202 RCP<Thyra::ModelEvaluator<Scalar> > modelI =
205 Teuchos::rcp_const_cast<Thyra::ModelEvaluator<Scalar> >(
206 stepperArray[i]->getModel()
208 RCP<ForwardSensitivityModelEvaluatorBase<Scalar> > fwdSensME =
209 Teuchos::rcp_dynamic_cast<ForwardSensitivityModelEvaluatorBase<Scalar> > (
212 if (Teuchos::nonnull(fwdSensME)) {
213 bool forceUpToDateW =
true;
214 fwdSensME->initializePointState( Teuchos::inOutArg(*baseStepper), forceUpToDateW);
216 stepperArray[i]->setStepControlData(*baseStepper);
221template<
class Scalar>
222Scalar ForwardSensitivityStackedStepperStepStrategy<Scalar>::evaluateStep(
223 Scalar local_dt_taken,
225 Array<RCP<StepperBase<Scalar> > > &stepperArray
229 dt_taken_ = local_dt_taken;
232 TEUCHOS_TEST_FOR_EXCEPTION( local_dt_taken != dt_taken_, std::logic_error,
233 "Error! sub-stepper["<<i<<
"] did not take the same "
234 "size step as the base stepper!"
241template<
class Scalar>
242void ForwardSensitivityStackedStepperStepStrategy<Scalar>::setStateBasePoint(
243 const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint
246 stateBasePoint_ = stateBasePoint;
250template<
class Scalar>
252 :
virtual public StepperBase<Scalar>,
253 virtual public Teuchos::ParameterListAcceptorDefaultBase
258 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType ScalarMag;
268 void addStepper(
const RCP<StepperBase<Scalar> >& stepper);
272 ArrayView<const RCP<StepperBase<Scalar> > > getNonconstSteppers()
const;
276 void setStackedStepperStepControlStrategy(
277 const RCP<StackedStepperStepStrategyBase<Scalar> >& stepControl
282 RCP<const StackedStepperStepStrategyBase<Scalar> >
283 getStackedStepperStepCongrolStrategy()
const;
291 void setParameterList(RCP<Teuchos::ParameterList>
const& paramList);
293 RCP<const Teuchos::ParameterList> getValidParameters()
const;
301 bool acceptsModel()
const;
305 const RCP<
const Thyra::ModelEvaluator<Scalar> >& model
309 void setNonconstModel(
310 const RCP<Thyra::ModelEvaluator<Scalar> >& model
319 RCP<const Thyra::ModelEvaluator<Scalar> > getModel()
const;
322 RCP<Thyra::ModelEvaluator<Scalar> > getNonconstModel();
336 void setInitialCondition(
337 const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stacked_ic
342 Thyra::ModelEvaluatorBase::InArgs<Scalar> getInitialCondition()
const;
345 Scalar takeStep( Scalar dt, StepSizeType stepType );
348 const StepStatus<Scalar> getStepStatus()
const;
361 RCP<const Thyra::VectorSpaceBase<Scalar> >
366 const Array<Scalar>& time_vec,
367 const Array<Teuchos::RCP<
const Thyra::VectorBase<Scalar> > >& x_vec,
368 const Array<Teuchos::RCP<
const Thyra::VectorBase<Scalar> > >& xdot_vec
372 TimeRange<Scalar> getTimeRange()
const;
376 const Array<Scalar>& time_vec,
377 Array<RCP<
const Thyra::VectorBase<Scalar> > >* x_vec,
378 Array<RCP<
const Thyra::VectorBase<Scalar> > >* xdot_vec,
379 Array<ScalarMag>* accuracy_vec
383 void getNodes(Array<Scalar>* time_vec)
const;
386 void removeNodes(Array<Scalar>& time_vec);
389 int getOrder()
const;
395 mutable bool isInitialized_;
397 Array<RCP<StepperBase<Scalar> > > stepperArray_;
399 mutable Array<RCP<const Thyra::VectorSpaceBase<Scalar> > > xSpaceArray_;
402 mutable Teuchos::RCP<const Thyra::ProductVectorSpaceBase<Scalar> >
407 RCP<StackedStepperStepStrategyBase<Scalar> > stackedStepperStepStrategy_;
410 void setupSpaces_()
const;
419template<
class Scalar>
421RCP<StackedStepper<Scalar> >
424 return Teuchos::rcp(
new StackedStepper<Scalar>());
431template<
class Scalar>
432StackedStepper<Scalar>::StackedStepper()
433 : isInitialized_(false)
436template<
class Scalar>
437void StackedStepper<Scalar>::setupSpaces_()
const
440 if (!isInitialized_) {
441 for (
int i=0 ; i < as<int>(stepperArray_.size()) ; ++i) {
442 xSpaceArray_.push_back(stepperArray_[i]->get_x_space());
446 x_bar_space_ = Thyra::productVectorSpace<Scalar>(xSpaceArray_);
449 isInitialized_ =
true;
453template<
class Scalar>
454void StackedStepper<Scalar>::addStepper(
455 const RCP<StepperBase<Scalar> >& stepper
458 TEUCHOS_ASSERT(!is_null(stepper));
459 stepperArray_.push_back(stepper);
460 isInitialized_ =
false;
465template<
class Scalar>
466ArrayView<const RCP<StepperBase<Scalar> > >
467StackedStepper<Scalar>::getNonconstSteppers()
const
469 return stepperArray_();
473template<
class Scalar>
474void StackedStepper<Scalar>::setStackedStepperStepControlStrategy(
475 const RCP<StackedStepperStepStrategyBase<Scalar> >& stepControl
478 TEUCHOS_ASSERT( Teuchos::nonnull(stepControl) );
479 stackedStepperStepStrategy_ = stepControl;
483template<
class Scalar>
484RCP<const StackedStepperStepStrategyBase<Scalar> >
485StackedStepper<Scalar>::getStackedStepperStepCongrolStrategy()
const
487 return stackedStepperStepStrategy_;
491template<
class Scalar>
492void StackedStepper<Scalar>::setParameterList(
493 RCP<Teuchos::ParameterList>
const& paramList
496 TEUCHOS_TEST_FOR_EXCEPT(is_null(paramList));
497 paramList->validateParameters(*getValidParameters());
498 this->setMyParamList(paramList);
499 Teuchos::readVerboseObjectSublist(&*paramList,
this);
503template<
class Scalar>
504RCP<const Teuchos::ParameterList>
505StackedStepper<Scalar>::getValidParameters()
const
507 static RCP<const ParameterList> validPL;
508 if (is_null(validPL) ) {
509 RCP<ParameterList> pl = Teuchos::parameterList();
510 Teuchos::setupVerboseObjectSublist(&*pl);
519template<
class Scalar>
520bool StackedStepper<Scalar>::acceptsModel()
const
525template<
class Scalar>
526void StackedStepper<Scalar>::setModel(
527 const RCP<
const Thyra::ModelEvaluator<Scalar> >& model
530 TEUCHOS_TEST_FOR_EXCEPT_MSG(
true,
531 "Error, this stepper subclass does not accept a model"
532 " as defined by the StepperBase interface!");
536template<
class Scalar>
537void StackedStepper<Scalar>::setNonconstModel(
538 const RCP<Thyra::ModelEvaluator<Scalar> >& model
541 TEUCHOS_TEST_FOR_EXCEPT_MSG(
true,
542 "Error, this stepper subclass does not accept a model"
543 " as defined by the StepperBase interface!");
547template<
class Scalar>
548RCP<const Thyra::ModelEvaluator<Scalar> >
549StackedStepper<Scalar>::getModel()
const
551 TEUCHOS_TEST_FOR_EXCEPT(
true);
552 return Teuchos::null;
556template<
class Scalar>
557RCP<Thyra::ModelEvaluator<Scalar> >
558StackedStepper<Scalar>::getNonconstModel()
560 TEUCHOS_TEST_FOR_EXCEPT(
true);
561 return Teuchos::null;
565template<
class Scalar>
566void StackedStepper<Scalar>::setInitialCondition(
567 const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stacked_ic
570 TEUCHOS_TEST_FOR_EXCEPT(
true);
574template<
class Scalar>
575Thyra::ModelEvaluatorBase::InArgs<Scalar>
576StackedStepper<Scalar>::getInitialCondition()
const
578 TEUCHOS_TEST_FOR_EXCEPT(
true);
579 Thyra::ModelEvaluatorBase::InArgs<Scalar> inArgs;
584template<
class Scalar>
586StackedStepper<Scalar>::takeStep(
587 Scalar dt, StepSizeType stepType
590 typedef Teuchos::ScalarTraits<Scalar> ST;
592 TEUCHOS_TEST_FOR_EXCEPTION( stepperArray_.size() == 0, std::logic_error,
593 "Error! Rythmos::StackedStepper::takeStep: "
594 "addStepper must be called at least once before takeStep!"
596 this->setupSpaces_();
598 RYTHMOS_FUNC_TIME_MONITOR(
"Rythmos:StackedStepper::takeStep");
600 if (Teuchos::is_null(stackedStepperStepStrategy_)) {
601 stackedStepperStepStrategy_ = defaultStackedStepperStepStrategy<Scalar>();
604 Scalar dt_taken = Scalar(-ST::one());
605 for (
int i=0 ; i<Teuchos::as<int>(stepperArray_.size()) ; ++i) {
606 stackedStepperStepStrategy_->setupNextStepper(i,stepperArray_);
607 Scalar local_dt_taken = stepperArray_[i]->takeStep(dt,stepType);
608 dt_taken = stackedStepperStepStrategy_->evaluateStep(
630template<
class Scalar>
631const StepStatus<Scalar>
632StackedStepper<Scalar>::getStepStatus()
const
634 TEUCHOS_TEST_FOR_EXCEPT(
true);
635 const StepStatus<Scalar> stepStatus;
644template<
class Scalar>
645RCP<const Thyra::VectorSpaceBase<Scalar> >
646StackedStepper<Scalar>::get_x_space()
const
648 this->setupSpaces_();
649 TEUCHOS_TEST_FOR_EXCEPTION( is_null(x_bar_space_), std::logic_error,
650 "Rythmos::StackedStepper::get_x_space(): "
651 "addStepper must be called at least once before get_x_space()!"
653 return(x_bar_space_);
657template<
class Scalar>
658void StackedStepper<Scalar>::addPoints(
659 const Array<Scalar>& time_vec,
660 const Array<Teuchos::RCP<
const Thyra::VectorBase<Scalar> > >& x_vec,
661 const Array<Teuchos::RCP<
const Thyra::VectorBase<Scalar> > >& xdot_vec
664 TEUCHOS_TEST_FOR_EXCEPT(
665 "Not implemented addPoints(...) yet but we could if we wanted!"
670template<
class Scalar>
672StackedStepper<Scalar>::getTimeRange()
const
674 TEUCHOS_TEST_FOR_EXCEPT(
true);
675 TimeRange<Scalar> tr;
680template<
class Scalar>
681void StackedStepper<Scalar>::getPoints(
682 const Array<Scalar>& time_vec,
683 Array<RCP<
const Thyra::VectorBase<Scalar> > >* x_bar_vec,
684 Array<RCP<
const Thyra::VectorBase<Scalar> > >* x_bar_dot_vec,
685 Array<ScalarMag>* accuracy_vec
689 this->setupSpaces_();
690 TEUCHOS_TEST_FOR_EXCEPTION( stepperArray_.size() == 0, std::logic_error,
691 "Rythmos::StackedStepper:getPoints: Error! "
692 "addStepper must be called at least once before getPoints!"
694 int numSteppers = as<int>(stepperArray_.size());
695 int numTimePoints = as<int>(time_vec.size());
698 Array<RCP<Thyra::VectorBase<Scalar> > > x_bar_vec_local;
699 if (x_bar_vec != NULL) {
700 x_bar_vec_local.clear();
702 for (
int n = 0 ; n < numTimePoints ; ++n) {
703 x_bar_vec_local.push_back(Thyra::createMember(this->get_x_space()));
704 x_bar_vec->push_back(x_bar_vec_local[n]);
709 Array<RCP<Thyra::VectorBase<Scalar> > > x_bar_dot_vec_local;
710 if (x_bar_dot_vec != NULL) {
711 x_bar_dot_vec_local.clear();
712 x_bar_dot_vec->clear();
713 for (
int n = 0 ; n < numTimePoints ; ++n) {
714 x_bar_dot_vec_local.push_back(Thyra::createMember(this->get_x_space()));
715 x_bar_dot_vec->push_back(x_bar_dot_vec_local[n]);
721 accuracy_vec->clear();
722 accuracy_vec->resize(numTimePoints);
725 for (
int i=0 ; i < numSteppers; ++i) {
726 Array<RCP<const Thyra::VectorBase<Scalar> > > sub_x_bar_vec;
727 Array<RCP<const Thyra::VectorBase<Scalar> > > sub_x_bar_dot_vec;
728 Array<ScalarMag> sub_accuracy_vec;
729 stepperArray_[i]->getPoints(
731 x_bar_vec ? &sub_x_bar_vec : 0,
732 x_bar_dot_vec ? &sub_x_bar_dot_vec : 0,
733 accuracy_vec ? &sub_accuracy_vec: 0
737 for (
int n=0; n < numTimePoints ; ++n ) {
739 RCP<Thyra::ProductVectorBase<Scalar> > x_bar_pv =
740 Teuchos::rcp_dynamic_cast<Thyra::ProductVectorBase<Scalar> >(
744 RCP<Thyra::VectorBase<Scalar> > xi =
745 x_bar_pv->getNonconstVectorBlock(i);
746 V_V(Teuchos::outArg(*xi),*sub_x_bar_vec[n]);
749 RCP<Thyra::ProductVectorBase<Scalar> > x_bar_dot_pv =
750 Teuchos::rcp_dynamic_cast<Thyra::ProductVectorBase<Scalar> >(
751 x_bar_dot_vec_local[n],
754 RCP<Thyra::VectorBase<Scalar> > xdoti =
755 x_bar_dot_pv->getNonconstVectorBlock(i);
756 V_V(Teuchos::outArg(*xdoti),*sub_x_bar_dot_vec[n]);
759 if ( (i == 0) && accuracy_vec) {
760 (*accuracy_vec)[n] = sub_accuracy_vec[n];
767template<
class Scalar>
768void StackedStepper<Scalar>::getNodes(
769 Array<Scalar>* time_vec
772 TEUCHOS_TEST_FOR_EXCEPT(
true);
776template<
class Scalar>
777void StackedStepper<Scalar>::removeNodes(
778 Array<Scalar>& time_vec
781 TEUCHOS_TEST_FOR_EXCEPT(
"Not implemented yet but we can!");
785template<
class Scalar>
786int StackedStepper<Scalar>::getOrder()
const
788 TEUCHOS_TEST_FOR_EXCEPT(
true);
virtual void addPoints(const Array< Scalar > &time_vec, const Array< RCP< const Thyra::VectorBase< Scalar > > > &x_vec, const Array< RCP< const Thyra::VectorBase< Scalar > > > &xdot_vec)=0
Add points to the buffer.