Rythmos_ForwardSensitivityModelEvaluator.hpp

//@HEADER
// ***********************************************************************
//
//                           Rythmos Package
//                 Copyright (2006) Sandia Corporation
//
// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
// license for use of this work by or on behalf of the U.S. Government.
//
// This library is free software; you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as
// published by the Free Software Foundation; either version 2.1 of the
// License, or (at your option) any later version.
//
// This library is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
// USA
// Questions? Contact Todd S. Coffey (tscoffe@sandia.gov)
//
// ***********************************************************************
//@HEADER

#ifndef RYTHMOS_FORWARD_SENSITIVITY_MODEL_EVALUATOR_HPP
#define RYTHMOS_FORWARD_SENSITIVITY_MODEL_EVALUATOR_HPP


#include "Thyra_ModelEvaluator.hpp" // Interface
#include "Thyra_StateFuncModelEvaluatorBase.hpp" // Implementation
#include "Thyra_DefaultProductVectorSpace.hpp"
#include "Thyra_PhysicallyBlockedLinearOpWithSolveBase.hpp" // Interface
#include "Thyra_DefaultBlockedTriangularLinearOpWithSolve.hpp" // Implementation
#include "Thyra_ModelEvaluatorDelegatorBase.hpp"
#include "Thyra_ModelEvaluatorHelpers.hpp"
#include "Thyra_DefaultMultiVectorProductVectorSpace.hpp"
#include "Thyra_DefaultMultiVectorProductVector.hpp"
#include "Thyra_DefaultMultiVectorLinearOpWithSolve.hpp"
#include "Teuchos_implicit_cast.hpp"
#include "Teuchos_Assert.hpp"


