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

Generated on Wed May 12 21:25:43 2010 for Rythmos - Transient Integration for Differential Equations by  doxygen 1.4.7