Tempus_IntegratorForwardSensitivity_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_IntegratorForwardSensitivity_impl_hpp
#define Tempus_IntegratorForwardSensitivity_impl_hpp

#include "Tempus_SolutionHistory_decl.hpp"
#include "Thyra_DefaultMultiVectorProductVector.hpp"
#include "Thyra_VectorStdOps.hpp"
#include "Thyra_MultiVectorStdOps.hpp"

#include "Tempus_CombinedForwardSensitivityModelEvaluator.hpp"
#include "Tempus_WrapCombinedFSAModelEvaluator.hpp"

namespace Tempus {

template <class Scalar>
IntegratorForwardSensitivity<Scalar>::IntegratorForwardSensitivity(
    const Teuchos::RCP<Thyra::ModelEvaluator<Scalar>> &model,
    const Teuchos::RCP<IntegratorBasic<Scalar>> &integrator,
    const Teuchos::RCP<SensitivityModelEvaluatorBase<Scalar>> &sens_model,
    const Teuchos::RCP<StepperStaggeredForwardSensitivity<Scalar>> &sens_stepper,
    bool use_combined_method)
    : model_(model)
    , integrator_(integrator)
    , sens_model_(sens_model)
    , sens_stepper_(sens_stepper)
    , use_combined_method_(use_combined_method)
{
  integrator_->initialize();
}

template<class Scalar>
IntegratorForwardSensitivity<Scalar>::
IntegratorForwardSensitivity()
{
  integrator_ = createIntegratorBasic<Scalar>();
  integrator_->initialize();
}

template<class Scalar>
void IntegratorForwardSensitivity<Scalar>::
setStepper(
  Teuchos::RCP<Thyra::ModelEvaluator<Scalar> > model)
{
  if (use_combined_method_)
    integrator_->setModel(sens_model_);
  else
    integrator_->setStepper(sens_stepper_);
}

template<class Scalar>
void IntegratorForwardSensitivity<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)
{
  using Teuchos::RCP;
  using Teuchos::rcp_dynamic_cast;
  using Thyra::VectorSpaceBase;
  using Thyra::assign;
  using Thyra::createMember;
  typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;

  //
  // Create and initialize product X, Xdot, Xdotdot

  RCP< const VectorSpaceBase<Scalar> > space;
  if (use_combined_method_)
    space = sens_model_->get_x_space();
  else
    space = sens_stepper_->get_x_space();
  RCP<DMVPV> X       = rcp_dynamic_cast<DMVPV>(createMember(space));
  RCP<DMVPV> Xdot    = rcp_dynamic_cast<DMVPV>(createMember(space));
  RCP<DMVPV> Xdotdot = rcp_dynamic_cast<DMVPV>(createMember(space));

  const int num_param = X->getNonconstMultiVector()->domain()->dim()-1;
  const Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
  const Teuchos::Range1D rng(1,num_param);

  // x
  assign(X->getNonconstMultiVector()->col(0).ptr(), *x0);
  if (DxDp0 == Teuchos::null)
    assign(X->getNonconstMultiVector()->subView(rng).ptr(), zero);
  else
    assign(X->getNonconstMultiVector()->subView(rng).ptr(), *DxDp0);

  // xdot
  if (xdot0 == Teuchos::null)
    assign(Xdot->getNonconstMultiVector()->col(0).ptr(), zero);
  else
    assign(Xdot->getNonconstMultiVector()->col(0).ptr(), *xdot0);
  if (DxdotDp0 == Teuchos::null)
    assign(Xdot->getNonconstMultiVector()->subView(rng).ptr(), zero);
  else
    assign(Xdot->getNonconstMultiVector()->subView(rng).ptr(), *DxdotDp0);

  // xdotdot
  if (xdotdot0 == Teuchos::null)
    assign(Xdotdot->getNonconstMultiVector()->col(0).ptr(), zero);
  else
    assign(Xdotdot->getNonconstMultiVector()->col(0).ptr(), *xdotdot0);
  if (DxdotDp0 == Teuchos::null)
    assign(Xdotdot->getNonconstMultiVector()->subView(rng).ptr(), zero);
  else
    assign(Xdotdot->getNonconstMultiVector()->subView(rng).ptr(), *DxdotdotDp0);

  integrator_->initializeSolutionHistory(t0, X, Xdot, Xdotdot);
}

template<class Scalar>
Teuchos::RCP<const Thyra::VectorBase<Scalar> >
IntegratorForwardSensitivity<Scalar>::
getX() const
{
  using Teuchos::RCP;
  using Teuchos::rcp_dynamic_cast;
  typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;

  RCP<const DMVPV> X = rcp_dynamic_cast<const DMVPV>(integrator_->getX());
  return X->getMultiVector()->col(0);
}

