Tempus_IntegratorAdjointSensitivity_impl.hpp

Go to the documentation of this file.

// @HEADER
// ****************************************************************************
//                Tempus: Copyright (2017) Sandia Corporation
//
// Distributed under BSD 3-clause license (See accompanying file Copyright.txt)
// ****************************************************************************
// @HEADER
 
#ifndef Tempus_IntegratorAdjointSensitivity_impl_hpp
#define Tempus_IntegratorAdjointSensitivity_impl_hpp
 
#include "Teuchos_ParameterList.hpp"
#include "Thyra_DefaultMultiVectorProductVectorSpace.hpp"
#include "Thyra_DefaultMultiVectorProductVector.hpp"
#include "Thyra_DefaultProductVectorSpace.hpp"
#include "Thyra_DefaultProductVector.hpp"
#include "Thyra_VectorStdOps.hpp"
#include "Thyra_MultiVectorStdOps.hpp"
#include "Thyra_ImplicitAdjointModelEvaluator.hpp"
#include "Thyra_LinearOpWithSolveFactoryHelpers.hpp"
 
namespace Tempus {
 
template <class Scalar>
IntegratorAdjointSensitivity<Scalar>::IntegratorAdjointSensitivity(
    const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> > &model,
    const Teuchos::RCP<IntegratorBasic<Scalar> > &state_integrator,
    const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> > &adjoint_model,
    const Teuchos::RCP<AdjointAuxSensitivityModelEvaluator<Scalar> > &adjoint_aux_model,
    const Teuchos::RCP<IntegratorBasic<Scalar> > &adjoint_integrator,
    const Teuchos::RCP<SolutionHistory<Scalar> > &solutionHistory,
    const int p_index,
    const int g_index,
    const bool g_depends_on_p,
    const bool f_depends_on_p,
    const bool ic_depends_on_p,
    const bool mass_matrix_is_identity)
    : model_(model)
    , state_integrator_(state_integrator)
    , adjoint_model_(adjoint_model)
    , adjoint_aux_model_(adjoint_aux_model)
    , adjoint_integrator_(adjoint_integrator)
    , solutionHistory_(solutionHistory)
    , p_index_(p_index)
    , g_index_(g_index)
    , g_depends_on_p_(g_depends_on_p)
    , f_depends_on_p_(f_depends_on_p)
    , ic_depends_on_p_(ic_depends_on_p)
    , mass_matrix_is_identity_(mass_matrix_is_identity)
    , stepMode_(SensitivityStepMode::Forward)
{
 
  TEUCHOS_TEST_FOR_EXCEPTION(
      getStepper()->getUseFSAL(), std::logic_error,
      "Error - IntegratorAdjointSensitivity(): Cannot use FSAL with\n"
      "        IntegratorAdjointSensitivity, because the state and adjoint\n"
      "        integrators require ModelEvaluator evaluation in the\n"
      "        constructor to make the initial conditions consistent.\n"
      "        For the adjoint integrator, this requires special construction\n"
      "        which has not been implemented yet.\n");
}
 
template<class Scalar>
IntegratorAdjointSensitivity<Scalar>::
IntegratorAdjointSensitivity()
{
  state_integrator_   = createIntegratorBasic<Scalar>();
  adjoint_integrator_ = createIntegratorBasic<Scalar>();
  stepMode_ = SensitivityStepMode::Forward;
}
 
template<class Scalar>
bool
IntegratorAdjointSensitivity<Scalar>::
advanceTime()
{
  const Scalar tfinal =
    state_integrator_->getTimeStepControl()->getFinalTime();
  return advanceTime(tfinal);
}
 
template<class Scalar>
bool
IntegratorAdjointSensitivity<Scalar>::
advanceTime(const Scalar timeFinal)
{
  TEMPUS_FUNC_TIME_MONITOR_DIFF("Tempus::IntegratorAdjointSensitivity::advanceTime()", TEMPUS_AS_AT);
 
  using Teuchos::RCP;
  using Teuchos::rcp_dynamic_cast;
  using Thyra::VectorBase;
  using Thyra::VectorSpaceBase;
  using Thyra::MultiVectorBase;
  using Thyra::LinearOpBase;
  using Thyra::LinearOpWithSolveBase;
  using Thyra::LinearOpWithSolveFactoryBase;
  using Thyra::PreconditionerBase;
  using Thyra::PreconditionerFactoryBase;
  using Thyra::createMember;
  using Thyra::createMembers;
  using Thyra::assign;
  typedef Thyra::ModelEvaluatorBase MEB;
  typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
  typedef Thyra::DefaultProductVector<Scalar> DPV;
 
  // Get initial state for later
  RCP<const SolutionHistory<Scalar> > state_solution_history =
    state_integrator_->getSolutionHistory();
  RCP<const SolutionState<Scalar> > initial_state =
    (*state_solution_history)[0];
 
  // Run state integrator and get solution
  bool state_status = true;
  {
    TEMPUS_FUNC_TIME_MONITOR_DIFF("Tempus::IntegratorAdjointSensitivity::advanceTime::state", TEMPUS_AS_AT_FWD);
    stepMode_ = SensitivityStepMode::Forward;
    state_status = state_integrator_->advanceTime(timeFinal);
  }
 
  // For at least some time-stepping methods, the time of the last time step
  // may not be timeFinal (e.g., it may be greater by at most delta_t).
  // But since the adjoint model requires timeFinal in its formulation, reset
  // it to the achieved final time.
  adjoint_aux_model_->setFinalTime(state_integrator_->getTime());
 
  // Set solution history in adjoint stepper
  adjoint_aux_model_->setForwardSolutionHistory(state_solution_history);
 
  // Compute dg/dx
  RCP<const VectorSpaceBase<Scalar> > g_space = model_->get_g_space(g_index_);
  RCP<const VectorSpaceBase<Scalar> > x_space = model_->get_x_space();
  const int num_g = g_space->dim();
  RCP<MultiVectorBase<Scalar> > dgdx = createMembers(x_space, num_g);
  MEB::InArgs<Scalar> inargs = model_->getNominalValues();
  RCP<const SolutionState<Scalar> > state =
    state_solution_history->getCurrentState();
  inargs.set_t(state->getTime());
  inargs.set_x(state->getX());
  inargs.set_x_dot(state->getXDot());
  MEB::OutArgs<Scalar> outargs = model_->createOutArgs();
  MEB::OutArgs<Scalar> adj_outargs = adjoint_model_->createOutArgs();
  outargs.set_DgDx(g_index_,
                   MEB::Derivative<Scalar>(dgdx, MEB::DERIV_MV_GRADIENT_FORM));
  model_->evalModel(inargs, outargs);
  outargs.set_DgDx(g_index_, MEB::Derivative<Scalar>());
 
  // Compute ICs == [ (df/dx_dot)^{-T} (dg/dx)^T; 0 ]
  // For explicit form, we are relying on the user to inform us the
  // the mass matrix is the identity.  It would be nice to be able to determine
  // somehow automatically that we are using an explicit stepper.
  RCP<DPV> adjoint_init =
    rcp_dynamic_cast<DPV>(Thyra::createMember(adjoint_aux_model_->get_x_space()));
  RCP<MultiVectorBase<Scalar> > adjoint_init_mv =
    rcp_dynamic_cast<DMVPV>(adjoint_init->getNonconstVectorBlock(0))->getNonconstMultiVector();
  assign(adjoint_init->getNonconstVectorBlock(1).ptr(),
         Teuchos::ScalarTraits<Scalar>::zero());
  if (mass_matrix_is_identity_)
    assign(adjoint_init_mv.ptr(), *dgdx);
  else {
    inargs.set_alpha(1.0);
    inargs.set_beta(0.0);
    RCP<LinearOpWithSolveBase<Scalar> > W;
    if (adj_outargs.supports(MEB::OUT_ARG_W)) {
      // Model supports W
      W = adjoint_model_->create_W();
      adj_outargs.set_W(W);
      adjoint_model_->evalModel(inargs, adj_outargs);
      adj_outargs.set_W(Teuchos::null);
    }
    else {
      // Otherwise model must support a W_op and W factory
      RCP<const LinearOpWithSolveFactoryBase<Scalar> > lowsfb =
        adjoint_model_->get_W_factory();
      TEUCHOS_TEST_FOR_EXCEPTION(
        lowsfb == Teuchos::null, std::logic_error,
        "Adjoint ME must support W out-arg or provide a W_factory for non-identity mass matrix");
 
      // Compute W_op (and W_prec if supported)
      RCP<LinearOpBase<Scalar> > W_op = adjoint_model_->create_W_op();
      adj_outargs.set_W_op(W_op);
      RCP<PreconditionerFactoryBase<Scalar> > prec_factory =
        lowsfb->getPreconditionerFactory();
      RCP<PreconditionerBase<Scalar> > W_prec;
      if (prec_factory != Teuchos::null)
        W_prec = prec_factory->createPrec();
      else if (adj_outargs.supports(MEB::OUT_ARG_W_prec)) {
        W_prec = adjoint_model_->create_W_prec();
        adj_outargs.set_W_prec(W_prec);
      }
      adjoint_model_->evalModel(inargs, adj_outargs);
      adj_outargs.set_W_op(Teuchos::null);
      if (adj_outargs.supports(MEB::OUT_ARG_W_prec))
        adj_outargs.set_W_prec(Teuchos::null);
 
      // Create and initialize W
      W = lowsfb->createOp();
      if (W_prec != Teuchos::null) {
        if (prec_factory != Teuchos::null)
          prec_factory->initializePrec(
            Thyra::defaultLinearOpSource<Scalar>(W_op), W_prec.get());
        Thyra::initializePreconditionedOp<Scalar>(
          *lowsfb, W_op, W_prec, W.ptr());
      }
      else
        Thyra::initializeOp<Scalar>(*lowsfb, W_op, W.ptr());
    }
    TEUCHOS_TEST_FOR_EXCEPTION(
      W == Teuchos::null, std::logic_error,
      "A null W has been encountered in Tempus::IntegratorAdjointSensitivity::advanceTime!\n");
    // Initialize adjoint_init_mv to zero before solve for linear solvers that
    // use what is passed in as the initial guess
    assign(adjoint_init_mv.ptr(), Teuchos::ScalarTraits<Scalar>::zero());
    W->solve(Thyra::NOTRANS, *dgdx, adjoint_init_mv.ptr());
  }
 
  // Run sensitivity integrator and get solution
  bool sens_status = true;
  {
    TEMPUS_FUNC_TIME_MONITOR_DIFF("Tempus::IntegratorAdjointSensitivity::advanceTime::adjoint", TEMPUS_AS_AT_ADJ);
    stepMode_ = SensitivityStepMode::Adjoint;
    const Scalar tinit = adjoint_integrator_->getTimeStepControl()->getInitTime();
    adjoint_integrator_->initializeSolutionHistory(tinit, adjoint_init);
    sens_status = adjoint_integrator_->advanceTime(timeFinal);
  }
  RCP<const SolutionHistory<Scalar> > adjoint_solution_history =
    adjoint_integrator_->getSolutionHistory();
 
  // Compute dg/dp at final time T
  RCP<const VectorSpaceBase<Scalar> > p_space = model_->get_p_space(p_index_);
  dgdp_ = createMembers(p_space, num_g);
  if (g_depends_on_p_) {
    MEB::DerivativeSupport dgdp_support =
      outargs.supports(MEB::OUT_ARG_DgDp, g_index_, p_index_);
    if (dgdp_support.supports(MEB::DERIV_MV_GRADIENT_FORM)) {
      outargs.set_DgDp(g_index_, p_index_,
                       MEB::Derivative<Scalar>(dgdp_,
                                               MEB::DERIV_MV_GRADIENT_FORM));
      model_->evalModel(inargs, outargs);
    }
    else if (dgdp_support.supports(MEB::DERIV_MV_JACOBIAN_FORM)) {
      const int num_p = p_space->dim();
      RCP<MultiVectorBase<Scalar> > dgdp_trans =
        createMembers(g_space, num_p);
      outargs.set_DgDp(g_index_, p_index_,
                       MEB::Derivative<Scalar>(dgdp_trans,
                                               MEB::DERIV_MV_JACOBIAN_FORM));
      model_->evalModel(inargs, outargs);
      Thyra::DetachedMultiVectorView<Scalar> dgdp_view(*dgdp_);
      Thyra::DetachedMultiVectorView<Scalar> dgdp_trans_view(*dgdp_trans);
      for (int i=0; i<num_p; ++i)
        for (int j=0; j<num_g; ++j)
          dgdp_view(i,j) = dgdp_trans_view(j,i);
    }
    else
      TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
                                 "Invalid dg/dp support");
    outargs.set_DgDp(g_index_, p_index_, MEB::Derivative<Scalar>());
  }
  else
    assign(dgdp_.ptr(), Scalar(0.0));
 
