Tempus_IntegratorAdjointSensitivity_decl.hpp

Go to the documentation of this file.

//@HEADER
// *****************************************************************************
//          Tempus: Time Integration and Sensitivity Analysis Package
//
// Copyright 2017 NTESS and the Tempus contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
//@HEADER
 
#ifndef Tempus_IntegratorAdjointSensitivity_decl_hpp
#define Tempus_IntegratorAdjointSensitivity_decl_hpp
 
#include <vector>
 
// Tempus
#include "Tempus_config.hpp"
#include "Tempus_IntegratorBasic.hpp"
#include "Tempus_AdjointAuxSensitivityModelEvaluator.hpp"
#include "Tempus_AdjointSensitivityUtils.hpp"
#include "Tempus_AdjointDAEJumpSolver.hpp"
#include "Tempus_DAEStructure.hpp"
 
#include "Tempus_StepperStaggeredForwardSensitivity.hpp"  // For SensitivityStepMode
 
namespace Tempus {

template <class Scalar>
class IntegratorAdjointSensitivity : virtual public Tempus::Integrator<Scalar> {
 public:

  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>> &solution_history,
      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,
      const AdjointResponseType response_type   = AdjointResponseType::FinalTime,
      const std::vector<Scalar> &response_times = std::vector<Scalar>(),
      const Teuchos::RCP<const DAEStructure<Scalar>> &dae_structure =
          Teuchos::null,
      const Teuchos::RCP<const Teuchos::ParameterList> &dae_jump_solver_pl =
          Teuchos::null);

  IntegratorAdjointSensitivity();

  virtual ~IntegratorAdjointSensitivity() {}



  virtual bool advanceTime();
  virtual bool advanceTime(const Scalar timeFinal) override;
  virtual Scalar getTime() const override;
  virtual int getIndex() const override;
  virtual Status getStatus() const override;
  virtual void setStatus(const Status st) override;
  virtual Teuchos::RCP<Stepper<Scalar>> getStepper() const override;
  virtual Teuchos::RCP<const SolutionHistory<Scalar>> getSolutionHistory()
      const override;
  Teuchos::RCP<const SolutionHistory<Scalar>> getStateSolutionHistory() const;
  Teuchos::RCP<const SolutionHistory<Scalar>> getSensSolutionHistory() const;
  virtual Teuchos::RCP<SolutionHistory<Scalar>> getNonConstSolutionHistory()
      override;
  virtual Teuchos::RCP<const TimeStepControl<Scalar>> getTimeStepControl()
      const override;
  virtual Teuchos::RCP<TimeStepControl<Scalar>> getNonConstTimeStepControl()
      override;
  Teuchos::RCP<TimeStepControl<Scalar>> getStateNonConstTimeStepControl();
  Teuchos::RCP<TimeStepControl<Scalar>> getSensNonConstTimeStepControl();
  virtual Teuchos::RCP<Teuchos::Time> getIntegratorTimer() const override
  {
    return state_integrator_->getIntegratorTimer();
  }
  virtual Teuchos::RCP<Teuchos::Time> getStepperTimer() const override
  {
    return state_integrator_->getStepperTimer();
  }


  virtual void initializeSolutionHistory(
      Scalar t0, Teuchos::RCP<const Thyra::VectorBase<Scalar>> x0,
      Teuchos::RCP<const Thyra::VectorBase<Scalar>> xdot0      = Teuchos::null,
      Teuchos::RCP<const Thyra::VectorBase<Scalar>> xdotdot0   = Teuchos::null,
      Teuchos::RCP<const Thyra::MultiVectorBase<Scalar>> DxDp0 = Teuchos::null,
      Teuchos::RCP<const Thyra::MultiVectorBase<Scalar>> DxdotDp0 =
          Teuchos::null,
      Teuchos::RCP<const Thyra::MultiVectorBase<Scalar>> DxdotdotDp0 =
          Teuchos::null);

  virtual Teuchos::RCP<IntegratorObserver<Scalar>> getObserver();
  virtual void setObserver(
      Teuchos::RCP<IntegratorObserver<Scalar>> obs = Teuchos::null);
  virtual void initialize();

  virtual Teuchos::RCP<const Thyra::VectorBase<Scalar>> getX() const;
  virtual Teuchos::RCP<const Thyra::VectorBase<Scalar>> getXDot() const;
  virtual Teuchos::RCP<const Thyra::VectorBase<Scalar>> getXDotDot() const;