template<class Scalar>
Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> >
IntegratorForwardSensitivity<Scalar>::
getDxDp() const
{
  using Teuchos::RCP;
  using Teuchos::rcp_dynamic_cast;
  typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;

  RCP<const DMVPV> X = rcp_dynamic_cast<const DMVPV>(integrator_->getX());
  const int num_param = X->getMultiVector()->domain()->dim()-1;
  const Teuchos::Range1D rng(1,num_param);
  return X->getMultiVector()->subView(rng);
}

template<class Scalar>
Teuchos::RCP<const Thyra::VectorBase<Scalar> >
IntegratorForwardSensitivity<Scalar>::
getXDot() const
{
  using Teuchos::RCP;
  using Teuchos::rcp_dynamic_cast;
  typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;

  RCP<const DMVPV> Xdot = rcp_dynamic_cast<const DMVPV>(integrator_->getXDot());
  return Xdot->getMultiVector()->col(0);
}

template<class Scalar>
Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> >
IntegratorForwardSensitivity<Scalar>::
getDXDotDp() const
{
  using Teuchos::RCP;
  using Teuchos::rcp_dynamic_cast;
  typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;

  RCP<const DMVPV> Xdot = rcp_dynamic_cast<const DMVPV>(integrator_->getXDot());
  const int num_param = Xdot->getMultiVector()->domain()->dim()-1;
  const Teuchos::Range1D rng(1,num_param);
  return Xdot->getMultiVector()->subView(rng);
}

template<class Scalar>
Teuchos::RCP<const Thyra::VectorBase<Scalar> >
IntegratorForwardSensitivity<Scalar>::
getXDotDot() const
{
  using Teuchos::RCP;
  using Teuchos::rcp_dynamic_cast;
  typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;

  RCP<const DMVPV> Xdotdot =
    rcp_dynamic_cast<const DMVPV>(integrator_->getXDotDot());
  return Xdotdot->getMultiVector()->col(0);
}

template<class Scalar>
Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> >
IntegratorForwardSensitivity<Scalar>::
getDXDotDotDp() const
{
  using Teuchos::RCP;
  using Teuchos::rcp_dynamic_cast;
  typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;

  RCP<const DMVPV> Xdotdot =
    rcp_dynamic_cast<const DMVPV>(integrator_->getXDotDot());
  const int num_param = Xdotdot->getMultiVector()->domain()->dim()-1;
  const Teuchos::Range1D rng(1,num_param);
  return Xdotdot->getMultiVector()->subView(rng);
}

template<class Scalar>
Teuchos::RCP<const Thyra::VectorBase<Scalar> >
IntegratorForwardSensitivity<Scalar>::
getG() const
{
  typedef Thyra::ModelEvaluatorBase MEB;

  // Compute g which is computed by response 1 of the
  // sensitivity model evaluator
  Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > smodel;
  if (use_combined_method_)
    smodel = sens_model_;
  else
    smodel = sens_stepper_->getModel();
  MEB::InArgs<Scalar> inargs = smodel->getNominalValues();
  MEB::OutArgs<Scalar> outargs = smodel->createOutArgs();
  inargs.set_t(integrator_->getTime());
  inargs.set_x(integrator_->getX());
  if (inargs.supports(MEB::IN_ARG_x_dot))
    inargs.set_x_dot(integrator_->getXDot());
  if (inargs.supports(MEB::IN_ARG_x_dot_dot))
    inargs.set_x_dot_dot(integrator_->getXDotDot());

  Teuchos::RCP<Thyra::VectorBase<Scalar> > g =
    Thyra::createMember(smodel->get_g_space(1));
  outargs.set_g(1, g);

  smodel->evalModel(inargs, outargs);
  return g;
}

template<class Scalar>
Teuchos::RCP<const Thyra::MultiVectorBase<Scalar> >
IntegratorForwardSensitivity<Scalar>::
getDgDp() const
{
  typedef Thyra::ModelEvaluatorBase MEB;
  typedef Thyra::DefaultMultiVectorProductVector<Scalar> DMVPV;

  // Compute final dg/dp  which is computed by response 0  of the
  // sensitivity model evaluator
    Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > smodel;
  if (use_combined_method_)
    smodel = sens_model_;
  else
    smodel = sens_stepper_->getModel();
  MEB::InArgs<Scalar> inargs = smodel->getNominalValues();
  MEB::OutArgs<Scalar> outargs = smodel->createOutArgs();
  inargs.set_t(integrator_->getTime());
  inargs.set_x(integrator_->getX());
  if (inargs.supports(MEB::IN_ARG_x_dot))
    inargs.set_x_dot(integrator_->getXDot());
  if (inargs.supports(MEB::IN_ARG_x_dot_dot))
    inargs.set_x_dot_dot(integrator_->getXDotDot());

  Teuchos::RCP<Thyra::VectorBase<Scalar> > G =
    Thyra::createMember(smodel->get_g_space(0));
  Teuchos::RCP<DMVPV> dgdp = Teuchos::rcp_dynamic_cast<DMVPV>(G);
  outargs.set_g(0, G);

  smodel->evalModel(inargs, outargs);
  return dgdp->getMultiVector();
}