  // Add in initial condition term = (dx/dp^T(0))*(df/dx_dot^T(0))*y(0)
  // If dxdp_init_ is null, assume it is zero
  if (ic_depends_on_p_ && dxdp_init_ != Teuchos::null) {
    RCP<const SolutionState<Scalar> > adjoint_state =
      adjoint_solution_history->getCurrentState();
    RCP<const VectorBase<Scalar> > adjoint_x =
      rcp_dynamic_cast<const DPV>(adjoint_state->getX())->getVectorBlock(0);
    RCP<const MultiVectorBase<Scalar> > adjoint_mv =
      rcp_dynamic_cast<const DMVPV>(adjoint_x)->getMultiVector();
    if (mass_matrix_is_identity_)
      dxdp_init_->apply(Thyra::CONJTRANS, *adjoint_mv, dgdp_.ptr(), Scalar(1.0),
                        Scalar(1.0));
    else {
      inargs.set_t(initial_state->getTime());
      inargs.set_x(initial_state->getX());
      inargs.set_x_dot(initial_state->getXDot());
      inargs.set_alpha(1.0);
      inargs.set_beta(0.0);
      RCP<LinearOpBase<Scalar> > W_op = adjoint_model_->create_W_op();
      adj_outargs.set_W_op(W_op);
      adjoint_model_->evalModel(inargs, adj_outargs);
      adj_outargs.set_W_op(Teuchos::null);
      RCP<MultiVectorBase<Scalar> > tmp = createMembers(x_space, num_g);
      W_op->apply(Thyra::NOTRANS, *adjoint_mv, tmp.ptr(), Scalar(1.0),
                  Scalar(0.0));
      dxdp_init_->apply(Thyra::CONJTRANS, *tmp, dgdp_.ptr(), Scalar(1.0),
                        Scalar(1.0));
    }
  }
 
