Rythmos_BackwardEulerStepper.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., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Todd S. Coffey (tscoffe@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 //@HEADER
00028 
00029 #ifndef Rythmos_BACKWARD_EULER_STEPPER_H
00030 #define Rythmos_BACKWARD_EULER_STEPPER_H
00031 
00032 #include "Rythmos_StepperBase.hpp"
00033 #include "Rythmos_DataStore.hpp"
00034 #include "Rythmos_LinearInterpolator.hpp"
00035 #include "Rythmos_SingleResidualModelEvaluator.hpp"
00036 #include "Rythmos_SolverAcceptingStepperBase.hpp"
00037 
00038 #include "Thyra_VectorBase.hpp"
00039 #include "Thyra_ModelEvaluator.hpp"
00040 #include "Thyra_ModelEvaluatorHelpers.hpp"
00041 #include "Thyra_AssertOp.hpp"
00042 #include "Thyra_NonlinearSolverBase.hpp"
00043 #include "Thyra_TestingTools.hpp"
00044 
00045 #include "Teuchos_VerboseObjectParameterListHelpers.hpp"
00046 #include "Teuchos_as.hpp"
00047 
00048 
00049 namespace Rythmos {
00050 
00051 
00062 template<class Scalar>
00063 class BackwardEulerStepper : virtual public SolverAcceptingStepperBase<Scalar>
00064 {
00065 public:
00066   
00068   typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType ScalarMag;
00069   
00072 
00074   BackwardEulerStepper();
00075   
00077   BackwardEulerStepper(
00078     const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >  &model
00079     ,const Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> >  &solver
00080     );
00081   
00083   void setInterpolator(Teuchos::RCP<InterpolatorBase<Scalar> > interpolator);
00084   
00086   Teuchos::RCP<InterpolatorBase<Scalar> > unsetInterpolator();
00087 
00089 
00092 
00094   void setSolver(
00095     const Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> > &solver
00096     );
00097 
00099   Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> >
00100   getNonconstSolver();
00101 
00103   Teuchos::RCP<const Thyra::NonlinearSolverBase<Scalar> >
00104   getSolver() const;
00105 
00107 
00110  
00112   bool supportsCloning() const;
00113 
00121   Teuchos::RCP<StepperBase<Scalar> > cloneStepperAlgorithm() const;
00122 
00124   void setModel(const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > &model);
00125   
00127   Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >
00128   getModel() const;
00129 
00131   void setInitialCondition(
00132     const Thyra::ModelEvaluatorBase::InArgs<Scalar> &initialCondition
00133     );
00134 
00136   Scalar takeStep(Scalar dt, StepSizeType flag);
00137   
00139   const StepStatus<Scalar> getStepStatus() const;
00140   
00142 
00145 
00147   Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
00148   get_x_space() const;
00149 
00151   void addPoints(
00152     const Array<Scalar>& time_vec,
00153     const Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >& x_vec,
00154     const Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >& xdot_vec
00155     );
00156   
00158   TimeRange<Scalar> getTimeRange() const;
00159   
00161   void getPoints(
00162     const Array<Scalar>& time_vec,
00163     Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >* x_vec,
00164     Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec,
00165     Array<ScalarMag>* accuracy_vec
00166     ) const;
00167   
00169   void getNodes(Array<Scalar>* time_vec) const;
00170   
00172   void removeNodes(Array<Scalar>& time_vec);
00173 
00175   int getOrder() const;
00176 
00178   
00181 
00183   void setParameterList(Teuchos::RCP<Teuchos::ParameterList> const& paramList);
00184   
00186   Teuchos::RCP<Teuchos::ParameterList> getParameterList();
00187   
00189   Teuchos::RCP<Teuchos::ParameterList> unsetParameterList();
00190   
00192   Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const;
00193  
00195 
00198   
00200   void describe(
00201     Teuchos::FancyOStream  &out,
00202     const Teuchos::EVerbosityLevel verbLevel
00203     ) const;
00204 
00206 
00207 private:
00208 
00209   // ///////////////////////
00210   // Private date members
00211 
00212   bool isInitialized_;
00213   Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > model_;
00214   Thyra::ModelEvaluatorBase::InArgs<Scalar> basePoint_;
00215   Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> > solver_;
00216   Teuchos::RCP<Thyra::VectorBase<Scalar> > x_;
00217   Teuchos::RCP<Thyra::VectorBase<Scalar> > x_dot_;
00218   Teuchos::RCP<Thyra::VectorBase<Scalar> > scaled_x_old_;
00219   Teuchos::RCP<Thyra::VectorBase<Scalar> > x_dot_old_;
00220   Scalar t_;
00221   Scalar t_old_;
00222   Scalar dt_;
00223   int numSteps_;
00224 
00225   Teuchos::RCP<Rythmos::SingleResidualModelEvaluator<Scalar> >  neModel_;
00226 
00227   Teuchos::RCP<Teuchos::ParameterList> parameterList_;
00228 
00229   Teuchos::RCP<InterpolatorBase<Scalar> > interpolator_;
00230 
00231 
00232   // //////////////////////////
00233   // Private member functions
00234 
00235   void initialize_();
00236 
00237 };
00238 
00239 
00240 // ////////////////////////////
00241 // Defintions
00242 
00243 
00244 // Constructors, intializers, Misc.
00245 
00246 
00247 template<class Scalar>
00248 BackwardEulerStepper<Scalar>::BackwardEulerStepper()
00249   :t_(-1.0)
00250   ,t_old_(0.0)
00251 {
00252   numSteps_ = 0;
00253   isInitialized_=false;
00254 }
00255 
00256 
00257 template<class Scalar>
00258 BackwardEulerStepper<Scalar>::BackwardEulerStepper(
00259   const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > &model
00260   ,const Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> > &solver
00261   )
00262   :t_(-1.0)
00263   ,t_old_(0.0)
00264 {
00265   setModel(model);
00266   setSolver(solver);
00267   numSteps_ = 0;
00268   isInitialized_ = false;
00269 }
00270 
00271 
00272 template<class Scalar>
00273 void BackwardEulerStepper<Scalar>::setInterpolator(
00274   Teuchos::RCP<InterpolatorBase<Scalar> > interpolator
00275   )
00276 {
00277 #ifdef TEUCHOS_DEBUG
00278   TEST_FOR_EXCEPT(is_null(interpolator));
00279 #endif
00280   interpolator_ = interpolator;
00281   isInitialized_ = false;
00282 }
00283 
00284 
00285 template<class Scalar>
00286 Teuchos::RCP<InterpolatorBase<Scalar> >
00287 BackwardEulerStepper<Scalar>::unsetInterpolator()
00288 {
00289   Teuchos::RCP<InterpolatorBase<Scalar> > temp_interpolator = interpolator_;
00290   interpolator_ = Teuchos::null;
00291   return(temp_interpolator);
00292   isInitialized_ = false;
00293 }
00294 
00295 
00296 // Overridden from SolverAcceptingStepperBase
00297 
00298 
00299 template<class Scalar>
00300 void BackwardEulerStepper<Scalar>::setSolver(
00301   const Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> > &solver
00302   )
00303 {
00304   using Teuchos::as;
00305 
00306   TEST_FOR_EXCEPT(solver == Teuchos::null)
00307 
00308   Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
00309   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00310   Teuchos::OSTab ostab(out,1,"BES::setSolver");
00311   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00312     *out << "solver = " << solver->description() << std::endl;
00313   }
00314 
00315   solver_ = solver;
00316 
00317   isInitialized_ = false;
00318 
00319 }
00320 
00321 
00322 template<class Scalar>
00323 Teuchos::RCP<Thyra::NonlinearSolverBase<Scalar> >
00324 BackwardEulerStepper<Scalar>::getNonconstSolver()
00325 {
00326   return solver_;
00327 }
00328 
00329 
00330 template<class Scalar>
00331 Teuchos::RCP<const Thyra::NonlinearSolverBase<Scalar> >
00332 BackwardEulerStepper<Scalar>::getSolver() const
00333 {
00334   return solver_;
00335 }
00336 
00337 
00338 // Overridden from StepperBase
00339  
00340 
00341 template<class Scalar>
00342 bool BackwardEulerStepper<Scalar>::supportsCloning() const
00343 {
00344   return true;
00345 }
00346 
00347 
00348 template<class Scalar>
00349 Teuchos::RCP<StepperBase<Scalar> >
00350 BackwardEulerStepper<Scalar>::cloneStepperAlgorithm() const
00351 {
00352   Teuchos::RCP<BackwardEulerStepper<Scalar> >
00353     stepper = Teuchos::rcp(new BackwardEulerStepper<Scalar>);
00354   stepper->isInitialized_ = isInitialized_;
00355   stepper->model_ = model_; // Model is stateless so shallow copy is okay!
00356   if (!is_null(solver_))
00357     stepper->solver_ = solver_->cloneNonlinearSolver().assert_not_null();
00358   if (!is_null(x_))
00359     stepper->x_ = x_->clone_v().assert_not_null();
00360   if (!is_null(scaled_x_old_))
00361     stepper->scaled_x_old_ = scaled_x_old_->clone_v().assert_not_null();
00362   stepper->t_ = t_;
00363   stepper->t_old_ = t_old_;
00364   stepper->dt_ = dt_;
00365   stepper->numSteps_ = numSteps_;
00366   if (!is_null(neModel_))
00367     stepper->neModel_
00368     = Teuchos::rcp(new Rythmos::SingleResidualModelEvaluator<Scalar>);
00369   if (!is_null(parameterList_))
00370     stepper->parameterList_ = parameterList(*parameterList_);
00371   if (!is_null(interpolator_))
00372     stepper->interpolator_
00373       = interpolator_->cloneInterpolator().assert_not_null(); // ToDo: Implement cloneInterpolator()
00374   return stepper;
00375 }
00376 
00377 
00378 template<class Scalar>
00379 void BackwardEulerStepper<Scalar>::setModel(
00380   const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > &model
00381   )
00382 {
00383 
00384   using Teuchos::as;
00385 
00386   TEST_FOR_EXCEPT( is_null(model) );
00387 
00388   Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
00389   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00390   Teuchos::OSTab ostab(out,1,"BES::setModel");
00391   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00392     *out << "model = " << model->description() << std::endl;
00393   }
00394   model_ = model;
00395 
00396   // Wipe out x.  This will either be set thorugh setInitialCondition(...) or
00397   // it will be taken from the model's nominal vlaues!
00398   x_ = Teuchos::null;
00399   scaled_x_old_ = Teuchos::null;
00400   x_dot_ = Teuchos::null;
00401   x_dot_old_ = Teuchos::null;
00402 
00403   isInitialized_ = false;
00404   
00405 }
00406 
00407 
00408 template<class Scalar>
00409 Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >
00410 BackwardEulerStepper<Scalar>::getModel() const
00411 {
00412   return model_;
00413 }
00414 
00415 
00416 template<class Scalar>
00417 void BackwardEulerStepper<Scalar>::setInitialCondition(
00418   const Thyra::ModelEvaluatorBase::InArgs<Scalar> &initialCondition
00419   )
00420 {
00421 
00422   typedef Teuchos::ScalarTraits<Scalar> ST;
00423   typedef Thyra::ModelEvaluatorBase MEB;
00424 
00425   TEST_FOR_EXCEPT( is_null(model_) );
00426 
00427   basePoint_ = initialCondition;
00428 
00429   // x
00430 
00431   Teuchos::RCP<const Thyra::VectorBase<Scalar> >
00432     x_init = initialCondition.get_x();
00433 
00434 #ifdef TEUCHOS_DEBUG
00435   TEST_FOR_EXCEPTION(
00436     is_null(x_init), std::logic_error,
00437     "Error, if the client passes in an intial condition to setInitialCondition(...),\n"
00438     "then x can not be null!" );
00439   THYRA_ASSERT_VEC_SPACES(
00440     "Rythmos::BackwardEulerStepper::setInitialCondition(...)",
00441     *x_init->space(), *model_->get_x_space() );
00442 #endif
00443 
00444   x_ = x_init->clone_v();
00445 
00446   // x_dot
00447 
00448   x_dot_ = createMember(model_->get_x_space());
00449 
00450   Teuchos::RCP<const Thyra::VectorBase<Scalar> >
00451     x_dot_init = initialCondition.get_x_dot();
00452 
00453   if (!is_null(x_dot_init))
00454     assign(&*x_dot_,*x_dot_init);
00455   else
00456     assign(&*x_dot_,ST::zero());
00457   
00458   // t
00459   
00460   t_ = (
00461     initialCondition.supports(MEB::IN_ARG_t)
00462     ? initialCondition.get_t()
00463     : ST::zero()
00464     );
00465   t_old_ = t_;
00466 
00467 }
00468 
00469 
00470 template<class Scalar>
00471 Scalar BackwardEulerStepper<Scalar>::takeStep(Scalar dt, StepSizeType flag)
00472 {
00473 
00474   using Teuchos::as;
00475   using Teuchos::incrVerbLevel;
00476   typedef Teuchos::ScalarTraits<Scalar> ST;
00477   typedef Thyra::NonlinearSolverBase<Scalar> NSB;
00478   typedef Teuchos::VerboseObjectTempState<NSB> VOTSNSB;
00479 
00480   Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
00481   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00482   Teuchos::OSTab ostab(out,1,"BES::takeStep");
00483   VOTSNSB solver_outputTempState(solver_,out,incrVerbLevel(verbLevel,-1));
00484 
00485   if ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
00486     *out
00487       << "\nEntering " << Teuchos::TypeNameTraits<BackwardEulerStepper<Scalar> >::name()
00488       << "::takeStep("<<dt<<","<<toString(flag)<<") ...\n"; 
00489   }
00490 
00491   if (!isInitialized_) {
00492     initialize_(); // This will set x_ if not already set!
00493   }
00494 
00495   dt_ = dt;
00496 
00497   if ((flag == STEP_TYPE_VARIABLE) || (dt == ST::zero())) {
00498     if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) )
00499       *out << "\nThe arguments to takeStep are not valid for BackwardEulerStepper at this time." << std::endl;
00500     // print something out about this method not supporting automatic variable step-size
00501     return(Scalar(-ST::one()));
00502   }
00503   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00504     *out << "\ndt = " << dt << std::endl;
00505   }
00506 
00507 
00508   //
00509   // Setup the nonlinear equations:
00510   //
00511   //   f( (1/dt)* x + (-1/dt)*x_old), x, t ) = 0
00512   //
00513 
00514   V_StV( &*scaled_x_old_, Scalar(-ST::one()/dt), *x_ );
00515   t_old_ = t_;
00516   if(!neModel_.get()) {
00517     neModel_ = Teuchos::rcp(new Rythmos::SingleResidualModelEvaluator<Scalar>());
00518   }
00519   neModel_->initializeSingleResidualModel(
00520     model_, basePoint_,
00521     Scalar(ST::one()/dt), scaled_x_old_,
00522     ST::one(), Teuchos::null,
00523     t_old_+dt,
00524     Teuchos::null
00525     );
00526   if( solver_->getModel().get() != neModel_.get() ) {
00527     solver_->setModel(neModel_);
00528   }
00529   // 2007/05/18: rabartl: ToDo: Above, set the stream and the verbosity level
00530   // on solver_ so that we an see what it is doing!
00531 
00532   //
00533   // Solve the implicit nonlinear system to a tolerance of ???
00534   //
00535   
00536   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
00537     *out << "\nSolving the implicit backward-Euler timestep equation ...\n";
00538   }
00539 
00540   Thyra::SolveStatus<Scalar>
00541     neSolveStatus = solver_->solve(&*x_);
00542 
00543   // In the above solve, on input *x_ is the old value of x for the previous
00544   // time step which is used as the initial guess for the solver.  On output,
00545   // *x_ is the converged timestep solution.
00546  
00547   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
00548     *out << "\nOutput status of nonlinear solve:\n" << neSolveStatus;
00549   }
00550 
00551   // 2007/05/18: rabartl: ToDo: Above, get the solve status from the above
00552   // solve and at least print warning message if the solve fails!  Actually,
00553   // you should most likely thrown an exception if the solve fails or return
00554   // false if appropriate
00555 
00556   //
00557   // Update the step
00558   //
00559 
00560   assign( &*x_dot_old_, *x_dot_ );
00561 
00562   // x_dot = (1/dt)*x - (1/dt)*x_old 
00563   V_StV( &*x_dot_, Scalar(ST::one()/dt), *x_ );
00564   Vp_StV( &*x_dot_, Scalar(ST::one()), *scaled_x_old_ );
00565 
00566   t_ += dt;
00567 
00568   numSteps_++;
00569 
00570   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00571     *out << "\nt_old_ = " << t_old_ << std::endl;
00572     *out << "\nt_ = " << t_ << std::endl;
00573   }
00574 
00575 #ifdef TEUCHOS_DEBUG
00576 
00577   if ( includesVerbLevel(verbLevel,Teuchos::VERB_LOW) )
00578     *out << "\nChecking to make sure that solution and the interpolated solution are the same! ...\n";
00579 
00580   {
00581 
00582     typedef ScalarTraits<Scalar> ST;
00583     typedef typename ST::magnitudeType ScalarMag;
00584     typedef ScalarTraits<ScalarMag> SMT;
00585     
00586     Teuchos::OSTab tab(out);
00587 
00588     const StepStatus<Scalar> stepStatus = this->getStepStatus();
00589 
00590     RCP<const Thyra::VectorBase<Scalar> >
00591       x = stepStatus.solution,
00592       xdot = stepStatus.solutionDot;
00593 
00594     Array<Scalar> time_vec = Teuchos::tuple(stepStatus.time);
00595     Array<RCP<const Thyra::VectorBase<Scalar> > > x_vec, xdot_vec;
00596     this->getPoints(time_vec,&x_vec,&xdot_vec,0);
00597 
00598     RCP<const Thyra::VectorBase<Scalar> >
00599       x_interp = x_vec[0],
00600       xdot_interp = xdot_vec[0];
00601 
00602     TEST_FOR_EXCEPT(
00603       !Thyra::testRelNormDiffErr(
00604         "x", *x, "x_interp", *x_interp,
00605         "2*epsilon", ScalarMag(100.0*SMT::eps()),
00606         "2*epsilon", ScalarMag(100.0*SMT::eps()),
00607         includesVerbLevel(verbLevel,Teuchos::VERB_HIGH) ? out.get() : 0
00608         )
00609       );
00610 
00611     TEST_FOR_EXCEPT(
00612       !Thyra::testRelNormDiffErr(
00613         "xdot", *xdot, "xdot_interp", *xdot_interp,
00614         "2*epsilon", ScalarMag(100.0*SMT::eps()),
00615         "2*epsilon", ScalarMag(100.0*SMT::eps()),
00616         includesVerbLevel(verbLevel,Teuchos::VERB_HIGH) ? out.get() : 0
00617         )
00618       );
00619 
00620   }
00621 
00622   // 2007/07/25: rabartl: ToDo: Move the above test into a helper function so
00623   // that it can be used from lots of different places!
00624 
00625 #endif
00626 
00627   if ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
00628     *out
00629       << "\nLeaving " << Teuchos::TypeNameTraits<BackwardEulerStepper<Scalar> >::name()
00630       << "::takeStep(...) ...\n"; 
00631   }
00632 
00633   return(dt);
00634 
00635 }
00636 
00637 
00638 template<class Scalar>
00639 const StepStatus<Scalar> BackwardEulerStepper<Scalar>::getStepStatus() const
00640 {
00641 
00642   typedef Teuchos::ScalarTraits<Scalar> ST;
00643 
00644   StepStatus<Scalar> stepStatus;
00645 
00646   if (!isInitialized_) {
00647     stepStatus.message = "This stepper is uninitialized.";
00648     stepStatus.stepStatus = STEP_STATUS_UNINITIALIZED;
00649     stepStatus.stepSize = Scalar(-ST::one());
00650     stepStatus.order = -1;
00651     stepStatus.time = Scalar(-ST::one());
00652     // 2007/05/21: rabartl: ToDo: Need to return something other than -1.0
00653     // since negative time might be just fine.
00654     stepStatus.solution = Teuchos::null;
00655     // 2007/05/18: rabartl: I changed the above to null since just returning a
00656     // vector fully of junk makes no sense.
00657     stepStatus.solutionDot = Teuchos::null;
00658     return(stepStatus);
00659   }
00660 
00661   if (numSteps_ > 0) {
00662     stepStatus.stepStatus = STEP_STATUS_CONVERGED; 
00663   }
00664 
00665   stepStatus.stepSize = dt_;
00666   stepStatus.order = 1;
00667   stepStatus.time = t_;
00668   stepStatus.solution = x_;
00669   stepStatus.solutionDot = x_dot_;
00670 
00671   return(stepStatus);
00672 }
00673 
00674 
00675 // Overridden from InterpolationBufferBase
00676 
00677 
00678 template<class Scalar>
00679 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
00680 BackwardEulerStepper<Scalar>::get_x_space() const
00681 {
00682   return ( !is_null(model_) ? model_->get_x_space() : Teuchos::null );
00683 }
00684 
00685 
00686 template<class Scalar>
00687 void BackwardEulerStepper<Scalar>::addPoints(
00688     const Array<Scalar>& time_vec
00689     ,const Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >& x_vec
00690     ,const Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >& xdot_vec
00691     )
00692 {
00693   typedef Teuchos::ScalarTraits<Scalar> ST;
00694 #ifdef TEUCHOS_DEBUG
00695   TEST_FOR_EXCEPTION(!isInitialized_,std::logic_error,"Error, attempting to call addPoints before the stepper is initialized!\n");
00696   TEST_FOR_EXCEPTION(time_vec.size() == 0, std::logic_error, "Error, addPoints called with an empty time_vec array!\n");
00697 #endif // TEUCHOS_DEBUG
00698   using Teuchos::as;
00699   Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
00700   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00701   Teuchos::OSTab ostab(out,1,"BES::setPoints");
00702   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00703     *out << "time_vec = " << std::endl;
00704     for (int i=0 ; i<Teuchos::as<int>(time_vec.size()) ; ++i) {
00705       *out << "time_vec[" << i << "] = " << time_vec[i] << std::endl;
00706     }
00707   } else if (time_vec.size() == 1) {
00708     int n = 0;
00709     t_ = time_vec[n];
00710     t_old_ = t_;
00711     Thyra::V_V(&*x_,*x_vec[n]);
00712     Thyra::V_V(&*scaled_x_old_,*x_);
00713   } else {
00714     int n = time_vec.size()-1;
00715     int nm1 = time_vec.size()-2;
00716     t_ = time_vec[n];
00717     t_old_ = time_vec[nm1];
00718     Thyra::V_V(&*x_,*x_vec[n]);
00719     Scalar dt = t_ - t_old_;
00720     Thyra::V_StV(&*scaled_x_old_,Scalar(-ST::one()/dt),*x_vec[nm1]);
00721   }
00722 }
00723 
00724 
00725 
00726 template<class Scalar>
00727 TimeRange<Scalar> BackwardEulerStepper<Scalar>::getTimeRange() const
00728 {
00729   return timeRange(t_old_,t_);
00730 }
00731 
00732 
00733 template<class Scalar>
00734 void BackwardEulerStepper<Scalar>::getPoints(
00735     const Array<Scalar>& time_vec
00736     ,Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >* x_vec
00737     ,Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec
00738     ,Array<ScalarMag>* accuracy_vec) const
00739 {
00740 
00741   using Teuchos::as;
00742   typedef Teuchos::ScalarTraits<Scalar> ST;
00743   typedef typename ST::magnitudeType ScalarMag;
00744   typedef Teuchos::ScalarTraits<ScalarMag> SMT;
00745   typename DataStore<Scalar>::DataStoreVector_t ds_nodes;
00746   typename DataStore<Scalar>::DataStoreVector_t ds_out;
00747 
00748 #ifdef TEUCHOS_DEBUG
00749   TEST_FOR_EXCEPT( 0 == x_vec );
00750   TEST_FOR_EXCEPTION(!isInitialized_,std::logic_error,"Error, attempting to call getPoints before the stepper is initialized!\n");
00751 #endif
00752   Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
00753   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00754   Teuchos::OSTab ostab(out,1,"BES::getPoints");
00755   if ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
00756     *out
00757       << "\nEntering " << Teuchos::TypeNameTraits<BackwardEulerStepper<Scalar> >::name()
00758       << "::getPoints(...) ...\n"; 
00759   }
00760   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00761     for (int i=0 ; i<Teuchos::as<int>(time_vec.size()) ; ++i) {
00762       *out << "time_vec[" << i << "] = " << time_vec[i] << std::endl;
00763     }
00764     *out << "I can interpolate in the interval [" << t_old_ << "," << t_ << "]." << std::endl;
00765   }
00766   if (t_old_ != t_) {
00767     if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00768       *out << "Passing two points to interpolator:  " << t_old_ << " and " << t_ << std::endl;
00769     }
00770     DataStore<Scalar> ds_temp;
00771     Scalar dt = t_ - t_old_;
00772 #ifdef TEUCHOS_DEBUG
00773     TEST_FOR_EXCEPT(
00774       !Thyra::testRelErr(
00775         "dt", dt, "dt_", dt_,
00776         "1e+4*epsilon", ScalarMag(1e+4*SMT::eps()),
00777         "1e+2*epsilon", ScalarMag(1e+2*SMT::eps()),
00778         as<int>(verbLevel) >= as<int>(Teuchos::VERB_MEDIUM) ? out.get() : 0
00779         )
00780       );
00781 #endif
00782     Teuchos::RCP<Thyra::VectorBase<Scalar> >
00783       x_temp = scaled_x_old_->clone_v();
00784     Thyra::Vt_S(&*x_temp,Scalar(-ST::one()*dt));  // undo the scaling
00785     ds_temp.time = t_old_;
00786     ds_temp.x = x_temp;
00787     ds_temp.xdot = x_dot_old_;
00788     ds_temp.accuracy = ScalarMag(dt);
00789     ds_nodes.push_back(ds_temp);
00790     ds_temp.time = t_;
00791     ds_temp.x = x_;
00792     ds_temp.xdot = x_dot_;
00793     ds_temp.accuracy = ScalarMag(dt);
00794     ds_nodes.push_back(ds_temp);
00795   }
00796   else {
00797     if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00798       *out << "Passing one point to interpolator:  " << t_ << std::endl;
00799     }
00800     DataStore<Scalar> ds_temp;
00801     ds_temp.time = t_;
00802     ds_temp.x = x_;
00803     ds_temp.xdot = x_dot_;
00804     ds_temp.accuracy = ScalarMag(ST::zero());
00805     ds_nodes.push_back(ds_temp);
00806   }
00807   interpolator_->interpolate(ds_nodes,time_vec,&ds_out);
00808   Array<Scalar> time_out;
00809   dataStoreVectorToVector(ds_out,&time_out,x_vec,xdot_vec,accuracy_vec);
00810   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00811     *out << "Passing out the interpolated values:" << std::endl;
00812     for (int i=0; i<Teuchos::as<int>(time_out.size()) ; ++i) {
00813       *out << "time[" << i << "] = " << time_out[i] << std::endl;
00814       *out << "x_vec[" << i << "] = " << std::endl;
00815       (*x_vec)[i]->describe(*out,Teuchos::VERB_EXTREME);
00816       if (xdot_vec) {
00817         if ( (*xdot_vec)[i] == Teuchos::null) {
00818           *out << "xdot_vec[" << i << "] = Teuchos::null" << std::endl;
00819         }
00820         else {
00821           *out << "xdot_vec[" << i << "] = " << std::endl;
00822           (*xdot_vec)[i]->describe(*out,Teuchos::VERB_EXTREME);
00823         }
00824       }
00825       if(accuracy_vec) {
00826         *out << "accuracy[" << i << "] = " << (*accuracy_vec)[i] << std::endl;
00827       }
00828     }
00829   }
00830   if ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
00831     *out
00832       << "Leaving " << Teuchos::TypeNameTraits<BackwardEulerStepper<Scalar> >::name()
00833       << "::getPoints(...) ...\n"; 
00834   }
00835 }
00836 
00837 
00838 template<class Scalar>
00839 void BackwardEulerStepper<Scalar>::getNodes(Array<Scalar>* time_vec) const
00840 {
00841   using Teuchos::as;
00842 #ifdef TEUCHOS_DEBUG
00843   TEST_FOR_EXCEPTION(!isInitialized_,std::logic_error,"Error, attempting to call getNodes before stepper is initialized!\n");
00844 #endif // TEUCHOS_DEBUG
00845   time_vec->clear();
00846   time_vec->push_back(t_old_);
00847   if (numSteps_ > 0) {
00848     time_vec->push_back(t_);
00849   }
00850   Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
00851   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00852   Teuchos::OSTab ostab(out,1,"BES::getNodes");
00853   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00854     *out << this->description() << std::endl;
00855     for (int i=0 ; i<Teuchos::as<int>(time_vec->size()) ; ++i) {
00856       *out << "time_vec[" << i << "] = " << (*time_vec)[i] << std::endl;
00857     }
00858   }
00859 }
00860 
00861 
00862 template<class Scalar>
00863 void BackwardEulerStepper<Scalar>::removeNodes(Array<Scalar>& time_vec) 
00864 {
00865 #ifdef TEUCHOS_DEBUG
00866   TEST_FOR_EXCEPTION(!isInitialized_,std::logic_error,"Error, attempting to call removeNodes before stepper is initialized!\n");
00867 #endif // TEUCHOS_DEBUG
00868   using Teuchos::as;
00869   Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
00870   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00871   Teuchos::OSTab ostab(out,1,"BES::removeNodes");
00872   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00873     *out << "time_vec = " << std::endl;
00874     for (int i=0 ; i<Teuchos::as<int>(time_vec.size()) ; ++i) {
00875       *out << "time_vec[" << i << "] = " << time_vec[i] << std::endl;
00876     }
00877   }
00878   TEST_FOR_EXCEPTION(true,std::logic_error,"Error, removeNodes is not implemented for BackwardEulerStepper at this time.\n");
00879   // TODO:
00880   // if any time in time_vec matches t_ or t_old_, then do the following:
00881   // remove t_old_:  set t_old_ = t_ and set scaled_x_old_ = x_
00882   // remove t_:  set t_ = t_old_ and set x_ = -dt*scaled_x_old_
00883 }
00884 
00885 
00886 template<class Scalar>
00887 int BackwardEulerStepper<Scalar>::getOrder() const
00888 {
00889   return(1);
00890 }
00891 
00892 
00893 // Overridden from Teuchos::ParameterListAcceptor
00894 
00895 
00896 template <class Scalar>
00897 void BackwardEulerStepper<Scalar>::setParameterList(
00898   Teuchos::RCP<Teuchos::ParameterList> const& paramList
00899   )
00900 {
00901   TEST_FOR_EXCEPT(is_null(paramList));
00902   paramList->validateParametersAndSetDefaults(*this->getValidParameters());
00903   parameterList_ = paramList;
00904   Teuchos::readVerboseObjectSublist(&*parameterList_,this);
00905 }
00906 
00907 
00908 template <class Scalar>
00909 Teuchos::RCP<Teuchos::ParameterList>
00910 BackwardEulerStepper<Scalar>::getParameterList()
00911 {
00912   return(parameterList_);
00913 }
00914 
00915 
00916 template <class Scalar>
00917 Teuchos::RCP<Teuchos::ParameterList>
00918 BackwardEulerStepper<Scalar>::unsetParameterList()
00919 {
00920   Teuchos::RCP<Teuchos::ParameterList>
00921     temp_param_list = parameterList_;
00922   parameterList_ = Teuchos::null;
00923   return(temp_param_list);
00924 }
00925 
00926 
00927 template<class Scalar>
00928 Teuchos::RCP<const Teuchos::ParameterList>
00929 BackwardEulerStepper<Scalar>::getValidParameters() const
00930 {
00931   using Teuchos::RCP;
00932   using Teuchos::ParameterList;
00933   static RCP<const ParameterList> validPL;
00934   if (is_null(validPL)) {
00935     RCP<ParameterList> pl = Teuchos::parameterList();
00936     Teuchos::setupVerboseObjectSublist(&*pl);
00937     validPL = pl;
00938   }
00939   return validPL;
00940 }
00941 
00942 
00943 // Overridden from Teuchos::Describable
00944 
00945 
00946 template<class Scalar>
00947 void BackwardEulerStepper<Scalar>::describe(
00948       Teuchos::FancyOStream                &out
00949       ,const Teuchos::EVerbosityLevel      verbLevel
00950       ) const
00951 {
00952   using Teuchos::as;
00953   Teuchos::OSTab tab(out);
00954   if (!isInitialized_) {
00955     out << this->description() << " : This stepper is not initialized yet" << std::endl;
00956     return;
00957   }
00958   if (
00959     as<int>(verbLevel) == as<int>(Teuchos::VERB_DEFAULT)
00960     ||
00961     as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW)
00962     )
00963   {
00964     out << this->description() << "::describe:" << std::endl;
00965     out << "model = " << model_->description() << std::endl;
00966     out << "solver = " << solver_->description() << std::endl;
00967     if (neModel_ == Teuchos::null) {
00968       out << "neModel = Teuchos::null" << std::endl;
00969     } else {
00970       out << "neModel = " << neModel_->description() << std::endl;
00971     }
00972   }
00973   else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW)) {
00974     out << "t_ = " << t_ << std::endl;
00975     out << "t_old_ = " << t_old_ << std::endl;
00976   }
00977   else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_MEDIUM)) {
00978   }
00979   else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH)) {
00980     out << "model_ = " << std::endl;
00981     model_->describe(out,verbLevel);
00982     out << "solver_ = " << std::endl;
00983     solver_->describe(out,verbLevel);
00984     if (neModel_ == Teuchos::null) {
00985       out << "neModel = Teuchos::null" << std::endl;
00986     } else {
00987       out << "neModel = " << std::endl;
00988       neModel_->describe(out,verbLevel);
00989     }
00990     out << "x_ = " << std::endl;
00991     x_->describe(out,verbLevel);
00992     out << "scaled_x_old_ = " << std::endl;
00993     scaled_x_old_->describe(out,verbLevel);
00994   }
00995 }
00996 
00997 
00998 // private
00999 
01000 
01001 template <class Scalar>
01002 void BackwardEulerStepper<Scalar>::initialize_()
01003 {
01004   TEST_FOR_EXCEPT(model_ == Teuchos::null);
01005   TEST_FOR_EXCEPT(solver_ == Teuchos::null);
01006 
01007   if (is_null(x_)) {
01008     // If x has not been set then setInitialCondition(...) was not
01009     // called by the client so the model had better have the 
01010     // initial condition!
01011     this->setInitialCondition(model_->getNominalValues());
01012   }
01013 
01014   if ( is_null(interpolator_) ) {
01015     // If an interpolator has not been explicitly set, then just create
01016     // a default linear interpolator.
01017     interpolator_ = Teuchos::rcp(new LinearInterpolator<Scalar> );
01018     // 2007/05/18: rabartl: ToDo: Replace this with a Hermete interplator
01019     // when it is implementated!
01020   }
01021 
01022   if (is_null(scaled_x_old_))
01023     scaled_x_old_ = createMember(model_->get_x_space());
01024 
01025   if (is_null(x_dot_))
01026     x_dot_ = createMember(model_->get_x_space());
01027 
01028   if (is_null(x_dot_old_))
01029     x_dot_old_ = createMember(model_->get_x_space());
01030 
01031   isInitialized_ = true;
01032 
01033 }
01034 
01035 
01036 } // namespace Rythmos
01037 
01038 #endif //Rythmos_BACKWARD_EULER_STEPPER_H

Generated on Tue Oct 20 12:46:07 2009 for Rythmos - Transient Integration for Differential Equations by doxygen 1.4.7