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., 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_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) &&
00340                               (desiredOrder <= maxOrder_)));
00341     TEUCHOS_TEST_FOR_EXCEPT(!(desiredOrder <= usedOrder_+1));
00342     currentOrder_ = desiredOrder;
00343     if (numberOfSteps_ == 0) {
00344       psi_[0] = hh_;
00345       if (nef_ == 0) {
00346         Vt_S(xHistory_[1].ptr(),hh_);
00347       } else {
00348         Vt_S(xHistory_[1].ptr(),hh_/hh_old);
00349       }
00350     }
00351     this->updateCoeffs_();
00352     // compute predictor
00353     obtainPredictor_();
00354     // solve nonlinear problem (as follows)
00355 
00356     //
00357     // Setup the nonlinear equations:
00358     //
00359     //   f_bar( x_dot_coeff * x_bar + x_dot_base,
00360     //          x_coeff * x_bar + x_base, t_base ) = 0
00361     //   x_dot_coeff = -alpha_s/dt
00362     //   x_dot_base = x_prime_pred + (alpha_s/dt) * x_pred
00363     //   x_coeff = 1
00364     //   x_base = 0
00365     //   t_base = tn+dt
00366     //
00367     Scalar coeff_x_dot = Scalar(-ST::one())*alpha_s_/hh_;
00368     V_StVpStV( x_dot_base_.ptr(), ST::one(), *xpn0_, alpha_s_/hh_, *xn0_ );
00369     if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_EXTREME) ) {
00370       *out << "model_ = " << std::endl;
00371       model_->describe(*out,verbLevel);
00372       *out << "basePoint_ = " << std::endl;
00373       basePoint_.describe(*out,verbLevel);
00374       *out << "coeff_x_dot = " << coeff_x_dot << std::endl;
00375       *out << "x_dot_base_ = " << std::endl;
00376       x_dot_base_->describe(*out,verbLevel);
00377       *out << "time_+hh_ = " << time_+hh_ << std::endl;
00378       *out << "xn0_ = " << std::endl;
00379       xn0_->describe(*out,verbLevel);
00380     }
00381     neModel_.initializeSingleResidualModel(
00382       model_, basePoint_,
00383       coeff_x_dot, x_dot_base_,
00384       ST::one(), Teuchos::null,
00385       time_+hh_,
00386       xn0_
00387       );
00388     //
00389     // Solve the implicit nonlinear system to a tolerance of ???
00390     //
00391     if (solver_->getModel().get() != &neModel_) {
00392       solver_->setModel( Teuchos::rcpFromRef(neModel_) );
00393     }
00394     // Rythmos::TimeStepNonlinearSolver uses a built in solveCriteria,
00395     // so you can't pass one in.
00396     // I believe this is the correct solveCriteria for IDA though.
00397     /*
00398     Thyra::SolveMeasureType nonlinear_solve_measure_type(
00399       Thyra::SOLVE_MEASURE_NORM_RESIDUAL,Thyra::SOLVE_MEASURE_ONE);
00400     // This should be changed to match the condition in IDA
00401     TScalarMag tolerance = relErrTol_/TScalarMag(10.0);
00402     Thyra::SolveCriteria<Scalar> nonlinearSolveCriteria(
00403       nonlinear_solve_measure_type, tolerance);
00404     Thyra::SolveStatus<Scalar> nonlinearSolveStatus =
00405       solver_->solve( &*xn0_, &nonlinearSolveCriteria, &*delta_ );
00406     */
00407     if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_EXTREME) ) {
00408       *out << "xn0 = " << std::endl;
00409       xn0_->describe(*out,verbLevel);
00410       *out << "ee = " << std::endl;
00411       ee_->describe(*out,verbLevel);
00412     }
00413     nonlinearSolveStatus_ = solver_->solve( &*xn0_, NULL, &*ee_ );
00414     if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_EXTREME) ) {
00415       *out << "xn0 = " << std::endl;
00416       xn0_->describe(*out,verbLevel);
00417       *out << "ee = " << std::endl;
00418       ee_->describe(*out,verbLevel);
00419     }
00420     // In the above solve, on input *xn0_ is the initial guess that comes from
00421     // the predictor.  On output, *xn0_ is the solved for solution and *ee_ is
00422     // the difference computed from the intial guess in *xn0_ to the final
00423     // solved value of *xn0_.  This is needed for basic numerical stability.
00424     if (nonlinearSolveStatus_.solveStatus == Thyra::SOLVE_STATUS_CONVERGED)  {
00425       newtonConvergenceStatus_ = 0;
00426     }
00427     else {
00428       newtonConvergenceStatus_ = -1;
00429     }
00430 
00431     // check error and evaluate LTE
00432     stepControl_->setCorrection(*this,xn0_,ee_,newtonConvergenceStatus_);
00433     bool stepPass = stepControl_->acceptStep(*this,&LETvalue_);
00434 
00435     if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00436       *out << "xn0_ = " << std::endl;
00437       xn0_->describe(*out,verbLevel);
00438       *out << "xpn0_ = " << std::endl;
00439       xpn0_->describe(*out,verbLevel);
00440       *out << "ee_ = " << std::endl;
00441       ee_->describe(*out,verbLevel);
00442       for (int i=0; i<std::max(2,maxOrder_); ++i) {
00443         *out << "xHistory_[" << i << "] = " << std::endl;
00444         xHistory_[i]->describe(*out,verbLevel);
00445       }
00446     }
00447 
00448     // Check LTE here:
00449     if (!stepPass) { // stepPass = false
00450       stepLETStatus_ = STEP_LET_STATUS_FAILED;
00451       status = stepControl_->rejectStep(*this);
00452       nef_++;
00453       if (status == CONTINUE_ANYWAY) {
00454         break;
00455       } else {
00456         restoreHistory_();
00457       }
00458     } else { // stepPass = true
00459       stepLETStatus_ = STEP_LET_STATUS_PASSED;
00460       break;
00461     }
00462   }
00463 
00464   // 08/22/07 the history array must be updated before
00465   // stepControl_->completeStep.
00466   completeStep_();
00467   stepControl_->completeStep(*this);
00468 
00469   if ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
00470     *out
00471       << "\nLeaving " << this->Teuchos::Describable::description()
00472       << "::takeStep("<<dt<<","<<toString(stepType)<<") ...\n";
00473   }
00474 
00475   return(usedStep_);
00476 }
00477 
00478 
00479 template<class Scalar>
00480 const StepStatus<Scalar> ImplicitBDFStepper<Scalar>::getStepStatus() const
00481 {
00482 
00483   // 2007/08/24: rabartl: We agreed that getStepStatus() would be free
00484   // so I have commented out removed all code that is not free
00485 
00486   typedef Teuchos::ScalarTraits<Scalar> ST;
00487   StepStatus<Scalar> stepStatus;
00488   if (!isInitialized_) {
00489     stepStatus.message = "This stepper is uninitialized.";
00490     stepStatus.stepStatus = STEP_STATUS_UNINITIALIZED;
00491     stepStatus.stepSize = Scalar(-ST::one());
00492     stepStatus.order = -1;
00493     stepStatus.time = Scalar(-ST::one());
00494 //    return(stepStatus);
00495   }
00496   else if (numberOfSteps_ > 0) {
00497     stepStatus.stepStatus = STEP_STATUS_CONVERGED;
00498   } else {
00499     stepStatus.stepStatus = STEP_STATUS_UNKNOWN;
00500   }
00501 
00502   stepStatus.stepLETStatus = stepLETStatus_;
00503   stepStatus.stepSize = usedStep_;
00504   stepStatus.order = usedOrder_;
00505   stepStatus.time = time_;
00506   stepStatus.stepLETValue = LETvalue_;
00507   if(xHistory_.size() > 0)
00508     stepStatus.solution = xHistory_[0];
00509   else
00510     stepStatus.solution = Teuchos::null;
00511   stepStatus.solutionDot = Teuchos::null; // This is *not* free!
00512   stepStatus.residual = Teuchos::null; // This is *not* free!
00513 
00514   return(stepStatus);
00515 
00516 }
00517 
00518 
00519 // Overridden from InterpolationBufferBase
00520 
00521 
00522 template<class Scalar>
00523 RCP<const Thyra::VectorSpaceBase<Scalar> >
00524 ImplicitBDFStepper<Scalar>::get_x_space() const
00525 {
00526   //TEUCHOS_TEST_FOR_EXCEPTION(!isInitialized_,std::logic_error,"Error, attempting to call get_x_space before initialization!\n");
00527   return ( !is_null(model_) ? model_->get_x_space() : Teuchos::null );
00528 }
00529 
00530 
00531 template<class Scalar>
00532 void ImplicitBDFStepper<Scalar>::addPoints(
00533   const Array<Scalar>& time_vec,
00534   const Array<RCP<const Thyra::VectorBase<Scalar> > >& x_vec,
00535   const Array<RCP<const Thyra::VectorBase<Scalar> > >& xdot_vec
00536   )
00537 {
00538   TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
00539     "Error, addPoints is not implemented for ImplicitBDFStepper.\n");
00540 }
00541 
00542 
00543 template<class Scalar>
00544 TimeRange<Scalar> ImplicitBDFStepper<Scalar>::getTimeRange() const
00545 {
00546   if ( !isInitialized_ && haveInitialCondition_ )
00547     return timeRange<Scalar>(time_,time_);
00548   if ( !isInitialized_ && !haveInitialCondition_ )
00549     return invalidTimeRange<Scalar>();
00550   return timeRange<Scalar>(time_-usedStep_,time_);
00551 }
00552 
00553 
00554 template<class Scalar>
00555 void ImplicitBDFStepper<Scalar>::getPoints(
00556   const Array<Scalar>& time_vec
00557   ,Array<RCP<const Thyra::VectorBase<Scalar> > >* x_vec
00558   ,Array<RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec
00559   ,Array<ScalarMag>* accuracy_vec) const
00560 {
00561   using Teuchos::as;
00562   using Teuchos::constOptInArg;
00563   using Teuchos::null;
00564   using Teuchos::ptr;
00565   typedef Teuchos::ScalarTraits<Scalar> ST;
00566   typedef typename ST::magnitudeType mScalarMag;
00567 
00568   TEUCHOS_ASSERT(haveInitialCondition_);
00569   // Only do this if we're being called pre-initialization to get the IC.
00570   if ( (numberOfSteps_ == -1) &&
00571        (time_vec.length() == 1) &&
00572        (compareTimeValues<Scalar>(time_vec[0],time_)==0) ) {
00573     defaultGetPoints<Scalar>(
00574         time_, constOptInArg(*xn0_), constOptInArg(*xpn0_),
00575         time_, constOptInArg(*xn0_), constOptInArg(*xpn0_),
00576         time_vec, ptr(x_vec), ptr(xdot_vec), ptr(accuracy_vec),
00577         Ptr<InterpolatorBase<Scalar> >(null)
00578         );
00579     return;
00580   }
00581   TEUCHOS_ASSERT(isInitialized_);
00582   RYTHMOS_FUNC_TIME_MONITOR("Rythmos::ImplicitBDFStepper::getPoints");
00583   if (x_vec)
00584     x_vec->clear();
00585   if (xdot_vec)
00586     xdot_vec->clear();
00587   for (Teuchos::Ordinal i=0 ; i<time_vec.size() ; ++i) {
00588     RCP<Thyra::VectorBase<Scalar> >
00589       x_temp = createMember(xn0_->space());
00590     RCP<Thyra::VectorBase<Scalar> >
00591       xdot_temp = createMember(xn0_->space());
00592     mScalarMag accuracy = -ST::zero();
00593     interpolateSolution_(
00594       time_vec[i], &*x_temp, &*xdot_temp,
00595       accuracy_vec ? &accuracy : 0
00596       );
00597     if (x_vec)
00598       x_vec->push_back(x_temp);
00599     if (xdot_vec)
00600       xdot_vec->push_back(xdot_temp);
00601     if (accuracy_vec)
00602       accuracy_vec->push_back(accuracy);
00603   }
00604   if ( as<int>(this->getVerbLevel()) >= as<int>(Teuchos::VERB_HIGH) ) {
00605     RCP<Teuchos::FancyOStream> out = this->getOStream();
00606     Teuchos::OSTab ostab(out,1,"getPoints");
00607     *out << "Passing out the interpolated values:" << std::endl;
00608     for (Teuchos::Ordinal i=0; i<time_vec.size() ; ++i) {
00609       *out << "time_[" << i << "] = " << time_vec[i] << std::endl;
00610       if (x_vec) {
00611         *out << "x_vec[" << i << "] = " << std::endl;
00612         (*x_vec)[i]->describe(*out,this->getVerbLevel());
00613       }
00614       if (xdot_vec) {
00615         *out << "xdot_vec[" << i << "] = ";
00616         if ( (*xdot_vec)[i] == Teuchos::null) {
00617           *out << "Teuchos::null" << std::endl;
00618         }
00619         else {
00620           *out << std::endl << Teuchos::describe(*(*xdot_vec)[i],this->getVerbLevel());
00621         }
00622       }
00623       if (accuracy_vec)
00624         *out << "accuracy[" << i << "] = " << (*accuracy_vec)[i] << std::endl;
00625     }
00626   }
00627 }
00628 
00629 
00630 template<class Scalar>
00631 void ImplicitBDFStepper<Scalar>::getNodes(Array<Scalar>* time_vec) const
00632 {
00633   TEUCHOS_ASSERT( time_vec != NULL );
00634   time_vec->clear();
00635   if (!haveInitialCondition_) {
00636     return;
00637   }
00638   if (numberOfSteps_ > 0) {
00639     time_vec->push_back(time_-usedStep_);
00640   }
00641   time_vec->push_back(time_);
00642 }
00643 
00644 
00645 template<class Scalar>
00646 void ImplicitBDFStepper<Scalar>::removeNodes(Array<Scalar>& time_vec)
00647 {
00648   TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
00649     "Error, removeNodes is not implemented for ImplicitBDFStepper.\n");
00650 }
00651 
00652 
00653 template<class Scalar>
00654 int ImplicitBDFStepper<Scalar>::getOrder() const
00655 {
00656   if (!isInitialized_) {
00657     return(-1);
00658   }
00659   return(usedOrder_);
00660 }
00661 
00662 
00663 // Overridden from Teuchos::ParameterListAcceptor
00664 
00665 
00666 template<class Scalar>
00667 void ImplicitBDFStepper<Scalar>::setParameterList(
00668   RCP<Teuchos::ParameterList> const& paramList
00669   )
00670 {
00671   TEUCHOS_TEST_FOR_EXCEPT(paramList == Teuchos::null);
00672   paramList->validateParameters(*this->getValidParameters(),0);
00673   parameterList_ = paramList;
00674   Teuchos::readVerboseObjectSublist(&*parameterList_,this);
00675 }
00676 
00677 
00678 template<class Scalar>
00679 RCP<Teuchos::ParameterList> ImplicitBDFStepper<Scalar>::getNonconstParameterList()
00680 {
00681   return(parameterList_);
00682 }
00683 
00684 
00685 template<class Scalar>
00686 RCP<Teuchos::ParameterList>
00687 ImplicitBDFStepper<Scalar>::unsetParameterList()
00688 {
00689   RCP<Teuchos::ParameterList> temp_param_list = parameterList_;
00690   parameterList_ = Teuchos::null;
00691   return(temp_param_list);
00692 }
00693 
00694 
00695 template<class Scalar>
00696 RCP<const Teuchos::ParameterList>
00697 ImplicitBDFStepper<Scalar>::getValidParameters() const
00698 {
00699   static RCP<Teuchos::ParameterList> validPL;
00700   if (is_null(validPL)) {
00701     RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
00702     // This line is required to pass StepperValidator UnitTest!
00703     pl->sublist(RythmosStepControlSettings_name);
00704     Teuchos::setupVerboseObjectSublist(&*pl);
00705     validPL = pl;
00706   }
00707 
00708   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00709   if (Teuchos::as<int>(verbLevel) == Teuchos::VERB_HIGH) {
00710     RCP<Teuchos::FancyOStream> out = this->getOStream();
00711     Teuchos::OSTab ostab(out,1,"getValidParameters");
00712     *out << "Setting up valid parameterlist." << std::endl;
00713     validPL->print(*out);
00714   }
00715 
00716   return (validPL);
00717 
00718 }
00719 
00720 
00721 // Overridden from Teuchos::Describable
00722 
00723 
00724 template<class Scalar>
00725 std::string ImplicitBDFStepper<Scalar>::description() const
00726 {
00727   std::ostringstream out;
00728   out << this->Teuchos::Describable::description();
00729   const TimeRange<Scalar> timeRange = this->getTimeRange();
00730   if (timeRange.isValid())
00731     out << " (timeRange="<<timeRange<<")";
00732   else
00733     out << " (This stepper is not initialized yet)";
00734   out << std::endl;
00735   return out.str();
00736 }
00737 
00738 
00739 template<class Scalar>
00740 void ImplicitBDFStepper<Scalar>::describe(
00741   Teuchos::FancyOStream &out,
00742   const Teuchos::EVerbosityLevel verbLevel
00743   ) const
00744 {
00745 
00746   using Teuchos::as;
00747 
00748   if (!isInitialized_) {
00749     out << this->description();
00750     return;
00751   }
00752 
00753   if ( (as<int>(verbLevel) == as<int>(Teuchos::VERB_DEFAULT) ) ||
00754     (as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW)     )
00755     )
00756   {
00757     out << this->description() << std::endl;
00758     out << "model_ = " << Teuchos::describe(*model_,verbLevel);
00759     out << "solver_ = " << Teuchos::describe(*solver_,verbLevel);
00760     out << "neModel_ = " << Teuchos::describe(neModel_,verbLevel);
00761   }
00762   if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW)) {
00763     out << "time_ = " << time_ << std::endl;
00764     out << "hh_ = " << hh_ << std::endl;
00765     out << "currentOrder_ = " << currentOrder_ << std::endl;
00766   }
00767   if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH)) {
00768     out << "xn0_ = " << Teuchos::describe(*xn0_,verbLevel);
00769     out << "xpn0_ = " << Teuchos::describe(*xpn0_,verbLevel);
00770     out << "x_dot_base_ = " << Teuchos::describe(*x_dot_base_,verbLevel);
00771     for (int i=0 ; i < std::max(2,maxOrder_) ; ++i) {
00772       out << "xHistory_[" << i << "] = "
00773           << Teuchos::describe(*xHistory_[i],verbLevel);
00774     }
00775     out << "ee_ = " << Teuchos::describe(*ee_,verbLevel);
00776   }
00777 }
00778 
00779 
00780 // private
00781 
00782 
00783 // 2007/08/24: rabartl: Belos: We really should initialize all of this data in
00784 // a member initialization list but since there are like three constructors
00785 // this would mean that we would have to duplicate code (which is error prone)
00786 // or use a macro (which is not easy to debug).  We really should remove all
00787 // but the default constructor which then would set this data once in the
00788 // initialization list.
00789 
00790 template<class Scalar>
00791 void ImplicitBDFStepper<Scalar>::defaultInitializeAll_()
00792 {
00793   typedef Teuchos::ScalarTraits<Scalar> ST;
00794   const Scalar nan = ST::nan(), one = ST::one(), zero = ST::zero();
00795   // Initialize some data members to their rightful default values
00796   haveInitialCondition_ = false;
00797   isInitialized_ = false;
00798   currentOrder_ = 1;
00799   usedOrder_ = 1;
00800   usedStep_ = zero;
00801   // Initialize the rest of the private data members to invalid values to
00802   // avoid uninitialed memory
00803   time_ = nan;
00804   hh_ = nan;
00805   maxOrder_ = -1;
00806   LETvalue_ = -one;
00807   stepLETStatus_ = STEP_LET_STATUS_UNKNOWN;
00808   alpha_s_ = -one;
00809   numberOfSteps_ = -1;
00810   nef_ = -1;
00811   nscsco_ = -1;
00812   newtonConvergenceStatus_ = -1;
00813 }
00814 
00815 
00816 template<class Scalar>
00817 void ImplicitBDFStepper<Scalar>::obtainPredictor_()
00818 {
00819 
00820   using Teuchos::as;
00821   typedef Teuchos::ScalarTraits<Scalar> ST;
00822 
00823   if (!isInitialized_) {
00824     return;
00825   }
00826 
00827   RCP<Teuchos::FancyOStream> out = this->getOStream();
00828   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00829   Teuchos::OSTab ostab(out,1,"obtainPredictor_");
00830   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00831     *out << "currentOrder_ = " << currentOrder_ << std::endl;
00832   }
00833 
00834   // prepare history array for prediction
00835   for (int i=nscsco_;i<=currentOrder_;++i) {
00836     Vt_S(xHistory_[i].ptr(),beta_[i]);
00837   }
00838 
00839   // evaluate predictor
00840   V_V(xn0_.ptr(),*xHistory_[0]);
00841   V_S(xpn0_.ptr(),ST::zero());
00842   for (int i=1;i<=currentOrder_;++i) {
00843     Vp_V(xn0_.ptr(),*xHistory_[i]);
00844     Vp_StV(xpn0_.ptr(),gamma_[i],*xHistory_[i]);
00845   }
00846   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00847     *out << "xn0_ = " << std::endl;
00848     xn0_->describe(*out,verbLevel);
00849     *out << "xpn0_ = " << std::endl;
00850     xpn0_->describe(*out,verbLevel);
00851   }
00852 }
00853 
00854 
00855 template<class Scalar>
00856 void ImplicitBDFStepper<Scalar>::interpolateSolution_(
00857   const Scalar& timepoint,
00858   Thyra::VectorBase<Scalar>* x_ptr,
00859   Thyra::VectorBase<Scalar>* xdot_ptr,
00860   ScalarMag* accuracy_ptr
00861   ) const
00862 {
00863 
00864   typedef std::numeric_limits<Scalar> NL;
00865   typedef Teuchos::ScalarTraits<Scalar> ST;
00866 
00867 #ifdef HAVE_RYTHMOS_DEBUG
00868   TEUCHOS_TEST_FOR_EXCEPTION(
00869     !isInitialized_,std::logic_error,
00870     "Error, attempting to call interpolateSolution before initialization!\n");
00871   const TimeRange<Scalar> currTimeRange = this->getTimeRange();
00872   TEUCHOS_TEST_FOR_EXCEPTION(
00873     !currTimeRange.isInRange(timepoint), std::logic_error,
00874     "Error, timepoint = " << timepoint << " is not in the time range "
00875     << currTimeRange << "!" );
00876 #endif
00877 
00878   const Scalar tn = time_;
00879   const int kused = usedOrder_;
00880 
00881   // order of interpolation
00882   int kord = kused;
00883   if ( (kused == 0) || (timepoint == tn) )  {
00884     kord = 1;
00885   }
00886 
00887   // Initialize interploation
00888   Thyra::V_V(ptr(x_ptr),*xHistory_[0]);
00889   Thyra::V_S(ptr(xdot_ptr),ST::zero());
00890 
00891   // Add history array contributions
00892   const Scalar delt = timepoint - tn;
00893   Scalar c = ST::one(); // coefficient for interpolation of x
00894   Scalar d = ST::zero(); // coefficient for interpolation of xdot
00895   Scalar gam = delt/psi_[0]; // coefficient for interpolation
00896   for (int j=1 ; j <= kord ; ++j) {
00897     d = d*gam + c/psi_[j-1];
00898     c = c*gam;
00899     gam = (delt + psi_[j-1])/psi_[j];
00900     Thyra::Vp_StV(ptr(x_ptr),c,*xHistory_[j]);
00901     Thyra::Vp_StV(ptr(xdot_ptr),d,*xHistory_[j]);
00902   }
00903 
00904   // Set approximate accuracy
00905   if (accuracy_ptr) {
00906     *accuracy_ptr = Teuchos::ScalarTraits<Scalar>::pow(usedStep_,kord);
00907   }
00908 
00909 }
00910 
00911 
00912 template<class Scalar>
00913 void ImplicitBDFStepper<Scalar>::updateHistory_()
00914 {
00915 
00916   using Teuchos::as;
00917 
00918   // Save Newton correction for potential order increase on next step.
00919   if (usedOrder_ < maxOrder_)  {
00920     assign( xHistory_[usedOrder_+1].ptr(), *ee_ );
00921   }
00922   // Update history arrays
00923   Vp_V( xHistory_[usedOrder_].ptr(), *ee_ );
00924   for (int j=usedOrder_-1;j>=0;j--) {
00925     Vp_V( xHistory_[j].ptr(), *xHistory_[j+1] );
00926   }
00927   RCP<Teuchos::FancyOStream> out = this->getOStream();
00928   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00929   Teuchos::OSTab ostab(out,1,"updateHistory_");
00930   if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00931     for (int i=0;i<std::max(2,maxOrder_);++i) {
00932       *out << "xHistory_[" << i << "] = " << std::endl;
00933       xHistory_[i]->describe(*out,verbLevel);
00934     }
00935   }
00936 
00937 }
00938 
00939 
00940 template<class Scalar>
00941 void ImplicitBDFStepper<Scalar>::restoreHistory_()
00942 {
00943 
00944   using Teuchos::as;
00945   typedef Teuchos::ScalarTraits<Scalar> ST;
00946 
00947   // undo preparation of history array for prediction
00948   for (int i=nscsco_;i<=currentOrder_;++i) {
00949     Vt_S( xHistory_[i].ptr(), ST::one()/beta_[i] );
00950   }
00951   for (int i=1;i<=currentOrder_;++i) {
00952     psi_[i-1] = psi_[i] - hh_;
00953   }
00954   RCP<Teuchos::FancyOStream> out = this->getOStream();
00955   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00956   Teuchos::OSTab ostab(out,1,"restoreHistory_");
00957   if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00958     for (int i=0;i<maxOrder_;++i) {
00959       *out << "psi_[" << i << "] = " << psi_[i] << std::endl;
00960     }
00961     for (int i=0;i<maxOrder_;++i) {
00962       *out << "xHistory_[" << i << "] = " << std::endl;
00963       xHistory_[i]->describe(*out,verbLevel);
00964     }
00965   }
00966 
00967 }
00968 
00969 
00970 template<class Scalar>
00971 void ImplicitBDFStepper<Scalar>::updateCoeffs_()
00972 {
00973 
00974   using Teuchos::as;
00975   typedef Teuchos::ScalarTraits<Scalar> ST;
00976 
00977   // If the number of steps taken with constant order and constant stepsize is
00978   // more than the current order + 1 then we don't bother to update the
00979   // coefficients because we've reached a constant step-size formula.  When
00980   // this is is not true, then we update the coefficients for the variable
00981   // step-sizes.
00982   if ((hh_ != usedStep_) || (currentOrder_ != usedOrder_)) {
00983     nscsco_ = 0;
00984   }
00985   nscsco_ = std::min(nscsco_+1,usedOrder_+2);
00986   if (currentOrder_+1 >= nscsco_) {
00987     beta_[0] = ST::one();
00988     alpha_[0] = ST::one();
00989     Scalar temp1 = hh_;
00990     gamma_[0] = ST::zero();
00991     for (int i=1;i<=currentOrder_;++i) {
00992       Scalar temp2 = psi_[i-1];
00993       psi_[i-1] = temp1;
00994       beta_[i] = beta_[i-1]*psi_[i-1]/temp2;
00995       temp1 = temp2 + hh_;
00996       alpha_[i] = hh_/temp1;
00997       gamma_[i] = gamma_[i-1]+alpha_[i-1]/hh_;
00998     }
00999     psi_[currentOrder_] = temp1;
01000   }
01001   alpha_s_ = ST::zero();
01002   for (int i=0;i<currentOrder_;++i) {
01003     alpha_s_ = alpha_s_ - Scalar(ST::one()/(i+ST::one()));
01004   }
01005   RCP<Teuchos::FancyOStream> out = this->getOStream();
01006   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
01007   Teuchos::OSTab ostab(out,1,"updateCoeffs_");
01008   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
01009     for (int i=0;i<=maxOrder_;++i) {
01010       *out << "alpha_[" << i << "] = " << alpha_[i] << std::endl;
01011       *out << "beta_[" << i << "] = " << beta_[i] << std::endl;
01012       *out << "gamma_[" << i << "] = " << gamma_[i] << std::endl;
01013       *out << "psi_[" << i << "] = " << psi_[i] << std::endl;
01014       *out << "alpha_s_ = " << alpha_s_ << std::endl;
01015     }
01016     //std::cout << "alpha_s_ = " << alpha_s_ << std::endl;
01017     //for (int i=0;i<=maxOrder_;++i) {
01018     //  std::cout << "  alpha_[" << i << "] = " << alpha_[i] << std::endl;
01019     //  std::cout << "  beta_[" << i << "] = " << beta_[i] << std::endl;
01020     //  std::cout << "  gamma_[" << i << "] = " << gamma_[i] << std::endl;
01021     //  std::cout << "  psi_[" << i << "] = " << psi_[i] << std::endl;
01022     //}
01023   }
01024 }
01025 
01026 
01027 template<class Scalar>
01028 void ImplicitBDFStepper<Scalar>::initialize_()
01029 {
01030 
01031   using Teuchos::as;
01032   typedef Teuchos::ScalarTraits<Scalar> ST;
01033   using Thyra::createMember;
01034 
01035   RCP<Teuchos::FancyOStream> out = this->getOStream();
01036   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
01037   const bool doTrace = (as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH));
01038   Teuchos::OSTab ostab(out,1,"initialize_");
01039 
01040   if (doTrace) {
01041     *out
01042       << "\nEntering " << this->Teuchos::Describable::description()
01043       << "::initialize_()...\n";
01044   }
01045 
01046   TEUCHOS_TEST_FOR_EXCEPT(model_ == Teuchos::null);
01047   TEUCHOS_TEST_FOR_EXCEPT(solver_ == Teuchos::null);
01048   TEUCHOS_ASSERT(haveInitialCondition_);
01049 
01050   // Initialize Parameter List if none provided.
01051   if (parameterList_ == Teuchos::null) {
01052     RCP<Teuchos::ParameterList> emptyParameterList =
01053       Teuchos::rcp(new Teuchos::ParameterList);
01054     this->setParameterList(emptyParameterList);
01055   }
01056 
01057   // Initialize StepControl
01058   if (stepControl_ == Teuchos::null) {
01059     RCP<ImplicitBDFStepperStepControl<Scalar> > implicitBDFStepperStepControl =
01060       Teuchos::rcp(new ImplicitBDFStepperStepControl<Scalar>());
01061     RCP<Teuchos::ParameterList> stepControlPL =
01062       Teuchos::sublist(parameterList_, RythmosStepControlSettings_name);
01063     implicitBDFStepperStepControl->setParameterList(stepControlPL);
01064     this->setStepControlStrategy(implicitBDFStepperStepControl);
01065   }
01066   stepControl_->initialize(*this);
01067 
01068   maxOrder_ = stepControl_->getMaxOrder(); // maximum order
01069   TEUCHOS_TEST_FOR_EXCEPTION(
01070       !((1 <= maxOrder_) && (maxOrder_ <= 5)), std::logic_error,
01071       "Error, maxOrder returned from stepControl_->getMaxOrder() = "
01072       << maxOrder_ << " is outside range of [1,5]!\n");
01073 
01074   Scalar zero = ST::zero();
01075 
01076   currentOrder_ = 1; // Current order of integration
01077   usedOrder_ = 1;    // order used in current step (used after currentOrder_
01078                      // is updated)
01079   usedStep_ = zero;
01080   nscsco_  =  0;
01081   LETvalue_ = zero;
01082 
01083   alpha_.clear();  // $\alpha_j(n)=h_n/\psi_j(n)$ coefficient used in
01084                    // local error test
01085   // note:   $h_n$ = current step size, n = current time step
01086   gamma_.clear();  // calculate time derivative of history array for predictor
01087   beta_.clear();   // coefficients used to evaluate predictor from history array
01088   psi_.clear();    // $\psi_j(n) = t_n-t_{n-j}$ intermediary variable used to
01089   // compute $\beta_j(n):$
01090   for (int i=0 ; i<=maxOrder_ ; ++i) {
01091     alpha_.push_back(zero);
01092     beta_.push_back(zero);
01093     gamma_.push_back(zero);
01094     psi_.push_back(zero);
01095   }
01096   alpha_s_=Scalar(-ST::one());  // $\alpha_s$ fixed-leading coefficient of
01097                                 // this BDF method
01098   hh_=zero;
01099   numberOfSteps_=0;   // number of total time integration steps taken
01100   nef_ = 0;
01101 
01102   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
01103     *out << "alpha_s_ = " << alpha_s_ << std::endl;
01104     for (int i=0 ; i<=maxOrder_ ; ++i) {
01105       *out << "alpha_[" << i << "] = " << alpha_[i] << std::endl;
01106       *out << "beta_[" << i << "] = " << beta_[i] << std::endl;
01107       *out << "gamma_[" << i << "] = " << gamma_[i] << std::endl;
01108       *out << "psi_[" << i << "] = " << psi_[i] << std::endl;
01109     }
01110     *out << "numberOfSteps_ = " << numberOfSteps_ << std::endl;
01111   }
01112 
01113   // setInitialCondition initialized xHistory with xn0, xpn0.
01114   // Now we add the rest of the vectors.  Store maxOrder_+1 vectors
01115   for (int i=2 ; i<=maxOrder_ ; ++i) {
01116     xHistory_.push_back(createMember(xn0_->space()));
01117     V_S(xHistory_[i].ptr(),zero);
01118   }
01119 
01120   isInitialized_ = true;
01121 
01122   if (doTrace) {
01123     *out
01124       << "\nLeaving " << this->Teuchos::Describable::description()
01125       << "::initialize_()...\n";
01126   }
01127 
01128 }
01129 
01130 
01131 template<class Scalar>
01132 void ImplicitBDFStepper<Scalar>::completeStep_()
01133 {
01134 
01135   using Teuchos::as;
01136   typedef Teuchos::ScalarTraits<Scalar> ST;
01137 
01138 #ifdef HAVE_RYTHMOS_DEBUG
01139   TEUCHOS_TEST_FOR_EXCEPT(ST::isnaninf(hh_));
01140 #endif
01141 
01142   numberOfSteps_ ++;
01143   nef_ = 0;
01144   usedStep_ = hh_;
01145   usedOrder_ = currentOrder_;
01146   time_ += hh_;
01147   RCP<Teuchos::FancyOStream> out = this->getOStream();
01148   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
01149   Teuchos::OSTab ostab(out,1,"completeStep_");
01150 
01151   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
01152     *out << "numberOfSteps_ = " << numberOfSteps_ << std::endl;
01153     *out << "time_ = " << time_ << std::endl;
01154   }
01155 
01156   // 03/22/04 tscoffe:  Note that updating the history has nothing to do with
01157   // the step-size and everything to do with the newton correction vectors.
01158   updateHistory_();
01159 
01160 }
01161 
01162 template<class Scalar>
01163 void ImplicitBDFStepper<Scalar>::setStepControlData(
01164   const StepperBase<Scalar> & stepper)
01165 {
01166   if (!isInitialized_) {
01167     initialize_();
01168   }
01169   stepControl_->setStepControlData(stepper);
01170 }
01171 
01172 template<class Scalar>
01173 const Thyra::SolveStatus<Scalar>& ImplicitBDFStepper<Scalar>::getNonlinearSolveStatus() const
01174 {
01175   return nonlinearSolveStatus_;
01176 }
01177 
01178 //
01179 // Explicit Instantiation macro
01180 //
01181 // Must be expanded from within the Rythmos namespace!
01182 //
01183 
01184 #define RYTHMOS_IMPLICITBDF_STEPPER_INSTANT(SCALAR) \
01185   \
01186   template class ImplicitBDFStepper< SCALAR >; \
01187   \
01188   template RCP< ImplicitBDFStepper< SCALAR > > \
01189   implicitBDFStepper();  \
01190   \
01191   template RCP< ImplicitBDFStepper< SCALAR > > \
01192   implicitBDFStepper( \
01193     const RCP<Thyra::ModelEvaluator< SCALAR > >& model, \
01194     const RCP<Thyra::NonlinearSolverBase< SCALAR > >& solver, \
01195     const RCP<Teuchos::ParameterList>& parameterList \
01196     ); \
01197 
01198 
01199 } // namespace Rythmos
01200 
01201 
01202 #endif //Rythmos_IMPLICITBDF_STEPPER_DEF_H
 All Classes Functions Variables Typedefs Friends