Rythmos_ImplicitBDFStepper.hpp

Go to the documentation of this file.
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_H
00030 #define Rythmos_IMPLICITBDF_STEPPER_H
00031 
00032 #include "Rythmos_Stepper.hpp"
00033 #include "Teuchos_RefCountPtr.hpp"
00034 #include "Thyra_VectorBase.hpp"
00035 #include "Thyra_ModelEvaluator.hpp"
00036 #include "Thyra_ModelEvaluatorHelpers.hpp"
00037 #include "Thyra_NonlinearSolverBase.hpp"
00038 #include "Thyra_SingleResidSSDAEModelEvaluator.hpp"
00039 #include "Thyra_SolveSupportTypes.hpp"
00040 
00041 namespace Rythmos {
00042 
00043 enum BDFactionFlag { ACTION_UNSET, ACTION_LOWER, ACTION_MAINTAIN, ACTION_RAISE };
00044 enum BDFstatusFlag { PREDICT_AGAIN, CONTINUE_ANYWAY, REP_ERR_FAIL, REP_CONV_FAIL };
00045 
00047 template<class Scalar>
00048 class ImplicitBDFStepper : virtual public Stepper<Scalar>
00049 {
00050   public:
00051 
00052     typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType ScalarMag;
00053 
00055     ImplicitBDFStepper(
00056       const Teuchos::RefCountPtr<const Thyra::ModelEvaluator<Scalar> >  &model
00057       ,const Teuchos::RefCountPtr<Thyra::NonlinearSolverBase<Scalar> >  &solver
00058       );
00059 
00061     ImplicitBDFStepper(
00062       const Teuchos::RefCountPtr<const Thyra::ModelEvaluator<Scalar> >  &model
00063       ,const Teuchos::RefCountPtr<Thyra::NonlinearSolverBase<Scalar> >  &solver
00064       ,Teuchos::ParameterList &parameterList
00065       );
00066 
00068     void InitializeStepper(Teuchos::ParameterList &implicitBDFParameters);
00069 
00071     void setModel(const Teuchos::RefCountPtr<const Thyra::ModelEvaluator<Scalar> > &model);
00072 
00074     void setSolver(const Teuchos::RefCountPtr<Thyra::NonlinearSolverBase<Scalar> > &solver);
00075 
00077     Scalar TakeStep(Scalar dt);
00078    
00080     Scalar TakeStep();
00081 
00083     Teuchos::RefCountPtr<const Thyra::VectorBase<Scalar> > get_solution() const;
00084 
00086     Teuchos::RefCountPtr<const Thyra::VectorBase<Scalar> > get_residual() const;
00087 
00089     std::string description() const;
00090 
00092     void describe(
00093       Teuchos::FancyOStream       &out
00094       ,const Teuchos::EVerbosityLevel      verbLevel
00095       ) const;
00096 
00097     bool ErrWtVecSet(
00098       Thyra::VectorBase<Scalar> *w, 
00099       const Thyra::VectorBase<Scalar> &y
00100       );
00101 
00102     Scalar WRMSNorm(
00103       const Thyra::VectorBase<Scalar> &w
00104       , const Thyra::VectorBase<Scalar> &y
00105       ) const;
00106     
00109     bool SetPoints(
00110       const std::vector<Scalar>& time_list
00111       ,const std::vector<Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > >& x_list
00112       ,const std::vector<Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > >& xdot_list);
00113     
00115     bool GetPoints(
00116       const std::vector<Scalar>& time_list
00117       ,std::vector<Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > >* x_list
00118       ,std::vector<Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > >* xdot_list
00119       ,std::vector<ScalarMag>* accuracy_list) const;
00120 
00122     bool SetRange(
00123       const Scalar& time_lower
00124       ,const Scalar& time_upper
00125       ,const InterpolationBuffer<Scalar> & IB);
00126 
00128     bool GetNodes(std::vector<Scalar>* time_list) const;
00129 
00131     bool RemoveNodes(std::vector<Scalar>& time_list) const;
00132 
00134     int GetOrder() const;
00135 
00136 
00137   private:
00138 
00139 
00140     void obtainPredictor();
00141     void obtainResidual();
00142     void obtainJacobian();
00143     // bool interpolateSolution(Scalar time, Thyra::VectorBase<Scalar> &tmpSolVector);
00144     void updateHistory();
00145     void restoreHistory();
00146     void updateCoeffs();
00147     void initialize();
00148     Scalar checkReduceOrder();
00149     BDFstatusFlag rejectStep();
00150     void completeStep();
00151 
00152     void setDefaultMagicNumbers(Teuchos::ParameterList &magicNumberList);
00153 
00154     // 05/05/06 tscoffe:  I hate the underscores for private variables!
00155     Teuchos::RefCountPtr<const Thyra::ModelEvaluator<Scalar> > model;
00156     Teuchos::RefCountPtr<Thyra::NonlinearSolverBase<Scalar> > solver;
00157     Thyra::SingleResidSSDAEModelEvaluator<Scalar>   neModel;
00158 
00159     Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > xn0;
00160     Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > xpn0;
00161     Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > x_dot_base;
00162     std::vector<Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > > xHistory;
00163     Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > ee;
00164     Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > delta;
00165     Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > residual;
00166     Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > errWtVec;
00167 
00168     Scalar time;
00169 
00170 
00171     //typedef typename Thyra::ModelEvaluatorBase::InArgs<Scalar>::ScalarMag ScalarMag;
00172     ScalarMag relErrTol; // relative error tolerance
00173     ScalarMag absErrTol; // absolute error tolerance
00174     Scalar hh;        // Current step-size
00175     int currentOrder; // Current order of integration
00176     int oldOrder;     // previous order of integration
00177     int maxOrder;     // maximum order = min(5,user option maxord) - see below.
00178     int usedOrder;    // order used in current step (used after currentOrder is updated)
00179     Scalar alpha_s;    // $\alpha_s$ fixed-leading coefficient of this BDF method
00180     vector<Scalar> alpha;    // $\alpha_j(n)=h_n/\psi_j(n)$ coefficient used in local error test
00181                       // note:   $h_n$ = current step size, n = current time step
00182     Scalar alpha_0;     // $-\sum_{j=1}^k \alpha_j(n)$ coefficient used in local error test
00183     Scalar cj ;        // $-\alpha_s/h_n$ coefficient used in local error test
00184     Scalar ck ;        // local error coefficient gamma[0] = 0; // $\gamma_j(n)=\sum_{l=1}^{j-1}1/\psi_l(n)$ coefficient used to
00185     vector<Scalar> gamma;    // calculate time derivative of history array for predictor 
00186     vector<Scalar> beta;     // coefficients used to evaluate predictor from history array
00187     vector<Scalar> psi;      // $\psi_j(n) = t_n-t_{n-j}$ intermediary variable used to 
00188                       // compute $\beta_j(n)$
00189     vector<Scalar> sigma;    // $\sigma_j(n) = \frac{h_n^j(j-1)!}{\psi_1(n)*\cdots *\psi_j(n)}$
00190     int numberOfSteps;// number of total time integration steps taken
00191     int  nef;
00192     Scalar usedStep;
00193     int nscsco;
00194     Scalar Ek;
00195     Scalar Ekm1;
00196     Scalar Ekm2;
00197     Scalar Ekp1;
00198     Scalar Est;
00199     Scalar Tk;
00200     Scalar Tkm1;
00201     Scalar Tkm2;
00202     Scalar Tkp1;
00203     int newOrder;
00204     bool initialPhase;
00205     Scalar stopTime;
00206     bool constantStepSize;
00207 
00208     // Magic Numbers:
00209     Scalar h0_safety;
00210     Scalar h0_max_factor;
00211     Scalar h_phase0_incr;
00212     Scalar h_max_inv;
00213     Scalar Tkm1_Tk_safety;
00214     Scalar Tkp1_Tk_safety;
00215     Scalar r_factor;
00216     Scalar r_safety;
00217     Scalar r_fudge;
00218     Scalar r_min;
00219     Scalar r_max;
00220     Scalar r_hincr_test;
00221     Scalar r_hincr;
00222     int    max_LET_fail;
00223     Scalar minTimeStep;
00224     Scalar maxTimeStep;
00225 
00226     int newtonConvergenceStatus;
00227 
00228 #ifdef Rythmos_DEBUG
00229     int debugLevel;
00230     Teuchos::RefCountPtr<Teuchos::FancyOStream> debug_out;
00231 #endif // Rythmos_DEBUG
00232 
00233 
00234 };
00235 
00236 // ////////////////////////////
00237 // Defintions
00238 
00239 template<class Scalar>
00240 void ImplicitBDFStepper<Scalar>::InitializeStepper(
00241     Teuchos::ParameterList &implicitBDFParameters
00242     )
00243 {
00244   typedef Teuchos::ScalarTraits<Scalar> ST;
00245   // Initialize algorithm coefficients
00246   maxOrder = implicitBDFParameters.get("maxOrder",int(5)); // maximum order
00247   maxOrder = max(1,min(maxOrder,5)); // 1 <= maxOrder <= 5
00248   currentOrder=1; // Current order of integration
00249   oldOrder=1; // previous order of integration
00250   usedOrder=1;  // order used in current step (used after currentOrder is updated)
00251   alpha_s=Scalar(-ST::one());  // $\alpha_s$ fixed-leading coefficient of this BDF method
00252   alpha.reserve(maxOrder+1);  // $\alpha_j(n)=h_n/\psi_j(n)$ coefficient used in local error test
00253                   // note:   $h_n$ = current step size, n = current time step
00254   gamma.reserve(maxOrder+1);  // calculate time derivative of history array for predictor 
00255   beta.reserve(maxOrder+1);   // coefficients used to evaluate predictor from history array
00256   psi.reserve(maxOrder+1);    // $\psi_j(n) = t_n-t_{n-j}$ intermediary variable used to 
00257                   // compute $\beta_j(n;$
00258   sigma.reserve(maxOrder+1);  // $\sigma_j(n) = \frac{h_n^j(j-1)!}{\psi_1(n)*\cdots *\psi_j(n)}$
00259   for (int i=0 ; i<maxOrder ; ++i)
00260   {
00261     alpha.push_back(ST::zero());
00262     beta.push_back(ST::zero());
00263     gamma.push_back(ST::zero());
00264     psi.push_back(ST::zero());
00265     sigma.push_back(ST::zero());
00266   }
00267   alpha_0=ST::zero();   // $-\sum_{j=1}^k \alpha_j(n)$ coefficient used in local error test
00268   cj=ST::zero();      // $-\alpha_s/h_n$ coefficient used in local error test
00269   ck=ST::zero();      // local error coefficient gamma_[0] = 0; // $\gamma_j(n)=\sum_{l=1}^{j-1}1/\psi_l(n)$ coefficient used to
00270   hh=ST::zero();
00271   numberOfSteps=0;   // number of total time integration steps taken
00272   nef=0;
00273   usedStep=ST::zero();
00274   nscsco=0;
00275   Ek=ST::zero();
00276   Ekm1=ST::zero();
00277   Ekm2=ST::zero();
00278   Ekp1=ST::zero();
00279   Est=ST::zero();
00280   Tk=ST::zero();
00281   Tkm1=ST::zero();
00282   Tkm2=ST::zero();
00283   Tkp1=ST::zero();
00284   newOrder=1;
00285   initialPhase=true;
00286 
00287   relErrTol = implicitBDFParameters.get( "relErrTol", Scalar(1.0e-4) );
00288   absErrTol = implicitBDFParameters.get( "absErrTol", Scalar(1.0e-6) );
00289   constantStepSize = implicitBDFParameters.get( "constantStepSize", false );
00290   stopTime = implicitBDFParameters.get( "stopTime", Scalar(10.0) );
00291 
00292 
00293 #ifdef Rythmos_DEBUG
00294   debugLevel = implicitBDFParameters.get( "debugLevel", int(1) );
00295   debug_out = Teuchos::VerboseObjectBase::getDefaultOStream();
00296   debug_out->precision(15);
00297   debug_out->setMaxLenLinePrefix(28);
00298   debug_out->pushLinePrefix("Rythmos::ImplicitBDFStepper");
00299   debug_out->setShowLinePrefix(true);
00300   debug_out->setTabIndentStr("    ");
00301 #endif // Rythmos_DEBUG
00302 
00303   setDefaultMagicNumbers(implicitBDFParameters.sublist("magicNumbers"));
00304 
00305 #ifdef Rythmos_DEBUG
00306   Teuchos::OSTab ostab(debug_out,1,"InitializeStepper");
00307   if (debugLevel > 1)
00308   {
00309     *debug_out << "maxOrder = " << maxOrder << endl;
00310     *debug_out << "currentOrder = " << currentOrder << endl;
00311     *debug_out << "oldOrder = " << oldOrder << endl;
00312     *debug_out << "usedOrder = " << usedOrder << endl;
00313     *debug_out << "alpha_s = " << alpha_s << endl;
00314     for (int i=0 ; i<maxOrder ; ++i)
00315     {
00316       *debug_out << "alpha[" << i << "] = " << alpha[i] << endl;
00317       *debug_out << "beta[" << i << "] = " << beta[i] << endl;
00318       *debug_out << "gamma[" << i << "] = " << gamma[i] << endl;
00319       *debug_out << "psi[" << i << "] = " << psi[i] << endl;
00320       *debug_out << "sigma[" << i << "] = " << sigma[i] << endl;
00321     }
00322     *debug_out << "alpha_0 = " << alpha_0 << endl;
00323     *debug_out << "cj = " << cj << endl;
00324     *debug_out << "ck = " << ck << endl;
00325     *debug_out << "numberOfSteps = " << numberOfSteps << endl;
00326     *debug_out << "nef = " << nef << endl;
00327     *debug_out << "usedStep = " << usedStep << endl;
00328     *debug_out << "nscsco = " << nscsco << endl;
00329     *debug_out << "Ek = " << Ek << endl;
00330     *debug_out << "Ekm1 = " << Ekm1 << endl;
00331     *debug_out << "Ekm2 = " << Ekm2 << endl;
00332     *debug_out << "Ekp1 = " << Ekp1 << endl;
00333     *debug_out << "Est = " << Est << endl;
00334     *debug_out << "Tk = " << Tk << endl;
00335     *debug_out << "Tkm1 = " << Tkm1 << endl;
00336     *debug_out << "Tkm2 = " << Tkm2 << endl;
00337     *debug_out << "Tkp1 = " << Tkp1 << endl;
00338     *debug_out << "newOrder = " << newOrder << endl;
00339     *debug_out << "initialPhase = " << initialPhase << endl;
00340     *debug_out << "relErrTol  = " << relErrTol  << endl;
00341     *debug_out << "absErrTol  = " << absErrTol  << endl;
00342     *debug_out << "constantStepSize  = " << constantStepSize  << endl;
00343     *debug_out << "stopTime  = " << stopTime  << endl;
00344   }
00345 #endif // Rythmos_DEBUG
00346 }
00347 
00348 template<class Scalar>
00349 ImplicitBDFStepper<Scalar>::ImplicitBDFStepper(
00350   const Teuchos::RefCountPtr<const Thyra::ModelEvaluator<Scalar> > &model
00351   ,const Teuchos::RefCountPtr<Thyra::NonlinearSolverBase<Scalar> > &solver
00352   ,Teuchos::ParameterList &parameterList
00353   )
00354 {
00355   InitializeStepper(parameterList);
00356 
00357   // Now we instantiate the model and the solver
00358   setModel(model);
00359   setSolver(solver);
00360 }
00361 
00362 template<class Scalar>
00363 ImplicitBDFStepper<Scalar>::ImplicitBDFStepper(
00364   const Teuchos::RefCountPtr<const Thyra::ModelEvaluator<Scalar> > &model
00365   ,const Teuchos::RefCountPtr<Thyra::NonlinearSolverBase<Scalar> > &solver
00366   )
00367 {
00368   Teuchos::ParameterList emptyList;
00369   InitializeStepper(emptyList);
00370 
00371   // Now we instantiate the model and the solver
00372   setModel(model);
00373   setSolver(solver);
00374 }
00375 
00376 template<class Scalar>
00377 void ImplicitBDFStepper<Scalar>::setDefaultMagicNumbers(
00378     Teuchos::ParameterList &magicNumberList)
00379 {
00380   // Magic Number Defaults:
00381   h0_safety      = magicNumberList.get( "h0_safety",      Scalar(2.0)     );
00382   h0_max_factor  = magicNumberList.get( "h0_max_factor",  Scalar(0.001)   );
00383   h_phase0_incr  = magicNumberList.get( "h_phase0_incr",  Scalar(2.0)     );
00384   h_max_inv      = magicNumberList.get( "h_max_inv",      Scalar(0.0)     );
00385   Tkm1_Tk_safety = magicNumberList.get( "Tkm1_Tk_safety", Scalar(2.0)     );
00386   Tkp1_Tk_safety = magicNumberList.get( "Tkp1_Tk_safety", Scalar(0.5)     );
00387   r_factor       = magicNumberList.get( "r_factor",       Scalar(0.9)     );
00388   r_safety       = magicNumberList.get( "r_safety",       Scalar(2.0)     );
00389   r_fudge        = magicNumberList.get( "r_fudge",        Scalar(0.0001)  );
00390   r_min          = magicNumberList.get( "r_min",          Scalar(0.125)   );
00391   r_max          = magicNumberList.get( "r_max",          Scalar(0.9)     );
00392   r_hincr_test   = magicNumberList.get( "r_hincr_test",   Scalar(2.0)     );
00393   r_hincr        = magicNumberList.get( "r_hincr",        Scalar(2.0)     );
00394   max_LET_fail   = magicNumberList.get( "max_LET_fail",   int(15)         );
00395   minTimeStep    = magicNumberList.get( "minTimeStep",    Scalar(0.0)     );
00396   maxTimeStep    = magicNumberList.get( "maxTimeStep",    Scalar(10.0)    ); 
00397 
00398 #ifdef Rythmos_DEBUG
00399   Teuchos::OSTab ostab(debug_out,1,"setDefaultMagicNumbers");
00400   if (debugLevel > 1)
00401   {
00402     *debug_out << "h0_safety = " << h0_safety << endl;
00403     *debug_out << "h0_max_factor = " << h0_max_factor << endl;
00404     *debug_out << "h_phase0_incr = " << h_phase0_incr << endl;
00405     *debug_out << "h_max_inv = " << h_max_inv << endl;
00406     *debug_out << "Tkm1_Tk_safety = " << Tkm1_Tk_safety  << endl;
00407     *debug_out << "Tkp1_Tk_safety = " << Tkp1_Tk_safety << endl;
00408     *debug_out << "r_factor = " << r_factor << endl;
00409     *debug_out << "r_safety = " << r_safety << endl;
00410     *debug_out << "r_fudge = " << r_fudge << endl;
00411     *debug_out << "r_min = " << r_min << endl;
00412     *debug_out << "r_max = " << r_max << endl;
00413     *debug_out << "r_hincr_test = " << r_hincr_test << endl;
00414     *debug_out << "r_hincr = " << r_hincr << endl;
00415     *debug_out << "max_LET_fail = " << max_LET_fail << endl;
00416     *debug_out << "minTimeStep = " << minTimeStep << endl;
00417     *debug_out << "maxTimeStep = " << maxTimeStep << endl;
00418   }
00419 #endif // Rythmos_DEBUG
00420 }
00421 
00422 template<class Scalar>
00423 void ImplicitBDFStepper<Scalar>::setModel(const Teuchos::RefCountPtr<const Thyra::ModelEvaluator<Scalar> > &model_)
00424 {
00425   typedef Teuchos::ScalarTraits<Scalar> ST;
00426   model = model_;
00427   time = ST::zero();
00428   xn0 = model->getNominalValues().get_x()->clone_v();
00429   xpn0 = model->getNominalValues().get_x_dot()->clone_v(); 
00430   x_dot_base = model->getNominalValues().get_x_dot()->clone_v();
00431   ee = model->getNominalValues().get_x()->clone_v();
00432   delta = model->getNominalValues().get_x()->clone_v();
00433   residual = Thyra::createMember(model->get_f_space());
00434   errWtVec = xn0->clone_v();
00435   ErrWtVecSet(&*errWtVec,*xn0);
00436   xHistory.push_back(xn0->clone_v());
00437   xHistory.push_back(xpn0->clone_v());
00438   for (int i=2 ; i<=maxOrder ; ++i) // Store maxOrder+1 vectors
00439   {
00440     xHistory.push_back(xn0->clone_v()); 
00441     V_S(&*xHistory[i],ST::zero());
00442   }
00443   initialize(); // Now that we have the model, we can do our initialization 
00444 }
00445 
00446 template<class Scalar>
00447 void ImplicitBDFStepper<Scalar>::setSolver(const Teuchos::RefCountPtr<Thyra::NonlinearSolverBase<Scalar> > &solver_)
00448 {
00449   solver = solver_;
00450 }
00451 
00452 template<class Scalar>
00453 Scalar ImplicitBDFStepper<Scalar>::TakeStep()
00454 {
00455   typedef Teuchos::ScalarTraits<Scalar> ST;
00456   typedef typename Thyra::ModelEvaluatorBase::InArgs<Scalar>::ScalarMag ScalarMag;
00457 #ifdef Rythmos_DEBUG
00458   Teuchos::OSTab ostab(debug_out,1,"TakeStep");
00459 #endif // Rythmos_DEBUG
00460   BDFstatusFlag status;
00461   while (1)
00462   {
00463     // Set up problem coefficients (and handle first step)
00464     updateCoeffs();
00465     // Set Error weight vector
00466     ErrWtVecSet(&*errWtVec,*xn0);
00467     // compute predictor
00468     obtainPredictor();
00469     // solve nonlinear problem (as follows)
00470     
00471     //
00472     // Setup the nonlinear equations:
00473     //
00474     //   f_bar( x_dot_coeff * x_bar + x_dot_base, x_coeff * x_bar + x_base, t_base ) = 0
00475     //   x_dot_coeff = -alpha_s/dt
00476     //   x_dot_base = x_prime_pred + (alpha_s/dt) * x_pred
00477     //   x_coeff = 1
00478     //   x_base = 0
00479     //   t_base = tn+dt
00480     //
00481     Scalar coeff_x_dot = Scalar(-ST::one())*alpha_s/hh;
00482     V_StVpStV( &*x_dot_base, ST::one(), *xpn0, alpha_s/hh, *xn0 );
00483     neModel.initialize(model,coeff_x_dot,x_dot_base,ST::one(),Teuchos::null,time+hh,xn0);
00484     //
00485     // Solve the implicit nonlinear system to a tolerance of ???
00486     // 
00487     // 05/08/06 tscoffe:  I really need to get the update, not the solution from
00488     // the nonlinear solver.
00489     if(solver->getModel().get()!=&neModel)
00490       solver->setModel( Teuchos::rcp(&neModel,false) );
00491     /* // Thyra::TimeStepNewtonNonlinearSolver uses a built in solveCriteria, so you can't pass one in.
00492        // I believe this is the correct solveCriteria for IDA though.
00493     Thyra::SolveMeasureType nonlinear_solve_measure_type(Thyra::SOLVE_MEASURE_NORM_RESIDUAL,Thyra::SOLVE_MEASURE_ONE); 
00494     ScalarMag tolerance = relErrTol/ScalarMag(10.0); // This should be changed to match the condition in IDA
00495     Thyra::SolveCriteria<Scalar> nonlinearSolveCriteria(nonlinear_solve_measure_type, tolerance);
00496     Thyra::SolveStatus<Scalar> nonlinearSolveStatus = solver->solve( &*xn0, &nonlinearSolveCriteria, &*delta ); 
00497     */
00498     //Thyra::assign(&*xn0,ST::zero()); // 08/10/06 tscoffe:  what is this doing here?  It hoses the solve.
00499     Thyra::SolveStatus<Scalar> nonlinearSolveStatus = solver->solve( &*xn0, NULL, &*ee ); 
00500     if (nonlinearSolveStatus.solveStatus == Thyra::SOLVE_STATUS_CONVERGED) 
00501       newtonConvergenceStatus = 0;
00502     else 
00503       newtonConvergenceStatus = -1;
00504 
00505     // check error and evaluate LTE
00506     Scalar enorm = checkReduceOrder();
00507     
00508 #ifdef Rythmos_DEBUG
00509     if (debugLevel > 1)
00510     {
00511       *debug_out << "xn0 = " << std::endl;
00512       xn0->describe(*debug_out,Teuchos::VERB_EXTREME);
00513       *debug_out << "ee = " << std::endl;
00514       ee->describe(*debug_out,Teuchos::VERB_EXTREME);
00515       *debug_out << "delta = " << std::endl;
00516       delta->describe(*debug_out,Teuchos::VERB_EXTREME);
00517       for (int i=0; i<max(2,maxOrder); ++i)
00518       {
00519         *debug_out << "xHistory[" << i << "] = " << std::endl;
00520         xHistory[i]->describe(*debug_out,Teuchos::VERB_EXTREME);
00521       }
00522       *debug_out << "ck = " << ck << endl;
00523       *debug_out << "enorm = " << enorm << endl;
00524       *debug_out << "Local Truncation Error Check: (ck*enorm) < 1:  (" << ck*enorm << ") <?= 1" << endl;
00525     }
00526 #endif // Rythmos_DEBUG
00527     // Check LTE here:
00528     if ((ck*enorm) > ST::one())
00529       status = rejectStep(); 
00530     else 
00531       break;
00532     if (status == CONTINUE_ANYWAY)
00533       break;
00534     if (status == REP_ERR_FAIL)
00535       return(Scalar(-ST::one()));
00536   }
00537 
00538   completeStep();  
00539   return(usedStep);
00540 }
00541 
00542 template<class Scalar>
00543 Scalar ImplicitBDFStepper<Scalar>::TakeStep(Scalar dt)
00544 {
00545   constantStepSize = true;
00546   hh = dt;
00547   dt = TakeStep();
00548   return(dt);
00549 }
00550 
00551 template<class Scalar>
00552 Teuchos::RefCountPtr<const Thyra::VectorBase<Scalar> > ImplicitBDFStepper<Scalar>::get_solution() const
00553 {
00554   return(xHistory[0]);
00555 }
00556 
00557 template<class Scalar>
00558 Teuchos::RefCountPtr<const Thyra::VectorBase<Scalar> > ImplicitBDFStepper<Scalar>::get_residual() const
00559 {
00560   return(residual);
00561 }
00562 
00563 template<class Scalar>
00564 std::string ImplicitBDFStepper<Scalar>::description() const
00565 {
00566   std::string name = "Rythmos::ImplicitBDFStepper";
00567   return(name);
00568 }
00569 
00570 template<class Scalar>
00571 void ImplicitBDFStepper<Scalar>::describe(
00572       Teuchos::FancyOStream                &out
00573       ,const Teuchos::EVerbosityLevel      verbLevel
00574       ) const
00575 {
00576   if (verbLevel == Teuchos::VERB_EXTREME)
00577   {
00578     out << description() << "::describe" << std::endl;
00579     out << "model = " << std::endl;
00580     model->describe(out,verbLevel);
00581     out << "solver = " << std::endl;
00582     solver->describe(out,verbLevel);
00583     out << "xn0 = " << std::endl;
00584     xn0->describe(out,verbLevel);
00585     out << "xpn0 = " << std::endl;
00586     xpn0->describe(out,verbLevel);
00587     out << "x_dot_base = " << std::endl;
00588     x_dot_base->describe(out,verbLevel);
00589     out << "xHistory = " << std::endl;
00590     for (int i=0 ; i < max(2,maxOrder) ; ++i)
00591     {
00592       out << "xHistory[" << i << "] = " << std::endl;
00593       xHistory[i]->describe(out,verbLevel);
00594     }
00595     out << "ee = " << std::endl;
00596     ee->describe(out,verbLevel);
00597     out << "delta = " << std::endl;
00598     delta->describe(out,verbLevel);
00599     out << "residual = " << std::endl;
00600     residual->describe(out,verbLevel);
00601     out << "errWtVec = " << std::endl;
00602     errWtVec->describe(out,verbLevel);
00603     out << "time = " << time << std::endl;
00604     out << "hh = " << hh << std::endl;
00605     out << "currentOrder = " << currentOrder << std::endl;
00606     out << "neModel = " << std::endl;
00607     neModel.describe(out,verbLevel);
00608   }
00609 }
00610 
00611 template<class Scalar>
00612 void ImplicitBDFStepper<Scalar>::obtainPredictor()
00613 {
00614   typedef Teuchos::ScalarTraits<Scalar> ST;
00615 
00616 #ifdef Rythmos_DEBUG
00617   Teuchos::OSTab ostab(debug_out,1,"obtainPredictor");
00618   if (debugLevel > 1)
00619   {
00620     *debug_out << "currentOrder = " << currentOrder << std::endl;
00621   }
00622 #endif // Rythmos_DEBUG
00623   
00624   // prepare history array for prediction
00625   for (int i=nscsco;i<=currentOrder;++i)
00626   {
00627     Vt_S(&*xHistory[i],beta[i]);
00628   }
00629   
00630   // evaluate predictor
00631   V_V(&*xn0,*xHistory[0]);
00632   V_S(&*xpn0,ST::zero());
00633   for (int i=1;i<=currentOrder;++i)
00634   {
00635     Vp_V(&*xn0,*xHistory[i]);
00636     Vp_StV(&*xpn0,gamma[i],*xHistory[i]);
00637   }
00638 #ifdef Rythmos_DEBUG
00639   if (debugLevel > 1)
00640   {
00641     *debug_out << "xn0 = " << std::endl;
00642     xn0->describe(*debug_out,Teuchos::VERB_EXTREME);
00643     *debug_out << "xpn0 = " << std::endl;
00644     xpn0->describe(*debug_out,Teuchos::VERB_EXTREME);
00645   }
00646 #endif // Rythmos_DEBUG
00647 }
00648 
00649 /* // This function is not ready yet.
00650 template<class Scalar>
00651 void ImplicitBDFStepper<Scalar>::interpolateSolution(Scalar timepoint, 
00652       Thyra::VectorBase<Scalar>  &tmpSolVector)
00653 {
00654   typedef Teuchos::ScalarTraits<Scalar> ST;
00655 
00656 // 03/23/04 tscoffe:  Currently this code is nearly identical to the IDA code
00657 // for interpolating to an output time.  Either we acknowledge the copyright,
00658 // the list of conditions in the license and the disclaimer or we rewrite this
00659 // function.  The IDA license is included after this routine.
00660   Scalar tfuzz;   // fuzz factor to check for valid output time
00661   Scalar tp;      // approximately t{n-1}
00662   Scalar delt;    // distance between timepoint and time
00663   Scalar c = Scalar(1.0); // coefficient for interpolation
00664   Scalar gam;     // coefficient for interpolation
00665   int kord;       // order of interpolation
00666   Scalar tn = time;
00667   Scalar hh = hh;
00668   Scalar hused = usedStep;
00669   int kused = usedOrder;
00670   Scalar uround = ST::zero();  // unit round-off (set to zero for now)
00671 
00672   tfuzz = 100 * uround * (tn + hh);
00673   tp = tn - hused - tfuzz;
00674   if ( (timepoint - tp)*hh < ST::zero() ) 
00675     return false;
00676 
00677   assign(&*tmpSolVector,*xHistory[0]),
00678   kord = kused;
00679   if ( (kused == 0) || (timepoint == tn) ) 
00680     kord = 1;
00681 
00682   delt = timepoint - tn;
00683   gam = delt/psi[0];
00684   for (int j=1 ; j <= kord ; ++j)
00685   {
00686     c = c*gam;
00687     gam = (delt + psi[j-1])/psi[j];
00688     Vp_StV(&*tmpSolVector,c,*xHistory[j]);
00689   }
00690   return true;
00691 }
00692 */
00693 
00694 template<class Scalar>
00695 void ImplicitBDFStepper<Scalar>::updateHistory()
00696 {
00697 
00698   // Save Newton correction for potential order increase on next step.
00699   if (usedOrder < maxOrder)  
00700   {
00701     assign( &*xHistory[usedOrder+1], *ee );
00702   }
00703   // Update history arrays
00704   Vp_V( &*xHistory[usedOrder], *ee );
00705   for (int j=usedOrder-1;j>=0;j--) 
00706   {
00707     Vp_V( &*xHistory[j], *xHistory[j+1] );
00708   }
00709 #ifdef Rythmos_DEBUG
00710   Teuchos::OSTab ostab(debug_out,1,"updateHistory");
00711   if (debugLevel > 1)
00712   {
00713     for (int i=0;i<max(2,maxOrder);++i)
00714     {
00715       *debug_out << "xHistory[" << i << "] = " << endl;
00716       xHistory[i]->describe(*debug_out,Teuchos::VERB_EXTREME);
00717     }
00718   }
00719 #endif // Rythmos_DEBUG
00720 
00721 }
00722 
00723 template<class Scalar>
00724 void ImplicitBDFStepper<Scalar>::restoreHistory()
00725 {
00726   typedef Teuchos::ScalarTraits<Scalar> ST;
00727 
00728   // undo preparation of history array for prediction
00729   for (int i=nscsco;i<=currentOrder;++i)
00730   {
00731     Vt_S( &*xHistory[i], ST::one()/beta[i] );
00732   }
00733   for (int i=1;i<=currentOrder;++i)
00734   {
00735     psi[i-1] = psi[i] - hh;
00736   }
00737 #ifdef Rythmos_DEBUG
00738   Teuchos::OSTab ostab(debug_out,1,"restoreHistory");
00739   if (debugLevel > 1)
00740   {
00741     for (int i=0;i<maxOrder;++i)
00742       *debug_out << "psi[" << i << "] = " << psi[i] << endl;
00743     for (int i=0;i<maxOrder;++i)
00744     {
00745       *debug_out << "xHistory[" << i << "] = " << endl;
00746       xHistory[i]->describe(*debug_out,Teuchos::VERB_EXTREME);
00747     }
00748   }
00749 #endif // Rythmos_DEBUG
00750 } 
00751 
00752 template<class Scalar>
00753 void ImplicitBDFStepper<Scalar>::updateCoeffs()
00754 {
00755   typedef Teuchos::ScalarTraits<Scalar> ST;
00756   // If the number of steps taken with constant order and constant stepsize is
00757   // more than the current order + 1 then we don't bother to update the
00758   // coefficients because we've reached a constant step-size formula.  When
00759   // this is is not true, then we update the coefficients for the variable
00760   // step-sizes. 
00761   if ((hh != usedStep) || (currentOrder != usedOrder))
00762     nscsco = 0;
00763   nscsco = min(nscsco+1,usedOrder+2);
00764   if (currentOrder+1 >= nscsco)
00765   {
00766     beta[0] = ST::one();
00767     alpha[0] = ST::one();
00768     Scalar temp1 = hh;
00769     sigma[0] = ST::one();
00770     gamma[0] = ST::zero();
00771     for (int i=1;i<=currentOrder;++i)
00772     {
00773       Scalar temp2 = psi[i-1];
00774       psi[i-1] = temp1;
00775       beta[i] = beta[i-1]*psi[i-1]/temp2;
00776       temp1 = temp2 + hh;
00777       alpha[i] = hh/temp1;
00778       sigma[i] = Scalar(i+1)*sigma[i-1]*alpha[i];
00779       gamma[i] = gamma[i-1]+alpha[i-1]/hh;
00780     }
00781     psi[currentOrder] = temp1;
00782     alpha_s = ST::zero();
00783     alpha_0 = ST::zero();
00784     for (int i=0;i<currentOrder;++i)
00785     {
00786       alpha_s = alpha_s - Scalar(ST::one()/(i+ST::one()));
00787       alpha_0 = alpha_0 - alpha[i];
00788     }
00789     cj = -alpha_s/hh;
00790     ck = abs(alpha[currentOrder]+alpha_s-alpha_0);
00791     ck = max(ck,alpha[currentOrder]);
00792   }
00793 #ifdef Rythmos_DEBUG
00794   Teuchos::OSTab ostab(debug_out,1,"updateCoeffs");
00795   if (debugLevel > 1)
00796   {
00797     for (int i=0;i<=maxOrder;++i)
00798     {
00799       *debug_out << "alpha[" << i << "] = " << alpha[i] << endl;
00800       *debug_out << "beta[" << i << "] = " << beta[i] << endl;
00801       *debug_out << "sigma[" << i << "] = " << sigma[i] << endl;
00802       *debug_out << "gamma[" << i << "] = " << gamma[i] << endl;
00803       *debug_out << "psi[" << i << "] = " << psi[i] << endl;
00804       *debug_out << "alpha_s = " << alpha_s << endl;
00805       *debug_out << "alpha_0 = " << alpha_0 << endl;
00806       *debug_out << "cj = " << cj << endl;
00807       *debug_out << "ck = " << ck << endl;
00808     }
00809   }
00810 #endif // Rythmos_DEBUG
00811 
00812 }
00813 
00814 template<class Scalar>
00815 void ImplicitBDFStepper<Scalar>::initialize()
00816 {
00817   typedef Teuchos::ScalarTraits<Scalar> ST;
00818 
00819   // Choose initial step-size
00820   Scalar time_to_stop = stopTime - time;
00821   Scalar currentTimeStep;
00822   if (constantStepSize)
00823   {
00824     currentTimeStep = hh;
00825     //currentTimeStep = 0.1 * time_to_stop;
00826     //currentTimeStep = min(hh, currentTimeStep);
00827   }
00828   else
00829   {
00830     // compute an initial step-size based on rate of change in the solution initially
00831     Scalar ypnorm = WRMSNorm(*errWtVec,*xHistory[1]);
00832     if (ypnorm > ST::zero())  // time-dependent DAE
00833     {
00834       currentTimeStep = min(h0_max_factor*abs(time_to_stop),sqrt(2.0)/(h0_safety*ypnorm));
00835     } 
00836     else  // non-time-dependent DAE
00837     {
00838       currentTimeStep = h0_max_factor*abs(time_to_stop);
00839     }
00840     // choose min of user specified value and our value:
00841     if (hh > ST::zero())
00842       currentTimeStep = min(hh, currentTimeStep);
00843     // check for maximum step-size:
00844     Scalar rh = abs(currentTimeStep)*h_max_inv; 
00845     if (rh>1.0) currentTimeStep = currentTimeStep/rh;
00846   }
00847   hh = currentTimeStep;
00848 #ifdef Rythmos_DEBUG
00849   Teuchos::OSTab ostab(debug_out,1,"initialize");
00850   if (debugLevel > 1)
00851   {
00852     *debug_out << "hh = " << hh << endl;
00853   }
00854 #endif // Rythmos_DEBUG
00855 
00856   // x history
00857   assign(&*xHistory[0],*xn0);
00858   V_S(&*xHistory[1],ST::zero());
00859 
00860   // Coefficient initialization 
00861   numberOfSteps = 0;    // number of total time integration steps taken
00862   currentOrder = 1;
00863   usedOrder = 1;
00864   psi[0] = hh;
00865   cj = 1/psi[0];
00866   nscsco = 0;
00867 }
00868 
00869 template<class Scalar>
00870 Scalar ImplicitBDFStepper<Scalar>::checkReduceOrder()
00871 {
00872 // This routine puts its output in newOrder_
00873 
00874 // This routine changes the following variables:
00875 //    Ek, Tk, Est, newOrder, dsDae.delta_x, dsDae.delta_q,
00876 //    Ekm1, Tkm1, Ekm2, Tkm2 
00877 
00878 // This routine reads but does not change the following variables:
00879 //    currentOrder, sigma, dsDae.newtonCorrectionPtr, dsDae.qNewtonCorrectionPtr,
00880 //    dsDae.errWtVecPtr, dsDae.qErrWtVecPtr, dsDae.xHistory, dsDae.qHistory
00881 
00882   Scalar enorm = WRMSNorm(*errWtVec,*ee);
00883   Ek = sigma[currentOrder]*enorm;
00884   Tk = Scalar(currentOrder+1)*Ek;
00885   Est = Ek;
00886   newOrder = currentOrder;
00887 #ifdef Rythmos_DEBUG
00888   Teuchos::OSTab ostab(debug_out,1,"checkReduceOrder");
00889   if (debugLevel > 1)
00890   {
00891     *debug_out << "currentOrder = " << currentOrder << std::endl;
00892     *debug_out << "Ek = " << Ek << std::endl;
00893     *debug_out << "Tk = " << Tk << std::endl;
00894     *debug_out << "enorm = " << enorm << std::endl;
00895   }
00896 #endif // Rythmos_DEBUG
00897   if (currentOrder>1)
00898   {
00899     V_VpV(&*delta,*xHistory[currentOrder],*ee);
00900     Ekm1 = sigma[currentOrder-1]*WRMSNorm(*errWtVec,*delta);
00901     Tkm1 = currentOrder*Ekm1;
00902 #ifdef Rythmos_DEBUG
00903     if (debugLevel > 1)
00904     {
00905       *debug_out << "Ekm1 = " << Ekm1 << endl;
00906       *debug_out << "Tkm1 = " << Tkm1 << endl;
00907     }
00908 #endif // Rythmos_DEBUG
00909     if (currentOrder>2)
00910     {
00911       Vp_V(&*delta,*xHistory[currentOrder-1]);
00912       Ekm2 = sigma[currentOrder-2]*WRMSNorm(*errWtVec,*delta);
00913       Tkm2 = (currentOrder-1)*Ekm2;
00914 #ifdef Rythmos_DEBUG
00915       if (debugLevel > 1)
00916       {
00917         *debug_out << "Ekm2 = " << Ekm2 << endl;
00918         *debug_out << "Tkm2 = " << Tkm2 << endl;
00919       }
00920 #endif // Rythmos_DEBUG
00921       if (max(Tkm1,Tkm2)<=Tk)
00922       {
00923         newOrder--;
00924         Est = Ekm1;
00925       }
00926     }
00927     else if (Tkm1 <= Tkm1_Tk_safety * Tk)
00928     {
00929       newOrder--;
00930       Est = Ekm1;
00931     }
00932   }
00933 #ifdef Rythmos_DEBUG
00934   if (debugLevel > 1)
00935   {
00936     *debug_out << "Est = " << Est << endl;
00937     *debug_out << "newOrder= " << newOrder << endl;
00938   }
00939 #endif // Rythmos_DEBUG
00940   return(enorm);
00941 }
00942 
00943 template<class Scalar>
00944 BDFstatusFlag ImplicitBDFStepper<Scalar>::rejectStep()
00945 {
00946 
00947 // This routine puts its output in newTimeStep and newOrder
00948 
00949 // This routine changes the following variables:
00950 //    initialPhase, nef, psi, newTimeStep,
00951 //    newOrder, currentOrder, currentTimeStep, dsDae.xHistory,
00952 //    dsDae.qHistory, nextTimePt, 
00953 //    currentTimeStepSum, nextTimePt
00954 
00955 // This routine reads but does not change the following variables:
00956 //    r_factor, r_safety, Est, r_fudge, r_min, r_max,
00957 //    minTimeStep, maxTimeStep, time, stopTime 
00958 
00959   // Only update the time step if we are NOT running constant stepsize.
00960   bool adjustStep = (!constantStepSize);
00961 
00962   Scalar newTimeStep = hh;
00963   Scalar rr = 1.0; // step size ratio = new step / old step
00964   nef++;
00965 #ifdef Rythmos_DEBUG
00966   Teuchos::OSTab ostab(debug_out,1,"rejectStep");
00967   if (debugLevel > 1)
00968   {
00969     *debug_out << "adjustStep = " << adjustStep << endl;
00970     *debug_out << "nef = " << nef << endl;
00971   }
00972 #endif // Rythmos_DEBUG
00973   if (nef >= max_LET_fail)  
00974   {
00975     cerr << "Rythmos_Stepper_ImplicitBDF::rejectStep:  " 
00976           << "  Maximum number of local error test failures.  " << endl;
00977     return(REP_ERR_FAIL);
00978   }
00979   if (adjustStep)
00980   {
00981     initialPhase = false;
00982 #ifdef Rythmos_DEBUG
00983     if (debugLevel > 1)
00984     {
00985       *debug_out << "initialPhase = " << initialPhase << endl;
00986     }
00987 #endif // Rythmos_DEBUG
00988     restoreHistory();
00989     // restore psi_
00990 //    for (int i=1;i<=currentOrder;++i)
00991 //      psi[i-1] = psi[i] - hh;
00992 
00993     if ((newtonConvergenceStatus < 0))
00994     {
00996       // rely on the error estimate - it may be full of Nan's.
00997       rr = r_min;
00998       newTimeStep = rr * hh;
00999 
01000       if (nef > 2) newOrder = 1;//consistent with block below.
01001     }
01002     else
01003     {
01004       // 03/11/04 tscoffe:  Here is the block for choosing order & 
01005       // step-size when the local error test FAILS (but Newton 
01006       // succeeded). 
01007       if (nef == 1) // first local error test failure
01008       {
01009         rr = r_factor*pow(r_safety*Est+r_fudge,-1.0/(newOrder+1.0));
01010         rr = max(r_min,min(r_max,rr));
01011         newTimeStep = rr * hh;
01012       }
01013       else if (nef == 2) // second failure
01014       {
01015         rr = r_min;
01016         newTimeStep = rr * hh;
01017       }
01018       else if (nef > 2) // third and later failures
01019       {
01020         newOrder = 1;
01021         rr = r_min;
01022         newTimeStep = rr * hh;
01023       }
01024     }
01025 #ifdef Rythmos_DEBUG
01026     if (debugLevel > 1)
01027     {
01028       *debug_out << "rr = " << rr << endl;
01029       *debug_out << "newOrder = " << newOrder << endl;
01030     }
01031 #endif // Rythmos_DEBUG
01032     currentOrder = newOrder;
01033     if (numberOfSteps == 0) // still first step
01034     {
01035       psi[0] = newTimeStep;
01036       Vt_S(&*xHistory[1],rr);
01037 #ifdef Rythmos_DEBUG
01038       if (debugLevel > 1)
01039       {
01040         *debug_out << "numberOfSteps == 0:" << endl;
01041         *debug_out << "psi[0] = " << psi[0] << endl;
01042         *debug_out << "xHistory[1] = " << std::endl;
01043         xHistory[1]->describe(*debug_out,Teuchos::VERB_EXTREME);
01044       }
01045 #endif // Rythmos_DEBUG
01046     }
01047   }
01048   else if (!adjustStep)
01049   {
01050     cerr << "Rythmos_Stepper_ImplicitBDF::rejectStep:  "
01051          << "Warning: Local error test failed with constant step-size." << endl;
01052   }
01053 
01054   BDFstatusFlag return_status = PREDICT_AGAIN;
01055 
01056   // If the step needs to be adjusted:
01057   if (adjustStep)
01058   {
01059     newTimeStep = max(newTimeStep, minTimeStep);
01060     newTimeStep = min(newTimeStep, maxTimeStep);
01061 
01062     Scalar nextTimePt = time + newTimeStep;
01063 
01064     if (nextTimePt > stopTime)
01065     {
01066       nextTimePt  = stopTime;
01067       newTimeStep = stopTime - time;
01068     }
01069 
01070     hh = newTimeStep;
01071   }
01072   else // if time step is constant for this step:
01073   {
01074     Scalar nextTimePt = time + hh;
01075 
01076     if (nextTimePt > stopTime)
01077     {
01078       nextTimePt      = stopTime;
01079       hh = stopTime - time;
01080     }
01081     return_status = CONTINUE_ANYWAY;
01082   }
01083 #ifdef Rythmos_DEBUG
01084   if (debugLevel > 1)
01085   {
01086     *debug_out << "hh = " << hh << endl;
01087   }
01088 #endif // Rythmos_DEBUG
01089   return(return_status);
01090 }
01091 
01092 template<class Scalar>
01093 void ImplicitBDFStepper<Scalar>::completeStep()
01094 {
01095   typedef Teuchos::ScalarTraits<Scalar> ST;
01096 
01097   numberOfSteps ++;
01098   nef = 0;
01099   time += hh;
01100 #ifdef Rythmos_DEBUG
01101   Teuchos::OSTab ostab(debug_out,1,"completeStep");
01102   if (debugLevel > 1)
01103   {
01104     *debug_out << "numberOfSteps = " << numberOfSteps << endl;
01105     *debug_out << "nef = " << nef << endl;
01106     *debug_out << "time = " << time << endl;
01107   }
01108 #endif // Rythmos_DEBUG
01109   
01110   // Only update the time step if we are NOT running constant stepsize.
01111   bool adjustStep = (!constantStepSize);
01112 
01113   Scalar newTimeStep = hh;
01114   Scalar rr = ST::one(); // step size ratio = new step / old step
01115   // 03/11/04 tscoffe:  Here is the block for choosing order & step-size when
01116   // the local error test PASSES (and Newton succeeded). 
01117   int orderDiff = currentOrder - usedOrder;
01118   usedOrder = currentOrder;
01119   usedStep = hh;
01120   if ((newOrder == currentOrder-1) || (currentOrder == maxOrder))
01121   {
01122     // If we reduced our order or reached max order then move to the next phase
01123     // of integration where we don't automatically double the step-size and
01124     // increase the order.
01125     initialPhase = false;
01126   }
01127 #ifdef Rythmos_DEBUG
01128   if (debugLevel > 1)
01129   {
01130     *debug_out << "initialPhase = " << initialPhase << endl;
01131   }
01132 #endif // Rythmos_DEBUG
01133   if (initialPhase)
01134   {
01135     currentOrder++;
01136     newTimeStep = h_phase0_incr * hh;
01137 #ifdef Rythmos_DEBUG
01138     if (debugLevel > 1)
01139     {
01140       *debug_out << "currentOrder = " << currentOrder << endl;
01141       *debug_out << "newTimeStep = " << newTimeStep << endl;
01142     }
01143 #endif // Rythmos_DEBUG
01144   }
01145   else // not in the initial phase of integration
01146   {
01147     BDFactionFlag action = ACTION_UNSET;
01148     if (newOrder == currentOrder-1)
01149       action = ACTION_LOWER;
01150     else if (newOrder == maxOrder)
01151       action = ACTION_MAINTAIN;
01152     else if ((currentOrder+1>=nscsco) || (orderDiff == 1))
01153     {
01154       // If we just raised the order last time then we won't raise it again
01155       // until we've taken currentOrder+1 steps at order currentOrder.
01156       action = ACTION_MAINTAIN;
01157     }
01158     else // consider changing the order 
01159     {
01160       V_StVpStV(&*delta,ST::one(),*ee,Scalar(-ST::one()),*xHistory[currentOrder+1]);
01161       Tkp1 = WRMSNorm(*errWtVec,*delta);
01162       Ekp1 = Tkp1/(currentOrder+2);
01163 #ifdef Rythmos_DEBUG
01164       if (debugLevel > 1)
01165       {
01166         *debug_out << "delta = " << endl;
01167         delta->describe(*debug_out,Teuchos::VERB_EXTREME);
01168         *debug_out << "Tkp1 = ||delta||_WRMS = " << Tkp1 << endl;
01169         *debug_out << "Ekp1 = " << Ekp1 << endl;
01170       }
01171 #endif // Rythmos_DEBUG
01172       if (currentOrder == 1)
01173       {
01174         if (Tkp1 >= Tkp1_Tk_safety * Tk)
01175           action = ACTION_MAINTAIN;
01176         else
01177           action = ACTION_RAISE;
01178       }
01179       else
01180       {
01181         if (Tkm1 <= min(Tk,Tkp1))
01182           action = ACTION_LOWER;
01183         else if (Tkp1 >= Tk)
01184           action = ACTION_MAINTAIN;
01185         else
01186           action = ACTION_RAISE;
01187       }
01188     }
01189     if (action == ACTION_RAISE)
01190     {
01191       currentOrder++;
01192       Est = Ekp1;
01193     }
01194     else if (action == ACTION_LOWER)
01195     {
01196       currentOrder--;
01197       Est = Ekm1;
01198     }
01199     newTimeStep = hh;
01200     rr = pow(r_safety*Est+r_fudge,-1.0/(currentOrder+1.0));
01201     if (rr >= r_hincr_test)
01202     {
01203       rr = r_hincr;
01204       newTimeStep = rr*hh;
01205     }
01206     else if (rr <= 1)
01207     {
01208       rr = max(r_min,min(r_max,rr));
01209       newTimeStep = rr*hh;
01210     }
01211 #ifdef Rythmos_DEBUG
01212     if (debugLevel > 1)
01213     {
01214       *debug_out << "Est = " << Est << endl;
01215       *debug_out << "currentOrder = " << currentOrder << endl;
01216       *debug_out << "rr  = " << rr << endl;
01217       *debug_out << "newTimeStep = " << newTimeStep << endl;
01218     }
01219 #endif // Rythmos_DEBUG
01220   }
01221   // 03/22/04 tscoffe:  Note that updating the history has nothing to do with
01222   // the step-size and everything to do with the newton correction vectors.
01223   updateHistory();
01224 
01225   // 12/01/05 tscoffe:  This is necessary to avoid currentTimeStep == 0 right
01226   // before a breakpoint.  So I'm checking to see if time is identically
01227   // equal to stopTime, in which case we are right before a breakpoint and we
01228   // should not adjust currentStepSize because that would result in
01229   // currentStepSize == 0.
01230   if (time < stopTime)
01231   {
01232     // If the step needs to be adjusted:
01233     if (adjustStep)
01234     {
01235       newTimeStep = max(newTimeStep, minTimeStep);
01236       newTimeStep = min(newTimeStep, maxTimeStep);
01237 
01238       Scalar nextTimePt = time + newTimeStep;
01239 
01240       if (nextTimePt > stopTime)
01241       {
01242         nextTimePt  = stopTime;
01243         newTimeStep = stopTime - time;
01244       }
01245 
01246       hh = newTimeStep;
01247     }
01248     else // if time step is constant for this step:
01249     {
01250       Scalar nextTimePt = time + hh;
01251 
01252       if (nextTimePt > stopTime)
01253       {
01254         nextTimePt      = stopTime;
01255         hh = stopTime - time;
01256       }
01257     }
01258   }
01259 #ifdef Rythmos_DEBUG
01260   if (debugLevel > 1)
01261   {
01262     *debug_out << "hh = " << hh << endl;
01263   }
01264 #endif // Rythmos_DEBUG
01265 }
01266 
01267 template<class Scalar>
01268 bool ImplicitBDFStepper<Scalar>::ErrWtVecSet(Thyra::VectorBase<Scalar> *w_in, const Thyra::VectorBase<Scalar> &y)
01269 {
01270   typedef Teuchos::ScalarTraits<Scalar> ST;
01271   Thyra::VectorBase<Scalar> &w = *w_in;
01272   abs(&w,y);
01273   Vt_S(&w,relErrTol);
01274   Vp_S(&w,absErrTol);
01275   reciprocal(&w,w);
01276   Vt_StV(&w,ST::one(),w); // We square w because of how weighted norm_2 is computed.
01277   // divide by N to get RMS norm
01278   int N = y.space()->dim();
01279   Vt_S(&w,Scalar(1.0/N));
01280   // Now you can compute WRMS norm as:
01281   // Scalar WRMSnorm = norm_2(w,y); // WRMS norm of y with respect to weights w.
01282   return true; // This should be updated to reflect success of reciprocal
01283 }
01284 
01285 template<class Scalar>
01286 Scalar ImplicitBDFStepper<Scalar>::WRMSNorm(const Thyra::VectorBase<Scalar> &w, const Thyra::VectorBase<Scalar> &y) const
01287 {
01288   return(norm_2(w,y));
01289 }
01290 
01291 template<class Scalar>
01292 bool ImplicitBDFStepper<Scalar>::SetPoints(
01293     const std::vector<Scalar>& time_list
01294     ,const std::vector<Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > >& x_list
01295     ,const std::vector<Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > >& xdot_list)
01296 {
01297   return(false);
01298 }
01299 
01300 template<class Scalar>
01301 bool ImplicitBDFStepper<Scalar>::GetPoints(
01302     const std::vector<Scalar>& time_list
01303     ,std::vector<Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > >* x_list
01304     ,std::vector<Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > >* xdot_list
01305     ,std::vector<ScalarMag>* accuracy_list) const
01306 {
01307   return(false);
01308 }
01309 
01310 template<class Scalar>
01311 bool ImplicitBDFStepper<Scalar>::SetRange(
01312     const Scalar& time_lower 
01313     ,const Scalar& time_upper
01314     ,const InterpolationBuffer<Scalar>& IB)
01315 {
01316   return(false);
01317 }
01318 
01319 template<class Scalar>
01320 bool ImplicitBDFStepper<Scalar>::GetNodes(std::vector<Scalar>* time_list) const
01321 {
01322   return(false);
01323 }
01324 
01325 template<class Scalar>
01326 bool ImplicitBDFStepper<Scalar>::RemoveNodes(std::vector<Scalar>& time_list) const
01327 {
01328   return(false);
01329 }
01330 
01331 template<class Scalar>
01332 int ImplicitBDFStepper<Scalar>::GetOrder() const
01333 {
01334   return(currentOrder);
01335 }
01336 
01337 
01338 } // namespace Rythmos
01339 
01340 #endif //Rythmos_IMPLICITBDF_STEPPER_H

Generated on Thu Sep 18 12:30:05 2008 for Rythmos - Transient Integration for Differential Equations by doxygen 1.3.9.1