Rythmos_ImplicitBDFStepper.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_H
00030 #define Rythmos_IMPLICITBDF_STEPPER_H
00031 
00032 #include "Rythmos_ImplicitBDFStepperDecl.hpp"
00033 #include "Rythmos_ImplicitBDFStepperStepControl.hpp"
00034 
00035 namespace {
00036   const std::string RythmosStepControlSettings_name = "Step Control Settings";
00037 }
00038 
00039 namespace Rythmos {
00040 
00041 // ////////////////////////////
00042 // Defintions
00043 
00044 
00045 // Constructors, intializers, Misc.
00046 
00047 
00048 template<class Scalar>
00049 ImplicitBDFStepper<Scalar>::ImplicitBDFStepper()
00050 {
00051   defaultInitializeAllData_();
00052   haveInitialCondition_ = false;
00053   isInitialized_=false;
00054 }
00055 
00056 
00057 template<class Scalar>
00058 ImplicitBDFStepper<Scalar>::ImplicitBDFStepper(
00059   const RCP<const Thyra::ModelEvaluator<Scalar> > &model
00060   ,const RCP<Thyra::NonlinearSolverBase<Scalar> > &solver
00061   ,RCP<Teuchos::ParameterList> &parameterList
00062   )
00063 {
00064   defaultInitializeAllData_();
00065   this->setParameterList(parameterList);
00066   // Now we instantiate the model and the solver
00067   setModel(model);
00068   setSolver(solver);
00069   haveInitialCondition_ = false;
00070   isInitialized_=false;
00071 }
00072 
00073 
00074 template<class Scalar>
00075 ImplicitBDFStepper<Scalar>::ImplicitBDFStepper(
00076   const RCP<const Thyra::ModelEvaluator<Scalar> > &model
00077   ,const RCP<Thyra::NonlinearSolverBase<Scalar> > &solver
00078   )
00079 {
00080   defaultInitializeAllData_();
00081   // Now we instantiate the model and the solver
00082   setModel(model);
00083   setSolver(solver);
00084   haveInitialCondition_ = false;
00085   isInitialized_=false;
00086 }
00087 
00088 template<class Scalar>
00089 const Thyra::VectorBase<Scalar>& 
00090   ImplicitBDFStepper<Scalar>::getxHistory(int index) const
00091 {
00092   TEST_FOR_EXCEPTION(!isInitialized_,std::logic_error,
00093       "Error, attempting to call getxHistory before initialization!\n");
00094   TEST_FOR_EXCEPT( !(( 0 <= index ) && ( index <= maxOrder_ )) );
00095   TEST_FOR_EXCEPT( !( index <= usedOrder_+1 ) );
00096   return(*(xHistory_[index]));
00097 }
00098 
00099 template<class Scalar>
00100 void ImplicitBDFStepper<Scalar>::setStepControlStrategy(const RCP<StepControlStrategyBase<Scalar> >& stepControl)
00101 {
00102   TEST_FOR_EXCEPTION(stepControl == Teuchos::null,std::logic_error,"Error, stepControl == Teuchos::null!\n");
00103   stepControl_ = stepControl;    
00104 }
00105 
00106 template<class Scalar>
00107 RCP<StepControlStrategyBase<Scalar> > ImplicitBDFStepper<Scalar>::getNonconstStepControlStrategy() 
00108 {
00109   return(stepControl_);
00110 }
00111 
00112 template<class Scalar>
00113 RCP<const StepControlStrategyBase<Scalar> > ImplicitBDFStepper<Scalar>::getStepControlStrategy() const
00114 {
00115   return(stepControl_);
00116 }
00117 
00118 
00119 // Overridden from SolverAcceptingStepperBase
00120 
00121 
00122 template<class Scalar>
00123 void ImplicitBDFStepper<Scalar>::setSolver(const RCP<Thyra::NonlinearSolverBase<Scalar> > &solver)
00124 {
00125   TEST_FOR_EXCEPT(solver == Teuchos::null)
00126     solver_ = solver;
00127 }
00128 
00129 
00130 template<class Scalar>
00131 RCP<Thyra::NonlinearSolverBase<Scalar> >
00132 ImplicitBDFStepper<Scalar>::getNonconstSolver()
00133 {
00134   return (solver_);
00135 }
00136 
00137 
00138 template<class Scalar>
00139 RCP<const Thyra::NonlinearSolverBase<Scalar> >
00140 ImplicitBDFStepper<Scalar>::getSolver() const
00141 {
00142   return (solver_);
00143 }
00144 
00145 
00146 // Overridden from StepperBase
00147 
00148 
00149 template<class Scalar>
00150 bool ImplicitBDFStepper<Scalar>::supportsCloning() const
00151 {
00152   return true;
00153 }
00154 
00155 
00156 template<class Scalar>
00157 RCP<StepperBase<Scalar> >
00158 ImplicitBDFStepper<Scalar>::cloneStepperAlgorithm() const
00159 {
00160 
00161   // Just use the interface to clone the algorithm in an basically
00162   // uninitialized state
00163 
00164   RCP<ImplicitBDFStepper<Scalar> >
00165     stepper = Teuchos::rcp(new ImplicitBDFStepper<Scalar>());
00166 
00167   if (!is_null(model_))
00168     stepper->setModel(model_); // Shallow copy is okay!
00169 
00170   if (!is_null(solver_))
00171     stepper->setSolver(solver_->cloneNonlinearSolver().assert_not_null());
00172   
00173   if (!is_null(parameterList_))
00174     stepper->setParameterList(Teuchos::parameterList(*parameterList_));
00175 
00176   if (!is_null(stepControl_)) {
00177     if (stepControl_->supportsCloning())
00178       stepper->setStepControlStrategy(
00179         stepControl_->cloneStepControlStrategyAlgorithm().assert_not_null());
00180   }
00181   
00182   // At this point, we should have a valid algorithm state.  What might be
00183   // missing is the initial condition if it was not given in *model_ but was
00184   // set explicitly.  However, the specification for this function does not
00185   // guarantee that the full state will be copied in any case!
00186 
00187   return stepper;
00188 
00189 }
00190 
00191 
00192 template<class Scalar>
00193 void ImplicitBDFStepper<Scalar>::setModel(
00194   const RCP<const Thyra::ModelEvaluator<Scalar> > &model
00195   )
00196 {
00197   typedef Teuchos::ScalarTraits<Scalar> ST;
00198   TEST_FOR_EXCEPT(model == Teuchos::null)
00199     model_= model;
00200 }
00201 
00202 
00203 template<class Scalar>
00204 RCP<const Thyra::ModelEvaluator<Scalar> >
00205 ImplicitBDFStepper<Scalar>::getModel() const
00206 {
00207   return model_;
00208 }
00209 
00210 
00211 template<class Scalar>
00212 void ImplicitBDFStepper<Scalar>::setInitialCondition(
00213   const Thyra::ModelEvaluatorBase::InArgs<Scalar> &initialCondition
00214   )
00215 {
00216   typedef Teuchos::ScalarTraits<Scalar> ST;
00217   typedef Thyra::ModelEvaluatorBase MEB;
00218   TEST_FOR_EXCEPT(initialCondition.get_x()==Teuchos::null);
00219   TEST_FOR_EXCEPT(initialCondition.get_x_dot()==Teuchos::null);
00220   basePoint_ = initialCondition;
00221   xn0_ = initialCondition.get_x()->clone_v();
00222   xpn0_ = initialCondition.get_x_dot()->clone_v(); 
00223   time_ = initialCondition.get_t();
00224   haveInitialCondition_ = true;
00225   if (isInitialized_) {
00226     initialize_();
00227   }
00228 }
00229 
00230 
00231 template<class Scalar>
00232 Scalar ImplicitBDFStepper<Scalar>::takeStep(Scalar dt, StepSizeType stepType)
00233 {
00234   
00235 #ifdef ENABLE_RYTHMOS_TIMERS
00236   TEUCHOS_FUNC_TIME_MONITOR("Rythmos::ImplicitBDFStepper::takeStep");
00237 #endif
00238   
00239   using Teuchos::as;
00240   using Teuchos::incrVerbLevel;
00241   typedef Teuchos::ScalarTraits<Scalar> ST;
00242   typedef typename Thyra::ModelEvaluatorBase::InArgs<Scalar>::ScalarMag ScalarMag;
00243   typedef Thyra::NonlinearSolverBase<Scalar> NSB;
00244   typedef Teuchos::VerboseObjectTempState<NSB> VOTSNSB;
00245 
00246   RCP<Teuchos::FancyOStream> out = this->getOStream();
00247   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00248   Teuchos::OSTab ostab(out,1,"takeStep");
00249   VOTSNSB solver_outputTempState(solver_,out,incrVerbLevel(verbLevel,-1));
00250 
00251   if ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
00252     *out
00253       << "\nEntering " << this->Teuchos::Describable::description()
00254       << "::takeStep("<<dt<<","<<toString(stepType)<<") ...\n"; 
00255   }
00256 
00257   if (!isInitialized_) {
00258     initialize_(); 
00259   }
00260 
00261   stepControl_->setOStream(out);
00262   stepControl_->setVerbLevel(verbLevel);
00263   stepControl_->setRequestedStepSize(*this,dt,stepType);
00264 
00265   AttemptedStepStatusFlag status;
00266   while (1) {
00267     // Set up problem coefficients (and handle first step)
00268     Scalar hh_old = hh_;
00269     int desiredOrder;
00270     stepControl_->nextStepSize(*this,&hh_,&stepType,&desiredOrder);
00271     TEST_FOR_EXCEPT(!((1 <= desiredOrder) && (desiredOrder <= maxOrder_)));
00272     TEST_FOR_EXCEPT(!(desiredOrder <= usedOrder_+1));
00273     currentOrder_ = desiredOrder;
00274     if (numberOfSteps_ == 0) {
00275       psi_[0] = hh_;
00276       if (nef_ == 0) {
00277         Vt_S(&*xHistory_[1],hh_);
00278       } else {
00279         Vt_S(&*xHistory_[1],hh_/hh_old);
00280       }
00281     }
00282     this->updateCoeffs_();
00283     // compute predictor
00284     obtainPredictor_();
00285     // solve nonlinear problem (as follows)
00286     
00287     //
00288     // Setup the nonlinear equations:
00289     //
00290     //   f_bar( x_dot_coeff * x_bar + x_dot_base, x_coeff * x_bar + x_base, t_base ) = 0
00291     //   x_dot_coeff = -alpha_s/dt
00292     //   x_dot_base = x_prime_pred + (alpha_s/dt) * x_pred
00293     //   x_coeff = 1
00294     //   x_base = 0
00295     //   t_base = tn+dt
00296     //
00297     Scalar coeff_x_dot = Scalar(-ST::one())*alpha_s_/hh_;
00298     V_StVpStV( &*x_dot_base_, ST::one(), *xpn0_, alpha_s_/hh_, *xn0_ );
00299     if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_EXTREME) ) {
00300       *out << "model_ = " << std::endl;
00301       model_->describe(*out,verbLevel);
00302       *out << "basePoint_ = " << std::endl;
00303       basePoint_.describe(*out,verbLevel);
00304       *out << "coeff_x_dot = " << coeff_x_dot << std::endl;
00305       *out << "x_dot_base_ = " << std::endl;
00306       x_dot_base_->describe(*out,verbLevel);
00307       *out << "time_+hh_ = " << time_+hh_ << std::endl;
00308       *out << "xn0_ = " << std::endl;
00309       xn0_->describe(*out,verbLevel);
00310     }
00311     neModel_.initializeSingleResidualModel(
00312       model_, basePoint_,
00313       coeff_x_dot, x_dot_base_,
00314       ST::one(), Teuchos::null,
00315       time_+hh_,
00316       xn0_
00317       );
00318     //
00319     // Solve the implicit nonlinear system to a tolerance of ???
00320     // 
00321     if(solver_->getModel().get()!=&neModel_) {
00322       solver_->setModel( Teuchos::rcp(&neModel_,false) );
00323     }
00324     /* // Rythmos::TimeStepNonlinearSolver uses a built in solveCriteria, so you can't pass one in.
00325     // I believe this is the correct solveCriteria for IDA though.
00326     Thyra::SolveMeasureType nonlinear_solve_measure_type(Thyra::SOLVE_MEASURE_NORM_RESIDUAL,Thyra::SOLVE_MEASURE_ONE); 
00327     ScalarMag tolerance = relErrTol_/ScalarMag(10.0); // This should be changed to match the condition in IDA
00328     Thyra::SolveCriteria<Scalar> nonlinearSolveCriteria(nonlinear_solve_measure_type, tolerance);
00329     Thyra::SolveStatus<Scalar> nonlinearSolveStatus = solver_->solve( &*xn0_, &nonlinearSolveCriteria, &*delta_ ); 
00330     */
00331     if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_EXTREME) ) {
00332       *out << "xn0 = " << std::endl;
00333       xn0_->describe(*out,verbLevel);
00334       *out << "ee = " << std::endl;
00335       ee_->describe(*out,verbLevel);
00336     }
00337     Thyra::SolveStatus<Scalar> nonlinearSolveStatus = solver_->solve( &*xn0_, NULL, &*ee_ ); 
00338     if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_EXTREME) ) {
00339       *out << "xn0 = " << std::endl;
00340       xn0_->describe(*out,verbLevel);
00341       *out << "ee = " << std::endl;
00342       ee_->describe(*out,verbLevel);
00343     }
00344     // In the above solve, on input *xn0_ is the initial guess that comes from
00345     // the predictor.  On output, *xn0_ is the solved for solution and *ee_ is
00346     // the difference computed from the intial guess in *xn0_ to the final
00347     // solved value of *xn0_.  This is needed for basic numerical stability.
00348     if (nonlinearSolveStatus.solveStatus == Thyra::SOLVE_STATUS_CONVERGED)  {
00349       newtonConvergenceStatus_ = 0;
00350     }
00351     else {
00352       newtonConvergenceStatus_ = -1;
00353     }
00354 
00355     // check error and evaluate LTE
00356     stepControl_->setCorrection(*this,xn0_,ee_,newtonConvergenceStatus_);
00357     bool stepPass = stepControl_->acceptStep(*this,&LETvalue_);
00358     
00359     if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00360       *out << "xn0_ = " << std::endl;
00361       xn0_->describe(*out,verbLevel);
00362       *out << "xpn0_ = " << std::endl;
00363       xpn0_->describe(*out,verbLevel);
00364       *out << "ee_ = " << std::endl;
00365       ee_->describe(*out,verbLevel);
00366       for (int i=0; i<std::max(2,maxOrder_); ++i) {
00367         *out << "xHistory_[" << i << "] = " << std::endl;
00368         xHistory_[i]->describe(*out,verbLevel);
00369       }
00370     }
00371 
00372     // Check LTE here:
00373     if (!stepPass) { // stepPass = false
00374       stepLETStatus_ = STEP_LET_STATUS_FAILED;
00375       status = stepControl_->rejectStep(*this);
00376       nef_++;
00377       if (status == CONTINUE_ANYWAY) {
00378         break;
00379       } else {
00380         restoreHistory_();
00381       }
00382     } else { // stepPass = true
00383       stepLETStatus_ = STEP_LET_STATUS_PASSED;
00384       break;
00385     }
00386   }
00387 
00388   // 08/22/07 the history array must be updated before stepControl_->completeStep.
00389   completeStep_();   
00390   stepControl_->completeStep(*this);
00391 
00392   if ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
00393     *out
00394       << "\nLeaving " << this->Teuchos::Describable::description()
00395       << "::takeStep("<<dt<<","<<toString(stepType)<<") ...\n"; 
00396   }
00397 
00398   return(usedStep_);
00399 
00400 }
00401 
00402 
00403 template<class Scalar>
00404 const StepStatus<Scalar> ImplicitBDFStepper<Scalar>::getStepStatus() const
00405 {
00406 
00407   // 2007/08/24: rabartl: We agreed that getStepStatus() would be free
00408   // so I have commented out removed all code that is not free
00409 
00410   typedef Teuchos::ScalarTraits<Scalar> ST;
00411   StepStatus<Scalar> stepStatus;
00412   if (!isInitialized_) {
00413     stepStatus.message = "This stepper is uninitialized.";
00414     stepStatus.stepStatus = STEP_STATUS_UNINITIALIZED;
00415     stepStatus.stepSize = Scalar(-ST::one());
00416     stepStatus.order = -1;
00417     stepStatus.time = Scalar(-ST::one());
00418     return(stepStatus);
00419   }
00420 
00421   if (numberOfSteps_ > 0) {
00422     stepStatus.stepStatus = STEP_STATUS_CONVERGED; 
00423   } else {
00424     stepStatus.stepStatus = STEP_STATUS_UNKNOWN;
00425   }
00426   stepStatus.stepLETStatus = stepLETStatus_;
00427   stepStatus.stepSize = usedStep_; 
00428   stepStatus.order = usedOrder_;
00429   stepStatus.time = time_;
00430   stepStatus.stepLETValue = LETvalue_; 
00431   stepStatus.solution = xHistory_[0];
00432   stepStatus.solutionDot = Teuchos::null; // This is *not* free!
00433   stepStatus.residual = Teuchos::null; // This is *not* free!
00434 
00435   return(stepStatus);
00436 
00437 }
00438 
00439 
00440 // Overridden from InterpolationBufferBase
00441 
00442 
00443 template<class Scalar>
00444 RCP<const Thyra::VectorSpaceBase<Scalar> >
00445 ImplicitBDFStepper<Scalar>::get_x_space() const
00446 {
00447   //TEST_FOR_EXCEPTION(!isInitialized_,std::logic_error,"Error, attempting to call get_x_space before initialization!\n");
00448   return ( !is_null(model_) ? model_->get_x_space() : Teuchos::null );
00449 }
00450 
00451 
00452 template<class Scalar>
00453 void ImplicitBDFStepper<Scalar>::addPoints(
00454   const Array<Scalar>& time_vec,
00455   const Array<RCP<const Thyra::VectorBase<Scalar> > >& x_vec,
00456   const Array<RCP<const Thyra::VectorBase<Scalar> > >& xdot_vec
00457   )
00458 {
00459   TEST_FOR_EXCEPTION(true,std::logic_error,
00460     "Error, addPoints is not implemented for ImplicitBDFStepper.\n");
00461 }
00462 
00463 
00464 template<class Scalar>
00465 TimeRange<Scalar> ImplicitBDFStepper<Scalar>::getTimeRange() const
00466 {
00467   if ( !isInitialized_ && haveInitialCondition_ )
00468     return timeRange<Scalar>(time_,time_);
00469   if ( !isInitialized_ && !haveInitialCondition_ )
00470     return invalidTimeRange<Scalar>();
00471   return timeRange<Scalar>(time_-usedStep_,time_);
00472 }
00473 
00474 
00475 template<class Scalar>
00476 void ImplicitBDFStepper<Scalar>::getPoints(
00477   const Array<Scalar>& time_vec
00478   ,Array<RCP<const Thyra::VectorBase<Scalar> > >* x_vec
00479   ,Array<RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec
00480   ,Array<ScalarMag>* accuracy_vec) const
00481 {
00482   using Teuchos::as;
00483   typedef Teuchos::ScalarTraits<Scalar> ST;
00484   typedef typename ST::magnitudeType ScalarMag;
00485 #ifdef ENABLE_RYTHMOS_TIMERS
00486   TEUCHOS_FUNC_TIME_MONITOR("Rythmos::ImplicitBDFStepper::getPoints");
00487 #endif
00488   if (x_vec)
00489     x_vec->clear();
00490   if (xdot_vec)
00491     xdot_vec->clear();
00492   for (unsigned int i=0 ; i<time_vec.size() ; ++i) {
00493     RCP<Thyra::VectorBase<Scalar> >
00494       x_temp = createMember(xn0_->space());
00495     RCP<Thyra::VectorBase<Scalar> >
00496       xdot_temp = createMember(xn0_->space());
00497     ScalarMag accuracy = -ST::zero();
00498     interpolateSolution_(
00499       time_vec[i], &*x_temp, &*xdot_temp,
00500       accuracy_vec ? &accuracy : 0
00501       );
00502     if (x_vec)
00503       x_vec->push_back(x_temp);
00504     if (xdot_vec)
00505       xdot_vec->push_back(xdot_temp);
00506     if (accuracy_vec)
00507       accuracy_vec->push_back(accuracy);
00508   }
00509   if ( as<int>(this->getVerbLevel()) >= as<int>(Teuchos::VERB_HIGH) ) {
00510     RCP<Teuchos::FancyOStream> out = this->getOStream();
00511     Teuchos::OSTab ostab(out,1,"getPoints");
00512     *out << "Passing out the interpolated values:" << std::endl;
00513     for (unsigned int i=0; i<time_vec.size() ; ++i) {
00514       *out << "time_[" << i << "] = " << time_vec[i] << std::endl;
00515       if (x_vec) {
00516         *out << "x_vec[" << i << "] = " << std::endl;
00517         (*x_vec)[i]->describe(*out,this->getVerbLevel());
00518       }
00519       if (xdot_vec) {
00520         *out << "xdot_vec[" << i << "] = ";
00521         if ( (*xdot_vec)[i] == Teuchos::null) {
00522           *out << "Teuchos::null" << std::endl;
00523         }
00524         else {
00525           *out << std::endl << Teuchos::describe(*(*xdot_vec)[i],this->getVerbLevel());
00526         }
00527       }
00528       if (accuracy_vec)
00529         *out << "accuracy[" << i << "] = " << (*accuracy_vec)[i] << std::endl;
00530     }
00531   }
00532 }
00533 
00534 
00535 template<class Scalar>
00536 void ImplicitBDFStepper<Scalar>::getNodes(Array<Scalar>* time_vec) const
00537 {
00538   TEST_FOR_EXCEPTION(!isInitialized_,std::logic_error,
00539     "Error, attempting to call getNodes before initialization!\n");
00540   if (numberOfSteps_ > 0) {
00541     time_vec->push_back(time_-usedStep_);
00542   }
00543   time_vec->push_back(time_);
00544 }
00545 
00546 
00547 template<class Scalar>
00548 void ImplicitBDFStepper<Scalar>::removeNodes(Array<Scalar>& time_vec) 
00549 {
00550   TEST_FOR_EXCEPTION(true,std::logic_error,
00551     "Error, removeNodes is not implemented for ImplicitBDFStepper.\n");
00552 }
00553 
00554 
00555 template<class Scalar>
00556 int ImplicitBDFStepper<Scalar>::getOrder() const
00557 {
00558   if (!isInitialized_) {
00559     return(-1);
00560   }
00561   return(usedOrder_);
00562 }
00563 
00564 
00565 // Overridden from Teuchos::ParameterListAcceptor
00566 
00567 
00568 template<class Scalar>
00569 void ImplicitBDFStepper<Scalar>::setParameterList(
00570   RCP<Teuchos::ParameterList> const& paramList
00571   )
00572 {
00573   TEST_FOR_EXCEPT(paramList == Teuchos::null);
00574   paramList->validateParameters(*this->getValidParameters(),0);
00575   parameterList_ = paramList;
00576   Teuchos::readVerboseObjectSublist(&*parameterList_,this);
00577 }
00578 
00579 
00580 template<class Scalar>
00581 RCP<Teuchos::ParameterList> ImplicitBDFStepper<Scalar>::getParameterList()
00582 {
00583   return(parameterList_);
00584 }
00585 
00586 
00587 template<class Scalar>
00588 RCP<Teuchos::ParameterList>
00589 ImplicitBDFStepper<Scalar>::unsetParameterList()
00590 {
00591   RCP<Teuchos::ParameterList> temp_param_list = parameterList_;
00592   parameterList_ = Teuchos::null;
00593   return(temp_param_list);
00594 }
00595 
00596 
00597 template<class Scalar>
00598 RCP<const Teuchos::ParameterList>
00599 ImplicitBDFStepper<Scalar>::getValidParameters() const
00600 {
00601 
00602   static RCP<Teuchos::ParameterList> validPL;
00603 
00604   if (is_null(validPL)) {
00605 
00606     RCP<Teuchos::ParameterList>
00607       pl = Teuchos::parameterList();
00608 
00609     pl->sublist(RythmosStepControlSettings_name);
00610 
00611     Teuchos::setupVerboseObjectSublist(&*pl);
00612 
00613     validPL = pl;
00614 
00615   }
00616 
00617   RCP<Teuchos::FancyOStream> out = this->getOStream();
00618   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00619   Teuchos::OSTab ostab(out,1,"getValidParameters");
00620   if (Teuchos::as<int>(verbLevel) == Teuchos::VERB_HIGH) {
00621     *out << "Setting up valid parameterlist." << std::endl;
00622     validPL->print(*out);
00623   }
00624 
00625   return (validPL);
00626   
00627 }
00628 
00629 
00630 // Overridden from Teuchos::Describable
00631 
00632 
00633 template<class Scalar>
00634 std::string ImplicitBDFStepper<Scalar>::description() const
00635 {
00636   std::ostringstream out;
00637   out << this->Teuchos::Describable::description();
00638   const TimeRange<Scalar> timeRange = this->getTimeRange();
00639   if (timeRange.isValid())
00640     out << " (timeRange="<<timeRange<<")";
00641   else
00642     out << " (This stepper is not initialized yet)";
00643   out << std::endl;
00644   return out.str();
00645 }
00646 
00647 
00648 template<class Scalar>
00649 void ImplicitBDFStepper<Scalar>::describe(
00650   Teuchos::FancyOStream &out,
00651   const Teuchos::EVerbosityLevel verbLevel
00652   ) const
00653 {
00654 
00655   using Teuchos::as;
00656 
00657   if (!isInitialized_) {
00658     out << this->description();
00659     return;
00660   }
00661   
00662   if ( (as<int>(verbLevel) == as<int>(Teuchos::VERB_DEFAULT) ) ||
00663     (as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW)     )
00664     )
00665   {
00666     out << this->description() << std::endl;
00667     out << "model_ = " << Teuchos::describe(*model_,verbLevel);
00668     out << "solver_ = " << Teuchos::describe(*solver_,verbLevel);
00669     out << "neModel_ = " << Teuchos::describe(neModel_,verbLevel);
00670   }
00671   if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW)) {
00672     out << "time_ = " << time_ << std::endl;
00673     out << "hh_ = " << hh_ << std::endl;
00674     out << "currentOrder_ = " << currentOrder_ << std::endl;
00675   }
00676   if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH)) {
00677     out << "xn0_ = " << Teuchos::describe(*xn0_,verbLevel);
00678     out << "xpn0_ = " << Teuchos::describe(*xpn0_,verbLevel);
00679     out << "x_dot_base_ = " << Teuchos::describe(*x_dot_base_,verbLevel);
00680     for (int i=0 ; i < std::max(2,maxOrder_) ; ++i) {
00681       out << "xHistory_[" << i << "] = "
00682           << Teuchos::describe(*xHistory_[i],verbLevel);
00683     }
00684     out << "ee_ = " << Teuchos::describe(*ee_,verbLevel);
00685   }
00686 }
00687 
00688 
00689 // private
00690 
00691 
00692 // 2007/08/24: rabartl: Belos: We really should initialize all of this data in
00693 // a member initialization list but since there are like three constructors
00694 // this would mean that we would have to duplicate code (which is error prone)
00695 // or use a macro (which is not easy to debug).  We really should remove all
00696 // but the default constructor which then would set this data once in the
00697 // initialization list.
00698 
00699 template<class Scalar>
00700 void ImplicitBDFStepper<Scalar>::defaultInitializeAllData_()
00701 {
00702   typedef Teuchos::ScalarTraits<Scalar> ST;
00703   const Scalar nan = ST::nan(), one = ST::one(), zero = ST::zero();
00704   // Initialize some data members to their rightful default values
00705   haveInitialCondition_ = false;
00706   isInitialized_ = false;
00707   currentOrder_ = 1;
00708   usedOrder_ = 1;
00709   usedStep_ = zero;
00710   // Initialize the rest of the private data members to invalid values to
00711   // avoid uninitialed memory
00712   time_ = nan;
00713   hh_ = nan;
00714   maxOrder_ = -1;
00715   LETvalue_ = -one;
00716   stepLETStatus_ = STEP_LET_STATUS_UNKNOWN;
00717   alpha_s_ = -one;
00718   numberOfSteps_ = -1;
00719   nef_ = -1;
00720   nscsco_ = -1;
00721   newtonConvergenceStatus_ = -1;
00722 }
00723 
00724 
00725 template<class Scalar>
00726 void ImplicitBDFStepper<Scalar>::getInitialCondition_()
00727 {
00728   typedef Teuchos::ScalarTraits<Scalar> ST;
00729   if (!haveInitialCondition_) {
00730     TEST_FOR_EXCEPT(model_->getNominalValues().get_x()==Teuchos::null);
00731     TEST_FOR_EXCEPT(model_->getNominalValues().get_x_dot()==Teuchos::null);
00732     basePoint_ = model_->getNominalValues();
00733     xn0_ = basePoint_.get_x()->clone_v();
00734     xpn0_ = basePoint_.get_x_dot()->clone_v(); 
00735     time_ = basePoint_.get_t();
00736     haveInitialCondition_ = true;
00737   }
00738 }
00739 
00740 
00741 template<class Scalar>
00742 void ImplicitBDFStepper<Scalar>::obtainPredictor_()
00743 {
00744 
00745   using Teuchos::as;
00746   typedef Teuchos::ScalarTraits<Scalar> ST;
00747 
00748   if (!isInitialized_) {
00749     return;
00750   }
00751 
00752   RCP<Teuchos::FancyOStream> out = this->getOStream();
00753   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00754   Teuchos::OSTab ostab(out,1,"obtainPredictor_");
00755   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00756     *out << "currentOrder_ = " << currentOrder_ << std::endl;
00757   }
00758   
00759   // prepare history array for prediction
00760   for (int i=nscsco_;i<=currentOrder_;++i) {
00761     Vt_S(&*xHistory_[i],beta_[i]);
00762   }
00763   
00764   // evaluate predictor
00765   V_V(&*xn0_,*xHistory_[0]);
00766   V_S(&*xpn0_,ST::zero());
00767   for (int i=1;i<=currentOrder_;++i) {
00768     Vp_V(&*xn0_,*xHistory_[i]);
00769     Vp_StV(&*xpn0_,gamma_[i],*xHistory_[i]);
00770   }
00771   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00772     *out << "xn0_ = " << std::endl;
00773     xn0_->describe(*out,verbLevel);
00774     *out << "xpn0_ = " << std::endl;
00775     xpn0_->describe(*out,verbLevel);
00776   }
00777 }
00778 
00779 
00780 template<class Scalar>
00781 void ImplicitBDFStepper<Scalar>::interpolateSolution_(
00782   const Scalar& timepoint,
00783   Thyra::VectorBase<Scalar>* x_ptr,
00784   Thyra::VectorBase<Scalar>* xdot_ptr,
00785   ScalarMag* accuracy_ptr
00786   ) const
00787 {
00788 
00789   typedef std::numeric_limits<Scalar> NL;
00790   typedef Teuchos::ScalarTraits<Scalar> ST;
00791 
00792 #ifdef TEUCHOS_DEBUG
00793   TEST_FOR_EXCEPTION(
00794     !isInitialized_,std::logic_error,
00795     "Error, attempting to call interpolateSolution before initialization!\n");
00796   const TimeRange<Scalar> currTimeRange = this->getTimeRange();
00797   TEST_FOR_EXCEPTION(
00798     !currTimeRange.isInRange(timepoint), std::logic_error,
00799     "Error, timepoint = " << timepoint << " is not in the time range "
00800     << currTimeRange << "!" );
00801 #endif
00802 
00803   const Scalar tn = time_;
00804   const int kused = usedOrder_;
00805 
00806   // order of interpolation
00807   int kord = kused;
00808   if ( (kused == 0) || (timepoint == tn) )  {
00809     kord = 1;
00810   }
00811 
00812   // Initialize interploation
00813   Thyra::V_V(x_ptr,*xHistory_[0]);
00814   Thyra::V_S(xdot_ptr,ST::zero());
00815   
00816   // Add history array contributions
00817   const Scalar delt = timepoint - tn;
00818   Scalar c = ST::one(); // coefficient for interpolation of x
00819   Scalar d = ST::zero(); // coefficient for interpolation of xdot
00820   Scalar gam = delt/psi_[0]; // coefficient for interpolation
00821   for (int j=1 ; j <= kord ; ++j) {
00822     d = d*gam + c/psi_[j-1];
00823     c = c*gam;
00824     gam = (delt + psi_[j-1])/psi_[j];
00825     Thyra::Vp_StV(x_ptr,c,*xHistory_[j]);
00826     Thyra::Vp_StV(xdot_ptr,d,*xHistory_[j]);
00827   }
00828 
00829   // Set approximate accuracy
00830   if (accuracy_ptr)
00831     *accuracy_ptr = pow(usedStep_,kord);
00832   
00833 }
00834 
00835 
00836 template<class Scalar>
00837 void ImplicitBDFStepper<Scalar>::updateHistory_()
00838 {
00839 
00840   using Teuchos::as;
00841 
00842   // Save Newton correction for potential order increase on next step.
00843   if (usedOrder_ < maxOrder_)  {
00844     assign( &*xHistory_[usedOrder_+1], *ee_ );
00845   }
00846   // Update history arrays
00847   Vp_V( &*xHistory_[usedOrder_], *ee_ );
00848   for (int j=usedOrder_-1;j>=0;j--) {
00849     Vp_V( &*xHistory_[j], *xHistory_[j+1] );
00850   }
00851   RCP<Teuchos::FancyOStream> out = this->getOStream();
00852   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00853   Teuchos::OSTab ostab(out,1,"updateHistory_");
00854   if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00855     for (int i=0;i<std::max(2,maxOrder_);++i) {
00856       *out << "xHistory_[" << i << "] = " << std::endl;
00857       xHistory_[i]->describe(*out,verbLevel);
00858     }
00859   }
00860 
00861 }
00862 
00863 
00864 template<class Scalar>
00865 void ImplicitBDFStepper<Scalar>::restoreHistory_()
00866 {
00867 
00868   using Teuchos::as;
00869   typedef Teuchos::ScalarTraits<Scalar> ST;
00870 
00871   // undo preparation of history array for prediction
00872   for (int i=nscsco_;i<=currentOrder_;++i) {
00873     Vt_S( &*xHistory_[i], ST::one()/beta_[i] );
00874   }
00875   for (int i=1;i<=currentOrder_;++i) {
00876     psi_[i-1] = psi_[i] - hh_;
00877   }
00878   RCP<Teuchos::FancyOStream> out = this->getOStream();
00879   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00880   Teuchos::OSTab ostab(out,1,"restoreHistory_");
00881   if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00882     for (int i=0;i<maxOrder_;++i) {
00883       *out << "psi_[" << i << "] = " << psi_[i] << std::endl;
00884     }
00885     for (int i=0;i<maxOrder_;++i) {
00886       *out << "xHistory_[" << i << "] = " << std::endl;
00887       xHistory_[i]->describe(*out,verbLevel);
00888     }
00889   }
00890 
00891 } 
00892 
00893 
00894 template<class Scalar>
00895 void ImplicitBDFStepper<Scalar>::updateCoeffs_()
00896 {
00897 
00898   using Teuchos::as;
00899   typedef Teuchos::ScalarTraits<Scalar> ST;
00900 
00901   // If the number of steps taken with constant order and constant stepsize is
00902   // more than the current order + 1 then we don't bother to update the
00903   // coefficients because we've reached a constant step-size formula.  When
00904   // this is is not true, then we update the coefficients for the variable
00905   // step-sizes. 
00906   if ((hh_ != usedStep_) || (currentOrder_ != usedOrder_)) {
00907     nscsco_ = 0;
00908   }
00909   nscsco_ = std::min(nscsco_+1,usedOrder_+2);
00910   if (currentOrder_+1 >= nscsco_) {
00911     beta_[0] = ST::one();
00912     alpha_[0] = ST::one();
00913     Scalar temp1 = hh_;
00914     gamma_[0] = ST::zero();
00915     for (int i=1;i<=currentOrder_;++i) {
00916       Scalar temp2 = psi_[i-1];
00917       psi_[i-1] = temp1;
00918       beta_[i] = beta_[i-1]*psi_[i-1]/temp2;
00919       temp1 = temp2 + hh_;
00920       alpha_[i] = hh_/temp1;
00921       gamma_[i] = gamma_[i-1]+alpha_[i-1]/hh_;
00922     }
00923     psi_[currentOrder_] = temp1;
00924   }
00925   alpha_s_ = ST::zero();
00926   for (int i=0;i<currentOrder_;++i) {
00927     alpha_s_ = alpha_s_ - Scalar(ST::one()/(i+ST::one()));
00928   }
00929   RCP<Teuchos::FancyOStream> out = this->getOStream();
00930   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00931   Teuchos::OSTab ostab(out,1,"updateCoeffs_");
00932   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00933     for (int i=0;i<=maxOrder_;++i) {
00934       *out << "alpha_[" << i << "] = " << alpha_[i] << std::endl;
00935       *out << "beta_[" << i << "] = " << beta_[i] << std::endl;
00936       *out << "gamma_[" << i << "] = " << gamma_[i] << std::endl;
00937       *out << "psi_[" << i << "] = " << psi_[i] << std::endl;
00938       *out << "alpha_s_ = " << alpha_s_ << std::endl;
00939     }
00940   }
00941 }
00942 
00943 
00944 template<class Scalar>
00945 void ImplicitBDFStepper<Scalar>::initialize_()
00946 {
00947 
00948   using Teuchos::as;
00949   typedef Teuchos::ScalarTraits<Scalar> ST;
00950   using Thyra::createMember;
00951 
00952   RCP<Teuchos::FancyOStream> out = this->getOStream();
00953   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00954   const bool doTrace = (as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH));
00955   Teuchos::OSTab ostab(out,1,"initialize_");
00956 
00957   if (doTrace) {
00958     *out
00959       << "\nEntering " << this->Teuchos::Describable::description()
00960       << "::initialize_()...\n";
00961   }
00962 
00963   TEST_FOR_EXCEPT(model_ == Teuchos::null);
00964   TEST_FOR_EXCEPT(solver_ == Teuchos::null);
00965 
00966   // 08/30/07 tscoffe:
00967   // This must be done before stepControl_->initialize() in order for
00968   // stepControl_ to get a valid TimeRange from the stepper.
00969   this->getInitialCondition_();
00970 
00971   // Initialize Parameter List if none provided.
00972   if (parameterList_ == Teuchos::null) {
00973     RCP<Teuchos::ParameterList> emptyParameterList = Teuchos::rcp(new Teuchos::ParameterList);
00974     this->setParameterList(emptyParameterList);
00975   }
00976 
00977   // Initialize StepControl
00978   if (stepControl_ == Teuchos::null) {
00979     RCP<ImplicitBDFStepperStepControl<Scalar> > implicitBDFStepperStepControl =
00980       Teuchos::rcp(new ImplicitBDFStepperStepControl<Scalar>());
00981     RCP<Teuchos::ParameterList> stepControlPL = 
00982       Teuchos::sublist(parameterList_, RythmosStepControlSettings_name);
00983     implicitBDFStepperStepControl->setParameterList(stepControlPL);
00984     this->setStepControlStrategy(implicitBDFStepperStepControl);
00985   }
00986   stepControl_->initialize(*this);
00987 
00988   maxOrder_ = stepControl_->getMaxOrder(); // maximum order
00989   TEST_FOR_EXCEPTION(
00990       !((1 <= maxOrder_) && (maxOrder_ <= 5)), std::logic_error,
00991       "Error, maxOrder returned from stepControl_->getMaxOrder() = " << maxOrder_ << " is outside range of [1,5]!\n"
00992       );
00993 
00994   Scalar zero = ST::zero();
00995 
00996   currentOrder_ = 1; // Current order of integration
00997   usedOrder_ = 1;  // order used in current step (used after currentOrder_ is updated)
00998   usedStep_ = zero;
00999   nscsco_  =  0;
01000   LETvalue_ = zero;
01001 
01002   alpha_.clear();  // $\alpha_j(n)=h_n/\psi_j(n)$ coefficient used in local error test
01003   // note:   $h_n$ = current step size, n = current time step
01004   gamma_.clear();  // calculate time derivative of history array for predictor 
01005   beta_.clear();   // coefficients used to evaluate predictor from history array
01006   psi_.clear();    // $\psi_j(n) = t_n-t_{n-j}$ intermediary variable used to 
01007   // compute $\beta_j(n):$
01008   for (int i=0 ; i<=maxOrder_ ; ++i) {
01009     alpha_.push_back(zero);
01010     beta_.push_back(zero);
01011     gamma_.push_back(zero);
01012     psi_.push_back(zero);
01013   }
01014   alpha_s_=Scalar(-ST::one());  // $\alpha_s$ fixed-leading coefficient of this BDF method
01015   hh_=zero;
01016   numberOfSteps_=0;   // number of total time integration steps taken
01017   nef_ = 0;
01018 
01019   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
01020     *out << "alpha_s_ = " << alpha_s_ << std::endl;
01021     for (int i=0 ; i<=maxOrder_ ; ++i) {
01022       *out << "alpha_[" << i << "] = " << alpha_[i] << std::endl;
01023       *out << "beta_[" << i << "] = " << beta_[i] << std::endl;
01024       *out << "gamma_[" << i << "] = " << gamma_[i] << std::endl;
01025       *out << "psi_[" << i << "] = " << psi_[i] << std::endl;
01026     }
01027     *out << "numberOfSteps_ = " << numberOfSteps_ << std::endl;
01028   }
01029 
01030   // Generate vectors for use in calculations
01031   RCP<const Thyra::VectorSpaceBase<Scalar> > 
01032     x_space = model_->get_x_space();
01033   RCP<const Thyra::VectorSpaceBase<Scalar> > 
01034     f_space = model_->get_f_space();
01035   x_dot_base_ = createMember(x_space);
01036   ee_ = createMember(x_space);
01037   // x history
01038   xHistory_.clear();
01039   xHistory_.push_back(xn0_->clone_v());
01040   xHistory_.push_back(xpn0_->clone_v());
01041   // Store maxOrder_+1 vectors
01042   for (int i=2 ; i<=maxOrder_ ; ++i) {
01043     xHistory_.push_back(createMember(x_space)); 
01044     V_S(&*xHistory_[i],zero);
01045   }
01046 
01047   isInitialized_ = true;
01048 
01049   if (doTrace) {
01050     *out
01051       << "\nLeaving " << this->Teuchos::Describable::description()
01052       << "::initialize_()...\n";
01053   }
01054 
01055 }
01056 
01057 
01058 template<class Scalar>
01059 void ImplicitBDFStepper<Scalar>::completeStep_()
01060 {
01061 
01062   using Teuchos::as;
01063   typedef Teuchos::ScalarTraits<Scalar> ST;
01064 
01065 #ifdef TEUCHOS_DEBUG
01066   TEST_FOR_EXCEPT(ST::isnaninf(hh_));
01067 #endif  
01068 
01069   numberOfSteps_ ++;
01070   nef_ = 0;
01071   usedStep_ = hh_;
01072   usedOrder_ = currentOrder_;
01073   time_ += hh_;
01074   RCP<Teuchos::FancyOStream> out = this->getOStream();
01075   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
01076   Teuchos::OSTab ostab(out,1,"completeStep_");
01077 
01078   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
01079     *out << "numberOfSteps_ = " << numberOfSteps_ << std::endl;
01080     *out << "time_ = " << time_ << std::endl;
01081   }
01082   
01083   // 03/22/04 tscoffe:  Note that updating the history has nothing to do with
01084   // the step-size and everything to do with the newton correction vectors.
01085   updateHistory_();
01086 
01087 }
01088 
01089 template<class Scalar>
01090 void ImplicitBDFStepper<Scalar>::setStepControlData(const StepperBase<Scalar> & stepper)
01091 {
01092   if (!isInitialized_) {
01093     initialize_();
01094   }
01095   stepControl_->setStepControlData(stepper);
01096 }
01097 
01098 
01099 
01100 } // namespace Rythmos
01101 
01102 
01103 #endif //Rythmos_IMPLICITBDF_STEPPER_H

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