  // Add in model parameter term = \int_0^T( (df/dp^T(t)*y(t) )dt which
  // is computed during the adjoint integration as an auxiliary integral
  // (2nd block of the solution vector)
  if (f_depends_on_p_) {
    RCP<const SolutionState<Scalar> > adjoint_state =
      adjoint_solution_history->getCurrentState();
    RCP<const VectorBase<Scalar> > z =
      rcp_dynamic_cast<const DPV>(adjoint_state->getX())->getVectorBlock(1);
    RCP<const MultiVectorBase<Scalar> > z_mv =
      rcp_dynamic_cast<const DMVPV>(z)->getMultiVector();
    Thyra::V_VmV(dgdp_.ptr(), *dgdp_, *z_mv);
  }
 
  buildSolutionHistory(state_solution_history, adjoint_solution_history);
 
  return state_status && sens_status;
}
 
template<class Scalar>
Scalar
IntegratorAdjointSensitivity<Scalar>::
getTime() const
{
  return solutionHistory_->getCurrentTime();
}
 
template<class Scalar>
int
IntegratorAdjointSensitivity<Scalar>::
getIndex() const
{
  return solutionHistory_->getCurrentIndex();
}
 
template<class Scalar>
Status
IntegratorAdjointSensitivity<Scalar>::
getStatus() const
{
  Status state_status = state_integrator_->getStatus();
  Status sens_status = adjoint_integrator_->getStatus();
  if (state_status == FAILED || sens_status == FAILED)
    return FAILED;
  if (state_status == WORKING || sens_status == WORKING)
    return WORKING;
  return PASSED;
}
 
