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<const 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<const 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   RCP<Rythmos::SingleResidualModelEvaluator<Scalar> >  neModel_ = Teuchos::null;
00115   RCP<Teuchos::ParameterList> parameterList_ = Teuchos::null;
00116   RCP<InterpolatorBase<Scalar> > 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 RYTHMOS_DEBUG
00127   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   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   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 
00280 template<class Scalar>
00281 RCP<const Thyra::ModelEvaluator<Scalar> >
00282 BackwardEulerStepper<Scalar>::getModel() const
00283 {
00284   return model_;
00285 }
00286 
00287 
00288 template<class Scalar>
00289 void BackwardEulerStepper<Scalar>::setInitialCondition(
00290   const Thyra::ModelEvaluatorBase::InArgs<Scalar> &initialCondition
00291   )
00292 {
00293 
00294   typedef Teuchos::ScalarTraits<Scalar> ST;
00295   typedef Thyra::ModelEvaluatorBase MEB;
00296 
00297   basePoint_ = initialCondition;
00298 
00299   // x
00300 
00301   RCP<const Thyra::VectorBase<Scalar> >
00302     x_init = initialCondition.get_x();
00303 
00304 #ifdef RYTHMOS_DEBUG
00305   TEST_FOR_EXCEPTION(
00306     is_null(x_init), std::logic_error,
00307     "Error, if the client passes in an intial condition to setInitialCondition(...),\n"
00308     "then x can not be null!" );
00309 #endif
00310 
00311   x_ = x_init->clone_v();
00312 
00313   // x_dot
00314 
00315   RCP<const Thyra::VectorBase<Scalar> >
00316     x_dot_init = initialCondition.get_x_dot();
00317 
00318   if (!is_null(x_dot_init)) {
00319     x_dot_ = x_dot_init->clone_v();
00320   }
00321   else {
00322     x_dot_ = createMember(x_->space());
00323     assign(&*x_dot_,ST::zero());
00324   }
00325   
00326   // t
00327   
00328   t_ = initialCondition.get_t();
00329 
00330   t_old_ = t_;
00331 
00332   // x_old 
00333 
00334   scaled_x_old_ = x_->clone_v();
00335 
00336   // x_dot_old
00337   
00338   x_dot_old_ = x_dot_->clone_v();
00339 
00340 
00341   haveInitialCondition_ = true;
00342 
00343 }
00344 
00345 
00346 template<class Scalar>
00347 Scalar BackwardEulerStepper<Scalar>::takeStep(Scalar dt, StepSizeType stepSizeType)
00348 {
00349 
00350   using Teuchos::as;
00351   using Teuchos::incrVerbLevel;
00352   typedef Teuchos::ScalarTraits<Scalar> ST;
00353   typedef Thyra::NonlinearSolverBase<Scalar> NSB;
00354   typedef Teuchos::VerboseObjectTempState<NSB> VOTSNSB;
00355 
00356   initialize();
00357 
00358   RCP<Teuchos::FancyOStream> out = this->getOStream();
00359   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00360   Teuchos::OSTab ostab(out,1,"BES::takeStep");
00361   VOTSNSB solver_outputTempState(solver_,out,incrVerbLevel(verbLevel,-1));
00362 
00363   if ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
00364     *out
00365       << "\nEntering " << Teuchos::TypeNameTraits<BackwardEulerStepper<Scalar> >::name()
00366       << "::takeStep("<<dt<<","<<toString(stepSizeType)<<") ...\n"; 
00367   }
00368 
00369   dt_ = dt;
00370 
00371   if ((stepSizeType == STEP_TYPE_VARIABLE) || (dt == ST::zero())) {
00372     if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) )
00373       *out << "\nThe arguments to takeStep are not valid for BackwardEulerStepper at this time." << std::endl;
00374     // print something out about this method not supporting automatic variable step-size
00375     return(Scalar(-ST::one()));
00376   }
00377   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00378     *out << "\ndt = " << dt << std::endl;
00379   }
00380 
00381 
00382   //
00383   // Setup the nonlinear equations:
00384   //
00385   //   f( (1/dt)* x + (-1/dt)*x_old), x, t ) = 0
00386   //
00387 
00388   V_StV( &*scaled_x_old_, Scalar(-ST::one()/dt), *x_ );
00389   t_old_ = t_;
00390   if(!neModel_.get()) {
00391     neModel_ = Teuchos::rcp(new Rythmos::SingleResidualModelEvaluator<Scalar>());
00392   }
00393   neModel_->initializeSingleResidualModel(
00394     model_, basePoint_,
00395     Scalar(ST::one()/dt), scaled_x_old_,
00396     ST::one(), Teuchos::null,
00397     t_old_+dt,
00398     Teuchos::null
00399     );
00400   if( solver_->getModel().get() != neModel_.get() ) {
00401     solver_->setModel(neModel_);
00402   }
00403   // 2007/05/18: rabartl: ToDo: Above, set the stream and the verbosity level
00404   // on solver_ so that we an see what it is doing!
00405 
00406   //
00407   // Solve the implicit nonlinear system to a tolerance of ???
00408   //
00409   
00410   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
00411     *out << "\nSolving the implicit backward-Euler timestep equation ...\n";
00412   }
00413 
00414   Thyra::SolveStatus<Scalar>
00415     neSolveStatus = solver_->solve(&*x_);
00416 
00417   // In the above solve, on input *x_ is the old value of x for the previous
00418   // time step which is used as the initial guess for the solver.  On output,
00419   // *x_ is the converged timestep solution.
00420  
00421   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
00422     *out << "\nOutput status of nonlinear solve:\n" << neSolveStatus;
00423   }
00424 
00425   // 2007/05/18: rabartl: ToDo: Above, get the solve status from the above
00426   // solve and at least print warning message if the solve fails!  Actually,
00427   // you should most likely thrown an exception if the solve fails or return
00428   // false if appropriate
00429 
00430   //
00431   // Update the step
00432   //
00433 
00434   assign( &*x_dot_old_, *x_dot_ );
00435 
00436   // x_dot = (1/dt)*x - (1/dt)*x_old 
00437   V_StV( &*x_dot_, Scalar(ST::one()/dt), *x_ );
00438   Vp_StV( &*x_dot_, Scalar(ST::one()), *scaled_x_old_ );
00439 
00440   t_ += dt;
00441 
00442   numSteps_++;
00443 
00444   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00445     *out << "\nt_old_ = " << t_old_ << std::endl;
00446     *out << "\nt_ = " << t_ << std::endl;
00447   }
00448 
00449 #ifdef RYTHMOS_DEBUG
00450   // 04/14/09 tscoffe: This code should be moved to StepperValidator
00451 
00452   if ( includesVerbLevel(verbLevel,Teuchos::VERB_LOW) )
00453     *out << "\nChecking to make sure that solution and the interpolated solution are the same! ...\n";
00454 
00455   {
00456 
00457     typedef ScalarTraits<Scalar> ST;
00458     typedef typename ST::magnitudeType ScalarMag;
00459     typedef ScalarTraits<ScalarMag> SMT;
00460     
00461     Teuchos::OSTab tab(out);
00462 
00463     const StepStatus<Scalar> stepStatus = this->getStepStatus();
00464 
00465     RCP<const Thyra::VectorBase<Scalar> >
00466       x = stepStatus.solution,
00467       xdot = stepStatus.solutionDot;
00468 
00469     Array<Scalar> time_vec = Teuchos::tuple(stepStatus.time);
00470     Array<RCP<const Thyra::VectorBase<Scalar> > > x_vec, xdot_vec;
00471     this->getPoints(time_vec,&x_vec,&xdot_vec,0);
00472 
00473     RCP<const Thyra::VectorBase<Scalar> >
00474       x_interp = x_vec[0],
00475       xdot_interp = xdot_vec[0];
00476 
00477     TEST_FOR_EXCEPT(
00478       !Thyra::testRelNormDiffErr(
00479         "x", *x, "x_interp", *x_interp,
00480         "2*epsilon", ScalarMag(100.0*SMT::eps()),
00481         "2*epsilon", ScalarMag(100.0*SMT::eps()),
00482         includesVerbLevel(verbLevel,Teuchos::VERB_HIGH) ? out.get() : 0
00483         )
00484       );
00485 
00486     TEST_FOR_EXCEPT(
00487       !Thyra::testRelNormDiffErr(
00488         "xdot", *xdot, "xdot_interp", *xdot_interp,
00489         "2*epsilon", ScalarMag(100.0*SMT::eps()),
00490         "2*epsilon", ScalarMag(100.0*SMT::eps()),
00491         includesVerbLevel(verbLevel,Teuchos::VERB_HIGH) ? out.get() : 0
00492         )
00493       );
00494 
00495   }
00496 
00497   // 2007/07/25: rabartl: ToDo: Move the above test into a helper function so
00498   // that it can be used from lots of different places!
00499 
00500 #endif // RYTHMOS_DEBUG
00501 
00502   if ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
00503     *out
00504       << "\nLeaving " << Teuchos::TypeNameTraits<BackwardEulerStepper<Scalar> >::name()
00505       << "::takeStep(...) ...\n"; 
00506   }
00507 
00508   return(dt);
00509 
00510 }
00511 
00512 
00513 template<class Scalar>
00514 const StepStatus<Scalar> BackwardEulerStepper<Scalar>::getStepStatus() const
00515 {
00516 
00517   typedef Teuchos::ScalarTraits<Scalar> ST;
00518 
00519   StepStatus<Scalar> stepStatus; // Defaults to unknown status
00520 
00521   if (!isInitialized_) {
00522     stepStatus.stepStatus = STEP_STATUS_UNINITIALIZED;
00523   }
00524   else if (numSteps_ > 0) {
00525     stepStatus.stepStatus = STEP_STATUS_CONVERGED; 
00526   }
00527   // else unknown
00528 
00529   stepStatus.stepSize = dt_;
00530   stepStatus.order = 1;
00531   stepStatus.time = t_;
00532   stepStatus.solution = x_;
00533   stepStatus.solutionDot = x_dot_;
00534 
00535   return(stepStatus);
00536 
00537 }
00538 
00539 
00540 // Overridden from InterpolationBufferBase
00541 
00542 
00543 template<class Scalar>
00544 RCP<const Thyra::VectorSpaceBase<Scalar> >
00545 BackwardEulerStepper<Scalar>::get_x_space() const
00546 {
00547   return ( !is_null(model_) ? model_->get_x_space() : Teuchos::null );
00548 }
00549 
00550 
00551 template<class Scalar>
00552 void BackwardEulerStepper<Scalar>::addPoints(
00553   const Array<Scalar>& time_vec,
00554   const Array<RCP<const Thyra::VectorBase<Scalar> > >& x_vec,
00555   const Array<RCP<const Thyra::VectorBase<Scalar> > >& xdot_vec
00556   )
00557 {
00558 
00559   typedef Teuchos::ScalarTraits<Scalar> ST;
00560   using Teuchos::as;
00561 
00562   initialize();
00563 
00564 #ifdef RYTHMOS_DEBUG
00565   TEST_FOR_EXCEPTION(
00566     time_vec.size() == 0, std::logic_error,
00567     "Error, addPoints called with an empty time_vec array!\n");
00568 #endif // RYTHMOS_DEBUG
00569 
00570   RCP<Teuchos::FancyOStream> out = this->getOStream();
00571   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00572   Teuchos::OSTab ostab(out,1,"BES::setPoints");
00573 
00574   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00575     *out << "time_vec = " << std::endl;
00576     for (int i=0 ; i<Teuchos::as<int>(time_vec.size()) ; ++i) {
00577       *out << "time_vec[" << i << "] = " << time_vec[i] << std::endl;
00578     }
00579   }
00580   else if (time_vec.size() == 1) {
00581     int n = 0;
00582     t_ = time_vec[n];
00583     t_old_ = t_;
00584     Thyra::V_V(&*x_,*x_vec[n]);
00585     Thyra::V_V(&*scaled_x_old_,*x_);
00586   }
00587   else {
00588     int n = time_vec.size()-1;
00589     int nm1 = time_vec.size()-2;
00590     t_ = time_vec[n];
00591     t_old_ = time_vec[nm1];
00592     Thyra::V_V(&*x_,*x_vec[n]);
00593     Scalar dt = t_ - t_old_;
00594     Thyra::V_StV(&*scaled_x_old_,Scalar(-ST::one()/dt),*x_vec[nm1]);
00595   }
00596 }
00597 
00598 
00599 template<class Scalar>
00600 TimeRange<Scalar> BackwardEulerStepper<Scalar>::getTimeRange() const
00601 {
00602   if ( !isInitialized_ && haveInitialCondition_ )
00603     return timeRange<Scalar>(t_,t_);
00604   if ( !isInitialized_ && !haveInitialCondition_ )
00605     return invalidTimeRange<Scalar>();
00606   return timeRange<Scalar>(t_old_,t_);
00607 }
00608 
00609 
00610 template<class Scalar>
00611 void BackwardEulerStepper<Scalar>::getPoints(
00612   const Array<Scalar>& time_vec,
00613   Array<RCP<const Thyra::VectorBase<Scalar> > >* x_vec,
00614   Array<RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec,
00615   Array<ScalarMag>* accuracy_vec
00616   ) const
00617 {
00618   typedef Teuchos::ScalarTraits<Scalar> ST;
00619   using Teuchos::constOptInArg;
00620   using Teuchos::ptr;
00621 #ifdef RYTHMOS_DEBUG
00622   TEUCHOS_ASSERT(haveInitialCondition_);
00623 #endif // RYTHMOS_DEBUG
00624   RCP<Thyra::VectorBase<Scalar> > x_temp = x_;
00625   if (compareTimeValues(t_old_,t_)!=0) {
00626     Scalar dt = t_ - t_old_;
00627     x_temp = scaled_x_old_->clone_v();
00628     Thyra::Vt_S(&*x_temp,Scalar(-ST::one()*dt));  // undo the scaling
00629   }
00630   defaultGetPoints<Scalar>(
00631       t_old_, constOptInArg(*x_temp), constOptInArg(*x_dot_old_),
00632       t_, constOptInArg(*x_), constOptInArg(*x_dot_),
00633       time_vec, ptr(x_vec), ptr(xdot_vec), ptr(accuracy_vec),
00634       ptr(interpolator_.get())
00635       );
00636 
00637   /*
00638   using Teuchos::as;
00639   typedef Teuchos::ScalarTraits<Scalar> ST;
00640   typedef typename ST::magnitudeType ScalarMag;
00641   typedef Teuchos::ScalarTraits<ScalarMag> SMT;
00642   typename DataStore<Scalar>::DataStoreVector_t ds_nodes;
00643   typename DataStore<Scalar>::DataStoreVector_t ds_out;
00644 
00645 #ifdef RYTHMOS_DEBUG
00646   TEST_FOR_EXCEPT(!haveInitialCondition_);
00647   TEST_FOR_EXCEPT( 0 == x_vec );
00648 #endif
00649 
00650   RCP<Teuchos::FancyOStream> out = this->getOStream();
00651   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00652   Teuchos::OSTab ostab(out,1,"BES::getPoints");
00653   if ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
00654     *out
00655       << "\nEntering " << Teuchos::TypeNameTraits<BackwardEulerStepper<Scalar> >::name()
00656       << "::getPoints(...) ...\n"; 
00657   }
00658   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00659     for (int i=0 ; i<Teuchos::as<int>(time_vec.size()) ; ++i) {
00660       *out << "time_vec[" << i << "] = " << time_vec[i] << std::endl;
00661     }
00662     *out << "I can interpolate in the interval [" << t_old_ << "," << t_ << "]." << std::endl;
00663   }
00664   if (t_old_ != t_) {
00665     if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00666       *out << "Passing two points to interpolator:  " << t_old_ << " and " << t_ << std::endl;
00667     }
00668     DataStore<Scalar> ds_temp;
00669     Scalar dt = t_ - t_old_;
00670 #ifdef RYTHMOS_DEBUG
00671     TEST_FOR_EXCEPT(
00672       !Thyra::testRelErr(
00673         "dt", dt, "dt_", dt_,
00674         "1e+4*epsilon", ScalarMag(1e+4*SMT::eps()),
00675         "1e+2*epsilon", ScalarMag(1e+2*SMT::eps()),
00676         as<int>(verbLevel) >= as<int>(Teuchos::VERB_MEDIUM) ? out.get() : 0
00677         )
00678       );
00679 #endif
00680     RCP<Thyra::VectorBase<Scalar> >
00681       x_temp = scaled_x_old_->clone_v();
00682     Thyra::Vt_S(&*x_temp,Scalar(-ST::one()*dt));  // undo the scaling
00683     ds_temp.time = t_old_;
00684     ds_temp.x = x_temp;
00685     ds_temp.xdot = x_dot_old_;
00686     ds_temp.accuracy = ScalarMag(dt);
00687     ds_nodes.push_back(ds_temp);
00688     ds_temp.time = t_;
00689     ds_temp.x = x_;
00690     ds_temp.xdot = x_dot_;
00691     ds_temp.accuracy = ScalarMag(dt);
00692     ds_nodes.push_back(ds_temp);
00693   }
00694   else {
00695     if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00696       *out << "Passing one point to interpolator:  " << t_ << std::endl;
00697     }
00698     DataStore<Scalar> ds_temp;
00699     ds_temp.time = t_;
00700     ds_temp.x = x_;
00701     ds_temp.xdot = x_dot_;
00702     ds_temp.accuracy = ScalarMag(ST::zero());
00703     ds_nodes.push_back(ds_temp);
00704   }
00705   interpolate<Scalar>(*interpolator_,rcp(&ds_nodes,false),time_vec,&ds_out);
00706   Array<Scalar> time_out;
00707   dataStoreVectorToVector(ds_out,&time_out,x_vec,xdot_vec,accuracy_vec);
00708   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00709     *out << "Passing out the interpolated values:" << std::endl;
00710     for (int i=0; i<Teuchos::as<int>(time_out.size()) ; ++i) {
00711       *out << "time[" << i << "] = " << time_out[i] << std::endl;
00712       *out << "x_vec[" << i << "] = " << std::endl;
00713       (*x_vec)[i]->describe(*out,Teuchos::VERB_EXTREME);
00714       if (xdot_vec) {
00715         if ( (*xdot_vec)[i] == Teuchos::null) {
00716           *out << "xdot_vec[" << i << "] = Teuchos::null" << std::endl;
00717         }
00718         else {
00719           *out << "xdot_vec[" << i << "] = " << std::endl;
00720           (*xdot_vec)[i]->describe(*out,Teuchos::VERB_EXTREME);
00721         }
00722       }
00723       if(accuracy_vec) {
00724         *out << "accuracy[" << i << "] = " << (*accuracy_vec)[i] << std::endl;
00725       }
00726     }
00727   }
00728   if ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
00729     *out
00730       << "Leaving " << Teuchos::TypeNameTraits<BackwardEulerStepper<Scalar> >::name()
00731       << "::getPoints(...) ...\n"; 
00732   }
00733   */
00734 
00735 }
00736 
00737 
00738 template<class Scalar>
00739 void BackwardEulerStepper<Scalar>::getNodes(Array<Scalar>* time_vec) const
00740 {
00741   using Teuchos::as;
00742 
00743   TEUCHOS_ASSERT( time_vec != NULL );
00744 
00745   time_vec->clear();
00746 
00747   if (!haveInitialCondition_) {
00748     return;
00749   }
00750 
00751   time_vec->push_back(t_old_);
00752 
00753   if (numSteps_ > 0) {
00754     time_vec->push_back(t_);
00755   }
00756 
00757   RCP<Teuchos::FancyOStream> out = this->getOStream();
00758   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00759   Teuchos::OSTab ostab(out,1,"BES::getNodes");
00760   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00761     *out << this->description() << std::endl;
00762     for (int i=0 ; i<Teuchos::as<int>(time_vec->size()) ; ++i) {
00763       *out << "time_vec[" << i << "] = " << (*time_vec)[i] << std::endl;
00764     }
00765   }
00766 }
00767 
00768 
00769 template<class Scalar>
00770 void BackwardEulerStepper<Scalar>::removeNodes(Array<Scalar>& time_vec) 
00771 {
00772   initialize();
00773   using Teuchos::as;
00774   RCP<Teuchos::FancyOStream> out = this->getOStream();
00775   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00776   Teuchos::OSTab ostab(out,1,"BES::removeNodes");
00777   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00778     *out << "time_vec = " << std::endl;
00779     for (int i=0 ; i<Teuchos::as<int>(time_vec.size()) ; ++i) {
00780       *out << "time_vec[" << i << "] = " << time_vec[i] << std::endl;
00781     }
00782   }
00783   TEST_FOR_EXCEPTION(true,std::logic_error,"Error, removeNodes is not implemented for BackwardEulerStepper at this time.\n");
00784   // TODO:
00785   // if any time in time_vec matches t_ or t_old_, then do the following:
00786   // remove t_old_:  set t_old_ = t_ and set scaled_x_old_ = x_
00787   // remove t_:  set t_ = t_old_ and set x_ = -dt*scaled_x_old_
00788 }
00789 
00790 
00791 template<class Scalar>
00792 int BackwardEulerStepper<Scalar>::getOrder() const
00793 {
00794   return(1);
00795 }
00796 
00797 
00798 // Overridden from Teuchos::ParameterListAcceptor
00799 
00800 
00801 template <class Scalar>
00802 void BackwardEulerStepper<Scalar>::setParameterList(
00803   RCP<Teuchos::ParameterList> const& paramList
00804   )
00805 {
00806   TEST_FOR_EXCEPT(is_null(paramList));
00807   paramList->validateParametersAndSetDefaults(*this->getValidParameters());
00808   parameterList_ = paramList;
00809   Teuchos::readVerboseObjectSublist(&*parameterList_,this);
00810 }
00811 
00812 
00813 template <class Scalar>
00814 RCP<Teuchos::ParameterList>
00815 BackwardEulerStepper<Scalar>::getNonconstParameterList()
00816 {
00817   return(parameterList_);
00818 }
00819 
00820 
00821 template <class Scalar>
00822 RCP<Teuchos::ParameterList>
00823 BackwardEulerStepper<Scalar>::unsetParameterList()
00824 {
00825   RCP<Teuchos::ParameterList>
00826     temp_param_list = parameterList_;
00827   parameterList_ = Teuchos::null;
00828   return(temp_param_list);
00829 }
00830 
00831 
00832 template<class Scalar>
00833 RCP<const Teuchos::ParameterList>
00834 BackwardEulerStepper<Scalar>::getValidParameters() const
00835 {
00836   using Teuchos::ParameterList;
00837   static RCP<const ParameterList> validPL;
00838   if (is_null(validPL)) {
00839     RCP<ParameterList> pl = Teuchos::parameterList();
00840     Teuchos::setupVerboseObjectSublist(&*pl);
00841     validPL = pl;
00842   }
00843   return validPL;
00844 }
00845 
00846 
00847 // Overridden from Teuchos::Describable
00848 
00849 
00850 template<class Scalar>
00851 void BackwardEulerStepper<Scalar>::describe(
00852   Teuchos::FancyOStream &out,
00853   const Teuchos::EVerbosityLevel verbLevel
00854   ) const
00855 {
00856   using Teuchos::as;
00857   Teuchos::OSTab tab(out);
00858   if (!isInitialized_) {
00859     out << this->description() << " : This stepper is not initialized yet" << std::endl;
00860     return;
00861   }
00862   if (
00863     as<int>(verbLevel) == as<int>(Teuchos::VERB_DEFAULT)
00864     ||
00865     as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW)
00866     )
00867   {
00868     out << this->description() << "::describe:" << std::endl;
00869     out << "model = " << model_->description() << std::endl;
00870     out << "solver = " << solver_->description() << std::endl;
00871     if (neModel_ == Teuchos::null) {
00872       out << "neModel = Teuchos::null" << std::endl;
00873     } else {
00874       out << "neModel = " << neModel_->description() << std::endl;
00875     }
00876   }
00877   else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW)) {
00878     out << "t_ = " << t_ << std::endl;
00879     out << "t_old_ = " << t_old_ << std::endl;
00880   }
00881   else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_MEDIUM)) {
00882   }
00883   else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH)) {
00884     out << "model_ = " << std::endl;
00885     model_->describe(out,verbLevel);
00886     out << "solver_ = " << std::endl;
00887     solver_->describe(out,verbLevel);
00888     if (neModel_ == Teuchos::null) {
00889       out << "neModel = Teuchos::null" << std::endl;
00890     } else {
00891       out << "neModel = " << std::endl;
00892       neModel_->describe(out,verbLevel);
00893     }
00894     out << "x_ = " << std::endl;
00895     x_->describe(out,verbLevel);
00896     out << "scaled_x_old_ = " << std::endl;
00897     scaled_x_old_->describe(out,verbLevel);
00898   }
00899 }
00900 
00901 
00902 // private
00903 
00904 
00905 template <class Scalar>
00906 void BackwardEulerStepper<Scalar>::initialize()
00907 {
00908 
00909   if (isInitialized_)
00910     return;
00911 
00912   TEST_FOR_EXCEPT(is_null(model_));
00913   TEST_FOR_EXCEPT(is_null(solver_));
00914   TEST_FOR_EXCEPT(!haveInitialCondition_);
00915 
00916 #ifdef RYTHMOS_DEBUG
00917   THYRA_ASSERT_VEC_SPACES(
00918     "Rythmos::BackwardEulerStepper::initialize(...)",
00919     *x_->space(), *model_->get_x_space() );
00920 #endif // RYTHMOS_DEBUG
00921 
00922   if ( is_null(interpolator_) ) {
00923     // If an interpolator has not been explicitly set, then just create
00924     // a default linear interpolator.
00925     interpolator_ = linearInterpolator<Scalar>();
00926     // 2007/05/18: rabartl: ToDo: Replace this with a Hermete interplator
00927     // when it is implementated!
00928   }
00929 
00930   isInitialized_ = true;
00931 
00932 }
00933 
00934 template<class Scalar>
00935 RCP<const MomentoBase<Scalar> >
00936 BackwardEulerStepper<Scalar>::getMomento() const
00937 {
00938   RCP<BackwardEulerStepperMomento<Scalar> > momento = Teuchos::rcp(new BackwardEulerStepperMomento<Scalar>());
00939   momento->set_scaled_x_old(scaled_x_old_);
00940   momento->set_x_dot_old(x_dot_old_);
00941   momento->set_x(x_);
00942   momento->set_x_dot(x_dot_);
00943   momento->set_t(t_);
00944   momento->set_t_old(t_old_);
00945   momento->set_dt(dt_);
00946   momento->set_numSteps(numSteps_);
00947   momento->set_isInitialized(isInitialized_);
00948   momento->set_haveInitialCondition(haveInitialCondition_);
00949   momento->set_parameterList(parameterList_);
00950   momento->set_basePoint(basePoint_);
00951   momento->set_neModel(neModel_);
00952   momento->set_interpolator(interpolator_);
00953   return momento;
00954 }
00955 
00956 template<class Scalar>
00957 void BackwardEulerStepper<Scalar>::setMomento(
00958     const Ptr<const MomentoBase<Scalar> >& momentoPtr,
00959     const RCP<const Thyra::ModelEvaluator<Scalar> >& model,
00960     const RCP<Thyra::NonlinearSolverBase<Scalar> >& solver
00961     ) 
00962 { 
00963   Ptr<const BackwardEulerStepperMomento<Scalar> > feMomentoPtr = 
00964     Teuchos::ptr_dynamic_cast<const BackwardEulerStepperMomento<Scalar> >(momentoPtr,true);
00965   const BackwardEulerStepperMomento<Scalar>& feMomento = *feMomentoPtr;
00966   model_ = model;
00967   solver_ = solver;
00968   scaled_x_old_ = feMomento.get_scaled_x_old();
00969   x_dot_old_ = feMomento.get_x_dot_old();
00970   x_ = feMomento.get_x();
00971   x_dot_ = feMomento.get_x_dot();
00972   t_ = feMomento.get_t();
00973   t_old_ = feMomento.get_t_old();
00974   dt_ = feMomento.get_dt();
00975   numSteps_ = feMomento.get_numSteps();
00976   isInitialized_ = feMomento.get_isInitialized();
00977   haveInitialCondition_ = feMomento.get_haveInitialCondition();
00978   parameterList_ = feMomento.get_parameterList();
00979   basePoint_ = feMomento.get_basePoint();
00980   neModel_ = feMomento.get_neModel();
00981   interpolator_ = feMomento.get_interpolator();
00982   this->checkConsistentState_();
00983 }
00984 
00985 template<class Scalar>
00986 void BackwardEulerStepper<Scalar>::checkConsistentState_()
00987 {
00988   if (isInitialized_) {
00989     TEUCHOS_ASSERT( !Teuchos::is_null(model_) );
00990     TEUCHOS_ASSERT( !Teuchos::is_null(solver_) );
00991     TEUCHOS_ASSERT( haveInitialCondition_ );
00992     TEUCHOS_ASSERT( !Teuchos::is_null(interpolator_) );
00993   }
00994   if (haveInitialCondition_) {
00995     // basePoint_ should be defined
00996     typedef Teuchos::ScalarTraits<Scalar> ST;
00997     TEUCHOS_ASSERT( !ST::isnaninf(t_) );
00998     TEUCHOS_ASSERT( !ST::isnaninf(t_old_) );
00999     TEUCHOS_ASSERT( !Teuchos::is_null(scaled_x_old_) );
01000     TEUCHOS_ASSERT( !Teuchos::is_null(x_dot_old_) );
01001     TEUCHOS_ASSERT( !Teuchos::is_null(x_) );
01002     TEUCHOS_ASSERT( !Teuchos::is_null(x_dot_) );
01003     TEUCHOS_ASSERT( t_ >= basePoint_.get_t() );
01004     TEUCHOS_ASSERT( t_old_ >= basePoint_.get_t() );
01005   }
01006   if (numSteps_ > 0) {
01007     TEUCHOS_ASSERT(isInitialized_);
01008   }
01009 }
01010 
01011 
01012 // 
01013 // Explicit Instantiation macro
01014 //
01015 // Must be expanded from within the Rythmos namespace!
01016 //
01017 
01018 #define RYTHMOS_BACKWARD_EULER_STEPPER_INSTANT(SCALAR) \
01019   \
01020   template class BackwardEulerStepper< SCALAR >; \
01021   \
01022   template RCP< BackwardEulerStepper< SCALAR > > \
01023   backwardEulerStepper( \
01024     const RCP<const Thyra::ModelEvaluator< SCALAR > > &model, \
01025     const RCP<Thyra::NonlinearSolverBase< SCALAR > > &solver \
01026       ); \
01027   template RCP< BackwardEulerStepper< SCALAR > > \
01028   backwardEulerStepper(); 
01029    
01030 
01031 
01032 
01033 } // namespace Rythmos
01034 
01035 #endif //Rythmos_BACKWARD_EULER_STEPPER_DEF_H

Generated on Wed May 12 21:25:42 2010 for Rythmos - Transient Integration for Differential Equations by  doxygen 1.4.7