template<class Scalar>
std::string
IntegratorForwardSensitivity<Scalar>::
description() const
{
  std::string name = "Tempus::IntegratorForwardSensitivity";
  return(name);
}

template<class Scalar>
void
IntegratorForwardSensitivity<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;
  integrator_->describe(*l_out, verbLevel);
}

template <class Scalar>
SensitivityStepMode
IntegratorForwardSensitivity<Scalar>::
getStepMode() const
{
  if (use_combined_method_)
    return SensitivityStepMode::Combined;
  return sens_stepper_->getStepMode(); // Staggered case
}

template<class Scalar>
Teuchos::RCP<IntegratorForwardSensitivity<Scalar> >
createIntegratorForwardSensitivity(
  Teuchos::RCP<Teuchos::ParameterList>                     pList,
  const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >&      model,
  const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >&      sens_residual_model,
  const Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >&      sens_solve_model)
{


  Teuchos::RCP<SensitivityModelEvaluatorBase<Scalar> > sens_model;
  Teuchos::RCP<StepperStaggeredForwardSensitivity<Scalar>> sens_stepper;

  //1. create integrator
  auto fwd_integrator = createIntegratorBasic<Scalar>(pList, model);

  //2. set parameter list
  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::rcp(new Teuchos::ParameterList);
  {
    Teuchos::RCP<const Teuchos::ParameterList> integrator_pl = fwd_integrator->getValidParameters();
    Teuchos::RCP<const Teuchos::ParameterList> sensitivity_pl =
      CombinedForwardSensitivityModelEvaluator<Scalar>::getValidParameters();
    pl->setParameters(*integrator_pl);
    Teuchos::ParameterList &spl = pl->sublist("Sensitivities");
    spl.setParameters(*sensitivity_pl);
    spl.set("Sensitivity Method", "Combined");
    spl.set("Reuse State Linear Solver", false);
  }
  pList->setParametersNotAlreadySet(*pl);

  auto sens_pl                    = Teuchos::sublist(pList, "Sensitivities", false);
  std::string integratorName       = pList->get<std::string>("Integrator Name", "Default Integrator");
  std::string stepperName          = pList->sublist(integratorName).get<std::string>("Stepper Name");
  auto stepper_pl                 = Teuchos::sublist(pList, stepperName, true);
  std::string sensitivity_method   = sens_pl->get<std::string>("Sensitivity Method");
  bool use_combined_method        = sensitivity_method == "Combined";

  //3. create sensitivity model and stepper
  // createSensitivityModelAndStepper
  {
    sens_pl->remove("Sensitivity Method");

    if (use_combined_method)
    {
      sens_pl->remove("Reuse State Linear Solver");
      sens_model     = wrapCombinedFSAModelEvaluator(
        model, sens_residual_model, sens_solve_model, sens_pl);
      fwd_integrator = createIntegratorBasic<Scalar>(pList, sens_model);
    }
    else
    {
      sens_stepper          = Teuchos::rcp(new StepperStaggeredForwardSensitivity<Scalar>(model, sens_residual_model, sens_solve_model, stepper_pl, sens_pl));
      auto fsa_staggered_me = Teuchos::rcp_const_cast<Thyra::ModelEvaluator<Scalar>>(sens_stepper->getModel());
      fwd_integrator        = createIntegratorBasic<Scalar>(pList, fsa_staggered_me);
      fwd_integrator->setStepper(sens_stepper);
    }
  }

  //4. initialize propoer integrator
  Teuchos::RCP<IntegratorForwardSensitivity<Scalar>> integrator =
    Teuchos::rcp(new IntegratorForwardSensitivity<Scalar>(model, fwd_integrator, sens_model, sens_stepper, use_combined_method));
  return (integrator);
}

template<class Scalar>
Teuchos::RCP<IntegratorForwardSensitivity<Scalar> >
createIntegratorForwardSensitivity()
{

  Teuchos::RCP<IntegratorForwardSensitivity<Scalar> > integrator =
    Teuchos::rcp(new IntegratorForwardSensitivity<Scalar>());
  return(integrator);
}

} // namespace Tempus
#endif // Tempus_IntegratorForwardSensitivity_impl_hpp