Rythmos - Transient Integration for Differential Equations Version of the Day
Rythmos_BackwardEulerStepper_def.hpp
00001 //@HEADER
00002 // ***********************************************************************
00003 //
00004 //                           Rythmos Package
00005 //                 Copyright (2006) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
00023 // USA
00024 // Questions? Contact Todd S. Coffey (tscoffe@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 //@HEADER
00028 
00029 #ifndef Rythmos_BACKWARD_EULER_STEPPER_DEF_H
00030 #define Rythmos_BACKWARD_EULER_STEPPER_DEF_H
00031 
00032 #include "Rythmos_BackwardEulerStepper_decl.hpp"
00033 #include "Rythmos_DataStore.hpp"
00034 #include "Rythmos_LinearInterpolator.hpp"
00035 #include "Rythmos_InterpolatorBaseHelpers.hpp"
00036 #include "Rythmos_StepperHelpers.hpp"
00037 #include "Rythmos_FixedStepControlStrategy.hpp"
00038 #include "Rythmos_SimpleStepControlStrategy.hpp"
00039 #include "Rythmos_FirstOrderErrorStepControlStrategy.hpp"
00040 
00041 #include "Thyra_ModelEvaluatorHelpers.hpp"
00042 #include "Thyra_AssertOp.hpp"
00043 #include "Thyra_TestingTools.hpp"
00044 
00045 #include "Teuchos_VerboseObjectParameterListHelpers.hpp"
00046 #include "Teuchos_as.hpp"
00047 
00048 namespace Rythmos {
00049 
00050 
00051 template<class Scalar>
00052 RCP<BackwardEulerStepper<Scalar> >
00053 backwardEulerStepper(
00054   const RCP<Thyra::ModelEvaluator<Scalar> >& model,
00055   const RCP<Thyra::NonlinearSolverBase<Scalar> >& solver
00056   )
00057 {
00058   return Teuchos::rcp(new BackwardEulerStepper<Scalar>(model, solver));
00059 }
00060 
00061 template<class Scalar>
00062 RCP<BackwardEulerStepper<Scalar> >
00063 backwardEulerStepper()
00064 {
00065   return Teuchos::rcp(new BackwardEulerStepper<Scalar>());
00066 }
00067 
00068 
00069 // ////////////////////////////
00070 // Defintions
00071 
00072 
00073 // Constructors, intializers, Misc.
00074 
00075 
00076 template<class Scalar>
00077 BackwardEulerStepper<Scalar>::BackwardEulerStepper()
00078 {
00079   typedef Teuchos::ScalarTraits<Scalar> ST;
00080   this->defaultInitializeAll_();
00081   numSteps_ = 0;
00082   dt_ = ST::zero();
00083 }
00084 
00085 
00086 template<class Scalar>
00087 BackwardEulerStepper<Scalar>::BackwardEulerStepper(
00088   const RCP<Thyra::ModelEvaluator<Scalar> >& model,
00089   const RCP<Thyra::NonlinearSolverBase<Scalar> >& solver
00090   )
00091 {
00092   typedef Teuchos::ScalarTraits<Scalar> ST;
00093   this->defaultInitializeAll_();
00094   dt_ = ST::zero();
00095   numSteps_ = 0;
00096   setModel(model);
00097   setSolver(solver);
00098 }
00099 
00100 template<class Scalar>
00101 void BackwardEulerStepper<Scalar>::defaultInitializeAll_()
00102 {
00103   typedef Teuchos::ScalarTraits<Scalar> ST;
00104   isInitialized_ = false;
00105   haveInitialCondition_ = false;
00106   model_ = Teuchos::null;
00107   solver_ = Teuchos::null;
00108   x_old_ = Teuchos::null;
00109   scaled_x_old_ = Teuchos::null;
00110   x_dot_old_ = Teuchos::null;
00111   // basePoint_;
00112   x_ = Teuchos::null;
00113   x_dot_ = Teuchos::null;
00114   dx_ = Teuchos::null;
00115   t_ = ST::nan();
00116   t_old_ = ST::nan();
00117   dt_ = ST::nan();
00118   numSteps_ = -1;
00119   neModel_ = Teuchos::null;
00120   parameterList_ = Teuchos::null;
00121   interpolator_ = Teuchos::null;
00122   stepControl_ = Teuchos::null;
00123   newtonConvergenceStatus_ = -1;
00124 }
00125 
00126 // Overridden from InterpolatorAcceptingObjectBase
00127 
00128 template<class Scalar>
00129 void BackwardEulerStepper<Scalar>::setInterpolator(
00130   const RCP<InterpolatorBase<Scalar> >& interpolator
00131   )
00132 {
00133 #ifdef HAVE_RYTHMOS_DEBUG
00134   TEUCHOS_TEST_FOR_EXCEPT(is_null(interpolator));
00135 #endif
00136   interpolator_ = interpolator;
00137   isInitialized_ = false;
00138 }
00139 
00140 template<class Scalar>
00141 RCP<InterpolatorBase<Scalar> >
00142 BackwardEulerStepper<Scalar>::getNonconstInterpolator()
00143 {
00144   return interpolator_;
00145 }
00146 
00147 template<class Scalar>
00148 RCP<const InterpolatorBase<Scalar> >
00149 BackwardEulerStepper<Scalar>::getInterpolator() const
00150 {
00151   return interpolator_;
00152 }
00153 
00154 template<class Scalar>
00155 RCP<InterpolatorBase<Scalar> >
00156 BackwardEulerStepper<Scalar>::unSetInterpolator()
00157 {
00158   RCP<InterpolatorBase<Scalar> > temp_interpolator = interpolator_;
00159   interpolator_ = Teuchos::null;
00160   return(temp_interpolator);
00161   isInitialized_ = false;
00162 }
00163 
00164 
00165 // Overridden from SolverAcceptingStepperBase
00166 
00167 
00168 template<class Scalar>
00169 void BackwardEulerStepper<Scalar>::setSolver(
00170   const RCP<Thyra::NonlinearSolverBase<Scalar> > &solver
00171   )
00172 {
00173   using Teuchos::as;
00174 
00175   TEUCHOS_TEST_FOR_EXCEPTION(solver == Teuchos::null, std::logic_error,
00176       "Error!  Thyra::NonlinearSolverBase RCP passed in through "
00177       "BackwardEulerStepper::setSolver is null!");
00178 
00179   RCP<Teuchos::FancyOStream> out = this->getOStream();
00180   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00181   Teuchos::OSTab ostab(out,1,"BES::setSolver");
00182   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00183     *out << "solver = " << solver->description() << std::endl;
00184   }
00185 
00186   solver_ = solver;
00187 
00188   isInitialized_ = false;
00189 
00190 }
00191 
00192 
00193 template<class Scalar>
00194 RCP<Thyra::NonlinearSolverBase<Scalar> >
00195 BackwardEulerStepper<Scalar>::getNonconstSolver()
00196 {
00197   return solver_;
00198 }
00199 
00200 
00201 template<class Scalar>
00202 RCP<const Thyra::NonlinearSolverBase<Scalar> >
00203 BackwardEulerStepper<Scalar>::getSolver() const
00204 {
00205   return solver_;
00206 }
00207 
00208 
00209 // Overridden from StepperBase
00210 
00211 
00212 template<class Scalar>
00213 bool BackwardEulerStepper<Scalar>::supportsCloning() const
00214 {
00215   return true;
00216 }
00217 
00218 
00219 template<class Scalar>
00220 RCP<StepperBase<Scalar> >
00221 BackwardEulerStepper<Scalar>::cloneStepperAlgorithm() const
00222 {
00223   RCP<BackwardEulerStepper<Scalar> >
00224     stepper = Teuchos::rcp(new BackwardEulerStepper<Scalar>);
00225   stepper->isInitialized_ = isInitialized_;
00226   stepper->model_ = model_; // Model is stateless so shallow copy is okay!
00227   if (!is_null(solver_))
00228     stepper->solver_ = solver_->cloneNonlinearSolver().assert_not_null();
00229   if (!is_null(x_))
00230     stepper->x_ = x_->clone_v().assert_not_null();
00231   if (!is_null(scaled_x_old_))
00232     stepper->scaled_x_old_ = scaled_x_old_->clone_v().assert_not_null();
00233   stepper->t_ = t_;
00234   stepper->t_old_ = t_old_;
00235   stepper->dt_ = dt_;
00236   stepper->numSteps_ = numSteps_;
00237   if (!is_null(neModel_))
00238     stepper->neModel_
00239     = Teuchos::rcp(new Rythmos::SingleResidualModelEvaluator<Scalar>);
00240   if (!is_null(parameterList_))
00241     stepper->parameterList_ = parameterList(*parameterList_);
00242   if (!is_null(interpolator_))
00243     stepper->interpolator_
00244       = interpolator_->cloneInterpolator().assert_not_null(); // ToDo: Implement cloneInterpolator()
00245   if (!is_null(stepControl_)) {
00246     if (stepControl_->supportsCloning())
00247       stepper->setStepControlStrategy(
00248         stepControl_->cloneStepControlStrategyAlgorithm().assert_not_null());
00249   }
00250   return stepper;
00251 }
00252 
00253 template<class Scalar>
00254 bool BackwardEulerStepper<Scalar>::isImplicit() const
00255 {
00256   return true;
00257 }
00258 
00259 template<class Scalar>
00260 void BackwardEulerStepper<Scalar>::setModel(
00261   const RCP<const Thyra::ModelEvaluator<Scalar> >& model)
00262 {
00263   using Teuchos::as;
00264 
00265   TEUCHOS_TEST_FOR_EXCEPT( is_null(model) );
00266   assertValidModel( *this, *model );
00267 
00268   RCP<Teuchos::FancyOStream> out = this->getOStream();
00269   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00270   Teuchos::OSTab ostab(out,1,"BES::setModel");
00271   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00272     *out << "model = " << model->description() << std::endl;
00273   }
00274   model_ = model;
00275 }
00276 
00277 template<class Scalar>
00278 void BackwardEulerStepper<Scalar>::setNonconstModel(
00279   const RCP<Thyra::ModelEvaluator<Scalar> >& model
00280   )
00281 {
00282   this->setModel(model); // TODO:  09/09/09 tscoffe:  Use ConstNonconstObjectContainer here!
00283 }
00284 
00285 template<class Scalar>
00286 RCP<const Thyra::ModelEvaluator<Scalar> >
00287 BackwardEulerStepper<Scalar>::getModel() const
00288 {
00289   return model_;
00290 }
00291 
00292 
00293 template<class Scalar>
00294 RCP<Thyra::ModelEvaluator<Scalar> >
00295 BackwardEulerStepper<Scalar>::getNonconstModel()
00296 {
00297   return Teuchos::null;
00298 }
00299 
00300 
00301 template<class Scalar>
00302 void BackwardEulerStepper<Scalar>::setInitialCondition(
00303   const Thyra::ModelEvaluatorBase::InArgs<Scalar> &initialCondition)
00304 {
00305 
00306   typedef Teuchos::ScalarTraits<Scalar> ST;
00307   typedef Thyra::ModelEvaluatorBase MEB;
00308 
00309   basePoint_ = initialCondition;
00310 
00311   // x
00312   RCP<const Thyra::VectorBase<Scalar> > x_init = initialCondition.get_x();
00313   TEUCHOS_TEST_FOR_EXCEPTION(
00314     is_null(x_init), std::logic_error,
00315     "Error, if the client passes in an intial condition to "
00316     "setInitialCondition(...), then x can not be null!" );
00317   x_ = x_init->clone_v();
00318 
00319   // x_dot
00320   RCP<const Thyra::VectorBase<Scalar> >
00321     x_dot_init = initialCondition.get_x_dot();
00322   if (!is_null(x_dot_init)) {
00323     x_dot_ = x_dot_init->clone_v();
00324   }
00325   else {
00326     x_dot_ = createMember(x_->space());
00327     V_S(x_dot_.ptr(),ST::zero());
00328   }
00329 
00330   // dx
00331   if (is_null(dx_)) dx_ = createMember(x_->space());
00332   V_S(dx_.ptr(), ST::zero());
00333 
00334   // t
00335   t_ = initialCondition.get_t();
00336   t_old_ = t_;
00337 
00338   // x_old
00339   if (is_null(x_old_)) x_old_ = createMember(x_->space());
00340   x_old_ = x_init->clone_v();
00341   scaled_x_old_ = x_->clone_v();
00342 
00343   // x_dot_old
00344   if (is_null(x_dot_old_)) x_dot_old_ = createMember(x_->space());
00345   x_dot_old_ = x_dot_->clone_v();
00346 
00347   haveInitialCondition_ = true;
00348 }
00349 
00350 
00351 template<class Scalar>
00352 Thyra::ModelEvaluatorBase::InArgs<Scalar>
00353 BackwardEulerStepper<Scalar>::getInitialCondition() const
00354 {
00355   return basePoint_;
00356 }
00357 
00358 template<class Scalar>
00359 void BackwardEulerStepper<Scalar>::setStepControlStrategy(
00360   const RCP<StepControlStrategyBase<Scalar> >& stepControl)
00361 {
00362   TEUCHOS_TEST_FOR_EXCEPTION(stepControl == Teuchos::null,std::logic_error,
00363     "Error, stepControl == Teuchos::null!\n");
00364   stepControl_ = stepControl;
00365 }
00366 
00367 template<class Scalar>
00368 RCP<StepControlStrategyBase<Scalar> >
00369 BackwardEulerStepper<Scalar>::getNonconstStepControlStrategy()
00370 {
00371   return(stepControl_);
00372 }
00373 
00374 template<class Scalar>
00375 RCP<const StepControlStrategyBase<Scalar> >
00376 BackwardEulerStepper<Scalar>::getStepControlStrategy() const
00377 {
00378   return(stepControl_);
00379 }
00380 
00381 template<class Scalar>
00382 Scalar BackwardEulerStepper<Scalar>::takeStep(Scalar dt,
00383                                               StepSizeType stepSizeType)
00384 {
00385 
00386   using Teuchos::as;
00387   using Teuchos::incrVerbLevel;
00388   typedef Teuchos::ScalarTraits<Scalar> ST;
00389   typedef Thyra::NonlinearSolverBase<Scalar> NSB;
00390   typedef Teuchos::VerboseObjectTempState<NSB> VOTSNSB;
00391 
00392   RCP<Teuchos::FancyOStream> out = this->getOStream();
00393   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00394   Teuchos::OSTab ostab(out,1,"BES::takeStep");
00395   VOTSNSB solver_outputTempState(solver_,out,incrVerbLevel(verbLevel,-1));
00396 
00397   if ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
00398     *out << "\nEntering "
00399          << Teuchos::TypeNameTraits<BackwardEulerStepper<Scalar> >::name()
00400          << "::takeStep("<<dt<<","<<toString(stepSizeType)<<") ...\n";
00401   }
00402 
00403   initialize_();
00404 
00405   if(!neModel_.get()) {
00406     neModel_ =Teuchos::rcp(new Rythmos::SingleResidualModelEvaluator<Scalar>());
00407   }
00408 
00409   dt_ = dt;
00410 
00411   if (dt <= ST::zero()) {
00412     if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) )
00413       *out << "\nThe arguments to takeStep are not valid for "
00414            << "BackwardEulerStepper at this time.\n"
00415            << "  dt = " << dt << "\n"
00416            << "BackwardEulerStepper requires positive dt.\n" << std::endl;
00417     return(Scalar(-ST::one()));
00418   }
00419   if ((stepSizeType == STEP_TYPE_VARIABLE) && (stepControl_ == Teuchos::null)) {
00420     if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) )
00421       *out << "\nFor 'variable' time stepping (step size adjustable by "
00422            << "BackwardEulerStepper), BackwardEulerStepper requires "
00423            << "Step-Control Strategy.\n"
00424            << "  stepType = " << toString(stepSizeType) << "\n" << std::endl;
00425     return(Scalar(-ST::one()));
00426   }
00427 
00428   stepControl_->setRequestedStepSize(*this,dt_,stepSizeType);
00429   AttemptedStepStatusFlag status;
00430   bool stepPass = false;
00431   while (1) {
00432 
00433     stepControl_->nextStepSize(*this,&dt_,&stepSizeType,NULL);
00434     if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00435       *out << "\nrequested dt = " << dt
00436            << "\ncurrent dt   = " << dt_ << "\n";
00437     }
00438 
00439     // Setup the nonlinear equations:
00440     //
00441     //   f( (1/dt)* x + (-1/dt)*x_old), x, t ) = 0
00442 
00443     V_StV( scaled_x_old_.ptr(), Scalar(-ST::one()/dt_), *x_old_ );
00444     t_old_ = t_;
00445     neModel_->initializeSingleResidualModel(
00446       model_, basePoint_,
00447       Scalar(ST::one()/dt_), scaled_x_old_,
00448       ST::one(), Teuchos::null,
00449       t_old_+dt_,
00450       Teuchos::null
00451       );
00452     if( solver_->getModel().get() != neModel_.get() ) {
00453       solver_->setModel(neModel_);
00454     }
00455     // 2007/05/18: rabartl: ToDo: Above, set the stream and the verbosity level
00456     // on solver_ so that we an see what it is doing!
00457 
00458     obtainPredictor_();
00459 
00460     // Solve the implicit nonlinear system to a tolerance of ???
00461 
00462     if ( as<int>(verbLevel) > as<int>(Teuchos::VERB_LOW) ) {
00463       *out << "\nSolving the implicit backward-Euler timestep equation ...\n";
00464     }
00465 
00466     Thyra::SolveStatus<Scalar> neSolveStatus =
00467       solver_->solve(&*x_, NULL, &*dx_);
00468 
00469     // In the above solve, on input *x_ is the initial guess that comes from
00470     // the predictor.  On output, *x_ is the converged timestep solution and
00471     // *dx_ is the difference computed from the intial guess in *x_ to the
00472     // final solved value of *x_.  This is needed for basic numerical stability.
00473 
00474     if ( as<int>(verbLevel) > as<int>(Teuchos::VERB_LOW) ) {
00475       *out << "\nOutput status of nonlinear solve:\n" << neSolveStatus;
00476     }
00477 
00478     // 2007/05/18: rabartl: ToDo: Above, get the solve status from the above
00479     // solve and at least print warning message if the solve fails!  Actually,
00480     // you should most likely thrown an exception if the solve fails or return
00481     // false if appropriate
00482 
00483     if (neSolveStatus.solveStatus == Thyra::SOLVE_STATUS_CONVERGED)  {
00484       newtonConvergenceStatus_ = 0;
00485     }
00486     else {
00487       newtonConvergenceStatus_ = -1;
00488     }
00489 
00490     stepControl_->setCorrection(*this,x_,dx_,newtonConvergenceStatus_);
00491 
00492     stepPass = stepControl_->acceptStep(*this,NULL);
00493 
00494     if (!stepPass) { // stepPass = false
00495       status = stepControl_->rejectStep(*this);
00496 
00497       if (status != PREDICT_AGAIN)
00498         break;
00499     } else { // stepPass = true
00500       break;
00501     }
00502   }
00503 
00504   // Update the step
00505 
00506   if (stepPass) {
00507     V_V( x_dot_old_.ptr(), *x_dot_ );
00508     // x_dot = (1/dt)*x - (1/dt)*x_old
00509     V_StVpStV( x_dot_.ptr(), Scalar(ST::one()/dt_), *x_,
00510                              Scalar(-ST::one()/dt_), *x_old_ );
00511     V_V( x_old_.ptr(), *x_ );
00512     t_ += dt_;
00513     numSteps_++;
00514     stepControl_->completeStep(*this);
00515   } else {
00516     // Complete failure.  Return to Integrator with bad step size.
00517     dt_ = Scalar(-ST::one());
00518   }
00519 
00520   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00521     *out << "\nt_old_ = " << t_old_ << std::endl;
00522     *out << "\nt_ = " << t_ << std::endl;
00523   }
00524 
00525 #ifdef HAVE_RYTHMOS_DEBUG
00526   // 04/14/09 tscoffe: This code should be moved to StepperValidator
00527 
00528   if ( includesVerbLevel(verbLevel,Teuchos::VERB_LOW) )
00529     *out << "\nChecking to make sure that solution and "
00530          << "the interpolated solution are the same! ...\n";
00531 
00532   {
00533 
00534     typedef ScalarTraits<Scalar> ST;
00535     typedef ScalarTraits<ScalarMag> SMT;
00536 
00537     Teuchos::OSTab tab(out);
00538 
00539     const StepStatus<Scalar> stepStatus = this->getStepStatus();
00540 
00541     RCP<const Thyra::VectorBase<Scalar> >
00542       x = stepStatus.solution,
00543       xdot = stepStatus.solutionDot;
00544 
00545     Array<Scalar> time_vec = Teuchos::tuple(stepStatus.time);
00546     Array<RCP<const Thyra::VectorBase<Scalar> > > x_vec, xdot_vec;
00547     this->getPoints(time_vec,&x_vec,&xdot_vec,0);
00548 
00549     RCP<const Thyra::VectorBase<Scalar> >
00550       x_interp = x_vec[0],
00551       xdot_interp = xdot_vec[0];
00552 
00553     TEUCHOS_TEST_FOR_EXCEPT(
00554       !Thyra::testRelNormDiffErr(
00555         "x", *x, "x_interp", *x_interp,
00556         "2*epsilon", ScalarMag(100.0*SMT::eps()),
00557         "2*epsilon", ScalarMag(100.0*SMT::eps()),
00558         includesVerbLevel(verbLevel,Teuchos::VERB_HIGH) ? out.get() : 0
00559         )
00560       );
00561 
00562     TEUCHOS_TEST_FOR_EXCEPT(
00563       !Thyra::testRelNormDiffErr(
00564         "xdot", *xdot, "xdot_interp", *xdot_interp,
00565         "2*epsilon", ScalarMag(100.0*SMT::eps()),
00566         "2*epsilon", ScalarMag(100.0*SMT::eps()),
00567         includesVerbLevel(verbLevel,Teuchos::VERB_HIGH) ? out.get() : 0
00568         )
00569       );
00570 
00571   }
00572 
00573   // 2007/07/25: rabartl: ToDo: Move the above test into a helper function so
00574   // that it can be used from lots of different places!
00575 
00576 #endif // HAVE_RYTHMOS_DEBUG
00577 
00578   if ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
00579     *out << "\nLeaving "
00580          << Teuchos::TypeNameTraits<BackwardEulerStepper<Scalar> >::name()
00581          << "::takeStep("<<dt_<<","<<toString(stepSizeType)<<") ...\n";
00582   }
00583 
00584   return(dt_);
00585 }
00586 
00587 
00588 template<class Scalar>
00589 const StepStatus<Scalar> BackwardEulerStepper<Scalar>::getStepStatus() const
00590 {
00591 
00592   typedef Teuchos::ScalarTraits<Scalar> ST;
00593 
00594   StepStatus<Scalar> stepStatus; // Defaults to unknown status
00595 
00596   if (!isInitialized_) {
00597     stepStatus.stepStatus = STEP_STATUS_UNINITIALIZED;
00598   }
00599   else if (numSteps_ > 0) {
00600     stepStatus.stepStatus = STEP_STATUS_CONVERGED;
00601   }
00602   // else unknown
00603 
00604   stepStatus.stepSize = dt_;
00605   stepStatus.order = 1;
00606   stepStatus.time = t_;
00607   stepStatus.solution = x_;
00608   stepStatus.solutionDot = x_dot_;
00609 
00610   return(stepStatus);
00611 
00612 }
00613 
00614 
00615 // Overridden from InterpolationBufferBase
00616 
00617 
00618 template<class Scalar>
00619 RCP<const Thyra::VectorSpaceBase<Scalar> >
00620 BackwardEulerStepper<Scalar>::get_x_space() const
00621 {
00622   return ( !is_null(model_) ? model_->get_x_space() : Teuchos::null );
00623 }
00624 
00625 
00626 template<class Scalar>
00627 void BackwardEulerStepper<Scalar>::addPoints(
00628   const Array<Scalar>& time_vec,
00629   const Array<RCP<const Thyra::VectorBase<Scalar> > >& x_vec,
00630   const Array<RCP<const Thyra::VectorBase<Scalar> > >& xdot_vec
00631   )
00632 {
00633 
00634   typedef Teuchos::ScalarTraits<Scalar> ST;
00635   using Teuchos::as;
00636 
00637   initialize_();
00638 
00639 #ifdef HAVE_RYTHMOS_DEBUG
00640   TEUCHOS_TEST_FOR_EXCEPTION(
00641     time_vec.size() == 0, std::logic_error,
00642     "Error, addPoints called with an empty time_vec array!\n");
00643 #endif // HAVE_RYTHMOS_DEBUG
00644 
00645   RCP<Teuchos::FancyOStream> out = this->getOStream();
00646   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00647   Teuchos::OSTab ostab(out,1,"BES::setPoints");
00648 
00649   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00650     *out << "time_vec = " << std::endl;
00651     for (int i=0 ; i<Teuchos::as<int>(time_vec.size()) ; ++i) {
00652       *out << "time_vec[" << i << "] = " << time_vec[i] << std::endl;
00653     }
00654   }
00655   else if (time_vec.size() == 1) {
00656     int n = 0;
00657     t_ = time_vec[n];
00658     t_old_ = t_;
00659     Thyra::V_V(x_.ptr(),*x_vec[n]);
00660     Thyra::V_V(scaled_x_old_.ptr(),*x_);
00661   }
00662   else {
00663     int n = time_vec.size()-1;
00664     int nm1 = time_vec.size()-2;
00665     t_ = time_vec[n];
00666     t_old_ = time_vec[nm1];
00667     Thyra::V_V(x_.ptr(),*x_vec[n]);
00668     Scalar dt = t_ - t_old_;
00669     Thyra::V_StV(scaled_x_old_.ptr(),Scalar(-ST::one()/dt),*x_vec[nm1]);
00670   }
00671 }
00672 
00673 
00674 template<class Scalar>
00675 TimeRange<Scalar> BackwardEulerStepper<Scalar>::getTimeRange() const
00676 {
00677   if ( !isInitialized_ && haveInitialCondition_ )
00678     return timeRange<Scalar>(t_,t_);
00679   if ( !isInitialized_ && !haveInitialCondition_ )
00680     return invalidTimeRange<Scalar>();
00681   return timeRange<Scalar>(t_old_,t_);
00682 }
00683 
00684 
00685 template<class Scalar>
00686 void BackwardEulerStepper<Scalar>::getPoints(
00687   const Array<Scalar>& time_vec,
00688   Array<RCP<const Thyra::VectorBase<Scalar> > >* x_vec,
00689   Array<RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec,
00690   Array<ScalarMag>* accuracy_vec
00691   ) const
00692 {
00693   typedef Teuchos::ScalarTraits<Scalar> ST;
00694   using Teuchos::constOptInArg;
00695   using Teuchos::ptr;
00696 #ifdef HAVE_RYTHMOS_DEBUG
00697   TEUCHOS_ASSERT(haveInitialCondition_);
00698 #endif // HAVE_RYTHMOS_DEBUG
00699   RCP<Thyra::VectorBase<Scalar> > x_temp = x_;
00700   if (compareTimeValues(t_old_,t_)!=0) {
00701     Scalar dt = t_ - t_old_;
00702     x_temp = scaled_x_old_->clone_v();
00703     Thyra::Vt_S(x_temp.ptr(),Scalar(-ST::one()*dt));  // undo the scaling
00704   }
00705   defaultGetPoints<Scalar>(
00706       t_old_, constOptInArg(*x_temp), constOptInArg(*x_dot_old_),
00707       t_, constOptInArg(*x_), constOptInArg(*x_dot_),
00708       time_vec, ptr(x_vec), ptr(xdot_vec), ptr(accuracy_vec),
00709       ptr(interpolator_.get())
00710       );
00711 
00712   /*
00713   using Teuchos::as;
00714   typedef Teuchos::ScalarTraits<Scalar> ST;
00715   typedef typename ST::magnitudeType ScalarMag;
00716   typedef Teuchos::ScalarTraits<ScalarMag> SMT;
00717   typename DataStore<Scalar>::DataStoreVector_t ds_nodes;
00718   typename DataStore<Scalar>::DataStoreVector_t ds_out;
00719 
00720 #ifdef HAVE_RYTHMOS_DEBUG
00721   TEUCHOS_TEST_FOR_EXCEPT(!haveInitialCondition_);
00722   TEUCHOS_TEST_FOR_EXCEPT( 0 == x_vec );
00723 #endif
00724 
00725   RCP<Teuchos::FancyOStream> out = this->getOStream();
00726   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00727   Teuchos::OSTab ostab(out,1,"BES::getPoints");
00728   if ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
00729     *out
00730       << "\nEntering " << Teuchos::TypeNameTraits<BackwardEulerStepper<Scalar> >::name()
00731       << "::getPoints(...) ...\n";
00732   }
00733   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00734     for (int i=0 ; i<Teuchos::as<int>(time_vec.size()) ; ++i) {
00735       *out << "time_vec[" << i << "] = " << time_vec[i] << std::endl;
00736     }
00737     *out << "I can interpolate in the interval [" << t_old_ << "," << t_ << "]." << std::endl;
00738   }
00739   if (t_old_ != t_) {
00740     if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00741       *out << "Passing two points to interpolator:  " << t_old_ << " and " << t_ << std::endl;
00742     }
00743     DataStore<Scalar> ds_temp;
00744     Scalar dt = t_ - t_old_;
00745 #ifdef HAVE_RYTHMOS_DEBUG
00746     TEUCHOS_TEST_FOR_EXCEPT(
00747       !Thyra::testRelErr(
00748         "dt", dt, "dt_", dt_,
00749         "1e+4*epsilon", ScalarMag(1e+4*SMT::eps()),
00750         "1e+2*epsilon", ScalarMag(1e+2*SMT::eps()),
00751         as<int>(verbLevel) >= as<int>(Teuchos::VERB_MEDIUM) ? out.get() : 0
00752         )
00753       );
00754 #endif
00755     RCP<Thyra::VectorBase<Scalar> >
00756       x_temp = scaled_x_old_->clone_v();
00757     Thyra::Vt_S(&*x_temp,Scalar(-ST::one()*dt));  // undo the scaling
00758     ds_temp.time = t_old_;
00759     ds_temp.x = x_temp;
00760     ds_temp.xdot = x_dot_old_;
00761     ds_temp.accuracy = ScalarMag(dt);
00762     ds_nodes.push_back(ds_temp);
00763     ds_temp.time = t_;
00764     ds_temp.x = x_;
00765     ds_temp.xdot = x_dot_;
00766     ds_temp.accuracy = ScalarMag(dt);
00767     ds_nodes.push_back(ds_temp);
00768   }
00769   else {
00770     if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00771       *out << "Passing one point to interpolator:  " << t_ << std::endl;
00772     }
00773     DataStore<Scalar> ds_temp;
00774     ds_temp.time = t_;
00775     ds_temp.x = x_;
00776     ds_temp.xdot = x_dot_;
00777     ds_temp.accuracy = ScalarMag(ST::zero());
00778     ds_nodes.push_back(ds_temp);
00779   }
00780   interpolate<Scalar>(*interpolator_,rcp(&ds_nodes,false),time_vec,&ds_out);
00781   Array<Scalar> time_out;
00782   dataStoreVectorToVector(ds_out,&time_out,x_vec,xdot_vec,accuracy_vec);
00783   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00784     *out << "Passing out the interpolated values:" << std::endl;
00785     for (int i=0; i<Teuchos::as<int>(time_out.size()) ; ++i) {
00786       *out << "time[" << i << "] = " << time_out[i] << std::endl;
00787       *out << "x_vec[" << i << "] = " << std::endl;
00788       (*x_vec)[i]->describe(*out,Teuchos::VERB_EXTREME);
00789       if (xdot_vec) {
00790         if ( (*xdot_vec)[i] == Teuchos::null) {
00791           *out << "xdot_vec[" << i << "] = Teuchos::null" << std::endl;
00792         }
00793         else {
00794           *out << "xdot_vec[" << i << "] = " << std::endl;
00795           (*xdot_vec)[i]->describe(*out,Teuchos::VERB_EXTREME);
00796         }
00797       }
00798       if(accuracy_vec) {
00799         *out << "accuracy[" << i << "] = " << (*accuracy_vec)[i] << std::endl;
00800       }
00801     }
00802   }
00803   if ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
00804     *out
00805       << "Leaving " << Teuchos::TypeNameTraits<BackwardEulerStepper<Scalar> >::name()
00806       << "::getPoints(...) ...\n";
00807   }
00808   */
00809 
00810 }
00811 
00812 
00813 template<class Scalar>
00814 void BackwardEulerStepper<Scalar>::getNodes(Array<Scalar>* time_vec) const
00815 {
00816   using Teuchos::as;
00817 
00818   TEUCHOS_ASSERT( time_vec != NULL );
00819 
00820   time_vec->clear();
00821 
00822   if (!haveInitialCondition_) {
00823     return;
00824   }
00825 
00826   time_vec->push_back(t_old_);
00827 
00828   if (numSteps_ > 0) {
00829     time_vec->push_back(t_);
00830   }
00831 
00832   RCP<Teuchos::FancyOStream> out = this->getOStream();
00833   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00834   Teuchos::OSTab ostab(out,1,"BES::getNodes");
00835   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00836     *out << this->description() << std::endl;
00837     for (int i=0 ; i<Teuchos::as<int>(time_vec->size()) ; ++i) {
00838       *out << "time_vec[" << i << "] = " << (*time_vec)[i] << std::endl;
00839     }
00840   }
00841 }
00842 
00843 
00844 template<class Scalar>
00845 void BackwardEulerStepper<Scalar>::removeNodes(Array<Scalar>& time_vec)
00846 {
00847   initialize_();
00848   using Teuchos::as;
00849   RCP<Teuchos::FancyOStream> out = this->getOStream();
00850   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00851   Teuchos::OSTab ostab(out,1,"BES::removeNodes");
00852   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00853     *out << "time_vec = " << std::endl;
00854     for (int i=0 ; i<Teuchos::as<int>(time_vec.size()) ; ++i) {
00855       *out << "time_vec[" << i << "] = " << time_vec[i] << std::endl;
00856     }
00857   }
00858   TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Error, removeNodes is not implemented for BackwardEulerStepper at this time.\n");
00859   // TODO:
00860   // if any time in time_vec matches t_ or t_old_, then do the following:
00861   // remove t_old_:  set t_old_ = t_ and set scaled_x_old_ = x_
00862   // remove t_:  set t_ = t_old_ and set x_ = -dt*scaled_x_old_
00863 }
00864 
00865 
00866 template<class Scalar>
00867 int BackwardEulerStepper<Scalar>::getOrder() const
00868 {
00869   return(1);
00870 }
00871 
00872 
00873 // Overridden from Teuchos::ParameterListAcceptor
00874 
00875 
00876 template <class Scalar>
00877 void BackwardEulerStepper<Scalar>::setParameterList(
00878   RCP<Teuchos::ParameterList> const& paramList
00879   )
00880 {
00881   TEUCHOS_TEST_FOR_EXCEPT(is_null(paramList));
00882   paramList->validateParametersAndSetDefaults(*this->getValidParameters());
00883   parameterList_ = paramList;
00884   Teuchos::readVerboseObjectSublist(&*parameterList_,this);
00885 }
00886 
00887 
00888 template <class Scalar>
00889 RCP<Teuchos::ParameterList>
00890 BackwardEulerStepper<Scalar>::getNonconstParameterList()
00891 {
00892   return(parameterList_);
00893 }
00894 
00895 
00896 template <class Scalar>
00897 RCP<Teuchos::ParameterList>
00898 BackwardEulerStepper<Scalar>::unsetParameterList()
00899 {
00900   RCP<Teuchos::ParameterList>
00901     temp_param_list = parameterList_;
00902   parameterList_ = Teuchos::null;
00903   return(temp_param_list);
00904 }
00905 
00906 
00907 template<class Scalar>
00908 RCP<const Teuchos::ParameterList>
00909 BackwardEulerStepper<Scalar>::getValidParameters() const
00910 {
00911   using Teuchos::ParameterList;
00912   static RCP<const ParameterList> validPL;
00913   if (is_null(validPL)) {
00914     RCP<ParameterList> pl = Teuchos::parameterList();
00915     // This line is required to pass StepperValidator UnitTest!
00916     pl->sublist(RythmosStepControlSettings_name);
00917     Teuchos::setupVerboseObjectSublist(&*pl);
00918     validPL = pl;
00919   }
00920   return validPL;
00921 }
00922 
00923 
00924 // Overridden from Teuchos::Describable
00925 
00926 
00927 template<class Scalar>
00928 void BackwardEulerStepper<Scalar>::describe(
00929   Teuchos::FancyOStream &out,
00930   const Teuchos::EVerbosityLevel verbLevel
00931   ) const
00932 {
00933   using Teuchos::as;
00934   Teuchos::OSTab tab(out);
00935   if (!isInitialized_) {
00936     out << this->description() << " : This stepper is not initialized yet"
00937         << std::endl;
00938     return;
00939   }
00940   if (
00941     as<int>(verbLevel) == as<int>(Teuchos::VERB_DEFAULT)
00942     ||
00943     as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW)
00944     )
00945   {
00946     out << this->description() << "::describe:" << std::endl;
00947     out << "model = " << model_->description() << std::endl;
00948     out << "solver = " << solver_->description() << std::endl;
00949     if (neModel_ == Teuchos::null) {
00950       out << "neModel = Teuchos::null" << std::endl;
00951     } else {
00952       out << "neModel = " << neModel_->description() << std::endl;
00953     }
00954   }
00955   else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW)) {
00956     out << "t_ = " << t_ << std::endl;
00957     out << "t_old_ = " << t_old_ << std::endl;
00958   }
00959   else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_MEDIUM)) {
00960   }
00961   else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH)) {
00962     out << "model_ = " << std::endl;
00963     model_->describe(out,verbLevel);
00964     out << "solver_ = " << std::endl;
00965     solver_->describe(out,verbLevel);
00966     if (neModel_ == Teuchos::null) {
00967       out << "neModel = Teuchos::null" << std::endl;
00968     } else {
00969       out << "neModel = " << std::endl;
00970       neModel_->describe(out,verbLevel);
00971     }
00972     out << "x_ = " << std::endl;
00973     x_->describe(out,verbLevel);
00974     out << "scaled_x_old_ = " << std::endl;
00975     scaled_x_old_->describe(out,verbLevel);
00976   }
00977 }
00978 
00979 
00980 // private
00981 
00982 
00983 template <class Scalar>
00984 void BackwardEulerStepper<Scalar>::initialize_()
00985 {
00986 
00987   if (isInitialized_) return;
00988 
00989   TEUCHOS_TEST_FOR_EXCEPT(is_null(model_));
00990   TEUCHOS_TEST_FOR_EXCEPT(is_null(solver_));
00991   TEUCHOS_TEST_FOR_EXCEPT(!haveInitialCondition_);
00992 
00993   // Initialize Parameter List if none provided.
00994   if (parameterList_ == Teuchos::null) {
00995     RCP<Teuchos::ParameterList> emptyParameterList =
00996       Teuchos::rcp(new Teuchos::ParameterList);
00997     this->setParameterList(emptyParameterList);
00998   }
00999 
01000   // Initialize StepControl
01001   if (stepControl_ == Teuchos::null) {
01002     RCP<StepControlStrategyBase<Scalar> > stepControlStrategy =
01003       Teuchos::rcp(new FixedStepControlStrategy<Scalar>());
01004     RCP<Teuchos::ParameterList> stepControlPL =
01005       Teuchos::sublist(parameterList_, RythmosStepControlSettings_name);
01006 
01007     stepControlStrategy->setParameterList(stepControlPL);
01008     this->setStepControlStrategy(stepControlStrategy);
01009   }
01010   stepControl_->initialize(*this);
01011   stepControl_->setOStream(this->getOStream());
01012   stepControl_->setVerbLevel(this->getVerbLevel());
01013 
01014   //maxOrder_ = stepControl_->getMaxOrder(); // maximum order
01015 
01016 #ifdef HAVE_RYTHMOS_DEBUG
01017   THYRA_ASSERT_VEC_SPACES(
01018     "Rythmos::BackwardEulerStepper::initialize_(...)",
01019     *x_->space(), *model_->get_x_space() );
01020 #endif // HAVE_RYTHMOS_DEBUG
01021 
01022   if ( is_null(interpolator_) ) {
01023     // If an interpolator has not been explicitly set, then just create
01024     // a default linear interpolator.
01025     interpolator_ = linearInterpolator<Scalar>();
01026     // 2007/05/18: rabartl: ToDo: Replace this with a Hermite interplator
01027     // when it is implementated!
01028   }
01029   isInitialized_ = true;
01030 }
01031 
01032 
01033 template<class Scalar>
01034 RCP<const MomentoBase<Scalar> >
01035 BackwardEulerStepper<Scalar>::getMomento() const
01036 {
01037   RCP<BackwardEulerStepperMomento<Scalar> > momento =
01038     Teuchos::rcp(new BackwardEulerStepperMomento<Scalar>());
01039   momento->set_scaled_x_old(scaled_x_old_);
01040   momento->set_x_old(x_old_);
01041   momento->set_x_dot_old(x_dot_old_);
01042   momento->set_x(x_);
01043   momento->set_x_dot(x_dot_);
01044   momento->set_dx(dx_);
01045   momento->set_t(t_);
01046   momento->set_t_old(t_old_);
01047   momento->set_dt(dt_);
01048   momento->set_numSteps(numSteps_);
01049   momento->set_newtonConvergenceStatus(newtonConvergenceStatus_);
01050   momento->set_isInitialized(isInitialized_);
01051   momento->set_haveInitialCondition(haveInitialCondition_);
01052   momento->set_parameterList(parameterList_);
01053   momento->set_basePoint(basePoint_);
01054   momento->set_neModel(neModel_);
01055   momento->set_interpolator(interpolator_);
01056   momento->set_stepControl(stepControl_);
01057   return momento;
01058 }
01059 
01060 template<class Scalar>
01061 void BackwardEulerStepper<Scalar>::setMomento(
01062     const Ptr<const MomentoBase<Scalar> >& momentoPtr,
01063     const RCP<Thyra::ModelEvaluator<Scalar> >& model,
01064     const RCP<Thyra::NonlinearSolverBase<Scalar> >& solver
01065     )
01066 {
01067   Ptr<const BackwardEulerStepperMomento<Scalar> > beMomentoPtr =
01068     Teuchos::ptr_dynamic_cast<const BackwardEulerStepperMomento<Scalar> >
01069       (momentoPtr,true);
01070   const BackwardEulerStepperMomento<Scalar>& beMomento = *beMomentoPtr;
01071   model_ = model;
01072   solver_ = solver;
01073   scaled_x_old_ = beMomento.get_scaled_x_old();
01074   x_old_ = beMomento.get_x_old();
01075   x_dot_old_ = beMomento.get_x_dot_old();
01076   x_ = beMomento.get_x();
01077   x_dot_ = beMomento.get_x_dot();
01078   dx_ = beMomento.get_dx();
01079   t_ = beMomento.get_t();
01080   t_old_ = beMomento.get_t_old();
01081   dt_ = beMomento.get_dt();
01082   numSteps_ = beMomento.get_numSteps();
01083   newtonConvergenceStatus_ = beMomento.get_newtonConvergenceStatus();
01084   isInitialized_ = beMomento.get_isInitialized();
01085   haveInitialCondition_ = beMomento.get_haveInitialCondition();
01086   parameterList_ = beMomento.get_parameterList();
01087   basePoint_ = beMomento.get_basePoint();
01088   neModel_ = beMomento.get_neModel();
01089   interpolator_ = beMomento.get_interpolator();
01090   stepControl_ = beMomento.get_stepControl();
01091   this->checkConsistentState_();
01092 }
01093 
01094 template<class Scalar>
01095 void BackwardEulerStepper<Scalar>::checkConsistentState_()
01096 {
01097   if (isInitialized_) {
01098     TEUCHOS_ASSERT( !Teuchos::is_null(model_) );
01099     TEUCHOS_ASSERT( !Teuchos::is_null(solver_) );
01100     TEUCHOS_ASSERT( haveInitialCondition_ );
01101     TEUCHOS_ASSERT( !Teuchos::is_null(interpolator_) );
01102     TEUCHOS_ASSERT( !Teuchos::is_null(stepControl_) );
01103   }
01104   if (haveInitialCondition_) {
01105     // basePoint_ should be defined
01106     typedef Teuchos::ScalarTraits<Scalar> ST;
01107     TEUCHOS_ASSERT( !ST::isnaninf(t_) );
01108     TEUCHOS_ASSERT( !ST::isnaninf(t_old_) );
01109     TEUCHOS_ASSERT( !Teuchos::is_null(scaled_x_old_) );
01110     TEUCHOS_ASSERT( !Teuchos::is_null(x_dot_old_) );
01111     TEUCHOS_ASSERT( !Teuchos::is_null(x_) );
01112     TEUCHOS_ASSERT( !Teuchos::is_null(x_dot_) );
01113     TEUCHOS_ASSERT( !Teuchos::is_null(x_old_) );
01114     TEUCHOS_ASSERT( !Teuchos::is_null(dx_) );
01115     TEUCHOS_ASSERT( t_ >= basePoint_.get_t() );
01116     TEUCHOS_ASSERT( t_old_ >= basePoint_.get_t() );
01117   }
01118   if (numSteps_ > 0) {
01119     TEUCHOS_ASSERT(isInitialized_);
01120   }
01121 }
01122 
01123 template<class Scalar>
01124 void BackwardEulerStepper<Scalar>::obtainPredictor_()
01125 {
01126   using Teuchos::as;
01127   typedef Teuchos::ScalarTraits<Scalar> ST;
01128 
01129   if (!isInitialized_) return;
01130 
01131   RCP<Teuchos::FancyOStream> out = this->getOStream();
01132   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
01133   Teuchos::OSTab ostab(out,1,"obtainPredictor_");
01134 
01135   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
01136     *out << "Before predictor x_ = " << std::endl;
01137     x_->describe(*out,verbLevel);
01138   }
01139   // evaluate predictor -- basic Forward Euler
01140   // x_ = dt_*x_dot_old_ + x_old_
01141   V_StVpV(x_.ptr(), dt_, *x_dot_old_, *x_old_);
01142 
01143   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
01144     *out << "After predictor x_ = " << std::endl;
01145     x_->describe(*out,verbLevel);
01146   }
01147 }
01148 
01149 //
01150 // Explicit Instantiation macro
01151 //
01152 // Must be expanded from within the Rythmos namespace!
01153 //
01154 
01155 #define RYTHMOS_BACKWARD_EULER_STEPPER_INSTANT(SCALAR) \
01156   \
01157   template class BackwardEulerStepper< SCALAR >; \
01158   \
01159   template RCP< BackwardEulerStepper< SCALAR > > \
01160   backwardEulerStepper( \
01161     const RCP<Thyra::ModelEvaluator< SCALAR > >& model, \
01162     const RCP<Thyra::NonlinearSolverBase< SCALAR > >& solver \
01163       ); \
01164   template RCP< BackwardEulerStepper< SCALAR > > \
01165   backwardEulerStepper();
01166 
01167 
01168 
01169 
01170 } // namespace Rythmos
01171 
01172 #endif //Rythmos_BACKWARD_EULER_STEPPER_DEF_H
 All Classes Functions Variables Typedefs Friends