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>
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>
118 template<
class Scalar>
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> >
134 return interpolator_;
137 template<
class Scalar>
138 RCP<const InterpolatorBase<Scalar> >
141 return interpolator_;
145 template<
class Scalar>
146 RCP<InterpolatorBase<Scalar> >
149 RCP<InterpolatorBase<Scalar> > temp_interpolator = interpolator_;
150 interpolator_ = Teuchos::null;
151 return(temp_interpolator);
152 isInitialized_ =
false;
159 template<
class Scalar>
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> >
192 template<
class Scalar>
193 RCP<const Thyra::NonlinearSolverBase<Scalar> >
203 template<
class Scalar>
209 template<
class Scalar>
210 RCP<StepperBase<Scalar> >
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>
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>
304 const RCP<Thyra::ModelEvaluator<Scalar> >& model
307 this->setModel(model);
311 template<
class Scalar>
312 RCP<const Thyra::ModelEvaluator<Scalar> >
319 template<
class Scalar>
320 RCP<Thyra::ModelEvaluator<Scalar> >
323 TEUCHOS_TEST_FOR_EXCEPT(
true);
324 return Teuchos::null;
328 template<
class Scalar>
330 const Thyra::ModelEvaluatorBase::InArgs<Scalar> &initialCondition
334 typedef Teuchos::ScalarTraits<Scalar> ST;
336 basePoint_ = initialCondition;
340 RCP<const Thyra::VectorBase<Scalar> >
341 x_init = initialCondition.get_x();
343 #ifdef HAVE_RYTHMOS_DEBUG
344 TEUCHOS_TEST_FOR_EXCEPTION(
345 is_null(x_init), std::logic_error,
346 "Error, if the client passes in an intial condition to setInitialCondition(...),\n"
347 "then x can not be null!" );
350 x_ = x_init->clone_v();
354 RCP<const Thyra::VectorBase<Scalar> >
355 x_dot_init = initialCondition.get_x_dot();
357 if (!is_null(x_dot_init)) {
358 x_dot_ = x_dot_init->clone_v();
361 x_dot_ = createMember(x_->space());
362 assign(x_dot_.ptr(),ST::zero());
367 t_ = initialCondition.get_t();
375 x_pre_ = x_->clone_v();
379 x_old_ = x_->clone_v();
383 x_dot_base_ = x_->clone_v();
387 x_dot_old_ = x_dot_->clone_v();
391 x_dot_really_old_ = x_dot_->clone_v();
393 haveInitialCondition_ =
true;
398 template<
class Scalar>
399 Thyra::ModelEvaluatorBase::InArgs<Scalar>
406 template<
class Scalar>
411 using Teuchos::incrVerbLevel;
412 typedef Teuchos::ScalarTraits<Scalar> ST;
413 typedef Thyra::NonlinearSolverBase<Scalar> NSB;
414 typedef Teuchos::VerboseObjectTempState<NSB> VOTSNSB;
422 RCP<Teuchos::FancyOStream> out = this->getOStream();
423 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
424 Teuchos::OSTab ostab(out,1,
"TS::takeStep");
425 VOTSNSB solver_outputTempState(solver_,out,incrVerbLevel(verbLevel,-1));
427 if ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
429 <<
"\nEntering " << Teuchos::TypeNameTraits<ThetaStepper<Scalar> >::name()
430 <<
"::takeStep("<<dt<<
","<<toString(stepSizeType)<<
") ...\n";
434 V_StV( x_old_.ptr(), Scalar(ST::one()), *x_ );
435 V_StV( x_dot_really_old_.ptr(), Scalar(ST::one()), *x_dot_old_ );
436 V_StV( x_dot_old_.ptr(), Scalar(ST::one()), *x_dot_ );
438 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_EXTREME) ) {
439 *out <<
"\nSetting dt_ and old data ..." << std::endl;
440 *out <<
"\ndt_ = " << dt_;
441 *out <<
"\nx_old_ = " << *x_old_;
442 *out <<
"\nx_dot_old_ = " << *x_dot_old_;
443 *out <<
"\nx_dot_really_old_ = " << *x_dot_really_old_;
446 if ((stepSizeType == STEP_TYPE_VARIABLE) || (dt == ST::zero())) {
447 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) )
448 *out <<
"\nThe arguments to takeStep are not valid for ThetaStepper at this time." << std::endl;
449 return(Scalar(-ST::one()));
451 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
452 *out <<
"\ndt = " << dt << std::endl;
453 *out <<
"\nstepSizeType = " << stepSizeType << std::endl;
469 const double theta = (numSteps_+1>=predictor_corrector_begin_after_step_) ? theta_ : 1.0;
471 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_EXTREME) ) {
472 *out <<
"\nSetting x_dot_base_ ..." << std::endl;
473 *out <<
"\ntheta = " << theta;
474 *out <<
"\nx_ = " << *x_;
475 *out <<
"\nx_dot_old_ = " << *x_dot_old_;
478 const Scalar coeff_x_dot = Scalar(ST::one()/(theta*dt));
479 const Scalar coeff_x = ST::one();
481 const Scalar x_coeff = Scalar(-coeff_x_dot);
482 const Scalar x_dot_old_coeff = Scalar( -(ST::one()-theta)/theta);
484 V_StV( x_dot_base_.ptr(), x_coeff, *x_old_ );
485 Vp_StV( x_dot_base_.ptr(), x_dot_old_coeff, *x_dot_old_);
487 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_EXTREME) ) {
488 *out <<
"\nx_dot_base_ = " << *x_dot_base_;
491 if(!neModel_.get()) {
495 neModel_->initializeSingleResidualModel(
505 if( solver_->getModel().get() != neModel_.get() ) {
506 solver_->setModel(neModel_);
509 solver_->setVerbLevel(this->getVerbLevel());
515 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
516 *out <<
"\nSolving the implicit theta-stepper timestep equation"
517 <<
" with theta = " << theta <<
"\n";
520 Thyra::SolveStatus<Scalar>
521 neSolveStatus = solver_->solve(&*x_);
527 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
528 *out <<
"\nOutput status of nonlinear solve:\n" << neSolveStatus;
542 V_StV( x_dot_.ptr(), Scalar( ST::one()/(theta*dt)), *x_ );
543 Vp_StV( x_dot_.ptr(), Scalar(-ST::one()/(theta*dt)), *x_old_ );
544 Vp_StV( x_dot_.ptr(), Scalar( -(ST::one()-theta)/theta), *x_dot_old_ );
546 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_EXTREME) ) {
547 *out <<
"\nUpdating x_dot_ ...\n";
548 *out <<
"\nx_dot_ = " << *x_dot_ << std::endl;
557 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
558 *out <<
"\nt_old_ = " << t_old_ << std::endl;
559 *out <<
"\nt_ = " << t_ << std::endl;
562 #ifdef HAVE_RYTHMOS_DEBUG
564 if ( includesVerbLevel(verbLevel,Teuchos::VERB_LOW) )
565 *out <<
"\nChecking to make sure that solution and the interpolated solution are the same! ...\n";
569 typedef typename ST::magnitudeType ScalarMag;
570 typedef ScalarTraits<ScalarMag> SMT;
572 Teuchos::OSTab tab(out);
574 const StepStatus<Scalar> stepStatus = this->getStepStatus();
576 RCP<const Thyra::VectorBase<Scalar> >
577 x = stepStatus.solution,
578 xdot = stepStatus.solutionDot;
580 Array<Scalar> time_vec = Teuchos::tuple(stepStatus.time);
581 Array<RCP<const Thyra::VectorBase<Scalar> > > x_vec, xdot_vec;
582 this->getPoints(time_vec,&x_vec,&xdot_vec,0);
584 RCP<const Thyra::VectorBase<Scalar> >
586 xdot_interp = xdot_vec[0];
588 TEUCHOS_TEST_FOR_EXCEPT(
589 !Thyra::testRelNormDiffErr(
590 "x", *x,
"x_interp", *x_interp,
591 "2*epsilon", ScalarMag(100.0*SMT::eps()),
592 "2*epsilon", ScalarMag(100.0*SMT::eps()),
593 includesVerbLevel(verbLevel,Teuchos::VERB_HIGH) ? out.get() : 0
597 TEUCHOS_TEST_FOR_EXCEPT(
598 !Thyra::testRelNormDiffErr(
599 "xdot", *xdot,
"xdot_interp", *xdot_interp,
600 "2*epsilon", ScalarMag(100.0*SMT::eps()),
601 "2*epsilon", ScalarMag(100.0*SMT::eps()),
602 includesVerbLevel(verbLevel,Teuchos::VERB_HIGH) ? out.get() : 0
611 #endif // HAVE_RYTHMOS_DEBUG
613 if ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
615 <<
"\nLeaving " << Teuchos::TypeNameTraits<ThetaStepper<Scalar> >::name()
616 <<
"::takeStep(...) ...\n";
624 template<
class Scalar>
628 StepStatus<Scalar> stepStatus;
630 if (!isInitialized_) {
631 stepStatus.stepStatus = STEP_STATUS_UNINITIALIZED;
633 else if (numSteps_ > 0) {
634 stepStatus.stepStatus = STEP_STATUS_CONVERGED;
638 stepStatus.stepSize = dt_;
639 stepStatus.order = 1;
640 stepStatus.time = t_;
641 stepStatus.solution = x_;
642 stepStatus.solutionDot = x_dot_;
652 template<
class Scalar>
653 RCP<const Thyra::VectorSpaceBase<Scalar> >
656 return ( !is_null(model_) ? model_->get_x_space() : Teuchos::null );
660 template<
class Scalar>
662 const Array<Scalar>& time_vec,
663 const Array<RCP<
const Thyra::VectorBase<Scalar> > >& x_vec,
664 const Array<RCP<
const Thyra::VectorBase<Scalar> > >&
668 typedef Teuchos::ScalarTraits<Scalar> ST;
673 #ifdef HAVE_RYTHMOS_DEBUG
674 TEUCHOS_TEST_FOR_EXCEPTION(
675 time_vec.size() == 0, std::logic_error,
676 "Error, addPoints called with an empty time_vec array!\n");
677 #endif // HAVE_RYTHMOS_DEBUG
679 RCP<Teuchos::FancyOStream> out = this->getOStream();
680 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
681 Teuchos::OSTab ostab(out,1,
"TS::setPoints");
683 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
684 *out <<
"time_vec = " << std::endl;
685 for (
int i=0 ; i<Teuchos::as<int>(time_vec.size()) ; ++i) {
686 *out <<
"time_vec[" << i <<
"] = " << time_vec[i] << std::endl;
689 else if (time_vec.size() == 1) {
693 Thyra::V_V(x_.ptr(),*x_vec[n]);
694 Thyra::V_V(x_dot_base_.ptr(),*x_);
697 int n = time_vec.size()-1;
698 int nm1 = time_vec.size()-2;
700 t_old_ = time_vec[nm1];
701 Thyra::V_V(x_.ptr(),*x_vec[n]);
702 Scalar dt = t_ - t_old_;
703 Thyra::V_StV(x_dot_base_.ptr(),Scalar(-ST::one()/dt),*x_vec[nm1]);
708 template<
class Scalar>
711 if ( !isInitialized_ && haveInitialCondition_ )
712 return timeRange<Scalar>(t_,t_);
713 if ( !isInitialized_ && !haveInitialCondition_ )
714 return invalidTimeRange<Scalar>();
715 return timeRange<Scalar>(t_old_,t_);
719 template<
class Scalar>
721 const Array<Scalar>& time_vec,
722 Array<RCP<
const Thyra::VectorBase<Scalar> > >* x_vec,
723 Array<RCP<
const Thyra::VectorBase<Scalar> > >* xdot_vec,
724 Array<ScalarMag>* accuracy_vec
727 using Teuchos::constOptInArg;
729 typedef Teuchos::ScalarTraits<Scalar> ST;
731 TEUCHOS_ASSERT(haveInitialCondition_);
733 RCP<Thyra::VectorBase<Scalar> > x_temp = x_;
734 if (compareTimeValues(t_old_,t_)!=0) {
735 Scalar dt = t_ - t_old_;
736 x_temp = x_dot_base_->clone_v();
737 Thyra::Vt_S(x_temp.ptr(),Scalar(-ST::one()*dt));
739 defaultGetPoints<Scalar>(
740 t_old_, constOptInArg(*x_temp), constOptInArg(*x_dot_old_),
741 t_, constOptInArg(*x_), constOptInArg(*x_dot_),
742 time_vec, ptr(x_vec), ptr(xdot_vec), ptr(accuracy_vec),
743 ptr(interpolator_.get())
852 template<
class Scalar>
857 TEUCHOS_ASSERT( time_vec != NULL );
860 if (!haveInitialCondition_) {
864 time_vec->push_back(t_old_);
866 time_vec->push_back(t_);
869 RCP<Teuchos::FancyOStream> out = this->getOStream();
870 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
871 Teuchos::OSTab ostab(out,1,
"TS::getNodes");
872 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
873 *out << this->description() << std::endl;
874 for (
int i=0 ; i<Teuchos::as<int>(time_vec->size()) ; ++i) {
875 *out <<
"time_vec[" << i <<
"] = " << (*time_vec)[i] << std::endl;
881 template<
class Scalar>
886 RCP<Teuchos::FancyOStream> out = this->getOStream();
887 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
888 Teuchos::OSTab ostab(out,1,
"TS::removeNodes");
889 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
890 *out <<
"time_vec = " << std::endl;
891 for (
int i=0 ; i<Teuchos::as<int>(time_vec.size()) ; ++i) {
892 *out <<
"time_vec[" << i <<
"] = " << time_vec[i] << std::endl;
895 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::logic_error,
"Error, removeNodes is not implemented for ThetaStepper at this time.\n");
903 template<
class Scalar>
906 return (thetaStepperType_==ImplicitEuler) ? 1 : 2;
913 template <
class Scalar>
915 RCP<Teuchos::ParameterList>
const& paramList
918 TEUCHOS_TEST_FOR_EXCEPT(is_null(paramList));
919 paramList->validateParametersAndSetDefaults(*this->getValidParameters());
920 parameterList_ = paramList;
921 Teuchos::readVerboseObjectSublist(&*parameterList_,
this);
923 RCP<ParameterList> pl_theta = Teuchos::sublist(parameterList_, RythmosStepControlSettings_name);
925 std::string thetaStepperTypeString =
926 Teuchos::getParameter<std::string>(*pl_theta, ThetaStepperType_name);
928 if (thetaStepperTypeString ==
"Implicit Euler")
929 thetaStepperType_ = ImplicitEuler;
930 else if (thetaStepperTypeString ==
"Trapezoid")
931 thetaStepperType_ = Trapezoid;
933 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
934 "Value of " << ThetaStepperType_name <<
" = " << thetaStepperTypeString
935 <<
" is invalid for Rythmos::ThetaStepper");
937 default_predictor_order_ =
938 Teuchos::getParameter<int>(*pl_theta, PredictorOrder_name);
940 Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
941 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
942 if ( Teuchos::as<int>(verbLevel) >= Teuchos::as<int>(Teuchos::VERB_HIGH) ) {
943 *out << ThetaStepperType_name <<
" = " << thetaStepperTypeString << std::endl;
948 template <
class Scalar>
949 RCP<Teuchos::ParameterList>
952 return(parameterList_);
956 template <
class Scalar>
957 RCP<Teuchos::ParameterList>
960 RCP<Teuchos::ParameterList>
961 temp_param_list = parameterList_;
962 parameterList_ = Teuchos::null;
963 return(temp_param_list);
967 template<
class Scalar>
968 RCP<const Teuchos::ParameterList>
971 using Teuchos::ParameterList;
973 static RCP<const ParameterList> validPL;
975 if (is_null(validPL)) {
977 RCP<ParameterList> pl_top_level = Teuchos::parameterList();
979 RCP<ParameterList> pl = Teuchos::sublist(pl_top_level, RythmosStepControlSettings_name);
981 pl->set<std::string> ( ThetaStepperType_name, ThetaStepperType_default,
982 "Name of Stepper Type in Theta Stepper"
985 pl->set<
int> ( PredictorOrder_name, PredictorOrder_default,
986 "Order of Predictor in Theta Stepper, can be 0, 1, 2"
989 Teuchos::setupVerboseObjectSublist(&*pl_top_level);
990 validPL = pl_top_level;
999 template<
class Scalar>
1001 Teuchos::FancyOStream &out,
1002 const Teuchos::EVerbosityLevel verbLevel
1006 Teuchos::OSTab tab(out);
1007 if (!isInitialized_) {
1008 out << this->description() <<
" : This stepper is not initialized yet" << std::endl;
1012 as<int>(verbLevel) == as<int>(Teuchos::VERB_DEFAULT)
1014 as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW)
1017 out << this->description() <<
"::describe:" << std::endl;
1018 out <<
"model = " << model_->description() << std::endl;
1019 out <<
"solver = " << solver_->description() << std::endl;
1020 if (neModel_ == Teuchos::null) {
1021 out <<
"neModel = Teuchos::null" << std::endl;
1023 out <<
"neModel = " << neModel_->description() << std::endl;
1026 else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW)) {
1027 out <<
"t_ = " << t_ << std::endl;
1028 out <<
"t_old_ = " << t_old_ << std::endl;
1030 else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_MEDIUM)) {
1032 else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH)) {
1033 out <<
"model_ = " << std::endl;
1034 model_->describe(out,verbLevel);
1035 out <<
"solver_ = " << std::endl;
1036 solver_->describe(out,verbLevel);
1037 if (neModel_ == Teuchos::null) {
1038 out <<
"neModel = Teuchos::null" << std::endl;
1040 out <<
"neModel = " << std::endl;
1041 neModel_->describe(out,verbLevel);
1043 out <<
"x_ = " << std::endl;
1044 x_->describe(out,verbLevel);
1045 out <<
"x_dot_base_ = " << std::endl;
1046 x_dot_base_->describe(out,verbLevel);
1054 template <
class Scalar>
1055 void ThetaStepper<Scalar>::initialize_()
1061 TEUCHOS_TEST_FOR_EXCEPT(is_null(model_));
1062 TEUCHOS_TEST_FOR_EXCEPT(is_null(solver_));
1063 TEUCHOS_TEST_FOR_EXCEPT(!haveInitialCondition_);
1065 #ifdef HAVE_RYTHMOS_DEBUG
1066 THYRA_ASSERT_VEC_SPACES(
1067 "Rythmos::ThetaStepper::setInitialCondition(...)",
1068 *x_->space(), *model_->get_x_space() );
1069 #endif // HAVE_RYTHMOS_DEBUG
1071 if ( is_null(interpolator_) ) {
1074 interpolator_ = Teuchos::rcp(
new LinearInterpolator<Scalar> );
1079 if (thetaStepperType_ == ImplicitEuler)
1082 predictor_corrector_begin_after_step_ = 2;
1087 predictor_corrector_begin_after_step_ = 3;
1090 isInitialized_ =
true;
1094 template<
class Scalar>
1095 void ThetaStepper<Scalar>::obtainPredictor_()
1099 typedef Teuchos::ScalarTraits<Scalar> ST;
1101 if (!isInitialized_) {
1105 RCP<Teuchos::FancyOStream> out = this->getOStream();
1106 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
1107 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
1108 *out <<
"Obtaining predictor..." << std::endl;
1111 const int preferred_predictor_order = std::min(default_predictor_order_, thetaStepperType_ + 1);
1112 const int max_predictor_order_at_this_timestep = std::max(0, numSteps_);
1114 const int predictor_order = std::min(preferred_predictor_order, max_predictor_order_at_this_timestep);
1116 switch (predictor_order)
1119 V_StV(x_pre_.ptr(), Scalar(ST::one()), *x_old_);
1123 V_StV(x_pre_.ptr(), Scalar(ST::one()), *x_old_);
1125 TEUCHOS_TEST_FOR_EXCEPT (dt_ <= 0.0);
1127 Vp_StV(x_pre_.ptr(), dt_, *x_dot_old_);
1132 V_StV(x_pre_.ptr(), Scalar(ST::one()), *x_old_);
1134 TEUCHOS_TEST_FOR_EXCEPT (dt_ <= 0.0);
1135 TEUCHOS_TEST_FOR_EXCEPT (dt_old_ <= 0.0);
1137 const Scalar coeff_x_dot_old = (0.5 * dt_) * (2.0 + dt_/dt_old_);
1138 const Scalar coeff_x_dot_really_old = - (0.5 * dt_) * (dt_/dt_old_);
1140 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
1141 *out <<
"x_dot_old_ = " << *x_dot_old_ << std::endl;
1142 *out <<
"x_dot_really_old_ = " << *x_dot_really_old_ << std::endl;
1145 Vp_StV( x_pre_.ptr(), coeff_x_dot_old, *x_dot_old_);
1146 Vp_StV( x_pre_.ptr(), coeff_x_dot_really_old, *x_dot_really_old_);
1150 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
1151 "Invalid predictor order " << predictor_order <<
". Valid values are 0, 1, and 2.");
1154 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
1155 *out <<
"x_pre_ = " << *x_pre_ << std::endl;
1159 V_StV(x_.ptr(), Scalar(ST::one()), *x_pre_);
1168 #define RYTHMOS_THETA_STEPPER_INSTANT(SCALAR) \
1170 template class ThetaStepper< SCALAR >; \
1172 template RCP< ThetaStepper< SCALAR > > \
1174 const RCP<Thyra::ModelEvaluator< SCALAR > >& model, \
1175 const RCP<Thyra::NonlinearSolverBase< SCALAR > >& solver, \
1176 RCP<Teuchos::ParameterList>& parameterList \
1182 #endif // HAVE_RYTHMOS_EXPERIMENTAL
1184 #endif //Rythmos_THETA_STEPPER_DEF_H
bool supportsCloning() const
Returns true.
RCP< const Teuchos::ParameterList > getValidParameters() const
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
RCP< InterpolatorBase< Scalar > > getNonconstInterpolator()
void removeNodes(Array< Scalar > &time_vec)
RCP< Teuchos::ParameterList > unsetParameterList()
void setInterpolator(const RCP< InterpolatorBase< Scalar > > &interpolator)
Redefined from InterpolatorAcceptingObjectBase.
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
RCP< const Thyra::NonlinearSolverBase< Scalar > > getSolver() const
RCP< const InterpolatorBase< Scalar > > getInterpolator() 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)
void setModel(const RCP< const Thyra::ModelEvaluator< Scalar > > &model)
RCP< const Thyra::VectorSpaceBase< Scalar > > get_x_space() const
void setParameterList(RCP< Teuchos::ParameterList > const ¶mList)
void setInitialCondition(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &initialCondition)
void setSolver(const RCP< Thyra::NonlinearSolverBase< Scalar > > &solver)
void getNodes(Array< Scalar > *time_vec) const
Decorator subclass for a steady-state version of a DAE for single-residual time stepper methods...
bool isImplicit() const
Return if this stepper is an implicit stepper.
RCP< Thyra::ModelEvaluator< Scalar > > getNonconstModel()
RCP< InterpolatorBase< Scalar > > unSetInterpolator()
TimeRange< Scalar > getTimeRange() const
RCP< Thyra::NonlinearSolverBase< Scalar > > getNonconstSolver()
Thyra::ModelEvaluatorBase::InArgs< Scalar > getInitialCondition() const
const StepStatus< Scalar > getStepStatus() const
Scalar takeStep(Scalar dt, StepSizeType flag)
RCP< StepperBase< Scalar > > cloneStepperAlgorithm() const
Creates copies of all internal data (including the parameter list) except the model which is assumed ...
RCP< const Thyra::ModelEvaluator< Scalar > > getModel() const
void setNonconstModel(const RCP< Thyra::ModelEvaluator< Scalar > > &model)