Rythmos - Transient Integration for Differential Equations Version of the Day
Rythmos_ImplicitBDFStepper_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_IMPLICITBDF_STEPPER_DEF_H
00030 #define Rythmos_IMPLICITBDF_STEPPER_DEF_H
00031 
00032 #include "Rythmos_ImplicitBDFStepper_decl.hpp"
00033 #include "Rythmos_StepperHelpers.hpp"
00034 #include "Rythmos_ImplicitBDFStepperStepControl.hpp"
00035 
00036 namespace Rythmos {
00037 
00038 // ////////////////////////////
00039 // Defintions
00040 
00041 // Nonmember constructor
00042 template<class Scalar>
00043 RCP<ImplicitBDFStepper<Scalar> > implicitBDFStepper() {
00044   RCP<ImplicitBDFStepper<Scalar> >
00045     stepper = rcp(new ImplicitBDFStepper<Scalar>() );
00046   return stepper;
00047 }
00048 
00049 template<class Scalar>
00050 RCP<ImplicitBDFStepper<Scalar> > implicitBDFStepper(
00051   const RCP<Thyra::ModelEvaluator<Scalar> >& model,
00052   const RCP<Thyra::NonlinearSolverBase<Scalar> >& solver
00053   )
00054 {
00055   RCP<ImplicitBDFStepper<Scalar> >
00056     stepper = Teuchos::rcp(new ImplicitBDFStepper<Scalar>(model, solver));
00057   return stepper;
00058 }
00059 
00060 template<class Scalar>
00061 RCP<ImplicitBDFStepper<Scalar> > implicitBDFStepper(
00062   const RCP<Thyra::ModelEvaluator<Scalar> >& model,
00063   const RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
00064   const RCP<Teuchos::ParameterList>& parameterList
00065   )
00066 {
00067   RCP<ImplicitBDFStepper<Scalar> >
00068     stepper = Teuchos::rcp(new ImplicitBDFStepper<Scalar>(model,
00069                                                           solver,
00070                                                           parameterList));
00071   return stepper;
00072 }
00073 
00074 // Constructors, intializers, Misc.
00075 
00076 
00077 template<class Scalar>
00078 ImplicitBDFStepper<Scalar>::ImplicitBDFStepper()
00079 {
00080   this->defaultInitializeAll_();
00081   haveInitialCondition_ = false;
00082   isInitialized_=false;
00083 }
00084 
00085 
00086 template<class Scalar>
00087 ImplicitBDFStepper<Scalar>::ImplicitBDFStepper(
00088   const RCP<Thyra::ModelEvaluator<Scalar> >& model
00089   ,const RCP<Thyra::NonlinearSolverBase<Scalar> >& solver
00090   ,const RCP<Teuchos::ParameterList>& parameterList
00091   )
00092 {
00093   this->defaultInitializeAll_();
00094   this->setParameterList(parameterList);
00095   // Now we instantiate the model and the solver
00096   setModel(model);
00097   setSolver(solver);
00098   haveInitialCondition_ = false;
00099   isInitialized_=false;
00100 }
00101 
00102 
00103 template<class Scalar>
00104 ImplicitBDFStepper<Scalar>::ImplicitBDFStepper(
00105   const RCP<Thyra::ModelEvaluator<Scalar> >& model
00106   ,const RCP<Thyra::NonlinearSolverBase<Scalar> >& solver
00107   )
00108 {
00109   this->defaultInitializeAll_();
00110   // Now we instantiate the model and the solver
00111   setModel(model);
00112   setSolver(solver);
00113   haveInitialCondition_ = false;
00114   isInitialized_=false;
00115 }
00116 
00117 template<class Scalar>
00118 const Thyra::VectorBase<Scalar>& 
00119   ImplicitBDFStepper<Scalar>::getxHistory(int index) const
00120 {
00121   TEUCHOS_TEST_FOR_EXCEPTION(!isInitialized_,std::logic_error,
00122       "Error, attempting to call getxHistory before initialization!\n");
00123   TEUCHOS_TEST_FOR_EXCEPT( !(( 0 <= index ) && ( index <= maxOrder_ )) );
00124   TEUCHOS_TEST_FOR_EXCEPT( !( index <= usedOrder_+1 ) );
00125   return(*(xHistory_[index]));
00126 }
00127 
00128 template<class Scalar>
00129 void ImplicitBDFStepper<Scalar>::setStepControlStrategy(const RCP<StepControlStrategyBase<Scalar> >& stepControl)
00130 {
00131   TEUCHOS_TEST_FOR_EXCEPTION(stepControl == Teuchos::null,std::logic_error,"Error, stepControl == Teuchos::null!\n");
00132   stepControl_ = stepControl;    
00133 }
00134 
00135 template<class Scalar>
00136 RCP<StepControlStrategyBase<Scalar> > ImplicitBDFStepper<Scalar>::getNonconstStepControlStrategy() 
00137 {
00138   return(stepControl_);
00139 }
00140 
00141 template<class Scalar>
00142 RCP<const StepControlStrategyBase<Scalar> > ImplicitBDFStepper<Scalar>::getStepControlStrategy() const
00143 {
00144   return(stepControl_);
00145 }
00146 
00147 
00148 // Overridden from SolverAcceptingStepperBase
00149 
00150 
00151 template<class Scalar>
00152 void ImplicitBDFStepper<Scalar>::setSolver(const RCP<Thyra::NonlinearSolverBase<Scalar> > &solver)
00153 {
00154   TEUCHOS_TEST_FOR_EXCEPT(solver == Teuchos::null)
00155     solver_ = solver;
00156 }
00157 
00158 
00159 template<class Scalar>
00160 RCP<Thyra::NonlinearSolverBase<Scalar> >
00161 ImplicitBDFStepper<Scalar>::getNonconstSolver()
00162 {
00163   return (solver_);
00164 }
00165 
00166 
00167 template<class Scalar>
00168 RCP<const Thyra::NonlinearSolverBase<Scalar> >
00169 ImplicitBDFStepper<Scalar>::getSolver() const
00170 {
00171   return (solver_);
00172 }
00173 
00174 
00175 // Overridden from StepperBase
00176 
00177 
00178 template<class Scalar>
00179 bool ImplicitBDFStepper<Scalar>::isImplicit() const
00180 {
00181   return true;
00182 }
00183 
00184 template<class Scalar>
00185 bool ImplicitBDFStepper<Scalar>::supportsCloning() const
00186 {
00187   return true;
00188 }
00189 
00190 
00191 template<class Scalar>
00192 RCP<StepperBase<Scalar> >
00193 ImplicitBDFStepper<Scalar>::cloneStepperAlgorithm() const
00194 {
00195 
00196   // Just use the interface to clone the algorithm in an basically
00197   // uninitialized state
00198 
00199   RCP<ImplicitBDFStepper<Scalar> >
00200     stepper = Teuchos::rcp(new ImplicitBDFStepper<Scalar>());
00201 
00202   if (!is_null(model_))
00203     stepper->setModel(model_); // Shallow copy is okay!
00204 
00205   if (!is_null(solver_))
00206     stepper->setSolver(solver_->cloneNonlinearSolver().assert_not_null());
00207   
00208   if (!is_null(parameterList_))
00209     stepper->setParameterList(Teuchos::parameterList(*parameterList_));
00210 
00211   if (!is_null(stepControl_)) {
00212     if (stepControl_->supportsCloning())
00213       stepper->setStepControlStrategy(
00214         stepControl_->cloneStepControlStrategyAlgorithm().assert_not_null());
00215   }
00216   
00217   // At this point, we should have a valid algorithm state.  What might be
00218   // missing is the initial condition if it was not given in *model_ but was
00219   // set explicitly.  However, the specification for this function does not
00220   // guarantee that the full state will be copied in any case!
00221 
00222   return stepper;
00223 
00224 }
00225 
00226 
00227 template<class Scalar>
00228 void ImplicitBDFStepper<Scalar>::setModel(
00229   const RCP<const Thyra::ModelEvaluator<Scalar> >& model
00230   )
00231 {
00232   typedef Teuchos::ScalarTraits<Scalar> ST;
00233   TEUCHOS_TEST_FOR_EXCEPT( is_null(model) );
00234   assertValidModel( *this, *model );
00235   model_ = model;
00236 }
00237 
00238 template<class Scalar>
00239 void ImplicitBDFStepper<Scalar>::setNonconstModel(
00240   const RCP<Thyra::ModelEvaluator<Scalar> >& model
00241   )
00242 {
00243   this->setModel(model); // TODO 09/09/09 tscoffe:  use ConstNonconstObjectContainer!
00244 }
00245 
00246 
00247 template<class Scalar>
00248 RCP<const Thyra::ModelEvaluator<Scalar> >
00249 ImplicitBDFStepper<Scalar>::getModel() const
00250 {
00251   return model_;
00252 }
00253 
00254 
00255 template<class Scalar>
00256 RCP<Thyra::ModelEvaluator<Scalar> >
00257 ImplicitBDFStepper<Scalar>::getNonconstModel()
00258 {
00259   return Teuchos::null;
00260 }
00261 
00262 
00263 template<class Scalar>
00264 void ImplicitBDFStepper<Scalar>::setInitialCondition(
00265   const Thyra::ModelEvaluatorBase::InArgs<Scalar> &initialCondition
00266   )
00267 {
00268   typedef Teuchos::ScalarTraits<Scalar> ST;
00269   typedef Thyra::ModelEvaluatorBase MEB;
00270   TEUCHOS_TEST_FOR_EXCEPT(is_null(initialCondition.get_x()));
00271   TEUCHOS_TEST_FOR_EXCEPT(is_null(initialCondition.get_x_dot()));
00272   basePoint_ = initialCondition;
00273   xn0_ = initialCondition.get_x()->clone_v();
00274   xpn0_ = initialCondition.get_x_dot()->clone_v(); 
00275   time_ = initialCondition.get_t();
00276   // Generate vectors for use in calculations
00277   x_dot_base_ = createMember(xpn0_->space());
00278   V_S(x_dot_base_.ptr(),ST::zero());
00279   ee_ = createMember(xn0_->space());
00280   V_S(ee_.ptr(),ST::zero());
00281   // x history
00282   xHistory_.clear();
00283   xHistory_.push_back(xn0_->clone_v());
00284   xHistory_.push_back(xpn0_->clone_v());
00285 
00286   haveInitialCondition_ = true;
00287   if (isInitialized_) {
00288     initialize_();
00289   }
00290 }
00291 
00292 
00293 template<class Scalar>
00294 Thyra::ModelEvaluatorBase::InArgs<Scalar> 
00295 ImplicitBDFStepper<Scalar>::getInitialCondition() const
00296 {
00297   return basePoint_;
00298 }
00299 
00300 
00301 template<class Scalar>
00302 Scalar ImplicitBDFStepper<Scalar>::takeStep(Scalar dt, StepSizeType stepType)
00303 {
00304   
00305   RYTHMOS_FUNC_TIME_MONITOR("Rythmos::ImplicitBDFStepper::takeStep");
00306   
00307   using Teuchos::as;
00308   using Teuchos::incrVerbLevel;
00309   typedef Teuchos::ScalarTraits<Scalar> ST;
00310   typedef typename Thyra::ModelEvaluatorBase::InArgs<Scalar>::ScalarMag TScalarMag;
00311   typedef Thyra::NonlinearSolverBase<Scalar> NSB;
00312   typedef Teuchos::VerboseObjectTempState<NSB> VOTSNSB;
00313 
00314   RCP<Teuchos::FancyOStream> out = this->getOStream();
00315   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00316   Teuchos::OSTab ostab(out,1,"takeStep");
00317   VOTSNSB solver_outputTempState(solver_,out,incrVerbLevel(verbLevel,-1));
00318 
00319   if ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
00320     *out
00321       << "\nEntering " << this->Teuchos::Describable::description()
00322       << "::takeStep("<<dt<<","<<toString(stepType)<<") ...\n"; 
00323   }
00324 
00325   if (!isInitialized_) {
00326     initialize_(); 
00327   }
00328 
00329   stepControl_->setOStream(out);
00330   stepControl_->setVerbLevel(verbLevel);
00331   stepControl_->setRequestedStepSize(*this,dt,stepType);
00332 
00333   AttemptedStepStatusFlag status;
00334   while (1) {
00335     // Set up problem coefficients (and handle first step)
00336     Scalar hh_old = hh_;
00337     int desiredOrder;
00338     stepControl_->nextStepSize(*this,&hh_,&stepType,&desiredOrder);
00339     TEUCHOS_TEST_FOR_EXCEPT(!((1 <= desiredOrder) && (desiredOrder <= maxOrder_)));
00340     TEUCHOS_TEST_FOR_EXCEPT(!(desiredOrder <= usedOrder_+1));
00341     currentOrder_ = desiredOrder;
00342     if (numberOfSteps_ == 0) {
00343       psi_[0] = hh_;
00344       if (nef_ == 0) {
00345         Vt_S(xHistory_[1].ptr(),hh_);
00346       } else {
00347         Vt_S(xHistory_[1].ptr(),hh_/hh_old);
00348       }
00349     }
00350     this->updateCoeffs_();
00351     // compute predictor
00352     obtainPredictor_();
00353     // solve nonlinear problem (as follows)
00354     
00355     //
00356     // Setup the nonlinear equations:
00357     //
00358     //   f_bar( x_dot_coeff * x_bar + x_dot_base, x_coeff * x_bar + x_base, t_base ) = 0
00359     //   x_dot_coeff = -alpha_s/dt
00360     //   x_dot_base = x_prime_pred + (alpha_s/dt) * x_pred
00361     //   x_coeff = 1
00362     //   x_base = 0
00363     //   t_base = tn+dt
00364     //
00365     Scalar coeff_x_dot = Scalar(-ST::one())*alpha_s_/hh_;
00366     V_StVpStV( x_dot_base_.ptr(), ST::one(), *xpn0_, alpha_s_/hh_, *xn0_ );
00367     if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_EXTREME) ) {
00368       *out << "model_ = " << std::endl;
00369       model_->describe(*out,verbLevel);
00370       *out << "basePoint_ = " << std::endl;
00371       basePoint_.describe(*out,verbLevel);
00372       *out << "coeff_x_dot = " << coeff_x_dot << std::endl;
00373       *out << "x_dot_base_ = " << std::endl;
00374       x_dot_base_->describe(*out,verbLevel);
00375       *out << "time_+hh_ = " << time_+hh_ << std::endl;
00376       *out << "xn0_ = " << std::endl;
00377       xn0_->describe(*out,verbLevel);
00378     }
00379     neModel_.initializeSingleResidualModel(
00380       model_, basePoint_,
00381       coeff_x_dot, x_dot_base_,
00382       ST::one(), Teuchos::null,
00383       time_+hh_,
00384       xn0_
00385       );
00386     //
00387     // Solve the implicit nonlinear system to a tolerance of ???
00388     // 
00389     if (solver_->getModel().get() != &neModel_) {
00390       solver_->setModel( Teuchos::rcpFromRef(neModel_) );
00391     }
00392     /* // Rythmos::TimeStepNonlinearSolver uses a built in solveCriteria, so you can't pass one in.
00393     // I believe this is the correct solveCriteria for IDA though.
00394     Thyra::SolveMeasureType nonlinear_solve_measure_type(Thyra::SOLVE_MEASURE_NORM_RESIDUAL,Thyra::SOLVE_MEASURE_ONE); 
00395     TScalarMag tolerance = relErrTol_/TScalarMag(10.0); // This should be changed to match the condition in IDA
00396     Thyra::SolveCriteria<Scalar> nonlinearSolveCriteria(nonlinear_solve_measure_type, tolerance);
00397     Thyra::SolveStatus<Scalar> nonlinearSolveStatus = solver_->solve( &*xn0_, &nonlinearSolveCriteria, &*delta_ ); 
00398     */
00399     if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_EXTREME) ) {
00400       *out << "xn0 = " << std::endl;
00401       xn0_->describe(*out,verbLevel);
00402       *out << "ee = " << std::endl;
00403       ee_->describe(*out,verbLevel);
00404     }
00405     Thyra::SolveStatus<Scalar> nonlinearSolveStatus = solver_->solve( &*xn0_, NULL, &*ee_ ); 
00406     if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_EXTREME) ) {
00407       *out << "xn0 = " << std::endl;
00408       xn0_->describe(*out,verbLevel);
00409       *out << "ee = " << std::endl;
00410       ee_->describe(*out,verbLevel);
00411     }
00412     // In the above solve, on input *xn0_ is the initial guess that comes from
00413     // the predictor.  On output, *xn0_ is the solved for solution and *ee_ is
00414     // the difference computed from the intial guess in *xn0_ to the final
00415     // solved value of *xn0_.  This is needed for basic numerical stability.
00416     if (nonlinearSolveStatus.solveStatus == Thyra::SOLVE_STATUS_CONVERGED)  {
00417       newtonConvergenceStatus_ = 0;
00418     }
00419     else {
00420       newtonConvergenceStatus_ = -1;
00421     }
00422 
00423     // check error and evaluate LTE
00424     stepControl_->setCorrection(*this,xn0_,ee_,newtonConvergenceStatus_);
00425     bool stepPass = stepControl_->acceptStep(*this,&LETvalue_);
00426     
00427     if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00428       *out << "xn0_ = " << std::endl;
00429       xn0_->describe(*out,verbLevel);
00430       *out << "xpn0_ = " << std::endl;
00431       xpn0_->describe(*out,verbLevel);
00432       *out << "ee_ = " << std::endl;
00433       ee_->describe(*out,verbLevel);
00434       for (int i=0; i<std::max(2,maxOrder_); ++i) {
00435         *out << "xHistory_[" << i << "] = " << std::endl;
00436         xHistory_[i]->describe(*out,verbLevel);
00437       }
00438     }
00439 
00440     // Check LTE here:
00441     if (!stepPass) { // stepPass = false
00442       stepLETStatus_ = STEP_LET_STATUS_FAILED;
00443       status = stepControl_->rejectStep(*this);
00444       nef_++;
00445       if (status == CONTINUE_ANYWAY) {
00446         break;
00447       } else {
00448         restoreHistory_();
00449       }
00450     } else { // stepPass = true
00451       stepLETStatus_ = STEP_LET_STATUS_PASSED;
00452       break;
00453     }
00454   }
00455 
00456   // 08/22/07 the history array must be updated before stepControl_->completeStep.
00457   completeStep_();   
00458   stepControl_->completeStep(*this);
00459 
00460   if ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
00461     *out
00462       << "\nLeaving " << this->Teuchos::Describable::description()
00463       << "::takeStep("<<dt<<","<<toString(stepType)<<") ...\n"; 
00464   }
00465 
00466   return(usedStep_);
00467 
00468 }
00469 
00470 
00471 template<class Scalar>
00472 const StepStatus<Scalar> ImplicitBDFStepper<Scalar>::getStepStatus() const
00473 {
00474 
00475   // 2007/08/24: rabartl: We agreed that getStepStatus() would be free
00476   // so I have commented out removed all code that is not free
00477 
00478   typedef Teuchos::ScalarTraits<Scalar> ST;
00479   StepStatus<Scalar> stepStatus;
00480   if (!isInitialized_) {
00481     stepStatus.message = "This stepper is uninitialized.";
00482     stepStatus.stepStatus = STEP_STATUS_UNINITIALIZED;
00483     stepStatus.stepSize = Scalar(-ST::one());
00484     stepStatus.order = -1;
00485     stepStatus.time = Scalar(-ST::one());
00486     return(stepStatus);
00487   }
00488 
00489   if (numberOfSteps_ > 0) {
00490     stepStatus.stepStatus = STEP_STATUS_CONVERGED; 
00491   } else {
00492     stepStatus.stepStatus = STEP_STATUS_UNKNOWN;
00493   }
00494   stepStatus.stepLETStatus = stepLETStatus_;
00495   stepStatus.stepSize = usedStep_; 
00496   stepStatus.order = usedOrder_;
00497   stepStatus.time = time_;
00498   stepStatus.stepLETValue = LETvalue_; 
00499   stepStatus.solution = xHistory_[0];
00500   stepStatus.solutionDot = Teuchos::null; // This is *not* free!
00501   stepStatus.residual = Teuchos::null; // This is *not* free!
00502 
00503   return(stepStatus);
00504 
00505 }
00506 
00507 
00508 // Overridden from InterpolationBufferBase
00509 
00510 
00511 template<class Scalar>
00512 RCP<const Thyra::VectorSpaceBase<Scalar> >
00513 ImplicitBDFStepper<Scalar>::get_x_space() const
00514 {
00515   //TEUCHOS_TEST_FOR_EXCEPTION(!isInitialized_,std::logic_error,"Error, attempting to call get_x_space before initialization!\n");
00516   return ( !is_null(model_) ? model_->get_x_space() : Teuchos::null );
00517 }
00518 
00519 
00520 template<class Scalar>
00521 void ImplicitBDFStepper<Scalar>::addPoints(
00522   const Array<Scalar>& time_vec,
00523   const Array<RCP<const Thyra::VectorBase<Scalar> > >& x_vec,
00524   const Array<RCP<const Thyra::VectorBase<Scalar> > >& xdot_vec
00525   )
00526 {
00527   TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
00528     "Error, addPoints is not implemented for ImplicitBDFStepper.\n");
00529 }
00530 
00531 
00532 template<class Scalar>
00533 TimeRange<Scalar> ImplicitBDFStepper<Scalar>::getTimeRange() const
00534 {
00535   if ( !isInitialized_ && haveInitialCondition_ )
00536     return timeRange<Scalar>(time_,time_);
00537   if ( !isInitialized_ && !haveInitialCondition_ )
00538     return invalidTimeRange<Scalar>();
00539   return timeRange<Scalar>(time_-usedStep_,time_);
00540 }
00541 
00542 
00543 template<class Scalar>
00544 void ImplicitBDFStepper<Scalar>::getPoints(
00545   const Array<Scalar>& time_vec
00546   ,Array<RCP<const Thyra::VectorBase<Scalar> > >* x_vec
00547   ,Array<RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec
00548   ,Array<ScalarMag>* accuracy_vec) const
00549 {
00550   using Teuchos::as;
00551   using Teuchos::constOptInArg;
00552   using Teuchos::null;
00553   using Teuchos::ptr;
00554   typedef Teuchos::ScalarTraits<Scalar> ST;
00555   typedef typename ST::magnitudeType mScalarMag;
00556 
00557   TEUCHOS_ASSERT(haveInitialCondition_);
00558   // Only do this if we're being called pre-initialization to get the IC.
00559   if ( (numberOfSteps_ == -1) && 
00560        (time_vec.length() == 1) &&
00561        (compareTimeValues<Scalar>(time_vec[0],time_)==0) ) {
00562     defaultGetPoints<Scalar>(
00563         time_, constOptInArg(*xn0_), constOptInArg(*xpn0_),
00564         time_, constOptInArg(*xn0_), constOptInArg(*xpn0_),
00565         time_vec, ptr(x_vec), ptr(xdot_vec), ptr(accuracy_vec),
00566         Ptr<InterpolatorBase<Scalar> >(null)
00567         );
00568     return;
00569   }
00570   TEUCHOS_ASSERT(isInitialized_);
00571   RYTHMOS_FUNC_TIME_MONITOR("Rythmos::ImplicitBDFStepper::getPoints");
00572   if (x_vec)
00573     x_vec->clear();
00574   if (xdot_vec)
00575     xdot_vec->clear();
00576   for (Teuchos::Ordinal i=0 ; i<time_vec.size() ; ++i) {
00577     RCP<Thyra::VectorBase<Scalar> >
00578       x_temp = createMember(xn0_->space());
00579     RCP<Thyra::VectorBase<Scalar> >
00580       xdot_temp = createMember(xn0_->space());
00581     mScalarMag accuracy = -ST::zero();
00582     interpolateSolution_(
00583       time_vec[i], &*x_temp, &*xdot_temp,
00584       accuracy_vec ? &accuracy : 0
00585       );
00586     if (x_vec)
00587       x_vec->push_back(x_temp);
00588     if (xdot_vec)
00589       xdot_vec->push_back(xdot_temp);
00590     if (accuracy_vec)
00591       accuracy_vec->push_back(accuracy);
00592   }
00593   if ( as<int>(this->getVerbLevel()) >= as<int>(Teuchos::VERB_HIGH) ) {
00594     RCP<Teuchos::FancyOStream> out = this->getOStream();
00595     Teuchos::OSTab ostab(out,1,"getPoints");
00596     *out << "Passing out the interpolated values:" << std::endl;
00597     for (Teuchos::Ordinal i=0; i<time_vec.size() ; ++i) {
00598       *out << "time_[" << i << "] = " << time_vec[i] << std::endl;
00599       if (x_vec) {
00600         *out << "x_vec[" << i << "] = " << std::endl;
00601         (*x_vec)[i]->describe(*out,this->getVerbLevel());
00602       }
00603       if (xdot_vec) {
00604         *out << "xdot_vec[" << i << "] = ";
00605         if ( (*xdot_vec)[i] == Teuchos::null) {
00606           *out << "Teuchos::null" << std::endl;
00607         }
00608         else {
00609           *out << std::endl << Teuchos::describe(*(*xdot_vec)[i],this->getVerbLevel());
00610         }
00611       }
00612       if (accuracy_vec)
00613         *out << "accuracy[" << i << "] = " << (*accuracy_vec)[i] << std::endl;
00614     }
00615   }
00616 }
00617 
00618 
00619 template<class Scalar>
00620 void ImplicitBDFStepper<Scalar>::getNodes(Array<Scalar>* time_vec) const
00621 {
00622   TEUCHOS_ASSERT( time_vec != NULL );
00623   time_vec->clear();
00624   if (!haveInitialCondition_) {
00625     return;
00626   }
00627   if (numberOfSteps_ > 0) {
00628     time_vec->push_back(time_-usedStep_);
00629   }
00630   time_vec->push_back(time_);
00631 }
00632 
00633 
00634 template<class Scalar>
00635 void ImplicitBDFStepper<Scalar>::removeNodes(Array<Scalar>& time_vec) 
00636 {
00637   TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
00638     "Error, removeNodes is not implemented for ImplicitBDFStepper.\n");
00639 }
00640 
00641 
00642 template<class Scalar>
00643 int ImplicitBDFStepper<Scalar>::getOrder() const
00644 {
00645   if (!isInitialized_) {
00646     return(-1);
00647   }
00648   return(usedOrder_);
00649 }
00650 
00651 
00652 // Overridden from Teuchos::ParameterListAcceptor
00653 
00654 
00655 template<class Scalar>
00656 void ImplicitBDFStepper<Scalar>::setParameterList(
00657   RCP<Teuchos::ParameterList> const& paramList
00658   )
00659 {
00660   TEUCHOS_TEST_FOR_EXCEPT(paramList == Teuchos::null);
00661   paramList->validateParameters(*this->getValidParameters(),0);
00662   parameterList_ = paramList;
00663   Teuchos::readVerboseObjectSublist(&*parameterList_,this);
00664 }
00665 
00666 
00667 template<class Scalar>
00668 RCP<Teuchos::ParameterList> ImplicitBDFStepper<Scalar>::getNonconstParameterList()
00669 {
00670   return(parameterList_);
00671 }
00672 
00673 
00674 template<class Scalar>
00675 RCP<Teuchos::ParameterList>
00676 ImplicitBDFStepper<Scalar>::unsetParameterList()
00677 {
00678   RCP<Teuchos::ParameterList> temp_param_list = parameterList_;
00679   parameterList_ = Teuchos::null;
00680   return(temp_param_list);
00681 }
00682 
00683 
00684 template<class Scalar>
00685 RCP<const Teuchos::ParameterList>
00686 ImplicitBDFStepper<Scalar>::getValidParameters() const
00687 {
00688 
00689   static RCP<Teuchos::ParameterList> validPL;
00690 
00691   if (is_null(validPL)) {
00692 
00693     RCP<Teuchos::ParameterList>
00694       pl = Teuchos::parameterList();
00695 
00696     pl->sublist(RythmosStepControlSettings_name);
00697 
00698     Teuchos::setupVerboseObjectSublist(&*pl);
00699 
00700     validPL = pl;
00701 
00702   }
00703 
00704   RCP<Teuchos::FancyOStream> out = this->getOStream();
00705   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00706   Teuchos::OSTab ostab(out,1,"getValidParameters");
00707   if (Teuchos::as<int>(verbLevel) == Teuchos::VERB_HIGH) {
00708     *out << "Setting up valid parameterlist." << std::endl;
00709     validPL->print(*out);
00710   }
00711 
00712   return (validPL);
00713   
00714 }
00715 
00716 
00717 // Overridden from Teuchos::Describable
00718 
00719 
00720 template<class Scalar>
00721 std::string ImplicitBDFStepper<Scalar>::description() const
00722 {
00723   std::ostringstream out;
00724   out << this->Teuchos::Describable::description();
00725   const TimeRange<Scalar> timeRange = this->getTimeRange();
00726   if (timeRange.isValid())
00727     out << " (timeRange="<<timeRange<<")";
00728   else
00729     out << " (This stepper is not initialized yet)";
00730   out << std::endl;
00731   return out.str();
00732 }
00733 
00734 
00735 template<class Scalar>
00736 void ImplicitBDFStepper<Scalar>::describe(
00737   Teuchos::FancyOStream &out,
00738   const Teuchos::EVerbosityLevel verbLevel
00739   ) const
00740 {
00741 
00742   using Teuchos::as;
00743 
00744   if (!isInitialized_) {
00745     out << this->description();
00746     return;
00747   }
00748   
00749   if ( (as<int>(verbLevel) == as<int>(Teuchos::VERB_DEFAULT) ) ||
00750     (as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW)     )
00751     )
00752   {
00753     out << this->description() << std::endl;
00754     out << "model_ = " << Teuchos::describe(*model_,verbLevel);
00755     out << "solver_ = " << Teuchos::describe(*solver_,verbLevel);
00756     out << "neModel_ = " << Teuchos::describe(neModel_,verbLevel);
00757   }
00758   if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW)) {
00759     out << "time_ = " << time_ << std::endl;
00760     out << "hh_ = " << hh_ << std::endl;
00761     out << "currentOrder_ = " << currentOrder_ << std::endl;
00762   }
00763   if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH)) {
00764     out << "xn0_ = " << Teuchos::describe(*xn0_,verbLevel);
00765     out << "xpn0_ = " << Teuchos::describe(*xpn0_,verbLevel);
00766     out << "x_dot_base_ = " << Teuchos::describe(*x_dot_base_,verbLevel);
00767     for (int i=0 ; i < std::max(2,maxOrder_) ; ++i) {
00768       out << "xHistory_[" << i << "] = "
00769           << Teuchos::describe(*xHistory_[i],verbLevel);
00770     }
00771     out << "ee_ = " << Teuchos::describe(*ee_,verbLevel);
00772   }
00773 }
00774 
00775 
00776 // private
00777 
00778 
00779 // 2007/08/24: rabartl: Belos: We really should initialize all of this data in
00780 // a member initialization list but since there are like three constructors
00781 // this would mean that we would have to duplicate code (which is error prone)
00782 // or use a macro (which is not easy to debug).  We really should remove all
00783 // but the default constructor which then would set this data once in the
00784 // initialization list.
00785 
00786 template<class Scalar>
00787 void ImplicitBDFStepper<Scalar>::defaultInitializeAll_()
00788 {
00789   typedef Teuchos::ScalarTraits<Scalar> ST;
00790   const Scalar nan = ST::nan(), one = ST::one(), zero = ST::zero();
00791   // Initialize some data members to their rightful default values
00792   haveInitialCondition_ = false;
00793   isInitialized_ = false;
00794   currentOrder_ = 1;
00795   usedOrder_ = 1;
00796   usedStep_ = zero;
00797   // Initialize the rest of the private data members to invalid values to
00798   // avoid uninitialed memory
00799   time_ = nan;
00800   hh_ = nan;
00801   maxOrder_ = -1;
00802   LETvalue_ = -one;
00803   stepLETStatus_ = STEP_LET_STATUS_UNKNOWN;
00804   alpha_s_ = -one;
00805   numberOfSteps_ = -1;
00806   nef_ = -1;
00807   nscsco_ = -1;
00808   newtonConvergenceStatus_ = -1;
00809 }
00810 
00811 
00812 template<class Scalar>
00813 void ImplicitBDFStepper<Scalar>::obtainPredictor_()
00814 {
00815 
00816   using Teuchos::as;
00817   typedef Teuchos::ScalarTraits<Scalar> ST;
00818 
00819   if (!isInitialized_) {
00820     return;
00821   }
00822 
00823   RCP<Teuchos::FancyOStream> out = this->getOStream();
00824   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00825   Teuchos::OSTab ostab(out,1,"obtainPredictor_");
00826   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00827     *out << "currentOrder_ = " << currentOrder_ << std::endl;
00828   }
00829   
00830   // prepare history array for prediction
00831   for (int i=nscsco_;i<=currentOrder_;++i) {
00832     Vt_S(xHistory_[i].ptr(),beta_[i]);
00833   }
00834   
00835   // evaluate predictor
00836   V_V(xn0_.ptr(),*xHistory_[0]);
00837   V_S(xpn0_.ptr(),ST::zero());
00838   for (int i=1;i<=currentOrder_;++i) {
00839     Vp_V(xn0_.ptr(),*xHistory_[i]);
00840     Vp_StV(xpn0_.ptr(),gamma_[i],*xHistory_[i]);
00841   }
00842   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00843     *out << "xn0_ = " << std::endl;
00844     xn0_->describe(*out,verbLevel);
00845     *out << "xpn0_ = " << std::endl;
00846     xpn0_->describe(*out,verbLevel);
00847   }
00848 }
00849 
00850 
00851 template<class Scalar>
00852 void ImplicitBDFStepper<Scalar>::interpolateSolution_(
00853   const Scalar& timepoint,
00854   Thyra::VectorBase<Scalar>* x_ptr,
00855   Thyra::VectorBase<Scalar>* xdot_ptr,
00856   ScalarMag* accuracy_ptr
00857   ) const
00858 {
00859 
00860   typedef std::numeric_limits<Scalar> NL;
00861   typedef Teuchos::ScalarTraits<Scalar> ST;
00862 
00863 #ifdef HAVE_RYTHMOS_DEBUG
00864   TEUCHOS_TEST_FOR_EXCEPTION(
00865     !isInitialized_,std::logic_error,
00866     "Error, attempting to call interpolateSolution before initialization!\n");
00867   const TimeRange<Scalar> currTimeRange = this->getTimeRange();
00868   TEUCHOS_TEST_FOR_EXCEPTION(
00869     !currTimeRange.isInRange(timepoint), std::logic_error,
00870     "Error, timepoint = " << timepoint << " is not in the time range "
00871     << currTimeRange << "!" );
00872 #endif
00873 
00874   const Scalar tn = time_;
00875   const int kused = usedOrder_;
00876 
00877   // order of interpolation
00878   int kord = kused;
00879   if ( (kused == 0) || (timepoint == tn) )  {
00880     kord = 1;
00881   }
00882 
00883   // Initialize interploation
00884   Thyra::V_V(ptr(x_ptr),*xHistory_[0]);
00885   Thyra::V_S(ptr(xdot_ptr),ST::zero());
00886   
00887   // Add history array contributions
00888   const Scalar delt = timepoint - tn;
00889   Scalar c = ST::one(); // coefficient for interpolation of x
00890   Scalar d = ST::zero(); // coefficient for interpolation of xdot
00891   Scalar gam = delt/psi_[0]; // coefficient for interpolation
00892   for (int j=1 ; j <= kord ; ++j) {
00893     d = d*gam + c/psi_[j-1];
00894     c = c*gam;
00895     gam = (delt + psi_[j-1])/psi_[j];
00896     Thyra::Vp_StV(ptr(x_ptr),c,*xHistory_[j]);
00897     Thyra::Vp_StV(ptr(xdot_ptr),d,*xHistory_[j]);
00898   }
00899 
00900   // Set approximate accuracy
00901   if (accuracy_ptr) {
00902     *accuracy_ptr = pow(usedStep_,kord);
00903   }
00904   
00905 }
00906 
00907 
00908 template<class Scalar>
00909 void ImplicitBDFStepper<Scalar>::updateHistory_()
00910 {
00911 
00912   using Teuchos::as;
00913 
00914   // Save Newton correction for potential order increase on next step.
00915   if (usedOrder_ < maxOrder_)  {
00916     assign( xHistory_[usedOrder_+1].ptr(), *ee_ );
00917   }
00918   // Update history arrays
00919   Vp_V( xHistory_[usedOrder_].ptr(), *ee_ );
00920   for (int j=usedOrder_-1;j>=0;j--) {
00921     Vp_V( xHistory_[j].ptr(), *xHistory_[j+1] );
00922   }
00923   RCP<Teuchos::FancyOStream> out = this->getOStream();
00924   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00925   Teuchos::OSTab ostab(out,1,"updateHistory_");
00926   if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00927     for (int i=0;i<std::max(2,maxOrder_);++i) {
00928       *out << "xHistory_[" << i << "] = " << std::endl;
00929       xHistory_[i]->describe(*out,verbLevel);
00930     }
00931   }
00932 
00933 }
00934 
00935 
00936 template<class Scalar>
00937 void ImplicitBDFStepper<Scalar>::restoreHistory_()
00938 {
00939 
00940   using Teuchos::as;
00941   typedef Teuchos::ScalarTraits<Scalar> ST;
00942 
00943   // undo preparation of history array for prediction
00944   for (int i=nscsco_;i<=currentOrder_;++i) {
00945     Vt_S( xHistory_[i].ptr(), ST::one()/beta_[i] );
00946   }
00947   for (int i=1;i<=currentOrder_;++i) {
00948     psi_[i-1] = psi_[i] - hh_;
00949   }
00950   RCP<Teuchos::FancyOStream> out = this->getOStream();
00951   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00952   Teuchos::OSTab ostab(out,1,"restoreHistory_");
00953   if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00954     for (int i=0;i<maxOrder_;++i) {
00955       *out << "psi_[" << i << "] = " << psi_[i] << std::endl;
00956     }
00957     for (int i=0;i<maxOrder_;++i) {
00958       *out << "xHistory_[" << i << "] = " << std::endl;
00959       xHistory_[i]->describe(*out,verbLevel);
00960     }
00961   }
00962 
00963 } 
00964 
00965 
00966 template<class Scalar>
00967 void ImplicitBDFStepper<Scalar>::updateCoeffs_()
00968 {
00969 
00970   using Teuchos::as;
00971   typedef Teuchos::ScalarTraits<Scalar> ST;
00972 
00973   // If the number of steps taken with constant order and constant stepsize is
00974   // more than the current order + 1 then we don't bother to update the
00975   // coefficients because we've reached a constant step-size formula.  When
00976   // this is is not true, then we update the coefficients for the variable
00977   // step-sizes. 
00978   if ((hh_ != usedStep_) || (currentOrder_ != usedOrder_)) {
00979     nscsco_ = 0;
00980   }
00981   nscsco_ = std::min(nscsco_+1,usedOrder_+2);
00982   if (currentOrder_+1 >= nscsco_) {
00983     beta_[0] = ST::one();
00984     alpha_[0] = ST::one();
00985     Scalar temp1 = hh_;
00986     gamma_[0] = ST::zero();
00987     for (int i=1;i<=currentOrder_;++i) {
00988       Scalar temp2 = psi_[i-1];
00989       psi_[i-1] = temp1;
00990       beta_[i] = beta_[i-1]*psi_[i-1]/temp2;
00991       temp1 = temp2 + hh_;
00992       alpha_[i] = hh_/temp1;
00993       gamma_[i] = gamma_[i-1]+alpha_[i-1]/hh_;
00994     }
00995     psi_[currentOrder_] = temp1;
00996   }
00997   alpha_s_ = ST::zero();
00998   for (int i=0;i<currentOrder_;++i) {
00999     alpha_s_ = alpha_s_ - Scalar(ST::one()/(i+ST::one()));
01000   }
01001   RCP<Teuchos::FancyOStream> out = this->getOStream();
01002   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
01003   Teuchos::OSTab ostab(out,1,"updateCoeffs_");
01004   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
01005     for (int i=0;i<=maxOrder_;++i) {
01006       *out << "alpha_[" << i << "] = " << alpha_[i] << std::endl;
01007       *out << "beta_[" << i << "] = " << beta_[i] << std::endl;
01008       *out << "gamma_[" << i << "] = " << gamma_[i] << std::endl;
01009       *out << "psi_[" << i << "] = " << psi_[i] << std::endl;
01010       *out << "alpha_s_ = " << alpha_s_ << std::endl;
01011     }
01012     //std::cout << "alpha_s_ = " << alpha_s_ << std::endl;
01013     //for (int i=0;i<=maxOrder_;++i) {
01014     //  std::cout << "  alpha_[" << i << "] = " << alpha_[i] << std::endl;
01015     //  std::cout << "  beta_[" << i << "] = " << beta_[i] << std::endl;
01016     //  std::cout << "  gamma_[" << i << "] = " << gamma_[i] << std::endl;
01017     //  std::cout << "  psi_[" << i << "] = " << psi_[i] << std::endl;
01018     //}
01019   }
01020 }
01021 
01022 
01023 template<class Scalar>
01024 void ImplicitBDFStepper<Scalar>::initialize_()
01025 {
01026 
01027   using Teuchos::as;
01028   typedef Teuchos::ScalarTraits<Scalar> ST;
01029   using Thyra::createMember;
01030 
01031   RCP<Teuchos::FancyOStream> out = this->getOStream();
01032   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
01033   const bool doTrace = (as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH));
01034   Teuchos::OSTab ostab(out,1,"initialize_");
01035 
01036   if (doTrace) {
01037     *out
01038       << "\nEntering " << this->Teuchos::Describable::description()
01039       << "::initialize_()...\n";
01040   }
01041 
01042   TEUCHOS_TEST_FOR_EXCEPT(model_ == Teuchos::null);
01043   TEUCHOS_TEST_FOR_EXCEPT(solver_ == Teuchos::null);
01044   TEUCHOS_ASSERT(haveInitialCondition_);
01045 
01046   // Initialize Parameter List if none provided.
01047   if (parameterList_ == Teuchos::null) {
01048     RCP<Teuchos::ParameterList> emptyParameterList = Teuchos::rcp(new Teuchos::ParameterList);
01049     this->setParameterList(emptyParameterList);
01050   }
01051 
01052   // Initialize StepControl
01053   if (stepControl_ == Teuchos::null) {
01054     RCP<ImplicitBDFStepperStepControl<Scalar> > implicitBDFStepperStepControl =
01055       Teuchos::rcp(new ImplicitBDFStepperStepControl<Scalar>());
01056     RCP<Teuchos::ParameterList> stepControlPL = 
01057       Teuchos::sublist(parameterList_, RythmosStepControlSettings_name);
01058     implicitBDFStepperStepControl->setParameterList(stepControlPL);
01059     this->setStepControlStrategy(implicitBDFStepperStepControl);
01060   }
01061   stepControl_->initialize(*this);
01062 
01063   maxOrder_ = stepControl_->getMaxOrder(); // maximum order
01064   TEUCHOS_TEST_FOR_EXCEPTION(
01065       !((1 <= maxOrder_) && (maxOrder_ <= 5)), std::logic_error,
01066       "Error, maxOrder returned from stepControl_->getMaxOrder() = " << maxOrder_ << " is outside range of [1,5]!\n"
01067       );
01068 
01069   Scalar zero = ST::zero();
01070 
01071   currentOrder_ = 1; // Current order of integration
01072   usedOrder_ = 1;  // order used in current step (used after currentOrder_ is updated)
01073   usedStep_ = zero;
01074   nscsco_  =  0;
01075   LETvalue_ = zero;
01076 
01077   alpha_.clear();  // $\alpha_j(n)=h_n/\psi_j(n)$ coefficient used in local error test
01078   // note:   $h_n$ = current step size, n = current time step
01079   gamma_.clear();  // calculate time derivative of history array for predictor 
01080   beta_.clear();   // coefficients used to evaluate predictor from history array
01081   psi_.clear();    // $\psi_j(n) = t_n-t_{n-j}$ intermediary variable used to 
01082   // compute $\beta_j(n):$
01083   for (int i=0 ; i<=maxOrder_ ; ++i) {
01084     alpha_.push_back(zero);
01085     beta_.push_back(zero);
01086     gamma_.push_back(zero);
01087     psi_.push_back(zero);
01088   }
01089   alpha_s_=Scalar(-ST::one());  // $\alpha_s$ fixed-leading coefficient of this BDF method
01090   hh_=zero;
01091   numberOfSteps_=0;   // number of total time integration steps taken
01092   nef_ = 0;
01093 
01094   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
01095     *out << "alpha_s_ = " << alpha_s_ << std::endl;
01096     for (int i=0 ; i<=maxOrder_ ; ++i) {
01097       *out << "alpha_[" << i << "] = " << alpha_[i] << std::endl;
01098       *out << "beta_[" << i << "] = " << beta_[i] << std::endl;
01099       *out << "gamma_[" << i << "] = " << gamma_[i] << std::endl;
01100       *out << "psi_[" << i << "] = " << psi_[i] << std::endl;
01101     }
01102     *out << "numberOfSteps_ = " << numberOfSteps_ << std::endl;
01103   }
01104 
01105   // setInitialCondition initialized xHistory with xn0, xpn0.  Now we add the rest of the vectors.
01106   // Store maxOrder_+1 vectors
01107   for (int i=2 ; i<=maxOrder_ ; ++i) {
01108     xHistory_.push_back(createMember(xn0_->space())); 
01109     V_S(xHistory_[i].ptr(),zero);
01110   }
01111 
01112   isInitialized_ = true;
01113 
01114   if (doTrace) {
01115     *out
01116       << "\nLeaving " << this->Teuchos::Describable::description()
01117       << "::initialize_()...\n";
01118   }
01119 
01120 }
01121 
01122 
01123 template<class Scalar>
01124 void ImplicitBDFStepper<Scalar>::completeStep_()
01125 {
01126 
01127   using Teuchos::as;
01128   typedef Teuchos::ScalarTraits<Scalar> ST;
01129 
01130 #ifdef HAVE_RYTHMOS_DEBUG
01131   TEUCHOS_TEST_FOR_EXCEPT(ST::isnaninf(hh_));
01132 #endif  
01133 
01134   numberOfSteps_ ++;
01135   nef_ = 0;
01136   usedStep_ = hh_;
01137   usedOrder_ = currentOrder_;
01138   time_ += hh_;
01139   RCP<Teuchos::FancyOStream> out = this->getOStream();
01140   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
01141   Teuchos::OSTab ostab(out,1,"completeStep_");
01142 
01143   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
01144     *out << "numberOfSteps_ = " << numberOfSteps_ << std::endl;
01145     *out << "time_ = " << time_ << std::endl;
01146   }
01147   
01148   // 03/22/04 tscoffe:  Note that updating the history has nothing to do with
01149   // the step-size and everything to do with the newton correction vectors.
01150   updateHistory_();
01151 
01152 }
01153 
01154 template<class Scalar>
01155 void ImplicitBDFStepper<Scalar>::setStepControlData(const StepperBase<Scalar> & stepper)
01156 {
01157   if (!isInitialized_) {
01158     initialize_();
01159   }
01160   stepControl_->setStepControlData(stepper);
01161 }
01162 
01163 // 
01164 // Explicit Instantiation macro
01165 //
01166 // Must be expanded from within the Rythmos namespace!
01167 //
01168 
01169 #define RYTHMOS_IMPLICITBDF_STEPPER_INSTANT(SCALAR) \
01170   \
01171   template class ImplicitBDFStepper< SCALAR >; \
01172   \
01173   template RCP< ImplicitBDFStepper< SCALAR > > \
01174   implicitBDFStepper();  \
01175   \
01176   template RCP< ImplicitBDFStepper< SCALAR > > \
01177   implicitBDFStepper( \
01178     const RCP<Thyra::ModelEvaluator< SCALAR > >& model, \
01179     const RCP<Thyra::NonlinearSolverBase< SCALAR > >& solver, \
01180     const RCP<Teuchos::ParameterList>& parameterList \
01181     ); \
01182 
01183 
01184 } // namespace Rythmos
01185 
01186 
01187 #endif //Rythmos_IMPLICITBDF_STEPPER_DEF_H
 All Classes Functions Variables Typedefs Friends