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 
00403   numberOfSteps_ ++;
00404   nef_ = 0;
00405   time_ += hh_;
00406   RCP<Teuchos::FancyOStream> out = this->getOStream();
00407   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00408   Teuchos::OSTab ostab(out,1,"completeStep_");
00409 
00410   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00411     *out << "numberOfSteps_ = " << numberOfSteps_ << std::endl;
00412     *out << "nef_ = " << nef_ << std::endl;
00413     *out << "time_ = " << time_ << std::endl;
00414   }
00415   
00416   // Only update the time step if we are NOT running constant stepsize.
00417   bool adjustStep = (stepSizeType_ == STEP_TYPE_VARIABLE);
00418 
00419   Scalar newTimeStep = hh_;
00420   Scalar rr = ST::one(); // step size ratio = new step / old step
00421   // 03/11/04 tscoffe:  Here is the block for choosing order & step-size when
00422   // the local error test PASSES (and Newton succeeded). 
00423   int orderDiff = currentOrder_ - usedOrder_;
00424   usedOrder_ = currentOrder_;
00425   usedStep_ = hh_;
00426   if ((newOrder_ == currentOrder_-1) || (currentOrder_ == maxOrder_)) {
00427     // If we reduced our order or reached std::max order then move to the next phase
00428     // of integration where we don't automatically double the step-size and
00429     // increase the order.
00430     initialPhase_ = false;
00431   }
00432   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00433     *out << "initialPhase_ = " << initialPhase_ << std::endl;
00434   }
00435   if (initialPhase_) {
00436     currentOrder_++;
00437     newTimeStep = h_phase0_incr_ * hh_;
00438     if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00439       *out << "currentOrder_ = " << currentOrder_ << std::endl;
00440       *out << "newTimeStep = " << newTimeStep << std::endl;
00441     }
00442   } else { // not in the initial phase of integration
00443     BDFactionFlag action = ACTION_UNSET;
00444     if (newOrder_ == currentOrder_-1) {
00445       action = ACTION_LOWER;
00446     } else if (newOrder_ == maxOrder_) {
00447       action = ACTION_MAINTAIN;
00448     } else if ((currentOrder_+1>=nscsco_) || (orderDiff == 1)) {
00449       // If we just raised the order last time then we won't raise it again
00450       // until we've taken currentOrder_+1 steps at order currentOrder_.
00451       action = ACTION_MAINTAIN;
00452     } else { // consider changing the order 
00453       const ImplicitBDFStepper<Scalar>& implicitBDFStepper = Teuchos::dyn_cast<const ImplicitBDFStepper<Scalar> >(stepper);
00454       const Thyra::VectorBase<Scalar>& xHistory = implicitBDFStepper.getxHistory(currentOrder_+1);
00455       V_StVpStV(delta_.ptr(),ST::one(),*ee_,Scalar(-ST::one()),xHistory);
00456       Tkp1_ = wRMSNorm_(*errWtVec_,*delta_);
00457       Ekp1_ = Tkp1_/(currentOrder_+2);
00458       if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00459         *out << "delta_ = " << std::endl;
00460         delta_->describe(*out,verbLevel);
00461         *out << "Tkp1_ = ||delta_||_WRMS = " << Tkp1_ << std::endl;
00462         *out << "Ekp1_ = " << Ekp1_ << std::endl;
00463       }
00464       if (currentOrder_ == 1) {
00465         if (Tkp1_ >= Tkp1_Tk_safety_ * Tk_) {
00466           action = ACTION_MAINTAIN;
00467         } else {
00468           action = ACTION_RAISE;
00469         }
00470       } else {
00471         if (Tkm1_ <= std::min(Tk_,Tkp1_)) {
00472           action = ACTION_LOWER;
00473         } else if (Tkp1_ >= Tk_) {
00474           action = ACTION_MAINTAIN;
00475         } else {
00476           action = ACTION_RAISE;
00477         }
00478       }
00479     }
00480     if (currentOrder_ < minOrder_) {
00481       action = ACTION_RAISE;
00482     } else if ( (currentOrder_ == minOrder_) && (action == ACTION_LOWER) ) {
00483       action = ACTION_MAINTAIN;
00484     }
00485     if (action == ACTION_RAISE) {
00486       currentOrder_++;
00487       Est_ = Ekp1_;
00488     } else if (action == ACTION_LOWER) {
00489       currentOrder_--;
00490       Est_ = Ekm1_;
00491     }
00492     newTimeStep = hh_;
00493     rr = pow(r_safety_*Est_+r_fudge_,-1.0/(currentOrder_+1.0));
00494     if (rr >= r_hincr_test_) {
00495       rr = r_hincr_;
00496       newTimeStep = rr*hh_;
00497     } else if (rr <= 1) {
00498       rr = std::max(r_min_,std::min(r_max_,rr));
00499       newTimeStep = rr*hh_;
00500     }
00501     if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00502       *out << "Est_ = " << Est_ << std::endl;
00503       *out << "rr  = " << rr << std::endl;
00504       *out << "newTimeStep = " << newTimeStep << std::endl;
00505     }
00506   }
00507   
00508   if (time_ < stopTime_) {
00509     // If the step needs to be adjusted:
00510     if (adjustStep) {
00511       newTimeStep = std::max(newTimeStep, minTimeStep_);
00512       newTimeStep = std::min(newTimeStep, maxTimeStep_);
00513 
00514       Scalar nextTimePt = time_ + newTimeStep;
00515 
00516       if (nextTimePt > stopTime_) {
00517         nextTimePt  = stopTime_;
00518         newTimeStep = stopTime_ - time_;
00519       }
00520 
00521       hh_ = newTimeStep;
00522 
00523     } else { // if time step is constant for this step:
00524       Scalar nextTimePt = time_ + hh_;
00525 
00526       if (nextTimePt > stopTime_) {
00527         nextTimePt      = stopTime_;
00528         hh_ = stopTime_ - time_;
00529       }
00530     }
00531   }
00532   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00533     *out << "hh_ = " << hh_ << std::endl;
00534     *out << "currentOrder_ = " << currentOrder_ << std::endl;
00535   }
00536   setStepControlState_(READY_FOR_NEXT_STEP);
00537 }
00538 
00539 template<class Scalar>
00540 AttemptedStepStatusFlag ImplicitBDFStepperStepControl<Scalar>::rejectStep(const StepperBase<Scalar>& stepper)
00541 {
00542   TEST_FOR_EXCEPT(stepControlState_ != AFTER_CORRECTION);
00543 
00544   using Teuchos::as;
00545 
00546   // This routine puts its output in newTimeStep and newOrder
00547 
00548   // This routine changes the following variables:
00549   //    initialPhase, nef, psi, newTimeStep,
00550   //    newOrder, currentOrder_, currentTimeStep, dsDae.xHistory,
00551   //    dsDae.qHistory, nextTimePt, 
00552   //    currentTimeStepSum, nextTimePt
00553 
00554   // This routine reads but does not change the following variables:
00555   //    r_factor, r_safety, Est_, r_fudge_, r_min_, r_max_,
00556   //    minTimeStep_, maxTimeStep_, time, stopTime_ 
00557 
00558   // Only update the time step if we are NOT running constant stepsize.
00559   bool adjustStep = (stepSizeType_ == STEP_TYPE_VARIABLE);
00560 
00561   Scalar newTimeStep = hh_;
00562   Scalar rr = 1.0; // step size ratio = new step / old step
00563   nef_++;
00564   RCP<Teuchos::FancyOStream> out = this->getOStream();
00565   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00566   Teuchos::OSTab ostab(out,1,"rejectStep_");
00567   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00568     *out << "adjustStep = " << adjustStep << std::endl;
00569     *out << "nef_ = " << nef_ << std::endl;
00570   }
00571   if (nef_ >= max_LET_fail_)  {
00572     TEST_FOR_EXCEPTION(nef_ >= max_LET_fail_, std::logic_error, "Error, maximum number of local error test failures.\n");
00573   }
00574   initialPhase_ = false;
00575   if (adjustStep) {
00576     if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00577       *out << "initialPhase_ = " << initialPhase_ << std::endl;
00578     }
00579     for (int i=1;i<=currentOrder_;++i) {
00580       psi_[i-1] = psi_[i] - hh_;
00581     }
00582 
00583     if ((newtonConvergenceStatus_ < 0)) {
00585       // rely on the error estimate - it may be full of Nan's.
00586       rr = r_min_;
00587       newTimeStep = rr * hh_;
00588 
00589       if (nef_ > 2) {
00590         newOrder_ = 1;//consistent with block below.
00591       }
00592     } else {
00593       // 03/11/04 tscoffe:  Here is the block for choosing order & 
00594       // step-size when the local error test FAILS (but Newton 
00595       // succeeded). 
00596       if (nef_ == 1) { // first local error test failure
00597         rr = r_factor_*pow(r_safety_*Est_+r_fudge_,-1.0/(newOrder_+1.0));
00598         rr = std::max(r_min_,std::min(r_max_,rr));
00599         newTimeStep = rr * hh_;
00600       } else if (nef_ == 2) { // second failure
00601         rr = r_min_;
00602         newTimeStep = rr * hh_;
00603       } else if (nef_ > 2) { // third and later failures
00604         newOrder_ = 1;
00605         rr = r_min_;
00606         newTimeStep = rr * hh_;
00607       }
00608     }
00609     if (newOrder_ >= minOrder_) {
00610       currentOrder_ = newOrder_;
00611     }
00612     if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00613       *out << "rr = " << rr << std::endl;
00614       *out << "newOrder_ = " << newOrder_ << std::endl;
00615       *out << "currentOrder_ = " << currentOrder_ << std::endl;
00616     }
00617     if (numberOfSteps_ == 0) { // still first step
00618       psi_[0] = newTimeStep;
00619       if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00620         *out << "numberOfSteps_ == 0:" << std::endl;
00621         *out << "psi_[0] = " << psi_[0] << std::endl;
00622       }
00623     }
00624   } else if (!adjustStep) {
00625     if ( as<int>(verbLevel) != as<int>(Teuchos::VERB_NONE) ) {
00626       *out << "Rythmos_ImplicitBDFStepperStepControl::rejectStep(...):  "
00627           << "Warning:  Local error test failed with constant step-size."
00628           << std::endl;
00629     }
00630   }
00631 
00632   AttemptedStepStatusFlag return_status = PREDICT_AGAIN;
00633 
00634   // If the step needs to be adjusted:
00635   if (adjustStep) {
00636     newTimeStep = std::max(newTimeStep, minTimeStep_);
00637     newTimeStep = std::min(newTimeStep, maxTimeStep_);
00638 
00639     Scalar nextTimePt = time_ + newTimeStep;
00640 
00641     if (nextTimePt > stopTime_) {
00642       nextTimePt  = stopTime_;
00643       newTimeStep = stopTime_ - time_;
00644     }
00645 
00646     hh_ = newTimeStep;
00647   
00648   } else { // if time step is constant for this step:
00649     Scalar nextTimePt = time_ + hh_;
00650 
00651     if (nextTimePt > stopTime_) {
00652       nextTimePt      = stopTime_;
00653       hh_ = stopTime_ - time_;
00654     }
00655     return_status = CONTINUE_ANYWAY;
00656   }
00657   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00658     *out << "hh_ = " << hh_ << std::endl;
00659   }
00660 
00661   if (return_status == PREDICT_AGAIN) {
00662     setStepControlState_(READY_FOR_NEXT_STEP);
00663   } else if (return_status == CONTINUE_ANYWAY) {
00664     // do nothing, as we'll call completeStep which must be in AFTER_CORRECTION state.
00665   }
00666 
00667   return(return_status);
00668 }
00669 
00670 template<class Scalar>
00671 Scalar ImplicitBDFStepperStepControl<Scalar>::checkReduceOrder_(const StepperBase<Scalar>& stepper)
00672 {
00673   TEST_FOR_EXCEPT(stepControlState_ != AFTER_CORRECTION);
00674   TEST_FOR_EXCEPT(checkReduceOrderCalled_ == true);
00675 
00676   using Teuchos::as;
00677 
00678   const ImplicitBDFStepper<Scalar>& implicitBDFStepper = Teuchos::dyn_cast<const ImplicitBDFStepper<Scalar> >(stepper);
00679 
00680   // This routine puts its output in newOrder_
00681   
00682   RCP<Teuchos::FancyOStream> out = this->getOStream();
00683   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00684   Teuchos::OSTab ostab(out,1,"checkReduceOrder_");
00685   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00686     *out << "sigma_[" << currentOrder_ << "] = " << sigma_[currentOrder_] << std::endl;
00687     *out << "ee_ = " << std::endl;
00688     ee_->describe(*out,verbLevel);
00689     *out << "errWtVec_ = " << std::endl;
00690     errWtVec_->describe(*out,verbLevel);
00691   }
00692 
00693   Scalar enorm = wRMSNorm_(*errWtVec_,*ee_);
00694   Ek_ = sigma_[currentOrder_]*enorm;
00695   Tk_ = Scalar(currentOrder_+1)*Ek_;
00696   Est_ = Ek_;
00697   newOrder_ = currentOrder_;
00698   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00699     *out << "currentOrder_ = " << currentOrder_ << std::endl;
00700     *out << "Ek_ = " << Ek_ << std::endl;
00701     *out << "Tk_ = " << Tk_ << std::endl;
00702     *out << "enorm = " << enorm << std::endl;
00703   }
00704   if (currentOrder_>1) {
00705     const Thyra::VectorBase<Scalar>& xHistoryCur = implicitBDFStepper.getxHistory(currentOrder_);
00706     V_VpV(delta_.ptr(),xHistoryCur,*ee_);
00707     Ekm1_ = sigma_[currentOrder_-1]*wRMSNorm_(*errWtVec_,*delta_);
00708     Tkm1_ = currentOrder_*Ekm1_;
00709     if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00710       *out << "Ekm1_ = " << Ekm1_ << std::endl;
00711       *out << "Tkm1_ = " << Tkm1_ << std::endl;
00712     }
00713     if (currentOrder_>2) {
00714       const Thyra::VectorBase<Scalar>& xHistoryPrev = implicitBDFStepper.getxHistory(currentOrder_-1);
00715       Vp_V(delta_.ptr(),xHistoryPrev);
00716       Ekm2_ = sigma_[currentOrder_-2]*wRMSNorm_(*errWtVec_,*delta_);
00717       Tkm2_ = (currentOrder_-1)*Ekm2_;
00718       if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00719         *out << "Ekm2_ = " << Ekm2_ << std::endl;
00720         *out << "Tkm2_ = " << Tkm2_ << std::endl;
00721       }
00722       if (std::max(Tkm1_,Tkm2_)<=Tk_) {
00723         newOrder_--;
00724         Est_ = Ekm1_;
00725       }
00726     }
00727     else if (Tkm1_ <= Tkm1_Tk_safety_ * Tk_) {
00728       newOrder_--;
00729       Est_ = Ekm1_;
00730     }
00731   }
00732   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00733     *out << "Est_ = " << Est_ << std::endl;
00734     *out << "newOrder_= " << newOrder_ << std::endl;
00735   }
00736   checkReduceOrderCalled_ = true;
00737   return(enorm);
00738 }
00739 
00740 template<class Scalar>
00741 bool ImplicitBDFStepperStepControl<Scalar>::acceptStep(const StepperBase<Scalar>& stepper, Scalar* LETValue)
00742 {
00743   TEST_FOR_EXCEPT(stepControlState_ != AFTER_CORRECTION);
00744   typedef Teuchos::ScalarTraits<Scalar> ST;
00745   RCP<Teuchos::FancyOStream> out = this->getOStream();
00746   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00747   Teuchos::OSTab ostab(out,1,"acceptStep");
00748 
00749   bool return_status = false;
00750   Scalar enorm = checkReduceOrder_(stepper);
00751   Scalar LET = ck_*enorm;
00752   if (LETValue) {
00753     *LETValue = LET;
00754   }
00755   if (LET < ST::one()) {
00756     return_status = true;
00757   }
00758   if ( Teuchos::as<int>(verbLevel) >= Teuchos::as<int>(Teuchos::VERB_HIGH) ) {
00759     *out << "ck_ = " << ck_ << std::endl;
00760     *out << "enorm = " << enorm << std::endl;
00761     *out << "Local Truncation Error Check: (ck*enorm) < 1:  (" << LET << ") <?= 1" << std::endl;
00762   }
00763   return(return_status);
00764 }
00765 
00766 template<class Scalar>
00767 void ImplicitBDFStepperStepControl<Scalar>::describe(
00768   Teuchos::FancyOStream &out,
00769   const Teuchos::EVerbosityLevel verbLevel
00770   ) const
00771 {
00772 
00773   using Teuchos::as;
00774 
00775   if ( (as<int>(verbLevel) == as<int>(Teuchos::VERB_DEFAULT) ) ||
00776     (as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW)     )
00777     ) {
00778     out << this->description() << "::describe" << std::endl;
00779   }
00780   else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW)) {
00781     out << "time_ = " << time_ << std::endl;
00782     out << "hh_ = " << hh_ << std::endl;
00783     out << "currentOrder_ = " << currentOrder_ << std::endl;
00784   }
00785   else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_MEDIUM)) {
00786   }
00787   else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH)) {
00788     out << "ee_ = "; 
00789     if (ee_ == Teuchos::null) {
00790       out << "Teuchos::null" << std::endl;
00791     } else {
00792       ee_->describe(out,verbLevel);
00793     }
00794     out << "delta_ = ";
00795     if (delta_ == Teuchos::null) {
00796       out << "Teuchos::null" << std::endl;
00797     } else {
00798       delta_->describe(out,verbLevel);
00799     }
00800     out << "errWtVec_ = ";
00801     if (errWtVec_ == Teuchos::null) {
00802       out << "Teuchos::null" << std::endl;
00803     } else {
00804       errWtVec_->describe(out,verbLevel);
00805     }
00806   }
00807 }
00808 
00809 template<class Scalar>
00810 void ImplicitBDFStepperStepControl<Scalar>::setParameterList(
00811   RCP<Teuchos::ParameterList> const& paramList
00812   )
00813 {
00814 
00815   using Teuchos::as;
00816   typedef Teuchos::ScalarTraits<Scalar> ST;
00817 
00818   TEST_FOR_EXCEPT(paramList == Teuchos::null);
00819   paramList->validateParameters(*this->getValidParameters(),0);
00820   parameterList_ = paramList;
00821   Teuchos::readVerboseObjectSublist(&*parameterList_,this);
00822 
00823   minOrder_ = parameterList_->get("minOrder",int(1)); // minimum order
00824   TEST_FOR_EXCEPTION(
00825       !((1 <= minOrder_) && (minOrder_ <= 5)), std::logic_error,
00826       "Error, minOrder_ = " << minOrder_ << " is not in range [1,5]!\n"
00827       );
00828   maxOrder_ = parameterList_->get("maxOrder",int(5)); // maximum order
00829   TEST_FOR_EXCEPTION(
00830       !((1 <= maxOrder_) && (maxOrder_ <= 5)), std::logic_error,
00831       "Error, maxOrder_ = " << maxOrder_ << " is not in range [1,5]!\n"
00832       );
00833 
00834   relErrTol_ = parameterList_->get( "relErrTol", Scalar(1.0e-4) );
00835   absErrTol_ = parameterList_->get( "absErrTol", Scalar(1.0e-6) );
00836   bool constantStepSize = parameterList_->get( "constantStepSize", false );
00837   stopTime_ = parameterList_->get( "stopTime", Scalar(1.0) );
00838   
00839   if (constantStepSize == true) {
00840     stepSizeType_ = STEP_TYPE_FIXED;
00841   } else {
00842     stepSizeType_ = STEP_TYPE_VARIABLE;
00843   }
00844 
00845   RCP<Teuchos::FancyOStream> out = this->getOStream();
00846   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00847   Teuchos::OSTab ostab(out,1,"setParameterList");
00848   out->precision(15);
00849 
00850   setDefaultMagicNumbers_(parameterList_->sublist("magicNumbers"));
00851 
00852   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00853     *out << "minOrder_ = " << minOrder_ << std::endl;
00854     *out << "maxOrder_ = " << maxOrder_ << std::endl;
00855     *out << "relErrTol  = " << relErrTol_  << std::endl;
00856     *out << "absErrTol  = " << absErrTol_  << std::endl;
00857     *out << "stepSizeType = " << stepSizeType_  << std::endl;
00858     *out << "stopTime_  = " << stopTime_  << std::endl;
00859   }
00860 
00861 }
00862 
00863 template<class Scalar>
00864 void ImplicitBDFStepperStepControl<Scalar>::setDefaultMagicNumbers_(
00865   Teuchos::ParameterList &magicNumberList)
00866 {
00867 
00868   using Teuchos::as;
00869 
00870   // Magic Number Defaults:
00871   h0_safety_      = magicNumberList.get( "h0_safety",      Scalar(2.0)     );
00872   h0_max_factor_  = magicNumberList.get( "h0_max_factor",  Scalar(0.001)   );
00873   h_phase0_incr_  = magicNumberList.get( "h_phase0_incr",  Scalar(2.0)     );
00874   h_max_inv_      = magicNumberList.get( "h_max_inv",      Scalar(0.0)     );
00875   Tkm1_Tk_safety_ = magicNumberList.get( "Tkm1_Tk_safety", Scalar(2.0)     );
00876   Tkp1_Tk_safety_ = magicNumberList.get( "Tkp1_Tk_safety", Scalar(0.5)     );
00877   r_factor_       = magicNumberList.get( "r_factor",       Scalar(0.9)     );
00878   r_safety_       = magicNumberList.get( "r_safety",       Scalar(2.0)     );
00879   r_fudge_        = magicNumberList.get( "r_fudge",        Scalar(0.0001)  );
00880   r_min_          = magicNumberList.get( "r_min",          Scalar(0.125)   );
00881   r_max_          = magicNumberList.get( "r_max",          Scalar(0.9)     );
00882   r_hincr_test_   = magicNumberList.get( "r_hincr_test",   Scalar(2.0)     );
00883   r_hincr_        = magicNumberList.get( "r_hincr",        Scalar(2.0)     );
00884   max_LET_fail_   = magicNumberList.get( "max_LET_fail",   int(15)         );
00885   minTimeStep_    = magicNumberList.get( "minTimeStep",    Scalar(0.0)     );
00886   maxTimeStep_    = magicNumberList.get( "maxTimeStep",    Scalar(10.0)    ); 
00887 
00888   RCP<Teuchos::FancyOStream> out = this->getOStream();
00889   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00890   Teuchos::OSTab ostab(out,1,"setDefaultMagicNumbers_");
00891   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00892     *out << "h0_safety_ = " << h0_safety_ << std::endl;
00893     *out << "h0_max_factor_ = " << h0_max_factor_ << std::endl;
00894     *out << "h_phase0_incr_ = " << h_phase0_incr_ << std::endl;
00895     *out << "h_max_inv_ = " << h_max_inv_ << std::endl;
00896     *out << "Tkm1_Tk_safety_ = " << Tkm1_Tk_safety_  << std::endl;
00897     *out << "Tkp1_Tk_safety_ = " << Tkp1_Tk_safety_ << std::endl;
00898     *out << "r_factor_ = " << r_factor_ << std::endl;
00899     *out << "r_safety_ = " << r_safety_ << std::endl;
00900     *out << "r_fudge_ = " << r_fudge_ << std::endl;
00901     *out << "r_min_ = " << r_min_ << std::endl;
00902     *out << "r_max_ = " << r_max_ << std::endl;
00903     *out << "r_hincr_test_ = " << r_hincr_test_ << std::endl;
00904     *out << "r_hincr_ = " << r_hincr_ << std::endl;
00905     *out << "max_LET_fail_ = " << max_LET_fail_ << std::endl;
00906     *out << "minTimeStep_ = " << minTimeStep_ << std::endl;
00907     *out << "maxTimeStep_ = " << maxTimeStep_ << std::endl;
00908   }
00909 
00910 }
00911 
00912 template<class Scalar>
00913 RCP<const Teuchos::ParameterList>
00914 ImplicitBDFStepperStepControl<Scalar>::getValidParameters() const
00915 {
00916 
00917   static RCP<Teuchos::ParameterList> validPL;
00918 
00919   if (is_null(validPL)) {
00920 
00921     RCP<Teuchos::ParameterList>
00922       pl = Teuchos::parameterList();
00923 
00924     pl->set<int>   ( "minOrder",         1,
00925         "lower limit of order selection, guaranteed"
00926         );
00927     pl->set<int>   ( "maxOrder",         5,
00928         "upper limit of order selection, does not guarantee this order"        
00929         );
00930     pl->set<Scalar>( "relErrTol",        Scalar(1.0e-4) );
00931     pl->set<Scalar>( "absErrTol",        Scalar(1.0e-6) );
00932     pl->set<bool>  ( "constantStepSize", false          );
00933     pl->set<Scalar>( "stopTime",         Scalar(10.0)   );
00934 
00935     Teuchos::ParameterList
00936       &magicNumberList = pl->sublist("magicNumbers", 
00937           false,
00938           "These are knobs in the algorithm that have been set to reasonable values using lots of testing and heuristics and some theory."
00939           );
00940     magicNumberList.set<Scalar>( "h0_safety",      Scalar(2.0)     );
00941     magicNumberList.set<Scalar>( "h0_max_factor",  Scalar(0.001)   );
00942     magicNumberList.set<Scalar>( "h_phase0_incr",  Scalar(2.0),
00943         "initial ramp-up in variable mode (stepSize multiplier) "     
00944         );
00945     magicNumberList.set<Scalar>( "h_max_inv",      Scalar(0.0)     );
00946     magicNumberList.set<Scalar>( "Tkm1_Tk_safety", Scalar(2.0)     );
00947     magicNumberList.set<Scalar>( "Tkp1_Tk_safety", Scalar(0.5)     );
00948     magicNumberList.set<Scalar>( "r_factor",       Scalar(0.9),
00949         "used in rejectStep:  time step ratio multiplier"
00950         );
00951     magicNumberList.set<Scalar>( "r_safety",       Scalar(2.0),
00952         "local error multiplier as part of time step ratio calculation"
00953         );
00954     magicNumberList.set<Scalar>( "r_fudge",        Scalar(0.0001),
00955         "local error addition as part of time step ratio calculation"
00956         );
00957     magicNumberList.set<Scalar>( "r_min",          Scalar(0.125),
00958         "used in rejectStep:  how much to cut step and lower bound for time step ratio"   
00959         );
00960     magicNumberList.set<Scalar>( "r_max",          Scalar(0.9),
00961         "upper bound for time step ratio"
00962         );
00963     magicNumberList.set<Scalar>( "r_hincr_test",   Scalar(2.0),     
00964         "used in completeStep:  if time step ratio > this then set time step ratio to r_hincr"
00965         );
00966     magicNumberList.set<Scalar>( "r_hincr",        Scalar(2.0),
00967         "used in completeStep:  limit on time step ratio increases, not checked by r_max"
00968         );
00969     magicNumberList.set<int>   ( "max_LET_fail",   int(15),
00970         "Max number of rejected steps"
00971         );
00972     magicNumberList.set<Scalar>( "minTimeStep",    Scalar(0.0),
00973         "bound on smallest time step in variable mode."     
00974         );
00975     magicNumberList.set<Scalar>( "maxTimeStep",    Scalar(10.0),
00976         "bound on largest time step in variable mode."    
00977         ); 
00978 
00979     Teuchos::setupVerboseObjectSublist(&*pl);
00980 
00981     validPL = pl;
00982 
00983   }
00984 
00985   return (validPL);
00986   
00987 }
00988 
00989 template<class Scalar>
00990 RCP<Teuchos::ParameterList>
00991 ImplicitBDFStepperStepControl<Scalar>::unsetParameterList()
00992 {
00993   RCP<Teuchos::ParameterList> temp_param_list = parameterList_;
00994   parameterList_ = Teuchos::null;
00995   return(temp_param_list);
00996 }
00997 
00998 template<class Scalar>
00999 RCP<Teuchos::ParameterList>
01000 ImplicitBDFStepperStepControl<Scalar>::getNonconstParameterList()
01001 {
01002   return(parameterList_);
01003 }
01004 
01005 template<class Scalar>
01006 void ImplicitBDFStepperStepControl<Scalar>::setStepControlData(const StepperBase<Scalar>& stepper)
01007 {
01008   if (stepControlState_ == UNINITIALIZED) {
01009     initialize(stepper);
01010   }
01011   const ImplicitBDFStepper<Scalar>& bdfstepper = Teuchos::dyn_cast<const ImplicitBDFStepper<Scalar> >(stepper);
01012   int desiredOrder = bdfstepper.getOrder();
01013   TEST_FOR_EXCEPT(!((1 <= desiredOrder) && (desiredOrder <= maxOrder_)));
01014   if (stepControlState_ == BEFORE_FIRST_STEP) {
01015     TEST_FOR_EXCEPTION(
01016         desiredOrder > 1, 
01017         std::logic_error, 
01018         "Error, this ImplicitBDF stepper has not taken a step yet, so it cannot take a step of order " << desiredOrder << " > 1!\n"
01019         );
01020   }
01021   TEST_FOR_EXCEPT(!(desiredOrder <= currentOrder_+1));
01022   currentOrder_ = desiredOrder;
01023 
01024   using Teuchos::as;
01025   RCP<Teuchos::FancyOStream> out = this->getOStream();
01026   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
01027   Teuchos::OSTab ostab(out,1,"setStepControlData");
01028 
01029   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_EXTREME) ) {
01030     *out << "currentOrder_ = " << currentOrder_ << std::endl;
01031   }
01032 }
01033 
01034 template<class Scalar>
01035 bool ImplicitBDFStepperStepControl<Scalar>::supportsCloning() const
01036 {
01037   return true;
01038 }
01039 
01040 
01041 template<class Scalar>
01042 RCP<StepControlStrategyBase<Scalar> >
01043 ImplicitBDFStepperStepControl<Scalar>::cloneStepControlStrategyAlgorithm() const
01044 {
01045 
01046   RCP<ImplicitBDFStepperStepControl<Scalar> > stepControl = rcp(new ImplicitBDFStepperStepControl<Scalar>());
01047 
01048   if (!is_null(parameterList_)) {
01049     stepControl->setParameterList(parameterList_);
01050   }
01051 
01052   return stepControl;
01053 }
01054 
01055 template<class Scalar>
01056 void ImplicitBDFStepperStepControl<Scalar>::setErrWtVecCalc(const RCP<ErrWtVecCalcBase<Scalar> >& errWtVecCalc)
01057 {
01058   TEST_FOR_EXCEPT(is_null(errWtVecCalc));
01059   errWtVecCalc_ = errWtVecCalc;
01060 }
01061 
01062 template<class Scalar>
01063 RCP<const ErrWtVecCalcBase<Scalar> > ImplicitBDFStepperStepControl<Scalar>::getErrWtVecCalc() const
01064 {
01065   return(errWtVecCalc_);
01066 }
01067 
01068 template<class Scalar>
01069 Scalar ImplicitBDFStepperStepControl<Scalar>::wRMSNorm_(
01070     const Thyra::VectorBase<Scalar>& weight, 
01071     const Thyra::VectorBase<Scalar>& vector) const
01072 {
01073   return(norm_2(weight,vector));
01074 }
01075 
01076 template<class Scalar>
01077 void ImplicitBDFStepperStepControl<Scalar>::defaultInitializeAllData_()
01078 {
01079   typedef Teuchos::ScalarTraits<Scalar> ST;
01080   Scalar zero = ST::zero();
01081   Scalar mone = Scalar(-ST::one());
01082 
01083   stepControlState_ = UNINITIALIZED;
01084   hh_ = zero;
01085   numberOfSteps_ = 0;
01086   stepSizeType_ = STEP_TYPE_VARIABLE;
01087 
01088   minOrder_ = -1;
01089   maxOrder_ = -1;
01090   nef_ = 0;
01091   midStep_ = false;
01092   checkReduceOrderCalled_ = false;
01093   time_ = -std::numeric_limits<Scalar>::max();
01094   relErrTol_ = mone;
01095   absErrTol_ = mone;
01096   usedStep_ = mone;
01097   currentOrder_ = 1;
01098   usedOrder_ = -1;
01099   nscsco_ = -1;
01100   alpha_s_ = mone;
01101   alpha_0_ = mone;
01102   cj_ = mone;
01103   ck_ = mone;
01104   ck_enorm_ = mone;
01105   constantStepSize_ = false;
01106   Ek_ = mone;
01107   Ekm1_ = mone;
01108   Ekm2_ = mone;
01109   Ekp1_ = mone;
01110   Est_ = mone;
01111   Tk_ = mone;
01112   Tkm1_ = mone;
01113   Tkm2_ = mone;
01114   Tkp1_ = mone;
01115   newOrder_ = -1;
01116   oldOrder_ = -1;
01117   initialPhase_ = false;
01118   stopTime_ = mone;
01119   h0_safety_ = mone;
01120   h0_max_factor_ = mone;
01121   h_phase0_incr_ = mone;
01122   h_max_inv_ = mone;
01123   Tkm1_Tk_safety_ = mone;
01124   Tkp1_Tk_safety_ = mone;
01125   r_factor_ = mone;
01126   r_safety_ = mone;
01127   r_fudge_ = mone;
01128   r_min_ = mone;
01129   r_max_ = mone;
01130   r_hincr_test_ = mone;
01131   r_hincr_ = mone;
01132   max_LET_fail_ = -1;
01133   minTimeStep_ = mone;
01134   maxTimeStep_ = mone;
01135   newtonConvergenceStatus_ = -1;
01136 }
01137 
01138 template<class Scalar>
01139 int ImplicitBDFStepperStepControl<Scalar>::getMinOrder() const
01140 {
01141   TEST_FOR_EXCEPTION(
01142       stepControlState_ == UNINITIALIZED, std::logic_error,
01143       "Error, attempting to call getMinOrder before intiialization!\n"
01144       );
01145   return(minOrder_);
01146 }
01147 
01148 template<class Scalar>
01149 int ImplicitBDFStepperStepControl<Scalar>::getMaxOrder() const
01150 {
01151   TEST_FOR_EXCEPTION(
01152       stepControlState_ == UNINITIALIZED, std::logic_error,
01153       "Error, attempting to call getMaxOrder before initialization!\n"
01154       );
01155   return(maxOrder_);
01156 }
01157 
01158 // 
01159 // Explicit Instantiation macro
01160 //
01161 // Must be expanded from within the Rythmos namespace!
01162 //
01163 
01164 #define RYTHMOS_IMPLICITBDF_STEPPER_STEPCONTROL_INSTANT(SCALAR) \
01165   template class ImplicitBDFStepperStepControl< SCALAR >; 
01166 
01167 
01168 } // namespace Rythmos
01169 #endif // Rythmos_IMPLICITBDF_STEPPER_STEP_CONTROL_DEF_H
01170 
 All Classes Functions Variables Typedefs Friends