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., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
00023 // USA
00024 // Questions? Contact Todd S. Coffey (tscoffe@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 //@HEADER
00028 
00029 #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 initializeState(...) 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   TEUCHOS_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   TEUCHOS_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 HAVE_RYTHMOS_DEBUG
00689   TEUCHOS_TEST_FOR_EXCEPTION(
00690     is_null(stateModel_), std::logic_error,
00691     "Error, you must call intializeStructure(...) before you call initializeState(...)"
00692     );
00693   TEUCHOS_TEST_FOR_EXCEPT( is_null(stateBasePoint_.get_x()) );
00694   TEUCHOS_TEST_FOR_EXCEPT( is_null(stateBasePoint_.get_x_dot()) );
00695   TEUCHOS_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 HAVE_RYTHMOS_DEBUG
00729   TEUCHOS_TEST_FOR_EXCEPTION(
00730     is_null(stateModel_), std::logic_error,
00731     "Error, you must call intializeStructure(...) before you call initializeState(...)"
00732     );
00733   TEUCHOS_TEST_FOR_EXCEPT( is_null(stateBasePoint.get_x()) );
00734   TEUCHOS_TEST_FOR_EXCEPT( is_null(stateBasePoint.get_x_dot()) );
00735   if (hasStateFuncParams()) {
00736     TEUCHOS_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::as;
00901   using Teuchos::rcp_dynamic_cast;
00902   typedef Teuchos::ScalarTraits<Scalar> ST;
00903   typedef Thyra::ModelEvaluatorBase MEB;
00904   typedef Teuchos::VerboseObjectTempState<Thyra::ModelEvaluatorBase> VOTSME;
00905 
00906   THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_GEN_BEGIN(
00907     "ForwardSensitivityImplicitModelEvaluator", inArgs, outArgs, Teuchos::null );
00908 
00909   //
00910   // Update the derivative matrices if they are not already updated for the
00911   // given time!.
00912   //
00913   
00914   {
00915     RYTHMOS_FUNC_TIME_MONITOR_DIFF(
00916       "Rythmos:ForwardSensitivityImplicitModelEvaluator::evalModel: computeMatrices",
00917       RythmosFSIMEmain
00918       );
00919     computeDerivativeMatrices(inArgs);
00920   }
00921 
00922   //
00923   // InArgs
00924   //
00925 
00926   RCP<const Thyra::DefaultMultiVectorProductVector<Scalar> >
00927     s_bar = rcp_dynamic_cast<const Thyra::DefaultMultiVectorProductVector<Scalar> >(
00928       inArgs.get_x().assert_not_null(), true
00929       );
00930   RCP<const Thyra::DefaultMultiVectorProductVector<Scalar> >
00931     s_bar_dot = rcp_dynamic_cast<const Thyra::DefaultMultiVectorProductVector<Scalar> >(
00932       inArgs.get_x_dot().assert_not_null(), true
00933       );
00934   const Scalar
00935     alpha = inArgs.get_alpha();
00936   const Scalar
00937     beta = inArgs.get_beta();
00938 
00939   RCP<const Thyra::MultiVectorBase<Scalar> >
00940     S = s_bar->getMultiVector();
00941   RCP<const Thyra::MultiVectorBase<Scalar> >
00942     S_dot = s_bar_dot->getMultiVector();
00943   
00944   //
00945   // OutArgs
00946   //
00947 
00948   RCP<Thyra::DefaultMultiVectorProductVector<Scalar> >
00949     f_sens = rcp_dynamic_cast<Thyra::DefaultMultiVectorProductVector<Scalar> >(
00950       outArgs.get_f(), true
00951       );
00952   RCP<Thyra::DefaultMultiVectorLinearOpWithSolve<Scalar> >
00953     W_sens = rcp_dynamic_cast<Thyra::DefaultMultiVectorLinearOpWithSolve<Scalar> >(
00954       outArgs.get_W(), true
00955       );
00956 
00957   RCP<Thyra::MultiVectorBase<Scalar> >
00958     F_sens = f_sens->getNonconstMultiVector().assert_not_null();
00959 
00960   //
00961   // Compute the requested functions
00962   //
00963 
00964   if(nonnull(F_sens)) {
00965 
00966     RYTHMOS_FUNC_TIME_MONITOR_DIFF(
00967       "Rythmos:ForwardSensitivityImplicitModelEvaluator::evalModel: computeSens",
00968       Rythmos_FSIME);
00969 
00970     // S_diff =  -(coeff_x_dot/coeff_x)*S + S_dot
00971     RCP<Thyra::MultiVectorBase<Scalar> >
00972       S_diff = createMembers( stateModel_->get_x_space(), np_ );
00973     V_StVpV( S_diff.ptr(), as<Scalar>(-coeff_x_dot_/coeff_x_), *S, *S_dot );
00974     // F_sens = (1/coeff_x) * W_tilde * S
00975     Thyra::apply(
00976       *W_tilde_, Thyra::NOTRANS,
00977       *S, F_sens.ptr(),
00978       as<Scalar>(1.0/coeff_x_), ST::zero()
00979       );
00980     // F_sens += d(f)/d(x_dot) * S_diff
00981     Thyra::apply(
00982       *DfDx_dot_, Thyra::NOTRANS,
00983       *S_diff, F_sens.ptr(),
00984       ST::one(), ST::one()
00985       );
00986     // F_sens += d(f)/d(p)
00987     if (hasStateFuncParams())
00988       Vp_V( F_sens.ptr(), *DfDp_ );
00989   }
00990   
00991   if(nonnull(W_sens)) {
00992     TEUCHOS_TEST_FOR_EXCEPTION(
00993       alpha != coeff_x_dot_, std::logic_error,
00994       "Error, alpha="<<alpha<<" != coeff_x_dot="<<coeff_x_dot_
00995       <<" with difference = "<<(alpha-coeff_x_dot_)<<"!" );
00996     TEUCHOS_TEST_FOR_EXCEPTION(
00997       beta != coeff_x_, std::logic_error,
00998       "Error, beta="<<beta<<" != coeff_x="<<coeff_x_
00999       <<" with difference = "<<(beta-coeff_x_)<<"!" );
01000     W_sens->initialize( W_tilde_, s_bar_space_, f_sens_space_ );
01001     
01002   }
01003 
01004   THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_END();
01005 
01006 }
01007 
01008 
01009 // private
01010 
01011 
01012 template<class Scalar>
01013 void ForwardSensitivityImplicitModelEvaluator<Scalar>::initializeStructureCommon(
01014   const RCP<const Thyra::ModelEvaluator<Scalar> >& stateModel,
01015   const int p_index,
01016   const RCP<const Thyra::VectorSpaceBase<Scalar> >& p_space
01017   )
01018 {
01019 
01020   typedef Thyra::ModelEvaluatorBase MEB;
01021 
01022   //
01023   // Validate input
01024   //
01025 
01026   TEUCHOS_ASSERT(nonnull(stateModel));
01027   TEUCHOS_ASSERT(p_index >= 0 || nonnull(p_space));
01028   if (p_index >= 0) {
01029     TEUCHOS_TEST_FOR_EXCEPTION(
01030       !( 0 <= p_index && p_index < stateModel->Np() ), std::logic_error,
01031       "Error, p_index does not fall in the range [0,"<<(stateModel->Np()-1)<<"]!" );
01032     // ToDo: Validate support for DfDp!
01033   }
01034   else {
01035     TEUCHOS_ASSERT_EQUALITY(p_index, -1);
01036   }
01037 
01038   //
01039   // Set the input objects
01040   //
01041 
01042   stateModel_ = stateModel;
01043 
01044   if (p_index >= 0) {
01045     p_index_ = p_index;
01046     p_space_ = stateModel_->get_p_space(p_index);
01047   }
01048   else {
01049     p_index_ = -1;
01050     p_space_ = p_space;
01051   }
01052 
01053   np_ = p_space_->dim();
01054 
01055   //
01056   // Create the structure of the model
01057   //
01058 
01059   s_bar_space_ = Thyra::multiVectorProductVectorSpace(
01060     stateModel_->get_x_space(), np_
01061     );
01062 
01063   f_sens_space_ = Thyra::multiVectorProductVectorSpace(
01064     stateModel_->get_f_space(), np_
01065     );
01066 
01067   nominalValues_ = this->createInArgs();
01068 
01069   //
01070   // Wipe out matrix storage
01071   //
01072 
01073   stateBasePoint_ = MEB::InArgs<Scalar>();
01074   W_tilde_ = Teuchos::null;
01075   coeff_x_dot_ = 0.0;
01076   coeff_x_ = 0.0;
01077   DfDx_dot_ = Teuchos::null;
01078   DfDp_ = Teuchos::null;
01079   W_tilde_compute_ = Teuchos::null;
01080   DfDx_dot_compute_ = Teuchos::null;
01081   DfDp_compute_ = Teuchos::null;
01082 
01083 }
01084 
01085 
01086 template<class Scalar>
01087 void ForwardSensitivityImplicitModelEvaluator<Scalar>::wrapNominalValuesAndBounds()
01088 {
01089 
01090   using Teuchos::rcp_dynamic_cast;
01091   typedef Thyra::ModelEvaluatorBase MEB;
01092 
01093   // nominalValues_.clear(); // ToDo: Implement this!
01094 
01095   nominalValues_.set_t(stateModel_->getNominalValues().get_t());
01096   
01097   // 2007/05/22: rabartl: Note: Currently there is not much of a reason to set
01098   // an initial condition here since the initial condition for the
01099   // sensitivities is really being set in the ForwardSensitivityStepper
01100   // object!  In the future, a more general use of this class might benefit
01101   // from setting the initial condition here.
01102 
01103 }
01104 
01105 
01106 template<class Scalar>
01107 void ForwardSensitivityImplicitModelEvaluator<Scalar>::computeDerivativeMatrices(
01108   const Thyra::ModelEvaluatorBase::InArgs<Scalar> &point
01109   ) const
01110 {
01111 
01112   typedef Thyra::ModelEvaluatorBase MEB;
01113   typedef Teuchos::VerboseObjectTempState<MEB> VOTSME;
01114 
01115   Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
01116   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
01117 
01118   const Scalar
01119     t_base = stateBasePoint_.get_t(),
01120     t = point.get_t();
01121 
01122   TEUCHOS_ASSERT_EQUALITY( t , t_base );
01123 
01124   if (is_null(W_tilde_)) {
01125     TEUCHOS_TEST_FOR_EXCEPT_MSG(true, "ToDo: compute W_tilde from scratch!");
01126   }
01127   
01128   if ( is_null(DfDx_dot_) || is_null(DfDp_) ) {
01129 
01130     MEB::InArgs<Scalar> inArgs = stateBasePoint_;
01131     MEB::OutArgs<Scalar> outArgs = stateModel_->createOutArgs();
01132     
01133     Teuchos::RCP<Thyra::LinearOpBase<Scalar> > DfDx_dot_compute;
01134     if (is_null(DfDx_dot_)) {
01135       DfDx_dot_compute = stateModel_->create_W_op();
01136       inArgs.set_alpha(1.0);
01137       inArgs.set_beta(0.0);
01138       outArgs.set_W_op(DfDx_dot_compute);
01139     }
01140 
01141     Teuchos::RCP<Thyra::MultiVectorBase<Scalar> > DfDp_compute;
01142     if (is_null(DfDp_) && hasStateFuncParams()) {
01143       DfDp_compute = Thyra::create_DfDp_mv(
01144         *stateModel_, p_index_,
01145         MEB::DERIV_MV_BY_COL
01146         ).getMultiVector();
01147       outArgs.set_DfDp(
01148         p_index_,
01149         MEB::Derivative<Scalar>(DfDp_compute,MEB::DERIV_MV_BY_COL)
01150         );
01151     }
01152     
01153     VOTSME stateModel_outputTempState(stateModel_,out,verbLevel);
01154     stateModel_->evalModel(inArgs,outArgs);
01155     if (nonnull(DfDx_dot_compute))
01156       DfDx_dot_ = DfDx_dot_compute;
01157     if (nonnull(DfDp_compute))
01158       DfDp_ = DfDp_compute;
01159   
01160   }
01161 
01162   TEUCHOS_TEST_FOR_EXCEPT_MSG( nonnull(stateIntegrator_),
01163     "ToDo: Update for using the stateIntegrator!" );
01164 
01165 /* 2007/12/11: rabartl: ToDo: Update the code below to work for the general
01166  * case.  It does not work in its current form!
01167  */
01168 
01169 /*
01170 
01171   typedef Thyra::ModelEvaluatorBase MEB;
01172   typedef Teuchos::VerboseObjectTempState<MEB> VOTSME;
01173 
01174   RCP<Teuchos::FancyOStream> out = this->getOStream();
01175   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
01176   
01177   const Scalar t = point.get_t();
01178 
01179   //
01180   // A) Set up the base point at time t and determine what needs to be
01181   // computed here.
01182   //
01183 
01184   bool update_W_tilde = false;
01185   bool update_DfDx_dot = false;
01186   bool update_DfDp = false;
01187 
01188   if (nonnull(stateIntegrator_)) {
01189     // Get x and x_dot at t and flag to recompute all matrices!
01190     RCP<const Thyra::VectorBase<Scalar> > x, x_dot;
01191     get_fwd_x_and_x_dot( *stateIntegrator_, t, &x, &x_dot );
01192     stateBasePoint_.set_t(t);
01193     stateBasePoint_.set_x(x);
01194     stateBasePoint_.set_x_dot(x_dot);
01195     update_W_tilde = true;
01196     update_DfDx_dot = true;
01197     update_DfDp = true;
01198   }
01199   else {
01200     // The time point should be the same so only compute matrices that have
01201     // not already been computed.
01202     TEUCHOS_ASSERT_EQUALITY( t , stateBasePoint_.get_t() );
01203     if (is_null(W_tilde_))
01204       update_W_tilde = true;
01205     if (is_null(DfDx_dot_))
01206       update_DfDx_dot = true;
01207     if (is_null(DfDp_))
01208       update_DfDx_dot = true;
01209   }
01210 
01211   //
01212   // B) Compute DfDx_dot and DfDp at the same time if needed
01213   //
01214 
01215   if ( update_DfDx_dot || update_DfDp ) {
01216 
01217     // B.1) Allocate storage for matrices.  Note: All of these matrices are
01218     // needed so any of these that is null must be coputed!
01219 
01220     if ( is_null(DfDx_dot_) || is_null(DfDx_dot_compute_) )
01221       DfDx_dot_compute_ = stateModel_->create_W_op();
01222 
01223     if ( is_null(DfDp_) || is_null(DfDp_compute_) )
01224       DfDp_compute_ = Thyra::create_DfDp_mv(
01225         *stateModel_,p_index_,
01226         MEB::DERIV_MV_BY_COL
01227         ).getMultiVector();
01228 
01229     // B.2) Setup the InArgs and OutArgs
01230 
01231     MEB::InArgs<Scalar> inArgs = stateBasePoint_;
01232     MEB::OutArgs<Scalar> outArgs = stateModel_->createOutArgs();
01233     
01234     if (update_DfDx_dot) {
01235       inArgs.set_alpha(1.0);
01236       inArgs.set_beta(0.0);
01237       outArgs.set_W_op(DfDx_dot_compute_);
01238     }
01239     // 2007/12/07: rabartl: ToDo: I need to change the structure of the
01240     // derivative properties in OutArgs to keep track of whether DfDx_dot is
01241     // constant or not separate from W.  Right now I can't do this.  This is
01242     // one reason to add a new DfDx_dot output object to OutArgs.  A better
01243     // solution is a more general data structure giving the dependence of f()
01244     // and g[j]() on x, x_dot, p[l], t, etc ...
01245 
01246     if (update_DfDp) {
01247       outArgs.set_DfDp(
01248         p_index_,
01249         MEB::Derivative<Scalar>(DfDp_compute_, MEB::DERIV_MV_BY_COL)
01250         );
01251     }
01252 
01253     // B.3) Evaluate the outputs
01254     
01255     VOTSME stateModel_outputTempState(stateModel_,out,verbLevel);
01256     stateModel_->evalModel(inArgs,outArgs);
01257 
01258     // B.4) Set the outputs
01259 
01260     if (nonnull(DfDx_dot_compute_))
01261       DfDx_dot_ = DfDx_dot_compute_;
01262 
01263     if (nonnull(DfDp_compute_))
01264       DfDp_ = DfDp_compute_;
01265   
01266   }
01267 
01268   //
01269   // C) Compute W_tilde separately if needed
01270   //
01271 
01272   if ( update_W_tilde ) {
01273 
01274     // C.1) Allocate storage for matrices.  Note: All of these matrices are
01275     // needed so any of these that is null must be coputed!
01276 
01277     if ( is_null(W_tilde_) || is_null(W_tilde_compute_) )
01278       W_tilde_compute_ = stateModel_->create_W();
01279 
01280     // C.2) Setup the InArgs and OutArgs
01281 
01282     MEB::InArgs<Scalar> inArgs = stateBasePoint_;
01283     MEB::OutArgs<Scalar> outArgs = stateModel_->createOutArgs();
01284     
01285     if (is_null(W_tilde_)) {
01286       coeff_x_dot_ = point.get_alpha();
01287       coeff_x_ = point.get_beta();
01288       inArgs.set_alpha(coeff_x_dot_);
01289       inArgs.set_beta(coeff_x_);
01290       outArgs.set_W(W_tilde_compute_);
01291     }
01292 
01293     // C.3) Evaluate the outputs
01294     
01295     VOTSME stateModel_outputTempState(stateModel_,out,verbLevel);
01296     stateModel_->evalModel(inArgs,outArgs);
01297 
01298     // B.4) Set the outputs
01299 
01300     if (nonnull(W_tilde_compute_))
01301       W_tilde_ = W_tilde_compute_;
01302 
01303   }
01304 
01305   // 2007/12/07: rabartl: Note: Above, we see an example of where it would be
01306   // good to have a separate OutArg for DfDx_dot (as a LOB object) so that we
01307   // can compute W and DfDx_dot (and any other output quantitiy) at the same
01308   // time.
01309 
01310 */
01311 
01312 
01313 }
01314 
01315 
01316 } // namespace Rythmos
01317 
01318 
01319 #endif // RYTHMOS_FORWARD_SENSITIVITY_IMPLICIT_MODEL_EVALUATOR_HPP
 All Classes Functions Variables Typedefs Friends