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