Rythmos - Transient Integration for Differential Equations Version of the Day
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 "Rythmos_SolverAcceptingStepperBase.hpp"
00036 #include "Rythmos_SingleResidualModelEvaluator.hpp"
00037 #include "Thyra_ModelEvaluator.hpp" // Interface
00038 #include "Thyra_StateFuncModelEvaluatorBase.hpp" // Implementation
00039 #include "Thyra_DefaultProductVectorSpace.hpp"
00040 #include "Thyra_PhysicallyBlockedLinearOpWithSolveBase.hpp" // Interface
00041 #include "Thyra_DefaultBlockedTriangularLinearOpWithSolve.hpp" // Implementation
00042 #include "Thyra_ModelEvaluatorDelegatorBase.hpp"
00043 #include "Thyra_ModelEvaluatorHelpers.hpp"
00044 #include "Thyra_DefaultMultiVectorProductVectorSpace.hpp"
00045 #include "Thyra_DefaultMultiVectorProductVector.hpp"
00046 #include "Thyra_DefaultMultiVectorLinearOpWithSolve.hpp"
00047 #include "Thyra_MultiVectorStdOps.hpp"
00048 #include "Teuchos_implicit_cast.hpp"
00049 #include "Teuchos_Assert.hpp"
00050 
00051 
00052 namespace Rythmos {
00053 
00054 
00301 template<class Scalar>
00302 class ForwardSensitivityImplicitModelEvaluator
00303   : virtual public ForwardSensitivityModelEvaluatorBase<Scalar>
00304 {
00305 public:
00306 
00309 
00311   ForwardSensitivityImplicitModelEvaluator();
00312 
00328   void initializeStructure(
00329     const RCP<const Thyra::ModelEvaluator<Scalar> > &stateModel,
00330     const int p_index
00331     );
00332   
00334   RCP<const Thyra::ModelEvaluator<Scalar> >
00335   getStateModel() const;
00336 
00338   RCP<Thyra::ModelEvaluator<Scalar> >
00339   getNonconstStateModel() const;
00340   
00342   int get_p_index() const;
00343 
00359   void setStateIntegrator(
00360     const RCP<IntegratorBase<Scalar> > &stateIntegrator,
00361     const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint
00362     );
00363 
00403   void initializePointState(
00404       Ptr<StepperBase<Scalar> > stateStepper,
00405       bool forceUpToDateW
00406       );
00407 
00408 //  /** \brief Initialize full state for a single point in time.
00409 //   *
00410 //   * \param stateBasePoint [in] The base point
00411 //   * <tt>(x_dot_tilde,x_tilde,t_tilde,...)</tt> for which the sensitivities
00412 //   * will be computed and represented at time <tt>t_tilde</tt>.  Any
00413 //   * sensitivities that are needed must be computed at this same time point.
00414 //   * The value of <tt>alpha</tt> and <tt>beta</tt> here are ignored.
00415 //   *
00416 //   * \param W_tilde [in,persisting] The composite state derivative matrix
00417 //   * computed at the base point <tt>stateBasePoint</tt> with
00418 //   * <tt>alpha=coeff_x_dot</tt> and <tt>beta=coeff_x</tt>.  This argument can
00419 //   * be null, in which case it can be computed internally if needed or not at
00420 //   * all.
00421 //   *
00422 //   * \param coeff_x_dot [in] The value of <tt>alpha</tt> for which
00423 //   * <tt>W_tilde</tt> was computed.
00424 //   *
00425 //   * \param coeff_x [in] The value of <tt>beta</tt> for which <tt>W_tilde</tt>
00426 //   * was computed.
00427 //   *
00428 //   * \param DfDx_dot [in] The matrix <tt>d(f)/d(x_dot)</tt> computed at
00429 //   * <tt>stateBasePoint</tt> if available.  If this argument is null, then it
00430 //   * will be computed internally if needed.  The default value is
00431 //   * <tt>Teuchos::null</tt>.
00432 //   *
00433 //   * \param DfDp [in] The matrix <tt>d(f)/d(p(p_index))</tt> computed at
00434 //   * <tt>stateBasePoint</tt> if available.  If this argument is null, then it
00435 //   * will be computed internally if needed.  The default value is
00436 //   * <tt>Teuchos::null</tt>.
00437 //   *
00438 //   * <b>Preconditions:</b><ul>
00439 //   * <li><tt>!is_null(W_tilde)</tt>
00440 //   * </ul>
00441 //   *
00442 //   * This function must be called after <tt>intializeStructure()</tt> and
00443 //   * before <tt>evalModel()</tt> is called.  After this function is called,
00444 //   * <tt>*this</tt> model is fully initialized and ready to compute the
00445 //   * requested outputs.
00446 //   */
00447 //  void initializePointState(
00448 //    const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint,
00449 //    const RCP<const Thyra::LinearOpWithSolveBase<Scalar> > &W_tilde,
00450 //    const Scalar &coeff_x_dot,
00451 //    const Scalar &coeff_x,
00452 //    const RCP<const Thyra::LinearOpBase<Scalar> > &DfDx_dot = Teuchos::null,
00453 //    const RCP<const Thyra::MultiVectorBase<Scalar> > &DfDp = Teuchos::null
00454 //    );
00455 
00456   // 2007/05/22: rabartl: ToDo: Add an InterpolationBufferBase
00457   // stateInterpBuffer object to the initailizeState(...) function that can be
00458   // used to get x and x_dot at different points in time t.  Then, modify the
00459   // logic to recompute all of the needed matrices if t != t_base (as passed
00460   // in through stateBasePoint).  The values of x(t) and xdot(t) can then be
00461   // gotten from the stateInterpBuffer object!
00462   
00464 
00467 
00469   void initializeStructureInitCondOnly(
00470     const RCP<const Thyra::ModelEvaluator<Scalar> >& stateModel,
00471     const RCP<const Thyra::VectorSpaceBase<Scalar> >& p_space
00472     );
00473   
00475   RCP<const Thyra::DefaultMultiVectorProductVectorSpace<Scalar> >
00476   get_s_bar_space() const;
00477   
00479   RCP<const Thyra::VectorSpaceBase<Scalar> > get_p_sens_space() const;
00480 
00482 
00485 
00487   RCP<const Thyra::VectorSpaceBase<Scalar> > get_x_space() const;
00489   RCP<const Thyra::VectorSpaceBase<Scalar> > get_f_space() const;
00491   Thyra::ModelEvaluatorBase::InArgs<Scalar> getNominalValues() const;
00493   RCP<Thyra::LinearOpWithSolveBase<Scalar> > create_W() const;
00495   Thyra::ModelEvaluatorBase::InArgs<Scalar> createInArgs() const;
00496 
00498 
00501 
00503   void initializeState(
00504     const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint,
00505     const RCP<const Thyra::LinearOpWithSolveBase<Scalar> > &W_tilde,
00506     const Scalar &coeff_x_dot,
00507     const Scalar &coeff_x,
00508     const RCP<const Thyra::LinearOpBase<Scalar> > &DfDx_dot = Teuchos::null,
00509     const RCP<const Thyra::MultiVectorBase<Scalar> > &DfDp = Teuchos::null
00510     );
00511 
00513 
00514 private:
00515 
00518 
00520   Thyra::ModelEvaluatorBase::OutArgs<Scalar> createOutArgsImpl() const;
00522   void evalModelImpl(
00523     const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
00524     const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs
00525     ) const;
00526 
00528 
00529 private:
00530 
00531   // /////////////////////////
00532   // Private data members
00533 
00534   RCP<const Thyra::ModelEvaluator<Scalar> > stateModel_;
00535   int p_index_;
00536   RCP<const Thyra::VectorSpaceBase<Scalar> > p_space_;
00537   int np_;
00538 
00539   RCP<IntegratorBase<Scalar> > stateIntegrator_;
00540 
00541   RCP<const Thyra::DefaultMultiVectorProductVectorSpace<Scalar> > s_bar_space_;
00542   RCP<const Thyra::DefaultMultiVectorProductVectorSpace<Scalar> > f_sens_space_;
00543 
00544   Thyra::ModelEvaluatorBase::InArgs<Scalar> nominalValues_;
00545 
00546   mutable Thyra::ModelEvaluatorBase::InArgs<Scalar> stateBasePoint_;
00547 
00548   mutable RCP<const Thyra::LinearOpWithSolveBase<Scalar> > W_tilde_;
00549   mutable Scalar coeff_x_dot_;
00550   mutable Scalar coeff_x_;
00551   mutable RCP<const Thyra::LinearOpBase<Scalar> > DfDx_dot_;
00552   mutable RCP<const Thyra::MultiVectorBase<Scalar> > DfDp_;
00553 
00554   mutable RCP<Thyra::LinearOpWithSolveBase<Scalar> > W_tilde_compute_;
00555   mutable RCP<Thyra::LinearOpBase<Scalar> > DfDx_dot_compute_;
00556   mutable RCP<Thyra::MultiVectorBase<Scalar> > DfDp_compute_;
00557 
00558   // /////////////////////////
00559   // Private member functions
00560 
00561   bool hasStateFuncParams() const { return p_index_ >= 0; }
00562   
00563   void initializeStructureCommon(
00564     const RCP<const Thyra::ModelEvaluator<Scalar> > &stateModel,
00565     const int p_index,
00566     const RCP<const Thyra::VectorSpaceBase<Scalar> > &p_space
00567     );
00568 
00569   void wrapNominalValuesAndBounds();
00570 
00571   void computeDerivativeMatrices(
00572     const Thyra::ModelEvaluatorBase::InArgs<Scalar> &point
00573     ) const;
00574   
00575 };
00576 
00577 
00578 // /////////////////////////////////
00579 // Implementations
00580 
00581 
00582 // Constructors/Intializers/Accessors
00583 
00584 template<class Scalar>
00585 RCP<ForwardSensitivityImplicitModelEvaluator<Scalar> > forwardSensitivityImplicitModelEvaluator()
00586 {
00587   return rcp(new ForwardSensitivityImplicitModelEvaluator<Scalar>());
00588 }
00589 
00590 template<class Scalar>
00591 ForwardSensitivityImplicitModelEvaluator<Scalar>::ForwardSensitivityImplicitModelEvaluator()
00592   : p_index_(0), np_(-1)
00593 {}
00594 
00595 
00596 template<class Scalar>
00597 RCP<const Thyra::ModelEvaluator<Scalar> >
00598 ForwardSensitivityImplicitModelEvaluator<Scalar>::getStateModel() const
00599 {
00600   return stateModel_;
00601 }
00602 
00603 
00604 template<class Scalar>
00605 RCP<Thyra::ModelEvaluator<Scalar> >
00606 ForwardSensitivityImplicitModelEvaluator<Scalar>::getNonconstStateModel() const
00607 {
00608   return Teuchos::null;
00609 }
00610 
00611 
00612 template<class Scalar>
00613 int ForwardSensitivityImplicitModelEvaluator<Scalar>::get_p_index() const
00614 {
00615   return p_index_;
00616 }
00617 
00618 
00619 template<class Scalar>
00620 void ForwardSensitivityImplicitModelEvaluator<Scalar>::setStateIntegrator(
00621   const RCP<IntegratorBase<Scalar> > &stateIntegrator,
00622   const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint
00623   )
00624 {
00625   stateIntegrator_ = stateIntegrator.assert_not_null();
00626   stateBasePoint_ = stateBasePoint;
00627 }
00628 
00629 template<class Scalar>
00630 void ForwardSensitivityImplicitModelEvaluator<Scalar>::initializePointState(
00631       Ptr<StepperBase<Scalar> > stateStepper,
00632       bool forceUpToDateW
00633       )
00634 {
00635   TEUCHOS_ASSERT( Teuchos::nonnull(stateStepper) );
00636   Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
00637   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00638 
00639   const StepStatus<Scalar> stateStepStatus = stateStepper->getStepStatus();
00640   TEST_FOR_EXCEPTION(
00641       stateStepStatus.stepStatus != STEP_STATUS_CONVERGED, std::logic_error,
00642       "Error, the status should be converged since a positive step size was returned!"
00643       );
00644   Scalar curr_t = stateStepStatus.time;
00645 
00646   // Get both x and x_dot since these are needed compute other derivative
00647   // objects at these points.
00648   RCP<const Thyra::VectorBase<Scalar> > x, x_dot;
00649   get_x_and_x_dot(*stateStepper,curr_t,&x,&x_dot);
00650       
00651   stateBasePoint_ = stateStepper->getInitialCondition();
00652   stateBasePoint_.set_x_dot( x_dot );
00653   stateBasePoint_.set_x( x );
00654   stateBasePoint_.set_t( curr_t );
00655 
00656   // Grab the SingleResidualModel that was used to compute the state timestep.
00657   // From this, we can get the constants that where used to compute W!
00658   RCP<SolverAcceptingStepperBase<Scalar> > 
00659     sasStepper = Teuchos::rcp_dynamic_cast<SolverAcceptingStepperBase<Scalar> >(
00660         Teuchos::rcpFromRef(*stateStepper),true
00661         );
00662   RCP<Thyra::NonlinearSolverBase<Scalar> > 
00663     stateTimeStepSolver = sasStepper->getNonconstSolver();
00664   RCP<const SingleResidualModelEvaluatorBase<Scalar> >
00665     singleResidualModel
00666     = Teuchos::rcp_dynamic_cast<const SingleResidualModelEvaluatorBase<Scalar> >(
00667         stateTimeStepSolver->getModel(), true
00668         );
00669   const Scalar
00670     coeff_x_dot = singleResidualModel->get_coeff_x_dot(),
00671     coeff_x = singleResidualModel->get_coeff_x();
00672 
00673   // Get W (and force an update if not up to date already)
00674 
00675   using Teuchos::as;
00676   if ((as<int>(verbLevel) >= as<int>(Teuchos::VERB_MEDIUM)) && forceUpToDateW) {
00677     *out << "\nForcing an update of W at the converged state timestep ...\n";
00678   }
00679 
00680   RCP<Thyra::LinearOpWithSolveBase<Scalar> >
00681     W_tilde = stateTimeStepSolver->get_nonconst_W(forceUpToDateW);
00682 
00683   TEST_FOR_EXCEPTION(
00684       is_null(W_tilde), std::logic_error,
00685       "Error, the W from the state time step must be non-null!"
00686       );
00687 
00688 #ifdef RYTHMOS_DEBUG
00689   TEST_FOR_EXCEPTION(
00690     is_null(stateModel_), std::logic_error,
00691     "Error, you must call intializeStructure(...) before you call initializeState(...)"
00692     );
00693   TEST_FOR_EXCEPT( is_null(stateBasePoint_.get_x()) );
00694   TEST_FOR_EXCEPT( is_null(stateBasePoint_.get_x_dot()) );
00695   TEST_FOR_EXCEPT( is_null(stateBasePoint_.get_p(p_index_)) );
00696   // What about the other parameter values?  We really can't say anything
00697   // about these and we can't check them.  They can be null just fine.
00698   if (!is_null(W_tilde)) {
00699     THYRA_ASSERT_VEC_SPACES("",*W_tilde->domain(),*stateModel_->get_x_space());
00700     THYRA_ASSERT_VEC_SPACES("",*W_tilde->range(),*stateModel_->get_f_space());
00701   }
00702 #endif
00703 
00704   W_tilde_ = W_tilde;
00705   coeff_x_dot_ = coeff_x_dot;
00706   coeff_x_ = coeff_x;
00707   DfDx_dot_ = Teuchos::null;
00708   DfDp_ = Teuchos::null;
00709 
00710   wrapNominalValuesAndBounds();
00711 
00712 }
00713 
00714 
00715 template<class Scalar>
00716 void ForwardSensitivityImplicitModelEvaluator<Scalar>::initializeState(
00717   const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint,
00718   const RCP<const Thyra::LinearOpWithSolveBase<Scalar> > &W_tilde,
00719   const Scalar &coeff_x_dot,
00720   const Scalar &coeff_x,
00721   const RCP<const Thyra::LinearOpBase<Scalar> > &DfDx_dot,
00722   const RCP<const Thyra::MultiVectorBase<Scalar> > &DfDp
00723   )
00724 {
00725 
00726   typedef Thyra::ModelEvaluatorBase MEB;
00727 
00728 #ifdef RYTHMOS_DEBUG
00729   TEST_FOR_EXCEPTION(
00730     is_null(stateModel_), std::logic_error,
00731     "Error, you must call intializeStructure(...) before you call initializeState(...)"
00732     );
00733   TEST_FOR_EXCEPT( is_null(stateBasePoint.get_x()) );
00734   TEST_FOR_EXCEPT( is_null(stateBasePoint.get_x_dot()) );
00735   if (hasStateFuncParams()) {
00736     TEST_FOR_EXCEPT( is_null(stateBasePoint.get_p(p_index_)) );
00737   }
00738   // What about the other parameter values?  We really can't say anything
00739   // about these and we can't check them.  They can be null just fine.
00740   if (nonnull(W_tilde)) {
00741     THYRA_ASSERT_VEC_SPACES("", *W_tilde->domain(), *stateModel_->get_x_space());
00742     THYRA_ASSERT_VEC_SPACES("", *W_tilde->range(), *stateModel_->get_f_space());
00743   }
00744   if (nonnull(DfDx_dot)) {
00745     THYRA_ASSERT_VEC_SPACES("", *DfDx_dot->domain(), *stateModel_->get_x_space());
00746     THYRA_ASSERT_VEC_SPACES("", *DfDx_dot->range(), *stateModel_->get_f_space());
00747   }
00748   if (nonnull(DfDp)) {
00749     TEUCHOS_ASSERT(hasStateFuncParams());
00750     THYRA_ASSERT_VEC_SPACES("", *DfDp->domain(), *p_space_);
00751     THYRA_ASSERT_VEC_SPACES("", *DfDp->range(), *stateModel_->get_f_space());
00752   }
00753 #endif
00754 
00755   stateBasePoint_ = stateBasePoint;
00756 
00757   // Set whatever derivatives where passed in.  If an input in null, then the
00758   // member will be null and the null linear operators will be computed later
00759   // just in time.
00760 
00761   W_tilde_ = W_tilde;
00762   coeff_x_dot_ = coeff_x_dot;
00763   coeff_x_ = coeff_x;
00764   DfDx_dot_ = DfDx_dot;
00765   DfDp_ = DfDp;
00766 
00767   wrapNominalValuesAndBounds();
00768 
00769 }
00770 
00771 
00772 // Public functions overridden from ForwardSensitivityModelEvaluatorBase
00773 
00774 
00775 template<class Scalar>
00776 void ForwardSensitivityImplicitModelEvaluator<Scalar>::initializeStructure(
00777   const RCP<const Thyra::ModelEvaluator<Scalar> >& stateModel,
00778   const int p_index
00779   )
00780 {
00781   initializeStructureCommon(stateModel, p_index, Teuchos::null);
00782 }
00783 
00784 
00785 template<class Scalar>
00786 void ForwardSensitivityImplicitModelEvaluator<Scalar>::initializeStructureInitCondOnly(
00787   const RCP<const Thyra::ModelEvaluator<Scalar> >& stateModel,
00788   const RCP<const Thyra::VectorSpaceBase<Scalar> >& p_space
00789   )
00790 {
00791   initializeStructureCommon(stateModel, -1, p_space);
00792 }
00793 
00794 
00795 template<class Scalar>
00796 RCP<const Thyra::DefaultMultiVectorProductVectorSpace<Scalar> >
00797 ForwardSensitivityImplicitModelEvaluator<Scalar>::get_s_bar_space() const
00798 {
00799   return s_bar_space_;
00800 }
00801 
00802 
00803 template<class Scalar>
00804 RCP<const Thyra::VectorSpaceBase<Scalar> >
00805 ForwardSensitivityImplicitModelEvaluator<Scalar>::get_p_sens_space() const
00806 {
00807   return p_space_;
00808 }
00809 
00810 
00811 // Public functions overridden from ModelEvaulator
00812 
00813 
00814 template<class Scalar>
00815 RCP<const Thyra::VectorSpaceBase<Scalar> >
00816 ForwardSensitivityImplicitModelEvaluator<Scalar>::get_x_space() const
00817 {
00818   return s_bar_space_;
00819 }
00820 
00821 
00822 template<class Scalar>
00823 RCP<const Thyra::VectorSpaceBase<Scalar> >
00824 ForwardSensitivityImplicitModelEvaluator<Scalar>::get_f_space() const
00825 {
00826   return f_sens_space_;
00827 }
00828 
00829 
00830 template<class Scalar>
00831 Thyra::ModelEvaluatorBase::InArgs<Scalar>
00832 ForwardSensitivityImplicitModelEvaluator<Scalar>::getNominalValues() const
00833 {
00834   return nominalValues_;
00835 }
00836 
00837 
00838 template<class Scalar>
00839 RCP<Thyra::LinearOpWithSolveBase<Scalar> >
00840 ForwardSensitivityImplicitModelEvaluator<Scalar>::create_W() const
00841 {
00842   return Thyra::multiVectorLinearOpWithSolve<Scalar>();
00843 }
00844 
00845 
00846 template<class Scalar>
00847 Thyra::ModelEvaluatorBase::InArgs<Scalar>
00848 ForwardSensitivityImplicitModelEvaluator<Scalar>::createInArgs() const
00849 {
00850   typedef Thyra::ModelEvaluatorBase MEB;
00851   MEB::InArgs<Scalar> stateModelInArgs = stateModel_->createInArgs();
00852   MEB::InArgsSetup<Scalar> inArgs;
00853   inArgs.setModelEvalDescription(this->description());
00854   inArgs.setSupports( MEB::IN_ARG_x_dot,
00855     stateModelInArgs.supports(MEB::IN_ARG_x_dot) );
00856   inArgs.setSupports( MEB::IN_ARG_x );
00857   inArgs.setSupports( MEB::IN_ARG_t );
00858   inArgs.setSupports( MEB::IN_ARG_alpha,
00859     stateModelInArgs.supports(MEB::IN_ARG_alpha) );
00860   inArgs.setSupports( MEB::IN_ARG_beta,
00861     stateModelInArgs.supports(MEB::IN_ARG_beta) );
00862   return inArgs;
00863 }
00864 
00865 
00866 // Private functions overridden from ModelEvaulatorDefaultBase
00867 
00868 
00869 template<class Scalar>
00870 Thyra::ModelEvaluatorBase::OutArgs<Scalar>
00871 ForwardSensitivityImplicitModelEvaluator<Scalar>::createOutArgsImpl() const
00872 {
00873 
00874   typedef Thyra::ModelEvaluatorBase MEB;
00875 
00876   MEB::OutArgs<Scalar> stateModelOutArgs = stateModel_->createOutArgs();
00877   MEB::OutArgsSetup<Scalar> outArgs;
00878 
00879   outArgs.setModelEvalDescription(this->description());
00880 
00881   outArgs.setSupports(MEB::OUT_ARG_f);
00882 
00883   if (stateModelOutArgs.supports(MEB::OUT_ARG_W) ) {
00884     outArgs.setSupports(MEB::OUT_ARG_W);
00885     outArgs.set_W_properties(stateModelOutArgs.get_W_properties());
00886   }
00887 
00888   return outArgs;
00889 
00890 }
00891 
00892 
00893 template<class Scalar>
00894 void ForwardSensitivityImplicitModelEvaluator<Scalar>::evalModelImpl(
00895   const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
00896   const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs
00897   ) const
00898 {
00899 
00900   using Teuchos::rcp_dynamic_cast;
00901   typedef Teuchos::ScalarTraits<Scalar> ST;
00902   typedef Thyra::ModelEvaluatorBase MEB;
00903   typedef Teuchos::VerboseObjectTempState<Thyra::ModelEvaluatorBase> VOTSME;
00904 
00905   THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_GEN_BEGIN(
00906     "ForwardSensitivityImplicitModelEvaluator", inArgs, outArgs, Teuchos::null );
00907 
00908   //
00909   // Update the derivative matrices if they are not already updated for the
00910   // given time!.
00911   //
00912   
00913   {
00914 #ifdef ENABLE_RYTHMOS_TIMERS
00915     TEUCHOS_FUNC_TIME_MONITOR_DIFF(
00916       "Rythmos:ForwardSensitivityImplicitModelEvaluator::evalModel: computeMatrices",
00917       RythmosFSIMEmain
00918       );
00919 #endif
00920     computeDerivativeMatrices(inArgs);
00921   }
00922 
00923   //
00924   // InArgs
00925   //
00926 
00927   RCP<const Thyra::DefaultMultiVectorProductVector<Scalar> >
00928     s_bar = rcp_dynamic_cast<const Thyra::DefaultMultiVectorProductVector<Scalar> >(
00929       inArgs.get_x().assert_not_null(), true
00930       );
00931   RCP<const Thyra::DefaultMultiVectorProductVector<Scalar> >
00932     s_bar_dot = rcp_dynamic_cast<const Thyra::DefaultMultiVectorProductVector<Scalar> >(
00933       inArgs.get_x_dot().assert_not_null(), true
00934       );
00935   const Scalar
00936     alpha = inArgs.get_alpha();
00937   const Scalar
00938     beta = inArgs.get_beta();
00939 
00940   RCP<const Thyra::MultiVectorBase<Scalar> >
00941     S = s_bar->getMultiVector();
00942   RCP<const Thyra::MultiVectorBase<Scalar> >
00943     S_dot = s_bar_dot->getMultiVector();
00944   
00945   //
00946   // OutArgs
00947   //
00948 
00949   RCP<Thyra::DefaultMultiVectorProductVector<Scalar> >
00950     f_sens = rcp_dynamic_cast<Thyra::DefaultMultiVectorProductVector<Scalar> >(
00951       outArgs.get_f(), true
00952       );
00953   RCP<Thyra::DefaultMultiVectorLinearOpWithSolve<Scalar> >
00954     W_sens = rcp_dynamic_cast<Thyra::DefaultMultiVectorLinearOpWithSolve<Scalar> >(
00955       outArgs.get_W(), true
00956       );
00957 
00958   RCP<Thyra::MultiVectorBase<Scalar> >
00959     F_sens = f_sens->getNonconstMultiVector().assert_not_null();
00960 
00961   //
00962   // Compute the requested functions
00963   //
00964 
00965   if(nonnull(F_sens)) {
00966 
00967 #ifdef ENABLE_RYTHMOS_TIMERS
00968     TEUCHOS_FUNC_TIME_MONITOR_DIFF(
00969       "Rythmos:ForwardSensitivityImplicitModelEvaluator::evalModel: computeSens",
00970       Rythmos_FSIME);
00971 #endif
00972 
00973     // S_diff =  -(coeff_x_dot/coeff_x)*S + S_dot
00974     RCP<Thyra::MultiVectorBase<Scalar> >
00975       S_diff = createMembers( stateModel_->get_x_space(), np_ );
00976     V_StVpV( &*S_diff, Scalar(-coeff_x_dot_/coeff_x_), *S, *S_dot );
00977     // F_sens = (1/coeff_x) * W_tilde * S
00978     Thyra::apply(
00979       *W_tilde_, Thyra::NOTRANS,
00980       *S, F_sens.ptr(),
00981       Scalar(1.0/coeff_x_), ST::zero()
00982       );
00983     // F_sens += d(f)/d(x_dot) * S_diff
00984     Thyra::apply(
00985       *DfDx_dot_, Thyra::NOTRANS,
00986       *S_diff, F_sens.ptr(),
00987       ST::one(), ST::one()
00988       );
00989     // F_sens += d(f)/d(p)
00990     if (hasStateFuncParams())
00991       Vp_V( &*F_sens, *DfDp_ );
00992   }
00993   
00994   if(nonnull(W_sens)) {
00995     TEST_FOR_EXCEPTION(
00996       alpha != coeff_x_dot_, std::logic_error,
00997       "Error, alpha="<<alpha<<" != coeff_x_dot="<<coeff_x_dot_
00998       <<" with difference = "<<(alpha-coeff_x_dot_)<<"!" );
00999     TEST_FOR_EXCEPTION(
01000       beta != coeff_x_, std::logic_error,
01001       "Error, beta="<<beta<<" != coeff_x="<<coeff_x_
01002       <<" with difference = "<<(beta-coeff_x_)<<"!" );
01003     W_sens->initialize( W_tilde_, s_bar_space_, f_sens_space_ );
01004     
01005   }
01006 
01007   THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_END();
01008 
01009 }
01010 
01011 
01012 // private
01013 
01014 
01015 template<class Scalar>
01016 void ForwardSensitivityImplicitModelEvaluator<Scalar>::initializeStructureCommon(
01017   const RCP<const Thyra::ModelEvaluator<Scalar> >& stateModel,
01018   const int p_index,
01019   const RCP<const Thyra::VectorSpaceBase<Scalar> >& p_space
01020   )
01021 {
01022 
01023   typedef Thyra::ModelEvaluatorBase MEB;
01024 
01025   //
01026   // Validate input
01027   //
01028 
01029   TEUCHOS_ASSERT(nonnull(stateModel));
01030   TEUCHOS_ASSERT(p_index >= 0 || nonnull(p_space));
01031   if (p_index >= 0) {
01032     TEST_FOR_EXCEPTION(
01033       !( 0 <= p_index && p_index < stateModel->Np() ), std::logic_error,
01034       "Error, p_index does not fall in the range [0,"<<(stateModel->Np()-1)<<"]!" );
01035     // ToDo: Validate support for DfDp!
01036   }
01037   else {
01038     TEUCHOS_ASSERT_EQUALITY(p_index, -1);
01039   }
01040 
01041   //
01042   // Set the input objects
01043   //
01044 
01045   stateModel_ = stateModel;
01046 
01047   if (p_index >= 0) {
01048     p_index_ = p_index;
01049     p_space_ = stateModel_->get_p_space(p_index);
01050   }
01051   else {
01052     p_index_ = -1;
01053     p_space_ = p_space;
01054   }
01055 
01056   np_ = p_space_->dim();
01057 
01058   //
01059   // Create the structure of the model
01060   //
01061 
01062   s_bar_space_ = Thyra::multiVectorProductVectorSpace(
01063     stateModel_->get_x_space(), np_
01064     );
01065 
01066   f_sens_space_ = Thyra::multiVectorProductVectorSpace(
01067     stateModel_->get_f_space(), np_
01068     );
01069 
01070   nominalValues_ = this->createInArgs();
01071 
01072   //
01073   // Wipe out matrix storage
01074   //
01075 
01076   stateBasePoint_ = MEB::InArgs<Scalar>();
01077   W_tilde_ = Teuchos::null;
01078   coeff_x_dot_ = 0.0;
01079   coeff_x_ = 0.0;
01080   DfDx_dot_ = Teuchos::null;
01081   DfDp_ = Teuchos::null;
01082   W_tilde_compute_ = Teuchos::null;
01083   DfDx_dot_compute_ = Teuchos::null;
01084   DfDp_compute_ = Teuchos::null;
01085 
01086 }
01087 
01088 
01089 template<class Scalar>
01090 void ForwardSensitivityImplicitModelEvaluator<Scalar>::wrapNominalValuesAndBounds()
01091 {
01092 
01093   using Teuchos::rcp_dynamic_cast;
01094   typedef Thyra::ModelEvaluatorBase MEB;
01095 
01096   // nominalValues_.clear(); // ToDo: Implement this!
01097 
01098   nominalValues_.set_t(stateModel_->getNominalValues().get_t());
01099   
01100   // 2007/05/22: rabartl: Note: Currently there is not much of a reason to set
01101   // an initial condition here since the initial condition for the
01102   // sensitivities is really being set in the ForwardSensitivityStepper
01103   // object!  In the future, a more general use of this class might benefit
01104   // from setting the initial condition here.
01105 
01106 }
01107 
01108 
01109 template<class Scalar>
01110 void ForwardSensitivityImplicitModelEvaluator<Scalar>::computeDerivativeMatrices(
01111   const Thyra::ModelEvaluatorBase::InArgs<Scalar> &point
01112   ) const
01113 {
01114 
01115   typedef Thyra::ModelEvaluatorBase MEB;
01116   typedef Teuchos::VerboseObjectTempState<MEB> VOTSME;
01117 
01118   Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
01119   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
01120 
01121   const Scalar
01122     t_base = stateBasePoint_.get_t(),
01123     t = point.get_t();
01124 
01125   TEUCHOS_ASSERT_EQUALITY( t , t_base );
01126 
01127   if (is_null(W_tilde_)) {
01128     TEST_FOR_EXCEPT_MSG(true, "ToDo: compute W_tilde from scratch!");
01129   }
01130   
01131   if ( is_null(DfDx_dot_) || is_null(DfDp_) ) {
01132 
01133     MEB::InArgs<Scalar> inArgs = stateBasePoint_;
01134     MEB::OutArgs<Scalar> outArgs = stateModel_->createOutArgs();
01135     
01136     Teuchos::RCP<Thyra::LinearOpBase<Scalar> > DfDx_dot_compute;
01137     if (is_null(DfDx_dot_)) {
01138       DfDx_dot_compute = stateModel_->create_W_op();
01139       inArgs.set_alpha(1.0);
01140       inArgs.set_beta(0.0);
01141       outArgs.set_W_op(DfDx_dot_compute);
01142     }
01143 
01144     Teuchos::RCP<Thyra::MultiVectorBase<Scalar> > DfDp_compute;
01145     if (is_null(DfDp_) && hasStateFuncParams()) {
01146       DfDp_compute = Thyra::create_DfDp_mv(
01147         *stateModel_, p_index_,
01148         MEB::DERIV_MV_BY_COL
01149         ).getMultiVector();
01150       outArgs.set_DfDp(
01151         p_index_,
01152         MEB::Derivative<Scalar>(DfDp_compute,MEB::DERIV_MV_BY_COL)
01153         );
01154     }
01155     
01156     VOTSME stateModel_outputTempState(stateModel_,out,verbLevel);
01157     stateModel_->evalModel(inArgs,outArgs);
01158     if (nonnull(DfDx_dot_compute))
01159       DfDx_dot_ = DfDx_dot_compute;
01160     if (nonnull(DfDp_compute))
01161       DfDp_ = DfDp_compute;
01162   
01163   }
01164 
01165   TEST_FOR_EXCEPT_MSG( nonnull(stateIntegrator_),
01166     "ToDo: Update for using the stateIntegrator!" );
01167 
01168 /* 2007/12/11: rabartl: ToDo: Update the code below to work for the general
01169  * case.  It does not work in its current form!
01170  */
01171 
01172 /*
01173 
01174   typedef Thyra::ModelEvaluatorBase MEB;
01175   typedef Teuchos::VerboseObjectTempState<MEB> VOTSME;
01176 
01177   RCP<Teuchos::FancyOStream> out = this->getOStream();
01178   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
01179   
01180   const Scalar t = point.get_t();
01181 
01182   //
01183   // A) Set up the base point at time t and determine what needs to be
01184   // computed here.
01185   //
01186 
01187   bool update_W_tilde = false;
01188   bool update_DfDx_dot = false;
01189   bool update_DfDp = false;
01190 
01191   if (nonnull(stateIntegrator_)) {
01192     // Get x and x_dot at t and flag to recompute all matrices!
01193     RCP<const Thyra::VectorBase<Scalar> > x, x_dot;
01194     get_fwd_x_and_x_dot( *stateIntegrator_, t, &x, &x_dot );
01195     stateBasePoint_.set_t(t);
01196     stateBasePoint_.set_x(x);
01197     stateBasePoint_.set_x_dot(x_dot);
01198     update_W_tilde = true;
01199     update_DfDx_dot = true;
01200     update_DfDp = true;
01201   }
01202   else {
01203     // The time point should be the same so only compute matrices that have
01204     // not already been computed.
01205     TEUCHOS_ASSERT_EQUALITY( t , stateBasePoint_.get_t() );
01206     if (is_null(W_tilde_))
01207       update_W_tilde = true;
01208     if (is_null(DfDx_dot_))
01209       update_DfDx_dot = true;
01210     if (is_null(DfDp_))
01211       update_DfDx_dot = true;
01212   }
01213 
01214   //
01215   // B) Compute DfDx_dot and DfDp at the same time if needed
01216   //
01217 
01218   if ( update_DfDx_dot || update_DfDp ) {
01219 
01220     // B.1) Allocate storage for matrices.  Note: All of these matrices are
01221     // needed so any of these that is null must be coputed!
01222 
01223     if ( is_null(DfDx_dot_) || is_null(DfDx_dot_compute_) )
01224       DfDx_dot_compute_ = stateModel_->create_W_op();
01225 
01226     if ( is_null(DfDp_) || is_null(DfDp_compute_) )
01227       DfDp_compute_ = Thyra::create_DfDp_mv(
01228         *stateModel_,p_index_,
01229         MEB::DERIV_MV_BY_COL
01230         ).getMultiVector();
01231 
01232     // B.2) Setup the InArgs and OutArgs
01233 
01234     MEB::InArgs<Scalar> inArgs = stateBasePoint_;
01235     MEB::OutArgs<Scalar> outArgs = stateModel_->createOutArgs();
01236     
01237     if (update_DfDx_dot) {
01238       inArgs.set_alpha(1.0);
01239       inArgs.set_beta(0.0);
01240       outArgs.set_W_op(DfDx_dot_compute_);
01241     }
01242     // 2007/12/07: rabartl: ToDo: I need to change the structure of the
01243     // derivative properties in OutArgs to keep track of whether DfDx_dot is
01244     // constant or not separate from W.  Right now I can't do this.  This is
01245     // one reason to add a new DfDx_dot output object to OutArgs.  A better
01246     // solution is a more general data structure giving the dependence of f()
01247     // and g[j]() on x, x_dot, p[l], t, etc ...
01248 
01249     if (update_DfDp) {
01250       outArgs.set_DfDp(
01251         p_index_,
01252         MEB::Derivative<Scalar>(DfDp_compute_, MEB::DERIV_MV_BY_COL)
01253         );
01254     }
01255 
01256     // B.3) Evaluate the outputs
01257     
01258     VOTSME stateModel_outputTempState(stateModel_,out,verbLevel);
01259     stateModel_->evalModel(inArgs,outArgs);
01260 
01261     // B.4) Set the outputs
01262 
01263     if (nonnull(DfDx_dot_compute_))
01264       DfDx_dot_ = DfDx_dot_compute_;
01265 
01266     if (nonnull(DfDp_compute_))
01267       DfDp_ = DfDp_compute_;
01268   
01269   }
01270 
01271   //
01272   // C) Compute W_tilde separately if needed
01273   //
01274 
01275   if ( update_W_tilde ) {
01276 
01277     // C.1) Allocate storage for matrices.  Note: All of these matrices are
01278     // needed so any of these that is null must be coputed!
01279 
01280     if ( is_null(W_tilde_) || is_null(W_tilde_compute_) )
01281       W_tilde_compute_ = stateModel_->create_W();
01282 
01283     // C.2) Setup the InArgs and OutArgs
01284 
01285     MEB::InArgs<Scalar> inArgs = stateBasePoint_;
01286     MEB::OutArgs<Scalar> outArgs = stateModel_->createOutArgs();
01287     
01288     if (is_null(W_tilde_)) {
01289       coeff_x_dot_ = point.get_alpha();
01290       coeff_x_ = point.get_beta();
01291       inArgs.set_alpha(coeff_x_dot_);
01292       inArgs.set_beta(coeff_x_);
01293       outArgs.set_W(W_tilde_compute_);
01294     }
01295 
01296     // C.3) Evaluate the outputs
01297     
01298     VOTSME stateModel_outputTempState(stateModel_,out,verbLevel);
01299     stateModel_->evalModel(inArgs,outArgs);
01300 
01301     // B.4) Set the outputs
01302 
01303     if (nonnull(W_tilde_compute_))
01304       W_tilde_ = W_tilde_compute_;
01305 
01306   }
01307 
01308   // 2007/12/07: rabartl: Note: Above, we see an example of where it would be
01309   // good to have a separate OutArg for DfDx_dot (as a LOB object) so that we
01310   // can compute W and DfDx_dot (and any other output quantitiy) at the same
01311   // time.
01312 
01313 */
01314 
01315 
01316 }
01317 
01318 
01319 } // namespace Rythmos
01320 
01321 
01322 #endif // RYTHMOS_FORWARD_SENSITIVITY_IMPLICIT_MODEL_EVALUATOR_HPP
 All Classes Functions Variables Typedefs Friends