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