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

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