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., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
00023 // USA
00024 // Questions? Contact Todd S. Coffey (tscoffe@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 //@HEADER
00028 
00029 #ifndef Rythmos_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
00591   // and after the first step
00592   local_success = true;
00593   try {
00594     this->validateGetIC_();
00595   }
00596   TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose,std::cerr,local_success);
00597   success_array.push_back(local_success);
00598 
00599   // Verify that getInitialCondition() returns the IC after
00600   // setInitialCondition(...)
00601   local_success = true;
00602   try {
00603     this->validateGetIC2_();
00604   }
00605   TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose,std::cerr,local_success);
00606   success_array.push_back(local_success);
00607 
00608   // Validate that the stepper supports getNodes, which is used by
00609   // the Trailing Interpolation Buffer feature of the Integrator.
00610   local_success = true;
00611   try {
00612     this->validateGetNodes_();
00613   }
00614   TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose,std::cerr,local_success);
00615   success_array.push_back(local_success);
00616 
00617   // Verify that getPoints(t) returns the same vectors as getStepStatus
00618   // TODO
00619 
00620   bool global_success = true;
00621   for (int i=0 ; i < Teuchos::as<int>(success_array.size()) ; ++i) {
00622     global_success = global_success && success_array[i];
00623   }
00624 
00625   TEUCHOS_TEST_FOR_EXCEPTION( !global_success, std::logic_error,
00626       "Error!  StepperValidator:  The stepper " << stepperName_
00627       << " did not pass stepper validation.");
00628 }
00629 
00630 template<class Scalar>
00631   void StepperValidator<Scalar>::setParameterList(RCP<Teuchos::ParameterList> const& paramList)
00632 {
00633   TEUCHOS_TEST_FOR_EXCEPT(true);
00634 }
00635 
00636 template<class Scalar>
00637   RCP<Teuchos::ParameterList> StepperValidator<Scalar>::getNonconstParameterList()
00638 {
00639   TEUCHOS_TEST_FOR_EXCEPT(true);
00640   return Teuchos::null;
00641 }
00642 
00643 template<class Scalar>
00644   RCP<Teuchos::ParameterList> StepperValidator<Scalar>::unsetParameterList()
00645 {
00646   TEUCHOS_TEST_FOR_EXCEPT(true);
00647   return Teuchos::null;
00648 }
00649 
00650 template<class Scalar>
00651   RCP<const Teuchos::ParameterList> StepperValidator<Scalar>::getValidParameters() const
00652 {
00653   TEUCHOS_TEST_FOR_EXCEPT(true);
00654   return Teuchos::null;
00655 }
00656 
00657 template<class Scalar>
00658   std::string StepperValidator<Scalar>::description() const
00659 {
00660   TEUCHOS_TEST_FOR_EXCEPT(true);
00661   return("");
00662 }
00663 
00664 template<class Scalar>
00665   void StepperValidator<Scalar>::describe(
00666     Teuchos::FancyOStream &out,
00667     const Teuchos::EVerbosityLevel verbLevel
00668     ) const
00669 {
00670   TEUCHOS_TEST_FOR_EXCEPT(true);
00671 }
00672 
00673 template<class Scalar>
00674 RCP<StepperBase<Scalar> > StepperValidator<Scalar>::getStepper_(
00675     const RCP<Thyra::ModelEvaluator<Scalar> >& model,
00676     const Thyra::ModelEvaluatorBase::InArgs<Scalar>& initialCondition,
00677     const RCP<Thyra::NonlinearSolverBase<Scalar> >& nlSolver
00678     ) const
00679 {
00680   RCP<IntegratorBase<Scalar> > integrator = integratorBuilder_->create(model,initialCondition,nlSolver);
00681   RCP<StepperBase<Scalar> > stepper = integrator->getNonconstStepper();
00682   return stepper;
00683 }
00684 
00685 template<class Scalar>
00686 bool StepperValidator<Scalar>::isImplicitStepper_() const
00687 {
00688   RCP<const StepperBuilder<Scalar> > sb = integratorBuilder_->getStepperBuilder();
00689   RCP<StepperBase<Scalar> > stepper = sb->create(stepperName_);
00690   return stepper->isImplicit();
00691 }
00692 
00693 template<class Scalar>
00694 Thyra::ModelEvaluatorBase::InArgs<Scalar> StepperValidator<Scalar>::getSomeIC_(
00695     const Thyra::ModelEvaluator<Scalar>& model
00696     ) const
00697 {
00698   // Set up some initial condition:
00699   Thyra::ModelEvaluatorBase::InArgs<Scalar> model_ic = model.createInArgs();
00700   if (model_ic.supports(Thyra::ModelEvaluatorBase::IN_ARG_t)) {
00701     Scalar time = Scalar(0.125);
00702     model_ic.set_t(time);
00703   }
00704   if (model_ic.supports(Thyra::ModelEvaluatorBase::IN_ARG_x)) {
00705     RCP<VectorBase<Scalar> > x_ic = Thyra::createMember(*(model.get_x_space()));
00706     {
00707       Thyra::DetachedVectorView<Scalar> x_ic_view( *x_ic );
00708       x_ic_view[0] = Scalar(10.0);
00709     }
00710     model_ic.set_x(x_ic);
00711   }
00712   for (int i=0 ; i<model_ic.Np() ; ++i) {
00713     RCP<VectorBase<Scalar> > p_ic = Thyra::createMember(*(model.get_p_space(i)));
00714     {
00715       Thyra::DetachedVectorView<Scalar> p_ic_view( *p_ic );
00716       p_ic_view[i] = Scalar(11.0+i);
00717     }
00718     model_ic.set_p(i,p_ic);
00719   }
00720   if (model_ic.supports(Thyra::ModelEvaluatorBase::IN_ARG_x_dot)) {
00721     RCP<VectorBase<Scalar> > xdot_ic = Thyra::createMember(*(model.get_x_space()));
00722     {
00723       Thyra::DetachedVectorView<Scalar> xdot_ic_view( *xdot_ic );
00724       xdot_ic_view[0] = Scalar(12.0);
00725     }
00726     model_ic.set_x_dot(xdot_ic);
00727   }
00728   return model_ic;
00729 }
00730 
00731 
00732 template<class Scalar>
00733 void StepperValidator<Scalar>::validateIC_() const
00734 {
00735   typedef Teuchos::ScalarTraits<Scalar> ST;
00736   // Determine if the stepper is implicit or not:
00737   bool isImplicit = this->isImplicitStepper_();
00738   RCP<StepperValidatorMockModel<Scalar> > model =
00739     stepperValidatorMockModel<Scalar>(isImplicit);
00740   // Set up some initial condition:
00741   Thyra::ModelEvaluatorBase::InArgs<Scalar> stepper_ic = this->getSomeIC_(*model);
00742   // Create nonlinear solver (if needed)
00743   RCP<Thyra::NonlinearSolverBase<Scalar> > nlSolver;
00744   if (isImplicit) {
00745     nlSolver = Rythmos::timeStepNonlinearSolver<Scalar>();
00746   }
00747   RCP<StepperBase<Scalar> > stepper = this->getStepper_(model,stepper_ic,nlSolver);
00748   Scalar dt = Scalar(0.1);
00749   stepper->takeStep(dt,STEP_TYPE_FIXED);
00750   // Verify that the parameter got set on the model by asking the model for the
00751   // inArgs passed to it through evalModel
00752   RCP<const std::vector<Thyra::ModelEvaluatorBase::InArgs<Scalar> > >
00753     passedInArgs_ptr = model->getPassedInArgs();
00754   const std::vector<Thyra::ModelEvaluatorBase::InArgs<Scalar> >&
00755     passedInArgs = *passedInArgs_ptr;
00756   bool valid_t = false;
00757   // Technically this will fail for any Runge Kutta Butcher Tableau where:
00758   //  c(0) != 0 and c(s) != 1, where s = # of stages.
00759   if ( (passedInArgs[0].get_t() == stepper_ic.get_t()   )   ||
00760        (passedInArgs[0].get_t() == stepper_ic.get_t()+dt) )
00761   {
00762     valid_t = true;
00763   }
00764   TEUCHOS_TEST_FOR_EXCEPTION( !valid_t, std::logic_error,
00765     "Error!  StepperValidator::validateIC:  Time did not get correctly set on"
00766     " the model through StepperBase::setInitialCondition!"
00767     );
00768   // If a stepper uses a predictor, then the x passed into the model will not
00769   // be the same as the IC, so we can't check it.
00770   RCP<const VectorBase<Scalar> > p_out = passedInArgs[0].get_p(0);
00771   TEUCHOS_TEST_FOR_EXCEPTION( is_null(p_out), std::logic_error,
00772     "Error!  StepperValidator::validateIC:  Parameter 0 did not get set on the"
00773     " model through StepperBase::setInitialCondition!"
00774     );
00775   {
00776     Thyra::ConstDetachedVectorView<Scalar> p_out_view( *p_out );
00777     TEUCHOS_TEST_FOR_EXCEPTION( p_out_view[0] != Scalar(11.0), std::logic_error,
00778         "Error!  StepperValidator::validateIC:  Parameter 0 did not get set correctly on the model through StepperBase::setInitialCondition!"
00779         );
00780   }
00781   if (isImplicit) {
00782     // Each time integration method will approximate xdot according to its own
00783     // algorithm, so we can't really test it here.
00784   }
00785 }
00786 
00787 template<class Scalar>
00788 void StepperValidator<Scalar>::validateStates_() const
00789 {
00790   typedef Teuchos::ScalarTraits<Scalar> ST;
00791   RCP<const StepperBuilder<Scalar> > sb =
00792     integratorBuilder_->getStepperBuilder();
00793   RCP<StepperBase<Scalar> > stepper = sb->create(stepperName_);
00794   {
00795     // Uninitialized Stepper:
00796     TimeRange<Scalar> tr = stepper->getTimeRange();
00797     TEUCHOS_TEST_FOR_EXCEPTION( tr.isValid(), std::logic_error,
00798         "Error!  StepperValidator::validateStates:  Uninitialized "
00799         "stepper returned a valid time range!"
00800         );
00801     const StepStatus<Scalar> ss = stepper->getStepStatus();
00802     TEUCHOS_TEST_FOR_EXCEPTION( ss.stepStatus != STEP_STATUS_UNINITIALIZED,
00803         std::logic_error,
00804         "Error!  StepperValidator::validateStates:  Uninitialized "
00805         "stepper returned a valid step status!"
00806         );
00807   }
00808   bool isImplicit = stepper->isImplicit();
00809   RCP<StepperValidatorMockModel<Scalar> > model =
00810     stepperValidatorMockModel<Scalar>(isImplicit);
00811   Thyra::ModelEvaluatorBase::InArgs<Scalar> stepper_ic =
00812     this->getSomeIC_(*model);
00813   stepper->setInitialCondition(stepper_ic);
00814   {
00815     // Has initial condition:
00816     TimeRange<Scalar> tr = stepper->getTimeRange();
00817     TEUCHOS_TEST_FOR_EXCEPTION( compareTimeValues(tr.lower(),tr.upper()) != 0,
00818         std::logic_error,
00819         "Error!  StepperValidator::validateStates:  Stepper with "
00820         "initial condition returned a non zero time range!"
00821         );
00822 //    const StepStatus<Scalar> ss = stepper->getStepStatus();
00823 //    TEUCHOS_TEST_FOR_EXCEPTION( ss.stepStatus != STEP_STATUS_UNKNOWN,
00824 //        std::logic_error,
00825 //        "Error!  StepperValidator::validateStates:  Stepper with "
00826 //        "initial condition did not return STEP_STATUS_UNKNOWN!"
00827 //        );
00828   }
00829   // 04/16/09 tscoffe:  Now we use the integratorBuilder to create a fully
00830   // initialized stepper which we can use for taking a step.  We can't just
00831   // continue setting them up because we don't know what they all require.
00832   RCP<Thyra::NonlinearSolverBase<Scalar> > nlSolver;
00833   if (isImplicit) {
00834     nlSolver = timeStepNonlinearSolver<Scalar>();
00835   }
00836   stepper = this->getStepper_(model,stepper_ic,nlSolver);
00837   {
00838     // Still has initial condition:
00839     TimeRange<Scalar> tr = stepper->getTimeRange();
00840     TEUCHOS_TEST_FOR_EXCEPTION( compareTimeValues(tr.lower(),tr.upper()) != 0,
00841         std::logic_error,
00842         "Error!  StepperValidator::validateStates:  Fully initialized "
00843         "stepper returned a non zero time range!"
00844         );
00845 //    const StepStatus<Scalar> ss = stepper->getStepStatus();
00846 //    TEUCHOS_TEST_FOR_EXCEPTION( ss.stepStatus != STEP_STATUS_UNKNOWN,
00847 //        std::logic_error,
00848 //        "Error!  StepperValidator::validateStates:  Fully initialized "
00849 //        "stepper, prior to taking a step, did not return STEP_STATUS_UNKNOWN!"
00850 //        );
00851   }
00852   Scalar dt = Scalar(0.1);
00853   stepper->takeStep(dt,STEP_TYPE_FIXED);
00854   {
00855     // Taken Step:
00856     TimeRange<Scalar> tr = stepper->getTimeRange();
00857     TEUCHOS_TEST_FOR_EXCEPTION( compareTimeValues(tr.lower(),tr.upper()) >= 0,
00858         std::logic_error,
00859         "Error!  StepperValidator::validateStates:  Stepper returned "
00860         "a zero (or invalid) time range after taking a step!"
00861         );
00862     const StepStatus<Scalar> ss = stepper->getStepStatus();
00863     TEUCHOS_TEST_FOR_EXCEPTION( ss.stepStatus != STEP_STATUS_CONVERGED,
00864         std::logic_error,
00865         "Error!  StepperValidator::validateStates:  Stepper did not "
00866         "return converged step status after taking a step!"
00867         );
00868   }
00869 }
00870 
00871 template<class Scalar>
00872 void StepperValidator<Scalar>::validateGetIC_() const
00873 {
00874   typedef Teuchos::ScalarTraits<Scalar> ST;
00875   typedef typename ScalarTraits<Scalar>::magnitudeType ScalarMag;
00876   // Determine if the stepper is implicit or not:
00877   bool isImplicit = this->isImplicitStepper_();
00878   RCP<StepperValidatorMockModel<Scalar> > model =
00879     stepperValidatorMockModel<Scalar>(isImplicit);
00880   // Set up some initial condition:
00881   Thyra::ModelEvaluatorBase::InArgs<Scalar> stepper_ic = this->getSomeIC_(*model);
00882   // Create nonlinear solver (if needed)
00883   RCP<Thyra::NonlinearSolverBase<Scalar> > nlSolver;
00884   if (isImplicit) {
00885     nlSolver = Rythmos::timeStepNonlinearSolver<Scalar>();
00886   }
00887   RCP<StepperBase<Scalar> > stepper = this->getStepper_(model,stepper_ic,nlSolver);
00888   // Verify we can get the IC through getPoints prior to taking a step:
00889   {
00890     Array<Scalar> time_vec;
00891     Array<RCP<const VectorBase<Scalar> > > x_vec;
00892     Array<RCP<const VectorBase<Scalar> > > xdot_vec;
00893     Array<ScalarMag> accuracy_vec;
00894     time_vec.push_back(stepper_ic.get_t());
00895     stepper->getPoints(time_vec,&x_vec,&xdot_vec,&accuracy_vec);
00896     {
00897       Thyra::ConstDetachedVectorView<Scalar> x_view( *x_vec[0] );
00898       TEUCHOS_TEST_FOR_EXCEPTION( compareTimeValues<Scalar>(x_view[0],Scalar(10.0)) != 0,
00899         std::logic_error,
00900         "Error!  StepperValidator::validateGetIC:  Stepper did not return the initial"
00901         " condition for X through getPoints prior to taking a step!"
00902         );
00903     }
00904     if (isImplicit && !is_null(xdot_vec[0])) {
00905       Thyra::ConstDetachedVectorView<Scalar> xdot_view( *xdot_vec[0] );
00906       TEUCHOS_TEST_FOR_EXCEPTION( compareTimeValues<Scalar>(xdot_view[0],Scalar(12.0)) != 0,
00907         std::logic_error,
00908         "Error!  StepperValidator::validateGetIC:  Stepper did not return the initial"
00909         " condition for XDOT through getPoints prior to taking a step!"
00910         );
00911     }
00912   }
00913   // Verify we can get the IC through getPoints after taking a step:
00914   {
00915     Scalar dt = Scalar(0.1);
00916     stepper->takeStep(dt,STEP_TYPE_FIXED);
00917     Array<Scalar> time_vec;
00918     Array<RCP<const VectorBase<Scalar> > > x_vec;
00919     Array<RCP<const VectorBase<Scalar> > > xdot_vec;
00920     Array<ScalarMag> accuracy_vec;
00921     time_vec.push_back(stepper_ic.get_t());
00922     stepper->getPoints(time_vec,&x_vec,&xdot_vec,&accuracy_vec);
00923     {
00924       Thyra::ConstDetachedVectorView<Scalar> x_view( *x_vec[0] );
00925       TEUCHOS_TEST_FOR_EXCEPTION( compareTimeValues<Scalar>(x_view[0],Scalar(10.0)) != 0,
00926         std::logic_error,
00927         "Error!  StepperValidator::validateGetIC:  Stepper did not return the initial"
00928         " condition for X through getPoints after taking a step!"
00929         );
00930     }
00931     if (isImplicit && !is_null(xdot_vec[0])) {
00932       Thyra::ConstDetachedVectorView<Scalar> xdot_view( *xdot_vec[0] );
00933       TEUCHOS_TEST_FOR_EXCEPTION( compareTimeValues<Scalar>(xdot_view[0],Scalar(12.0)) != 0,
00934         std::logic_error,
00935         "Error!  StepperValidator::validateGetIC:  Stepper did not return the initial"
00936         " condition for XDOT through getPoints after taking a step!"
00937         );
00938     }
00939   }
00940 }
00941 
00942 
00943 template<class Scalar>
00944 void StepperValidator<Scalar>::validateGetIC2_() const
00945 {
00946   typedef Teuchos::ScalarTraits<Scalar> ST;
00947   typedef typename ScalarTraits<Scalar>::magnitudeType ScalarMag;
00948   // Determine if the stepper is implicit or not:
00949   bool isImplicit = this->isImplicitStepper_();
00950   RCP<StepperValidatorMockModel<Scalar> > model =
00951     stepperValidatorMockModel<Scalar>(isImplicit);
00952   // Set up some initial condition:
00953   Thyra::ModelEvaluatorBase::InArgs<Scalar> stepper_ic = this->getSomeIC_(*model);
00954   // Create nonlinear solver (if needed)
00955   RCP<Thyra::NonlinearSolverBase<Scalar> > nlSolver;
00956   if (isImplicit) {
00957     nlSolver = Rythmos::timeStepNonlinearSolver<Scalar>();
00958   }
00959   RCP<StepperBase<Scalar> > stepper = this->getStepper_(model,stepper_ic,nlSolver);
00960   // Verify we can get the IC back through getInitialCondition:
00961   Thyra::ModelEvaluatorBase::InArgs<Scalar> new_ic = stepper->getInitialCondition();
00962   //TEUCHOS_ASSERT( new_ic == stepper_ic );
00963   // Verify new_ic == stepper_ic
00964   {
00965     TEUCHOS_ASSERT( new_ic.get_t() == stepper_ic.get_t() );
00966     TEUCHOS_ASSERT( new_ic.get_x() == stepper_ic.get_x() );
00967     for (int i=0 ; i<stepper_ic.Np() ; ++i) {
00968       TEUCHOS_ASSERT( new_ic.get_p(i) == stepper_ic.get_p(i) );
00969     }
00970     if (isImplicit) {
00971       TEUCHOS_ASSERT( new_ic.get_x_dot() == stepper_ic.get_x_dot() );
00972     }
00973   }
00974 }
00975 
00976 
00977 template<class Scalar>
00978 void StepperValidator<Scalar>::validateGetNodes_() const
00979 {
00980   typedef Teuchos::ScalarTraits<Scalar> ST;
00981   // Create uninitialized stepper and verify we get no nodes back
00982   {
00983     RCP<const StepperBuilder<Scalar> > sb = integratorBuilder_->getStepperBuilder();
00984     RCP<StepperBase<Scalar> > stepper = sb->create(stepperName_);
00985     Array<Scalar> nodes;
00986     stepper->getNodes(&nodes);
00987     TEUCHOS_TEST_FOR_EXCEPTION( nodes.size() != 0, std::logic_error,
00988         "Error!  StepperValidator::validateGetNodes:  Uninitialized stepper returned non-empty node list!"
00989         );
00990   }
00991   // Create fully initialize stepper and verify we get back one node for IC
00992   bool isImplicit = this->isImplicitStepper_();
00993   RCP<StepperValidatorMockModel<Scalar> > model =
00994     stepperValidatorMockModel<Scalar>(isImplicit);
00995   Thyra::ModelEvaluatorBase::InArgs<Scalar> stepper_ic = this->getSomeIC_(*model);
00996   RCP<Thyra::NonlinearSolverBase<Scalar> > nlSolver;
00997   if (isImplicit) {
00998     nlSolver = Rythmos::timeStepNonlinearSolver<Scalar>();
00999   }
01000   RCP<StepperBase<Scalar> > stepper = this->getStepper_(model,stepper_ic,nlSolver);
01001   {
01002     Array<Scalar> nodes;
01003     stepper->getNodes(&nodes);
01004     TEUCHOS_TEST_FOR_EXCEPTION( nodes.size() == 0, std::logic_error,
01005         "Error!  StepperValidator::validateGetNodes:  Initialized stepper returned empty node list!"
01006         );
01007     TEUCHOS_TEST_FOR_EXCEPTION( nodes.size() > 1, std::logic_error,
01008         "Error!  StepperValidator::validateGetNodes:  Initialized stepper returned node list with more than one node!"
01009         );
01010   }
01011   // Take a step with the stepper and verify we get back two nodes
01012   Scalar dt = Scalar(0.1);
01013   stepper->takeStep(dt,STEP_TYPE_FIXED);
01014   {
01015     Array<Scalar> nodes;
01016     stepper->getNodes(&nodes);
01017     TEUCHOS_TEST_FOR_EXCEPTION( nodes.size() == 0, std::logic_error,
01018         "Error!  StepperValidator::validateGetNodes:  After taking a step, stepper returned empty node list!"
01019         );
01020     TEUCHOS_TEST_FOR_EXCEPTION( nodes.size() == 1, std::logic_error,
01021         "Error!  StepperValidator::validateGetNodes:  After taking a step, stepper returned node list with only one node!"
01022         );
01023     TEUCHOS_TEST_FOR_EXCEPTION( nodes.size() > 2, std::logic_error,
01024         "Error!  StepperValidator::validateGetNodes:  After taking a step, stepper returned node list with more than two nodes!"
01025         );
01026   }
01027 }
01028 
01029 } // namespace Rythmos
01030 
01031 #endif //Rythmos_STEPPER_VALIDATOR_H
 All Classes Functions Variables Typedefs Friends