namespace Rythmos {


template<class Scalar>
class ForwardSensitivityModelEvaluator
  : virtual public Thyra::StateFuncModelEvaluatorBase<Scalar>
{
public:


  ForwardSensitivityModelEvaluator();

  void initializeStructure(
    const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > &stateModel,
    const int p_index
    );
  
  Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >
  getStateModel() const;
  
  int get_p_index() const;

  void initializeState(
    const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint,
    const Teuchos::RCP<const Thyra::LinearOpWithSolveBase<Scalar> > &W_tilde,
    const Scalar &coeff_x_dot,
    const Scalar &coeff_x,
    const Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > &DfDx_dot = Teuchos::null,
    const Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> > &DfDp = Teuchos::null
    );

  // 2007/05/22: rabartl: ToDo: Add an InterpolationBufferBase
  // stateInterpBuffer object to the initailizeState(...) function that can be
  // used to get x and x_dot at different points in time t.  Then, modify the
  // logic to recompute all of the needed matrices if t != t_base (as passed
  // in through stateBasePoint).  The values of x(t) and xdot(t) can then be
  // gotten from the stateInterpBuffer object!
  


  Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > get_x_space() const;
  Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > get_f_space() const;
  Thyra::ModelEvaluatorBase::InArgs<Scalar> getNominalValues() const;
  Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> > create_W() const;
  Thyra::ModelEvaluatorBase::InArgs<Scalar> createInArgs() const;


private:


  Thyra::ModelEvaluatorBase::OutArgs<Scalar> createOutArgsImpl() const;
  void evalModelImpl(
    const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
    const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs
    ) const;


private:

  // /////////////////////////
  // Private data members

  Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > stateModel_;
  int p_index_;
  int np_;
  Teuchos::RCP<const Thyra::DefaultMultiVectorProductVectorSpace<Scalar> > s_bar_space_;
  Teuchos::RCP<const Thyra::DefaultMultiVectorProductVectorSpace<Scalar> > f_sens_space_;
  Thyra::ModelEvaluatorBase::InArgs<Scalar> nominalValues_;

  Thyra::ModelEvaluatorBase::InArgs<Scalar> stateBasePoint_;

  mutable Teuchos::RCP<const Thyra::LinearOpWithSolveBase<Scalar> > W_tilde_;
  mutable Scalar coeff_x_dot_;
  mutable Scalar coeff_x_;
  mutable Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > DfDx_dot_;
  mutable Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> > DfDp_;


  // /////////////////////////
  // Private member functions

  void wrapNominalValuesAndBounds();

  void computeDerivativeMatrices(
    const Thyra::ModelEvaluatorBase::InArgs<Scalar> &point
    ) const;
  
};


// /////////////////////////////////
// Implementations


// Constructors/Intializers/Accessors


template<class Scalar>
ForwardSensitivityModelEvaluator<Scalar>::ForwardSensitivityModelEvaluator()
  : p_index_(0), np_(-1)
{}


template<class Scalar>
void ForwardSensitivityModelEvaluator<Scalar>::initializeStructure(
  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > &stateModel,
  const int p_index
  )
{

  typedef Thyra::ModelEvaluatorBase MEB;

  //
  // Validate input
  //

  TEST_FOR_EXCEPT( is_null(stateModel) );
  TEST_FOR_EXCEPTION(
    !( 0 <= p_index && p_index < stateModel->Np() ), std::logic_error,
    "Error, p_index does not fall in the range [0,"<<(stateModel->Np()-1)<<"]!" );
  // ToDo: Validate support for DfDp!

  //
  // Set the input objects
  //

  stateModel_ = stateModel;
  p_index_ = p_index;
  np_ = stateModel_->get_p_space(p_index)->dim();

  //
  // Create the structure of the model
  //

  s_bar_space_ = Thyra::multiVectorProductVectorSpace(
    stateModel_->get_x_space(), np_
    );

  f_sens_space_ = Thyra::multiVectorProductVectorSpace(
    stateModel_->get_f_space(), np_
    );

  nominalValues_ = this->createInArgs();

}


template<class Scalar>
Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >
ForwardSensitivityModelEvaluator<Scalar>::getStateModel() const
{
  return stateModel_;
}


template<class Scalar>
int ForwardSensitivityModelEvaluator<Scalar>::get_p_index() const
{
  return p_index_;
}


template<class Scalar>
void ForwardSensitivityModelEvaluator<Scalar>::initializeState(
  const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint,
  const Teuchos::RCP<const Thyra::LinearOpWithSolveBase<Scalar> > &W_tilde,
  const Scalar &coeff_x_dot,
  const Scalar &coeff_x,
  const Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > &DfDx_dot,
  const Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> > &DfDp
  )
{

  typedef Thyra::ModelEvaluatorBase MEB;

#ifdef TEUCHOS_DEBUG
  TEST_FOR_EXCEPTION(
    is_null(stateModel_), std::logic_error,
    "Error, you must call intializeStructure(...) before you call initializeState(...)"
    );
  TEST_FOR_EXCEPT( is_null(stateBasePoint.get_x()) );
  TEST_FOR_EXCEPT( is_null(stateBasePoint.get_x_dot()) );
  TEST_FOR_EXCEPT( is_null(stateBasePoint.get_p(p_index_)) );
  // What about the other parameter values?  We really can't say anything
  // about these and we can't check them.  They can be null just fine.
  if (!is_null(W_tilde)) {
    THYRA_ASSERT_VEC_SPACES("",*W_tilde->domain(),*stateModel_->get_x_space());
    THYRA_ASSERT_VEC_SPACES("",*W_tilde->range(),*stateModel_->get_f_space());
  }
  if (!is_null(DfDx_dot)) {
    THYRA_ASSERT_VEC_SPACES("",*DfDx_dot->domain(),*stateModel_->get_x_space());
    THYRA_ASSERT_VEC_SPACES("",*DfDx_dot->range(),*stateModel_->get_f_space());
  }
  if (!is_null(DfDp)) {
    THYRA_ASSERT_VEC_SPACES("",*DfDp->domain(),*stateModel_->get_p_space(p_index_));
    THYRA_ASSERT_VEC_SPACES("",*DfDp->range(),*stateModel_->get_f_space());
  }
#endif

  stateBasePoint_ = stateBasePoint;

  // Set whatever derivatives where passed in.  If an input in null, then the
  // member will be null and the null linear operators will be computed later
  // just in time.

  W_tilde_ = W_tilde;
  coeff_x_dot_ = coeff_x_dot;
  coeff_x_ = coeff_x;
  DfDx_dot_ = DfDx_dot;
  DfDp_ = DfDp;

  wrapNominalValuesAndBounds();

}


// Public functions overridden from ModelEvaulator


template<class Scalar>
Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
ForwardSensitivityModelEvaluator<Scalar>::get_x_space() const
{
  return s_bar_space_;
}


template<class Scalar>
Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
ForwardSensitivityModelEvaluator<Scalar>::get_f_space() const
{
  return f_sens_space_;
}


template<class Scalar>
Thyra::ModelEvaluatorBase::InArgs<Scalar>
ForwardSensitivityModelEvaluator<Scalar>::getNominalValues() const
{
  return nominalValues_;
}


template<class Scalar>
Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> >
ForwardSensitivityModelEvaluator<Scalar>::create_W() const
{
  return Thyra::multiVectorLinearOpWithSolve<Scalar>();
}


template<class Scalar>
Thyra::ModelEvaluatorBase::InArgs<Scalar>
ForwardSensitivityModelEvaluator<Scalar>::createInArgs() const
{
  typedef Thyra::ModelEvaluatorBase MEB;
  MEB::InArgs<Scalar> stateModelInArgs = stateModel_->createInArgs();
  MEB::InArgsSetup<Scalar> inArgs;
  inArgs.setModelEvalDescription(this->description());
  inArgs.setSupports( MEB::IN_ARG_x_dot,
    stateModelInArgs.supports(MEB::IN_ARG_x_dot) );
  inArgs.setSupports( MEB::IN_ARG_x );
  inArgs.setSupports( MEB::IN_ARG_t );
  inArgs.setSupports( MEB::IN_ARG_alpha,
    stateModelInArgs.supports(MEB::IN_ARG_alpha) );
  inArgs.setSupports( MEB::IN_ARG_beta,
    stateModelInArgs.supports(MEB::IN_ARG_beta) );
  return inArgs;
}


// Private functions overridden from ModelEvaulatorDefaultBase


template<class Scalar>
Thyra::ModelEvaluatorBase::OutArgs<Scalar>
ForwardSensitivityModelEvaluator<Scalar>::createOutArgsImpl() const
{

  typedef Thyra::ModelEvaluatorBase MEB;

  MEB::OutArgs<Scalar> stateModelOutArgs = stateModel_->createOutArgs();
  MEB::OutArgsSetup<Scalar> outArgs;

  outArgs.setModelEvalDescription(this->description());

  outArgs.setSupports(MEB::OUT_ARG_f);

  if (stateModelOutArgs.supports(MEB::OUT_ARG_W) ) {
    outArgs.setSupports(MEB::OUT_ARG_W);
    outArgs.set_W_properties(stateModelOutArgs.get_W_properties());
  }

  return outArgs;

}


template<class Scalar>
void ForwardSensitivityModelEvaluator<Scalar>::evalModelImpl(
  const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
  const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs
  ) const
{

  using Teuchos::rcp_dynamic_cast;
  typedef Teuchos::ScalarTraits<Scalar> ST;
  typedef Thyra::ModelEvaluatorBase MEB;
  typedef Teuchos::VerboseObjectTempState<Thyra::ModelEvaluatorBase> VOTSME;

  THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_GEN_BEGIN(
    "ForwardSensitivityModelEvaluator", inArgs, outArgs, Teuchos::null );

  //
  // Update the derivative matrices if they are not already updated for the
  // given time!.
  //
  
  {
#ifdef ENABLE_RYTHMOS_TIMERS
    TEUCHOS_FUNC_TIME_MONITOR("Rythmos:ForwardSensitivityModelEvaluator::evalModel: computeMatrices");
#endif
    computeDerivativeMatrices(inArgs);
  }

  //
  // InArgs
  //

  RCP<const Thyra::DefaultMultiVectorProductVector<Scalar> >
    s_bar = rcp_dynamic_cast<const Thyra::DefaultMultiVectorProductVector<Scalar> >(
      inArgs.get_x().assert_not_null(), true
      );
  RCP<const Thyra::DefaultMultiVectorProductVector<Scalar> >
    s_bar_dot = rcp_dynamic_cast<const Thyra::DefaultMultiVectorProductVector<Scalar> >(
      inArgs.get_x_dot().assert_not_null(), true
      );
  const Scalar
    alpha = inArgs.get_alpha();
  const Scalar
    beta = inArgs.get_beta();

  RCP<const Thyra::MultiVectorBase<Scalar> >
    S = s_bar->getMultiVector();
  RCP<const Thyra::MultiVectorBase<Scalar> >
    S_dot = s_bar_dot->getMultiVector();
  
  //
  // OutArgs
  //

  RCP<Thyra::DefaultMultiVectorProductVector<Scalar> >
    f_sens = rcp_dynamic_cast<Thyra::DefaultMultiVectorProductVector<Scalar> >(
      outArgs.get_f(), true
      );
  RCP<Thyra::DefaultMultiVectorLinearOpWithSolve<Scalar> >
    W_sens = rcp_dynamic_cast<Thyra::DefaultMultiVectorLinearOpWithSolve<Scalar> >(
      outArgs.get_W(), true
      );

  RCP<Thyra::MultiVectorBase<Scalar> >
    F_sens = f_sens->getNonconstMultiVector().assert_not_null();

  //
  // Compute the requested functions
  //

  if(!is_null(F_sens)) {

#ifdef ENABLE_RYTHMOS_TIMERS
    TEUCHOS_FUNC_TIME_MONITOR("Rythmos:ForwardSensitivityModelEvaluator::evalModel: computeSens");
#endif

    // S_diff =  -(coeff_x_dot/coeff_x)*S + S_dot
    Teuchos::RCP<Thyra::MultiVectorBase<Scalar> >
      S_diff = createMembers( stateModel_->get_x_space(), np_ );
    V_StVpV( &*S_diff, Scalar(-coeff_x_dot_/coeff_x_), *S, *S_dot );
    // F_sens = (1/coeff_x) * W_tilde * S
    Thyra::apply(
      *W_tilde_, Thyra::NOTRANS,
      *S, &*F_sens,
      Scalar(1.0/coeff_x_), ST::zero()
      );
    // F_sens += d(f)/d(x_dot) * S_diff
    Thyra::apply(
      *DfDx_dot_, Thyra::NOTRANS,
      *S_diff, &*F_sens,
      ST::one(), ST::one()
      );
    // F_sens += d(f)/d(p)
    Vp_V( &*F_sens, *DfDp_ );
  }
  
  if(!is_null(W_sens)) {
    TEST_FOR_EXCEPTION(
      alpha != coeff_x_dot_, std::logic_error,
      "Error, alpha="<<alpha<<" != coeff_x_dot="<<coeff_x_dot_
      <<" with difference = "<<(alpha-coeff_x_dot_)<<"!" );
    TEST_FOR_EXCEPTION(
      beta != coeff_x_, std::logic_error,
      "Error, beta="<<beta<<" != coeff_x="<<coeff_x_
      <<" with difference = "<<(beta-coeff_x_)<<"!" );
    W_sens->initialize( W_tilde_, s_bar_space_, f_sens_space_ );
    
  }

  THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_END();

}


// private


template<class Scalar>
void ForwardSensitivityModelEvaluator<Scalar>::wrapNominalValuesAndBounds()
{

  using Teuchos::rcp_dynamic_cast;
  typedef Thyra::ModelEvaluatorBase MEB;

  // nominalValues_.clear(); // ToDo: Implement this!

  nominalValues_.set_t(stateModel_->getNominalValues().get_t());
  
  // 2007/05/22: rabartl: Note: Currently there is not much of a reason to set
  // an initial condition here since the initial condition for the
  // sensitivities is really being set in the ForwardSensitivityStepper
  // object!  In the future, a more general use of this class might benefit
  // from setting the initial condition here.

}


template<class Scalar>
void ForwardSensitivityModelEvaluator<Scalar>::computeDerivativeMatrices(
  const Thyra::ModelEvaluatorBase::InArgs<Scalar> &point
  ) const
{

  typedef Thyra::ModelEvaluatorBase MEB;
  typedef Teuchos::VerboseObjectTempState<MEB> VOTSME;

  Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
  Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();

  const Scalar
    t_base = stateBasePoint_.get_t(),
    t = point.get_t();

  TEUCHOS_ASSERT_EQUALITY( t , t_base );

  if (is_null(W_tilde_)) {
    TEST_FOR_EXCEPT("ToDo: compute W_tilde from scratch!");
  }
  
  if ( is_null(DfDx_dot_) || is_null(DfDp_) ) {

    MEB::InArgs<Scalar> inArgs = stateBasePoint_;
    MEB::OutArgs<Scalar> outArgs = stateModel_->createOutArgs();
    
    Teuchos::RCP<Thyra::LinearOpBase<Scalar> > DfDx_dot_compute;
    if (is_null(DfDx_dot_)) {
      DfDx_dot_compute = stateModel_->create_W_op();
      inArgs.set_alpha(1.0);
      inArgs.set_beta(0.0);
      outArgs.set_W_op(DfDx_dot_compute);
    }

    Teuchos::RCP<Thyra::MultiVectorBase<Scalar> > DfDp_compute;
    if (is_null(DfDp_)) {
      DfDp_compute = Thyra::create_DfDp_mv(
        *stateModel_,p_index_,
        MEB::DERIV_MV_BY_COL
        ).getMultiVector();
      outArgs.set_DfDp(
        p_index_,
        MEB::Derivative<Scalar>(DfDp_compute,MEB::DERIV_MV_BY_COL)
        );
    }
    
    VOTSME stateModel_outputTempState(stateModel_,out,verbLevel);
    stateModel_->evalModel(inArgs,outArgs);
    if (!is_null(DfDx_dot_compute))
      DfDx_dot_ = DfDx_dot_compute;
    if (!is_null(DfDp_compute))
      DfDp_ = DfDp_compute;
  
  }

}


} // namespace Rythmos


#endif // RYTHMOS_FORWARD_SENSITIVITY_MODEL_EVALUATOR_HPP

---