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 "Rythmos_IntegratorBase.hpp"
00034 #include "Thyra_ModelEvaluator.hpp" // Interface
00035 #include "Thyra_StateFuncModelEvaluatorBase.hpp" // Implementation
00036 #include "Thyra_DefaultProductVectorSpace.hpp"
00037 #include "Thyra_PhysicallyBlockedLinearOpWithSolveBase.hpp" // Interface
00038 #include "Thyra_DefaultBlockedTriangularLinearOpWithSolve.hpp" // Implementation
00039 #include "Thyra_ModelEvaluatorDelegatorBase.hpp"
00040 #include "Thyra_ModelEvaluatorHelpers.hpp"
00041 #include "Thyra_DefaultMultiVectorProductVectorSpace.hpp"
00042 #include "Thyra_DefaultMultiVectorProductVector.hpp"
00043 #include "Thyra_DefaultMultiVectorLinearOpWithSolve.hpp"
00044 #include "Teuchos_implicit_cast.hpp"
00045 #include "Teuchos_Assert.hpp"
00046 
00047 
00048 namespace Rythmos {
00049 
00050 
00297 template<class Scalar>
00298 class ForwardSensitivityModelEvaluator
00299   : virtual public Thyra::StateFuncModelEvaluatorBase<Scalar>
00300 {
00301 public:
00302 
00305 
00307   ForwardSensitivityModelEvaluator();
00308 
00324   void initializeStructure(
00325     const RCP<const Thyra::ModelEvaluator<Scalar> > &stateModel,
00326     const int p_index
00327     );
00328   
00330   RCP<const Thyra::ModelEvaluator<Scalar> >
00331   getStateModel() const;
00332   
00334   int get_p_index() const;
00335 
00351   void setStateIntegrator(
00352     const RCP<IntegratorBase<Scalar> > &stateIntegrator,
00353     const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint
00354     );
00355 
00395   void initializePointState(
00396     const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint,
00397     const RCP<const Thyra::LinearOpWithSolveBase<Scalar> > &W_tilde,
00398     const Scalar &coeff_x_dot,
00399     const Scalar &coeff_x,
00400     const RCP<const Thyra::LinearOpBase<Scalar> > &DfDx_dot = Teuchos::null,
00401     const RCP<const Thyra::MultiVectorBase<Scalar> > &DfDp = Teuchos::null
00402     );
00403 
00404   // 2007/05/22: rabartl: ToDo: Add an InterpolationBufferBase
00405   // stateInterpBuffer object to the initailizeState(...) function that can be
00406   // used to get x and x_dot at different points in time t.  Then, modify the
00407   // logic to recompute all of the needed matrices if t != t_base (as passed
00408   // in through stateBasePoint).  The values of x(t) and xdot(t) can then be
00409   // gotten from the stateInterpBuffer object!
00410   
00412 
00415 
00417   RCP<const Thyra::VectorSpaceBase<Scalar> > get_x_space() const;
00419   RCP<const Thyra::VectorSpaceBase<Scalar> > get_f_space() const;
00421   Thyra::ModelEvaluatorBase::InArgs<Scalar> getNominalValues() const;
00423   RCP<Thyra::LinearOpWithSolveBase<Scalar> > create_W() const;
00425   Thyra::ModelEvaluatorBase::InArgs<Scalar> createInArgs() const;
00426 
00428 
00431 
00433   void initializeState(
00434     const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint,
00435     const RCP<const Thyra::LinearOpWithSolveBase<Scalar> > &W_tilde,
00436     const Scalar &coeff_x_dot,
00437     const Scalar &coeff_x,
00438     const RCP<const Thyra::LinearOpBase<Scalar> > &DfDx_dot = Teuchos::null,
00439     const RCP<const Thyra::MultiVectorBase<Scalar> > &DfDp = Teuchos::null
00440     )
00441     {
00442       initializePointState( stateBasePoint, W_tilde, coeff_x_dot, coeff_x,
00443         DfDx_dot, DfDp );
00444     }
00445 
00447 
00448 private:
00449 
00452 
00454   Thyra::ModelEvaluatorBase::OutArgs<Scalar> createOutArgsImpl() const;
00456   void evalModelImpl(
00457     const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
00458     const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs
00459     ) const;
00460 
00462 
00463 private:
00464 
00465   // /////////////////////////
00466   // Private data members
00467 
00468   RCP<const Thyra::ModelEvaluator<Scalar> > stateModel_;
00469   int p_index_;
00470   int np_;
00471 
00472   RCP<IntegratorBase<Scalar> > stateIntegrator_;
00473 
00474   RCP<const Thyra::DefaultMultiVectorProductVectorSpace<Scalar> > s_bar_space_;
00475   RCP<const Thyra::DefaultMultiVectorProductVectorSpace<Scalar> > f_sens_space_;
00476 
00477   Thyra::ModelEvaluatorBase::InArgs<Scalar> nominalValues_;
00478 
00479   mutable Thyra::ModelEvaluatorBase::InArgs<Scalar> stateBasePoint_;
00480 
00481   mutable RCP<const Thyra::LinearOpWithSolveBase<Scalar> > W_tilde_;
00482   mutable Scalar coeff_x_dot_;
00483   mutable Scalar coeff_x_;
00484   mutable RCP<const Thyra::LinearOpBase<Scalar> > DfDx_dot_;
00485   mutable RCP<const Thyra::MultiVectorBase<Scalar> > DfDp_;
00486 
00487   mutable RCP<Thyra::LinearOpWithSolveBase<Scalar> > W_tilde_compute_;
00488   mutable RCP<Thyra::LinearOpBase<Scalar> > DfDx_dot_compute_;
00489   mutable RCP<Thyra::MultiVectorBase<Scalar> > DfDp_compute_;
00490 
00491   // /////////////////////////
00492   // Private member functions
00493 
00494   void wrapNominalValuesAndBounds();
00495 
00496   void computeDerivativeMatrices(
00497     const Thyra::ModelEvaluatorBase::InArgs<Scalar> &point
00498     ) const;
00499   
00500 };
00501 
00502 
00503 // /////////////////////////////////
00504 // Implementations
00505 
00506 
00507 // Constructors/Intializers/Accessors
00508 
00509 
00510 template<class Scalar>
00511 ForwardSensitivityModelEvaluator<Scalar>::ForwardSensitivityModelEvaluator()
00512   : p_index_(0), np_(-1)
00513 {}
00514 
00515 
00516 template<class Scalar>
00517 void ForwardSensitivityModelEvaluator<Scalar>::initializeStructure(
00518   const RCP<const Thyra::ModelEvaluator<Scalar> > &stateModel,
00519   const int p_index
00520   )
00521 {
00522 
00523   typedef Thyra::ModelEvaluatorBase MEB;
00524 
00525   //
00526   // Validate input
00527   //
00528 
00529   TEST_FOR_EXCEPT( is_null(stateModel) );
00530   TEST_FOR_EXCEPTION(
00531     !( 0 <= p_index && p_index < stateModel->Np() ), std::logic_error,
00532     "Error, p_index does not fall in the range [0,"<<(stateModel->Np()-1)<<"]!" );
00533   // ToDo: Validate support for DfDp!
00534 
00535   //
00536   // Set the input objects
00537   //
00538 
00539   stateModel_ = stateModel;
00540   p_index_ = p_index;
00541   np_ = stateModel_->get_p_space(p_index)->dim();
00542 
00543   //
00544   // Create the structure of the model
00545   //
00546 
00547   s_bar_space_ = Thyra::multiVectorProductVectorSpace(
00548     stateModel_->get_x_space(), np_
00549     );
00550 
00551   f_sens_space_ = Thyra::multiVectorProductVectorSpace(
00552     stateModel_->get_f_space(), np_
00553     );
00554 
00555   nominalValues_ = this->createInArgs();
00556 
00557   //
00558   // Wipe out matrix storage
00559   //
00560 
00561   stateBasePoint_ = MEB::InArgs<Scalar>();
00562   W_tilde_ = Teuchos::null;
00563   coeff_x_dot_ = 0.0;
00564   coeff_x_ = 0.0;
00565   DfDx_dot_ = Teuchos::null;
00566   DfDp_ = Teuchos::null;
00567   W_tilde_compute_ = Teuchos::null;
00568   DfDx_dot_compute_ = Teuchos::null;
00569   DfDp_compute_ = Teuchos::null;
00570 
00571 }
00572 
00573 
00574 template<class Scalar>
00575 RCP<const Thyra::ModelEvaluator<Scalar> >
00576 ForwardSensitivityModelEvaluator<Scalar>::getStateModel() const
00577 {
00578   return stateModel_;
00579 }
00580 
00581 
00582 template<class Scalar>
00583 int ForwardSensitivityModelEvaluator<Scalar>::get_p_index() const
00584 {
00585   return p_index_;
00586 }
00587 
00588 
00589 template<class Scalar>
00590 void ForwardSensitivityModelEvaluator<Scalar>::setStateIntegrator(
00591   const RCP<IntegratorBase<Scalar> > &stateIntegrator,
00592   const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint
00593   )
00594 {
00595   stateIntegrator_ = stateIntegrator.assert_not_null();
00596   stateBasePoint_ = stateBasePoint;
00597 }
00598 
00599 
00600 template<class Scalar>
00601 void ForwardSensitivityModelEvaluator<Scalar>::initializePointState(
00602   const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint,
00603   const RCP<const Thyra::LinearOpWithSolveBase<Scalar> > &W_tilde,
00604   const Scalar &coeff_x_dot,
00605   const Scalar &coeff_x,
00606   const RCP<const Thyra::LinearOpBase<Scalar> > &DfDx_dot,
00607   const RCP<const Thyra::MultiVectorBase<Scalar> > &DfDp
00608   )
00609 {
00610 
00611   typedef Thyra::ModelEvaluatorBase MEB;
00612 
00613 #ifdef RYTHMOS_DEBUG
00614   TEST_FOR_EXCEPTION(
00615     is_null(stateModel_), std::logic_error,
00616     "Error, you must call intializeStructure(...) before you call initializeState(...)"
00617     );
00618   TEST_FOR_EXCEPT( is_null(stateBasePoint.get_x()) );
00619   TEST_FOR_EXCEPT( is_null(stateBasePoint.get_x_dot()) );
00620   TEST_FOR_EXCEPT( is_null(stateBasePoint.get_p(p_index_)) );
00621   // What about the other parameter values?  We really can't say anything
00622   // about these and we can't check them.  They can be null just fine.
00623   if (!is_null(W_tilde)) {
00624     THYRA_ASSERT_VEC_SPACES("",*W_tilde->domain(),*stateModel_->get_x_space());
00625     THYRA_ASSERT_VEC_SPACES("",*W_tilde->range(),*stateModel_->get_f_space());
00626   }
00627   if (!is_null(DfDx_dot)) {
00628     THYRA_ASSERT_VEC_SPACES("",*DfDx_dot->domain(),*stateModel_->get_x_space());
00629     THYRA_ASSERT_VEC_SPACES("",*DfDx_dot->range(),*stateModel_->get_f_space());
00630   }
00631   if (!is_null(DfDp)) {
00632     THYRA_ASSERT_VEC_SPACES("",*DfDp->domain(),*stateModel_->get_p_space(p_index_));
00633     THYRA_ASSERT_VEC_SPACES("",*DfDp->range(),*stateModel_->get_f_space());
00634   }
00635 #endif
00636 
00637   stateBasePoint_ = stateBasePoint;
00638 
00639   // Set whatever derivatives where passed in.  If an input in null, then the
00640   // member will be null and the null linear operators will be computed later
00641   // just in time.
00642 
00643   W_tilde_ = W_tilde;
00644   coeff_x_dot_ = coeff_x_dot;
00645   coeff_x_ = coeff_x;
00646   DfDx_dot_ = DfDx_dot;
00647   DfDp_ = DfDp;
00648 
00649   wrapNominalValuesAndBounds();
00650 
00651 }
00652 
00653 
00654 // Public functions overridden from ModelEvaulator
00655 
00656 
00657 template<class Scalar>
00658 RCP<const Thyra::VectorSpaceBase<Scalar> >
00659 ForwardSensitivityModelEvaluator<Scalar>::get_x_space() const
00660 {
00661   return s_bar_space_;
00662 }
00663 
00664 
00665 template<class Scalar>
00666 RCP<const Thyra::VectorSpaceBase<Scalar> >
00667 ForwardSensitivityModelEvaluator<Scalar>::get_f_space() const
00668 {
00669   return f_sens_space_;
00670 }
00671 
00672 
00673 template<class Scalar>
00674 Thyra::ModelEvaluatorBase::InArgs<Scalar>
00675 ForwardSensitivityModelEvaluator<Scalar>::getNominalValues() const
00676 {
00677   return nominalValues_;
00678 }
00679 
00680 
00681 template<class Scalar>
00682 RCP<Thyra::LinearOpWithSolveBase<Scalar> >
00683 ForwardSensitivityModelEvaluator<Scalar>::create_W() const
00684 {
00685   return Thyra::multiVectorLinearOpWithSolve<Scalar>();
00686 }
00687 
00688 
00689 template<class Scalar>
00690 Thyra::ModelEvaluatorBase::InArgs<Scalar>
00691 ForwardSensitivityModelEvaluator<Scalar>::createInArgs() const
00692 {
00693   typedef Thyra::ModelEvaluatorBase MEB;
00694   MEB::InArgs<Scalar> stateModelInArgs = stateModel_->createInArgs();
00695   MEB::InArgsSetup<Scalar> inArgs;
00696   inArgs.setModelEvalDescription(this->description());
00697   inArgs.setSupports( MEB::IN_ARG_x_dot,
00698     stateModelInArgs.supports(MEB::IN_ARG_x_dot) );
00699   inArgs.setSupports( MEB::IN_ARG_x );
00700   inArgs.setSupports( MEB::IN_ARG_t );
00701   inArgs.setSupports( MEB::IN_ARG_alpha,
00702     stateModelInArgs.supports(MEB::IN_ARG_alpha) );
00703   inArgs.setSupports( MEB::IN_ARG_beta,
00704     stateModelInArgs.supports(MEB::IN_ARG_beta) );
00705   return inArgs;
00706 }
00707 
00708 
00709 // Private functions overridden from ModelEvaulatorDefaultBase
00710 
00711 
00712 template<class Scalar>
00713 Thyra::ModelEvaluatorBase::OutArgs<Scalar>
00714 ForwardSensitivityModelEvaluator<Scalar>::createOutArgsImpl() const
00715 {
00716 
00717   typedef Thyra::ModelEvaluatorBase MEB;
00718 
00719   MEB::OutArgs<Scalar> stateModelOutArgs = stateModel_->createOutArgs();
00720   MEB::OutArgsSetup<Scalar> outArgs;
00721 
00722   outArgs.setModelEvalDescription(this->description());
00723 
00724   outArgs.setSupports(MEB::OUT_ARG_f);
00725 
00726   if (stateModelOutArgs.supports(MEB::OUT_ARG_W) ) {
00727     outArgs.setSupports(MEB::OUT_ARG_W);
00728     outArgs.set_W_properties(stateModelOutArgs.get_W_properties());
00729   }
00730 
00731   return outArgs;
00732 
00733 }
00734 
00735 
00736 template<class Scalar>
00737 void ForwardSensitivityModelEvaluator<Scalar>::evalModelImpl(
00738   const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
00739   const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs
00740   ) const
00741 {
00742 
00743   using Teuchos::rcp_dynamic_cast;
00744   typedef Teuchos::ScalarTraits<Scalar> ST;
00745   typedef Thyra::ModelEvaluatorBase MEB;
00746   typedef Teuchos::VerboseObjectTempState<Thyra::ModelEvaluatorBase> VOTSME;
00747 
00748   THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_GEN_BEGIN(
00749     "ForwardSensitivityModelEvaluator", inArgs, outArgs, Teuchos::null );
00750 
00751   //
00752   // Update the derivative matrices if they are not already updated for the
00753   // given time!.
00754   //
00755   
00756   {
00757 #ifdef ENABLE_RYTHMOS_TIMERS
00758     TEUCHOS_FUNC_TIME_MONITOR("Rythmos:ForwardSensitivityModelEvaluator::evalModel: computeMatrices");
00759 #endif
00760     computeDerivativeMatrices(inArgs);
00761   }
00762 
00763   //
00764   // InArgs
00765   //
00766 
00767   RCP<const Thyra::DefaultMultiVectorProductVector<Scalar> >
00768     s_bar = rcp_dynamic_cast<const Thyra::DefaultMultiVectorProductVector<Scalar> >(
00769       inArgs.get_x().assert_not_null(), true
00770       );
00771   RCP<const Thyra::DefaultMultiVectorProductVector<Scalar> >
00772     s_bar_dot = rcp_dynamic_cast<const Thyra::DefaultMultiVectorProductVector<Scalar> >(
00773       inArgs.get_x_dot().assert_not_null(), true
00774       );
00775   const Scalar
00776     alpha = inArgs.get_alpha();
00777   const Scalar
00778     beta = inArgs.get_beta();
00779 
00780   RCP<const Thyra::MultiVectorBase<Scalar> >
00781     S = s_bar->getMultiVector();
00782   RCP<const Thyra::MultiVectorBase<Scalar> >
00783     S_dot = s_bar_dot->getMultiVector();
00784   
00785   //
00786   // OutArgs
00787   //
00788 
00789   RCP<Thyra::DefaultMultiVectorProductVector<Scalar> >
00790     f_sens = rcp_dynamic_cast<Thyra::DefaultMultiVectorProductVector<Scalar> >(
00791       outArgs.get_f(), true
00792       );
00793   RCP<Thyra::DefaultMultiVectorLinearOpWithSolve<Scalar> >
00794     W_sens = rcp_dynamic_cast<Thyra::DefaultMultiVectorLinearOpWithSolve<Scalar> >(
00795       outArgs.get_W(), true
00796       );
00797 
00798   RCP<Thyra::MultiVectorBase<Scalar> >
00799     F_sens = f_sens->getNonconstMultiVector().assert_not_null();
00800 
00801   //
00802   // Compute the requested functions
00803   //
00804 
00805   if(!is_null(F_sens)) {
00806 
00807 #ifdef ENABLE_RYTHMOS_TIMERS
00808     TEUCHOS_FUNC_TIME_MONITOR("Rythmos:ForwardSensitivityModelEvaluator::evalModel: computeSens");
00809 #endif
00810 
00811     // S_diff =  -(coeff_x_dot/coeff_x)*S + S_dot
00812     RCP<Thyra::MultiVectorBase<Scalar> >
00813       S_diff = createMembers( stateModel_->get_x_space(), np_ );
00814     V_StVpV( &*S_diff, Scalar(-coeff_x_dot_/coeff_x_), *S, *S_dot );
00815     // F_sens = (1/coeff_x) * W_tilde * S
00816     Thyra::apply(
00817       *W_tilde_, Thyra::NOTRANS,
00818       *S, &*F_sens,
00819       Scalar(1.0/coeff_x_), ST::zero()
00820       );
00821     // F_sens += d(f)/d(x_dot) * S_diff
00822     Thyra::apply(
00823       *DfDx_dot_, Thyra::NOTRANS,
00824       *S_diff, &*F_sens,
00825       ST::one(), ST::one()
00826       );
00827     // F_sens += d(f)/d(p)
00828     Vp_V( &*F_sens, *DfDp_ );
00829   }
00830   
00831   if(!is_null(W_sens)) {
00832     TEST_FOR_EXCEPTION(
00833       alpha != coeff_x_dot_, std::logic_error,
00834       "Error, alpha="<<alpha<<" != coeff_x_dot="<<coeff_x_dot_
00835       <<" with difference = "<<(alpha-coeff_x_dot_)<<"!" );
00836     TEST_FOR_EXCEPTION(
00837       beta != coeff_x_, std::logic_error,
00838       "Error, beta="<<beta<<" != coeff_x="<<coeff_x_
00839       <<" with difference = "<<(beta-coeff_x_)<<"!" );
00840     W_sens->initialize( W_tilde_, s_bar_space_, f_sens_space_ );
00841     
00842   }
00843 
00844   THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_END();
00845 
00846 }
00847 
00848 
00849 // private
00850 
00851 
00852 template<class Scalar>
00853 void ForwardSensitivityModelEvaluator<Scalar>::wrapNominalValuesAndBounds()
00854 {
00855 
00856   using Teuchos::rcp_dynamic_cast;
00857   typedef Thyra::ModelEvaluatorBase MEB;
00858 
00859   // nominalValues_.clear(); // ToDo: Implement this!
00860 
00861   nominalValues_.set_t(stateModel_->getNominalValues().get_t());
00862   
00863   // 2007/05/22: rabartl: Note: Currently there is not much of a reason to set
00864   // an initial condition here since the initial condition for the
00865   // sensitivities is really being set in the ForwardSensitivityStepper
00866   // object!  In the future, a more general use of this class might benefit
00867   // from setting the initial condition here.
00868 
00869 }
00870 
00871 
00872 template<class Scalar>
00873 void ForwardSensitivityModelEvaluator<Scalar>::computeDerivativeMatrices(
00874   const Thyra::ModelEvaluatorBase::InArgs<Scalar> &point
00875   ) const
00876 {
00877 
00878   typedef Thyra::ModelEvaluatorBase MEB;
00879   typedef Teuchos::VerboseObjectTempState<MEB> VOTSME;
00880 
00881   Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
00882   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00883 
00884   const Scalar
00885     t_base = stateBasePoint_.get_t(),
00886     t = point.get_t();
00887 
00888   TEUCHOS_ASSERT_EQUALITY( t , t_base );
00889 
00890   if (is_null(W_tilde_)) {
00891     TEST_FOR_EXCEPT_MSG(true, "ToDo: compute W_tilde from scratch!");
00892   }
00893   
00894   if ( is_null(DfDx_dot_) || is_null(DfDp_) ) {
00895 
00896     MEB::InArgs<Scalar> inArgs = stateBasePoint_;
00897     MEB::OutArgs<Scalar> outArgs = stateModel_->createOutArgs();
00898     
00899     Teuchos::RCP<Thyra::LinearOpBase<Scalar> > DfDx_dot_compute;
00900     if (is_null(DfDx_dot_)) {
00901       DfDx_dot_compute = stateModel_->create_W_op();
00902       inArgs.set_alpha(1.0);
00903       inArgs.set_beta(0.0);
00904       outArgs.set_W_op(DfDx_dot_compute);
00905     }
00906 
00907     Teuchos::RCP<Thyra::MultiVectorBase<Scalar> > DfDp_compute;
00908     if (is_null(DfDp_)) {
00909       DfDp_compute = Thyra::create_DfDp_mv(
00910         *stateModel_,p_index_,
00911         MEB::DERIV_MV_BY_COL
00912         ).getMultiVector();
00913       outArgs.set_DfDp(
00914         p_index_,
00915         MEB::Derivative<Scalar>(DfDp_compute,MEB::DERIV_MV_BY_COL)
00916         );
00917     }
00918     
00919     VOTSME stateModel_outputTempState(stateModel_,out,verbLevel);
00920     stateModel_->evalModel(inArgs,outArgs);
00921     if (!is_null(DfDx_dot_compute))
00922       DfDx_dot_ = DfDx_dot_compute;
00923     if (!is_null(DfDp_compute))
00924       DfDp_ = DfDp_compute;
00925   
00926   }
00927 
00928   TEST_FOR_EXCEPT_MSG( !is_null(stateIntegrator_),
00929     "ToDo: Update for using the stateIntegrator!" );
00930 
00931 /* 2007/12/11: rabartl: ToDo: Update the code below to work for the general
00932  * case.  It does not work in its current form!
00933  */
00934 
00935 /*
00936 
00937   typedef Thyra::ModelEvaluatorBase MEB;
00938   typedef Teuchos::VerboseObjectTempState<MEB> VOTSME;
00939 
00940   RCP<Teuchos::FancyOStream> out = this->getOStream();
00941   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00942   
00943   const Scalar t = point.get_t();
00944 
00945   //
00946   // A) Set up the base point at time t and determine what needs to be
00947   // computed here.
00948   //
00949 
00950   bool update_W_tilde = false;
00951   bool update_DfDx_dot = false;
00952   bool update_DfDp = false;
00953 
00954   if (!is_null(stateIntegrator_)) {
00955     // Get x and x_dot at t and flag to recompute all matrices!
00956     RCP<const Thyra::VectorBase<Scalar> > x, x_dot;
00957     get_fwd_x_and_x_dot( *stateIntegrator_, t, &x, &x_dot );
00958     stateBasePoint_.set_t(t);
00959     stateBasePoint_.set_x(x);
00960     stateBasePoint_.set_x_dot(x_dot);
00961     update_W_tilde = true;
00962     update_DfDx_dot = true;
00963     update_DfDp = true;
00964   }
00965   else {
00966     // The time point should be the same so only compute matrices that have
00967     // not already been computed.
00968     TEUCHOS_ASSERT_EQUALITY( t , stateBasePoint_.get_t() );
00969     if (is_null(W_tilde_))
00970       update_W_tilde = true;
00971     if (is_null(DfDx_dot_))
00972       update_DfDx_dot = true;
00973     if (is_null(DfDp_))
00974       update_DfDx_dot = true;
00975   }
00976 
00977   //
00978   // B) Compute DfDx_dot and DfDp at the same time if needed
00979   //
00980 
00981   if ( update_DfDx_dot || update_DfDp ) {
00982 
00983     // B.1) Allocate storage for matrices.  Note: All of these matrices are
00984     // needed so any of these that is null must be coputed!
00985 
00986     if ( is_null(DfDx_dot_) || is_null(DfDx_dot_compute_) )
00987       DfDx_dot_compute_ = stateModel_->create_W_op();
00988 
00989     if ( is_null(DfDp_) || is_null(DfDp_compute_) )
00990       DfDp_compute_ = Thyra::create_DfDp_mv(
00991         *stateModel_,p_index_,
00992         MEB::DERIV_MV_BY_COL
00993         ).getMultiVector();
00994 
00995     // B.2) Setup the InArgs and OutArgs
00996 
00997     MEB::InArgs<Scalar> inArgs = stateBasePoint_;
00998     MEB::OutArgs<Scalar> outArgs = stateModel_->createOutArgs();
00999     
01000     if (update_DfDx_dot) {
01001       inArgs.set_alpha(1.0);
01002       inArgs.set_beta(0.0);
01003       outArgs.set_W_op(DfDx_dot_compute_);
01004     }
01005     // 2007/12/07: rabartl: ToDo: I need to change the structure of the
01006     // derivative properties in OutArgs to keep track of whether DfDx_dot is
01007     // constant or not separate from W.  Right now I can't do this.  This is
01008     // one reason to add a new DfDx_dot output object to OutArgs.  A better
01009     // solution is a more general data structure giving the dependence of f()
01010     // and g[j]() on x, x_dot, p[l], t, etc ...
01011 
01012     if (update_DfDp) {
01013       outArgs.set_DfDp(
01014         p_index_,
01015         MEB::Derivative<Scalar>(DfDp_compute_, MEB::DERIV_MV_BY_COL)
01016         );
01017     }
01018 
01019     // B.3) Evaluate the outputs
01020     
01021     VOTSME stateModel_outputTempState(stateModel_,out,verbLevel);
01022     stateModel_->evalModel(inArgs,outArgs);
01023 
01024     // B.4) Set the outputs
01025 
01026     if (!is_null(DfDx_dot_compute_))
01027       DfDx_dot_ = DfDx_dot_compute_;
01028 
01029     if (!is_null(DfDp_compute_))
01030       DfDp_ = DfDp_compute_;
01031   
01032   }
01033 
01034   //
01035   // C) Compute W_tilde separately if needed
01036   //
01037 
01038   if ( update_W_tilde ) {
01039 
01040     // C.1) Allocate storage for matrices.  Note: All of these matrices are
01041     // needed so any of these that is null must be coputed!
01042 
01043     if ( is_null(W_tilde_) || is_null(W_tilde_compute_) )
01044       W_tilde_compute_ = stateModel_->create_W();
01045 
01046     // C.2) Setup the InArgs and OutArgs
01047 
01048     MEB::InArgs<Scalar> inArgs = stateBasePoint_;
01049     MEB::OutArgs<Scalar> outArgs = stateModel_->createOutArgs();
01050     
01051     if (is_null(W_tilde_)) {
01052       coeff_x_dot_ = point.get_alpha();
01053       coeff_x_ = point.get_beta();
01054       inArgs.set_alpha(coeff_x_dot_);
01055       inArgs.set_beta(coeff_x_);
01056       outArgs.set_W(W_tilde_compute_);
01057     }
01058 
01059     // C.3) Evaluate the outputs
01060     
01061     VOTSME stateModel_outputTempState(stateModel_,out,verbLevel);
01062     stateModel_->evalModel(inArgs,outArgs);
01063 
01064     // B.4) Set the outputs
01065 
01066     if (!is_null(W_tilde_compute_))
01067       W_tilde_ = W_tilde_compute_;
01068 
01069   }
01070 
01071   // 2007/12/07: rabartl: Note: Above, we see an example of where it would be
01072   // good to have a separate OutArg for DfDx_dot (as a LOB object) so that we
01073   // can compute W and DfDx_dot (and any other output quantitiy) at the same
01074   // time.
01075 
01076 */
01077 
01078 
01079 }
01080 
01081 
01082 } // namespace Rythmos
01083 
01084 
01085 #endif // RYTHMOS_FORWARD_SENSITIVITY_MODEL_EVALUATOR_HPP

Generated on Wed May 12 21:25:43 2010 for Rythmos - Transient Integration for Differential Equations by  doxygen 1.4.7