Rythmos - Transient Integration for Differential Equations Version of the Day
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., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
00023 // USA
00024 // Questions? Contact Todd S. Coffey (tscoffe@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 //@HEADER
00028 
00029 #ifndef Rythmos_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<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 HAVE_RYTHMOS_DEBUG
00124   TEUCHOS_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   TEUCHOS_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   TEUCHOS_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 template<class Scalar>
00303 void ThetaStepper<Scalar>::setNonconstModel(
00304   const RCP<Thyra::ModelEvaluator<Scalar> >& model
00305   )
00306 {
00307   this->setModel(model); // TODO 09/09/09 tscoffe:  use ConstNonconstObjectContainer!
00308 }
00309 
00310 
00311 template<class Scalar>
00312 RCP<const Thyra::ModelEvaluator<Scalar> >
00313 ThetaStepper<Scalar>::getModel() const
00314 {
00315   return model_;
00316 }
00317 
00318 
00319 template<class Scalar>
00320 RCP<Thyra::ModelEvaluator<Scalar> >
00321 ThetaStepper<Scalar>::getNonconstModel() 
00322 {
00323   TEUCHOS_TEST_FOR_EXCEPT(true);
00324   return Teuchos::null;
00325 }
00326 
00327 
00328 template<class Scalar>
00329 void ThetaStepper<Scalar>::setInitialCondition(
00330   const Thyra::ModelEvaluatorBase::InArgs<Scalar> &initialCondition
00331   )
00332 {
00333 
00334   typedef Teuchos::ScalarTraits<Scalar> ST;
00335   typedef Thyra::ModelEvaluatorBase MEB;
00336 
00337   basePoint_ = initialCondition;
00338 
00339   // x
00340 
00341   RCP<const Thyra::VectorBase<Scalar> >
00342     x_init = initialCondition.get_x();
00343 
00344 #ifdef HAVE_RYTHMOS_DEBUG
00345   TEUCHOS_TEST_FOR_EXCEPTION(
00346     is_null(x_init), std::logic_error,
00347     "Error, if the client passes in an intial condition to setInitialCondition(...),\n"
00348     "then x can not be null!" );
00349 #endif
00350 
00351   x_ = x_init->clone_v();
00352 
00353   // x_dot
00354 
00355   RCP<const Thyra::VectorBase<Scalar> >
00356     x_dot_init = initialCondition.get_x_dot();
00357 
00358   if (!is_null(x_dot_init)) {
00359     x_dot_ = x_dot_init->clone_v();
00360   }
00361   else {
00362     x_dot_ = createMember(x_->space());
00363     assign(x_dot_.ptr(),ST::zero());
00364   }
00365 
00366   // t
00367   
00368   t_ = initialCondition.get_t();
00369 
00370   t_old_ = t_;
00371 
00372   dt_old_ = 0.0;
00373 
00374   // x pre
00375   
00376   x_pre_ = x_->clone_v();
00377 
00378   // x old
00379 
00380   x_old_ = x_->clone_v();
00381 
00382   // x dot base
00383 
00384   x_dot_base_ = x_->clone_v();
00385 
00386   // x dot old
00387 
00388   x_dot_old_ = x_dot_->clone_v();
00389 
00390   // x dot really old
00391 
00392   x_dot_really_old_ = x_dot_->clone_v();
00393 
00394   haveInitialCondition_ = true;
00395 
00396 }
00397 
00398 
00399 template<class Scalar>
00400 Thyra::ModelEvaluatorBase::InArgs<Scalar> 
00401 ThetaStepper<Scalar>::getInitialCondition() const
00402 {
00403   return basePoint_;
00404 }
00405 
00406 
00407 template<class Scalar>
00408 Scalar ThetaStepper<Scalar>::takeStep(Scalar dt, StepSizeType stepSizeType)
00409 {
00410 
00411   using Teuchos::as;
00412   using Teuchos::incrVerbLevel;
00413   typedef Teuchos::ScalarTraits<Scalar> ST;
00414   typedef Thyra::NonlinearSolverBase<Scalar> NSB;
00415   typedef Teuchos::VerboseObjectTempState<NSB> VOTSNSB;
00416 
00417   initialize_();
00418 
00419   // DEBUG
00420   //this->setOverridingVerbLevel(Teuchos::VERB_EXTREME);
00421   //this->setVerbLevel(Teuchos::VERB_EXTREME);
00422 
00423   RCP<Teuchos::FancyOStream> out = this->getOStream();
00424   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00425   Teuchos::OSTab ostab(out,1,"TS::takeStep");
00426   VOTSNSB solver_outputTempState(solver_,out,incrVerbLevel(verbLevel,-1));
00427 
00428   if ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
00429     *out
00430       << "\nEntering " << Teuchos::TypeNameTraits<ThetaStepper<Scalar> >::name()
00431       << "::takeStep("<<dt<<","<<toString(stepSizeType)<<") ...\n"; 
00432   }
00433 
00434   dt_ = dt;
00435   V_StV( x_old_.ptr(), Scalar(ST::one()), *x_ );
00436   V_StV( x_dot_really_old_.ptr(), Scalar(ST::one()), *x_dot_old_ );
00437   V_StV( x_dot_old_.ptr(),        Scalar(ST::one()), *x_dot_ );
00438 
00439   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_EXTREME) ) {
00440     *out << "\nSetting dt_ and old data ..." << std::endl;
00441     *out << "\ndt_ = " << dt_;
00442     *out << "\nx_old_ = " << *x_old_;
00443     *out << "\nx_dot_old_ = " << *x_dot_old_;
00444     *out << "\nx_dot_really_old_ = " << *x_dot_really_old_;
00445   }
00446 
00447   if ((stepSizeType == STEP_TYPE_VARIABLE) || (dt == ST::zero())) {
00448     if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) )
00449       *out << "\nThe arguments to takeStep are not valid for ThetaStepper at this time." << std::endl;
00450     return(Scalar(-ST::one()));
00451   }
00452   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00453     *out << "\ndt = " << dt << std::endl;
00454     *out << "\nstepSizeType = " << stepSizeType << std::endl;
00455   }
00456 
00457   // compute predictor
00458   obtainPredictor_();
00459 
00460   //
00461   // Setup the nonlinear equations:
00462   //
00463   //   substitute:
00464   // 
00465   //   x_dot = ( 1/(theta*dt) )*x + ( -1/(theta*dt) )*x_old + ( -(1-theta)/theta )*x_dot_old
00466   //   
00467   //   f( x_dot, x, t ) = 0
00468   //
00469 
00470   const double theta = (numSteps_+1>=predictor_corrector_begin_after_step_) ? theta_ : 1.0;
00471 
00472   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_EXTREME) ) {
00473     *out << "\nSetting x_dot_base_ ..." << std::endl;
00474     *out << "\ntheta = " << theta;
00475     *out << "\nx_ = " << *x_;
00476     *out << "\nx_dot_old_ = " << *x_dot_old_;
00477   }
00478 
00479   const Scalar coeff_x_dot = Scalar(ST::one()/(theta*dt)); 
00480   const Scalar coeff_x = ST::one();
00481 
00482   const Scalar x_coeff = Scalar(-coeff_x_dot);
00483   const Scalar x_dot_old_coeff = Scalar( -(ST::one()-theta)/theta);
00484 
00485   V_StV( x_dot_base_.ptr(), x_coeff, *x_old_ );
00486   Vp_StV( x_dot_base_.ptr(), x_dot_old_coeff, *x_dot_old_);
00487 
00488   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_EXTREME) ) {
00489     *out << "\nx_dot_base_ = " << *x_dot_base_;
00490   }
00491 
00492   if(!neModel_.get()) {
00493     neModel_ = Teuchos::rcp(new Rythmos::SingleResidualModelEvaluator<Scalar>());
00494   }
00495 
00496   neModel_->initializeSingleResidualModel(
00497     model_, 
00498     basePoint_,
00499     coeff_x_dot,
00500     x_dot_base_,
00501     coeff_x,
00502     Teuchos::null, // x_base
00503     t_+dt, // t_base
00504     Teuchos::null // x_bar_init
00505     );
00506   if( solver_->getModel().get() != neModel_.get() ) {
00507     solver_->setModel(neModel_);
00508   }
00509   
00510   solver_->setVerbLevel(this->getVerbLevel());
00511 
00512   //
00513   // Solve the implicit nonlinear system to a tolerance of ???
00514   //
00515   
00516   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
00517     *out << "\nSolving the implicit theta-stepper timestep equation"
00518    << " with theta = " << theta << "\n";
00519   }
00520 
00521   Thyra::SolveStatus<Scalar>
00522     neSolveStatus = solver_->solve(&*x_);
00523 
00524   // In the above solve, on input *x_ is the old value of x for the previous
00525   // time step which is used as the initial guess for the solver.  On output,
00526   // *x_ is the converged timestep solution.
00527  
00528   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
00529     *out << "\nOutput status of nonlinear solve:\n" << neSolveStatus;
00530   }
00531 
00532   // 2007/05/18: rabartl: ToDo: Above, get the solve status from the above
00533   // solve and at least print warning message if the solve fails!  Actually,
00534   // you should most likely thrown an exception if the solve fails or return
00535   // false if appropriate
00536 
00537   //
00538   // Update the step data
00539   //
00540 
00541   // x_dot = ( 1/(theta*dt) )*x + ( -1/(theta*dt) )*x_old + ( -(1-theta)/theta )*x_dot_old
00542 
00543   V_StV(  x_dot_.ptr(), Scalar( ST::one()/(theta*dt)), *x_ );
00544   Vp_StV( x_dot_.ptr(), Scalar(-ST::one()/(theta*dt)), *x_old_ );
00545   Vp_StV( x_dot_.ptr(), Scalar( -(ST::one()-theta)/theta), *x_dot_old_ );
00546 
00547   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_EXTREME) ) {
00548     *out << "\nUpdating x_dot_ ...\n";
00549     *out << "\nx_dot_ = " << *x_dot_ << std::endl;
00550   }
00551 
00552   t_old_ = t_;
00553   dt_old_ = dt;
00554   t_ += dt;
00555 
00556   numSteps_++;
00557 
00558   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00559     *out << "\nt_old_ = " << t_old_ << std::endl;
00560     *out << "\nt_ = " << t_ << std::endl;
00561   }
00562 
00563 #ifdef HAVE_RYTHMOS_DEBUG
00564 
00565   if ( includesVerbLevel(verbLevel,Teuchos::VERB_LOW) )
00566     *out << "\nChecking to make sure that solution and the interpolated solution are the same! ...\n";
00567 
00568   {
00569 
00570     typedef ScalarTraits<Scalar> ST;
00571     typedef typename ST::magnitudeType ScalarMag;
00572     typedef ScalarTraits<ScalarMag> SMT;
00573     
00574     Teuchos::OSTab tab(out);
00575 
00576     const StepStatus<Scalar> stepStatus = this->getStepStatus();
00577 
00578     RCP<const Thyra::VectorBase<Scalar> >
00579       x = stepStatus.solution,
00580       xdot = stepStatus.solutionDot;
00581 
00582     Array<Scalar> time_vec = Teuchos::tuple(stepStatus.time);
00583     Array<RCP<const Thyra::VectorBase<Scalar> > > x_vec, xdot_vec;
00584     this->getPoints(time_vec,&x_vec,&xdot_vec,0);
00585 
00586     RCP<const Thyra::VectorBase<Scalar> >
00587       x_interp = x_vec[0],
00588       xdot_interp = xdot_vec[0];
00589 
00590     TEUCHOS_TEST_FOR_EXCEPT(
00591       !Thyra::testRelNormDiffErr(
00592         "x", *x, "x_interp", *x_interp,
00593         "2*epsilon", ScalarMag(100.0*SMT::eps()),
00594         "2*epsilon", ScalarMag(100.0*SMT::eps()),
00595         includesVerbLevel(verbLevel,Teuchos::VERB_HIGH) ? out.get() : 0
00596         )
00597       );
00598 
00599     TEUCHOS_TEST_FOR_EXCEPT(
00600       !Thyra::testRelNormDiffErr(
00601         "xdot", *xdot, "xdot_interp", *xdot_interp,
00602         "2*epsilon", ScalarMag(100.0*SMT::eps()),
00603         "2*epsilon", ScalarMag(100.0*SMT::eps()),
00604         includesVerbLevel(verbLevel,Teuchos::VERB_HIGH) ? out.get() : 0
00605         )
00606       );
00607 
00608   }
00609 
00610   // 2007/07/25: rabartl: ToDo: Move the above test into a helper function so
00611   // that it can be used from lots of different places!
00612 
00613 #endif // HAVE_RYTHMOS_DEBUG
00614 
00615   if ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
00616     *out
00617       << "\nLeaving " << Teuchos::TypeNameTraits<ThetaStepper<Scalar> >::name()
00618       << "::takeStep(...) ...\n"; 
00619   }
00620 
00621   return(dt);
00622 
00623 }
00624 
00625 
00626 template<class Scalar>
00627 const StepStatus<Scalar> ThetaStepper<Scalar>::getStepStatus() const
00628 {
00629 
00630   typedef Teuchos::ScalarTraits<Scalar> ST;
00631 
00632   StepStatus<Scalar> stepStatus; // Defaults to unknown status
00633 
00634   if (!isInitialized_) {
00635     stepStatus.stepStatus = STEP_STATUS_UNINITIALIZED;
00636   }
00637   else if (numSteps_ > 0) {
00638     stepStatus.stepStatus = STEP_STATUS_CONVERGED; 
00639   }
00640   // else unknown
00641 
00642   stepStatus.stepSize = dt_;
00643   stepStatus.order = 1;
00644   stepStatus.time = t_;
00645   stepStatus.solution = x_;
00646   stepStatus.solutionDot = x_dot_;
00647 
00648   return(stepStatus);
00649 
00650 }
00651 
00652 
00653 // Overridden from InterpolationBufferBase
00654 
00655 
00656 template<class Scalar>
00657 RCP<const Thyra::VectorSpaceBase<Scalar> >
00658 ThetaStepper<Scalar>::get_x_space() const
00659 {
00660   return ( !is_null(model_) ? model_->get_x_space() : Teuchos::null );
00661 }
00662 
00663 
00664 template<class Scalar>
00665 void ThetaStepper<Scalar>::addPoints(
00666   const Array<Scalar>& time_vec,
00667   const Array<RCP<const Thyra::VectorBase<Scalar> > >& x_vec,
00668   const Array<RCP<const Thyra::VectorBase<Scalar> > >& xdot_vec
00669   )
00670 {
00671 
00672   typedef Teuchos::ScalarTraits<Scalar> ST;
00673   using Teuchos::as;
00674 
00675   initialize_();
00676 
00677 #ifdef HAVE_RYTHMOS_DEBUG
00678   TEUCHOS_TEST_FOR_EXCEPTION(
00679     time_vec.size() == 0, std::logic_error,
00680     "Error, addPoints called with an empty time_vec array!\n");
00681 #endif // HAVE_RYTHMOS_DEBUG
00682 
00683   RCP<Teuchos::FancyOStream> out = this->getOStream();
00684   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00685   Teuchos::OSTab ostab(out,1,"TS::setPoints");
00686 
00687   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00688     *out << "time_vec = " << std::endl;
00689     for (int i=0 ; i<Teuchos::as<int>(time_vec.size()) ; ++i) {
00690       *out << "time_vec[" << i << "] = " << time_vec[i] << std::endl;
00691     }
00692   }
00693   else if (time_vec.size() == 1) {
00694     int n = 0;
00695     t_ = time_vec[n];
00696     t_old_ = t_;
00697     Thyra::V_V(x_.ptr(),*x_vec[n]);
00698     Thyra::V_V(x_dot_base_.ptr(),*x_);
00699   }
00700   else {
00701     int n = time_vec.size()-1;
00702     int nm1 = time_vec.size()-2;
00703     t_ = time_vec[n];
00704     t_old_ = time_vec[nm1];
00705     Thyra::V_V(x_.ptr(),*x_vec[n]);
00706     Scalar dt = t_ - t_old_;
00707     Thyra::V_StV(x_dot_base_.ptr(),Scalar(-ST::one()/dt),*x_vec[nm1]);
00708   }
00709 }
00710 
00711 
00712 template<class Scalar>
00713 TimeRange<Scalar> ThetaStepper<Scalar>::getTimeRange() const
00714 {
00715   if ( !isInitialized_ && haveInitialCondition_ )
00716     return timeRange<Scalar>(t_,t_);
00717   if ( !isInitialized_ && !haveInitialCondition_ )
00718     return invalidTimeRange<Scalar>();
00719   return timeRange<Scalar>(t_old_,t_);
00720 }
00721 
00722 
00723 template<class Scalar>
00724 void ThetaStepper<Scalar>::getPoints(
00725   const Array<Scalar>& time_vec,
00726   Array<RCP<const Thyra::VectorBase<Scalar> > >* x_vec,
00727   Array<RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec,
00728   Array<ScalarMag>* accuracy_vec
00729   ) const
00730 {
00731   using Teuchos::constOptInArg;
00732   using Teuchos::ptr;
00733   typedef Teuchos::ScalarTraits<Scalar> ST;
00734 
00735   TEUCHOS_ASSERT(haveInitialCondition_);
00736 
00737   RCP<Thyra::VectorBase<Scalar> > x_temp = x_;
00738   if (compareTimeValues(t_old_,t_)!=0) {
00739     Scalar dt = t_ - t_old_;
00740     x_temp = x_dot_base_->clone_v();
00741     Thyra::Vt_S(x_temp.ptr(),Scalar(-ST::one()*dt));  // undo the scaling
00742   }
00743   defaultGetPoints<Scalar>(
00744       t_old_, constOptInArg(*x_temp), constOptInArg(*x_dot_old_),
00745       t_, constOptInArg(*x_), constOptInArg(*x_dot_),
00746       time_vec, ptr(x_vec), ptr(xdot_vec), ptr(accuracy_vec),
00747       ptr(interpolator_.get())
00748       );
00749 
00750   /*
00751   using Teuchos::as;
00752   typedef Teuchos::ScalarTraits<Scalar> ST;
00753   typedef typename ST::magnitudeType ScalarMag;
00754   typedef Teuchos::ScalarTraits<ScalarMag> SMT;
00755   typename DataStore<Scalar>::DataStoreVector_t ds_nodes;
00756   typename DataStore<Scalar>::DataStoreVector_t ds_out;
00757 
00758 #ifdef HAVE_RYTHMOS_DEBUG
00759   TEUCHOS_TEST_FOR_EXCEPT(!haveInitialCondition_);
00760   TEUCHOS_TEST_FOR_EXCEPT( 0 == x_vec );
00761 #endif
00762 
00763   RCP<Teuchos::FancyOStream> out = this->getOStream();
00764   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00765   Teuchos::OSTab ostab(out,1,"TS::getPoints");
00766   if ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
00767     *out
00768       << "\nEntering " << Teuchos::TypeNameTraits<ThetaStepper<Scalar> >::name()
00769       << "::getPoints(...) ...\n"; 
00770   }
00771   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00772     for (int i=0 ; i<Teuchos::as<int>(time_vec.size()) ; ++i) {
00773       *out << "time_vec[" << i << "] = " << time_vec[i] << std::endl;
00774     }
00775     *out << "I can interpolate in the interval [" << t_old_ << "," << t_ << "]." << std::endl;
00776   }
00777 
00778   if (t_old_ != t_) {
00779     if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00780       *out << "Passing two points to interpolator:  " << t_old_ << " and " << t_ << std::endl;
00781     }
00782     DataStore<Scalar> ds_temp;
00783     Scalar dt = t_ - t_old_;
00784 #ifdef HAVE_RYTHMOS_DEBUG
00785     TEUCHOS_TEST_FOR_EXCEPT(
00786       !Thyra::testRelErr(
00787         "dt", dt, "dt_", dt_,
00788         "1e+4*epsilon", ScalarMag(1e+4*SMT::eps()),
00789         "1e+2*epsilon", ScalarMag(1e+2*SMT::eps()),
00790         as<int>(verbLevel) >= as<int>(Teuchos::VERB_MEDIUM) ? out.get() : 0
00791         )
00792       );
00793 #endif
00794     ds_temp.time = t_old_;
00795     ds_temp.x = x_old_;
00796     ds_temp.xdot = x_dot_old_;
00797     ds_temp.accuracy = ScalarMag(dt);
00798     ds_nodes.push_back(ds_temp);
00799     ds_temp.time = t_;
00800     ds_temp.x = x_;
00801     ds_temp.xdot = x_dot_;
00802     ds_temp.accuracy = ScalarMag(dt);
00803     ds_nodes.push_back(ds_temp);
00804   }
00805   else {
00806     if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00807       *out << "Passing one point to interpolator:  " << t_ << std::endl;
00808     }
00809     DataStore<Scalar> ds_temp;
00810     ds_temp.time = t_;
00811     ds_temp.x = x_;
00812     ds_temp.xdot = x_dot_;
00813     ds_temp.accuracy = ScalarMag(ST::zero());
00814     ds_nodes.push_back(ds_temp);
00815   }
00816   interpolate<Scalar>(*interpolator_,rcp(&ds_nodes,false),time_vec,&ds_out);
00817   Array<Scalar> time_out;
00818   dataStoreVectorToVector(ds_out,&time_out,x_vec,xdot_vec,accuracy_vec);
00819   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00820     *out << "Passing out the interpolated values:" << std::endl;
00821     for (int i=0; i<Teuchos::as<int>(time_out.size()) ; ++i) {
00822       if (x_vec) {
00823         if ( (*x_vec)[i] == Teuchos::null) {
00824           *out << "x_vec[" << i << "] = Teuchos::null" << std::endl;
00825   }
00826   else {
00827     *out << "time[" << i << "] = " << time_out[i] << std::endl;
00828     *out << "x_vec[" << i << "] = " << std::endl;
00829     (*x_vec)[i]->describe(*out,Teuchos::VERB_EXTREME);
00830   }
00831       }
00832       if (xdot_vec) {
00833         if ( (*xdot_vec)[i] == Teuchos::null) {
00834           *out << "xdot_vec[" << i << "] = Teuchos::null" << std::endl;
00835         }
00836         else {
00837           *out << "xdot_vec[" << i << "] = " << std::endl;
00838           (*xdot_vec)[i]->describe(*out,Teuchos::VERB_EXTREME);
00839         }
00840       }
00841       if(accuracy_vec) {
00842         *out << "accuracy[" << i << "] = " << (*accuracy_vec)[i] << std::endl;
00843       }
00844     }
00845   }
00846   if ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
00847     *out
00848       << "Leaving " << Teuchos::TypeNameTraits<ThetaStepper<Scalar> >::name()
00849       << "::getPoints(...) ...\n"; 
00850   }
00851   */
00852 
00853 }
00854 
00855 
00856 template<class Scalar>
00857 void ThetaStepper<Scalar>::getNodes(Array<Scalar>* time_vec) const
00858 {
00859   using Teuchos::as;
00860 
00861   TEUCHOS_ASSERT( time_vec != NULL );
00862 
00863   time_vec->clear();
00864   if (!haveInitialCondition_) {
00865     return;
00866   }
00867 
00868   time_vec->push_back(t_old_);
00869   if (numSteps_ > 0) {
00870     time_vec->push_back(t_);
00871   }
00872 
00873   RCP<Teuchos::FancyOStream> out = this->getOStream();
00874   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00875   Teuchos::OSTab ostab(out,1,"TS::getNodes");
00876   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00877     *out << this->description() << std::endl;
00878     for (int i=0 ; i<Teuchos::as<int>(time_vec->size()) ; ++i) {
00879       *out << "time_vec[" << i << "] = " << (*time_vec)[i] << std::endl;
00880     }
00881   }
00882 }
00883 
00884 
00885 template<class Scalar>
00886 void ThetaStepper<Scalar>::removeNodes(Array<Scalar>& time_vec) 
00887 {
00888   initialize_();
00889   using Teuchos::as;
00890   RCP<Teuchos::FancyOStream> out = this->getOStream();
00891   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00892   Teuchos::OSTab ostab(out,1,"TS::removeNodes");
00893   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00894     *out << "time_vec = " << std::endl;
00895     for (int i=0 ; i<Teuchos::as<int>(time_vec.size()) ; ++i) {
00896       *out << "time_vec[" << i << "] = " << time_vec[i] << std::endl;
00897     }
00898   }
00899   TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Error, removeNodes is not implemented for ThetaStepper at this time.\n");
00900   // TODO:
00901   // if any time in time_vec matches t_ or t_old_, then do the following:
00902   // remove t_old_:  set t_old_ = t_ and set x_dot_base_ = x_
00903   // remove t_:  set t_ = t_old_ and set x_ = -dt*x_dot_base_
00904 }
00905 
00906 
00907 template<class Scalar>
00908 int ThetaStepper<Scalar>::getOrder() const
00909 {
00910   return (thetaStepperType_==ImplicitEuler) ? 1 : 2;
00911 }
00912 
00913 
00914 // Overridden from Teuchos::ParameterListAcceptor
00915 
00916 
00917 template <class Scalar>
00918 void ThetaStepper<Scalar>::setParameterList(
00919   RCP<Teuchos::ParameterList> const& paramList
00920   )
00921 {
00922   TEUCHOS_TEST_FOR_EXCEPT(is_null(paramList));
00923   paramList->validateParametersAndSetDefaults(*this->getValidParameters());
00924   parameterList_ = paramList;
00925   Teuchos::readVerboseObjectSublist(&*parameterList_,this);
00926 
00927   RCP<ParameterList> pl_theta = Teuchos::sublist(parameterList_, RythmosStepControlSettings_name);
00928 
00929   std::string thetaStepperTypeString = 
00930     Teuchos::getParameter<std::string>(*pl_theta, ThetaStepperType_name);
00931   
00932   if (thetaStepperTypeString == "Implicit Euler")
00933     thetaStepperType_ = ImplicitEuler;
00934   else if (thetaStepperTypeString == "Trapezoid")
00935     thetaStepperType_ = Trapezoid;
00936   else
00937     TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, 
00938            "Value of " << ThetaStepperType_name << " = " << thetaStepperTypeString 
00939            << " is invalid for Rythmos::ThetaStepper");
00940 
00941   default_predictor_order_ = 
00942     Teuchos::getParameter<int>(*pl_theta, PredictorOrder_name);
00943   
00944   Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
00945   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00946   if ( Teuchos::as<int>(verbLevel) >= Teuchos::as<int>(Teuchos::VERB_HIGH) ) {
00947     *out << ThetaStepperType_name << " = " << thetaStepperTypeString << std::endl;
00948   }
00949 }
00950 
00951 
00952 template <class Scalar>
00953 RCP<Teuchos::ParameterList>
00954 ThetaStepper<Scalar>::getNonconstParameterList()
00955 {
00956   return(parameterList_);
00957 }
00958 
00959 
00960 template <class Scalar>
00961 RCP<Teuchos::ParameterList>
00962 ThetaStepper<Scalar>::unsetParameterList()
00963 {
00964   RCP<Teuchos::ParameterList>
00965     temp_param_list = parameterList_;
00966   parameterList_ = Teuchos::null;
00967   return(temp_param_list);
00968 }
00969 
00970 
00971 template<class Scalar>
00972 RCP<const Teuchos::ParameterList>
00973 ThetaStepper<Scalar>::getValidParameters() const
00974 {
00975   using Teuchos::ParameterList;
00976 
00977   static RCP<const ParameterList> validPL;
00978 
00979   if (is_null(validPL)) {
00980 
00981     RCP<ParameterList> pl_top_level = Teuchos::parameterList();
00982 
00983     RCP<ParameterList> pl = Teuchos::sublist(pl_top_level, RythmosStepControlSettings_name);
00984 
00985     pl->set<std::string> ( ThetaStepperType_name, ThetaStepperType_default,
00986         "Name of Stepper Type in Theta Stepper"
00987         );
00988 
00989     pl->set<int> ( PredictorOrder_name, PredictorOrder_default,
00990         "Order of Predictor in Theta Stepper, can be 0, 1, 2"
00991         );
00992 
00993     Teuchos::setupVerboseObjectSublist(&*pl_top_level);
00994     validPL = pl_top_level;
00995   }
00996   return validPL;
00997 }
00998 
00999 
01000 // Overridden from Teuchos::Describable
01001 
01002 
01003 template<class Scalar>
01004 void ThetaStepper<Scalar>::describe(
01005   Teuchos::FancyOStream &out,
01006   const Teuchos::EVerbosityLevel verbLevel
01007   ) const
01008 {
01009   using Teuchos::as;
01010   Teuchos::OSTab tab(out);
01011   if (!isInitialized_) {
01012     out << this->description() << " : This stepper is not initialized yet" << std::endl;
01013     return;
01014   }
01015   if (
01016     as<int>(verbLevel) == as<int>(Teuchos::VERB_DEFAULT)
01017     ||
01018     as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW)
01019     )
01020   {
01021     out << this->description() << "::describe:" << std::endl;
01022     out << "model = " << model_->description() << std::endl;
01023     out << "solver = " << solver_->description() << std::endl;
01024     if (neModel_ == Teuchos::null) {
01025       out << "neModel = Teuchos::null" << std::endl;
01026     } else {
01027       out << "neModel = " << neModel_->description() << std::endl;
01028     }
01029   }
01030   else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW)) {
01031     out << "t_ = " << t_ << std::endl;
01032     out << "t_old_ = " << t_old_ << std::endl;
01033   }
01034   else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_MEDIUM)) {
01035   }
01036   else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH)) {
01037     out << "model_ = " << std::endl;
01038     model_->describe(out,verbLevel);
01039     out << "solver_ = " << std::endl;
01040     solver_->describe(out,verbLevel);
01041     if (neModel_ == Teuchos::null) {
01042       out << "neModel = Teuchos::null" << std::endl;
01043     } else {
01044       out << "neModel = " << std::endl;
01045       neModel_->describe(out,verbLevel);
01046     }
01047     out << "x_ = " << std::endl;
01048     x_->describe(out,verbLevel);
01049     out << "x_dot_base_ = " << std::endl;
01050     x_dot_base_->describe(out,verbLevel);
01051   }
01052 }
01053 
01054 
01055 // private
01056 
01057 
01058 template <class Scalar>
01059 void ThetaStepper<Scalar>::initialize_()
01060 {
01061 
01062   typedef Teuchos::ScalarTraits<Scalar> ST;
01063 
01064   if (isInitialized_)
01065     return;
01066 
01067   TEUCHOS_TEST_FOR_EXCEPT(is_null(model_));
01068   TEUCHOS_TEST_FOR_EXCEPT(is_null(solver_));
01069   TEUCHOS_TEST_FOR_EXCEPT(!haveInitialCondition_);
01070 
01071 #ifdef HAVE_RYTHMOS_DEBUG
01072   THYRA_ASSERT_VEC_SPACES(
01073     "Rythmos::ThetaStepper::setInitialCondition(...)",
01074     *x_->space(), *model_->get_x_space() );
01075 #endif // HAVE_RYTHMOS_DEBUG
01076 
01077   if ( is_null(interpolator_) ) {
01078     // If an interpolator has not been explicitly set, then just create
01079     // a default linear interpolator.
01080     interpolator_ = Teuchos::rcp(new LinearInterpolator<Scalar> );
01081     // 2007/05/18: rabartl: ToDo: Replace this with a Hermete interplator
01082     // when it is implementated!
01083   }
01084 
01085   if (thetaStepperType_ == ImplicitEuler)
01086   {
01087     theta_ = 1.0;
01088     predictor_corrector_begin_after_step_ = 2;
01089   }
01090   else
01091   {
01092     theta_ = 0.5;
01093     predictor_corrector_begin_after_step_ = 3;
01094   }
01095 
01096   isInitialized_ = true;
01097 
01098 }
01099 
01100 template<class Scalar>
01101 void ThetaStepper<Scalar>::obtainPredictor_()
01102 {
01103 
01104   using Teuchos::as;
01105   typedef Teuchos::ScalarTraits<Scalar> ST;
01106 
01107   if (!isInitialized_) {
01108     return;
01109   }
01110 
01111   RCP<Teuchos::FancyOStream> out = this->getOStream();
01112   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
01113   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
01114     *out << "Obtaining predictor..." << std::endl;
01115   }
01116 
01117   const int preferred_predictor_order = std::min(default_predictor_order_, thetaStepperType_ + 1);
01118   const int max_predictor_order_at_this_timestep = std::max(0, numSteps_);
01119 
01120   const int predictor_order = std::min(preferred_predictor_order, max_predictor_order_at_this_timestep);
01121 
01122   switch (predictor_order) 
01123   {
01124     case 0:
01125       V_StV(x_pre_.ptr(), Scalar(ST::one()), *x_old_);
01126       break;
01127     case 1:
01128     {
01129       V_StV(x_pre_.ptr(), Scalar(ST::one()), *x_old_);
01130 
01131       TEUCHOS_TEST_FOR_EXCEPT (dt_ <= 0.0);
01132 
01133       Vp_StV(x_pre_.ptr(), dt_, *x_dot_old_);
01134     }
01135     break;
01136     case 2:
01137     {
01138       V_StV(x_pre_.ptr(), Scalar(ST::one()), *x_old_);
01139 
01140       TEUCHOS_TEST_FOR_EXCEPT (dt_ <= 0.0);
01141       TEUCHOS_TEST_FOR_EXCEPT (dt_old_ <= 0.0);
01142 
01143       const Scalar coeff_x_dot_old = (0.5 * dt_) * (2.0 + dt_/dt_old_);
01144       const Scalar coeff_x_dot_really_old = - (0.5 * dt_) * (dt_/dt_old_);
01145       
01146       if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
01147   *out << "x_dot_old_ = " << *x_dot_old_ << std::endl;
01148   *out << "x_dot_really_old_ = " << *x_dot_really_old_ << std::endl;
01149       }
01150       
01151       Vp_StV( x_pre_.ptr(), coeff_x_dot_old, *x_dot_old_);
01152       Vp_StV( x_pre_.ptr(), coeff_x_dot_really_old, *x_dot_really_old_);
01153     }
01154     break;
01155     default:
01156       TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
01157        "Invalid predictor order " << predictor_order << ". Valid values are 0, 1, and 2.");
01158   }
01159   
01160   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
01161     *out << "x_pre_ = " << *x_pre_ << std::endl;
01162   }
01163 
01164   // copy to current solution
01165   V_StV(x_.ptr(), Scalar(ST::one()), *x_pre_);
01166 }
01167 
01168 // 
01169 // Explicit Instantiation macro
01170 //
01171 // Must be expanded from within the Rythmos namespace!
01172 //
01173 
01174 #define RYTHMOS_THETA_STEPPER_INSTANT(SCALAR) \
01175   \
01176   template class ThetaStepper< SCALAR >; \
01177   \
01178   template RCP< ThetaStepper< SCALAR > > \
01179   thetaStepper( \
01180     const RCP<Thyra::ModelEvaluator< SCALAR > >& model, \
01181     const RCP<Thyra::NonlinearSolverBase< SCALAR > >& solver, \
01182     RCP<Teuchos::ParameterList>& parameterList \
01183       );
01184    
01185 
01186 } // namespace Rythmos
01187 
01188 #endif // HAVE_RYTHMOS_EXPERIMENTAL
01189 
01190 #endif //Rythmos_THETA_STEPPER_DEF_H
 All Classes Functions Variables Typedefs Friends