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