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

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