template <class Scalar>
void IntegratorAdjointSensitivity<Scalar>::setStatus(const Status st) {
  state_integrator_->setStatus(st);
  adjoint_integrator_->setStatus(st);
}
 
template<class Scalar>
Teuchos::RCP<Stepper<Scalar> >
IntegratorAdjointSensitivity<Scalar>::
getStepper() const
{
  return state_integrator_->getStepper();
}
 
template<class Scalar>
Teuchos::RCP<const SolutionHistory<Scalar> >
IntegratorAdjointSensitivity<Scalar>::
getSolutionHistory() const
{
  return solutionHistory_;
}
 
template<class Scalar>
Teuchos::RCP<const SolutionHistory<Scalar> >
IntegratorAdjointSensitivity<Scalar>::
getStateSolutionHistory() const
{
  return state_integrator_->getSolutionHistory();
}
 
template<class Scalar>
Teuchos::RCP<const SolutionHistory<Scalar> >
IntegratorAdjointSensitivity<Scalar>::
getSensSolutionHistory() const
{
  return adjoint_integrator_->getSolutionHistory();
}
 
template<class Scalar>
Teuchos::RCP<SolutionHistory<Scalar> >
IntegratorAdjointSensitivity<Scalar>::
getNonConstSolutionHistory()
{
  return solutionHistory_;
}
 
