Rythmos - Transient Integration for Differential Equations Version of the Day
Rythmos_ForwardEulerStepper_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., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
00023 // USA
00024 // Questions? Contact Todd S. Coffey (tscoffe@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 //@HEADER
00028 
00029 #ifndef Rythmos_FORWARDEULER_STEPPER_DEF_H
00030 #define Rythmos_FORWARDEULER_STEPPER_DEF_H
00031 
00032 #include "Rythmos_ForwardEulerStepper_decl.hpp"
00033 #include "Rythmos_StepperHelpers.hpp"
00034 #include "Thyra_ModelEvaluatorHelpers.hpp"
00035 #include "Thyra_MultiVectorStdOps.hpp"
00036 #include "Thyra_VectorStdOps.hpp"
00037 #include "Teuchos_VerboseObjectParameterListHelpers.hpp"
00038 
00039 
00040 namespace Rythmos {
00041 
00042 
00043 // --------------------------------------------------------------------
00044 // ForwardEulerStepperMomento definitions:
00045 // --------------------------------------------------------------------
00046 
00047 template<class Scalar>
00048 ForwardEulerStepperMomento<Scalar>::ForwardEulerStepperMomento() 
00049 {}
00050 
00051 template<class Scalar>
00052 ForwardEulerStepperMomento<Scalar>::~ForwardEulerStepperMomento() 
00053 {}
00054 
00055 template<class Scalar>
00056 void ForwardEulerStepperMomento<Scalar>::serialize(
00057         const StateSerializerStrategy<Scalar>& stateSerializer,
00058         std::ostream& oStream
00059         ) const
00060 {
00061   using Teuchos::is_null;
00062   TEUCHOS_ASSERT( !is_null(model_) );
00063   RCP<VectorBase<Scalar> > sol_vec = solution_vector_;
00064   if (is_null(sol_vec)) {
00065     sol_vec = Thyra::createMember(model_->get_x_space());
00066   }
00067   RCP<VectorBase<Scalar> > res_vec = residual_vector_;
00068   if (is_null(res_vec)) {
00069     res_vec = Thyra::createMember(model_->get_f_space());
00070   }
00071   RCP<VectorBase<Scalar> > sol_vec_old = solution_vector_old_;
00072   if (is_null(sol_vec_old)) {
00073     sol_vec_old = Thyra::createMember(model_->get_x_space());
00074   }
00075   stateSerializer.serializeVectorBase(*sol_vec,oStream);
00076   stateSerializer.serializeVectorBase(*res_vec,oStream);
00077   stateSerializer.serializeVectorBase(*sol_vec_old,oStream);
00078   stateSerializer.serializeScalar(t_,oStream);
00079   stateSerializer.serializeScalar(t_old_,oStream);
00080   stateSerializer.serializeScalar(dt_,oStream);
00081   stateSerializer.serializeInt(numSteps_,oStream);
00082   stateSerializer.serializeBool(isInitialized_,oStream);
00083   stateSerializer.serializeBool(haveInitialCondition_,oStream);
00084   RCP<ParameterList> pl = parameterList_;
00085   if (Teuchos::is_null(pl)) {
00086     pl = Teuchos::parameterList();
00087   }
00088   stateSerializer.serializeParameterList(*pl,oStream);
00089 }
00090 
00091 template<class Scalar>
00092 void ForwardEulerStepperMomento<Scalar>::deSerialize(
00093         const StateSerializerStrategy<Scalar>& stateSerializer,
00094         std::istream& iStream
00095         )
00096 {
00097   using Teuchos::outArg;
00098   using Teuchos::is_null;
00099   TEUCHOS_ASSERT( !is_null(model_) );
00100   if (is_null(solution_vector_)) {
00101     solution_vector_ = Thyra::createMember(*model_->get_x_space());
00102   }
00103   if (is_null(residual_vector_)) {
00104     residual_vector_ = Thyra::createMember(*model_->get_f_space());
00105   }
00106   if (is_null(solution_vector_old_)) {
00107     solution_vector_old_ = Thyra::createMember(*model_->get_x_space());
00108   }
00109   stateSerializer.deSerializeVectorBase(outArg(*solution_vector_),iStream);
00110   stateSerializer.deSerializeVectorBase(outArg(*residual_vector_),iStream);
00111   stateSerializer.deSerializeVectorBase(outArg(*solution_vector_old_),iStream);
00112   stateSerializer.deSerializeScalar(outArg(t_),iStream);
00113   stateSerializer.deSerializeScalar(outArg(t_old_),iStream);
00114   stateSerializer.deSerializeScalar(outArg(dt_),iStream);
00115   stateSerializer.deSerializeInt(outArg(numSteps_),iStream);
00116   stateSerializer.deSerializeBool(outArg(isInitialized_),iStream);
00117   stateSerializer.deSerializeBool(outArg(haveInitialCondition_),iStream);
00118   if (is_null(parameterList_)) {
00119     parameterList_ = Teuchos::parameterList();
00120   }
00121   stateSerializer.deSerializeParameterList(outArg(*parameterList_),iStream);
00122 }
00123 
00124 template<class Scalar>
00125 RCP<MomentoBase<Scalar> > ForwardEulerStepperMomento<Scalar>::clone() const
00126 {
00127   RCP<ForwardEulerStepperMomento<Scalar> > m = rcp(new ForwardEulerStepperMomento<Scalar>());
00128   m->set_solution_vector(solution_vector_);
00129   m->set_residual_vector(residual_vector_);
00130   m->set_solution_vector_old(solution_vector_old_);
00131   m->set_t(t_);
00132   m->set_t_old(t_old_);
00133   m->set_dt(dt_);
00134   m->set_numSteps(numSteps_);
00135   m->set_isInitialized(isInitialized_);
00136   m->set_haveInitialCondition(haveInitialCondition_);
00137   m->set_parameterList(parameterList_);
00138   if (!Teuchos::is_null(this->getMyParamList())) {
00139     m->setParameterList(Teuchos::parameterList(*(this->getMyParamList())));
00140   }
00141   m->set_model(model_);
00142   m->set_basePoint(basePoint_);
00143   // How do I copy the VerboseObject data?  
00144   // 07/10/09 tscoffe:  Its not set up in Teuchos to do this yet
00145   return m;
00146 }
00147 
00148 template<class Scalar>
00149 void ForwardEulerStepperMomento<Scalar>::set_solution_vector(const RCP<const VectorBase<Scalar> >& solution_vector )
00150 { 
00151   solution_vector_ = Teuchos::null;
00152   if (!Teuchos::is_null(solution_vector)) {
00153     solution_vector_ = solution_vector->clone_v(); 
00154   }
00155 }
00156 
00157 template<class Scalar>
00158 RCP<VectorBase<Scalar> > ForwardEulerStepperMomento<Scalar>::get_solution_vector() const
00159 { 
00160   return solution_vector_; 
00161 }
00162 
00163 template<class Scalar>
00164 void ForwardEulerStepperMomento<Scalar>::set_residual_vector(const RCP<const VectorBase<Scalar> >& residual_vector)
00165 { 
00166   residual_vector_ = Teuchos::null;
00167   if (!Teuchos::is_null(residual_vector)) {
00168     residual_vector_ = residual_vector->clone_v(); 
00169   }
00170 }
00171 
00172 template<class Scalar>
00173 RCP<VectorBase<Scalar> > ForwardEulerStepperMomento<Scalar>::get_residual_vector() const
00174 { 
00175   return residual_vector_; 
00176 }
00177 
00178 template<class Scalar>
00179 void ForwardEulerStepperMomento<Scalar>::set_solution_vector_old(const RCP<const VectorBase<Scalar> >& solution_vector_old )
00180 { 
00181   solution_vector_old_ = Teuchos::null;
00182   if (!Teuchos::is_null(solution_vector_old)) {
00183     solution_vector_old_ = solution_vector_old->clone_v(); 
00184   }
00185 }
00186 
00187 template<class Scalar>
00188 RCP<VectorBase<Scalar> > ForwardEulerStepperMomento<Scalar>::get_solution_vector_old() const
00189 { 
00190   return solution_vector_old_; 
00191 }
00192 
00193 template<class Scalar>
00194 void ForwardEulerStepperMomento<Scalar>::set_t(const Scalar & t)
00195 { 
00196   t_ = t; 
00197 }
00198 
00199 template<class Scalar>
00200 Scalar ForwardEulerStepperMomento<Scalar>::get_t() const
00201 { 
00202   return t_; 
00203 }
00204 
00205 template<class Scalar>
00206 void ForwardEulerStepperMomento<Scalar>::set_t_old(const Scalar & t_old)
00207 { 
00208   t_old_ = t_old; 
00209 }
00210 
00211 template<class Scalar>
00212 Scalar ForwardEulerStepperMomento<Scalar>::get_t_old() const
00213 { 
00214   return t_old_; 
00215 }
00216 
00217 template<class Scalar>
00218 void ForwardEulerStepperMomento<Scalar>::set_dt(const Scalar & dt)
00219 { 
00220   dt_ = dt; 
00221 }
00222 
00223 template<class Scalar>
00224 Scalar ForwardEulerStepperMomento<Scalar>::get_dt() const
00225 { 
00226   return dt_; 
00227 }
00228 
00229 template<class Scalar>
00230 void ForwardEulerStepperMomento<Scalar>::set_numSteps(const int & numSteps)
00231 { 
00232   numSteps_ = numSteps; 
00233 }
00234 
00235 template<class Scalar>
00236 int ForwardEulerStepperMomento<Scalar>::get_numSteps() const
00237 { 
00238   return numSteps_; 
00239 }
00240 
00241 template<class Scalar>
00242 void ForwardEulerStepperMomento<Scalar>::set_isInitialized(const bool & isInitialized)
00243 { 
00244   isInitialized_ = isInitialized; 
00245 }
00246 
00247 template<class Scalar>
00248 bool ForwardEulerStepperMomento<Scalar>::get_isInitialized() const
00249 { 
00250   return isInitialized_; 
00251 }
00252 
00253 template<class Scalar>
00254 void ForwardEulerStepperMomento<Scalar>::set_haveInitialCondition(const bool & haveInitialCondition)
00255 { 
00256   haveInitialCondition_ = haveInitialCondition; 
00257 }
00258 
00259 template<class Scalar>
00260 bool ForwardEulerStepperMomento<Scalar>::get_haveInitialCondition() const
00261 { 
00262   return haveInitialCondition_; 
00263 }
00264 
00265 template<class Scalar>
00266 void ForwardEulerStepperMomento<Scalar>::set_parameterList(const RCP<const ParameterList>& pl)
00267 { 
00268   parameterList_ = Teuchos::null;
00269   if (!Teuchos::is_null(pl)) {
00270     parameterList_ = Teuchos::parameterList(*pl); 
00271   }
00272 }
00273 
00274 template<class Scalar>
00275 RCP<ParameterList> ForwardEulerStepperMomento<Scalar>::get_parameterList() const
00276 { 
00277   return parameterList_; 
00278 }
00279 
00280 template<class Scalar>
00281 void ForwardEulerStepperMomento<Scalar>::setParameterList(const RCP<ParameterList>& paramList)
00282 { 
00283   this->setMyParamList(paramList); 
00284 }
00285 
00286 template<class Scalar>
00287 RCP<const ParameterList> ForwardEulerStepperMomento<Scalar>::getValidParameters() const
00288 { 
00289   return Teuchos::null; 
00290 }
00291 
00292 template<class Scalar>
00293 void ForwardEulerStepperMomento<Scalar>::set_model(const RCP<const Thyra::ModelEvaluator<Scalar> >& model)
00294 { 
00295   model_ = model; 
00296 }
00297 
00298 template<class Scalar>
00299 RCP<const Thyra::ModelEvaluator<Scalar> > ForwardEulerStepperMomento<Scalar>::get_model() const
00300 { 
00301   return model_; 
00302 }
00303 
00304 template<class Scalar>
00305 void ForwardEulerStepperMomento<Scalar>::set_basePoint(const RCP<const Thyra::ModelEvaluatorBase::InArgs<Scalar> >& basePoint)
00306 { 
00307   basePoint_ = basePoint; 
00308 }
00309 
00310 template<class Scalar>
00311 RCP<const Thyra::ModelEvaluatorBase::InArgs<Scalar> > ForwardEulerStepperMomento<Scalar>::get_basePoint() const
00312 { 
00313   return basePoint_; 
00314 }
00315 
00316 // --------------------------------------------------------------------
00317 // ForwardEulerStepper definitions:
00318 // --------------------------------------------------------------------
00319 
00320 // Nonmember constructor
00321 template<class Scalar>
00322 RCP<ForwardEulerStepper<Scalar> > forwardEulerStepper()
00323 {
00324   RCP<ForwardEulerStepper<Scalar> > stepper = rcp(new ForwardEulerStepper<Scalar>());
00325   return stepper;
00326 }
00327 
00328 // Nonmember constructor
00329 template<class Scalar>
00330 RCP<ForwardEulerStepper<Scalar> > forwardEulerStepper(
00331     const RCP<Thyra::ModelEvaluator<Scalar> >& model
00332     )
00333 {
00334   RCP<ForwardEulerStepper<Scalar> > stepper = forwardEulerStepper<Scalar>();
00335   {
00336     RCP<StepperBase<Scalar> > stepperBase = 
00337       Teuchos::rcp_dynamic_cast<StepperBase<Scalar> >(stepper,true);
00338     setStepperModel(Teuchos::inOutArg(*stepperBase),model);
00339   }
00340   return stepper;
00341 }
00342 
00343 template<class Scalar>
00344 ForwardEulerStepper<Scalar>::ForwardEulerStepper()
00345 {
00346   this->defaultInitializAll_();
00347   dt_ = Teuchos::ScalarTraits<Scalar>::zero();
00348   numSteps_ = 0;
00349 }
00350 
00351 template<class Scalar>
00352 void ForwardEulerStepper<Scalar>::defaultInitializAll_()
00353 {
00354   model_ = Teuchos::null;
00355   solution_vector_ = Teuchos::null;
00356   residual_vector_ = Teuchos::null;
00357   t_ = ST::nan();
00358   dt_ = ST::nan();
00359   t_old_ = ST::nan();
00360   solution_vector_old_ = Teuchos::null;
00361   //basePoint_;
00362   numSteps_ = -1;
00363   haveInitialCondition_ = false;
00364   parameterList_ = Teuchos::null;
00365   isInitialized_ = false;
00366 }
00367 
00368 template<class Scalar>
00369 void ForwardEulerStepper<Scalar>::initialize_()
00370 {
00371   if (!isInitialized_) {
00372     TEUCHOS_TEST_FOR_EXCEPTION( is_null(model_), std::logic_error,
00373        "Error!  Please set a model on the stepper.\n" 
00374        );
00375     residual_vector_ = Thyra::createMember(model_->get_f_space());
00376     isInitialized_ = true;
00377   }
00378 }
00379 
00380 template<class Scalar>
00381 ForwardEulerStepper<Scalar>::~ForwardEulerStepper()
00382 {
00383 }
00384 
00385 template<class Scalar>
00386 RCP<const Thyra::VectorSpaceBase<Scalar> > ForwardEulerStepper<Scalar>::get_x_space() const
00387 {
00388   TEUCHOS_TEST_FOR_EXCEPTION(!haveInitialCondition_,std::logic_error,"Error, attempting to call get_x_space before setting an initial condition!\n");
00389   return(solution_vector_->space());
00390 }
00391 
00392 template<class Scalar>
00393 Scalar ForwardEulerStepper<Scalar>::takeStep(Scalar dt, StepSizeType flag)
00394 {
00395   TEUCHOS_TEST_FOR_EXCEPTION( !haveInitialCondition_, std::logic_error,
00396      "Error!  Attempting to call takeStep before setting an initial condition!\n" 
00397      );
00398   this->initialize_();
00399   if (flag == STEP_TYPE_VARIABLE) { 
00400     // print something out about this method not supporting automatic variable step-size
00401     return(-ST::one());
00402   }
00403   //Thyra::eval_f<Scalar>(*model_,*solution_vector_,t_+dt,&*residual_vector_);
00404   eval_model_explicit<Scalar>(*model_,basePoint_,*solution_vector_,t_+dt,Teuchos::outArg(*residual_vector_));
00405 
00406   // solution_vector_old_ = solution_vector_
00407   Thyra::V_V(Teuchos::outArg(*solution_vector_old_),*solution_vector_);
00408   // solution_vector = solution_vector + dt*residual_vector
00409   Thyra::Vp_StV(solution_vector_.ptr(),dt,*residual_vector_); 
00410   t_old_ = t_;
00411   t_ += dt;
00412   dt_ = dt;
00413   numSteps_++;
00414 
00415   return(dt);
00416 }
00417 
00418 template<class Scalar>
00419 const StepStatus<Scalar> ForwardEulerStepper<Scalar>::getStepStatus() const
00420 {
00421   StepStatus<Scalar> stepStatus;
00422 
00423   if (!haveInitialCondition_) {
00424     stepStatus.stepStatus = STEP_STATUS_UNINITIALIZED;
00425   } 
00426   else if (numSteps_ == 0) {
00427     stepStatus.stepStatus = STEP_STATUS_UNKNOWN;
00428     stepStatus.order = this->getOrder();
00429     stepStatus.time = t_;
00430     stepStatus.solution = solution_vector_;
00431   } 
00432   else {
00433     stepStatus.stepStatus = STEP_STATUS_CONVERGED;
00434     stepStatus.stepSize = dt_; 
00435     stepStatus.order = this->getOrder();
00436     stepStatus.time = t_;
00437     stepStatus.stepLETValue = Scalar(-ST::one()); 
00438     stepStatus.solution = solution_vector_;
00439     stepStatus.residual = residual_vector_;
00440   }
00441 
00442   return(stepStatus);
00443 }
00444 
00445 template<class Scalar>
00446 std::string ForwardEulerStepper<Scalar>::description() const
00447 {
00448   std::string name = "Rythmos::ForwardEulerStepper";
00449   return(name);
00450 }
00451 
00452 template<class Scalar>
00453 void ForwardEulerStepper<Scalar>::describe(
00454       Teuchos::FancyOStream            &out,
00455       const Teuchos::EVerbosityLevel   verbLevel
00456       ) const
00457 {
00458   if ( (static_cast<int>(verbLevel) == static_cast<int>(Teuchos::VERB_DEFAULT) ) ||
00459        (static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW)     )
00460      ) {
00461     out << description() << "::describe" << std::endl;
00462     out << "model = " << model_->description() << std::endl;
00463   } else if (static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW)) {
00464   } else if (static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_MEDIUM)) {
00465   } else if (static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_HIGH)) {
00466     out << "model = " << std::endl;
00467     model_->describe(out,verbLevel);
00468     out << "solution_vector = " << std::endl;
00469     solution_vector_->describe(out,verbLevel);
00470     out << "residual_vector = " << std::endl;
00471     residual_vector_->describe(out,verbLevel);
00472   }
00473 }
00474 
00475 template<class Scalar>
00476 void ForwardEulerStepper<Scalar>::addPoints(
00477     const Array<Scalar>& time_vec
00478     ,const Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >& x_vec
00479     ,const Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >& xdot_vec
00480     )
00481 {
00482   TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Error, addPoints is not implemented for ForwardEulerStepper.\n");
00483 }
00484 
00485 template<class Scalar>
00486 void ForwardEulerStepper<Scalar>::getPoints(
00487     const Array<Scalar>& time_vec
00488     ,Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >* x_vec
00489     ,Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec
00490     ,Array<ScalarMag>* accuracy_vec) const
00491 {
00492   TEUCHOS_ASSERT( haveInitialCondition_ );
00493   using Teuchos::constOptInArg;
00494   using Teuchos::null;
00495   defaultGetPoints<Scalar>(
00496       t_old_,
00497       constOptInArg(*solution_vector_old_),
00498       Ptr<const VectorBase<Scalar> >(null),
00499       t_,
00500       constOptInArg(*solution_vector_),
00501       Ptr<const VectorBase<Scalar> >(null),
00502       time_vec,
00503       ptr(x_vec),
00504       ptr(xdot_vec),
00505       ptr(accuracy_vec),
00506       Ptr<InterpolatorBase<Scalar> >(null)
00507       );
00508 }
00509 
00510 template<class Scalar>
00511 TimeRange<Scalar> ForwardEulerStepper<Scalar>::getTimeRange() const
00512 {
00513   if (!haveInitialCondition_) {
00514     return(invalidTimeRange<Scalar>());
00515   } else {
00516     return(TimeRange<Scalar>(t_old_,t_));
00517   }
00518 }
00519 
00520 template<class Scalar>
00521 void ForwardEulerStepper<Scalar>::getNodes(Array<Scalar>* time_vec) const
00522 {
00523   TEUCHOS_ASSERT( time_vec != NULL );
00524   time_vec->clear();
00525   if (!haveInitialCondition_) {
00526     return; 
00527   } else {
00528     time_vec->push_back(t_old_);
00529   }
00530   if (numSteps_ > 0) {
00531     time_vec->push_back(t_);
00532   }
00533 }
00534 
00535 template<class Scalar>
00536 void ForwardEulerStepper<Scalar>::removeNodes(Array<Scalar>& time_vec) 
00537 {
00538   TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Error, removeNodes is not implemented for ForwardEulerStepper.\n");
00539 }
00540 
00541 template<class Scalar>
00542 int ForwardEulerStepper<Scalar>::getOrder() const
00543 {
00544   return(1);
00545 }
00546 
00547 template <class Scalar>
00548 void ForwardEulerStepper<Scalar>::setParameterList(Teuchos::RCP<Teuchos::ParameterList> const& paramList)
00549 {
00550   TEUCHOS_TEST_FOR_EXCEPT(is_null(paramList));
00551   paramList->validateParametersAndSetDefaults(*this->getValidParameters());
00552   parameterList_ = paramList;
00553   Teuchos::readVerboseObjectSublist(&*parameterList_,this);
00554 }
00555 
00556 template <class Scalar>
00557 Teuchos::RCP<Teuchos::ParameterList> ForwardEulerStepper<Scalar>::getNonconstParameterList()
00558 {
00559   return(parameterList_);
00560 }
00561 
00562 template <class Scalar>
00563 Teuchos::RCP<Teuchos::ParameterList> ForwardEulerStepper<Scalar>::unsetParameterList()
00564 {
00565   Teuchos::RCP<Teuchos::ParameterList> temp_param_list = parameterList_;
00566   parameterList_ = Teuchos::null;
00567   return(temp_param_list);
00568 }
00569 
00570 
00571 template<class Scalar>
00572 RCP<const Teuchos::ParameterList>
00573 ForwardEulerStepper<Scalar>::getValidParameters() const
00574 {
00575   using Teuchos::ParameterList;
00576   static RCP<const ParameterList> validPL;
00577   if (is_null(validPL)) {
00578     RCP<ParameterList> pl = Teuchos::parameterList();
00579     Teuchos::setupVerboseObjectSublist(&*pl);
00580     validPL = pl;
00581   }
00582   return validPL;
00583 }
00584 
00585 
00586 template<class Scalar>
00587 void ForwardEulerStepper<Scalar>::setModel(const RCP<const Thyra::ModelEvaluator<Scalar> >& model)
00588 {
00589   TEUCHOS_TEST_FOR_EXCEPT( is_null(model) );
00590   assertValidModel( *this, *model );
00591   model_ = model;
00592 }
00593 
00594 
00595 template<class Scalar>
00596 void ForwardEulerStepper<Scalar>::setNonconstModel(const RCP<Thyra::ModelEvaluator<Scalar> >& model)
00597 {
00598   this->setModel(model); // TODO 09/09/09 tscoffe:  use ConstNonconstObjectContainer!
00599 }
00600 
00601 
00602 template<class Scalar>
00603 Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >
00604 ForwardEulerStepper<Scalar>::getModel() const
00605 {
00606   return model_;
00607 }
00608 
00609 template<class Scalar>
00610 Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >
00611 ForwardEulerStepper<Scalar>::getNonconstModel() 
00612 {
00613   return Teuchos::null;
00614 }
00615 
00616 template<class Scalar>
00617 void ForwardEulerStepper<Scalar>::setInitialCondition(
00618     const Thyra::ModelEvaluatorBase::InArgs<Scalar> &initialCondition
00619     )
00620 {
00621   basePoint_ = initialCondition;
00622   t_ = initialCondition.get_t();
00623   t_old_ = t_;
00624   solution_vector_ = initialCondition.get_x()->clone_v();
00625   solution_vector_old_ = solution_vector_->clone_v();
00626   haveInitialCondition_ = true;
00627 }
00628 
00629 
00630 template<class Scalar>
00631 Thyra::ModelEvaluatorBase::InArgs<Scalar> 
00632 ForwardEulerStepper<Scalar>::getInitialCondition() const
00633 {
00634   return basePoint_;
00635 }
00636 
00637 
00638 template<class Scalar>
00639 bool ForwardEulerStepper<Scalar>::supportsCloning() const
00640 {
00641   return true;
00642 }
00643 
00644 template<class Scalar>
00645 RCP<StepperBase<Scalar> >
00646 ForwardEulerStepper<Scalar>::cloneStepperAlgorithm() const
00647 {
00648 
00649   // Just use the interface to clone the algorithm in a basically
00650   // uninitialized state
00651 
00652   RCP<StepperBase<Scalar> >
00653     stepper = Teuchos::rcp(new ForwardEulerStepper<Scalar>());
00654 
00655   if (!is_null(model_)) {
00656     setStepperModel(Teuchos::inOutArg(*stepper),model_); // Shallow copy is okay!
00657   }
00658 
00659   if (!is_null(parameterList_)) {
00660     stepper->setParameterList(Teuchos::parameterList(*parameterList_));
00661   }
00662 
00663   return stepper;
00664 
00665 }
00666 
00667 template<class Scalar>
00668 RCP<const MomentoBase<Scalar> >
00669 ForwardEulerStepper<Scalar>::getMomento() const
00670 {
00671   RCP<ForwardEulerStepperMomento<Scalar> > momento = Teuchos::rcp(new ForwardEulerStepperMomento<Scalar>());
00672   momento->set_solution_vector(solution_vector_);
00673   momento->set_solution_vector_old(solution_vector_old_);
00674   momento->set_residual_vector(residual_vector_);
00675   momento->set_t(t_);
00676   momento->set_t_old(t_old_);
00677   momento->set_dt(dt_);
00678   momento->set_numSteps(numSteps_);
00679   momento->set_isInitialized(isInitialized_);
00680   momento->set_haveInitialCondition(haveInitialCondition_);
00681   momento->set_parameterList(parameterList_);
00682   momento->set_model(model_);
00683   RCP<Thyra::ModelEvaluatorBase::InArgs<Scalar> > bp = Teuchos::rcp(new Thyra::ModelEvaluatorBase::InArgs<Scalar>(basePoint_));
00684   momento->set_basePoint(bp);
00685   return momento;
00686 }
00687 
00688 template<class Scalar>
00689 void ForwardEulerStepper<Scalar>::setMomento( const Ptr<const MomentoBase<Scalar> >& momentoPtr ) 
00690 { 
00691   Ptr<const ForwardEulerStepperMomento<Scalar> > feMomentoPtr = 
00692     Teuchos::ptr_dynamic_cast<const ForwardEulerStepperMomento<Scalar> >(momentoPtr,true);
00693   TEUCHOS_TEST_FOR_EXCEPTION( Teuchos::is_null(feMomentoPtr->get_model()), std::logic_error,
00694       "Error!  Rythmos::ForwardEulerStepper::setMomento:  The momento must have a valid model through set_model(...) prior to calling ForwardEulerStepper::setMomento(...)."
00695       );
00696   TEUCHOS_TEST_FOR_EXCEPTION( Teuchos::is_null(feMomentoPtr->get_basePoint()), std::logic_error,
00697       "Error!  Rythmos::ForwardEulerStepper::setMomento:  The momento must have a valid base point through set_basePoint(...) prior to calling ForwardEulerStepper::setMomento(...)."
00698       );
00699   model_ = feMomentoPtr->get_model();
00700   basePoint_ = *(feMomentoPtr->get_basePoint());
00701   const ForwardEulerStepperMomento<Scalar>& feMomento = *feMomentoPtr;
00702   solution_vector_ = feMomento.get_solution_vector();
00703   solution_vector_old_ = feMomento.get_solution_vector_old();
00704   residual_vector_ = feMomento.get_residual_vector();
00705   t_ = feMomento.get_t();
00706   t_old_ = feMomento.get_t_old();
00707   dt_ = feMomento.get_dt();
00708   numSteps_ = feMomento.get_numSteps();
00709   isInitialized_ = feMomento.get_isInitialized();
00710   haveInitialCondition_ = feMomento.get_haveInitialCondition();
00711   parameterList_ = feMomento.get_parameterList();
00712   this->checkConsistentState_();
00713 }
00714 
00715 template<class Scalar>
00716 void ForwardEulerStepper<Scalar>::checkConsistentState_()
00717 {
00718   if (isInitialized_) {
00719     TEUCHOS_ASSERT( !Teuchos::is_null(model_) );
00720     TEUCHOS_ASSERT( !Teuchos::is_null(residual_vector_) );
00721   }
00722   if (haveInitialCondition_) {
00723     TEUCHOS_ASSERT( !ST::isnaninf(t_) );
00724     TEUCHOS_ASSERT( !ST::isnaninf(t_old_) );
00725     TEUCHOS_ASSERT( !Teuchos::is_null(solution_vector_) );
00726     TEUCHOS_ASSERT( !Teuchos::is_null(solution_vector_old_) );
00727     TEUCHOS_ASSERT( t_ >= basePoint_.get_t() );
00728     TEUCHOS_ASSERT( t_old_ >= basePoint_.get_t() );
00729   }
00730   if (numSteps_ > 0) {
00731     TEUCHOS_ASSERT(isInitialized_);
00732     TEUCHOS_ASSERT(haveInitialCondition_);
00733   }
00734 }
00735 
00736 // 
00737 // Explicit Instantiation macro
00738 //
00739 // Must be expanded from within the Rythmos namespace!
00740 //
00741 
00742 #define RYTHMOS_FORWARD_EULER_STEPPER_INSTANT(SCALAR) \
00743   \
00744   template class ForwardEulerStepperMomento< SCALAR >; \
00745   template class ForwardEulerStepper< SCALAR >; \
00746   \
00747   template RCP<ForwardEulerStepper< SCALAR > > forwardEulerStepper(); \
00748   template RCP<ForwardEulerStepper< SCALAR > > forwardEulerStepper( \
00749       const RCP<Thyra::ModelEvaluator< SCALAR > >& model \
00750       ); 
00751 
00752 
00753 } // namespace Rythmos
00754 
00755 #endif //Rythmos_FORWARDEULER_STEPPER_DEF_H
 All Classes Functions Variables Typedefs Friends