Rythmos_ForwardSensitivityExplicitModelEvaluator.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_FORWARD_SENSITIVITY_EXPLICIT_MODEL_EVALUATOR_HPP
00030 #define RYTHMOS_FORWARD_SENSITIVITY_EXPLICIT_MODEL_EVALUATOR_HPP
00031 
00032 
00033 #include "Rythmos_IntegratorBase.hpp"
00034 #include "Rythmos_ForwardSensitivityModelEvaluatorBase.hpp"
00035 #include "Thyra_ModelEvaluator.hpp" // Interface
00036 #include "Thyra_StateFuncModelEvaluatorBase.hpp" // Implementation
00037 #include "Thyra_ModelEvaluatorDelegatorBase.hpp"
00038 #include "Thyra_DefaultMultiVectorProductVectorSpace.hpp"
00039 #include "Thyra_DefaultMultiVectorProductVector.hpp"
00040 #include "Thyra_DefaultMultiVectorLinearOpWithSolve.hpp"
00041 
00042 
00043 namespace Rythmos {
00044 
00045 
00121 template<class Scalar>
00122 class ForwardSensitivityExplicitModelEvaluator
00123   : virtual public ForwardSensitivityModelEvaluatorBase<Scalar>
00124 {
00125 public:
00126 
00129 
00131   ForwardSensitivityExplicitModelEvaluator();
00132 
00148   void initializeStructure(
00149     const RCP<const Thyra::ModelEvaluator<Scalar> > &stateModel,
00150     const int p_index
00151     );
00152   
00154   RCP<const Thyra::ModelEvaluator<Scalar> >
00155   getStateModel() const;
00156   
00158   int get_p_index() const;
00159 
00173   void initializePointState(
00174     const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint
00175     );
00176 
00178 
00181 
00183   RCP<const Thyra::VectorSpaceBase<Scalar> > get_x_space() const;
00185   RCP<const Thyra::VectorSpaceBase<Scalar> > get_f_space() const;
00187   Thyra::ModelEvaluatorBase::InArgs<Scalar> getNominalValues() const;
00189   RCP<Thyra::LinearOpWithSolveBase<Scalar> > create_W() const;
00191   Thyra::ModelEvaluatorBase::InArgs<Scalar> createInArgs() const;
00192 
00194 
00195 private:
00196 
00199 
00201   Thyra::ModelEvaluatorBase::OutArgs<Scalar> createOutArgsImpl() const;
00203   void evalModelImpl(
00204     const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
00205     const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs
00206     ) const;
00207 
00209 
00210 private:
00211 
00212   // /////////////////////////
00213   // Private data members
00214 
00215   RCP<const Thyra::ModelEvaluator<Scalar> > stateModel_;
00216   int p_index_;
00217   int np_;
00218 
00219   RCP<const Thyra::DefaultMultiVectorProductVectorSpace<Scalar> > s_bar_space_;
00220   RCP<const Thyra::DefaultMultiVectorProductVectorSpace<Scalar> > f_sens_space_;
00221 
00222   Thyra::ModelEvaluatorBase::InArgs<Scalar> nominalValues_;
00223 
00224   mutable Thyra::ModelEvaluatorBase::InArgs<Scalar> stateBasePoint_;
00225 
00226   mutable RCP<Thyra::LinearOpBase<Scalar> > DfDx_;
00227   mutable RCP<Thyra::MultiVectorBase<Scalar> > DfDp_;
00228 
00229   mutable RCP<Thyra::LinearOpBase<Scalar> > DfDx_compute_;
00230   mutable RCP<Thyra::MultiVectorBase<Scalar> > DfDp_compute_;
00231 
00232   // /////////////////////////
00233   // Private member functions
00234 
00235   void wrapNominalValuesAndBounds();
00236 
00237   void computeDerivativeMatrices(
00238     const Thyra::ModelEvaluatorBase::InArgs<Scalar> &point
00239     ) const;
00240   
00241 };
00242 
00243 
00244 // /////////////////////////////////
00245 // Implementations
00246 
00247 
00248 // Constructors/Intializers/Accessors
00249 
00250 // Nonmember constructor
00251 template<class Scalar>
00252 RCP<ForwardSensitivityExplicitModelEvaluator<Scalar> > forwardSensitivityExplicitModelEvaluator()
00253 {
00254   RCP<ForwardSensitivityExplicitModelEvaluator<Scalar> > fseme = rcp(new ForwardSensitivityExplicitModelEvaluator<Scalar> );
00255   return fseme;
00256 }
00257 
00258 
00259 template<class Scalar>
00260 ForwardSensitivityExplicitModelEvaluator<Scalar>::ForwardSensitivityExplicitModelEvaluator()
00261   : p_index_(0), np_(-1)
00262 {}
00263 
00264 
00265 template<class Scalar>
00266 void ForwardSensitivityExplicitModelEvaluator<Scalar>::initializeStructure(
00267   const RCP<const Thyra::ModelEvaluator<Scalar> > &stateModel,
00268   const int p_index
00269   )
00270 {
00271 
00272   typedef Thyra::ModelEvaluatorBase MEB;
00273 
00274   //
00275   // Validate input
00276   //
00277 
00278   TEST_FOR_EXCEPT( is_null(stateModel) );
00279   TEST_FOR_EXCEPTION(
00280     !( 0 <= p_index && p_index < stateModel->Np() ), std::logic_error,
00281     "Error, p_index does not fall in the range [0,"<<(stateModel->Np()-1)<<"]!" );
00282   // ToDo: Validate support for DfDp!
00283 
00284   //
00285   // Set the input objects
00286   //
00287 
00288   stateModel_ = stateModel;
00289   p_index_ = p_index;
00290   np_ = stateModel_->get_p_space(p_index)->dim();
00291 
00292   //
00293   // Create the structure of the model
00294   //
00295 
00296   s_bar_space_ = Thyra::multiVectorProductVectorSpace(
00297     stateModel_->get_x_space(), np_
00298     );
00299 
00300   f_sens_space_ = Thyra::multiVectorProductVectorSpace(
00301     stateModel_->get_f_space(), np_
00302     );
00303 
00304   nominalValues_ = this->createInArgs();
00305 
00306   //
00307   // Wipe out matrix storage
00308   //
00309 
00310   stateBasePoint_ = MEB::InArgs<Scalar>();
00311   DfDx_ = Teuchos::null;
00312   DfDp_ = Teuchos::null;
00313   DfDx_compute_ = Teuchos::null;
00314   DfDp_compute_ = Teuchos::null;
00315 
00316 }
00317 
00318 
00319 template<class Scalar>
00320 RCP<const Thyra::ModelEvaluator<Scalar> >
00321 ForwardSensitivityExplicitModelEvaluator<Scalar>::getStateModel() const
00322 {
00323   return stateModel_;
00324 }
00325 
00326 
00327 template<class Scalar>
00328 int ForwardSensitivityExplicitModelEvaluator<Scalar>::get_p_index() const
00329 {
00330   return p_index_;
00331 }
00332 
00333 
00334 template<class Scalar>
00335 void ForwardSensitivityExplicitModelEvaluator<Scalar>::initializePointState(
00336   const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint
00337   )
00338 {
00339 
00340   typedef Thyra::ModelEvaluatorBase MEB;
00341 
00342 #ifdef RYTHMOS_DEBUG
00343   TEST_FOR_EXCEPTION(
00344     is_null(stateModel_), std::logic_error,
00345     "Error, you must call intializeStructure(...) before you call initializePointState(...)"
00346     );
00347   TEST_FOR_EXCEPT( is_null(stateBasePoint.get_x()) );
00348   TEST_FOR_EXCEPT( is_null(stateBasePoint.get_p(p_index_)) );
00349   // What about the other parameter values?  We really can't say anything
00350   // about these and we can't check them.  They can be null just fine.
00351 #endif
00352 
00353   stateBasePoint_ = stateBasePoint;
00354 
00355   // Set whatever derivatives where passed in.  If an input in null, then the
00356   // member will be null and the null linear operators will be computed later
00357   // just in time.
00358 
00359   wrapNominalValuesAndBounds();
00360 
00361 }
00362 
00363 
00364 // Public functions overridden from ModelEvaulator
00365 
00366 
00367 template<class Scalar>
00368 RCP<const Thyra::VectorSpaceBase<Scalar> >
00369 ForwardSensitivityExplicitModelEvaluator<Scalar>::get_x_space() const
00370 {
00371   return s_bar_space_;
00372 }
00373 
00374 
00375 template<class Scalar>
00376 RCP<const Thyra::VectorSpaceBase<Scalar> >
00377 ForwardSensitivityExplicitModelEvaluator<Scalar>::get_f_space() const
00378 {
00379   return f_sens_space_;
00380 }
00381 
00382 
00383 template<class Scalar>
00384 Thyra::ModelEvaluatorBase::InArgs<Scalar>
00385 ForwardSensitivityExplicitModelEvaluator<Scalar>::getNominalValues() const
00386 {
00387   return nominalValues_;
00388 }
00389 
00390 
00391 template<class Scalar>
00392 RCP<Thyra::LinearOpWithSolveBase<Scalar> >
00393 ForwardSensitivityExplicitModelEvaluator<Scalar>::create_W() const
00394 {
00395   return Thyra::multiVectorLinearOpWithSolve<Scalar>();
00396 }
00397 
00398 
00399 template<class Scalar>
00400 Thyra::ModelEvaluatorBase::InArgs<Scalar>
00401 ForwardSensitivityExplicitModelEvaluator<Scalar>::createInArgs() const
00402 {
00403   TEUCHOS_ASSERT( !is_null(stateModel_) );
00404   typedef Thyra::ModelEvaluatorBase MEB;
00405   MEB::InArgs<Scalar> stateModelInArgs = stateModel_->createInArgs();
00406   MEB::InArgsSetup<Scalar> inArgs;
00407   inArgs.setModelEvalDescription(this->description());
00408   inArgs.setSupports( MEB::IN_ARG_x );
00409   inArgs.setSupports( MEB::IN_ARG_t );
00410   //inArgs.setSupports( MEB::IN_ARG_beta,
00411   //  stateModelInArgs.supports(MEB::IN_ARG_beta) );
00412   return inArgs;
00413 }
00414 
00415 
00416 // Private functions overridden from ModelEvaulatorDefaultBase
00417 
00418 
00419 template<class Scalar>
00420 Thyra::ModelEvaluatorBase::OutArgs<Scalar>
00421 ForwardSensitivityExplicitModelEvaluator<Scalar>::createOutArgsImpl() const
00422 {
00423   TEUCHOS_ASSERT( !is_null(stateModel_) );
00424   typedef Thyra::ModelEvaluatorBase MEB;
00425 
00426   MEB::OutArgs<Scalar> stateModelOutArgs = stateModel_->createOutArgs();
00427   MEB::OutArgsSetup<Scalar> outArgs;
00428 
00429   outArgs.setModelEvalDescription(this->description());
00430 
00431   outArgs.setSupports(MEB::OUT_ARG_f);
00432 
00433   return outArgs;
00434 
00435 }
00436 
00437 
00438 template<class Scalar>
00439 void ForwardSensitivityExplicitModelEvaluator<Scalar>::evalModelImpl(
00440   const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
00441   const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs
00442   ) const
00443 {
00444 
00445   using Teuchos::rcp_dynamic_cast;
00446   typedef Teuchos::ScalarTraits<Scalar> ST;
00447   typedef Thyra::ModelEvaluatorBase MEB;
00448   typedef Teuchos::VerboseObjectTempState<Thyra::ModelEvaluatorBase> VOTSME;
00449 
00450   THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_GEN_BEGIN(
00451     "ForwardSensitivityExplicitModelEvaluator", inArgs, outArgs, Teuchos::null );
00452 
00453   //
00454   // Update the derivative matrices if they are not already updated for the
00455   // given time!.
00456   //
00457   
00458   {
00459 #ifdef ENABLE_RYTHMOS_TIMERS
00460     TEUCHOS_FUNC_TIME_MONITOR("Rythmos:ForwardSensitivityExplicitModelEvaluator::evalModel: computeMatrices");
00461 #endif
00462     computeDerivativeMatrices(inArgs);
00463   }
00464 
00465   //
00466   // InArgs
00467   //
00468 
00469   RCP<const Thyra::DefaultMultiVectorProductVector<Scalar> >
00470     s_bar = rcp_dynamic_cast<const Thyra::DefaultMultiVectorProductVector<Scalar> >(
00471       inArgs.get_x().assert_not_null(), true
00472       );
00473   RCP<const Thyra::MultiVectorBase<Scalar> >
00474     S = s_bar->getMultiVector();
00475   
00476   //
00477   // OutArgs
00478   //
00479 
00480   RCP<Thyra::DefaultMultiVectorProductVector<Scalar> >
00481     f_sens = rcp_dynamic_cast<Thyra::DefaultMultiVectorProductVector<Scalar> >(
00482       outArgs.get_f(), true
00483       );
00484 
00485   RCP<Thyra::MultiVectorBase<Scalar> >
00486     F_sens = f_sens->getNonconstMultiVector().assert_not_null();
00487 
00488   //
00489   // Compute the requested functions
00490   //
00491 
00492   if(!is_null(F_sens)) {
00493 
00494 #ifdef ENABLE_RYTHMOS_TIMERS
00495     TEUCHOS_FUNC_TIME_MONITOR("Rythmos:ForwardSensitivityExplicitModelEvaluator::evalModel: computeSens");
00496 #endif
00497     
00498     // Form the residual:  df/dx * S + df/dp
00499     // F_sens = df/dx * S
00500     Thyra::apply(
00501       *DfDx_, Thyra::NOTRANS,
00502       *S, &*F_sens,
00503       ST::one(), ST::zero()
00504       );
00505     // F_sens += d(f)/d(p)
00506     Vp_V( &*F_sens, *DfDp_ );
00507   }
00508   
00509   THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_END();
00510 
00511 }
00512 
00513 
00514 // private
00515 
00516 
00517 template<class Scalar>
00518 void ForwardSensitivityExplicitModelEvaluator<Scalar>::wrapNominalValuesAndBounds()
00519 {
00520   TEUCHOS_ASSERT( !is_null(stateModel_) );
00521 
00522   using Teuchos::rcp_dynamic_cast;
00523   typedef Thyra::ModelEvaluatorBase MEB;
00524 
00525   // nominalValues_.clear(); // ToDo: Implement this!
00526 
00527   nominalValues_.set_t(stateModel_->getNominalValues().get_t());
00528   
00529   // 2007/05/22: rabartl: Note: Currently there is not much of a reason to set
00530   // an initial condition here since the initial condition for the
00531   // sensitivities is really being set in the ForwardSensitivityStepper
00532   // object!  In the future, a more general use of this class might benefit
00533   // from setting the initial condition here.
00534 
00535 }
00536 
00537 
00538 template<class Scalar>
00539 void ForwardSensitivityExplicitModelEvaluator<Scalar>::computeDerivativeMatrices(
00540   const Thyra::ModelEvaluatorBase::InArgs<Scalar> &point
00541   ) const
00542 {
00543   TEUCHOS_ASSERT( !is_null(stateModel_) );
00544 
00545   typedef Thyra::ModelEvaluatorBase MEB;
00546   typedef Teuchos::VerboseObjectTempState<MEB> VOTSME;
00547 
00548   Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
00549   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00550 
00551   MEB::InArgs<Scalar> inArgs = stateBasePoint_;
00552   MEB::OutArgs<Scalar> outArgs = stateModel_->createOutArgs();
00553   
00554   if (is_null(DfDx_)) {
00555     DfDx_ = stateModel_->create_W_op();
00556   }
00557   inArgs.set_beta(1.0);
00558   outArgs.set_W_op(DfDx_);
00559 
00560   if (is_null(DfDp_)) {
00561     DfDp_ = Thyra::create_DfDp_mv(
00562       *stateModel_,p_index_,
00563       MEB::DERIV_MV_BY_COL
00564       ).getMultiVector();
00565   }
00566   outArgs.set_DfDp(
00567     p_index_,
00568     MEB::Derivative<Scalar>(DfDp_,MEB::DERIV_MV_BY_COL)
00569     );
00570   
00571   VOTSME stateModel_outputTempState(stateModel_,out,verbLevel);
00572   stateModel_->evalModel(inArgs,outArgs);
00573   
00574 
00575 }
00576 
00577 
00578 } // namespace Rythmos
00579 
00580 
00581 #endif // RYTHMOS_FORWARD_SENSITIVITY_EXPLICIT_MODEL_EVALUATOR_HPP

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