29 #ifndef Rythmos_THETA_STEPPER_DEF_H 30 #define Rythmos_THETA_STEPPER_DEF_H 32 #include "Rythmos_ConfigDefs.h" 33 #ifdef HAVE_RYTHMOS_EXPERIMENTAL 35 #include "Rythmos_ThetaStepper_decl.hpp" 44 template<
class Scalar>
45 RCP<ThetaStepper<Scalar> >
47 const RCP<Thyra::ModelEvaluator<Scalar> >& model,
48 const RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
49 RCP<Teuchos::ParameterList>& parameterList
52 Teuchos::RCP<ThetaStepper<Scalar> > stepper =
53 Teuchos::rcp(
new ThetaStepper<Scalar>());
54 stepper->setParameterList(parameterList);
55 stepper->setModel(model);
56 stepper->setSolver(solver);
68 template<
class Scalar>
69 ThetaStepper<Scalar>::ThetaStepper()
71 typedef Teuchos::ScalarTraits<Scalar> ST;
72 this->defaultInitializeAll_();
73 Scalar zero = ST::zero();
81 template<
class Scalar>
82 void ThetaStepper<Scalar>::defaultInitializeAll_()
84 typedef Teuchos::ScalarTraits<Scalar> ST;
85 isInitialized_ =
false;
86 haveInitialCondition_ =
false;
87 model_ = Teuchos::null;
88 solver_ = Teuchos::null;
91 x_old_ = Teuchos::null;
92 x_pre_ = Teuchos::null;
93 x_dot_ = Teuchos::null;
94 x_dot_old_ = Teuchos::null;
95 x_dot_really_old_ = Teuchos::null;
96 x_dot_base_ = Teuchos::null;
102 thetaStepperType_ = INVALID_THETA_STEPPER_TYPE;
104 neModel_ = Teuchos::null;
105 parameterList_ = Teuchos::null;
106 interpolator_ = Teuchos::null;
107 predictor_corrector_begin_after_step_ = -1;
108 default_predictor_order_ = -1;
111 template<
class Scalar>
112 bool ThetaStepper<Scalar>::isImplicit()
const 118 template<
class Scalar>
119 void ThetaStepper<Scalar>::setInterpolator(
120 const RCP<InterpolatorBase<Scalar> >& interpolator
123 #ifdef HAVE_RYTHMOS_DEBUG 124 TEUCHOS_TEST_FOR_EXCEPT(is_null(interpolator));
126 interpolator_ = interpolator;
127 isInitialized_ =
false;
130 template<
class Scalar>
131 RCP<InterpolatorBase<Scalar> >
132 ThetaStepper<Scalar>::getNonconstInterpolator()
134 return interpolator_;
137 template<
class Scalar>
138 RCP<const InterpolatorBase<Scalar> >
139 ThetaStepper<Scalar>::getInterpolator()
const 141 return interpolator_;
145 template<
class Scalar>
146 RCP<InterpolatorBase<Scalar> >
147 ThetaStepper<Scalar>::unSetInterpolator()
149 RCP<InterpolatorBase<Scalar> > temp_interpolator = interpolator_;
150 interpolator_ = Teuchos::null;
151 return(temp_interpolator);
152 isInitialized_ =
false;
159 template<
class Scalar>
160 void ThetaStepper<Scalar>::setSolver(
161 const RCP<Thyra::NonlinearSolverBase<Scalar> > &solver
166 TEUCHOS_TEST_FOR_EXCEPTION(solver == Teuchos::null, std::logic_error,
167 "Error! Thyra::NonlinearSolverBase RCP passed in through ThetaStepper::setSolver is null!" 170 RCP<Teuchos::FancyOStream> out = this->getOStream();
171 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
172 Teuchos::OSTab ostab(out,1,
"TS::setSolver");
173 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
174 *out <<
"solver = " << solver->description() << std::endl;
179 isInitialized_ =
false;
184 template<
class Scalar>
185 RCP<Thyra::NonlinearSolverBase<Scalar> >
186 ThetaStepper<Scalar>::getNonconstSolver()
192 template<
class Scalar>
193 RCP<const Thyra::NonlinearSolverBase<Scalar> >
194 ThetaStepper<Scalar>::getSolver()
const 203 template<
class Scalar>
204 bool ThetaStepper<Scalar>::supportsCloning()
const 209 template<
class Scalar>
210 RCP<StepperBase<Scalar> >
211 ThetaStepper<Scalar>::cloneStepperAlgorithm()
const 213 RCP<ThetaStepper<Scalar> >
214 stepper = Teuchos::rcp(
new ThetaStepper<Scalar>);
215 stepper->isInitialized_ = isInitialized_;
216 stepper->model_ = model_;
218 if (!is_null(solver_))
219 stepper->solver_ = solver_->cloneNonlinearSolver().assert_not_null();
221 stepper->basePoint_ = basePoint_;
224 stepper->x_ = x_->clone_v().assert_not_null();
225 if (!is_null(x_old_))
226 stepper->x_old_ = x_old_->clone_v().assert_not_null();
227 if (!is_null(x_pre_))
228 stepper->x_pre_ = x_pre_->clone_v().assert_not_null();
230 if (!is_null(x_dot_))
231 stepper->x_dot_ = x_dot_->clone_v().assert_not_null();
232 if (!is_null(x_dot_old_))
233 stepper->x_dot_old_ = x_dot_old_->clone_v().assert_not_null();
234 if (!is_null(x_dot_really_old_))
235 stepper->x_dot_really_old_ = x_dot_really_old_->clone_v().assert_not_null();
236 if (!is_null(x_dot_base_))
237 stepper->x_dot_base_ = x_dot_base_->clone_v().assert_not_null();
240 stepper->t_old_ = t_old_;
243 stepper->dt_old_ = dt_old_;
245 stepper->numSteps_ = numSteps_;
247 stepper->thetaStepperType_ = thetaStepperType_;
248 stepper->theta_ = theta_;
249 stepper->predictor_corrector_begin_after_step_ = predictor_corrector_begin_after_step_;
250 stepper->default_predictor_order_ = default_predictor_order_;
252 if (!is_null(neModel_))
256 if (!is_null(parameterList_))
257 stepper->parameterList_ = parameterList(*parameterList_);
259 if (!is_null(interpolator_))
260 stepper->interpolator_
261 = interpolator_->cloneInterpolator().assert_not_null();
266 template<
class Scalar>
267 void ThetaStepper<Scalar>::setModel(
268 const RCP<
const Thyra::ModelEvaluator<Scalar> >& model
274 TEUCHOS_TEST_FOR_EXCEPT( is_null(model) );
275 assertValidModel( *
this, *model );
277 RCP<Teuchos::FancyOStream> out = this->getOStream();
278 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
279 Teuchos::OSTab ostab(out,1,
"TS::setModel");
280 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
281 *out <<
"model = " << model->description() << std::endl;
288 x_old_ = Teuchos::null;
289 x_pre_ = Teuchos::null;
291 x_dot_ = Teuchos::null;
292 x_dot_old_ = Teuchos::null;
293 x_dot_really_old_ = Teuchos::null;
294 x_dot_base_ = Teuchos::null;
296 isInitialized_ =
false;
297 haveInitialCondition_ = setDefaultInitialConditionFromNominalValues<Scalar>(
298 *model_, Teuchos::ptr(
this) );
302 template<
class Scalar>
303 void ThetaStepper<Scalar>::setNonconstModel(
304 const RCP<Thyra::ModelEvaluator<Scalar> >& model
307 this->setModel(model);
311 template<
class Scalar>
312 RCP<const Thyra::ModelEvaluator<Scalar> >
313 ThetaStepper<Scalar>::getModel()
const 319 template<
class Scalar>
320 RCP<Thyra::ModelEvaluator<Scalar> >
321 ThetaStepper<Scalar>::getNonconstModel()
323 TEUCHOS_TEST_FOR_EXCEPT(
true);
324 return Teuchos::null;
328 template<
class Scalar>
329 void ThetaStepper<Scalar>::setInitialCondition(
330 const Thyra::ModelEvaluatorBase::InArgs<Scalar> &initialCondition
334 typedef Teuchos::ScalarTraits<Scalar> ST;
335 typedef Thyra::ModelEvaluatorBase MEB;
337 basePoint_ = initialCondition;
341 RCP<const Thyra::VectorBase<Scalar> >
342 x_init = initialCondition.get_x();
344 #ifdef HAVE_RYTHMOS_DEBUG 345 TEUCHOS_TEST_FOR_EXCEPTION(
346 is_null(x_init), std::logic_error,
347 "Error, if the client passes in an intial condition to setInitialCondition(...),\n" 348 "then x can not be null!" );
351 x_ = x_init->clone_v();
355 RCP<const Thyra::VectorBase<Scalar> >
356 x_dot_init = initialCondition.get_x_dot();
358 if (!is_null(x_dot_init)) {
359 x_dot_ = x_dot_init->clone_v();
362 x_dot_ = createMember(x_->space());
363 assign(x_dot_.ptr(),ST::zero());
368 t_ = initialCondition.get_t();
376 x_pre_ = x_->clone_v();
380 x_old_ = x_->clone_v();
384 x_dot_base_ = x_->clone_v();
388 x_dot_old_ = x_dot_->clone_v();
392 x_dot_really_old_ = x_dot_->clone_v();
394 haveInitialCondition_ =
true;
399 template<
class Scalar>
400 Thyra::ModelEvaluatorBase::InArgs<Scalar>
401 ThetaStepper<Scalar>::getInitialCondition()
const 407 template<
class Scalar>
408 Scalar ThetaStepper<Scalar>::takeStep(Scalar dt, StepSizeType stepSizeType)
412 using Teuchos::incrVerbLevel;
413 typedef Teuchos::ScalarTraits<Scalar> ST;
414 typedef Thyra::NonlinearSolverBase<Scalar> NSB;
415 typedef Teuchos::VerboseObjectTempState<NSB> VOTSNSB;
423 RCP<Teuchos::FancyOStream> out = this->getOStream();
424 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
425 Teuchos::OSTab ostab(out,1,
"TS::takeStep");
426 VOTSNSB solver_outputTempState(solver_,out,incrVerbLevel(verbLevel,-1));
428 if ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
430 <<
"\nEntering " << Teuchos::TypeNameTraits<ThetaStepper<Scalar> >::name()
431 <<
"::takeStep("<<dt<<
","<<toString(stepSizeType)<<
") ...\n";
435 V_StV( x_old_.ptr(), Scalar(ST::one()), *x_ );
436 V_StV( x_dot_really_old_.ptr(), Scalar(ST::one()), *x_dot_old_ );
437 V_StV( x_dot_old_.ptr(), Scalar(ST::one()), *x_dot_ );
439 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_EXTREME) ) {
440 *out <<
"\nSetting dt_ and old data ..." << std::endl;
441 *out <<
"\ndt_ = " << dt_;
442 *out <<
"\nx_old_ = " << *x_old_;
443 *out <<
"\nx_dot_old_ = " << *x_dot_old_;
444 *out <<
"\nx_dot_really_old_ = " << *x_dot_really_old_;
447 if ((stepSizeType == STEP_TYPE_VARIABLE) || (dt == ST::zero())) {
448 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) )
449 *out <<
"\nThe arguments to takeStep are not valid for ThetaStepper at this time." << std::endl;
450 return(Scalar(-ST::one()));
452 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
453 *out <<
"\ndt = " << dt << std::endl;
454 *out <<
"\nstepSizeType = " << stepSizeType << std::endl;
470 const double theta = (numSteps_+1>=predictor_corrector_begin_after_step_) ? theta_ : 1.0;
472 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_EXTREME) ) {
473 *out <<
"\nSetting x_dot_base_ ..." << std::endl;
474 *out <<
"\ntheta = " << theta;
475 *out <<
"\nx_ = " << *x_;
476 *out <<
"\nx_dot_old_ = " << *x_dot_old_;
479 const Scalar coeff_x_dot = Scalar(ST::one()/(theta*dt));
480 const Scalar coeff_x = ST::one();
482 const Scalar x_coeff = Scalar(-coeff_x_dot);
483 const Scalar x_dot_old_coeff = Scalar( -(ST::one()-theta)/theta);
485 V_StV( x_dot_base_.ptr(), x_coeff, *x_old_ );
486 Vp_StV( x_dot_base_.ptr(), x_dot_old_coeff, *x_dot_old_);
488 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_EXTREME) ) {
489 *out <<
"\nx_dot_base_ = " << *x_dot_base_;
492 if(!neModel_.get()) {
496 neModel_->initializeSingleResidualModel(
506 if( solver_->getModel().get() != neModel_.get() ) {
507 solver_->setModel(neModel_);
510 solver_->setVerbLevel(this->getVerbLevel());
516 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
517 *out <<
"\nSolving the implicit theta-stepper timestep equation" 518 <<
" with theta = " << theta <<
"\n";
521 Thyra::SolveStatus<Scalar>
522 neSolveStatus = solver_->solve(&*x_);
528 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
529 *out <<
"\nOutput status of nonlinear solve:\n" << neSolveStatus;
543 V_StV( x_dot_.ptr(), Scalar( ST::one()/(theta*dt)), *x_ );
544 Vp_StV( x_dot_.ptr(), Scalar(-ST::one()/(theta*dt)), *x_old_ );
545 Vp_StV( x_dot_.ptr(), Scalar( -(ST::one()-theta)/theta), *x_dot_old_ );
547 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_EXTREME) ) {
548 *out <<
"\nUpdating x_dot_ ...\n";
549 *out <<
"\nx_dot_ = " << *x_dot_ << std::endl;
558 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
559 *out <<
"\nt_old_ = " << t_old_ << std::endl;
560 *out <<
"\nt_ = " << t_ << std::endl;
563 #ifdef HAVE_RYTHMOS_DEBUG 565 if ( includesVerbLevel(verbLevel,Teuchos::VERB_LOW) )
566 *out <<
"\nChecking to make sure that solution and the interpolated solution are the same! ...\n";
570 typedef ScalarTraits<Scalar> ST;
571 typedef typename ST::magnitudeType ScalarMag;
572 typedef ScalarTraits<ScalarMag> SMT;
574 Teuchos::OSTab tab(out);
576 const StepStatus<Scalar> stepStatus = this->getStepStatus();
578 RCP<const Thyra::VectorBase<Scalar> >
579 x = stepStatus.solution,
580 xdot = stepStatus.solutionDot;
582 Array<Scalar> time_vec = Teuchos::tuple(stepStatus.time);
583 Array<RCP<const Thyra::VectorBase<Scalar> > > x_vec, xdot_vec;
584 this->getPoints(time_vec,&x_vec,&xdot_vec,0);
586 RCP<const Thyra::VectorBase<Scalar> >
588 xdot_interp = xdot_vec[0];
590 TEUCHOS_TEST_FOR_EXCEPT(
591 !Thyra::testRelNormDiffErr(
592 "x", *x,
"x_interp", *x_interp,
593 "2*epsilon", ScalarMag(100.0*SMT::eps()),
594 "2*epsilon", ScalarMag(100.0*SMT::eps()),
595 includesVerbLevel(verbLevel,Teuchos::VERB_HIGH) ? out.get() : 0
599 TEUCHOS_TEST_FOR_EXCEPT(
600 !Thyra::testRelNormDiffErr(
601 "xdot", *xdot,
"xdot_interp", *xdot_interp,
602 "2*epsilon", ScalarMag(100.0*SMT::eps()),
603 "2*epsilon", ScalarMag(100.0*SMT::eps()),
604 includesVerbLevel(verbLevel,Teuchos::VERB_HIGH) ? out.get() : 0
613 #endif // HAVE_RYTHMOS_DEBUG 615 if ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
617 <<
"\nLeaving " << Teuchos::TypeNameTraits<ThetaStepper<Scalar> >::name()
618 <<
"::takeStep(...) ...\n";
626 template<
class Scalar>
627 const StepStatus<Scalar> ThetaStepper<Scalar>::getStepStatus()
const 630 typedef Teuchos::ScalarTraits<Scalar> ST;
632 StepStatus<Scalar> stepStatus;
634 if (!isInitialized_) {
635 stepStatus.stepStatus = STEP_STATUS_UNINITIALIZED;
637 else if (numSteps_ > 0) {
638 stepStatus.stepStatus = STEP_STATUS_CONVERGED;
642 stepStatus.stepSize = dt_;
643 stepStatus.order = 1;
644 stepStatus.time = t_;
645 stepStatus.solution = x_;
646 stepStatus.solutionDot = x_dot_;
656 template<
class Scalar>
657 RCP<const Thyra::VectorSpaceBase<Scalar> >
658 ThetaStepper<Scalar>::get_x_space()
const 660 return ( !is_null(model_) ? model_->get_x_space() : Teuchos::null );
664 template<
class Scalar>
665 void ThetaStepper<Scalar>::addPoints(
666 const Array<Scalar>& time_vec,
667 const Array<RCP<
const Thyra::VectorBase<Scalar> > >& x_vec,
668 const Array<RCP<
const Thyra::VectorBase<Scalar> > >& xdot_vec
672 typedef Teuchos::ScalarTraits<Scalar> ST;
677 #ifdef HAVE_RYTHMOS_DEBUG 678 TEUCHOS_TEST_FOR_EXCEPTION(
679 time_vec.size() == 0, std::logic_error,
680 "Error, addPoints called with an empty time_vec array!\n");
681 #endif // HAVE_RYTHMOS_DEBUG 683 RCP<Teuchos::FancyOStream> out = this->getOStream();
684 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
685 Teuchos::OSTab ostab(out,1,
"TS::setPoints");
687 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
688 *out <<
"time_vec = " << std::endl;
689 for (
int i=0 ; i<Teuchos::as<int>(time_vec.size()) ; ++i) {
690 *out <<
"time_vec[" << i <<
"] = " << time_vec[i] << std::endl;
693 else if (time_vec.size() == 1) {
697 Thyra::V_V(x_.ptr(),*x_vec[n]);
698 Thyra::V_V(x_dot_base_.ptr(),*x_);
701 int n = time_vec.size()-1;
702 int nm1 = time_vec.size()-2;
704 t_old_ = time_vec[nm1];
705 Thyra::V_V(x_.ptr(),*x_vec[n]);
706 Scalar dt = t_ - t_old_;
707 Thyra::V_StV(x_dot_base_.ptr(),Scalar(-ST::one()/dt),*x_vec[nm1]);
712 template<
class Scalar>
713 TimeRange<Scalar> ThetaStepper<Scalar>::getTimeRange()
const 715 if ( !isInitialized_ && haveInitialCondition_ )
716 return timeRange<Scalar>(t_,t_);
717 if ( !isInitialized_ && !haveInitialCondition_ )
718 return invalidTimeRange<Scalar>();
719 return timeRange<Scalar>(t_old_,t_);
723 template<
class Scalar>
724 void ThetaStepper<Scalar>::getPoints(
725 const Array<Scalar>& time_vec,
726 Array<RCP<
const Thyra::VectorBase<Scalar> > >* x_vec,
727 Array<RCP<
const Thyra::VectorBase<Scalar> > >* xdot_vec,
728 Array<ScalarMag>* accuracy_vec
731 using Teuchos::constOptInArg;
733 typedef Teuchos::ScalarTraits<Scalar> ST;
735 TEUCHOS_ASSERT(haveInitialCondition_);
737 RCP<Thyra::VectorBase<Scalar> > x_temp = x_;
738 if (compareTimeValues(t_old_,t_)!=0) {
739 Scalar dt = t_ - t_old_;
740 x_temp = x_dot_base_->clone_v();
741 Thyra::Vt_S(x_temp.ptr(),Scalar(-ST::one()*dt));
743 defaultGetPoints<Scalar>(
744 t_old_, constOptInArg(*x_temp), constOptInArg(*x_dot_old_),
745 t_, constOptInArg(*x_), constOptInArg(*x_dot_),
746 time_vec, ptr(x_vec), ptr(xdot_vec), ptr(accuracy_vec),
747 ptr(interpolator_.get())
856 template<
class Scalar>
857 void ThetaStepper<Scalar>::getNodes(Array<Scalar>* time_vec)
const 861 TEUCHOS_ASSERT( time_vec != NULL );
864 if (!haveInitialCondition_) {
868 time_vec->push_back(t_old_);
870 time_vec->push_back(t_);
873 RCP<Teuchos::FancyOStream> out = this->getOStream();
874 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
875 Teuchos::OSTab ostab(out,1,
"TS::getNodes");
876 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
877 *out << this->description() << std::endl;
878 for (
int i=0 ; i<Teuchos::as<int>(time_vec->size()) ; ++i) {
879 *out <<
"time_vec[" << i <<
"] = " << (*time_vec)[i] << std::endl;
885 template<
class Scalar>
886 void ThetaStepper<Scalar>::removeNodes(Array<Scalar>& time_vec)
890 RCP<Teuchos::FancyOStream> out = this->getOStream();
891 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
892 Teuchos::OSTab ostab(out,1,
"TS::removeNodes");
893 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
894 *out <<
"time_vec = " << std::endl;
895 for (
int i=0 ; i<Teuchos::as<int>(time_vec.size()) ; ++i) {
896 *out <<
"time_vec[" << i <<
"] = " << time_vec[i] << std::endl;
899 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::logic_error,
"Error, removeNodes is not implemented for ThetaStepper at this time.\n");
907 template<
class Scalar>
908 int ThetaStepper<Scalar>::getOrder()
const 910 return (thetaStepperType_==ImplicitEuler) ? 1 : 2;
917 template <
class Scalar>
918 void ThetaStepper<Scalar>::setParameterList(
919 RCP<Teuchos::ParameterList>
const& paramList
922 TEUCHOS_TEST_FOR_EXCEPT(is_null(paramList));
923 paramList->validateParametersAndSetDefaults(*this->getValidParameters());
924 parameterList_ = paramList;
925 Teuchos::readVerboseObjectSublist(&*parameterList_,
this);
927 RCP<ParameterList> pl_theta = Teuchos::sublist(parameterList_, RythmosStepControlSettings_name);
929 std::string thetaStepperTypeString =
930 Teuchos::getParameter<std::string>(*pl_theta, ThetaStepperType_name);
932 if (thetaStepperTypeString ==
"Implicit Euler")
933 thetaStepperType_ = ImplicitEuler;
934 else if (thetaStepperTypeString ==
"Trapezoid")
935 thetaStepperType_ = Trapezoid;
937 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
938 "Value of " << ThetaStepperType_name <<
" = " << thetaStepperTypeString
939 <<
" is invalid for Rythmos::ThetaStepper");
941 default_predictor_order_ =
942 Teuchos::getParameter<int>(*pl_theta, PredictorOrder_name);
944 Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
945 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
946 if ( Teuchos::as<int>(verbLevel) >= Teuchos::as<int>(Teuchos::VERB_HIGH) ) {
947 *out << ThetaStepperType_name <<
" = " << thetaStepperTypeString << std::endl;
952 template <
class Scalar>
953 RCP<Teuchos::ParameterList>
954 ThetaStepper<Scalar>::getNonconstParameterList()
956 return(parameterList_);
960 template <
class Scalar>
961 RCP<Teuchos::ParameterList>
962 ThetaStepper<Scalar>::unsetParameterList()
964 RCP<Teuchos::ParameterList>
965 temp_param_list = parameterList_;
966 parameterList_ = Teuchos::null;
967 return(temp_param_list);
971 template<
class Scalar>
972 RCP<const Teuchos::ParameterList>
973 ThetaStepper<Scalar>::getValidParameters()
const 975 using Teuchos::ParameterList;
977 static RCP<const ParameterList> validPL;
979 if (is_null(validPL)) {
981 RCP<ParameterList> pl_top_level = Teuchos::parameterList();
983 RCP<ParameterList> pl = Teuchos::sublist(pl_top_level, RythmosStepControlSettings_name);
985 pl->set<std::string> ( ThetaStepperType_name, ThetaStepperType_default,
986 "Name of Stepper Type in Theta Stepper" 989 pl->set<
int> ( PredictorOrder_name, PredictorOrder_default,
990 "Order of Predictor in Theta Stepper, can be 0, 1, 2" 993 Teuchos::setupVerboseObjectSublist(&*pl_top_level);
994 validPL = pl_top_level;
1003 template<
class Scalar>
1004 void ThetaStepper<Scalar>::describe(
1005 Teuchos::FancyOStream &out,
1006 const Teuchos::EVerbosityLevel verbLevel
1010 Teuchos::OSTab tab(out);
1011 if (!isInitialized_) {
1012 out << this->description() <<
" : This stepper is not initialized yet" << std::endl;
1016 as<int>(verbLevel) == as<int>(Teuchos::VERB_DEFAULT)
1018 as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW)
1021 out << this->description() <<
"::describe:" << std::endl;
1022 out <<
"model = " << model_->description() << std::endl;
1023 out <<
"solver = " << solver_->description() << std::endl;
1024 if (neModel_ == Teuchos::null) {
1025 out <<
"neModel = Teuchos::null" << std::endl;
1027 out <<
"neModel = " << neModel_->description() << std::endl;
1030 else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW)) {
1031 out <<
"t_ = " << t_ << std::endl;
1032 out <<
"t_old_ = " << t_old_ << std::endl;
1034 else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_MEDIUM)) {
1036 else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH)) {
1037 out <<
"model_ = " << std::endl;
1038 model_->describe(out,verbLevel);
1039 out <<
"solver_ = " << std::endl;
1040 solver_->describe(out,verbLevel);
1041 if (neModel_ == Teuchos::null) {
1042 out <<
"neModel = Teuchos::null" << std::endl;
1044 out <<
"neModel = " << std::endl;
1045 neModel_->describe(out,verbLevel);
1047 out <<
"x_ = " << std::endl;
1048 x_->describe(out,verbLevel);
1049 out <<
"x_dot_base_ = " << std::endl;
1050 x_dot_base_->describe(out,verbLevel);
1058 template <
class Scalar>
1059 void ThetaStepper<Scalar>::initialize_()
1062 typedef Teuchos::ScalarTraits<Scalar> ST;
1067 TEUCHOS_TEST_FOR_EXCEPT(is_null(model_));
1068 TEUCHOS_TEST_FOR_EXCEPT(is_null(solver_));
1069 TEUCHOS_TEST_FOR_EXCEPT(!haveInitialCondition_);
1071 #ifdef HAVE_RYTHMOS_DEBUG 1072 THYRA_ASSERT_VEC_SPACES(
1073 "Rythmos::ThetaStepper::setInitialCondition(...)",
1074 *x_->space(), *model_->get_x_space() );
1075 #endif // HAVE_RYTHMOS_DEBUG 1077 if ( is_null(interpolator_) ) {
1080 interpolator_ = Teuchos::rcp(
new LinearInterpolator<Scalar> );
1085 if (thetaStepperType_ == ImplicitEuler)
1088 predictor_corrector_begin_after_step_ = 2;
1093 predictor_corrector_begin_after_step_ = 3;
1096 isInitialized_ =
true;
1100 template<
class Scalar>
1101 void ThetaStepper<Scalar>::obtainPredictor_()
1105 typedef Teuchos::ScalarTraits<Scalar> ST;
1107 if (!isInitialized_) {
1111 RCP<Teuchos::FancyOStream> out = this->getOStream();
1112 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
1113 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
1114 *out <<
"Obtaining predictor..." << std::endl;
1117 const int preferred_predictor_order = std::min(default_predictor_order_, thetaStepperType_ + 1);
1118 const int max_predictor_order_at_this_timestep = std::max(0, numSteps_);
1120 const int predictor_order = std::min(preferred_predictor_order, max_predictor_order_at_this_timestep);
1122 switch (predictor_order)
1125 V_StV(x_pre_.ptr(), Scalar(ST::one()), *x_old_);
1129 V_StV(x_pre_.ptr(), Scalar(ST::one()), *x_old_);
1131 TEUCHOS_TEST_FOR_EXCEPT (dt_ <= 0.0);
1133 Vp_StV(x_pre_.ptr(), dt_, *x_dot_old_);
1138 V_StV(x_pre_.ptr(), Scalar(ST::one()), *x_old_);
1140 TEUCHOS_TEST_FOR_EXCEPT (dt_ <= 0.0);
1141 TEUCHOS_TEST_FOR_EXCEPT (dt_old_ <= 0.0);
1143 const Scalar coeff_x_dot_old = (0.5 * dt_) * (2.0 + dt_/dt_old_);
1144 const Scalar coeff_x_dot_really_old = - (0.5 * dt_) * (dt_/dt_old_);
1146 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
1147 *out <<
"x_dot_old_ = " << *x_dot_old_ << std::endl;
1148 *out <<
"x_dot_really_old_ = " << *x_dot_really_old_ << std::endl;
1151 Vp_StV( x_pre_.ptr(), coeff_x_dot_old, *x_dot_old_);
1152 Vp_StV( x_pre_.ptr(), coeff_x_dot_really_old, *x_dot_really_old_);
1156 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
1157 "Invalid predictor order " << predictor_order <<
". Valid values are 0, 1, and 2.");
1160 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
1161 *out <<
"x_pre_ = " << *x_pre_ << std::endl;
1165 V_StV(x_.ptr(), Scalar(ST::one()), *x_pre_);
1174 #define RYTHMOS_THETA_STEPPER_INSTANT(SCALAR) \ 1176 template class ThetaStepper< SCALAR >; \ 1178 template RCP< ThetaStepper< SCALAR > > \ 1180 const RCP<Thyra::ModelEvaluator< SCALAR > >& model, \ 1181 const RCP<Thyra::NonlinearSolverBase< SCALAR > >& solver, \ 1182 RCP<Teuchos::ParameterList>& parameterList \ 1188 #endif // HAVE_RYTHMOS_EXPERIMENTAL 1190 #endif //Rythmos_THETA_STEPPER_DEF_H Decorator subclass for a steady-state version of a DAE for single-residual time stepper methods...