Rythmos - Transient Integration for Differential Equations Version of the Day
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., 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_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 #include "Thyra_VectorStdOps.hpp"
00042 #include "Thyra_MultiVectorStdOps.hpp"
00043 
00044 
00045 namespace Rythmos {
00046 
00047 
00125 template<class Scalar>
00126 class ForwardSensitivityExplicitModelEvaluator
00127   : virtual public ForwardSensitivityModelEvaluatorBase<Scalar>
00128 {
00129 public:
00130 
00133 
00135   ForwardSensitivityExplicitModelEvaluator();
00136 
00138 
00141 
00143   void initializeStructure(
00144     const RCP<const Thyra::ModelEvaluator<Scalar> > &stateModel,
00145     const int p_index
00146     );
00147   
00149   void initializeStructureInitCondOnly(
00150     const RCP<const Thyra::ModelEvaluator<Scalar> >& stateModel,
00151     const RCP<const Thyra::VectorSpaceBase<Scalar> >& p_space
00152     );
00153 
00155   RCP<const Thyra::ModelEvaluator<Scalar> >
00156   getStateModel() const;
00157 
00159   RCP<Thyra::ModelEvaluator<Scalar> >
00160   getNonconstStateModel() const;
00161   
00163   int get_p_index() const;
00164   
00166   RCP<const Thyra::DefaultMultiVectorProductVectorSpace<Scalar> >
00167   get_s_bar_space() const;
00168   
00170   RCP<const Thyra::VectorSpaceBase<Scalar> > get_p_sens_space() const;
00171 
00174   void initializePointState(
00175       Ptr<StepperBase<Scalar> > stateStepper,
00176       bool forceUpToDateW
00177       );
00178 
00180 
00183 
00185   RCP<const Thyra::VectorSpaceBase<Scalar> > get_x_space() const;
00187   RCP<const Thyra::VectorSpaceBase<Scalar> > get_f_space() const;
00189   Thyra::ModelEvaluatorBase::InArgs<Scalar> getNominalValues() const;
00191   RCP<Thyra::LinearOpWithSolveBase<Scalar> > create_W() const;
00193   Thyra::ModelEvaluatorBase::InArgs<Scalar> createInArgs() const;
00194 
00196 
00197 private:
00198 
00201 
00203   Thyra::ModelEvaluatorBase::OutArgs<Scalar> createOutArgsImpl() const;
00205   void evalModelImpl(
00206     const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
00207     const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs
00208     ) const;
00209 
00211 
00212 private:
00213 
00214   // /////////////////////////
00215   // Private data members
00216 
00217   RCP<const Thyra::ModelEvaluator<Scalar> > stateModel_;
00218   int p_index_;
00219   int np_;
00220 
00221   RCP<const Thyra::DefaultMultiVectorProductVectorSpace<Scalar> > s_bar_space_;
00222   RCP<const Thyra::DefaultMultiVectorProductVectorSpace<Scalar> > f_sens_space_;
00223 
00224   Thyra::ModelEvaluatorBase::InArgs<Scalar> nominalValues_;
00225 
00226   mutable Thyra::ModelEvaluatorBase::InArgs<Scalar> stateBasePoint_;
00227 
00228   mutable RCP<Thyra::LinearOpBase<Scalar> > DfDx_;
00229   mutable RCP<Thyra::MultiVectorBase<Scalar> > DfDp_;
00230 
00231   mutable RCP<Thyra::LinearOpBase<Scalar> > DfDx_compute_;
00232   mutable RCP<Thyra::MultiVectorBase<Scalar> > DfDp_compute_;
00233 
00234   // /////////////////////////
00235   // Private member functions
00236 
00237   void wrapNominalValuesAndBounds();
00238 
00239   void computeDerivativeMatrices(
00240     const Thyra::ModelEvaluatorBase::InArgs<Scalar> &point
00241     ) const;
00242   
00243 };
00244 
00245 
00246 // /////////////////////////////////
00247 // Implementations
00248 
00249 
00250 // Constructors/Intializers/Accessors
00251 
00252 
00253 template<class Scalar>
00254 RCP<ForwardSensitivityExplicitModelEvaluator<Scalar> >
00255 forwardSensitivityExplicitModelEvaluator()
00256 {
00257   RCP<ForwardSensitivityExplicitModelEvaluator<Scalar> > fseme =
00258     rcp(new ForwardSensitivityExplicitModelEvaluator<Scalar> );
00259   return fseme;
00260 }
00261 
00262 
00263 template<class Scalar>
00264 ForwardSensitivityExplicitModelEvaluator<Scalar>::ForwardSensitivityExplicitModelEvaluator()
00265   : p_index_(0), np_(-1)
00266 {}
00267 
00268 
00269 
00270 // Public functions overridden from ForwardSensitivityModelEvaluatorBase
00271 
00272 
00273 template<class Scalar>
00274 void ForwardSensitivityExplicitModelEvaluator<Scalar>::initializeStructure(
00275   const RCP<const Thyra::ModelEvaluator<Scalar> > &stateModel,
00276   const int p_index
00277   )
00278 {
00279 
00280   typedef Thyra::ModelEvaluatorBase MEB;
00281 
00282   //
00283   // Validate input
00284   //
00285 
00286   TEUCHOS_TEST_FOR_EXCEPT( is_null(stateModel) );
00287   TEUCHOS_TEST_FOR_EXCEPTION(
00288     !( 0 <= p_index && p_index < stateModel->Np() ), std::logic_error,
00289     "Error, p_index does not fall in the range [0,"<<(stateModel->Np()-1)<<"]!" );
00290   // ToDo: Validate support for DfDp!
00291 
00292   //
00293   // Set the input objects
00294   //
00295 
00296   stateModel_ = stateModel;
00297   p_index_ = p_index;
00298   np_ = stateModel_->get_p_space(p_index)->dim();
00299 
00300   //
00301   // Create the structure of the model
00302   //
00303 
00304   s_bar_space_ = Thyra::multiVectorProductVectorSpace(
00305     stateModel_->get_x_space(), np_
00306     );
00307 
00308   f_sens_space_ = Thyra::multiVectorProductVectorSpace(
00309     stateModel_->get_f_space(), np_
00310     );
00311 
00312   nominalValues_ = this->createInArgs();
00313 
00314   this->wrapNominalValuesAndBounds();
00315 
00316   //
00317   // Wipe out matrix storage
00318   //
00319 
00320   stateBasePoint_ = MEB::InArgs<Scalar>();
00321   DfDx_ = Teuchos::null;
00322   DfDp_ = Teuchos::null;
00323   DfDx_compute_ = Teuchos::null;
00324   DfDp_compute_ = Teuchos::null;
00325 
00326 }
00327 
00328 
00329 template<class Scalar>
00330 void ForwardSensitivityExplicitModelEvaluator<Scalar>::initializeStructureInitCondOnly(
00331   const RCP<const Thyra::ModelEvaluator<Scalar> >& stateModel,
00332   const RCP<const Thyra::VectorSpaceBase<Scalar> >& p_space
00333   )
00334 {
00335   TEUCHOS_TEST_FOR_EXCEPT_MSG(true, "ToDo: Implement initializeStructureInitCondOnly()!" );
00336 }
00337 
00338 
00339 template<class Scalar>
00340 RCP<const Thyra::ModelEvaluator<Scalar> >
00341 ForwardSensitivityExplicitModelEvaluator<Scalar>::getStateModel() const
00342 {
00343   return stateModel_;
00344 }
00345 
00346 
00347 template<class Scalar>
00348 RCP<Thyra::ModelEvaluator<Scalar> >
00349 ForwardSensitivityExplicitModelEvaluator<Scalar>::getNonconstStateModel() const
00350 {
00351   return Teuchos::null;
00352 }
00353 
00354 
00355 template<class Scalar>
00356 int ForwardSensitivityExplicitModelEvaluator<Scalar>::get_p_index() const
00357 {
00358   return p_index_;
00359 }
00360 
00361 
00362 template<class Scalar>
00363 RCP<const Thyra::DefaultMultiVectorProductVectorSpace<Scalar> >
00364 ForwardSensitivityExplicitModelEvaluator<Scalar>::get_s_bar_space() const
00365 {
00366   return s_bar_space_;
00367 }
00368 
00369 
00370 template<class Scalar>
00371 RCP<const Thyra::VectorSpaceBase<Scalar> >
00372 ForwardSensitivityExplicitModelEvaluator<Scalar>::get_p_sens_space() const
00373 {
00374   return stateModel_->get_p_space(p_index_);
00375 }
00376 
00377 
00378 template<class Scalar>
00379 void ForwardSensitivityExplicitModelEvaluator<Scalar>::initializePointState(
00380     Ptr<StepperBase<Scalar> > stateStepper,
00381     bool forceUpToDateW
00382     )
00383 {
00384   TEUCHOS_ASSERT( Teuchos::nonnull(stateStepper) );
00385 #ifdef HAVE_RYTHMOS_DEBUG
00386   TEUCHOS_TEST_FOR_EXCEPTION(
00387     is_null(stateModel_), std::logic_error,
00388     "Error, you must call intializeStructure(...) before you call initializePointState(...)"
00389     );
00390 #endif // HAVE_RYTHMOS_DEBUG
00391 
00392   Scalar curr_t = stateStepper->getStepStatus().time;
00393   RCP<const Thyra::VectorBase<Scalar> > x;
00394   x = get_x(*stateStepper,curr_t);
00395 #ifdef HAVE_RYTHMOS_DEBUG
00396   TEUCHOS_TEST_FOR_EXCEPT( Teuchos::is_null(x) );
00397 #endif // HAVE_RYTHMOS_DEBUG
00398       
00399   stateBasePoint_ = stateStepper->getInitialCondition(); // set parameters
00400   stateBasePoint_.set_x( x );
00401   stateBasePoint_.set_t( curr_t );
00402 
00403   // Set whatever derivatives where passed in.  If an input in null, then the
00404   // member will be null and the null linear operators will be computed later
00405   // just in time.
00406 
00407   wrapNominalValuesAndBounds();
00408   
00409 }
00410 
00411 
00412 // Public functions overridden from ModelEvaulator
00413 
00414 
00415 template<class Scalar>
00416 RCP<const Thyra::VectorSpaceBase<Scalar> >
00417 ForwardSensitivityExplicitModelEvaluator<Scalar>::get_x_space() const
00418 {
00419   return s_bar_space_;
00420 }
00421 
00422 
00423 template<class Scalar>
00424 RCP<const Thyra::VectorSpaceBase<Scalar> >
00425 ForwardSensitivityExplicitModelEvaluator<Scalar>::get_f_space() const
00426 {
00427   return f_sens_space_;
00428 }
00429 
00430 
00431 template<class Scalar>
00432 Thyra::ModelEvaluatorBase::InArgs<Scalar>
00433 ForwardSensitivityExplicitModelEvaluator<Scalar>::getNominalValues() const
00434 {
00435   return nominalValues_;
00436 }
00437 
00438 
00439 template<class Scalar>
00440 RCP<Thyra::LinearOpWithSolveBase<Scalar> >
00441 ForwardSensitivityExplicitModelEvaluator<Scalar>::create_W() const
00442 {
00443   return Thyra::multiVectorLinearOpWithSolve<Scalar>();
00444 }
00445 
00446 
00447 template<class Scalar>
00448 Thyra::ModelEvaluatorBase::InArgs<Scalar>
00449 ForwardSensitivityExplicitModelEvaluator<Scalar>::createInArgs() const
00450 {
00451   TEUCHOS_ASSERT( !is_null(stateModel_) );
00452   typedef Thyra::ModelEvaluatorBase MEB;
00453   MEB::InArgs<Scalar> stateModelInArgs = stateModel_->createInArgs();
00454   MEB::InArgsSetup<Scalar> inArgs;
00455   inArgs.setModelEvalDescription(this->description());
00456   inArgs.setSupports( MEB::IN_ARG_x );
00457   inArgs.setSupports( MEB::IN_ARG_t );
00458   inArgs.setSupports( MEB::IN_ARG_beta,
00459     stateModelInArgs.supports(MEB::IN_ARG_beta) );
00460   return inArgs;
00461 }
00462 
00463 
00464 // Private functions overridden from ModelEvaulatorDefaultBase
00465 
00466 
00467 template<class Scalar>
00468 Thyra::ModelEvaluatorBase::OutArgs<Scalar>
00469 ForwardSensitivityExplicitModelEvaluator<Scalar>::createOutArgsImpl() const
00470 {
00471   TEUCHOS_ASSERT( !is_null(stateModel_) );
00472   typedef Thyra::ModelEvaluatorBase MEB;
00473 
00474   MEB::OutArgs<Scalar> stateModelOutArgs = stateModel_->createOutArgs();
00475   MEB::OutArgsSetup<Scalar> outArgs;
00476 
00477   outArgs.setModelEvalDescription(this->description());
00478 
00479   outArgs.setSupports(MEB::OUT_ARG_f);
00480 
00481   return outArgs;
00482 
00483 }
00484 
00485 
00486 template<class Scalar>
00487 void ForwardSensitivityExplicitModelEvaluator<Scalar>::evalModelImpl(
00488   const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
00489   const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs
00490   ) const
00491 {
00492 
00493   using Teuchos::rcp_dynamic_cast;
00494   typedef Teuchos::ScalarTraits<Scalar> ST;
00495   typedef Thyra::ModelEvaluatorBase MEB;
00496   typedef Teuchos::VerboseObjectTempState<Thyra::ModelEvaluatorBase> VOTSME;
00497 
00498   THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_GEN_BEGIN(
00499     "ForwardSensitivityExplicitModelEvaluator", inArgs, outArgs, Teuchos::null );
00500 
00501   //
00502   // Update the derivative matrices if they are not already updated for the
00503   // given time!.
00504   //
00505   
00506   {
00507     RYTHMOS_FUNC_TIME_MONITOR_DIFF(
00508         "Rythmos:ForwardSensitivityExplicitModelEvaluator::evalModel: computeMatrices",
00509         Rythmos_FSEME
00510         );
00511     computeDerivativeMatrices(inArgs);
00512   }
00513 
00514   //
00515   // InArgs
00516   //
00517 
00518   RCP<const Thyra::DefaultMultiVectorProductVector<Scalar> >
00519     s_bar = rcp_dynamic_cast<const Thyra::DefaultMultiVectorProductVector<Scalar> >(
00520       inArgs.get_x().assert_not_null(), true
00521       );
00522   RCP<const Thyra::MultiVectorBase<Scalar> >
00523     S = s_bar->getMultiVector();
00524   
00525   //
00526   // OutArgs
00527   //
00528 
00529   RCP<Thyra::DefaultMultiVectorProductVector<Scalar> >
00530     f_sens = rcp_dynamic_cast<Thyra::DefaultMultiVectorProductVector<Scalar> >(
00531       outArgs.get_f(), true
00532       );
00533 
00534   RCP<Thyra::MultiVectorBase<Scalar> >
00535     F_sens = f_sens->getNonconstMultiVector().assert_not_null();
00536 
00537   //
00538   // Compute the requested functions
00539   //
00540 
00541   if(!is_null(F_sens)) {
00542 
00543     RYTHMOS_FUNC_TIME_MONITOR_DIFF(
00544         "Rythmos:ForwardSensitivityExplicitModelEvaluator::evalModel: computeSens",
00545         Rythmos_FSEME
00546         );
00547     
00548     // Form the residual:  df/dx * S + df/dp
00549     // F_sens = df/dx * S
00550     Thyra::apply(
00551       *DfDx_, Thyra::NOTRANS,
00552       *S, F_sens.ptr(),
00553       ST::one(), ST::zero()
00554       );
00555     // F_sens += d(f)/d(p)
00556     Vp_V( F_sens.ptr(), *DfDp_ );
00557   }
00558   
00559   THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_END();
00560 
00561 }
00562 
00563 
00564 // private
00565 
00566 
00567 template<class Scalar>
00568 void ForwardSensitivityExplicitModelEvaluator<Scalar>::wrapNominalValuesAndBounds()
00569 {
00570   TEUCHOS_ASSERT( !is_null(stateModel_) );
00571 
00572   using Teuchos::rcp_dynamic_cast;
00573   typedef Thyra::ModelEvaluatorBase MEB;
00574   typedef Teuchos::ScalarTraits<Scalar> ST;
00575 
00576   // nominalValues_.clear(); // ToDo: Implement this!
00577 
00578   nominalValues_.set_t(stateModel_->getNominalValues().get_t());
00579   
00580   // 2007/05/22: rabartl: Note: Currently there is not much of a reason to set
00581   // an initial condition here since the initial condition for the
00582   // sensitivities is really being set in the ForwardSensitivityStepper
00583   // object!  In the future, a more general use of this class might benefit
00584   // from setting the initial condition here.
00585 
00586   // 2009/07/20: tscoffe/ccober:  This is the future.  We're going to use this
00587   // in a more general way, so we need legitimate nominal values now.
00588   RCP<VectorBase<Scalar> > s_bar_ic = Thyra::createMember(this->get_x_space());
00589   Thyra::V_S(s_bar_ic.ptr(),ST::zero());
00590   nominalValues_.set_x(s_bar_ic);
00591 }
00592 
00593 
00594 template<class Scalar>
00595 void ForwardSensitivityExplicitModelEvaluator<Scalar>::computeDerivativeMatrices(
00596   const Thyra::ModelEvaluatorBase::InArgs<Scalar> &point
00597   ) const
00598 {
00599   TEUCHOS_ASSERT( !is_null(stateModel_) );
00600 
00601   typedef Thyra::ModelEvaluatorBase MEB;
00602   typedef Teuchos::VerboseObjectTempState<MEB> VOTSME;
00603 
00604   Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
00605   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00606 
00607   MEB::InArgs<Scalar> inArgs = stateBasePoint_;
00608   MEB::OutArgs<Scalar> outArgs = stateModel_->createOutArgs();
00609   
00610   if (is_null(DfDx_)) {
00611     DfDx_ = stateModel_->create_W_op();
00612   }
00613   if (inArgs.supports(MEB::IN_ARG_beta)) {
00614     inArgs.set_beta(1.0);
00615   }
00616   outArgs.set_W_op(DfDx_);
00617 
00618   if (is_null(DfDp_)) {
00619     DfDp_ = Thyra::create_DfDp_mv(
00620       *stateModel_,p_index_,
00621       MEB::DERIV_MV_BY_COL
00622       ).getMultiVector();
00623   }
00624   outArgs.set_DfDp(
00625     p_index_,
00626     MEB::Derivative<Scalar>(DfDp_,MEB::DERIV_MV_BY_COL)
00627     );
00628   
00629   VOTSME stateModel_outputTempState(stateModel_,out,verbLevel);
00630   stateModel_->evalModel(inArgs,outArgs);
00631   
00632 
00633 }
00634 
00635 
00636 } // namespace Rythmos
00637 
00638 
00639 #endif // RYTHMOS_FORWARD_SENSITIVITY_EXPLICIT_MODEL_EVALUATOR_HPP
 All Classes Functions Variables Typedefs Friends