Rythmos - Transient Integration for Differential Equations Version of the Day
Rythmos_ForwardSensitivityIntegratorAsModelEvaluator.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 
00030 #ifndef RYTHMOS_FORWARD_SENSITIVITY_INTEGRATOR_AS_MODEL_EVALUATOR_HPP
00031 #define RYTHMOS_FORWARD_SENSITIVITY_INTEGRATOR_AS_MODEL_EVALUATOR_HPP
00032 
00033 
00034 #include "Thyra_ResponseOnlyModelEvaluatorBase.hpp"
00035 
00036 #include "Rythmos_StepperBase.hpp"
00037 #include "Rythmos_IntegratorBase.hpp"
00038 #include "Rythmos_InterpolationBufferHelpers.hpp"
00039 #include "Rythmos_ForwardSensitivityStepper.hpp"
00040 #include "Rythmos_ForwardResponseSensitivityComputer.hpp"
00041 #include "Rythmos_extractStateAndSens.hpp"
00042 #include "Thyra_ModelEvaluatorDelegatorBase.hpp"
00043 #include "Thyra_DefaultMultiVectorProductVectorSpace.hpp"
00044 #include "Thyra_DefaultProductVectorSpace.hpp"
00045 #include "Teuchos_ParameterListAcceptor.hpp"
00046 #include "Teuchos_VerboseObjectParameterListHelpers.hpp"
00047 #include "Teuchos_Assert.hpp"
00048 
00049 
00050 namespace Rythmos {
00051 
00052 
00053 namespace ForwardSensitivityIntegratorAsModelEvaluatorTypes {
00055 enum EResponseType {
00057   RESPONSE_TYPE_SUM,
00059   RESPONSE_TYPE_BLOCK
00060 };
00061 } // namespace ForwardSensitivityIntegratorAsModelEvaluatorTypes
00062 
00063 
00138 template<class Scalar>
00139 class ForwardSensitivityIntegratorAsModelEvaluator
00140   : virtual public Thyra::ResponseOnlyModelEvaluatorBase<Scalar>
00141   , virtual public Teuchos::ParameterListAcceptor
00142 {
00143 public:
00144 
00146 
00149 
00151   ForwardSensitivityIntegratorAsModelEvaluator();
00152 
00206   void initialize(
00207     const RCP<StepperBase<Scalar> > &stateStepper,
00208     const RCP<IntegratorBase<Scalar> > &stateIntegrator,
00209     const RCP<ForwardSensitivityStepper<Scalar> > &stateAndSensStepper,
00210     const RCP<IntegratorBase<Scalar> > &stateAndSensIntegrator,
00211     const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateAndSensInitCond,
00212     const Array<Scalar> &responseTimes,
00213     const Array<RCP<const Thyra::ModelEvaluator<Scalar> > > &responseFuncs,
00214     const Array<Thyra::ModelEvaluatorBase::InArgs<Scalar> > &responseFuncBasePoints,
00215     const ForwardSensitivityIntegratorAsModelEvaluatorTypes::EResponseType responseType
00216     );
00217 
00219   const Array<Scalar>& getResponseTimes() const;
00220 
00222 
00223 
00226 
00234   void setParameterList(RCP<Teuchos::ParameterList> const& paramList);
00236   RCP<Teuchos::ParameterList> getNonconstParameterList();
00238   RCP<Teuchos::ParameterList> unsetParameterList();
00240   RCP<const Teuchos::ParameterList> getParameterList() const;
00249   RCP<const Teuchos::ParameterList> getValidParameters() const;
00250 
00252 
00255 
00257   int Np() const;
00259   int Ng() const;
00261   RCP<const Thyra::VectorSpaceBase<Scalar> > get_p_space(int l) const;
00263   RCP<const Thyra::VectorSpaceBase<Scalar> > get_g_space(int j) const;
00265   Thyra::ModelEvaluatorBase::InArgs<Scalar> createInArgs() const;
00266 
00268 
00269 private:
00270 
00273 
00275   Thyra::ModelEvaluatorBase::OutArgs<Scalar> createOutArgsImpl() const;
00277   void evalModelImpl(
00278     const Thyra::ModelEvaluatorBase::InArgs<Scalar>& inArgs,
00279     const Thyra::ModelEvaluatorBase::OutArgs<Scalar>& outArgs
00280     ) const;
00281 
00283 
00284 private:
00285 
00286   // //////////////////////
00287   // Private types
00288 
00289   typedef Teuchos::Array<RCP<const Thyra::VectorSpaceBase<Scalar> > > SpaceArray_t;
00290   
00291 
00292   // //////////////////////
00293   // Private data members
00294 
00295   RCP<Teuchos::ParameterList>  paramList_;
00296   bool dumpSensitivities_;
00297 
00298   RCP<StepperBase<Scalar> > stateStepper_;
00299   RCP<IntegratorBase<Scalar> > stateIntegrator_;
00300   mutable Thyra::ModelEvaluatorBase::InArgs<Scalar> stateInitCond_;
00301 
00302   RCP<ForwardSensitivityStepper<Scalar> > stateAndSensStepper_;
00303   RCP<IntegratorBase<Scalar> > stateAndSensIntegrator_;
00304   mutable Thyra::ModelEvaluatorBase::InArgs<Scalar> stateAndSensInitCond_;
00305 
00306   Array<Scalar> responseTimes_;
00307   Array<RCP<const Thyra::ModelEvaluator<Scalar> > > responseFuncs_;
00308   mutable Array<Thyra::ModelEvaluatorBase::InArgs<Scalar> > responseFuncInArgs_;
00309   ForwardSensitivityIntegratorAsModelEvaluatorTypes::EResponseType responseType_;
00310   bool response_func_supports_x_dot_;
00311   bool response_func_supports_D_x_dot_;
00312   bool response_func_supports_D_p_;
00313 
00314   int p_index_;
00315   int g_index_;
00316   Scalar finalTime_;
00317 
00318   int Np_;
00319   int Ng_;
00320 
00321   SpaceArray_t g_space_;
00322   SpaceArray_t p_space_;
00323 
00324   static const std::string dumpSensitivities_name_;
00325   static const bool dumpSensitivities_default_;
00326 
00327 };
00328 
00329 
00331 template<class Scalar>
00332 RCP<ForwardSensitivityIntegratorAsModelEvaluator<Scalar> >
00333 forwardSensitivityIntegratorAsModelEvaluator(
00334   const RCP<StepperBase<Scalar> > &stateStepper,
00335   const RCP<IntegratorBase<Scalar> > &stateIntegrator,
00336   const RCP<ForwardSensitivityStepper<Scalar> > &stateAndSensStepper,
00337   const RCP<IntegratorBase<Scalar> > &stateAndSensIntegrator,
00338   const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateAndSensInitCond,
00339   const Array<Scalar> &responseTimes,
00340   const Array<RCP<const Thyra::ModelEvaluator<Scalar> > > &responseFuncs,
00341   const Array<Thyra::ModelEvaluatorBase::InArgs<Scalar> > &responseFuncBasePoints,
00342   const ForwardSensitivityIntegratorAsModelEvaluatorTypes::EResponseType responseType
00343   )
00344 {
00345   using Teuchos::rcp;
00346   RCP<ForwardSensitivityIntegratorAsModelEvaluator<Scalar> >
00347     sensIntegratorAsModelEvaluator = rcp(new ForwardSensitivityIntegratorAsModelEvaluator<Scalar>());
00348   sensIntegratorAsModelEvaluator->initialize(
00349     stateStepper, stateIntegrator,
00350     stateAndSensStepper, stateAndSensIntegrator,
00351     stateAndSensInitCond,
00352     responseTimes, responseFuncs,
00353     responseFuncBasePoints, responseType
00354     );
00355   return sensIntegratorAsModelEvaluator;
00356 }
00357 
00358 
00359 // ////////////////////////
00360 // Definitions
00361 
00362 
00363 // Static members
00364 
00365 
00366 template<class Scalar>
00367 const std::string
00368 ForwardSensitivityIntegratorAsModelEvaluator<Scalar>::dumpSensitivities_name_
00369 = "Dump Sensitivities";
00370 
00371 template<class Scalar>
00372 const bool
00373 ForwardSensitivityIntegratorAsModelEvaluator<Scalar>::dumpSensitivities_default_
00374 = false;
00375 
00376 
00377 // Constructors, Initialization, Misc.
00378 
00379 
00380 template<class Scalar>
00381 ForwardSensitivityIntegratorAsModelEvaluator<Scalar>::ForwardSensitivityIntegratorAsModelEvaluator()
00382   :dumpSensitivities_(dumpSensitivities_default_),
00383    responseType_(ForwardSensitivityIntegratorAsModelEvaluatorTypes::RESPONSE_TYPE_SUM),
00384    response_func_supports_x_dot_(false),
00385    response_func_supports_D_x_dot_(false),
00386    response_func_supports_D_p_(false),
00387    p_index_(-1),
00388    g_index_(-1),
00389    finalTime_(-std::numeric_limits<Scalar>::max()),
00390    Np_(0),
00391    Ng_(0)
00392 {}
00393 
00394 
00395 template<class Scalar>
00396 void ForwardSensitivityIntegratorAsModelEvaluator<Scalar>::initialize(
00397   const RCP<StepperBase<Scalar> > &stateStepper,
00398   const RCP<IntegratorBase<Scalar> > &stateIntegrator,
00399   const RCP<ForwardSensitivityStepper<Scalar> > &stateAndSensStepper,
00400   const RCP<IntegratorBase<Scalar> > &stateAndSensIntegrator,
00401   const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateAndSensInitCond,
00402   const Array<Scalar> &responseTimes,
00403   const Array<RCP<const Thyra::ModelEvaluator<Scalar> > > &responseFuncs,
00404   const Array<Thyra::ModelEvaluatorBase::InArgs<Scalar> > &responseFuncBasePoints,
00405   const ForwardSensitivityIntegratorAsModelEvaluatorTypes::EResponseType responseType
00406   )
00407 {
00408 
00409   using Teuchos::as;
00410   typedef Thyra::ModelEvaluatorBase MEB;
00411   namespace FSIAMET = ForwardSensitivityIntegratorAsModelEvaluatorTypes;
00412 
00413   //
00414   // A) Validate and set input
00415   //
00416 
00417 #ifdef HAVE_RYTHMOS_DEBUG
00418   const int numResponseTimes = responseTimes.size();
00419 
00420   TEUCHOS_TEST_FOR_EXCEPT(is_null(stateStepper));
00421   TEUCHOS_TEST_FOR_EXCEPT(is_null(stateIntegrator));
00422   TEUCHOS_TEST_FOR_EXCEPT(is_null(stateAndSensStepper));
00423   TEUCHOS_TEST_FOR_EXCEPT(is_null(stateAndSensInitCond.get_x()));
00424   TEUCHOS_TEST_FOR_EXCEPT(is_null(stateAndSensInitCond.get_x_dot()));
00425   TEUCHOS_TEST_FOR_EXCEPT( !( numResponseTimes > 0 ) );
00426   assertTimePointsAreSorted(responseTimes);
00427   TEUCHOS_TEST_FOR_EXCEPT( as<int>(responseFuncs.size()) != numResponseTimes );
00428   TEUCHOS_TEST_FOR_EXCEPT( as<int>(responseFuncBasePoints.size()) != numResponseTimes );
00429   // ToDo: Assert that all of the observation models have the same response
00430   // function spaces so that they can be added together!
00431 #endif // HAVE_RYTHMOS_DEBUG
00432 
00433   stateStepper_ = stateStepper;
00434   stateAndSensStepper_ = stateAndSensStepper;
00435   stateAndSensInitCond_ = stateAndSensInitCond;
00436   stateIntegrator_ = stateIntegrator;
00437   if (!is_null(stateAndSensIntegrator))
00438     stateAndSensIntegrator_ = stateAndSensIntegrator;
00439   else
00440     stateAndSensIntegrator_ = stateIntegrator_->cloneIntegrator().assert_not_null();
00441   responseTimes_ = responseTimes;
00442   responseFuncs_ = responseFuncs;
00443   responseFuncInArgs_ = responseFuncBasePoints;
00444   responseType_ = responseType;
00445   
00446   finalTime_ = responseTimes_.back();
00447   
00448   //
00449   // B) Set the initial condition for the state-only problem
00450   //
00451 
00452   const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >
00453     stateModel = stateStepper_->getModel();
00454 
00455   // Get base point pased in with no x or x_dot
00456   MEB::InArgs<Scalar>
00457     basePoint_no_x = stateAndSensInitCond_;
00458   basePoint_no_x.set_x(Teuchos::null);
00459   basePoint_no_x.set_x_dot(Teuchos::null);
00460 
00461   // Create an empty InArgs for the state model/stepper
00462   stateInitCond_ = stateModel->createInArgs();
00463 
00464   // Set the base point (except x, x_dot).
00465   stateInitCond_.setArgs(basePoint_no_x);
00466 
00467   // Set x and x_dot
00468   const RCP<const Thyra::ProductVectorBase<Scalar> >
00469     x_bar_init = Thyra::productVectorBase<Scalar>(
00470       stateAndSensInitCond_.get_x()
00471       ),
00472     x_bar_dot_init = Thyra::productVectorBase<Scalar>(
00473       stateAndSensInitCond_.get_x_dot()
00474       );
00475   stateInitCond_.set_x(x_bar_init->getVectorBlock(0));
00476   stateInitCond_.set_x_dot(x_bar_dot_init->getVectorBlock(0));
00477 
00478   //
00479   // C) Set up the info for this model evaluators interface
00480   //
00481 
00482   Np_ = 1;
00483   p_index_ = getParameterIndex(*stateAndSensStepper_);
00484   p_space_.clear();
00485   p_space_.push_back(stateModel->get_p_space(p_index_));
00486   
00487   Ng_ = 1;
00488   g_index_ = 0; // ToDo: Accept this from input!
00489   g_space_.clear();
00490 
00491   if (responseType_ == FSIAMET::RESPONSE_TYPE_SUM) {
00492     g_space_.push_back(responseFuncs[0]->get_g_space(0));
00493   }
00494   else if (responseType_ == FSIAMET::RESPONSE_TYPE_BLOCK) {
00495     g_space_.push_back(
00496       Thyra::productVectorSpace(
00497         responseFuncs[0]->get_g_space(g_index_), responseFuncs.size()
00498         )
00499       );
00500   }
00501 
00502   MEB::InArgs<Scalar>
00503     responseInArgs = responseFuncs[0]->createInArgs();
00504   response_func_supports_x_dot_ =
00505     responseInArgs.supports(MEB::IN_ARG_x_dot);
00506   MEB::OutArgs<Scalar>
00507     responseOutArgs = responseFuncs[0]->createOutArgs();
00508   response_func_supports_D_x_dot_ =
00509     !responseOutArgs.supports(MEB::OUT_ARG_DgDx_dot,g_index_).none();
00510   response_func_supports_D_p_ =
00511     !responseOutArgs.supports(MEB::OUT_ARG_DgDp,g_index_,p_index_).none();
00512   
00513 }
00514 
00515 
00516 template<class Scalar>
00517 const Array<Scalar>&
00518 ForwardSensitivityIntegratorAsModelEvaluator<Scalar>::getResponseTimes() const
00519 {
00520   return responseTimes_;
00521 }
00522 
00523 
00524 // Overridden from Teuchos::ParameterListAcceptor
00525 
00526 
00527 template<class Scalar>
00528 void ForwardSensitivityIntegratorAsModelEvaluator<Scalar>::setParameterList(
00529   RCP<Teuchos::ParameterList> const& paramList
00530   )
00531 {
00532   TEUCHOS_TEST_FOR_EXCEPT(0==paramList.get());
00533   paramList->validateParameters(*getValidParameters());
00534   paramList_ = paramList;
00535   dumpSensitivities_ = paramList_->get(
00536     dumpSensitivities_name_, dumpSensitivities_default_);
00537   Teuchos::readVerboseObjectSublist(&*paramList_,this);
00538 #ifdef HAVE_RYTHMOS_DEBUG
00539   paramList_->validateParameters(*getValidParameters());
00540 #endif // HAVE_RYTHMOS_DEBUG
00541 
00542 }
00543 
00544 
00545 template<class Scalar>
00546 RCP<Teuchos::ParameterList>
00547 ForwardSensitivityIntegratorAsModelEvaluator<Scalar>::getNonconstParameterList()
00548 {
00549   return paramList_;
00550 }
00551 
00552 
00553 template<class Scalar>
00554 RCP<Teuchos::ParameterList>
00555 ForwardSensitivityIntegratorAsModelEvaluator<Scalar>::unsetParameterList()
00556 {
00557   RCP<Teuchos::ParameterList> _paramList = paramList_;
00558   paramList_ = Teuchos::null;
00559   return _paramList;
00560 }
00561 
00562 
00563 template<class Scalar>
00564 RCP<const Teuchos::ParameterList>
00565 ForwardSensitivityIntegratorAsModelEvaluator<Scalar>::getParameterList() const
00566 {
00567   return paramList_;
00568 }
00569 
00570 
00571 template<class Scalar>
00572 RCP<const Teuchos::ParameterList>
00573 ForwardSensitivityIntegratorAsModelEvaluator<Scalar>::getValidParameters() const
00574 {
00575   static RCP<ParameterList> validPL;
00576   if (is_null(validPL)) {
00577     RCP<ParameterList> pl = Teuchos::parameterList();
00578     pl->set(dumpSensitivities_name_, dumpSensitivities_default_,
00579       "Set to true to have the individual sensitivities printed to\n"
00580       "the output stream."
00581       );
00582     Teuchos::setupVerboseObjectSublist(&*pl);
00583     validPL = pl;
00584   }
00585   return validPL;
00586 }
00587 
00588 
00589 // Public functions overridden from ModelEvaulator
00590 
00591 
00592 template<class Scalar>
00593 int ForwardSensitivityIntegratorAsModelEvaluator<Scalar>::Np() const
00594 {
00595   return Np_;
00596 }
00597 
00598 
00599 template<class Scalar>
00600 int ForwardSensitivityIntegratorAsModelEvaluator<Scalar>::Ng() const
00601 {
00602   return Ng_;
00603 }
00604 
00605 
00606 template<class Scalar>
00607 RCP<const Thyra::VectorSpaceBase<Scalar> >
00608 ForwardSensitivityIntegratorAsModelEvaluator<Scalar>::get_p_space(int l) const
00609 {
00610 #ifdef HAVE_RYTHMOS_DEBUG
00611   TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( l, 0, Np_ );
00612 #endif
00613   return p_space_[l];
00614 }
00615 
00616 
00617 template<class Scalar>
00618 RCP<const Thyra::VectorSpaceBase<Scalar> >
00619 ForwardSensitivityIntegratorAsModelEvaluator<Scalar>::get_g_space(int j) const
00620 {
00621 #ifdef HAVE_RYTHMOS_DEBUG
00622   TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( j, 0, Ng_ );
00623 #endif
00624   return g_space_[j];
00625 }
00626 
00627 
00628 template<class Scalar>
00629 Thyra::ModelEvaluatorBase::InArgs<Scalar>
00630 ForwardSensitivityIntegratorAsModelEvaluator<Scalar>::createInArgs() const
00631 {
00632   typedef Thyra::ModelEvaluatorBase MEB;
00633   MEB::InArgsSetup<Scalar> inArgs;
00634   inArgs.setModelEvalDescription(this->description());
00635   inArgs.set_Np(Np_);
00636   return inArgs;
00637 }
00638 
00639 
00640 // Private functions overridden from ModelEvaulatorDefaultBase
00641 
00642 
00643 template<class Scalar>
00644 Thyra::ModelEvaluatorBase::OutArgs<Scalar>
00645 ForwardSensitivityIntegratorAsModelEvaluator<Scalar>::createOutArgsImpl() const
00646 {
00647   typedef Thyra::ModelEvaluatorBase MEB;
00648   typedef MEB::DerivativeSupport DS;
00649   namespace FSIAMET = ForwardSensitivityIntegratorAsModelEvaluatorTypes;
00650   MEB::OutArgsSetup<Scalar> outArgs;
00651   outArgs.setModelEvalDescription(this->description());
00652   outArgs.set_Np_Ng(Np_,Ng_);
00653   outArgs.setSupports(MEB::OUT_ARG_DgDp, 0 ,0, MEB::DERIV_MV_BY_COL);
00654   return outArgs;
00655 }
00656 
00657 
00658 template<class Scalar>
00659 void ForwardSensitivityIntegratorAsModelEvaluator<Scalar>::evalModelImpl(
00660   const Thyra::ModelEvaluatorBase::InArgs<Scalar>& inArgs,
00661   const Thyra::ModelEvaluatorBase::OutArgs<Scalar>& outArgs
00662   ) const
00663 {
00664 
00665   using Teuchos::as;
00666   using Teuchos::null;
00667   using Teuchos::describe;
00668   using Teuchos::OSTab;
00669   using Teuchos::rcp_dynamic_cast;
00670   typedef Teuchos::ScalarTraits<Scalar> ST;
00671   typedef Teuchos::VerboseObjectTempState<InterpolationBufferBase<Scalar> > VOTSSB;
00672   typedef Thyra::ModelEvaluatorBase MEB;
00673   namespace FSIAMET = ForwardSensitivityIntegratorAsModelEvaluatorTypes;
00674 
00675   //
00676   // Stream output and other basic setup
00677   //
00678   
00679   THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_GEN_BEGIN(
00680     "Rythmos::ForwardSensitivityIntegratorAsModelEvaluator", inArgs, outArgs, null
00681     );
00682   VOTSSB stateIntegrator_outputTempState(
00683     stateIntegrator_,out,incrVerbLevel(verbLevel,-1));
00684   VOTSSB stateAndSensIntegrator_outputTempState(
00685     stateAndSensIntegrator_,out,incrVerbLevel(verbLevel,-1));
00686 
00687   const bool trace =
00688     out.get() && includesVerbLevel(localVerbLevel,Teuchos::VERB_LOW);
00689   const bool print_norms =
00690     out.get() && includesVerbLevel(localVerbLevel,Teuchos::VERB_MEDIUM);
00691   const bool print_x =
00692     out.get() && includesVerbLevel(localVerbLevel,Teuchos::VERB_EXTREME);
00693 
00694   const RCP<const Thyra::ModelEvaluator<Scalar> >
00695     stateModel = stateStepper_->getModel();
00696   
00697   //
00698   // A) Process OutArgs first to see what functions we will be computing
00699   //
00700   
00701   RCP<Thyra::VectorBase<Scalar> >
00702     d_hat = outArgs.get_g(0);
00703 
00704   MEB::Derivative<Scalar>
00705     D_d_hat_D_p = outArgs.get_DgDp(0,p_index_);
00706     
00707   // Integrate state+sens or just state?
00708   const bool integrateStateAndSens = !D_d_hat_D_p.isEmpty();
00709 
00710   // Shortcut exit if no output was requested
00711   if ( is_null(d_hat) && D_d_hat_D_p.isEmpty() ) {
00712     if (trace)
00713       *out << "\nSkipping evaluation since no outputs were requested!\n";
00714     return;
00715   }
00716     
00717   //
00718   // B) Process the InArgs knowing if we will integrate just the state or the
00719   // state+sens.
00720   //
00721     
00722   const RCP<const Thyra::VectorBase<Scalar> >
00723     p = inArgs.get_p(0).assert_not_null();
00724     
00725   if (integrateStateAndSens) {
00726     if (trace)
00727       *out << "\nComputing D_d_hat_d_p by integrating state+sens ...\n";
00728     stateAndSensInitCond_.set_p(p_index_,p);
00729   }
00730   else {
00731     if (trace)
00732       *out << "\nComputing just d_hat by integrating the state ...\n";
00733     stateInitCond_.set_p(p_index_,p);
00734   }
00735     
00736   //
00737   // C) Setup the stepper and the integrator for state+sens or only state
00738   //
00739 
00740   RCP<IntegratorBase<Scalar> > integrator;
00741   if (integrateStateAndSens) {
00742     stateAndSensStepper_->setInitialCondition(stateAndSensInitCond_);
00743     stateAndSensIntegrator_->setStepper(stateAndSensStepper_,finalTime_);
00744     integrator = stateAndSensIntegrator_;
00745   }
00746   else {
00747     stateStepper_->setInitialCondition(stateInitCond_);
00748     stateIntegrator_->setStepper(stateStepper_,finalTime_);
00749     integrator = stateIntegrator_;
00750   }
00751     
00752   //
00753   // D) Setup for computing and processing the individual response functions
00754   //
00755 
00756   ForwardResponseSensitivityComputer<Scalar>
00757     forwardResponseSensitivityComputer;
00758   forwardResponseSensitivityComputer.setOStream(out);
00759   forwardResponseSensitivityComputer.setVerbLevel(localVerbLevel);
00760   forwardResponseSensitivityComputer.dumpSensitivities(dumpSensitivities_);
00761   forwardResponseSensitivityComputer.setResponseFunction(
00762     responseFuncs_[0],
00763     responseFuncs_[0]->createInArgs(), // Will replace this for each point!
00764     p_index_, g_index_
00765     );
00766 
00767   // D.a) Create storage for the individual response function ouptuts g_k
00768   // and its derivaitves
00769 
00770   RCP<Thyra::VectorBase<Scalar> > g_k;
00771   RCP<Thyra::MultiVectorBase<Scalar> > D_g_hat_D_p_k;
00772 
00773   if (!is_null(d_hat)) {
00774     g_k = forwardResponseSensitivityComputer.create_g_hat();
00775   }
00776   if (!D_d_hat_D_p.isEmpty()) {
00777     D_g_hat_D_p_k = forwardResponseSensitivityComputer.create_D_g_hat_D_p();
00778   }
00779 
00780   // D.b) Zero out d_hat and D_d_hat_D_p if we are doing a summation type of
00781   // evaluation
00782   if (responseType_ == FSIAMET::RESPONSE_TYPE_SUM) {
00783     if (!is_null(d_hat)) {
00784       assign( d_hat.ptr(), ST::zero() );
00785     }
00786     if (!D_d_hat_D_p.isEmpty()) {
00787       assign( D_d_hat_D_p.getMultiVector().ptr(), ST::zero() );
00788     }
00789   }
00790 
00791   // D.c) Get product vector and multi-vector interfaces if
00792   // we are just blocking up response functions
00793   RCP<Thyra::ProductVectorBase<Scalar> > prod_d_hat;
00794   RCP<Thyra::ProductMultiVectorBase<Scalar> > prod_D_d_hat_D_p_trans;
00795   if (responseType_ == FSIAMET::RESPONSE_TYPE_BLOCK) {
00796     prod_d_hat = rcp_dynamic_cast<Thyra::ProductVectorBase<Scalar> >(
00797       d_hat, true);
00798     prod_D_d_hat_D_p_trans = rcp_dynamic_cast<Thyra::ProductMultiVectorBase<Scalar> >(
00799       D_d_hat_D_p.getMultiVector(), true);
00800   }
00801     
00802   //
00803   // E) Run the integrator at the response time points and assemble
00804   //    the response function and/or the response function
00805   //    derivative
00806   //
00807 
00808   if (trace) *out << "\nStarting integration and assembly loop ...\n";
00809     
00810   const int numResponseFunctions = responseTimes_.size();
00811 
00812   for (int k = 0; k < numResponseFunctions; ++k ) {
00813       
00814     OSTab os_tab(out);
00815       
00816     const Scalar t = responseTimes_[k];
00817       
00818     //
00819     // E.1) Integrate up to t and get x_bar and x_bar_dot
00820     //
00821     // Note, x_bar and x_bar_dot may be the state or the state+sens!
00822       
00823     if (trace)
00824       *out << "\nIntegrating state (or state+sens) for response k = "
00825            << k << " for t = " << t << " ...\n";
00826       
00827     RCP<const Thyra::VectorBase<Scalar> > x_bar, x_bar_dot;
00828 
00829     {
00830       RYTHMOS_FUNC_TIME_MONITOR(
00831         "Rythmos:ForwardSensitivityIntegratorAsModelEvaluator::evalModel: integrate"
00832         );
00833       get_fwd_x_and_x_dot( *integrator, t, &x_bar, &x_bar_dot );
00834     }
00835       
00836     if (print_norms) {
00837       *out << "\n||x_bar||inf = " << norms_inf(*x_bar) << std::endl;
00838       *out << "\n||x_bar_dot||inf = " << norms_inf(*x_bar_dot) << std::endl;
00839     }
00840     if (print_x) {
00841       *out << "\nx_bar = " << *x_bar << std::endl;
00842       *out << "\nx_bar_dot = " << *x_bar_dot << std::endl;
00843     }
00844       
00845     // Extra underlying state and sensitivities
00846     RCP<const Thyra::VectorBase<Scalar> > x, x_dot;
00847     RCP<const Thyra::MultiVectorBase<Scalar> > S, S_dot;
00848     if (integrateStateAndSens) {
00849       extractStateAndSens( x_bar, x_bar_dot, &x, &S, &x_dot, &S_dot );
00850     }
00851     else {
00852       x = x_bar;
00853       x_dot = x_bar_dot;
00854     }
00855       
00856     //
00857     // E.2) Evaluate the response function g_k and its derivatives at this
00858     // time point
00859     //
00860       
00861     if (trace)
00862       *out << "\nEvaluating response function k = " << k << " at time t = " << t << " ...\n";
00863       
00864     RCP<const Thyra::ModelEvaluator<Scalar> > responseFunc = responseFuncs_[k];
00865       
00866     MEB::InArgs<Scalar> responseInArgs = responseFunc->createInArgs();
00867     responseInArgs.setArgs(responseFuncInArgs_[k]);
00868     responseInArgs.set_p(p_index_,p);
00869 
00870     forwardResponseSensitivityComputer.resetResponseFunction(
00871       responseFunc, responseInArgs
00872       );
00873 
00874     if (!is_null(D_g_hat_D_p_k)) {
00875       forwardResponseSensitivityComputer.computeResponseAndSensitivity(
00876         x_dot.get(), S_dot.get(), *x, *S, t, g_k.get(), &*D_g_hat_D_p_k
00877         );
00878     }
00879     else {
00880       forwardResponseSensitivityComputer.computeResponse(
00881         x_dot.get(), *x, t, g_k.get()
00882         );
00883     }
00884 
00885     //
00886     // E.3) Assemble the final response functions and derivatives
00887     //
00888 
00889     // E.3.a) Assemble d_hat from g_k
00890     if (!is_null(d_hat)) {
00891       if (responseType_ == FSIAMET::RESPONSE_TYPE_SUM) {
00892         if (trace) *out << "\nd_hat += g_k ...\n";
00893         Vp_V( d_hat.ptr(), *g_k );
00894       }
00895       else if (responseType_ == FSIAMET::RESPONSE_TYPE_BLOCK) {
00896         if (trace) *out << "\nd_hat["<<k<<"] = g_k ...\n";
00897         assign( prod_d_hat->getNonconstVectorBlock(k).ptr(), *g_k );
00898       }
00899     }
00900 
00901     // E.3.b) Assemble D_d_hat_Dp from D_g_hat_D_p_k
00902     if (!D_d_hat_D_p.isEmpty()) {
00903       if (responseType_ == FSIAMET::RESPONSE_TYPE_SUM) {
00904         if (trace) *out << "\nD_d_hat_D_p += D_g_hat_D_p_k ...\n";
00905         Vp_V( D_d_hat_D_p.getMultiVector().ptr(), *D_g_hat_D_p_k );
00906         if (dumpSensitivities_ || print_x)
00907           *out << "\nD_d_hat_D_p = "
00908                << Teuchos::describe(*D_d_hat_D_p.getMultiVector(),Teuchos::VERB_EXTREME);
00909       }
00910       else if (responseType_ == FSIAMET::RESPONSE_TYPE_BLOCK) {
00911         if (trace) *out << "\nD_d_hat_D_p["<<k<<"] = D_g_hat_D_p_k ...\n";
00912         assign(
00913           prod_D_d_hat_D_p_trans->getNonconstMultiVectorBlock(k).ptr(),
00914           *D_g_hat_D_p_k
00915           );
00916       }
00917     }
00918 
00919   }
00920     
00921   THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_END();
00922 
00923 }
00924 
00925 
00926 } // namespace Rythmos
00927 
00928 
00929 #endif // RYTHMOS_FORWARD_SENSITIVITY_INTEGRATOR_AS_MODEL_EVALUATOR_HPP
 All Classes Functions Variables Typedefs Friends