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

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