Rythmos_ForwardSensitivityModelEvaluator.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_MODEL_EVALUATOR_HPP
00030 #define RYTHMOS_FORWARD_SENSITIVITY_MODEL_EVALUATOR_HPP
00031 
00032 
00033 #include "Thyra_ModelEvaluator.hpp" // Interface
00034 #include "Thyra_StateFuncModelEvaluatorBase.hpp" // Implementation
00035 #include "Thyra_DefaultProductVectorSpace.hpp"
00036 #include "Thyra_PhysicallyBlockedLinearOpWithSolveBase.hpp" // Interface
00037 #include "Thyra_DefaultBlockedTriangularLinearOpWithSolve.hpp" // Implementation
00038 #include "Thyra_ModelEvaluatorDelegatorBase.hpp"
00039 #include "Thyra_ModelEvaluatorHelpers.hpp"
00040 #include "Thyra_DefaultMultiVectorProductVectorSpace.hpp"
00041 #include "Thyra_DefaultMultiVectorProductVector.hpp"
00042 #include "Thyra_DefaultMultiVectorLinearOpWithSolve.hpp"
00043 #include "Teuchos_implicit_cast.hpp"
00044 #include "Teuchos_Assert.hpp"
00045 
00046 
00047 namespace Rythmos {
00048 
00049 
00296 template<class Scalar>
00297 class ForwardSensitivityModelEvaluator
00298   : virtual public Thyra::StateFuncModelEvaluatorBase<Scalar>
00299 {
00300 public:
00301 
00304 
00306   ForwardSensitivityModelEvaluator();
00307 
00323   void initializeStructure(
00324     const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > &stateModel,
00325     const int p_index
00326     );
00327   
00329   Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >
00330   getStateModel() const;
00331   
00333   int get_p_index() const;
00334 
00375   void initializeState(
00376     const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint,
00377     const Teuchos::RCP<const Thyra::LinearOpWithSolveBase<Scalar> > &W_tilde,
00378     const Scalar &coeff_x_dot,
00379     const Scalar &coeff_x,
00380     const Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > &DfDx_dot = Teuchos::null,
00381     const Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> > &DfDp = Teuchos::null
00382     );
00383 
00384   // 2007/05/22: rabartl: ToDo: Add an InterpolationBufferBase
00385   // stateInterpBuffer object to the initailizeState(...) function that can be
00386   // used to get x and x_dot at different points in time t.  Then, modify the
00387   // logic to recompute all of the needed matrices if t != t_base (as passed
00388   // in through stateBasePoint).  The values of x(t) and xdot(t) can then be
00389   // gotten from the stateInterpBuffer object!
00390   
00392 
00395 
00397   Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > get_x_space() const;
00399   Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > get_f_space() const;
00401   Thyra::ModelEvaluatorBase::InArgs<Scalar> getNominalValues() const;
00403   Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> > create_W() const;
00405   Thyra::ModelEvaluatorBase::InArgs<Scalar> createInArgs() const;
00406 
00408 
00409 private:
00410 
00413 
00415   Thyra::ModelEvaluatorBase::OutArgs<Scalar> createOutArgsImpl() const;
00417   void evalModelImpl(
00418     const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
00419     const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs
00420     ) const;
00421 
00423 
00424 private:
00425 
00426   // /////////////////////////
00427   // Private data members
00428 
00429   Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > stateModel_;
00430   int p_index_;
00431   int np_;
00432   Teuchos::RCP<const Thyra::DefaultMultiVectorProductVectorSpace<Scalar> > s_bar_space_;
00433   Teuchos::RCP<const Thyra::DefaultMultiVectorProductVectorSpace<Scalar> > f_sens_space_;
00434   Thyra::ModelEvaluatorBase::InArgs<Scalar> nominalValues_;
00435 
00436   Thyra::ModelEvaluatorBase::InArgs<Scalar> stateBasePoint_;
00437 
00438   mutable Teuchos::RCP<const Thyra::LinearOpWithSolveBase<Scalar> > W_tilde_;
00439   mutable Scalar coeff_x_dot_;
00440   mutable Scalar coeff_x_;
00441   mutable Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > DfDx_dot_;
00442   mutable Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> > DfDp_;
00443 
00444 
00445   // /////////////////////////
00446   // Private member functions
00447 
00448   void wrapNominalValuesAndBounds();
00449 
00450   void computeDerivativeMatrices(
00451     const Thyra::ModelEvaluatorBase::InArgs<Scalar> &point
00452     ) const;
00453   
00454 };
00455 
00456 
00457 // /////////////////////////////////
00458 // Implementations
00459 
00460 
00461 // Constructors/Intializers/Accessors
00462 
00463 
00464 template<class Scalar>
00465 ForwardSensitivityModelEvaluator<Scalar>::ForwardSensitivityModelEvaluator()
00466   : p_index_(0), np_(-1)
00467 {}
00468 
00469 
00470 template<class Scalar>
00471 void ForwardSensitivityModelEvaluator<Scalar>::initializeStructure(
00472   const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > &stateModel,
00473   const int p_index
00474   )
00475 {
00476 
00477   typedef Thyra::ModelEvaluatorBase MEB;
00478 
00479   //
00480   // Validate input
00481   //
00482 
00483   TEST_FOR_EXCEPT( is_null(stateModel) );
00484   TEST_FOR_EXCEPTION(
00485     !( 0 <= p_index && p_index < stateModel->Np() ), std::logic_error,
00486     "Error, p_index does not fall in the range [0,"<<(stateModel->Np()-1)<<"]!" );
00487   // ToDo: Validate support for DfDp!
00488 
00489   //
00490   // Set the input objects
00491   //
00492 
00493   stateModel_ = stateModel;
00494   p_index_ = p_index;
00495   np_ = stateModel_->get_p_space(p_index)->dim();
00496 
00497   //
00498   // Create the structure of the model
00499   //
00500 
00501   s_bar_space_ = Thyra::multiVectorProductVectorSpace(
00502     stateModel_->get_x_space(), np_
00503     );
00504 
00505   f_sens_space_ = Thyra::multiVectorProductVectorSpace(
00506     stateModel_->get_f_space(), np_
00507     );
00508 
00509   nominalValues_ = this->createInArgs();
00510 
00511 }
00512 
00513 
00514 template<class Scalar>
00515 Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >
00516 ForwardSensitivityModelEvaluator<Scalar>::getStateModel() const
00517 {
00518   return stateModel_;
00519 }
00520 
00521 
00522 template<class Scalar>
00523 int ForwardSensitivityModelEvaluator<Scalar>::get_p_index() const
00524 {
00525   return p_index_;
00526 }
00527 
00528 
00529 template<class Scalar>
00530 void ForwardSensitivityModelEvaluator<Scalar>::initializeState(
00531   const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint,
00532   const Teuchos::RCP<const Thyra::LinearOpWithSolveBase<Scalar> > &W_tilde,
00533   const Scalar &coeff_x_dot,
00534   const Scalar &coeff_x,
00535   const Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > &DfDx_dot,
00536   const Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> > &DfDp
00537   )
00538 {
00539 
00540   typedef Thyra::ModelEvaluatorBase MEB;
00541 
00542 #ifdef TEUCHOS_DEBUG
00543   TEST_FOR_EXCEPTION(
00544     is_null(stateModel_), std::logic_error,
00545     "Error, you must call intializeStructure(...) before you call initializeState(...)"
00546     );
00547   TEST_FOR_EXCEPT( is_null(stateBasePoint.get_x()) );
00548   TEST_FOR_EXCEPT( is_null(stateBasePoint.get_x_dot()) );
00549   TEST_FOR_EXCEPT( is_null(stateBasePoint.get_p(p_index_)) );
00550   // What about the other parameter values?  We really can't say anything
00551   // about these and we can't check them.  They can be null just fine.
00552   if (!is_null(W_tilde)) {
00553     THYRA_ASSERT_VEC_SPACES("",*W_tilde->domain(),*stateModel_->get_x_space());
00554     THYRA_ASSERT_VEC_SPACES("",*W_tilde->range(),*stateModel_->get_f_space());
00555   }
00556   if (!is_null(DfDx_dot)) {
00557     THYRA_ASSERT_VEC_SPACES("",*DfDx_dot->domain(),*stateModel_->get_x_space());
00558     THYRA_ASSERT_VEC_SPACES("",*DfDx_dot->range(),*stateModel_->get_f_space());
00559   }
00560   if (!is_null(DfDp)) {
00561     THYRA_ASSERT_VEC_SPACES("",*DfDp->domain(),*stateModel_->get_p_space(p_index_));
00562     THYRA_ASSERT_VEC_SPACES("",*DfDp->range(),*stateModel_->get_f_space());
00563   }
00564 #endif
00565 
00566   stateBasePoint_ = stateBasePoint;
00567 
00568   // Set whatever derivatives where passed in.  If an input in null, then the
00569   // member will be null and the null linear operators will be computed later
00570   // just in time.
00571 
00572   W_tilde_ = W_tilde;
00573   coeff_x_dot_ = coeff_x_dot;
00574   coeff_x_ = coeff_x;
00575   DfDx_dot_ = DfDx_dot;
00576   DfDp_ = DfDp;
00577 
00578   wrapNominalValuesAndBounds();
00579 
00580 }
00581 
00582 
00583 // Public functions overridden from ModelEvaulator
00584 
00585 
00586 template<class Scalar>
00587 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
00588 ForwardSensitivityModelEvaluator<Scalar>::get_x_space() const
00589 {
00590   return s_bar_space_;
00591 }
00592 
00593 
00594 template<class Scalar>
00595 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
00596 ForwardSensitivityModelEvaluator<Scalar>::get_f_space() const
00597 {
00598   return f_sens_space_;
00599 }
00600 
00601 
00602 template<class Scalar>
00603 Thyra::ModelEvaluatorBase::InArgs<Scalar>
00604 ForwardSensitivityModelEvaluator<Scalar>::getNominalValues() const
00605 {
00606   return nominalValues_;
00607 }
00608 
00609 
00610 template<class Scalar>
00611 Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> >
00612 ForwardSensitivityModelEvaluator<Scalar>::create_W() const
00613 {
00614   return Thyra::multiVectorLinearOpWithSolve<Scalar>();
00615 }
00616 
00617 
00618 template<class Scalar>
00619 Thyra::ModelEvaluatorBase::InArgs<Scalar>
00620 ForwardSensitivityModelEvaluator<Scalar>::createInArgs() const
00621 {
00622   typedef Thyra::ModelEvaluatorBase MEB;
00623   MEB::InArgs<Scalar> stateModelInArgs = stateModel_->createInArgs();
00624   MEB::InArgsSetup<Scalar> inArgs;
00625   inArgs.setModelEvalDescription(this->description());
00626   inArgs.setSupports( MEB::IN_ARG_x_dot,
00627     stateModelInArgs.supports(MEB::IN_ARG_x_dot) );
00628   inArgs.setSupports( MEB::IN_ARG_x );
00629   inArgs.setSupports( MEB::IN_ARG_t );
00630   inArgs.setSupports( MEB::IN_ARG_alpha,
00631     stateModelInArgs.supports(MEB::IN_ARG_alpha) );
00632   inArgs.setSupports( MEB::IN_ARG_beta,
00633     stateModelInArgs.supports(MEB::IN_ARG_beta) );
00634   return inArgs;
00635 }
00636 
00637 
00638 // Private functions overridden from ModelEvaulatorDefaultBase
00639 
00640 
00641 template<class Scalar>
00642 Thyra::ModelEvaluatorBase::OutArgs<Scalar>
00643 ForwardSensitivityModelEvaluator<Scalar>::createOutArgsImpl() const
00644 {
00645 
00646   typedef Thyra::ModelEvaluatorBase MEB;
00647 
00648   MEB::OutArgs<Scalar> stateModelOutArgs = stateModel_->createOutArgs();
00649   MEB::OutArgsSetup<Scalar> outArgs;
00650 
00651   outArgs.setModelEvalDescription(this->description());
00652 
00653   outArgs.setSupports(MEB::OUT_ARG_f);
00654 
00655   if (stateModelOutArgs.supports(MEB::OUT_ARG_W) ) {
00656     outArgs.setSupports(MEB::OUT_ARG_W);
00657     outArgs.set_W_properties(stateModelOutArgs.get_W_properties());
00658   }
00659 
00660   return outArgs;
00661 
00662 }
00663 
00664 
00665 template<class Scalar>
00666 void ForwardSensitivityModelEvaluator<Scalar>::evalModelImpl(
00667   const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
00668   const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs
00669   ) const
00670 {
00671 
00672   using Teuchos::rcp_dynamic_cast;
00673   typedef Teuchos::ScalarTraits<Scalar> ST;
00674   typedef Thyra::ModelEvaluatorBase MEB;
00675   typedef Teuchos::VerboseObjectTempState<Thyra::ModelEvaluatorBase> VOTSME;
00676 
00677   THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_GEN_BEGIN(
00678     "ForwardSensitivityModelEvaluator", inArgs, outArgs, Teuchos::null );
00679 
00680   //
00681   // Update the derivative matrices if they are not already updated for the
00682   // given time!.
00683   //
00684   
00685   {
00686 #ifdef ENABLE_RYTHMOS_TIMERS
00687     TEUCHOS_FUNC_TIME_MONITOR("Rythmos:ForwardSensitivityModelEvaluator::evalModel: computeMatrices");
00688 #endif
00689     computeDerivativeMatrices(inArgs);
00690   }
00691 
00692   //
00693   // InArgs
00694   //
00695 
00696   RCP<const Thyra::DefaultMultiVectorProductVector<Scalar> >
00697     s_bar = rcp_dynamic_cast<const Thyra::DefaultMultiVectorProductVector<Scalar> >(
00698       inArgs.get_x().assert_not_null(), true
00699       );
00700   RCP<const Thyra::DefaultMultiVectorProductVector<Scalar> >
00701     s_bar_dot = rcp_dynamic_cast<const Thyra::DefaultMultiVectorProductVector<Scalar> >(
00702       inArgs.get_x_dot().assert_not_null(), true
00703       );
00704   const Scalar
00705     alpha = inArgs.get_alpha();
00706   const Scalar
00707     beta = inArgs.get_beta();
00708 
00709   RCP<const Thyra::MultiVectorBase<Scalar> >
00710     S = s_bar->getMultiVector();
00711   RCP<const Thyra::MultiVectorBase<Scalar> >
00712     S_dot = s_bar_dot->getMultiVector();
00713   
00714   //
00715   // OutArgs
00716   //
00717 
00718   RCP<Thyra::DefaultMultiVectorProductVector<Scalar> >
00719     f_sens = rcp_dynamic_cast<Thyra::DefaultMultiVectorProductVector<Scalar> >(
00720       outArgs.get_f(), true
00721       );
00722   RCP<Thyra::DefaultMultiVectorLinearOpWithSolve<Scalar> >
00723     W_sens = rcp_dynamic_cast<Thyra::DefaultMultiVectorLinearOpWithSolve<Scalar> >(
00724       outArgs.get_W(), true
00725       );
00726 
00727   RCP<Thyra::MultiVectorBase<Scalar> >
00728     F_sens = f_sens->getNonconstMultiVector().assert_not_null();
00729 
00730   //
00731   // Compute the requested functions
00732   //
00733 
00734   if(!is_null(F_sens)) {
00735 
00736 #ifdef ENABLE_RYTHMOS_TIMERS
00737     TEUCHOS_FUNC_TIME_MONITOR("Rythmos:ForwardSensitivityModelEvaluator::evalModel: computeSens");
00738 #endif
00739 
00740     // S_diff =  -(coeff_x_dot/coeff_x)*S + S_dot
00741     Teuchos::RCP<Thyra::MultiVectorBase<Scalar> >
00742       S_diff = createMembers( stateModel_->get_x_space(), np_ );
00743     V_StVpV( &*S_diff, Scalar(-coeff_x_dot_/coeff_x_), *S, *S_dot );
00744     // F_sens = (1/coeff_x) * W_tilde * S
00745     Thyra::apply(
00746       *W_tilde_, Thyra::NOTRANS,
00747       *S, &*F_sens,
00748       Scalar(1.0/coeff_x_), ST::zero()
00749       );
00750     // F_sens += d(f)/d(x_dot) * S_diff
00751     Thyra::apply(
00752       *DfDx_dot_, Thyra::NOTRANS,
00753       *S_diff, &*F_sens,
00754       ST::one(), ST::one()
00755       );
00756     // F_sens += d(f)/d(p)
00757     Vp_V( &*F_sens, *DfDp_ );
00758   }
00759   
00760   if(!is_null(W_sens)) {
00761     TEST_FOR_EXCEPTION(
00762       alpha != coeff_x_dot_, std::logic_error,
00763       "Error, alpha="<<alpha<<" != coeff_x_dot="<<coeff_x_dot_
00764       <<" with difference = "<<(alpha-coeff_x_dot_)<<"!" );
00765     TEST_FOR_EXCEPTION(
00766       beta != coeff_x_, std::logic_error,
00767       "Error, beta="<<beta<<" != coeff_x="<<coeff_x_
00768       <<" with difference = "<<(beta-coeff_x_)<<"!" );
00769     W_sens->initialize( W_tilde_, s_bar_space_, f_sens_space_ );
00770     
00771   }
00772 
00773   THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_END();
00774 
00775 }
00776 
00777 
00778 // private
00779 
00780 
00781 template<class Scalar>
00782 void ForwardSensitivityModelEvaluator<Scalar>::wrapNominalValuesAndBounds()
00783 {
00784 
00785   using Teuchos::rcp_dynamic_cast;
00786   typedef Thyra::ModelEvaluatorBase MEB;
00787 
00788   // nominalValues_.clear(); // ToDo: Implement this!
00789 
00790   nominalValues_.set_t(stateModel_->getNominalValues().get_t());
00791   
00792   // 2007/05/22: rabartl: Note: Currently there is not much of a reason to set
00793   // an initial condition here since the initial condition for the
00794   // sensitivities is really being set in the ForwardSensitivityStepper
00795   // object!  In the future, a more general use of this class might benefit
00796   // from setting the initial condition here.
00797 
00798 }
00799 
00800 
00801 template<class Scalar>
00802 void ForwardSensitivityModelEvaluator<Scalar>::computeDerivativeMatrices(
00803   const Thyra::ModelEvaluatorBase::InArgs<Scalar> &point
00804   ) const
00805 {
00806 
00807   typedef Thyra::ModelEvaluatorBase MEB;
00808   typedef Teuchos::VerboseObjectTempState<MEB> VOTSME;
00809 
00810   Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
00811   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00812 
00813   const Scalar
00814     t_base = stateBasePoint_.get_t(),
00815     t = point.get_t();
00816 
00817   TEUCHOS_ASSERT_EQUALITY( t , t_base );
00818 
00819   if (is_null(W_tilde_)) {
00820     TEST_FOR_EXCEPT("ToDo: compute W_tilde from scratch!");
00821   }
00822   
00823   if ( is_null(DfDx_dot_) || is_null(DfDp_) ) {
00824 
00825     MEB::InArgs<Scalar> inArgs = stateBasePoint_;
00826     MEB::OutArgs<Scalar> outArgs = stateModel_->createOutArgs();
00827     
00828     Teuchos::RCP<Thyra::LinearOpBase<Scalar> > DfDx_dot_compute;
00829     if (is_null(DfDx_dot_)) {
00830       DfDx_dot_compute = stateModel_->create_W_op();
00831       inArgs.set_alpha(1.0);
00832       inArgs.set_beta(0.0);
00833       outArgs.set_W_op(DfDx_dot_compute);
00834     }
00835 
00836     Teuchos::RCP<Thyra::MultiVectorBase<Scalar> > DfDp_compute;
00837     if (is_null(DfDp_)) {
00838       DfDp_compute = Thyra::create_DfDp_mv(
00839         *stateModel_,p_index_,
00840         MEB::DERIV_MV_BY_COL
00841         ).getMultiVector();
00842       outArgs.set_DfDp(
00843         p_index_,
00844         MEB::Derivative<Scalar>(DfDp_compute,MEB::DERIV_MV_BY_COL)
00845         );
00846     }
00847     
00848     VOTSME stateModel_outputTempState(stateModel_,out,verbLevel);
00849     stateModel_->evalModel(inArgs,outArgs);
00850     if (!is_null(DfDx_dot_compute))
00851       DfDx_dot_ = DfDx_dot_compute;
00852     if (!is_null(DfDp_compute))
00853       DfDp_ = DfDp_compute;
00854   
00855   }
00856 
00857 }
00858 
00859 
00860 } // namespace Rythmos
00861 
00862 
00863 #endif // RYTHMOS_FORWARD_SENSITIVITY_MODEL_EVALUATOR_HPP

Generated on Tue Oct 20 12:46:07 2009 for Rythmos - Transient Integration for Differential Equations by doxygen 1.4.7