  virtual Teuchos::RCP<const Thyra::MultiVectorBase<Scalar>> getY() const;
  virtual Teuchos::RCP<const Thyra::MultiVectorBase<Scalar>> getYDot() const;
  virtual Teuchos::RCP<const Thyra::MultiVectorBase<Scalar>> getYDotDot() const;

  virtual Teuchos::RCP<const Thyra::MultiVectorBase<Scalar>> getDgDp() const;
 
  Teuchos::RCP<Thyra::ModelEvaluator<Scalar>> getAdjointModel() const
  {
    return adjoint_model_;
  };

  Teuchos::RCP<const DAEStructure<Scalar>> getDAEStructure() const
  {
    return dae_structure_;
  }


  std::string description() const override;
  void describe(Teuchos::FancyOStream &out,
                const Teuchos::EVerbosityLevel verbLevel) const override;

  SensitivityStepMode getStepMode() const;
 
 protected:
  // Create sensitivity model evaluator from application model
  Teuchos::RCP<AdjointAuxSensitivityModelEvaluator<Scalar>> createAdjointModel(
      const Teuchos::RCP<Thyra::ModelEvaluator<Scalar>> &model,
      const Teuchos::RCP<Thyra::ModelEvaluator<Scalar>> &adjoint_model,
      const Teuchos::RCP<Teuchos::ParameterList> &inputPL);
 
  void buildSolutionHistory(
      const Teuchos::RCP<const SolutionHistory<Scalar>> &state_solution_history,
      const Teuchos::RCP<const SolutionHistory<Scalar>>
          &adjoint_solution_history);

  Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar>>
  createAdjointMassMatrixSolver(
      const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs) const;

  Teuchos::RCP<const SolutionState<Scalar>> findForwardState(
      const SolutionHistory<Scalar> &state_solution_history, const Scalar t,
      const Scalar tol) const;

  bool integrateAdjoint(const Scalar tau_begin, const Scalar tau_end,
                        Thyra::VectorBase<Scalar> &adjoint_state, int &index);

  void appendAdjointSolutionHistory(const Scalar tol);
 
  Teuchos::RCP<Thyra::ModelEvaluator<Scalar>> model_;
  Teuchos::RCP<IntegratorBasic<Scalar>> state_integrator_;
  Teuchos::RCP<Thyra::ModelEvaluator<Scalar>> adjoint_model_;
  Teuchos::RCP<AdjointAuxSensitivityModelEvaluator<Scalar>> adjoint_aux_model_;
  Teuchos::RCP<IntegratorBasic<Scalar>> adjoint_integrator_;
  Teuchos::RCP<SolutionHistory<Scalar>> solutionHistory_;
  int p_index_;
  int g_index_;
  bool g_depends_on_p_;
  bool f_depends_on_p_;
  bool ic_depends_on_p_;
  bool mass_matrix_is_identity_;
  AdjointResponseType response_type_;
  std::vector<Scalar> response_times_;
  Teuchos::RCP<SolutionHistory<Scalar>> adjoint_solution_history_;
  Teuchos::RCP<const Thyra::MultiVectorBase<Scalar>> dxdp_init_;
  Teuchos::RCP<Thyra::MultiVectorBase<Scalar>> dgdp_;
  SensitivityStepMode stepMode_;
  Teuchos::RCP<const DAEStructure<Scalar>> dae_structure_;
  Teuchos::RCP<const AdjointDAEJumpSolver<Scalar>> dae_jump_solver_;
  Teuchos::RCP<Thyra::LinearOpBase<Scalar>> dfdp_op_;
  Teuchos::RCP<Thyra::MultiVectorBase<Scalar>> dfdp_mv_;
};


template <class Scalar>
Teuchos::RCP<IntegratorAdjointSensitivity<Scalar>>
createIntegratorAdjointSensitivity(
    Teuchos::RCP<Teuchos::ParameterList> pList,
    const Teuchos::RCP<Thyra::ModelEvaluator<Scalar>> &model,
    const Teuchos::RCP<Thyra::ModelEvaluator<Scalar>> &adjoint_model =
        Teuchos::null);

template <class Scalar>
Teuchos::RCP<IntegratorAdjointSensitivity<Scalar>>
createIntegratorAdjointSensitivity();
 
}  // namespace Tempus
 
#endif  // Tempus_IntegratorAdjointSensitivity_decl_hpp