29#ifndef Rythmos_FORWARDEULER_STEPPER_DEF_H
30#define Rythmos_FORWARDEULER_STEPPER_DEF_H
32#include "Rythmos_ForwardEulerStepper_decl.hpp"
33#include "Rythmos_StepperHelpers.hpp"
34#include "Thyra_ModelEvaluatorHelpers.hpp"
35#include "Thyra_MultiVectorStdOps.hpp"
36#include "Thyra_VectorStdOps.hpp"
37#include "Teuchos_VerboseObjectParameterListHelpers.hpp"
48ForwardEulerStepperMomento<Scalar>::ForwardEulerStepperMomento()
52ForwardEulerStepperMomento<Scalar>::~ForwardEulerStepperMomento()
56void ForwardEulerStepperMomento<Scalar>::serialize(
57 const StateSerializerStrategy<Scalar>& stateSerializer,
61 using Teuchos::is_null;
62 TEUCHOS_ASSERT( !is_null(model_) );
63 RCP<VectorBase<Scalar> > sol_vec = solution_vector_;
64 if (is_null(sol_vec)) {
65 sol_vec = Thyra::createMember(model_->get_x_space());
67 RCP<VectorBase<Scalar> > res_vec = residual_vector_;
68 if (is_null(res_vec)) {
69 res_vec = Thyra::createMember(model_->get_f_space());
71 RCP<VectorBase<Scalar> > sol_vec_old = solution_vector_old_;
72 if (is_null(sol_vec_old)) {
73 sol_vec_old = Thyra::createMember(model_->get_x_space());
75 stateSerializer.serializeVectorBase(*sol_vec,oStream);
76 stateSerializer.serializeVectorBase(*res_vec,oStream);
77 stateSerializer.serializeVectorBase(*sol_vec_old,oStream);
78 stateSerializer.serializeScalar(t_,oStream);
79 stateSerializer.serializeScalar(t_old_,oStream);
80 stateSerializer.serializeScalar(dt_,oStream);
81 stateSerializer.serializeInt(numSteps_,oStream);
82 stateSerializer.serializeBool(isInitialized_,oStream);
83 stateSerializer.serializeBool(haveInitialCondition_,oStream);
84 RCP<ParameterList> pl = parameterList_;
85 if (Teuchos::is_null(pl)) {
86 pl = Teuchos::parameterList();
88 stateSerializer.serializeParameterList(*pl,oStream);
92void ForwardEulerStepperMomento<Scalar>::deSerialize(
93 const StateSerializerStrategy<Scalar>& stateSerializer,
97 using Teuchos::outArg;
98 using Teuchos::is_null;
99 TEUCHOS_ASSERT( !is_null(model_) );
100 if (is_null(solution_vector_)) {
101 solution_vector_ = Thyra::createMember(*model_->get_x_space());
103 if (is_null(residual_vector_)) {
104 residual_vector_ = Thyra::createMember(*model_->get_f_space());
106 if (is_null(solution_vector_old_)) {
107 solution_vector_old_ = Thyra::createMember(*model_->get_x_space());
109 stateSerializer.deSerializeVectorBase(outArg(*solution_vector_),iStream);
110 stateSerializer.deSerializeVectorBase(outArg(*residual_vector_),iStream);
111 stateSerializer.deSerializeVectorBase(outArg(*solution_vector_old_),iStream);
112 stateSerializer.deSerializeScalar(outArg(t_),iStream);
113 stateSerializer.deSerializeScalar(outArg(t_old_),iStream);
114 stateSerializer.deSerializeScalar(outArg(dt_),iStream);
115 stateSerializer.deSerializeInt(outArg(numSteps_),iStream);
116 stateSerializer.deSerializeBool(outArg(isInitialized_),iStream);
117 stateSerializer.deSerializeBool(outArg(haveInitialCondition_),iStream);
118 if (is_null(parameterList_)) {
119 parameterList_ = Teuchos::parameterList();
121 stateSerializer.deSerializeParameterList(outArg(*parameterList_),iStream);
124template<
class Scalar>
125RCP<MomentoBase<Scalar> > ForwardEulerStepperMomento<Scalar>::clone()
const
127 RCP<ForwardEulerStepperMomento<Scalar> > m = rcp(
new ForwardEulerStepperMomento<Scalar>());
128 m->set_solution_vector(solution_vector_);
129 m->set_residual_vector(residual_vector_);
130 m->set_solution_vector_old(solution_vector_old_);
132 m->set_t_old(t_old_);
134 m->set_numSteps(numSteps_);
135 m->set_isInitialized(isInitialized_);
136 m->set_haveInitialCondition(haveInitialCondition_);
137 m->set_parameterList(parameterList_);
138 if (!Teuchos::is_null(this->getMyParamList())) {
139 m->setParameterList(Teuchos::parameterList(*(this->getMyParamList())));
141 m->set_model(model_);
142 m->set_basePoint(basePoint_);
148template<
class Scalar>
149void ForwardEulerStepperMomento<Scalar>::set_solution_vector(
const RCP<
const VectorBase<Scalar> >& solution_vector )
151 solution_vector_ = Teuchos::null;
152 if (!Teuchos::is_null(solution_vector)) {
153 solution_vector_ = solution_vector->clone_v();
157template<
class Scalar>
158RCP<VectorBase<Scalar> > ForwardEulerStepperMomento<Scalar>::get_solution_vector()
const
160 return solution_vector_;
163template<
class Scalar>
164void ForwardEulerStepperMomento<Scalar>::set_residual_vector(
const RCP<
const VectorBase<Scalar> >& residual_vector)
166 residual_vector_ = Teuchos::null;
167 if (!Teuchos::is_null(residual_vector)) {
168 residual_vector_ = residual_vector->clone_v();
172template<
class Scalar>
173RCP<VectorBase<Scalar> > ForwardEulerStepperMomento<Scalar>::get_residual_vector()
const
175 return residual_vector_;
178template<
class Scalar>
179void ForwardEulerStepperMomento<Scalar>::set_solution_vector_old(
const RCP<
const VectorBase<Scalar> >& solution_vector_old )
181 solution_vector_old_ = Teuchos::null;
182 if (!Teuchos::is_null(solution_vector_old)) {
183 solution_vector_old_ = solution_vector_old->clone_v();
187template<
class Scalar>
188RCP<VectorBase<Scalar> > ForwardEulerStepperMomento<Scalar>::get_solution_vector_old()
const
190 return solution_vector_old_;
193template<
class Scalar>
194void ForwardEulerStepperMomento<Scalar>::set_t(
const Scalar & t)
199template<
class Scalar>
200Scalar ForwardEulerStepperMomento<Scalar>::get_t()
const
205template<
class Scalar>
206void ForwardEulerStepperMomento<Scalar>::set_t_old(
const Scalar & t_old)
211template<
class Scalar>
212Scalar ForwardEulerStepperMomento<Scalar>::get_t_old()
const
217template<
class Scalar>
218void ForwardEulerStepperMomento<Scalar>::set_dt(
const Scalar & dt)
223template<
class Scalar>
224Scalar ForwardEulerStepperMomento<Scalar>::get_dt()
const
229template<
class Scalar>
230void ForwardEulerStepperMomento<Scalar>::set_numSteps(
const int & numSteps)
232 numSteps_ = numSteps;
235template<
class Scalar>
236int ForwardEulerStepperMomento<Scalar>::get_numSteps()
const
241template<
class Scalar>
242void ForwardEulerStepperMomento<Scalar>::set_isInitialized(
const bool & isInitialized)
244 isInitialized_ = isInitialized;
247template<
class Scalar>
248bool ForwardEulerStepperMomento<Scalar>::get_isInitialized()
const
250 return isInitialized_;
253template<
class Scalar>
254void ForwardEulerStepperMomento<Scalar>::set_haveInitialCondition(
const bool & haveInitialCondition)
256 haveInitialCondition_ = haveInitialCondition;
259template<
class Scalar>
260bool ForwardEulerStepperMomento<Scalar>::get_haveInitialCondition()
const
262 return haveInitialCondition_;
265template<
class Scalar>
266void ForwardEulerStepperMomento<Scalar>::set_parameterList(
const RCP<const ParameterList>& pl)
268 parameterList_ = Teuchos::null;
269 if (!Teuchos::is_null(pl)) {
270 parameterList_ = Teuchos::parameterList(*pl);
274template<
class Scalar>
275RCP<ParameterList> ForwardEulerStepperMomento<Scalar>::get_parameterList()
const
277 return parameterList_;
280template<
class Scalar>
281void ForwardEulerStepperMomento<Scalar>::setParameterList(
const RCP<ParameterList>& paramList)
283 this->setMyParamList(paramList);
286template<
class Scalar>
287RCP<const ParameterList> ForwardEulerStepperMomento<Scalar>::getValidParameters()
const
289 return Teuchos::null;
292template<
class Scalar>
293void ForwardEulerStepperMomento<Scalar>::set_model(
const RCP<
const Thyra::ModelEvaluator<Scalar> >& model)
298template<
class Scalar>
299RCP<const Thyra::ModelEvaluator<Scalar> > ForwardEulerStepperMomento<Scalar>::get_model()
const
304template<
class Scalar>
305void ForwardEulerStepperMomento<Scalar>::set_basePoint(
const RCP<
const Thyra::ModelEvaluatorBase::InArgs<Scalar> >& basePoint)
307 basePoint_ = basePoint;
310template<
class Scalar>
311RCP<const Thyra::ModelEvaluatorBase::InArgs<Scalar> > ForwardEulerStepperMomento<Scalar>::get_basePoint()
const
321template<
class Scalar>
322RCP<ForwardEulerStepper<Scalar> > forwardEulerStepper()
324 RCP<ForwardEulerStepper<Scalar> > stepper = rcp(
new ForwardEulerStepper<Scalar>());
329template<
class Scalar>
330RCP<ForwardEulerStepper<Scalar> > forwardEulerStepper(
331 const RCP<Thyra::ModelEvaluator<Scalar> >& model
334 RCP<ForwardEulerStepper<Scalar> > stepper = forwardEulerStepper<Scalar>();
336 RCP<StepperBase<Scalar> > stepperBase =
337 Teuchos::rcp_dynamic_cast<StepperBase<Scalar> >(stepper,
true);
338 setStepperModel(Teuchos::inOutArg(*stepperBase),model);
343template<
class Scalar>
346 this->defaultInitializAll_();
347 dt_ = Teuchos::ScalarTraits<Scalar>::zero();
351template<
class Scalar>
354 model_ = Teuchos::null;
355 solution_vector_ = Teuchos::null;
356 residual_vector_ = Teuchos::null;
360 solution_vector_old_ = Teuchos::null;
363 haveInitialCondition_ =
false;
364 parameterList_ = Teuchos::null;
365 isInitialized_ =
false;
368template<
class Scalar>
369void ForwardEulerStepper<Scalar>::initialize_()
371 if (!isInitialized_) {
372 TEUCHOS_TEST_FOR_EXCEPTION( is_null(model_), std::logic_error,
373 "Error! Please set a model on the stepper.\n"
375 residual_vector_ = Thyra::createMember(model_->get_f_space());
376 isInitialized_ =
true;
380template<
class Scalar>
385template<
class Scalar>
388 TEUCHOS_TEST_FOR_EXCEPTION(!haveInitialCondition_,std::logic_error,
"Error, attempting to call get_x_space before setting an initial condition!\n");
389 return(solution_vector_->space());
392template<
class Scalar>
395 TEUCHOS_TEST_FOR_EXCEPTION( !haveInitialCondition_, std::logic_error,
396 "Error! Attempting to call takeStep before setting an initial condition!\n"
399 if (flag == STEP_TYPE_VARIABLE) {
404 eval_model_explicit<Scalar>(*model_,basePoint_,*solution_vector_,t_+dt,Teuchos::outArg(*residual_vector_));
407 Thyra::V_V(Teuchos::outArg(*solution_vector_old_),*solution_vector_);
409 Thyra::Vp_StV(solution_vector_.ptr(),dt,*residual_vector_);
418template<
class Scalar>
423 if (!haveInitialCondition_) {
424 stepStatus.
stepStatus = STEP_STATUS_UNINITIALIZED;
426 else if (numSteps_ == 0) {
428 stepStatus.
order = this->getOrder();
429 stepStatus.
time = t_;
430 stepStatus.
solution = solution_vector_;
433 stepStatus.
stepStatus = STEP_STATUS_CONVERGED;
435 stepStatus.
order = this->getOrder();
436 stepStatus.
time = t_;
438 stepStatus.
solution = solution_vector_;
439 stepStatus.
residual = residual_vector_;
445template<
class Scalar>
448 std::string name =
"Rythmos::ForwardEulerStepper";
452template<
class Scalar>
454 Teuchos::FancyOStream &out,
455 const Teuchos::EVerbosityLevel verbLevel
458 if ( (
static_cast<int>(verbLevel) ==
static_cast<int>(Teuchos::VERB_DEFAULT) ) ||
459 (
static_cast<int>(verbLevel) >=
static_cast<int>(Teuchos::VERB_LOW) )
461 out << description() <<
"::describe" << std::endl;
462 out <<
"model = " << model_->description() << std::endl;
463 }
else if (
static_cast<int>(verbLevel) >=
static_cast<int>(Teuchos::VERB_LOW)) {
464 }
else if (
static_cast<int>(verbLevel) >=
static_cast<int>(Teuchos::VERB_MEDIUM)) {
465 }
else if (
static_cast<int>(verbLevel) >=
static_cast<int>(Teuchos::VERB_HIGH)) {
466 out <<
"model = " << std::endl;
467 model_->describe(out,verbLevel);
468 out <<
"solution_vector = " << std::endl;
469 solution_vector_->describe(out,verbLevel);
470 out <<
"residual_vector = " << std::endl;
471 residual_vector_->describe(out,verbLevel);
475template<
class Scalar>
478 ,
const Array<Teuchos::RCP<
const Thyra::VectorBase<Scalar> > >&
479 ,
const Array<Teuchos::RCP<
const Thyra::VectorBase<Scalar> > >&
482 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::logic_error,
"Error, addPoints is not implemented for ForwardEulerStepper.\n");
485template<
class Scalar>
487 const Array<Scalar>& time_vec
488 ,Array<Teuchos::RCP<
const Thyra::VectorBase<Scalar> > >* x_vec
489 ,Array<Teuchos::RCP<
const Thyra::VectorBase<Scalar> > >* xdot_vec
490 ,Array<ScalarMag>* accuracy_vec)
const
492 TEUCHOS_ASSERT( haveInitialCondition_ );
493 using Teuchos::constOptInArg;
495 defaultGetPoints<Scalar>(
497 constOptInArg(*solution_vector_old_),
498 Ptr<
const VectorBase<Scalar> >(null),
500 constOptInArg(*solution_vector_),
501 Ptr<
const VectorBase<Scalar> >(null),
510template<
class Scalar>
513 if (!haveInitialCondition_) {
514 return(invalidTimeRange<Scalar>());
520template<
class Scalar>
523 TEUCHOS_ASSERT( time_vec != NULL );
525 if (!haveInitialCondition_) {
528 time_vec->push_back(t_old_);
531 time_vec->push_back(t_);
535template<
class Scalar>
538 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::logic_error,
"Error, removeNodes is not implemented for ForwardEulerStepper.\n");
541template<
class Scalar>
547template <
class Scalar>
550 TEUCHOS_TEST_FOR_EXCEPT(is_null(paramList));
551 paramList->validateParametersAndSetDefaults(*this->getValidParameters());
552 parameterList_ = paramList;
553 Teuchos::readVerboseObjectSublist(&*parameterList_,
this);
556template <
class Scalar>
559 return(parameterList_);
562template <
class Scalar>
565 Teuchos::RCP<Teuchos::ParameterList> temp_param_list = parameterList_;
566 parameterList_ = Teuchos::null;
567 return(temp_param_list);
571template<
class Scalar>
572RCP<const Teuchos::ParameterList>
575 using Teuchos::ParameterList;
576 static RCP<const ParameterList> validPL;
577 if (is_null(validPL)) {
578 RCP<ParameterList> pl = Teuchos::parameterList();
579 Teuchos::setupVerboseObjectSublist(&*pl);
586template<
class Scalar>
589 TEUCHOS_TEST_FOR_EXCEPT( is_null(model) );
590 assertValidModel( *
this, *model );
595template<
class Scalar>
598 this->setModel(model);
602template<
class Scalar>
603Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >
609template<
class Scalar>
610Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >
613 return Teuchos::null;
616template<
class Scalar>
618 const Thyra::ModelEvaluatorBase::InArgs<Scalar> &initialCondition
621 basePoint_ = initialCondition;
622 t_ = initialCondition.get_t();
624 solution_vector_ = initialCondition.get_x()->clone_v();
625 solution_vector_old_ = solution_vector_->clone_v();
626 haveInitialCondition_ =
true;
630template<
class Scalar>
631Thyra::ModelEvaluatorBase::InArgs<Scalar>
638template<
class Scalar>
644template<
class Scalar>
645RCP<StepperBase<Scalar> >
652 RCP<StepperBase<Scalar> >
655 if (!is_null(model_)) {
656 setStepperModel(Teuchos::inOutArg(*stepper),model_);
659 if (!is_null(parameterList_)) {
660 stepper->setParameterList(Teuchos::parameterList(*parameterList_));
667template<
class Scalar>
668RCP<const MomentoBase<Scalar> >
672 momento->set_solution_vector(solution_vector_);
673 momento->set_solution_vector_old(solution_vector_old_);
674 momento->set_residual_vector(residual_vector_);
676 momento->set_t_old(t_old_);
677 momento->set_dt(dt_);
678 momento->set_numSteps(numSteps_);
679 momento->set_isInitialized(isInitialized_);
680 momento->set_haveInitialCondition(haveInitialCondition_);
681 momento->set_parameterList(parameterList_);
682 momento->set_model(model_);
683 RCP<Thyra::ModelEvaluatorBase::InArgs<Scalar> > bp = Teuchos::rcp(
new Thyra::ModelEvaluatorBase::InArgs<Scalar>(basePoint_));
684 momento->set_basePoint(bp);
688template<
class Scalar>
691 Ptr<const ForwardEulerStepperMomento<Scalar> > feMomentoPtr =
692 Teuchos::ptr_dynamic_cast<const ForwardEulerStepperMomento<Scalar> >(momentoPtr,
true);
693 TEUCHOS_TEST_FOR_EXCEPTION( Teuchos::is_null(feMomentoPtr->get_model()), std::logic_error,
694 "Error! Rythmos::ForwardEulerStepper::setMomento: The momento must have a valid model through set_model(...) prior to calling ForwardEulerStepper::setMomento(...)."
696 TEUCHOS_TEST_FOR_EXCEPTION( Teuchos::is_null(feMomentoPtr->get_basePoint()), std::logic_error,
697 "Error! Rythmos::ForwardEulerStepper::setMomento: The momento must have a valid base point through set_basePoint(...) prior to calling ForwardEulerStepper::setMomento(...)."
699 model_ = feMomentoPtr->get_model();
700 basePoint_ = *(feMomentoPtr->get_basePoint());
702 solution_vector_ = feMomento.get_solution_vector();
703 solution_vector_old_ = feMomento.get_solution_vector_old();
704 residual_vector_ = feMomento.get_residual_vector();
705 t_ = feMomento.get_t();
706 t_old_ = feMomento.get_t_old();
707 dt_ = feMomento.get_dt();
708 numSteps_ = feMomento.get_numSteps();
709 isInitialized_ = feMomento.get_isInitialized();
710 haveInitialCondition_ = feMomento.get_haveInitialCondition();
711 parameterList_ = feMomento.get_parameterList();
712 this->checkConsistentState_();
715template<
class Scalar>
718 if (isInitialized_) {
719 TEUCHOS_ASSERT( !Teuchos::is_null(model_) );
720 TEUCHOS_ASSERT( !Teuchos::is_null(residual_vector_) );
722 if (haveInitialCondition_) {
723 TEUCHOS_ASSERT( !ST::isnaninf(t_) );
724 TEUCHOS_ASSERT( !ST::isnaninf(t_old_) );
725 TEUCHOS_ASSERT( !Teuchos::is_null(solution_vector_) );
726 TEUCHOS_ASSERT( !Teuchos::is_null(solution_vector_old_) );
727 TEUCHOS_ASSERT( t_ >= basePoint_.get_t() );
728 TEUCHOS_ASSERT( t_old_ >= basePoint_.get_t() );
731 TEUCHOS_ASSERT(isInitialized_);
732 TEUCHOS_ASSERT(haveInitialCondition_);
742#define RYTHMOS_FORWARD_EULER_STEPPER_INSTANT(SCALAR) \
744 template class ForwardEulerStepperMomento< SCALAR >; \
745 template class ForwardEulerStepper< SCALAR >; \
747 template RCP<ForwardEulerStepper< SCALAR > > forwardEulerStepper(); \
748 template RCP<ForwardEulerStepper< SCALAR > > forwardEulerStepper( \
749 const RCP<Thyra::ModelEvaluator< SCALAR > >& model \
Concrete momento class for the ForwardEulerStepper.
std::string description() const
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)
RCP< const Thyra::ModelEvaluator< Scalar > > getModel() const
RCP< const Thyra::VectorSpaceBase< Scalar > > get_x_space() const
void setNonconstModel(const RCP< Thyra::ModelEvaluator< Scalar > > &model)
void setInitialCondition(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &initialCondition)
RCP< Teuchos::ParameterList > unsetParameterList()
void getNodes(Array< Scalar > *time_vec) const
Get interpolation nodes.
void setMomento(const Ptr< const MomentoBase< Scalar > > &momentoPtr)
Set momento object for use in restarts.
RCP< Teuchos::ParameterList > getNonconstParameterList()
void getPoints(const Array< Scalar > &time_vec, Array< RCP< const Thyra::VectorBase< Scalar > > > *x_vec, Array< RCP< const Thyra::VectorBase< Scalar > > > *xdot_vec, Array< ScalarMag > *accuracy_vec) const
Get values from buffer.
bool supportsCloning() const
RCP< const Teuchos::ParameterList > getValidParameters() const
int getOrder() const
Get order of interpolation.
Thyra::ModelEvaluatorBase::InArgs< Scalar > getInitialCondition() const
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
RCP< Thyra::ModelEvaluator< Scalar > > getNonconstModel()
void removeNodes(Array< Scalar > &time_vec)
Remove interpolation nodes.
void setParameterList(RCP< Teuchos::ParameterList > const ¶mList)
Redefined from Teuchos::ParameterListAcceptor.
TimeRange< Scalar > getTimeRange() const
const StepStatus< Scalar > getStepStatus() const
RCP< StepperBase< Scalar > > cloneStepperAlgorithm() const
void setModel(const RCP< const Thyra::ModelEvaluator< Scalar > > &model)
Scalar takeStep(Scalar dt, StepSizeType flag)
RCP< const MomentoBase< Scalar > > getMomento() const
Get momento object for use in restarts.
Base strategy class for interpolation functionality.
Base class for a momento object.
RCP< const Thyra::VectorBase< Scalar > > solution
RCP< const Thyra::VectorBase< Scalar > > residual