template<class Scalar>
Teuchos::RCP<const TimeStepControl<Scalar> >
IntegratorAdjointSensitivity<Scalar>::
getTimeStepControl() const
{
  return state_integrator_->getTimeStepControl();
}
 
template<class Scalar>
Teuchos::RCP<TimeStepControl<Scalar> >
IntegratorAdjointSensitivity<Scalar>::
getNonConstTimeStepControl()
{
  return state_integrator_->getNonConstTimeStepControl();
}
 
template<class Scalar>
Teuchos::RCP<TimeStepControl<Scalar> >
IntegratorAdjointSensitivity<Scalar>::
getStateNonConstTimeStepControl()
{
  return state_integrator_->getNonConstTimeStepControl();
}
 
template<class Scalar>
Teuchos::RCP<TimeStepControl<Scalar> >
IntegratorAdjointSensitivity<Scalar>::
getSensNonConstTimeStepControl()
{
  return adjoint_integrator_->getNonConstTimeStepControl();
}
 
template<class Scalar>
void IntegratorAdjointSensitivity<Scalar>::
initializeSolutionHistory(Scalar t0,
  Teuchos::RCP<const Thyra::VectorBase<Scalar> > x0,
  Teuchos::RCP<const Thyra::VectorBase<Scalar> > xdot0,
  Teuchos::RCP<const Thyra::VectorBase<Scalar> > xdotdot0,
  Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> > DxDp0,
  Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> > /* DxdotDp0 */,
  Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> > /* DxdotdotDp0 */)
{
  state_integrator_->initializeSolutionHistory(t0, x0, xdot0, xdotdot0);
  dxdp_init_ = DxDp0;
}
 
template<class Scalar>
Teuchos::RCP<IntegratorObserver<Scalar> >
IntegratorAdjointSensitivity<Scalar>::
getObserver()
{
  return state_integrator_->getObserver();
}
 
template<class Scalar>
void IntegratorAdjointSensitivity<Scalar>::
setObserver(Teuchos::RCP<IntegratorObserver<Scalar> > obs)
{
  state_integrator_->setObserver(obs);
  // ETP 1/12/22 Disabling passing of the observer to the adjoint
  // integrator to work around issues in Piro
  //adjoint_integrator_->setObserver(obs);
}
 
template<class Scalar>
void IntegratorAdjointSensitivity<Scalar>::
initialize()
{
  state_integrator_->initialize();
  adjoint_integrator_->initialize();
}
 
template<class Scalar>
Teuchos::RCP<const Thyra::VectorBase<Scalar> >
IntegratorAdjointSensitivity<Scalar>::
getX() const
{
  return state_integrator_->getX();
}
 
