223 RCP<BackwardEulerStepper<Scalar> >
224 stepper = Teuchos::rcp(
new BackwardEulerStepper<Scalar>);
225 stepper->isInitialized_ = isInitialized_;
226 stepper->model_ = model_;
227 if (!is_null(solver_))
228 stepper->solver_ = solver_->cloneNonlinearSolver().assert_not_null();
230 stepper->x_ = x_->clone_v().assert_not_null();
231 if (!is_null(scaled_x_old_))
232 stepper->scaled_x_old_ = scaled_x_old_->clone_v().assert_not_null();
234 stepper->t_old_ = t_old_;
236 stepper->numSteps_ = numSteps_;
237 if (!is_null(neModel_))
239 = Teuchos::rcp(
new Rythmos::SingleResidualModelEvaluator<Scalar>);
240 if (!is_null(parameterList_))
241 stepper->parameterList_ = parameterList(*parameterList_);
242 if (!is_null(interpolator_))
243 stepper->interpolator_
244 = interpolator_->cloneInterpolator().assert_not_null();
245 if (!is_null(stepControl_)) {
246 if (stepControl_->supportsCloning())
247 stepper->setStepControlStrategy(
248 stepControl_->cloneStepControlStrategyAlgorithm().assert_not_null());
303 const Thyra::ModelEvaluatorBase::InArgs<Scalar> &initialCondition)
306 typedef Teuchos::ScalarTraits<Scalar> ST;
309 basePoint_ = initialCondition;
312 RCP<const Thyra::VectorBase<Scalar> > x_init = initialCondition.
get_x();
313 TEUCHOS_TEST_FOR_EXCEPTION(
314 is_null(x_init), std::logic_error,
315 "Error, if the client passes in an intial condition to "
316 "setInitialCondition(...), then x can not be null!" );
317 x_ = x_init->clone_v();
320 RCP<const Thyra::VectorBase<Scalar> >
321 x_dot_init = initialCondition.get_x_dot();
322 if (!is_null(x_dot_init)) {
323 x_dot_ = x_dot_init->clone_v();
326 x_dot_ = createMember(x_->space());
327 V_S(x_dot_.ptr(),ST::zero());
331 if (is_null(dx_)) dx_ = createMember(x_->space());
332 V_S(dx_.ptr(), ST::zero());
335 t_ = initialCondition.get_t();
339 if (is_null(x_old_)) x_old_ = createMember(x_->space());
340 x_old_ = x_init->clone_v();
341 scaled_x_old_ = x_->clone_v();
344 if (is_null(x_dot_old_)) x_dot_old_ = createMember(x_->space());
345 x_dot_old_ = x_dot_->clone_v();
347 haveInitialCondition_ =
true;
383 StepSizeType stepSizeType)
387 using Teuchos::incrVerbLevel;
388 typedef Teuchos::ScalarTraits<Scalar> ST;
389 typedef Thyra::NonlinearSolverBase<Scalar> NSB;
390 typedef Teuchos::VerboseObjectTempState<NSB> VOTSNSB;
392 RCP<Teuchos::FancyOStream> out = this->getOStream();
393 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
394 Teuchos::OSTab ostab(out,1,
"BES::takeStep");
395 VOTSNSB solver_outputTempState(solver_,out,incrVerbLevel(verbLevel,-1));
397 if ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
398 *out <<
"\nEntering "
399 << Teuchos::TypeNameTraits<BackwardEulerStepper<Scalar> >::name()
400 <<
"::takeStep("<<dt<<
","<<toString(stepSizeType)<<
") ...\n";
405 if(!neModel_.get()) {
406 neModel_ =Teuchos::rcp(
new Rythmos::SingleResidualModelEvaluator<Scalar>());
411 if (dt <= ST::zero()) {
412 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) )
413 *out <<
"\nThe arguments to takeStep are not valid for "
414 <<
"BackwardEulerStepper at this time.\n"
415 <<
" dt = " << dt <<
"\n"
416 <<
"BackwardEulerStepper requires positive dt.\n" << std::endl;
417 return(Scalar(-ST::one()));
419 if ((stepSizeType == STEP_TYPE_VARIABLE) && (stepControl_ == Teuchos::null)) {
420 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) )
421 *out <<
"\nFor 'variable' time stepping (step size adjustable by "
422 <<
"BackwardEulerStepper), BackwardEulerStepper requires "
423 <<
"Step-Control Strategy.\n"
424 <<
" stepType = " << toString(stepSizeType) <<
"\n" << std::endl;
425 return(Scalar(-ST::one()));
428 stepControl_->setRequestedStepSize(*
this,dt_,stepSizeType);
429 AttemptedStepStatusFlag status;
430 bool stepPass =
false;
433 stepControl_->nextStepSize(*
this,&dt_,&stepSizeType,NULL);
434 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
435 *out <<
"\nrequested dt = " << dt
436 <<
"\ncurrent dt = " << dt_ <<
"\n";
443 V_StV( scaled_x_old_.ptr(), Scalar(-ST::one()/dt_), *x_old_ );
445 neModel_->initializeSingleResidualModel(
447 Scalar(ST::one()/dt_), scaled_x_old_,
448 ST::one(), Teuchos::null,
452 if( solver_->getModel().get() != neModel_.get() ) {
453 solver_->setModel(neModel_);
462 if ( as<int>(verbLevel) > as<int>(Teuchos::VERB_LOW) ) {
463 *out <<
"\nSolving the implicit backward-Euler timestep equation ...\n";
466 Thyra::SolveStatus<Scalar> neSolveStatus =
467 solver_->solve(&*x_, NULL, &*dx_);
474 if ( as<int>(verbLevel) > as<int>(Teuchos::VERB_LOW) ) {
475 *out <<
"\nOutput status of nonlinear solve:\n" << neSolveStatus;
483 if (neSolveStatus.solveStatus == Thyra::SOLVE_STATUS_CONVERGED) {
484 newtonConvergenceStatus_ = 0;
487 newtonConvergenceStatus_ = -1;
490 stepControl_->setCorrection(*
this,x_,dx_,newtonConvergenceStatus_);
492 stepPass = stepControl_->acceptStep(*
this,NULL);
495 status = stepControl_->rejectStep(*
this);
497 if (status != PREDICT_AGAIN)
507 V_V( x_dot_old_.ptr(), *x_dot_ );
509 V_StVpStV( x_dot_.ptr(), Scalar(ST::one()/dt_), *x_,
510 Scalar(-ST::one()/dt_), *x_old_ );
511 V_V( x_old_.ptr(), *x_ );
514 stepControl_->completeStep(*
this);
517 dt_ = Scalar(-ST::one());
520 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
521 *out <<
"\nt_old_ = " << t_old_ << std::endl;
522 *out <<
"\nt_ = " << t_ << std::endl;
525#ifdef HAVE_RYTHMOS_DEBUG
528 if ( includesVerbLevel(verbLevel,Teuchos::VERB_LOW) )
529 *out <<
"\nChecking to make sure that solution and "
530 <<
"the interpolated solution are the same! ...\n";
534 typedef ScalarTraits<ScalarMag> SMT;
536 Teuchos::OSTab tab(out);
538 const StepStatus<Scalar> stepStatus = this->getStepStatus();
540 RCP<const Thyra::VectorBase<Scalar> >
541 x = stepStatus.solution,
542 xdot = stepStatus.solutionDot;
544 Array<Scalar> time_vec = Teuchos::tuple(stepStatus.time);
545 Array<RCP<const Thyra::VectorBase<Scalar> > > x_vec, xdot_vec;
546 this->getPoints(time_vec,&x_vec,&xdot_vec,0);
548 RCP<const Thyra::VectorBase<Scalar> >
550 xdot_interp = xdot_vec[0];
552 TEUCHOS_TEST_FOR_EXCEPT(
553 !Thyra::testRelNormDiffErr(
554 "x", *x,
"x_interp", *x_interp,
555 "2*epsilon",
ScalarMag(100.0*SMT::eps()),
556 "2*epsilon",
ScalarMag(100.0*SMT::eps()),
557 includesVerbLevel(verbLevel,Teuchos::VERB_HIGH) ? out.get() : 0
561 TEUCHOS_TEST_FOR_EXCEPT(
562 !Thyra::testRelNormDiffErr(
563 "xdot", *xdot,
"xdot_interp", *xdot_interp,
564 "2*epsilon",
ScalarMag(100.0*SMT::eps()),
565 "2*epsilon",
ScalarMag(100.0*SMT::eps()),
566 includesVerbLevel(verbLevel,Teuchos::VERB_HIGH) ? out.get() : 0
577 if ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
579 << Teuchos::TypeNameTraits<BackwardEulerStepper<Scalar> >::name()
580 <<
"::takeStep("<<dt_<<
","<<toString(stepSizeType)<<
") ...\n";
627 const Array<Scalar>& time_vec,
628 const Array<RCP<
const Thyra::VectorBase<Scalar> > >& x_vec,
629 const Array<RCP<
const Thyra::VectorBase<Scalar> > >&
633 typedef Teuchos::ScalarTraits<Scalar> ST;
638#ifdef HAVE_RYTHMOS_DEBUG
639 TEUCHOS_TEST_FOR_EXCEPTION(
640 time_vec.size() == 0, std::logic_error,
641 "Error, addPoints called with an empty time_vec array!\n");
644 RCP<Teuchos::FancyOStream> out = this->getOStream();
645 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
646 Teuchos::OSTab ostab(out,1,
"BES::setPoints");
648 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
649 *out <<
"time_vec = " << std::endl;
650 for (
int i=0 ; i<Teuchos::as<int>(time_vec.size()) ; ++i) {
651 *out <<
"time_vec[" << i <<
"] = " << time_vec[i] << std::endl;
654 else if (time_vec.size() == 1) {
658 Thyra::V_V(x_.ptr(),*x_vec[n]);
659 Thyra::V_V(scaled_x_old_.ptr(),*x_);
662 int n = time_vec.size()-1;
663 int nm1 = time_vec.size()-2;
665 t_old_ = time_vec[nm1];
666 Thyra::V_V(x_.ptr(),*x_vec[n]);
667 Scalar dt = t_ - t_old_;
668 Thyra::V_StV(scaled_x_old_.ptr(),Scalar(-ST::one()/dt),*x_vec[nm1]);
686 const Array<Scalar>& time_vec,
687 Array<RCP<
const Thyra::VectorBase<Scalar> > >* x_vec,
688 Array<RCP<
const Thyra::VectorBase<Scalar> > >* xdot_vec,
689 Array<ScalarMag>* accuracy_vec
692 typedef Teuchos::ScalarTraits<Scalar> ST;
693 using Teuchos::constOptInArg;
695#ifdef HAVE_RYTHMOS_DEBUG
696 TEUCHOS_ASSERT(haveInitialCondition_);
698 RCP<Thyra::VectorBase<Scalar> > x_temp = x_;
699 if (compareTimeValues(t_old_,t_)!=0) {
700 Scalar dt = t_ - t_old_;
701 x_temp = scaled_x_old_->clone_v();
702 Thyra::Vt_S(x_temp.ptr(),Scalar(-ST::one()*dt));
704 defaultGetPoints<Scalar>(
705 t_old_, constOptInArg(*x_temp), constOptInArg(*x_dot_old_),
706 t_, constOptInArg(*x_), constOptInArg(*x_dot_),
707 time_vec, ptr(x_vec), ptr(xdot_vec), ptr(accuracy_vec),
708 ptr(interpolator_.get())
928 Teuchos::FancyOStream &out,
929 const Teuchos::EVerbosityLevel verbLevel
933 Teuchos::OSTab tab(out);
934 if (!isInitialized_) {
935 out << this->description() <<
" : This stepper is not initialized yet"
940 as<int>(verbLevel) == as<int>(Teuchos::VERB_DEFAULT)
942 as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW)
945 out << this->description() <<
"::describe:" << std::endl;
946 out <<
"model = " << model_->description() << std::endl;
947 out <<
"solver = " << solver_->description() << std::endl;
948 if (neModel_ == Teuchos::null) {
949 out <<
"neModel = Teuchos::null" << std::endl;
951 out <<
"neModel = " << neModel_->description() << std::endl;
954 else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW)) {
955 out <<
"t_ = " << t_ << std::endl;
956 out <<
"t_old_ = " << t_old_ << std::endl;
958 else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_MEDIUM)) {
960 else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH)) {
961 out <<
"model_ = " << std::endl;
962 model_->describe(out,verbLevel);
963 out <<
"solver_ = " << std::endl;
964 solver_->describe(out,verbLevel);
965 if (neModel_ == Teuchos::null) {
966 out <<
"neModel = Teuchos::null" << std::endl;
968 out <<
"neModel = " << std::endl;
969 neModel_->describe(out,verbLevel);
971 out <<
"x_ = " << std::endl;
972 x_->describe(out,verbLevel);
973 out <<
"scaled_x_old_ = " << std::endl;
974 scaled_x_old_->describe(out,verbLevel);
1036 RCP<BackwardEulerStepperMomento<Scalar> > momento =
1037 Teuchos::rcp(
new BackwardEulerStepperMomento<Scalar>());
1038 momento->set_scaled_x_old(scaled_x_old_);
1039 momento->set_x_old(x_old_);
1040 momento->set_x_dot_old(x_dot_old_);
1042 momento->set_x_dot(x_dot_);
1043 momento->set_dx(dx_);
1045 momento->set_t_old(t_old_);
1046 momento->set_dt(dt_);
1047 momento->set_numSteps(numSteps_);
1048 momento->set_newtonConvergenceStatus(newtonConvergenceStatus_);
1049 momento->set_isInitialized(isInitialized_);
1050 momento->set_haveInitialCondition(haveInitialCondition_);
1051 momento->set_parameterList(parameterList_);
1052 momento->set_basePoint(basePoint_);
1053 momento->set_neModel(neModel_);
1054 momento->set_interpolator(interpolator_);
1055 momento->set_stepControl(stepControl_);
1061 const Ptr<
const MomentoBase<Scalar> >& momentoPtr,
1062 const RCP<Thyra::ModelEvaluator<Scalar> >& model,
1063 const RCP<Thyra::NonlinearSolverBase<Scalar> >& solver
1066 Ptr<const BackwardEulerStepperMomento<Scalar> > beMomentoPtr =
1067 Teuchos::ptr_dynamic_cast<const BackwardEulerStepperMomento<Scalar> >
1069 const BackwardEulerStepperMomento<Scalar>& beMomento = *beMomentoPtr;
1072 scaled_x_old_ = beMomento.get_scaled_x_old();
1073 x_old_ = beMomento.get_x_old();
1074 x_dot_old_ = beMomento.get_x_dot_old();
1075 x_ = beMomento.
get_x();
1076 x_dot_ = beMomento.get_x_dot();
1077 dx_ = beMomento.get_dx();
1078 t_ = beMomento.get_t();
1079 t_old_ = beMomento.get_t_old();
1080 dt_ = beMomento.get_dt();
1081 numSteps_ = beMomento.get_numSteps();
1082 newtonConvergenceStatus_ = beMomento.get_newtonConvergenceStatus();
1083 isInitialized_ = beMomento.get_isInitialized();
1084 haveInitialCondition_ = beMomento.get_haveInitialCondition();
1085 parameterList_ = beMomento.get_parameterList();
1086 basePoint_ = beMomento.get_basePoint();
1087 neModel_ = beMomento.get_neModel();
1088 interpolator_ = beMomento.get_interpolator();
1089 stepControl_ = beMomento.get_stepControl();
1090 this->checkConsistentState_();