Rythmos_StepperValidator.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_STEPPER_VALIDATOR_H
00030 #define Rythmos_STEPPER_VALIDATOR_H
00031 
00032 #include "Teuchos_ParameterList.hpp"
00033 #include "Teuchos_ParameterListAcceptor.hpp"
00034 #include "Teuchos_Describable.hpp"
00035 #include "Teuchos_VerboseObject.hpp"
00036 #include "Rythmos_IntegratorBuilder.hpp"
00037 #include "Thyra_StateFuncModelEvaluatorBase.hpp" 
00038 #include "Teuchos_ParameterListAcceptorDefaultBase.hpp"
00039 #include "Thyra_ModelEvaluator.hpp" 
00040 
00041 #include "Rythmos_StepperBase.hpp"
00042 #include "Rythmos_Types.hpp"
00043 #include "Teuchos_VerboseObjectParameterListHelpers.hpp"
00044 #include "Thyra_DetachedVectorView.hpp"
00045 #include "Thyra_DefaultSpmdVectorSpace.hpp"
00046 #include "Thyra_DefaultSerialDenseLinearOpWithSolveFactory.hpp"
00047 #include "Rythmos_TimeStepNonlinearSolver.hpp"
00048 
00049 #include "Teuchos_StandardCatchMacros.hpp"
00050 
00051 
00052 namespace Rythmos {
00053 
00054 
00062 template<class Scalar> 
00063 class StepperValidator 
00064   : virtual public Teuchos::Describable
00065   , virtual public Teuchos::ParameterListAcceptor
00066   , virtual public Teuchos::VerboseObject<StepperValidator<Scalar> >
00067 {
00068 public:
00069 
00070   StepperValidator();
00071 
00072   virtual ~StepperValidator();
00073 
00075   void setIntegratorBuilder(
00076     const RCP<IntegratorBuilder<Scalar> > &integratorBuilder
00077     );
00078 
00080   void validateStepper() const;
00081 
00084 
00086   void setParameterList(RCP<Teuchos::ParameterList> const& paramList);
00087 
00089   RCP<Teuchos::ParameterList> getNonconstParameterList();
00090 
00092   RCP<Teuchos::ParameterList> unsetParameterList();
00093 
00095   RCP<const Teuchos::ParameterList> getValidParameters() const;
00096 
00098 
00101 
00103   std::string description() const;
00104 
00106   void describe(
00107     Teuchos::FancyOStream &out,
00108     const Teuchos::EVerbosityLevel verbLevel
00109     ) const;
00110 
00112 
00113 private:
00114 
00115   // Private member functions:
00116   void defaultInitializeAll_();
00117   RCP<StepperBase<Scalar> > getStepper_(
00118       const RCP<const Thyra::ModelEvaluator<Scalar> > model,
00119       const Thyra::ModelEvaluatorBase::InArgs<Scalar> &initialCondition,
00120       const RCP<Thyra::NonlinearSolverBase<Scalar> > nlSolver
00121       ) const;
00122   bool isImplicitStepper_() const;
00123   Thyra::ModelEvaluatorBase::InArgs<Scalar> getSomeIC_(
00124     const Thyra::ModelEvaluator<Scalar>& model) const;
00125 
00126   // Validate that parameters set through setInitialCondition actually get set
00127   // on the model when evalModel is called.
00128   void validateIC_() const;
00129   // Validate that getTimeRange and getStepStatus return the correct states of
00130   // initialization.
00131   void validateStates_() const;
00132   // Validate that we can get the initial condition through getPoints after
00133   // setInitialCondition has been set and after the first step.
00134   void validateGetIC_() const;
00135   // Validate that the stepper supports getNodes, which is used by the
00136   // Trailing Interpolation Buffer feature of the Integrator.
00137   void validateGetNodes_() const;
00138 
00139   // Private member data:
00140   RCP<IntegratorBuilder<Scalar> > integratorBuilder_;
00141   RCP<ParameterList> paramList_;
00142   mutable std::string stepperName_;
00143 
00144 };
00145 
00146 //
00147 // StepperValidator nonmember constructor:
00148 //
00149 template<class Scalar>
00150 RCP<StepperValidator<Scalar> > stepperValidator()
00151 {
00152   return rcp(new StepperValidator<Scalar>() );
00153 }
00154 
00155 //
00156 //  Mock Model class for validating a stepper
00157 //
00158 template<class Scalar>
00159 class StepperValidatorMockModel
00160   : public Thyra::StateFuncModelEvaluatorBase<Scalar>,
00161     public Teuchos::ParameterListAcceptorDefaultBase
00162 {
00163 public:
00164 
00165   StepperValidatorMockModel();
00166 
00167   virtual ~StepperValidatorMockModel();
00168 
00169   RCP<const std::vector<Thyra::ModelEvaluatorBase::InArgs<Scalar> > > getPassedInArgs() const;
00170 
00173 
00175   RCP<const Thyra::VectorSpaceBase<Scalar> > get_x_space() const;
00177   RCP<const Thyra::VectorSpaceBase<Scalar> > get_f_space() const;
00179   Thyra::ModelEvaluatorBase::InArgs<Scalar> getNominalValues() const;
00181   RCP<Thyra::LinearOpWithSolveBase<Scalar> > create_W() const;
00183   RCP<Thyra::LinearOpBase<Scalar> > create_W_op() const;
00185   RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar> > get_W_factory() const;
00187   Thyra::ModelEvaluatorBase::InArgs<Scalar> createInArgs() const;
00188 
00190   RCP<const Thyra::VectorSpaceBase<Scalar> > get_p_space(int l) const;
00192   RCP<const Teuchos::Array<std::string> > get_p_names(int l) const;
00194   RCP<const Thyra::VectorSpaceBase<Scalar> > get_g_space(int j) const;
00195 
00197   
00200   
00202   void setParameterList(RCP<ParameterList> const& paramList);
00203 
00205   RCP<const ParameterList> getValidParameters() const;
00206 
00208   
00209 private:
00210 
00213 
00215   Thyra::ModelEvaluatorBase::OutArgs<Scalar> createOutArgsImpl() const;
00216 
00218   void evalModelImpl(
00219     const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs_bar,
00220     const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs_bar
00221     ) const;
00222 
00224   
00225   void defaultInitializeAll_();
00226   void initialize_();
00227 
00228   mutable RCP<std::vector<Thyra::ModelEvaluatorBase::InArgs<Scalar> > > passedInArgs_;
00229   bool isImplicit_;
00230   bool isInitialized_;
00231   Thyra::ModelEvaluatorBase::InArgs<Scalar> inArgs_;
00232   Thyra::ModelEvaluatorBase::OutArgs<Scalar> outArgs_;
00233   Thyra::ModelEvaluatorBase::InArgs<Scalar> nominalValues_;
00234   RCP<const Thyra::VectorSpaceBase<Scalar> > x_space_;
00235   RCP<const Thyra::VectorSpaceBase<Scalar> > f_space_;
00236   RCP<const Thyra::VectorSpaceBase<Scalar> > p_space_;
00237   int Np_; // number of parameter vectors (1)
00238   int Ng_; // number of observation functions (0)
00239   int np_; // length of parameter vector (1)
00240   int dim_; // length of solution vector
00241 
00242 };
00243 
00244 //
00245 // StepperValidatorMockModel nonmember constructors:
00246 //
00247 template<class Scalar>
00248 RCP<StepperValidatorMockModel<Scalar> > stepperValidatorMockModel() 
00249 {
00250   return(rcp(new StepperValidatorMockModel<Scalar>()));
00251 }
00252 
00253 template<class Scalar>
00254 RCP<StepperValidatorMockModel<Scalar> > stepperValidatorMockModel(
00255     bool isImplicit
00256     ) 
00257 {
00258   RCP<StepperValidatorMockModel<Scalar> > model = rcp(new StepperValidatorMockModel<Scalar>());
00259   RCP<ParameterList> pl = Teuchos::parameterList();
00260   pl->set("Is Implicit",isImplicit);
00261   model->setParameterList(pl);
00262   return model;
00263 }
00264 
00265 //
00266 // StepperValidatorMockModel Definitions:
00267 //
00268 template<class Scalar>
00269 StepperValidatorMockModel<Scalar>::StepperValidatorMockModel()
00270 { 
00271   this->defaultInitializeAll_();
00272   Np_ = 1;
00273   Ng_ = 0;
00274   np_ = 1;
00275   dim_ = 1;
00276   // Create x_space and f_space
00277   x_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(dim_);
00278   f_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(dim_);
00279   // Create p_space 
00280   p_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(np_);
00281   passedInArgs_ = rcp(new std::vector<Thyra::ModelEvaluatorBase::InArgs<Scalar> >);
00282   this->initialize_();
00283 }
00284 
00285 template<class Scalar>
00286 void StepperValidatorMockModel<Scalar>::defaultInitializeAll_()
00287 {
00288   passedInArgs_ = Teuchos::null;
00289   isImplicit_ = false;
00290   isInitialized_ = false;
00291   //inArgs_;
00292   //outArgs_;
00293   //nominalValues_;
00294   x_space_ = Teuchos::null;
00295   f_space_ = Teuchos::null;
00296   p_space_ = Teuchos::null;
00297   Np_ = -1;
00298   Ng_ = -1;
00299   np_ = -1;
00300   dim_ = -1;
00301 }
00302 
00303 template<class Scalar>
00304 StepperValidatorMockModel<Scalar>::~StepperValidatorMockModel()
00305 { }
00306 
00307 template<class Scalar>
00308 void StepperValidatorMockModel<Scalar>::initialize_()
00309 {
00310   if (!isInitialized_) {
00311     // Set up prototypical InArgs
00312     Thyra::ModelEvaluatorBase::InArgsSetup<Scalar> inArgs;
00313     inArgs.setModelEvalDescription(this->description());
00314     inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_t );
00315     inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_x );
00316     inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_beta );
00317     // For ExplicitTaylorPolynomialStepper
00318     inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_x_poly );
00319     inArgs.set_Np(1); 
00320     if (isImplicit_) {
00321       inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_x_dot );
00322       inArgs.setSupports( Thyra::ModelEvaluatorBase::IN_ARG_alpha );
00323     }
00324     inArgs_ = inArgs;
00325     // Set up prototypical OutArgs
00326     Thyra::ModelEvaluatorBase::OutArgsSetup<Scalar> outArgs;
00327     outArgs.setModelEvalDescription(this->description());
00328     outArgs.setSupports( Thyra::ModelEvaluatorBase::OUT_ARG_f );
00329     outArgs.setSupports( Thyra::ModelEvaluatorBase::OUT_ARG_W_op );
00330     // For ExplicitTaylorPolynomialStepper
00331     outArgs.setSupports( Thyra::ModelEvaluatorBase::OUT_ARG_f_poly );
00332     outArgs.set_Np_Ng(Np_,Ng_);
00333     outArgs_ = outArgs;
00334     // Set up nominal values
00335     nominalValues_ = inArgs_;
00336     nominalValues_.set_t(Scalar(1.0));
00337     const RCP<VectorBase<Scalar> > x_ic = Thyra::createMember(x_space_);
00338     Thyra::V_S(Teuchos::outArg(*x_ic),Scalar(2.0));
00339     nominalValues_.set_x(x_ic);
00340     const RCP<VectorBase<Scalar> > p_ic = Thyra::createMember(p_space_);
00341     Thyra::V_S(Teuchos::outArg(*p_ic),Scalar(3.0));
00342     nominalValues_.set_p(0,p_ic);
00343     isInitialized_ = true;
00344   }
00345 }
00346 
00347 
00348 template<class Scalar>
00349 RCP<const std::vector<Thyra::ModelEvaluatorBase::InArgs<Scalar> > > StepperValidatorMockModel<Scalar>::getPassedInArgs() const
00350 {
00351   return passedInArgs_;
00352 }
00353 
00354 
00355 template<class Scalar>
00356 RCP<const Thyra::VectorSpaceBase<Scalar> > StepperValidatorMockModel<Scalar>::get_x_space() const
00357 {
00358   TEST_FOR_EXCEPT(!isInitialized_);
00359   return x_space_;
00360 }
00361 
00362 
00363 template<class Scalar>
00364 RCP<const Thyra::VectorSpaceBase<Scalar> > StepperValidatorMockModel<Scalar>::get_f_space() const
00365 {
00366   TEST_FOR_EXCEPT(!isInitialized_);
00367   return f_space_;
00368 }
00369 
00370 
00371 template<class Scalar>
00372 Thyra::ModelEvaluatorBase::InArgs<Scalar> StepperValidatorMockModel<Scalar>::getNominalValues() const
00373 {
00374   TEST_FOR_EXCEPT(!isInitialized_);
00375   return nominalValues_;
00376 }
00377 
00378 
00379 template<class Scalar>
00380 RCP<Thyra::LinearOpWithSolveBase<Scalar> > StepperValidatorMockModel<Scalar>::create_W() const
00381 {
00382   RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar> > W_factory = this->get_W_factory();
00383   RCP<Thyra::LinearOpBase<Scalar> > matrix = this->create_W_op();
00384   {
00385     // 01/20/09 tscoffe:  This is a total hack to provide a full rank matrix to
00386     // linearOpWithSolve because it ends up factoring the matrix during
00387     // initialization, which it really shouldn't do, or I'm doing something
00388     // wrong here.   The net effect is that I get exceptions thrown in
00389     // optimized mode due to the matrix being rank deficient unless I do this.
00390     RCP<Thyra::MultiVectorBase<Scalar> > multivec = Teuchos::rcp_dynamic_cast<Thyra::MultiVectorBase<Scalar> >(matrix,true);
00391     {
00392       RCP<Thyra::VectorBase<Scalar> > vec = Thyra::createMember(x_space_);
00393       {
00394         Thyra::DetachedVectorView<Scalar> vec_view( *vec );
00395         vec_view[0] = 1.0;
00396       }
00397       V_V(&*(multivec->col(0)),*vec);
00398     }
00399   }
00400   RCP<Thyra::LinearOpWithSolveBase<Scalar> > W = 
00401     Thyra::linearOpWithSolve<Scalar>(
00402       *W_factory,
00403       matrix
00404       );
00405   return W;
00406 }
00407 
00408 
00409 template<class Scalar>
00410 RCP<Thyra::LinearOpBase<Scalar> > StepperValidatorMockModel<Scalar>::create_W_op() const
00411 {
00412   RCP<Thyra::MultiVectorBase<Scalar> > matrix = Thyra::createMembers(x_space_, dim_);
00413   return(matrix);
00414 }
00415 
00416 
00417 template<class Scalar>
00418 RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar> > StepperValidatorMockModel<Scalar>::get_W_factory() const
00419 {
00420   RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> > W_factory = 
00421     Thyra::defaultSerialDenseLinearOpWithSolveFactory<Scalar>();
00422   return W_factory;
00423 }
00424 
00425 
00426 template<class Scalar>
00427 Thyra::ModelEvaluatorBase::InArgs<Scalar> StepperValidatorMockModel<Scalar>::createInArgs() const
00428 {
00429   TEST_FOR_EXCEPT(!isInitialized_);
00430   return inArgs_;
00431 }
00432 
00433 
00434 template<class Scalar>
00435 RCP<const Thyra::VectorSpaceBase<Scalar> > StepperValidatorMockModel<Scalar>::get_p_space(int l) const
00436 {
00437   TEST_FOR_EXCEPT(!isInitialized_);
00438   return p_space_;
00439 }
00440 
00441 
00442 template<class Scalar>
00443 RCP<const Teuchos::Array<std::string> > StepperValidatorMockModel<Scalar>::get_p_names(int l) const
00444 {
00445   RCP<Teuchos::Array<std::string> > p_strings = 
00446     Teuchos::rcp(new Teuchos::Array<std::string>());
00447   p_strings->push_back("Model Coefficient");
00448   return p_strings;
00449 }
00450 
00451 
00452 template<class Scalar>
00453 RCP<const Thyra::VectorSpaceBase<Scalar> > StepperValidatorMockModel<Scalar>::get_g_space(int j) const
00454 {
00455   return Teuchos::null;
00456 }
00457 
00458 
00459 template<class Scalar>
00460 void StepperValidatorMockModel<Scalar>::setParameterList(RCP<ParameterList> const& paramList)
00461 {
00462   TEST_FOR_EXCEPT( is_null(paramList) );
00463   paramList->validateParameters(*this->getValidParameters());
00464   Teuchos::readVerboseObjectSublist(&*paramList,this);
00465   bool test_isImplicit = paramList->get("Is Implicit",false);
00466   if (isImplicit_ != test_isImplicit) {
00467     isImplicit_ = test_isImplicit;
00468     isInitialized_ = false;
00469     this->initialize_();
00470   }
00471   setMyParamList(paramList);
00472 }
00473 template<class Scalar>
00474 RCP<const ParameterList> StepperValidatorMockModel<Scalar>::getValidParameters() const
00475 {
00476   static RCP<const ParameterList> validPL;
00477   if (is_null(validPL)) {
00478     RCP<ParameterList> pl = Teuchos::parameterList();
00479     pl->set("Is Implicit",false);
00480     Teuchos::setupVerboseObjectSublist(&*pl);
00481     validPL = pl;
00482   }
00483   return validPL;
00484 }
00485 template<class Scalar>
00486 Thyra::ModelEvaluatorBase::OutArgs<Scalar> StepperValidatorMockModel<Scalar>::createOutArgsImpl() const
00487 {
00488   TEST_FOR_EXCEPT(!isInitialized_);
00489   return outArgs_;
00490 }
00491 
00492 template<class Scalar>
00493 void StepperValidatorMockModel<Scalar>::evalModelImpl(
00494   const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
00495   const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs
00496   ) const
00497 {
00498   typedef Teuchos::ScalarTraits<Scalar> ST;
00499   passedInArgs_->push_back(inArgs);
00500   // Fill f with zeros.
00501   RCP<VectorBase<Scalar> > f_out = outArgs.get_f();
00502   if (!is_null(f_out)) {
00503     Thyra::V_S(Teuchos::outArg(*f_out),ST::zero());
00504   }
00505   if (outArgs.supports(Thyra::ModelEvaluatorBase::OUT_ARG_f_poly)) {
00506     RCP<Teuchos::Polynomial<VectorBase<Scalar> > > f_poly_out = outArgs.get_f_poly();
00507     if (!is_null(f_poly_out)) {
00508       //Thyra::V_S(Teuchos::outArg(*f_poly_out),ST::zero());
00509     }
00510   }
00511 
00512 }
00513 
00514 
00515 //
00516 // StepperValidator Definitions:
00517 //
00518 
00519 template<class Scalar>
00520 StepperValidator<Scalar>::StepperValidator()
00521 { 
00522   this->defaultInitializeAll_();
00523 }
00524 
00525 template<class Scalar>
00526 void StepperValidator<Scalar>::defaultInitializeAll_()
00527 {
00528   integratorBuilder_ = Teuchos::null;
00529   paramList_ = Teuchos::null;
00530   stepperName_ = "";
00531 }
00532 
00533 template<class Scalar>
00534   StepperValidator<Scalar>::~StepperValidator()
00535 { }
00536 
00537 template<class Scalar>
00538   void StepperValidator<Scalar>::setIntegratorBuilder(
00539     const RCP<IntegratorBuilder<Scalar> > &integratorBuilder
00540     )
00541 {
00542   TEUCHOS_ASSERT( !is_null(integratorBuilder) );
00543   integratorBuilder_ = integratorBuilder;
00544 }
00545 
00546 template<class Scalar>
00547   void StepperValidator<Scalar>::validateStepper() const
00548 {
00549   // Extract the name of the stepper for later
00550   {
00551     RCP<const ParameterList> pl = integratorBuilder_->getParameterList();
00552     if (is_null(pl)) {
00553       pl = integratorBuilder_->getValidParameters();
00554     }
00555     const Teuchos::ParameterList& stepperSelectionPL = pl->sublist("Stepper Settings").sublist("Stepper Selection");
00556     stepperName_ = stepperSelectionPL.get<std::string>("Stepper Type");
00557   }
00558   bool verbose = true;
00559   Array<bool> success_array;
00560   bool local_success = true;
00561 
00562   // Verify that the stepper passes parameters to the model in evalModel:
00563   local_success = true;
00564   try {
00565     this->validateIC_();
00566   }
00567   TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose,std::cerr,local_success);
00568   success_array.push_back(local_success);
00569 
00570   // Verify that the stepper states are correct 
00571   //   uninitialized => getTimeRange == invalidTimeRange
00572   //   initialized, but no step => getTimeRange.length() == 0, [t_ic,t_ic]
00573   //   initialized, step taken => getTimeRange.length() > 0
00574   local_success = true;
00575   try {
00576     this->validateStates_();
00577   }
00578   TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose,std::cerr,local_success);
00579   success_array.push_back(local_success);
00580 
00581   // Verify that getPoints(t_ic) returns the IC after initialization and after the first step
00582   local_success = true;
00583   try {
00584     this->validateGetIC_();
00585   }
00586   TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose,std::cerr,local_success);
00587   success_array.push_back(local_success);
00588 
00589   // Validate that the stepper supports getNodes, which is used by the Trailing Interpolation Buffer feature of the Integrator.
00590   local_success = true;
00591   try {
00592     this->validateGetNodes_();
00593   }
00594   TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose,std::cerr,local_success);
00595   success_array.push_back(local_success);
00596 
00597   // Verify that getPoints(t) returns the same vectors as getStepStatus
00598   // TODO
00599 
00600   bool global_success = true;
00601   for (int i=0 ; i < Teuchos::as<int>(success_array.size()) ; ++i) {
00602     global_success = global_success && success_array[i];
00603   }
00604 
00605   TEST_FOR_EXCEPTION( !global_success, std::logic_error,
00606       "Error!  StepperValidator:  The stepper " << stepperName_ << " did not pass stepper validation."
00607       );
00608 }
00609 
00610 template<class Scalar>
00611   void StepperValidator<Scalar>::setParameterList(RCP<Teuchos::ParameterList> const& paramList)
00612 {
00613   TEST_FOR_EXCEPT(true);
00614 }
00615 
00616 template<class Scalar>
00617   RCP<Teuchos::ParameterList> StepperValidator<Scalar>::getNonconstParameterList()
00618 {
00619   TEST_FOR_EXCEPT(true);
00620   return Teuchos::null;
00621 }
00622 
00623 template<class Scalar>
00624   RCP<Teuchos::ParameterList> StepperValidator<Scalar>::unsetParameterList()
00625 {
00626   TEST_FOR_EXCEPT(true);
00627   return Teuchos::null;
00628 }
00629 
00630 template<class Scalar>
00631   RCP<const Teuchos::ParameterList> StepperValidator<Scalar>::getValidParameters() const
00632 {
00633   TEST_FOR_EXCEPT(true);
00634   return Teuchos::null;
00635 }
00636 
00637 template<class Scalar>
00638   std::string StepperValidator<Scalar>::description() const
00639 {
00640   TEST_FOR_EXCEPT(true);
00641   return("");
00642 }
00643 
00644 template<class Scalar>
00645   void StepperValidator<Scalar>::describe(
00646     Teuchos::FancyOStream &out,
00647     const Teuchos::EVerbosityLevel verbLevel
00648     ) const
00649 {
00650   TEST_FOR_EXCEPT(true);
00651 }
00652 
00653 template<class Scalar>
00654 RCP<StepperBase<Scalar> > StepperValidator<Scalar>::getStepper_(
00655     const RCP<const Thyra::ModelEvaluator<Scalar> > model,
00656     const Thyra::ModelEvaluatorBase::InArgs<Scalar> &initialCondition,
00657     const RCP<Thyra::NonlinearSolverBase<Scalar> > nlSolver
00658     ) const
00659 {
00660   RCP<IntegratorBase<Scalar> > integrator = integratorBuilder_->create(model,initialCondition,nlSolver);
00661   RCP<StepperBase<Scalar> > stepper = integrator->getNonconstStepper();
00662   return stepper;
00663 }
00664 
00665 template<class Scalar>
00666 bool StepperValidator<Scalar>::isImplicitStepper_() const
00667 {
00668   RCP<const StepperBuilder<Scalar> > sb = integratorBuilder_->getStepperBuilder();
00669   RCP<StepperBase<Scalar> > stepper = sb->create(stepperName_);
00670   return stepper->isImplicit();
00671 }
00672 
00673 template<class Scalar>
00674 Thyra::ModelEvaluatorBase::InArgs<Scalar> StepperValidator<Scalar>::getSomeIC_(
00675     const Thyra::ModelEvaluator<Scalar>& model 
00676     ) const
00677 {
00678   // Set up some initial condition:
00679   Thyra::ModelEvaluatorBase::InArgs<Scalar> model_ic = model.createInArgs();
00680   if (model_ic.supports(Thyra::ModelEvaluatorBase::IN_ARG_t)) {
00681     Scalar time = Scalar(0.125);
00682     model_ic.set_t(time);
00683   }
00684   if (model_ic.supports(Thyra::ModelEvaluatorBase::IN_ARG_x)) {
00685     RCP<VectorBase<Scalar> > x_ic = Thyra::createMember(*(model.get_x_space()));
00686     {
00687       Thyra::DetachedVectorView<Scalar> x_ic_view( *x_ic );
00688       x_ic_view[0] = Scalar(10.0);
00689     }
00690     model_ic.set_x(x_ic);
00691   }
00692   for (int i=0 ; i<model_ic.Np() ; ++i) {
00693     RCP<VectorBase<Scalar> > p_ic = Thyra::createMember(*(model.get_p_space(i)));
00694     {
00695       Thyra::DetachedVectorView<Scalar> p_ic_view( *p_ic );
00696       p_ic_view[i] = Scalar(11.0+i);
00697     }
00698     model_ic.set_p(i,p_ic);
00699   }
00700   if (model_ic.supports(Thyra::ModelEvaluatorBase::IN_ARG_x_dot)) {
00701     RCP<VectorBase<Scalar> > xdot_ic = Thyra::createMember(*(model.get_x_space()));
00702     {
00703       Thyra::DetachedVectorView<Scalar> xdot_ic_view( *xdot_ic );
00704       xdot_ic_view[0] = Scalar(12.0);
00705     }
00706     model_ic.set_x_dot(xdot_ic);
00707   }
00708   return model_ic;
00709 }
00710 
00711 
00712 template<class Scalar>
00713 void StepperValidator<Scalar>::validateIC_() const
00714 {
00715   typedef Teuchos::ScalarTraits<Scalar> ST;
00716   // Determine if the stepper is implicit or not:
00717   bool isImplicit = this->isImplicitStepper_();
00718   RCP<StepperValidatorMockModel<Scalar> > model = 
00719     stepperValidatorMockModel<Scalar>(isImplicit);
00720   // Set up some initial condition:
00721   Thyra::ModelEvaluatorBase::InArgs<Scalar> stepper_ic = this->getSomeIC_(*model);
00722   // Create nonlinear solver (if needed)
00723   RCP<Thyra::NonlinearSolverBase<Scalar> > nlSolver;
00724   if (isImplicit) {
00725     nlSolver = Rythmos::timeStepNonlinearSolver<Scalar>();
00726   }
00727   RCP<StepperBase<Scalar> > stepper = this->getStepper_(model,stepper_ic,nlSolver);
00728   Scalar dt = Scalar(0.1);
00729   Scalar dt_taken = ST::nan();
00730   dt_taken = stepper->takeStep(dt,STEP_TYPE_FIXED);
00731   // Verify that the parameter got set on the model by asking the model for the
00732   // inArgs passed to it through evalModel
00733   RCP<const std::vector<Thyra::ModelEvaluatorBase::InArgs<Scalar> > >
00734     passedInArgs_ptr = model->getPassedInArgs();
00735   const std::vector<Thyra::ModelEvaluatorBase::InArgs<Scalar> >&
00736     passedInArgs = *passedInArgs_ptr;
00737   bool valid_t = false;
00738   // Technically this will fail for any Runge Kutta Butcher Tableau where:
00739   //  c(0) != 0 and c(s) != 1, where s = # of stages.
00740   if ( (passedInArgs[0].get_t() == stepper_ic.get_t()   )   || 
00741        (passedInArgs[0].get_t() == stepper_ic.get_t()+dt) )
00742   {
00743     valid_t = true;
00744   }
00745   TEST_FOR_EXCEPTION( !valid_t, std::logic_error,
00746     "Error!  StepperValidator::validateIC:  Time did not get correctly set on"
00747     " the model through StepperBase::setInitialCondition!"
00748     );
00749   // If a stepper uses a predictor, then the x passed into the model will not
00750   // be the same as the IC, so we can't check it.
00751   RCP<const VectorBase<Scalar> > p_out = passedInArgs[0].get_p(0);
00752   TEST_FOR_EXCEPTION( is_null(p_out), std::logic_error,
00753     "Error!  StepperValidator::validateIC:  Parameter 0 did not get set on the"
00754     " model through StepperBase::setInitialCondition!"
00755     );
00756   {
00757     Thyra::ConstDetachedVectorView<Scalar> p_out_view( *p_out );
00758     TEST_FOR_EXCEPTION( p_out_view[0] != Scalar(11.0), std::logic_error,
00759         "Error!  StepperValidator::validateIC:  Parameter 0 did not get set correctly on the model through StepperBase::setInitialCondition!"
00760         );
00761   }
00762   if (isImplicit) {
00763     // Each time integration method will approximate xdot according to its own
00764     // algorithm, so we can't really test it here.
00765   }
00766 }
00767 
00768 template<class Scalar>
00769 void StepperValidator<Scalar>::validateStates_() const
00770 {
00771   typedef Teuchos::ScalarTraits<Scalar> ST;
00772   RCP<const StepperBuilder<Scalar> > sb = integratorBuilder_->getStepperBuilder();
00773   RCP<StepperBase<Scalar> > stepper = sb->create(stepperName_);
00774   {
00775     // Uninitialized Stepper:
00776     TimeRange<Scalar> tr = stepper->getTimeRange();
00777     TEST_FOR_EXCEPTION( tr.isValid(), std::logic_error,
00778         "Error!  StepperValidator::validateStates:  Uninitialized stepper returned a valid time range!"
00779         );
00780     const StepStatus<Scalar> ss = stepper->getStepStatus();
00781     TEST_FOR_EXCEPTION( ss.stepStatus != STEP_STATUS_UNINITIALIZED, std::logic_error,
00782         "Error!  StepperValidator::validateStates:  Uninitialized stepper returned a valid step status!"
00783         );
00784   }
00785   bool isImplicit = stepper->isImplicit();
00786   RCP<StepperValidatorMockModel<Scalar> > model = 
00787     stepperValidatorMockModel<Scalar>(isImplicit);
00788   Thyra::ModelEvaluatorBase::InArgs<Scalar> stepper_ic = this->getSomeIC_(*model);
00789   stepper->setInitialCondition(stepper_ic);
00790   {
00791     // Has initial condition:
00792     TimeRange<Scalar> tr = stepper->getTimeRange();
00793     TEST_FOR_EXCEPTION( compareTimeValues(tr.lower(),tr.upper()) != 0, std::logic_error,
00794         "Error!  StepperValidator::validateStates:  Stepper with initial condition returned a non zero time range!"
00795         );
00796 //    const StepStatus<Scalar> ss = stepper->getStepStatus();
00797 //    TEST_FOR_EXCEPTION( ss.stepStatus != STEP_STATUS_UNKNOWN, std::logic_error,
00798 //        "Error!  StepperValidator::validateStates:  Stepper with initial condition did not return STEP_STATUS_UNKNOWN!"
00799 //        );
00800   }
00801   // 04/16/09 tscoffe:  Now we use the integratorBuilder to create a fully
00802   // initialized stepper which we can use for taking a step.  We can't just
00803   // continue setting them up because we don't know what they all require.
00804   RCP<Thyra::NonlinearSolverBase<Scalar> > nlSolver;
00805   if (isImplicit) {
00806     nlSolver = timeStepNonlinearSolver<Scalar>();
00807   }
00808   stepper = this->getStepper_(model,stepper_ic,nlSolver);
00809   {
00810     // Still has initial condition:
00811     TimeRange<Scalar> tr = stepper->getTimeRange();
00812     TEST_FOR_EXCEPTION( compareTimeValues(tr.lower(),tr.upper()) != 0, std::logic_error,
00813         "Error!  StepperValidator::validateStates:  Fully initialized stepper returned a non zero time range!"
00814         );
00815 //    const StepStatus<Scalar> ss = stepper->getStepStatus();
00816 //    TEST_FOR_EXCEPTION( ss.stepStatus != STEP_STATUS_UNKNOWN, std::logic_error,
00817 //        "Error!  StepperValidator::validateStates:  Fully initialized stepper, prior to taking a step, did not return STEP_STATUS_UNKNOWN!"
00818 //        );
00819   }
00820   Scalar dt = Scalar(0.1);
00821   Scalar dt_taken = ST::nan();
00822   dt_taken = stepper->takeStep(dt,STEP_TYPE_FIXED);
00823   {
00824     // Taken Step:
00825     TimeRange<Scalar> tr = stepper->getTimeRange();
00826     TEST_FOR_EXCEPTION( compareTimeValues(tr.lower(),tr.upper()) >= 0, std::logic_error,
00827         "Error!  StepperValidator::validateStates:  Stepper returned a zero (or invalid) time range after taking a step!"
00828         );
00829     const StepStatus<Scalar> ss = stepper->getStepStatus();
00830     TEST_FOR_EXCEPTION( ss.stepStatus != STEP_STATUS_CONVERGED, std::logic_error,
00831         "Error!  StepperValidator::validateStates:  Stepper did not return converged step status after taking a step!"
00832         );
00833   }
00834 }
00835 
00836 template<class Scalar>
00837 void StepperValidator<Scalar>::validateGetIC_() const
00838 {
00839   typedef Teuchos::ScalarTraits<Scalar> ST;
00840   typedef typename ScalarTraits<Scalar>::magnitudeType ScalarMag;
00841   // Determine if the stepper is implicit or not:
00842   bool isImplicit = this->isImplicitStepper_();
00843   RCP<StepperValidatorMockModel<Scalar> > model = 
00844     stepperValidatorMockModel<Scalar>(isImplicit);
00845   // Set up some initial condition:
00846   Thyra::ModelEvaluatorBase::InArgs<Scalar> stepper_ic = this->getSomeIC_(*model);
00847   // Create nonlinear solver (if needed)
00848   RCP<Thyra::NonlinearSolverBase<Scalar> > nlSolver;
00849   if (isImplicit) {
00850     nlSolver = Rythmos::timeStepNonlinearSolver<Scalar>();
00851   }
00852   RCP<StepperBase<Scalar> > stepper = this->getStepper_(model,stepper_ic,nlSolver);
00853   // Verify we can get the IC through getPoints prior to taking a step:
00854   {
00855     Array<Scalar> time_vec;
00856     Array<RCP<const VectorBase<Scalar> > > x_vec;
00857     Array<RCP<const VectorBase<Scalar> > > xdot_vec;
00858     Array<ScalarMag> accuracy_vec;
00859     time_vec.push_back(stepper_ic.get_t());
00860     stepper->getPoints(time_vec,&x_vec,&xdot_vec,&accuracy_vec);
00861     {
00862       Thyra::ConstDetachedVectorView<Scalar> x_view( *x_vec[0] );
00863       TEST_FOR_EXCEPTION( compareTimeValues<Scalar>(x_view[0],Scalar(10.0)) != 0,
00864         std::logic_error,
00865         "Error!  StepperValidator::validateGetIC:  Stepper did not return the initial"
00866         " condition for X through getPoints prior to taking a step!"
00867         );
00868     }
00869     if (isImplicit && !is_null(xdot_vec[0])) {
00870       Thyra::ConstDetachedVectorView<Scalar> xdot_view( *xdot_vec[0] );
00871       TEST_FOR_EXCEPTION( compareTimeValues<Scalar>(xdot_view[0],Scalar(12.0)) != 0,
00872         std::logic_error,
00873         "Error!  StepperValidator::validateGetIC:  Stepper did not return the initial"
00874         " condition for XDOT through getPoints prior to taking a step!"
00875         );
00876     }
00877   }
00878   // Verify we can get the IC through getPoints after taking a step:
00879   {
00880     Scalar dt = Scalar(0.1);
00881     Scalar dt_taken = ST::nan();
00882     dt_taken = stepper->takeStep(dt,STEP_TYPE_FIXED);
00883     Array<Scalar> time_vec;
00884     Array<RCP<const VectorBase<Scalar> > > x_vec;
00885     Array<RCP<const VectorBase<Scalar> > > xdot_vec;
00886     Array<ScalarMag> accuracy_vec;
00887     time_vec.push_back(stepper_ic.get_t());
00888     stepper->getPoints(time_vec,&x_vec,&xdot_vec,&accuracy_vec);
00889     {
00890       Thyra::ConstDetachedVectorView<Scalar> x_view( *x_vec[0] );
00891       TEST_FOR_EXCEPTION( compareTimeValues<Scalar>(x_view[0],Scalar(10.0)) != 0,
00892         std::logic_error,
00893         "Error!  StepperValidator::validateGetIC:  Stepper did not return the initial"
00894         " condition for X through getPoints after taking a step!"
00895         );
00896     }
00897     if (isImplicit && !is_null(xdot_vec[0])) {
00898       Thyra::ConstDetachedVectorView<Scalar> xdot_view( *xdot_vec[0] );
00899       TEST_FOR_EXCEPTION( compareTimeValues<Scalar>(xdot_view[0],Scalar(12.0)) != 0,
00900         std::logic_error,
00901         "Error!  StepperValidator::validateGetIC:  Stepper did not return the initial"
00902         " condition for XDOT through getPoints after taking a step!"
00903         );
00904     }
00905   }
00906 }
00907 
00908 
00909 template<class Scalar>
00910 void StepperValidator<Scalar>::validateGetNodes_() const
00911 {
00912   typedef Teuchos::ScalarTraits<Scalar> ST;
00913   // Create uninitialized stepper and verify we get no nodes back
00914   {
00915     RCP<const StepperBuilder<Scalar> > sb = integratorBuilder_->getStepperBuilder();
00916     RCP<StepperBase<Scalar> > stepper = sb->create(stepperName_);
00917     Array<Scalar> nodes;
00918     stepper->getNodes(&nodes);
00919     TEST_FOR_EXCEPTION( nodes.size() != 0, std::logic_error,
00920         "Error!  StepperValidator::validateGetNodes:  Uninitialized stepper returned non-empty node list!"
00921         );
00922   }
00923   // Create fully initialize stepper and verify we get back one node for IC
00924   bool isImplicit = this->isImplicitStepper_();
00925   RCP<StepperValidatorMockModel<Scalar> > model = 
00926     stepperValidatorMockModel<Scalar>(isImplicit);
00927   Thyra::ModelEvaluatorBase::InArgs<Scalar> stepper_ic = this->getSomeIC_(*model);
00928   RCP<Thyra::NonlinearSolverBase<Scalar> > nlSolver;
00929   if (isImplicit) {
00930     nlSolver = Rythmos::timeStepNonlinearSolver<Scalar>();
00931   }
00932   RCP<StepperBase<Scalar> > stepper = this->getStepper_(model,stepper_ic,nlSolver);
00933   {
00934     Array<Scalar> nodes;
00935     stepper->getNodes(&nodes);
00936     TEST_FOR_EXCEPTION( nodes.size() == 0, std::logic_error,
00937         "Error!  StepperValidator::validateGetNodes:  Initialized stepper returned empty node list!"
00938         );
00939     TEST_FOR_EXCEPTION( nodes.size() > 1, std::logic_error,
00940         "Error!  StepperValidator::validateGetNodes:  Initialized stepper returned node list with more than one node!"
00941         );
00942   }
00943   // Take a step with the stepper and verify we get back two nodes
00944   Scalar dt = Scalar(0.1);
00945   Scalar dt_taken = ST::nan();
00946   dt_taken = stepper->takeStep(dt,STEP_TYPE_FIXED);
00947   {
00948     Array<Scalar> nodes;
00949     stepper->getNodes(&nodes);
00950     TEST_FOR_EXCEPTION( nodes.size() == 0, std::logic_error,
00951         "Error!  StepperValidator::validateGetNodes:  After taking a step, stepper returned empty node list!"
00952         );
00953     TEST_FOR_EXCEPTION( nodes.size() == 1, std::logic_error,
00954         "Error!  StepperValidator::validateGetNodes:  After taking a step, stepper returned node list with only one node!"
00955         );
00956     TEST_FOR_EXCEPTION( nodes.size() > 2, std::logic_error,
00957         "Error!  StepperValidator::validateGetNodes:  After taking a step, stepper returned node list with more than two nodes!"
00958         );
00959   }
00960 }
00961 
00962 } // namespace Rythmos
00963 
00964 #endif //Rythmos_STEPPER_VALIDATOR_H

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