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