Rythmos_ImplicitBDFStepperStepControl.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_H
00030 #define Rythmos_IMPLICITBDF_STEPPER_STEP_CONTROL_H
00031 
00032 #include "Rythmos_ImplicitBDFStepperStepControlDecl.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_,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>& xHistory = implicitBDFStepper.getxHistory(0);
00268   errWtVecCalc_->errWtVecSet(&*errWtVec_,xHistory,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>& xHistory = implicitBDFStepper.getxHistory(1);
00280     Scalar ypnorm = wRMSNorm_(*errWtVec_,xHistory);
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 TEUCHOS_DEBUG
00292       TEST_FOR_EXCEPT(ST::isnaninf(currentTimeStep));
00293 #endif
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_,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 (action == ACTION_RAISE) {
00481       currentOrder_++;
00482       Est_ = Ekp1_;
00483     } else if (action == ACTION_LOWER) {
00484       currentOrder_--;
00485       Est_ = Ekm1_;
00486     }
00487     newTimeStep = hh_;
00488     rr = pow(r_safety_*Est_+r_fudge_,-1.0/(currentOrder_+1.0));
00489     if (rr >= r_hincr_test_) {
00490       rr = r_hincr_;
00491       newTimeStep = rr*hh_;
00492     } else if (rr <= 1) {
00493       rr = std::max(r_min_,std::min(r_max_,rr));
00494       newTimeStep = rr*hh_;
00495     }
00496     if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00497       *out << "Est_ = " << Est_ << std::endl;
00498       *out << "rr  = " << rr << std::endl;
00499       *out << "newTimeStep = " << newTimeStep << std::endl;
00500     }
00501   }
00502   
00503   if (time_ < stopTime_) {
00504     // If the step needs to be adjusted:
00505     if (adjustStep) {
00506       newTimeStep = std::max(newTimeStep, minTimeStep_);
00507       newTimeStep = std::min(newTimeStep, maxTimeStep_);
00508 
00509       Scalar nextTimePt = time_ + newTimeStep;
00510 
00511       if (nextTimePt > stopTime_) {
00512         nextTimePt  = stopTime_;
00513         newTimeStep = stopTime_ - time_;
00514       }
00515 
00516       hh_ = newTimeStep;
00517 
00518     } else { // if time step is constant for this step:
00519       Scalar nextTimePt = time_ + hh_;
00520 
00521       if (nextTimePt > stopTime_) {
00522         nextTimePt      = stopTime_;
00523         hh_ = stopTime_ - time_;
00524       }
00525     }
00526   }
00527   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00528     *out << "hh_ = " << hh_ << std::endl;
00529     *out << "currentOrder_ = " << currentOrder_ << std::endl;
00530   }
00531   setStepControlState_(READY_FOR_NEXT_STEP);
00532 }
00533 
00534 template<class Scalar>
00535 AttemptedStepStatusFlag ImplicitBDFStepperStepControl<Scalar>::rejectStep(const StepperBase<Scalar>& stepper)
00536 {
00537   TEST_FOR_EXCEPT(stepControlState_ != AFTER_CORRECTION);
00538 
00539   using Teuchos::as;
00540 
00541   // This routine puts its output in newTimeStep and newOrder
00542 
00543   // This routine changes the following variables:
00544   //    initialPhase, nef, psi, newTimeStep,
00545   //    newOrder, currentOrder_, currentTimeStep, dsDae.xHistory,
00546   //    dsDae.qHistory, nextTimePt, 
00547   //    currentTimeStepSum, nextTimePt
00548 
00549   // This routine reads but does not change the following variables:
00550   //    r_factor, r_safety, Est_, r_fudge_, r_min_, r_max_,
00551   //    minTimeStep_, maxTimeStep_, time, stopTime_ 
00552 
00553   // Only update the time step if we are NOT running constant stepsize.
00554   bool adjustStep = (stepSizeType_ == STEP_TYPE_VARIABLE);
00555 
00556   Scalar newTimeStep = hh_;
00557   Scalar rr = 1.0; // step size ratio = new step / old step
00558   nef_++;
00559   RCP<Teuchos::FancyOStream> out = this->getOStream();
00560   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00561   Teuchos::OSTab ostab(out,1,"rejectStep_");
00562   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00563     *out << "adjustStep = " << adjustStep << std::endl;
00564     *out << "nef_ = " << nef_ << std::endl;
00565   }
00566   if (nef_ >= max_LET_fail_)  {
00567     TEST_FOR_EXCEPTION(nef_ >= max_LET_fail_, std::logic_error, "Error, maximum number of local error test failures.\n");
00568   }
00569   initialPhase_ = false;
00570   if (adjustStep) {
00571     if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00572       *out << "initialPhase_ = " << initialPhase_ << std::endl;
00573     }
00574     for (int i=1;i<=currentOrder_;++i) {
00575       psi_[i-1] = psi_[i] - hh_;
00576     }
00577 
00578     if ((newtonConvergenceStatus_ < 0)) {
00580       // rely on the error estimate - it may be full of Nan's.
00581       rr = r_min_;
00582       newTimeStep = rr * hh_;
00583 
00584       if (nef_ > 2) {
00585         newOrder_ = 1;//consistent with block below.
00586       }
00587     } else {
00588       // 03/11/04 tscoffe:  Here is the block for choosing order & 
00589       // step-size when the local error test FAILS (but Newton 
00590       // succeeded). 
00591       if (nef_ == 1) { // first local error test failure
00592         rr = r_factor_*pow(r_safety_*Est_+r_fudge_,-1.0/(newOrder_+1.0));
00593         rr = std::max(r_min_,std::min(r_max_,rr));
00594         newTimeStep = rr * hh_;
00595       } else if (nef_ == 2) { // second failure
00596         rr = r_min_;
00597         newTimeStep = rr * hh_;
00598       } else if (nef_ > 2) { // third and later failures
00599         newOrder_ = 1;
00600         rr = r_min_;
00601         newTimeStep = rr * hh_;
00602       }
00603     }
00604     if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00605       *out << "rr = " << rr << std::endl;
00606       *out << "newOrder_ = " << newOrder_ << std::endl;
00607     }
00608     currentOrder_ = newOrder_;
00609     if (numberOfSteps_ == 0) { // still first step
00610       psi_[0] = newTimeStep;
00611       if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00612         *out << "numberOfSteps_ == 0:" << std::endl;
00613         *out << "psi_[0] = " << psi_[0] << std::endl;
00614       }
00615     }
00616   } else if (!adjustStep) {
00617     *out << "Rythmos_ImplicitBDFStepperStepControl::rejectStep(...):  "
00618          << "Warning:  Local error test failed with constant step-size."
00619          << std::endl;
00620   }
00621 
00622   AttemptedStepStatusFlag return_status = PREDICT_AGAIN;
00623 
00624   // If the step needs to be adjusted:
00625   if (adjustStep) {
00626     newTimeStep = std::max(newTimeStep, minTimeStep_);
00627     newTimeStep = std::min(newTimeStep, maxTimeStep_);
00628 
00629     Scalar nextTimePt = time_ + newTimeStep;
00630 
00631     if (nextTimePt > stopTime_) {
00632       nextTimePt  = stopTime_;
00633       newTimeStep = stopTime_ - time_;
00634     }
00635 
00636     hh_ = newTimeStep;
00637   
00638   } else { // if time step is constant for this step:
00639     Scalar nextTimePt = time_ + hh_;
00640 
00641     if (nextTimePt > stopTime_) {
00642       nextTimePt      = stopTime_;
00643       hh_ = stopTime_ - time_;
00644     }
00645     return_status = CONTINUE_ANYWAY;
00646   }
00647   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00648     *out << "hh_ = " << hh_ << std::endl;
00649   }
00650 
00651   if (return_status == PREDICT_AGAIN) {
00652     setStepControlState_(READY_FOR_NEXT_STEP);
00653   } else if (return_status == CONTINUE_ANYWAY) {
00654     // do nothing, as we'll call completeStep which must be in AFTER_CORRECTION state.
00655   }
00656 
00657   return(return_status);
00658 }
00659 
00660 template<class Scalar>
00661 Scalar ImplicitBDFStepperStepControl<Scalar>::checkReduceOrder_(const StepperBase<Scalar>& stepper)
00662 {
00663   TEST_FOR_EXCEPT(stepControlState_ != AFTER_CORRECTION);
00664   TEST_FOR_EXCEPT(checkReduceOrderCalled_ == true);
00665 
00666   using Teuchos::as;
00667 
00668   const ImplicitBDFStepper<Scalar>& implicitBDFStepper = Teuchos::dyn_cast<const ImplicitBDFStepper<Scalar> >(stepper);
00669 
00670   // This routine puts its output in newOrder_
00671   
00672   RCP<Teuchos::FancyOStream> out = this->getOStream();
00673   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00674   Teuchos::OSTab ostab(out,1,"checkReduceOrder_");
00675   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00676     *out << "sigma_[" << currentOrder_ << "] = " << sigma_[currentOrder_] << std::endl;
00677     *out << "ee_ = " << std::endl;
00678     ee_->describe(*out,verbLevel);
00679     *out << "errWtVec_ = " << std::endl;
00680     errWtVec_->describe(*out,verbLevel);
00681   }
00682 
00683   Scalar enorm = wRMSNorm_(*errWtVec_,*ee_);
00684   Ek_ = sigma_[currentOrder_]*enorm;
00685   Tk_ = Scalar(currentOrder_+1)*Ek_;
00686   Est_ = Ek_;
00687   newOrder_ = currentOrder_;
00688   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00689     *out << "currentOrder_ = " << currentOrder_ << std::endl;
00690     *out << "Ek_ = " << Ek_ << std::endl;
00691     *out << "Tk_ = " << Tk_ << std::endl;
00692     *out << "enorm = " << enorm << std::endl;
00693   }
00694   if (currentOrder_>1) {
00695     const Thyra::VectorBase<Scalar>& xHistory = implicitBDFStepper.getxHistory(currentOrder_);
00696     V_VpV(&*delta_,xHistory,*ee_);
00697     Ekm1_ = sigma_[currentOrder_-1]*wRMSNorm_(*errWtVec_,*delta_);
00698     Tkm1_ = currentOrder_*Ekm1_;
00699     if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00700       *out << "Ekm1_ = " << Ekm1_ << std::endl;
00701       *out << "Tkm1_ = " << Tkm1_ << std::endl;
00702     }
00703     if (currentOrder_>2) {
00704       const Thyra::VectorBase<Scalar>& xHistory = implicitBDFStepper.getxHistory(currentOrder_-1);
00705       Vp_V(&*delta_,xHistory);
00706       Ekm2_ = sigma_[currentOrder_-2]*wRMSNorm_(*errWtVec_,*delta_);
00707       Tkm2_ = (currentOrder_-1)*Ekm2_;
00708       if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00709         *out << "Ekm2_ = " << Ekm2_ << std::endl;
00710         *out << "Tkm2_ = " << Tkm2_ << std::endl;
00711       }
00712       if (std::max(Tkm1_,Tkm2_)<=Tk_) {
00713         newOrder_--;
00714         Est_ = Ekm1_;
00715       }
00716     }
00717     else if (Tkm1_ <= Tkm1_Tk_safety_ * Tk_) {
00718       newOrder_--;
00719       Est_ = Ekm1_;
00720     }
00721   }
00722   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00723     *out << "Est_ = " << Est_ << std::endl;
00724     *out << "newOrder_= " << newOrder_ << std::endl;
00725   }
00726   checkReduceOrderCalled_ = true;
00727   return(enorm);
00728 }
00729 
00730 template<class Scalar>
00731 bool ImplicitBDFStepperStepControl<Scalar>::acceptStep(const StepperBase<Scalar>& stepper, Scalar* LETValue)
00732 {
00733   TEST_FOR_EXCEPT(stepControlState_ != AFTER_CORRECTION);
00734   typedef Teuchos::ScalarTraits<Scalar> ST;
00735   RCP<Teuchos::FancyOStream> out = this->getOStream();
00736   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00737   Teuchos::OSTab ostab(out,1,"acceptStep");
00738 
00739   bool return_status = false;
00740   Scalar enorm = checkReduceOrder_(stepper);
00741   Scalar LET = ck_*enorm;
00742   if (LETValue) {
00743     *LETValue = LET;
00744   }
00745   if (LET < ST::one()) {
00746     return_status = true;
00747   }
00748   if ( Teuchos::as<int>(verbLevel) >= Teuchos::as<int>(Teuchos::VERB_HIGH) ) {
00749     *out << "ck_ = " << ck_ << std::endl;
00750     *out << "enorm = " << enorm << std::endl;
00751     *out << "Local Truncation Error Check: (ck*enorm) < 1:  (" << LET << ") <?= 1" << std::endl;
00752   }
00753   return(return_status);
00754 }
00755 
00756 template<class Scalar>
00757 void ImplicitBDFStepperStepControl<Scalar>::describe(
00758   Teuchos::FancyOStream &out,
00759   const Teuchos::EVerbosityLevel verbLevel
00760   ) const
00761 {
00762 
00763   using Teuchos::as;
00764 
00765   if ( (as<int>(verbLevel) == as<int>(Teuchos::VERB_DEFAULT) ) ||
00766     (as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW)     )
00767     ) {
00768     out << this->description() << "::describe" << std::endl;
00769   }
00770   else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW)) {
00771     out << "time_ = " << time_ << std::endl;
00772     out << "hh_ = " << hh_ << std::endl;
00773     out << "currentOrder_ = " << currentOrder_ << std::endl;
00774   }
00775   else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_MEDIUM)) {
00776   }
00777   else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH)) {
00778     out << "ee_ = "; 
00779     if (ee_ == Teuchos::null) {
00780       out << "Teuchos::null" << std::endl;
00781     } else {
00782       ee_->describe(out,verbLevel);
00783     }
00784     out << "delta_ = ";
00785     if (delta_ == Teuchos::null) {
00786       out << "Teuchos::null" << std::endl;
00787     } else {
00788       delta_->describe(out,verbLevel);
00789     }
00790     out << "errWtVec_ = ";
00791     if (errWtVec_ == Teuchos::null) {
00792       out << "Teuchos::null" << std::endl;
00793     } else {
00794       errWtVec_->describe(out,verbLevel);
00795     }
00796   }
00797 }
00798 
00799 template<class Scalar>
00800 void ImplicitBDFStepperStepControl<Scalar>::setParameterList(
00801   RCP<Teuchos::ParameterList> const& paramList
00802   )
00803 {
00804 
00805   using Teuchos::as;
00806   typedef Teuchos::ScalarTraits<Scalar> ST;
00807 
00808   TEST_FOR_EXCEPT(paramList == Teuchos::null);
00809   paramList->validateParameters(*this->getValidParameters(),0);
00810   parameterList_ = paramList;
00811   Teuchos::readVerboseObjectSublist(&*parameterList_,this);
00812 
00813   maxOrder_ = parameterList_->get("maxOrder",int(5)); // maximum order
00814   TEST_FOR_EXCEPTION(
00815       !((1 <= maxOrder_) && (maxOrder_ <= 5)), std::logic_error,
00816       "Error, maxOrder_ = " << maxOrder_ << " is not in range [1,5]!\n"
00817       );
00818 
00819   relErrTol_ = parameterList_->get( "relErrTol", Scalar(1.0e-4) );
00820   absErrTol_ = parameterList_->get( "absErrTol", Scalar(1.0e-6) );
00821   bool constantStepSize = parameterList_->get( "constantStepSize", false );
00822   stopTime_ = parameterList_->get( "stopTime", Scalar(1.0) );
00823   
00824   if (constantStepSize == true) {
00825     stepSizeType_ = STEP_TYPE_FIXED;
00826   } else {
00827     stepSizeType_ = STEP_TYPE_VARIABLE;
00828   }
00829 
00830   RCP<Teuchos::FancyOStream> out = this->getOStream();
00831   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00832   Teuchos::OSTab ostab(out,1,"setParameterList");
00833   out->precision(15);
00834 
00835   setDefaultMagicNumbers_(parameterList_->sublist("magicNumbers"));
00836 
00837   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00838     *out << "maxOrder_ = " << maxOrder_ << std::endl;
00839     *out << "relErrTol  = " << relErrTol_  << std::endl;
00840     *out << "absErrTol  = " << absErrTol_  << std::endl;
00841     *out << "stepSizeType = " << stepSizeType_  << std::endl;
00842     *out << "stopTime_  = " << stopTime_  << std::endl;
00843   }
00844 
00845 }
00846 
00847 template<class Scalar>
00848 void ImplicitBDFStepperStepControl<Scalar>::setDefaultMagicNumbers_(
00849   Teuchos::ParameterList &magicNumberList)
00850 {
00851 
00852   using Teuchos::as;
00853 
00854   // Magic Number Defaults:
00855   h0_safety_      = magicNumberList.get( "h0_safety",      Scalar(2.0)     );
00856   h0_max_factor_  = magicNumberList.get( "h0_max_factor",  Scalar(0.001)   );
00857   h_phase0_incr_  = magicNumberList.get( "h_phase0_incr",  Scalar(2.0)     );
00858   h_max_inv_      = magicNumberList.get( "h_max_inv",      Scalar(0.0)     );
00859   Tkm1_Tk_safety_ = magicNumberList.get( "Tkm1_Tk_safety", Scalar(2.0)     );
00860   Tkp1_Tk_safety_ = magicNumberList.get( "Tkp1_Tk_safety", Scalar(0.5)     );
00861   r_factor_       = magicNumberList.get( "r_factor",       Scalar(0.9)     );
00862   r_safety_       = magicNumberList.get( "r_safety",       Scalar(2.0)     );
00863   r_fudge_        = magicNumberList.get( "r_fudge",        Scalar(0.0001)  );
00864   r_min_          = magicNumberList.get( "r_min",          Scalar(0.125)   );
00865   r_max_          = magicNumberList.get( "r_max",          Scalar(0.9)     );
00866   r_hincr_test_   = magicNumberList.get( "r_hincr_test",   Scalar(2.0)     );
00867   r_hincr_        = magicNumberList.get( "r_hincr",        Scalar(2.0)     );
00868   max_LET_fail_   = magicNumberList.get( "max_LET_fail",   int(15)         );
00869   minTimeStep_    = magicNumberList.get( "minTimeStep",    Scalar(0.0)     );
00870   maxTimeStep_    = magicNumberList.get( "maxTimeStep",    Scalar(10.0)    ); 
00871 
00872   RCP<Teuchos::FancyOStream> out = this->getOStream();
00873   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00874   Teuchos::OSTab ostab(out,1,"setDefaultMagicNumbers_");
00875   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00876     *out << "h0_safety_ = " << h0_safety_ << std::endl;
00877     *out << "h0_max_factor_ = " << h0_max_factor_ << std::endl;
00878     *out << "h_phase0_incr_ = " << h_phase0_incr_ << std::endl;
00879     *out << "h_max_inv_ = " << h_max_inv_ << std::endl;
00880     *out << "Tkm1_Tk_safety_ = " << Tkm1_Tk_safety_  << std::endl;
00881     *out << "Tkp1_Tk_safety_ = " << Tkp1_Tk_safety_ << std::endl;
00882     *out << "r_factor_ = " << r_factor_ << std::endl;
00883     *out << "r_safety_ = " << r_safety_ << std::endl;
00884     *out << "r_fudge_ = " << r_fudge_ << std::endl;
00885     *out << "r_min_ = " << r_min_ << std::endl;
00886     *out << "r_max_ = " << r_max_ << std::endl;
00887     *out << "r_hincr_test_ = " << r_hincr_test_ << std::endl;
00888     *out << "r_hincr_ = " << r_hincr_ << std::endl;
00889     *out << "max_LET_fail_ = " << max_LET_fail_ << std::endl;
00890     *out << "minTimeStep_ = " << minTimeStep_ << std::endl;
00891     *out << "maxTimeStep_ = " << maxTimeStep_ << std::endl;
00892   }
00893 
00894 }
00895 
00896 template<class Scalar>
00897 RCP<const Teuchos::ParameterList>
00898 ImplicitBDFStepperStepControl<Scalar>::getValidParameters() const
00899 {
00900 
00901   static RCP<Teuchos::ParameterList> validPL;
00902 
00903   if (is_null(validPL)) {
00904 
00905     RCP<Teuchos::ParameterList>
00906       pl = Teuchos::parameterList();
00907 
00908     pl->set<int>   ( "maxOrder",         5              );
00909     pl->set<Scalar>( "relErrTol",        Scalar(1.0e-4) );
00910     pl->set<Scalar>( "absErrTol",        Scalar(1.0e-6) );
00911     pl->set<bool>  ( "constantStepSize", false          );
00912     pl->set<Scalar>( "stopTime",         Scalar(10.0)   );
00913 
00914     Teuchos::ParameterList
00915       &magicNumberList = pl->sublist("magicNumbers");
00916     magicNumberList.set<Scalar>( "h0_safety",      Scalar(2.0)     );
00917     magicNumberList.set<Scalar>( "h0_max_factor",  Scalar(0.001)   );
00918     magicNumberList.set<Scalar>( "h_phase0_incr",  Scalar(2.0)     );
00919     magicNumberList.set<Scalar>( "h_max_inv",      Scalar(0.0)     );
00920     magicNumberList.set<Scalar>( "Tkm1_Tk_safety", Scalar(2.0)     );
00921     magicNumberList.set<Scalar>( "Tkp1_Tk_safety", Scalar(0.5)     );
00922     magicNumberList.set<Scalar>( "r_factor",       Scalar(0.9)     );
00923     magicNumberList.set<Scalar>( "r_safety",       Scalar(2.0)     );
00924     magicNumberList.set<Scalar>( "r_fudge",        Scalar(0.0001)  );
00925     magicNumberList.set<Scalar>( "r_min",          Scalar(0.125)   );
00926     magicNumberList.set<Scalar>( "r_max",          Scalar(0.9)     );
00927     magicNumberList.set<Scalar>( "r_hincr_test",   Scalar(2.0)     );
00928     magicNumberList.set<Scalar>( "r_hincr",        Scalar(2.0)     );
00929     magicNumberList.set<int>   ( "max_LET_fail",   int(15)         );
00930     magicNumberList.set<Scalar>( "minTimeStep",    Scalar(0.0)     );
00931     magicNumberList.set<Scalar>( "maxTimeStep",    Scalar(10.0)    ); 
00932 
00933     Teuchos::setupVerboseObjectSublist(&*pl);
00934 
00935     validPL = pl;
00936 
00937   }
00938 
00939   return (validPL);
00940   
00941 }
00942 
00943 template<class Scalar>
00944 RCP<Teuchos::ParameterList>
00945 ImplicitBDFStepperStepControl<Scalar>::unsetParameterList()
00946 {
00947   RCP<Teuchos::ParameterList> temp_param_list = parameterList_;
00948   parameterList_ = Teuchos::null;
00949   return(temp_param_list);
00950 }
00951 
00952 template<class Scalar>
00953 RCP<Teuchos::ParameterList> ImplicitBDFStepperStepControl<Scalar>::getParameterList()
00954 {
00955   return(parameterList_);
00956 }
00957 
00958 template<class Scalar>
00959 void ImplicitBDFStepperStepControl<Scalar>::setStepControlData(const StepperBase<Scalar>& stepper)
00960 {
00961   if (stepControlState_ == UNINITIALIZED) {
00962     initialize(stepper);
00963   }
00964   const ImplicitBDFStepper<Scalar>& bdfstepper = Teuchos::dyn_cast<const ImplicitBDFStepper<Scalar> >(stepper);
00965   int desiredOrder = bdfstepper.getOrder();
00966   TEST_FOR_EXCEPT(!((1 <= desiredOrder) && (desiredOrder <= maxOrder_)));
00967   if (stepControlState_ == BEFORE_FIRST_STEP) {
00968     TEST_FOR_EXCEPTION(
00969         desiredOrder > 1, 
00970         std::logic_error, 
00971         "Error, this ImplicitBDF stepper has not taken a step yet, so it cannot take a step of order " << desiredOrder << " > 1!\n"
00972         );
00973   }
00974   TEST_FOR_EXCEPT(!(desiredOrder <= currentOrder_+1));
00975   currentOrder_ = desiredOrder;
00976 
00977   using Teuchos::as;
00978   RCP<Teuchos::FancyOStream> out = this->getOStream();
00979   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00980   Teuchos::OSTab ostab(out,1,"setStepControlData");
00981 
00982   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_EXTREME) ) {
00983     *out << "currentOrder_ = " << currentOrder_ << std::endl;
00984   }
00985 }
00986 
00987 template<class Scalar>
00988 bool ImplicitBDFStepperStepControl<Scalar>::supportsCloning() const
00989 {
00990   return true;
00991 }
00992 
00993 
00994 template<class Scalar>
00995 RCP<StepControlStrategyBase<Scalar> >
00996 ImplicitBDFStepperStepControl<Scalar>::cloneStepControlStrategyAlgorithm() const
00997 {
00998 
00999   RCP<ImplicitBDFStepperStepControl<Scalar> > stepControl = rcp(new ImplicitBDFStepperStepControl<Scalar>());
01000 
01001   if (!is_null(parameterList_)) {
01002     stepControl->setParameterList(parameterList_);
01003   }
01004 
01005   return stepControl;
01006 }
01007 
01008 template<class Scalar>
01009 void ImplicitBDFStepperStepControl<Scalar>::setErrWtVecCalc(const RCP<ErrWtVecCalcBase<Scalar> >& errWtVecCalc)
01010 {
01011   TEST_FOR_EXCEPT(is_null(errWtVecCalc));
01012   errWtVecCalc_ = errWtVecCalc;
01013 }
01014 
01015 template<class Scalar>
01016 RCP<const ErrWtVecCalcBase<Scalar> > ImplicitBDFStepperStepControl<Scalar>::getErrWtVecCalc() const
01017 {
01018   return(errWtVecCalc_);
01019 }
01020 
01021 template<class Scalar>
01022 Scalar ImplicitBDFStepperStepControl<Scalar>::wRMSNorm_(
01023     const Thyra::VectorBase<Scalar>& weight, 
01024     const Thyra::VectorBase<Scalar>& vector) const
01025 {
01026   return(norm_2(weight,vector));
01027 }
01028 
01029 template<class Scalar>
01030 void ImplicitBDFStepperStepControl<Scalar>::defaultInitializeAllData_()
01031 {
01032   typedef Teuchos::ScalarTraits<Scalar> ST;
01033   Scalar zero = ST::zero();
01034   Scalar mone = Scalar(-ST::one());
01035 
01036   stepControlState_ = UNINITIALIZED;
01037   hh_ = zero;
01038   numberOfSteps_ = 0;
01039   stepSizeType_ = STEP_TYPE_VARIABLE;
01040 
01041   maxOrder_ = -1;
01042   nef_ = 0;
01043   midStep_ = false;
01044   checkReduceOrderCalled_ = false;
01045   time_ = -std::numeric_limits<Scalar>::max();
01046   relErrTol_ = mone;
01047   absErrTol_ = mone;
01048   usedStep_ = mone;
01049   currentOrder_ = 1;
01050   usedOrder_ = -1;
01051   nscsco_ = -1;
01052   alpha_s_ = mone;
01053   alpha_0_ = mone;
01054   cj_ = mone;
01055   ck_ = mone;
01056   ck_enorm_ = mone;
01057   constantStepSize_ = false;
01058   Ek_ = mone;
01059   Ekm1_ = mone;
01060   Ekm2_ = mone;
01061   Ekp1_ = mone;
01062   Est_ = mone;
01063   Tk_ = mone;
01064   Tkm1_ = mone;
01065   Tkm2_ = mone;
01066   Tkp1_ = mone;
01067   newOrder_ = -1;
01068   oldOrder_ = -1;
01069   initialPhase_ = false;
01070   stopTime_ = mone;
01071   h0_safety_ = mone;
01072   h0_max_factor_ = mone;
01073   h_phase0_incr_ = mone;
01074   h_max_inv_ = mone;
01075   Tkm1_Tk_safety_ = mone;
01076   Tkp1_Tk_safety_ = mone;
01077   r_factor_ = mone;
01078   r_safety_ = mone;
01079   r_fudge_ = mone;
01080   r_min_ = mone;
01081   r_max_ = mone;
01082   r_hincr_test_ = mone;
01083   r_hincr_ = mone;
01084   max_LET_fail_ = -1;
01085   minTimeStep_ = mone;
01086   maxTimeStep_ = mone;
01087   newtonConvergenceStatus_ = -1;
01088 }
01089 
01090 template<class Scalar>
01091 int ImplicitBDFStepperStepControl<Scalar>::getMaxOrder() const
01092 {
01093   TEST_FOR_EXCEPTION(
01094       stepControlState_ == UNINITIALIZED, std::logic_error,
01095       "Error, attempting to call getMaxOrder before intiialization!\n"
01096       );
01097   return(maxOrder_);
01098 }
01099 
01100 } // namespace Rythmos
01101 #endif // Rythmos_IMPLICITBDF_STEPPER_STEP_CONTROL_H
01102 

Generated on Tue Oct 20 12:46:08 2009 for Rythmos - Transient Integration for Differential Equations by doxygen 1.4.7