template<class Scalar>
Teuchos::RCP<const Thyra::VectorBase<Scalar> >
IntegratorAdjointSensitivity<Scalar>::
getXDot() const
{
  return state_integrator_->getXDot();
}
 
template<class Scalar>
Teuchos::RCP<const Thyra::VectorBase<Scalar> >
IntegratorAdjointSensitivity<Scalar>::
getXDotDot() const
{
  return state_integrator_->getXDotDot();
}
 
template<class Scalar>
Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> >
IntegratorAdjointSensitivity<Scalar>::
getY() const
{
  using Teuchos::RCP;
  using Teuchos::rcp_dynamic_cast;
  typedef Thyra::DefaultProductVector<Scalar> DPV;
  typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
  RCP<const DPV> pv =
    rcp_dynamic_cast<const DPV>(adjoint_integrator_->getX());
  RCP<const DMVPV> mvpv =
    rcp_dynamic_cast<const DMVPV>(pv->getVectorBlock(0));
  return mvpv->getMultiVector();
}
 
template<class Scalar>
Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> >
IntegratorAdjointSensitivity<Scalar>::
getYDot() const
{
  using Teuchos::RCP;
  using Teuchos::rcp_dynamic_cast;
  typedef Thyra::DefaultProductVector<Scalar> DPV;
  typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
  RCP<const DPV> pv =
    rcp_dynamic_cast<const DPV>(adjoint_integrator_->getXDot());
  RCP<const DMVPV> mvpv =
    rcp_dynamic_cast<const DMVPV>(pv->getVectorBlock(0));
  return mvpv->getMultiVector();
}
 
template<class Scalar>
Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> >
IntegratorAdjointSensitivity<Scalar>::
getYDotDot() const
{
  using Teuchos::RCP;
  using Teuchos::rcp_dynamic_cast;
  typedef Thyra::DefaultProductVector<Scalar> DPV;
  typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;
  RCP<const DPV> pv =
    rcp_dynamic_cast<const DPV>(adjoint_integrator_->getXDotDot());
  RCP<const DMVPV> mvpv =
    rcp_dynamic_cast<const DMVPV>(pv->getVectorBlock(0));
  return mvpv->getMultiVector();
}
 
template<class Scalar>
Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> >
IntegratorAdjointSensitivity<Scalar>::
getDgDp() const
{
  return dgdp_;
}
 
template<class Scalar>
std::string
IntegratorAdjointSensitivity<Scalar>::
description() const
{
  std::string name = "Tempus::IntegratorAdjointSensitivity";
  return(name);
}
 
template<class Scalar>
void
IntegratorAdjointSensitivity<Scalar>::
describe(
  Teuchos::FancyOStream          &out,
  const Teuchos::EVerbosityLevel verbLevel) const
{
  auto l_out = Teuchos::fancyOStream( out.getOStream() );
  Teuchos::OSTab ostab(*l_out, 2, this->description());
  l_out->setOutputToRootOnly(0);
 
  *l_out << description() << "::describe" << std::endl;
  state_integrator_->describe(*l_out, verbLevel);
  adjoint_integrator_->describe(*l_out, verbLevel);
}
 
template<class Scalar>
SensitivityStepMode
IntegratorAdjointSensitivity<Scalar>::
getStepMode() const
{
  return stepMode_;
}
 
template <class Scalar>
Teuchos::RCP<AdjointAuxSensitivityModelEvaluator<Scalar> >
IntegratorAdjointSensitivity<Scalar>::
createAdjointModel(
  const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& model,
  const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& adjoint_model,
  const Teuchos::RCP<Teuchos::ParameterList>& inputPL)
{
  using Teuchos::ParameterList;
  using Teuchos::RCP;
  using Teuchos::rcp;
 
  RCP<ParameterList> spl = Teuchos::parameterList();
  if (inputPL != Teuchos::null)
  {
    *spl = inputPL->sublist("Sensitivities");
  }
  if (spl->isParameter("Response Depends on Parameters"))
    spl->remove("Response Depends on Parameters");
  if (spl->isParameter("Residual Depends on Parameters"))
    spl->remove("Residual Depends on Parameters");
  if (spl->isParameter("IC Depends on Parameters"))
    spl->remove("IC Depends on Parameters");
 
  const Scalar tinit = state_integrator_->getTimeStepControl()->getInitTime();
  const Scalar tfinal = state_integrator_->getTimeStepControl()->getFinalTime();
  return rcp(new AdjointAuxSensitivityModelEvaluator<Scalar>(
               model, adjoint_model, tinit, tfinal, spl));
}
 
