Rythmos - Transient Integration for Differential Equations Version of the Day
Rythmos_ForwardResponseSensitivityComputer.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_RESPONSE_SENSITIVITY_COMPUTER_HPP
00030 #define RYTHMOS_FORWARD_RESPONSE_SENSITIVITY_COMPUTER_HPP
00031 
00032 
00033 #include "Rythmos_Types.hpp"
00034 #include "Thyra_ModelEvaluator.hpp"
00035 #include "Teuchos_VerboseObject.hpp"
00036 #include "Teuchos_StandardMemberCompositionMacros.hpp"
00037 
00038 
00039 namespace Rythmos {
00040 
00041 
00048 template<class Scalar>
00049 class ForwardResponseSensitivityComputer
00050   : public Teuchos::VerboseObject<ForwardResponseSensitivityComputer<Scalar> >
00051 {
00052 public:
00053   
00055   ForwardResponseSensitivityComputer();
00056 
00058   STANDARD_MEMBER_COMPOSITION_MEMBERS( bool, dumpSensitivities );
00059 
00076   void setResponseFunction(
00077     const RCP<const Thyra::ModelEvaluator<Scalar> > &responseFunc,
00078     const Thyra::ModelEvaluatorBase::InArgs<Scalar> &basePoint,
00079     const int p_index,
00080     const int g_index
00081     );
00082 
00087   void resetResponseFunction(
00088     const RCP<const Thyra::ModelEvaluator<Scalar> > &responseFunc,
00089     const Thyra::ModelEvaluatorBase::InArgs<Scalar> &basePoint
00090     );
00091 
00093   const RCP<Thyra::VectorBase<Scalar> > create_g_hat() const;
00094 
00096   const RCP<Thyra::MultiVectorBase<Scalar> > create_D_g_hat_D_p() const;
00097 
00110   void computeResponse(
00111     const Thyra::VectorBase<Scalar> *x_dot,
00112     const Thyra::VectorBase<Scalar> &x,
00113     const Scalar t,
00114     Thyra::VectorBase<Scalar> *g_hat
00115     ) const;
00116 
00136   void computeResponseAndSensitivity(
00137     const Thyra::VectorBase<Scalar> *x_dot,
00138     const Thyra::MultiVectorBase<Scalar> *S_dot,
00139     const Thyra::VectorBase<Scalar> &x,
00140     const Thyra::MultiVectorBase<Scalar> &S,
00141     const Scalar t,
00142     Thyra::VectorBase<Scalar> *g_hat,
00143     Thyra::MultiVectorBase<Scalar> *D_g_hat_D_p
00144     ) const;
00145 
00146 private: // Data members
00147 
00148   RCP<const Thyra::ModelEvaluator<Scalar> > responseFunc_;
00149   Thyra::ModelEvaluatorBase::InArgs<Scalar> basePoint_;
00150   int p_index_;
00151   int g_index_;
00152 
00153   RCP<const Thyra::VectorSpaceBase<Scalar> > p_space_;
00154   RCP<const Thyra::VectorSpaceBase<Scalar> > g_space_;
00155 
00156   bool response_func_supports_x_dot_;
00157   bool response_func_supports_D_x_dot_;
00158   bool response_func_supports_D_p_;
00159 
00160   mutable Thyra::ModelEvaluatorBase::InArgs<Scalar> responseInArgs_;
00161   mutable Thyra::ModelEvaluatorBase::OutArgs<Scalar> responseOutArgs_;
00162 
00163   mutable RCP<Thyra::LinearOpBase<Scalar> > D_g_D_x_dot_;
00164   mutable RCP<Thyra::LinearOpBase<Scalar> > D_g_D_x_;
00165   mutable RCP<Thyra::MultiVectorBase<Scalar> > D_g_D_p_;
00166 
00167 private: // Functions
00168 
00169   void clearCache();
00170 
00171   void createCache(const bool computeSens) const;
00172 
00173   void computeResponseAndSensitivityImpl(
00174     const Thyra::VectorBase<Scalar> *x_dot,
00175     const Thyra::MultiVectorBase<Scalar> *S_dot,
00176     const Thyra::VectorBase<Scalar> &x,
00177     const Thyra::MultiVectorBase<Scalar> *S,
00178     const Scalar t,
00179     Thyra::VectorBase<Scalar> *g_hat,
00180     Thyra::MultiVectorBase<Scalar> *D_g_hat_D_p
00181     ) const;
00182 
00183 };
00184 
00185 
00186 //
00187 // Implementations
00188 //
00189 
00190 
00191 template<class Scalar>
00192 ForwardResponseSensitivityComputer<Scalar>::ForwardResponseSensitivityComputer()
00193   :dumpSensitivities_(false),
00194    p_index_(-1),
00195    g_index_(-1),
00196    response_func_supports_x_dot_(false),
00197    response_func_supports_D_x_dot_(false),
00198    response_func_supports_D_p_(false)
00199 {}
00200 
00201 
00202 template<class Scalar>
00203 void ForwardResponseSensitivityComputer<Scalar>::setResponseFunction(
00204   const RCP<const Thyra::ModelEvaluator<Scalar> > &responseFunc,
00205   const Thyra::ModelEvaluatorBase::InArgs<Scalar> &basePoint,
00206   const int p_index,
00207   const int g_index
00208   )
00209 {
00210 
00211   typedef Thyra::ModelEvaluatorBase MEB;
00212 
00213   // ToDo: Validate input!
00214 
00215   responseFunc_ = responseFunc;
00216   basePoint_ = basePoint;
00217   p_index_ = p_index;
00218   g_index_ = g_index;
00219 
00220   p_space_ = responseFunc_->get_p_space(p_index_);
00221   g_space_ = responseFunc_->get_g_space(g_index_);
00222 
00223 
00224   MEB::InArgs<Scalar>
00225     responseInArgs = responseFunc_->createInArgs();
00226   response_func_supports_x_dot_ =
00227     responseInArgs.supports(MEB::IN_ARG_x_dot);
00228   MEB::OutArgs<Scalar>
00229     responseOutArgs = responseFunc_->createOutArgs();
00230   response_func_supports_D_x_dot_ =
00231     !responseOutArgs.supports(MEB::OUT_ARG_DgDx_dot,g_index_).none();
00232   response_func_supports_D_p_ =
00233     !responseOutArgs.supports(MEB::OUT_ARG_DgDp,g_index_,p_index_).none();
00234 
00235   clearCache();
00236 
00237 }
00238 
00239 
00240 template<class Scalar>
00241 void ForwardResponseSensitivityComputer<Scalar>::resetResponseFunction(
00242   const RCP<const Thyra::ModelEvaluator<Scalar> > &responseFunc,
00243   const Thyra::ModelEvaluatorBase::InArgs<Scalar> &basePoint
00244   )
00245 {
00246   // ToDo: Validate that responseFunc is the same structure as the one already
00247   // set!
00248   responseFunc_ = responseFunc;
00249   basePoint_ = basePoint;
00250 }
00251 
00252 
00253 template<class Scalar>
00254 const RCP<Thyra::VectorBase<Scalar> >
00255 ForwardResponseSensitivityComputer<Scalar>::create_g_hat() const
00256 {
00257   return Thyra::createMember(g_space_);
00258 }
00259 
00260 
00261 template<class Scalar>
00262 const RCP<Thyra::MultiVectorBase<Scalar> >
00263 ForwardResponseSensitivityComputer<Scalar>::create_D_g_hat_D_p() const
00264 {
00265   return Thyra::createMembers(g_space_,p_space_->dim());
00266 }
00267 
00268 
00269 template<class Scalar>
00270 void ForwardResponseSensitivityComputer<Scalar>::computeResponse(
00271   const Thyra::VectorBase<Scalar> *x_dot,
00272   const Thyra::VectorBase<Scalar> &x,
00273   const Scalar t,
00274   Thyra::VectorBase<Scalar> *g_hat
00275   ) const
00276 {
00277   computeResponseAndSensitivityImpl(x_dot,0,x,0,t,g_hat,0);
00278 }
00279 
00280 
00281 template<class Scalar>
00282 void ForwardResponseSensitivityComputer<Scalar>::computeResponseAndSensitivity(
00283   const Thyra::VectorBase<Scalar> *x_dot,
00284   const Thyra::MultiVectorBase<Scalar> *S_dot,
00285   const Thyra::VectorBase<Scalar> &x,
00286   const Thyra::MultiVectorBase<Scalar> &S,
00287   const Scalar t,
00288   Thyra::VectorBase<Scalar> *g_hat,
00289   Thyra::MultiVectorBase<Scalar> *D_g_hat_D_p
00290   ) const
00291 {
00292   computeResponseAndSensitivityImpl(x_dot,S_dot,x,&S,t,g_hat,D_g_hat_D_p);
00293 }
00294 
00295 
00296 // private
00297 
00298 
00299 template<class Scalar>
00300 void ForwardResponseSensitivityComputer<Scalar>::clearCache()
00301 {
00302   D_g_D_x_dot_ = Teuchos::null;
00303   D_g_D_x_ = Teuchos::null;
00304   D_g_D_p_ = Teuchos::null;
00305 }
00306 
00307 
00308 template<class Scalar>
00309 void ForwardResponseSensitivityComputer<Scalar>::createCache(
00310   const bool computeSens
00311   ) const
00312 {
00313   if (computeSens) {
00314     if (response_func_supports_D_x_dot_ && is_null(D_g_D_x_dot_))
00315       D_g_D_x_dot_ = responseFunc_->create_DgDx_dot_op(g_index_);
00316     D_g_D_x_ = responseFunc_->create_DgDx_op(g_index_);
00317     if (response_func_supports_D_p_ && is_null(D_g_D_p_))
00318       D_g_D_p_ = createMembers(g_space_,p_space_->dim());
00319   }
00320 }
00321 
00322 
00323 template<class Scalar>
00324 void ForwardResponseSensitivityComputer<Scalar>::computeResponseAndSensitivityImpl(
00325   const Thyra::VectorBase<Scalar> *x_dot,
00326   const Thyra::MultiVectorBase<Scalar> *S_dot,
00327   const Thyra::VectorBase<Scalar> &x,
00328   const Thyra::MultiVectorBase<Scalar> *S,
00329   const Scalar t,
00330   Thyra::VectorBase<Scalar> *g_hat,
00331   Thyra::MultiVectorBase<Scalar> *D_g_hat_D_p
00332   ) const
00333 {
00334 
00335   using Teuchos::rcp;
00336   using Teuchos::ptr;
00337   using Teuchos::includesVerbLevel;
00338   typedef ScalarTraits<Scalar> ST;
00339   using Thyra::apply;
00340   using Thyra::Vp_V;
00341   typedef Thyra::ModelEvaluatorBase MEB;
00342 
00343   //
00344   // A) Setup for output
00345   //
00346 
00347   const RCP<Teuchos::FancyOStream> out = this->getOStream();
00348   const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00349 
00350   const bool trace =
00351     out.get() && includesVerbLevel(verbLevel,Teuchos::VERB_LOW);
00352   const bool print_norms =
00353     out.get() && includesVerbLevel(verbLevel,Teuchos::VERB_MEDIUM);
00354   const bool print_x =
00355     out.get() && includesVerbLevel(verbLevel,Teuchos::VERB_EXTREME);
00356 
00357   //
00358   // B) Initialize storage
00359   //
00360 
00361   const bool computeSens = ( D_g_hat_D_p != 0 );
00362   createCache(computeSens);
00363 
00364   //
00365   // C) Evaluate the response function
00366   //
00367 
00368   //
00369   // C.1) Setup input/output and evaluate the response function
00370   //
00371 
00372   if (trace)
00373     *out << "\nEvaluating response function at time t = " << t << " ...\n";
00374    
00375   // C.1.a) Setup the input arguments
00376    
00377   MEB::InArgs<Scalar> responseInArgs = responseFunc_->createInArgs();
00378   responseInArgs.setArgs(basePoint_);
00379   responseInArgs.set_x(rcp(&x,false));
00380   if (response_func_supports_x_dot_)
00381     responseInArgs.set_x_dot(rcp(x_dot,false));
00382       
00383   // C.1.b) Setup output arguments
00384   
00385   MEB::OutArgs<Scalar> responseOutArgs = responseFunc_->createOutArgs();
00386   
00387   if (g_hat)
00388     responseOutArgs.set_g(g_index_,rcp(g_hat,false));
00389   
00390   if (computeSens) {
00391         
00392     // D_g_D_x_dot
00393     if (response_func_supports_D_x_dot_) {
00394       responseOutArgs.set_DgDx_dot(
00395         g_index_,
00396         MEB::Derivative<Scalar>(D_g_D_x_dot_)
00397         );
00398     }
00399     
00400     // D_g_D_x
00401     responseOutArgs.set_DgDx(
00402       g_index_,
00403       MEB::Derivative<Scalar>(D_g_D_x_)
00404       );
00405     
00406     // D_g_D_p
00407     if (response_func_supports_D_p_) {
00408       responseOutArgs.set_DgDp(
00409         g_index_, p_index_,
00410         MEB::Derivative<Scalar>(D_g_D_p_,MEB::DERIV_MV_BY_COL)
00411         );
00412     }
00413     
00414   }
00415       
00416   // C.1.c) Evaluate the response function k
00417 
00418   {
00419     RYTHMOS_FUNC_TIME_MONITOR("Rythmos:ForwardResponseSensitivityComputer::evalModel: evalResponse");
00420     responseFunc_->evalModel( responseInArgs, responseOutArgs );
00421   }
00422   
00423   // C.1.d) Print the outputs just coputed
00424   
00425   if (print_norms) {
00426     if (g_hat)
00427       *out << "\n||g_hat||inf = " << norm_inf(*g_hat) << endl;
00428     if (computeSens && !is_null(D_g_D_p_))
00429       *out << "\n||D_g_D_p_||inf = " << norms_inf(*D_g_D_p_) << endl;
00430   }
00431   
00432   if ( g_hat && (dumpSensitivities_ || print_x) )
00433     *out << "\ng_hat = " << *g_hat;
00434       
00435   if (computeSens && print_x) {
00436     if (!is_null(D_g_D_x_dot_))
00437       *out << "\nD_g_D_x_dot = " << *D_g_D_x_ << endl;
00438     if (!is_null(D_g_D_x_))
00439       *out << "\nD_g_D_x = " << *D_g_D_x_ << endl;
00440     if (!is_null(D_g_D_p_))
00441       *out << "\nD_g_D_p = " << *D_g_D_p_ << endl;
00442   }
00443   
00444   //
00445   // C.2) Assemble the output response function sensitivity D_d_hat_D_p
00446   //
00447 
00448   // D_g_hat_D_p = DgDx_dot * S_dot + DgDx * S + DgDp 
00449   
00450   if (computeSens) {
00451 
00452     RYTHMOS_FUNC_TIME_MONITOR("Rythmos:ForwardResponseSensitivityComputer::evalModel: computeSens");
00453     
00454     if (trace)
00455       *out << "\nD_g_hat_D_p = DgDx_dot * S_dot + DgDx * S + DgDp ...\n";
00456     
00457     assign( ptr(D_g_hat_D_p), ST::zero() );
00458       
00459     // D_g_hat_D_p += DgDx_dot * S_dot
00460     if (response_func_supports_D_x_dot_) {
00461       apply( *D_g_D_x_dot_, Thyra::NOTRANS, *S_dot,
00462         ptr(D_g_hat_D_p), ST::one(), ST::one() );
00463     }
00464     
00465     // D_g_hat_D_p += DgDx * S
00466     apply( *D_g_D_x_, Thyra::NOTRANS, *S,
00467       ptr(D_g_hat_D_p), ST::one(), ST::one() );
00468     
00469     // D_g_hat_D_p += DgDp
00470     if (response_func_supports_D_p_) {
00471       Vp_V( ptr(D_g_hat_D_p), *D_g_D_p_ );
00472     }
00473       
00474     if (dumpSensitivities_ || print_x)
00475       *out << "\nD_g_hat_D_p = "
00476            << Teuchos::describe(*D_g_hat_D_p,Teuchos::VERB_EXTREME);
00477     
00478   }
00479 
00480 }
00481 
00482 
00483 } // namespace Rythmos
00484 
00485 
00486 #endif // RYTHMOS_FORWARD_RESPONSE_SENSITIVITY_COMPUTER_HPP
 All Classes Functions Variables Typedefs Friends