Rythmos - Transient Integration for Differential Equations Version of the Day
Rythmos_ImplicitBDFStepperStepControl_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_STEP_CONTROL_DEF_H
00030 #define Rythmos_IMPLICITBDF_STEPPER_STEP_CONTROL_DEF_H
00031 
00032 #include "Rythmos_ImplicitBDFStepperStepControl_decl.hpp"
00033 #include "Rythmos_ImplicitBDFStepper.hpp"
00034 #include "Rythmos_ImplicitBDFStepperErrWtVecCalc.hpp"
00035 
00036 namespace Rythmos {
00037 
00038 template<class Scalar>
00039 void ImplicitBDFStepperStepControl<Scalar>::setStepControlState_(StepControlStrategyState newState)
00040 {
00041   if (stepControlState_ == UNINITIALIZED) {
00042     TEST_FOR_EXCEPT(newState != BEFORE_FIRST_STEP);
00043   } else if (stepControlState_ == BEFORE_FIRST_STEP) {
00044     TEST_FOR_EXCEPT(newState != MID_STEP);
00045   } else if (stepControlState_ == MID_STEP) {
00046     TEST_FOR_EXCEPT(newState != AFTER_CORRECTION);
00047   } else if (stepControlState_ == AFTER_CORRECTION) {
00048     TEST_FOR_EXCEPT(newState != READY_FOR_NEXT_STEP);
00049     checkReduceOrderCalled_ = false;
00050   } else if (stepControlState_ == READY_FOR_NEXT_STEP) {
00051     TEST_FOR_EXCEPT(newState != MID_STEP);
00052   }
00053   stepControlState_ = newState;
00054 }
00055 
00056 template<class Scalar>
00057 StepControlStrategyState ImplicitBDFStepperStepControl<Scalar>::getCurrentState()
00058 {
00059   return(stepControlState_);
00060 }
00061 
00062 template<class Scalar>
00063 void ImplicitBDFStepperStepControl<Scalar>::updateCoeffs_() 
00064 {
00065   TEST_FOR_EXCEPT(!((stepControlState_ == BEFORE_FIRST_STEP) || (stepControlState_ == READY_FOR_NEXT_STEP)));
00066   using Teuchos::as;
00067   typedef Teuchos::ScalarTraits<Scalar> ST;
00068 
00069   // If the number of steps taken with constant order and constant stepsize is
00070   // more than the current order + 1 then we don't bother to update the
00071   // coefficients because we've reached a constant step-size formula.  When
00072   // this is is not true, then we update the coefficients for the variable
00073   // step-sizes. 
00074   if ((hh_ != usedStep_) || (currentOrder_ != usedOrder_)) {
00075     nscsco_ = 0;
00076   }
00077   nscsco_ = std::min(nscsco_+1,usedOrder_+2);
00078   if (currentOrder_+1 >= nscsco_) {
00079     alpha_[0] = ST::one();
00080     Scalar temp1 = hh_;
00081     sigma_[0] = ST::one();
00082     gamma_[0] = ST::zero();
00083     for (int i=1;i<=currentOrder_;++i) {
00084       Scalar temp2 = psi_[i-1];
00085       psi_[i-1] = temp1;
00086       temp1 = temp2 + hh_;
00087       alpha_[i] = hh_/temp1;
00088       sigma_[i] = Scalar(i+1)*sigma_[i-1]*alpha_[i];
00089       gamma_[i] = gamma_[i-1]+alpha_[i-1]/hh_;
00090     }
00091     psi_[currentOrder_] = temp1;
00092     alpha_s_ = ST::zero();
00093     alpha_0_ = ST::zero();
00094     for (int i=0;i<currentOrder_;++i) {
00095       alpha_s_ = alpha_s_ - Scalar(ST::one()/(i+ST::one()));
00096       alpha_0_ = alpha_0_ - alpha_[i];
00097     }
00098     cj_ = -alpha_s_/hh_;
00099     ck_ = std::abs(alpha_[currentOrder_]+alpha_s_-alpha_0_);
00100     ck_ = std::max(ck_,alpha_[currentOrder_]);
00101   }
00102   RCP<Teuchos::FancyOStream> out = this->getOStream();
00103   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00104   Teuchos::OSTab ostab(out,1,"updateCoeffs_");
00105   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00106     for (int i=0;i<=maxOrder_;++i) {
00107       *out << "alpha_[" << i << "] = " << alpha_[i] << std::endl;
00108       *out << "sigma_[" << i << "] = " << sigma_[i] << std::endl;
00109       *out << "gamma_[" << i << "] = " << gamma_[i] << std::endl;
00110       *out << "psi_[" << i << "] = " << psi_[i] << std::endl;
00111       *out << "alpha_s_ = " << alpha_s_ << std::endl;
00112       *out << "alpha_0_ = " << alpha_0_ << std::endl;
00113       *out << "cj_ = " << cj_ << std::endl;
00114       *out << "ck_ = " << ck_ << std::endl;
00115     }
00116   }
00117 }
00118 
00119 template<class Scalar>
00120 ImplicitBDFStepperStepControl<Scalar>::ImplicitBDFStepperStepControl()
00121 {
00122   defaultInitializeAllData_();
00123 }
00124 
00125 template<class Scalar>
00126 void ImplicitBDFStepperStepControl<Scalar>::initialize(const StepperBase<Scalar>& stepper)
00127 {
00128   using Teuchos::as;
00129   typedef Teuchos::ScalarTraits<Scalar> ST;
00130   using Thyra::createMember;
00131 
00132   RCP<Teuchos::FancyOStream> out = this->getOStream();
00133   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00134   const bool doTrace = (as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH));
00135   Teuchos::OSTab ostab(out,1,"initialize");
00136 
00137   if (doTrace) {
00138     *out
00139       << "\nEntering " << this->Teuchos::Describable::description()
00140       << "::initialize()...\n";
00141   }
00142 
00143   // Set initial time:
00144   TimeRange<Scalar> stepperRange = stepper.getTimeRange();
00145   TEST_FOR_EXCEPTION(
00146       !stepperRange.isValid(),
00147       std::logic_error,
00148       "Error, Stepper does not have valid time range for initialization of ImplicitBDFStepperStepControl!\n"
00149       );
00150   time_ = stepperRange.upper();
00151 
00152   if (parameterList_ == Teuchos::null) {
00153     RCP<Teuchos::ParameterList> emptyParameterList = Teuchos::rcp(new Teuchos::ParameterList);
00154     this->setParameterList(emptyParameterList);
00155   }
00156 
00157   if (is_null(errWtVecCalc_)) {
00158     RCP<ImplicitBDFStepperErrWtVecCalc<Scalar> > IBDFErrWtVecCalc = rcp(new ImplicitBDFStepperErrWtVecCalc<Scalar>());
00159     errWtVecCalc_ = IBDFErrWtVecCalc;
00160   }
00161 
00162   // 08/22/07 initialize can be called from the stepper when setInitialCondition is called.
00163   checkReduceOrderCalled_ = false;
00164   stepControlState_ = UNINITIALIZED;
00165 
00166   currentOrder_=1; // Current order of integration
00167   oldOrder_=1; // previous order of integration
00168   usedOrder_=1;  // order used in current step (used after currentOrder_ is updated)
00169   alpha_s_=Scalar(-ST::one());  // $\alpha_s$ fixed-leading coefficient of this BDF method
00170   alpha_.clear();  // $\alpha_j(n)=h_n/\psi_j(n)$ coefficient used in local error test
00171   // note:   $h_n$ = current step size, n = current time step
00172   gamma_.clear();  // calculate time derivative of history array for predictor 
00173   psi_.clear();    // $\psi_j(n) = t_n-t_{n-j}$ intermediary variable used to 
00174   sigma_.clear();  // $\sigma_j(n) = \frac{h_n^j(j-1)!}{\psi_1(n)*\cdots *\psi_j(n)}$
00175   Scalar zero = ST::zero();
00176   for (int i=0 ; i<=maxOrder_ ; ++i) {
00177     alpha_.push_back(zero);
00178     gamma_.push_back(zero);
00179     psi_.push_back(zero);
00180     sigma_.push_back(zero);
00181   }
00182   alpha_0_=zero;   // $-\sum_{j=1}^k \alpha_j(n)$ coefficient used in local error test
00183   cj_=zero;      // $-\alpha_s/h_n$ coefficient used in local error test
00184   ck_=zero;      // local error coefficient gamma_[0] = 0; // $\gamma_j(n)=\sum_{l=1}^{j-1}1/\psi_l(n)$ coefficient used to
00185   hh_=zero;
00186   numberOfSteps_=0;   // number of total time integration steps taken
00187   nef_=0;
00188   usedStep_=zero;
00189   nscsco_=0;
00190   Ek_=zero;
00191   Ekm1_=zero;
00192   Ekm2_=zero;
00193   Ekp1_=zero;
00194   Est_=zero;
00195   Tk_=zero;
00196   Tkm1_=zero;
00197   Tkm2_=zero;
00198   Tkp1_=zero;
00199   newOrder_=currentOrder_;
00200   initialPhase_=true;
00201 
00202   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00203     *out << "currentOrder_ = " << currentOrder_ << std::endl;
00204     *out << "oldOrder_ = " << oldOrder_ << std::endl;
00205     *out << "usedOrder_ = " << usedOrder_ << std::endl;
00206     *out << "alpha_s_ = " << alpha_s_ << std::endl;
00207     for (int i=0 ; i<=maxOrder_ ; ++i) {
00208       *out << "alpha_[" << i << "] = " << alpha_[i] << std::endl;
00209       *out << "gamma_[" << i << "] = " << gamma_[i] << std::endl;
00210       *out << "psi_[" << i << "] = " << psi_[i] << std::endl;
00211       *out << "sigma_[" << i << "] = " << sigma_[i] << std::endl;
00212     }
00213     *out << "alpha_0_ = " << alpha_0_ << std::endl;
00214     *out << "cj_ = " << cj_ << std::endl;
00215     *out << "ck_ = " << ck_ << std::endl;
00216     *out << "numberOfSteps_ = " << numberOfSteps_ << std::endl;
00217     *out << "nef_ = " << nef_ << std::endl;
00218     *out << "usedStep_ = " << usedStep_ << std::endl;
00219     *out << "nscsco_ = " << nscsco_ << std::endl;
00220     *out << "Ek_ = " << Ek_ << std::endl;
00221     *out << "Ekm1_ = " << Ekm1_ << std::endl;
00222     *out << "Ekm2_ = " << Ekm2_ << std::endl;
00223     *out << "Ekp1_ = " << Ekp1_ << std::endl;
00224     *out << "Est_ = " << Est_ << std::endl;
00225     *out << "Tk_ = " << Tk_ << std::endl;
00226     *out << "Tkm1_ = " << Tkm1_ << std::endl;
00227     *out << "Tkm2_ = " << Tkm2_ << std::endl;
00228     *out << "Tkp1_ = " << Tkp1_ << std::endl;
00229     *out << "newOrder_ = " << newOrder_ << std::endl;
00230     *out << "initialPhase_ = " << initialPhase_ << std::endl;
00231   }
00232 
00233 
00234   if (is_null(delta_)) {
00235     delta_ = createMember(stepper.get_x_space());
00236   }
00237   if (is_null(errWtVec_)) {
00238     errWtVec_ = createMember(stepper.get_x_space());
00239   }
00240   V_S(delta_.ptr(),zero); 
00241 
00242   setStepControlState_(BEFORE_FIRST_STEP);
00243 
00244   if (doTrace) {
00245     *out
00246       << "\nLeaving " << this->Teuchos::Describable::description()
00247       << "::initialize()...\n";
00248   }
00249 }
00250 
00251 template<class Scalar>
00252 void ImplicitBDFStepperStepControl<Scalar>::getFirstTimeStep_(const StepperBase<Scalar>& stepper)
00253 {
00254   
00255   TEST_FOR_EXCEPT(!(stepControlState_ == BEFORE_FIRST_STEP));
00256 
00257   using Teuchos::as;
00258   typedef Teuchos::ScalarTraits<Scalar> ST;
00259   Scalar zero = ST::zero();
00260 
00261   RCP<Teuchos::FancyOStream> out = this->getOStream();
00262   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00263   const bool doTrace = (as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH));
00264   Teuchos::OSTab ostab(out,1,"getFirstTimeStep_");
00265 
00266   const ImplicitBDFStepper<Scalar>& implicitBDFStepper = Teuchos::dyn_cast<const ImplicitBDFStepper<Scalar> >(stepper);
00267   const Thyra::VectorBase<Scalar>& xHistory0 = implicitBDFStepper.getxHistory(0);
00268   errWtVecCalc_->errWtVecSet(&*errWtVec_,xHistory0,relErrTol_,absErrTol_);
00269 
00270   // Choose initial step-size
00271   Scalar time_to_stop = stopTime_ - time_;
00272   Scalar currentTimeStep = ST::nan();
00273   if (stepSizeType_ == STEP_TYPE_FIXED) {
00274     currentTimeStep = hh_;
00275     //currentTimeStep = 0.1 * time_to_stop;
00276     //currentTimeStep = std::min(hh_, currentTimeStep);
00277   } else {
00278     // compute an initial step-size based on rate of change in the solution initially
00279     const Thyra::VectorBase<Scalar>& xHistory1 = implicitBDFStepper.getxHistory(1);
00280     Scalar ypnorm = wRMSNorm_(*errWtVec_,xHistory1);
00281     if (ypnorm > zero) { // time-dependent DAE
00282       currentTimeStep = std::min(h0_max_factor_*std::abs(time_to_stop),std::sqrt(2.0)/(h0_safety_*ypnorm));
00283     } else { // non-time-dependent DAE
00284       currentTimeStep = h0_max_factor_*std::abs(time_to_stop);
00285     }
00286     // choose std::min of user specified value and our value:
00287     if (hh_ > zero) {
00288       currentTimeStep = std::min(hh_, currentTimeStep);
00289     }
00290     // check for maximum step-size:
00291 #ifdef RYTHMOS_DEBUG
00292       TEST_FOR_EXCEPT(ST::isnaninf(currentTimeStep));
00293 #endif // RYTHMOS_DEBUG
00294     Scalar rh = std::abs(currentTimeStep)*h_max_inv_; 
00295     if (rh>1.0) {
00296       currentTimeStep = currentTimeStep/rh;
00297     }
00298   }
00299   hh_ = currentTimeStep;
00300   
00301   // Coefficient initialization 
00302   psi_[0] = hh_;
00303   cj_ = 1/psi_[0];
00304 
00305   if (doTrace) {
00306     *out << "\nhh_ = " << hh_ << std::endl;
00307     *out << "psi_[0] = " << psi_[0] << std::endl;
00308     *out << "cj_ = " << cj_ << std::endl;
00309   }
00310   
00311 }
00312 
00313 template<class Scalar>
00314 void ImplicitBDFStepperStepControl<Scalar>::setRequestedStepSize(
00315     const StepperBase<Scalar>& stepper
00316     ,const Scalar& stepSize
00317     ,const StepSizeType& stepSizeType
00318     )
00319 {
00320   typedef Teuchos::ScalarTraits<Scalar> ST;
00321   TEST_FOR_EXCEPT(!((stepControlState_ == UNINITIALIZED) || (stepControlState_ == BEFORE_FIRST_STEP) || (stepControlState_ == READY_FOR_NEXT_STEP) || (stepControlState_ == MID_STEP)));
00322   TEST_FOR_EXCEPTION(
00323       ((stepSizeType == STEP_TYPE_FIXED) && (stepSize == ST::zero())), 
00324       std::logic_error, 
00325       "Error, step size type == STEP_TYPE_FIXED, but requested step size == 0!\n"
00326       );
00327   bool didInitialization = false; 
00328   if (stepControlState_ == UNINITIALIZED) {
00329     initialize(stepper);
00330     didInitialization = true;
00331   }
00332 
00333   // errWtVecSet_ is called during initialize
00334   if (!didInitialization) {
00335     const ImplicitBDFStepper<Scalar>& implicitBDFStepper = Teuchos::dyn_cast<const ImplicitBDFStepper<Scalar> >(stepper);
00336     const Thyra::VectorBase<Scalar>& xHistory = implicitBDFStepper.getxHistory(0);
00337     errWtVecCalc_->errWtVecSet(&*errWtVec_,xHistory,relErrTol_,absErrTol_);
00338   }
00339 
00340   stepSizeType_ = stepSizeType;
00341   if (stepSizeType_ == STEP_TYPE_FIXED) {
00342     hh_ = stepSize;
00343     if (numberOfSteps_ == 0) {
00344       psi_[0] = hh_;
00345       cj_ = 1/psi_[0];
00346     }
00347   } else { // STEP_TYPE_VARIABLE
00348     if (stepSize != ST::zero()) {
00349       maxTimeStep_ = stepSize;
00350       h_max_inv_ = Scalar(ST::one()/maxTimeStep_);
00351     }
00352   }
00353 }
00354 
00355 template<class Scalar>
00356 void ImplicitBDFStepperStepControl<Scalar>::nextStepSize(const StepperBase<Scalar>& stepper, Scalar* stepSize, StepSizeType* stepSizeType, int* order)
00357 {
00358   TEST_FOR_EXCEPT(!((stepControlState_ == BEFORE_FIRST_STEP) || 
00359          (stepControlState_ == MID_STEP) ||  
00360          (stepControlState_ == READY_FOR_NEXT_STEP) )
00361         );
00362   if (stepControlState_ == BEFORE_FIRST_STEP) {
00363     getFirstTimeStep_(stepper);
00364   }
00365   if (stepControlState_ != MID_STEP) {
00366     this->updateCoeffs_();
00367   }
00368   if (stepSizeType_ == STEP_TYPE_VARIABLE) {
00369     if (hh_ > maxTimeStep_) {
00370       hh_ = maxTimeStep_;
00371     }
00372   }
00373   *stepSize = hh_;
00374   *stepSizeType = stepSizeType_;
00375   *order = currentOrder_;
00376   if (stepControlState_ != MID_STEP) {
00377     setStepControlState_(MID_STEP);
00378   }
00379 }
00380 
00381 template<class Scalar>
00382 void ImplicitBDFStepperStepControl<Scalar>::setCorrection(
00383      const StepperBase<Scalar>& stepper
00384     ,const RCP<const Thyra::VectorBase<Scalar> >& soln
00385     ,const RCP<const Thyra::VectorBase<Scalar> >& ee
00386     ,int solveStatus)
00387 {
00388   TEST_FOR_EXCEPT(stepControlState_ != MID_STEP);
00389   TEST_FOR_EXCEPTION(is_null(ee), std::logic_error, "Error, ee == Teuchos::null!\n");
00390   ee_ = ee;
00391   newtonConvergenceStatus_ = solveStatus;
00392   setStepControlState_(AFTER_CORRECTION);
00393 } 
00394 
00395 template<class Scalar>
00396 void ImplicitBDFStepperStepControl<Scalar>::completeStep(const StepperBase<Scalar>& stepper)
00397 {
00398   TEST_FOR_EXCEPT(stepControlState_ != AFTER_CORRECTION);
00399   using Teuchos::as;
00400   typedef Teuchos::ScalarTraits<Scalar> ST;
00401 
00402   numberOfSteps_ ++;
00403   nef_ = 0;
00404   time_ += hh_;
00405   RCP<Teuchos::FancyOStream> out = this->getOStream();
00406   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00407   Teuchos::OSTab ostab(out,1,"completeStep_");
00408 
00409   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00410     *out << "numberOfSteps_ = " << numberOfSteps_ << std::endl;
00411     *out << "nef_ = " << nef_ << std::endl;
00412     *out << "time_ = " << time_ << std::endl;
00413   }
00414   
00415   // Only update the time step if we are NOT running constant stepsize.
00416   bool adjustStep = (stepSizeType_ == STEP_TYPE_VARIABLE);
00417 
00418   Scalar newTimeStep = hh_;
00419   Scalar rr = ST::one(); // step size ratio = new step / old step
00420   // 03/11/04 tscoffe:  Here is the block for choosing order & step-size when
00421   // the local error test PASSES (and Newton succeeded). 
00422   int orderDiff = currentOrder_ - usedOrder_;
00423   usedOrder_ = currentOrder_;
00424   usedStep_ = hh_;
00425   if ((newOrder_ == currentOrder_-1) || (currentOrder_ == maxOrder_)) {
00426     // If we reduced our order or reached std::max order then move to the next phase
00427     // of integration where we don't automatically double the step-size and
00428     // increase the order.
00429     initialPhase_ = false;
00430   }
00431   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00432     *out << "initialPhase_ = " << initialPhase_ << std::endl;
00433   }
00434   if (initialPhase_) {
00435     currentOrder_++;
00436     newTimeStep = h_phase0_incr_ * hh_;
00437     if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00438       *out << "currentOrder_ = " << currentOrder_ << std::endl;
00439       *out << "newTimeStep = " << newTimeStep << std::endl;
00440     }
00441   } else { // not in the initial phase of integration
00442     BDFactionFlag action = ACTION_UNSET;
00443     if (newOrder_ == currentOrder_-1) {
00444       action = ACTION_LOWER;
00445     } else if (newOrder_ == maxOrder_) {
00446       action = ACTION_MAINTAIN;
00447     } else if ((currentOrder_+1>=nscsco_) || (orderDiff == 1)) {
00448       // If we just raised the order last time then we won't raise it again
00449       // until we've taken currentOrder_+1 steps at order currentOrder_.
00450       action = ACTION_MAINTAIN;
00451     } else { // consider changing the order 
00452       const ImplicitBDFStepper<Scalar>& implicitBDFStepper = Teuchos::dyn_cast<const ImplicitBDFStepper<Scalar> >(stepper);
00453       const Thyra::VectorBase<Scalar>& xHistory = implicitBDFStepper.getxHistory(currentOrder_+1);
00454       V_StVpStV(delta_.ptr(),ST::one(),*ee_,Scalar(-ST::one()),xHistory);
00455       Tkp1_ = wRMSNorm_(*errWtVec_,*delta_);
00456       Ekp1_ = Tkp1_/(currentOrder_+2);
00457       if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00458         *out << "delta_ = " << std::endl;
00459         delta_->describe(*out,verbLevel);
00460         *out << "Tkp1_ = ||delta_||_WRMS = " << Tkp1_ << std::endl;
00461         *out << "Ekp1_ = " << Ekp1_ << std::endl;
00462       }
00463       if (currentOrder_ == 1) {
00464         if (Tkp1_ >= Tkp1_Tk_safety_ * Tk_) {
00465           action = ACTION_MAINTAIN;
00466         } else {
00467           action = ACTION_RAISE;
00468         }
00469       } else {
00470         if (Tkm1_ <= std::min(Tk_,Tkp1_)) {
00471           action = ACTION_LOWER;
00472         } else if (Tkp1_ >= Tk_) {
00473           action = ACTION_MAINTAIN;
00474         } else {
00475           action = ACTION_RAISE;
00476         }
00477       }
00478     }
00479     if (currentOrder_ < minOrder_) {
00480       action = ACTION_RAISE;
00481     } else if ( (currentOrder_ == minOrder_) && (action == ACTION_LOWER) ) {
00482       action = ACTION_MAINTAIN;
00483     }
00484     if (action == ACTION_RAISE) {
00485       currentOrder_++;
00486       Est_ = Ekp1_;
00487     } else if (action == ACTION_LOWER) {
00488       currentOrder_--;
00489       Est_ = Ekm1_;
00490     }
00491     newTimeStep = hh_;
00492     rr = pow(r_safety_*Est_+r_fudge_,-1.0/(currentOrder_+1.0));
00493     if (rr >= r_hincr_test_) {
00494       rr = r_hincr_;
00495       newTimeStep = rr*hh_;
00496     } else if (rr <= 1) {
00497       rr = std::max(r_min_,std::min(r_max_,rr));
00498       newTimeStep = rr*hh_;
00499     }
00500     if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00501       *out << "Est_ = " << Est_ << std::endl;
00502       *out << "rr  = " << rr << std::endl;
00503       *out << "newTimeStep = " << newTimeStep << std::endl;
00504     }
00505   }
00506   
00507   if (time_ < stopTime_) {
00508     // If the step needs to be adjusted:
00509     if (adjustStep) {
00510       newTimeStep = std::max(newTimeStep, minTimeStep_);
00511       newTimeStep = std::min(newTimeStep, maxTimeStep_);
00512 
00513       Scalar nextTimePt = time_ + newTimeStep;
00514 
00515       if (nextTimePt > stopTime_) {
00516         nextTimePt  = stopTime_;
00517         newTimeStep = stopTime_ - time_;
00518       }
00519 
00520       hh_ = newTimeStep;
00521 
00522     } else { // if time step is constant for this step:
00523       Scalar nextTimePt = time_ + hh_;
00524 
00525       if (nextTimePt > stopTime_) {
00526         nextTimePt      = stopTime_;
00527         hh_ = stopTime_ - time_;
00528       }
00529     }
00530   }
00531   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00532     *out << "hh_ = " << hh_ << std::endl;
00533     *out << "currentOrder_ = " << currentOrder_ << std::endl;
00534   }
00535   setStepControlState_(READY_FOR_NEXT_STEP);
00536 }
00537 
00538 template<class Scalar>
00539 AttemptedStepStatusFlag ImplicitBDFStepperStepControl<Scalar>::rejectStep(const StepperBase<Scalar>& stepper)
00540 {
00541   TEST_FOR_EXCEPT(stepControlState_ != AFTER_CORRECTION);
00542 
00543   using Teuchos::as;
00544 
00545   // This routine puts its output in newTimeStep and newOrder
00546 
00547   // This routine changes the following variables:
00548   //    initialPhase, nef, psi, newTimeStep,
00549   //    newOrder, currentOrder_, currentTimeStep, dsDae.xHistory,
00550   //    dsDae.qHistory, nextTimePt, 
00551   //    currentTimeStepSum, nextTimePt
00552 
00553   // This routine reads but does not change the following variables:
00554   //    r_factor, r_safety, Est_, r_fudge_, r_min_, r_max_,
00555   //    minTimeStep_, maxTimeStep_, time, stopTime_ 
00556 
00557   // Only update the time step if we are NOT running constant stepsize.
00558   bool adjustStep = (stepSizeType_ == STEP_TYPE_VARIABLE);
00559 
00560   Scalar newTimeStep = hh_;
00561   Scalar rr = 1.0; // step size ratio = new step / old step
00562   nef_++;
00563   RCP<Teuchos::FancyOStream> out = this->getOStream();
00564   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00565   Teuchos::OSTab ostab(out,1,"rejectStep_");
00566   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00567     *out << "adjustStep = " << adjustStep << std::endl;
00568     *out << "nef_ = " << nef_ << std::endl;
00569   }
00570   if (nef_ >= max_LET_fail_)  {
00571     TEST_FOR_EXCEPTION(nef_ >= max_LET_fail_, std::logic_error, "Error, maximum number of local error test failures.\n");
00572   }
00573   initialPhase_ = false;
00574   if (adjustStep) {
00575     if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00576       *out << "initialPhase_ = " << initialPhase_ << std::endl;
00577     }
00578     for (int i=1;i<=currentOrder_;++i) {
00579       psi_[i-1] = psi_[i] - hh_;
00580     }
00581 
00582     if ((newtonConvergenceStatus_ < 0)) {
00584       // rely on the error estimate - it may be full of Nan's.
00585       rr = r_min_;
00586       newTimeStep = rr * hh_;
00587 
00588       if (nef_ > 2) {
00589         newOrder_ = 1;//consistent with block below.
00590       }
00591     } else {
00592       // 03/11/04 tscoffe:  Here is the block for choosing order & 
00593       // step-size when the local error test FAILS (but Newton 
00594       // succeeded). 
00595       if (nef_ == 1) { // first local error test failure
00596         rr = r_factor_*pow(r_safety_*Est_+r_fudge_,-1.0/(newOrder_+1.0));
00597         rr = std::max(r_min_,std::min(r_max_,rr));
00598         newTimeStep = rr * hh_;
00599       } else if (nef_ == 2) { // second failure
00600         rr = r_min_;
00601         newTimeStep = rr * hh_;
00602       } else if (nef_ > 2) { // third and later failures
00603         newOrder_ = 1;
00604         rr = r_min_;
00605         newTimeStep = rr * hh_;
00606       }
00607     }
00608     if (newOrder_ >= minOrder_) {
00609       currentOrder_ = newOrder_;
00610     }
00611     if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00612       *out << "rr = " << rr << std::endl;
00613       *out << "newOrder_ = " << newOrder_ << std::endl;
00614       *out << "currentOrder_ = " << currentOrder_ << std::endl;
00615     }
00616     if (numberOfSteps_ == 0) { // still first step
00617       psi_[0] = newTimeStep;
00618       if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00619         *out << "numberOfSteps_ == 0:" << std::endl;
00620         *out << "psi_[0] = " << psi_[0] << std::endl;
00621       }
00622     }
00623   } else if (!adjustStep) {
00624     if ( as<int>(verbLevel) != as<int>(Teuchos::VERB_NONE) ) {
00625       *out << "Rythmos_ImplicitBDFStepperStepControl::rejectStep(...):  "
00626           << "Warning:  Local error test failed with constant step-size."
00627           << std::endl;
00628     }
00629   }
00630 
00631   AttemptedStepStatusFlag return_status = PREDICT_AGAIN;
00632 
00633   // If the step needs to be adjusted:
00634   if (adjustStep) {
00635     newTimeStep = std::max(newTimeStep, minTimeStep_);
00636     newTimeStep = std::min(newTimeStep, maxTimeStep_);
00637 
00638     Scalar nextTimePt = time_ + newTimeStep;
00639 
00640     if (nextTimePt > stopTime_) {
00641       nextTimePt  = stopTime_;
00642       newTimeStep = stopTime_ - time_;
00643     }
00644 
00645     hh_ = newTimeStep;
00646   
00647   } else { // if time step is constant for this step:
00648     Scalar nextTimePt = time_ + hh_;
00649 
00650     if (nextTimePt > stopTime_) {
00651       nextTimePt      = stopTime_;
00652       hh_ = stopTime_ - time_;
00653     }
00654     return_status = CONTINUE_ANYWAY;
00655   }
00656   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00657     *out << "hh_ = " << hh_ << std::endl;
00658   }
00659 
00660   if (return_status == PREDICT_AGAIN) {
00661     setStepControlState_(READY_FOR_NEXT_STEP);
00662   } else if (return_status == CONTINUE_ANYWAY) {
00663     // do nothing, as we'll call completeStep which must be in AFTER_CORRECTION state.
00664   }
00665 
00666   return(return_status);
00667 }
00668 
00669 template<class Scalar>
00670 Scalar ImplicitBDFStepperStepControl<Scalar>::checkReduceOrder_(const StepperBase<Scalar>& stepper)
00671 {
00672   TEST_FOR_EXCEPT(stepControlState_ != AFTER_CORRECTION);
00673   TEST_FOR_EXCEPT(checkReduceOrderCalled_ == true);
00674 
00675   using Teuchos::as;
00676 
00677   const ImplicitBDFStepper<Scalar>& implicitBDFStepper = Teuchos::dyn_cast<const ImplicitBDFStepper<Scalar> >(stepper);
00678 
00679   // This routine puts its output in newOrder_
00680   
00681   RCP<Teuchos::FancyOStream> out = this->getOStream();
00682   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00683   Teuchos::OSTab ostab(out,1,"checkReduceOrder_");
00684   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00685     *out << "sigma_[" << currentOrder_ << "] = " << sigma_[currentOrder_] << std::endl;
00686     *out << "ee_ = " << std::endl;
00687     ee_->describe(*out,verbLevel);
00688     *out << "errWtVec_ = " << std::endl;
00689     errWtVec_->describe(*out,verbLevel);
00690   }
00691 
00692   Scalar enorm = wRMSNorm_(*errWtVec_,*ee_);
00693   Ek_ = sigma_[currentOrder_]*enorm;
00694   Tk_ = Scalar(currentOrder_+1)*Ek_;
00695   Est_ = Ek_;
00696   newOrder_ = currentOrder_;
00697   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00698     *out << "currentOrder_ = " << currentOrder_ << std::endl;
00699     *out << "Ek_ = " << Ek_ << std::endl;
00700     *out << "Tk_ = " << Tk_ << std::endl;
00701     *out << "enorm = " << enorm << std::endl;
00702   }
00703   if (currentOrder_>1) {
00704     const Thyra::VectorBase<Scalar>& xHistoryCur = implicitBDFStepper.getxHistory(currentOrder_);
00705     V_VpV(delta_.ptr(),xHistoryCur,*ee_);
00706     Ekm1_ = sigma_[currentOrder_-1]*wRMSNorm_(*errWtVec_,*delta_);
00707     Tkm1_ = currentOrder_*Ekm1_;
00708     if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00709       *out << "Ekm1_ = " << Ekm1_ << std::endl;
00710       *out << "Tkm1_ = " << Tkm1_ << std::endl;
00711     }
00712     if (currentOrder_>2) {
00713       const Thyra::VectorBase<Scalar>& xHistoryPrev = implicitBDFStepper.getxHistory(currentOrder_-1);
00714       Vp_V(delta_.ptr(),xHistoryPrev);
00715       Ekm2_ = sigma_[currentOrder_-2]*wRMSNorm_(*errWtVec_,*delta_);
00716       Tkm2_ = (currentOrder_-1)*Ekm2_;
00717       if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00718         *out << "Ekm2_ = " << Ekm2_ << std::endl;
00719         *out << "Tkm2_ = " << Tkm2_ << std::endl;
00720       }
00721       if (std::max(Tkm1_,Tkm2_)<=Tk_) {
00722         newOrder_--;
00723         Est_ = Ekm1_;
00724       }
00725     }
00726     else if (Tkm1_ <= Tkm1_Tk_safety_ * Tk_) {
00727       newOrder_--;
00728       Est_ = Ekm1_;
00729     }
00730   }
00731   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00732     *out << "Est_ = " << Est_ << std::endl;
00733     *out << "newOrder_= " << newOrder_ << std::endl;
00734   }
00735   checkReduceOrderCalled_ = true;
00736   return(enorm);
00737 }
00738 
00739 template<class Scalar>
00740 bool ImplicitBDFStepperStepControl<Scalar>::acceptStep(const StepperBase<Scalar>& stepper, Scalar* LETValue)
00741 {
00742   TEST_FOR_EXCEPT(stepControlState_ != AFTER_CORRECTION);
00743   typedef Teuchos::ScalarTraits<Scalar> ST;
00744   RCP<Teuchos::FancyOStream> out = this->getOStream();
00745   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00746   Teuchos::OSTab ostab(out,1,"acceptStep");
00747 
00748   bool return_status = false;
00749   Scalar enorm = checkReduceOrder_(stepper);
00750   Scalar LET = ck_*enorm;
00751 
00752   if (failStepIfNonlinearSolveFails_ && (newtonConvergenceStatus_ < 0) )
00753     return false;
00754 
00755   if (LETValue) {
00756     *LETValue = LET;
00757   }
00758   if (LET < ST::one()) {
00759     return_status = true;
00760   }
00761   if ( Teuchos::as<int>(verbLevel) >= Teuchos::as<int>(Teuchos::VERB_HIGH) ) {
00762     *out << "ck_ = " << ck_ << std::endl;
00763     *out << "enorm = " << enorm << std::endl;
00764     *out << "Local Truncation Error Check: (ck*enorm) < 1:  (" << LET << ") <?= 1" << std::endl;
00765   }
00766   return(return_status);
00767 }
00768 
00769 template<class Scalar>
00770 void ImplicitBDFStepperStepControl<Scalar>::describe(
00771   Teuchos::FancyOStream &out,
00772   const Teuchos::EVerbosityLevel verbLevel
00773   ) const
00774 {
00775 
00776   using Teuchos::as;
00777 
00778   if ( (as<int>(verbLevel) == as<int>(Teuchos::VERB_DEFAULT) ) ||
00779     (as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW)     )
00780     ) {
00781     out << this->description() << "::describe" << std::endl;
00782   }
00783   else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW)) {
00784     out << "time_ = " << time_ << std::endl;
00785     out << "hh_ = " << hh_ << std::endl;
00786     out << "currentOrder_ = " << currentOrder_ << std::endl;
00787   }
00788   else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_MEDIUM)) {
00789   }
00790   else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH)) {
00791     out << "ee_ = "; 
00792     if (ee_ == Teuchos::null) {
00793       out << "Teuchos::null" << std::endl;
00794     } else {
00795       ee_->describe(out,verbLevel);
00796     }
00797     out << "delta_ = ";
00798     if (delta_ == Teuchos::null) {
00799       out << "Teuchos::null" << std::endl;
00800     } else {
00801       delta_->describe(out,verbLevel);
00802     }
00803     out << "errWtVec_ = ";
00804     if (errWtVec_ == Teuchos::null) {
00805       out << "Teuchos::null" << std::endl;
00806     } else {
00807       errWtVec_->describe(out,verbLevel);
00808     }
00809   }
00810 }
00811 
00812 template<class Scalar>
00813 void ImplicitBDFStepperStepControl<Scalar>::setParameterList(
00814   RCP<Teuchos::ParameterList> const& paramList
00815   )
00816 {
00817 
00818   using Teuchos::as;
00819   typedef Teuchos::ScalarTraits<Scalar> ST;
00820 
00821   TEST_FOR_EXCEPT(paramList == Teuchos::null);
00822   paramList->validateParameters(*this->getValidParameters(),0);
00823   parameterList_ = paramList;
00824   Teuchos::readVerboseObjectSublist(&*parameterList_,this);
00825 
00826   minOrder_ = parameterList_->get("minOrder",int(1)); // minimum order
00827   TEST_FOR_EXCEPTION(
00828       !((1 <= minOrder_) && (minOrder_ <= 5)), std::logic_error,
00829       "Error, minOrder_ = " << minOrder_ << " is not in range [1,5]!\n"
00830       );
00831   maxOrder_ = parameterList_->get("maxOrder",int(5)); // maximum order
00832   TEST_FOR_EXCEPTION(
00833       !((1 <= maxOrder_) && (maxOrder_ <= 5)), std::logic_error,
00834       "Error, maxOrder_ = " << maxOrder_ << " is not in range [1,5]!\n"
00835       );
00836 
00837   relErrTol_ = parameterList_->get( "relErrTol", Scalar(1.0e-4) );
00838   absErrTol_ = parameterList_->get( "absErrTol", Scalar(1.0e-6) );
00839   bool constantStepSize = parameterList_->get( "constantStepSize", false );
00840   stopTime_ = parameterList_->get( "stopTime", Scalar(1.0) );
00841   
00842   if (constantStepSize == true) {
00843     stepSizeType_ = STEP_TYPE_FIXED;
00844   } else {
00845     stepSizeType_ = STEP_TYPE_VARIABLE;
00846   }
00847 
00848   failStepIfNonlinearSolveFails_ = 
00849     parameterList_->get( "failStepIfNonlinearSolveFails", false );
00850 
00851   RCP<Teuchos::FancyOStream> out = this->getOStream();
00852   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00853   Teuchos::OSTab ostab(out,1,"setParameterList");
00854   out->precision(15);
00855 
00856   setDefaultMagicNumbers_(parameterList_->sublist("magicNumbers"));
00857 
00858   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00859     *out << "minOrder_ = " << minOrder_ << std::endl;
00860     *out << "maxOrder_ = " << maxOrder_ << std::endl;
00861     *out << "relErrTol  = " << relErrTol_  << std::endl;
00862     *out << "absErrTol  = " << absErrTol_  << std::endl;
00863     *out << "stepSizeType = " << stepSizeType_  << std::endl;
00864     *out << "stopTime_  = " << stopTime_  << std::endl;
00865     *out << "failStepIfNonlinearSolveFails_ = " 
00866    << failStepIfNonlinearSolveFails_  << std::endl;
00867   }
00868 
00869 }
00870 
00871 template<class Scalar>
00872 void ImplicitBDFStepperStepControl<Scalar>::setDefaultMagicNumbers_(
00873   Teuchos::ParameterList &magicNumberList)
00874 {
00875 
00876   using Teuchos::as;
00877 
00878   // Magic Number Defaults:
00879   h0_safety_      = magicNumberList.get( "h0_safety",      Scalar(2.0)     );
00880   h0_max_factor_  = magicNumberList.get( "h0_max_factor",  Scalar(0.001)   );
00881   h_phase0_incr_  = magicNumberList.get( "h_phase0_incr",  Scalar(2.0)     );
00882   h_max_inv_      = magicNumberList.get( "h_max_inv",      Scalar(0.0)     );
00883   Tkm1_Tk_safety_ = magicNumberList.get( "Tkm1_Tk_safety", Scalar(2.0)     );
00884   Tkp1_Tk_safety_ = magicNumberList.get( "Tkp1_Tk_safety", Scalar(0.5)     );
00885   r_factor_       = magicNumberList.get( "r_factor",       Scalar(0.9)     );
00886   r_safety_       = magicNumberList.get( "r_safety",       Scalar(2.0)     );
00887   r_fudge_        = magicNumberList.get( "r_fudge",        Scalar(0.0001)  );
00888   r_min_          = magicNumberList.get( "r_min",          Scalar(0.125)   );
00889   r_max_          = magicNumberList.get( "r_max",          Scalar(0.9)     );
00890   r_hincr_test_   = magicNumberList.get( "r_hincr_test",   Scalar(2.0)     );
00891   r_hincr_        = magicNumberList.get( "r_hincr",        Scalar(2.0)     );
00892   max_LET_fail_   = magicNumberList.get( "max_LET_fail",   int(15)         );
00893   minTimeStep_    = magicNumberList.get( "minTimeStep",    Scalar(0.0)     );
00894   maxTimeStep_    = magicNumberList.get( "maxTimeStep",    Scalar(10.0)    ); 
00895 
00896   RCP<Teuchos::FancyOStream> out = this->getOStream();
00897   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00898   Teuchos::OSTab ostab(out,1,"setDefaultMagicNumbers_");
00899   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00900     *out << "h0_safety_ = " << h0_safety_ << std::endl;
00901     *out << "h0_max_factor_ = " << h0_max_factor_ << std::endl;
00902     *out << "h_phase0_incr_ = " << h_phase0_incr_ << std::endl;
00903     *out << "h_max_inv_ = " << h_max_inv_ << std::endl;
00904     *out << "Tkm1_Tk_safety_ = " << Tkm1_Tk_safety_  << std::endl;
00905     *out << "Tkp1_Tk_safety_ = " << Tkp1_Tk_safety_ << std::endl;
00906     *out << "r_factor_ = " << r_factor_ << std::endl;
00907     *out << "r_safety_ = " << r_safety_ << std::endl;
00908     *out << "r_fudge_ = " << r_fudge_ << std::endl;
00909     *out << "r_min_ = " << r_min_ << std::endl;
00910     *out << "r_max_ = " << r_max_ << std::endl;
00911     *out << "r_hincr_test_ = " << r_hincr_test_ << std::endl;
00912     *out << "r_hincr_ = " << r_hincr_ << std::endl;
00913     *out << "max_LET_fail_ = " << max_LET_fail_ << std::endl;
00914     *out << "minTimeStep_ = " << minTimeStep_ << std::endl;
00915     *out << "maxTimeStep_ = " << maxTimeStep_ << std::endl;
00916   }
00917 
00918 }
00919 
00920 template<class Scalar>
00921 RCP<const Teuchos::ParameterList>
00922 ImplicitBDFStepperStepControl<Scalar>::getValidParameters() const
00923 {
00924 
00925   static RCP<Teuchos::ParameterList> validPL;
00926 
00927   if (is_null(validPL)) {
00928 
00929     RCP<Teuchos::ParameterList>
00930       pl = Teuchos::parameterList();
00931 
00932     pl->set<int>   ( "minOrder",         1,
00933         "lower limit of order selection, guaranteed"
00934         );
00935     pl->set<int>   ( "maxOrder",         5,
00936         "upper limit of order selection, does not guarantee this order"        
00937         );
00938     pl->set<Scalar>( "relErrTol",        Scalar(1.0e-4) );
00939     pl->set<Scalar>( "absErrTol",        Scalar(1.0e-6) );
00940     pl->set<bool>  ( "constantStepSize", false          );
00941     pl->set<Scalar>( "stopTime",         Scalar(10.0)   );
00942     pl->set<bool>("failStepIfNonlinearSolveFails", false,
00943       "Power user command. Will force the function acceptStep() to return false ieven if the LET is acceptable.  Used to run with loose tolerances but enforce a correct nonlinear solution to the step.");
00944 
00945     Teuchos::ParameterList
00946       &magicNumberList = pl->sublist("magicNumbers", 
00947           false,
00948           "These are knobs in the algorithm that have been set to reasonable values using lots of testing and heuristics and some theory."
00949           );
00950     magicNumberList.set<Scalar>( "h0_safety",      Scalar(2.0)     );
00951     magicNumberList.set<Scalar>( "h0_max_factor",  Scalar(0.001)   );
00952     magicNumberList.set<Scalar>( "h_phase0_incr",  Scalar(2.0),
00953         "initial ramp-up in variable mode (stepSize multiplier) "     
00954         );
00955     magicNumberList.set<Scalar>( "h_max_inv",      Scalar(0.0)     );
00956     magicNumberList.set<Scalar>( "Tkm1_Tk_safety", Scalar(2.0)     );
00957     magicNumberList.set<Scalar>( "Tkp1_Tk_safety", Scalar(0.5)     );
00958     magicNumberList.set<Scalar>( "r_factor",       Scalar(0.9),
00959         "used in rejectStep:  time step ratio multiplier"
00960         );
00961     magicNumberList.set<Scalar>( "r_safety",       Scalar(2.0),
00962         "local error multiplier as part of time step ratio calculation"
00963         );
00964     magicNumberList.set<Scalar>( "r_fudge",        Scalar(0.0001),
00965         "local error addition as part of time step ratio calculation"
00966         );
00967     magicNumberList.set<Scalar>( "r_min",          Scalar(0.125),
00968         "used in rejectStep:  how much to cut step and lower bound for time step ratio"   
00969         );
00970     magicNumberList.set<Scalar>( "r_max",          Scalar(0.9),
00971         "upper bound for time step ratio"
00972         );
00973     magicNumberList.set<Scalar>( "r_hincr_test",   Scalar(2.0),     
00974         "used in completeStep:  if time step ratio > this then set time step ratio to r_hincr"
00975         );
00976     magicNumberList.set<Scalar>( "r_hincr",        Scalar(2.0),
00977         "used in completeStep:  limit on time step ratio increases, not checked by r_max"
00978         );
00979     magicNumberList.set<int>   ( "max_LET_fail",   int(15),
00980         "Max number of rejected steps"
00981         );
00982     magicNumberList.set<Scalar>( "minTimeStep",    Scalar(0.0),
00983         "bound on smallest time step in variable mode."     
00984         );
00985     magicNumberList.set<Scalar>( "maxTimeStep",    Scalar(10.0),
00986         "bound on largest time step in variable mode."    
00987         ); 
00988 
00989     Teuchos::setupVerboseObjectSublist(&*pl);
00990 
00991     validPL = pl;
00992 
00993   }
00994 
00995   return (validPL);
00996   
00997 }
00998 
00999 template<class Scalar>
01000 RCP<Teuchos::ParameterList>
01001 ImplicitBDFStepperStepControl<Scalar>::unsetParameterList()
01002 {
01003   RCP<Teuchos::ParameterList> temp_param_list = parameterList_;
01004   parameterList_ = Teuchos::null;
01005   return(temp_param_list);
01006 }
01007 
01008 template<class Scalar>
01009 RCP<Teuchos::ParameterList>
01010 ImplicitBDFStepperStepControl<Scalar>::getNonconstParameterList()
01011 {
01012   return(parameterList_);
01013 }
01014 
01015 template<class Scalar>
01016 void ImplicitBDFStepperStepControl<Scalar>::setStepControlData(const StepperBase<Scalar>& stepper)
01017 {
01018   if (stepControlState_ == UNINITIALIZED) {
01019     initialize(stepper);
01020   }
01021   const ImplicitBDFStepper<Scalar>& bdfstepper = Teuchos::dyn_cast<const ImplicitBDFStepper<Scalar> >(stepper);
01022   int desiredOrder = bdfstepper.getOrder();
01023   TEST_FOR_EXCEPT(!((1 <= desiredOrder) && (desiredOrder <= maxOrder_)));
01024   if (stepControlState_ == BEFORE_FIRST_STEP) {
01025     TEST_FOR_EXCEPTION(
01026         desiredOrder > 1, 
01027         std::logic_error, 
01028         "Error, this ImplicitBDF stepper has not taken a step yet, so it cannot take a step of order " << desiredOrder << " > 1!\n"
01029         );
01030   }
01031   TEST_FOR_EXCEPT(!(desiredOrder <= currentOrder_+1));
01032   currentOrder_ = desiredOrder;
01033 
01034   using Teuchos::as;
01035   RCP<Teuchos::FancyOStream> out = this->getOStream();
01036   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
01037   Teuchos::OSTab ostab(out,1,"setStepControlData");
01038 
01039   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_EXTREME) ) {
01040     *out << "currentOrder_ = " << currentOrder_ << std::endl;
01041   }
01042 }
01043 
01044 template<class Scalar>
01045 bool ImplicitBDFStepperStepControl<Scalar>::supportsCloning() const
01046 {
01047   return true;
01048 }
01049 
01050 
01051 template<class Scalar>
01052 RCP<StepControlStrategyBase<Scalar> >
01053 ImplicitBDFStepperStepControl<Scalar>::cloneStepControlStrategyAlgorithm() const
01054 {
01055 
01056   RCP<ImplicitBDFStepperStepControl<Scalar> > stepControl = rcp(new ImplicitBDFStepperStepControl<Scalar>());
01057 
01058   if (!is_null(parameterList_)) {
01059     stepControl->setParameterList(parameterList_);
01060   }
01061 
01062   return stepControl;
01063 }
01064 
01065 template<class Scalar>
01066 void ImplicitBDFStepperStepControl<Scalar>::setErrWtVecCalc(const RCP<ErrWtVecCalcBase<Scalar> >& errWtVecCalc)
01067 {
01068   TEST_FOR_EXCEPT(is_null(errWtVecCalc));
01069   errWtVecCalc_ = errWtVecCalc;
01070 }
01071 
01072 template<class Scalar>
01073 RCP<const ErrWtVecCalcBase<Scalar> > ImplicitBDFStepperStepControl<Scalar>::getErrWtVecCalc() const
01074 {
01075   return(errWtVecCalc_);
01076 }
01077 
01078 template<class Scalar>
01079 Scalar ImplicitBDFStepperStepControl<Scalar>::wRMSNorm_(
01080     const Thyra::VectorBase<Scalar>& weight, 
01081     const Thyra::VectorBase<Scalar>& vector) const
01082 {
01083   return(norm_2(weight,vector));
01084 }
01085 
01086 template<class Scalar>
01087 void ImplicitBDFStepperStepControl<Scalar>::defaultInitializeAllData_()
01088 {
01089   typedef Teuchos::ScalarTraits<Scalar> ST;
01090   Scalar zero = ST::zero();
01091   Scalar mone = Scalar(-ST::one());
01092 
01093   stepControlState_ = UNINITIALIZED;
01094   hh_ = zero;
01095   numberOfSteps_ = 0;
01096   stepSizeType_ = STEP_TYPE_VARIABLE;
01097 
01098   minOrder_ = -1;
01099   maxOrder_ = -1;
01100   nef_ = 0;
01101   midStep_ = false;
01102   checkReduceOrderCalled_ = false;
01103   time_ = -std::numeric_limits<Scalar>::max();
01104   relErrTol_ = mone;
01105   absErrTol_ = mone;
01106   usedStep_ = mone;
01107   currentOrder_ = 1;
01108   usedOrder_ = -1;
01109   nscsco_ = -1;
01110   alpha_s_ = mone;
01111   alpha_0_ = mone;
01112   cj_ = mone;
01113   ck_ = mone;
01114   ck_enorm_ = mone;
01115   constantStepSize_ = false;
01116   Ek_ = mone;
01117   Ekm1_ = mone;
01118   Ekm2_ = mone;
01119   Ekp1_ = mone;
01120   Est_ = mone;
01121   Tk_ = mone;
01122   Tkm1_ = mone;
01123   Tkm2_ = mone;
01124   Tkp1_ = mone;
01125   newOrder_ = -1;
01126   oldOrder_ = -1;
01127   initialPhase_ = false;
01128   stopTime_ = mone;
01129   h0_safety_ = mone;
01130   h0_max_factor_ = mone;
01131   h_phase0_incr_ = mone;
01132   h_max_inv_ = mone;
01133   Tkm1_Tk_safety_ = mone;
01134   Tkp1_Tk_safety_ = mone;
01135   r_factor_ = mone;
01136   r_safety_ = mone;
01137   r_fudge_ = mone;
01138   r_min_ = mone;
01139   r_max_ = mone;
01140   r_hincr_test_ = mone;
01141   r_hincr_ = mone;
01142   max_LET_fail_ = -1;
01143   minTimeStep_ = mone;
01144   maxTimeStep_ = mone;
01145   newtonConvergenceStatus_ = -1;
01146 }
01147 
01148 template<class Scalar>
01149 int ImplicitBDFStepperStepControl<Scalar>::getMinOrder() const
01150 {
01151   TEST_FOR_EXCEPTION(
01152       stepControlState_ == UNINITIALIZED, std::logic_error,
01153       "Error, attempting to call getMinOrder before intiialization!\n"
01154       );
01155   return(minOrder_);
01156 }
01157 
01158 template<class Scalar>
01159 int ImplicitBDFStepperStepControl<Scalar>::getMaxOrder() const
01160 {
01161   TEST_FOR_EXCEPTION(
01162       stepControlState_ == UNINITIALIZED, std::logic_error,
01163       "Error, attempting to call getMaxOrder before initialization!\n"
01164       );
01165   return(maxOrder_);
01166 }
01167 
01168 // 
01169 // Explicit Instantiation macro
01170 //
01171 // Must be expanded from within the Rythmos namespace!
01172 //
01173 
01174 #define RYTHMOS_IMPLICITBDF_STEPPER_STEPCONTROL_INSTANT(SCALAR) \
01175   template class ImplicitBDFStepperStepControl< SCALAR >; 
01176 
01177 
01178 } // namespace Rythmos
01179 #endif // Rythmos_IMPLICITBDF_STEPPER_STEP_CONTROL_DEF_H
01180 
 All Classes Functions Variables Typedefs Friends