template<class Scalar>
void
IntegratorAdjointSensitivity<Scalar>::
buildSolutionHistory(
  const Teuchos::RCP<const SolutionHistory<Scalar> >& state_solution_history,
  const Teuchos::RCP<const SolutionHistory<Scalar> >& adjoint_solution_history)
{
  using Teuchos::RCP;
  using Teuchos::rcp;
  using Teuchos::rcp_dynamic_cast;
  using Teuchos::ParameterList;
  using Thyra::VectorBase;
  using Thyra::MultiVectorBase;
  using Thyra::VectorSpaceBase;
  using Thyra::createMembers;
  using Thyra::multiVectorProductVector;
  using Thyra::assign;
  typedef Thyra::DefaultProductVectorSpace<Scalar> DPVS;
  typedef Thyra::DefaultProductVector<Scalar> DPV;
 
  RCP<const VectorSpaceBase<Scalar> > x_space = model_->get_x_space();
  RCP<const VectorSpaceBase<Scalar> > adjoint_space =
    rcp_dynamic_cast<const DPVS>(adjoint_aux_model_->get_x_space())->getBlock(0);
  Teuchos::Array< RCP<const VectorSpaceBase<Scalar> > > spaces(2);
  spaces[0] = x_space;
  spaces[1] = adjoint_space;
  RCP<const DPVS > prod_space = Thyra::productVectorSpace(spaces());
 
  int num_states = state_solution_history->getNumStates();
  const Scalar t_init = state_integrator_->getTimeStepControl()->getInitTime();
  const Scalar t_final = state_integrator_->getTime();
  for (int i=0; i<num_states; ++i) {
    RCP<const SolutionState<Scalar> > forward_state =
      (*state_solution_history)[i];
    RCP<const SolutionState<Scalar> > adjoint_state =
      adjoint_solution_history->findState(t_final+t_init-forward_state->getTime());
 
    // X
    RCP<DPV> x = Thyra::defaultProductVector(prod_space);
    RCP<const VectorBase<Scalar> > adjoint_x =
      rcp_dynamic_cast<const DPV>(adjoint_state->getX())->getVectorBlock(0);
    assign(x->getNonconstVectorBlock(0).ptr(), *(forward_state->getX()));
    assign(x->getNonconstVectorBlock(1).ptr(), *(adjoint_x));
    RCP<VectorBase<Scalar> > x_b = x;
 
    // X-Dot
    RCP<DPV> x_dot = Thyra::defaultProductVector(prod_space);
    RCP<const VectorBase<Scalar> > adjoint_x_dot =
      rcp_dynamic_cast<const DPV>(adjoint_state->getXDot())->getVectorBlock(0);
    assign(x_dot->getNonconstVectorBlock(0).ptr(), *(forward_state->getXDot()));
    assign(x_dot->getNonconstVectorBlock(1).ptr(), *(adjoint_x_dot));
    RCP<VectorBase<Scalar> > x_dot_b = x_dot;
 
    // X-Dot-Dot
    RCP<DPV> x_dot_dot;
    if (forward_state->getXDotDot() != Teuchos::null) {
       x_dot_dot = Thyra::defaultProductVector(prod_space);
       RCP<const VectorBase<Scalar> > adjoint_x_dot_dot =
         rcp_dynamic_cast<const DPV>(
           adjoint_state->getXDotDot())->getVectorBlock(0);
       assign(x_dot_dot->getNonconstVectorBlock(0).ptr(),
              *(forward_state->getXDotDot()));
       assign(x_dot_dot->getNonconstVectorBlock(1).ptr(),
              *(adjoint_x_dot_dot));
    }
    RCP<VectorBase<Scalar> > x_dot_dot_b = x_dot_dot;
 
    RCP<SolutionState<Scalar> > prod_state = forward_state->clone();
    prod_state->setX(x_b);
    prod_state->setXDot(x_dot_b);
    prod_state->setXDotDot(x_dot_dot_b);
    prod_state->setPhysicsState(Teuchos::null);
    solutionHistory_->addState(prod_state);
  }
}

