Tempus: Add combined, staggered, and pseudotransient FSA methods

This commit adds integrators supporting the combined, staggered, and
pseudotransient forward sensitivity analysis methods where the
sensitivity equations are solved alongside the forward equations.  The
different methods refer to how the solves between states/sensitivities
are partitioned:

  * combined:  solved together (best for explicit)
  * staggered: solve state equations and then sensitivity equations for
    each time step (best for implicit)
  * pseudotransient: do full state equation integration followed by full
    sensitivity equation integration (for when time integrator is used
    as a steady-state solver)

In all cases there is an option, "Use DfDp as Tangent", that controls
whether the DfDp out-arg is interpreted as the tangent df/dp +
df/dx*dx/dp + df/dxdot*dxdot/dp, or not (which is not regular model
evaluator behavior).  If the model supports it, it is a much simpler and
more efficient way to assemble the sensitivity equations.

It also supports options to reuse the linear solver from the state
equations when solving the sensitivity equations.  This saves Jacobian
assembly time and preconditioner generation time.  Supporting it
required one small change to the NOX-Thyra interface to not throw when
accessing the preconditioner RCP if it is null (because there is no way
to know beforehand if it is null or not until you access it).

Tests were added for Backward Euler, Explicit RK, and BDF2 time
integration methods for all of the above methods and options.

packages/nox/src-thyra/NOX_Thyra_Group.C

@@ -380,15 +380,13 @@ NOX::Thyra::Group::getJacobian() const
 Teuchos::RCP< ::Thyra::PreconditionerBase<double> >
 NOX::Thyra::Group::getNonconstPreconditioner()
 {
-  NOX_ASSERT(nonnull(prec_));
   return prec_;
 }
 
 
 Teuchos::RCP<const ::Thyra::PreconditionerBase<double> >
 NOX::Thyra::Group::getPreconditioner() const
 {
-  NOX_ASSERT(nonnull(prec_));
   return prec_;
 }
 

packages/tempus/src/CMakeLists.txt

