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

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