template<class Scalar>
Teuchos::RCP<IntegratorAdjointSensitivity<Scalar> >
createIntegratorAdjointSensitivity(
    Teuchos::RCP<Teuchos::ParameterList>                     inputPL,
    const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >&      model
    ,const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >& adjoint_model)
{
  // set the parameters
  Teuchos::RCP<Teuchos::ParameterList> spl = Teuchos::parameterList();
  if (inputPL != Teuchos::null)
    *spl = inputPL->sublist("Sensitivities");
 
  int p_index                  = spl->get<int>("Sensitivity Parameter Index", 0);
  int g_index                  = spl->get<int>("Response Function Index", 0);
  bool g_depends_on_p          = spl->get<bool>("Response Depends on Parameters", true);
  bool f_depends_on_p          = spl->get<bool>("Residual Depends on Parameters", true);
  bool ic_depends_on_p         = spl->get<bool>("IC Depends on Parameters", true);
  bool mass_matrix_is_identity = spl->get<bool>("Mass Matrix Is Identity", false);
 
  auto state_integrator = createIntegratorBasic<Scalar>(inputPL, model);
 
  // createAdjointModel
  if (spl->isParameter("Response Depends on Parameters"))
    spl->remove("Response Depends on Parameters");
  if (spl->isParameter("Residual Depends on Parameters"))
    spl->remove("Residual Depends on Parameters");
  if (spl->isParameter("IC Depends on Parameters"))
    spl->remove("IC Depends on Parameters");
 
  const Scalar tinit = state_integrator->getTimeStepControl()->getInitTime();
  const Scalar tfinal = state_integrator->getTimeStepControl()->getFinalTime();
  //auto adjoint_model  = Teuchos::rcp(new AdjointAuxSensitivityModelEvaluator<Scalar>(model, tfinal, spl));
  //TODO: where is the adjoint ME coming from?
 
  Teuchos::RCP<Thyra::ModelEvaluator<Scalar>> adjt_model = adjoint_model;
  if (adjoint_model == Teuchos::null)
    adjt_model = Thyra::implicitAdjointModelEvaluator(model);
 
  auto adjoint_aux_model  = Teuchos::rcp(new AdjointAuxSensitivityModelEvaluator<Scalar>(model, adjt_model, tinit, tfinal, spl));
 
  // Create combined solution histories combining the forward and adjoint
  // solutions.  We do not include the auxiliary part from the adjoint solution.
  auto integrator_name           = inputPL->get<std::string>("Integrator Name");
  auto integratorPL              = Teuchos::sublist(inputPL, integrator_name, true);
  auto shPL                      = Teuchos::sublist(integratorPL, "Solution History", true);
  auto combined_solution_History = createSolutionHistoryPL<Scalar>(shPL);
 
  auto adjoint_integrator = createIntegratorBasic<Scalar>(inputPL, adjoint_aux_model);
 
  Teuchos::RCP<IntegratorAdjointSensitivity<Scalar>> integrator = Teuchos::rcp(new IntegratorAdjointSensitivity<Scalar>(
      model, state_integrator, adjt_model, adjoint_aux_model, adjoint_integrator, combined_solution_History, p_index, g_index, g_depends_on_p,
      f_depends_on_p, ic_depends_on_p, mass_matrix_is_identity));
 
  return (integrator);
}

template<class Scalar>
Teuchos::RCP<IntegratorAdjointSensitivity<Scalar> >
createIntegratorAdjointSensitivity()
{
  Teuchos::RCP<IntegratorAdjointSensitivity<Scalar> > integrator =
    Teuchos::rcp(new IntegratorAdjointSensitivity<Scalar>());
  return(integrator);
}
 
} // namespace Tempus
#endif // Tempus_IntegratorAdjointSensitivity_impl_hpp