@@ -10,7 +10,7 @@ SET_AND_INC_DIRS(DIR ${CMAKE_CURRENT_SOURCE_DIR})
 APPEND_GLOB(HEADERS ${DIR}/*.hpp)
 APPEND_GLOB(SOURCES ${DIR}/*.cpp)
 
-# Must glob the binary dir last to get all of the 
+# Must glob the binary dir last to get all of the
 # auto-generated ETI headers
 TRILINOS_CREATE_CLIENT_TEMPLATE_HEADERS(${DIR})
 SET_AND_INC_DIRS(DIR ${CMAKE_CURRENT_BINARY_DIR})
@@ -21,3 +21,4 @@ TRIBITS_ADD_LIBRARY(
   HEADERS ${HEADERS}
   SOURCES ${SOURCES}
   )
+

packages/tempus/src/Tempus_CombinedForwardSensitivityModelEvaluator.cpp

@@ -0,0 +1,19 @@
+// @HEADER
+// ****************************************************************************
+//                Tempus: Copyright (2017) Sandia Corporation
+//
+// Distributed under BSD 3-clause license (See accompanying file Copyright.txt)
+// ****************************************************************************
+// @HEADER
+
+#include "Tempus_ExplicitTemplateInstantiation.hpp"
+
+#ifdef HAVE_TEMPUS_EXPLICIT_INSTANTIATION
+#include "Tempus_CombinedForwardSensitivityModelEvaluator_decl.hpp"
+#include "Tempus_CombinedForwardSensitivityModelEvaluator_impl.hpp"
+
+namespace Tempus {
+  TEMPUS_INSTANTIATE_TEMPLATE_CLASS(CombinedForwardSensitivityModelEvaluator)
+}
+
+#endif

packages/tempus/src/Tempus_CombinedForwardSensitivityModelEvaluator_decl.hpp

@@ -0,0 +1,136 @@
+// @HEADER
+// ****************************************************************************
+//                Tempus: Copyright (2017) Sandia Corporation
+//
+// Distributed under BSD 3-clause license (See accompanying file Copyright.txt)
+// ****************************************************************************
+// @HEADER
+
+#ifndef Tempus_CombinedForwardSensitivityModelEvaluator_decl_hpp
+#define Tempus_CombinedForwardSensitivityModelEvaluator_decl_hpp
+
+#include "Thyra_StateFuncModelEvaluatorBase.hpp"
+#include "Thyra_DefaultMultiVectorProductVectorSpace.hpp"
+
+namespace Tempus {
+
+/** \brief Transform a ModelEvaluator's sensitivity equations to its residual
+ *
+ * This class wraps a given ModelEvalutor encapsulating f(x,p) and creates
+ * a new "residual" for it and the forward sensitivity equations:
+ *        F(X) = [                f(x,p)             ] = 0
+ *               [ (df/dx)(x,p) * dx/dp + df/dp(x,p) ]
+ * where X = [ x; dx/dp ] (transient terms supressed for simplicity). This model
+ * evaluator can then be handed to a regular (non)linear solver to compute X.
+ * Note that even though these equations are linear in X, it is not necessarily
+ * the case that the underlying model evaluator accurately computes df/dx in its
+ * evaluation of W.  Therefore this model evaluator can optionally reinterpret
+ * the model's df/dp out-arg as (df/dx)(x,p) * dx/dp + df/dp(x,p) where dx/dp is
+ * passed as another parameter (product) vector (encapsulated in the
+ * Thyra::DefaultMultiVectorProductVector).  This is not standard model
+ * evaluator behavior, but is useful for models where W is only an approximation
+ * to df/dx and/or the model is capable of directly computing
+ * (df/dx)(x,p) * dx/dp + df/dp(x,p).
+ */
+template <typename Scalar>
+class CombinedForwardSensitivityModelEvaluator :
+    public Thyra::StateFuncModelEvaluatorBase<Scalar> {
+public:
+  typedef Thyra::VectorBase<Scalar>  Vector;
+  typedef Thyra::MultiVectorBase<Scalar>  MultiVector;
+
+  //! Constructor
+  /*!
+   * The optionally supplied parameter list supports the following options:
+   * <ul>
+   *   <li> "Use DfDp as Tangent" (default:  false) Reinterpret the df/dp
+   *        out-arg as the tangent vector (df/dx)(x,p) * dx/dp + df/dp(x,p)
+   *        as described above.  If it is false, this implementation will
+   *        compute the tangent through several calls to the models
+   *        evalModel().
+   *   <li> "Sensitivity Parameter Index" (default: 0) Model evaluator
+   *        parameter index for which sensitivities will be computed.
+   *   <li> "Sensitivity X Tangent Index" (default: 1) If "Use DfDp as Tangent"
+   *        is true, the model evaluator parameter index for passing dx/dp
+   *        as a Thyra::DefaultMultiVectorProductVector.
+   *   <li> "Sensitivity X-Dot Tangent Index" (default: 2) If
+   *        "Use DfDp as Tangent" is true, the model evaluator parameter index
+   *        for passing dx_dot/dp as a Thyra::DefaultMultiVectorProductVector.
+   *   <li> "Sensitivity X-Dot-Dot Tangent Index" (default: 3) If
+   *        "Use DfDp as Tangent" is true, the model evaluator parameter index
+   *        for passing dx_dot_dot/dp as a
+   *        Thyra::DefaultMultiVectorProductVector (if the model supports
+   *        x_dot_dot).
+   * </ul>
+   */
+  CombinedForwardSensitivityModelEvaluator(
+    const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > & model,
+    const Teuchos::RCP<const Teuchos::ParameterList>& pList = Teuchos::null,
+    const Teuchos::RCP<MultiVector>& dxdp_init = Teuchos::null,
+    const Teuchos::RCP<MultiVector>& dx_dotdp_init = Teuchos::null,
+    const Teuchos::RCP<MultiVector>& dx_dotdot_dp_init = Teuchos::null);
+
+  //! Get the underlying model 'f'
+  Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > getModel() const
+  { return model_; }
+
+  /** \name Public functions overridden from ModelEvaulator. */
+  //@{
+
+  Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > get_p_space(int p) const;
+
+  Teuchos::RCP<const Teuchos::Array<std::string> > get_p_names(int p) const;
+
+  Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > get_x_space() const;
+
+  Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> > get_f_space() const;
+
+  Teuchos::RCP<Thyra::LinearOpBase<Scalar> > create_W_op() const;
+
+  Teuchos::RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar> >
+  get_W_factory() const;
+
+  Thyra::ModelEvaluatorBase::InArgs<Scalar> createInArgs() const;
+
+  Thyra::ModelEvaluatorBase::InArgs<Scalar> getNominalValues() const;
+
+  //@}
+
+  static Teuchos::RCP<const Teuchos::ParameterList> getValidParameters();
+
+private:
+
+  Thyra::ModelEvaluatorBase::OutArgs<Scalar> createOutArgsImpl() const;
+
+  void evalModelImpl(
+    const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
+    const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs) const;
+
+
+  Thyra::ModelEvaluatorBase::InArgs<Scalar> prototypeInArgs_;
+  Thyra::ModelEvaluatorBase::OutArgs<Scalar> prototypeOutArgs_;
+
+  Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> > model_;
+  Teuchos::RCP<MultiVector> dxdp_init_;
+  Teuchos::RCP<MultiVector> dx_dotdp_init_;
+  Teuchos::RCP<MultiVector> dx_dotdotdp_init_;
+  int p_index_;
+  int x_tangent_index_;
+  int xdot_tangent_index_;
+  int xdotdot_tangent_index_;
+  bool use_dfdp_as_tangent_;
+
+  int num_param_;
+  Teuchos::RCP<const Thyra::DefaultMultiVectorProductVectorSpace<Scalar> > dxdp_space_;
+  Teuchos::RCP<const Thyra::DefaultMultiVectorProductVectorSpace<Scalar> > x_dxdp_space_;
+  Teuchos::RCP<const Thyra::DefaultMultiVectorProductVectorSpace<Scalar> > dfdp_space_;
+  Teuchos::RCP<const Thyra::DefaultMultiVectorProductVectorSpace<Scalar> > f_dfdp_space_;
+
+  mutable Teuchos::RCP<Thyra::LinearOpBase<Scalar> > my_dfdx_;
+  mutable Teuchos::RCP<Thyra::LinearOpBase<Scalar> > my_dfdxdot_;
+  mutable Teuchos::RCP<Thyra::LinearOpBase<Scalar> > my_dfdxdotdot_;
+};
+
+} // namespace Tempus
+
+#endif